[petsc-users] Inquiry about resetting a petscsection for a dmplex

[neil liu]

Dear Pestc community,

I am currently working on a 3D adaptive vector FEM solver. In my case, I
need to solve two systems: one for the primal equation using a low-order
discretization and another for the adjoint equation using a high-order
discretization.

Afterward, I need to reset the section associated with the DMPlex.
Whichever is set first—20 DOFs (second-order) or 6 DOFs (first-order)—the
final mapping always follows that of the first-defined configuration.

Did I miss something?

Thanks,


Xiaodong

PetscErrorCode DMManage::SetupSection(CaseInfo &objCaseInfo){
PetscSection s;
PetscInt edgeStart, edgeEnd, pStart, pEnd;
PetscInt cellStart, cellEnd;
PetscInt faceStart, faceEnd;

PetscFunctionBeginUser;
DMPlexGetChart(dm, &pStart, &pEnd);
DMPlexGetHeightStratum(dm, 0, &cellStart, &cellEnd);
DMPlexGetHeightStratum(dm, 1, &faceStart, &faceEnd);
DMPlexGetHeightStratum(dm, 2, &edgeStart, &edgeEnd); /* edges */;
PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s);
PetscSectionSetNumFields(s, 1);
PetscSectionSetFieldComponents(s, 0, 1);
if (objCaseInfo.getnumberDof_local() == 6){
PetscSectionSetChart(s, edgeStart, edgeEnd);
for (PetscInt edgeIndex = edgeStart; edgeIndex < edgeEnd; ++edgeIndex) {
PetscSectionSetDof(s, edgeIndex, objCaseInfo.numdofPerEdge);
PetscSectionSetFieldDof(s, edgeIndex, 0, 1);
}
}
else if(objCaseInfo.getnumberDof_local() == 20){
PetscSectionSetChart(s, faceStart, edgeEnd);
for (PetscInt faceIndex = faceStart; faceIndex < faceEnd; ++faceIndex) {
PetscSectionSetDof(s, faceIndex, objCaseInfo.numdofPerFace);
PetscSectionSetFieldDof(s, faceIndex, 0, 1);
}
//Test
for (PetscInt edgeIndex = edgeStart; edgeIndex < edgeEnd; ++edgeIndex) {
PetscSectionSetDof(s, edgeIndex, objCaseInfo.numdofPerEdge);
PetscSectionSetFieldDof(s, edgeIndex, 0, 1);
}
}
//
PetscSectionSetUp(s);
DMSetLocalSection(dm, s);
PetscSectionDestroy(&s);

//Output map for check
ISLocalToGlobalMapping ltogm;
const PetscInt *g_idx;
DMGetLocalToGlobalMapping(dm, &ltogm);
ISLocalToGlobalMappingView(ltogm, PETSC_VIEWER_STDOUT_WORLD);
ISLocalToGlobalMappingGetIndices(ltogm, &g_idx);

PetscFunctionReturn(PETSC_SUCCESS);
}