On Sun, Jul 24, 2022 at 8:46 AM Duan Junming via petsc-users <
petsc-users@mcs.anl.gov> wrote:
> Dear all,
>
>
> I want to add more dofs to the coordinate DM to represent a curved mesh.
>
> I first create a simple box mesh with one cell:
>
>
> PetscCall(DMCreate(comm, &dm));
> PetscCall(DMSetType(dm, DMPLEX));
> dim = 2;
> PetscInt n_faces[2] = {1, 1};
> DMBoundaryType periodicity[2] = {DM_BOUNDARY_GHOSTED,
> DM_BOUNDARY_GHOSTED};
> DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, n_faces, NULL, NULL,
> periodicity, PETSC_TRUE, &dm);
>
> Then I add more dofs to cdm using section:
>
> PetscCall(DMGetCoordinateDM(dm, &cdm));
> PetscCall(DMGetLocalSection(cdm, &cs));
> PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
> PetscCall(PetscSectionSetChart(cs, pStart, pEnd));
> degree = 3;
> for(depth = 0; depth <= dim; ++ depth) {
> PetscCall(DMPlexGetDepthStratum(cdm, depth, &pStart, &pEnd));
> for(p = pStart; p < pEnd; ++p) {
> *PetscCall(PetscSectionSetDof(cs, p, dim*PetscPowInt(degree-1,
> depth)));*
> * PetscCall(PetscSectionSetFieldDof(cs, p, 0,
> dim*PetscPowInt(degree-1, depth)));*
> }
> }
> PetscCall(PetscSectionSetUp(cs));
> PetscCall(DMSetUp(cdm));
>
> Finally, I wish to get the vec to set the coordinates:
>
> PetscCall(DMCreateLocalVector(cdm, &coordinates));
> PetscScalar *pCoords = NULL;
> PetscCall(DMPlexVecGetClosure(cdm, NULL, coordinates, 0, &pSize,
> &pCoords));
> printf("pSize: %d\n", pSize);
>
> However, if I only use *PetscSectionSetDof*, the code reports "Section
> size 32 does not match Vec closure size 0".
> Instead, if I only use * PetscSectionSetFieldDof*, the code crashes with
> "Checking
> the memory for corruption."
> It runs without error when I use both * PetscSectionSetDof *and *
> PetscSectionSetFieldDof*.
>
> Maybe I have missed something, but could you help to explain why I should
> both functions *PetscSectionSetDof* and * PetscSectionSetFieldDof?*
> It is not straightforward that I should again set the dofs of the section
> after each field of the section has been set.
>
This is the intended behavior, and how we use it internally. The use of
fields is optional in order to allow Section to be as lightweight
as possible. We do not enforce consistency between the field and overall
sizes.
Thanks,
Matt
> I can also see that *PetscSectionSetFieldDof *call *PetscSectionSetDof *on
> the specific field of the section,
> and cs->field[0] doesn't share the same address as cs itself.
>
> Thanks in advance!
> Junming
>
>
>
--
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/>
