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 Level: intermediate
13
14 Note:
15 This should greatly improve the performance of the closure operations, at the cost of additional memory.
16
17 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSection`, `DMPlexVecGetClosure()`, `DMPlexVecRestoreClosure()`, `DMPlexVecSetClosure()`, `DMPlexMatSetClosure()`
18 @*/
DMPlexCreateClosureIndex(DM dm,PetscSection section)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(PETSC_SUCCESS);
74 }
75