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;
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(MPI_IN_PLACE, &anyNew, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)secNew)));
614: if (!anyNew) {
615: PetscCall(PetscSectionDestroy(&secNew));
616: *sectionNew = NULL;
617: *isNew = NULL;
618: } else {
619: if (compress) {
620: PetscSection secComp;
621: PetscInt *valsComp = NULL;
623: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), &secComp));
624: PetscCall(PetscSectionSetChart(secComp, pStart, pEnd));
625: for (p = pStart; p < pEnd; p++) {
626: PetscInt dof;
628: PetscCall(PetscSectionGetDof(secNew, p, &dof));
629: PetscCall(PetscSectionSetDof(secComp, p, dof));
630: }
631: PetscCall(PetscSectionSetUp(secComp));
632: PetscCall(PetscSectionGetStorageSize(secComp, &sizeNew));
633: PetscCall(PetscMalloc1(sizeNew, &valsComp));
634: for (p = pStart; p < pEnd; p++) {
635: PetscInt dof, off, offNew, j;
637: PetscCall(PetscSectionGetDof(secNew, p, &dof));
638: PetscCall(PetscSectionGetOffset(secNew, p, &off));
639: PetscCall(PetscSectionGetOffset(secComp, p, &offNew));
640: for (j = 0; j < dof; j++) valsComp[offNew + j] = valsNew[off + j];
641: }
642: PetscCall(PetscSectionDestroy(&secNew));
643: secNew = secComp;
644: PetscCall(PetscFree(valsNew));
645: valsNew = valsComp;
646: }
647: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), sizeNew, valsNew, PETSC_OWN_POINTER, isNew));
648: }
649: PetscFunctionReturn(PETSC_SUCCESS);
650: }
652: static PetscErrorCode DMPlexCreateAnchors_Tree(DM dm)
653: {
654: PetscInt p, pStart, pEnd, *anchors, size;
655: PetscInt aMin = PETSC_INT_MAX, aMax = PETSC_INT_MIN;
656: PetscSection aSec;
657: DMLabel canonLabel;
658: IS aIS;
660: PetscFunctionBegin;
662: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
663: PetscCall(DMGetLabel(dm, "canonical", &canonLabel));
664: for (p = pStart; p < pEnd; p++) {
665: PetscInt parent;
667: if (canonLabel) {
668: PetscInt canon;
670: PetscCall(DMLabelGetValue(canonLabel, p, &canon));
671: if (p != canon) continue;
672: }
673: PetscCall(DMPlexGetTreeParent(dm, p, &parent, NULL));
674: if (parent != p) {
675: aMin = PetscMin(aMin, p);
676: aMax = PetscMax(aMax, p + 1);
677: }
678: }
679: if (aMin > aMax) {
680: aMin = -1;
681: aMax = -1;
682: }
683: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &aSec));
684: PetscCall(PetscSectionSetChart(aSec, aMin, aMax));
685: for (p = aMin; p < aMax; p++) {
686: PetscInt parent, ancestor = p;
688: if (canonLabel) {
689: PetscInt canon;
691: PetscCall(DMLabelGetValue(canonLabel, p, &canon));
692: if (p != canon) continue;
693: }
694: PetscCall(DMPlexGetTreeParent(dm, p, &parent, NULL));
695: while (parent != ancestor) {
696: ancestor = parent;
697: PetscCall(DMPlexGetTreeParent(dm, ancestor, &parent, NULL));
698: }
699: if (ancestor != p) {
700: PetscInt closureSize, *closure = NULL;
702: PetscCall(DMPlexGetTransitiveClosure(dm, ancestor, PETSC_TRUE, &closureSize, &closure));
703: PetscCall(PetscSectionSetDof(aSec, p, closureSize));
704: PetscCall(DMPlexRestoreTransitiveClosure(dm, ancestor, PETSC_TRUE, &closureSize, &closure));
705: }
706: }
707: PetscCall(PetscSectionSetUp(aSec));
708: PetscCall(PetscSectionGetStorageSize(aSec, &size));
709: PetscCall(PetscMalloc1(size, &anchors));
710: for (p = aMin; p < aMax; p++) {
711: PetscInt parent, ancestor = p;
713: if (canonLabel) {
714: PetscInt canon;
716: PetscCall(DMLabelGetValue(canonLabel, p, &canon));
717: if (p != canon) continue;
718: }
719: PetscCall(DMPlexGetTreeParent(dm, p, &parent, NULL));
720: while (parent != ancestor) {
721: ancestor = parent;
722: PetscCall(DMPlexGetTreeParent(dm, ancestor, &parent, NULL));
723: }
724: if (ancestor != p) {
725: PetscInt j, closureSize, *closure = NULL, aOff;
727: PetscCall(PetscSectionGetOffset(aSec, p, &aOff));
729: PetscCall(DMPlexGetTransitiveClosure(dm, ancestor, PETSC_TRUE, &closureSize, &closure));
730: for (j = 0; j < closureSize; j++) anchors[aOff + j] = closure[2 * j];
731: PetscCall(DMPlexRestoreTransitiveClosure(dm, ancestor, PETSC_TRUE, &closureSize, &closure));
732: }
733: }
734: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size, anchors, PETSC_OWN_POINTER, &aIS));
735: {
736: PetscSection aSecNew = aSec;
737: IS aISNew = aIS;
739: PetscCall(PetscObjectReference((PetscObject)aSec));
740: PetscCall(PetscObjectReference((PetscObject)aIS));
741: while (aSecNew) {
742: PetscCall(PetscSectionDestroy(&aSec));
743: PetscCall(ISDestroy(&aIS));
744: aSec = aSecNew;
745: aIS = aISNew;
746: aSecNew = NULL;
747: aISNew = NULL;
748: PetscCall(AnchorsFlatten(aSec, aIS, &aSecNew, &aISNew));
749: }
750: }
751: PetscCall(DMPlexSetAnchors(dm, aSec, aIS));
752: PetscCall(PetscSectionDestroy(&aSec));
753: PetscCall(ISDestroy(&aIS));
754: PetscFunctionReturn(PETSC_SUCCESS);
755: }
757: static PetscErrorCode DMPlexGetTrueSupportSize(DM dm, PetscInt p, PetscInt *dof, PetscInt *numTrueSupp)
758: {
759: PetscFunctionBegin;
760: if (numTrueSupp[p] == -1) {
761: PetscInt i, alldof;
762: const PetscInt *supp;
763: PetscInt count = 0;
765: PetscCall(DMPlexGetSupportSize(dm, p, &alldof));
766: PetscCall(DMPlexGetSupport(dm, p, &supp));
767: for (i = 0; i < alldof; i++) {
768: PetscInt q = supp[i], numCones, j;
769: const PetscInt *cone;
771: PetscCall(DMPlexGetConeSize(dm, q, &numCones));
772: PetscCall(DMPlexGetCone(dm, q, &cone));
773: for (j = 0; j < numCones; j++) {
774: if (cone[j] == p) break;
775: }
776: if (j < numCones) count++;
777: }
778: numTrueSupp[p] = count;
779: }
780: *dof = numTrueSupp[p];
781: PetscFunctionReturn(PETSC_SUCCESS);
782: }
784: static PetscErrorCode DMPlexTreeExchangeSupports(DM dm)
785: {
786: DM_Plex *mesh = (DM_Plex *)dm->data;
787: PetscSection newSupportSection;
788: PetscInt newSize, *newSupports, pStart, pEnd, p, d, depth;
789: PetscInt *numTrueSupp;
790: PetscInt *offsets;
792: PetscFunctionBegin;
794: /* symmetrize the hierarchy */
795: PetscCall(DMPlexGetDepth(dm, &depth));
796: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)mesh->supportSection), &newSupportSection));
797: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
798: PetscCall(PetscSectionSetChart(newSupportSection, pStart, pEnd));
799: PetscCall(PetscCalloc1(pEnd, &offsets));
800: PetscCall(PetscMalloc1(pEnd, &numTrueSupp));
801: for (p = 0; p < pEnd; p++) numTrueSupp[p] = -1;
802: /* if a point is in the (true) support of q, it should be in the support of
803: * parent(q) */
804: for (d = 0; d <= depth; d++) {
805: PetscCall(DMPlexGetHeightStratum(dm, d, &pStart, &pEnd));
806: for (p = pStart; p < pEnd; ++p) {
807: PetscInt dof, q, qdof, parent;
809: PetscCall(DMPlexGetTrueSupportSize(dm, p, &dof, numTrueSupp));
810: PetscCall(PetscSectionAddDof(newSupportSection, p, dof));
811: q = p;
812: PetscCall(DMPlexGetTreeParent(dm, q, &parent, NULL));
813: while (parent != q && parent >= pStart && parent < pEnd) {
814: q = parent;
816: PetscCall(DMPlexGetTrueSupportSize(dm, q, &qdof, numTrueSupp));
817: PetscCall(PetscSectionAddDof(newSupportSection, p, qdof));
818: PetscCall(PetscSectionAddDof(newSupportSection, q, dof));
819: PetscCall(DMPlexGetTreeParent(dm, q, &parent, NULL));
820: }
821: }
822: }
823: PetscCall(PetscSectionSetUp(newSupportSection));
824: PetscCall(PetscSectionGetStorageSize(newSupportSection, &newSize));
825: PetscCall(PetscMalloc1(newSize, &newSupports));
826: for (d = 0; d <= depth; d++) {
827: PetscCall(DMPlexGetHeightStratum(dm, d, &pStart, &pEnd));
828: for (p = pStart; p < pEnd; p++) {
829: PetscInt dof, off, q, qdof, qoff, newDof, newOff, newqOff, i, parent;
831: PetscCall(PetscSectionGetDof(mesh->supportSection, p, &dof));
832: PetscCall(PetscSectionGetOffset(mesh->supportSection, p, &off));
833: PetscCall(PetscSectionGetDof(newSupportSection, p, &newDof));
834: PetscCall(PetscSectionGetOffset(newSupportSection, p, &newOff));
835: for (i = 0; i < dof; i++) {
836: PetscInt numCones, j;
837: const PetscInt *cone;
838: PetscInt q = mesh->supports[off + i];
840: PetscCall(DMPlexGetConeSize(dm, q, &numCones));
841: PetscCall(DMPlexGetCone(dm, q, &cone));
842: for (j = 0; j < numCones; j++) {
843: if (cone[j] == p) break;
844: }
845: if (j < numCones) newSupports[newOff + offsets[p]++] = q;
846: }
848: q = p;
849: PetscCall(DMPlexGetTreeParent(dm, q, &parent, NULL));
850: while (parent != q && parent >= pStart && parent < pEnd) {
851: q = parent;
852: PetscCall(PetscSectionGetDof(mesh->supportSection, q, &qdof));
853: PetscCall(PetscSectionGetOffset(mesh->supportSection, q, &qoff));
854: PetscCall(PetscSectionGetOffset(newSupportSection, q, &newqOff));
855: for (i = 0; i < qdof; i++) {
856: PetscInt numCones, j;
857: const PetscInt *cone;
858: PetscInt r = mesh->supports[qoff + i];
860: PetscCall(DMPlexGetConeSize(dm, r, &numCones));
861: PetscCall(DMPlexGetCone(dm, r, &cone));
862: for (j = 0; j < numCones; j++) {
863: if (cone[j] == q) break;
864: }
865: if (j < numCones) newSupports[newOff + offsets[p]++] = r;
866: }
867: for (i = 0; i < dof; i++) {
868: PetscInt numCones, j;
869: const PetscInt *cone;
870: PetscInt r = mesh->supports[off + i];
872: PetscCall(DMPlexGetConeSize(dm, r, &numCones));
873: PetscCall(DMPlexGetCone(dm, r, &cone));
874: for (j = 0; j < numCones; j++) {
875: if (cone[j] == p) break;
876: }
877: if (j < numCones) newSupports[newqOff + offsets[q]++] = r;
878: }
879: PetscCall(DMPlexGetTreeParent(dm, q, &parent, NULL));
880: }
881: }
882: }
883: PetscCall(PetscSectionDestroy(&mesh->supportSection));
884: mesh->supportSection = newSupportSection;
885: PetscCall(PetscFree(mesh->supports));
886: mesh->supports = newSupports;
887: PetscCall(PetscFree(offsets));
888: PetscCall(PetscFree(numTrueSupp));
889: PetscFunctionReturn(PETSC_SUCCESS);
890: }
892: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM, PetscSection, PetscSection, Mat);
893: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM, PetscSection, PetscSection, Mat);
895: static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical, PetscBool exchangeSupports)
896: {
897: DM_Plex *mesh = (DM_Plex *)dm->data;
898: DM refTree;
899: PetscInt size;
901: PetscFunctionBegin;
904: PetscCall(PetscObjectReference((PetscObject)parentSection));
905: PetscCall(PetscSectionDestroy(&mesh->parentSection));
906: mesh->parentSection = parentSection;
907: PetscCall(PetscSectionGetStorageSize(parentSection, &size));
908: if (parents != mesh->parents) {
909: PetscCall(PetscFree(mesh->parents));
910: PetscCall(PetscMalloc1(size, &mesh->parents));
911: PetscCall(PetscArraycpy(mesh->parents, parents, size));
912: }
913: if (childIDs != mesh->childIDs) {
914: PetscCall(PetscFree(mesh->childIDs));
915: PetscCall(PetscMalloc1(size, &mesh->childIDs));
916: PetscCall(PetscArraycpy(mesh->childIDs, childIDs, size));
917: }
918: PetscCall(DMPlexGetReferenceTree(dm, &refTree));
919: if (refTree) {
920: DMLabel canonLabel;
922: PetscCall(DMGetLabel(refTree, "canonical", &canonLabel));
923: if (canonLabel) {
924: PetscInt i;
926: for (i = 0; i < size; i++) {
927: PetscInt canon;
928: PetscCall(DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon));
929: if (canon >= 0) mesh->childIDs[i] = canon;
930: }
931: }
932: mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_FromReference;
933: } else {
934: mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_Direct;
935: }
936: PetscCall(DMPlexTreeSymmetrize(dm));
937: if (computeCanonical) {
938: PetscInt dim;
940: /* add the canonical label */
941: PetscCall(DMGetDimension(dm, &dim));
942: PetscCall(DMCreateLabel(dm, "canonical"));
943: for (PetscInt d = 0; d <= dim; d++) {
944: PetscInt p, dStart, dEnd, canon = -1, cNumChildren;
945: const PetscInt *cChildren;
947: PetscCall(DMPlexGetDepthStratum(dm, d, &dStart, &dEnd));
948: for (p = dStart; p < dEnd; p++) {
949: PetscCall(DMPlexGetTreeChildren(dm, p, &cNumChildren, &cChildren));
950: if (cNumChildren) {
951: canon = p;
952: break;
953: }
954: }
955: if (canon == -1) continue;
956: for (p = dStart; p < dEnd; p++) {
957: PetscInt numChildren;
958: const PetscInt *children;
960: PetscCall(DMPlexGetTreeChildren(dm, p, &numChildren, &children));
961: if (numChildren) {
962: 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);
963: PetscCall(DMSetLabelValue(dm, "canonical", p, canon));
964: for (PetscInt i = 0; i < numChildren; i++) PetscCall(DMSetLabelValue(dm, "canonical", children[i], cChildren[i]));
965: }
966: }
967: }
968: }
969: if (exchangeSupports) PetscCall(DMPlexTreeExchangeSupports(dm));
970: mesh->createanchors = DMPlexCreateAnchors_Tree;
971: /* reset anchors */
972: PetscCall(DMPlexSetAnchors(dm, NULL, NULL));
973: PetscFunctionReturn(PETSC_SUCCESS);
974: }
976: /*@
977: DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points. This routine also creates
978: the point-to-point constraints determined by the tree: a point is constrained to the points in the closure of its
979: tree root.
981: Collective
983: Input Parameters:
984: + dm - the `DMPLEX` object
985: . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
986: offset indexes the parent and childID list; the reference count of parentSection is incremented
987: . parents - a list of the point parents; copied, can be destroyed
988: - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
989: the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed
991: Level: intermediate
993: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetTree()`, `DMPlexSetReferenceTree()`, `DMPlexSetAnchors()`, `DMPlexGetTreeParent()`, `DMPlexGetTreeChildren()`
994: @*/
995: PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
996: {
997: PetscFunctionBegin;
998: PetscCall(DMPlexSetTree_Internal(dm, parentSection, parents, childIDs, PETSC_FALSE, PETSC_TRUE));
999: PetscFunctionReturn(PETSC_SUCCESS);
1000: }
1002: /*@
1003: DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points.
1004: Collective
1006: Input Parameter:
1007: . dm - the `DMPLEX` object
1009: Output Parameters:
1010: + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
1011: offset indexes the parent and childID list
1012: . parents - a list of the point parents
1013: . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
1014: the child corresponds to the point in the reference tree with index childID
1015: . childSection - the inverse of the parent section
1016: - children - a list of the point children
1018: Level: intermediate
1020: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetTree()`, `DMPlexSetReferenceTree()`, `DMPlexSetAnchors()`, `DMPlexGetTreeParent()`, `DMPlexGetTreeChildren()`
1021: @*/
1022: PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[])
1023: {
1024: DM_Plex *mesh = (DM_Plex *)dm->data;
1026: PetscFunctionBegin;
1028: if (parentSection) *parentSection = mesh->parentSection;
1029: if (parents) *parents = mesh->parents;
1030: if (childIDs) *childIDs = mesh->childIDs;
1031: if (childSection) *childSection = mesh->childSection;
1032: if (children) *children = mesh->children;
1033: PetscFunctionReturn(PETSC_SUCCESS);
1034: }
1036: /*@
1037: DMPlexGetTreeParent - get the parent of a point in the tree describing the point hierarchy (not the DAG)
1039: Input Parameters:
1040: + dm - the `DMPLEX` object
1041: - point - the query point
1043: Output Parameters:
1044: + parent - if not `NULL`, set to the parent of the point, or the point itself if the point does not have a parent
1045: - childID - if not `NULL`, set to the child ID of the point with respect to its parent, or 0 if the point
1046: does not have a parent
1048: Level: intermediate
1050: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetTree()`, `DMPlexGetTree()`, `DMPlexGetTreeChildren()`
1051: @*/
1052: PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID)
1053: {
1054: DM_Plex *mesh = (DM_Plex *)dm->data;
1055: PetscSection pSec;
1057: PetscFunctionBegin;
1059: pSec = mesh->parentSection;
1060: if (pSec && point >= pSec->pStart && point < pSec->pEnd) {
1061: PetscInt dof;
1063: PetscCall(PetscSectionGetDof(pSec, point, &dof));
1064: if (dof) {
1065: PetscInt off;
1067: PetscCall(PetscSectionGetOffset(pSec, point, &off));
1068: if (parent) *parent = mesh->parents[off];
1069: if (childID) *childID = mesh->childIDs[off];
1070: PetscFunctionReturn(PETSC_SUCCESS);
1071: }
1072: }
1073: if (parent) *parent = point;
1074: if (childID) *childID = 0;
1075: PetscFunctionReturn(PETSC_SUCCESS);
1076: }
1078: /*@C
1079: DMPlexGetTreeChildren - get the children of a point in the tree describing the point hierarchy (not the DAG)
1081: Input Parameters:
1082: + dm - the `DMPLEX` object
1083: - point - the query point
1085: Output Parameters:
1086: + numChildren - if not `NULL`, set to the number of children
1087: - children - if not `NULL`, set to a list children, or set to `NULL` if the point has no children
1089: Level: intermediate
1091: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetTree()`, `DMPlexGetTree()`, `DMPlexGetTreeParent()`
1092: @*/
1093: PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[])
1094: {
1095: DM_Plex *mesh = (DM_Plex *)dm->data;
1096: PetscSection childSec;
1097: PetscInt dof = 0;
1099: PetscFunctionBegin;
1101: childSec = mesh->childSection;
1102: if (childSec && point >= childSec->pStart && point < childSec->pEnd) PetscCall(PetscSectionGetDof(childSec, point, &dof));
1103: if (numChildren) *numChildren = dof;
1104: if (children) {
1105: if (dof) {
1106: PetscInt off;
1108: PetscCall(PetscSectionGetOffset(childSec, point, &off));
1109: *children = &mesh->children[off];
1110: } else {
1111: *children = NULL;
1112: }
1113: }
1114: PetscFunctionReturn(PETSC_SUCCESS);
1115: }
1117: 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)
1118: {
1119: PetscInt f, b, p, c, offset, qPoints;
1121: PetscFunctionBegin;
1122: PetscCall(PetscSpaceEvaluate(space, nPoints, points, work, NULL, NULL));
1123: for (f = 0, offset = 0; f < nFunctionals; f++) {
1124: qPoints = pointsPerFn[f];
1125: for (b = 0; b < nBasis; b++) {
1126: PetscScalar val = 0.;
1128: for (p = 0; p < qPoints; p++) {
1129: for (c = 0; c < nComps; c++) val += work[((offset + p) * nBasis + b) * nComps + c] * weights[(offset + p) * nComps + c];
1130: }
1131: PetscCall(MatSetValue(basisAtPoints, b, f, val, INSERT_VALUES));
1132: }
1133: offset += qPoints;
1134: }
1135: PetscCall(MatAssemblyBegin(basisAtPoints, MAT_FINAL_ASSEMBLY));
1136: PetscCall(MatAssemblyEnd(basisAtPoints, MAT_FINAL_ASSEMBLY));
1137: PetscFunctionReturn(PETSC_SUCCESS);
1138: }
1140: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM dm, PetscSection section, PetscSection cSec, Mat cMat)
1141: {
1142: PetscDS ds;
1143: PetscInt spdim;
1144: PetscInt numFields, f, c, cStart, cEnd, pStart, pEnd, conStart, conEnd;
1145: const PetscInt *anchors;
1146: PetscSection aSec;
1147: PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJparent, detJ, detJparent;
1148: IS aIS;
1150: PetscFunctionBegin;
1151: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1152: PetscCall(DMGetDS(dm, &ds));
1153: PetscCall(PetscDSGetNumFields(ds, &numFields));
1154: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1155: PetscCall(DMPlexGetAnchors(dm, &aSec, &aIS));
1156: PetscCall(ISGetIndices(aIS, &anchors));
1157: PetscCall(PetscSectionGetChart(cSec, &conStart, &conEnd));
1158: PetscCall(DMGetDimension(dm, &spdim));
1159: PetscCall(PetscMalloc6(spdim, &v0, spdim, &v0parent, spdim, &vtmp, spdim * spdim, &J, spdim * spdim, &Jparent, spdim * spdim, &invJparent));
1161: for (f = 0; f < numFields; f++) {
1162: PetscObject disc;
1163: PetscClassId id;
1164: PetscSpace bspace;
1165: PetscDualSpace dspace;
1166: PetscInt i, j, k, nPoints, Nc, offset;
1167: PetscInt fSize, maxDof;
1168: PetscReal *weights, *pointsRef, *pointsReal, *work;
1169: PetscScalar *scwork;
1170: const PetscScalar *X;
1171: PetscInt *sizes, *workIndRow, *workIndCol;
1172: Mat Amat, Bmat, Xmat;
1173: const PetscInt *numDof = NULL;
1174: const PetscInt ***perms = NULL;
1175: const PetscScalar ***flips = NULL;
1177: PetscCall(PetscDSGetDiscretization(ds, f, &disc));
1178: PetscCall(PetscObjectGetClassId(disc, &id));
1179: if (id == PETSCFE_CLASSID) {
1180: PetscFE fe = (PetscFE)disc;
1182: PetscCall(PetscFEGetBasisSpace(fe, &bspace));
1183: PetscCall(PetscFEGetDualSpace(fe, &dspace));
1184: PetscCall(PetscDualSpaceGetDimension(dspace, &fSize));
1185: PetscCall(PetscFEGetNumComponents(fe, &Nc));
1186: } else if (id == PETSCFV_CLASSID) {
1187: PetscFV fv = (PetscFV)disc;
1189: PetscCall(PetscFVGetNumComponents(fv, &Nc));
1190: PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)fv), &bspace));
1191: PetscCall(PetscSpaceSetType(bspace, PETSCSPACEPOLYNOMIAL));
1192: PetscCall(PetscSpaceSetDegree(bspace, 0, PETSC_DETERMINE));
1193: PetscCall(PetscSpaceSetNumComponents(bspace, Nc));
1194: PetscCall(PetscSpaceSetNumVariables(bspace, spdim));
1195: PetscCall(PetscSpaceSetUp(bspace));
1196: PetscCall(PetscFVGetDualSpace(fv, &dspace));
1197: PetscCall(PetscDualSpaceGetDimension(dspace, &fSize));
1198: } else SETERRQ(PetscObjectComm(disc), PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id);
1199: PetscCall(PetscDualSpaceGetNumDof(dspace, &numDof));
1200: for (i = 0, maxDof = 0; i <= spdim; i++) maxDof = PetscMax(maxDof, numDof[i]);
1201: PetscCall(PetscDualSpaceGetSymmetries(dspace, &perms, &flips));
1203: PetscCall(MatCreate(PETSC_COMM_SELF, &Amat));
1204: PetscCall(MatSetSizes(Amat, fSize, fSize, fSize, fSize));
1205: PetscCall(MatSetType(Amat, MATSEQDENSE));
1206: PetscCall(MatSetUp(Amat));
1207: PetscCall(MatDuplicate(Amat, MAT_DO_NOT_COPY_VALUES, &Bmat));
1208: PetscCall(MatDuplicate(Amat, MAT_DO_NOT_COPY_VALUES, &Xmat));
1209: nPoints = 0;
1210: for (i = 0; i < fSize; i++) {
1211: PetscInt qPoints, thisNc;
1212: PetscQuadrature quad;
1214: PetscCall(PetscDualSpaceGetFunctional(dspace, i, &quad));
1215: PetscCall(PetscQuadratureGetData(quad, NULL, &thisNc, &qPoints, NULL, NULL));
1216: PetscCheck(thisNc == Nc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Functional dim %" PetscInt_FMT " does not much basis dim %" PetscInt_FMT, thisNc, Nc);
1217: nPoints += qPoints;
1218: }
1219: PetscCall(PetscMalloc7(fSize, &sizes, nPoints * Nc, &weights, spdim * nPoints, &pointsRef, spdim * nPoints, &pointsReal, nPoints * fSize * Nc, &work, maxDof, &workIndRow, maxDof, &workIndCol));
1220: PetscCall(PetscMalloc1(maxDof * maxDof, &scwork));
1221: offset = 0;
1222: for (i = 0; i < fSize; i++) {
1223: PetscInt qPoints;
1224: const PetscReal *p, *w;
1225: PetscQuadrature quad;
1227: PetscCall(PetscDualSpaceGetFunctional(dspace, i, &quad));
1228: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &qPoints, &p, &w));
1229: PetscCall(PetscArraycpy(weights + Nc * offset, w, Nc * qPoints));
1230: PetscCall(PetscArraycpy(pointsRef + spdim * offset, p, spdim * qPoints));
1231: sizes[i] = qPoints;
1232: offset += qPoints;
1233: }
1234: PetscCall(EvaluateBasis(bspace, fSize, fSize, Nc, nPoints, sizes, pointsRef, weights, work, Amat));
1235: PetscCall(MatLUFactor(Amat, NULL, NULL, NULL));
1236: for (c = cStart; c < cEnd; c++) {
1237: PetscInt parent;
1238: PetscInt closureSize, closureSizeP, *closure = NULL, *closureP = NULL;
1239: PetscInt *childOffsets, *parentOffsets;
1241: PetscCall(DMPlexGetTreeParent(dm, c, &parent, NULL));
1242: if (parent == c) continue;
1243: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1244: for (i = 0; i < closureSize; i++) {
1245: PetscInt p = closure[2 * i];
1246: PetscInt conDof;
1248: if (p < conStart || p >= conEnd) continue;
1249: if (numFields) PetscCall(PetscSectionGetFieldDof(cSec, p, f, &conDof));
1250: else PetscCall(PetscSectionGetDof(cSec, p, &conDof));
1251: if (conDof) break;
1252: }
1253: if (i == closureSize) {
1254: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1255: continue;
1256: }
1258: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, NULL, &detJ));
1259: PetscCall(DMPlexComputeCellGeometryFEM(dm, parent, NULL, v0parent, Jparent, invJparent, &detJparent));
1260: for (i = 0; i < nPoints; i++) {
1261: const PetscReal xi0[3] = {-1., -1., -1.};
1263: CoordinatesRefToReal(spdim, spdim, xi0, v0, J, &pointsRef[i * spdim], vtmp);
1264: CoordinatesRealToRef(spdim, spdim, xi0, v0parent, invJparent, vtmp, &pointsReal[i * spdim]);
1265: }
1266: PetscCall(EvaluateBasis(bspace, fSize, fSize, Nc, nPoints, sizes, pointsReal, weights, work, Bmat));
1267: PetscCall(MatMatSolve(Amat, Bmat, Xmat));
1268: PetscCall(MatDenseGetArrayRead(Xmat, &X));
1269: PetscCall(DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSizeP, &closureP));
1270: PetscCall(PetscMalloc2(closureSize + 1, &childOffsets, closureSizeP + 1, &parentOffsets));
1271: childOffsets[0] = 0;
1272: for (i = 0; i < closureSize; i++) {
1273: PetscInt p = closure[2 * i];
1274: PetscInt dof;
1276: if (numFields) PetscCall(PetscSectionGetFieldDof(section, p, f, &dof));
1277: else PetscCall(PetscSectionGetDof(section, p, &dof));
1278: childOffsets[i + 1] = childOffsets[i] + dof;
1279: }
1280: parentOffsets[0] = 0;
1281: for (i = 0; i < closureSizeP; i++) {
1282: PetscInt p = closureP[2 * i];
1283: PetscInt dof;
1285: if (numFields) PetscCall(PetscSectionGetFieldDof(section, p, f, &dof));
1286: else PetscCall(PetscSectionGetDof(section, p, &dof));
1287: parentOffsets[i + 1] = parentOffsets[i] + dof;
1288: }
1289: for (i = 0; i < closureSize; i++) {
1290: PetscInt conDof, conOff, aDof, aOff, nWork;
1291: PetscInt p = closure[2 * i];
1292: PetscInt o = closure[2 * i + 1];
1293: const PetscInt *perm;
1294: const PetscScalar *flip;
1296: if (p < conStart || p >= conEnd) continue;
1297: if (numFields) {
1298: PetscCall(PetscSectionGetFieldDof(cSec, p, f, &conDof));
1299: PetscCall(PetscSectionGetFieldOffset(cSec, p, f, &conOff));
1300: } else {
1301: PetscCall(PetscSectionGetDof(cSec, p, &conDof));
1302: PetscCall(PetscSectionGetOffset(cSec, p, &conOff));
1303: }
1304: if (!conDof) continue;
1305: perm = (perms && perms[i]) ? perms[i][o] : NULL;
1306: flip = (flips && flips[i]) ? flips[i][o] : NULL;
1307: PetscCall(PetscSectionGetDof(aSec, p, &aDof));
1308: PetscCall(PetscSectionGetOffset(aSec, p, &aOff));
1309: nWork = childOffsets[i + 1] - childOffsets[i];
1310: for (k = 0; k < aDof; k++) {
1311: PetscInt a = anchors[aOff + k];
1312: PetscInt aSecDof, aSecOff;
1314: if (numFields) {
1315: PetscCall(PetscSectionGetFieldDof(section, a, f, &aSecDof));
1316: PetscCall(PetscSectionGetFieldOffset(section, a, f, &aSecOff));
1317: } else {
1318: PetscCall(PetscSectionGetDof(section, a, &aSecDof));
1319: PetscCall(PetscSectionGetOffset(section, a, &aSecOff));
1320: }
1321: if (!aSecDof) continue;
1323: for (j = 0; j < closureSizeP; j++) {
1324: PetscInt q = closureP[2 * j];
1325: PetscInt oq = closureP[2 * j + 1];
1327: if (q == a) {
1328: PetscInt r, s, nWorkP;
1329: const PetscInt *permP;
1330: const PetscScalar *flipP;
1332: permP = (perms && perms[j]) ? perms[j][oq] : NULL;
1333: flipP = (flips && flips[j]) ? flips[j][oq] : NULL;
1334: nWorkP = parentOffsets[j + 1] - parentOffsets[j];
1335: /* get a copy of the child-to-anchor portion of the matrix, and transpose so that rows correspond to the
1336: * child and columns correspond to the anchor: BUT the maxrix returned by MatDenseGetArrayRead() is
1337: * column-major, so transpose-transpose = do nothing */
1338: for (r = 0; r < nWork; r++) {
1339: for (s = 0; s < nWorkP; s++) scwork[r * nWorkP + s] = X[fSize * (r + childOffsets[i]) + (s + parentOffsets[j])];
1340: }
1341: for (r = 0; r < nWork; r++) workIndRow[perm ? perm[r] : r] = conOff + r;
1342: for (s = 0; s < nWorkP; s++) workIndCol[permP ? permP[s] : s] = aSecOff + s;
1343: if (flip) {
1344: for (r = 0; r < nWork; r++) {
1345: for (s = 0; s < nWorkP; s++) scwork[r * nWorkP + s] *= flip[r];
1346: }
1347: }
1348: if (flipP) {
1349: for (r = 0; r < nWork; r++) {
1350: for (s = 0; s < nWorkP; s++) scwork[r * nWorkP + s] *= flipP[s];
1351: }
1352: }
1353: PetscCall(MatSetValues(cMat, nWork, workIndRow, nWorkP, workIndCol, scwork, INSERT_VALUES));
1354: break;
1355: }
1356: }
1357: }
1358: }
1359: PetscCall(MatDenseRestoreArrayRead(Xmat, &X));
1360: PetscCall(PetscFree2(childOffsets, parentOffsets));
1361: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
1362: PetscCall(DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSizeP, &closureP));
1363: }
1364: PetscCall(MatDestroy(&Amat));
1365: PetscCall(MatDestroy(&Bmat));
1366: PetscCall(MatDestroy(&Xmat));
1367: PetscCall(PetscFree(scwork));
1368: PetscCall(PetscFree7(sizes, weights, pointsRef, pointsReal, work, workIndRow, workIndCol));
1369: if (id == PETSCFV_CLASSID) PetscCall(PetscSpaceDestroy(&bspace));
1370: }
1371: PetscCall(MatAssemblyBegin(cMat, MAT_FINAL_ASSEMBLY));
1372: PetscCall(MatAssemblyEnd(cMat, MAT_FINAL_ASSEMBLY));
1373: PetscCall(PetscFree6(v0, v0parent, vtmp, J, Jparent, invJparent));
1374: PetscCall(ISRestoreIndices(aIS, &anchors));
1375: PetscFunctionReturn(PETSC_SUCCESS);
1376: }
1378: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1379: {
1380: Mat refCmat;
1381: PetscDS ds;
1382: PetscInt numFields, maxFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, **refPointFieldN;
1383: PetscScalar ***refPointFieldMats;
1384: PetscSection refConSec, refAnSec, refSection;
1385: IS refAnIS;
1386: const PetscInt *refAnchors;
1387: const PetscInt **perms;
1388: const PetscScalar **flips;
1390: PetscFunctionBegin;
1391: PetscCall(DMGetDS(refTree, &ds));
1392: PetscCall(PetscDSGetNumFields(ds, &numFields));
1393: maxFields = PetscMax(1, numFields);
1394: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, &refCmat, NULL));
1395: PetscCall(DMPlexGetAnchors(refTree, &refAnSec, &refAnIS));
1396: PetscCall(ISGetIndices(refAnIS, &refAnchors));
1397: PetscCall(DMGetLocalSection(refTree, &refSection));
1398: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
1399: PetscCall(PetscMalloc1(pRefEnd - pRefStart, &refPointFieldMats));
1400: PetscCall(PetscMalloc1(pRefEnd - pRefStart, &refPointFieldN));
1401: PetscCall(PetscSectionGetMaxDof(refConSec, &maxDof));
1402: PetscCall(PetscSectionGetMaxDof(refAnSec, &maxAnDof));
1403: PetscCall(PetscMalloc1(maxDof, &rows));
1404: PetscCall(PetscMalloc1(maxDof * maxAnDof, &cols));
1405: for (p = pRefStart; p < pRefEnd; p++) {
1406: PetscInt parent, closureSize, *closure = NULL, pDof;
1408: PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL));
1409: PetscCall(PetscSectionGetDof(refConSec, p, &pDof));
1410: if (!pDof || parent == p) continue;
1412: PetscCall(PetscMalloc1(maxFields, &refPointFieldMats[p - pRefStart]));
1413: PetscCall(PetscCalloc1(maxFields, &refPointFieldN[p - pRefStart]));
1414: PetscCall(DMPlexGetTransitiveClosure(refTree, parent, PETSC_TRUE, &closureSize, &closure));
1415: for (f = 0; f < maxFields; f++) {
1416: PetscInt cDof, cOff, numCols, r, i;
1418: if (f < numFields) {
1419: PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof));
1420: PetscCall(PetscSectionGetFieldOffset(refConSec, p, f, &cOff));
1421: PetscCall(PetscSectionGetFieldPointSyms(refSection, f, closureSize, closure, &perms, &flips));
1422: } else {
1423: PetscCall(PetscSectionGetDof(refConSec, p, &cDof));
1424: PetscCall(PetscSectionGetOffset(refConSec, p, &cOff));
1425: PetscCall(PetscSectionGetPointSyms(refSection, closureSize, closure, &perms, &flips));
1426: }
1428: for (r = 0; r < cDof; r++) rows[r] = cOff + r;
1429: numCols = 0;
1430: for (i = 0; i < closureSize; i++) {
1431: PetscInt q = closure[2 * i];
1432: PetscInt aDof, aOff, j;
1433: const PetscInt *perm = perms ? perms[i] : NULL;
1435: if (numFields) {
1436: PetscCall(PetscSectionGetFieldDof(refSection, q, f, &aDof));
1437: PetscCall(PetscSectionGetFieldOffset(refSection, q, f, &aOff));
1438: } else {
1439: PetscCall(PetscSectionGetDof(refSection, q, &aDof));
1440: PetscCall(PetscSectionGetOffset(refSection, q, &aOff));
1441: }
1443: for (j = 0; j < aDof; j++) cols[numCols++] = aOff + (perm ? perm[j] : j);
1444: }
1445: refPointFieldN[p - pRefStart][f] = numCols;
1446: PetscCall(PetscMalloc1(cDof * numCols, &refPointFieldMats[p - pRefStart][f]));
1447: PetscCall(MatGetValues(refCmat, cDof, rows, numCols, cols, refPointFieldMats[p - pRefStart][f]));
1448: if (flips) {
1449: PetscInt colOff = 0;
1451: for (i = 0; i < closureSize; i++) {
1452: PetscInt q = closure[2 * i];
1453: PetscInt aDof, aOff, j;
1454: const PetscScalar *flip = flips ? flips[i] : NULL;
1456: if (numFields) {
1457: PetscCall(PetscSectionGetFieldDof(refSection, q, f, &aDof));
1458: PetscCall(PetscSectionGetFieldOffset(refSection, q, f, &aOff));
1459: } else {
1460: PetscCall(PetscSectionGetDof(refSection, q, &aDof));
1461: PetscCall(PetscSectionGetOffset(refSection, q, &aOff));
1462: }
1463: if (flip) {
1464: for (PetscInt k = 0; k < cDof; k++) {
1465: for (j = 0; j < aDof; j++) refPointFieldMats[p - pRefStart][f][k * numCols + colOff + j] *= flip[j];
1466: }
1467: }
1468: colOff += aDof;
1469: }
1470: }
1471: if (numFields) {
1472: PetscCall(PetscSectionRestoreFieldPointSyms(refSection, f, closureSize, closure, &perms, &flips));
1473: } else {
1474: PetscCall(PetscSectionRestorePointSyms(refSection, closureSize, closure, &perms, &flips));
1475: }
1476: }
1477: PetscCall(DMPlexRestoreTransitiveClosure(refTree, parent, PETSC_TRUE, &closureSize, &closure));
1478: }
1479: *childrenMats = refPointFieldMats;
1480: *childrenN = refPointFieldN;
1481: PetscCall(ISRestoreIndices(refAnIS, &refAnchors));
1482: PetscCall(PetscFree(rows));
1483: PetscCall(PetscFree(cols));
1484: PetscFunctionReturn(PETSC_SUCCESS);
1485: }
1487: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1488: {
1489: PetscDS ds;
1490: PetscInt **refPointFieldN;
1491: PetscScalar ***refPointFieldMats;
1492: PetscInt numFields, maxFields, pRefStart, pRefEnd, p, f;
1493: PetscSection refConSec;
1495: PetscFunctionBegin;
1496: refPointFieldN = *childrenN;
1497: *childrenN = NULL;
1498: refPointFieldMats = *childrenMats;
1499: *childrenMats = NULL;
1500: PetscCall(DMGetDS(refTree, &ds));
1501: PetscCall(PetscDSGetNumFields(ds, &numFields));
1502: maxFields = PetscMax(1, numFields);
1503: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
1504: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
1505: for (p = pRefStart; p < pRefEnd; p++) {
1506: PetscInt parent, pDof;
1508: PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL));
1509: PetscCall(PetscSectionGetDof(refConSec, p, &pDof));
1510: if (!pDof || parent == p) continue;
1512: for (f = 0; f < maxFields; f++) {
1513: PetscInt cDof;
1515: if (numFields) {
1516: PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof));
1517: } else {
1518: PetscCall(PetscSectionGetDof(refConSec, p, &cDof));
1519: }
1521: PetscCall(PetscFree(refPointFieldMats[p - pRefStart][f]));
1522: }
1523: PetscCall(PetscFree(refPointFieldMats[p - pRefStart]));
1524: PetscCall(PetscFree(refPointFieldN[p - pRefStart]));
1525: }
1526: PetscCall(PetscFree(refPointFieldMats));
1527: PetscCall(PetscFree(refPointFieldN));
1528: PetscFunctionReturn(PETSC_SUCCESS);
1529: }
1531: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm, PetscSection section, PetscSection conSec, Mat cMat)
1532: {
1533: DM refTree;
1534: PetscDS ds;
1535: Mat refCmat;
1536: PetscInt numFields, maxFields, f, pRefStart, pRefEnd, p, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN;
1537: PetscScalar ***refPointFieldMats, *pointWork;
1538: PetscSection refConSec, refAnSec, anSec;
1539: IS refAnIS, anIS;
1540: const PetscInt *anchors;
1542: PetscFunctionBegin;
1544: PetscCall(DMGetDS(dm, &ds));
1545: PetscCall(PetscDSGetNumFields(ds, &numFields));
1546: maxFields = PetscMax(1, numFields);
1547: PetscCall(DMPlexGetReferenceTree(dm, &refTree));
1548: PetscCall(DMCopyDisc(dm, refTree));
1549: PetscCall(DMSetLocalSection(refTree, NULL));
1550: PetscCall(DMSetDefaultConstraints(refTree, NULL, NULL, NULL));
1551: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, &refCmat, NULL));
1552: PetscCall(DMPlexGetAnchors(refTree, &refAnSec, &refAnIS));
1553: PetscCall(DMPlexGetAnchors(dm, &anSec, &anIS));
1554: PetscCall(ISGetIndices(anIS, &anchors));
1555: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
1556: PetscCall(PetscSectionGetChart(conSec, &conStart, &conEnd));
1557: PetscCall(PetscSectionGetMaxDof(refConSec, &maxDof));
1558: PetscCall(PetscSectionGetMaxDof(refAnSec, &maxAnDof));
1559: PetscCall(PetscMalloc1(maxDof * maxDof * maxAnDof, &pointWork));
1561: /* step 1: get submats for every constrained point in the reference tree */
1562: PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
1564: /* step 2: compute the preorder */
1565: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1566: PetscCall(PetscMalloc2(pEnd - pStart, &perm, pEnd - pStart, &iperm));
1567: for (p = pStart; p < pEnd; p++) {
1568: perm[p - pStart] = p;
1569: iperm[p - pStart] = p - pStart;
1570: }
1571: for (p = 0; p < pEnd - pStart;) {
1572: PetscInt point = perm[p];
1573: PetscInt parent;
1575: PetscCall(DMPlexGetTreeParent(dm, point, &parent, NULL));
1576: if (parent == point) {
1577: p++;
1578: } else {
1579: PetscInt size, closureSize, *closure = NULL, i;
1581: PetscCall(DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
1582: for (i = 0; i < closureSize; i++) {
1583: PetscInt q = closure[2 * i];
1584: if (iperm[q - pStart] > iperm[point - pStart]) {
1585: /* swap */
1586: perm[p] = q;
1587: perm[iperm[q - pStart]] = point;
1588: iperm[point - pStart] = iperm[q - pStart];
1589: iperm[q - pStart] = p;
1590: break;
1591: }
1592: }
1593: size = closureSize;
1594: PetscCall(DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
1595: if (i == size) p++;
1596: }
1597: }
1599: /* step 3: fill the constraint matrix */
1600: /* we are going to use a preorder progressive fill strategy. Mat doesn't
1601: * allow progressive fill without assembly, so we are going to set up the
1602: * values outside of the Mat first.
1603: */
1604: {
1605: PetscInt nRows, row, nnz;
1606: PetscBool done;
1607: PetscInt secStart, secEnd;
1608: const PetscInt *ia, *ja;
1609: PetscScalar *vals;
1611: PetscCall(PetscSectionGetChart(section, &secStart, &secEnd));
1612: PetscCall(MatGetRowIJ(cMat, 0, PETSC_FALSE, PETSC_FALSE, &nRows, &ia, &ja, &done));
1613: PetscCheck(done, PetscObjectComm((PetscObject)cMat), PETSC_ERR_PLIB, "Could not get RowIJ of constraint matrix");
1614: nnz = ia[nRows];
1615: /* malloc and then zero rows right before we fill them: this way valgrind
1616: * can tell if we are doing progressive fill in the wrong order */
1617: PetscCall(PetscMalloc1(nnz, &vals));
1618: for (p = 0; p < pEnd - pStart; p++) {
1619: PetscInt parent, childid, closureSize, *closure = NULL;
1620: PetscInt point = perm[p], pointDof;
1622: PetscCall(DMPlexGetTreeParent(dm, point, &parent, &childid));
1623: if ((point < conStart) || (point >= conEnd) || (parent == point)) continue;
1624: PetscCall(PetscSectionGetDof(conSec, point, &pointDof));
1625: if (!pointDof) continue;
1626: PetscCall(DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
1627: for (f = 0; f < maxFields; f++) {
1628: PetscInt cDof, cOff, numCols, numFillCols, i, r, matOffset, offset;
1629: PetscScalar *pointMat;
1630: const PetscInt **perms;
1631: const PetscScalar **flips;
1633: if (numFields) {
1634: PetscCall(PetscSectionGetFieldDof(conSec, point, f, &cDof));
1635: PetscCall(PetscSectionGetFieldOffset(conSec, point, f, &cOff));
1636: } else {
1637: PetscCall(PetscSectionGetDof(conSec, point, &cDof));
1638: PetscCall(PetscSectionGetOffset(conSec, point, &cOff));
1639: }
1640: if (!cDof) continue;
1641: if (numFields) PetscCall(PetscSectionGetFieldPointSyms(section, f, closureSize, closure, &perms, &flips));
1642: else PetscCall(PetscSectionGetPointSyms(section, closureSize, closure, &perms, &flips));
1644: /* make sure that every row for this point is the same size */
1645: if (PetscDefined(USE_DEBUG)) {
1646: for (r = 0; r < cDof; r++) {
1647: if (cDof > 1 && r) {
1648: 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]);
1649: }
1650: }
1651: }
1652: /* zero rows */
1653: for (i = ia[cOff]; i < ia[cOff + cDof]; i++) vals[i] = 0.;
1654: matOffset = ia[cOff];
1655: numFillCols = ia[cOff + 1] - matOffset;
1656: pointMat = refPointFieldMats[childid - pRefStart][f];
1657: numCols = refPointFieldN[childid - pRefStart][f];
1658: offset = 0;
1659: for (i = 0; i < closureSize; i++) {
1660: PetscInt q = closure[2 * i];
1661: PetscInt aDof, aOff, j, k, qConDof, qConOff;
1662: const PetscInt *perm = perms ? perms[i] : NULL;
1663: const PetscScalar *flip = flips ? flips[i] : NULL;
1665: qConDof = qConOff = 0;
1666: if (q < secStart || q >= secEnd) continue;
1667: if (numFields) {
1668: PetscCall(PetscSectionGetFieldDof(section, q, f, &aDof));
1669: PetscCall(PetscSectionGetFieldOffset(section, q, f, &aOff));
1670: if (q >= conStart && q < conEnd) {
1671: PetscCall(PetscSectionGetFieldDof(conSec, q, f, &qConDof));
1672: PetscCall(PetscSectionGetFieldOffset(conSec, q, f, &qConOff));
1673: }
1674: } else {
1675: PetscCall(PetscSectionGetDof(section, q, &aDof));
1676: PetscCall(PetscSectionGetOffset(section, q, &aOff));
1677: if (q >= conStart && q < conEnd) {
1678: PetscCall(PetscSectionGetDof(conSec, q, &qConDof));
1679: PetscCall(PetscSectionGetOffset(conSec, q, &qConOff));
1680: }
1681: }
1682: if (!aDof) continue;
1683: if (qConDof) {
1684: /* this point has anchors: its rows of the matrix should already
1685: * be filled, thanks to preordering */
1686: /* first multiply into pointWork, then set in matrix */
1687: PetscInt aMatOffset = ia[qConOff];
1688: PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset;
1689: for (r = 0; r < cDof; r++) {
1690: for (j = 0; j < aNumFillCols; j++) {
1691: PetscScalar inVal = 0;
1692: for (k = 0; k < aDof; k++) {
1693: PetscInt col = perm ? perm[k] : k;
1695: inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j] * (flip ? flip[col] : 1.);
1696: }
1697: pointWork[r * aNumFillCols + j] = inVal;
1698: }
1699: }
1700: /* assume that the columns are sorted, spend less time searching */
1701: for (j = 0, k = 0; j < aNumFillCols; j++) {
1702: PetscInt col = ja[aMatOffset + j];
1703: for (; k < numFillCols; k++) {
1704: if (ja[matOffset + k] == col) break;
1705: }
1706: PetscCheck(k != numFillCols, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No nonzero space for (%" PetscInt_FMT ", %" PetscInt_FMT ")", cOff, col);
1707: for (r = 0; r < cDof; r++) vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j];
1708: }
1709: } else {
1710: /* find where to put this portion of pointMat into the matrix */
1711: for (k = 0; k < numFillCols; k++) {
1712: if (ja[matOffset + k] == aOff) break;
1713: }
1714: PetscCheck(k != numFillCols, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No nonzero space for (%" PetscInt_FMT ", %" PetscInt_FMT ")", cOff, aOff);
1715: for (r = 0; r < cDof; r++) {
1716: for (j = 0; j < aDof; j++) {
1717: PetscInt col = perm ? perm[j] : j;
1719: vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + j] * (flip ? flip[col] : 1.);
1720: }
1721: }
1722: }
1723: offset += aDof;
1724: }
1725: if (numFields) {
1726: PetscCall(PetscSectionRestoreFieldPointSyms(section, f, closureSize, closure, &perms, &flips));
1727: } else {
1728: PetscCall(PetscSectionRestorePointSyms(section, closureSize, closure, &perms, &flips));
1729: }
1730: }
1731: PetscCall(DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
1732: }
1733: for (row = 0; row < nRows; row++) PetscCall(MatSetValues(cMat, 1, &row, ia[row + 1] - ia[row], &ja[ia[row]], &vals[ia[row]], INSERT_VALUES));
1734: PetscCall(MatRestoreRowIJ(cMat, 0, PETSC_FALSE, PETSC_FALSE, &nRows, &ia, &ja, &done));
1735: PetscCheck(done, PetscObjectComm((PetscObject)cMat), PETSC_ERR_PLIB, "Could not restore RowIJ of constraint matrix");
1736: PetscCall(MatAssemblyBegin(cMat, MAT_FINAL_ASSEMBLY));
1737: PetscCall(MatAssemblyEnd(cMat, MAT_FINAL_ASSEMBLY));
1738: PetscCall(PetscFree(vals));
1739: }
1741: /* clean up */
1742: PetscCall(ISRestoreIndices(anIS, &anchors));
1743: PetscCall(PetscFree2(perm, iperm));
1744: PetscCall(PetscFree(pointWork));
1745: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
1746: PetscFunctionReturn(PETSC_SUCCESS);
1747: }
1749: /*@
1750: DMPlexTreeRefineCell - Refine a single cell on rank 0 using the `DM`'s reference tree, producing a non-conforming mesh
1752: Collective
1754: Input Parameters:
1755: + dm - The `DM` with an attached reference tree (see `DMPlexSetReferenceTree()`)
1756: - cell - The cell to be refined
1758: Output Parameter:
1759: . ncdm - A new `DM` in which `cell` has been split according to the reference tree
1761: Level: developer
1763: Note:
1764: This routine is intended for testing and demonstration; it produces one example of a non-conforming
1765: mesh but is not a general local-refinement facility. Only rank 0 performs the refinement.
1767: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetReferenceTree()`, `DMPlexGetReferenceTree()`, `DMPlexSetTree()`
1768: @*/
1769: /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of
1770: * a non-conforming mesh. Local refinement comes later */
1771: PetscErrorCode DMPlexTreeRefineCell(DM dm, PetscInt cell, DM *ncdm)
1772: {
1773: DM K;
1774: PetscMPIInt rank;
1775: PetscInt dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd;
1776: PetscInt numNewCones, *newConeSizes, *newCones, *newOrientations;
1777: PetscInt *Kembedding;
1778: PetscInt *cellClosure = NULL, nc;
1779: PetscScalar *newVertexCoords;
1780: PetscInt numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset;
1781: PetscSection parentSection;
1783: PetscFunctionBegin;
1784: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1785: PetscCall(DMGetDimension(dm, &dim));
1786: PetscCall(DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm));
1787: PetscCall(DMSetDimension(*ncdm, dim));
1789: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1790: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &parentSection));
1791: PetscCall(DMPlexGetReferenceTree(dm, &K));
1792: PetscCall(DMGetCoordinatesLocalSetUp(dm));
1793: if (rank == 0) {
1794: /* compute the new charts */
1795: PetscCall(PetscMalloc5(dim + 1, &pNewCount, dim + 1, &pNewStart, dim + 1, &pNewEnd, dim + 1, &pOldStart, dim + 1, &pOldEnd));
1796: offset = 0;
1797: for (d = 0; d <= dim; d++) {
1798: PetscInt pOldCount, kStart, kEnd, k;
1800: pNewStart[d] = offset;
1801: PetscCall(DMPlexGetHeightStratum(dm, d, &pOldStart[d], &pOldEnd[d]));
1802: PetscCall(DMPlexGetHeightStratum(K, d, &kStart, &kEnd));
1803: pOldCount = pOldEnd[d] - pOldStart[d];
1804: /* adding the new points */
1805: pNewCount[d] = pOldCount + kEnd - kStart;
1806: if (!d) {
1807: /* removing the cell */
1808: pNewCount[d]--;
1809: }
1810: for (k = kStart; k < kEnd; k++) {
1811: PetscInt parent;
1812: PetscCall(DMPlexGetTreeParent(K, k, &parent, NULL));
1813: if (parent == k) {
1814: /* avoid double counting points that won't actually be new */
1815: pNewCount[d]--;
1816: }
1817: }
1818: pNewEnd[d] = pNewStart[d] + pNewCount[d];
1819: offset = pNewEnd[d];
1820: }
1821: 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]);
1822: /* get the current closure of the cell that we are removing */
1823: PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &nc, &cellClosure));
1825: PetscCall(PetscMalloc1(pNewEnd[dim], &newConeSizes));
1826: {
1827: DMPolytopeType pct, qct;
1828: PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j;
1830: PetscCall(DMPlexGetChart(K, &kStart, &kEnd));
1831: PetscCall(PetscMalloc4(kEnd - kStart, &Kembedding, kEnd - kStart, &perm, kEnd - kStart, &iperm, kEnd - kStart, &preOrient));
1833: for (k = kStart; k < kEnd; k++) {
1834: perm[k - kStart] = k;
1835: iperm[k - kStart] = k - kStart;
1836: preOrient[k - kStart] = 0;
1837: }
1839: PetscCall(DMPlexGetTransitiveClosure(K, 0, PETSC_TRUE, &closureSizeK, &closureK));
1840: for (j = 1; j < closureSizeK; j++) {
1841: PetscInt parentOrientA = closureK[2 * j + 1];
1842: PetscInt parentOrientB = cellClosure[2 * j + 1];
1843: PetscInt p, q;
1845: p = closureK[2 * j];
1846: q = cellClosure[2 * j];
1847: PetscCall(DMPlexGetCellType(K, p, &pct));
1848: PetscCall(DMPlexGetCellType(dm, q, &qct));
1849: for (d = 0; d <= dim; d++) {
1850: if (q >= pOldStart[d] && q < pOldEnd[d]) Kembedding[p] = (q - pOldStart[d]) + pNewStart[d];
1851: }
1852: parentOrientA = DMPolytopeConvertNewOrientation_Internal(pct, parentOrientA);
1853: parentOrientB = DMPolytopeConvertNewOrientation_Internal(qct, parentOrientB);
1854: if (parentOrientA != parentOrientB) {
1855: PetscInt numChildren;
1856: const PetscInt *children;
1858: PetscCall(DMPlexGetTreeChildren(K, p, &numChildren, &children));
1859: for (PetscInt i = 0; i < numChildren; i++) {
1860: PetscInt kPerm, oPerm;
1862: k = children[i];
1863: PetscCall(DMPlexReferenceTreeGetChildSymmetry(K, p, parentOrientA, 0, k, parentOrientB, &oPerm, &kPerm));
1864: /* perm = what refTree position I'm in */
1865: perm[kPerm - kStart] = k;
1866: /* iperm = who is at this position */
1867: iperm[k - kStart] = kPerm - kStart;
1868: preOrient[kPerm - kStart] = oPerm;
1869: }
1870: }
1871: }
1872: PetscCall(DMPlexRestoreTransitiveClosure(K, 0, PETSC_TRUE, &closureSizeK, &closureK));
1873: }
1874: PetscCall(PetscSectionSetChart(parentSection, 0, pNewEnd[dim]));
1875: offset = 0;
1876: numNewCones = 0;
1877: for (d = 0; d <= dim; d++) {
1878: PetscInt kStart, kEnd, k;
1879: PetscInt p;
1880: PetscInt size;
1882: for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1883: /* skip cell 0 */
1884: if (p == cell) continue;
1885: /* old cones to new cones */
1886: PetscCall(DMPlexGetConeSize(dm, p, &size));
1887: newConeSizes[offset++] = size;
1888: numNewCones += size;
1889: }
1891: PetscCall(DMPlexGetHeightStratum(K, d, &kStart, &kEnd));
1892: for (k = kStart; k < kEnd; k++) {
1893: PetscInt kParent;
1895: PetscCall(DMPlexGetTreeParent(K, k, &kParent, NULL));
1896: if (kParent != k) {
1897: Kembedding[k] = offset;
1898: PetscCall(DMPlexGetConeSize(K, k, &size));
1899: newConeSizes[offset++] = size;
1900: numNewCones += size;
1901: if (kParent != 0) PetscCall(PetscSectionSetDof(parentSection, Kembedding[k], 1));
1902: }
1903: }
1904: }
1906: PetscCall(PetscSectionSetUp(parentSection));
1907: PetscCall(PetscSectionGetStorageSize(parentSection, &numPointsWithParents));
1908: PetscCall(PetscMalloc2(numNewCones, &newCones, numNewCones, &newOrientations));
1909: PetscCall(PetscMalloc2(numPointsWithParents, &parents, numPointsWithParents, &childIDs));
1911: /* fill new cones */
1912: offset = 0;
1913: for (d = 0; d <= dim; d++) {
1914: PetscInt kStart, kEnd, k, l;
1915: PetscInt p;
1916: PetscInt size;
1917: const PetscInt *cone, *orientation;
1919: for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1920: /* skip cell 0 */
1921: if (p == cell) continue;
1922: /* old cones to new cones */
1923: PetscCall(DMPlexGetConeSize(dm, p, &size));
1924: PetscCall(DMPlexGetCone(dm, p, &cone));
1925: PetscCall(DMPlexGetConeOrientation(dm, p, &orientation));
1926: for (l = 0; l < size; l++) {
1927: newCones[offset] = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1];
1928: newOrientations[offset++] = orientation[l];
1929: }
1930: }
1932: PetscCall(DMPlexGetHeightStratum(K, d, &kStart, &kEnd));
1933: for (k = kStart; k < kEnd; k++) {
1934: PetscInt kPerm = perm[k], kParent;
1935: PetscInt preO = preOrient[k];
1937: PetscCall(DMPlexGetTreeParent(K, k, &kParent, NULL));
1938: if (kParent != k) {
1939: /* embed new cones */
1940: PetscCall(DMPlexGetConeSize(K, k, &size));
1941: PetscCall(DMPlexGetCone(K, kPerm, &cone));
1942: PetscCall(DMPlexGetConeOrientation(K, kPerm, &orientation));
1943: for (l = 0; l < size; l++) {
1944: PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size - (preO + 1) - l) % size);
1945: PetscInt newO, lSize, oTrue;
1946: DMPolytopeType ct = DM_NUM_POLYTOPES;
1948: q = iperm[cone[m]];
1949: newCones[offset] = Kembedding[q];
1950: PetscCall(DMPlexGetConeSize(K, q, &lSize));
1951: if (lSize == 2) ct = DM_POLYTOPE_SEGMENT;
1952: else if (lSize == 4) ct = DM_POLYTOPE_QUADRILATERAL;
1953: oTrue = DMPolytopeConvertNewOrientation_Internal(ct, orientation[m]);
1954: oTrue = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2);
1955: newO = DihedralCompose(lSize, oTrue, preOrient[q]);
1956: newOrientations[offset++] = DMPolytopeConvertOldOrientation_Internal(ct, newO);
1957: }
1958: if (kParent != 0) {
1959: PetscInt newPoint = Kembedding[kParent];
1960: PetscCall(PetscSectionGetOffset(parentSection, Kembedding[k], &pOffset));
1961: parents[pOffset] = newPoint;
1962: childIDs[pOffset] = k;
1963: }
1964: }
1965: }
1966: }
1968: PetscCall(PetscMalloc1(dim * (pNewEnd[dim] - pNewStart[dim]), &newVertexCoords));
1970: /* fill coordinates */
1971: offset = 0;
1972: {
1973: PetscInt kStart, kEnd, l;
1974: PetscSection vSection;
1975: Vec coords;
1976: PetscScalar *coordvals;
1977: PetscInt dof, off;
1978: PetscReal v0[3], J[9], detJ;
1980: if (PetscDefined(USE_DEBUG)) {
1981: PetscCall(DMPlexGetHeightStratum(K, 0, &kStart, &kEnd));
1982: for (PetscInt k = kStart; k < kEnd; k++) {
1983: PetscCall(DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ));
1984: PetscCheck(detJ > 0., PETSC_COMM_SELF, PETSC_ERR_PLIB, "reference tree cell %" PetscInt_FMT " has bad determinant", k);
1985: }
1986: }
1987: PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ));
1988: PetscCall(DMGetCoordinateSection(dm, &vSection));
1989: PetscCall(DMGetCoordinatesLocal(dm, &coords));
1990: PetscCall(VecGetArray(coords, &coordvals));
1991: for (PetscInt v = pOldStart[dim]; v < pOldEnd[dim]; v++) {
1992: PetscCall(PetscSectionGetDof(vSection, v, &dof));
1993: PetscCall(PetscSectionGetOffset(vSection, v, &off));
1994: for (l = 0; l < dof; l++) newVertexCoords[offset++] = coordvals[off + l];
1995: }
1996: PetscCall(VecRestoreArray(coords, &coordvals));
1998: PetscCall(DMGetCoordinateSection(K, &vSection));
1999: PetscCall(DMGetCoordinatesLocal(K, &coords));
2000: PetscCall(VecGetArray(coords, &coordvals));
2001: PetscCall(DMPlexGetDepthStratum(K, 0, &kStart, &kEnd));
2002: for (PetscInt v = kStart; v < kEnd; v++) {
2003: PetscReal coord[3], newCoord[3];
2004: PetscInt vPerm = perm[v];
2005: PetscInt kParent;
2006: const PetscReal xi0[3] = {-1., -1., -1.};
2008: PetscCall(DMPlexGetTreeParent(K, v, &kParent, NULL));
2009: if (kParent != v) {
2010: /* this is a new vertex */
2011: PetscCall(PetscSectionGetOffset(vSection, vPerm, &off));
2012: for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off + l]);
2013: CoordinatesRefToReal(dim, dim, xi0, v0, J, coord, newCoord);
2014: for (l = 0; l < dim; ++l) newVertexCoords[offset + l] = newCoord[l];
2015: offset += dim;
2016: }
2017: }
2018: PetscCall(VecRestoreArray(coords, &coordvals));
2019: }
2021: /* need to reverse the order of pNewCount: vertices first, cells last */
2022: for (d = 0; d < (dim + 1) / 2; d++) {
2023: PetscInt tmp;
2025: tmp = pNewCount[d];
2026: pNewCount[d] = pNewCount[dim - d];
2027: pNewCount[dim - d] = tmp;
2028: }
2030: PetscCall(DMPlexCreateFromDAG(*ncdm, dim, pNewCount, newConeSizes, newCones, newOrientations, newVertexCoords));
2031: PetscCall(DMPlexSetReferenceTree(*ncdm, K));
2032: PetscCall(DMPlexSetTree(*ncdm, parentSection, parents, childIDs));
2034: /* clean up */
2035: PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &nc, &cellClosure));
2036: PetscCall(PetscFree5(pNewCount, pNewStart, pNewEnd, pOldStart, pOldEnd));
2037: PetscCall(PetscFree(newConeSizes));
2038: PetscCall(PetscFree2(newCones, newOrientations));
2039: PetscCall(PetscFree(newVertexCoords));
2040: PetscCall(PetscFree2(parents, childIDs));
2041: PetscCall(PetscFree4(Kembedding, perm, iperm, preOrient));
2042: } else {
2043: PetscInt p, counts[4];
2044: PetscInt *coneSizes, *cones, *orientations;
2045: Vec coordVec;
2046: PetscScalar *coords;
2048: for (d = 0; d <= dim; d++) {
2049: PetscInt dStart, dEnd;
2051: PetscCall(DMPlexGetDepthStratum(dm, d, &dStart, &dEnd));
2052: counts[d] = dEnd - dStart;
2053: }
2054: PetscCall(PetscMalloc1(pEnd - pStart, &coneSizes));
2055: for (p = pStart; p < pEnd; p++) PetscCall(DMPlexGetConeSize(dm, p, &coneSizes[p - pStart]));
2056: PetscCall(DMPlexGetCones(dm, &cones));
2057: PetscCall(DMPlexGetConeOrientations(dm, &orientations));
2058: PetscCall(DMGetCoordinatesLocal(dm, &coordVec));
2059: PetscCall(VecGetArray(coordVec, &coords));
2061: PetscCall(PetscSectionSetChart(parentSection, pStart, pEnd));
2062: PetscCall(PetscSectionSetUp(parentSection));
2063: PetscCall(DMPlexCreateFromDAG(*ncdm, dim, counts, coneSizes, cones, orientations, NULL));
2064: PetscCall(DMPlexSetReferenceTree(*ncdm, K));
2065: PetscCall(DMPlexSetTree(*ncdm, parentSection, NULL, NULL));
2066: PetscCall(VecRestoreArray(coordVec, &coords));
2067: }
2068: PetscCall(PetscSectionDestroy(&parentSection));
2069: PetscFunctionReturn(PETSC_SUCCESS);
2070: }
2072: PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
2073: {
2074: PetscSF coarseToFineEmbedded;
2075: PetscSection globalCoarse, globalFine;
2076: PetscSection localCoarse, localFine;
2077: PetscSection aSec, cSec;
2078: PetscSection rootIndicesSec, rootMatricesSec;
2079: PetscSection leafIndicesSec, leafMatricesSec;
2080: PetscInt *rootIndices, *leafIndices;
2081: PetscScalar *rootMatrices, *leafMatrices;
2082: IS aIS;
2083: const PetscInt *anchors;
2084: Mat cMat;
2085: PetscInt numFields, maxFields;
2086: PetscInt pStartC, pEndC, pStartF, pEndF, p;
2087: PetscInt aStart, aEnd, cStart, cEnd;
2088: PetscInt *maxChildIds;
2089: PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
2090: const PetscInt ***perms;
2091: const PetscScalar ***flips;
2093: PetscFunctionBegin;
2094: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
2095: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
2096: PetscCall(DMGetGlobalSection(fine, &globalFine));
2097: { /* winnow fine points that don't have global dofs out of the sf */
2098: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, nleaves, l;
2099: const PetscInt *leaves;
2101: PetscCall(PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL));
2102: for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
2103: p = leaves ? leaves[l] : l;
2104: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
2105: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
2106: if ((dof - cdof) > 0) numPointsWithDofs++;
2107: }
2108: PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
2109: for (l = 0, offset = 0; l < nleaves; l++) {
2110: p = leaves ? leaves[l] : l;
2111: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
2112: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
2113: if ((dof - cdof) > 0) pointsWithDofs[offset++] = l;
2114: }
2115: PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded));
2116: PetscCall(PetscFree(pointsWithDofs));
2117: }
2118: /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
2119: PetscCall(PetscMalloc1(pEndC - pStartC, &maxChildIds));
2120: for (p = pStartC; p < pEndC; p++) maxChildIds[p - pStartC] = -2;
2121: PetscCall(PetscSFReduceBegin(coarseToFineEmbedded, MPIU_INT, childIds, maxChildIds, MPI_MAX));
2122: PetscCall(PetscSFReduceEnd(coarseToFineEmbedded, MPIU_INT, childIds, maxChildIds, MPI_MAX));
2124: PetscCall(DMGetLocalSection(coarse, &localCoarse));
2125: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
2127: PetscCall(DMPlexGetAnchors(coarse, &aSec, &aIS));
2128: PetscCall(ISGetIndices(aIS, &anchors));
2129: PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd));
2131: PetscCall(DMGetDefaultConstraints(coarse, &cSec, &cMat, NULL));
2132: PetscCall(PetscSectionGetChart(cSec, &cStart, &cEnd));
2134: /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
2135: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootIndicesSec));
2136: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootMatricesSec));
2137: PetscCall(PetscSectionSetChart(rootIndicesSec, pStartC, pEndC));
2138: PetscCall(PetscSectionSetChart(rootMatricesSec, pStartC, pEndC));
2139: PetscCall(PetscSectionGetNumFields(localCoarse, &numFields));
2140: maxFields = PetscMax(1, numFields);
2141: PetscCall(PetscMalloc7(maxFields + 1, &offsets, maxFields + 1, &offsetsCopy, maxFields + 1, &newOffsets, maxFields + 1, &newOffsetsCopy, maxFields + 1, &rowOffsets, maxFields + 1, &numD, maxFields + 1, &numO));
2142: PetscCall(PetscMalloc2(maxFields + 1, (PetscInt ****)&perms, maxFields + 1, (PetscScalar ****)&flips));
2143: PetscCall(PetscMemzero((void *)perms, (maxFields + 1) * sizeof(const PetscInt **)));
2144: PetscCall(PetscMemzero((void *)flips, (maxFields + 1) * sizeof(const PetscScalar **)));
2146: for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
2147: PetscInt dof, matSize = 0;
2148: PetscInt aDof = 0;
2149: PetscInt cDof = 0;
2150: PetscInt maxChildId = maxChildIds[p - pStartC];
2151: PetscInt numRowIndices = 0;
2152: PetscInt numColIndices = 0;
2153: PetscInt f;
2155: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
2156: if (dof < 0) dof = -(dof + 1);
2157: if (p >= aStart && p < aEnd) PetscCall(PetscSectionGetDof(aSec, p, &aDof));
2158: if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(cSec, p, &cDof));
2159: for (f = 0; f <= numFields; f++) offsets[f] = 0;
2160: for (f = 0; f <= numFields; f++) newOffsets[f] = 0;
2161: if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
2162: PetscInt *closure = NULL, closureSize, cl;
2164: PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2165: for (cl = 0; cl < closureSize; cl++) { /* get the closure */
2166: PetscInt c = closure[2 * cl], clDof;
2168: PetscCall(PetscSectionGetDof(localCoarse, c, &clDof));
2169: numRowIndices += clDof;
2170: for (f = 0; f < numFields; f++) {
2171: PetscCall(PetscSectionGetFieldDof(localCoarse, c, f, &clDof));
2172: offsets[f + 1] += clDof;
2173: }
2174: }
2175: for (f = 0; f < numFields; f++) {
2176: offsets[f + 1] += offsets[f];
2177: newOffsets[f + 1] = offsets[f + 1];
2178: }
2179: /* get the number of indices needed and their field offsets */
2180: PetscCall(DMPlexAnchorsModifyMat(coarse, localCoarse, closureSize, numRowIndices, closure, NULL, NULL, NULL, &numColIndices, NULL, NULL, newOffsets, PETSC_FALSE));
2181: PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2182: if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */
2183: numColIndices = numRowIndices;
2184: matSize = 0;
2185: } else if (numFields) { /* we send one submat for each field: sum their sizes */
2186: matSize = 0;
2187: for (f = 0; f < numFields; f++) {
2188: PetscInt numRow, numCol;
2190: numRow = offsets[f + 1] - offsets[f];
2191: numCol = newOffsets[f + 1] - newOffsets[f];
2192: matSize += numRow * numCol;
2193: }
2194: } else {
2195: matSize = numRowIndices * numColIndices;
2196: }
2197: } else if (maxChildId == -1) {
2198: if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */
2199: PetscInt aOff;
2201: PetscCall(PetscSectionGetOffset(aSec, p, &aOff));
2202: for (f = 0; f < numFields; f++) {
2203: PetscInt fDof;
2205: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
2206: offsets[f + 1] = fDof;
2207: }
2208: for (PetscInt a = 0; a < aDof; a++) {
2209: PetscInt anchor = anchors[a + aOff], aLocalDof;
2211: PetscCall(PetscSectionGetDof(localCoarse, anchor, &aLocalDof));
2212: numColIndices += aLocalDof;
2213: for (f = 0; f < numFields; f++) {
2214: PetscInt fDof;
2216: PetscCall(PetscSectionGetFieldDof(localCoarse, anchor, f, &fDof));
2217: newOffsets[f + 1] += fDof;
2218: }
2219: }
2220: if (numFields) {
2221: matSize = 0;
2222: for (f = 0; f < numFields; f++) matSize += offsets[f + 1] * newOffsets[f + 1];
2223: } else {
2224: matSize = numColIndices * dof;
2225: }
2226: } else { /* no children, and no constraints on dofs: just get the global indices */
2227: numColIndices = dof;
2228: matSize = 0;
2229: }
2230: }
2231: /* we will pack the column indices with the field offsets */
2232: PetscCall(PetscSectionSetDof(rootIndicesSec, p, numColIndices ? numColIndices + 2 * numFields : 0));
2233: PetscCall(PetscSectionSetDof(rootMatricesSec, p, matSize));
2234: }
2235: PetscCall(PetscSectionSetUp(rootIndicesSec));
2236: PetscCall(PetscSectionSetUp(rootMatricesSec));
2237: {
2238: PetscInt numRootIndices, numRootMatrices;
2240: PetscCall(PetscSectionGetStorageSize(rootIndicesSec, &numRootIndices));
2241: PetscCall(PetscSectionGetStorageSize(rootMatricesSec, &numRootMatrices));
2242: PetscCall(PetscMalloc2(numRootIndices, &rootIndices, numRootMatrices, &rootMatrices));
2243: for (p = pStartC; p < pEndC; p++) {
2244: PetscInt numRowIndices = 0, numColIndices, matSize, dof;
2245: PetscInt pIndOff, pMatOff, f;
2246: PetscInt *pInd;
2247: PetscInt maxChildId = maxChildIds[p - pStartC];
2248: PetscScalar *pMat = NULL;
2250: PetscCall(PetscSectionGetDof(rootIndicesSec, p, &numColIndices));
2251: if (!numColIndices) continue;
2252: for (f = 0; f <= numFields; f++) {
2253: offsets[f] = 0;
2254: newOffsets[f] = 0;
2255: offsetsCopy[f] = 0;
2256: newOffsetsCopy[f] = 0;
2257: }
2258: numColIndices -= 2 * numFields;
2259: PetscCall(PetscSectionGetOffset(rootIndicesSec, p, &pIndOff));
2260: pInd = &rootIndices[pIndOff];
2261: PetscCall(PetscSectionGetDof(rootMatricesSec, p, &matSize));
2262: if (matSize) {
2263: PetscCall(PetscSectionGetOffset(rootMatricesSec, p, &pMatOff));
2264: pMat = &rootMatrices[pMatOff];
2265: }
2266: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
2267: if (dof < 0) dof = -(dof + 1);
2268: if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
2269: PetscInt j;
2271: if (matSize == 0) { /* don't need to calculate the mat, just the indices */
2272: PetscInt numIndices, *indices;
2273: PetscCall(DMPlexGetClosureIndices(coarse, localCoarse, globalCoarse, p, PETSC_TRUE, &numIndices, &indices, offsets, NULL));
2274: PetscCheck(numIndices == numColIndices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "mismatching constraint indices calculations");
2275: for (PetscInt i = 0; i < numColIndices; i++) pInd[i] = indices[i];
2276: for (PetscInt i = 0; i < numFields; i++) {
2277: pInd[numColIndices + i] = offsets[i + 1];
2278: pInd[numColIndices + numFields + i] = offsets[i + 1];
2279: }
2280: PetscCall(DMPlexRestoreClosureIndices(coarse, localCoarse, globalCoarse, p, PETSC_TRUE, &numIndices, &indices, offsets, NULL));
2281: } else {
2282: PetscInt closureSize, *closure = NULL, cl;
2283: PetscScalar *pMatIn, *pMatModified;
2284: PetscInt numPoints, *points;
2286: {
2287: PetscInt *closure = NULL, closureSize, cl;
2289: PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2290: for (cl = 0; cl < closureSize; cl++) { /* get the closure */
2291: PetscInt c = closure[2 * cl], clDof;
2293: PetscCall(PetscSectionGetDof(localCoarse, c, &clDof));
2294: numRowIndices += clDof;
2295: }
2296: PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2297: }
2299: PetscCall(DMGetWorkArray(coarse, numRowIndices * numRowIndices, MPIU_SCALAR, &pMatIn));
2300: for (PetscInt i = 0; i < numRowIndices; i++) { /* initialize to the identity */
2301: for (j = 0; j < numRowIndices; j++) pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.;
2302: }
2303: PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2304: for (f = 0; f < maxFields; f++) {
2305: if (numFields) PetscCall(PetscSectionGetFieldPointSyms(localCoarse, f, closureSize, closure, &perms[f], &flips[f]));
2306: else PetscCall(PetscSectionGetPointSyms(localCoarse, closureSize, closure, &perms[f], &flips[f]));
2307: }
2308: if (numFields) {
2309: for (cl = 0; cl < closureSize; cl++) {
2310: PetscInt c = closure[2 * cl];
2312: for (f = 0; f < numFields; f++) {
2313: PetscInt fDof;
2315: PetscCall(PetscSectionGetFieldDof(localCoarse, c, f, &fDof));
2316: offsets[f + 1] += fDof;
2317: }
2318: }
2319: for (f = 0; f < numFields; f++) {
2320: offsets[f + 1] += offsets[f];
2321: newOffsets[f + 1] = offsets[f + 1];
2322: }
2323: }
2324: /* TODO : flips here ? */
2325: /* apply hanging node constraints on the right, get the new points and the new offsets */
2326: PetscCall(DMPlexAnchorsModifyMat(coarse, localCoarse, closureSize, numRowIndices, closure, perms, pMatIn, &numPoints, NULL, &points, &pMatModified, newOffsets, PETSC_FALSE));
2327: for (f = 0; f < maxFields; f++) {
2328: if (numFields) PetscCall(PetscSectionRestoreFieldPointSyms(localCoarse, f, closureSize, closure, &perms[f], &flips[f]));
2329: else PetscCall(PetscSectionRestorePointSyms(localCoarse, closureSize, closure, &perms[f], &flips[f]));
2330: }
2331: for (f = 0; f < maxFields; f++) {
2332: if (numFields) PetscCall(PetscSectionGetFieldPointSyms(localCoarse, f, numPoints, points, &perms[f], &flips[f]));
2333: else PetscCall(PetscSectionGetPointSyms(localCoarse, numPoints, points, &perms[f], &flips[f]));
2334: }
2335: if (!numFields) {
2336: for (PetscInt i = 0; i < numRowIndices * numColIndices; i++) pMat[i] = pMatModified[i];
2337: } else {
2338: PetscInt count;
2339: for (f = 0, count = 0; f < numFields; f++) {
2340: for (PetscInt i = offsets[f]; i < offsets[f + 1]; i++) {
2341: for (PetscInt j = newOffsets[f]; j < newOffsets[f + 1]; j++, count++) pMat[count] = pMatModified[i * numColIndices + j];
2342: }
2343: }
2344: }
2345: PetscCall(DMRestoreWorkArray(coarse, numRowIndices * numColIndices, MPIU_SCALAR, &pMatModified));
2346: PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2347: PetscCall(DMRestoreWorkArray(coarse, numRowIndices * numColIndices, MPIU_SCALAR, &pMatIn));
2348: if (numFields) {
2349: for (f = 0; f < numFields; f++) {
2350: pInd[numColIndices + f] = offsets[f + 1];
2351: pInd[numColIndices + numFields + f] = newOffsets[f + 1];
2352: }
2353: for (cl = 0; cl < numPoints; cl++) {
2354: PetscInt globalOff, c = points[2 * cl];
2355: PetscCall(PetscSectionGetOffset(globalCoarse, c, &globalOff));
2356: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff + 1) : globalOff, newOffsets, PETSC_FALSE, perms, cl, NULL, pInd));
2357: }
2358: } else {
2359: for (cl = 0; cl < numPoints; cl++) {
2360: PetscInt c = points[2 * cl], globalOff;
2361: const PetscInt *perm = perms[0] ? perms[0][cl] : NULL;
2363: PetscCall(PetscSectionGetOffset(globalCoarse, c, &globalOff));
2364: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff + 1) : globalOff, newOffsets, PETSC_FALSE, perm, NULL, pInd));
2365: }
2366: }
2367: for (f = 0; f < maxFields; f++) {
2368: if (numFields) PetscCall(PetscSectionRestoreFieldPointSyms(localCoarse, f, numPoints, points, &perms[f], &flips[f]));
2369: else PetscCall(PetscSectionRestorePointSyms(localCoarse, numPoints, points, &perms[f], &flips[f]));
2370: }
2371: PetscCall(DMRestoreWorkArray(coarse, numPoints, MPIU_SCALAR, &points));
2372: }
2373: } else if (matSize) {
2374: PetscInt cOff;
2375: PetscInt *rowIndices, *colIndices, a, aDof = 0, aOff;
2377: numRowIndices = dof;
2378: PetscCall(DMGetWorkArray(coarse, numRowIndices, MPIU_INT, &rowIndices));
2379: PetscCall(DMGetWorkArray(coarse, numColIndices, MPIU_INT, &colIndices));
2380: PetscCall(PetscSectionGetOffset(cSec, p, &cOff));
2381: PetscCall(PetscSectionGetDof(aSec, p, &aDof));
2382: PetscCall(PetscSectionGetOffset(aSec, p, &aOff));
2383: if (numFields) {
2384: for (f = 0; f < numFields; f++) {
2385: PetscInt fDof;
2387: PetscCall(PetscSectionGetFieldDof(cSec, p, f, &fDof));
2388: offsets[f + 1] = fDof;
2389: for (a = 0; a < aDof; a++) {
2390: PetscInt anchor = anchors[a + aOff];
2391: PetscCall(PetscSectionGetFieldDof(localCoarse, anchor, f, &fDof));
2392: newOffsets[f + 1] += fDof;
2393: }
2394: }
2395: for (f = 0; f < numFields; f++) {
2396: offsets[f + 1] += offsets[f];
2397: offsetsCopy[f + 1] = offsets[f + 1];
2398: newOffsets[f + 1] += newOffsets[f];
2399: newOffsetsCopy[f + 1] = newOffsets[f + 1];
2400: }
2401: PetscCall(DMPlexGetIndicesPointFields_Internal(cSec, PETSC_TRUE, p, cOff, offsetsCopy, PETSC_TRUE, NULL, -1, NULL, rowIndices));
2402: for (a = 0; a < aDof; a++) {
2403: PetscInt anchor = anchors[a + aOff], lOff;
2404: PetscCall(PetscSectionGetOffset(localCoarse, anchor, &lOff));
2405: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_TRUE, anchor, lOff, newOffsetsCopy, PETSC_TRUE, NULL, -1, NULL, colIndices));
2406: }
2407: } else {
2408: PetscCall(DMPlexGetIndicesPoint_Internal(cSec, PETSC_TRUE, p, cOff, offsetsCopy, PETSC_TRUE, NULL, NULL, rowIndices));
2409: for (a = 0; a < aDof; a++) {
2410: PetscInt anchor = anchors[a + aOff], lOff;
2411: PetscCall(PetscSectionGetOffset(localCoarse, anchor, &lOff));
2412: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_TRUE, anchor, lOff, newOffsetsCopy, PETSC_TRUE, NULL, NULL, colIndices));
2413: }
2414: }
2415: if (numFields) {
2416: PetscInt count, a;
2418: for (f = 0, count = 0; f < numFields; f++) {
2419: PetscInt iSize = offsets[f + 1] - offsets[f];
2420: PetscInt jSize = newOffsets[f + 1] - newOffsets[f];
2421: PetscCall(MatGetValues(cMat, iSize, &rowIndices[offsets[f]], jSize, &colIndices[newOffsets[f]], &pMat[count]));
2422: count += iSize * jSize;
2423: pInd[numColIndices + f] = offsets[f + 1];
2424: pInd[numColIndices + numFields + f] = newOffsets[f + 1];
2425: }
2426: for (a = 0; a < aDof; a++) {
2427: PetscInt anchor = anchors[a + aOff];
2428: PetscInt gOff;
2429: PetscCall(PetscSectionGetOffset(globalCoarse, anchor, &gOff));
2430: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, anchor, gOff < 0 ? -(gOff + 1) : gOff, newOffsets, PETSC_FALSE, NULL, -1, NULL, pInd));
2431: }
2432: } else {
2433: PetscCall(MatGetValues(cMat, numRowIndices, rowIndices, numColIndices, colIndices, pMat));
2434: for (PetscInt a = 0; a < aDof; a++) {
2435: PetscInt anchor = anchors[a + aOff];
2436: PetscInt gOff;
2437: PetscCall(PetscSectionGetOffset(globalCoarse, anchor, &gOff));
2438: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, anchor, gOff < 0 ? -(gOff + 1) : gOff, newOffsets, PETSC_FALSE, NULL, NULL, pInd));
2439: }
2440: }
2441: PetscCall(DMRestoreWorkArray(coarse, numColIndices, MPIU_INT, &colIndices));
2442: PetscCall(DMRestoreWorkArray(coarse, numRowIndices, MPIU_INT, &rowIndices));
2443: } else {
2444: PetscInt gOff;
2446: PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));
2447: if (numFields) {
2448: for (f = 0; f < numFields; f++) {
2449: PetscInt fDof;
2450: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
2451: offsets[f + 1] = fDof + offsets[f];
2452: }
2453: for (f = 0; f < numFields; f++) {
2454: pInd[numColIndices + f] = offsets[f + 1];
2455: pInd[numColIndices + numFields + f] = offsets[f + 1];
2456: }
2457: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, -1, NULL, pInd));
2458: } else {
2459: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, NULL, pInd));
2460: }
2461: }
2462: }
2463: PetscCall(PetscFree(maxChildIds));
2464: }
2465: {
2466: PetscSF indicesSF, matricesSF;
2467: PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices;
2469: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafIndicesSec));
2470: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafMatricesSec));
2471: PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootIndicesSec, &remoteOffsetsIndices, leafIndicesSec));
2472: PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootMatricesSec, &remoteOffsetsMatrices, leafMatricesSec));
2473: PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootIndicesSec, remoteOffsetsIndices, leafIndicesSec, &indicesSF));
2474: PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootMatricesSec, remoteOffsetsMatrices, leafMatricesSec, &matricesSF));
2475: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
2476: PetscCall(PetscFree(remoteOffsetsIndices));
2477: PetscCall(PetscFree(remoteOffsetsMatrices));
2478: PetscCall(PetscSectionGetStorageSize(leafIndicesSec, &numLeafIndices));
2479: PetscCall(PetscSectionGetStorageSize(leafMatricesSec, &numLeafMatrices));
2480: PetscCall(PetscMalloc2(numLeafIndices, &leafIndices, numLeafMatrices, &leafMatrices));
2481: PetscCall(PetscSFBcastBegin(indicesSF, MPIU_INT, rootIndices, leafIndices, MPI_REPLACE));
2482: PetscCall(PetscSFBcastBegin(matricesSF, MPIU_SCALAR, rootMatrices, leafMatrices, MPI_REPLACE));
2483: PetscCall(PetscSFBcastEnd(indicesSF, MPIU_INT, rootIndices, leafIndices, MPI_REPLACE));
2484: PetscCall(PetscSFBcastEnd(matricesSF, MPIU_SCALAR, rootMatrices, leafMatrices, MPI_REPLACE));
2485: PetscCall(PetscSFDestroy(&matricesSF));
2486: PetscCall(PetscSFDestroy(&indicesSF));
2487: PetscCall(PetscFree2(rootIndices, rootMatrices));
2488: PetscCall(PetscSectionDestroy(&rootIndicesSec));
2489: PetscCall(PetscSectionDestroy(&rootMatricesSec));
2490: }
2491: /* count to preallocate */
2492: PetscCall(DMGetLocalSection(fine, &localFine));
2493: {
2494: PetscInt nGlobal;
2495: PetscInt *dnnz, *onnz;
2496: PetscLayout rowMap, colMap;
2497: PetscInt rowStart, rowEnd, colStart, colEnd;
2498: PetscInt maxDof;
2499: PetscInt *rowIndices;
2500: DM refTree;
2501: PetscInt **refPointFieldN;
2502: PetscScalar ***refPointFieldMats;
2503: PetscSection refConSec, refAnSec;
2504: PetscInt pRefStart, pRefEnd, maxConDof, maxColumns, leafStart, leafEnd;
2505: PetscScalar *pointWork;
2507: PetscCall(PetscSectionGetConstrainedStorageSize(globalFine, &nGlobal));
2508: PetscCall(PetscCalloc2(nGlobal, &dnnz, nGlobal, &onnz));
2509: PetscCall(MatGetLayouts(mat, &rowMap, &colMap));
2510: PetscCall(PetscLayoutSetUp(rowMap));
2511: PetscCall(PetscLayoutSetUp(colMap));
2512: PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd));
2513: PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd));
2514: PetscCall(PetscSectionGetMaxDof(localFine, &maxDof));
2515: PetscCall(PetscSectionGetChart(leafIndicesSec, &leafStart, &leafEnd));
2516: PetscCall(DMGetWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
2517: for (p = leafStart; p < leafEnd; p++) {
2518: PetscInt gDof, gcDof, gOff;
2519: PetscInt numColIndices, pIndOff, *pInd;
2520: PetscInt matSize;
2521: PetscInt i;
2523: PetscCall(PetscSectionGetDof(globalFine, p, &gDof));
2524: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &gcDof));
2525: if ((gDof - gcDof) <= 0) continue;
2526: PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
2527: PetscCheck(gOff >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "I though having global dofs meant a non-negative offset");
2528: PetscCheck(gOff >= rowStart && (gOff + gDof - gcDof) <= rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "I thought the row map would constrain the global dofs");
2529: PetscCall(PetscSectionGetDof(leafIndicesSec, p, &numColIndices));
2530: PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &pIndOff));
2531: numColIndices -= 2 * numFields;
2532: PetscCheck(numColIndices > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "global fine dof with no dofs to interpolate from");
2533: pInd = &leafIndices[pIndOff];
2534: offsets[0] = 0;
2535: offsetsCopy[0] = 0;
2536: newOffsets[0] = 0;
2537: newOffsetsCopy[0] = 0;
2538: if (numFields) {
2539: PetscInt f;
2540: for (f = 0; f < numFields; f++) {
2541: PetscInt rowDof;
2543: PetscCall(PetscSectionGetFieldDof(localFine, p, f, &rowDof));
2544: offsets[f + 1] = offsets[f] + rowDof;
2545: offsetsCopy[f + 1] = offsets[f + 1];
2546: newOffsets[f + 1] = pInd[numColIndices + numFields + f];
2547: numD[f] = 0;
2548: numO[f] = 0;
2549: }
2550: PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices));
2551: for (f = 0; f < numFields; f++) {
2552: PetscInt colOffset = newOffsets[f];
2553: PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f];
2555: for (i = 0; i < numFieldCols; i++) {
2556: PetscInt gInd = pInd[i + colOffset];
2558: if (gInd >= colStart && gInd < colEnd) {
2559: numD[f]++;
2560: } else if (gInd >= 0) { /* negative means non-entry */
2561: numO[f]++;
2562: }
2563: }
2564: }
2565: } else {
2566: PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices));
2567: numD[0] = 0;
2568: numO[0] = 0;
2569: for (i = 0; i < numColIndices; i++) {
2570: PetscInt gInd = pInd[i];
2572: if (gInd >= colStart && gInd < colEnd) {
2573: numD[0]++;
2574: } else if (gInd >= 0) { /* negative means non-entry */
2575: numO[0]++;
2576: }
2577: }
2578: }
2579: PetscCall(PetscSectionGetDof(leafMatricesSec, p, &matSize));
2580: if (!matSize) { /* incoming matrix is identity */
2581: PetscInt childId;
2583: childId = childIds[p - pStartF];
2584: if (childId < 0) { /* no child interpolation: one nnz per */
2585: if (numFields) {
2586: PetscInt f;
2587: for (f = 0; f < numFields; f++) {
2588: PetscInt numRows = offsets[f + 1] - offsets[f], row;
2589: for (row = 0; row < numRows; row++) {
2590: PetscInt gIndCoarse = pInd[newOffsets[f] + row];
2591: PetscInt gIndFine = rowIndices[offsets[f] + row];
2592: if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2593: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2594: dnnz[gIndFine - rowStart] = 1;
2595: } else if (gIndCoarse >= 0) { /* remote */
2596: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2597: onnz[gIndFine - rowStart] = 1;
2598: } else { /* constrained */
2599: PetscCheck(gIndFine < 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2600: }
2601: }
2602: }
2603: } else {
2604: PetscInt i;
2605: for (i = 0; i < gDof; i++) {
2606: PetscInt gIndCoarse = pInd[i];
2607: PetscInt gIndFine = rowIndices[i];
2608: if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2609: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2610: dnnz[gIndFine - rowStart] = 1;
2611: } else if (gIndCoarse >= 0) { /* remote */
2612: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2613: onnz[gIndFine - rowStart] = 1;
2614: } else { /* constrained */
2615: PetscCheck(gIndFine < 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2616: }
2617: }
2618: }
2619: } else { /* interpolate from all */
2620: if (numFields) {
2621: for (PetscInt f = 0; f < numFields; f++) {
2622: PetscInt numRows = offsets[f + 1] - offsets[f], row;
2623: for (row = 0; row < numRows; row++) {
2624: PetscInt gIndFine = rowIndices[offsets[f] + row];
2625: if (gIndFine >= 0) {
2626: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2627: dnnz[gIndFine - rowStart] = numD[f];
2628: onnz[gIndFine - rowStart] = numO[f];
2629: }
2630: }
2631: }
2632: } else {
2633: for (PetscInt i = 0; i < gDof; i++) {
2634: PetscInt gIndFine = rowIndices[i];
2635: if (gIndFine >= 0) {
2636: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2637: dnnz[gIndFine - rowStart] = numD[0];
2638: onnz[gIndFine - rowStart] = numO[0];
2639: }
2640: }
2641: }
2642: }
2643: } else { /* interpolate from all */
2644: if (numFields) {
2645: for (PetscInt f = 0; f < numFields; f++) {
2646: PetscInt numRows = offsets[f + 1] - offsets[f], row;
2647: for (row = 0; row < numRows; row++) {
2648: PetscInt gIndFine = rowIndices[offsets[f] + row];
2649: if (gIndFine >= 0) {
2650: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2651: dnnz[gIndFine - rowStart] = numD[f];
2652: onnz[gIndFine - rowStart] = numO[f];
2653: }
2654: }
2655: }
2656: } else { /* every dof get a full row */
2657: for (PetscInt 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: for (PetscInt f = 0; f < numFields; f++) {
2733: PetscInt numRows = offsets[f + 1] - offsets[f];
2734: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2735: PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], refPointFieldMats[childId - pRefStart][f], INSERT_VALUES));
2736: }
2737: } else {
2738: PetscCall(MatSetValues(mat, gDof, rowIndices, numColIndices, pInd, refPointFieldMats[childId - pRefStart][0], INSERT_VALUES));
2739: }
2740: }
2741: } else { /* interpolate from all */
2742: PetscInt pMatOff;
2743: PetscScalar *pMat;
2745: PetscCall(PetscSectionGetOffset(leafMatricesSec, p, &pMatOff));
2746: pMat = &leafMatrices[pMatOff];
2747: if (childId < 0) { /* copy the incoming matrix */
2748: if (numFields) {
2749: PetscInt f, count;
2750: for (f = 0, count = 0; f < numFields; f++) {
2751: PetscInt numRows = offsets[f + 1] - offsets[f];
2752: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2753: PetscInt numInRows = rowOffsets[f + 1] - rowOffsets[f];
2754: PetscScalar *inMat = &pMat[count];
2756: PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], inMat, INSERT_VALUES));
2757: count += numCols * numInRows;
2758: }
2759: } else {
2760: PetscCall(MatSetValues(mat, gDof, rowIndices, numColIndices, pInd, pMat, INSERT_VALUES));
2761: }
2762: } else { /* multiply the incoming matrix by the child interpolation */
2763: if (numFields) {
2764: PetscInt f, count;
2765: for (f = 0, count = 0; f < numFields; f++) {
2766: PetscInt numRows = offsets[f + 1] - offsets[f];
2767: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2768: PetscInt numInRows = rowOffsets[f + 1] - rowOffsets[f];
2769: PetscScalar *inMat = &pMat[count];
2770: PetscInt i, j, k;
2771: PetscCheck(refPointFieldN[childId - pRefStart][f] == numInRows, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point constraint matrix multiply dimension mismatch");
2772: for (i = 0; i < numRows; i++) {
2773: for (j = 0; j < numCols; j++) {
2774: PetscScalar val = 0.;
2775: for (k = 0; k < numInRows; k++) val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j];
2776: pointWork[i * numCols + j] = val;
2777: }
2778: }
2779: PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], pointWork, INSERT_VALUES));
2780: count += numCols * numInRows;
2781: }
2782: } else { /* every dof gets a full row */
2783: PetscInt numRows = gDof;
2784: PetscInt numCols = numColIndices;
2785: PetscInt numInRows = matSize / numColIndices;
2786: PetscInt i, j, k;
2787: PetscCheck(refPointFieldN[childId - pRefStart][0] == numInRows, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point constraint matrix multiply dimension mismatch");
2788: for (i = 0; i < numRows; i++) {
2789: for (j = 0; j < numCols; j++) {
2790: PetscScalar val = 0.;
2791: for (k = 0; k < numInRows; k++) val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j];
2792: pointWork[i * numCols + j] = val;
2793: }
2794: }
2795: PetscCall(MatSetValues(mat, numRows, rowIndices, numCols, pInd, pointWork, INSERT_VALUES));
2796: }
2797: }
2798: }
2799: }
2800: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
2801: PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
2802: PetscCall(PetscFree(pointWork));
2803: }
2804: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
2805: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
2806: PetscCall(PetscSectionDestroy(&leafIndicesSec));
2807: PetscCall(PetscSectionDestroy(&leafMatricesSec));
2808: PetscCall(PetscFree2(leafIndices, leafMatrices));
2809: PetscCall(PetscFree2(*(PetscInt ****)&perms, *(PetscScalar ****)&flips));
2810: PetscCall(PetscFree7(offsets, offsetsCopy, newOffsets, newOffsetsCopy, rowOffsets, numD, numO));
2811: PetscCall(ISRestoreIndices(aIS, &anchors));
2812: PetscFunctionReturn(PETSC_SUCCESS);
2813: }
2815: /*
2816: * Assuming a nodal basis (w.r.t. the dual basis) basis:
2817: *
2818: * for each coarse dof \phi^c_i:
2819: * for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i:
2820: * for each fine dof \phi^f_j;
2821: * a_{i,j} = 0;
2822: * for each fine dof \phi^f_k:
2823: * a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l
2824: * [^^^ this is = \phi^c_i ^^^]
2825: */
2826: /*@
2827: DMPlexComputeInjectorReferenceTree - Compute the injection matrix from fine to coarse degrees of freedom on the reference tree
2829: Collective
2831: Input Parameter:
2832: . refTree - The reference-tree `DMPLEX` (see `DMPlexCreateDefaultReferenceTree()`)
2834: Output Parameter:
2835: . inj - The newly created injection `Mat` mapping fine-space coefficients on the reference tree to their coarse-space counterparts
2837: Level: developer
2839: Note:
2840: For a nodal basis, the injection is derived from the constraint matrix attached to the reference
2841: tree; the returned matrix is used internally by `DMPlexComputeInjectorTree()` to construct the
2842: global injection between refined and coarse meshes.
2844: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetReferenceTree()`, `DMPlexCreateDefaultReferenceTree()`, `DMPlexComputeInjectorTree()`, `DMPlexComputeInterpolatorTree()`
2845: @*/
2846: PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj)
2847: {
2848: PetscDS ds;
2849: PetscSection section, cSection;
2850: DMLabel canonical, depth;
2851: Mat cMat, mat;
2852: PetscInt *nnz;
2853: PetscInt f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd;
2854: PetscInt m, n;
2855: PetscScalar *pointScalar;
2856: PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent;
2858: PetscFunctionBegin;
2859: PetscCall(DMGetLocalSection(refTree, §ion));
2860: PetscCall(DMGetDimension(refTree, &dim));
2861: PetscCall(PetscMalloc6(dim, &v0, dim, &v0parent, dim, &vtmp, dim * dim, &J, dim * dim, &Jparent, dim * dim, &invJ));
2862: PetscCall(PetscMalloc2(dim, &pointScalar, dim, &pointRef));
2863: PetscCall(DMGetDS(refTree, &ds));
2864: PetscCall(PetscDSGetNumFields(ds, &numFields));
2865: PetscCall(PetscSectionGetNumFields(section, &numSecFields));
2866: PetscCall(DMGetLabel(refTree, "canonical", &canonical));
2867: PetscCall(DMGetLabel(refTree, "depth", &depth));
2868: PetscCall(DMGetDefaultConstraints(refTree, &cSection, &cMat, NULL));
2869: PetscCall(DMPlexGetChart(refTree, &pStart, &pEnd));
2870: PetscCall(DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd));
2871: PetscCall(MatGetSize(cMat, &n, &m)); /* the injector has transpose sizes from the constraint matrix */
2872: /* Step 1: compute non-zero pattern. A proper subset of constraint matrix non-zero */
2873: PetscCall(PetscCalloc1(m, &nnz));
2874: for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */
2875: const PetscInt *children;
2876: PetscInt numChildren;
2877: PetscInt i, numChildDof, numSelfDof;
2879: if (canonical) {
2880: PetscInt pCanonical;
2881: PetscCall(DMLabelGetValue(canonical, p, &pCanonical));
2882: if (p != pCanonical) continue;
2883: }
2884: PetscCall(DMPlexGetTreeChildren(refTree, p, &numChildren, &children));
2885: if (!numChildren) continue;
2886: for (i = 0, numChildDof = 0; i < numChildren; i++) {
2887: PetscInt child = children[i];
2888: PetscInt dof;
2890: PetscCall(PetscSectionGetDof(section, child, &dof));
2891: numChildDof += dof;
2892: }
2893: PetscCall(PetscSectionGetDof(section, p, &numSelfDof));
2894: if (!numChildDof || !numSelfDof) continue;
2895: for (f = 0; f < numFields; f++) {
2896: PetscInt selfOff;
2898: if (numSecFields) { /* count the dofs for just this field */
2899: for (i = 0, numChildDof = 0; i < numChildren; i++) {
2900: PetscInt child = children[i];
2901: PetscInt dof;
2903: PetscCall(PetscSectionGetFieldDof(section, child, f, &dof));
2904: numChildDof += dof;
2905: }
2906: PetscCall(PetscSectionGetFieldDof(section, p, f, &numSelfDof));
2907: PetscCall(PetscSectionGetFieldOffset(section, p, f, &selfOff));
2908: } else {
2909: PetscCall(PetscSectionGetOffset(section, p, &selfOff));
2910: }
2911: for (i = 0; i < numSelfDof; i++) nnz[selfOff + i] = numChildDof;
2912: }
2913: }
2914: PetscCall(MatCreateAIJ(PETSC_COMM_SELF, m, n, m, n, -1, nnz, -1, NULL, &mat));
2915: PetscCall(PetscFree(nnz));
2916: /* Setp 2: compute entries */
2917: for (p = pStart; p < pEnd; p++) {
2918: const PetscInt *children;
2919: PetscInt numChildren;
2920: PetscInt i, numChildDof, numSelfDof;
2922: /* same conditions about when entries occur */
2923: if (canonical) {
2924: PetscInt pCanonical;
2925: PetscCall(DMLabelGetValue(canonical, p, &pCanonical));
2926: if (p != pCanonical) continue;
2927: }
2928: PetscCall(DMPlexGetTreeChildren(refTree, p, &numChildren, &children));
2929: if (!numChildren) continue;
2930: for (i = 0, numChildDof = 0; i < numChildren; i++) {
2931: PetscInt child = children[i];
2932: PetscInt dof;
2934: PetscCall(PetscSectionGetDof(section, child, &dof));
2935: numChildDof += dof;
2936: }
2937: PetscCall(PetscSectionGetDof(section, p, &numSelfDof));
2938: if (!numChildDof || !numSelfDof) continue;
2940: for (f = 0; f < numFields; f++) {
2941: PetscInt pI = -1, cI = -1;
2942: PetscInt selfOff, Nc, parentCell;
2943: PetscInt cellShapeOff;
2944: PetscObject disc;
2945: PetscDualSpace dsp;
2946: PetscClassId classId;
2947: PetscScalar *pointMat;
2948: PetscInt *matRows, *matCols;
2949: PetscInt pO = PETSC_INT_MIN;
2950: const PetscInt *depthNumDof;
2952: if (numSecFields) {
2953: for (i = 0, numChildDof = 0; i < numChildren; i++) {
2954: PetscInt child = children[i];
2955: PetscInt dof;
2957: PetscCall(PetscSectionGetFieldDof(section, child, f, &dof));
2958: numChildDof += dof;
2959: }
2960: PetscCall(PetscSectionGetFieldDof(section, p, f, &numSelfDof));
2961: PetscCall(PetscSectionGetFieldOffset(section, p, f, &selfOff));
2962: } else {
2963: PetscCall(PetscSectionGetOffset(section, p, &selfOff));
2964: }
2966: /* find a cell whose closure contains p */
2967: if (p >= cStart && p < cEnd) {
2968: parentCell = p;
2969: } else {
2970: PetscInt *star = NULL;
2971: PetscInt numStar;
2973: parentCell = -1;
2974: PetscCall(DMPlexGetTransitiveClosure(refTree, p, PETSC_FALSE, &numStar, &star));
2975: for (i = numStar - 1; i >= 0; i--) {
2976: PetscInt c = star[2 * i];
2978: if (c >= cStart && c < cEnd) {
2979: parentCell = c;
2980: break;
2981: }
2982: }
2983: PetscCall(DMPlexRestoreTransitiveClosure(refTree, p, PETSC_FALSE, &numStar, &star));
2984: }
2985: /* determine the offset of p's shape functions within parentCell's shape functions */
2986: PetscCall(PetscDSGetDiscretization(ds, f, &disc));
2987: PetscCall(PetscObjectGetClassId(disc, &classId));
2988: if (classId == PETSCFE_CLASSID) PetscCall(PetscFEGetDualSpace((PetscFE)disc, &dsp));
2989: else {
2990: PetscCheck(classId == PETSCFV_CLASSID, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported discretization object");
2991: PetscCall(PetscFVGetDualSpace((PetscFV)disc, &dsp));
2992: }
2993: PetscCall(PetscDualSpaceGetNumDof(dsp, &depthNumDof));
2994: PetscCall(PetscDualSpaceGetNumComponents(dsp, &Nc));
2995: {
2996: PetscInt *closure = NULL;
2997: PetscInt numClosure;
2999: PetscCall(DMPlexGetTransitiveClosure(refTree, parentCell, PETSC_TRUE, &numClosure, &closure));
3000: for (i = 0, pI = -1, cellShapeOff = 0; i < numClosure; i++) {
3001: PetscInt point = closure[2 * i], pointDepth;
3003: pO = closure[2 * i + 1];
3004: if (point == p) {
3005: pI = i;
3006: break;
3007: }
3008: PetscCall(DMLabelGetValue(depth, point, &pointDepth));
3009: cellShapeOff += depthNumDof[pointDepth];
3010: }
3011: PetscCall(DMPlexRestoreTransitiveClosure(refTree, parentCell, PETSC_TRUE, &numClosure, &closure));
3012: }
3014: PetscCall(DMGetWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR, &pointMat));
3015: PetscCall(DMGetWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT, &matRows));
3016: matCols = matRows + numSelfDof;
3017: for (i = 0; i < numSelfDof; i++) matRows[i] = selfOff + i;
3018: for (i = 0; i < numSelfDof * numChildDof; i++) pointMat[i] = 0.;
3019: {
3020: PetscInt colOff = 0;
3022: for (i = 0; i < numChildren; i++) {
3023: PetscInt child = children[i];
3024: PetscInt dof, off, j;
3026: if (numSecFields) {
3027: PetscCall(PetscSectionGetFieldDof(cSection, child, f, &dof));
3028: PetscCall(PetscSectionGetFieldOffset(cSection, child, f, &off));
3029: } else {
3030: PetscCall(PetscSectionGetDof(cSection, child, &dof));
3031: PetscCall(PetscSectionGetOffset(cSection, child, &off));
3032: }
3034: for (j = 0; j < dof; j++) matCols[colOff++] = off + j;
3035: }
3036: }
3037: if (classId == PETSCFE_CLASSID) {
3038: PetscFE fe = (PetscFE)disc;
3039: PetscInt fSize;
3040: const PetscInt ***perms;
3041: const PetscScalar ***flips;
3042: const PetscInt *pperms;
3044: PetscCall(PetscFEGetDualSpace(fe, &dsp));
3045: PetscCall(PetscDualSpaceGetDimension(dsp, &fSize));
3046: PetscCall(PetscDualSpaceGetSymmetries(dsp, &perms, &flips));
3047: pperms = perms ? perms[pI] ? perms[pI][pO] : NULL : NULL;
3048: for (i = 0; i < numSelfDof; i++) { /* for every shape function */
3049: PetscQuadrature q;
3050: PetscInt dim, thisNc, numPoints, j, k;
3051: const PetscReal *points;
3052: const PetscReal *weights;
3053: PetscInt *closure = NULL;
3054: PetscInt numClosure;
3055: PetscInt iCell = pperms ? pperms[i] : i;
3056: PetscInt parentCellShapeDof = cellShapeOff + iCell;
3057: PetscTabulation Tparent;
3059: PetscCall(PetscDualSpaceGetFunctional(dsp, parentCellShapeDof, &q));
3060: PetscCall(PetscQuadratureGetData(q, &dim, &thisNc, &numPoints, &points, &weights));
3061: PetscCheck(thisNc == Nc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Functional dim %" PetscInt_FMT " does not much basis dim %" PetscInt_FMT, thisNc, Nc);
3062: PetscCall(PetscFECreateTabulation(fe, 1, numPoints, points, 0, &Tparent)); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */
3063: for (j = 0; j < numPoints; j++) {
3064: PetscInt childCell = -1;
3065: PetscReal *parentValAtPoint;
3066: const PetscReal xi0[3] = {-1., -1., -1.};
3067: const PetscReal *pointReal = &points[dim * j];
3068: const PetscScalar *point;
3069: PetscTabulation Tchild;
3070: PetscInt childCellShapeOff, pointMatOff;
3071: #if PetscDefined(USE_COMPLEX)
3072: for (PetscInt d = 0; d < dim; d++) pointScalar[d] = points[dim * j + d];
3073: point = pointScalar;
3074: #else
3075: point = pointReal;
3076: #endif
3078: parentValAtPoint = &Tparent->T[0][(fSize * j + parentCellShapeDof) * Nc];
3080: for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/
3081: PetscInt child = children[k];
3082: PetscInt *star = NULL;
3083: PetscInt numStar;
3085: PetscCall(DMPlexGetTransitiveClosure(refTree, child, PETSC_FALSE, &numStar, &star));
3086: for (PetscInt s = numStar - 1; s >= 0; s--) {
3087: PetscInt c = star[2 * s];
3089: if (c < cStart || c >= cEnd) continue;
3090: PetscCall(DMPlexLocatePoint_Internal(refTree, dim, point, c, &childCell));
3091: if (childCell >= 0) break;
3092: }
3093: PetscCall(DMPlexRestoreTransitiveClosure(refTree, child, PETSC_FALSE, &numStar, &star));
3094: if (childCell >= 0) break;
3095: }
3096: PetscCheck(childCell >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not locate quadrature point");
3097: PetscCall(DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ));
3098: PetscCall(DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent));
3099: CoordinatesRefToReal(dim, dim, xi0, v0parent, Jparent, pointReal, vtmp);
3100: CoordinatesRealToRef(dim, dim, xi0, v0, invJ, vtmp, pointRef);
3102: PetscCall(PetscFECreateTabulation(fe, 1, 1, pointRef, 0, &Tchild));
3103: PetscCall(DMPlexGetTransitiveClosure(refTree, childCell, PETSC_TRUE, &numClosure, &closure));
3104: for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */
3105: PetscInt child = children[k], childDepth, childDof, childO = PETSC_INT_MIN;
3106: PetscInt l;
3107: const PetscInt *cperms;
3109: PetscCall(DMLabelGetValue(depth, child, &childDepth));
3110: childDof = depthNumDof[childDepth];
3111: for (l = 0, cI = -1, childCellShapeOff = 0; l < numClosure; l++) {
3112: PetscInt point = closure[2 * l];
3113: PetscInt pointDepth;
3115: childO = closure[2 * l + 1];
3116: if (point == child) {
3117: cI = l;
3118: break;
3119: }
3120: PetscCall(DMLabelGetValue(depth, point, &pointDepth));
3121: childCellShapeOff += depthNumDof[pointDepth];
3122: }
3123: if (l == numClosure) {
3124: pointMatOff += childDof;
3125: continue; /* child is not in the closure of the cell: has nothing to contribute to this point */
3126: }
3127: cperms = perms ? perms[cI] ? perms[cI][childO] : NULL : NULL;
3128: for (l = 0; l < childDof; l++) {
3129: PetscInt lCell = cperms ? cperms[l] : l;
3130: PetscInt childCellDof = childCellShapeOff + lCell;
3131: PetscReal *childValAtPoint;
3132: PetscReal val = 0.;
3134: childValAtPoint = &Tchild->T[0][childCellDof * Nc];
3135: for (m = 0; m < Nc; m++) val += weights[j * Nc + m] * parentValAtPoint[m] * childValAtPoint[m];
3137: pointMat[i * numChildDof + pointMatOff + l] += val;
3138: }
3139: pointMatOff += childDof;
3140: }
3141: PetscCall(DMPlexRestoreTransitiveClosure(refTree, childCell, PETSC_TRUE, &numClosure, &closure));
3142: PetscCall(PetscTabulationDestroy(&Tchild));
3143: }
3144: PetscCall(PetscTabulationDestroy(&Tparent));
3145: }
3146: } else { /* just the volume-weighted averages of the children */
3147: PetscReal parentVol;
3148: PetscInt childCell;
3150: PetscCall(DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL));
3151: for (i = 0, childCell = 0; i < numChildren; i++) {
3152: PetscInt child = children[i], j;
3153: PetscReal childVol;
3155: if (child < cStart || child >= cEnd) continue;
3156: PetscCall(DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL));
3157: for (j = 0; j < Nc; j++) pointMat[j * numChildDof + Nc * childCell + j] = childVol / parentVol;
3158: childCell++;
3159: }
3160: }
3161: /* Insert pointMat into mat */
3162: PetscCall(MatSetValues(mat, numSelfDof, matRows, numChildDof, matCols, pointMat, INSERT_VALUES));
3163: PetscCall(DMRestoreWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT, &matRows));
3164: PetscCall(DMRestoreWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR, &pointMat));
3165: }
3166: }
3167: PetscCall(PetscFree6(v0, v0parent, vtmp, J, Jparent, invJ));
3168: PetscCall(PetscFree2(pointScalar, pointRef));
3169: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
3170: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
3171: *inj = mat;
3172: PetscFunctionReturn(PETSC_SUCCESS);
3173: }
3175: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3176: {
3177: PetscDS ds;
3178: PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof;
3179: PetscScalar ***refPointFieldMats;
3180: PetscSection refConSec, refSection;
3182: PetscFunctionBegin;
3183: PetscCall(DMGetDS(refTree, &ds));
3184: PetscCall(PetscDSGetNumFields(ds, &numFields));
3185: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
3186: PetscCall(DMGetLocalSection(refTree, &refSection));
3187: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
3188: PetscCall(PetscMalloc1(pRefEnd - pRefStart, &refPointFieldMats));
3189: PetscCall(PetscSectionGetMaxDof(refConSec, &maxDof));
3190: PetscCall(PetscMalloc1(maxDof, &rows));
3191: PetscCall(PetscMalloc1(maxDof * maxDof, &cols));
3192: for (p = pRefStart; p < pRefEnd; p++) {
3193: PetscInt parent, pDof, parentDof;
3195: PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL));
3196: PetscCall(PetscSectionGetDof(refConSec, p, &pDof));
3197: PetscCall(PetscSectionGetDof(refSection, parent, &parentDof));
3198: if (!pDof || !parentDof || parent == p) continue;
3200: PetscCall(PetscMalloc1(numFields, &refPointFieldMats[p - pRefStart]));
3201: for (f = 0; f < numFields; f++) {
3202: PetscInt cDof, cOff, numCols, r;
3204: if (numFields > 1) {
3205: PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof));
3206: PetscCall(PetscSectionGetFieldOffset(refConSec, p, f, &cOff));
3207: } else {
3208: PetscCall(PetscSectionGetDof(refConSec, p, &cDof));
3209: PetscCall(PetscSectionGetOffset(refConSec, p, &cOff));
3210: }
3212: for (r = 0; r < cDof; r++) rows[r] = cOff + r;
3213: numCols = 0;
3214: {
3215: PetscInt aDof, aOff, j;
3217: if (numFields > 1) {
3218: PetscCall(PetscSectionGetFieldDof(refSection, parent, f, &aDof));
3219: PetscCall(PetscSectionGetFieldOffset(refSection, parent, f, &aOff));
3220: } else {
3221: PetscCall(PetscSectionGetDof(refSection, parent, &aDof));
3222: PetscCall(PetscSectionGetOffset(refSection, parent, &aOff));
3223: }
3225: for (j = 0; j < aDof; j++) cols[numCols++] = aOff + j;
3226: }
3227: PetscCall(PetscMalloc1(cDof * numCols, &refPointFieldMats[p - pRefStart][f]));
3228: /* transpose of constraint matrix */
3229: PetscCall(MatGetValues(inj, numCols, cols, cDof, rows, refPointFieldMats[p - pRefStart][f]));
3230: }
3231: }
3232: *childrenMats = refPointFieldMats;
3233: PetscCall(PetscFree(rows));
3234: PetscCall(PetscFree(cols));
3235: PetscFunctionReturn(PETSC_SUCCESS);
3236: }
3238: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3239: {
3240: PetscDS ds;
3241: PetscScalar ***refPointFieldMats;
3242: PetscInt numFields, pRefStart, pRefEnd, p, f;
3243: PetscSection refConSec, refSection;
3245: PetscFunctionBegin;
3246: refPointFieldMats = *childrenMats;
3247: *childrenMats = NULL;
3248: PetscCall(DMGetDS(refTree, &ds));
3249: PetscCall(DMGetLocalSection(refTree, &refSection));
3250: PetscCall(PetscDSGetNumFields(ds, &numFields));
3251: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
3252: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
3253: for (p = pRefStart; p < pRefEnd; p++) {
3254: PetscInt parent, pDof, parentDof;
3256: PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL));
3257: PetscCall(PetscSectionGetDof(refConSec, p, &pDof));
3258: PetscCall(PetscSectionGetDof(refSection, parent, &parentDof));
3259: if (!pDof || !parentDof || parent == p) continue;
3261: for (f = 0; f < numFields; f++) {
3262: PetscInt cDof;
3264: if (numFields > 1) {
3265: PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof));
3266: } else {
3267: PetscCall(PetscSectionGetDof(refConSec, p, &cDof));
3268: }
3270: PetscCall(PetscFree(refPointFieldMats[p - pRefStart][f]));
3271: }
3272: PetscCall(PetscFree(refPointFieldMats[p - pRefStart]));
3273: }
3274: PetscCall(PetscFree(refPointFieldMats));
3275: PetscFunctionReturn(PETSC_SUCCESS);
3276: }
3278: static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree, Mat *injRef)
3279: {
3280: Mat cMatRef;
3281: PetscObject injRefObj;
3283: PetscFunctionBegin;
3284: PetscCall(DMGetDefaultConstraints(refTree, NULL, &cMatRef, NULL));
3285: PetscCall(PetscObjectQuery((PetscObject)cMatRef, "DMPlexComputeInjectorTree_refTree", &injRefObj));
3286: *injRef = (Mat)injRefObj;
3287: if (!*injRef) {
3288: PetscCall(DMPlexComputeInjectorReferenceTree(refTree, injRef));
3289: PetscCall(PetscObjectCompose((PetscObject)cMatRef, "DMPlexComputeInjectorTree_refTree", (PetscObject)*injRef));
3290: /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */
3291: PetscCall(PetscObjectDereference((PetscObject)*injRef));
3292: }
3293: PetscFunctionReturn(PETSC_SUCCESS);
3294: }
3296: 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)
3297: {
3298: PetscInt pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti;
3299: PetscSection globalCoarse, globalFine;
3300: PetscSection localCoarse, localFine, leafIndicesSec;
3301: PetscSection multiRootSec, rootIndicesSec;
3302: PetscInt *leafInds, *rootInds = NULL;
3303: const PetscInt *rootDegrees;
3304: PetscScalar *leafVals = NULL, *rootVals = NULL;
3305: PetscSF coarseToFineEmbedded;
3307: PetscFunctionBegin;
3308: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
3309: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
3310: PetscCall(DMGetLocalSection(fine, &localFine));
3311: PetscCall(DMGetGlobalSection(fine, &globalFine));
3312: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafIndicesSec));
3313: PetscCall(PetscSectionSetChart(leafIndicesSec, pStartF, pEndF));
3314: PetscCall(PetscSectionGetMaxDof(localFine, &maxDof));
3315: { /* winnow fine points that don't have global dofs out of the sf */
3316: PetscInt l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices;
3317: const PetscInt *leaves;
3319: PetscCall(PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL));
3320: for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
3321: p = leaves ? leaves[l] : l;
3322: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
3323: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
3324: if ((dof - cdof) > 0) {
3325: numPointsWithDofs++;
3327: PetscCall(PetscSectionGetDof(localFine, p, &dof));
3328: PetscCall(PetscSectionSetDof(leafIndicesSec, p, dof + 1));
3329: }
3330: }
3331: PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
3332: PetscCall(PetscSectionSetUp(leafIndicesSec));
3333: PetscCall(PetscSectionGetStorageSize(leafIndicesSec, &numIndices));
3334: PetscCall(PetscMalloc1((gatheredIndices ? numIndices : (maxDof + 1)), &leafInds));
3335: if (gatheredValues) PetscCall(PetscMalloc1(numIndices, &leafVals));
3336: for (l = 0, offset = 0; l < nleaves; l++) {
3337: p = leaves ? leaves[l] : l;
3338: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
3339: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
3340: if ((dof - cdof) > 0) {
3341: PetscInt off, gOff;
3342: PetscInt *pInd;
3343: PetscScalar *pVal = NULL;
3345: pointsWithDofs[offset++] = l;
3347: PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &off));
3349: pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds;
3350: if (gatheredValues) {
3351: pVal = &leafVals[off + 1];
3352: for (PetscInt i = 0; i < dof; i++) pVal[i] = 0.;
3353: }
3354: PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
3356: offsets[0] = 0;
3357: if (numFields) {
3358: for (PetscInt f = 0; f < numFields; f++) {
3359: PetscInt fDof;
3360: PetscCall(PetscSectionGetFieldDof(localFine, p, f, &fDof));
3361: offsets[f + 1] = fDof + offsets[f];
3362: }
3363: PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, -1, NULL, pInd));
3364: } else {
3365: PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, NULL, pInd));
3366: }
3367: if (gatheredValues) PetscCall(VecGetValues(fineVec, dof, pInd, pVal));
3368: }
3369: }
3370: PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded));
3371: PetscCall(PetscFree(pointsWithDofs));
3372: }
3374: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
3375: PetscCall(DMGetLocalSection(coarse, &localCoarse));
3376: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
3378: { /* there may be the case where an sf root has a parent: broadcast parents back to children */
3379: MPI_Datatype threeInt;
3380: PetscMPIInt rank;
3381: PetscInt (*parentNodeAndIdCoarse)[3];
3382: PetscInt (*parentNodeAndIdFine)[3];
3383: PetscInt p, nleaves, nleavesToParents;
3384: PetscSF pointSF, sfToParents;
3385: const PetscInt *ilocal;
3386: const PetscSFNode *iremote;
3387: PetscSFNode *iremoteToParents;
3388: PetscInt *ilocalToParents;
3390: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank));
3391: PetscCallMPI(MPI_Type_contiguous(3, MPIU_INT, &threeInt));
3392: PetscCallMPI(MPI_Type_commit(&threeInt));
3393: PetscCall(PetscMalloc2(pEndC - pStartC, &parentNodeAndIdCoarse, pEndF - pStartF, &parentNodeAndIdFine));
3394: PetscCall(DMGetPointSF(coarse, &pointSF));
3395: PetscCall(PetscSFGetGraph(pointSF, NULL, &nleaves, &ilocal, &iremote));
3396: for (p = pStartC; p < pEndC; p++) {
3397: PetscInt parent, childId;
3398: PetscCall(DMPlexGetTreeParent(coarse, p, &parent, &childId));
3399: parentNodeAndIdCoarse[p - pStartC][0] = rank;
3400: parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC;
3401: parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId;
3402: if (nleaves > 0) {
3403: PetscInt leaf = -1;
3405: if (ilocal) {
3406: PetscCall(PetscFindInt(parent, nleaves, ilocal, &leaf));
3407: } else {
3408: leaf = p - pStartC;
3409: }
3410: if (leaf >= 0) {
3411: parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank;
3412: parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index;
3413: }
3414: }
3415: }
3416: for (p = pStartF; p < pEndF; p++) {
3417: parentNodeAndIdFine[p - pStartF][0] = -1;
3418: parentNodeAndIdFine[p - pStartF][1] = -1;
3419: parentNodeAndIdFine[p - pStartF][2] = -1;
3420: }
3421: PetscCall(PetscSFBcastBegin(coarseToFineEmbedded, threeInt, parentNodeAndIdCoarse, parentNodeAndIdFine, MPI_REPLACE));
3422: PetscCall(PetscSFBcastEnd(coarseToFineEmbedded, threeInt, parentNodeAndIdCoarse, parentNodeAndIdFine, MPI_REPLACE));
3423: for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3424: PetscInt dof;
3426: PetscCall(PetscSectionGetDof(leafIndicesSec, p, &dof));
3427: if (dof) {
3428: PetscInt off;
3430: PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &off));
3431: if (gatheredIndices) {
3432: leafInds[off] = PetscMax(childIds[p - pStartF], parentNodeAndIdFine[p - pStartF][2]);
3433: } else if (gatheredValues) {
3434: leafVals[off] = (PetscScalar)PetscMax(childIds[p - pStartF], parentNodeAndIdFine[p - pStartF][2]);
3435: }
3436: }
3437: if (parentNodeAndIdFine[p - pStartF][0] >= 0) nleavesToParents++;
3438: }
3439: PetscCall(PetscMalloc1(nleavesToParents, &ilocalToParents));
3440: PetscCall(PetscMalloc1(nleavesToParents, &iremoteToParents));
3441: for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3442: if (parentNodeAndIdFine[p - pStartF][0] >= 0) {
3443: ilocalToParents[nleavesToParents] = p - pStartF;
3444: // FIXME PetscCall(PetscMPIIntCast(parentNodeAndIdFine[p - pStartF][0],&iremoteToParents[nleavesToParents].rank));
3445: iremoteToParents[nleavesToParents].rank = parentNodeAndIdFine[p - pStartF][0];
3446: iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p - pStartF][1];
3447: nleavesToParents++;
3448: }
3449: }
3450: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)coarse), &sfToParents));
3451: PetscCall(PetscSFSetGraph(sfToParents, pEndC - pStartC, nleavesToParents, ilocalToParents, PETSC_OWN_POINTER, iremoteToParents, PETSC_OWN_POINTER));
3452: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
3454: coarseToFineEmbedded = sfToParents;
3456: PetscCall(PetscFree2(parentNodeAndIdCoarse, parentNodeAndIdFine));
3457: PetscCallMPI(MPI_Type_free(&threeInt));
3458: }
3460: { /* winnow out coarse points that don't have dofs */
3461: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;
3462: PetscSF sfDofsOnly;
3464: for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) {
3465: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3466: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
3467: if ((dof - cdof) > 0) numPointsWithDofs++;
3468: }
3469: PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
3470: for (p = pStartC, offset = 0; p < pEndC; p++) {
3471: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3472: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
3473: if ((dof - cdof) > 0) pointsWithDofs[offset++] = p - pStartC;
3474: }
3475: PetscCall(PetscSFCreateEmbeddedRootSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly));
3476: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
3477: PetscCall(PetscFree(pointsWithDofs));
3478: coarseToFineEmbedded = sfDofsOnly;
3479: }
3481: /* communicate back to the coarse mesh which coarse points have children (that may require injection) */
3482: PetscCall(PetscSFComputeDegreeBegin(coarseToFineEmbedded, &rootDegrees));
3483: PetscCall(PetscSFComputeDegreeEnd(coarseToFineEmbedded, &rootDegrees));
3484: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &multiRootSec));
3485: PetscCall(PetscSectionSetChart(multiRootSec, pStartC, pEndC));
3486: for (p = pStartC; p < pEndC; p++) PetscCall(PetscSectionSetDof(multiRootSec, p, rootDegrees[p - pStartC]));
3487: PetscCall(PetscSectionSetUp(multiRootSec));
3488: PetscCall(PetscSectionGetStorageSize(multiRootSec, &numMulti));
3489: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootIndicesSec));
3490: { /* distribute the leaf section */
3491: PetscSF multi, multiInv, indicesSF;
3492: PetscInt *remoteOffsets, numRootIndices;
3494: PetscCall(PetscSFGetMultiSF(coarseToFineEmbedded, &multi));
3495: PetscCall(PetscSFCreateInverseSF(multi, &multiInv));
3496: PetscCall(PetscSFDistributeSection(multiInv, leafIndicesSec, &remoteOffsets, rootIndicesSec));
3497: PetscCall(PetscSFCreateSectionSF(multiInv, leafIndicesSec, remoteOffsets, rootIndicesSec, &indicesSF));
3498: PetscCall(PetscFree(remoteOffsets));
3499: PetscCall(PetscSFDestroy(&multiInv));
3500: PetscCall(PetscSectionGetStorageSize(rootIndicesSec, &numRootIndices));
3501: if (gatheredIndices) {
3502: PetscCall(PetscMalloc1(numRootIndices, &rootInds));
3503: PetscCall(PetscSFBcastBegin(indicesSF, MPIU_INT, leafInds, rootInds, MPI_REPLACE));
3504: PetscCall(PetscSFBcastEnd(indicesSF, MPIU_INT, leafInds, rootInds, MPI_REPLACE));
3505: }
3506: if (gatheredValues) {
3507: PetscCall(PetscMalloc1(numRootIndices, &rootVals));
3508: PetscCall(PetscSFBcastBegin(indicesSF, MPIU_SCALAR, leafVals, rootVals, MPI_REPLACE));
3509: PetscCall(PetscSFBcastEnd(indicesSF, MPIU_SCALAR, leafVals, rootVals, MPI_REPLACE));
3510: }
3511: PetscCall(PetscSFDestroy(&indicesSF));
3512: }
3513: PetscCall(PetscSectionDestroy(&leafIndicesSec));
3514: PetscCall(PetscFree(leafInds));
3515: PetscCall(PetscFree(leafVals));
3516: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
3517: *rootMultiSec = multiRootSec;
3518: *multiLeafSec = rootIndicesSec;
3519: if (gatheredIndices) *gatheredIndices = rootInds;
3520: if (gatheredValues) *gatheredValues = rootVals;
3521: PetscFunctionReturn(PETSC_SUCCESS);
3522: }
3524: PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
3525: {
3526: DM refTree;
3527: PetscSection multiRootSec, rootIndicesSec;
3528: PetscSection globalCoarse, globalFine;
3529: PetscSection localCoarse, localFine;
3530: PetscSection cSecRef;
3531: PetscInt *rootIndices = NULL, *parentIndices, pRefStart, pRefEnd;
3532: Mat injRef;
3533: PetscInt numFields, maxDof;
3534: PetscInt pStartC, pEndC, pStartF, pEndF, p;
3535: PetscInt *offsets, *offsetsCopy, *rowOffsets;
3536: PetscLayout rowMap, colMap;
3537: PetscInt rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO;
3538: PetscScalar ***childrenMats = NULL; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */
3540: PetscFunctionBegin;
3541: /* get the templates for the fine-to-coarse injection from the reference tree */
3542: PetscCall(DMPlexGetReferenceTree(coarse, &refTree));
3543: PetscCall(DMCopyDisc(coarse, refTree));
3544: PetscCall(DMSetLocalSection(refTree, NULL));
3545: PetscCall(DMSetDefaultConstraints(refTree, NULL, NULL, NULL));
3546: PetscCall(DMGetDefaultConstraints(refTree, &cSecRef, NULL, NULL));
3547: PetscCall(PetscSectionGetChart(cSecRef, &pRefStart, &pRefEnd));
3548: PetscCall(DMPlexReferenceTreeGetInjector(refTree, &injRef));
3550: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
3551: PetscCall(DMGetLocalSection(fine, &localFine));
3552: PetscCall(DMGetGlobalSection(fine, &globalFine));
3553: PetscCall(PetscSectionGetNumFields(localFine, &numFields));
3554: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
3555: PetscCall(DMGetLocalSection(coarse, &localCoarse));
3556: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
3557: PetscCall(PetscSectionGetMaxDof(localCoarse, &maxDof));
3558: {
3559: PetscInt maxFields = PetscMax(1, numFields) + 1;
3560: PetscCall(PetscMalloc3(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &rowOffsets));
3561: }
3563: PetscCall(DMPlexTransferInjectorTree(coarse, fine, coarseToFine, childIds, NULL, numFields, offsets, &multiRootSec, &rootIndicesSec, &rootIndices, NULL));
3565: PetscCall(PetscMalloc1(maxDof, &parentIndices));
3567: /* count indices */
3568: PetscCall(MatGetLayouts(mat, &rowMap, &colMap));
3569: PetscCall(PetscLayoutSetUp(rowMap));
3570: PetscCall(PetscLayoutSetUp(colMap));
3571: PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd));
3572: PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd));
3573: PetscCall(PetscCalloc2(rowEnd - rowStart, &nnzD, rowEnd - rowStart, &nnzO));
3574: for (p = pStartC; p < pEndC; p++) {
3575: PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
3577: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3578: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
3579: if ((dof - cdof) <= 0) continue;
3580: PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));
3582: rowOffsets[0] = 0;
3583: offsetsCopy[0] = 0;
3584: if (numFields) {
3585: PetscInt f;
3587: for (f = 0; f < numFields; f++) {
3588: PetscInt fDof;
3589: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
3590: rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3591: }
3592: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices));
3593: } else {
3594: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices));
3595: rowOffsets[1] = offsetsCopy[0];
3596: }
3598: PetscCall(PetscSectionGetDof(multiRootSec, p, &numLeaves));
3599: PetscCall(PetscSectionGetOffset(multiRootSec, p, &leafStart));
3600: leafEnd = leafStart + numLeaves;
3601: for (l = leafStart; l < leafEnd; l++) {
3602: PetscInt numIndices, childId, offset;
3603: const PetscInt *childIndices;
3605: PetscCall(PetscSectionGetDof(rootIndicesSec, l, &numIndices));
3606: PetscCall(PetscSectionGetOffset(rootIndicesSec, l, &offset));
3607: childId = rootIndices[offset++];
3608: childIndices = &rootIndices[offset];
3609: numIndices--;
3611: if (childId == -1) { /* equivalent points: scatter */
3612: PetscInt i;
3614: for (i = 0; i < numIndices; i++) {
3615: PetscInt colIndex = childIndices[i];
3616: PetscInt rowIndex = parentIndices[i];
3617: if (rowIndex < 0) continue;
3618: PetscCheck(colIndex >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unconstrained fine and constrained coarse");
3619: if (colIndex >= colStart && colIndex < colEnd) {
3620: nnzD[rowIndex - rowStart] = 1;
3621: } else {
3622: nnzO[rowIndex - rowStart] = 1;
3623: }
3624: }
3625: } else {
3626: PetscInt parentId, lim;
3628: PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL));
3630: lim = PetscMax(1, numFields);
3631: offsets[0] = 0;
3632: if (numFields) {
3633: for (PetscInt f = 0; f < numFields; f++) {
3634: PetscInt fDof;
3635: PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof));
3637: offsets[f + 1] = fDof + offsets[f];
3638: }
3639: } else {
3640: PetscInt cDof;
3642: PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof));
3643: offsets[1] = cDof;
3644: }
3645: for (PetscInt f = 0; f < lim; f++) {
3646: PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1];
3647: PetscInt childStart = offsets[f], childEnd = offsets[f + 1];
3648: PetscInt i, numD = 0, numO = 0;
3650: for (i = childStart; i < childEnd; i++) {
3651: PetscInt colIndex = childIndices[i];
3653: if (colIndex < 0) continue;
3654: if (colIndex >= colStart && colIndex < colEnd) {
3655: numD++;
3656: } else {
3657: numO++;
3658: }
3659: }
3660: for (i = parentStart; i < parentEnd; i++) {
3661: PetscInt rowIndex = parentIndices[i];
3663: if (rowIndex < 0) continue;
3664: nnzD[rowIndex - rowStart] += numD;
3665: nnzO[rowIndex - rowStart] += numO;
3666: }
3667: }
3668: }
3669: }
3670: }
3671: /* preallocate */
3672: PetscCall(MatXAIJSetPreallocation(mat, 1, nnzD, nnzO, NULL, NULL));
3673: PetscCall(PetscFree2(nnzD, nnzO));
3674: /* insert values */
3675: PetscCall(DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree, injRef, &childrenMats));
3676: for (p = pStartC; p < pEndC; p++) {
3677: PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
3679: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3680: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
3681: if ((dof - cdof) <= 0) continue;
3682: PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));
3684: rowOffsets[0] = 0;
3685: offsetsCopy[0] = 0;
3686: if (numFields) {
3687: PetscInt f;
3689: for (f = 0; f < numFields; f++) {
3690: PetscInt fDof;
3691: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
3692: rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3693: }
3694: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices));
3695: } else {
3696: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices));
3697: rowOffsets[1] = offsetsCopy[0];
3698: }
3700: PetscCall(PetscSectionGetDof(multiRootSec, p, &numLeaves));
3701: PetscCall(PetscSectionGetOffset(multiRootSec, p, &leafStart));
3702: leafEnd = leafStart + numLeaves;
3703: for (l = leafStart; l < leafEnd; l++) {
3704: PetscInt numIndices, childId, offset;
3705: const PetscInt *childIndices;
3707: PetscCall(PetscSectionGetDof(rootIndicesSec, l, &numIndices));
3708: PetscCall(PetscSectionGetOffset(rootIndicesSec, l, &offset));
3709: childId = rootIndices[offset++];
3710: childIndices = &rootIndices[offset];
3711: numIndices--;
3713: if (childId == -1) { /* equivalent points: scatter */
3714: for (PetscInt i = 0; i < numIndices; i++) PetscCall(MatSetValue(mat, parentIndices[i], childIndices[i], 1., INSERT_VALUES));
3715: } else {
3716: PetscInt parentId, lim;
3718: PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL));
3720: lim = PetscMax(1, numFields);
3721: offsets[0] = 0;
3722: if (numFields) {
3723: for (PetscInt f = 0; f < numFields; f++) {
3724: PetscInt fDof;
3725: PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof));
3727: offsets[f + 1] = fDof + offsets[f];
3728: }
3729: } else {
3730: PetscInt cDof;
3732: PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof));
3733: offsets[1] = cDof;
3734: }
3735: for (PetscInt f = 0; f < lim; f++) {
3736: PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0];
3737: PetscInt *rowIndices = &parentIndices[rowOffsets[f]];
3738: const PetscInt *colIndices = &childIndices[offsets[f]];
3740: PetscCall(MatSetValues(mat, rowOffsets[f + 1] - rowOffsets[f], rowIndices, offsets[f + 1] - offsets[f], colIndices, childMat, INSERT_VALUES));
3741: }
3742: }
3743: }
3744: }
3745: PetscCall(PetscSectionDestroy(&multiRootSec));
3746: PetscCall(PetscSectionDestroy(&rootIndicesSec));
3747: PetscCall(PetscFree(parentIndices));
3748: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree, injRef, &childrenMats));
3749: PetscCall(PetscFree(rootIndices));
3750: PetscCall(PetscFree3(offsets, offsetsCopy, rowOffsets));
3752: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
3753: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
3754: PetscFunctionReturn(PETSC_SUCCESS);
3755: }
3757: static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom)
3758: {
3759: PetscSF coarseToFineEmbedded;
3760: PetscSection globalCoarse, globalFine;
3761: PetscSection localCoarse, localFine;
3762: PetscSection aSec, cSec;
3763: PetscSection rootValuesSec;
3764: PetscSection leafValuesSec;
3765: PetscScalar *rootValues, *leafValues;
3766: IS aIS;
3767: const PetscInt *anchors;
3768: Mat cMat;
3769: PetscInt numFields;
3770: PetscInt pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd;
3771: PetscInt aStart, aEnd, cStart, cEnd;
3772: PetscInt *maxChildIds;
3773: PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
3774: PetscFV fv = NULL;
3775: PetscInt dim, numFVcomps = -1, fvField = -1;
3776: DM cellDM = NULL, gradDM = NULL;
3777: const PetscScalar *cellGeomArray = NULL;
3778: const PetscScalar *gradArray = NULL;
3780: PetscFunctionBegin;
3781: PetscCall(VecSetOption(vecFine, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE));
3782: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
3783: PetscCall(DMPlexGetSimplexOrBoxCells(coarse, 0, &cellStart, &cellEnd));
3784: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
3785: PetscCall(DMGetGlobalSection(fine, &globalFine));
3786: PetscCall(DMGetCoordinateDim(coarse, &dim));
3787: { /* winnow fine points that don't have global dofs out of the sf */
3788: PetscInt nleaves, l;
3789: const PetscInt *leaves;
3790: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;
3792: PetscCall(PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL));
3794: for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
3795: PetscInt p = leaves ? leaves[l] : l;
3797: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
3798: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
3799: if ((dof - cdof) > 0) numPointsWithDofs++;
3800: }
3801: PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
3802: for (l = 0, offset = 0; l < nleaves; l++) {
3803: PetscInt p = leaves ? leaves[l] : l;
3805: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
3806: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
3807: if ((dof - cdof) > 0) pointsWithDofs[offset++] = l;
3808: }
3809: PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded));
3810: PetscCall(PetscFree(pointsWithDofs));
3811: }
3812: /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
3813: PetscCall(PetscMalloc1(pEndC - pStartC, &maxChildIds));
3814: for (p = pStartC; p < pEndC; p++) maxChildIds[p - pStartC] = -2;
3815: PetscCall(PetscSFReduceBegin(coarseToFineEmbedded, MPIU_INT, cids, maxChildIds, MPIU_MAX));
3816: PetscCall(PetscSFReduceEnd(coarseToFineEmbedded, MPIU_INT, cids, maxChildIds, MPIU_MAX));
3818: PetscCall(DMGetLocalSection(coarse, &localCoarse));
3819: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
3821: PetscCall(DMPlexGetAnchors(coarse, &aSec, &aIS));
3822: PetscCall(ISGetIndices(aIS, &anchors));
3823: PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd));
3825: PetscCall(DMGetDefaultConstraints(coarse, &cSec, &cMat, NULL));
3826: PetscCall(PetscSectionGetChart(cSec, &cStart, &cEnd));
3828: /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
3829: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootValuesSec));
3830: PetscCall(PetscSectionSetChart(rootValuesSec, pStartC, pEndC));
3831: PetscCall(PetscSectionGetNumFields(localCoarse, &numFields));
3832: {
3833: PetscInt maxFields = PetscMax(1, numFields) + 1;
3834: PetscCall(PetscMalloc7(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &newOffsets, maxFields, &newOffsetsCopy, maxFields, &rowOffsets, maxFields, &numD, maxFields, &numO));
3835: }
3836: if (grad) {
3837: PetscInt i;
3839: PetscCall(VecGetDM(cellGeom, &cellDM));
3840: PetscCall(VecGetArrayRead(cellGeom, &cellGeomArray));
3841: PetscCall(VecGetDM(grad, &gradDM));
3842: PetscCall(VecGetArrayRead(grad, &gradArray));
3843: for (i = 0; i < PetscMax(1, numFields); i++) {
3844: PetscObject obj;
3845: PetscClassId id;
3847: PetscCall(DMGetField(coarse, i, NULL, &obj));
3848: PetscCall(PetscObjectGetClassId(obj, &id));
3849: if (id == PETSCFV_CLASSID) {
3850: fv = (PetscFV)obj;
3851: PetscCall(PetscFVGetNumComponents(fv, &numFVcomps));
3852: fvField = i;
3853: break;
3854: }
3855: }
3856: }
3858: for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
3859: PetscInt dof;
3860: PetscInt maxChildId = maxChildIds[p - pStartC];
3861: PetscInt numValues = 0;
3863: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3864: if (dof < 0) dof = -(dof + 1);
3865: offsets[0] = 0;
3866: newOffsets[0] = 0;
3867: if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
3868: PetscInt *closure = NULL, closureSize, cl;
3870: PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
3871: for (cl = 0; cl < closureSize; cl++) { /* get the closure */
3872: PetscInt c = closure[2 * cl], clDof;
3874: PetscCall(PetscSectionGetDof(localCoarse, c, &clDof));
3875: numValues += clDof;
3876: }
3877: PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
3878: } else if (maxChildId == -1) {
3879: PetscCall(PetscSectionGetDof(localCoarse, p, &numValues));
3880: }
3881: /* we will pack the column indices with the field offsets */
3882: if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) {
3883: /* also send the centroid, and the gradient */
3884: numValues += dim * (1 + numFVcomps);
3885: }
3886: PetscCall(PetscSectionSetDof(rootValuesSec, p, numValues));
3887: }
3888: PetscCall(PetscSectionSetUp(rootValuesSec));
3889: {
3890: PetscInt numRootValues;
3891: const PetscScalar *coarseArray;
3893: PetscCall(PetscSectionGetStorageSize(rootValuesSec, &numRootValues));
3894: PetscCall(PetscMalloc1(numRootValues, &rootValues));
3895: PetscCall(VecGetArrayRead(vecCoarseLocal, &coarseArray));
3896: for (p = pStartC; p < pEndC; p++) {
3897: PetscInt numValues;
3898: PetscInt pValOff;
3899: PetscScalar *pVal;
3900: PetscInt maxChildId = maxChildIds[p - pStartC];
3902: PetscCall(PetscSectionGetDof(rootValuesSec, p, &numValues));
3903: if (!numValues) continue;
3904: PetscCall(PetscSectionGetOffset(rootValuesSec, p, &pValOff));
3905: pVal = &rootValues[pValOff];
3906: if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
3907: PetscInt closureSize = numValues;
3908: PetscCall(DMPlexVecGetClosure(coarse, NULL, vecCoarseLocal, p, &closureSize, &pVal));
3909: if (grad && p >= cellStart && p < cellEnd) {
3910: PetscFVCellGeom *cg;
3911: PetscScalar *gradVals = NULL;
3912: PetscInt i;
3914: pVal += (numValues - dim * (1 + numFVcomps));
3916: PetscCall(DMPlexPointLocalRead(cellDM, p, cellGeomArray, (void *)&cg));
3917: for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i];
3918: pVal += dim;
3919: PetscCall(DMPlexPointGlobalRead(gradDM, p, gradArray, (void *)&gradVals));
3920: for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i];
3921: }
3922: } else if (maxChildId == -1) {
3923: PetscInt lDof, lOff, i;
3925: PetscCall(PetscSectionGetDof(localCoarse, p, &lDof));
3926: PetscCall(PetscSectionGetOffset(localCoarse, p, &lOff));
3927: for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i];
3928: }
3929: }
3930: PetscCall(VecRestoreArrayRead(vecCoarseLocal, &coarseArray));
3931: PetscCall(PetscFree(maxChildIds));
3932: }
3933: {
3934: PetscSF valuesSF;
3935: PetscInt *remoteOffsetsValues, numLeafValues;
3937: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafValuesSec));
3938: PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootValuesSec, &remoteOffsetsValues, leafValuesSec));
3939: PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootValuesSec, remoteOffsetsValues, leafValuesSec, &valuesSF));
3940: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
3941: PetscCall(PetscFree(remoteOffsetsValues));
3942: PetscCall(PetscSectionGetStorageSize(leafValuesSec, &numLeafValues));
3943: PetscCall(PetscMalloc1(numLeafValues, &leafValues));
3944: PetscCall(PetscSFBcastBegin(valuesSF, MPIU_SCALAR, rootValues, leafValues, MPI_REPLACE));
3945: PetscCall(PetscSFBcastEnd(valuesSF, MPIU_SCALAR, rootValues, leafValues, MPI_REPLACE));
3946: PetscCall(PetscSFDestroy(&valuesSF));
3947: PetscCall(PetscFree(rootValues));
3948: PetscCall(PetscSectionDestroy(&rootValuesSec));
3949: }
3950: PetscCall(DMGetLocalSection(fine, &localFine));
3951: {
3952: PetscInt maxDof;
3953: PetscInt *rowIndices;
3954: DM refTree;
3955: PetscInt **refPointFieldN;
3956: PetscScalar ***refPointFieldMats;
3957: PetscSection refConSec, refAnSec;
3958: PetscInt pRefStart, pRefEnd, leafStart, leafEnd;
3959: PetscScalar *pointWork;
3961: PetscCall(PetscSectionGetMaxDof(localFine, &maxDof));
3962: PetscCall(DMGetWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
3963: PetscCall(DMGetWorkArray(fine, maxDof, MPIU_SCALAR, &pointWork));
3964: PetscCall(DMPlexGetReferenceTree(fine, &refTree));
3965: PetscCall(DMCopyDisc(fine, refTree));
3966: PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
3967: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
3968: PetscCall(DMPlexGetAnchors(refTree, &refAnSec, NULL));
3969: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
3970: PetscCall(PetscSectionGetChart(leafValuesSec, &leafStart, &leafEnd));
3971: PetscCall(DMPlexGetSimplexOrBoxCells(fine, 0, &cellStart, &cellEnd));
3972: for (p = leafStart; p < leafEnd; p++) {
3973: PetscInt gDof, gcDof, gOff, lDof;
3974: PetscInt numValues, pValOff;
3975: PetscInt childId;
3976: const PetscScalar *pVal;
3977: const PetscScalar *fvGradData = NULL;
3979: PetscCall(PetscSectionGetDof(globalFine, p, &gDof));
3980: PetscCall(PetscSectionGetDof(localFine, p, &lDof));
3981: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &gcDof));
3982: if ((gDof - gcDof) <= 0) continue;
3983: PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
3984: PetscCall(PetscSectionGetDof(leafValuesSec, p, &numValues));
3985: if (!numValues) continue;
3986: PetscCall(PetscSectionGetOffset(leafValuesSec, p, &pValOff));
3987: pVal = &leafValues[pValOff];
3988: offsets[0] = 0;
3989: offsetsCopy[0] = 0;
3990: newOffsets[0] = 0;
3991: newOffsetsCopy[0] = 0;
3992: childId = cids[p - pStartF];
3993: if (numFields) {
3994: PetscInt f;
3995: for (f = 0; f < numFields; f++) {
3996: PetscInt rowDof;
3998: PetscCall(PetscSectionGetFieldDof(localFine, p, f, &rowDof));
3999: offsets[f + 1] = offsets[f] + rowDof;
4000: offsetsCopy[f + 1] = offsets[f + 1];
4001: /* TODO: closure indices */
4002: newOffsets[f + 1] = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]);
4003: }
4004: PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices));
4005: } else {
4006: offsets[0] = 0;
4007: offsets[1] = lDof;
4008: newOffsets[0] = 0;
4009: newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0];
4010: PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices));
4011: }
4012: if (childId == -1) { /* no child interpolation: one nnz per */
4013: PetscCall(VecSetValues(vecFine, numValues, rowIndices, pVal, INSERT_VALUES));
4014: } else {
4015: if (grad && p >= cellStart && p < cellEnd) {
4016: numValues -= (dim * (1 + numFVcomps));
4017: fvGradData = &pVal[numValues];
4018: }
4019: for (PetscInt f = 0; f < PetscMax(1, numFields); f++) {
4020: const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f];
4021: PetscInt numRows = offsets[f + 1] - offsets[f];
4022: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
4023: const PetscScalar *cVal = &pVal[newOffsets[f]];
4024: PetscScalar *rVal = &pointWork[offsets[f]];
4026: #if 0
4027: 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));
4028: #endif
4029: for (PetscInt i = 0; i < numRows; i++) {
4030: PetscScalar val = 0.;
4031: for (PetscInt j = 0; j < numCols; j++) val += childMat[i * numCols + j] * cVal[j];
4032: rVal[i] = val;
4033: }
4034: if (f == fvField && p >= cellStart && p < cellEnd) {
4035: PetscReal centroid[3];
4036: PetscScalar diff[3];
4037: const PetscScalar *parentCentroid = &fvGradData[0];
4038: const PetscScalar *gradient = &fvGradData[dim];
4040: PetscCall(DMPlexComputeCellGeometryFVM(fine, p, NULL, centroid, NULL));
4041: for (PetscInt i = 0; i < dim; i++) diff[i] = centroid[i] - parentCentroid[i];
4042: for (PetscInt i = 0; i < numFVcomps; i++) {
4043: PetscScalar val = 0.;
4045: for (PetscInt j = 0; j < dim; j++) val += gradient[dim * i + j] * diff[j];
4046: rVal[i] += val;
4047: }
4048: }
4049: PetscCall(VecSetValues(vecFine, numRows, &rowIndices[offsets[f]], rVal, INSERT_VALUES));
4050: }
4051: }
4052: }
4053: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
4054: PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_SCALAR, &pointWork));
4055: PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
4056: }
4057: PetscCall(PetscFree(leafValues));
4058: PetscCall(PetscSectionDestroy(&leafValuesSec));
4059: PetscCall(PetscFree7(offsets, offsetsCopy, newOffsets, newOffsetsCopy, rowOffsets, numD, numO));
4060: PetscCall(ISRestoreIndices(aIS, &anchors));
4061: PetscFunctionReturn(PETSC_SUCCESS);
4062: }
4064: static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids)
4065: {
4066: DM refTree;
4067: PetscSection multiRootSec, rootIndicesSec;
4068: PetscSection globalCoarse, globalFine;
4069: PetscSection localCoarse, localFine;
4070: PetscSection cSecRef;
4071: PetscInt *parentIndices, pRefStart, pRefEnd;
4072: PetscScalar *rootValues, *parentValues;
4073: Mat injRef;
4074: PetscInt numFields, maxDof;
4075: PetscInt pStartC, pEndC, pStartF, pEndF, p;
4076: PetscInt *offsets, *offsetsCopy, *rowOffsets;
4077: PetscLayout rowMap, colMap;
4078: PetscInt rowStart, rowEnd, colStart, colEnd;
4079: PetscScalar ***childrenMats = NULL; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */
4081: PetscFunctionBegin;
4082: /* get the templates for the fine-to-coarse injection from the reference tree */
4083: PetscCall(VecSetOption(vecFine, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE));
4084: PetscCall(VecSetOption(vecCoarse, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE));
4085: PetscCall(DMPlexGetReferenceTree(coarse, &refTree));
4086: PetscCall(DMCopyDisc(coarse, refTree));
4087: PetscCall(DMGetDefaultConstraints(refTree, &cSecRef, NULL, NULL));
4088: PetscCall(PetscSectionGetChart(cSecRef, &pRefStart, &pRefEnd));
4089: PetscCall(DMPlexReferenceTreeGetInjector(refTree, &injRef));
4091: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
4092: PetscCall(DMGetLocalSection(fine, &localFine));
4093: PetscCall(DMGetGlobalSection(fine, &globalFine));
4094: PetscCall(PetscSectionGetNumFields(localFine, &numFields));
4095: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
4096: PetscCall(DMGetLocalSection(coarse, &localCoarse));
4097: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
4098: PetscCall(PetscSectionGetMaxDof(localCoarse, &maxDof));
4099: {
4100: PetscInt maxFields = PetscMax(1, numFields) + 1;
4101: PetscCall(PetscMalloc3(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &rowOffsets));
4102: }
4104: PetscCall(DMPlexTransferInjectorTree(coarse, fine, coarseToFine, cids, vecFine, numFields, offsets, &multiRootSec, &rootIndicesSec, NULL, &rootValues));
4106: PetscCall(PetscMalloc2(maxDof, &parentIndices, maxDof, &parentValues));
4108: /* count indices */
4109: PetscCall(VecGetLayout(vecFine, &colMap));
4110: PetscCall(VecGetLayout(vecCoarse, &rowMap));
4111: PetscCall(PetscLayoutSetUp(rowMap));
4112: PetscCall(PetscLayoutSetUp(colMap));
4113: PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd));
4114: PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd));
4115: /* insert values */
4116: PetscCall(DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree, injRef, &childrenMats));
4117: for (p = pStartC; p < pEndC; p++) {
4118: PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
4119: PetscBool contribute = PETSC_FALSE;
4121: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
4122: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
4123: if ((dof - cdof) <= 0) continue;
4124: PetscCall(PetscSectionGetDof(localCoarse, p, &dof));
4125: PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));
4127: rowOffsets[0] = 0;
4128: offsetsCopy[0] = 0;
4129: if (numFields) {
4130: for (PetscInt 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: contribute = PETSC_TRUE;
4159: for (PetscInt m = 0; m < numIndices; m++) parentValues[m] = childValues[m];
4160: } else { /* contributions from children: sum with injectors from reference tree */
4161: PetscInt parentId, lim;
4163: contribute = PETSC_TRUE;
4164: PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL));
4166: lim = PetscMax(1, numFields);
4167: offsets[0] = 0;
4168: if (numFields) {
4169: for (PetscInt f = 0; f < numFields; f++) {
4170: PetscInt fDof;
4171: PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof));
4173: offsets[f + 1] = fDof + offsets[f];
4174: }
4175: } else {
4176: PetscInt cDof;
4178: PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof));
4179: offsets[1] = cDof;
4180: }
4181: for (PetscInt f = 0; f < lim; f++) {
4182: PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0];
4183: PetscInt n = offsets[f + 1] - offsets[f];
4184: PetscInt m = rowOffsets[f + 1] - rowOffsets[f];
4185: const PetscScalar *colValues = &childValues[offsets[f]];
4187: for (PetscInt i = 0; i < m; i++) {
4188: PetscScalar val = 0.;
4189: for (PetscInt j = 0; j < n; j++) val += childMat[n * i + j] * colValues[j];
4190: parentValues[rowOffsets[f] + i] += val;
4191: }
4192: }
4193: }
4194: }
4195: if (contribute) PetscCall(VecSetValues(vecCoarse, dof, parentIndices, parentValues, INSERT_VALUES));
4196: }
4197: PetscCall(PetscSectionDestroy(&multiRootSec));
4198: PetscCall(PetscSectionDestroy(&rootIndicesSec));
4199: PetscCall(PetscFree2(parentIndices, parentValues));
4200: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree, injRef, &childrenMats));
4201: PetscCall(PetscFree(rootValues));
4202: PetscCall(PetscFree3(offsets, offsetsCopy, rowOffsets));
4203: PetscFunctionReturn(PETSC_SUCCESS);
4204: }
4206: /*@
4207: DMPlexTransferVecTree - transfer a vector between two meshes that differ from each other by refinement/coarsening
4208: that can be represented by a common reference tree used by both. This routine can be used for a combination of
4209: coarsening and refinement at the same time.
4211: Collective
4213: Input Parameters:
4214: + dmIn - The `DMPLEX` mesh for the input vector
4215: . dmOut - The second `DMPLEX` mesh
4216: . vecIn - The input vector
4217: . sfRefine - A star forest indicating points in the mesh `dmIn` (roots in the star forest) that are parents to points in
4218: the mesh `dmOut` (leaves in the star forest), i.e. where `dmOut` is more refined than `dmIn`
4219: . sfCoarsen - A star forest indicating points in the mesh `dmOut` (roots in the star forest) that are parents to points in
4220: the mesh `dmIn` (leaves in the star forest), i.e. where `dmOut` is more coarsened than `dmIn`
4221: . cidsRefine - The childIds of the points in `dmOut`. These childIds relate back to the reference tree: childid[j] = k implies
4222: that mesh point j of `dmOut` was refined from a point in `dmIn` just as the mesh point k in the reference
4223: tree was refined from its parent. childid[j] = -1 indicates that the point j in `dmOut` is exactly
4224: equivalent to its root in `dmIn`, so no interpolation is necessary. childid[j] = -2 indicates that this
4225: point j in `dmOut` is not a leaf of `sfRefine`.
4226: . cidsCoarsen - The childIds of the points in `dmIn`. These childIds relate back to the reference tree: childid[j] = k implies
4227: that mesh point j of dmIn coarsens to a point in `dmOut` just as the mesh point k in the reference
4228: tree coarsens to its parent. childid[j] = -2 indicates that point j in `dmOut` is not a leaf in `sfCoarsen`.
4229: . useBCs - `PETSC_TRUE` indicates that boundary values should be inserted into `vecIn` before transfer.
4230: - time - Used if boundary values are time dependent.
4232: Output Parameter:
4233: . vecOut - Using interpolation and injection operators calculated on the reference tree, the transferred
4234: projection of `vecIn` from `dmIn` to `dmOut`. Note that any field discretized with a `PetscFV` finite volume
4235: method that uses gradient reconstruction will use reconstructed gradients when interpolating from
4236: coarse points to fine points.
4238: Level: developer
4240: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `Vec`, `PetscFV`, `DMPlexSetReferenceTree()`, `DMPlexGetReferenceTree()`, `PetscFVGetComputeGradients()`
4241: @*/
4242: PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time)
4243: {
4244: PetscFunctionBegin;
4245: PetscCall(VecSet(vecOut, 0.0));
4246: if (sfRefine) {
4247: Vec vecInLocal;
4248: DM dmGrad = NULL;
4249: Vec faceGeom = NULL, cellGeom = NULL, grad = NULL;
4251: PetscCall(DMGetLocalVector(dmIn, &vecInLocal));
4252: PetscCall(VecSet(vecInLocal, 0.0));
4253: {
4254: PetscInt numFields;
4256: PetscCall(DMGetNumFields(dmIn, &numFields));
4257: for (PetscInt i = 0; i < numFields; i++) {
4258: PetscObject obj;
4259: PetscClassId classid;
4261: PetscCall(DMGetField(dmIn, i, NULL, &obj));
4262: PetscCall(PetscObjectGetClassId(obj, &classid));
4263: if (classid == PETSCFV_CLASSID) {
4264: PetscCall(DMPlexGetDataFVM(dmIn, (PetscFV)obj, &cellGeom, &faceGeom, &dmGrad));
4265: break;
4266: }
4267: }
4268: }
4269: if (useBCs) PetscCall(DMPlexInsertBoundaryValues(dmIn, PETSC_TRUE, vecInLocal, time, faceGeom, cellGeom, NULL));
4270: PetscCall(DMGlobalToLocalBegin(dmIn, vecIn, INSERT_VALUES, vecInLocal));
4271: PetscCall(DMGlobalToLocalEnd(dmIn, vecIn, INSERT_VALUES, vecInLocal));
4272: if (dmGrad) {
4273: PetscCall(DMGetGlobalVector(dmGrad, &grad));
4274: PetscCall(DMPlexReconstructGradientsFVM(dmIn, vecInLocal, grad));
4275: }
4276: PetscCall(DMPlexTransferVecTree_Interpolate(dmIn, vecInLocal, dmOut, vecOut, sfRefine, cidsRefine, grad, cellGeom));
4277: PetscCall(DMRestoreLocalVector(dmIn, &vecInLocal));
4278: if (dmGrad) PetscCall(DMRestoreGlobalVector(dmGrad, &grad));
4279: }
4280: if (sfCoarsen) PetscCall(DMPlexTransferVecTree_Inject(dmIn, vecIn, dmOut, vecOut, sfCoarsen, cidsCoarsen));
4281: PetscCall(VecAssemblyBegin(vecOut));
4282: PetscCall(VecAssemblyEnd(vecOut));
4283: PetscFunctionReturn(PETSC_SUCCESS);
4284: }