1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 /*@ 4 DMPlexCreateClosureIndex - Calculate an index for the given PetscSection for the closure operation on the DM 5 6 Not collective 7 8 Input Parameters: 9 + dm - The DM 10 - section - The section describing the layout in the local vector, or NULL to use the default section 11 12 Note: 13 This should greatly improve the performance of the closure operations, at the cost of additional memory. 14 15 Level: intermediate 16 17 .seealso `DMPlexVecGetClosure()`, `DMPlexVecRestoreClosure()`, `DMPlexVecSetClosure()`, `DMPlexMatSetClosure()` 18 @*/ 19 PetscErrorCode DMPlexCreateClosureIndex(DM dm, PetscSection section) 20 { 21 PetscSection closureSection; 22 IS closureIS; 23 PetscInt *clPoints; 24 PetscInt pStart, pEnd, sStart, sEnd, point, clSize; 25 26 PetscFunctionBegin; 27 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 28 if (!section) PetscCall(DMGetLocalSection(dm, §ion)); 29 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 2); 30 PetscCall(PetscSectionGetChart(section, &sStart, &sEnd)); 31 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 32 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject) section), &closureSection)); 33 PetscCall(PetscSectionSetChart(closureSection, pStart, pEnd)); 34 for (point = pStart; point < pEnd; ++point) { 35 PetscInt *points = NULL, numPoints, p, dof, cldof = 0; 36 37 PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points)); 38 for (p = 0; p < numPoints*2; p += 2) { 39 if ((points[p] >= sStart) && (points[p] < sEnd)) { 40 PetscCall(PetscSectionGetDof(section, points[p], &dof)); 41 if (dof) cldof += 2; 42 } 43 } 44 PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points)); 45 PetscCall(PetscSectionSetDof(closureSection, point, cldof)); 46 } 47 PetscCall(PetscSectionSetUp(closureSection)); 48 PetscCall(PetscSectionGetStorageSize(closureSection, &clSize)); 49 PetscCall(PetscMalloc1(clSize, &clPoints)); 50 for (point = pStart; point < pEnd; ++point) { 51 PetscInt *points = NULL, numPoints, p, q, dof, cldof, cloff; 52 53 PetscCall(PetscSectionGetDof(closureSection, point, &cldof)); 54 PetscCall(PetscSectionGetOffset(closureSection, point, &cloff)); 55 PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points)); 56 for (p = 0, q = 0; p < numPoints*2; p += 2) { 57 if ((points[p] >= sStart) && (points[p] < sEnd)) { 58 PetscCall(PetscSectionGetDof(section, points[p], &dof)); 59 if (dof) { 60 clPoints[cloff+q*2] = points[p]; 61 clPoints[cloff+q*2+1] = points[p+1]; 62 ++q; 63 } 64 } 65 } 66 PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points)); 67 PetscCheck(q*2 == cldof,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid size for closure %" PetscInt_FMT " should be %" PetscInt_FMT, q*2, cldof); 68 } 69 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPoints, PETSC_OWN_POINTER, &closureIS)); 70 PetscCall(PetscSectionSetClosureIndex(section, (PetscObject) dm, closureSection, closureIS)); 71 PetscCall(PetscSectionDestroy(&closureSection)); 72 PetscCall(ISDestroy(&closureIS)); 73 PetscFunctionReturn(0); 74 } 75