Thank you very much! ________________________________ From: [email protected] <[email protected]> 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 <[email protected]<mailto:[email protected]>> 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/>
