[petsc-users] Adding more dofs to the coordinate DM crashes with either PetscSectionSetDof or PetscSectionFieldDof

[Duan Junming]

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.

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


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220724/7e4dbb12/attachment.html>