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

[Duan Junming]

Thank you very much!

________________________________
From: knepley at gmail.com <knepley at gmail.com>
Sent: Sunday, July 24, 2022 5:32:34 PM
To: Duan Junming
Cc: PETSc
Subject: Re: [petsc-users] Adding more dofs to the coordinate DM crashes with either PetscSectionSetDof or PetscSectionFieldDof

On Sun, Jul 24, 2022 at 8:46 AM Duan Junming via petsc-users <petsc-users at mcs.anl.gov<mailto:petsc-users at 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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220725/f04c45bc/attachment.html>