Actual source code: plextree.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/isimpl.h>
3: #include <petsc/private/petscfeimpl.h>
4: #include <petscsf.h>
5: #include <petscds.h>
7: /* hierarchy routines */
9: /*@
10: DMPlexSetReferenceTree - set the reference tree for hierarchically non-conforming meshes.
12: Not Collective
14: Input Parameters:
15: + dm - The `DMPLEX` object
16: - ref - The reference tree `DMPLEX` object
18: Level: intermediate
20: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetReferenceTree()`, `DMPlexCreateDefaultReferenceTree()`
21: @*/
22: PetscErrorCode DMPlexSetReferenceTree(DM dm, DM ref)
23: {
24: DM_Plex *mesh = (DM_Plex *)dm->data;
26: PetscFunctionBegin;
29: PetscCall(PetscObjectReference((PetscObject)ref));
30: PetscCall(DMDestroy(&mesh->referenceTree));
31: mesh->referenceTree = ref;
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: /*@
36: DMPlexGetReferenceTree - get the reference tree for hierarchically non-conforming meshes.
38: Not Collective
40: Input Parameter:
41: . dm - The `DMPLEX` object
43: Output Parameter:
44: . ref - The reference tree `DMPLEX` object
46: Level: intermediate
48: Developer Notes:
49: The reference tree is shallow copied during `DMClone()`, thus it is may be shared by different `DM`s.
50: It is not a topological-only object, since some parts of the library use its local section to compute
51: interpolation and injection matrices. This may lead to unexpected failures during those calls.
53: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetReferenceTree()`, `DMPlexCreateDefaultReferenceTree()`
54: @*/
55: PetscErrorCode DMPlexGetReferenceTree(DM dm, DM *ref)
56: {
57: DM_Plex *mesh = (DM_Plex *)dm->data;
59: PetscFunctionBegin;
61: PetscAssertPointer(ref, 2);
62: *ref = mesh->referenceTree;
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: static PetscErrorCode DMPlexReferenceTreeGetChildSymmetry_Default(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
67: {
68: PetscInt coneSize, dStart, dEnd, dim, ABswap, oAvert, oBvert, ABswapVert;
70: PetscFunctionBegin;
71: if (parentOrientA == parentOrientB) {
72: if (childOrientB) *childOrientB = childOrientA;
73: if (childB) *childB = childA;
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
76: for (dim = 0; dim < 3; dim++) {
77: PetscCall(DMPlexGetDepthStratum(dm, dim, &dStart, &dEnd));
78: if (parent >= dStart && parent <= dEnd) break;
79: }
80: PetscCheck(dim <= 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot perform child symmetry for %" PetscInt_FMT "-cells", dim);
81: PetscCheck(dim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "A vertex has no children");
82: if (childA < dStart || childA >= dEnd) {
83: /* this is a lower-dimensional child: bootstrap */
84: PetscInt size, i, sA = -1, sB, sOrientB, sConeSize;
85: const PetscInt *supp, *coneA, *coneB, *oA, *oB;
87: PetscCall(DMPlexGetSupportSize(dm, childA, &size));
88: PetscCall(DMPlexGetSupport(dm, childA, &supp));
90: /* find a point sA in supp(childA) that has the same parent */
91: for (i = 0; i < size; i++) {
92: PetscInt sParent;
94: sA = supp[i];
95: if (sA == parent) continue;
96: PetscCall(DMPlexGetTreeParent(dm, sA, &sParent, NULL));
97: if (sParent == parent) break;
98: }
99: PetscCheck(i != size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "could not find support in children");
100: /* find out which point sB is in an equivalent position to sA under
101: * parentOrientB */
102: PetscCall(DMPlexReferenceTreeGetChildSymmetry_Default(dm, parent, parentOrientA, 0, sA, parentOrientB, &sOrientB, &sB));
103: PetscCall(DMPlexGetConeSize(dm, sA, &sConeSize));
104: PetscCall(DMPlexGetCone(dm, sA, &coneA));
105: PetscCall(DMPlexGetCone(dm, sB, &coneB));
106: PetscCall(DMPlexGetConeOrientation(dm, sA, &oA));
107: PetscCall(DMPlexGetConeOrientation(dm, sB, &oB));
108: /* step through the cone of sA in natural order */
109: for (i = 0; i < sConeSize; i++) {
110: if (coneA[i] == childA) {
111: /* if childA is at position i in coneA,
112: * then we want the point that is at sOrientB*i in coneB */
113: PetscInt j = (sOrientB >= 0) ? ((sOrientB + i) % sConeSize) : ((sConeSize - (sOrientB + 1) - i) % sConeSize);
114: if (childB) *childB = coneB[j];
115: if (childOrientB) {
116: DMPolytopeType ct;
117: PetscInt oBtrue;
119: PetscCall(DMPlexGetConeSize(dm, childA, &coneSize));
120: /* compose sOrientB and oB[j] */
121: PetscCheck(coneSize == 0 || coneSize == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected a vertex or an edge");
122: ct = coneSize ? DM_POLYTOPE_SEGMENT : DM_POLYTOPE_POINT;
123: /* we may have to flip an edge */
124: oBtrue = (sOrientB >= 0) ? oB[j] : DMPolytopeTypeComposeOrientation(ct, -1, oB[j]);
125: oBtrue = DMPolytopeConvertNewOrientation_Internal(ct, oBtrue);
126: ABswap = DihedralSwap(coneSize, DMPolytopeConvertNewOrientation_Internal(ct, oA[i]), oBtrue);
127: *childOrientB = DihedralCompose(coneSize, childOrientA, ABswap);
128: }
129: break;
130: }
131: }
132: PetscCheck(i != sConeSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "support cone mismatch");
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }
135: /* get the cone size and symmetry swap */
136: PetscCall(DMPlexGetConeSize(dm, parent, &coneSize));
137: ABswap = DihedralSwap(coneSize, parentOrientA, parentOrientB);
138: if (dim == 2) {
139: /* orientations refer to cones: we want them to refer to vertices:
140: * if it's a rotation, they are the same, but if the order is reversed, a
141: * permutation that puts side i first does *not* put vertex i first */
142: oAvert = (parentOrientA >= 0) ? parentOrientA : -((-parentOrientA % coneSize) + 1);
143: oBvert = (parentOrientB >= 0) ? parentOrientB : -((-parentOrientB % coneSize) + 1);
144: ABswapVert = DihedralSwap(coneSize, oAvert, oBvert);
145: } else {
146: ABswapVert = ABswap;
147: }
148: if (childB) {
149: /* assume that each child corresponds to a vertex, in the same order */
150: PetscInt p, posA = -1, numChildren, i;
151: const PetscInt *children;
153: /* count which position the child is in */
154: PetscCall(DMPlexGetTreeChildren(dm, parent, &numChildren, &children));
155: for (i = 0; i < numChildren; i++) {
156: p = children[i];
157: if (p == childA) {
158: posA = i;
159: break;
160: }
161: }
162: if (posA >= coneSize) {
163: /* this is the triangle in the middle of a uniformly refined triangle: it
164: * is invariant */
165: PetscCheck(dim == 2 && posA == 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Expected a middle triangle, got something else");
166: *childB = childA;
167: } else {
168: /* figure out position B by applying ABswapVert */
169: PetscInt posB;
171: posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize - (ABswapVert + 1) - posA) % coneSize);
172: if (childB) *childB = children[posB];
173: }
174: }
175: if (childOrientB) *childOrientB = DihedralCompose(coneSize, childOrientA, ABswap);
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: /*@
180: DMPlexReferenceTreeGetChildSymmetry - Given a reference tree, transform a childid and orientation from one parent frame to another
182: Input Parameters:
183: + dm - the reference tree `DMPLEX` object
184: . parent - the parent point
185: . parentOrientA - the reference orientation for describing the parent
186: . childOrientA - the reference orientation for describing the child
187: . childA - the reference childID for describing the child
188: - parentOrientB - the new orientation for describing the parent
190: Output Parameters:
191: + childOrientB - if not `NULL`, set to the new orientation for describing the child
192: - childB - if not `NULL`, the new childID for describing the child
194: Level: developer
196: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetReferenceTree()`, `DMPlexSetReferenceTree()`, `DMPlexSetTree()`
197: @*/
198: PetscErrorCode DMPlexReferenceTreeGetChildSymmetry(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
199: {
200: DM_Plex *mesh = (DM_Plex *)dm->data;
202: PetscFunctionBegin;
204: PetscCheck(mesh->getchildsymmetry, PETSC_COMM_SELF, PETSC_ERR_SUP, "DMPlexReferenceTreeGetChildSymmetry not implemented");
205: PetscCall(mesh->getchildsymmetry(dm, parent, parentOrientA, childOrientA, childA, parentOrientB, childOrientB, childB));
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: static PetscErrorCode DMPlexSetTree_Internal(DM, PetscSection, PetscInt *, PetscInt *, PetscBool, PetscBool);
211: PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
212: {
213: PetscFunctionBegin;
214: PetscCall(DMPlexSetTree_Internal(dm, parentSection, parents, childIDs, PETSC_TRUE, PETSC_FALSE));
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: PetscErrorCode DMPlexCreateReferenceTree_Union(DM K, DM Kref, const char *labelName, DM *ref)
219: {
220: MPI_Comm comm;
221: PetscInt dim, p, pStart, pEnd, pRefStart, pRefEnd, d, offset, parentSize, *parents, *childIDs;
222: PetscInt *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts;
223: DMLabel identity, identityRef;
224: PetscSection unionSection, unionConeSection, parentSection;
225: PetscScalar *unionCoords;
226: IS perm;
228: PetscFunctionBegin;
229: comm = PetscObjectComm((PetscObject)K);
230: PetscCall(DMGetDimension(K, &dim));
231: PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
232: PetscCall(DMGetLabel(K, labelName, &identity));
233: PetscCall(DMGetLabel(Kref, labelName, &identityRef));
234: PetscCall(DMPlexGetChart(Kref, &pRefStart, &pRefEnd));
235: PetscCall(PetscSectionCreate(comm, &unionSection));
236: PetscCall(PetscSectionSetChart(unionSection, 0, (pEnd - pStart) + (pRefEnd - pRefStart)));
237: /* count points that will go in the union */
238: for (p = pStart; p < pEnd; p++) PetscCall(PetscSectionSetDof(unionSection, p - pStart, 1));
239: for (p = pRefStart; p < pRefEnd; p++) {
240: PetscInt q, qSize;
241: PetscCall(DMLabelGetValue(identityRef, p, &q));
242: PetscCall(DMLabelGetStratumSize(identityRef, q, &qSize));
243: if (qSize > 1) PetscCall(PetscSectionSetDof(unionSection, p - pRefStart + (pEnd - pStart), 1));
244: }
245: PetscCall(PetscMalloc1(pEnd - pStart + pRefEnd - pRefStart, &permvals));
246: offset = 0;
247: /* stratify points in the union by topological dimension */
248: for (d = 0; d <= dim; d++) {
249: PetscInt cStart, cEnd, c;
251: PetscCall(DMPlexGetHeightStratum(K, d, &cStart, &cEnd));
252: for (c = cStart; c < cEnd; c++) permvals[offset++] = c;
254: PetscCall(DMPlexGetHeightStratum(Kref, d, &cStart, &cEnd));
255: for (c = cStart; c < cEnd; c++) permvals[offset++] = c + (pEnd - pStart);
256: }
257: PetscCall(ISCreateGeneral(comm, (pEnd - pStart) + (pRefEnd - pRefStart), permvals, PETSC_OWN_POINTER, &perm));
258: PetscCall(PetscSectionSetPermutation(unionSection, perm));
259: PetscCall(PetscSectionSetUp(unionSection));
260: PetscCall(PetscSectionGetStorageSize(unionSection, &numUnionPoints));
261: PetscCall(PetscMalloc2(numUnionPoints, &coneSizes, dim + 1, &numDimPoints));
262: /* count dimension points */
263: for (d = 0; d <= dim; d++) {
264: PetscInt cStart, cOff, cOff2;
265: PetscCall(DMPlexGetHeightStratum(K, d, &cStart, NULL));
266: PetscCall(PetscSectionGetOffset(unionSection, cStart - pStart, &cOff));
267: if (d < dim) {
268: PetscCall(DMPlexGetHeightStratum(K, d + 1, &cStart, NULL));
269: PetscCall(PetscSectionGetOffset(unionSection, cStart - pStart, &cOff2));
270: } else {
271: cOff2 = numUnionPoints;
272: }
273: numDimPoints[dim - d] = cOff2 - cOff;
274: }
275: PetscCall(PetscSectionCreate(comm, &unionConeSection));
276: PetscCall(PetscSectionSetChart(unionConeSection, 0, numUnionPoints));
277: /* count the cones in the union */
278: for (p = pStart; p < pEnd; p++) {
279: PetscInt dof, uOff;
281: PetscCall(DMPlexGetConeSize(K, p, &dof));
282: PetscCall(PetscSectionGetOffset(unionSection, p - pStart, &uOff));
283: PetscCall(PetscSectionSetDof(unionConeSection, uOff, dof));
284: coneSizes[uOff] = dof;
285: }
286: for (p = pRefStart; p < pRefEnd; p++) {
287: PetscInt dof, uDof, uOff;
289: PetscCall(DMPlexGetConeSize(Kref, p, &dof));
290: PetscCall(PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart), &uDof));
291: PetscCall(PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart), &uOff));
292: if (uDof) {
293: PetscCall(PetscSectionSetDof(unionConeSection, uOff, dof));
294: coneSizes[uOff] = dof;
295: }
296: }
297: PetscCall(PetscSectionSetUp(unionConeSection));
298: PetscCall(PetscSectionGetStorageSize(unionConeSection, &numCones));
299: PetscCall(PetscMalloc2(numCones, &unionCones, numCones, &unionOrientations));
300: /* write the cones in the union */
301: for (p = pStart; p < pEnd; p++) {
302: PetscInt dof, uOff, c, cOff;
303: const PetscInt *cone, *orientation;
305: PetscCall(DMPlexGetConeSize(K, p, &dof));
306: PetscCall(DMPlexGetCone(K, p, &cone));
307: PetscCall(DMPlexGetConeOrientation(K, p, &orientation));
308: PetscCall(PetscSectionGetOffset(unionSection, p - pStart, &uOff));
309: PetscCall(PetscSectionGetOffset(unionConeSection, uOff, &cOff));
310: for (c = 0; c < dof; c++) {
311: PetscInt e, eOff;
312: e = cone[c];
313: PetscCall(PetscSectionGetOffset(unionSection, e - pStart, &eOff));
314: unionCones[cOff + c] = eOff;
315: unionOrientations[cOff + c] = orientation[c];
316: }
317: }
318: for (p = pRefStart; p < pRefEnd; p++) {
319: PetscInt dof, uDof, uOff, c, cOff;
320: const PetscInt *cone, *orientation;
322: PetscCall(DMPlexGetConeSize(Kref, p, &dof));
323: PetscCall(DMPlexGetCone(Kref, p, &cone));
324: PetscCall(DMPlexGetConeOrientation(Kref, p, &orientation));
325: PetscCall(PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart), &uDof));
326: PetscCall(PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart), &uOff));
327: if (uDof) {
328: PetscCall(PetscSectionGetOffset(unionConeSection, uOff, &cOff));
329: for (c = 0; c < dof; c++) {
330: PetscInt e, eOff, eDof;
332: e = cone[c];
333: PetscCall(PetscSectionGetDof(unionSection, e - pRefStart + (pEnd - pStart), &eDof));
334: if (eDof) {
335: PetscCall(PetscSectionGetOffset(unionSection, e - pRefStart + (pEnd - pStart), &eOff));
336: } else {
337: PetscCall(DMLabelGetValue(identityRef, e, &e));
338: PetscCall(PetscSectionGetOffset(unionSection, e - pStart, &eOff));
339: }
340: unionCones[cOff + c] = eOff;
341: unionOrientations[cOff + c] = orientation[c];
342: }
343: }
344: }
345: /* get the coordinates */
346: {
347: PetscInt vStart, vEnd, vRefStart, vRefEnd, v, vDof, vOff;
348: PetscSection KcoordsSec, KrefCoordsSec;
349: Vec KcoordsVec, KrefCoordsVec;
350: PetscScalar *Kcoords;
352: PetscCall(DMGetCoordinateSection(K, &KcoordsSec));
353: PetscCall(DMGetCoordinatesLocal(K, &KcoordsVec));
354: PetscCall(DMGetCoordinateSection(Kref, &KrefCoordsSec));
355: PetscCall(DMGetCoordinatesLocal(Kref, &KrefCoordsVec));
357: numVerts = numDimPoints[0];
358: PetscCall(PetscMalloc1(numVerts * dim, &unionCoords));
359: PetscCall(DMPlexGetDepthStratum(K, 0, &vStart, &vEnd));
361: offset = 0;
362: for (v = vStart; v < vEnd; v++) {
363: PetscCall(PetscSectionGetOffset(unionSection, v - pStart, &vOff));
364: PetscCall(VecGetValuesSection(KcoordsVec, KcoordsSec, v, &Kcoords));
365: for (d = 0; d < dim; d++) unionCoords[offset * dim + d] = Kcoords[d];
366: offset++;
367: }
368: PetscCall(DMPlexGetDepthStratum(Kref, 0, &vRefStart, &vRefEnd));
369: for (v = vRefStart; v < vRefEnd; v++) {
370: PetscCall(PetscSectionGetDof(unionSection, v - pRefStart + (pEnd - pStart), &vDof));
371: PetscCall(PetscSectionGetOffset(unionSection, v - pRefStart + (pEnd - pStart), &vOff));
372: PetscCall(VecGetValuesSection(KrefCoordsVec, KrefCoordsSec, v, &Kcoords));
373: if (vDof) {
374: for (d = 0; d < dim; d++) unionCoords[offset * dim + d] = Kcoords[d];
375: offset++;
376: }
377: }
378: }
379: PetscCall(DMCreate(comm, ref));
380: PetscCall(DMSetType(*ref, DMPLEX));
381: PetscCall(DMSetDimension(*ref, dim));
382: PetscCall(DMPlexCreateFromDAG(*ref, dim, numDimPoints, coneSizes, unionCones, unionOrientations, unionCoords));
383: /* set the tree */
384: PetscCall(PetscSectionCreate(comm, &parentSection));
385: PetscCall(PetscSectionSetChart(parentSection, 0, numUnionPoints));
386: for (p = pRefStart; p < pRefEnd; p++) {
387: PetscInt uDof, uOff;
389: PetscCall(PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart), &uDof));
390: PetscCall(PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart), &uOff));
391: if (uDof) PetscCall(PetscSectionSetDof(parentSection, uOff, 1));
392: }
393: PetscCall(PetscSectionSetUp(parentSection));
394: PetscCall(PetscSectionGetStorageSize(parentSection, &parentSize));
395: PetscCall(PetscMalloc2(parentSize, &parents, parentSize, &childIDs));
396: for (p = pRefStart; p < pRefEnd; p++) {
397: PetscInt uDof, uOff;
399: PetscCall(PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart), &uDof));
400: PetscCall(PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart), &uOff));
401: if (uDof) {
402: PetscInt pOff, parent, parentU;
403: PetscCall(PetscSectionGetOffset(parentSection, uOff, &pOff));
404: PetscCall(DMLabelGetValue(identityRef, p, &parent));
405: PetscCall(PetscSectionGetOffset(unionSection, parent - pStart, &parentU));
406: parents[pOff] = parentU;
407: childIDs[pOff] = uOff;
408: }
409: }
410: PetscCall(DMPlexCreateReferenceTree_SetTree(*ref, parentSection, parents, childIDs));
411: PetscCall(PetscSectionDestroy(&parentSection));
412: PetscCall(PetscFree2(parents, childIDs));
414: /* clean up */
415: PetscCall(PetscSectionDestroy(&unionSection));
416: PetscCall(PetscSectionDestroy(&unionConeSection));
417: PetscCall(ISDestroy(&perm));
418: PetscCall(PetscFree(unionCoords));
419: PetscCall(PetscFree2(unionCones, unionOrientations));
420: PetscCall(PetscFree2(coneSizes, numDimPoints));
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: /*@
425: DMPlexCreateDefaultReferenceTree - create a reference tree for isotropic hierarchical mesh refinement.
427: Collective
429: Input Parameters:
430: + comm - the MPI communicator
431: . dim - the spatial dimension
432: - simplex - Flag for simplex, otherwise use a tensor-product cell
434: Output Parameter:
435: . ref - the reference tree `DMPLEX` object
437: Level: intermediate
439: .seealso: `DMPlexSetReferenceTree()`, `DMPlexGetReferenceTree()`
440: @*/
441: PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref)
442: {
443: DM_Plex *mesh;
444: DM K, Kref;
445: PetscInt p, pStart, pEnd;
446: DMLabel identity;
448: PetscFunctionBegin;
449: #if 1
450: comm = PETSC_COMM_SELF;
451: #endif
452: /* create a reference element */
453: PetscCall(DMPlexCreateReferenceCell(comm, DMPolytopeTypeSimpleShape(dim, simplex), &K));
454: PetscCall(DMCreateLabel(K, "identity"));
455: PetscCall(DMGetLabel(K, "identity", &identity));
456: PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
457: for (p = pStart; p < pEnd; p++) PetscCall(DMLabelSetValue(identity, p, p));
458: /* refine it */
459: PetscCall(DMRefine(K, comm, &Kref));
461: /* the reference tree is the union of these two, without duplicating
462: * points that appear in both */
463: PetscCall(DMPlexCreateReferenceTree_Union(K, Kref, "identity", ref));
464: mesh = (DM_Plex *)(*ref)->data;
465: mesh->getchildsymmetry = DMPlexReferenceTreeGetChildSymmetry_Default;
466: PetscCall(DMDestroy(&K));
467: PetscCall(DMDestroy(&Kref));
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }
471: static PetscErrorCode DMPlexTreeSymmetrize(DM dm)
472: {
473: DM_Plex *mesh = (DM_Plex *)dm->data;
474: PetscSection childSec, pSec;
475: PetscInt p, pSize, cSize, parMax = PETSC_INT_MIN, parMin = PETSC_INT_MAX;
476: PetscInt *offsets, *children, pStart, pEnd;
478: PetscFunctionBegin;
480: PetscCall(PetscSectionDestroy(&mesh->childSection));
481: PetscCall(PetscFree(mesh->children));
482: pSec = mesh->parentSection;
483: if (!pSec) PetscFunctionReturn(PETSC_SUCCESS);
484: PetscCall(PetscSectionGetStorageSize(pSec, &pSize));
485: for (p = 0; p < pSize; p++) {
486: PetscInt par = mesh->parents[p];
488: parMax = PetscMax(parMax, par + 1);
489: parMin = PetscMin(parMin, par);
490: }
491: if (parMin > parMax) {
492: parMin = -1;
493: parMax = -1;
494: }
495: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)pSec), &childSec));
496: PetscCall(PetscSectionSetChart(childSec, parMin, parMax));
497: for (p = 0; p < pSize; p++) {
498: PetscInt par = mesh->parents[p];
500: PetscCall(PetscSectionAddDof(childSec, par, 1));
501: }
502: PetscCall(PetscSectionSetUp(childSec));
503: PetscCall(PetscSectionGetStorageSize(childSec, &cSize));
504: PetscCall(PetscMalloc1(cSize, &children));
505: PetscCall(PetscCalloc1(parMax - parMin, &offsets));
506: PetscCall(PetscSectionGetChart(pSec, &pStart, &pEnd));
507: for (p = pStart; p < pEnd; p++) {
508: PetscInt dof, off, i;
510: PetscCall(PetscSectionGetDof(pSec, p, &dof));
511: PetscCall(PetscSectionGetOffset(pSec, p, &off));
512: for (i = 0; i < dof; i++) {
513: PetscInt par = mesh->parents[off + i], cOff;
515: PetscCall(PetscSectionGetOffset(childSec, par, &cOff));
516: children[cOff + offsets[par - parMin]++] = p;
517: }
518: }
519: mesh->childSection = childSec;
520: mesh->children = children;
521: PetscCall(PetscFree(offsets));
522: PetscFunctionReturn(PETSC_SUCCESS);
523: }
525: static PetscErrorCode AnchorsFlatten(PetscSection section, IS is, PetscSection *sectionNew, IS *isNew)
526: {
527: PetscInt pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL;
528: const PetscInt *vals;
529: PetscSection secNew;
530: PetscBool anyNew, globalAnyNew;
531: PetscBool compress;
533: PetscFunctionBegin;
534: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
535: PetscCall(ISGetLocalSize(is, &size));
536: PetscCall(ISGetIndices(is, &vals));
537: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), &secNew));
538: PetscCall(PetscSectionSetChart(secNew, pStart, pEnd));
539: for (i = 0; i < size; i++) {
540: PetscInt dof;
542: p = vals[i];
543: if (p < pStart || p >= pEnd) continue;
544: PetscCall(PetscSectionGetDof(section, p, &dof));
545: if (dof) break;
546: }
547: if (i == size) {
548: PetscCall(PetscSectionSetUp(secNew));
549: anyNew = PETSC_FALSE;
550: compress = PETSC_FALSE;
551: sizeNew = 0;
552: } else {
553: anyNew = PETSC_TRUE;
554: for (p = pStart; p < pEnd; p++) {
555: PetscInt dof, off;
557: PetscCall(PetscSectionGetDof(section, p, &dof));
558: PetscCall(PetscSectionGetOffset(section, p, &off));
559: for (i = 0; i < dof; i++) {
560: PetscInt q = vals[off + i], qDof = 0;
562: if (q >= pStart && q < pEnd) PetscCall(PetscSectionGetDof(section, q, &qDof));
563: if (qDof) PetscCall(PetscSectionAddDof(secNew, p, qDof));
564: else PetscCall(PetscSectionAddDof(secNew, p, 1));
565: }
566: }
567: PetscCall(PetscSectionSetUp(secNew));
568: PetscCall(PetscSectionGetStorageSize(secNew, &sizeNew));
569: PetscCall(PetscMalloc1(sizeNew, &valsNew));
570: compress = PETSC_FALSE;
571: for (p = pStart; p < pEnd; p++) {
572: PetscInt dof, off, count, offNew, dofNew;
574: PetscCall(PetscSectionGetDof(section, p, &dof));
575: PetscCall(PetscSectionGetOffset(section, p, &off));
576: PetscCall(PetscSectionGetDof(secNew, p, &dofNew));
577: PetscCall(PetscSectionGetOffset(secNew, p, &offNew));
578: count = 0;
579: for (i = 0; i < dof; i++) {
580: PetscInt q = vals[off + i], qDof = 0, qOff = 0, j;
582: if (q >= pStart && q < pEnd) {
583: PetscCall(PetscSectionGetDof(section, q, &qDof));
584: PetscCall(PetscSectionGetOffset(section, q, &qOff));
585: }
586: if (qDof) {
587: PetscInt oldCount = count;
589: for (j = 0; j < qDof; j++) {
590: PetscInt k, r = vals[qOff + j];
592: for (k = 0; k < oldCount; k++) {
593: if (valsNew[offNew + k] == r) break;
594: }
595: if (k == oldCount) valsNew[offNew + count++] = r;
596: }
597: } else {
598: PetscInt k, oldCount = count;
600: for (k = 0; k < oldCount; k++) {
601: if (valsNew[offNew + k] == q) break;
602: }
603: if (k == oldCount) valsNew[offNew + count++] = q;
604: }
605: }
606: if (count < dofNew) {
607: PetscCall(PetscSectionSetDof(secNew, p, count));
608: compress = PETSC_TRUE;
609: }
610: }
611: }
612: PetscCall(ISRestoreIndices(is, &vals));
613: PetscCallMPI(MPIU_Allreduce(&anyNew, &globalAnyNew, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)secNew)));
614: if (!globalAnyNew) {
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 d, dim;
940: /* add the canonical label */
941: PetscCall(DMGetDimension(dm, &dim));
942: PetscCall(DMCreateLabel(dm, "canonical"));
943: for (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, i;
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 (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: PetscInt k;
1465: for (k = 0; k < cDof; k++) {
1466: for (j = 0; j < aDof; j++) refPointFieldMats[p - pRefStart][f][k * numCols + colOff + j] *= flip[j];
1467: }
1468: }
1469: colOff += aDof;
1470: }
1471: }
1472: if (numFields) {
1473: PetscCall(PetscSectionRestoreFieldPointSyms(refSection, f, closureSize, closure, &perms, &flips));
1474: } else {
1475: PetscCall(PetscSectionRestorePointSyms(refSection, closureSize, closure, &perms, &flips));
1476: }
1477: }
1478: PetscCall(DMPlexRestoreTransitiveClosure(refTree, parent, PETSC_TRUE, &closureSize, &closure));
1479: }
1480: *childrenMats = refPointFieldMats;
1481: *childrenN = refPointFieldN;
1482: PetscCall(ISRestoreIndices(refAnIS, &refAnchors));
1483: PetscCall(PetscFree(rows));
1484: PetscCall(PetscFree(cols));
1485: PetscFunctionReturn(PETSC_SUCCESS);
1486: }
1488: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1489: {
1490: PetscDS ds;
1491: PetscInt **refPointFieldN;
1492: PetscScalar ***refPointFieldMats;
1493: PetscInt numFields, maxFields, pRefStart, pRefEnd, p, f;
1494: PetscSection refConSec;
1496: PetscFunctionBegin;
1497: refPointFieldN = *childrenN;
1498: *childrenN = NULL;
1499: refPointFieldMats = *childrenMats;
1500: *childrenMats = NULL;
1501: PetscCall(DMGetDS(refTree, &ds));
1502: PetscCall(PetscDSGetNumFields(ds, &numFields));
1503: maxFields = PetscMax(1, numFields);
1504: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
1505: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
1506: for (p = pRefStart; p < pRefEnd; p++) {
1507: PetscInt parent, pDof;
1509: PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL));
1510: PetscCall(PetscSectionGetDof(refConSec, p, &pDof));
1511: if (!pDof || parent == p) continue;
1513: for (f = 0; f < maxFields; f++) {
1514: PetscInt cDof;
1516: if (numFields) {
1517: PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof));
1518: } else {
1519: PetscCall(PetscSectionGetDof(refConSec, p, &cDof));
1520: }
1522: PetscCall(PetscFree(refPointFieldMats[p - pRefStart][f]));
1523: }
1524: PetscCall(PetscFree(refPointFieldMats[p - pRefStart]));
1525: PetscCall(PetscFree(refPointFieldN[p - pRefStart]));
1526: }
1527: PetscCall(PetscFree(refPointFieldMats));
1528: PetscCall(PetscFree(refPointFieldN));
1529: PetscFunctionReturn(PETSC_SUCCESS);
1530: }
1532: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm, PetscSection section, PetscSection conSec, Mat cMat)
1533: {
1534: DM refTree;
1535: PetscDS ds;
1536: Mat refCmat;
1537: PetscInt numFields, maxFields, f, pRefStart, pRefEnd, p, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN;
1538: PetscScalar ***refPointFieldMats, *pointWork;
1539: PetscSection refConSec, refAnSec, anSec;
1540: IS refAnIS, anIS;
1541: const PetscInt *anchors;
1543: PetscFunctionBegin;
1545: PetscCall(DMGetDS(dm, &ds));
1546: PetscCall(PetscDSGetNumFields(ds, &numFields));
1547: maxFields = PetscMax(1, numFields);
1548: PetscCall(DMPlexGetReferenceTree(dm, &refTree));
1549: PetscCall(DMCopyDisc(dm, refTree));
1550: PetscCall(DMSetLocalSection(refTree, NULL));
1551: PetscCall(DMSetDefaultConstraints(refTree, NULL, NULL, NULL));
1552: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, &refCmat, NULL));
1553: PetscCall(DMPlexGetAnchors(refTree, &refAnSec, &refAnIS));
1554: PetscCall(DMPlexGetAnchors(dm, &anSec, &anIS));
1555: PetscCall(ISGetIndices(anIS, &anchors));
1556: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
1557: PetscCall(PetscSectionGetChart(conSec, &conStart, &conEnd));
1558: PetscCall(PetscSectionGetMaxDof(refConSec, &maxDof));
1559: PetscCall(PetscSectionGetMaxDof(refAnSec, &maxAnDof));
1560: PetscCall(PetscMalloc1(maxDof * maxDof * maxAnDof, &pointWork));
1562: /* step 1: get submats for every constrained point in the reference tree */
1563: PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
1565: /* step 2: compute the preorder */
1566: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1567: PetscCall(PetscMalloc2(pEnd - pStart, &perm, pEnd - pStart, &iperm));
1568: for (p = pStart; p < pEnd; p++) {
1569: perm[p - pStart] = p;
1570: iperm[p - pStart] = p - pStart;
1571: }
1572: for (p = 0; p < pEnd - pStart;) {
1573: PetscInt point = perm[p];
1574: PetscInt parent;
1576: PetscCall(DMPlexGetTreeParent(dm, point, &parent, NULL));
1577: if (parent == point) {
1578: p++;
1579: } else {
1580: PetscInt size, closureSize, *closure = NULL, i;
1582: PetscCall(DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
1583: for (i = 0; i < closureSize; i++) {
1584: PetscInt q = closure[2 * i];
1585: if (iperm[q - pStart] > iperm[point - pStart]) {
1586: /* swap */
1587: perm[p] = q;
1588: perm[iperm[q - pStart]] = point;
1589: iperm[point - pStart] = iperm[q - pStart];
1590: iperm[q - pStart] = p;
1591: break;
1592: }
1593: }
1594: size = closureSize;
1595: PetscCall(DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
1596: if (i == size) p++;
1597: }
1598: }
1600: /* step 3: fill the constraint matrix */
1601: /* we are going to use a preorder progressive fill strategy. Mat doesn't
1602: * allow progressive fill without assembly, so we are going to set up the
1603: * values outside of the Mat first.
1604: */
1605: {
1606: PetscInt nRows, row, nnz;
1607: PetscBool done;
1608: PetscInt secStart, secEnd;
1609: const PetscInt *ia, *ja;
1610: PetscScalar *vals;
1612: PetscCall(PetscSectionGetChart(section, &secStart, &secEnd));
1613: PetscCall(MatGetRowIJ(cMat, 0, PETSC_FALSE, PETSC_FALSE, &nRows, &ia, &ja, &done));
1614: PetscCheck(done, PetscObjectComm((PetscObject)cMat), PETSC_ERR_PLIB, "Could not get RowIJ of constraint matrix");
1615: nnz = ia[nRows];
1616: /* malloc and then zero rows right before we fill them: this way valgrind
1617: * can tell if we are doing progressive fill in the wrong order */
1618: PetscCall(PetscMalloc1(nnz, &vals));
1619: for (p = 0; p < pEnd - pStart; p++) {
1620: PetscInt parent, childid, closureSize, *closure = NULL;
1621: PetscInt point = perm[p], pointDof;
1623: PetscCall(DMPlexGetTreeParent(dm, point, &parent, &childid));
1624: if ((point < conStart) || (point >= conEnd) || (parent == point)) continue;
1625: PetscCall(PetscSectionGetDof(conSec, point, &pointDof));
1626: if (!pointDof) continue;
1627: PetscCall(DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
1628: for (f = 0; f < maxFields; f++) {
1629: PetscInt cDof, cOff, numCols, numFillCols, i, r, matOffset, offset;
1630: PetscScalar *pointMat;
1631: const PetscInt **perms;
1632: const PetscScalar **flips;
1634: if (numFields) {
1635: PetscCall(PetscSectionGetFieldDof(conSec, point, f, &cDof));
1636: PetscCall(PetscSectionGetFieldOffset(conSec, point, f, &cOff));
1637: } else {
1638: PetscCall(PetscSectionGetDof(conSec, point, &cDof));
1639: PetscCall(PetscSectionGetOffset(conSec, point, &cOff));
1640: }
1641: if (!cDof) continue;
1642: if (numFields) PetscCall(PetscSectionGetFieldPointSyms(section, f, closureSize, closure, &perms, &flips));
1643: else PetscCall(PetscSectionGetPointSyms(section, closureSize, closure, &perms, &flips));
1645: /* make sure that every row for this point is the same size */
1646: if (PetscDefined(USE_DEBUG)) {
1647: for (r = 0; r < cDof; r++) {
1648: if (cDof > 1 && r) {
1649: 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]);
1650: }
1651: }
1652: }
1653: /* zero rows */
1654: for (i = ia[cOff]; i < ia[cOff + cDof]; i++) vals[i] = 0.;
1655: matOffset = ia[cOff];
1656: numFillCols = ia[cOff + 1] - matOffset;
1657: pointMat = refPointFieldMats[childid - pRefStart][f];
1658: numCols = refPointFieldN[childid - pRefStart][f];
1659: offset = 0;
1660: for (i = 0; i < closureSize; i++) {
1661: PetscInt q = closure[2 * i];
1662: PetscInt aDof, aOff, j, k, qConDof, qConOff;
1663: const PetscInt *perm = perms ? perms[i] : NULL;
1664: const PetscScalar *flip = flips ? flips[i] : NULL;
1666: qConDof = qConOff = 0;
1667: if (q < secStart || q >= secEnd) continue;
1668: if (numFields) {
1669: PetscCall(PetscSectionGetFieldDof(section, q, f, &aDof));
1670: PetscCall(PetscSectionGetFieldOffset(section, q, f, &aOff));
1671: if (q >= conStart && q < conEnd) {
1672: PetscCall(PetscSectionGetFieldDof(conSec, q, f, &qConDof));
1673: PetscCall(PetscSectionGetFieldOffset(conSec, q, f, &qConOff));
1674: }
1675: } else {
1676: PetscCall(PetscSectionGetDof(section, q, &aDof));
1677: PetscCall(PetscSectionGetOffset(section, q, &aOff));
1678: if (q >= conStart && q < conEnd) {
1679: PetscCall(PetscSectionGetDof(conSec, q, &qConDof));
1680: PetscCall(PetscSectionGetOffset(conSec, q, &qConOff));
1681: }
1682: }
1683: if (!aDof) continue;
1684: if (qConDof) {
1685: /* this point has anchors: its rows of the matrix should already
1686: * be filled, thanks to preordering */
1687: /* first multiply into pointWork, then set in matrix */
1688: PetscInt aMatOffset = ia[qConOff];
1689: PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset;
1690: for (r = 0; r < cDof; r++) {
1691: for (j = 0; j < aNumFillCols; j++) {
1692: PetscScalar inVal = 0;
1693: for (k = 0; k < aDof; k++) {
1694: PetscInt col = perm ? perm[k] : k;
1696: inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j] * (flip ? flip[col] : 1.);
1697: }
1698: pointWork[r * aNumFillCols + j] = inVal;
1699: }
1700: }
1701: /* assume that the columns are sorted, spend less time searching */
1702: for (j = 0, k = 0; j < aNumFillCols; j++) {
1703: PetscInt col = ja[aMatOffset + j];
1704: for (; k < numFillCols; k++) {
1705: if (ja[matOffset + k] == col) break;
1706: }
1707: PetscCheck(k != numFillCols, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No nonzero space for (%" PetscInt_FMT ", %" PetscInt_FMT ")", cOff, col);
1708: for (r = 0; r < cDof; r++) vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j];
1709: }
1710: } else {
1711: /* find where to put this portion of pointMat into the matrix */
1712: for (k = 0; k < numFillCols; k++) {
1713: if (ja[matOffset + k] == aOff) break;
1714: }
1715: PetscCheck(k != numFillCols, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No nonzero space for (%" PetscInt_FMT ", %" PetscInt_FMT ")", cOff, aOff);
1716: for (r = 0; r < cDof; r++) {
1717: for (j = 0; j < aDof; j++) {
1718: PetscInt col = perm ? perm[j] : j;
1720: vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + j] * (flip ? flip[col] : 1.);
1721: }
1722: }
1723: }
1724: offset += aDof;
1725: }
1726: if (numFields) {
1727: PetscCall(PetscSectionRestoreFieldPointSyms(section, f, closureSize, closure, &perms, &flips));
1728: } else {
1729: PetscCall(PetscSectionRestorePointSyms(section, closureSize, closure, &perms, &flips));
1730: }
1731: }
1732: PetscCall(DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
1733: }
1734: for (row = 0; row < nRows; row++) PetscCall(MatSetValues(cMat, 1, &row, ia[row + 1] - ia[row], &ja[ia[row]], &vals[ia[row]], INSERT_VALUES));
1735: PetscCall(MatRestoreRowIJ(cMat, 0, PETSC_FALSE, PETSC_FALSE, &nRows, &ia, &ja, &done));
1736: PetscCheck(done, PetscObjectComm((PetscObject)cMat), PETSC_ERR_PLIB, "Could not restore RowIJ of constraint matrix");
1737: PetscCall(MatAssemblyBegin(cMat, MAT_FINAL_ASSEMBLY));
1738: PetscCall(MatAssemblyEnd(cMat, MAT_FINAL_ASSEMBLY));
1739: PetscCall(PetscFree(vals));
1740: }
1742: /* clean up */
1743: PetscCall(ISRestoreIndices(anIS, &anchors));
1744: PetscCall(PetscFree2(perm, iperm));
1745: PetscCall(PetscFree(pointWork));
1746: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
1747: PetscFunctionReturn(PETSC_SUCCESS);
1748: }
1750: /*@
1751: DMPlexTreeRefineCell - Refine a single cell on rank 0 using the `DM`'s reference tree, producing a non-conforming mesh
1753: Collective
1755: Input Parameters:
1756: + dm - The `DM` with an attached reference tree (see `DMPlexSetReferenceTree()`)
1757: - cell - The cell to be refined
1759: Output Parameter:
1760: . ncdm - A new `DM` in which `cell` has been split according to the reference tree
1762: Level: developer
1764: Note:
1765: This routine is intended for testing and demonstration; it produces one example of a non-conforming
1766: mesh but is not a general local-refinement facility. Only rank 0 performs the refinement.
1768: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetReferenceTree()`, `DMPlexGetReferenceTree()`, `DMPlexSetTree()`
1769: @*/
1770: /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of
1771: * a non-conforming mesh. Local refinement comes later */
1772: PetscErrorCode DMPlexTreeRefineCell(DM dm, PetscInt cell, DM *ncdm)
1773: {
1774: DM K;
1775: PetscMPIInt rank;
1776: PetscInt dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd;
1777: PetscInt numNewCones, *newConeSizes, *newCones, *newOrientations;
1778: PetscInt *Kembedding;
1779: PetscInt *cellClosure = NULL, nc;
1780: PetscScalar *newVertexCoords;
1781: PetscInt numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset;
1782: PetscSection parentSection;
1784: PetscFunctionBegin;
1785: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1786: PetscCall(DMGetDimension(dm, &dim));
1787: PetscCall(DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm));
1788: PetscCall(DMSetDimension(*ncdm, dim));
1790: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1791: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &parentSection));
1792: PetscCall(DMPlexGetReferenceTree(dm, &K));
1793: PetscCall(DMGetCoordinatesLocalSetUp(dm));
1794: if (rank == 0) {
1795: /* compute the new charts */
1796: PetscCall(PetscMalloc5(dim + 1, &pNewCount, dim + 1, &pNewStart, dim + 1, &pNewEnd, dim + 1, &pOldStart, dim + 1, &pOldEnd));
1797: offset = 0;
1798: for (d = 0; d <= dim; d++) {
1799: PetscInt pOldCount, kStart, kEnd, k;
1801: pNewStart[d] = offset;
1802: PetscCall(DMPlexGetHeightStratum(dm, d, &pOldStart[d], &pOldEnd[d]));
1803: PetscCall(DMPlexGetHeightStratum(K, d, &kStart, &kEnd));
1804: pOldCount = pOldEnd[d] - pOldStart[d];
1805: /* adding the new points */
1806: pNewCount[d] = pOldCount + kEnd - kStart;
1807: if (!d) {
1808: /* removing the cell */
1809: pNewCount[d]--;
1810: }
1811: for (k = kStart; k < kEnd; k++) {
1812: PetscInt parent;
1813: PetscCall(DMPlexGetTreeParent(K, k, &parent, NULL));
1814: if (parent == k) {
1815: /* avoid double counting points that won't actually be new */
1816: pNewCount[d]--;
1817: }
1818: }
1819: pNewEnd[d] = pNewStart[d] + pNewCount[d];
1820: offset = pNewEnd[d];
1821: }
1822: 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]);
1823: /* get the current closure of the cell that we are removing */
1824: PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &nc, &cellClosure));
1826: PetscCall(PetscMalloc1(pNewEnd[dim], &newConeSizes));
1827: {
1828: DMPolytopeType pct, qct;
1829: PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j;
1831: PetscCall(DMPlexGetChart(K, &kStart, &kEnd));
1832: PetscCall(PetscMalloc4(kEnd - kStart, &Kembedding, kEnd - kStart, &perm, kEnd - kStart, &iperm, kEnd - kStart, &preOrient));
1834: for (k = kStart; k < kEnd; k++) {
1835: perm[k - kStart] = k;
1836: iperm[k - kStart] = k - kStart;
1837: preOrient[k - kStart] = 0;
1838: }
1840: PetscCall(DMPlexGetTransitiveClosure(K, 0, PETSC_TRUE, &closureSizeK, &closureK));
1841: for (j = 1; j < closureSizeK; j++) {
1842: PetscInt parentOrientA = closureK[2 * j + 1];
1843: PetscInt parentOrientB = cellClosure[2 * j + 1];
1844: PetscInt p, q;
1846: p = closureK[2 * j];
1847: q = cellClosure[2 * j];
1848: PetscCall(DMPlexGetCellType(K, p, &pct));
1849: PetscCall(DMPlexGetCellType(dm, q, &qct));
1850: for (d = 0; d <= dim; d++) {
1851: if (q >= pOldStart[d] && q < pOldEnd[d]) Kembedding[p] = (q - pOldStart[d]) + pNewStart[d];
1852: }
1853: parentOrientA = DMPolytopeConvertNewOrientation_Internal(pct, parentOrientA);
1854: parentOrientB = DMPolytopeConvertNewOrientation_Internal(qct, parentOrientB);
1855: if (parentOrientA != parentOrientB) {
1856: PetscInt numChildren, i;
1857: const PetscInt *children;
1859: PetscCall(DMPlexGetTreeChildren(K, p, &numChildren, &children));
1860: for (i = 0; i < numChildren; i++) {
1861: PetscInt kPerm, oPerm;
1863: k = children[i];
1864: PetscCall(DMPlexReferenceTreeGetChildSymmetry(K, p, parentOrientA, 0, k, parentOrientB, &oPerm, &kPerm));
1865: /* perm = what refTree position I'm in */
1866: perm[kPerm - kStart] = k;
1867: /* iperm = who is at this position */
1868: iperm[k - kStart] = kPerm - kStart;
1869: preOrient[kPerm - kStart] = oPerm;
1870: }
1871: }
1872: }
1873: PetscCall(DMPlexRestoreTransitiveClosure(K, 0, PETSC_TRUE, &closureSizeK, &closureK));
1874: }
1875: PetscCall(PetscSectionSetChart(parentSection, 0, pNewEnd[dim]));
1876: offset = 0;
1877: numNewCones = 0;
1878: for (d = 0; d <= dim; d++) {
1879: PetscInt kStart, kEnd, k;
1880: PetscInt p;
1881: PetscInt size;
1883: for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1884: /* skip cell 0 */
1885: if (p == cell) continue;
1886: /* old cones to new cones */
1887: PetscCall(DMPlexGetConeSize(dm, p, &size));
1888: newConeSizes[offset++] = size;
1889: numNewCones += size;
1890: }
1892: PetscCall(DMPlexGetHeightStratum(K, d, &kStart, &kEnd));
1893: for (k = kStart; k < kEnd; k++) {
1894: PetscInt kParent;
1896: PetscCall(DMPlexGetTreeParent(K, k, &kParent, NULL));
1897: if (kParent != k) {
1898: Kembedding[k] = offset;
1899: PetscCall(DMPlexGetConeSize(K, k, &size));
1900: newConeSizes[offset++] = size;
1901: numNewCones += size;
1902: if (kParent != 0) PetscCall(PetscSectionSetDof(parentSection, Kembedding[k], 1));
1903: }
1904: }
1905: }
1907: PetscCall(PetscSectionSetUp(parentSection));
1908: PetscCall(PetscSectionGetStorageSize(parentSection, &numPointsWithParents));
1909: PetscCall(PetscMalloc2(numNewCones, &newCones, numNewCones, &newOrientations));
1910: PetscCall(PetscMalloc2(numPointsWithParents, &parents, numPointsWithParents, &childIDs));
1912: /* fill new cones */
1913: offset = 0;
1914: for (d = 0; d <= dim; d++) {
1915: PetscInt kStart, kEnd, k, l;
1916: PetscInt p;
1917: PetscInt size;
1918: const PetscInt *cone, *orientation;
1920: for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1921: /* skip cell 0 */
1922: if (p == cell) continue;
1923: /* old cones to new cones */
1924: PetscCall(DMPlexGetConeSize(dm, p, &size));
1925: PetscCall(DMPlexGetCone(dm, p, &cone));
1926: PetscCall(DMPlexGetConeOrientation(dm, p, &orientation));
1927: for (l = 0; l < size; l++) {
1928: newCones[offset] = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1];
1929: newOrientations[offset++] = orientation[l];
1930: }
1931: }
1933: PetscCall(DMPlexGetHeightStratum(K, d, &kStart, &kEnd));
1934: for (k = kStart; k < kEnd; k++) {
1935: PetscInt kPerm = perm[k], kParent;
1936: PetscInt preO = preOrient[k];
1938: PetscCall(DMPlexGetTreeParent(K, k, &kParent, NULL));
1939: if (kParent != k) {
1940: /* embed new cones */
1941: PetscCall(DMPlexGetConeSize(K, k, &size));
1942: PetscCall(DMPlexGetCone(K, kPerm, &cone));
1943: PetscCall(DMPlexGetConeOrientation(K, kPerm, &orientation));
1944: for (l = 0; l < size; l++) {
1945: PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size - (preO + 1) - l) % size);
1946: PetscInt newO, lSize, oTrue;
1947: DMPolytopeType ct = DM_NUM_POLYTOPES;
1949: q = iperm[cone[m]];
1950: newCones[offset] = Kembedding[q];
1951: PetscCall(DMPlexGetConeSize(K, q, &lSize));
1952: if (lSize == 2) ct = DM_POLYTOPE_SEGMENT;
1953: else if (lSize == 4) ct = DM_POLYTOPE_QUADRILATERAL;
1954: oTrue = DMPolytopeConvertNewOrientation_Internal(ct, orientation[m]);
1955: oTrue = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2);
1956: newO = DihedralCompose(lSize, oTrue, preOrient[q]);
1957: newOrientations[offset++] = DMPolytopeConvertOldOrientation_Internal(ct, newO);
1958: }
1959: if (kParent != 0) {
1960: PetscInt newPoint = Kembedding[kParent];
1961: PetscCall(PetscSectionGetOffset(parentSection, Kembedding[k], &pOffset));
1962: parents[pOffset] = newPoint;
1963: childIDs[pOffset] = k;
1964: }
1965: }
1966: }
1967: }
1969: PetscCall(PetscMalloc1(dim * (pNewEnd[dim] - pNewStart[dim]), &newVertexCoords));
1971: /* fill coordinates */
1972: offset = 0;
1973: {
1974: PetscInt kStart, kEnd, l;
1975: PetscSection vSection;
1976: PetscInt v;
1977: Vec coords;
1978: PetscScalar *coordvals;
1979: PetscInt dof, off;
1980: PetscReal v0[3], J[9], detJ;
1982: if (PetscDefined(USE_DEBUG)) {
1983: PetscInt k;
1984: PetscCall(DMPlexGetHeightStratum(K, 0, &kStart, &kEnd));
1985: for (k = kStart; k < kEnd; k++) {
1986: PetscCall(DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ));
1987: PetscCheck(detJ > 0., PETSC_COMM_SELF, PETSC_ERR_PLIB, "reference tree cell %" PetscInt_FMT " has bad determinant", k);
1988: }
1989: }
1990: PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ));
1991: PetscCall(DMGetCoordinateSection(dm, &vSection));
1992: PetscCall(DMGetCoordinatesLocal(dm, &coords));
1993: PetscCall(VecGetArray(coords, &coordvals));
1994: for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) {
1995: PetscCall(PetscSectionGetDof(vSection, v, &dof));
1996: PetscCall(PetscSectionGetOffset(vSection, v, &off));
1997: for (l = 0; l < dof; l++) newVertexCoords[offset++] = coordvals[off + l];
1998: }
1999: PetscCall(VecRestoreArray(coords, &coordvals));
2001: PetscCall(DMGetCoordinateSection(K, &vSection));
2002: PetscCall(DMGetCoordinatesLocal(K, &coords));
2003: PetscCall(VecGetArray(coords, &coordvals));
2004: PetscCall(DMPlexGetDepthStratum(K, 0, &kStart, &kEnd));
2005: for (v = kStart; v < kEnd; v++) {
2006: PetscReal coord[3], newCoord[3];
2007: PetscInt vPerm = perm[v];
2008: PetscInt kParent;
2009: const PetscReal xi0[3] = {-1., -1., -1.};
2011: PetscCall(DMPlexGetTreeParent(K, v, &kParent, NULL));
2012: if (kParent != v) {
2013: /* this is a new vertex */
2014: PetscCall(PetscSectionGetOffset(vSection, vPerm, &off));
2015: for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off + l]);
2016: CoordinatesRefToReal(dim, dim, xi0, v0, J, coord, newCoord);
2017: for (l = 0; l < dim; ++l) newVertexCoords[offset + l] = newCoord[l];
2018: offset += dim;
2019: }
2020: }
2021: PetscCall(VecRestoreArray(coords, &coordvals));
2022: }
2024: /* need to reverse the order of pNewCount: vertices first, cells last */
2025: for (d = 0; d < (dim + 1) / 2; d++) {
2026: PetscInt tmp;
2028: tmp = pNewCount[d];
2029: pNewCount[d] = pNewCount[dim - d];
2030: pNewCount[dim - d] = tmp;
2031: }
2033: PetscCall(DMPlexCreateFromDAG(*ncdm, dim, pNewCount, newConeSizes, newCones, newOrientations, newVertexCoords));
2034: PetscCall(DMPlexSetReferenceTree(*ncdm, K));
2035: PetscCall(DMPlexSetTree(*ncdm, parentSection, parents, childIDs));
2037: /* clean up */
2038: PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &nc, &cellClosure));
2039: PetscCall(PetscFree5(pNewCount, pNewStart, pNewEnd, pOldStart, pOldEnd));
2040: PetscCall(PetscFree(newConeSizes));
2041: PetscCall(PetscFree2(newCones, newOrientations));
2042: PetscCall(PetscFree(newVertexCoords));
2043: PetscCall(PetscFree2(parents, childIDs));
2044: PetscCall(PetscFree4(Kembedding, perm, iperm, preOrient));
2045: } else {
2046: PetscInt p, counts[4];
2047: PetscInt *coneSizes, *cones, *orientations;
2048: Vec coordVec;
2049: PetscScalar *coords;
2051: for (d = 0; d <= dim; d++) {
2052: PetscInt dStart, dEnd;
2054: PetscCall(DMPlexGetDepthStratum(dm, d, &dStart, &dEnd));
2055: counts[d] = dEnd - dStart;
2056: }
2057: PetscCall(PetscMalloc1(pEnd - pStart, &coneSizes));
2058: for (p = pStart; p < pEnd; p++) PetscCall(DMPlexGetConeSize(dm, p, &coneSizes[p - pStart]));
2059: PetscCall(DMPlexGetCones(dm, &cones));
2060: PetscCall(DMPlexGetConeOrientations(dm, &orientations));
2061: PetscCall(DMGetCoordinatesLocal(dm, &coordVec));
2062: PetscCall(VecGetArray(coordVec, &coords));
2064: PetscCall(PetscSectionSetChart(parentSection, pStart, pEnd));
2065: PetscCall(PetscSectionSetUp(parentSection));
2066: PetscCall(DMPlexCreateFromDAG(*ncdm, dim, counts, coneSizes, cones, orientations, NULL));
2067: PetscCall(DMPlexSetReferenceTree(*ncdm, K));
2068: PetscCall(DMPlexSetTree(*ncdm, parentSection, NULL, NULL));
2069: PetscCall(VecRestoreArray(coordVec, &coords));
2070: }
2071: PetscCall(PetscSectionDestroy(&parentSection));
2072: PetscFunctionReturn(PETSC_SUCCESS);
2073: }
2075: PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
2076: {
2077: PetscSF coarseToFineEmbedded;
2078: PetscSection globalCoarse, globalFine;
2079: PetscSection localCoarse, localFine;
2080: PetscSection aSec, cSec;
2081: PetscSection rootIndicesSec, rootMatricesSec;
2082: PetscSection leafIndicesSec, leafMatricesSec;
2083: PetscInt *rootIndices, *leafIndices;
2084: PetscScalar *rootMatrices, *leafMatrices;
2085: IS aIS;
2086: const PetscInt *anchors;
2087: Mat cMat;
2088: PetscInt numFields, maxFields;
2089: PetscInt pStartC, pEndC, pStartF, pEndF, p;
2090: PetscInt aStart, aEnd, cStart, cEnd;
2091: PetscInt *maxChildIds;
2092: PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
2093: const PetscInt ***perms;
2094: const PetscScalar ***flips;
2096: PetscFunctionBegin;
2097: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
2098: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
2099: PetscCall(DMGetGlobalSection(fine, &globalFine));
2100: { /* winnow fine points that don't have global dofs out of the sf */
2101: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, nleaves, l;
2102: const PetscInt *leaves;
2104: PetscCall(PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL));
2105: for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
2106: p = leaves ? leaves[l] : l;
2107: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
2108: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
2109: if ((dof - cdof) > 0) numPointsWithDofs++;
2110: }
2111: PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
2112: for (l = 0, offset = 0; l < nleaves; l++) {
2113: p = leaves ? leaves[l] : l;
2114: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
2115: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
2116: if ((dof - cdof) > 0) pointsWithDofs[offset++] = l;
2117: }
2118: PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded));
2119: PetscCall(PetscFree(pointsWithDofs));
2120: }
2121: /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
2122: PetscCall(PetscMalloc1(pEndC - pStartC, &maxChildIds));
2123: for (p = pStartC; p < pEndC; p++) maxChildIds[p - pStartC] = -2;
2124: PetscCall(PetscSFReduceBegin(coarseToFineEmbedded, MPIU_INT, childIds, maxChildIds, MPI_MAX));
2125: PetscCall(PetscSFReduceEnd(coarseToFineEmbedded, MPIU_INT, childIds, maxChildIds, MPI_MAX));
2127: PetscCall(DMGetLocalSection(coarse, &localCoarse));
2128: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
2130: PetscCall(DMPlexGetAnchors(coarse, &aSec, &aIS));
2131: PetscCall(ISGetIndices(aIS, &anchors));
2132: PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd));
2134: PetscCall(DMGetDefaultConstraints(coarse, &cSec, &cMat, NULL));
2135: PetscCall(PetscSectionGetChart(cSec, &cStart, &cEnd));
2137: /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
2138: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootIndicesSec));
2139: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootMatricesSec));
2140: PetscCall(PetscSectionSetChart(rootIndicesSec, pStartC, pEndC));
2141: PetscCall(PetscSectionSetChart(rootMatricesSec, pStartC, pEndC));
2142: PetscCall(PetscSectionGetNumFields(localCoarse, &numFields));
2143: maxFields = PetscMax(1, numFields);
2144: PetscCall(PetscMalloc7(maxFields + 1, &offsets, maxFields + 1, &offsetsCopy, maxFields + 1, &newOffsets, maxFields + 1, &newOffsetsCopy, maxFields + 1, &rowOffsets, maxFields + 1, &numD, maxFields + 1, &numO));
2145: PetscCall(PetscMalloc2(maxFields + 1, (PetscInt ****)&perms, maxFields + 1, (PetscScalar ****)&flips));
2146: PetscCall(PetscMemzero((void *)perms, (maxFields + 1) * sizeof(const PetscInt **)));
2147: PetscCall(PetscMemzero((void *)flips, (maxFields + 1) * sizeof(const PetscScalar **)));
2149: for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
2150: PetscInt dof, matSize = 0;
2151: PetscInt aDof = 0;
2152: PetscInt cDof = 0;
2153: PetscInt maxChildId = maxChildIds[p - pStartC];
2154: PetscInt numRowIndices = 0;
2155: PetscInt numColIndices = 0;
2156: PetscInt f;
2158: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
2159: if (dof < 0) dof = -(dof + 1);
2160: if (p >= aStart && p < aEnd) PetscCall(PetscSectionGetDof(aSec, p, &aDof));
2161: if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(cSec, p, &cDof));
2162: for (f = 0; f <= numFields; f++) offsets[f] = 0;
2163: for (f = 0; f <= numFields; f++) newOffsets[f] = 0;
2164: if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
2165: PetscInt *closure = NULL, closureSize, cl;
2167: PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2168: for (cl = 0; cl < closureSize; cl++) { /* get the closure */
2169: PetscInt c = closure[2 * cl], clDof;
2171: PetscCall(PetscSectionGetDof(localCoarse, c, &clDof));
2172: numRowIndices += clDof;
2173: for (f = 0; f < numFields; f++) {
2174: PetscCall(PetscSectionGetFieldDof(localCoarse, c, f, &clDof));
2175: offsets[f + 1] += clDof;
2176: }
2177: }
2178: for (f = 0; f < numFields; f++) {
2179: offsets[f + 1] += offsets[f];
2180: newOffsets[f + 1] = offsets[f + 1];
2181: }
2182: /* get the number of indices needed and their field offsets */
2183: PetscCall(DMPlexAnchorsModifyMat(coarse, localCoarse, closureSize, numRowIndices, closure, NULL, NULL, NULL, &numColIndices, NULL, NULL, newOffsets, PETSC_FALSE));
2184: PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2185: if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */
2186: numColIndices = numRowIndices;
2187: matSize = 0;
2188: } else if (numFields) { /* we send one submat for each field: sum their sizes */
2189: matSize = 0;
2190: for (f = 0; f < numFields; f++) {
2191: PetscInt numRow, numCol;
2193: numRow = offsets[f + 1] - offsets[f];
2194: numCol = newOffsets[f + 1] - newOffsets[f];
2195: matSize += numRow * numCol;
2196: }
2197: } else {
2198: matSize = numRowIndices * numColIndices;
2199: }
2200: } else if (maxChildId == -1) {
2201: if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */
2202: PetscInt aOff, a;
2204: PetscCall(PetscSectionGetOffset(aSec, p, &aOff));
2205: for (f = 0; f < numFields; f++) {
2206: PetscInt fDof;
2208: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
2209: offsets[f + 1] = fDof;
2210: }
2211: for (a = 0; a < aDof; a++) {
2212: PetscInt anchor = anchors[a + aOff], aLocalDof;
2214: PetscCall(PetscSectionGetDof(localCoarse, anchor, &aLocalDof));
2215: numColIndices += aLocalDof;
2216: for (f = 0; f < numFields; f++) {
2217: PetscInt fDof;
2219: PetscCall(PetscSectionGetFieldDof(localCoarse, anchor, f, &fDof));
2220: newOffsets[f + 1] += fDof;
2221: }
2222: }
2223: if (numFields) {
2224: matSize = 0;
2225: for (f = 0; f < numFields; f++) matSize += offsets[f + 1] * newOffsets[f + 1];
2226: } else {
2227: matSize = numColIndices * dof;
2228: }
2229: } else { /* no children, and no constraints on dofs: just get the global indices */
2230: numColIndices = dof;
2231: matSize = 0;
2232: }
2233: }
2234: /* we will pack the column indices with the field offsets */
2235: PetscCall(PetscSectionSetDof(rootIndicesSec, p, numColIndices ? numColIndices + 2 * numFields : 0));
2236: PetscCall(PetscSectionSetDof(rootMatricesSec, p, matSize));
2237: }
2238: PetscCall(PetscSectionSetUp(rootIndicesSec));
2239: PetscCall(PetscSectionSetUp(rootMatricesSec));
2240: {
2241: PetscInt numRootIndices, numRootMatrices;
2243: PetscCall(PetscSectionGetStorageSize(rootIndicesSec, &numRootIndices));
2244: PetscCall(PetscSectionGetStorageSize(rootMatricesSec, &numRootMatrices));
2245: PetscCall(PetscMalloc2(numRootIndices, &rootIndices, numRootMatrices, &rootMatrices));
2246: for (p = pStartC; p < pEndC; p++) {
2247: PetscInt numRowIndices = 0, numColIndices, matSize, dof;
2248: PetscInt pIndOff, pMatOff, f;
2249: PetscInt *pInd;
2250: PetscInt maxChildId = maxChildIds[p - pStartC];
2251: PetscScalar *pMat = NULL;
2253: PetscCall(PetscSectionGetDof(rootIndicesSec, p, &numColIndices));
2254: if (!numColIndices) continue;
2255: for (f = 0; f <= numFields; f++) {
2256: offsets[f] = 0;
2257: newOffsets[f] = 0;
2258: offsetsCopy[f] = 0;
2259: newOffsetsCopy[f] = 0;
2260: }
2261: numColIndices -= 2 * numFields;
2262: PetscCall(PetscSectionGetOffset(rootIndicesSec, p, &pIndOff));
2263: pInd = &rootIndices[pIndOff];
2264: PetscCall(PetscSectionGetDof(rootMatricesSec, p, &matSize));
2265: if (matSize) {
2266: PetscCall(PetscSectionGetOffset(rootMatricesSec, p, &pMatOff));
2267: pMat = &rootMatrices[pMatOff];
2268: }
2269: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
2270: if (dof < 0) dof = -(dof + 1);
2271: if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
2272: PetscInt i, j;
2274: if (matSize == 0) { /* don't need to calculate the mat, just the indices */
2275: PetscInt numIndices, *indices;
2276: PetscCall(DMPlexGetClosureIndices(coarse, localCoarse, globalCoarse, p, PETSC_TRUE, &numIndices, &indices, offsets, NULL));
2277: PetscCheck(numIndices == numColIndices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "mismatching constraint indices calculations");
2278: for (i = 0; i < numColIndices; i++) pInd[i] = indices[i];
2279: for (i = 0; i < numFields; i++) {
2280: pInd[numColIndices + i] = offsets[i + 1];
2281: pInd[numColIndices + numFields + i] = offsets[i + 1];
2282: }
2283: PetscCall(DMPlexRestoreClosureIndices(coarse, localCoarse, globalCoarse, p, PETSC_TRUE, &numIndices, &indices, offsets, NULL));
2284: } else {
2285: PetscInt closureSize, *closure = NULL, cl;
2286: PetscScalar *pMatIn, *pMatModified;
2287: PetscInt numPoints, *points;
2289: {
2290: PetscInt *closure = NULL, closureSize, cl;
2292: PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2293: for (cl = 0; cl < closureSize; cl++) { /* get the closure */
2294: PetscInt c = closure[2 * cl], clDof;
2296: PetscCall(PetscSectionGetDof(localCoarse, c, &clDof));
2297: numRowIndices += clDof;
2298: }
2299: PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2300: }
2302: PetscCall(DMGetWorkArray(coarse, numRowIndices * numRowIndices, MPIU_SCALAR, &pMatIn));
2303: for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */
2304: for (j = 0; j < numRowIndices; j++) pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.;
2305: }
2306: PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2307: for (f = 0; f < maxFields; f++) {
2308: if (numFields) PetscCall(PetscSectionGetFieldPointSyms(localCoarse, f, closureSize, closure, &perms[f], &flips[f]));
2309: else PetscCall(PetscSectionGetPointSyms(localCoarse, closureSize, closure, &perms[f], &flips[f]));
2310: }
2311: if (numFields) {
2312: for (cl = 0; cl < closureSize; cl++) {
2313: PetscInt c = closure[2 * cl];
2315: for (f = 0; f < numFields; f++) {
2316: PetscInt fDof;
2318: PetscCall(PetscSectionGetFieldDof(localCoarse, c, f, &fDof));
2319: offsets[f + 1] += fDof;
2320: }
2321: }
2322: for (f = 0; f < numFields; f++) {
2323: offsets[f + 1] += offsets[f];
2324: newOffsets[f + 1] = offsets[f + 1];
2325: }
2326: }
2327: /* TODO : flips here ? */
2328: /* apply hanging node constraints on the right, get the new points and the new offsets */
2329: PetscCall(DMPlexAnchorsModifyMat(coarse, localCoarse, closureSize, numRowIndices, closure, perms, pMatIn, &numPoints, NULL, &points, &pMatModified, newOffsets, PETSC_FALSE));
2330: for (f = 0; f < maxFields; f++) {
2331: if (numFields) PetscCall(PetscSectionRestoreFieldPointSyms(localCoarse, f, closureSize, closure, &perms[f], &flips[f]));
2332: else PetscCall(PetscSectionRestorePointSyms(localCoarse, closureSize, closure, &perms[f], &flips[f]));
2333: }
2334: for (f = 0; f < maxFields; f++) {
2335: if (numFields) PetscCall(PetscSectionGetFieldPointSyms(localCoarse, f, numPoints, points, &perms[f], &flips[f]));
2336: else PetscCall(PetscSectionGetPointSyms(localCoarse, numPoints, points, &perms[f], &flips[f]));
2337: }
2338: if (!numFields) {
2339: for (i = 0; i < numRowIndices * numColIndices; i++) pMat[i] = pMatModified[i];
2340: } else {
2341: PetscInt i, j, count;
2342: for (f = 0, count = 0; f < numFields; f++) {
2343: for (i = offsets[f]; i < offsets[f + 1]; i++) {
2344: for (j = newOffsets[f]; j < newOffsets[f + 1]; j++, count++) pMat[count] = pMatModified[i * numColIndices + j];
2345: }
2346: }
2347: }
2348: PetscCall(DMRestoreWorkArray(coarse, numRowIndices * numColIndices, MPIU_SCALAR, &pMatModified));
2349: PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2350: PetscCall(DMRestoreWorkArray(coarse, numRowIndices * numColIndices, MPIU_SCALAR, &pMatIn));
2351: if (numFields) {
2352: for (f = 0; f < numFields; f++) {
2353: pInd[numColIndices + f] = offsets[f + 1];
2354: pInd[numColIndices + numFields + f] = newOffsets[f + 1];
2355: }
2356: for (cl = 0; cl < numPoints; cl++) {
2357: PetscInt globalOff, c = points[2 * cl];
2358: PetscCall(PetscSectionGetOffset(globalCoarse, c, &globalOff));
2359: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff + 1) : globalOff, newOffsets, PETSC_FALSE, perms, cl, NULL, pInd));
2360: }
2361: } else {
2362: for (cl = 0; cl < numPoints; cl++) {
2363: PetscInt c = points[2 * cl], globalOff;
2364: const PetscInt *perm = perms[0] ? perms[0][cl] : NULL;
2366: PetscCall(PetscSectionGetOffset(globalCoarse, c, &globalOff));
2367: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff + 1) : globalOff, newOffsets, PETSC_FALSE, perm, NULL, pInd));
2368: }
2369: }
2370: for (f = 0; f < maxFields; f++) {
2371: if (numFields) PetscCall(PetscSectionRestoreFieldPointSyms(localCoarse, f, numPoints, points, &perms[f], &flips[f]));
2372: else PetscCall(PetscSectionRestorePointSyms(localCoarse, numPoints, points, &perms[f], &flips[f]));
2373: }
2374: PetscCall(DMRestoreWorkArray(coarse, numPoints, MPIU_SCALAR, &points));
2375: }
2376: } else if (matSize) {
2377: PetscInt cOff;
2378: PetscInt *rowIndices, *colIndices, a, aDof = 0, aOff;
2380: numRowIndices = dof;
2381: PetscCall(DMGetWorkArray(coarse, numRowIndices, MPIU_INT, &rowIndices));
2382: PetscCall(DMGetWorkArray(coarse, numColIndices, MPIU_INT, &colIndices));
2383: PetscCall(PetscSectionGetOffset(cSec, p, &cOff));
2384: PetscCall(PetscSectionGetDof(aSec, p, &aDof));
2385: PetscCall(PetscSectionGetOffset(aSec, p, &aOff));
2386: if (numFields) {
2387: for (f = 0; f < numFields; f++) {
2388: PetscInt fDof;
2390: PetscCall(PetscSectionGetFieldDof(cSec, p, f, &fDof));
2391: offsets[f + 1] = fDof;
2392: for (a = 0; a < aDof; a++) {
2393: PetscInt anchor = anchors[a + aOff];
2394: PetscCall(PetscSectionGetFieldDof(localCoarse, anchor, f, &fDof));
2395: newOffsets[f + 1] += fDof;
2396: }
2397: }
2398: for (f = 0; f < numFields; f++) {
2399: offsets[f + 1] += offsets[f];
2400: offsetsCopy[f + 1] = offsets[f + 1];
2401: newOffsets[f + 1] += newOffsets[f];
2402: newOffsetsCopy[f + 1] = newOffsets[f + 1];
2403: }
2404: PetscCall(DMPlexGetIndicesPointFields_Internal(cSec, PETSC_TRUE, p, cOff, offsetsCopy, PETSC_TRUE, NULL, -1, NULL, rowIndices));
2405: for (a = 0; a < aDof; a++) {
2406: PetscInt anchor = anchors[a + aOff], lOff;
2407: PetscCall(PetscSectionGetOffset(localCoarse, anchor, &lOff));
2408: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_TRUE, anchor, lOff, newOffsetsCopy, PETSC_TRUE, NULL, -1, NULL, colIndices));
2409: }
2410: } else {
2411: PetscCall(DMPlexGetIndicesPoint_Internal(cSec, PETSC_TRUE, p, cOff, offsetsCopy, PETSC_TRUE, NULL, NULL, rowIndices));
2412: for (a = 0; a < aDof; a++) {
2413: PetscInt anchor = anchors[a + aOff], lOff;
2414: PetscCall(PetscSectionGetOffset(localCoarse, anchor, &lOff));
2415: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_TRUE, anchor, lOff, newOffsetsCopy, PETSC_TRUE, NULL, NULL, colIndices));
2416: }
2417: }
2418: if (numFields) {
2419: PetscInt count, a;
2421: for (f = 0, count = 0; f < numFields; f++) {
2422: PetscInt iSize = offsets[f + 1] - offsets[f];
2423: PetscInt jSize = newOffsets[f + 1] - newOffsets[f];
2424: PetscCall(MatGetValues(cMat, iSize, &rowIndices[offsets[f]], jSize, &colIndices[newOffsets[f]], &pMat[count]));
2425: count += iSize * jSize;
2426: pInd[numColIndices + f] = offsets[f + 1];
2427: pInd[numColIndices + numFields + f] = newOffsets[f + 1];
2428: }
2429: for (a = 0; a < aDof; a++) {
2430: PetscInt anchor = anchors[a + aOff];
2431: PetscInt gOff;
2432: PetscCall(PetscSectionGetOffset(globalCoarse, anchor, &gOff));
2433: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, anchor, gOff < 0 ? -(gOff + 1) : gOff, newOffsets, PETSC_FALSE, NULL, -1, NULL, pInd));
2434: }
2435: } else {
2436: PetscInt a;
2437: PetscCall(MatGetValues(cMat, numRowIndices, rowIndices, numColIndices, colIndices, pMat));
2438: for (a = 0; a < aDof; a++) {
2439: PetscInt anchor = anchors[a + aOff];
2440: PetscInt gOff;
2441: PetscCall(PetscSectionGetOffset(globalCoarse, anchor, &gOff));
2442: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, anchor, gOff < 0 ? -(gOff + 1) : gOff, newOffsets, PETSC_FALSE, NULL, NULL, pInd));
2443: }
2444: }
2445: PetscCall(DMRestoreWorkArray(coarse, numColIndices, MPIU_INT, &colIndices));
2446: PetscCall(DMRestoreWorkArray(coarse, numRowIndices, MPIU_INT, &rowIndices));
2447: } else {
2448: PetscInt gOff;
2450: PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));
2451: if (numFields) {
2452: for (f = 0; f < numFields; f++) {
2453: PetscInt fDof;
2454: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
2455: offsets[f + 1] = fDof + offsets[f];
2456: }
2457: for (f = 0; f < numFields; f++) {
2458: pInd[numColIndices + f] = offsets[f + 1];
2459: pInd[numColIndices + numFields + f] = offsets[f + 1];
2460: }
2461: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, -1, NULL, pInd));
2462: } else {
2463: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, NULL, pInd));
2464: }
2465: }
2466: }
2467: PetscCall(PetscFree(maxChildIds));
2468: }
2469: {
2470: PetscSF indicesSF, matricesSF;
2471: PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices;
2473: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafIndicesSec));
2474: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafMatricesSec));
2475: PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootIndicesSec, &remoteOffsetsIndices, leafIndicesSec));
2476: PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootMatricesSec, &remoteOffsetsMatrices, leafMatricesSec));
2477: PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootIndicesSec, remoteOffsetsIndices, leafIndicesSec, &indicesSF));
2478: PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootMatricesSec, remoteOffsetsMatrices, leafMatricesSec, &matricesSF));
2479: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
2480: PetscCall(PetscFree(remoteOffsetsIndices));
2481: PetscCall(PetscFree(remoteOffsetsMatrices));
2482: PetscCall(PetscSectionGetStorageSize(leafIndicesSec, &numLeafIndices));
2483: PetscCall(PetscSectionGetStorageSize(leafMatricesSec, &numLeafMatrices));
2484: PetscCall(PetscMalloc2(numLeafIndices, &leafIndices, numLeafMatrices, &leafMatrices));
2485: PetscCall(PetscSFBcastBegin(indicesSF, MPIU_INT, rootIndices, leafIndices, MPI_REPLACE));
2486: PetscCall(PetscSFBcastBegin(matricesSF, MPIU_SCALAR, rootMatrices, leafMatrices, MPI_REPLACE));
2487: PetscCall(PetscSFBcastEnd(indicesSF, MPIU_INT, rootIndices, leafIndices, MPI_REPLACE));
2488: PetscCall(PetscSFBcastEnd(matricesSF, MPIU_SCALAR, rootMatrices, leafMatrices, MPI_REPLACE));
2489: PetscCall(PetscSFDestroy(&matricesSF));
2490: PetscCall(PetscSFDestroy(&indicesSF));
2491: PetscCall(PetscFree2(rootIndices, rootMatrices));
2492: PetscCall(PetscSectionDestroy(&rootIndicesSec));
2493: PetscCall(PetscSectionDestroy(&rootMatricesSec));
2494: }
2495: /* count to preallocate */
2496: PetscCall(DMGetLocalSection(fine, &localFine));
2497: {
2498: PetscInt nGlobal;
2499: PetscInt *dnnz, *onnz;
2500: PetscLayout rowMap, colMap;
2501: PetscInt rowStart, rowEnd, colStart, colEnd;
2502: PetscInt maxDof;
2503: PetscInt *rowIndices;
2504: DM refTree;
2505: PetscInt **refPointFieldN;
2506: PetscScalar ***refPointFieldMats;
2507: PetscSection refConSec, refAnSec;
2508: PetscInt pRefStart, pRefEnd, maxConDof, maxColumns, leafStart, leafEnd;
2509: PetscScalar *pointWork;
2511: PetscCall(PetscSectionGetConstrainedStorageSize(globalFine, &nGlobal));
2512: PetscCall(PetscCalloc2(nGlobal, &dnnz, nGlobal, &onnz));
2513: PetscCall(MatGetLayouts(mat, &rowMap, &colMap));
2514: PetscCall(PetscLayoutSetUp(rowMap));
2515: PetscCall(PetscLayoutSetUp(colMap));
2516: PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd));
2517: PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd));
2518: PetscCall(PetscSectionGetMaxDof(localFine, &maxDof));
2519: PetscCall(PetscSectionGetChart(leafIndicesSec, &leafStart, &leafEnd));
2520: PetscCall(DMGetWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
2521: for (p = leafStart; p < leafEnd; p++) {
2522: PetscInt gDof, gcDof, gOff;
2523: PetscInt numColIndices, pIndOff, *pInd;
2524: PetscInt matSize;
2525: PetscInt i;
2527: PetscCall(PetscSectionGetDof(globalFine, p, &gDof));
2528: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &gcDof));
2529: if ((gDof - gcDof) <= 0) continue;
2530: PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
2531: PetscCheck(gOff >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "I though having global dofs meant a non-negative offset");
2532: PetscCheck(gOff >= rowStart && (gOff + gDof - gcDof) <= rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "I thought the row map would constrain the global dofs");
2533: PetscCall(PetscSectionGetDof(leafIndicesSec, p, &numColIndices));
2534: PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &pIndOff));
2535: numColIndices -= 2 * numFields;
2536: PetscCheck(numColIndices > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "global fine dof with no dofs to interpolate from");
2537: pInd = &leafIndices[pIndOff];
2538: offsets[0] = 0;
2539: offsetsCopy[0] = 0;
2540: newOffsets[0] = 0;
2541: newOffsetsCopy[0] = 0;
2542: if (numFields) {
2543: PetscInt f;
2544: for (f = 0; f < numFields; f++) {
2545: PetscInt rowDof;
2547: PetscCall(PetscSectionGetFieldDof(localFine, p, f, &rowDof));
2548: offsets[f + 1] = offsets[f] + rowDof;
2549: offsetsCopy[f + 1] = offsets[f + 1];
2550: newOffsets[f + 1] = pInd[numColIndices + numFields + f];
2551: numD[f] = 0;
2552: numO[f] = 0;
2553: }
2554: PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices));
2555: for (f = 0; f < numFields; f++) {
2556: PetscInt colOffset = newOffsets[f];
2557: PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f];
2559: for (i = 0; i < numFieldCols; i++) {
2560: PetscInt gInd = pInd[i + colOffset];
2562: if (gInd >= colStart && gInd < colEnd) {
2563: numD[f]++;
2564: } else if (gInd >= 0) { /* negative means non-entry */
2565: numO[f]++;
2566: }
2567: }
2568: }
2569: } else {
2570: PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices));
2571: numD[0] = 0;
2572: numO[0] = 0;
2573: for (i = 0; i < numColIndices; i++) {
2574: PetscInt gInd = pInd[i];
2576: if (gInd >= colStart && gInd < colEnd) {
2577: numD[0]++;
2578: } else if (gInd >= 0) { /* negative means non-entry */
2579: numO[0]++;
2580: }
2581: }
2582: }
2583: PetscCall(PetscSectionGetDof(leafMatricesSec, p, &matSize));
2584: if (!matSize) { /* incoming matrix is identity */
2585: PetscInt childId;
2587: childId = childIds[p - pStartF];
2588: if (childId < 0) { /* no child interpolation: one nnz per */
2589: if (numFields) {
2590: PetscInt f;
2591: for (f = 0; f < numFields; f++) {
2592: PetscInt numRows = offsets[f + 1] - offsets[f], row;
2593: for (row = 0; row < numRows; row++) {
2594: PetscInt gIndCoarse = pInd[newOffsets[f] + row];
2595: PetscInt gIndFine = rowIndices[offsets[f] + row];
2596: if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2597: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2598: dnnz[gIndFine - rowStart] = 1;
2599: } else if (gIndCoarse >= 0) { /* remote */
2600: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2601: onnz[gIndFine - rowStart] = 1;
2602: } else { /* constrained */
2603: PetscCheck(gIndFine < 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2604: }
2605: }
2606: }
2607: } else {
2608: PetscInt i;
2609: for (i = 0; i < gDof; i++) {
2610: PetscInt gIndCoarse = pInd[i];
2611: PetscInt gIndFine = rowIndices[i];
2612: if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2613: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2614: dnnz[gIndFine - rowStart] = 1;
2615: } else if (gIndCoarse >= 0) { /* remote */
2616: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2617: onnz[gIndFine - rowStart] = 1;
2618: } else { /* constrained */
2619: PetscCheck(gIndFine < 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2620: }
2621: }
2622: }
2623: } else { /* interpolate from all */
2624: if (numFields) {
2625: PetscInt f;
2626: for (f = 0; f < numFields; f++) {
2627: PetscInt numRows = offsets[f + 1] - offsets[f], row;
2628: for (row = 0; row < numRows; row++) {
2629: PetscInt gIndFine = rowIndices[offsets[f] + row];
2630: if (gIndFine >= 0) {
2631: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2632: dnnz[gIndFine - rowStart] = numD[f];
2633: onnz[gIndFine - rowStart] = numO[f];
2634: }
2635: }
2636: }
2637: } else {
2638: PetscInt i;
2639: for (i = 0; i < gDof; i++) {
2640: PetscInt gIndFine = rowIndices[i];
2641: if (gIndFine >= 0) {
2642: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2643: dnnz[gIndFine - rowStart] = numD[0];
2644: onnz[gIndFine - rowStart] = numO[0];
2645: }
2646: }
2647: }
2648: }
2649: } else { /* interpolate from all */
2650: if (numFields) {
2651: PetscInt f;
2652: for (f = 0; f < numFields; f++) {
2653: PetscInt numRows = offsets[f + 1] - offsets[f], row;
2654: for (row = 0; row < numRows; row++) {
2655: PetscInt gIndFine = rowIndices[offsets[f] + row];
2656: if (gIndFine >= 0) {
2657: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2658: dnnz[gIndFine - rowStart] = numD[f];
2659: onnz[gIndFine - rowStart] = numO[f];
2660: }
2661: }
2662: }
2663: } else { /* every dof get a full row */
2664: PetscInt i;
2665: for (i = 0; i < gDof; i++) {
2666: PetscInt gIndFine = rowIndices[i];
2667: if (gIndFine >= 0) {
2668: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2669: dnnz[gIndFine - rowStart] = numD[0];
2670: onnz[gIndFine - rowStart] = numO[0];
2671: }
2672: }
2673: }
2674: }
2675: }
2676: PetscCall(MatXAIJSetPreallocation(mat, 1, dnnz, onnz, NULL, NULL));
2677: PetscCall(PetscFree2(dnnz, onnz));
2679: PetscCall(DMPlexGetReferenceTree(fine, &refTree));
2680: PetscCall(DMCopyDisc(fine, refTree));
2681: PetscCall(DMSetLocalSection(refTree, NULL));
2682: PetscCall(DMSetDefaultConstraints(refTree, NULL, NULL, NULL));
2683: PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
2684: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
2685: PetscCall(DMPlexGetAnchors(refTree, &refAnSec, NULL));
2686: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
2687: PetscCall(PetscSectionGetMaxDof(refConSec, &maxConDof));
2688: PetscCall(PetscSectionGetMaxDof(leafIndicesSec, &maxColumns));
2689: PetscCall(PetscMalloc1(maxConDof * maxColumns, &pointWork));
2690: for (p = leafStart; p < leafEnd; p++) {
2691: PetscInt gDof, gcDof, gOff;
2692: PetscInt numColIndices, pIndOff, *pInd;
2693: PetscInt matSize;
2694: PetscInt childId;
2696: PetscCall(PetscSectionGetDof(globalFine, p, &gDof));
2697: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &gcDof));
2698: if ((gDof - gcDof) <= 0) continue;
2699: childId = childIds[p - pStartF];
2700: PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
2701: PetscCall(PetscSectionGetDof(leafIndicesSec, p, &numColIndices));
2702: PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &pIndOff));
2703: numColIndices -= 2 * numFields;
2704: pInd = &leafIndices[pIndOff];
2705: offsets[0] = 0;
2706: offsetsCopy[0] = 0;
2707: newOffsets[0] = 0;
2708: newOffsetsCopy[0] = 0;
2709: rowOffsets[0] = 0;
2710: if (numFields) {
2711: PetscInt f;
2712: for (f = 0; f < numFields; f++) {
2713: PetscInt rowDof;
2715: PetscCall(PetscSectionGetFieldDof(localFine, p, f, &rowDof));
2716: offsets[f + 1] = offsets[f] + rowDof;
2717: offsetsCopy[f + 1] = offsets[f + 1];
2718: rowOffsets[f + 1] = pInd[numColIndices + f];
2719: newOffsets[f + 1] = pInd[numColIndices + numFields + f];
2720: }
2721: PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices));
2722: } else {
2723: PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices));
2724: }
2725: PetscCall(PetscSectionGetDof(leafMatricesSec, p, &matSize));
2726: if (!matSize) { /* incoming matrix is identity */
2727: if (childId < 0) { /* no child interpolation: scatter */
2728: if (numFields) {
2729: PetscInt f;
2730: for (f = 0; f < numFields; f++) {
2731: PetscInt numRows = offsets[f + 1] - offsets[f], row;
2732: for (row = 0; row < numRows; row++) PetscCall(MatSetValue(mat, rowIndices[offsets[f] + row], pInd[newOffsets[f] + row], 1., INSERT_VALUES));
2733: }
2734: } else {
2735: PetscInt numRows = gDof, row;
2736: for (row = 0; row < numRows; row++) PetscCall(MatSetValue(mat, rowIndices[row], pInd[row], 1., INSERT_VALUES));
2737: }
2738: } else { /* interpolate from all */
2739: if (numFields) {
2740: PetscInt f;
2741: for (f = 0; f < numFields; f++) {
2742: PetscInt numRows = offsets[f + 1] - offsets[f];
2743: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2744: PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], refPointFieldMats[childId - pRefStart][f], INSERT_VALUES));
2745: }
2746: } else {
2747: PetscCall(MatSetValues(mat, gDof, rowIndices, numColIndices, pInd, refPointFieldMats[childId - pRefStart][0], INSERT_VALUES));
2748: }
2749: }
2750: } else { /* interpolate from all */
2751: PetscInt pMatOff;
2752: PetscScalar *pMat;
2754: PetscCall(PetscSectionGetOffset(leafMatricesSec, p, &pMatOff));
2755: pMat = &leafMatrices[pMatOff];
2756: if (childId < 0) { /* copy the incoming matrix */
2757: if (numFields) {
2758: PetscInt f, count;
2759: for (f = 0, count = 0; f < numFields; f++) {
2760: PetscInt numRows = offsets[f + 1] - offsets[f];
2761: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2762: PetscInt numInRows = rowOffsets[f + 1] - rowOffsets[f];
2763: PetscScalar *inMat = &pMat[count];
2765: PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], inMat, INSERT_VALUES));
2766: count += numCols * numInRows;
2767: }
2768: } else {
2769: PetscCall(MatSetValues(mat, gDof, rowIndices, numColIndices, pInd, pMat, INSERT_VALUES));
2770: }
2771: } else { /* multiply the incoming matrix by the child interpolation */
2772: if (numFields) {
2773: PetscInt f, count;
2774: for (f = 0, count = 0; f < numFields; f++) {
2775: PetscInt numRows = offsets[f + 1] - offsets[f];
2776: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2777: PetscInt numInRows = rowOffsets[f + 1] - rowOffsets[f];
2778: PetscScalar *inMat = &pMat[count];
2779: PetscInt i, j, k;
2780: PetscCheck(refPointFieldN[childId - pRefStart][f] == numInRows, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point constraint matrix multiply dimension mismatch");
2781: for (i = 0; i < numRows; i++) {
2782: for (j = 0; j < numCols; j++) {
2783: PetscScalar val = 0.;
2784: for (k = 0; k < numInRows; k++) val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j];
2785: pointWork[i * numCols + j] = val;
2786: }
2787: }
2788: PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], pointWork, INSERT_VALUES));
2789: count += numCols * numInRows;
2790: }
2791: } else { /* every dof gets a full row */
2792: PetscInt numRows = gDof;
2793: PetscInt numCols = numColIndices;
2794: PetscInt numInRows = matSize / numColIndices;
2795: PetscInt i, j, k;
2796: PetscCheck(refPointFieldN[childId - pRefStart][0] == numInRows, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point constraint matrix multiply dimension mismatch");
2797: for (i = 0; i < numRows; i++) {
2798: for (j = 0; j < numCols; j++) {
2799: PetscScalar val = 0.;
2800: for (k = 0; k < numInRows; k++) val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j];
2801: pointWork[i * numCols + j] = val;
2802: }
2803: }
2804: PetscCall(MatSetValues(mat, numRows, rowIndices, numCols, pInd, pointWork, INSERT_VALUES));
2805: }
2806: }
2807: }
2808: }
2809: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
2810: PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
2811: PetscCall(PetscFree(pointWork));
2812: }
2813: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
2814: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
2815: PetscCall(PetscSectionDestroy(&leafIndicesSec));
2816: PetscCall(PetscSectionDestroy(&leafMatricesSec));
2817: PetscCall(PetscFree2(leafIndices, leafMatrices));
2818: PetscCall(PetscFree2(*(PetscInt ****)&perms, *(PetscScalar ****)&flips));
2819: PetscCall(PetscFree7(offsets, offsetsCopy, newOffsets, newOffsetsCopy, rowOffsets, numD, numO));
2820: PetscCall(ISRestoreIndices(aIS, &anchors));
2821: PetscFunctionReturn(PETSC_SUCCESS);
2822: }
2824: /*
2825: * Assuming a nodal basis (w.r.t. the dual basis) basis:
2826: *
2827: * for each coarse dof \phi^c_i:
2828: * for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i:
2829: * for each fine dof \phi^f_j;
2830: * a_{i,j} = 0;
2831: * for each fine dof \phi^f_k:
2832: * a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l
2833: * [^^^ this is = \phi^c_i ^^^]
2834: */
2835: /*@
2836: DMPlexComputeInjectorReferenceTree - Compute the injection matrix from fine to coarse degrees of freedom on the reference tree
2838: Collective
2840: Input Parameter:
2841: . refTree - The reference-tree `DMPLEX` (see `DMPlexCreateDefaultReferenceTree()`)
2843: Output Parameter:
2844: . inj - The newly created injection `Mat` mapping fine-space coefficients on the reference tree to their coarse-space counterparts
2846: Level: developer
2848: Note:
2849: For a nodal basis, the injection is derived from the constraint matrix attached to the reference
2850: tree; the returned matrix is used internally by `DMPlexComputeInjectorTree()` to construct the
2851: global injection between refined and coarse meshes.
2853: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSetReferenceTree()`, `DMPlexCreateDefaultReferenceTree()`, `DMPlexComputeInjectorTree()`, `DMPlexComputeInterpolatorTree()`
2854: @*/
2855: PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj)
2856: {
2857: PetscDS ds;
2858: PetscSection section, cSection;
2859: DMLabel canonical, depth;
2860: Mat cMat, mat;
2861: PetscInt *nnz;
2862: PetscInt f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd;
2863: PetscInt m, n;
2864: PetscScalar *pointScalar;
2865: PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent;
2867: PetscFunctionBegin;
2868: PetscCall(DMGetLocalSection(refTree, §ion));
2869: PetscCall(DMGetDimension(refTree, &dim));
2870: PetscCall(PetscMalloc6(dim, &v0, dim, &v0parent, dim, &vtmp, dim * dim, &J, dim * dim, &Jparent, dim * dim, &invJ));
2871: PetscCall(PetscMalloc2(dim, &pointScalar, dim, &pointRef));
2872: PetscCall(DMGetDS(refTree, &ds));
2873: PetscCall(PetscDSGetNumFields(ds, &numFields));
2874: PetscCall(PetscSectionGetNumFields(section, &numSecFields));
2875: PetscCall(DMGetLabel(refTree, "canonical", &canonical));
2876: PetscCall(DMGetLabel(refTree, "depth", &depth));
2877: PetscCall(DMGetDefaultConstraints(refTree, &cSection, &cMat, NULL));
2878: PetscCall(DMPlexGetChart(refTree, &pStart, &pEnd));
2879: PetscCall(DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd));
2880: PetscCall(MatGetSize(cMat, &n, &m)); /* the injector has transpose sizes from the constraint matrix */
2881: /* Step 1: compute non-zero pattern. A proper subset of constraint matrix non-zero */
2882: PetscCall(PetscCalloc1(m, &nnz));
2883: for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */
2884: const PetscInt *children;
2885: PetscInt numChildren;
2886: PetscInt i, numChildDof, numSelfDof;
2888: if (canonical) {
2889: PetscInt pCanonical;
2890: PetscCall(DMLabelGetValue(canonical, p, &pCanonical));
2891: if (p != pCanonical) continue;
2892: }
2893: PetscCall(DMPlexGetTreeChildren(refTree, p, &numChildren, &children));
2894: if (!numChildren) continue;
2895: for (i = 0, numChildDof = 0; i < numChildren; i++) {
2896: PetscInt child = children[i];
2897: PetscInt dof;
2899: PetscCall(PetscSectionGetDof(section, child, &dof));
2900: numChildDof += dof;
2901: }
2902: PetscCall(PetscSectionGetDof(section, p, &numSelfDof));
2903: if (!numChildDof || !numSelfDof) continue;
2904: for (f = 0; f < numFields; f++) {
2905: PetscInt selfOff;
2907: if (numSecFields) { /* count the dofs for just this field */
2908: for (i = 0, numChildDof = 0; i < numChildren; i++) {
2909: PetscInt child = children[i];
2910: PetscInt dof;
2912: PetscCall(PetscSectionGetFieldDof(section, child, f, &dof));
2913: numChildDof += dof;
2914: }
2915: PetscCall(PetscSectionGetFieldDof(section, p, f, &numSelfDof));
2916: PetscCall(PetscSectionGetFieldOffset(section, p, f, &selfOff));
2917: } else {
2918: PetscCall(PetscSectionGetOffset(section, p, &selfOff));
2919: }
2920: for (i = 0; i < numSelfDof; i++) nnz[selfOff + i] = numChildDof;
2921: }
2922: }
2923: PetscCall(MatCreateAIJ(PETSC_COMM_SELF, m, n, m, n, -1, nnz, -1, NULL, &mat));
2924: PetscCall(PetscFree(nnz));
2925: /* Setp 2: compute entries */
2926: for (p = pStart; p < pEnd; p++) {
2927: const PetscInt *children;
2928: PetscInt numChildren;
2929: PetscInt i, numChildDof, numSelfDof;
2931: /* same conditions about when entries occur */
2932: if (canonical) {
2933: PetscInt pCanonical;
2934: PetscCall(DMLabelGetValue(canonical, p, &pCanonical));
2935: if (p != pCanonical) continue;
2936: }
2937: PetscCall(DMPlexGetTreeChildren(refTree, p, &numChildren, &children));
2938: if (!numChildren) continue;
2939: for (i = 0, numChildDof = 0; i < numChildren; i++) {
2940: PetscInt child = children[i];
2941: PetscInt dof;
2943: PetscCall(PetscSectionGetDof(section, child, &dof));
2944: numChildDof += dof;
2945: }
2946: PetscCall(PetscSectionGetDof(section, p, &numSelfDof));
2947: if (!numChildDof || !numSelfDof) continue;
2949: for (f = 0; f < numFields; f++) {
2950: PetscInt pI = -1, cI = -1;
2951: PetscInt selfOff, Nc, parentCell;
2952: PetscInt cellShapeOff;
2953: PetscObject disc;
2954: PetscDualSpace dsp;
2955: PetscClassId classId;
2956: PetscScalar *pointMat;
2957: PetscInt *matRows, *matCols;
2958: PetscInt pO = PETSC_INT_MIN;
2959: const PetscInt *depthNumDof;
2961: if (numSecFields) {
2962: for (i = 0, numChildDof = 0; i < numChildren; i++) {
2963: PetscInt child = children[i];
2964: PetscInt dof;
2966: PetscCall(PetscSectionGetFieldDof(section, child, f, &dof));
2967: numChildDof += dof;
2968: }
2969: PetscCall(PetscSectionGetFieldDof(section, p, f, &numSelfDof));
2970: PetscCall(PetscSectionGetFieldOffset(section, p, f, &selfOff));
2971: } else {
2972: PetscCall(PetscSectionGetOffset(section, p, &selfOff));
2973: }
2975: /* find a cell whose closure contains p */
2976: if (p >= cStart && p < cEnd) {
2977: parentCell = p;
2978: } else {
2979: PetscInt *star = NULL;
2980: PetscInt numStar;
2982: parentCell = -1;
2983: PetscCall(DMPlexGetTransitiveClosure(refTree, p, PETSC_FALSE, &numStar, &star));
2984: for (i = numStar - 1; i >= 0; i--) {
2985: PetscInt c = star[2 * i];
2987: if (c >= cStart && c < cEnd) {
2988: parentCell = c;
2989: break;
2990: }
2991: }
2992: PetscCall(DMPlexRestoreTransitiveClosure(refTree, p, PETSC_FALSE, &numStar, &star));
2993: }
2994: /* determine the offset of p's shape functions within parentCell's shape functions */
2995: PetscCall(PetscDSGetDiscretization(ds, f, &disc));
2996: PetscCall(PetscObjectGetClassId(disc, &classId));
2997: if (classId == PETSCFE_CLASSID) PetscCall(PetscFEGetDualSpace((PetscFE)disc, &dsp));
2998: else {
2999: PetscCheck(classId == PETSCFV_CLASSID, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported discretization object");
3000: PetscCall(PetscFVGetDualSpace((PetscFV)disc, &dsp));
3001: }
3002: PetscCall(PetscDualSpaceGetNumDof(dsp, &depthNumDof));
3003: PetscCall(PetscDualSpaceGetNumComponents(dsp, &Nc));
3004: {
3005: PetscInt *closure = NULL;
3006: PetscInt numClosure;
3008: PetscCall(DMPlexGetTransitiveClosure(refTree, parentCell, PETSC_TRUE, &numClosure, &closure));
3009: for (i = 0, pI = -1, cellShapeOff = 0; i < numClosure; i++) {
3010: PetscInt point = closure[2 * i], pointDepth;
3012: pO = closure[2 * i + 1];
3013: if (point == p) {
3014: pI = i;
3015: break;
3016: }
3017: PetscCall(DMLabelGetValue(depth, point, &pointDepth));
3018: cellShapeOff += depthNumDof[pointDepth];
3019: }
3020: PetscCall(DMPlexRestoreTransitiveClosure(refTree, parentCell, PETSC_TRUE, &numClosure, &closure));
3021: }
3023: PetscCall(DMGetWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR, &pointMat));
3024: PetscCall(DMGetWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT, &matRows));
3025: matCols = matRows + numSelfDof;
3026: for (i = 0; i < numSelfDof; i++) matRows[i] = selfOff + i;
3027: for (i = 0; i < numSelfDof * numChildDof; i++) pointMat[i] = 0.;
3028: {
3029: PetscInt colOff = 0;
3031: for (i = 0; i < numChildren; i++) {
3032: PetscInt child = children[i];
3033: PetscInt dof, off, j;
3035: if (numSecFields) {
3036: PetscCall(PetscSectionGetFieldDof(cSection, child, f, &dof));
3037: PetscCall(PetscSectionGetFieldOffset(cSection, child, f, &off));
3038: } else {
3039: PetscCall(PetscSectionGetDof(cSection, child, &dof));
3040: PetscCall(PetscSectionGetOffset(cSection, child, &off));
3041: }
3043: for (j = 0; j < dof; j++) matCols[colOff++] = off + j;
3044: }
3045: }
3046: if (classId == PETSCFE_CLASSID) {
3047: PetscFE fe = (PetscFE)disc;
3048: PetscInt fSize;
3049: const PetscInt ***perms;
3050: const PetscScalar ***flips;
3051: const PetscInt *pperms;
3053: PetscCall(PetscFEGetDualSpace(fe, &dsp));
3054: PetscCall(PetscDualSpaceGetDimension(dsp, &fSize));
3055: PetscCall(PetscDualSpaceGetSymmetries(dsp, &perms, &flips));
3056: pperms = perms ? perms[pI] ? perms[pI][pO] : NULL : NULL;
3057: for (i = 0; i < numSelfDof; i++) { /* for every shape function */
3058: PetscQuadrature q;
3059: PetscInt dim, thisNc, numPoints, j, k;
3060: const PetscReal *points;
3061: const PetscReal *weights;
3062: PetscInt *closure = NULL;
3063: PetscInt numClosure;
3064: PetscInt iCell = pperms ? pperms[i] : i;
3065: PetscInt parentCellShapeDof = cellShapeOff + iCell;
3066: PetscTabulation Tparent;
3068: PetscCall(PetscDualSpaceGetFunctional(dsp, parentCellShapeDof, &q));
3069: PetscCall(PetscQuadratureGetData(q, &dim, &thisNc, &numPoints, &points, &weights));
3070: PetscCheck(thisNc == Nc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Functional dim %" PetscInt_FMT " does not much basis dim %" PetscInt_FMT, thisNc, Nc);
3071: PetscCall(PetscFECreateTabulation(fe, 1, numPoints, points, 0, &Tparent)); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */
3072: for (j = 0; j < numPoints; j++) {
3073: PetscInt childCell = -1;
3074: PetscReal *parentValAtPoint;
3075: const PetscReal xi0[3] = {-1., -1., -1.};
3076: const PetscReal *pointReal = &points[dim * j];
3077: const PetscScalar *point;
3078: PetscTabulation Tchild;
3079: PetscInt childCellShapeOff, pointMatOff;
3080: #if defined(PETSC_USE_COMPLEX)
3081: PetscInt d;
3083: for (d = 0; d < dim; d++) pointScalar[d] = points[dim * j + d];
3084: point = pointScalar;
3085: #else
3086: point = pointReal;
3087: #endif
3089: parentValAtPoint = &Tparent->T[0][(fSize * j + parentCellShapeDof) * Nc];
3091: for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/
3092: PetscInt child = children[k];
3093: PetscInt *star = NULL;
3094: PetscInt numStar, s;
3096: PetscCall(DMPlexGetTransitiveClosure(refTree, child, PETSC_FALSE, &numStar, &star));
3097: for (s = numStar - 1; s >= 0; s--) {
3098: PetscInt c = star[2 * s];
3100: if (c < cStart || c >= cEnd) continue;
3101: PetscCall(DMPlexLocatePoint_Internal(refTree, dim, point, c, &childCell));
3102: if (childCell >= 0) break;
3103: }
3104: PetscCall(DMPlexRestoreTransitiveClosure(refTree, child, PETSC_FALSE, &numStar, &star));
3105: if (childCell >= 0) break;
3106: }
3107: PetscCheck(childCell >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not locate quadrature point");
3108: PetscCall(DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ));
3109: PetscCall(DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent));
3110: CoordinatesRefToReal(dim, dim, xi0, v0parent, Jparent, pointReal, vtmp);
3111: CoordinatesRealToRef(dim, dim, xi0, v0, invJ, vtmp, pointRef);
3113: PetscCall(PetscFECreateTabulation(fe, 1, 1, pointRef, 0, &Tchild));
3114: PetscCall(DMPlexGetTransitiveClosure(refTree, childCell, PETSC_TRUE, &numClosure, &closure));
3115: for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */
3116: PetscInt child = children[k], childDepth, childDof, childO = PETSC_INT_MIN;
3117: PetscInt l;
3118: const PetscInt *cperms;
3120: PetscCall(DMLabelGetValue(depth, child, &childDepth));
3121: childDof = depthNumDof[childDepth];
3122: for (l = 0, cI = -1, childCellShapeOff = 0; l < numClosure; l++) {
3123: PetscInt point = closure[2 * l];
3124: PetscInt pointDepth;
3126: childO = closure[2 * l + 1];
3127: if (point == child) {
3128: cI = l;
3129: break;
3130: }
3131: PetscCall(DMLabelGetValue(depth, point, &pointDepth));
3132: childCellShapeOff += depthNumDof[pointDepth];
3133: }
3134: if (l == numClosure) {
3135: pointMatOff += childDof;
3136: continue; /* child is not in the closure of the cell: has nothing to contribute to this point */
3137: }
3138: cperms = perms ? perms[cI] ? perms[cI][childO] : NULL : NULL;
3139: for (l = 0; l < childDof; l++) {
3140: PetscInt lCell = cperms ? cperms[l] : l;
3141: PetscInt childCellDof = childCellShapeOff + lCell;
3142: PetscReal *childValAtPoint;
3143: PetscReal val = 0.;
3145: childValAtPoint = &Tchild->T[0][childCellDof * Nc];
3146: for (m = 0; m < Nc; m++) val += weights[j * Nc + m] * parentValAtPoint[m] * childValAtPoint[m];
3148: pointMat[i * numChildDof + pointMatOff + l] += val;
3149: }
3150: pointMatOff += childDof;
3151: }
3152: PetscCall(DMPlexRestoreTransitiveClosure(refTree, childCell, PETSC_TRUE, &numClosure, &closure));
3153: PetscCall(PetscTabulationDestroy(&Tchild));
3154: }
3155: PetscCall(PetscTabulationDestroy(&Tparent));
3156: }
3157: } else { /* just the volume-weighted averages of the children */
3158: PetscReal parentVol;
3159: PetscInt childCell;
3161: PetscCall(DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL));
3162: for (i = 0, childCell = 0; i < numChildren; i++) {
3163: PetscInt child = children[i], j;
3164: PetscReal childVol;
3166: if (child < cStart || child >= cEnd) continue;
3167: PetscCall(DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL));
3168: for (j = 0; j < Nc; j++) pointMat[j * numChildDof + Nc * childCell + j] = childVol / parentVol;
3169: childCell++;
3170: }
3171: }
3172: /* Insert pointMat into mat */
3173: PetscCall(MatSetValues(mat, numSelfDof, matRows, numChildDof, matCols, pointMat, INSERT_VALUES));
3174: PetscCall(DMRestoreWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT, &matRows));
3175: PetscCall(DMRestoreWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR, &pointMat));
3176: }
3177: }
3178: PetscCall(PetscFree6(v0, v0parent, vtmp, J, Jparent, invJ));
3179: PetscCall(PetscFree2(pointScalar, pointRef));
3180: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
3181: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
3182: *inj = mat;
3183: PetscFunctionReturn(PETSC_SUCCESS);
3184: }
3186: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3187: {
3188: PetscDS ds;
3189: PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof;
3190: PetscScalar ***refPointFieldMats;
3191: PetscSection refConSec, refSection;
3193: PetscFunctionBegin;
3194: PetscCall(DMGetDS(refTree, &ds));
3195: PetscCall(PetscDSGetNumFields(ds, &numFields));
3196: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
3197: PetscCall(DMGetLocalSection(refTree, &refSection));
3198: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
3199: PetscCall(PetscMalloc1(pRefEnd - pRefStart, &refPointFieldMats));
3200: PetscCall(PetscSectionGetMaxDof(refConSec, &maxDof));
3201: PetscCall(PetscMalloc1(maxDof, &rows));
3202: PetscCall(PetscMalloc1(maxDof * maxDof, &cols));
3203: for (p = pRefStart; p < pRefEnd; p++) {
3204: PetscInt parent, pDof, parentDof;
3206: PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL));
3207: PetscCall(PetscSectionGetDof(refConSec, p, &pDof));
3208: PetscCall(PetscSectionGetDof(refSection, parent, &parentDof));
3209: if (!pDof || !parentDof || parent == p) continue;
3211: PetscCall(PetscMalloc1(numFields, &refPointFieldMats[p - pRefStart]));
3212: for (f = 0; f < numFields; f++) {
3213: PetscInt cDof, cOff, numCols, r;
3215: if (numFields > 1) {
3216: PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof));
3217: PetscCall(PetscSectionGetFieldOffset(refConSec, p, f, &cOff));
3218: } else {
3219: PetscCall(PetscSectionGetDof(refConSec, p, &cDof));
3220: PetscCall(PetscSectionGetOffset(refConSec, p, &cOff));
3221: }
3223: for (r = 0; r < cDof; r++) rows[r] = cOff + r;
3224: numCols = 0;
3225: {
3226: PetscInt aDof, aOff, j;
3228: if (numFields > 1) {
3229: PetscCall(PetscSectionGetFieldDof(refSection, parent, f, &aDof));
3230: PetscCall(PetscSectionGetFieldOffset(refSection, parent, f, &aOff));
3231: } else {
3232: PetscCall(PetscSectionGetDof(refSection, parent, &aDof));
3233: PetscCall(PetscSectionGetOffset(refSection, parent, &aOff));
3234: }
3236: for (j = 0; j < aDof; j++) cols[numCols++] = aOff + j;
3237: }
3238: PetscCall(PetscMalloc1(cDof * numCols, &refPointFieldMats[p - pRefStart][f]));
3239: /* transpose of constraint matrix */
3240: PetscCall(MatGetValues(inj, numCols, cols, cDof, rows, refPointFieldMats[p - pRefStart][f]));
3241: }
3242: }
3243: *childrenMats = refPointFieldMats;
3244: PetscCall(PetscFree(rows));
3245: PetscCall(PetscFree(cols));
3246: PetscFunctionReturn(PETSC_SUCCESS);
3247: }
3249: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3250: {
3251: PetscDS ds;
3252: PetscScalar ***refPointFieldMats;
3253: PetscInt numFields, pRefStart, pRefEnd, p, f;
3254: PetscSection refConSec, refSection;
3256: PetscFunctionBegin;
3257: refPointFieldMats = *childrenMats;
3258: *childrenMats = NULL;
3259: PetscCall(DMGetDS(refTree, &ds));
3260: PetscCall(DMGetLocalSection(refTree, &refSection));
3261: PetscCall(PetscDSGetNumFields(ds, &numFields));
3262: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
3263: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
3264: for (p = pRefStart; p < pRefEnd; p++) {
3265: PetscInt parent, pDof, parentDof;
3267: PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL));
3268: PetscCall(PetscSectionGetDof(refConSec, p, &pDof));
3269: PetscCall(PetscSectionGetDof(refSection, parent, &parentDof));
3270: if (!pDof || !parentDof || parent == p) continue;
3272: for (f = 0; f < numFields; f++) {
3273: PetscInt cDof;
3275: if (numFields > 1) {
3276: PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof));
3277: } else {
3278: PetscCall(PetscSectionGetDof(refConSec, p, &cDof));
3279: }
3281: PetscCall(PetscFree(refPointFieldMats[p - pRefStart][f]));
3282: }
3283: PetscCall(PetscFree(refPointFieldMats[p - pRefStart]));
3284: }
3285: PetscCall(PetscFree(refPointFieldMats));
3286: PetscFunctionReturn(PETSC_SUCCESS);
3287: }
3289: static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree, Mat *injRef)
3290: {
3291: Mat cMatRef;
3292: PetscObject injRefObj;
3294: PetscFunctionBegin;
3295: PetscCall(DMGetDefaultConstraints(refTree, NULL, &cMatRef, NULL));
3296: PetscCall(PetscObjectQuery((PetscObject)cMatRef, "DMPlexComputeInjectorTree_refTree", &injRefObj));
3297: *injRef = (Mat)injRefObj;
3298: if (!*injRef) {
3299: PetscCall(DMPlexComputeInjectorReferenceTree(refTree, injRef));
3300: PetscCall(PetscObjectCompose((PetscObject)cMatRef, "DMPlexComputeInjectorTree_refTree", (PetscObject)*injRef));
3301: /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */
3302: PetscCall(PetscObjectDereference((PetscObject)*injRef));
3303: }
3304: PetscFunctionReturn(PETSC_SUCCESS);
3305: }
3307: 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)
3308: {
3309: PetscInt pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti;
3310: PetscSection globalCoarse, globalFine;
3311: PetscSection localCoarse, localFine, leafIndicesSec;
3312: PetscSection multiRootSec, rootIndicesSec;
3313: PetscInt *leafInds, *rootInds = NULL;
3314: const PetscInt *rootDegrees;
3315: PetscScalar *leafVals = NULL, *rootVals = NULL;
3316: PetscSF coarseToFineEmbedded;
3318: PetscFunctionBegin;
3319: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
3320: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
3321: PetscCall(DMGetLocalSection(fine, &localFine));
3322: PetscCall(DMGetGlobalSection(fine, &globalFine));
3323: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafIndicesSec));
3324: PetscCall(PetscSectionSetChart(leafIndicesSec, pStartF, pEndF));
3325: PetscCall(PetscSectionGetMaxDof(localFine, &maxDof));
3326: { /* winnow fine points that don't have global dofs out of the sf */
3327: PetscInt l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices;
3328: const PetscInt *leaves;
3330: PetscCall(PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL));
3331: for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
3332: p = leaves ? leaves[l] : l;
3333: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
3334: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
3335: if ((dof - cdof) > 0) {
3336: numPointsWithDofs++;
3338: PetscCall(PetscSectionGetDof(localFine, p, &dof));
3339: PetscCall(PetscSectionSetDof(leafIndicesSec, p, dof + 1));
3340: }
3341: }
3342: PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
3343: PetscCall(PetscSectionSetUp(leafIndicesSec));
3344: PetscCall(PetscSectionGetStorageSize(leafIndicesSec, &numIndices));
3345: PetscCall(PetscMalloc1((gatheredIndices ? numIndices : (maxDof + 1)), &leafInds));
3346: if (gatheredValues) PetscCall(PetscMalloc1(numIndices, &leafVals));
3347: for (l = 0, offset = 0; l < nleaves; l++) {
3348: p = leaves ? leaves[l] : l;
3349: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
3350: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
3351: if ((dof - cdof) > 0) {
3352: PetscInt off, gOff;
3353: PetscInt *pInd;
3354: PetscScalar *pVal = NULL;
3356: pointsWithDofs[offset++] = l;
3358: PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &off));
3360: pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds;
3361: if (gatheredValues) {
3362: PetscInt i;
3364: pVal = &leafVals[off + 1];
3365: for (i = 0; i < dof; i++) pVal[i] = 0.;
3366: }
3367: PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
3369: offsets[0] = 0;
3370: if (numFields) {
3371: PetscInt f;
3373: for (f = 0; f < numFields; f++) {
3374: PetscInt fDof;
3375: PetscCall(PetscSectionGetFieldDof(localFine, p, f, &fDof));
3376: offsets[f + 1] = fDof + offsets[f];
3377: }
3378: PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, -1, NULL, pInd));
3379: } else {
3380: PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, NULL, pInd));
3381: }
3382: if (gatheredValues) PetscCall(VecGetValues(fineVec, dof, pInd, pVal));
3383: }
3384: }
3385: PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded));
3386: PetscCall(PetscFree(pointsWithDofs));
3387: }
3389: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
3390: PetscCall(DMGetLocalSection(coarse, &localCoarse));
3391: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
3393: { /* there may be the case where an sf root has a parent: broadcast parents back to children */
3394: MPI_Datatype threeInt;
3395: PetscMPIInt rank;
3396: PetscInt (*parentNodeAndIdCoarse)[3];
3397: PetscInt (*parentNodeAndIdFine)[3];
3398: PetscInt p, nleaves, nleavesToParents;
3399: PetscSF pointSF, sfToParents;
3400: const PetscInt *ilocal;
3401: const PetscSFNode *iremote;
3402: PetscSFNode *iremoteToParents;
3403: PetscInt *ilocalToParents;
3405: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank));
3406: PetscCallMPI(MPI_Type_contiguous(3, MPIU_INT, &threeInt));
3407: PetscCallMPI(MPI_Type_commit(&threeInt));
3408: PetscCall(PetscMalloc2(pEndC - pStartC, &parentNodeAndIdCoarse, pEndF - pStartF, &parentNodeAndIdFine));
3409: PetscCall(DMGetPointSF(coarse, &pointSF));
3410: PetscCall(PetscSFGetGraph(pointSF, NULL, &nleaves, &ilocal, &iremote));
3411: for (p = pStartC; p < pEndC; p++) {
3412: PetscInt parent, childId;
3413: PetscCall(DMPlexGetTreeParent(coarse, p, &parent, &childId));
3414: parentNodeAndIdCoarse[p - pStartC][0] = rank;
3415: parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC;
3416: parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId;
3417: if (nleaves > 0) {
3418: PetscInt leaf = -1;
3420: if (ilocal) {
3421: PetscCall(PetscFindInt(parent, nleaves, ilocal, &leaf));
3422: } else {
3423: leaf = p - pStartC;
3424: }
3425: if (leaf >= 0) {
3426: parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank;
3427: parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index;
3428: }
3429: }
3430: }
3431: for (p = pStartF; p < pEndF; p++) {
3432: parentNodeAndIdFine[p - pStartF][0] = -1;
3433: parentNodeAndIdFine[p - pStartF][1] = -1;
3434: parentNodeAndIdFine[p - pStartF][2] = -1;
3435: }
3436: PetscCall(PetscSFBcastBegin(coarseToFineEmbedded, threeInt, parentNodeAndIdCoarse, parentNodeAndIdFine, MPI_REPLACE));
3437: PetscCall(PetscSFBcastEnd(coarseToFineEmbedded, threeInt, parentNodeAndIdCoarse, parentNodeAndIdFine, MPI_REPLACE));
3438: for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3439: PetscInt dof;
3441: PetscCall(PetscSectionGetDof(leafIndicesSec, p, &dof));
3442: if (dof) {
3443: PetscInt off;
3445: PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &off));
3446: if (gatheredIndices) {
3447: leafInds[off] = PetscMax(childIds[p - pStartF], parentNodeAndIdFine[p - pStartF][2]);
3448: } else if (gatheredValues) {
3449: leafVals[off] = (PetscScalar)PetscMax(childIds[p - pStartF], parentNodeAndIdFine[p - pStartF][2]);
3450: }
3451: }
3452: if (parentNodeAndIdFine[p - pStartF][0] >= 0) nleavesToParents++;
3453: }
3454: PetscCall(PetscMalloc1(nleavesToParents, &ilocalToParents));
3455: PetscCall(PetscMalloc1(nleavesToParents, &iremoteToParents));
3456: for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3457: if (parentNodeAndIdFine[p - pStartF][0] >= 0) {
3458: ilocalToParents[nleavesToParents] = p - pStartF;
3459: // FIXME PetscCall(PetscMPIIntCast(parentNodeAndIdFine[p - pStartF][0],&iremoteToParents[nleavesToParents].rank));
3460: iremoteToParents[nleavesToParents].rank = parentNodeAndIdFine[p - pStartF][0];
3461: iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p - pStartF][1];
3462: nleavesToParents++;
3463: }
3464: }
3465: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)coarse), &sfToParents));
3466: PetscCall(PetscSFSetGraph(sfToParents, pEndC - pStartC, nleavesToParents, ilocalToParents, PETSC_OWN_POINTER, iremoteToParents, PETSC_OWN_POINTER));
3467: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
3469: coarseToFineEmbedded = sfToParents;
3471: PetscCall(PetscFree2(parentNodeAndIdCoarse, parentNodeAndIdFine));
3472: PetscCallMPI(MPI_Type_free(&threeInt));
3473: }
3475: { /* winnow out coarse points that don't have dofs */
3476: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;
3477: PetscSF sfDofsOnly;
3479: for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) {
3480: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3481: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
3482: if ((dof - cdof) > 0) numPointsWithDofs++;
3483: }
3484: PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
3485: for (p = pStartC, offset = 0; p < pEndC; p++) {
3486: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3487: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
3488: if ((dof - cdof) > 0) pointsWithDofs[offset++] = p - pStartC;
3489: }
3490: PetscCall(PetscSFCreateEmbeddedRootSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly));
3491: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
3492: PetscCall(PetscFree(pointsWithDofs));
3493: coarseToFineEmbedded = sfDofsOnly;
3494: }
3496: /* communicate back to the coarse mesh which coarse points have children (that may require injection) */
3497: PetscCall(PetscSFComputeDegreeBegin(coarseToFineEmbedded, &rootDegrees));
3498: PetscCall(PetscSFComputeDegreeEnd(coarseToFineEmbedded, &rootDegrees));
3499: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &multiRootSec));
3500: PetscCall(PetscSectionSetChart(multiRootSec, pStartC, pEndC));
3501: for (p = pStartC; p < pEndC; p++) PetscCall(PetscSectionSetDof(multiRootSec, p, rootDegrees[p - pStartC]));
3502: PetscCall(PetscSectionSetUp(multiRootSec));
3503: PetscCall(PetscSectionGetStorageSize(multiRootSec, &numMulti));
3504: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootIndicesSec));
3505: { /* distribute the leaf section */
3506: PetscSF multi, multiInv, indicesSF;
3507: PetscInt *remoteOffsets, numRootIndices;
3509: PetscCall(PetscSFGetMultiSF(coarseToFineEmbedded, &multi));
3510: PetscCall(PetscSFCreateInverseSF(multi, &multiInv));
3511: PetscCall(PetscSFDistributeSection(multiInv, leafIndicesSec, &remoteOffsets, rootIndicesSec));
3512: PetscCall(PetscSFCreateSectionSF(multiInv, leafIndicesSec, remoteOffsets, rootIndicesSec, &indicesSF));
3513: PetscCall(PetscFree(remoteOffsets));
3514: PetscCall(PetscSFDestroy(&multiInv));
3515: PetscCall(PetscSectionGetStorageSize(rootIndicesSec, &numRootIndices));
3516: if (gatheredIndices) {
3517: PetscCall(PetscMalloc1(numRootIndices, &rootInds));
3518: PetscCall(PetscSFBcastBegin(indicesSF, MPIU_INT, leafInds, rootInds, MPI_REPLACE));
3519: PetscCall(PetscSFBcastEnd(indicesSF, MPIU_INT, leafInds, rootInds, MPI_REPLACE));
3520: }
3521: if (gatheredValues) {
3522: PetscCall(PetscMalloc1(numRootIndices, &rootVals));
3523: PetscCall(PetscSFBcastBegin(indicesSF, MPIU_SCALAR, leafVals, rootVals, MPI_REPLACE));
3524: PetscCall(PetscSFBcastEnd(indicesSF, MPIU_SCALAR, leafVals, rootVals, MPI_REPLACE));
3525: }
3526: PetscCall(PetscSFDestroy(&indicesSF));
3527: }
3528: PetscCall(PetscSectionDestroy(&leafIndicesSec));
3529: PetscCall(PetscFree(leafInds));
3530: PetscCall(PetscFree(leafVals));
3531: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
3532: *rootMultiSec = multiRootSec;
3533: *multiLeafSec = rootIndicesSec;
3534: if (gatheredIndices) *gatheredIndices = rootInds;
3535: if (gatheredValues) *gatheredValues = rootVals;
3536: PetscFunctionReturn(PETSC_SUCCESS);
3537: }
3539: PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
3540: {
3541: DM refTree;
3542: PetscSection multiRootSec, rootIndicesSec;
3543: PetscSection globalCoarse, globalFine;
3544: PetscSection localCoarse, localFine;
3545: PetscSection cSecRef;
3546: PetscInt *rootIndices = NULL, *parentIndices, pRefStart, pRefEnd;
3547: Mat injRef;
3548: PetscInt numFields, maxDof;
3549: PetscInt pStartC, pEndC, pStartF, pEndF, p;
3550: PetscInt *offsets, *offsetsCopy, *rowOffsets;
3551: PetscLayout rowMap, colMap;
3552: PetscInt rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO;
3553: PetscScalar ***childrenMats = NULL; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */
3555: PetscFunctionBegin;
3556: /* get the templates for the fine-to-coarse injection from the reference tree */
3557: PetscCall(DMPlexGetReferenceTree(coarse, &refTree));
3558: PetscCall(DMCopyDisc(coarse, refTree));
3559: PetscCall(DMSetLocalSection(refTree, NULL));
3560: PetscCall(DMSetDefaultConstraints(refTree, NULL, NULL, NULL));
3561: PetscCall(DMGetDefaultConstraints(refTree, &cSecRef, NULL, NULL));
3562: PetscCall(PetscSectionGetChart(cSecRef, &pRefStart, &pRefEnd));
3563: PetscCall(DMPlexReferenceTreeGetInjector(refTree, &injRef));
3565: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
3566: PetscCall(DMGetLocalSection(fine, &localFine));
3567: PetscCall(DMGetGlobalSection(fine, &globalFine));
3568: PetscCall(PetscSectionGetNumFields(localFine, &numFields));
3569: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
3570: PetscCall(DMGetLocalSection(coarse, &localCoarse));
3571: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
3572: PetscCall(PetscSectionGetMaxDof(localCoarse, &maxDof));
3573: {
3574: PetscInt maxFields = PetscMax(1, numFields) + 1;
3575: PetscCall(PetscMalloc3(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &rowOffsets));
3576: }
3578: PetscCall(DMPlexTransferInjectorTree(coarse, fine, coarseToFine, childIds, NULL, numFields, offsets, &multiRootSec, &rootIndicesSec, &rootIndices, NULL));
3580: PetscCall(PetscMalloc1(maxDof, &parentIndices));
3582: /* count indices */
3583: PetscCall(MatGetLayouts(mat, &rowMap, &colMap));
3584: PetscCall(PetscLayoutSetUp(rowMap));
3585: PetscCall(PetscLayoutSetUp(colMap));
3586: PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd));
3587: PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd));
3588: PetscCall(PetscCalloc2(rowEnd - rowStart, &nnzD, rowEnd - rowStart, &nnzO));
3589: for (p = pStartC; p < pEndC; p++) {
3590: PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
3592: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3593: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
3594: if ((dof - cdof) <= 0) continue;
3595: PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));
3597: rowOffsets[0] = 0;
3598: offsetsCopy[0] = 0;
3599: if (numFields) {
3600: PetscInt f;
3602: for (f = 0; f < numFields; f++) {
3603: PetscInt fDof;
3604: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
3605: rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3606: }
3607: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices));
3608: } else {
3609: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices));
3610: rowOffsets[1] = offsetsCopy[0];
3611: }
3613: PetscCall(PetscSectionGetDof(multiRootSec, p, &numLeaves));
3614: PetscCall(PetscSectionGetOffset(multiRootSec, p, &leafStart));
3615: leafEnd = leafStart + numLeaves;
3616: for (l = leafStart; l < leafEnd; l++) {
3617: PetscInt numIndices, childId, offset;
3618: const PetscInt *childIndices;
3620: PetscCall(PetscSectionGetDof(rootIndicesSec, l, &numIndices));
3621: PetscCall(PetscSectionGetOffset(rootIndicesSec, l, &offset));
3622: childId = rootIndices[offset++];
3623: childIndices = &rootIndices[offset];
3624: numIndices--;
3626: if (childId == -1) { /* equivalent points: scatter */
3627: PetscInt i;
3629: for (i = 0; i < numIndices; i++) {
3630: PetscInt colIndex = childIndices[i];
3631: PetscInt rowIndex = parentIndices[i];
3632: if (rowIndex < 0) continue;
3633: PetscCheck(colIndex >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unconstrained fine and constrained coarse");
3634: if (colIndex >= colStart && colIndex < colEnd) {
3635: nnzD[rowIndex - rowStart] = 1;
3636: } else {
3637: nnzO[rowIndex - rowStart] = 1;
3638: }
3639: }
3640: } else {
3641: PetscInt parentId, f, lim;
3643: PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL));
3645: lim = PetscMax(1, numFields);
3646: offsets[0] = 0;
3647: if (numFields) {
3648: PetscInt f;
3650: for (f = 0; f < numFields; f++) {
3651: PetscInt fDof;
3652: PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof));
3654: offsets[f + 1] = fDof + offsets[f];
3655: }
3656: } else {
3657: PetscInt cDof;
3659: PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof));
3660: offsets[1] = cDof;
3661: }
3662: for (f = 0; f < lim; f++) {
3663: PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1];
3664: PetscInt childStart = offsets[f], childEnd = offsets[f + 1];
3665: PetscInt i, numD = 0, numO = 0;
3667: for (i = childStart; i < childEnd; i++) {
3668: PetscInt colIndex = childIndices[i];
3670: if (colIndex < 0) continue;
3671: if (colIndex >= colStart && colIndex < colEnd) {
3672: numD++;
3673: } else {
3674: numO++;
3675: }
3676: }
3677: for (i = parentStart; i < parentEnd; i++) {
3678: PetscInt rowIndex = parentIndices[i];
3680: if (rowIndex < 0) continue;
3681: nnzD[rowIndex - rowStart] += numD;
3682: nnzO[rowIndex - rowStart] += numO;
3683: }
3684: }
3685: }
3686: }
3687: }
3688: /* preallocate */
3689: PetscCall(MatXAIJSetPreallocation(mat, 1, nnzD, nnzO, NULL, NULL));
3690: PetscCall(PetscFree2(nnzD, nnzO));
3691: /* insert values */
3692: PetscCall(DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree, injRef, &childrenMats));
3693: for (p = pStartC; p < pEndC; p++) {
3694: PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
3696: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3697: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
3698: if ((dof - cdof) <= 0) continue;
3699: PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));
3701: rowOffsets[0] = 0;
3702: offsetsCopy[0] = 0;
3703: if (numFields) {
3704: PetscInt f;
3706: for (f = 0; f < numFields; f++) {
3707: PetscInt fDof;
3708: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
3709: rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3710: }
3711: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices));
3712: } else {
3713: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices));
3714: rowOffsets[1] = offsetsCopy[0];
3715: }
3717: PetscCall(PetscSectionGetDof(multiRootSec, p, &numLeaves));
3718: PetscCall(PetscSectionGetOffset(multiRootSec, p, &leafStart));
3719: leafEnd = leafStart + numLeaves;
3720: for (l = leafStart; l < leafEnd; l++) {
3721: PetscInt numIndices, childId, offset;
3722: const PetscInt *childIndices;
3724: PetscCall(PetscSectionGetDof(rootIndicesSec, l, &numIndices));
3725: PetscCall(PetscSectionGetOffset(rootIndicesSec, l, &offset));
3726: childId = rootIndices[offset++];
3727: childIndices = &rootIndices[offset];
3728: numIndices--;
3730: if (childId == -1) { /* equivalent points: scatter */
3731: PetscInt i;
3733: for (i = 0; i < numIndices; i++) PetscCall(MatSetValue(mat, parentIndices[i], childIndices[i], 1., INSERT_VALUES));
3734: } else {
3735: PetscInt parentId, f, lim;
3737: PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL));
3739: lim = PetscMax(1, numFields);
3740: offsets[0] = 0;
3741: if (numFields) {
3742: PetscInt f;
3744: for (f = 0; f < numFields; f++) {
3745: PetscInt fDof;
3746: PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof));
3748: offsets[f + 1] = fDof + offsets[f];
3749: }
3750: } else {
3751: PetscInt cDof;
3753: PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof));
3754: offsets[1] = cDof;
3755: }
3756: for (f = 0; f < lim; f++) {
3757: PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0];
3758: PetscInt *rowIndices = &parentIndices[rowOffsets[f]];
3759: const PetscInt *colIndices = &childIndices[offsets[f]];
3761: PetscCall(MatSetValues(mat, rowOffsets[f + 1] - rowOffsets[f], rowIndices, offsets[f + 1] - offsets[f], colIndices, childMat, INSERT_VALUES));
3762: }
3763: }
3764: }
3765: }
3766: PetscCall(PetscSectionDestroy(&multiRootSec));
3767: PetscCall(PetscSectionDestroy(&rootIndicesSec));
3768: PetscCall(PetscFree(parentIndices));
3769: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree, injRef, &childrenMats));
3770: PetscCall(PetscFree(rootIndices));
3771: PetscCall(PetscFree3(offsets, offsetsCopy, rowOffsets));
3773: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
3774: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
3775: PetscFunctionReturn(PETSC_SUCCESS);
3776: }
3778: static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom)
3779: {
3780: PetscSF coarseToFineEmbedded;
3781: PetscSection globalCoarse, globalFine;
3782: PetscSection localCoarse, localFine;
3783: PetscSection aSec, cSec;
3784: PetscSection rootValuesSec;
3785: PetscSection leafValuesSec;
3786: PetscScalar *rootValues, *leafValues;
3787: IS aIS;
3788: const PetscInt *anchors;
3789: Mat cMat;
3790: PetscInt numFields;
3791: PetscInt pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd;
3792: PetscInt aStart, aEnd, cStart, cEnd;
3793: PetscInt *maxChildIds;
3794: PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
3795: PetscFV fv = NULL;
3796: PetscInt dim, numFVcomps = -1, fvField = -1;
3797: DM cellDM = NULL, gradDM = NULL;
3798: const PetscScalar *cellGeomArray = NULL;
3799: const PetscScalar *gradArray = NULL;
3801: PetscFunctionBegin;
3802: PetscCall(VecSetOption(vecFine, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE));
3803: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
3804: PetscCall(DMPlexGetSimplexOrBoxCells(coarse, 0, &cellStart, &cellEnd));
3805: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
3806: PetscCall(DMGetGlobalSection(fine, &globalFine));
3807: PetscCall(DMGetCoordinateDim(coarse, &dim));
3808: { /* winnow fine points that don't have global dofs out of the sf */
3809: PetscInt nleaves, l;
3810: const PetscInt *leaves;
3811: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;
3813: PetscCall(PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL));
3815: for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
3816: PetscInt p = leaves ? leaves[l] : l;
3818: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
3819: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
3820: if ((dof - cdof) > 0) numPointsWithDofs++;
3821: }
3822: PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
3823: for (l = 0, offset = 0; l < nleaves; l++) {
3824: PetscInt p = leaves ? leaves[l] : l;
3826: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
3827: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
3828: if ((dof - cdof) > 0) pointsWithDofs[offset++] = l;
3829: }
3830: PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded));
3831: PetscCall(PetscFree(pointsWithDofs));
3832: }
3833: /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
3834: PetscCall(PetscMalloc1(pEndC - pStartC, &maxChildIds));
3835: for (p = pStartC; p < pEndC; p++) maxChildIds[p - pStartC] = -2;
3836: PetscCall(PetscSFReduceBegin(coarseToFineEmbedded, MPIU_INT, cids, maxChildIds, MPIU_MAX));
3837: PetscCall(PetscSFReduceEnd(coarseToFineEmbedded, MPIU_INT, cids, maxChildIds, MPIU_MAX));
3839: PetscCall(DMGetLocalSection(coarse, &localCoarse));
3840: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
3842: PetscCall(DMPlexGetAnchors(coarse, &aSec, &aIS));
3843: PetscCall(ISGetIndices(aIS, &anchors));
3844: PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd));
3846: PetscCall(DMGetDefaultConstraints(coarse, &cSec, &cMat, NULL));
3847: PetscCall(PetscSectionGetChart(cSec, &cStart, &cEnd));
3849: /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
3850: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootValuesSec));
3851: PetscCall(PetscSectionSetChart(rootValuesSec, pStartC, pEndC));
3852: PetscCall(PetscSectionGetNumFields(localCoarse, &numFields));
3853: {
3854: PetscInt maxFields = PetscMax(1, numFields) + 1;
3855: PetscCall(PetscMalloc7(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &newOffsets, maxFields, &newOffsetsCopy, maxFields, &rowOffsets, maxFields, &numD, maxFields, &numO));
3856: }
3857: if (grad) {
3858: PetscInt i;
3860: PetscCall(VecGetDM(cellGeom, &cellDM));
3861: PetscCall(VecGetArrayRead(cellGeom, &cellGeomArray));
3862: PetscCall(VecGetDM(grad, &gradDM));
3863: PetscCall(VecGetArrayRead(grad, &gradArray));
3864: for (i = 0; i < PetscMax(1, numFields); i++) {
3865: PetscObject obj;
3866: PetscClassId id;
3868: PetscCall(DMGetField(coarse, i, NULL, &obj));
3869: PetscCall(PetscObjectGetClassId(obj, &id));
3870: if (id == PETSCFV_CLASSID) {
3871: fv = (PetscFV)obj;
3872: PetscCall(PetscFVGetNumComponents(fv, &numFVcomps));
3873: fvField = i;
3874: break;
3875: }
3876: }
3877: }
3879: for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
3880: PetscInt dof;
3881: PetscInt maxChildId = maxChildIds[p - pStartC];
3882: PetscInt numValues = 0;
3884: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3885: if (dof < 0) dof = -(dof + 1);
3886: offsets[0] = 0;
3887: newOffsets[0] = 0;
3888: if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
3889: PetscInt *closure = NULL, closureSize, cl;
3891: PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
3892: for (cl = 0; cl < closureSize; cl++) { /* get the closure */
3893: PetscInt c = closure[2 * cl], clDof;
3895: PetscCall(PetscSectionGetDof(localCoarse, c, &clDof));
3896: numValues += clDof;
3897: }
3898: PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
3899: } else if (maxChildId == -1) {
3900: PetscCall(PetscSectionGetDof(localCoarse, p, &numValues));
3901: }
3902: /* we will pack the column indices with the field offsets */
3903: if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) {
3904: /* also send the centroid, and the gradient */
3905: numValues += dim * (1 + numFVcomps);
3906: }
3907: PetscCall(PetscSectionSetDof(rootValuesSec, p, numValues));
3908: }
3909: PetscCall(PetscSectionSetUp(rootValuesSec));
3910: {
3911: PetscInt numRootValues;
3912: const PetscScalar *coarseArray;
3914: PetscCall(PetscSectionGetStorageSize(rootValuesSec, &numRootValues));
3915: PetscCall(PetscMalloc1(numRootValues, &rootValues));
3916: PetscCall(VecGetArrayRead(vecCoarseLocal, &coarseArray));
3917: for (p = pStartC; p < pEndC; p++) {
3918: PetscInt numValues;
3919: PetscInt pValOff;
3920: PetscScalar *pVal;
3921: PetscInt maxChildId = maxChildIds[p - pStartC];
3923: PetscCall(PetscSectionGetDof(rootValuesSec, p, &numValues));
3924: if (!numValues) continue;
3925: PetscCall(PetscSectionGetOffset(rootValuesSec, p, &pValOff));
3926: pVal = &rootValues[pValOff];
3927: if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
3928: PetscInt closureSize = numValues;
3929: PetscCall(DMPlexVecGetClosure(coarse, NULL, vecCoarseLocal, p, &closureSize, &pVal));
3930: if (grad && p >= cellStart && p < cellEnd) {
3931: PetscFVCellGeom *cg;
3932: PetscScalar *gradVals = NULL;
3933: PetscInt i;
3935: pVal += (numValues - dim * (1 + numFVcomps));
3937: PetscCall(DMPlexPointLocalRead(cellDM, p, cellGeomArray, (void *)&cg));
3938: for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i];
3939: pVal += dim;
3940: PetscCall(DMPlexPointGlobalRead(gradDM, p, gradArray, (void *)&gradVals));
3941: for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i];
3942: }
3943: } else if (maxChildId == -1) {
3944: PetscInt lDof, lOff, i;
3946: PetscCall(PetscSectionGetDof(localCoarse, p, &lDof));
3947: PetscCall(PetscSectionGetOffset(localCoarse, p, &lOff));
3948: for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i];
3949: }
3950: }
3951: PetscCall(VecRestoreArrayRead(vecCoarseLocal, &coarseArray));
3952: PetscCall(PetscFree(maxChildIds));
3953: }
3954: {
3955: PetscSF valuesSF;
3956: PetscInt *remoteOffsetsValues, numLeafValues;
3958: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafValuesSec));
3959: PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootValuesSec, &remoteOffsetsValues, leafValuesSec));
3960: PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootValuesSec, remoteOffsetsValues, leafValuesSec, &valuesSF));
3961: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
3962: PetscCall(PetscFree(remoteOffsetsValues));
3963: PetscCall(PetscSectionGetStorageSize(leafValuesSec, &numLeafValues));
3964: PetscCall(PetscMalloc1(numLeafValues, &leafValues));
3965: PetscCall(PetscSFBcastBegin(valuesSF, MPIU_SCALAR, rootValues, leafValues, MPI_REPLACE));
3966: PetscCall(PetscSFBcastEnd(valuesSF, MPIU_SCALAR, rootValues, leafValues, MPI_REPLACE));
3967: PetscCall(PetscSFDestroy(&valuesSF));
3968: PetscCall(PetscFree(rootValues));
3969: PetscCall(PetscSectionDestroy(&rootValuesSec));
3970: }
3971: PetscCall(DMGetLocalSection(fine, &localFine));
3972: {
3973: PetscInt maxDof;
3974: PetscInt *rowIndices;
3975: DM refTree;
3976: PetscInt **refPointFieldN;
3977: PetscScalar ***refPointFieldMats;
3978: PetscSection refConSec, refAnSec;
3979: PetscInt pRefStart, pRefEnd, leafStart, leafEnd;
3980: PetscScalar *pointWork;
3982: PetscCall(PetscSectionGetMaxDof(localFine, &maxDof));
3983: PetscCall(DMGetWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
3984: PetscCall(DMGetWorkArray(fine, maxDof, MPIU_SCALAR, &pointWork));
3985: PetscCall(DMPlexGetReferenceTree(fine, &refTree));
3986: PetscCall(DMCopyDisc(fine, refTree));
3987: PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
3988: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
3989: PetscCall(DMPlexGetAnchors(refTree, &refAnSec, NULL));
3990: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
3991: PetscCall(PetscSectionGetChart(leafValuesSec, &leafStart, &leafEnd));
3992: PetscCall(DMPlexGetSimplexOrBoxCells(fine, 0, &cellStart, &cellEnd));
3993: for (p = leafStart; p < leafEnd; p++) {
3994: PetscInt gDof, gcDof, gOff, lDof;
3995: PetscInt numValues, pValOff;
3996: PetscInt childId;
3997: const PetscScalar *pVal;
3998: const PetscScalar *fvGradData = NULL;
4000: PetscCall(PetscSectionGetDof(globalFine, p, &gDof));
4001: PetscCall(PetscSectionGetDof(localFine, p, &lDof));
4002: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &gcDof));
4003: if ((gDof - gcDof) <= 0) continue;
4004: PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
4005: PetscCall(PetscSectionGetDof(leafValuesSec, p, &numValues));
4006: if (!numValues) continue;
4007: PetscCall(PetscSectionGetOffset(leafValuesSec, p, &pValOff));
4008: pVal = &leafValues[pValOff];
4009: offsets[0] = 0;
4010: offsetsCopy[0] = 0;
4011: newOffsets[0] = 0;
4012: newOffsetsCopy[0] = 0;
4013: childId = cids[p - pStartF];
4014: if (numFields) {
4015: PetscInt f;
4016: for (f = 0; f < numFields; f++) {
4017: PetscInt rowDof;
4019: PetscCall(PetscSectionGetFieldDof(localFine, p, f, &rowDof));
4020: offsets[f + 1] = offsets[f] + rowDof;
4021: offsetsCopy[f + 1] = offsets[f + 1];
4022: /* TODO: closure indices */
4023: newOffsets[f + 1] = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]);
4024: }
4025: PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices));
4026: } else {
4027: offsets[0] = 0;
4028: offsets[1] = lDof;
4029: newOffsets[0] = 0;
4030: newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0];
4031: PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices));
4032: }
4033: if (childId == -1) { /* no child interpolation: one nnz per */
4034: PetscCall(VecSetValues(vecFine, numValues, rowIndices, pVal, INSERT_VALUES));
4035: } else {
4036: PetscInt f;
4038: if (grad && p >= cellStart && p < cellEnd) {
4039: numValues -= (dim * (1 + numFVcomps));
4040: fvGradData = &pVal[numValues];
4041: }
4042: for (f = 0; f < PetscMax(1, numFields); f++) {
4043: const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f];
4044: PetscInt numRows = offsets[f + 1] - offsets[f];
4045: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
4046: const PetscScalar *cVal = &pVal[newOffsets[f]];
4047: PetscScalar *rVal = &pointWork[offsets[f]];
4048: PetscInt i, j;
4050: #if 0
4051: 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));
4052: #endif
4053: for (i = 0; i < numRows; i++) {
4054: PetscScalar val = 0.;
4055: for (j = 0; j < numCols; j++) val += childMat[i * numCols + j] * cVal[j];
4056: rVal[i] = val;
4057: }
4058: if (f == fvField && p >= cellStart && p < cellEnd) {
4059: PetscReal centroid[3];
4060: PetscScalar diff[3];
4061: const PetscScalar *parentCentroid = &fvGradData[0];
4062: const PetscScalar *gradient = &fvGradData[dim];
4064: PetscCall(DMPlexComputeCellGeometryFVM(fine, p, NULL, centroid, NULL));
4065: for (i = 0; i < dim; i++) diff[i] = centroid[i] - parentCentroid[i];
4066: for (i = 0; i < numFVcomps; i++) {
4067: PetscScalar val = 0.;
4069: for (j = 0; j < dim; j++) val += gradient[dim * i + j] * diff[j];
4070: rVal[i] += val;
4071: }
4072: }
4073: PetscCall(VecSetValues(vecFine, numRows, &rowIndices[offsets[f]], rVal, INSERT_VALUES));
4074: }
4075: }
4076: }
4077: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
4078: PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_SCALAR, &pointWork));
4079: PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
4080: }
4081: PetscCall(PetscFree(leafValues));
4082: PetscCall(PetscSectionDestroy(&leafValuesSec));
4083: PetscCall(PetscFree7(offsets, offsetsCopy, newOffsets, newOffsetsCopy, rowOffsets, numD, numO));
4084: PetscCall(ISRestoreIndices(aIS, &anchors));
4085: PetscFunctionReturn(PETSC_SUCCESS);
4086: }
4088: static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids)
4089: {
4090: DM refTree;
4091: PetscSection multiRootSec, rootIndicesSec;
4092: PetscSection globalCoarse, globalFine;
4093: PetscSection localCoarse, localFine;
4094: PetscSection cSecRef;
4095: PetscInt *parentIndices, pRefStart, pRefEnd;
4096: PetscScalar *rootValues, *parentValues;
4097: Mat injRef;
4098: PetscInt numFields, maxDof;
4099: PetscInt pStartC, pEndC, pStartF, pEndF, p;
4100: PetscInt *offsets, *offsetsCopy, *rowOffsets;
4101: PetscLayout rowMap, colMap;
4102: PetscInt rowStart, rowEnd, colStart, colEnd;
4103: PetscScalar ***childrenMats = NULL; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */
4105: PetscFunctionBegin;
4106: /* get the templates for the fine-to-coarse injection from the reference tree */
4107: PetscCall(VecSetOption(vecFine, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE));
4108: PetscCall(VecSetOption(vecCoarse, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE));
4109: PetscCall(DMPlexGetReferenceTree(coarse, &refTree));
4110: PetscCall(DMCopyDisc(coarse, refTree));
4111: PetscCall(DMGetDefaultConstraints(refTree, &cSecRef, NULL, NULL));
4112: PetscCall(PetscSectionGetChart(cSecRef, &pRefStart, &pRefEnd));
4113: PetscCall(DMPlexReferenceTreeGetInjector(refTree, &injRef));
4115: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
4116: PetscCall(DMGetLocalSection(fine, &localFine));
4117: PetscCall(DMGetGlobalSection(fine, &globalFine));
4118: PetscCall(PetscSectionGetNumFields(localFine, &numFields));
4119: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
4120: PetscCall(DMGetLocalSection(coarse, &localCoarse));
4121: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
4122: PetscCall(PetscSectionGetMaxDof(localCoarse, &maxDof));
4123: {
4124: PetscInt maxFields = PetscMax(1, numFields) + 1;
4125: PetscCall(PetscMalloc3(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &rowOffsets));
4126: }
4128: PetscCall(DMPlexTransferInjectorTree(coarse, fine, coarseToFine, cids, vecFine, numFields, offsets, &multiRootSec, &rootIndicesSec, NULL, &rootValues));
4130: PetscCall(PetscMalloc2(maxDof, &parentIndices, maxDof, &parentValues));
4132: /* count indices */
4133: PetscCall(VecGetLayout(vecFine, &colMap));
4134: PetscCall(VecGetLayout(vecCoarse, &rowMap));
4135: PetscCall(PetscLayoutSetUp(rowMap));
4136: PetscCall(PetscLayoutSetUp(colMap));
4137: PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd));
4138: PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd));
4139: /* insert values */
4140: PetscCall(DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree, injRef, &childrenMats));
4141: for (p = pStartC; p < pEndC; p++) {
4142: PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
4143: PetscBool contribute = PETSC_FALSE;
4145: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
4146: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
4147: if ((dof - cdof) <= 0) continue;
4148: PetscCall(PetscSectionGetDof(localCoarse, p, &dof));
4149: PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));
4151: rowOffsets[0] = 0;
4152: offsetsCopy[0] = 0;
4153: if (numFields) {
4154: PetscInt f;
4156: for (f = 0; f < numFields; f++) {
4157: PetscInt fDof;
4158: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
4159: rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
4160: }
4161: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices));
4162: } else {
4163: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices));
4164: rowOffsets[1] = offsetsCopy[0];
4165: }
4167: PetscCall(PetscSectionGetDof(multiRootSec, p, &numLeaves));
4168: PetscCall(PetscSectionGetOffset(multiRootSec, p, &leafStart));
4169: leafEnd = leafStart + numLeaves;
4170: for (l = 0; l < dof; l++) parentValues[l] = 0.;
4171: for (l = leafStart; l < leafEnd; l++) {
4172: PetscInt numIndices, childId, offset;
4173: const PetscScalar *childValues;
4175: PetscCall(PetscSectionGetDof(rootIndicesSec, l, &numIndices));
4176: PetscCall(PetscSectionGetOffset(rootIndicesSec, l, &offset));
4177: childId = (PetscInt)PetscRealPart(rootValues[offset++]);
4178: childValues = &rootValues[offset];
4179: numIndices--;
4181: if (childId == -2) { /* skip */
4182: continue;
4183: } else if (childId == -1) { /* equivalent points: scatter */
4184: PetscInt m;
4186: contribute = PETSC_TRUE;
4187: for (m = 0; m < numIndices; m++) parentValues[m] = childValues[m];
4188: } else { /* contributions from children: sum with injectors from reference tree */
4189: PetscInt parentId, f, lim;
4191: contribute = PETSC_TRUE;
4192: PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL));
4194: lim = PetscMax(1, numFields);
4195: offsets[0] = 0;
4196: if (numFields) {
4197: PetscInt f;
4199: for (f = 0; f < numFields; f++) {
4200: PetscInt fDof;
4201: PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof));
4203: offsets[f + 1] = fDof + offsets[f];
4204: }
4205: } else {
4206: PetscInt cDof;
4208: PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof));
4209: offsets[1] = cDof;
4210: }
4211: for (f = 0; f < lim; f++) {
4212: PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0];
4213: PetscInt n = offsets[f + 1] - offsets[f];
4214: PetscInt m = rowOffsets[f + 1] - rowOffsets[f];
4215: PetscInt i, j;
4216: const PetscScalar *colValues = &childValues[offsets[f]];
4218: for (i = 0; i < m; i++) {
4219: PetscScalar val = 0.;
4220: for (j = 0; j < n; j++) val += childMat[n * i + j] * colValues[j];
4221: parentValues[rowOffsets[f] + i] += val;
4222: }
4223: }
4224: }
4225: }
4226: if (contribute) PetscCall(VecSetValues(vecCoarse, dof, parentIndices, parentValues, INSERT_VALUES));
4227: }
4228: PetscCall(PetscSectionDestroy(&multiRootSec));
4229: PetscCall(PetscSectionDestroy(&rootIndicesSec));
4230: PetscCall(PetscFree2(parentIndices, parentValues));
4231: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree, injRef, &childrenMats));
4232: PetscCall(PetscFree(rootValues));
4233: PetscCall(PetscFree3(offsets, offsetsCopy, rowOffsets));
4234: PetscFunctionReturn(PETSC_SUCCESS);
4235: }
4237: /*@
4238: DMPlexTransferVecTree - transfer a vector between two meshes that differ from each other by refinement/coarsening
4239: that can be represented by a common reference tree used by both. This routine can be used for a combination of
4240: coarsening and refinement at the same time.
4242: Collective
4244: Input Parameters:
4245: + dmIn - The `DMPLEX` mesh for the input vector
4246: . dmOut - The second `DMPLEX` mesh
4247: . vecIn - The input vector
4248: . sfRefine - A star forest indicating points in the mesh `dmIn` (roots in the star forest) that are parents to points in
4249: the mesh `dmOut` (leaves in the star forest), i.e. where `dmOut` is more refined than `dmIn`
4250: . sfCoarsen - A star forest indicating points in the mesh `dmOut` (roots in the star forest) that are parents to points in
4251: the mesh `dmIn` (leaves in the star forest), i.e. where `dmOut` is more coarsened than `dmIn`
4252: . cidsRefine - The childIds of the points in `dmOut`. These childIds relate back to the reference tree: childid[j] = k implies
4253: that mesh point j of `dmOut` was refined from a point in `dmIn` just as the mesh point k in the reference
4254: tree was refined from its parent. childid[j] = -1 indicates that the point j in `dmOut` is exactly
4255: equivalent to its root in `dmIn`, so no interpolation is necessary. childid[j] = -2 indicates that this
4256: point j in `dmOut` is not a leaf of `sfRefine`.
4257: . cidsCoarsen - The childIds of the points in `dmIn`. These childIds relate back to the reference tree: childid[j] = k implies
4258: that mesh point j of dmIn coarsens to a point in `dmOut` just as the mesh point k in the reference
4259: tree coarsens to its parent. childid[j] = -2 indicates that point j in `dmOut` is not a leaf in `sfCoarsen`.
4260: . useBCs - `PETSC_TRUE` indicates that boundary values should be inserted into `vecIn` before transfer.
4261: - time - Used if boundary values are time dependent.
4263: Output Parameter:
4264: . vecOut - Using interpolation and injection operators calculated on the reference tree, the transferred
4265: projection of `vecIn` from `dmIn` to `dmOut`. Note that any field discretized with a `PetscFV` finite volume
4266: method that uses gradient reconstruction will use reconstructed gradients when interpolating from
4267: coarse points to fine points.
4269: Level: developer
4271: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `Vec`, `PetscFV`, `DMPlexSetReferenceTree()`, `DMPlexGetReferenceTree()`, `PetscFVGetComputeGradients()`
4272: @*/
4273: PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time)
4274: {
4275: PetscFunctionBegin;
4276: PetscCall(VecSet(vecOut, 0.0));
4277: if (sfRefine) {
4278: Vec vecInLocal;
4279: DM dmGrad = NULL;
4280: Vec faceGeom = NULL, cellGeom = NULL, grad = NULL;
4282: PetscCall(DMGetLocalVector(dmIn, &vecInLocal));
4283: PetscCall(VecSet(vecInLocal, 0.0));
4284: {
4285: PetscInt numFields, i;
4287: PetscCall(DMGetNumFields(dmIn, &numFields));
4288: for (i = 0; i < numFields; i++) {
4289: PetscObject obj;
4290: PetscClassId classid;
4292: PetscCall(DMGetField(dmIn, i, NULL, &obj));
4293: PetscCall(PetscObjectGetClassId(obj, &classid));
4294: if (classid == PETSCFV_CLASSID) {
4295: PetscCall(DMPlexGetDataFVM(dmIn, (PetscFV)obj, &cellGeom, &faceGeom, &dmGrad));
4296: break;
4297: }
4298: }
4299: }
4300: if (useBCs) PetscCall(DMPlexInsertBoundaryValues(dmIn, PETSC_TRUE, vecInLocal, time, faceGeom, cellGeom, NULL));
4301: PetscCall(DMGlobalToLocalBegin(dmIn, vecIn, INSERT_VALUES, vecInLocal));
4302: PetscCall(DMGlobalToLocalEnd(dmIn, vecIn, INSERT_VALUES, vecInLocal));
4303: if (dmGrad) {
4304: PetscCall(DMGetGlobalVector(dmGrad, &grad));
4305: PetscCall(DMPlexReconstructGradientsFVM(dmIn, vecInLocal, grad));
4306: }
4307: PetscCall(DMPlexTransferVecTree_Interpolate(dmIn, vecInLocal, dmOut, vecOut, sfRefine, cidsRefine, grad, cellGeom));
4308: PetscCall(DMRestoreLocalVector(dmIn, &vecInLocal));
4309: if (dmGrad) PetscCall(DMRestoreGlobalVector(dmGrad, &grad));
4310: }
4311: if (sfCoarsen) PetscCall(DMPlexTransferVecTree_Inject(dmIn, vecIn, dmOut, vecOut, sfCoarsen, cidsCoarsen));
4312: PetscCall(VecAssemblyBegin(vecOut));
4313: PetscCall(VecAssemblyEnd(vecOut));
4314: PetscFunctionReturn(PETSC_SUCCESS);
4315: }