Hi matt,
Many thanks for the help !!
Le 3/26/2020 à 11:38 PM, Matthew Knepley a écrit :
17 18
7-----8-----9-----14----15
| | |
16,3 4 5,6 11 12,13
| | |
1-----2-----3-----9-----10
19 20
so each face gets 2 dofs, one for each cell.When doing a cell integral,
you only use the dof that is for that cell.
I think i understood. I put the beginning of the code for dof management
in attachment. i defined :
numDof[0*(user.dim+1)+user.dim] = internalDof; /* internalDof
defined on cells */
numDof[0*(user.dim+1)+user.dim-1] = nbDofFaceEle; /* nbDofFaceEle
defined on faces */
This way, the dof faces are not duplicated. They are defined on the
face. (this part is commented in the source file)
Am i right ?
> The local-to-global would update the face dofs, so you would get each
side.
>
This way, a local-to-global will update the face dofs ?
There is a reordering when you extract the closure. I have written one
for spectral elements. We would need
another here that ordered all the "other" face dofs to the end.
You mean using PetscSectionSetPermutation ?
Or rewrite DMPlexSetClosurePermutationTensor in order to put all the
face dofs at the end ?
I don't understand this part. I'm looking at the documentation for
DMPlexSetClosurePermutationTensor.
https://www.mcs.anl.gov/petsc/petsc-current/src/dm/impls/plex/plex.c.html#DMPlexSetClosurePermutationTensor
It says that
----------------------
The closure in BFS ordering works through height strata (cells, edges,
vertices) to produce the ordering
.vb
0 1 2 3 8 9 14 15 11 10 13 12 4 5 7 6
.ve
--------------------------
which is what we want no ?
This seems a little complicated to me. Do you know how Andreas Klockner
does it in Hedge? Or Tim Warburton?
I'm reading his book (nodal discontinuous galerkin methods)
I just want to make sure I am not missing an elegant way to handle this.
Here field 1 is synchronised with field 0, locally.
But the external Face dof, field 2, have to be synchronised with the
values of the adjacent cell.
Is it possible to use something like DMPlexGetFaceFields ?
Is there an example of such use of PetscSection and synchronisation
process ?
For the parallel part, should i use PetscSF object ?
In parallel, integrals would be summed into the global vector, so each
side has a 0 for the other face dof and the right contribution
for its face dof. Then both sides get both solution dofs. It seems to
work in my head.
I read your article "Mesh Algorithms for PDE with Sieve I: Mesh
Distribution". But it's refereeing to Matthew G. Knepley and Dmitry A.
Karpeev. Sieve implementation.
Technical Report ANL/MCS to appear, Argonne National Laboratory,
January 2008.
I couldn't find it. It is freely available ?
Don't bother reading that. There are later ones:
There are two pretty good sources:
https://arxiv.org/abs/1505.04633
https://arxiv.org/abs/1506.06194
The last one is a follow-on to this paper
https://arxiv.org/abs/0908.4427
Thanks,
Matt
>
> I don't think you need extra vertices, > or coordinates, and for
output I
> recommend using DMPlexProject() to get
> the solution in some space that can be plotted like P1, or
anything else
> supported by your visualization.
I would like to use DMplex as much as i can, as i would in the future
refine locally the mesh.
I hope you're good in this difficult situation (covid19),
Best regards,
Yann
>
> Thanks,
>
> Matt
>
> >
> > We have an implementation of spectral element ordering
> >
>
(https://gitlab.com/petsc/petsc/-/blob/master/src/dm/impls/plex/examples/tutorials/ex6.c).
>
> > Those share
> > the whole element boundary.
> >
> > 2) What ghosts do you need?
> In order to compute the numerical fluxes of one element, i
need the
> values of the surrounding nodes connected to the adjacent
elements.
> >
> > 3) You want to store real space coordinates for a
quadrature?
> It should be basically the same as PetscFE of higher order.
> I add some vertex needed to compute a polynomal solution of
the desired
> order. That means that if i have a N, order of the local
approximation,
> i need 0.5*(N+1)*(N+2) vertex to store in the DMPlex (in 2D), in
> order to :
> 1) have the correct number of dof
> 2) use ghost nodes to sync the values of the
vertex/edge/facet for
> 1D/2D/3D problem
> 2) save correctly the solution
>
> Does it make sense to you ?
>
> Maybe like
>
https://www.mcs.anl.gov/petsc/petsc-current/src/ts/examples/tutorials/ex11.c.html
> With the use of the function SplitFaces, which i didn't fully
> understood
> so far.
>
> Thanks,
>
> Yann
>
> >
> > We usually define a quadrature on the reference
element once.
> >
> > Thanks,
> >
> > Matt
> >
> > I found elements of answers in those threads :
> >
>
https://lists.mcs.anl.gov/pipermail/petsc-users/2016-August/029985.html
> >
>
https://lists.mcs.anl.gov/mailman/htdig/petsc-users/2019-October/039581.html
> >
> > However, it's not clear for me where to begin.
> >
> > Quoting Matt, i should :
> > " DMGetCoordinateDM(dm, &cdm);
> > <Set field information into cdm instead of dm>
> > DMCreateLocalVector(cdm, &coordinatesLocal);
> > <Fill in higher order coordinate values>
> > DMSetCoordinatesLocal(dm, coordinatesLocal);"
> >
> > However, i will not create ghost nodes this way. And
i'm not
> sure to
> > keep the good ordering.
> > This part should be implemented in the PetscFE
interface, for
> high
> > order
> > discrete solutions.
> > I did not succeed in finding the correct part of the
source
> doing it.
> >
> > Could you please give me some hint to begin correctly
thoses
> tasks ?
> >
> > Thanks,
> >
> > Yann
> >
> >
> >
> > --
> > What most experimenters take for granted before they begin
their
> > experiments is infinitely more interesting than any
results to which
> > their experiments lead.
> > -- Norbert Wiener
> >
> > https://www.cse.buffalo.edu/~knepley/
> <http://www.cse.buffalo.edu/~knepley/>
>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which
> their experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
<http://www.cse.buffalo.edu/~knepley/>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which
their experiments lead.
-- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
static char help[] = "dof management with dmplex for DG problems, for simplices
\n\n";
#include <petscdmplex.h>
typedef struct {
PetscInt dim; /* Topological problem dimension */
PetscInt Nf; /* Number of fields */
PetscInt Nc[1]; /* Number of components per field */
PetscInt N; /* Order of polymomials used for approximation */
} AppCtx;
static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
{
PetscBool flg;
PetscErrorCode ierr;
PetscFunctionBeginUser;
options->dim = 2;
options->Nf = 1;
options->Nc[0] = 1;
options->N = 4;
ierr = PetscOptionsBegin(comm, "", "dof management with dmplex for DG
problems (simplices) Options", "DMPLEX");CHKERRQ(ierr);
ierr = PetscOptionsRangeInt("-dim", "Problem dimension", "dof_dg.c",
options->dim, &options->dim, NULL,1,3);CHKERRQ(ierr);
ierr = PetscOptionsInt("-order", "Order of polymomials used for
approximation", "dof_dg.c", options->N, &options->N, &flg);CHKERRQ(ierr);
ierr = PetscOptionsEnd();
PetscFunctionReturn(0);
}
int main(int argc, char **argv) {
DM dm;
PetscSection s;
Vec u;
PetscViewer viewer;
AppCtx user;
PetscInt cells[3] = {2, 2, 2};
PetscReal lower[3] = {-1,-1,-1};
PetscReal upper[3] = {1,1,1};
/* DG specification of the PetscSection */
/* order of poly approx is stored in AppCtx user : user.N */
PetscInt dofPerEle=0; /* Number of dof per cell */
PetscInt internalDof; /* dof defined on the internal cell */
/* so without the faces dof */
PetscInt nbFace=0; /* Number of faces per cell */
PetscInt nbDofFaceEle=0; /* Number of dof per face */
PetscErrorCode ierr;
ierr = PetscInitialize(&argc, &argv, NULL, help); if (ierr) return ierr;
ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, user.dim, PETSC_FALSE, cells,
lower, upper, NULL, PETSC_TRUE, &dm);CHKERRQ(ierr);
ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
/* Total number of dof per element, number of faces per element and number of
dof per face */
if (user.dim == 2) {
dofPerEle = (user.N+1)*(user.N+2)/2;
nbFace = 3;
nbDofFaceEle = user.N + 1;
internalDof = dofPerEle - nbFace*nbDofFaceEle + 3;
} else if (user.dim == 3) {
dofPerEle = (user.N+1)*(user.N+2)*(user.N+3)/6;
nbFace = 4;
nbDofFaceEle = (user.N+1)*(user.N+2)/2;
internalDof = dofPerEle - nbFace*nbDofFaceEle + (user.N + 1)*6 - 4;
}
PetscPrintf(MPI_COMM_WORLD,"\ndof per cell : %d\n",dofPerEle);
PetscPrintf(MPI_COMM_WORLD,"dof per face : %d\n",nbDofFaceEle);
PetscPrintf(MPI_COMM_WORLD,"internal dof : %d\n",internalDof);
{
PetscInt *numDof, d;
/* numDof[f*(dim+1)+d] gives the number of dof for field f on */
/* points of dimension d. For instance, numDof[1] is the number */
/* of dof for field 0 on each edge. */
ierr = PetscMalloc1(user.Nf*(user.dim+1), &numDof);CHKERRQ(ierr);
/* fields 0, written generally to add easily more fields */
for (d = 0; d < user.Nf*(user.dim+1); ++d) numDof[d] = 0;
/* we put only the interior dof for not duplicating the face nodes */
numDof[0*(user.dim+1)+user.dim] = internalDof; /* internalDof defined
on cells */
numDof[0*(user.dim+1)+user.dim-1] = nbDofFaceEle; /* nbDofFaceEle defined
on faces */
ierr = DMSetNumFields(dm, user.Nf);CHKERRQ(ierr);
ierr = DMPlexCreateSection(dm, NULL, user.Nc, numDof, 0, NULL, NULL, NULL,
NULL, &s);CHKERRQ(ierr);
ierr = PetscFree(numDof);CHKERRQ(ierr);
}
/* Name the Field variable */
ierr = PetscSectionSetFieldName(s, 0, "u");CHKERRQ(ierr);
/* Tell the DM to use this data layout */
ierr = DMSetLocalSection(dm, s);CHKERRQ(ierr);
/* Create a Vec with this layout and view it */
ierr = DMGetGlobalVector(dm, &u);CHKERRQ(ierr);
ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer);CHKERRQ(ierr);
ierr = PetscViewerSetType(viewer, PETSCVIEWERVTK);CHKERRQ(ierr);
ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);
ierr = PetscViewerFileSetName(viewer, "sol.vtk");CHKERRQ(ierr);
ierr = VecView(u, viewer);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
ierr = DMRestoreGlobalVector(dm, &u);CHKERRQ(ierr);
/* Cleanup */
ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);
ierr = DMDestroy(&dm);CHKERRQ(ierr);
ierr = PetscFinalize();
return ierr;
}
