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