Actual source code: plexindices.c
1: #include <petsc/private/dmpleximpl.h>
3: /*@
4: DMPlexCreateClosureIndex - Calculate an index for the given `PetscSection` for the closure operation on the `DM`
6: Not Collective
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
12: Level: intermediate
14: Note:
15: This should greatly improve the performance of the closure operations, at the cost of additional memory.
17: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `PetscSection`, `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;
27: if (!section) DMGetLocalSection(dm, §ion);
29: PetscSectionGetChart(section, &sStart, &sEnd);
30: DMPlexGetChart(dm, &pStart, &pEnd);
31: PetscSectionCreate(PetscObjectComm((PetscObject)section), &closureSection);
32: PetscSectionSetChart(closureSection, pStart, pEnd);
33: for (point = pStart; point < pEnd; ++point) {
34: PetscInt *points = NULL, numPoints, p, dof, cldof = 0;
36: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
37: for (p = 0; p < numPoints * 2; p += 2) {
38: if ((points[p] >= sStart) && (points[p] < sEnd)) {
39: PetscSectionGetDof(section, points[p], &dof);
40: if (dof) cldof += 2;
41: }
42: }
43: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
44: PetscSectionSetDof(closureSection, point, cldof);
45: }
46: PetscSectionSetUp(closureSection);
47: PetscSectionGetStorageSize(closureSection, &clSize);
48: PetscMalloc1(clSize, &clPoints);
49: for (point = pStart; point < pEnd; ++point) {
50: PetscInt *points = NULL, numPoints, p, q, dof, cldof, cloff;
52: PetscSectionGetDof(closureSection, point, &cldof);
53: PetscSectionGetOffset(closureSection, point, &cloff);
54: DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
55: for (p = 0, q = 0; p < numPoints * 2; p += 2) {
56: if ((points[p] >= sStart) && (points[p] < sEnd)) {
57: PetscSectionGetDof(section, points[p], &dof);
58: if (dof) {
59: clPoints[cloff + q * 2] = points[p];
60: clPoints[cloff + q * 2 + 1] = points[p + 1];
61: ++q;
62: }
63: }
64: }
65: DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
67: }
68: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPoints, PETSC_OWN_POINTER, &closureIS);
69: PetscSectionSetClosureIndex(section, (PetscObject)dm, closureSection, closureIS);
70: PetscSectionDestroy(&closureSection);
71: ISDestroy(&closureIS);
72: return 0;
73: }