Actual source code: plextree.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/isimpl.h>
3: #include <petsc/private/petscfeimpl.h>
4: #include <petscsf.h>
5: #include <petscds.h>
7: /* hierarchy routines */
9: /*@
10: DMPlexSetReferenceTree - set the reference tree for hierarchically non-conforming meshes.
12: Not Collective
14: Input Parameters:
15: + dm - The `DMPLEX` object
16: - ref - The reference tree `DMPLEX` object
18: Level: intermediate
20: .seealso: [](ch_unstructured), `DM`, `DMPLEX`,`DMPlexGetReferenceTree()`, `DMPlexCreateDefaultReferenceTree()`
21: @*/
22: PetscErrorCode DMPlexSetReferenceTree(DM dm, DM ref)
23: {
24: DM_Plex *mesh = (DM_Plex *)dm->data;
26: PetscFunctionBegin;
29: PetscCall(PetscObjectReference((PetscObject)ref));
30: PetscCall(DMDestroy(&mesh->referenceTree));
31: mesh->referenceTree = ref;
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: /*@
36: DMPlexGetReferenceTree - get the reference tree for hierarchically non-conforming meshes.
38: Not Collective
40: Input Parameter:
41: . dm - The `DMPLEX` object
43: Output Parameter:
44: . ref - The reference tree `DMPLEX` object
46: Level: intermediate
48: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetReferenceTree()`, `DMPlexCreateDefaultReferenceTree()`
49: @*/
50: PetscErrorCode DMPlexGetReferenceTree(DM dm, DM *ref)
51: {
52: DM_Plex *mesh = (DM_Plex *)dm->data;
54: PetscFunctionBegin;
56: PetscAssertPointer(ref, 2);
57: *ref = mesh->referenceTree;
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: static PetscErrorCode DMPlexReferenceTreeGetChildSymmetry_Default(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
62: {
63: PetscInt coneSize, dStart, dEnd, dim, ABswap, oAvert, oBvert, ABswapVert;
65: PetscFunctionBegin;
66: if (parentOrientA == parentOrientB) {
67: if (childOrientB) *childOrientB = childOrientA;
68: if (childB) *childB = childA;
69: PetscFunctionReturn(PETSC_SUCCESS);
70: }
71: for (dim = 0; dim < 3; dim++) {
72: PetscCall(DMPlexGetDepthStratum(dm, dim, &dStart, &dEnd));
73: if (parent >= dStart && parent <= dEnd) break;
74: }
75: PetscCheck(dim <= 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot perform child symmetry for %" PetscInt_FMT "-cells", dim);
76: PetscCheck(dim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "A vertex has no children");
77: if (childA < dStart || childA >= dEnd) {
78: /* this is a lower-dimensional child: bootstrap */
79: PetscInt size, i, sA = -1, sB, sOrientB, sConeSize;
80: const PetscInt *supp, *coneA, *coneB, *oA, *oB;
82: PetscCall(DMPlexGetSupportSize(dm, childA, &size));
83: PetscCall(DMPlexGetSupport(dm, childA, &supp));
85: /* find a point sA in supp(childA) that has the same parent */
86: for (i = 0; i < size; i++) {
87: PetscInt sParent;
89: sA = supp[i];
90: if (sA == parent) continue;
91: PetscCall(DMPlexGetTreeParent(dm, sA, &sParent, NULL));
92: if (sParent == parent) break;
93: }
94: PetscCheck(i != size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "could not find support in children");
95: /* find out which point sB is in an equivalent position to sA under
96: * parentOrientB */
97: PetscCall(DMPlexReferenceTreeGetChildSymmetry_Default(dm, parent, parentOrientA, 0, sA, parentOrientB, &sOrientB, &sB));
98: PetscCall(DMPlexGetConeSize(dm, sA, &sConeSize));
99: PetscCall(DMPlexGetCone(dm, sA, &coneA));
100: PetscCall(DMPlexGetCone(dm, sB, &coneB));
101: PetscCall(DMPlexGetConeOrientation(dm, sA, &oA));
102: PetscCall(DMPlexGetConeOrientation(dm, sB, &oB));
103: /* step through the cone of sA in natural order */
104: for (i = 0; i < sConeSize; i++) {
105: if (coneA[i] == childA) {
106: /* if childA is at position i in coneA,
107: * then we want the point that is at sOrientB*i in coneB */
108: PetscInt j = (sOrientB >= 0) ? ((sOrientB + i) % sConeSize) : ((sConeSize - (sOrientB + 1) - i) % sConeSize);
109: if (childB) *childB = coneB[j];
110: if (childOrientB) {
111: DMPolytopeType ct;
112: PetscInt oBtrue;
114: PetscCall(DMPlexGetConeSize(dm, childA, &coneSize));
115: /* compose sOrientB and oB[j] */
116: PetscCheck(coneSize == 0 || coneSize == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected a vertex or an edge");
117: ct = coneSize ? DM_POLYTOPE_SEGMENT : DM_POLYTOPE_POINT;
118: /* we may have to flip an edge */
119: oBtrue = (sOrientB >= 0) ? oB[j] : DMPolytopeTypeComposeOrientation(ct, -1, oB[j]);
120: oBtrue = DMPolytopeConvertNewOrientation_Internal(ct, oBtrue);
121: ABswap = DihedralSwap(coneSize, DMPolytopeConvertNewOrientation_Internal(ct, oA[i]), oBtrue);
122: *childOrientB = DihedralCompose(coneSize, childOrientA, ABswap);
123: }
124: break;
125: }
126: }
127: PetscCheck(i != sConeSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "support cone mismatch");
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
130: /* get the cone size and symmetry swap */
131: PetscCall(DMPlexGetConeSize(dm, parent, &coneSize));
132: ABswap = DihedralSwap(coneSize, parentOrientA, parentOrientB);
133: if (dim == 2) {
134: /* orientations refer to cones: we want them to refer to vertices:
135: * if it's a rotation, they are the same, but if the order is reversed, a
136: * permutation that puts side i first does *not* put vertex i first */
137: oAvert = (parentOrientA >= 0) ? parentOrientA : -((-parentOrientA % coneSize) + 1);
138: oBvert = (parentOrientB >= 0) ? parentOrientB : -((-parentOrientB % coneSize) + 1);
139: ABswapVert = DihedralSwap(coneSize, oAvert, oBvert);
140: } else {
141: ABswapVert = ABswap;
142: }
143: if (childB) {
144: /* assume that each child corresponds to a vertex, in the same order */
145: PetscInt p, posA = -1, numChildren, i;
146: const PetscInt *children;
148: /* count which position the child is in */
149: PetscCall(DMPlexGetTreeChildren(dm, parent, &numChildren, &children));
150: for (i = 0; i < numChildren; i++) {
151: p = children[i];
152: if (p == childA) {
153: posA = i;
154: break;
155: }
156: }
157: if (posA >= coneSize) {
158: /* this is the triangle in the middle of a uniformly refined triangle: it
159: * is invariant */
160: PetscCheck(dim == 2 && posA == 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Expected a middle triangle, got something else");
161: *childB = childA;
162: } else {
163: /* figure out position B by applying ABswapVert */
164: PetscInt posB;
166: posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize - (ABswapVert + 1) - posA) % coneSize);
167: if (childB) *childB = children[posB];
168: }
169: }
170: if (childOrientB) *childOrientB = DihedralCompose(coneSize, childOrientA, ABswap);
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: /*@
175: DMPlexReferenceTreeGetChildSymmetry - Given a reference tree, transform a childid and orientation from one parent frame to another
177: Input Parameters:
178: + dm - the reference tree `DMPLEX` object
179: . parent - the parent point
180: . parentOrientA - the reference orientation for describing the parent
181: . childOrientA - the reference orientation for describing the child
182: . childA - the reference childID for describing the child
183: - parentOrientB - the new orientation for describing the parent
185: Output Parameters:
186: + childOrientB - if not `NULL`, set to the new orientation for describing the child
187: - childB - if not `NULL`, the new childID for describing the child
189: Level: developer
191: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetReferenceTree()`, `DMPlexSetReferenceTree()`, `DMPlexSetTree()`
192: @*/
193: PetscErrorCode DMPlexReferenceTreeGetChildSymmetry(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
194: {
195: DM_Plex *mesh = (DM_Plex *)dm->data;
197: PetscFunctionBegin;
199: PetscCheck(mesh->getchildsymmetry, PETSC_COMM_SELF, PETSC_ERR_SUP, "DMPlexReferenceTreeGetChildSymmetry not implemented");
200: PetscCall(mesh->getchildsymmetry(dm, parent, parentOrientA, childOrientA, childA, parentOrientB, childOrientB, childB));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: static PetscErrorCode DMPlexSetTree_Internal(DM, PetscSection, PetscInt *, PetscInt *, PetscBool, PetscBool);
206: PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
207: {
208: PetscFunctionBegin;
209: PetscCall(DMPlexSetTree_Internal(dm, parentSection, parents, childIDs, PETSC_TRUE, PETSC_FALSE));
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: PetscErrorCode DMPlexCreateReferenceTree_Union(DM K, DM Kref, const char *labelName, DM *ref)
214: {
215: MPI_Comm comm;
216: PetscInt dim, p, pStart, pEnd, pRefStart, pRefEnd, d, offset, parentSize, *parents, *childIDs;
217: PetscInt *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts;
218: DMLabel identity, identityRef;
219: PetscSection unionSection, unionConeSection, parentSection;
220: PetscScalar *unionCoords;
221: IS perm;
223: PetscFunctionBegin;
224: comm = PetscObjectComm((PetscObject)K);
225: PetscCall(DMGetDimension(K, &dim));
226: PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
227: PetscCall(DMGetLabel(K, labelName, &identity));
228: PetscCall(DMGetLabel(Kref, labelName, &identityRef));
229: PetscCall(DMPlexGetChart(Kref, &pRefStart, &pRefEnd));
230: PetscCall(PetscSectionCreate(comm, &unionSection));
231: PetscCall(PetscSectionSetChart(unionSection, 0, (pEnd - pStart) + (pRefEnd - pRefStart)));
232: /* count points that will go in the union */
233: for (p = pStart; p < pEnd; p++) PetscCall(PetscSectionSetDof(unionSection, p - pStart, 1));
234: for (p = pRefStart; p < pRefEnd; p++) {
235: PetscInt q, qSize;
236: PetscCall(DMLabelGetValue(identityRef, p, &q));
237: PetscCall(DMLabelGetStratumSize(identityRef, q, &qSize));
238: if (qSize > 1) PetscCall(PetscSectionSetDof(unionSection, p - pRefStart + (pEnd - pStart), 1));
239: }
240: PetscCall(PetscMalloc1(pEnd - pStart + pRefEnd - pRefStart, &permvals));
241: offset = 0;
242: /* stratify points in the union by topological dimension */
243: for (d = 0; d <= dim; d++) {
244: PetscInt cStart, cEnd, c;
246: PetscCall(DMPlexGetHeightStratum(K, d, &cStart, &cEnd));
247: for (c = cStart; c < cEnd; c++) permvals[offset++] = c;
249: PetscCall(DMPlexGetHeightStratum(Kref, d, &cStart, &cEnd));
250: for (c = cStart; c < cEnd; c++) permvals[offset++] = c + (pEnd - pStart);
251: }
252: PetscCall(ISCreateGeneral(comm, (pEnd - pStart) + (pRefEnd - pRefStart), permvals, PETSC_OWN_POINTER, &perm));
253: PetscCall(PetscSectionSetPermutation(unionSection, perm));
254: PetscCall(PetscSectionSetUp(unionSection));
255: PetscCall(PetscSectionGetStorageSize(unionSection, &numUnionPoints));
256: PetscCall(PetscMalloc2(numUnionPoints, &coneSizes, dim + 1, &numDimPoints));
257: /* count dimension points */
258: for (d = 0; d <= dim; d++) {
259: PetscInt cStart, cOff, cOff2;
260: PetscCall(DMPlexGetHeightStratum(K, d, &cStart, NULL));
261: PetscCall(PetscSectionGetOffset(unionSection, cStart - pStart, &cOff));
262: if (d < dim) {
263: PetscCall(DMPlexGetHeightStratum(K, d + 1, &cStart, NULL));
264: PetscCall(PetscSectionGetOffset(unionSection, cStart - pStart, &cOff2));
265: } else {
266: cOff2 = numUnionPoints;
267: }
268: numDimPoints[dim - d] = cOff2 - cOff;
269: }
270: PetscCall(PetscSectionCreate(comm, &unionConeSection));
271: PetscCall(PetscSectionSetChart(unionConeSection, 0, numUnionPoints));
272: /* count the cones in the union */
273: for (p = pStart; p < pEnd; p++) {
274: PetscInt dof, uOff;
276: PetscCall(DMPlexGetConeSize(K, p, &dof));
277: PetscCall(PetscSectionGetOffset(unionSection, p - pStart, &uOff));
278: PetscCall(PetscSectionSetDof(unionConeSection, uOff, dof));
279: coneSizes[uOff] = dof;
280: }
281: for (p = pRefStart; p < pRefEnd; p++) {
282: PetscInt dof, uDof, uOff;
284: PetscCall(DMPlexGetConeSize(Kref, p, &dof));
285: PetscCall(PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart), &uDof));
286: PetscCall(PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart), &uOff));
287: if (uDof) {
288: PetscCall(PetscSectionSetDof(unionConeSection, uOff, dof));
289: coneSizes[uOff] = dof;
290: }
291: }
292: PetscCall(PetscSectionSetUp(unionConeSection));
293: PetscCall(PetscSectionGetStorageSize(unionConeSection, &numCones));
294: PetscCall(PetscMalloc2(numCones, &unionCones, numCones, &unionOrientations));
295: /* write the cones in the union */
296: for (p = pStart; p < pEnd; p++) {
297: PetscInt dof, uOff, c, cOff;
298: const PetscInt *cone, *orientation;
300: PetscCall(DMPlexGetConeSize(K, p, &dof));
301: PetscCall(DMPlexGetCone(K, p, &cone));
302: PetscCall(DMPlexGetConeOrientation(K, p, &orientation));
303: PetscCall(PetscSectionGetOffset(unionSection, p - pStart, &uOff));
304: PetscCall(PetscSectionGetOffset(unionConeSection, uOff, &cOff));
305: for (c = 0; c < dof; c++) {
306: PetscInt e, eOff;
307: e = cone[c];
308: PetscCall(PetscSectionGetOffset(unionSection, e - pStart, &eOff));
309: unionCones[cOff + c] = eOff;
310: unionOrientations[cOff + c] = orientation[c];
311: }
312: }
313: for (p = pRefStart; p < pRefEnd; p++) {
314: PetscInt dof, uDof, uOff, c, cOff;
315: const PetscInt *cone, *orientation;
317: PetscCall(DMPlexGetConeSize(Kref, p, &dof));
318: PetscCall(DMPlexGetCone(Kref, p, &cone));
319: PetscCall(DMPlexGetConeOrientation(Kref, p, &orientation));
320: PetscCall(PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart), &uDof));
321: PetscCall(PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart), &uOff));
322: if (uDof) {
323: PetscCall(PetscSectionGetOffset(unionConeSection, uOff, &cOff));
324: for (c = 0; c < dof; c++) {
325: PetscInt e, eOff, eDof;
327: e = cone[c];
328: PetscCall(PetscSectionGetDof(unionSection, e - pRefStart + (pEnd - pStart), &eDof));
329: if (eDof) {
330: PetscCall(PetscSectionGetOffset(unionSection, e - pRefStart + (pEnd - pStart), &eOff));
331: } else {
332: PetscCall(DMLabelGetValue(identityRef, e, &e));
333: PetscCall(PetscSectionGetOffset(unionSection, e - pStart, &eOff));
334: }
335: unionCones[cOff + c] = eOff;
336: unionOrientations[cOff + c] = orientation[c];
337: }
338: }
339: }
340: /* get the coordinates */
341: {
342: PetscInt vStart, vEnd, vRefStart, vRefEnd, v, vDof, vOff;
343: PetscSection KcoordsSec, KrefCoordsSec;
344: Vec KcoordsVec, KrefCoordsVec;
345: PetscScalar *Kcoords;
347: PetscCall(DMGetCoordinateSection(K, &KcoordsSec));
348: PetscCall(DMGetCoordinatesLocal(K, &KcoordsVec));
349: PetscCall(DMGetCoordinateSection(Kref, &KrefCoordsSec));
350: PetscCall(DMGetCoordinatesLocal(Kref, &KrefCoordsVec));
352: numVerts = numDimPoints[0];
353: PetscCall(PetscMalloc1(numVerts * dim, &unionCoords));
354: PetscCall(DMPlexGetDepthStratum(K, 0, &vStart, &vEnd));
356: offset = 0;
357: for (v = vStart; v < vEnd; v++) {
358: PetscCall(PetscSectionGetOffset(unionSection, v - pStart, &vOff));
359: PetscCall(VecGetValuesSection(KcoordsVec, KcoordsSec, v, &Kcoords));
360: for (d = 0; d < dim; d++) unionCoords[offset * dim + d] = Kcoords[d];
361: offset++;
362: }
363: PetscCall(DMPlexGetDepthStratum(Kref, 0, &vRefStart, &vRefEnd));
364: for (v = vRefStart; v < vRefEnd; v++) {
365: PetscCall(PetscSectionGetDof(unionSection, v - pRefStart + (pEnd - pStart), &vDof));
366: PetscCall(PetscSectionGetOffset(unionSection, v - pRefStart + (pEnd - pStart), &vOff));
367: PetscCall(VecGetValuesSection(KrefCoordsVec, KrefCoordsSec, v, &Kcoords));
368: if (vDof) {
369: for (d = 0; d < dim; d++) unionCoords[offset * dim + d] = Kcoords[d];
370: offset++;
371: }
372: }
373: }
374: PetscCall(DMCreate(comm, ref));
375: PetscCall(DMSetType(*ref, DMPLEX));
376: PetscCall(DMSetDimension(*ref, dim));
377: PetscCall(DMPlexCreateFromDAG(*ref, dim, numDimPoints, coneSizes, unionCones, unionOrientations, unionCoords));
378: /* set the tree */
379: PetscCall(PetscSectionCreate(comm, &parentSection));
380: PetscCall(PetscSectionSetChart(parentSection, 0, numUnionPoints));
381: for (p = pRefStart; p < pRefEnd; p++) {
382: PetscInt uDof, uOff;
384: PetscCall(PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart), &uDof));
385: PetscCall(PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart), &uOff));
386: if (uDof) PetscCall(PetscSectionSetDof(parentSection, uOff, 1));
387: }
388: PetscCall(PetscSectionSetUp(parentSection));
389: PetscCall(PetscSectionGetStorageSize(parentSection, &parentSize));
390: PetscCall(PetscMalloc2(parentSize, &parents, parentSize, &childIDs));
391: for (p = pRefStart; p < pRefEnd; p++) {
392: PetscInt uDof, uOff;
394: PetscCall(PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart), &uDof));
395: PetscCall(PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart), &uOff));
396: if (uDof) {
397: PetscInt pOff, parent, parentU;
398: PetscCall(PetscSectionGetOffset(parentSection, uOff, &pOff));
399: PetscCall(DMLabelGetValue(identityRef, p, &parent));
400: PetscCall(PetscSectionGetOffset(unionSection, parent - pStart, &parentU));
401: parents[pOff] = parentU;
402: childIDs[pOff] = uOff;
403: }
404: }
405: PetscCall(DMPlexCreateReferenceTree_SetTree(*ref, parentSection, parents, childIDs));
406: PetscCall(PetscSectionDestroy(&parentSection));
407: PetscCall(PetscFree2(parents, childIDs));
409: /* clean up */
410: PetscCall(PetscSectionDestroy(&unionSection));
411: PetscCall(PetscSectionDestroy(&unionConeSection));
412: PetscCall(ISDestroy(&perm));
413: PetscCall(PetscFree(unionCoords));
414: PetscCall(PetscFree2(unionCones, unionOrientations));
415: PetscCall(PetscFree2(coneSizes, numDimPoints));
416: PetscFunctionReturn(PETSC_SUCCESS);
417: }
419: /*@
420: DMPlexCreateDefaultReferenceTree - create a reference tree for isotropic hierarchical mesh refinement.
422: Collective
424: Input Parameters:
425: + comm - the MPI communicator
426: . dim - the spatial dimension
427: - simplex - Flag for simplex, otherwise use a tensor-product cell
429: Output Parameter:
430: . ref - the reference tree `DMPLEX` object
432: Level: intermediate
434: .seealso: `DMPlexSetReferenceTree()`, `DMPlexGetReferenceTree()`
435: @*/
436: PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref)
437: {
438: DM_Plex *mesh;
439: DM K, Kref;
440: PetscInt p, pStart, pEnd;
441: DMLabel identity;
443: PetscFunctionBegin;
444: #if 1
445: comm = PETSC_COMM_SELF;
446: #endif
447: /* create a reference element */
448: PetscCall(DMPlexCreateReferenceCell(comm, DMPolytopeTypeSimpleShape(dim, simplex), &K));
449: PetscCall(DMCreateLabel(K, "identity"));
450: PetscCall(DMGetLabel(K, "identity", &identity));
451: PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
452: for (p = pStart; p < pEnd; p++) PetscCall(DMLabelSetValue(identity, p, p));
453: /* refine it */
454: PetscCall(DMRefine(K, comm, &Kref));
456: /* the reference tree is the union of these two, without duplicating
457: * points that appear in both */
458: PetscCall(DMPlexCreateReferenceTree_Union(K, Kref, "identity", ref));
459: mesh = (DM_Plex *)(*ref)->data;
460: mesh->getchildsymmetry = DMPlexReferenceTreeGetChildSymmetry_Default;
461: PetscCall(DMDestroy(&K));
462: PetscCall(DMDestroy(&Kref));
463: PetscFunctionReturn(PETSC_SUCCESS);
464: }
466: static PetscErrorCode DMPlexTreeSymmetrize(DM dm)
467: {
468: DM_Plex *mesh = (DM_Plex *)dm->data;
469: PetscSection childSec, pSec;
470: PetscInt p, pSize, cSize, parMax = PETSC_MIN_INT, parMin = PETSC_MAX_INT;
471: PetscInt *offsets, *children, pStart, pEnd;
473: PetscFunctionBegin;
475: PetscCall(PetscSectionDestroy(&mesh->childSection));
476: PetscCall(PetscFree(mesh->children));
477: pSec = mesh->parentSection;
478: if (!pSec) PetscFunctionReturn(PETSC_SUCCESS);
479: PetscCall(PetscSectionGetStorageSize(pSec, &pSize));
480: for (p = 0; p < pSize; p++) {
481: PetscInt par = mesh->parents[p];
483: parMax = PetscMax(parMax, par + 1);
484: parMin = PetscMin(parMin, par);
485: }
486: if (parMin > parMax) {
487: parMin = -1;
488: parMax = -1;
489: }
490: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)pSec), &childSec));
491: PetscCall(PetscSectionSetChart(childSec, parMin, parMax));
492: for (p = 0; p < pSize; p++) {
493: PetscInt par = mesh->parents[p];
495: PetscCall(PetscSectionAddDof(childSec, par, 1));
496: }
497: PetscCall(PetscSectionSetUp(childSec));
498: PetscCall(PetscSectionGetStorageSize(childSec, &cSize));
499: PetscCall(PetscMalloc1(cSize, &children));
500: PetscCall(PetscCalloc1(parMax - parMin, &offsets));
501: PetscCall(PetscSectionGetChart(pSec, &pStart, &pEnd));
502: for (p = pStart; p < pEnd; p++) {
503: PetscInt dof, off, i;
505: PetscCall(PetscSectionGetDof(pSec, p, &dof));
506: PetscCall(PetscSectionGetOffset(pSec, p, &off));
507: for (i = 0; i < dof; i++) {
508: PetscInt par = mesh->parents[off + i], cOff;
510: PetscCall(PetscSectionGetOffset(childSec, par, &cOff));
511: children[cOff + offsets[par - parMin]++] = p;
512: }
513: }
514: mesh->childSection = childSec;
515: mesh->children = children;
516: PetscCall(PetscFree(offsets));
517: PetscFunctionReturn(PETSC_SUCCESS);
518: }
520: static PetscErrorCode AnchorsFlatten(PetscSection section, IS is, PetscSection *sectionNew, IS *isNew)
521: {
522: PetscInt pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL;
523: const PetscInt *vals;
524: PetscSection secNew;
525: PetscBool anyNew, globalAnyNew;
526: PetscBool compress;
528: PetscFunctionBegin;
529: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
530: PetscCall(ISGetLocalSize(is, &size));
531: PetscCall(ISGetIndices(is, &vals));
532: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), &secNew));
533: PetscCall(PetscSectionSetChart(secNew, pStart, pEnd));
534: for (i = 0; i < size; i++) {
535: PetscInt dof;
537: p = vals[i];
538: if (p < pStart || p >= pEnd) continue;
539: PetscCall(PetscSectionGetDof(section, p, &dof));
540: if (dof) break;
541: }
542: if (i == size) {
543: PetscCall(PetscSectionSetUp(secNew));
544: anyNew = PETSC_FALSE;
545: compress = PETSC_FALSE;
546: sizeNew = 0;
547: } else {
548: anyNew = PETSC_TRUE;
549: for (p = pStart; p < pEnd; p++) {
550: PetscInt dof, off;
552: PetscCall(PetscSectionGetDof(section, p, &dof));
553: PetscCall(PetscSectionGetOffset(section, p, &off));
554: for (i = 0; i < dof; i++) {
555: PetscInt q = vals[off + i], qDof = 0;
557: if (q >= pStart && q < pEnd) PetscCall(PetscSectionGetDof(section, q, &qDof));
558: if (qDof) PetscCall(PetscSectionAddDof(secNew, p, qDof));
559: else PetscCall(PetscSectionAddDof(secNew, p, 1));
560: }
561: }
562: PetscCall(PetscSectionSetUp(secNew));
563: PetscCall(PetscSectionGetStorageSize(secNew, &sizeNew));
564: PetscCall(PetscMalloc1(sizeNew, &valsNew));
565: compress = PETSC_FALSE;
566: for (p = pStart; p < pEnd; p++) {
567: PetscInt dof, off, count, offNew, dofNew;
569: PetscCall(PetscSectionGetDof(section, p, &dof));
570: PetscCall(PetscSectionGetOffset(section, p, &off));
571: PetscCall(PetscSectionGetDof(secNew, p, &dofNew));
572: PetscCall(PetscSectionGetOffset(secNew, p, &offNew));
573: count = 0;
574: for (i = 0; i < dof; i++) {
575: PetscInt q = vals[off + i], qDof = 0, qOff = 0, j;
577: if (q >= pStart && q < pEnd) {
578: PetscCall(PetscSectionGetDof(section, q, &qDof));
579: PetscCall(PetscSectionGetOffset(section, q, &qOff));
580: }
581: if (qDof) {
582: PetscInt oldCount = count;
584: for (j = 0; j < qDof; j++) {
585: PetscInt k, r = vals[qOff + j];
587: for (k = 0; k < oldCount; k++) {
588: if (valsNew[offNew + k] == r) break;
589: }
590: if (k == oldCount) valsNew[offNew + count++] = r;
591: }
592: } else {
593: PetscInt k, oldCount = count;
595: for (k = 0; k < oldCount; k++) {
596: if (valsNew[offNew + k] == q) break;
597: }
598: if (k == oldCount) valsNew[offNew + count++] = q;
599: }
600: }
601: if (count < dofNew) {
602: PetscCall(PetscSectionSetDof(secNew, p, count));
603: compress = PETSC_TRUE;
604: }
605: }
606: }
607: PetscCall(ISRestoreIndices(is, &vals));
608: PetscCall(MPIU_Allreduce(&anyNew, &globalAnyNew, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)secNew)));
609: if (!globalAnyNew) {
610: PetscCall(PetscSectionDestroy(&secNew));
611: *sectionNew = NULL;
612: *isNew = NULL;
613: } else {
614: PetscBool globalCompress;
616: PetscCall(MPIU_Allreduce(&compress, &globalCompress, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)secNew)));
617: if (compress) {
618: PetscSection secComp;
619: PetscInt *valsComp = NULL;
621: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), &secComp));
622: PetscCall(PetscSectionSetChart(secComp, pStart, pEnd));
623: for (p = pStart; p < pEnd; p++) {
624: PetscInt dof;
626: PetscCall(PetscSectionGetDof(secNew, p, &dof));
627: PetscCall(PetscSectionSetDof(secComp, p, dof));
628: }
629: PetscCall(PetscSectionSetUp(secComp));
630: PetscCall(PetscSectionGetStorageSize(secComp, &sizeNew));
631: PetscCall(PetscMalloc1(sizeNew, &valsComp));
632: for (p = pStart; p < pEnd; p++) {
633: PetscInt dof, off, offNew, j;
635: PetscCall(PetscSectionGetDof(secNew, p, &dof));
636: PetscCall(PetscSectionGetOffset(secNew, p, &off));
637: PetscCall(PetscSectionGetOffset(secComp, p, &offNew));
638: for (j = 0; j < dof; j++) valsComp[offNew + j] = valsNew[off + j];
639: }
640: PetscCall(PetscSectionDestroy(&secNew));
641: secNew = secComp;
642: PetscCall(PetscFree(valsNew));
643: valsNew = valsComp;
644: }
645: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), sizeNew, valsNew, PETSC_OWN_POINTER, isNew));
646: }
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
650: static PetscErrorCode DMPlexCreateAnchors_Tree(DM dm)
651: {
652: PetscInt p, pStart, pEnd, *anchors, size;
653: PetscInt aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT;
654: PetscSection aSec;
655: DMLabel canonLabel;
656: IS aIS;
658: PetscFunctionBegin;
660: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
661: PetscCall(DMGetLabel(dm, "canonical", &canonLabel));
662: for (p = pStart; p < pEnd; p++) {
663: PetscInt parent;
665: if (canonLabel) {
666: PetscInt canon;
668: PetscCall(DMLabelGetValue(canonLabel, p, &canon));
669: if (p != canon) continue;
670: }
671: PetscCall(DMPlexGetTreeParent(dm, p, &parent, NULL));
672: if (parent != p) {
673: aMin = PetscMin(aMin, p);
674: aMax = PetscMax(aMax, p + 1);
675: }
676: }
677: if (aMin > aMax) {
678: aMin = -1;
679: aMax = -1;
680: }
681: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &aSec));
682: PetscCall(PetscSectionSetChart(aSec, aMin, aMax));
683: for (p = aMin; p < aMax; p++) {
684: PetscInt parent, ancestor = p;
686: if (canonLabel) {
687: PetscInt canon;
689: PetscCall(DMLabelGetValue(canonLabel, p, &canon));
690: if (p != canon) continue;
691: }
692: PetscCall(DMPlexGetTreeParent(dm, p, &parent, NULL));
693: while (parent != ancestor) {
694: ancestor = parent;
695: PetscCall(DMPlexGetTreeParent(dm, ancestor, &parent, NULL));
696: }
697: if (ancestor != p) {
698: PetscInt closureSize, *closure = NULL;
700: PetscCall(DMPlexGetTransitiveClosure(dm, ancestor, PETSC_TRUE, &closureSize, &closure));
701: PetscCall(PetscSectionSetDof(aSec, p, closureSize));
702: PetscCall(DMPlexRestoreTransitiveClosure(dm, ancestor, PETSC_TRUE, &closureSize, &closure));
703: }
704: }
705: PetscCall(PetscSectionSetUp(aSec));
706: PetscCall(PetscSectionGetStorageSize(aSec, &size));
707: PetscCall(PetscMalloc1(size, &anchors));
708: for (p = aMin; p < aMax; p++) {
709: PetscInt parent, ancestor = p;
711: if (canonLabel) {
712: PetscInt canon;
714: PetscCall(DMLabelGetValue(canonLabel, p, &canon));
715: if (p != canon) continue;
716: }
717: PetscCall(DMPlexGetTreeParent(dm, p, &parent, NULL));
718: while (parent != ancestor) {
719: ancestor = parent;
720: PetscCall(DMPlexGetTreeParent(dm, ancestor, &parent, NULL));
721: }
722: if (ancestor != p) {
723: PetscInt j, closureSize, *closure = NULL, aOff;
725: PetscCall(PetscSectionGetOffset(aSec, p, &aOff));
727: PetscCall(DMPlexGetTransitiveClosure(dm, ancestor, PETSC_TRUE, &closureSize, &closure));
728: for (j = 0; j < closureSize; j++) anchors[aOff + j] = closure[2 * j];
729: PetscCall(DMPlexRestoreTransitiveClosure(dm, ancestor, PETSC_TRUE, &closureSize, &closure));
730: }
731: }
732: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size, anchors, PETSC_OWN_POINTER, &aIS));
733: {
734: PetscSection aSecNew = aSec;
735: IS aISNew = aIS;
737: PetscCall(PetscObjectReference((PetscObject)aSec));
738: PetscCall(PetscObjectReference((PetscObject)aIS));
739: while (aSecNew) {
740: PetscCall(PetscSectionDestroy(&aSec));
741: PetscCall(ISDestroy(&aIS));
742: aSec = aSecNew;
743: aIS = aISNew;
744: aSecNew = NULL;
745: aISNew = NULL;
746: PetscCall(AnchorsFlatten(aSec, aIS, &aSecNew, &aISNew));
747: }
748: }
749: PetscCall(DMPlexSetAnchors(dm, aSec, aIS));
750: PetscCall(PetscSectionDestroy(&aSec));
751: PetscCall(ISDestroy(&aIS));
752: PetscFunctionReturn(PETSC_SUCCESS);
753: }
755: static PetscErrorCode DMPlexGetTrueSupportSize(DM dm, PetscInt p, PetscInt *dof, PetscInt *numTrueSupp)
756: {
757: PetscFunctionBegin;
758: if (numTrueSupp[p] == -1) {
759: PetscInt i, alldof;
760: const PetscInt *supp;
761: PetscInt count = 0;
763: PetscCall(DMPlexGetSupportSize(dm, p, &alldof));
764: PetscCall(DMPlexGetSupport(dm, p, &supp));
765: for (i = 0; i < alldof; i++) {
766: PetscInt q = supp[i], numCones, j;
767: const PetscInt *cone;
769: PetscCall(DMPlexGetConeSize(dm, q, &numCones));
770: PetscCall(DMPlexGetCone(dm, q, &cone));
771: for (j = 0; j < numCones; j++) {
772: if (cone[j] == p) break;
773: }
774: if (j < numCones) count++;
775: }
776: numTrueSupp[p] = count;
777: }
778: *dof = numTrueSupp[p];
779: PetscFunctionReturn(PETSC_SUCCESS);
780: }
782: static PetscErrorCode DMPlexTreeExchangeSupports(DM dm)
783: {
784: DM_Plex *mesh = (DM_Plex *)dm->data;
785: PetscSection newSupportSection;
786: PetscInt newSize, *newSupports, pStart, pEnd, p, d, depth;
787: PetscInt *numTrueSupp;
788: PetscInt *offsets;
790: PetscFunctionBegin;
792: /* symmetrize the hierarchy */
793: PetscCall(DMPlexGetDepth(dm, &depth));
794: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)(mesh->supportSection)), &newSupportSection));
795: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
796: PetscCall(PetscSectionSetChart(newSupportSection, pStart, pEnd));
797: PetscCall(PetscCalloc1(pEnd, &offsets));
798: PetscCall(PetscMalloc1(pEnd, &numTrueSupp));
799: for (p = 0; p < pEnd; p++) numTrueSupp[p] = -1;
800: /* if a point is in the (true) support of q, it should be in the support of
801: * parent(q) */
802: for (d = 0; d <= depth; d++) {
803: PetscCall(DMPlexGetHeightStratum(dm, d, &pStart, &pEnd));
804: for (p = pStart; p < pEnd; ++p) {
805: PetscInt dof, q, qdof, parent;
807: PetscCall(DMPlexGetTrueSupportSize(dm, p, &dof, numTrueSupp));
808: PetscCall(PetscSectionAddDof(newSupportSection, p, dof));
809: q = p;
810: PetscCall(DMPlexGetTreeParent(dm, q, &parent, NULL));
811: while (parent != q && parent >= pStart && parent < pEnd) {
812: q = parent;
814: PetscCall(DMPlexGetTrueSupportSize(dm, q, &qdof, numTrueSupp));
815: PetscCall(PetscSectionAddDof(newSupportSection, p, qdof));
816: PetscCall(PetscSectionAddDof(newSupportSection, q, dof));
817: PetscCall(DMPlexGetTreeParent(dm, q, &parent, NULL));
818: }
819: }
820: }
821: PetscCall(PetscSectionSetUp(newSupportSection));
822: PetscCall(PetscSectionGetStorageSize(newSupportSection, &newSize));
823: PetscCall(PetscMalloc1(newSize, &newSupports));
824: for (d = 0; d <= depth; d++) {
825: PetscCall(DMPlexGetHeightStratum(dm, d, &pStart, &pEnd));
826: for (p = pStart; p < pEnd; p++) {
827: PetscInt dof, off, q, qdof, qoff, newDof, newOff, newqOff, i, parent;
829: PetscCall(PetscSectionGetDof(mesh->supportSection, p, &dof));
830: PetscCall(PetscSectionGetOffset(mesh->supportSection, p, &off));
831: PetscCall(PetscSectionGetDof(newSupportSection, p, &newDof));
832: PetscCall(PetscSectionGetOffset(newSupportSection, p, &newOff));
833: for (i = 0; i < dof; i++) {
834: PetscInt numCones, j;
835: const PetscInt *cone;
836: PetscInt q = mesh->supports[off + i];
838: PetscCall(DMPlexGetConeSize(dm, q, &numCones));
839: PetscCall(DMPlexGetCone(dm, q, &cone));
840: for (j = 0; j < numCones; j++) {
841: if (cone[j] == p) break;
842: }
843: if (j < numCones) newSupports[newOff + offsets[p]++] = q;
844: }
846: q = p;
847: PetscCall(DMPlexGetTreeParent(dm, q, &parent, NULL));
848: while (parent != q && parent >= pStart && parent < pEnd) {
849: q = parent;
850: PetscCall(PetscSectionGetDof(mesh->supportSection, q, &qdof));
851: PetscCall(PetscSectionGetOffset(mesh->supportSection, q, &qoff));
852: PetscCall(PetscSectionGetOffset(newSupportSection, q, &newqOff));
853: for (i = 0; i < qdof; i++) {
854: PetscInt numCones, j;
855: const PetscInt *cone;
856: PetscInt r = mesh->supports[qoff + i];
858: PetscCall(DMPlexGetConeSize(dm, r, &numCones));
859: PetscCall(DMPlexGetCone(dm, r, &cone));
860: for (j = 0; j < numCones; j++) {
861: if (cone[j] == q) break;
862: }
863: if (j < numCones) newSupports[newOff + offsets[p]++] = r;
864: }
865: for (i = 0; i < dof; i++) {
866: PetscInt numCones, j;
867: const PetscInt *cone;
868: PetscInt r = mesh->supports[off + i];
870: PetscCall(DMPlexGetConeSize(dm, r, &numCones));
871: PetscCall(DMPlexGetCone(dm, r, &cone));
872: for (j = 0; j < numCones; j++) {
873: if (cone[j] == p) break;
874: }
875: if (j < numCones) newSupports[newqOff + offsets[q]++] = r;
876: }
877: PetscCall(DMPlexGetTreeParent(dm, q, &parent, NULL));
878: }
879: }
880: }
881: PetscCall(PetscSectionDestroy(&mesh->supportSection));
882: mesh->supportSection = newSupportSection;
883: PetscCall(PetscFree(mesh->supports));
884: mesh->supports = newSupports;
885: PetscCall(PetscFree(offsets));
886: PetscCall(PetscFree(numTrueSupp));
888: PetscFunctionReturn(PETSC_SUCCESS);
889: }
891: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM, PetscSection, PetscSection, Mat);
892: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM, PetscSection, PetscSection, Mat);
894: static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical, PetscBool exchangeSupports)
895: {
896: DM_Plex *mesh = (DM_Plex *)dm->data;
897: DM refTree;
898: PetscInt size;
900: PetscFunctionBegin;
903: PetscCall(PetscObjectReference((PetscObject)parentSection));
904: PetscCall(PetscSectionDestroy(&mesh->parentSection));
905: mesh->parentSection = parentSection;
906: PetscCall(PetscSectionGetStorageSize(parentSection, &size));
907: if (parents != mesh->parents) {
908: PetscCall(PetscFree(mesh->parents));
909: PetscCall(PetscMalloc1(size, &mesh->parents));
910: PetscCall(PetscArraycpy(mesh->parents, parents, size));
911: }
912: if (childIDs != mesh->childIDs) {
913: PetscCall(PetscFree(mesh->childIDs));
914: PetscCall(PetscMalloc1(size, &mesh->childIDs));
915: PetscCall(PetscArraycpy(mesh->childIDs, childIDs, size));
916: }
917: PetscCall(DMPlexGetReferenceTree(dm, &refTree));
918: if (refTree) {
919: DMLabel canonLabel;
921: PetscCall(DMGetLabel(refTree, "canonical", &canonLabel));
922: if (canonLabel) {
923: PetscInt i;
925: for (i = 0; i < size; i++) {
926: PetscInt canon;
927: PetscCall(DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon));
928: if (canon >= 0) mesh->childIDs[i] = canon;
929: }
930: }
931: mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_FromReference;
932: } else {
933: mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_Direct;
934: }
935: PetscCall(DMPlexTreeSymmetrize(dm));
936: if (computeCanonical) {
937: PetscInt d, dim;
939: /* add the canonical label */
940: PetscCall(DMGetDimension(dm, &dim));
941: PetscCall(DMCreateLabel(dm, "canonical"));
942: for (d = 0; d <= dim; d++) {
943: PetscInt p, dStart, dEnd, canon = -1, cNumChildren;
944: const PetscInt *cChildren;
946: PetscCall(DMPlexGetDepthStratum(dm, d, &dStart, &dEnd));
947: for (p = dStart; p < dEnd; p++) {
948: PetscCall(DMPlexGetTreeChildren(dm, p, &cNumChildren, &cChildren));
949: if (cNumChildren) {
950: canon = p;
951: break;
952: }
953: }
954: if (canon == -1) continue;
955: for (p = dStart; p < dEnd; p++) {
956: PetscInt numChildren, i;
957: const PetscInt *children;
959: PetscCall(DMPlexGetTreeChildren(dm, p, &numChildren, &children));
960: if (numChildren) {
961: PetscCheck(numChildren == cNumChildren, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "All parent points in a stratum should have the same number of children: %" PetscInt_FMT " != %" PetscInt_FMT, numChildren, cNumChildren);
962: PetscCall(DMSetLabelValue(dm, "canonical", p, canon));
963: for (i = 0; i < numChildren; i++) PetscCall(DMSetLabelValue(dm, "canonical", children[i], cChildren[i]));
964: }
965: }
966: }
967: }
968: if (exchangeSupports) PetscCall(DMPlexTreeExchangeSupports(dm));
969: mesh->createanchors = DMPlexCreateAnchors_Tree;
970: /* reset anchors */
971: PetscCall(DMPlexSetAnchors(dm, NULL, NULL));
972: PetscFunctionReturn(PETSC_SUCCESS);
973: }
975: /*@
976: DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points. This routine also creates
977: the point-to-point constraints determined by the tree: a point is constrained to the points in the closure of its
978: tree root.
980: Collective
982: Input Parameters:
983: + dm - the `DMPLEX` object
984: . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
985: offset indexes the parent and childID list; the reference count of parentSection is incremented
986: . parents - a list of the point parents; copied, can be destroyed
987: - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
988: the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed
990: Level: intermediate
992: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetTree()`, `DMPlexSetReferenceTree()`, `DMPlexSetAnchors()`, `DMPlexGetTreeParent()`, `DMPlexGetTreeChildren()`
993: @*/
994: PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
995: {
996: PetscFunctionBegin;
997: PetscCall(DMPlexSetTree_Internal(dm, parentSection, parents, childIDs, PETSC_FALSE, PETSC_TRUE));
998: PetscFunctionReturn(PETSC_SUCCESS);
999: }
1001: /*@
1002: DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points.
1003: Collective
1005: Input Parameter:
1006: . dm - the `DMPLEX` object
1008: Output Parameters:
1009: + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
1010: offset indexes the parent and childID list
1011: . parents - a list of the point parents
1012: . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
1013: the child corresponds to the point in the reference tree with index childID
1014: . childSection - the inverse of the parent section
1015: - children - a list of the point children
1017: Level: intermediate
1019: .seealso: [](ch_unstructured), `DM`, `DMPLEX`,`DMPlexSetTree()`, `DMPlexSetReferenceTree()`, `DMPlexSetAnchors()`, `DMPlexGetTreeParent()`, `DMPlexGetTreeChildren()`
1020: @*/
1021: PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[])
1022: {
1023: DM_Plex *mesh = (DM_Plex *)dm->data;
1025: PetscFunctionBegin;
1027: if (parentSection) *parentSection = mesh->parentSection;
1028: if (parents) *parents = mesh->parents;
1029: if (childIDs) *childIDs = mesh->childIDs;
1030: if (childSection) *childSection = mesh->childSection;
1031: if (children) *children = mesh->children;
1032: PetscFunctionReturn(PETSC_SUCCESS);
1033: }
1035: /*@
1036: DMPlexGetTreeParent - get the parent of a point in the tree describing the point hierarchy (not the DAG)
1038: Input Parameters:
1039: + dm - the `DMPLEX` object
1040: - point - the query point
1042: Output Parameters:
1043: + parent - if not `NULL`, set to the parent of the point, or the point itself if the point does not have a parent
1044: - childID - if not `NULL`, set to the child ID of the point with respect to its parent, or 0 if the point
1045: does not have a parent
1047: Level: intermediate
1049: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetTree()`, `DMPlexGetTree()`, `DMPlexGetTreeChildren()`
1050: @*/
1051: PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID)
1052: {
1053: DM_Plex *mesh = (DM_Plex *)dm->data;
1054: PetscSection pSec;
1056: PetscFunctionBegin;
1058: pSec = mesh->parentSection;
1059: if (pSec && point >= pSec->pStart && point < pSec->pEnd) {
1060: PetscInt dof;
1062: PetscCall(PetscSectionGetDof(pSec, point, &dof));
1063: if (dof) {
1064: PetscInt off;
1066: PetscCall(PetscSectionGetOffset(pSec, point, &off));
1067: if (parent) *parent = mesh->parents[off];
1068: if (childID) *childID = mesh->childIDs[off];
1069: PetscFunctionReturn(PETSC_SUCCESS);
1070: }
1071: }
1072: if (parent) *parent = point;
1073: if (childID) *childID = 0;
1074: PetscFunctionReturn(PETSC_SUCCESS);
1075: }
1077: /*@C
1078: DMPlexGetTreeChildren - get the children of a point in the tree describing the point hierarchy (not the DAG)
1080: Input Parameters:
1081: + dm - the `DMPLEX` object
1082: - point - the query point
1084: Output Parameters:
1085: + numChildren - if not `NULL`, set to the number of children
1086: - children - if not `NULL`, set to a list children, or set to `NULL` if the point has no children
1088: Level: intermediate
1090: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetTree()`, `DMPlexGetTree()`, `DMPlexGetTreeParent()`
1091: @*/
1092: PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[])
1093: {
1094: DM_Plex *mesh = (DM_Plex *)dm->data;
1095: PetscSection childSec;
1096: PetscInt dof = 0;
1098: PetscFunctionBegin;
1100: childSec = mesh->childSection;
1101: if (childSec && point >= childSec->pStart && point < childSec->pEnd) PetscCall(PetscSectionGetDof(childSec, point, &dof));
1102: if (numChildren) *numChildren = dof;
1103: if (children) {
1104: if (dof) {
1105: PetscInt off;
1107: PetscCall(PetscSectionGetOffset(childSec, point, &off));
1108: *children = &mesh->children[off];
1109: } else {
1110: *children = NULL;
1111: }
1112: }
1113: PetscFunctionReturn(PETSC_SUCCESS);
1114: }
1116: static PetscErrorCode EvaluateBasis(PetscSpace space, PetscInt nBasis, PetscInt nFunctionals, PetscInt nComps, PetscInt nPoints, const PetscInt *pointsPerFn, const PetscReal *points, const PetscReal *weights, PetscReal *work, Mat basisAtPoints)
1117: {
1118: PetscInt f, b, p, c, offset, qPoints;
1120: PetscFunctionBegin;
1121: PetscCall(PetscSpaceEvaluate(space, nPoints, points, work, NULL, NULL));
1122: for (f = 0, offset = 0; f < nFunctionals; f++) {
1123: qPoints = pointsPerFn[f];
1124: for (b = 0; b < nBasis; b++) {
1125: PetscScalar val = 0.;
1127: for (p = 0; p < qPoints; p++) {
1128: for (c = 0; c < nComps; c++) val += work[((offset + p) * nBasis + b) * nComps + c] * weights[(offset + p) * nComps + c];
1129: }
1130: PetscCall(MatSetValue(basisAtPoints, b, f, val, INSERT_VALUES));
1131: }
1132: offset += qPoints;
1133: }
1134: PetscCall(MatAssemblyBegin(basisAtPoints, MAT_FINAL_ASSEMBLY));
1135: PetscCall(MatAssemblyEnd(basisAtPoints, MAT_FINAL_ASSEMBLY));
1136: PetscFunctionReturn(PETSC_SUCCESS);
1137: }
1139: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM dm, PetscSection section, PetscSection cSec, Mat cMat)
1140: {
1141: PetscDS ds;
1142: PetscInt spdim;
1143: PetscInt numFields, f, c, cStart, cEnd, pStart, pEnd, conStart, conEnd;
1144: const PetscInt *anchors;
1145: PetscSection aSec;
1146: PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJparent, detJ, detJparent;
1147: IS aIS;
1149: PetscFunctionBegin;
1150: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1151: PetscCall(DMGetDS(dm, &ds));
1152: PetscCall(PetscDSGetNumFields(ds, &numFields));
1153: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1154: PetscCall(DMPlexGetAnchors(dm, &aSec, &aIS));
1155: PetscCall(ISGetIndices(aIS, &anchors));
1156: PetscCall(PetscSectionGetChart(cSec, &conStart, &conEnd));
1157: PetscCall(DMGetDimension(dm, &spdim));
1158: PetscCall(PetscMalloc6(spdim, &v0, spdim, &v0parent, spdim, &vtmp, spdim * spdim, &J, spdim * spdim, &Jparent, spdim * spdim, &invJparent));
1160: for (f = 0; f < numFields; f++) {
1161: PetscObject disc;
1162: PetscClassId id;
1163: PetscSpace bspace;
1164: PetscDualSpace dspace;
1165: PetscInt i, j, k, nPoints, Nc, offset;
1166: PetscInt fSize, maxDof;
1167: PetscReal *weights, *pointsRef, *pointsReal, *work;
1168: PetscScalar *scwork;
1169: const PetscScalar *X;
1170: PetscInt *sizes, *workIndRow, *workIndCol;
1171: Mat Amat, Bmat, Xmat;
1172: const PetscInt *numDof = NULL;
1173: const PetscInt ***perms = NULL;
1174: const PetscScalar ***flips = NULL;
1176: PetscCall(PetscDSGetDiscretization(ds, f, &disc));
1177: PetscCall(PetscObjectGetClassId(disc, &id));
1178: if (id == PETSCFE_CLASSID) {
1179: PetscFE fe = (PetscFE)disc;
1181: PetscCall(PetscFEGetBasisSpace(fe, &bspace));
1182: PetscCall(PetscFEGetDualSpace(fe, &dspace));
1183: PetscCall(PetscDualSpaceGetDimension(dspace, &fSize));
1184: PetscCall(PetscFEGetNumComponents(fe, &Nc));
1185: } else if (id == PETSCFV_CLASSID) {
1186: PetscFV fv = (PetscFV)disc;
1188: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1189: PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)fv), &bspace));
1190: PetscCall(PetscSpaceSetType(bspace, PETSCSPACEPOLYNOMIAL));
1191: PetscCall(PetscSpaceSetDegree(bspace, 0, PETSC_DETERMINE));
1192: PetscCall(PetscSpaceSetNumComponents(bspace, Nc));
1193: PetscCall(PetscSpaceSetNumVariables(bspace, spdim));
1194: PetscCall(PetscSpaceSetUp(bspace));
1195: PetscCall(PetscFVGetDualSpace(fv, &dspace));
1196: PetscCall(PetscDualSpaceGetDimension(dspace, &fSize));
1197: } else SETERRQ(PetscObjectComm(disc), PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id);
1198: PetscCall(PetscDualSpaceGetNumDof(dspace, &numDof));
1199: for (i = 0, maxDof = 0; i <= spdim; i++) maxDof = PetscMax(maxDof, numDof[i]);
1200: PetscCall(PetscDualSpaceGetSymmetries(dspace, &perms, &flips));
1202: PetscCall(MatCreate(PETSC_COMM_SELF, &Amat));
1203: PetscCall(MatSetSizes(Amat, fSize, fSize, fSize, fSize));
1204: PetscCall(MatSetType(Amat, MATSEQDENSE));
1205: PetscCall(MatSetUp(Amat));
1206: PetscCall(MatDuplicate(Amat, MAT_DO_NOT_COPY_VALUES, &Bmat));
1207: PetscCall(MatDuplicate(Amat, MAT_DO_NOT_COPY_VALUES, &Xmat));
1208: nPoints = 0;
1209: for (i = 0; i < fSize; i++) {
1210: PetscInt qPoints, thisNc;
1211: PetscQuadrature quad;
1213: PetscCall(PetscDualSpaceGetFunctional(dspace, i, &quad));
1214: PetscCall(PetscQuadratureGetData(quad, NULL, &thisNc, &qPoints, NULL, NULL));
1215: PetscCheck(thisNc == Nc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Functional dim %" PetscInt_FMT " does not much basis dim %" PetscInt_FMT, thisNc, Nc);
1216: nPoints += qPoints;
1217: }
1218: PetscCall(PetscMalloc7(fSize, &sizes, nPoints * Nc, &weights, spdim * nPoints, &pointsRef, spdim * nPoints, &pointsReal, nPoints * fSize * Nc, &work, maxDof, &workIndRow, maxDof, &workIndCol));
1219: PetscCall(PetscMalloc1(maxDof * maxDof, &scwork));
1220: offset = 0;
1221: for (i = 0; i < fSize; i++) {
1222: PetscInt qPoints;
1223: const PetscReal *p, *w;
1224: PetscQuadrature quad;
1226: PetscCall(PetscDualSpaceGetFunctional(dspace, i, &quad));
1227: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &qPoints, &p, &w));
1228: PetscCall(PetscArraycpy(weights + Nc * offset, w, Nc * qPoints));
1229: PetscCall(PetscArraycpy(pointsRef + spdim * offset, p, spdim * qPoints));
1230: sizes[i] = qPoints;
1231: offset += qPoints;
1232: }
1233: PetscCall(EvaluateBasis(bspace, fSize, fSize, Nc, nPoints, sizes, pointsRef, weights, work, Amat));
1234: PetscCall(MatLUFactor(Amat, NULL, NULL, NULL));
1235: for (c = cStart; c < cEnd; c++) {
1236: PetscInt parent;
1237: PetscInt closureSize, closureSizeP, *closure = NULL, *closureP = NULL;
1238: PetscInt *childOffsets, *parentOffsets;
1240: PetscCall(DMPlexGetTreeParent(dm, c, &parent, NULL));
1241: if (parent == c) continue;
1242: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1243: for (i = 0; i < closureSize; i++) {
1244: PetscInt p = closure[2 * i];
1245: PetscInt conDof;
1247: if (p < conStart || p >= conEnd) continue;
1248: if (numFields) {
1249: PetscCall(PetscSectionGetFieldDof(cSec, p, f, &conDof));
1250: } else {
1251: PetscCall(PetscSectionGetDof(cSec, p, &conDof));
1252: }
1253: if (conDof) break;
1254: }
1255: if (i == closureSize) {
1256: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1257: continue;
1258: }
1260: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, NULL, &detJ));
1261: PetscCall(DMPlexComputeCellGeometryFEM(dm, parent, NULL, v0parent, Jparent, invJparent, &detJparent));
1262: for (i = 0; i < nPoints; i++) {
1263: const PetscReal xi0[3] = {-1., -1., -1.};
1265: CoordinatesRefToReal(spdim, spdim, xi0, v0, J, &pointsRef[i * spdim], vtmp);
1266: CoordinatesRealToRef(spdim, spdim, xi0, v0parent, invJparent, vtmp, &pointsReal[i * spdim]);
1267: }
1268: PetscCall(EvaluateBasis(bspace, fSize, fSize, Nc, nPoints, sizes, pointsReal, weights, work, Bmat));
1269: PetscCall(MatMatSolve(Amat, Bmat, Xmat));
1270: PetscCall(MatDenseGetArrayRead(Xmat, &X));
1271: PetscCall(DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSizeP, &closureP));
1272: PetscCall(PetscMalloc2(closureSize + 1, &childOffsets, closureSizeP + 1, &parentOffsets));
1273: childOffsets[0] = 0;
1274: for (i = 0; i < closureSize; i++) {
1275: PetscInt p = closure[2 * i];
1276: PetscInt dof;
1278: if (numFields) {
1279: PetscCall(PetscSectionGetFieldDof(section, p, f, &dof));
1280: } else {
1281: PetscCall(PetscSectionGetDof(section, p, &dof));
1282: }
1283: childOffsets[i + 1] = childOffsets[i] + dof;
1284: }
1285: parentOffsets[0] = 0;
1286: for (i = 0; i < closureSizeP; i++) {
1287: PetscInt p = closureP[2 * i];
1288: PetscInt dof;
1290: if (numFields) {
1291: PetscCall(PetscSectionGetFieldDof(section, p, f, &dof));
1292: } else {
1293: PetscCall(PetscSectionGetDof(section, p, &dof));
1294: }
1295: parentOffsets[i + 1] = parentOffsets[i] + dof;
1296: }
1297: for (i = 0; i < closureSize; i++) {
1298: PetscInt conDof, conOff, aDof, aOff, nWork;
1299: PetscInt p = closure[2 * i];
1300: PetscInt o = closure[2 * i + 1];
1301: const PetscInt *perm;
1302: const PetscScalar *flip;
1304: if (p < conStart || p >= conEnd) continue;
1305: if (numFields) {
1306: PetscCall(PetscSectionGetFieldDof(cSec, p, f, &conDof));
1307: PetscCall(PetscSectionGetFieldOffset(cSec, p, f, &conOff));
1308: } else {
1309: PetscCall(PetscSectionGetDof(cSec, p, &conDof));
1310: PetscCall(PetscSectionGetOffset(cSec, p, &conOff));
1311: }
1312: if (!conDof) continue;
1313: perm = (perms && perms[i]) ? perms[i][o] : NULL;
1314: flip = (flips && flips[i]) ? flips[i][o] : NULL;
1315: PetscCall(PetscSectionGetDof(aSec, p, &aDof));
1316: PetscCall(PetscSectionGetOffset(aSec, p, &aOff));
1317: nWork = childOffsets[i + 1] - childOffsets[i];
1318: for (k = 0; k < aDof; k++) {
1319: PetscInt a = anchors[aOff + k];
1320: PetscInt aSecDof, aSecOff;
1322: if (numFields) {
1323: PetscCall(PetscSectionGetFieldDof(section, a, f, &aSecDof));
1324: PetscCall(PetscSectionGetFieldOffset(section, a, f, &aSecOff));
1325: } else {
1326: PetscCall(PetscSectionGetDof(section, a, &aSecDof));
1327: PetscCall(PetscSectionGetOffset(section, a, &aSecOff));
1328: }
1329: if (!aSecDof) continue;
1331: for (j = 0; j < closureSizeP; j++) {
1332: PetscInt q = closureP[2 * j];
1333: PetscInt oq = closureP[2 * j + 1];
1335: if (q == a) {
1336: PetscInt r, s, nWorkP;
1337: const PetscInt *permP;
1338: const PetscScalar *flipP;
1340: permP = (perms && perms[j]) ? perms[j][oq] : NULL;
1341: flipP = (flips && flips[j]) ? flips[j][oq] : NULL;
1342: nWorkP = parentOffsets[j + 1] - parentOffsets[j];
1343: /* get a copy of the child-to-anchor portion of the matrix, and transpose so that rows correspond to the
1344: * child and columns correspond to the anchor: BUT the maxrix returned by MatDenseGetArrayRead() is
1345: * column-major, so transpose-transpose = do nothing */
1346: for (r = 0; r < nWork; r++) {
1347: for (s = 0; s < nWorkP; s++) scwork[r * nWorkP + s] = X[fSize * (r + childOffsets[i]) + (s + parentOffsets[j])];
1348: }
1349: for (r = 0; r < nWork; r++) workIndRow[perm ? perm[r] : r] = conOff + r;
1350: for (s = 0; s < nWorkP; s++) workIndCol[permP ? permP[s] : s] = aSecOff + s;
1351: if (flip) {
1352: for (r = 0; r < nWork; r++) {
1353: for (s = 0; s < nWorkP; s++) scwork[r * nWorkP + s] *= flip[r];
1354: }
1355: }
1356: if (flipP) {
1357: for (r = 0; r < nWork; r++) {
1358: for (s = 0; s < nWorkP; s++) scwork[r * nWorkP + s] *= flipP[s];
1359: }
1360: }
1361: PetscCall(MatSetValues(cMat, nWork, workIndRow, nWorkP, workIndCol, scwork, INSERT_VALUES));
1362: break;
1363: }
1364: }
1365: }
1366: }
1367: PetscCall(MatDenseRestoreArrayRead(Xmat, &X));
1368: PetscCall(PetscFree2(childOffsets, parentOffsets));
1369: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1370: PetscCall(DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSizeP, &closureP));
1371: }
1372: PetscCall(MatDestroy(&Amat));
1373: PetscCall(MatDestroy(&Bmat));
1374: PetscCall(MatDestroy(&Xmat));
1375: PetscCall(PetscFree(scwork));
1376: PetscCall(PetscFree7(sizes, weights, pointsRef, pointsReal, work, workIndRow, workIndCol));
1377: if (id == PETSCFV_CLASSID) PetscCall(PetscSpaceDestroy(&bspace));
1378: }
1379: PetscCall(MatAssemblyBegin(cMat, MAT_FINAL_ASSEMBLY));
1380: PetscCall(MatAssemblyEnd(cMat, MAT_FINAL_ASSEMBLY));
1381: PetscCall(PetscFree6(v0, v0parent, vtmp, J, Jparent, invJparent));
1382: PetscCall(ISRestoreIndices(aIS, &anchors));
1384: PetscFunctionReturn(PETSC_SUCCESS);
1385: }
1387: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1388: {
1389: Mat refCmat;
1390: PetscDS ds;
1391: PetscInt numFields, maxFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, **refPointFieldN;
1392: PetscScalar ***refPointFieldMats;
1393: PetscSection refConSec, refAnSec, refSection;
1394: IS refAnIS;
1395: const PetscInt *refAnchors;
1396: const PetscInt **perms;
1397: const PetscScalar **flips;
1399: PetscFunctionBegin;
1400: PetscCall(DMGetDS(refTree, &ds));
1401: PetscCall(PetscDSGetNumFields(ds, &numFields));
1402: maxFields = PetscMax(1, numFields);
1403: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, &refCmat, NULL));
1404: PetscCall(DMPlexGetAnchors(refTree, &refAnSec, &refAnIS));
1405: PetscCall(ISGetIndices(refAnIS, &refAnchors));
1406: PetscCall(DMGetLocalSection(refTree, &refSection));
1407: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
1408: PetscCall(PetscMalloc1(pRefEnd - pRefStart, &refPointFieldMats));
1409: PetscCall(PetscMalloc1(pRefEnd - pRefStart, &refPointFieldN));
1410: PetscCall(PetscSectionGetMaxDof(refConSec, &maxDof));
1411: PetscCall(PetscSectionGetMaxDof(refAnSec, &maxAnDof));
1412: PetscCall(PetscMalloc1(maxDof, &rows));
1413: PetscCall(PetscMalloc1(maxDof * maxAnDof, &cols));
1414: for (p = pRefStart; p < pRefEnd; p++) {
1415: PetscInt parent, closureSize, *closure = NULL, pDof;
1417: PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL));
1418: PetscCall(PetscSectionGetDof(refConSec, p, &pDof));
1419: if (!pDof || parent == p) continue;
1421: PetscCall(PetscMalloc1(maxFields, &refPointFieldMats[p - pRefStart]));
1422: PetscCall(PetscCalloc1(maxFields, &refPointFieldN[p - pRefStart]));
1423: PetscCall(DMPlexGetTransitiveClosure(refTree, parent, PETSC_TRUE, &closureSize, &closure));
1424: for (f = 0; f < maxFields; f++) {
1425: PetscInt cDof, cOff, numCols, r, i;
1427: if (f < numFields) {
1428: PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof));
1429: PetscCall(PetscSectionGetFieldOffset(refConSec, p, f, &cOff));
1430: PetscCall(PetscSectionGetFieldPointSyms(refSection, f, closureSize, closure, &perms, &flips));
1431: } else {
1432: PetscCall(PetscSectionGetDof(refConSec, p, &cDof));
1433: PetscCall(PetscSectionGetOffset(refConSec, p, &cOff));
1434: PetscCall(PetscSectionGetPointSyms(refSection, closureSize, closure, &perms, &flips));
1435: }
1437: for (r = 0; r < cDof; r++) rows[r] = cOff + r;
1438: numCols = 0;
1439: for (i = 0; i < closureSize; i++) {
1440: PetscInt q = closure[2 * i];
1441: PetscInt aDof, aOff, j;
1442: const PetscInt *perm = perms ? perms[i] : NULL;
1444: if (numFields) {
1445: PetscCall(PetscSectionGetFieldDof(refSection, q, f, &aDof));
1446: PetscCall(PetscSectionGetFieldOffset(refSection, q, f, &aOff));
1447: } else {
1448: PetscCall(PetscSectionGetDof(refSection, q, &aDof));
1449: PetscCall(PetscSectionGetOffset(refSection, q, &aOff));
1450: }
1452: for (j = 0; j < aDof; j++) cols[numCols++] = aOff + (perm ? perm[j] : j);
1453: }
1454: refPointFieldN[p - pRefStart][f] = numCols;
1455: PetscCall(PetscMalloc1(cDof * numCols, &refPointFieldMats[p - pRefStart][f]));
1456: PetscCall(MatGetValues(refCmat, cDof, rows, numCols, cols, refPointFieldMats[p - pRefStart][f]));
1457: if (flips) {
1458: PetscInt colOff = 0;
1460: for (i = 0; i < closureSize; i++) {
1461: PetscInt q = closure[2 * i];
1462: PetscInt aDof, aOff, j;
1463: const PetscScalar *flip = flips ? flips[i] : NULL;
1465: if (numFields) {
1466: PetscCall(PetscSectionGetFieldDof(refSection, q, f, &aDof));
1467: PetscCall(PetscSectionGetFieldOffset(refSection, q, f, &aOff));
1468: } else {
1469: PetscCall(PetscSectionGetDof(refSection, q, &aDof));
1470: PetscCall(PetscSectionGetOffset(refSection, q, &aOff));
1471: }
1472: if (flip) {
1473: PetscInt k;
1474: for (k = 0; k < cDof; k++) {
1475: for (j = 0; j < aDof; j++) refPointFieldMats[p - pRefStart][f][k * numCols + colOff + j] *= flip[j];
1476: }
1477: }
1478: colOff += aDof;
1479: }
1480: }
1481: if (numFields) {
1482: PetscCall(PetscSectionRestoreFieldPointSyms(refSection, f, closureSize, closure, &perms, &flips));
1483: } else {
1484: PetscCall(PetscSectionRestorePointSyms(refSection, closureSize, closure, &perms, &flips));
1485: }
1486: }
1487: PetscCall(DMPlexRestoreTransitiveClosure(refTree, parent, PETSC_TRUE, &closureSize, &closure));
1488: }
1489: *childrenMats = refPointFieldMats;
1490: *childrenN = refPointFieldN;
1491: PetscCall(ISRestoreIndices(refAnIS, &refAnchors));
1492: PetscCall(PetscFree(rows));
1493: PetscCall(PetscFree(cols));
1494: PetscFunctionReturn(PETSC_SUCCESS);
1495: }
1497: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1498: {
1499: PetscDS ds;
1500: PetscInt **refPointFieldN;
1501: PetscScalar ***refPointFieldMats;
1502: PetscInt numFields, maxFields, pRefStart, pRefEnd, p, f;
1503: PetscSection refConSec;
1505: PetscFunctionBegin;
1506: refPointFieldN = *childrenN;
1507: *childrenN = NULL;
1508: refPointFieldMats = *childrenMats;
1509: *childrenMats = NULL;
1510: PetscCall(DMGetDS(refTree, &ds));
1511: PetscCall(PetscDSGetNumFields(ds, &numFields));
1512: maxFields = PetscMax(1, numFields);
1513: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
1514: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
1515: for (p = pRefStart; p < pRefEnd; p++) {
1516: PetscInt parent, pDof;
1518: PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL));
1519: PetscCall(PetscSectionGetDof(refConSec, p, &pDof));
1520: if (!pDof || parent == p) continue;
1522: for (f = 0; f < maxFields; f++) {
1523: PetscInt cDof;
1525: if (numFields) {
1526: PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof));
1527: } else {
1528: PetscCall(PetscSectionGetDof(refConSec, p, &cDof));
1529: }
1531: PetscCall(PetscFree(refPointFieldMats[p - pRefStart][f]));
1532: }
1533: PetscCall(PetscFree(refPointFieldMats[p - pRefStart]));
1534: PetscCall(PetscFree(refPointFieldN[p - pRefStart]));
1535: }
1536: PetscCall(PetscFree(refPointFieldMats));
1537: PetscCall(PetscFree(refPointFieldN));
1538: PetscFunctionReturn(PETSC_SUCCESS);
1539: }
1541: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm, PetscSection section, PetscSection conSec, Mat cMat)
1542: {
1543: DM refTree;
1544: PetscDS ds;
1545: Mat refCmat;
1546: PetscInt numFields, maxFields, f, pRefStart, pRefEnd, p, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN;
1547: PetscScalar ***refPointFieldMats, *pointWork;
1548: PetscSection refConSec, refAnSec, anSec;
1549: IS refAnIS, anIS;
1550: const PetscInt *anchors;
1552: PetscFunctionBegin;
1554: PetscCall(DMGetDS(dm, &ds));
1555: PetscCall(PetscDSGetNumFields(ds, &numFields));
1556: maxFields = PetscMax(1, numFields);
1557: PetscCall(DMPlexGetReferenceTree(dm, &refTree));
1558: PetscCall(DMCopyDisc(dm, refTree));
1559: PetscCall(DMSetLocalSection(refTree, NULL));
1560: PetscCall(DMSetDefaultConstraints(refTree, NULL, NULL, NULL));
1561: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, &refCmat, NULL));
1562: PetscCall(DMPlexGetAnchors(refTree, &refAnSec, &refAnIS));
1563: PetscCall(DMPlexGetAnchors(dm, &anSec, &anIS));
1564: PetscCall(ISGetIndices(anIS, &anchors));
1565: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
1566: PetscCall(PetscSectionGetChart(conSec, &conStart, &conEnd));
1567: PetscCall(PetscSectionGetMaxDof(refConSec, &maxDof));
1568: PetscCall(PetscSectionGetMaxDof(refAnSec, &maxAnDof));
1569: PetscCall(PetscMalloc1(maxDof * maxDof * maxAnDof, &pointWork));
1571: /* step 1: get submats for every constrained point in the reference tree */
1572: PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
1574: /* step 2: compute the preorder */
1575: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1576: PetscCall(PetscMalloc2(pEnd - pStart, &perm, pEnd - pStart, &iperm));
1577: for (p = pStart; p < pEnd; p++) {
1578: perm[p - pStart] = p;
1579: iperm[p - pStart] = p - pStart;
1580: }
1581: for (p = 0; p < pEnd - pStart;) {
1582: PetscInt point = perm[p];
1583: PetscInt parent;
1585: PetscCall(DMPlexGetTreeParent(dm, point, &parent, NULL));
1586: if (parent == point) {
1587: p++;
1588: } else {
1589: PetscInt size, closureSize, *closure = NULL, i;
1591: PetscCall(DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
1592: for (i = 0; i < closureSize; i++) {
1593: PetscInt q = closure[2 * i];
1594: if (iperm[q - pStart] > iperm[point - pStart]) {
1595: /* swap */
1596: perm[p] = q;
1597: perm[iperm[q - pStart]] = point;
1598: iperm[point - pStart] = iperm[q - pStart];
1599: iperm[q - pStart] = p;
1600: break;
1601: }
1602: }
1603: size = closureSize;
1604: PetscCall(DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
1605: if (i == size) p++;
1606: }
1607: }
1609: /* step 3: fill the constraint matrix */
1610: /* we are going to use a preorder progressive fill strategy. Mat doesn't
1611: * allow progressive fill without assembly, so we are going to set up the
1612: * values outside of the Mat first.
1613: */
1614: {
1615: PetscInt nRows, row, nnz;
1616: PetscBool done;
1617: PetscInt secStart, secEnd;
1618: const PetscInt *ia, *ja;
1619: PetscScalar *vals;
1621: PetscCall(PetscSectionGetChart(section, &secStart, &secEnd));
1622: PetscCall(MatGetRowIJ(cMat, 0, PETSC_FALSE, PETSC_FALSE, &nRows, &ia, &ja, &done));
1623: PetscCheck(done, PetscObjectComm((PetscObject)cMat), PETSC_ERR_PLIB, "Could not get RowIJ of constraint matrix");
1624: nnz = ia[nRows];
1625: /* malloc and then zero rows right before we fill them: this way valgrind
1626: * can tell if we are doing progressive fill in the wrong order */
1627: PetscCall(PetscMalloc1(nnz, &vals));
1628: for (p = 0; p < pEnd - pStart; p++) {
1629: PetscInt parent, childid, closureSize, *closure = NULL;
1630: PetscInt point = perm[p], pointDof;
1632: PetscCall(DMPlexGetTreeParent(dm, point, &parent, &childid));
1633: if ((point < conStart) || (point >= conEnd) || (parent == point)) continue;
1634: PetscCall(PetscSectionGetDof(conSec, point, &pointDof));
1635: if (!pointDof) continue;
1636: PetscCall(DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
1637: for (f = 0; f < maxFields; f++) {
1638: PetscInt cDof, cOff, numCols, numFillCols, i, r, matOffset, offset;
1639: PetscScalar *pointMat;
1640: const PetscInt **perms;
1641: const PetscScalar **flips;
1643: if (numFields) {
1644: PetscCall(PetscSectionGetFieldDof(conSec, point, f, &cDof));
1645: PetscCall(PetscSectionGetFieldOffset(conSec, point, f, &cOff));
1646: } else {
1647: PetscCall(PetscSectionGetDof(conSec, point, &cDof));
1648: PetscCall(PetscSectionGetOffset(conSec, point, &cOff));
1649: }
1650: if (!cDof) continue;
1651: if (numFields) PetscCall(PetscSectionGetFieldPointSyms(section, f, closureSize, closure, &perms, &flips));
1652: else PetscCall(PetscSectionGetPointSyms(section, closureSize, closure, &perms, &flips));
1654: /* make sure that every row for this point is the same size */
1655: if (PetscDefined(USE_DEBUG)) {
1656: for (r = 0; r < cDof; r++) {
1657: if (cDof > 1 && r) {
1658: PetscCheck((ia[cOff + r + 1] - ia[cOff + r]) == (ia[cOff + r] - ia[cOff + r - 1]), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two point rows have different nnz: %" PetscInt_FMT " vs. %" PetscInt_FMT, (ia[cOff + r + 1] - ia[cOff + r]), (ia[cOff + r] - ia[cOff + r - 1]));
1659: }
1660: }
1661: }
1662: /* zero rows */
1663: for (i = ia[cOff]; i < ia[cOff + cDof]; i++) vals[i] = 0.;
1664: matOffset = ia[cOff];
1665: numFillCols = ia[cOff + 1] - matOffset;
1666: pointMat = refPointFieldMats[childid - pRefStart][f];
1667: numCols = refPointFieldN[childid - pRefStart][f];
1668: offset = 0;
1669: for (i = 0; i < closureSize; i++) {
1670: PetscInt q = closure[2 * i];
1671: PetscInt aDof, aOff, j, k, qConDof, qConOff;
1672: const PetscInt *perm = perms ? perms[i] : NULL;
1673: const PetscScalar *flip = flips ? flips[i] : NULL;
1675: qConDof = qConOff = 0;
1676: if (q < secStart || q >= secEnd) continue;
1677: if (numFields) {
1678: PetscCall(PetscSectionGetFieldDof(section, q, f, &aDof));
1679: PetscCall(PetscSectionGetFieldOffset(section, q, f, &aOff));
1680: if (q >= conStart && q < conEnd) {
1681: PetscCall(PetscSectionGetFieldDof(conSec, q, f, &qConDof));
1682: PetscCall(PetscSectionGetFieldOffset(conSec, q, f, &qConOff));
1683: }
1684: } else {
1685: PetscCall(PetscSectionGetDof(section, q, &aDof));
1686: PetscCall(PetscSectionGetOffset(section, q, &aOff));
1687: if (q >= conStart && q < conEnd) {
1688: PetscCall(PetscSectionGetDof(conSec, q, &qConDof));
1689: PetscCall(PetscSectionGetOffset(conSec, q, &qConOff));
1690: }
1691: }
1692: if (!aDof) continue;
1693: if (qConDof) {
1694: /* this point has anchors: its rows of the matrix should already
1695: * be filled, thanks to preordering */
1696: /* first multiply into pointWork, then set in matrix */
1697: PetscInt aMatOffset = ia[qConOff];
1698: PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset;
1699: for (r = 0; r < cDof; r++) {
1700: for (j = 0; j < aNumFillCols; j++) {
1701: PetscScalar inVal = 0;
1702: for (k = 0; k < aDof; k++) {
1703: PetscInt col = perm ? perm[k] : k;
1705: inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j] * (flip ? flip[col] : 1.);
1706: }
1707: pointWork[r * aNumFillCols + j] = inVal;
1708: }
1709: }
1710: /* assume that the columns are sorted, spend less time searching */
1711: for (j = 0, k = 0; j < aNumFillCols; j++) {
1712: PetscInt col = ja[aMatOffset + j];
1713: for (; k < numFillCols; k++) {
1714: if (ja[matOffset + k] == col) break;
1715: }
1716: PetscCheck(k != numFillCols, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No nonzero space for (%" PetscInt_FMT ", %" PetscInt_FMT ")", cOff, col);
1717: for (r = 0; r < cDof; r++) vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j];
1718: }
1719: } else {
1720: /* find where to put this portion of pointMat into the matrix */
1721: for (k = 0; k < numFillCols; k++) {
1722: if (ja[matOffset + k] == aOff) break;
1723: }
1724: PetscCheck(k != numFillCols, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No nonzero space for (%" PetscInt_FMT ", %" PetscInt_FMT ")", cOff, aOff);
1725: for (r = 0; r < cDof; r++) {
1726: for (j = 0; j < aDof; j++) {
1727: PetscInt col = perm ? perm[j] : j;
1729: vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + j] * (flip ? flip[col] : 1.);
1730: }
1731: }
1732: }
1733: offset += aDof;
1734: }
1735: if (numFields) {
1736: PetscCall(PetscSectionRestoreFieldPointSyms(section, f, closureSize, closure, &perms, &flips));
1737: } else {
1738: PetscCall(PetscSectionRestorePointSyms(section, closureSize, closure, &perms, &flips));
1739: }
1740: }
1741: PetscCall(DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
1742: }
1743: for (row = 0; row < nRows; row++) PetscCall(MatSetValues(cMat, 1, &row, ia[row + 1] - ia[row], &ja[ia[row]], &vals[ia[row]], INSERT_VALUES));
1744: PetscCall(MatRestoreRowIJ(cMat, 0, PETSC_FALSE, PETSC_FALSE, &nRows, &ia, &ja, &done));
1745: PetscCheck(done, PetscObjectComm((PetscObject)cMat), PETSC_ERR_PLIB, "Could not restore RowIJ of constraint matrix");
1746: PetscCall(MatAssemblyBegin(cMat, MAT_FINAL_ASSEMBLY));
1747: PetscCall(MatAssemblyEnd(cMat, MAT_FINAL_ASSEMBLY));
1748: PetscCall(PetscFree(vals));
1749: }
1751: /* clean up */
1752: PetscCall(ISRestoreIndices(anIS, &anchors));
1753: PetscCall(PetscFree2(perm, iperm));
1754: PetscCall(PetscFree(pointWork));
1755: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
1756: PetscFunctionReturn(PETSC_SUCCESS);
1757: }
1759: /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of
1760: * a non-conforming mesh. Local refinement comes later */
1761: PetscErrorCode DMPlexTreeRefineCell(DM dm, PetscInt cell, DM *ncdm)
1762: {
1763: DM K;
1764: PetscMPIInt rank;
1765: PetscInt dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd;
1766: PetscInt numNewCones, *newConeSizes, *newCones, *newOrientations;
1767: PetscInt *Kembedding;
1768: PetscInt *cellClosure = NULL, nc;
1769: PetscScalar *newVertexCoords;
1770: PetscInt numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset;
1771: PetscSection parentSection;
1773: PetscFunctionBegin;
1774: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1775: PetscCall(DMGetDimension(dm, &dim));
1776: PetscCall(DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm));
1777: PetscCall(DMSetDimension(*ncdm, dim));
1779: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1780: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &parentSection));
1781: PetscCall(DMPlexGetReferenceTree(dm, &K));
1782: PetscCall(DMGetCoordinatesLocalSetUp(dm));
1783: if (rank == 0) {
1784: /* compute the new charts */
1785: PetscCall(PetscMalloc5(dim + 1, &pNewCount, dim + 1, &pNewStart, dim + 1, &pNewEnd, dim + 1, &pOldStart, dim + 1, &pOldEnd));
1786: offset = 0;
1787: for (d = 0; d <= dim; d++) {
1788: PetscInt pOldCount, kStart, kEnd, k;
1790: pNewStart[d] = offset;
1791: PetscCall(DMPlexGetHeightStratum(dm, d, &pOldStart[d], &pOldEnd[d]));
1792: PetscCall(DMPlexGetHeightStratum(K, d, &kStart, &kEnd));
1793: pOldCount = pOldEnd[d] - pOldStart[d];
1794: /* adding the new points */
1795: pNewCount[d] = pOldCount + kEnd - kStart;
1796: if (!d) {
1797: /* removing the cell */
1798: pNewCount[d]--;
1799: }
1800: for (k = kStart; k < kEnd; k++) {
1801: PetscInt parent;
1802: PetscCall(DMPlexGetTreeParent(K, k, &parent, NULL));
1803: if (parent == k) {
1804: /* avoid double counting points that won't actually be new */
1805: pNewCount[d]--;
1806: }
1807: }
1808: pNewEnd[d] = pNewStart[d] + pNewCount[d];
1809: offset = pNewEnd[d];
1810: }
1811: PetscCheck(cell >= pOldStart[0] && cell < pOldEnd[0], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "%" PetscInt_FMT " not in cell range [%" PetscInt_FMT ", %" PetscInt_FMT ")", cell, pOldStart[0], pOldEnd[0]);
1812: /* get the current closure of the cell that we are removing */
1813: PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &nc, &cellClosure));
1815: PetscCall(PetscMalloc1(pNewEnd[dim], &newConeSizes));
1816: {
1817: DMPolytopeType pct, qct;
1818: PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j;
1820: PetscCall(DMPlexGetChart(K, &kStart, &kEnd));
1821: PetscCall(PetscMalloc4(kEnd - kStart, &Kembedding, kEnd - kStart, &perm, kEnd - kStart, &iperm, kEnd - kStart, &preOrient));
1823: for (k = kStart; k < kEnd; k++) {
1824: perm[k - kStart] = k;
1825: iperm[k - kStart] = k - kStart;
1826: preOrient[k - kStart] = 0;
1827: }
1829: PetscCall(DMPlexGetTransitiveClosure(K, 0, PETSC_TRUE, &closureSizeK, &closureK));
1830: for (j = 1; j < closureSizeK; j++) {
1831: PetscInt parentOrientA = closureK[2 * j + 1];
1832: PetscInt parentOrientB = cellClosure[2 * j + 1];
1833: PetscInt p, q;
1835: p = closureK[2 * j];
1836: q = cellClosure[2 * j];
1837: PetscCall(DMPlexGetCellType(K, p, &pct));
1838: PetscCall(DMPlexGetCellType(dm, q, &qct));
1839: for (d = 0; d <= dim; d++) {
1840: if (q >= pOldStart[d] && q < pOldEnd[d]) Kembedding[p] = (q - pOldStart[d]) + pNewStart[d];
1841: }
1842: parentOrientA = DMPolytopeConvertNewOrientation_Internal(pct, parentOrientA);
1843: parentOrientB = DMPolytopeConvertNewOrientation_Internal(qct, parentOrientB);
1844: if (parentOrientA != parentOrientB) {
1845: PetscInt numChildren, i;
1846: const PetscInt *children;
1848: PetscCall(DMPlexGetTreeChildren(K, p, &numChildren, &children));
1849: for (i = 0; i < numChildren; i++) {
1850: PetscInt kPerm, oPerm;
1852: k = children[i];
1853: PetscCall(DMPlexReferenceTreeGetChildSymmetry(K, p, parentOrientA, 0, k, parentOrientB, &oPerm, &kPerm));
1854: /* perm = what refTree position I'm in */
1855: perm[kPerm - kStart] = k;
1856: /* iperm = who is at this position */
1857: iperm[k - kStart] = kPerm - kStart;
1858: preOrient[kPerm - kStart] = oPerm;
1859: }
1860: }
1861: }
1862: PetscCall(DMPlexRestoreTransitiveClosure(K, 0, PETSC_TRUE, &closureSizeK, &closureK));
1863: }
1864: PetscCall(PetscSectionSetChart(parentSection, 0, pNewEnd[dim]));
1865: offset = 0;
1866: numNewCones = 0;
1867: for (d = 0; d <= dim; d++) {
1868: PetscInt kStart, kEnd, k;
1869: PetscInt p;
1870: PetscInt size;
1872: for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1873: /* skip cell 0 */
1874: if (p == cell) continue;
1875: /* old cones to new cones */
1876: PetscCall(DMPlexGetConeSize(dm, p, &size));
1877: newConeSizes[offset++] = size;
1878: numNewCones += size;
1879: }
1881: PetscCall(DMPlexGetHeightStratum(K, d, &kStart, &kEnd));
1882: for (k = kStart; k < kEnd; k++) {
1883: PetscInt kParent;
1885: PetscCall(DMPlexGetTreeParent(K, k, &kParent, NULL));
1886: if (kParent != k) {
1887: Kembedding[k] = offset;
1888: PetscCall(DMPlexGetConeSize(K, k, &size));
1889: newConeSizes[offset++] = size;
1890: numNewCones += size;
1891: if (kParent != 0) PetscCall(PetscSectionSetDof(parentSection, Kembedding[k], 1));
1892: }
1893: }
1894: }
1896: PetscCall(PetscSectionSetUp(parentSection));
1897: PetscCall(PetscSectionGetStorageSize(parentSection, &numPointsWithParents));
1898: PetscCall(PetscMalloc2(numNewCones, &newCones, numNewCones, &newOrientations));
1899: PetscCall(PetscMalloc2(numPointsWithParents, &parents, numPointsWithParents, &childIDs));
1901: /* fill new cones */
1902: offset = 0;
1903: for (d = 0; d <= dim; d++) {
1904: PetscInt kStart, kEnd, k, l;
1905: PetscInt p;
1906: PetscInt size;
1907: const PetscInt *cone, *orientation;
1909: for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1910: /* skip cell 0 */
1911: if (p == cell) continue;
1912: /* old cones to new cones */
1913: PetscCall(DMPlexGetConeSize(dm, p, &size));
1914: PetscCall(DMPlexGetCone(dm, p, &cone));
1915: PetscCall(DMPlexGetConeOrientation(dm, p, &orientation));
1916: for (l = 0; l < size; l++) {
1917: newCones[offset] = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1];
1918: newOrientations[offset++] = orientation[l];
1919: }
1920: }
1922: PetscCall(DMPlexGetHeightStratum(K, d, &kStart, &kEnd));
1923: for (k = kStart; k < kEnd; k++) {
1924: PetscInt kPerm = perm[k], kParent;
1925: PetscInt preO = preOrient[k];
1927: PetscCall(DMPlexGetTreeParent(K, k, &kParent, NULL));
1928: if (kParent != k) {
1929: /* embed new cones */
1930: PetscCall(DMPlexGetConeSize(K, k, &size));
1931: PetscCall(DMPlexGetCone(K, kPerm, &cone));
1932: PetscCall(DMPlexGetConeOrientation(K, kPerm, &orientation));
1933: for (l = 0; l < size; l++) {
1934: PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size - (preO + 1) - l) % size);
1935: PetscInt newO, lSize, oTrue;
1936: DMPolytopeType ct = DM_NUM_POLYTOPES;
1938: q = iperm[cone[m]];
1939: newCones[offset] = Kembedding[q];
1940: PetscCall(DMPlexGetConeSize(K, q, &lSize));
1941: if (lSize == 2) ct = DM_POLYTOPE_SEGMENT;
1942: else if (lSize == 4) ct = DM_POLYTOPE_QUADRILATERAL;
1943: oTrue = DMPolytopeConvertNewOrientation_Internal(ct, orientation[m]);
1944: oTrue = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2);
1945: newO = DihedralCompose(lSize, oTrue, preOrient[q]);
1946: newOrientations[offset++] = DMPolytopeConvertOldOrientation_Internal(ct, newO);
1947: }
1948: if (kParent != 0) {
1949: PetscInt newPoint = Kembedding[kParent];
1950: PetscCall(PetscSectionGetOffset(parentSection, Kembedding[k], &pOffset));
1951: parents[pOffset] = newPoint;
1952: childIDs[pOffset] = k;
1953: }
1954: }
1955: }
1956: }
1958: PetscCall(PetscMalloc1(dim * (pNewEnd[dim] - pNewStart[dim]), &newVertexCoords));
1960: /* fill coordinates */
1961: offset = 0;
1962: {
1963: PetscInt kStart, kEnd, l;
1964: PetscSection vSection;
1965: PetscInt v;
1966: Vec coords;
1967: PetscScalar *coordvals;
1968: PetscInt dof, off;
1969: PetscReal v0[3], J[9], detJ;
1971: if (PetscDefined(USE_DEBUG)) {
1972: PetscInt k;
1973: PetscCall(DMPlexGetHeightStratum(K, 0, &kStart, &kEnd));
1974: for (k = kStart; k < kEnd; k++) {
1975: PetscCall(DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ));
1976: PetscCheck(detJ > 0., PETSC_COMM_SELF, PETSC_ERR_PLIB, "reference tree cell %" PetscInt_FMT " has bad determinant", k);
1977: }
1978: }
1979: PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ));
1980: PetscCall(DMGetCoordinateSection(dm, &vSection));
1981: PetscCall(DMGetCoordinatesLocal(dm, &coords));
1982: PetscCall(VecGetArray(coords, &coordvals));
1983: for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) {
1984: PetscCall(PetscSectionGetDof(vSection, v, &dof));
1985: PetscCall(PetscSectionGetOffset(vSection, v, &off));
1986: for (l = 0; l < dof; l++) newVertexCoords[offset++] = coordvals[off + l];
1987: }
1988: PetscCall(VecRestoreArray(coords, &coordvals));
1990: PetscCall(DMGetCoordinateSection(K, &vSection));
1991: PetscCall(DMGetCoordinatesLocal(K, &coords));
1992: PetscCall(VecGetArray(coords, &coordvals));
1993: PetscCall(DMPlexGetDepthStratum(K, 0, &kStart, &kEnd));
1994: for (v = kStart; v < kEnd; v++) {
1995: PetscReal coord[3], newCoord[3];
1996: PetscInt vPerm = perm[v];
1997: PetscInt kParent;
1998: const PetscReal xi0[3] = {-1., -1., -1.};
2000: PetscCall(DMPlexGetTreeParent(K, v, &kParent, NULL));
2001: if (kParent != v) {
2002: /* this is a new vertex */
2003: PetscCall(PetscSectionGetOffset(vSection, vPerm, &off));
2004: for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off + l]);
2005: CoordinatesRefToReal(dim, dim, xi0, v0, J, coord, newCoord);
2006: for (l = 0; l < dim; ++l) newVertexCoords[offset + l] = newCoord[l];
2007: offset += dim;
2008: }
2009: }
2010: PetscCall(VecRestoreArray(coords, &coordvals));
2011: }
2013: /* need to reverse the order of pNewCount: vertices first, cells last */
2014: for (d = 0; d < (dim + 1) / 2; d++) {
2015: PetscInt tmp;
2017: tmp = pNewCount[d];
2018: pNewCount[d] = pNewCount[dim - d];
2019: pNewCount[dim - d] = tmp;
2020: }
2022: PetscCall(DMPlexCreateFromDAG(*ncdm, dim, pNewCount, newConeSizes, newCones, newOrientations, newVertexCoords));
2023: PetscCall(DMPlexSetReferenceTree(*ncdm, K));
2024: PetscCall(DMPlexSetTree(*ncdm, parentSection, parents, childIDs));
2026: /* clean up */
2027: PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &nc, &cellClosure));
2028: PetscCall(PetscFree5(pNewCount, pNewStart, pNewEnd, pOldStart, pOldEnd));
2029: PetscCall(PetscFree(newConeSizes));
2030: PetscCall(PetscFree2(newCones, newOrientations));
2031: PetscCall(PetscFree(newVertexCoords));
2032: PetscCall(PetscFree2(parents, childIDs));
2033: PetscCall(PetscFree4(Kembedding, perm, iperm, preOrient));
2034: } else {
2035: PetscInt p, counts[4];
2036: PetscInt *coneSizes, *cones, *orientations;
2037: Vec coordVec;
2038: PetscScalar *coords;
2040: for (d = 0; d <= dim; d++) {
2041: PetscInt dStart, dEnd;
2043: PetscCall(DMPlexGetDepthStratum(dm, d, &dStart, &dEnd));
2044: counts[d] = dEnd - dStart;
2045: }
2046: PetscCall(PetscMalloc1(pEnd - pStart, &coneSizes));
2047: for (p = pStart; p < pEnd; p++) PetscCall(DMPlexGetConeSize(dm, p, &coneSizes[p - pStart]));
2048: PetscCall(DMPlexGetCones(dm, &cones));
2049: PetscCall(DMPlexGetConeOrientations(dm, &orientations));
2050: PetscCall(DMGetCoordinatesLocal(dm, &coordVec));
2051: PetscCall(VecGetArray(coordVec, &coords));
2053: PetscCall(PetscSectionSetChart(parentSection, pStart, pEnd));
2054: PetscCall(PetscSectionSetUp(parentSection));
2055: PetscCall(DMPlexCreateFromDAG(*ncdm, dim, counts, coneSizes, cones, orientations, NULL));
2056: PetscCall(DMPlexSetReferenceTree(*ncdm, K));
2057: PetscCall(DMPlexSetTree(*ncdm, parentSection, NULL, NULL));
2058: PetscCall(VecRestoreArray(coordVec, &coords));
2059: }
2060: PetscCall(PetscSectionDestroy(&parentSection));
2062: PetscFunctionReturn(PETSC_SUCCESS);
2063: }
2065: PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
2066: {
2067: PetscSF coarseToFineEmbedded;
2068: PetscSection globalCoarse, globalFine;
2069: PetscSection localCoarse, localFine;
2070: PetscSection aSec, cSec;
2071: PetscSection rootIndicesSec, rootMatricesSec;
2072: PetscSection leafIndicesSec, leafMatricesSec;
2073: PetscInt *rootIndices, *leafIndices;
2074: PetscScalar *rootMatrices, *leafMatrices;
2075: IS aIS;
2076: const PetscInt *anchors;
2077: Mat cMat;
2078: PetscInt numFields, maxFields;
2079: PetscInt pStartC, pEndC, pStartF, pEndF, p;
2080: PetscInt aStart, aEnd, cStart, cEnd;
2081: PetscInt *maxChildIds;
2082: PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
2083: const PetscInt ***perms;
2084: const PetscScalar ***flips;
2086: PetscFunctionBegin;
2087: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
2088: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
2089: PetscCall(DMGetGlobalSection(fine, &globalFine));
2090: { /* winnow fine points that don't have global dofs out of the sf */
2091: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, nleaves, l;
2092: const PetscInt *leaves;
2094: PetscCall(PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL));
2095: for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
2096: p = leaves ? leaves[l] : l;
2097: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
2098: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
2099: if ((dof - cdof) > 0) numPointsWithDofs++;
2100: }
2101: PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
2102: for (l = 0, offset = 0; l < nleaves; l++) {
2103: p = leaves ? leaves[l] : l;
2104: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
2105: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
2106: if ((dof - cdof) > 0) pointsWithDofs[offset++] = l;
2107: }
2108: PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded));
2109: PetscCall(PetscFree(pointsWithDofs));
2110: }
2111: /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
2112: PetscCall(PetscMalloc1(pEndC - pStartC, &maxChildIds));
2113: for (p = pStartC; p < pEndC; p++) maxChildIds[p - pStartC] = -2;
2114: PetscCall(PetscSFReduceBegin(coarseToFineEmbedded, MPIU_INT, childIds, maxChildIds, MPI_MAX));
2115: PetscCall(PetscSFReduceEnd(coarseToFineEmbedded, MPIU_INT, childIds, maxChildIds, MPI_MAX));
2117: PetscCall(DMGetLocalSection(coarse, &localCoarse));
2118: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
2120: PetscCall(DMPlexGetAnchors(coarse, &aSec, &aIS));
2121: PetscCall(ISGetIndices(aIS, &anchors));
2122: PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd));
2124: PetscCall(DMGetDefaultConstraints(coarse, &cSec, &cMat, NULL));
2125: PetscCall(PetscSectionGetChart(cSec, &cStart, &cEnd));
2127: /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
2128: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootIndicesSec));
2129: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootMatricesSec));
2130: PetscCall(PetscSectionSetChart(rootIndicesSec, pStartC, pEndC));
2131: PetscCall(PetscSectionSetChart(rootMatricesSec, pStartC, pEndC));
2132: PetscCall(PetscSectionGetNumFields(localCoarse, &numFields));
2133: maxFields = PetscMax(1, numFields);
2134: PetscCall(PetscMalloc7(maxFields + 1, &offsets, maxFields + 1, &offsetsCopy, maxFields + 1, &newOffsets, maxFields + 1, &newOffsetsCopy, maxFields + 1, &rowOffsets, maxFields + 1, &numD, maxFields + 1, &numO));
2135: PetscCall(PetscMalloc2(maxFields + 1, (PetscInt ****)&perms, maxFields + 1, (PetscScalar ****)&flips));
2136: PetscCall(PetscMemzero((void *)perms, (maxFields + 1) * sizeof(const PetscInt **)));
2137: PetscCall(PetscMemzero((void *)flips, (maxFields + 1) * sizeof(const PetscScalar **)));
2139: for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
2140: PetscInt dof, matSize = 0;
2141: PetscInt aDof = 0;
2142: PetscInt cDof = 0;
2143: PetscInt maxChildId = maxChildIds[p - pStartC];
2144: PetscInt numRowIndices = 0;
2145: PetscInt numColIndices = 0;
2146: PetscInt f;
2148: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
2149: if (dof < 0) dof = -(dof + 1);
2150: if (p >= aStart && p < aEnd) PetscCall(PetscSectionGetDof(aSec, p, &aDof));
2151: if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(cSec, p, &cDof));
2152: for (f = 0; f <= numFields; f++) offsets[f] = 0;
2153: for (f = 0; f <= numFields; f++) newOffsets[f] = 0;
2154: if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
2155: PetscInt *closure = NULL, closureSize, cl;
2157: PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2158: for (cl = 0; cl < closureSize; cl++) { /* get the closure */
2159: PetscInt c = closure[2 * cl], clDof;
2161: PetscCall(PetscSectionGetDof(localCoarse, c, &clDof));
2162: numRowIndices += clDof;
2163: for (f = 0; f < numFields; f++) {
2164: PetscCall(PetscSectionGetFieldDof(localCoarse, c, f, &clDof));
2165: offsets[f + 1] += clDof;
2166: }
2167: }
2168: for (f = 0; f < numFields; f++) {
2169: offsets[f + 1] += offsets[f];
2170: newOffsets[f + 1] = offsets[f + 1];
2171: }
2172: /* get the number of indices needed and their field offsets */
2173: PetscCall(DMPlexAnchorsModifyMat(coarse, localCoarse, closureSize, numRowIndices, closure, NULL, NULL, NULL, &numColIndices, NULL, NULL, newOffsets, PETSC_FALSE));
2174: PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2175: if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */
2176: numColIndices = numRowIndices;
2177: matSize = 0;
2178: } else if (numFields) { /* we send one submat for each field: sum their sizes */
2179: matSize = 0;
2180: for (f = 0; f < numFields; f++) {
2181: PetscInt numRow, numCol;
2183: numRow = offsets[f + 1] - offsets[f];
2184: numCol = newOffsets[f + 1] - newOffsets[f];
2185: matSize += numRow * numCol;
2186: }
2187: } else {
2188: matSize = numRowIndices * numColIndices;
2189: }
2190: } else if (maxChildId == -1) {
2191: if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */
2192: PetscInt aOff, a;
2194: PetscCall(PetscSectionGetOffset(aSec, p, &aOff));
2195: for (f = 0; f < numFields; f++) {
2196: PetscInt fDof;
2198: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
2199: offsets[f + 1] = fDof;
2200: }
2201: for (a = 0; a < aDof; a++) {
2202: PetscInt anchor = anchors[a + aOff], aLocalDof;
2204: PetscCall(PetscSectionGetDof(localCoarse, anchor, &aLocalDof));
2205: numColIndices += aLocalDof;
2206: for (f = 0; f < numFields; f++) {
2207: PetscInt fDof;
2209: PetscCall(PetscSectionGetFieldDof(localCoarse, anchor, f, &fDof));
2210: newOffsets[f + 1] += fDof;
2211: }
2212: }
2213: if (numFields) {
2214: matSize = 0;
2215: for (f = 0; f < numFields; f++) matSize += offsets[f + 1] * newOffsets[f + 1];
2216: } else {
2217: matSize = numColIndices * dof;
2218: }
2219: } else { /* no children, and no constraints on dofs: just get the global indices */
2220: numColIndices = dof;
2221: matSize = 0;
2222: }
2223: }
2224: /* we will pack the column indices with the field offsets */
2225: PetscCall(PetscSectionSetDof(rootIndicesSec, p, numColIndices ? numColIndices + 2 * numFields : 0));
2226: PetscCall(PetscSectionSetDof(rootMatricesSec, p, matSize));
2227: }
2228: PetscCall(PetscSectionSetUp(rootIndicesSec));
2229: PetscCall(PetscSectionSetUp(rootMatricesSec));
2230: {
2231: PetscInt numRootIndices, numRootMatrices;
2233: PetscCall(PetscSectionGetStorageSize(rootIndicesSec, &numRootIndices));
2234: PetscCall(PetscSectionGetStorageSize(rootMatricesSec, &numRootMatrices));
2235: PetscCall(PetscMalloc2(numRootIndices, &rootIndices, numRootMatrices, &rootMatrices));
2236: for (p = pStartC; p < pEndC; p++) {
2237: PetscInt numRowIndices, numColIndices, matSize, dof;
2238: PetscInt pIndOff, pMatOff, f;
2239: PetscInt *pInd;
2240: PetscInt maxChildId = maxChildIds[p - pStartC];
2241: PetscScalar *pMat = NULL;
2243: PetscCall(PetscSectionGetDof(rootIndicesSec, p, &numColIndices));
2244: if (!numColIndices) continue;
2245: for (f = 0; f <= numFields; f++) {
2246: offsets[f] = 0;
2247: newOffsets[f] = 0;
2248: offsetsCopy[f] = 0;
2249: newOffsetsCopy[f] = 0;
2250: }
2251: numColIndices -= 2 * numFields;
2252: PetscCall(PetscSectionGetOffset(rootIndicesSec, p, &pIndOff));
2253: pInd = &(rootIndices[pIndOff]);
2254: PetscCall(PetscSectionGetDof(rootMatricesSec, p, &matSize));
2255: if (matSize) {
2256: PetscCall(PetscSectionGetOffset(rootMatricesSec, p, &pMatOff));
2257: pMat = &rootMatrices[pMatOff];
2258: }
2259: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
2260: if (dof < 0) dof = -(dof + 1);
2261: if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
2262: PetscInt i, j;
2263: PetscInt numRowIndices = matSize / numColIndices;
2265: if (!numRowIndices) { /* don't need to calculate the mat, just the indices */
2266: PetscInt numIndices, *indices;
2267: PetscCall(DMPlexGetClosureIndices(coarse, localCoarse, globalCoarse, p, PETSC_TRUE, &numIndices, &indices, offsets, NULL));
2268: PetscCheck(numIndices == numColIndices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "mismatching constraint indices calculations");
2269: for (i = 0; i < numColIndices; i++) pInd[i] = indices[i];
2270: for (i = 0; i < numFields; i++) {
2271: pInd[numColIndices + i] = offsets[i + 1];
2272: pInd[numColIndices + numFields + i] = offsets[i + 1];
2273: }
2274: PetscCall(DMPlexRestoreClosureIndices(coarse, localCoarse, globalCoarse, p, PETSC_TRUE, &numIndices, &indices, offsets, NULL));
2275: } else {
2276: PetscInt closureSize, *closure = NULL, cl;
2277: PetscScalar *pMatIn, *pMatModified;
2278: PetscInt numPoints, *points;
2280: PetscCall(DMGetWorkArray(coarse, numRowIndices * numRowIndices, MPIU_SCALAR, &pMatIn));
2281: for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */
2282: for (j = 0; j < numRowIndices; j++) pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.;
2283: }
2284: PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2285: for (f = 0; f < maxFields; f++) {
2286: if (numFields) PetscCall(PetscSectionGetFieldPointSyms(localCoarse, f, closureSize, closure, &perms[f], &flips[f]));
2287: else PetscCall(PetscSectionGetPointSyms(localCoarse, closureSize, closure, &perms[f], &flips[f]));
2288: }
2289: if (numFields) {
2290: for (cl = 0; cl < closureSize; cl++) {
2291: PetscInt c = closure[2 * cl];
2293: for (f = 0; f < numFields; f++) {
2294: PetscInt fDof;
2296: PetscCall(PetscSectionGetFieldDof(localCoarse, c, f, &fDof));
2297: offsets[f + 1] += fDof;
2298: }
2299: }
2300: for (f = 0; f < numFields; f++) {
2301: offsets[f + 1] += offsets[f];
2302: newOffsets[f + 1] = offsets[f + 1];
2303: }
2304: }
2305: /* TODO : flips here ? */
2306: /* apply hanging node constraints on the right, get the new points and the new offsets */
2307: PetscCall(DMPlexAnchorsModifyMat(coarse, localCoarse, closureSize, numRowIndices, closure, perms, pMatIn, &numPoints, NULL, &points, &pMatModified, newOffsets, PETSC_FALSE));
2308: for (f = 0; f < maxFields; f++) {
2309: if (numFields) PetscCall(PetscSectionRestoreFieldPointSyms(localCoarse, f, closureSize, closure, &perms[f], &flips[f]));
2310: else PetscCall(PetscSectionRestorePointSyms(localCoarse, closureSize, closure, &perms[f], &flips[f]));
2311: }
2312: for (f = 0; f < maxFields; f++) {
2313: if (numFields) PetscCall(PetscSectionGetFieldPointSyms(localCoarse, f, numPoints, points, &perms[f], &flips[f]));
2314: else PetscCall(PetscSectionGetPointSyms(localCoarse, numPoints, points, &perms[f], &flips[f]));
2315: }
2316: if (!numFields) {
2317: for (i = 0; i < numRowIndices * numColIndices; i++) pMat[i] = pMatModified[i];
2318: } else {
2319: PetscInt i, j, count;
2320: for (f = 0, count = 0; f < numFields; f++) {
2321: for (i = offsets[f]; i < offsets[f + 1]; i++) {
2322: for (j = newOffsets[f]; j < newOffsets[f + 1]; j++, count++) pMat[count] = pMatModified[i * numColIndices + j];
2323: }
2324: }
2325: }
2326: PetscCall(DMRestoreWorkArray(coarse, numRowIndices * numColIndices, MPIU_SCALAR, &pMatModified));
2327: PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2328: PetscCall(DMRestoreWorkArray(coarse, numRowIndices * numColIndices, MPIU_SCALAR, &pMatIn));
2329: if (numFields) {
2330: for (f = 0; f < numFields; f++) {
2331: pInd[numColIndices + f] = offsets[f + 1];
2332: pInd[numColIndices + numFields + f] = newOffsets[f + 1];
2333: }
2334: for (cl = 0; cl < numPoints; cl++) {
2335: PetscInt globalOff, c = points[2 * cl];
2336: PetscCall(PetscSectionGetOffset(globalCoarse, c, &globalOff));
2337: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff + 1) : globalOff, newOffsets, PETSC_FALSE, perms, cl, NULL, pInd));
2338: }
2339: } else {
2340: for (cl = 0; cl < numPoints; cl++) {
2341: PetscInt c = points[2 * cl], globalOff;
2342: const PetscInt *perm = perms[0] ? perms[0][cl] : NULL;
2344: PetscCall(PetscSectionGetOffset(globalCoarse, c, &globalOff));
2345: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff + 1) : globalOff, newOffsets, PETSC_FALSE, perm, NULL, pInd));
2346: }
2347: }
2348: for (f = 0; f < maxFields; f++) {
2349: if (numFields) PetscCall(PetscSectionRestoreFieldPointSyms(localCoarse, f, numPoints, points, &perms[f], &flips[f]));
2350: else PetscCall(PetscSectionRestorePointSyms(localCoarse, numPoints, points, &perms[f], &flips[f]));
2351: }
2352: PetscCall(DMRestoreWorkArray(coarse, numPoints, MPIU_SCALAR, &points));
2353: }
2354: } else if (matSize) {
2355: PetscInt cOff;
2356: PetscInt *rowIndices, *colIndices, a, aDof, aOff;
2358: numRowIndices = matSize / numColIndices;
2359: PetscCheck(numRowIndices == dof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Miscounted dofs");
2360: PetscCall(DMGetWorkArray(coarse, numRowIndices, MPIU_INT, &rowIndices));
2361: PetscCall(DMGetWorkArray(coarse, numColIndices, MPIU_INT, &colIndices));
2362: PetscCall(PetscSectionGetOffset(cSec, p, &cOff));
2363: PetscCall(PetscSectionGetDof(aSec, p, &aDof));
2364: PetscCall(PetscSectionGetOffset(aSec, p, &aOff));
2365: if (numFields) {
2366: for (f = 0; f < numFields; f++) {
2367: PetscInt fDof;
2369: PetscCall(PetscSectionGetFieldDof(cSec, p, f, &fDof));
2370: offsets[f + 1] = fDof;
2371: for (a = 0; a < aDof; a++) {
2372: PetscInt anchor = anchors[a + aOff];
2373: PetscCall(PetscSectionGetFieldDof(localCoarse, anchor, f, &fDof));
2374: newOffsets[f + 1] += fDof;
2375: }
2376: }
2377: for (f = 0; f < numFields; f++) {
2378: offsets[f + 1] += offsets[f];
2379: offsetsCopy[f + 1] = offsets[f + 1];
2380: newOffsets[f + 1] += newOffsets[f];
2381: newOffsetsCopy[f + 1] = newOffsets[f + 1];
2382: }
2383: PetscCall(DMPlexGetIndicesPointFields_Internal(cSec, PETSC_TRUE, p, cOff, offsetsCopy, PETSC_TRUE, NULL, -1, NULL, rowIndices));
2384: for (a = 0; a < aDof; a++) {
2385: PetscInt anchor = anchors[a + aOff], lOff;
2386: PetscCall(PetscSectionGetOffset(localCoarse, anchor, &lOff));
2387: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_TRUE, anchor, lOff, newOffsetsCopy, PETSC_TRUE, NULL, -1, NULL, colIndices));
2388: }
2389: } else {
2390: PetscCall(DMPlexGetIndicesPoint_Internal(cSec, PETSC_TRUE, p, cOff, offsetsCopy, PETSC_TRUE, NULL, NULL, rowIndices));
2391: for (a = 0; a < aDof; a++) {
2392: PetscInt anchor = anchors[a + aOff], lOff;
2393: PetscCall(PetscSectionGetOffset(localCoarse, anchor, &lOff));
2394: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_TRUE, anchor, lOff, newOffsetsCopy, PETSC_TRUE, NULL, NULL, colIndices));
2395: }
2396: }
2397: if (numFields) {
2398: PetscInt count, a;
2400: for (f = 0, count = 0; f < numFields; f++) {
2401: PetscInt iSize = offsets[f + 1] - offsets[f];
2402: PetscInt jSize = newOffsets[f + 1] - newOffsets[f];
2403: PetscCall(MatGetValues(cMat, iSize, &rowIndices[offsets[f]], jSize, &colIndices[newOffsets[f]], &pMat[count]));
2404: count += iSize * jSize;
2405: pInd[numColIndices + f] = offsets[f + 1];
2406: pInd[numColIndices + numFields + f] = newOffsets[f + 1];
2407: }
2408: for (a = 0; a < aDof; a++) {
2409: PetscInt anchor = anchors[a + aOff];
2410: PetscInt gOff;
2411: PetscCall(PetscSectionGetOffset(globalCoarse, anchor, &gOff));
2412: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, anchor, gOff < 0 ? -(gOff + 1) : gOff, newOffsets, PETSC_FALSE, NULL, -1, NULL, pInd));
2413: }
2414: } else {
2415: PetscInt a;
2416: PetscCall(MatGetValues(cMat, numRowIndices, rowIndices, numColIndices, colIndices, pMat));
2417: for (a = 0; a < aDof; a++) {
2418: PetscInt anchor = anchors[a + aOff];
2419: PetscInt gOff;
2420: PetscCall(PetscSectionGetOffset(globalCoarse, anchor, &gOff));
2421: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, anchor, gOff < 0 ? -(gOff + 1) : gOff, newOffsets, PETSC_FALSE, NULL, NULL, pInd));
2422: }
2423: }
2424: PetscCall(DMRestoreWorkArray(coarse, numColIndices, MPIU_INT, &colIndices));
2425: PetscCall(DMRestoreWorkArray(coarse, numRowIndices, MPIU_INT, &rowIndices));
2426: } else {
2427: PetscInt gOff;
2429: PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));
2430: if (numFields) {
2431: for (f = 0; f < numFields; f++) {
2432: PetscInt fDof;
2433: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
2434: offsets[f + 1] = fDof + offsets[f];
2435: }
2436: for (f = 0; f < numFields; f++) {
2437: pInd[numColIndices + f] = offsets[f + 1];
2438: pInd[numColIndices + numFields + f] = offsets[f + 1];
2439: }
2440: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, -1, NULL, pInd));
2441: } else {
2442: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, NULL, pInd));
2443: }
2444: }
2445: }
2446: PetscCall(PetscFree(maxChildIds));
2447: }
2448: {
2449: PetscSF indicesSF, matricesSF;
2450: PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices;
2452: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafIndicesSec));
2453: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafMatricesSec));
2454: PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootIndicesSec, &remoteOffsetsIndices, leafIndicesSec));
2455: PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootMatricesSec, &remoteOffsetsMatrices, leafMatricesSec));
2456: PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootIndicesSec, remoteOffsetsIndices, leafIndicesSec, &indicesSF));
2457: PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootMatricesSec, remoteOffsetsMatrices, leafMatricesSec, &matricesSF));
2458: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
2459: PetscCall(PetscFree(remoteOffsetsIndices));
2460: PetscCall(PetscFree(remoteOffsetsMatrices));
2461: PetscCall(PetscSectionGetStorageSize(leafIndicesSec, &numLeafIndices));
2462: PetscCall(PetscSectionGetStorageSize(leafMatricesSec, &numLeafMatrices));
2463: PetscCall(PetscMalloc2(numLeafIndices, &leafIndices, numLeafMatrices, &leafMatrices));
2464: PetscCall(PetscSFBcastBegin(indicesSF, MPIU_INT, rootIndices, leafIndices, MPI_REPLACE));
2465: PetscCall(PetscSFBcastBegin(matricesSF, MPIU_SCALAR, rootMatrices, leafMatrices, MPI_REPLACE));
2466: PetscCall(PetscSFBcastEnd(indicesSF, MPIU_INT, rootIndices, leafIndices, MPI_REPLACE));
2467: PetscCall(PetscSFBcastEnd(matricesSF, MPIU_SCALAR, rootMatrices, leafMatrices, MPI_REPLACE));
2468: PetscCall(PetscSFDestroy(&matricesSF));
2469: PetscCall(PetscSFDestroy(&indicesSF));
2470: PetscCall(PetscFree2(rootIndices, rootMatrices));
2471: PetscCall(PetscSectionDestroy(&rootIndicesSec));
2472: PetscCall(PetscSectionDestroy(&rootMatricesSec));
2473: }
2474: /* count to preallocate */
2475: PetscCall(DMGetLocalSection(fine, &localFine));
2476: {
2477: PetscInt nGlobal;
2478: PetscInt *dnnz, *onnz;
2479: PetscLayout rowMap, colMap;
2480: PetscInt rowStart, rowEnd, colStart, colEnd;
2481: PetscInt maxDof;
2482: PetscInt *rowIndices;
2483: DM refTree;
2484: PetscInt **refPointFieldN;
2485: PetscScalar ***refPointFieldMats;
2486: PetscSection refConSec, refAnSec;
2487: PetscInt pRefStart, pRefEnd, maxConDof, maxColumns, leafStart, leafEnd;
2488: PetscScalar *pointWork;
2490: PetscCall(PetscSectionGetConstrainedStorageSize(globalFine, &nGlobal));
2491: PetscCall(PetscCalloc2(nGlobal, &dnnz, nGlobal, &onnz));
2492: PetscCall(MatGetLayouts(mat, &rowMap, &colMap));
2493: PetscCall(PetscLayoutSetUp(rowMap));
2494: PetscCall(PetscLayoutSetUp(colMap));
2495: PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd));
2496: PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd));
2497: PetscCall(PetscSectionGetMaxDof(localFine, &maxDof));
2498: PetscCall(PetscSectionGetChart(leafIndicesSec, &leafStart, &leafEnd));
2499: PetscCall(DMGetWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
2500: for (p = leafStart; p < leafEnd; p++) {
2501: PetscInt gDof, gcDof, gOff;
2502: PetscInt numColIndices, pIndOff, *pInd;
2503: PetscInt matSize;
2504: PetscInt i;
2506: PetscCall(PetscSectionGetDof(globalFine, p, &gDof));
2507: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &gcDof));
2508: if ((gDof - gcDof) <= 0) continue;
2509: PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
2510: PetscCheck(gOff >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "I though having global dofs meant a non-negative offset");
2511: PetscCheck(gOff >= rowStart && (gOff + gDof - gcDof) <= rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "I thought the row map would constrain the global dofs");
2512: PetscCall(PetscSectionGetDof(leafIndicesSec, p, &numColIndices));
2513: PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &pIndOff));
2514: numColIndices -= 2 * numFields;
2515: PetscCheck(numColIndices > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "global fine dof with no dofs to interpolate from");
2516: pInd = &leafIndices[pIndOff];
2517: offsets[0] = 0;
2518: offsetsCopy[0] = 0;
2519: newOffsets[0] = 0;
2520: newOffsetsCopy[0] = 0;
2521: if (numFields) {
2522: PetscInt f;
2523: for (f = 0; f < numFields; f++) {
2524: PetscInt rowDof;
2526: PetscCall(PetscSectionGetFieldDof(localFine, p, f, &rowDof));
2527: offsets[f + 1] = offsets[f] + rowDof;
2528: offsetsCopy[f + 1] = offsets[f + 1];
2529: newOffsets[f + 1] = pInd[numColIndices + numFields + f];
2530: numD[f] = 0;
2531: numO[f] = 0;
2532: }
2533: PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices));
2534: for (f = 0; f < numFields; f++) {
2535: PetscInt colOffset = newOffsets[f];
2536: PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f];
2538: for (i = 0; i < numFieldCols; i++) {
2539: PetscInt gInd = pInd[i + colOffset];
2541: if (gInd >= colStart && gInd < colEnd) {
2542: numD[f]++;
2543: } else if (gInd >= 0) { /* negative means non-entry */
2544: numO[f]++;
2545: }
2546: }
2547: }
2548: } else {
2549: PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices));
2550: numD[0] = 0;
2551: numO[0] = 0;
2552: for (i = 0; i < numColIndices; i++) {
2553: PetscInt gInd = pInd[i];
2555: if (gInd >= colStart && gInd < colEnd) {
2556: numD[0]++;
2557: } else if (gInd >= 0) { /* negative means non-entry */
2558: numO[0]++;
2559: }
2560: }
2561: }
2562: PetscCall(PetscSectionGetDof(leafMatricesSec, p, &matSize));
2563: if (!matSize) { /* incoming matrix is identity */
2564: PetscInt childId;
2566: childId = childIds[p - pStartF];
2567: if (childId < 0) { /* no child interpolation: one nnz per */
2568: if (numFields) {
2569: PetscInt f;
2570: for (f = 0; f < numFields; f++) {
2571: PetscInt numRows = offsets[f + 1] - offsets[f], row;
2572: for (row = 0; row < numRows; row++) {
2573: PetscInt gIndCoarse = pInd[newOffsets[f] + row];
2574: PetscInt gIndFine = rowIndices[offsets[f] + row];
2575: if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2576: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2577: dnnz[gIndFine - rowStart] = 1;
2578: } else if (gIndCoarse >= 0) { /* remote */
2579: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2580: onnz[gIndFine - rowStart] = 1;
2581: } else { /* constrained */
2582: PetscCheck(gIndFine < 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2583: }
2584: }
2585: }
2586: } else {
2587: PetscInt i;
2588: for (i = 0; i < gDof; i++) {
2589: PetscInt gIndCoarse = pInd[i];
2590: PetscInt gIndFine = rowIndices[i];
2591: if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2592: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2593: dnnz[gIndFine - rowStart] = 1;
2594: } else if (gIndCoarse >= 0) { /* remote */
2595: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2596: onnz[gIndFine - rowStart] = 1;
2597: } else { /* constrained */
2598: PetscCheck(gIndFine < 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2599: }
2600: }
2601: }
2602: } else { /* interpolate from all */
2603: if (numFields) {
2604: PetscInt f;
2605: for (f = 0; f < numFields; f++) {
2606: PetscInt numRows = offsets[f + 1] - offsets[f], row;
2607: for (row = 0; row < numRows; row++) {
2608: PetscInt gIndFine = rowIndices[offsets[f] + row];
2609: if (gIndFine >= 0) {
2610: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2611: dnnz[gIndFine - rowStart] = numD[f];
2612: onnz[gIndFine - rowStart] = numO[f];
2613: }
2614: }
2615: }
2616: } else {
2617: PetscInt i;
2618: for (i = 0; i < gDof; i++) {
2619: PetscInt gIndFine = rowIndices[i];
2620: if (gIndFine >= 0) {
2621: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2622: dnnz[gIndFine - rowStart] = numD[0];
2623: onnz[gIndFine - rowStart] = numO[0];
2624: }
2625: }
2626: }
2627: }
2628: } else { /* interpolate from all */
2629: if (numFields) {
2630: PetscInt f;
2631: for (f = 0; f < numFields; f++) {
2632: PetscInt numRows = offsets[f + 1] - offsets[f], row;
2633: for (row = 0; row < numRows; row++) {
2634: PetscInt gIndFine = rowIndices[offsets[f] + row];
2635: if (gIndFine >= 0) {
2636: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2637: dnnz[gIndFine - rowStart] = numD[f];
2638: onnz[gIndFine - rowStart] = numO[f];
2639: }
2640: }
2641: }
2642: } else { /* every dof get a full row */
2643: PetscInt i;
2644: for (i = 0; i < gDof; i++) {
2645: PetscInt gIndFine = rowIndices[i];
2646: if (gIndFine >= 0) {
2647: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2648: dnnz[gIndFine - rowStart] = numD[0];
2649: onnz[gIndFine - rowStart] = numO[0];
2650: }
2651: }
2652: }
2653: }
2654: }
2655: PetscCall(MatXAIJSetPreallocation(mat, 1, dnnz, onnz, NULL, NULL));
2656: PetscCall(PetscFree2(dnnz, onnz));
2658: PetscCall(DMPlexGetReferenceTree(fine, &refTree));
2659: PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
2660: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
2661: PetscCall(DMPlexGetAnchors(refTree, &refAnSec, NULL));
2662: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
2663: PetscCall(PetscSectionGetMaxDof(refConSec, &maxConDof));
2664: PetscCall(PetscSectionGetMaxDof(leafIndicesSec, &maxColumns));
2665: PetscCall(PetscMalloc1(maxConDof * maxColumns, &pointWork));
2666: for (p = leafStart; p < leafEnd; p++) {
2667: PetscInt gDof, gcDof, gOff;
2668: PetscInt numColIndices, pIndOff, *pInd;
2669: PetscInt matSize;
2670: PetscInt childId;
2672: PetscCall(PetscSectionGetDof(globalFine, p, &gDof));
2673: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &gcDof));
2674: if ((gDof - gcDof) <= 0) continue;
2675: childId = childIds[p - pStartF];
2676: PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
2677: PetscCall(PetscSectionGetDof(leafIndicesSec, p, &numColIndices));
2678: PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &pIndOff));
2679: numColIndices -= 2 * numFields;
2680: pInd = &leafIndices[pIndOff];
2681: offsets[0] = 0;
2682: offsetsCopy[0] = 0;
2683: newOffsets[0] = 0;
2684: newOffsetsCopy[0] = 0;
2685: rowOffsets[0] = 0;
2686: if (numFields) {
2687: PetscInt f;
2688: for (f = 0; f < numFields; f++) {
2689: PetscInt rowDof;
2691: PetscCall(PetscSectionGetFieldDof(localFine, p, f, &rowDof));
2692: offsets[f + 1] = offsets[f] + rowDof;
2693: offsetsCopy[f + 1] = offsets[f + 1];
2694: rowOffsets[f + 1] = pInd[numColIndices + f];
2695: newOffsets[f + 1] = pInd[numColIndices + numFields + f];
2696: }
2697: PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices));
2698: } else {
2699: PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices));
2700: }
2701: PetscCall(PetscSectionGetDof(leafMatricesSec, p, &matSize));
2702: if (!matSize) { /* incoming matrix is identity */
2703: if (childId < 0) { /* no child interpolation: scatter */
2704: if (numFields) {
2705: PetscInt f;
2706: for (f = 0; f < numFields; f++) {
2707: PetscInt numRows = offsets[f + 1] - offsets[f], row;
2708: for (row = 0; row < numRows; row++) PetscCall(MatSetValue(mat, rowIndices[offsets[f] + row], pInd[newOffsets[f] + row], 1., INSERT_VALUES));
2709: }
2710: } else {
2711: PetscInt numRows = gDof, row;
2712: for (row = 0; row < numRows; row++) PetscCall(MatSetValue(mat, rowIndices[row], pInd[row], 1., INSERT_VALUES));
2713: }
2714: } else { /* interpolate from all */
2715: if (numFields) {
2716: PetscInt f;
2717: for (f = 0; f < numFields; f++) {
2718: PetscInt numRows = offsets[f + 1] - offsets[f];
2719: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2720: PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], refPointFieldMats[childId - pRefStart][f], INSERT_VALUES));
2721: }
2722: } else {
2723: PetscCall(MatSetValues(mat, gDof, rowIndices, numColIndices, pInd, refPointFieldMats[childId - pRefStart][0], INSERT_VALUES));
2724: }
2725: }
2726: } else { /* interpolate from all */
2727: PetscInt pMatOff;
2728: PetscScalar *pMat;
2730: PetscCall(PetscSectionGetOffset(leafMatricesSec, p, &pMatOff));
2731: pMat = &leafMatrices[pMatOff];
2732: if (childId < 0) { /* copy the incoming matrix */
2733: if (numFields) {
2734: PetscInt f, count;
2735: for (f = 0, count = 0; f < numFields; f++) {
2736: PetscInt numRows = offsets[f + 1] - offsets[f];
2737: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2738: PetscInt numInRows = rowOffsets[f + 1] - rowOffsets[f];
2739: PetscScalar *inMat = &pMat[count];
2741: PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], inMat, INSERT_VALUES));
2742: count += numCols * numInRows;
2743: }
2744: } else {
2745: PetscCall(MatSetValues(mat, gDof, rowIndices, numColIndices, pInd, pMat, INSERT_VALUES));
2746: }
2747: } else { /* multiply the incoming matrix by the child interpolation */
2748: if (numFields) {
2749: PetscInt f, count;
2750: for (f = 0, count = 0; f < numFields; f++) {
2751: PetscInt numRows = offsets[f + 1] - offsets[f];
2752: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2753: PetscInt numInRows = rowOffsets[f + 1] - rowOffsets[f];
2754: PetscScalar *inMat = &pMat[count];
2755: PetscInt i, j, k;
2756: PetscCheck(refPointFieldN[childId - pRefStart][f] == numInRows, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point constraint matrix multiply dimension mismatch");
2757: for (i = 0; i < numRows; i++) {
2758: for (j = 0; j < numCols; j++) {
2759: PetscScalar val = 0.;
2760: for (k = 0; k < numInRows; k++) val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j];
2761: pointWork[i * numCols + j] = val;
2762: }
2763: }
2764: PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], pointWork, INSERT_VALUES));
2765: count += numCols * numInRows;
2766: }
2767: } else { /* every dof gets a full row */
2768: PetscInt numRows = gDof;
2769: PetscInt numCols = numColIndices;
2770: PetscInt numInRows = matSize / numColIndices;
2771: PetscInt i, j, k;
2772: PetscCheck(refPointFieldN[childId - pRefStart][0] == numInRows, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point constraint matrix multiply dimension mismatch");
2773: for (i = 0; i < numRows; i++) {
2774: for (j = 0; j < numCols; j++) {
2775: PetscScalar val = 0.;
2776: for (k = 0; k < numInRows; k++) val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j];
2777: pointWork[i * numCols + j] = val;
2778: }
2779: }
2780: PetscCall(MatSetValues(mat, numRows, rowIndices, numCols, pInd, pointWork, INSERT_VALUES));
2781: }
2782: }
2783: }
2784: }
2785: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
2786: PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
2787: PetscCall(PetscFree(pointWork));
2788: }
2789: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
2790: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
2791: PetscCall(PetscSectionDestroy(&leafIndicesSec));
2792: PetscCall(PetscSectionDestroy(&leafMatricesSec));
2793: PetscCall(PetscFree2(leafIndices, leafMatrices));
2794: PetscCall(PetscFree2(*(PetscInt ****)&perms, *(PetscScalar ****)&flips));
2795: PetscCall(PetscFree7(offsets, offsetsCopy, newOffsets, newOffsetsCopy, rowOffsets, numD, numO));
2796: PetscCall(ISRestoreIndices(aIS, &anchors));
2797: PetscFunctionReturn(PETSC_SUCCESS);
2798: }
2800: /*
2801: * Assuming a nodal basis (w.r.t. the dual basis) basis:
2802: *
2803: * for each coarse dof \phi^c_i:
2804: * for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i:
2805: * for each fine dof \phi^f_j;
2806: * a_{i,j} = 0;
2807: * for each fine dof \phi^f_k:
2808: * a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l
2809: * [^^^ this is = \phi^c_i ^^^]
2810: */
2811: PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj)
2812: {
2813: PetscDS ds;
2814: PetscSection section, cSection;
2815: DMLabel canonical, depth;
2816: Mat cMat, mat;
2817: PetscInt *nnz;
2818: PetscInt f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd;
2819: PetscInt m, n;
2820: PetscScalar *pointScalar;
2821: PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent;
2823: PetscFunctionBegin;
2824: PetscCall(DMGetLocalSection(refTree, §ion));
2825: PetscCall(DMGetDimension(refTree, &dim));
2826: PetscCall(PetscMalloc6(dim, &v0, dim, &v0parent, dim, &vtmp, dim * dim, &J, dim * dim, &Jparent, dim * dim, &invJ));
2827: PetscCall(PetscMalloc2(dim, &pointScalar, dim, &pointRef));
2828: PetscCall(DMGetDS(refTree, &ds));
2829: PetscCall(PetscDSGetNumFields(ds, &numFields));
2830: PetscCall(PetscSectionGetNumFields(section, &numSecFields));
2831: PetscCall(DMGetLabel(refTree, "canonical", &canonical));
2832: PetscCall(DMGetLabel(refTree, "depth", &depth));
2833: PetscCall(DMGetDefaultConstraints(refTree, &cSection, &cMat, NULL));
2834: PetscCall(DMPlexGetChart(refTree, &pStart, &pEnd));
2835: PetscCall(DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd));
2836: PetscCall(MatGetSize(cMat, &n, &m)); /* the injector has transpose sizes from the constraint matrix */
2837: /* Step 1: compute non-zero pattern. A proper subset of constraint matrix non-zero */
2838: PetscCall(PetscCalloc1(m, &nnz));
2839: for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */
2840: const PetscInt *children;
2841: PetscInt numChildren;
2842: PetscInt i, numChildDof, numSelfDof;
2844: if (canonical) {
2845: PetscInt pCanonical;
2846: PetscCall(DMLabelGetValue(canonical, p, &pCanonical));
2847: if (p != pCanonical) continue;
2848: }
2849: PetscCall(DMPlexGetTreeChildren(refTree, p, &numChildren, &children));
2850: if (!numChildren) continue;
2851: for (i = 0, numChildDof = 0; i < numChildren; i++) {
2852: PetscInt child = children[i];
2853: PetscInt dof;
2855: PetscCall(PetscSectionGetDof(section, child, &dof));
2856: numChildDof += dof;
2857: }
2858: PetscCall(PetscSectionGetDof(section, p, &numSelfDof));
2859: if (!numChildDof || !numSelfDof) continue;
2860: for (f = 0; f < numFields; f++) {
2861: PetscInt selfOff;
2863: if (numSecFields) { /* count the dofs for just this field */
2864: for (i = 0, numChildDof = 0; i < numChildren; i++) {
2865: PetscInt child = children[i];
2866: PetscInt dof;
2868: PetscCall(PetscSectionGetFieldDof(section, child, f, &dof));
2869: numChildDof += dof;
2870: }
2871: PetscCall(PetscSectionGetFieldDof(section, p, f, &numSelfDof));
2872: PetscCall(PetscSectionGetFieldOffset(section, p, f, &selfOff));
2873: } else {
2874: PetscCall(PetscSectionGetOffset(section, p, &selfOff));
2875: }
2876: for (i = 0; i < numSelfDof; i++) nnz[selfOff + i] = numChildDof;
2877: }
2878: }
2879: PetscCall(MatCreateAIJ(PETSC_COMM_SELF, m, n, m, n, -1, nnz, -1, NULL, &mat));
2880: PetscCall(PetscFree(nnz));
2881: /* Setp 2: compute entries */
2882: for (p = pStart; p < pEnd; p++) {
2883: const PetscInt *children;
2884: PetscInt numChildren;
2885: PetscInt i, numChildDof, numSelfDof;
2887: /* same conditions about when entries occur */
2888: if (canonical) {
2889: PetscInt pCanonical;
2890: PetscCall(DMLabelGetValue(canonical, p, &pCanonical));
2891: if (p != pCanonical) continue;
2892: }
2893: PetscCall(DMPlexGetTreeChildren(refTree, p, &numChildren, &children));
2894: if (!numChildren) continue;
2895: for (i = 0, numChildDof = 0; i < numChildren; i++) {
2896: PetscInt child = children[i];
2897: PetscInt dof;
2899: PetscCall(PetscSectionGetDof(section, child, &dof));
2900: numChildDof += dof;
2901: }
2902: PetscCall(PetscSectionGetDof(section, p, &numSelfDof));
2903: if (!numChildDof || !numSelfDof) continue;
2905: for (f = 0; f < numFields; f++) {
2906: PetscInt pI = -1, cI = -1;
2907: PetscInt selfOff, Nc, parentCell;
2908: PetscInt cellShapeOff;
2909: PetscObject disc;
2910: PetscDualSpace dsp;
2911: PetscClassId classId;
2912: PetscScalar *pointMat;
2913: PetscInt *matRows, *matCols;
2914: PetscInt pO = PETSC_MIN_INT;
2915: const PetscInt *depthNumDof;
2917: if (numSecFields) {
2918: for (i = 0, numChildDof = 0; i < numChildren; i++) {
2919: PetscInt child = children[i];
2920: PetscInt dof;
2922: PetscCall(PetscSectionGetFieldDof(section, child, f, &dof));
2923: numChildDof += dof;
2924: }
2925: PetscCall(PetscSectionGetFieldDof(section, p, f, &numSelfDof));
2926: PetscCall(PetscSectionGetFieldOffset(section, p, f, &selfOff));
2927: } else {
2928: PetscCall(PetscSectionGetOffset(section, p, &selfOff));
2929: }
2931: /* find a cell whose closure contains p */
2932: if (p >= cStart && p < cEnd) {
2933: parentCell = p;
2934: } else {
2935: PetscInt *star = NULL;
2936: PetscInt numStar;
2938: parentCell = -1;
2939: PetscCall(DMPlexGetTransitiveClosure(refTree, p, PETSC_FALSE, &numStar, &star));
2940: for (i = numStar - 1; i >= 0; i--) {
2941: PetscInt c = star[2 * i];
2943: if (c >= cStart && c < cEnd) {
2944: parentCell = c;
2945: break;
2946: }
2947: }
2948: PetscCall(DMPlexRestoreTransitiveClosure(refTree, p, PETSC_FALSE, &numStar, &star));
2949: }
2950: /* determine the offset of p's shape functions within parentCell's shape functions */
2951: PetscCall(PetscDSGetDiscretization(ds, f, &disc));
2952: PetscCall(PetscObjectGetClassId(disc, &classId));
2953: if (classId == PETSCFE_CLASSID) {
2954: PetscCall(PetscFEGetDualSpace((PetscFE)disc, &dsp));
2955: } else if (classId == PETSCFV_CLASSID) {
2956: PetscCall(PetscFVGetDualSpace((PetscFV)disc, &dsp));
2957: } else {
2958: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported discretization object");
2959: }
2960: PetscCall(PetscDualSpaceGetNumDof(dsp, &depthNumDof));
2961: PetscCall(PetscDualSpaceGetNumComponents(dsp, &Nc));
2962: {
2963: PetscInt *closure = NULL;
2964: PetscInt numClosure;
2966: PetscCall(DMPlexGetTransitiveClosure(refTree, parentCell, PETSC_TRUE, &numClosure, &closure));
2967: for (i = 0, pI = -1, cellShapeOff = 0; i < numClosure; i++) {
2968: PetscInt point = closure[2 * i], pointDepth;
2970: pO = closure[2 * i + 1];
2971: if (point == p) {
2972: pI = i;
2973: break;
2974: }
2975: PetscCall(DMLabelGetValue(depth, point, &pointDepth));
2976: cellShapeOff += depthNumDof[pointDepth];
2977: }
2978: PetscCall(DMPlexRestoreTransitiveClosure(refTree, parentCell, PETSC_TRUE, &numClosure, &closure));
2979: }
2981: PetscCall(DMGetWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR, &pointMat));
2982: PetscCall(DMGetWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT, &matRows));
2983: matCols = matRows + numSelfDof;
2984: for (i = 0; i < numSelfDof; i++) matRows[i] = selfOff + i;
2985: for (i = 0; i < numSelfDof * numChildDof; i++) pointMat[i] = 0.;
2986: {
2987: PetscInt colOff = 0;
2989: for (i = 0; i < numChildren; i++) {
2990: PetscInt child = children[i];
2991: PetscInt dof, off, j;
2993: if (numSecFields) {
2994: PetscCall(PetscSectionGetFieldDof(cSection, child, f, &dof));
2995: PetscCall(PetscSectionGetFieldOffset(cSection, child, f, &off));
2996: } else {
2997: PetscCall(PetscSectionGetDof(cSection, child, &dof));
2998: PetscCall(PetscSectionGetOffset(cSection, child, &off));
2999: }
3001: for (j = 0; j < dof; j++) matCols[colOff++] = off + j;
3002: }
3003: }
3004: if (classId == PETSCFE_CLASSID) {
3005: PetscFE fe = (PetscFE)disc;
3006: PetscInt fSize;
3007: const PetscInt ***perms;
3008: const PetscScalar ***flips;
3009: const PetscInt *pperms;
3011: PetscCall(PetscFEGetDualSpace(fe, &dsp));
3012: PetscCall(PetscDualSpaceGetDimension(dsp, &fSize));
3013: PetscCall(PetscDualSpaceGetSymmetries(dsp, &perms, &flips));
3014: pperms = perms ? perms[pI] ? perms[pI][pO] : NULL : NULL;
3015: for (i = 0; i < numSelfDof; i++) { /* for every shape function */
3016: PetscQuadrature q;
3017: PetscInt dim, thisNc, numPoints, j, k;
3018: const PetscReal *points;
3019: const PetscReal *weights;
3020: PetscInt *closure = NULL;
3021: PetscInt numClosure;
3022: PetscInt iCell = pperms ? pperms[i] : i;
3023: PetscInt parentCellShapeDof = cellShapeOff + iCell;
3024: PetscTabulation Tparent;
3026: PetscCall(PetscDualSpaceGetFunctional(dsp, parentCellShapeDof, &q));
3027: PetscCall(PetscQuadratureGetData(q, &dim, &thisNc, &numPoints, &points, &weights));
3028: PetscCheck(thisNc == Nc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Functional dim %" PetscInt_FMT " does not much basis dim %" PetscInt_FMT, thisNc, Nc);
3029: PetscCall(PetscFECreateTabulation(fe, 1, numPoints, points, 0, &Tparent)); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */
3030: for (j = 0; j < numPoints; j++) {
3031: PetscInt childCell = -1;
3032: PetscReal *parentValAtPoint;
3033: const PetscReal xi0[3] = {-1., -1., -1.};
3034: const PetscReal *pointReal = &points[dim * j];
3035: const PetscScalar *point;
3036: PetscTabulation Tchild;
3037: PetscInt childCellShapeOff, pointMatOff;
3038: #if defined(PETSC_USE_COMPLEX)
3039: PetscInt d;
3041: for (d = 0; d < dim; d++) pointScalar[d] = points[dim * j + d];
3042: point = pointScalar;
3043: #else
3044: point = pointReal;
3045: #endif
3047: parentValAtPoint = &Tparent->T[0][(fSize * j + parentCellShapeDof) * Nc];
3049: for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/
3050: PetscInt child = children[k];
3051: PetscInt *star = NULL;
3052: PetscInt numStar, s;
3054: PetscCall(DMPlexGetTransitiveClosure(refTree, child, PETSC_FALSE, &numStar, &star));
3055: for (s = numStar - 1; s >= 0; s--) {
3056: PetscInt c = star[2 * s];
3058: if (c < cStart || c >= cEnd) continue;
3059: PetscCall(DMPlexLocatePoint_Internal(refTree, dim, point, c, &childCell));
3060: if (childCell >= 0) break;
3061: }
3062: PetscCall(DMPlexRestoreTransitiveClosure(refTree, child, PETSC_FALSE, &numStar, &star));
3063: if (childCell >= 0) break;
3064: }
3065: PetscCheck(childCell >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not locate quadrature point");
3066: PetscCall(DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ));
3067: PetscCall(DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent));
3068: CoordinatesRefToReal(dim, dim, xi0, v0parent, Jparent, pointReal, vtmp);
3069: CoordinatesRealToRef(dim, dim, xi0, v0, invJ, vtmp, pointRef);
3071: PetscCall(PetscFECreateTabulation(fe, 1, 1, pointRef, 0, &Tchild));
3072: PetscCall(DMPlexGetTransitiveClosure(refTree, childCell, PETSC_TRUE, &numClosure, &closure));
3073: for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */
3074: PetscInt child = children[k], childDepth, childDof, childO = PETSC_MIN_INT;
3075: PetscInt l;
3076: const PetscInt *cperms;
3078: PetscCall(DMLabelGetValue(depth, child, &childDepth));
3079: childDof = depthNumDof[childDepth];
3080: for (l = 0, cI = -1, childCellShapeOff = 0; l < numClosure; l++) {
3081: PetscInt point = closure[2 * l];
3082: PetscInt pointDepth;
3084: childO = closure[2 * l + 1];
3085: if (point == child) {
3086: cI = l;
3087: break;
3088: }
3089: PetscCall(DMLabelGetValue(depth, point, &pointDepth));
3090: childCellShapeOff += depthNumDof[pointDepth];
3091: }
3092: if (l == numClosure) {
3093: pointMatOff += childDof;
3094: continue; /* child is not in the closure of the cell: has nothing to contribute to this point */
3095: }
3096: cperms = perms ? perms[cI] ? perms[cI][childO] : NULL : NULL;
3097: for (l = 0; l < childDof; l++) {
3098: PetscInt lCell = cperms ? cperms[l] : l;
3099: PetscInt childCellDof = childCellShapeOff + lCell;
3100: PetscReal *childValAtPoint;
3101: PetscReal val = 0.;
3103: childValAtPoint = &Tchild->T[0][childCellDof * Nc];
3104: for (m = 0; m < Nc; m++) val += weights[j * Nc + m] * parentValAtPoint[m] * childValAtPoint[m];
3106: pointMat[i * numChildDof + pointMatOff + l] += val;
3107: }
3108: pointMatOff += childDof;
3109: }
3110: PetscCall(DMPlexRestoreTransitiveClosure(refTree, childCell, PETSC_TRUE, &numClosure, &closure));
3111: PetscCall(PetscTabulationDestroy(&Tchild));
3112: }
3113: PetscCall(PetscTabulationDestroy(&Tparent));
3114: }
3115: } else { /* just the volume-weighted averages of the children */
3116: PetscReal parentVol;
3117: PetscInt childCell;
3119: PetscCall(DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL));
3120: for (i = 0, childCell = 0; i < numChildren; i++) {
3121: PetscInt child = children[i], j;
3122: PetscReal childVol;
3124: if (child < cStart || child >= cEnd) continue;
3125: PetscCall(DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL));
3126: for (j = 0; j < Nc; j++) pointMat[j * numChildDof + Nc * childCell + j] = childVol / parentVol;
3127: childCell++;
3128: }
3129: }
3130: /* Insert pointMat into mat */
3131: PetscCall(MatSetValues(mat, numSelfDof, matRows, numChildDof, matCols, pointMat, INSERT_VALUES));
3132: PetscCall(DMRestoreWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT, &matRows));
3133: PetscCall(DMRestoreWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR, &pointMat));
3134: }
3135: }
3136: PetscCall(PetscFree6(v0, v0parent, vtmp, J, Jparent, invJ));
3137: PetscCall(PetscFree2(pointScalar, pointRef));
3138: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
3139: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
3140: *inj = mat;
3141: PetscFunctionReturn(PETSC_SUCCESS);
3142: }
3144: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3145: {
3146: PetscDS ds;
3147: PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof;
3148: PetscScalar ***refPointFieldMats;
3149: PetscSection refConSec, refSection;
3151: PetscFunctionBegin;
3152: PetscCall(DMGetDS(refTree, &ds));
3153: PetscCall(PetscDSGetNumFields(ds, &numFields));
3154: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
3155: PetscCall(DMGetLocalSection(refTree, &refSection));
3156: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
3157: PetscCall(PetscMalloc1(pRefEnd - pRefStart, &refPointFieldMats));
3158: PetscCall(PetscSectionGetMaxDof(refConSec, &maxDof));
3159: PetscCall(PetscMalloc1(maxDof, &rows));
3160: PetscCall(PetscMalloc1(maxDof * maxDof, &cols));
3161: for (p = pRefStart; p < pRefEnd; p++) {
3162: PetscInt parent, pDof, parentDof;
3164: PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL));
3165: PetscCall(PetscSectionGetDof(refConSec, p, &pDof));
3166: PetscCall(PetscSectionGetDof(refSection, parent, &parentDof));
3167: if (!pDof || !parentDof || parent == p) continue;
3169: PetscCall(PetscMalloc1(numFields, &refPointFieldMats[p - pRefStart]));
3170: for (f = 0; f < numFields; f++) {
3171: PetscInt cDof, cOff, numCols, r;
3173: if (numFields > 1) {
3174: PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof));
3175: PetscCall(PetscSectionGetFieldOffset(refConSec, p, f, &cOff));
3176: } else {
3177: PetscCall(PetscSectionGetDof(refConSec, p, &cDof));
3178: PetscCall(PetscSectionGetOffset(refConSec, p, &cOff));
3179: }
3181: for (r = 0; r < cDof; r++) rows[r] = cOff + r;
3182: numCols = 0;
3183: {
3184: PetscInt aDof, aOff, j;
3186: if (numFields > 1) {
3187: PetscCall(PetscSectionGetFieldDof(refSection, parent, f, &aDof));
3188: PetscCall(PetscSectionGetFieldOffset(refSection, parent, f, &aOff));
3189: } else {
3190: PetscCall(PetscSectionGetDof(refSection, parent, &aDof));
3191: PetscCall(PetscSectionGetOffset(refSection, parent, &aOff));
3192: }
3194: for (j = 0; j < aDof; j++) cols[numCols++] = aOff + j;
3195: }
3196: PetscCall(PetscMalloc1(cDof * numCols, &refPointFieldMats[p - pRefStart][f]));
3197: /* transpose of constraint matrix */
3198: PetscCall(MatGetValues(inj, numCols, cols, cDof, rows, refPointFieldMats[p - pRefStart][f]));
3199: }
3200: }
3201: *childrenMats = refPointFieldMats;
3202: PetscCall(PetscFree(rows));
3203: PetscCall(PetscFree(cols));
3204: PetscFunctionReturn(PETSC_SUCCESS);
3205: }
3207: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3208: {
3209: PetscDS ds;
3210: PetscScalar ***refPointFieldMats;
3211: PetscInt numFields, pRefStart, pRefEnd, p, f;
3212: PetscSection refConSec, refSection;
3214: PetscFunctionBegin;
3215: refPointFieldMats = *childrenMats;
3216: *childrenMats = NULL;
3217: PetscCall(DMGetDS(refTree, &ds));
3218: PetscCall(DMGetLocalSection(refTree, &refSection));
3219: PetscCall(PetscDSGetNumFields(ds, &numFields));
3220: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
3221: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
3222: for (p = pRefStart; p < pRefEnd; p++) {
3223: PetscInt parent, pDof, parentDof;
3225: PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL));
3226: PetscCall(PetscSectionGetDof(refConSec, p, &pDof));
3227: PetscCall(PetscSectionGetDof(refSection, parent, &parentDof));
3228: if (!pDof || !parentDof || parent == p) continue;
3230: for (f = 0; f < numFields; f++) {
3231: PetscInt cDof;
3233: if (numFields > 1) {
3234: PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof));
3235: } else {
3236: PetscCall(PetscSectionGetDof(refConSec, p, &cDof));
3237: }
3239: PetscCall(PetscFree(refPointFieldMats[p - pRefStart][f]));
3240: }
3241: PetscCall(PetscFree(refPointFieldMats[p - pRefStart]));
3242: }
3243: PetscCall(PetscFree(refPointFieldMats));
3244: PetscFunctionReturn(PETSC_SUCCESS);
3245: }
3247: static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree, Mat *injRef)
3248: {
3249: Mat cMatRef;
3250: PetscObject injRefObj;
3252: PetscFunctionBegin;
3253: PetscCall(DMGetDefaultConstraints(refTree, NULL, &cMatRef, NULL));
3254: PetscCall(PetscObjectQuery((PetscObject)cMatRef, "DMPlexComputeInjectorTree_refTree", &injRefObj));
3255: *injRef = (Mat)injRefObj;
3256: if (!*injRef) {
3257: PetscCall(DMPlexComputeInjectorReferenceTree(refTree, injRef));
3258: PetscCall(PetscObjectCompose((PetscObject)cMatRef, "DMPlexComputeInjectorTree_refTree", (PetscObject)*injRef));
3259: /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */
3260: PetscCall(PetscObjectDereference((PetscObject)*injRef));
3261: }
3262: PetscFunctionReturn(PETSC_SUCCESS);
3263: }
3265: static PetscErrorCode DMPlexTransferInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, const PetscInt *childIds, Vec fineVec, PetscInt numFields, PetscInt *offsets, PetscSection *rootMultiSec, PetscSection *multiLeafSec, PetscInt **gatheredIndices, PetscScalar **gatheredValues)
3266: {
3267: PetscInt pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti;
3268: PetscSection globalCoarse, globalFine;
3269: PetscSection localCoarse, localFine, leafIndicesSec;
3270: PetscSection multiRootSec, rootIndicesSec;
3271: PetscInt *leafInds, *rootInds = NULL;
3272: const PetscInt *rootDegrees;
3273: PetscScalar *leafVals = NULL, *rootVals = NULL;
3274: PetscSF coarseToFineEmbedded;
3276: PetscFunctionBegin;
3277: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
3278: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
3279: PetscCall(DMGetLocalSection(fine, &localFine));
3280: PetscCall(DMGetGlobalSection(fine, &globalFine));
3281: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafIndicesSec));
3282: PetscCall(PetscSectionSetChart(leafIndicesSec, pStartF, pEndF));
3283: PetscCall(PetscSectionGetMaxDof(localFine, &maxDof));
3284: { /* winnow fine points that don't have global dofs out of the sf */
3285: PetscInt l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices;
3286: const PetscInt *leaves;
3288: PetscCall(PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL));
3289: for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
3290: p = leaves ? leaves[l] : l;
3291: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
3292: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
3293: if ((dof - cdof) > 0) {
3294: numPointsWithDofs++;
3296: PetscCall(PetscSectionGetDof(localFine, p, &dof));
3297: PetscCall(PetscSectionSetDof(leafIndicesSec, p, dof + 1));
3298: }
3299: }
3300: PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
3301: PetscCall(PetscSectionSetUp(leafIndicesSec));
3302: PetscCall(PetscSectionGetStorageSize(leafIndicesSec, &numIndices));
3303: PetscCall(PetscMalloc1(gatheredIndices ? numIndices : (maxDof + 1), &leafInds));
3304: if (gatheredValues) PetscCall(PetscMalloc1(numIndices, &leafVals));
3305: for (l = 0, offset = 0; l < nleaves; l++) {
3306: p = leaves ? leaves[l] : l;
3307: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
3308: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
3309: if ((dof - cdof) > 0) {
3310: PetscInt off, gOff;
3311: PetscInt *pInd;
3312: PetscScalar *pVal = NULL;
3314: pointsWithDofs[offset++] = l;
3316: PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &off));
3318: pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds;
3319: if (gatheredValues) {
3320: PetscInt i;
3322: pVal = &leafVals[off + 1];
3323: for (i = 0; i < dof; i++) pVal[i] = 0.;
3324: }
3325: PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
3327: offsets[0] = 0;
3328: if (numFields) {
3329: PetscInt f;
3331: for (f = 0; f < numFields; f++) {
3332: PetscInt fDof;
3333: PetscCall(PetscSectionGetFieldDof(localFine, p, f, &fDof));
3334: offsets[f + 1] = fDof + offsets[f];
3335: }
3336: PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, -1, NULL, pInd));
3337: } else {
3338: PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, NULL, pInd));
3339: }
3340: if (gatheredValues) PetscCall(VecGetValues(fineVec, dof, pInd, pVal));
3341: }
3342: }
3343: PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded));
3344: PetscCall(PetscFree(pointsWithDofs));
3345: }
3347: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
3348: PetscCall(DMGetLocalSection(coarse, &localCoarse));
3349: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
3351: { /* there may be the case where an sf root has a parent: broadcast parents back to children */
3352: MPI_Datatype threeInt;
3353: PetscMPIInt rank;
3354: PetscInt(*parentNodeAndIdCoarse)[3];
3355: PetscInt(*parentNodeAndIdFine)[3];
3356: PetscInt p, nleaves, nleavesToParents;
3357: PetscSF pointSF, sfToParents;
3358: const PetscInt *ilocal;
3359: const PetscSFNode *iremote;
3360: PetscSFNode *iremoteToParents;
3361: PetscInt *ilocalToParents;
3363: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank));
3364: PetscCallMPI(MPI_Type_contiguous(3, MPIU_INT, &threeInt));
3365: PetscCallMPI(MPI_Type_commit(&threeInt));
3366: PetscCall(PetscMalloc2(pEndC - pStartC, &parentNodeAndIdCoarse, pEndF - pStartF, &parentNodeAndIdFine));
3367: PetscCall(DMGetPointSF(coarse, &pointSF));
3368: PetscCall(PetscSFGetGraph(pointSF, NULL, &nleaves, &ilocal, &iremote));
3369: for (p = pStartC; p < pEndC; p++) {
3370: PetscInt parent, childId;
3371: PetscCall(DMPlexGetTreeParent(coarse, p, &parent, &childId));
3372: parentNodeAndIdCoarse[p - pStartC][0] = rank;
3373: parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC;
3374: parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId;
3375: if (nleaves > 0) {
3376: PetscInt leaf = -1;
3378: if (ilocal) {
3379: PetscCall(PetscFindInt(parent, nleaves, ilocal, &leaf));
3380: } else {
3381: leaf = p - pStartC;
3382: }
3383: if (leaf >= 0) {
3384: parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank;
3385: parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index;
3386: }
3387: }
3388: }
3389: for (p = pStartF; p < pEndF; p++) {
3390: parentNodeAndIdFine[p - pStartF][0] = -1;
3391: parentNodeAndIdFine[p - pStartF][1] = -1;
3392: parentNodeAndIdFine[p - pStartF][2] = -1;
3393: }
3394: PetscCall(PetscSFBcastBegin(coarseToFineEmbedded, threeInt, parentNodeAndIdCoarse, parentNodeAndIdFine, MPI_REPLACE));
3395: PetscCall(PetscSFBcastEnd(coarseToFineEmbedded, threeInt, parentNodeAndIdCoarse, parentNodeAndIdFine, MPI_REPLACE));
3396: for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3397: PetscInt dof;
3399: PetscCall(PetscSectionGetDof(leafIndicesSec, p, &dof));
3400: if (dof) {
3401: PetscInt off;
3403: PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &off));
3404: if (gatheredIndices) {
3405: leafInds[off] = PetscMax(childIds[p - pStartF], parentNodeAndIdFine[p - pStartF][2]);
3406: } else if (gatheredValues) {
3407: leafVals[off] = (PetscScalar)PetscMax(childIds[p - pStartF], parentNodeAndIdFine[p - pStartF][2]);
3408: }
3409: }
3410: if (parentNodeAndIdFine[p - pStartF][0] >= 0) nleavesToParents++;
3411: }
3412: PetscCall(PetscMalloc1(nleavesToParents, &ilocalToParents));
3413: PetscCall(PetscMalloc1(nleavesToParents, &iremoteToParents));
3414: for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3415: if (parentNodeAndIdFine[p - pStartF][0] >= 0) {
3416: ilocalToParents[nleavesToParents] = p - pStartF;
3417: iremoteToParents[nleavesToParents].rank = parentNodeAndIdFine[p - pStartF][0];
3418: iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p - pStartF][1];
3419: nleavesToParents++;
3420: }
3421: }
3422: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)coarse), &sfToParents));
3423: PetscCall(PetscSFSetGraph(sfToParents, pEndC - pStartC, nleavesToParents, ilocalToParents, PETSC_OWN_POINTER, iremoteToParents, PETSC_OWN_POINTER));
3424: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
3426: coarseToFineEmbedded = sfToParents;
3428: PetscCall(PetscFree2(parentNodeAndIdCoarse, parentNodeAndIdFine));
3429: PetscCallMPI(MPI_Type_free(&threeInt));
3430: }
3432: { /* winnow out coarse points that don't have dofs */
3433: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;
3434: PetscSF sfDofsOnly;
3436: for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) {
3437: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3438: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
3439: if ((dof - cdof) > 0) numPointsWithDofs++;
3440: }
3441: PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
3442: for (p = pStartC, offset = 0; p < pEndC; p++) {
3443: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3444: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
3445: if ((dof - cdof) > 0) pointsWithDofs[offset++] = p - pStartC;
3446: }
3447: PetscCall(PetscSFCreateEmbeddedRootSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly));
3448: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
3449: PetscCall(PetscFree(pointsWithDofs));
3450: coarseToFineEmbedded = sfDofsOnly;
3451: }
3453: /* communicate back to the coarse mesh which coarse points have children (that may require injection) */
3454: PetscCall(PetscSFComputeDegreeBegin(coarseToFineEmbedded, &rootDegrees));
3455: PetscCall(PetscSFComputeDegreeEnd(coarseToFineEmbedded, &rootDegrees));
3456: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &multiRootSec));
3457: PetscCall(PetscSectionSetChart(multiRootSec, pStartC, pEndC));
3458: for (p = pStartC; p < pEndC; p++) PetscCall(PetscSectionSetDof(multiRootSec, p, rootDegrees[p - pStartC]));
3459: PetscCall(PetscSectionSetUp(multiRootSec));
3460: PetscCall(PetscSectionGetStorageSize(multiRootSec, &numMulti));
3461: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootIndicesSec));
3462: { /* distribute the leaf section */
3463: PetscSF multi, multiInv, indicesSF;
3464: PetscInt *remoteOffsets, numRootIndices;
3466: PetscCall(PetscSFGetMultiSF(coarseToFineEmbedded, &multi));
3467: PetscCall(PetscSFCreateInverseSF(multi, &multiInv));
3468: PetscCall(PetscSFDistributeSection(multiInv, leafIndicesSec, &remoteOffsets, rootIndicesSec));
3469: PetscCall(PetscSFCreateSectionSF(multiInv, leafIndicesSec, remoteOffsets, rootIndicesSec, &indicesSF));
3470: PetscCall(PetscFree(remoteOffsets));
3471: PetscCall(PetscSFDestroy(&multiInv));
3472: PetscCall(PetscSectionGetStorageSize(rootIndicesSec, &numRootIndices));
3473: if (gatheredIndices) {
3474: PetscCall(PetscMalloc1(numRootIndices, &rootInds));
3475: PetscCall(PetscSFBcastBegin(indicesSF, MPIU_INT, leafInds, rootInds, MPI_REPLACE));
3476: PetscCall(PetscSFBcastEnd(indicesSF, MPIU_INT, leafInds, rootInds, MPI_REPLACE));
3477: }
3478: if (gatheredValues) {
3479: PetscCall(PetscMalloc1(numRootIndices, &rootVals));
3480: PetscCall(PetscSFBcastBegin(indicesSF, MPIU_SCALAR, leafVals, rootVals, MPI_REPLACE));
3481: PetscCall(PetscSFBcastEnd(indicesSF, MPIU_SCALAR, leafVals, rootVals, MPI_REPLACE));
3482: }
3483: PetscCall(PetscSFDestroy(&indicesSF));
3484: }
3485: PetscCall(PetscSectionDestroy(&leafIndicesSec));
3486: PetscCall(PetscFree(leafInds));
3487: PetscCall(PetscFree(leafVals));
3488: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
3489: *rootMultiSec = multiRootSec;
3490: *multiLeafSec = rootIndicesSec;
3491: if (gatheredIndices) *gatheredIndices = rootInds;
3492: if (gatheredValues) *gatheredValues = rootVals;
3493: PetscFunctionReturn(PETSC_SUCCESS);
3494: }
3496: PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
3497: {
3498: DM refTree;
3499: PetscSection multiRootSec, rootIndicesSec;
3500: PetscSection globalCoarse, globalFine;
3501: PetscSection localCoarse, localFine;
3502: PetscSection cSecRef;
3503: PetscInt *rootIndices = NULL, *parentIndices, pRefStart, pRefEnd;
3504: Mat injRef;
3505: PetscInt numFields, maxDof;
3506: PetscInt pStartC, pEndC, pStartF, pEndF, p;
3507: PetscInt *offsets, *offsetsCopy, *rowOffsets;
3508: PetscLayout rowMap, colMap;
3509: PetscInt rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO;
3510: PetscScalar ***childrenMats = NULL; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */
3512: PetscFunctionBegin;
3514: /* get the templates for the fine-to-coarse injection from the reference tree */
3515: PetscCall(DMPlexGetReferenceTree(coarse, &refTree));
3516: PetscCall(DMGetDefaultConstraints(refTree, &cSecRef, NULL, NULL));
3517: PetscCall(PetscSectionGetChart(cSecRef, &pRefStart, &pRefEnd));
3518: PetscCall(DMPlexReferenceTreeGetInjector(refTree, &injRef));
3520: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
3521: PetscCall(DMGetLocalSection(fine, &localFine));
3522: PetscCall(DMGetGlobalSection(fine, &globalFine));
3523: PetscCall(PetscSectionGetNumFields(localFine, &numFields));
3524: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
3525: PetscCall(DMGetLocalSection(coarse, &localCoarse));
3526: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
3527: PetscCall(PetscSectionGetMaxDof(localCoarse, &maxDof));
3528: {
3529: PetscInt maxFields = PetscMax(1, numFields) + 1;
3530: PetscCall(PetscMalloc3(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &rowOffsets));
3531: }
3533: PetscCall(DMPlexTransferInjectorTree(coarse, fine, coarseToFine, childIds, NULL, numFields, offsets, &multiRootSec, &rootIndicesSec, &rootIndices, NULL));
3535: PetscCall(PetscMalloc1(maxDof, &parentIndices));
3537: /* count indices */
3538: PetscCall(MatGetLayouts(mat, &rowMap, &colMap));
3539: PetscCall(PetscLayoutSetUp(rowMap));
3540: PetscCall(PetscLayoutSetUp(colMap));
3541: PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd));
3542: PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd));
3543: PetscCall(PetscCalloc2(rowEnd - rowStart, &nnzD, rowEnd - rowStart, &nnzO));
3544: for (p = pStartC; p < pEndC; p++) {
3545: PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
3547: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3548: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
3549: if ((dof - cdof) <= 0) continue;
3550: PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));
3552: rowOffsets[0] = 0;
3553: offsetsCopy[0] = 0;
3554: if (numFields) {
3555: PetscInt f;
3557: for (f = 0; f < numFields; f++) {
3558: PetscInt fDof;
3559: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
3560: rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3561: }
3562: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices));
3563: } else {
3564: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices));
3565: rowOffsets[1] = offsetsCopy[0];
3566: }
3568: PetscCall(PetscSectionGetDof(multiRootSec, p, &numLeaves));
3569: PetscCall(PetscSectionGetOffset(multiRootSec, p, &leafStart));
3570: leafEnd = leafStart + numLeaves;
3571: for (l = leafStart; l < leafEnd; l++) {
3572: PetscInt numIndices, childId, offset;
3573: const PetscInt *childIndices;
3575: PetscCall(PetscSectionGetDof(rootIndicesSec, l, &numIndices));
3576: PetscCall(PetscSectionGetOffset(rootIndicesSec, l, &offset));
3577: childId = rootIndices[offset++];
3578: childIndices = &rootIndices[offset];
3579: numIndices--;
3581: if (childId == -1) { /* equivalent points: scatter */
3582: PetscInt i;
3584: for (i = 0; i < numIndices; i++) {
3585: PetscInt colIndex = childIndices[i];
3586: PetscInt rowIndex = parentIndices[i];
3587: if (rowIndex < 0) continue;
3588: PetscCheck(colIndex >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unconstrained fine and constrained coarse");
3589: if (colIndex >= colStart && colIndex < colEnd) {
3590: nnzD[rowIndex - rowStart] = 1;
3591: } else {
3592: nnzO[rowIndex - rowStart] = 1;
3593: }
3594: }
3595: } else {
3596: PetscInt parentId, f, lim;
3598: PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL));
3600: lim = PetscMax(1, numFields);
3601: offsets[0] = 0;
3602: if (numFields) {
3603: PetscInt f;
3605: for (f = 0; f < numFields; f++) {
3606: PetscInt fDof;
3607: PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof));
3609: offsets[f + 1] = fDof + offsets[f];
3610: }
3611: } else {
3612: PetscInt cDof;
3614: PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof));
3615: offsets[1] = cDof;
3616: }
3617: for (f = 0; f < lim; f++) {
3618: PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1];
3619: PetscInt childStart = offsets[f], childEnd = offsets[f + 1];
3620: PetscInt i, numD = 0, numO = 0;
3622: for (i = childStart; i < childEnd; i++) {
3623: PetscInt colIndex = childIndices[i];
3625: if (colIndex < 0) continue;
3626: if (colIndex >= colStart && colIndex < colEnd) {
3627: numD++;
3628: } else {
3629: numO++;
3630: }
3631: }
3632: for (i = parentStart; i < parentEnd; i++) {
3633: PetscInt rowIndex = parentIndices[i];
3635: if (rowIndex < 0) continue;
3636: nnzD[rowIndex - rowStart] += numD;
3637: nnzO[rowIndex - rowStart] += numO;
3638: }
3639: }
3640: }
3641: }
3642: }
3643: /* preallocate */
3644: PetscCall(MatXAIJSetPreallocation(mat, 1, nnzD, nnzO, NULL, NULL));
3645: PetscCall(PetscFree2(nnzD, nnzO));
3646: /* insert values */
3647: PetscCall(DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree, injRef, &childrenMats));
3648: for (p = pStartC; p < pEndC; p++) {
3649: PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
3651: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3652: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
3653: if ((dof - cdof) <= 0) continue;
3654: PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));
3656: rowOffsets[0] = 0;
3657: offsetsCopy[0] = 0;
3658: if (numFields) {
3659: PetscInt f;
3661: for (f = 0; f < numFields; f++) {
3662: PetscInt fDof;
3663: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
3664: rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3665: }
3666: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices));
3667: } else {
3668: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices));
3669: rowOffsets[1] = offsetsCopy[0];
3670: }
3672: PetscCall(PetscSectionGetDof(multiRootSec, p, &numLeaves));
3673: PetscCall(PetscSectionGetOffset(multiRootSec, p, &leafStart));
3674: leafEnd = leafStart + numLeaves;
3675: for (l = leafStart; l < leafEnd; l++) {
3676: PetscInt numIndices, childId, offset;
3677: const PetscInt *childIndices;
3679: PetscCall(PetscSectionGetDof(rootIndicesSec, l, &numIndices));
3680: PetscCall(PetscSectionGetOffset(rootIndicesSec, l, &offset));
3681: childId = rootIndices[offset++];
3682: childIndices = &rootIndices[offset];
3683: numIndices--;
3685: if (childId == -1) { /* equivalent points: scatter */
3686: PetscInt i;
3688: for (i = 0; i < numIndices; i++) PetscCall(MatSetValue(mat, parentIndices[i], childIndices[i], 1., INSERT_VALUES));
3689: } else {
3690: PetscInt parentId, f, lim;
3692: PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL));
3694: lim = PetscMax(1, numFields);
3695: offsets[0] = 0;
3696: if (numFields) {
3697: PetscInt f;
3699: for (f = 0; f < numFields; f++) {
3700: PetscInt fDof;
3701: PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof));
3703: offsets[f + 1] = fDof + offsets[f];
3704: }
3705: } else {
3706: PetscInt cDof;
3708: PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof));
3709: offsets[1] = cDof;
3710: }
3711: for (f = 0; f < lim; f++) {
3712: PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0];
3713: PetscInt *rowIndices = &parentIndices[rowOffsets[f]];
3714: const PetscInt *colIndices = &childIndices[offsets[f]];
3716: PetscCall(MatSetValues(mat, rowOffsets[f + 1] - rowOffsets[f], rowIndices, offsets[f + 1] - offsets[f], colIndices, childMat, INSERT_VALUES));
3717: }
3718: }
3719: }
3720: }
3721: PetscCall(PetscSectionDestroy(&multiRootSec));
3722: PetscCall(PetscSectionDestroy(&rootIndicesSec));
3723: PetscCall(PetscFree(parentIndices));
3724: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree, injRef, &childrenMats));
3725: PetscCall(PetscFree(rootIndices));
3726: PetscCall(PetscFree3(offsets, offsetsCopy, rowOffsets));
3728: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
3729: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
3730: PetscFunctionReturn(PETSC_SUCCESS);
3731: }
3733: static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom)
3734: {
3735: PetscSF coarseToFineEmbedded;
3736: PetscSection globalCoarse, globalFine;
3737: PetscSection localCoarse, localFine;
3738: PetscSection aSec, cSec;
3739: PetscSection rootValuesSec;
3740: PetscSection leafValuesSec;
3741: PetscScalar *rootValues, *leafValues;
3742: IS aIS;
3743: const PetscInt *anchors;
3744: Mat cMat;
3745: PetscInt numFields;
3746: PetscInt pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd;
3747: PetscInt aStart, aEnd, cStart, cEnd;
3748: PetscInt *maxChildIds;
3749: PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
3750: PetscFV fv = NULL;
3751: PetscInt dim, numFVcomps = -1, fvField = -1;
3752: DM cellDM = NULL, gradDM = NULL;
3753: const PetscScalar *cellGeomArray = NULL;
3754: const PetscScalar *gradArray = NULL;
3756: PetscFunctionBegin;
3757: PetscCall(VecSetOption(vecFine, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE));
3758: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
3759: PetscCall(DMPlexGetSimplexOrBoxCells(coarse, 0, &cellStart, &cellEnd));
3760: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
3761: PetscCall(DMGetGlobalSection(fine, &globalFine));
3762: PetscCall(DMGetCoordinateDim(coarse, &dim));
3763: { /* winnow fine points that don't have global dofs out of the sf */
3764: PetscInt nleaves, l;
3765: const PetscInt *leaves;
3766: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;
3768: PetscCall(PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL));
3770: for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
3771: PetscInt p = leaves ? leaves[l] : l;
3773: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
3774: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
3775: if ((dof - cdof) > 0) numPointsWithDofs++;
3776: }
3777: PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
3778: for (l = 0, offset = 0; l < nleaves; l++) {
3779: PetscInt p = leaves ? leaves[l] : l;
3781: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
3782: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
3783: if ((dof - cdof) > 0) pointsWithDofs[offset++] = l;
3784: }
3785: PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded));
3786: PetscCall(PetscFree(pointsWithDofs));
3787: }
3788: /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
3789: PetscCall(PetscMalloc1(pEndC - pStartC, &maxChildIds));
3790: for (p = pStartC; p < pEndC; p++) maxChildIds[p - pStartC] = -2;
3791: PetscCall(PetscSFReduceBegin(coarseToFineEmbedded, MPIU_INT, cids, maxChildIds, MPIU_MAX));
3792: PetscCall(PetscSFReduceEnd(coarseToFineEmbedded, MPIU_INT, cids, maxChildIds, MPIU_MAX));
3794: PetscCall(DMGetLocalSection(coarse, &localCoarse));
3795: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
3797: PetscCall(DMPlexGetAnchors(coarse, &aSec, &aIS));
3798: PetscCall(ISGetIndices(aIS, &anchors));
3799: PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd));
3801: PetscCall(DMGetDefaultConstraints(coarse, &cSec, &cMat, NULL));
3802: PetscCall(PetscSectionGetChart(cSec, &cStart, &cEnd));
3804: /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
3805: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootValuesSec));
3806: PetscCall(PetscSectionSetChart(rootValuesSec, pStartC, pEndC));
3807: PetscCall(PetscSectionGetNumFields(localCoarse, &numFields));
3808: {
3809: PetscInt maxFields = PetscMax(1, numFields) + 1;
3810: PetscCall(PetscMalloc7(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &newOffsets, maxFields, &newOffsetsCopy, maxFields, &rowOffsets, maxFields, &numD, maxFields, &numO));
3811: }
3812: if (grad) {
3813: PetscInt i;
3815: PetscCall(VecGetDM(cellGeom, &cellDM));
3816: PetscCall(VecGetArrayRead(cellGeom, &cellGeomArray));
3817: PetscCall(VecGetDM(grad, &gradDM));
3818: PetscCall(VecGetArrayRead(grad, &gradArray));
3819: for (i = 0; i < PetscMax(1, numFields); i++) {
3820: PetscObject obj;
3821: PetscClassId id;
3823: PetscCall(DMGetField(coarse, i, NULL, &obj));
3824: PetscCall(PetscObjectGetClassId(obj, &id));
3825: if (id == PETSCFV_CLASSID) {
3826: fv = (PetscFV)obj;
3827: PetscCall(PetscFVGetNumComponents(fv, &numFVcomps));
3828: fvField = i;
3829: break;
3830: }
3831: }
3832: }
3834: for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
3835: PetscInt dof;
3836: PetscInt maxChildId = maxChildIds[p - pStartC];
3837: PetscInt numValues = 0;
3839: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3840: if (dof < 0) dof = -(dof + 1);
3841: offsets[0] = 0;
3842: newOffsets[0] = 0;
3843: if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
3844: PetscInt *closure = NULL, closureSize, cl;
3846: PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
3847: for (cl = 0; cl < closureSize; cl++) { /* get the closure */
3848: PetscInt c = closure[2 * cl], clDof;
3850: PetscCall(PetscSectionGetDof(localCoarse, c, &clDof));
3851: numValues += clDof;
3852: }
3853: PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
3854: } else if (maxChildId == -1) {
3855: PetscCall(PetscSectionGetDof(localCoarse, p, &numValues));
3856: }
3857: /* we will pack the column indices with the field offsets */
3858: if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) {
3859: /* also send the centroid, and the gradient */
3860: numValues += dim * (1 + numFVcomps);
3861: }
3862: PetscCall(PetscSectionSetDof(rootValuesSec, p, numValues));
3863: }
3864: PetscCall(PetscSectionSetUp(rootValuesSec));
3865: {
3866: PetscInt numRootValues;
3867: const PetscScalar *coarseArray;
3869: PetscCall(PetscSectionGetStorageSize(rootValuesSec, &numRootValues));
3870: PetscCall(PetscMalloc1(numRootValues, &rootValues));
3871: PetscCall(VecGetArrayRead(vecCoarseLocal, &coarseArray));
3872: for (p = pStartC; p < pEndC; p++) {
3873: PetscInt numValues;
3874: PetscInt pValOff;
3875: PetscScalar *pVal;
3876: PetscInt maxChildId = maxChildIds[p - pStartC];
3878: PetscCall(PetscSectionGetDof(rootValuesSec, p, &numValues));
3879: if (!numValues) continue;
3880: PetscCall(PetscSectionGetOffset(rootValuesSec, p, &pValOff));
3881: pVal = &(rootValues[pValOff]);
3882: if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
3883: PetscInt closureSize = numValues;
3884: PetscCall(DMPlexVecGetClosure(coarse, NULL, vecCoarseLocal, p, &closureSize, &pVal));
3885: if (grad && p >= cellStart && p < cellEnd) {
3886: PetscFVCellGeom *cg;
3887: PetscScalar *gradVals = NULL;
3888: PetscInt i;
3890: pVal += (numValues - dim * (1 + numFVcomps));
3892: PetscCall(DMPlexPointLocalRead(cellDM, p, cellGeomArray, (void *)&cg));
3893: for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i];
3894: pVal += dim;
3895: PetscCall(DMPlexPointGlobalRead(gradDM, p, gradArray, (void *)&gradVals));
3896: for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i];
3897: }
3898: } else if (maxChildId == -1) {
3899: PetscInt lDof, lOff, i;
3901: PetscCall(PetscSectionGetDof(localCoarse, p, &lDof));
3902: PetscCall(PetscSectionGetOffset(localCoarse, p, &lOff));
3903: for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i];
3904: }
3905: }
3906: PetscCall(VecRestoreArrayRead(vecCoarseLocal, &coarseArray));
3907: PetscCall(PetscFree(maxChildIds));
3908: }
3909: {
3910: PetscSF valuesSF;
3911: PetscInt *remoteOffsetsValues, numLeafValues;
3913: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafValuesSec));
3914: PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootValuesSec, &remoteOffsetsValues, leafValuesSec));
3915: PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootValuesSec, remoteOffsetsValues, leafValuesSec, &valuesSF));
3916: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
3917: PetscCall(PetscFree(remoteOffsetsValues));
3918: PetscCall(PetscSectionGetStorageSize(leafValuesSec, &numLeafValues));
3919: PetscCall(PetscMalloc1(numLeafValues, &leafValues));
3920: PetscCall(PetscSFBcastBegin(valuesSF, MPIU_SCALAR, rootValues, leafValues, MPI_REPLACE));
3921: PetscCall(PetscSFBcastEnd(valuesSF, MPIU_SCALAR, rootValues, leafValues, MPI_REPLACE));
3922: PetscCall(PetscSFDestroy(&valuesSF));
3923: PetscCall(PetscFree(rootValues));
3924: PetscCall(PetscSectionDestroy(&rootValuesSec));
3925: }
3926: PetscCall(DMGetLocalSection(fine, &localFine));
3927: {
3928: PetscInt maxDof;
3929: PetscInt *rowIndices;
3930: DM refTree;
3931: PetscInt **refPointFieldN;
3932: PetscScalar ***refPointFieldMats;
3933: PetscSection refConSec, refAnSec;
3934: PetscInt pRefStart, pRefEnd, leafStart, leafEnd;
3935: PetscScalar *pointWork;
3937: PetscCall(PetscSectionGetMaxDof(localFine, &maxDof));
3938: PetscCall(DMGetWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
3939: PetscCall(DMGetWorkArray(fine, maxDof, MPIU_SCALAR, &pointWork));
3940: PetscCall(DMPlexGetReferenceTree(fine, &refTree));
3941: PetscCall(DMCopyDisc(fine, refTree));
3942: PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
3943: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
3944: PetscCall(DMPlexGetAnchors(refTree, &refAnSec, NULL));
3945: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
3946: PetscCall(PetscSectionGetChart(leafValuesSec, &leafStart, &leafEnd));
3947: PetscCall(DMPlexGetSimplexOrBoxCells(fine, 0, &cellStart, &cellEnd));
3948: for (p = leafStart; p < leafEnd; p++) {
3949: PetscInt gDof, gcDof, gOff, lDof;
3950: PetscInt numValues, pValOff;
3951: PetscInt childId;
3952: const PetscScalar *pVal;
3953: const PetscScalar *fvGradData = NULL;
3955: PetscCall(PetscSectionGetDof(globalFine, p, &gDof));
3956: PetscCall(PetscSectionGetDof(localFine, p, &lDof));
3957: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &gcDof));
3958: if ((gDof - gcDof) <= 0) continue;
3959: PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
3960: PetscCall(PetscSectionGetDof(leafValuesSec, p, &numValues));
3961: if (!numValues) continue;
3962: PetscCall(PetscSectionGetOffset(leafValuesSec, p, &pValOff));
3963: pVal = &leafValues[pValOff];
3964: offsets[0] = 0;
3965: offsetsCopy[0] = 0;
3966: newOffsets[0] = 0;
3967: newOffsetsCopy[0] = 0;
3968: childId = cids[p - pStartF];
3969: if (numFields) {
3970: PetscInt f;
3971: for (f = 0; f < numFields; f++) {
3972: PetscInt rowDof;
3974: PetscCall(PetscSectionGetFieldDof(localFine, p, f, &rowDof));
3975: offsets[f + 1] = offsets[f] + rowDof;
3976: offsetsCopy[f + 1] = offsets[f + 1];
3977: /* TODO: closure indices */
3978: newOffsets[f + 1] = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]);
3979: }
3980: PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices));
3981: } else {
3982: offsets[0] = 0;
3983: offsets[1] = lDof;
3984: newOffsets[0] = 0;
3985: newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0];
3986: PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices));
3987: }
3988: if (childId == -1) { /* no child interpolation: one nnz per */
3989: PetscCall(VecSetValues(vecFine, numValues, rowIndices, pVal, INSERT_VALUES));
3990: } else {
3991: PetscInt f;
3993: if (grad && p >= cellStart && p < cellEnd) {
3994: numValues -= (dim * (1 + numFVcomps));
3995: fvGradData = &pVal[numValues];
3996: }
3997: for (f = 0; f < PetscMax(1, numFields); f++) {
3998: const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f];
3999: PetscInt numRows = offsets[f + 1] - offsets[f];
4000: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
4001: const PetscScalar *cVal = &pVal[newOffsets[f]];
4002: PetscScalar *rVal = &pointWork[offsets[f]];
4003: PetscInt i, j;
4005: #if 0
4006: PetscCall(PetscInfo(coarse,"childId %" PetscInt_FMT ", numRows %" PetscInt_FMT ", numCols %" PetscInt_FMT ", refPointFieldN %" PetscInt_FMT " maxDof %" PetscInt_FMT "\n",childId,numRows,numCols,refPointFieldN[childId - pRefStart][f], maxDof));
4007: #endif
4008: for (i = 0; i < numRows; i++) {
4009: PetscScalar val = 0.;
4010: for (j = 0; j < numCols; j++) val += childMat[i * numCols + j] * cVal[j];
4011: rVal[i] = val;
4012: }
4013: if (f == fvField && p >= cellStart && p < cellEnd) {
4014: PetscReal centroid[3];
4015: PetscScalar diff[3];
4016: const PetscScalar *parentCentroid = &fvGradData[0];
4017: const PetscScalar *gradient = &fvGradData[dim];
4019: PetscCall(DMPlexComputeCellGeometryFVM(fine, p, NULL, centroid, NULL));
4020: for (i = 0; i < dim; i++) diff[i] = centroid[i] - parentCentroid[i];
4021: for (i = 0; i < numFVcomps; i++) {
4022: PetscScalar val = 0.;
4024: for (j = 0; j < dim; j++) val += gradient[dim * i + j] * diff[j];
4025: rVal[i] += val;
4026: }
4027: }
4028: PetscCall(VecSetValues(vecFine, numRows, &rowIndices[offsets[f]], rVal, INSERT_VALUES));
4029: }
4030: }
4031: }
4032: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
4033: PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_SCALAR, &pointWork));
4034: PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
4035: }
4036: PetscCall(PetscFree(leafValues));
4037: PetscCall(PetscSectionDestroy(&leafValuesSec));
4038: PetscCall(PetscFree7(offsets, offsetsCopy, newOffsets, newOffsetsCopy, rowOffsets, numD, numO));
4039: PetscCall(ISRestoreIndices(aIS, &anchors));
4040: PetscFunctionReturn(PETSC_SUCCESS);
4041: }
4043: static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids)
4044: {
4045: DM refTree;
4046: PetscSection multiRootSec, rootIndicesSec;
4047: PetscSection globalCoarse, globalFine;
4048: PetscSection localCoarse, localFine;
4049: PetscSection cSecRef;
4050: PetscInt *parentIndices, pRefStart, pRefEnd;
4051: PetscScalar *rootValues, *parentValues;
4052: Mat injRef;
4053: PetscInt numFields, maxDof;
4054: PetscInt pStartC, pEndC, pStartF, pEndF, p;
4055: PetscInt *offsets, *offsetsCopy, *rowOffsets;
4056: PetscLayout rowMap, colMap;
4057: PetscInt rowStart, rowEnd, colStart, colEnd;
4058: PetscScalar ***childrenMats = NULL; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */
4060: PetscFunctionBegin;
4062: /* get the templates for the fine-to-coarse injection from the reference tree */
4063: PetscCall(VecSetOption(vecFine, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE));
4064: PetscCall(VecSetOption(vecCoarse, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE));
4065: PetscCall(DMPlexGetReferenceTree(coarse, &refTree));
4066: PetscCall(DMCopyDisc(coarse, refTree));
4067: PetscCall(DMGetDefaultConstraints(refTree, &cSecRef, NULL, NULL));
4068: PetscCall(PetscSectionGetChart(cSecRef, &pRefStart, &pRefEnd));
4069: PetscCall(DMPlexReferenceTreeGetInjector(refTree, &injRef));
4071: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
4072: PetscCall(DMGetLocalSection(fine, &localFine));
4073: PetscCall(DMGetGlobalSection(fine, &globalFine));
4074: PetscCall(PetscSectionGetNumFields(localFine, &numFields));
4075: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
4076: PetscCall(DMGetLocalSection(coarse, &localCoarse));
4077: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
4078: PetscCall(PetscSectionGetMaxDof(localCoarse, &maxDof));
4079: {
4080: PetscInt maxFields = PetscMax(1, numFields) + 1;
4081: PetscCall(PetscMalloc3(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &rowOffsets));
4082: }
4084: PetscCall(DMPlexTransferInjectorTree(coarse, fine, coarseToFine, cids, vecFine, numFields, offsets, &multiRootSec, &rootIndicesSec, NULL, &rootValues));
4086: PetscCall(PetscMalloc2(maxDof, &parentIndices, maxDof, &parentValues));
4088: /* count indices */
4089: PetscCall(VecGetLayout(vecFine, &colMap));
4090: PetscCall(VecGetLayout(vecCoarse, &rowMap));
4091: PetscCall(PetscLayoutSetUp(rowMap));
4092: PetscCall(PetscLayoutSetUp(colMap));
4093: PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd));
4094: PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd));
4095: /* insert values */
4096: PetscCall(DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree, injRef, &childrenMats));
4097: for (p = pStartC; p < pEndC; p++) {
4098: PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
4099: PetscBool contribute = PETSC_FALSE;
4101: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
4102: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
4103: if ((dof - cdof) <= 0) continue;
4104: PetscCall(PetscSectionGetDof(localCoarse, p, &dof));
4105: PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));
4107: rowOffsets[0] = 0;
4108: offsetsCopy[0] = 0;
4109: if (numFields) {
4110: PetscInt f;
4112: for (f = 0; f < numFields; f++) {
4113: PetscInt fDof;
4114: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
4115: rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
4116: }
4117: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices));
4118: } else {
4119: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices));
4120: rowOffsets[1] = offsetsCopy[0];
4121: }
4123: PetscCall(PetscSectionGetDof(multiRootSec, p, &numLeaves));
4124: PetscCall(PetscSectionGetOffset(multiRootSec, p, &leafStart));
4125: leafEnd = leafStart + numLeaves;
4126: for (l = 0; l < dof; l++) parentValues[l] = 0.;
4127: for (l = leafStart; l < leafEnd; l++) {
4128: PetscInt numIndices, childId, offset;
4129: const PetscScalar *childValues;
4131: PetscCall(PetscSectionGetDof(rootIndicesSec, l, &numIndices));
4132: PetscCall(PetscSectionGetOffset(rootIndicesSec, l, &offset));
4133: childId = (PetscInt)PetscRealPart(rootValues[offset++]);
4134: childValues = &rootValues[offset];
4135: numIndices--;
4137: if (childId == -2) { /* skip */
4138: continue;
4139: } else if (childId == -1) { /* equivalent points: scatter */
4140: PetscInt m;
4142: contribute = PETSC_TRUE;
4143: for (m = 0; m < numIndices; m++) parentValues[m] = childValues[m];
4144: } else { /* contributions from children: sum with injectors from reference tree */
4145: PetscInt parentId, f, lim;
4147: contribute = PETSC_TRUE;
4148: PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL));
4150: lim = PetscMax(1, numFields);
4151: offsets[0] = 0;
4152: if (numFields) {
4153: PetscInt f;
4155: for (f = 0; f < numFields; f++) {
4156: PetscInt fDof;
4157: PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof));
4159: offsets[f + 1] = fDof + offsets[f];
4160: }
4161: } else {
4162: PetscInt cDof;
4164: PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof));
4165: offsets[1] = cDof;
4166: }
4167: for (f = 0; f < lim; f++) {
4168: PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0];
4169: PetscInt n = offsets[f + 1] - offsets[f];
4170: PetscInt m = rowOffsets[f + 1] - rowOffsets[f];
4171: PetscInt i, j;
4172: const PetscScalar *colValues = &childValues[offsets[f]];
4174: for (i = 0; i < m; i++) {
4175: PetscScalar val = 0.;
4176: for (j = 0; j < n; j++) val += childMat[n * i + j] * colValues[j];
4177: parentValues[rowOffsets[f] + i] += val;
4178: }
4179: }
4180: }
4181: }
4182: if (contribute) PetscCall(VecSetValues(vecCoarse, dof, parentIndices, parentValues, INSERT_VALUES));
4183: }
4184: PetscCall(PetscSectionDestroy(&multiRootSec));
4185: PetscCall(PetscSectionDestroy(&rootIndicesSec));
4186: PetscCall(PetscFree2(parentIndices, parentValues));
4187: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree, injRef, &childrenMats));
4188: PetscCall(PetscFree(rootValues));
4189: PetscCall(PetscFree3(offsets, offsetsCopy, rowOffsets));
4190: PetscFunctionReturn(PETSC_SUCCESS);
4191: }
4193: /*@
4194: DMPlexTransferVecTree - transfer a vector between two meshes that differ from each other by refinement/coarsening
4195: that can be represented by a common reference tree used by both. This routine can be used for a combination of
4196: coarsening and refinement at the same time.
4198: Collective
4200: Input Parameters:
4201: + dmIn - The `DMPLEX` mesh for the input vector
4202: . dmOut - The second `DMPLEX` mesh
4203: . vecIn - The input vector
4204: . sfRefine - A star forest indicating points in the mesh `dmIn` (roots in the star forest) that are parents to points in
4205: the mesh `dmOut` (leaves in the star forest), i.e. where `dmOut` is more refined than `dmIn`
4206: . sfCoarsen - A star forest indicating points in the mesh `dmOut` (roots in the star forest) that are parents to points in
4207: the mesh `dmIn` (leaves in the star forest), i.e. where `dmOut` is more coarsened than `dmIn`
4208: . cidsRefine - The childIds of the points in `dmOut`. These childIds relate back to the reference tree: childid[j] = k implies
4209: that mesh point j of `dmOut` was refined from a point in `dmIn` just as the mesh point k in the reference
4210: tree was refined from its parent. childid[j] = -1 indicates that the point j in `dmOut` is exactly
4211: equivalent to its root in `dmIn`, so no interpolation is necessary. childid[j] = -2 indicates that this
4212: point j in `dmOut` is not a leaf of `sfRefine`.
4213: . cidsCoarsen - The childIds of the points in `dmIn`. These childIds relate back to the reference tree: childid[j] = k implies
4214: that mesh point j of dmIn coarsens to a point in `dmOut` just as the mesh point k in the reference
4215: tree coarsens to its parent. childid[j] = -2 indicates that point j in `dmOut` is not a leaf in `sfCoarsen`.
4216: . useBCs - `PETSC_TRUE` indicates that boundary values should be inserted into `vecIn` before transfer.
4217: - time - Used if boundary values are time dependent.
4219: Output Parameter:
4220: . vecOut - Using interpolation and injection operators calculated on the reference tree, the transferred
4221: projection of `vecIn` from `dmIn` to `dmOut`. Note that any field discretized with a `PetscFV` finite volume
4222: method that uses gradient reconstruction will use reconstructed gradients when interpolating from
4223: coarse points to fine points.
4225: Level: developer
4227: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `Vec`, `PetscFV`, `DMPlexSetReferenceTree()`, `DMPlexGetReferenceTree()`, `PetscFVGetComputeGradients()`
4228: @*/
4229: PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time)
4230: {
4231: PetscFunctionBegin;
4232: PetscCall(VecSet(vecOut, 0.0));
4233: if (sfRefine) {
4234: Vec vecInLocal;
4235: DM dmGrad = NULL;
4236: Vec faceGeom = NULL, cellGeom = NULL, grad = NULL;
4238: PetscCall(DMGetLocalVector(dmIn, &vecInLocal));
4239: PetscCall(VecSet(vecInLocal, 0.0));
4240: {
4241: PetscInt numFields, i;
4243: PetscCall(DMGetNumFields(dmIn, &numFields));
4244: for (i = 0; i < numFields; i++) {
4245: PetscObject obj;
4246: PetscClassId classid;
4248: PetscCall(DMGetField(dmIn, i, NULL, &obj));
4249: PetscCall(PetscObjectGetClassId(obj, &classid));
4250: if (classid == PETSCFV_CLASSID) {
4251: PetscCall(DMPlexGetDataFVM(dmIn, (PetscFV)obj, &cellGeom, &faceGeom, &dmGrad));
4252: break;
4253: }
4254: }
4255: }
4256: if (useBCs) PetscCall(DMPlexInsertBoundaryValues(dmIn, PETSC_TRUE, vecInLocal, time, faceGeom, cellGeom, NULL));
4257: PetscCall(DMGlobalToLocalBegin(dmIn, vecIn, INSERT_VALUES, vecInLocal));
4258: PetscCall(DMGlobalToLocalEnd(dmIn, vecIn, INSERT_VALUES, vecInLocal));
4259: if (dmGrad) {
4260: PetscCall(DMGetGlobalVector(dmGrad, &grad));
4261: PetscCall(DMPlexReconstructGradientsFVM(dmIn, vecInLocal, grad));
4262: }
4263: PetscCall(DMPlexTransferVecTree_Interpolate(dmIn, vecInLocal, dmOut, vecOut, sfRefine, cidsRefine, grad, cellGeom));
4264: PetscCall(DMRestoreLocalVector(dmIn, &vecInLocal));
4265: if (dmGrad) PetscCall(DMRestoreGlobalVector(dmGrad, &grad));
4266: }
4267: if (sfCoarsen) PetscCall(DMPlexTransferVecTree_Inject(dmIn, vecIn, dmOut, vecOut, sfCoarsen, cidsCoarsen));
4268: PetscCall(VecAssemblyBegin(vecOut));
4269: PetscCall(VecAssemblyEnd(vecOut));
4270: PetscFunctionReturn(PETSC_SUCCESS);
4271: }