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: /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of
1751: * a non-conforming mesh. Local refinement comes later */
1752: PetscErrorCode DMPlexTreeRefineCell(DM dm, PetscInt cell, DM *ncdm)
1753: {
1754: DM K;
1755: PetscMPIInt rank;
1756: PetscInt dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd;
1757: PetscInt numNewCones, *newConeSizes, *newCones, *newOrientations;
1758: PetscInt *Kembedding;
1759: PetscInt *cellClosure = NULL, nc;
1760: PetscScalar *newVertexCoords;
1761: PetscInt numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset;
1762: PetscSection parentSection;
1764: PetscFunctionBegin;
1765: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1766: PetscCall(DMGetDimension(dm, &dim));
1767: PetscCall(DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm));
1768: PetscCall(DMSetDimension(*ncdm, dim));
1770: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1771: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &parentSection));
1772: PetscCall(DMPlexGetReferenceTree(dm, &K));
1773: PetscCall(DMGetCoordinatesLocalSetUp(dm));
1774: if (rank == 0) {
1775: /* compute the new charts */
1776: PetscCall(PetscMalloc5(dim + 1, &pNewCount, dim + 1, &pNewStart, dim + 1, &pNewEnd, dim + 1, &pOldStart, dim + 1, &pOldEnd));
1777: offset = 0;
1778: for (d = 0; d <= dim; d++) {
1779: PetscInt pOldCount, kStart, kEnd, k;
1781: pNewStart[d] = offset;
1782: PetscCall(DMPlexGetHeightStratum(dm, d, &pOldStart[d], &pOldEnd[d]));
1783: PetscCall(DMPlexGetHeightStratum(K, d, &kStart, &kEnd));
1784: pOldCount = pOldEnd[d] - pOldStart[d];
1785: /* adding the new points */
1786: pNewCount[d] = pOldCount + kEnd - kStart;
1787: if (!d) {
1788: /* removing the cell */
1789: pNewCount[d]--;
1790: }
1791: for (k = kStart; k < kEnd; k++) {
1792: PetscInt parent;
1793: PetscCall(DMPlexGetTreeParent(K, k, &parent, NULL));
1794: if (parent == k) {
1795: /* avoid double counting points that won't actually be new */
1796: pNewCount[d]--;
1797: }
1798: }
1799: pNewEnd[d] = pNewStart[d] + pNewCount[d];
1800: offset = pNewEnd[d];
1801: }
1802: 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]);
1803: /* get the current closure of the cell that we are removing */
1804: PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &nc, &cellClosure));
1806: PetscCall(PetscMalloc1(pNewEnd[dim], &newConeSizes));
1807: {
1808: DMPolytopeType pct, qct;
1809: PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j;
1811: PetscCall(DMPlexGetChart(K, &kStart, &kEnd));
1812: PetscCall(PetscMalloc4(kEnd - kStart, &Kembedding, kEnd - kStart, &perm, kEnd - kStart, &iperm, kEnd - kStart, &preOrient));
1814: for (k = kStart; k < kEnd; k++) {
1815: perm[k - kStart] = k;
1816: iperm[k - kStart] = k - kStart;
1817: preOrient[k - kStart] = 0;
1818: }
1820: PetscCall(DMPlexGetTransitiveClosure(K, 0, PETSC_TRUE, &closureSizeK, &closureK));
1821: for (j = 1; j < closureSizeK; j++) {
1822: PetscInt parentOrientA = closureK[2 * j + 1];
1823: PetscInt parentOrientB = cellClosure[2 * j + 1];
1824: PetscInt p, q;
1826: p = closureK[2 * j];
1827: q = cellClosure[2 * j];
1828: PetscCall(DMPlexGetCellType(K, p, &pct));
1829: PetscCall(DMPlexGetCellType(dm, q, &qct));
1830: for (d = 0; d <= dim; d++) {
1831: if (q >= pOldStart[d] && q < pOldEnd[d]) Kembedding[p] = (q - pOldStart[d]) + pNewStart[d];
1832: }
1833: parentOrientA = DMPolytopeConvertNewOrientation_Internal(pct, parentOrientA);
1834: parentOrientB = DMPolytopeConvertNewOrientation_Internal(qct, parentOrientB);
1835: if (parentOrientA != parentOrientB) {
1836: PetscInt numChildren, i;
1837: const PetscInt *children;
1839: PetscCall(DMPlexGetTreeChildren(K, p, &numChildren, &children));
1840: for (i = 0; i < numChildren; i++) {
1841: PetscInt kPerm, oPerm;
1843: k = children[i];
1844: PetscCall(DMPlexReferenceTreeGetChildSymmetry(K, p, parentOrientA, 0, k, parentOrientB, &oPerm, &kPerm));
1845: /* perm = what refTree position I'm in */
1846: perm[kPerm - kStart] = k;
1847: /* iperm = who is at this position */
1848: iperm[k - kStart] = kPerm - kStart;
1849: preOrient[kPerm - kStart] = oPerm;
1850: }
1851: }
1852: }
1853: PetscCall(DMPlexRestoreTransitiveClosure(K, 0, PETSC_TRUE, &closureSizeK, &closureK));
1854: }
1855: PetscCall(PetscSectionSetChart(parentSection, 0, pNewEnd[dim]));
1856: offset = 0;
1857: numNewCones = 0;
1858: for (d = 0; d <= dim; d++) {
1859: PetscInt kStart, kEnd, k;
1860: PetscInt p;
1861: PetscInt size;
1863: for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1864: /* skip cell 0 */
1865: if (p == cell) continue;
1866: /* old cones to new cones */
1867: PetscCall(DMPlexGetConeSize(dm, p, &size));
1868: newConeSizes[offset++] = size;
1869: numNewCones += size;
1870: }
1872: PetscCall(DMPlexGetHeightStratum(K, d, &kStart, &kEnd));
1873: for (k = kStart; k < kEnd; k++) {
1874: PetscInt kParent;
1876: PetscCall(DMPlexGetTreeParent(K, k, &kParent, NULL));
1877: if (kParent != k) {
1878: Kembedding[k] = offset;
1879: PetscCall(DMPlexGetConeSize(K, k, &size));
1880: newConeSizes[offset++] = size;
1881: numNewCones += size;
1882: if (kParent != 0) PetscCall(PetscSectionSetDof(parentSection, Kembedding[k], 1));
1883: }
1884: }
1885: }
1887: PetscCall(PetscSectionSetUp(parentSection));
1888: PetscCall(PetscSectionGetStorageSize(parentSection, &numPointsWithParents));
1889: PetscCall(PetscMalloc2(numNewCones, &newCones, numNewCones, &newOrientations));
1890: PetscCall(PetscMalloc2(numPointsWithParents, &parents, numPointsWithParents, &childIDs));
1892: /* fill new cones */
1893: offset = 0;
1894: for (d = 0; d <= dim; d++) {
1895: PetscInt kStart, kEnd, k, l;
1896: PetscInt p;
1897: PetscInt size;
1898: const PetscInt *cone, *orientation;
1900: for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1901: /* skip cell 0 */
1902: if (p == cell) continue;
1903: /* old cones to new cones */
1904: PetscCall(DMPlexGetConeSize(dm, p, &size));
1905: PetscCall(DMPlexGetCone(dm, p, &cone));
1906: PetscCall(DMPlexGetConeOrientation(dm, p, &orientation));
1907: for (l = 0; l < size; l++) {
1908: newCones[offset] = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1];
1909: newOrientations[offset++] = orientation[l];
1910: }
1911: }
1913: PetscCall(DMPlexGetHeightStratum(K, d, &kStart, &kEnd));
1914: for (k = kStart; k < kEnd; k++) {
1915: PetscInt kPerm = perm[k], kParent;
1916: PetscInt preO = preOrient[k];
1918: PetscCall(DMPlexGetTreeParent(K, k, &kParent, NULL));
1919: if (kParent != k) {
1920: /* embed new cones */
1921: PetscCall(DMPlexGetConeSize(K, k, &size));
1922: PetscCall(DMPlexGetCone(K, kPerm, &cone));
1923: PetscCall(DMPlexGetConeOrientation(K, kPerm, &orientation));
1924: for (l = 0; l < size; l++) {
1925: PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size - (preO + 1) - l) % size);
1926: PetscInt newO, lSize, oTrue;
1927: DMPolytopeType ct = DM_NUM_POLYTOPES;
1929: q = iperm[cone[m]];
1930: newCones[offset] = Kembedding[q];
1931: PetscCall(DMPlexGetConeSize(K, q, &lSize));
1932: if (lSize == 2) ct = DM_POLYTOPE_SEGMENT;
1933: else if (lSize == 4) ct = DM_POLYTOPE_QUADRILATERAL;
1934: oTrue = DMPolytopeConvertNewOrientation_Internal(ct, orientation[m]);
1935: oTrue = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2);
1936: newO = DihedralCompose(lSize, oTrue, preOrient[q]);
1937: newOrientations[offset++] = DMPolytopeConvertOldOrientation_Internal(ct, newO);
1938: }
1939: if (kParent != 0) {
1940: PetscInt newPoint = Kembedding[kParent];
1941: PetscCall(PetscSectionGetOffset(parentSection, Kembedding[k], &pOffset));
1942: parents[pOffset] = newPoint;
1943: childIDs[pOffset] = k;
1944: }
1945: }
1946: }
1947: }
1949: PetscCall(PetscMalloc1(dim * (pNewEnd[dim] - pNewStart[dim]), &newVertexCoords));
1951: /* fill coordinates */
1952: offset = 0;
1953: {
1954: PetscInt kStart, kEnd, l;
1955: PetscSection vSection;
1956: PetscInt v;
1957: Vec coords;
1958: PetscScalar *coordvals;
1959: PetscInt dof, off;
1960: PetscReal v0[3], J[9], detJ;
1962: if (PetscDefined(USE_DEBUG)) {
1963: PetscInt k;
1964: PetscCall(DMPlexGetHeightStratum(K, 0, &kStart, &kEnd));
1965: for (k = kStart; k < kEnd; k++) {
1966: PetscCall(DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ));
1967: PetscCheck(detJ > 0., PETSC_COMM_SELF, PETSC_ERR_PLIB, "reference tree cell %" PetscInt_FMT " has bad determinant", k);
1968: }
1969: }
1970: PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ));
1971: PetscCall(DMGetCoordinateSection(dm, &vSection));
1972: PetscCall(DMGetCoordinatesLocal(dm, &coords));
1973: PetscCall(VecGetArray(coords, &coordvals));
1974: for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) {
1975: PetscCall(PetscSectionGetDof(vSection, v, &dof));
1976: PetscCall(PetscSectionGetOffset(vSection, v, &off));
1977: for (l = 0; l < dof; l++) newVertexCoords[offset++] = coordvals[off + l];
1978: }
1979: PetscCall(VecRestoreArray(coords, &coordvals));
1981: PetscCall(DMGetCoordinateSection(K, &vSection));
1982: PetscCall(DMGetCoordinatesLocal(K, &coords));
1983: PetscCall(VecGetArray(coords, &coordvals));
1984: PetscCall(DMPlexGetDepthStratum(K, 0, &kStart, &kEnd));
1985: for (v = kStart; v < kEnd; v++) {
1986: PetscReal coord[3], newCoord[3];
1987: PetscInt vPerm = perm[v];
1988: PetscInt kParent;
1989: const PetscReal xi0[3] = {-1., -1., -1.};
1991: PetscCall(DMPlexGetTreeParent(K, v, &kParent, NULL));
1992: if (kParent != v) {
1993: /* this is a new vertex */
1994: PetscCall(PetscSectionGetOffset(vSection, vPerm, &off));
1995: for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off + l]);
1996: CoordinatesRefToReal(dim, dim, xi0, v0, J, coord, newCoord);
1997: for (l = 0; l < dim; ++l) newVertexCoords[offset + l] = newCoord[l];
1998: offset += dim;
1999: }
2000: }
2001: PetscCall(VecRestoreArray(coords, &coordvals));
2002: }
2004: /* need to reverse the order of pNewCount: vertices first, cells last */
2005: for (d = 0; d < (dim + 1) / 2; d++) {
2006: PetscInt tmp;
2008: tmp = pNewCount[d];
2009: pNewCount[d] = pNewCount[dim - d];
2010: pNewCount[dim - d] = tmp;
2011: }
2013: PetscCall(DMPlexCreateFromDAG(*ncdm, dim, pNewCount, newConeSizes, newCones, newOrientations, newVertexCoords));
2014: PetscCall(DMPlexSetReferenceTree(*ncdm, K));
2015: PetscCall(DMPlexSetTree(*ncdm, parentSection, parents, childIDs));
2017: /* clean up */
2018: PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &nc, &cellClosure));
2019: PetscCall(PetscFree5(pNewCount, pNewStart, pNewEnd, pOldStart, pOldEnd));
2020: PetscCall(PetscFree(newConeSizes));
2021: PetscCall(PetscFree2(newCones, newOrientations));
2022: PetscCall(PetscFree(newVertexCoords));
2023: PetscCall(PetscFree2(parents, childIDs));
2024: PetscCall(PetscFree4(Kembedding, perm, iperm, preOrient));
2025: } else {
2026: PetscInt p, counts[4];
2027: PetscInt *coneSizes, *cones, *orientations;
2028: Vec coordVec;
2029: PetscScalar *coords;
2031: for (d = 0; d <= dim; d++) {
2032: PetscInt dStart, dEnd;
2034: PetscCall(DMPlexGetDepthStratum(dm, d, &dStart, &dEnd));
2035: counts[d] = dEnd - dStart;
2036: }
2037: PetscCall(PetscMalloc1(pEnd - pStart, &coneSizes));
2038: for (p = pStart; p < pEnd; p++) PetscCall(DMPlexGetConeSize(dm, p, &coneSizes[p - pStart]));
2039: PetscCall(DMPlexGetCones(dm, &cones));
2040: PetscCall(DMPlexGetConeOrientations(dm, &orientations));
2041: PetscCall(DMGetCoordinatesLocal(dm, &coordVec));
2042: PetscCall(VecGetArray(coordVec, &coords));
2044: PetscCall(PetscSectionSetChart(parentSection, pStart, pEnd));
2045: PetscCall(PetscSectionSetUp(parentSection));
2046: PetscCall(DMPlexCreateFromDAG(*ncdm, dim, counts, coneSizes, cones, orientations, NULL));
2047: PetscCall(DMPlexSetReferenceTree(*ncdm, K));
2048: PetscCall(DMPlexSetTree(*ncdm, parentSection, NULL, NULL));
2049: PetscCall(VecRestoreArray(coordVec, &coords));
2050: }
2051: PetscCall(PetscSectionDestroy(&parentSection));
2052: PetscFunctionReturn(PETSC_SUCCESS);
2053: }
2055: PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
2056: {
2057: PetscSF coarseToFineEmbedded;
2058: PetscSection globalCoarse, globalFine;
2059: PetscSection localCoarse, localFine;
2060: PetscSection aSec, cSec;
2061: PetscSection rootIndicesSec, rootMatricesSec;
2062: PetscSection leafIndicesSec, leafMatricesSec;
2063: PetscInt *rootIndices, *leafIndices;
2064: PetscScalar *rootMatrices, *leafMatrices;
2065: IS aIS;
2066: const PetscInt *anchors;
2067: Mat cMat;
2068: PetscInt numFields, maxFields;
2069: PetscInt pStartC, pEndC, pStartF, pEndF, p;
2070: PetscInt aStart, aEnd, cStart, cEnd;
2071: PetscInt *maxChildIds;
2072: PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
2073: const PetscInt ***perms;
2074: const PetscScalar ***flips;
2076: PetscFunctionBegin;
2077: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
2078: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
2079: PetscCall(DMGetGlobalSection(fine, &globalFine));
2080: { /* winnow fine points that don't have global dofs out of the sf */
2081: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, nleaves, l;
2082: const PetscInt *leaves;
2084: PetscCall(PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL));
2085: for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
2086: p = leaves ? leaves[l] : l;
2087: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
2088: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
2089: if ((dof - cdof) > 0) numPointsWithDofs++;
2090: }
2091: PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
2092: for (l = 0, offset = 0; l < nleaves; l++) {
2093: p = leaves ? leaves[l] : l;
2094: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
2095: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
2096: if ((dof - cdof) > 0) pointsWithDofs[offset++] = l;
2097: }
2098: PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded));
2099: PetscCall(PetscFree(pointsWithDofs));
2100: }
2101: /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
2102: PetscCall(PetscMalloc1(pEndC - pStartC, &maxChildIds));
2103: for (p = pStartC; p < pEndC; p++) maxChildIds[p - pStartC] = -2;
2104: PetscCall(PetscSFReduceBegin(coarseToFineEmbedded, MPIU_INT, childIds, maxChildIds, MPI_MAX));
2105: PetscCall(PetscSFReduceEnd(coarseToFineEmbedded, MPIU_INT, childIds, maxChildIds, MPI_MAX));
2107: PetscCall(DMGetLocalSection(coarse, &localCoarse));
2108: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
2110: PetscCall(DMPlexGetAnchors(coarse, &aSec, &aIS));
2111: PetscCall(ISGetIndices(aIS, &anchors));
2112: PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd));
2114: PetscCall(DMGetDefaultConstraints(coarse, &cSec, &cMat, NULL));
2115: PetscCall(PetscSectionGetChart(cSec, &cStart, &cEnd));
2117: /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
2118: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootIndicesSec));
2119: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootMatricesSec));
2120: PetscCall(PetscSectionSetChart(rootIndicesSec, pStartC, pEndC));
2121: PetscCall(PetscSectionSetChart(rootMatricesSec, pStartC, pEndC));
2122: PetscCall(PetscSectionGetNumFields(localCoarse, &numFields));
2123: maxFields = PetscMax(1, numFields);
2124: PetscCall(PetscMalloc7(maxFields + 1, &offsets, maxFields + 1, &offsetsCopy, maxFields + 1, &newOffsets, maxFields + 1, &newOffsetsCopy, maxFields + 1, &rowOffsets, maxFields + 1, &numD, maxFields + 1, &numO));
2125: PetscCall(PetscMalloc2(maxFields + 1, (PetscInt ****)&perms, maxFields + 1, (PetscScalar ****)&flips));
2126: PetscCall(PetscMemzero((void *)perms, (maxFields + 1) * sizeof(const PetscInt **)));
2127: PetscCall(PetscMemzero((void *)flips, (maxFields + 1) * sizeof(const PetscScalar **)));
2129: for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
2130: PetscInt dof, matSize = 0;
2131: PetscInt aDof = 0;
2132: PetscInt cDof = 0;
2133: PetscInt maxChildId = maxChildIds[p - pStartC];
2134: PetscInt numRowIndices = 0;
2135: PetscInt numColIndices = 0;
2136: PetscInt f;
2138: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
2139: if (dof < 0) dof = -(dof + 1);
2140: if (p >= aStart && p < aEnd) PetscCall(PetscSectionGetDof(aSec, p, &aDof));
2141: if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(cSec, p, &cDof));
2142: for (f = 0; f <= numFields; f++) offsets[f] = 0;
2143: for (f = 0; f <= numFields; f++) newOffsets[f] = 0;
2144: if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
2145: PetscInt *closure = NULL, closureSize, cl;
2147: PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2148: for (cl = 0; cl < closureSize; cl++) { /* get the closure */
2149: PetscInt c = closure[2 * cl], clDof;
2151: PetscCall(PetscSectionGetDof(localCoarse, c, &clDof));
2152: numRowIndices += clDof;
2153: for (f = 0; f < numFields; f++) {
2154: PetscCall(PetscSectionGetFieldDof(localCoarse, c, f, &clDof));
2155: offsets[f + 1] += clDof;
2156: }
2157: }
2158: for (f = 0; f < numFields; f++) {
2159: offsets[f + 1] += offsets[f];
2160: newOffsets[f + 1] = offsets[f + 1];
2161: }
2162: /* get the number of indices needed and their field offsets */
2163: PetscCall(DMPlexAnchorsModifyMat(coarse, localCoarse, closureSize, numRowIndices, closure, NULL, NULL, NULL, &numColIndices, NULL, NULL, newOffsets, PETSC_FALSE));
2164: PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2165: if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */
2166: numColIndices = numRowIndices;
2167: matSize = 0;
2168: } else if (numFields) { /* we send one submat for each field: sum their sizes */
2169: matSize = 0;
2170: for (f = 0; f < numFields; f++) {
2171: PetscInt numRow, numCol;
2173: numRow = offsets[f + 1] - offsets[f];
2174: numCol = newOffsets[f + 1] - newOffsets[f];
2175: matSize += numRow * numCol;
2176: }
2177: } else {
2178: matSize = numRowIndices * numColIndices;
2179: }
2180: } else if (maxChildId == -1) {
2181: if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */
2182: PetscInt aOff, a;
2184: PetscCall(PetscSectionGetOffset(aSec, p, &aOff));
2185: for (f = 0; f < numFields; f++) {
2186: PetscInt fDof;
2188: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
2189: offsets[f + 1] = fDof;
2190: }
2191: for (a = 0; a < aDof; a++) {
2192: PetscInt anchor = anchors[a + aOff], aLocalDof;
2194: PetscCall(PetscSectionGetDof(localCoarse, anchor, &aLocalDof));
2195: numColIndices += aLocalDof;
2196: for (f = 0; f < numFields; f++) {
2197: PetscInt fDof;
2199: PetscCall(PetscSectionGetFieldDof(localCoarse, anchor, f, &fDof));
2200: newOffsets[f + 1] += fDof;
2201: }
2202: }
2203: if (numFields) {
2204: matSize = 0;
2205: for (f = 0; f < numFields; f++) matSize += offsets[f + 1] * newOffsets[f + 1];
2206: } else {
2207: matSize = numColIndices * dof;
2208: }
2209: } else { /* no children, and no constraints on dofs: just get the global indices */
2210: numColIndices = dof;
2211: matSize = 0;
2212: }
2213: }
2214: /* we will pack the column indices with the field offsets */
2215: PetscCall(PetscSectionSetDof(rootIndicesSec, p, numColIndices ? numColIndices + 2 * numFields : 0));
2216: PetscCall(PetscSectionSetDof(rootMatricesSec, p, matSize));
2217: }
2218: PetscCall(PetscSectionSetUp(rootIndicesSec));
2219: PetscCall(PetscSectionSetUp(rootMatricesSec));
2220: {
2221: PetscInt numRootIndices, numRootMatrices;
2223: PetscCall(PetscSectionGetStorageSize(rootIndicesSec, &numRootIndices));
2224: PetscCall(PetscSectionGetStorageSize(rootMatricesSec, &numRootMatrices));
2225: PetscCall(PetscMalloc2(numRootIndices, &rootIndices, numRootMatrices, &rootMatrices));
2226: for (p = pStartC; p < pEndC; p++) {
2227: PetscInt numRowIndices = 0, numColIndices, matSize, dof;
2228: PetscInt pIndOff, pMatOff, f;
2229: PetscInt *pInd;
2230: PetscInt maxChildId = maxChildIds[p - pStartC];
2231: PetscScalar *pMat = NULL;
2233: PetscCall(PetscSectionGetDof(rootIndicesSec, p, &numColIndices));
2234: if (!numColIndices) continue;
2235: for (f = 0; f <= numFields; f++) {
2236: offsets[f] = 0;
2237: newOffsets[f] = 0;
2238: offsetsCopy[f] = 0;
2239: newOffsetsCopy[f] = 0;
2240: }
2241: numColIndices -= 2 * numFields;
2242: PetscCall(PetscSectionGetOffset(rootIndicesSec, p, &pIndOff));
2243: pInd = &rootIndices[pIndOff];
2244: PetscCall(PetscSectionGetDof(rootMatricesSec, p, &matSize));
2245: if (matSize) {
2246: PetscCall(PetscSectionGetOffset(rootMatricesSec, p, &pMatOff));
2247: pMat = &rootMatrices[pMatOff];
2248: }
2249: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
2250: if (dof < 0) dof = -(dof + 1);
2251: if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
2252: PetscInt i, j;
2254: if (matSize == 0) { /* don't need to calculate the mat, just the indices */
2255: PetscInt numIndices, *indices;
2256: PetscCall(DMPlexGetClosureIndices(coarse, localCoarse, globalCoarse, p, PETSC_TRUE, &numIndices, &indices, offsets, NULL));
2257: PetscCheck(numIndices == numColIndices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "mismatching constraint indices calculations");
2258: for (i = 0; i < numColIndices; i++) pInd[i] = indices[i];
2259: for (i = 0; i < numFields; i++) {
2260: pInd[numColIndices + i] = offsets[i + 1];
2261: pInd[numColIndices + numFields + i] = offsets[i + 1];
2262: }
2263: PetscCall(DMPlexRestoreClosureIndices(coarse, localCoarse, globalCoarse, p, PETSC_TRUE, &numIndices, &indices, offsets, NULL));
2264: } else {
2265: PetscInt closureSize, *closure = NULL, cl;
2266: PetscScalar *pMatIn, *pMatModified;
2267: PetscInt numPoints, *points;
2269: {
2270: PetscInt *closure = NULL, closureSize, cl;
2272: PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2273: for (cl = 0; cl < closureSize; cl++) { /* get the closure */
2274: PetscInt c = closure[2 * cl], clDof;
2276: PetscCall(PetscSectionGetDof(localCoarse, c, &clDof));
2277: numRowIndices += clDof;
2278: }
2279: PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2280: }
2282: PetscCall(DMGetWorkArray(coarse, numRowIndices * numRowIndices, MPIU_SCALAR, &pMatIn));
2283: for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */
2284: for (j = 0; j < numRowIndices; j++) pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.;
2285: }
2286: PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2287: for (f = 0; f < maxFields; f++) {
2288: if (numFields) PetscCall(PetscSectionGetFieldPointSyms(localCoarse, f, closureSize, closure, &perms[f], &flips[f]));
2289: else PetscCall(PetscSectionGetPointSyms(localCoarse, closureSize, closure, &perms[f], &flips[f]));
2290: }
2291: if (numFields) {
2292: for (cl = 0; cl < closureSize; cl++) {
2293: PetscInt c = closure[2 * cl];
2295: for (f = 0; f < numFields; f++) {
2296: PetscInt fDof;
2298: PetscCall(PetscSectionGetFieldDof(localCoarse, c, f, &fDof));
2299: offsets[f + 1] += fDof;
2300: }
2301: }
2302: for (f = 0; f < numFields; f++) {
2303: offsets[f + 1] += offsets[f];
2304: newOffsets[f + 1] = offsets[f + 1];
2305: }
2306: }
2307: /* TODO : flips here ? */
2308: /* apply hanging node constraints on the right, get the new points and the new offsets */
2309: PetscCall(DMPlexAnchorsModifyMat(coarse, localCoarse, closureSize, numRowIndices, closure, perms, pMatIn, &numPoints, NULL, &points, &pMatModified, newOffsets, PETSC_FALSE));
2310: for (f = 0; f < maxFields; f++) {
2311: if (numFields) PetscCall(PetscSectionRestoreFieldPointSyms(localCoarse, f, closureSize, closure, &perms[f], &flips[f]));
2312: else PetscCall(PetscSectionRestorePointSyms(localCoarse, closureSize, closure, &perms[f], &flips[f]));
2313: }
2314: for (f = 0; f < maxFields; f++) {
2315: if (numFields) PetscCall(PetscSectionGetFieldPointSyms(localCoarse, f, numPoints, points, &perms[f], &flips[f]));
2316: else PetscCall(PetscSectionGetPointSyms(localCoarse, numPoints, points, &perms[f], &flips[f]));
2317: }
2318: if (!numFields) {
2319: for (i = 0; i < numRowIndices * numColIndices; i++) pMat[i] = pMatModified[i];
2320: } else {
2321: PetscInt i, j, count;
2322: for (f = 0, count = 0; f < numFields; f++) {
2323: for (i = offsets[f]; i < offsets[f + 1]; i++) {
2324: for (j = newOffsets[f]; j < newOffsets[f + 1]; j++, count++) pMat[count] = pMatModified[i * numColIndices + j];
2325: }
2326: }
2327: }
2328: PetscCall(DMRestoreWorkArray(coarse, numRowIndices * numColIndices, MPIU_SCALAR, &pMatModified));
2329: PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
2330: PetscCall(DMRestoreWorkArray(coarse, numRowIndices * numColIndices, MPIU_SCALAR, &pMatIn));
2331: if (numFields) {
2332: for (f = 0; f < numFields; f++) {
2333: pInd[numColIndices + f] = offsets[f + 1];
2334: pInd[numColIndices + numFields + f] = newOffsets[f + 1];
2335: }
2336: for (cl = 0; cl < numPoints; cl++) {
2337: PetscInt globalOff, c = points[2 * cl];
2338: PetscCall(PetscSectionGetOffset(globalCoarse, c, &globalOff));
2339: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff + 1) : globalOff, newOffsets, PETSC_FALSE, perms, cl, NULL, pInd));
2340: }
2341: } else {
2342: for (cl = 0; cl < numPoints; cl++) {
2343: PetscInt c = points[2 * cl], globalOff;
2344: const PetscInt *perm = perms[0] ? perms[0][cl] : NULL;
2346: PetscCall(PetscSectionGetOffset(globalCoarse, c, &globalOff));
2347: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff + 1) : globalOff, newOffsets, PETSC_FALSE, perm, NULL, pInd));
2348: }
2349: }
2350: for (f = 0; f < maxFields; f++) {
2351: if (numFields) PetscCall(PetscSectionRestoreFieldPointSyms(localCoarse, f, numPoints, points, &perms[f], &flips[f]));
2352: else PetscCall(PetscSectionRestorePointSyms(localCoarse, numPoints, points, &perms[f], &flips[f]));
2353: }
2354: PetscCall(DMRestoreWorkArray(coarse, numPoints, MPIU_SCALAR, &points));
2355: }
2356: } else if (matSize) {
2357: PetscInt cOff;
2358: PetscInt *rowIndices, *colIndices, a, aDof = 0, aOff;
2360: numRowIndices = dof;
2361: PetscCall(DMGetWorkArray(coarse, numRowIndices, MPIU_INT, &rowIndices));
2362: PetscCall(DMGetWorkArray(coarse, numColIndices, MPIU_INT, &colIndices));
2363: PetscCall(PetscSectionGetOffset(cSec, p, &cOff));
2364: PetscCall(PetscSectionGetDof(aSec, p, &aDof));
2365: PetscCall(PetscSectionGetOffset(aSec, p, &aOff));
2366: if (numFields) {
2367: for (f = 0; f < numFields; f++) {
2368: PetscInt fDof;
2370: PetscCall(PetscSectionGetFieldDof(cSec, p, f, &fDof));
2371: offsets[f + 1] = fDof;
2372: for (a = 0; a < aDof; a++) {
2373: PetscInt anchor = anchors[a + aOff];
2374: PetscCall(PetscSectionGetFieldDof(localCoarse, anchor, f, &fDof));
2375: newOffsets[f + 1] += fDof;
2376: }
2377: }
2378: for (f = 0; f < numFields; f++) {
2379: offsets[f + 1] += offsets[f];
2380: offsetsCopy[f + 1] = offsets[f + 1];
2381: newOffsets[f + 1] += newOffsets[f];
2382: newOffsetsCopy[f + 1] = newOffsets[f + 1];
2383: }
2384: PetscCall(DMPlexGetIndicesPointFields_Internal(cSec, PETSC_TRUE, p, cOff, offsetsCopy, PETSC_TRUE, NULL, -1, NULL, rowIndices));
2385: for (a = 0; a < aDof; a++) {
2386: PetscInt anchor = anchors[a + aOff], lOff;
2387: PetscCall(PetscSectionGetOffset(localCoarse, anchor, &lOff));
2388: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_TRUE, anchor, lOff, newOffsetsCopy, PETSC_TRUE, NULL, -1, NULL, colIndices));
2389: }
2390: } else {
2391: PetscCall(DMPlexGetIndicesPoint_Internal(cSec, PETSC_TRUE, p, cOff, offsetsCopy, PETSC_TRUE, NULL, NULL, rowIndices));
2392: for (a = 0; a < aDof; a++) {
2393: PetscInt anchor = anchors[a + aOff], lOff;
2394: PetscCall(PetscSectionGetOffset(localCoarse, anchor, &lOff));
2395: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_TRUE, anchor, lOff, newOffsetsCopy, PETSC_TRUE, NULL, NULL, colIndices));
2396: }
2397: }
2398: if (numFields) {
2399: PetscInt count, a;
2401: for (f = 0, count = 0; f < numFields; f++) {
2402: PetscInt iSize = offsets[f + 1] - offsets[f];
2403: PetscInt jSize = newOffsets[f + 1] - newOffsets[f];
2404: PetscCall(MatGetValues(cMat, iSize, &rowIndices[offsets[f]], jSize, &colIndices[newOffsets[f]], &pMat[count]));
2405: count += iSize * jSize;
2406: pInd[numColIndices + f] = offsets[f + 1];
2407: pInd[numColIndices + numFields + f] = newOffsets[f + 1];
2408: }
2409: for (a = 0; a < aDof; a++) {
2410: PetscInt anchor = anchors[a + aOff];
2411: PetscInt gOff;
2412: PetscCall(PetscSectionGetOffset(globalCoarse, anchor, &gOff));
2413: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, anchor, gOff < 0 ? -(gOff + 1) : gOff, newOffsets, PETSC_FALSE, NULL, -1, NULL, pInd));
2414: }
2415: } else {
2416: PetscInt a;
2417: PetscCall(MatGetValues(cMat, numRowIndices, rowIndices, numColIndices, colIndices, pMat));
2418: for (a = 0; a < aDof; a++) {
2419: PetscInt anchor = anchors[a + aOff];
2420: PetscInt gOff;
2421: PetscCall(PetscSectionGetOffset(globalCoarse, anchor, &gOff));
2422: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, anchor, gOff < 0 ? -(gOff + 1) : gOff, newOffsets, PETSC_FALSE, NULL, NULL, pInd));
2423: }
2424: }
2425: PetscCall(DMRestoreWorkArray(coarse, numColIndices, MPIU_INT, &colIndices));
2426: PetscCall(DMRestoreWorkArray(coarse, numRowIndices, MPIU_INT, &rowIndices));
2427: } else {
2428: PetscInt gOff;
2430: PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));
2431: if (numFields) {
2432: for (f = 0; f < numFields; f++) {
2433: PetscInt fDof;
2434: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
2435: offsets[f + 1] = fDof + offsets[f];
2436: }
2437: for (f = 0; f < numFields; f++) {
2438: pInd[numColIndices + f] = offsets[f + 1];
2439: pInd[numColIndices + numFields + f] = offsets[f + 1];
2440: }
2441: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, -1, NULL, pInd));
2442: } else {
2443: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, NULL, pInd));
2444: }
2445: }
2446: }
2447: PetscCall(PetscFree(maxChildIds));
2448: }
2449: {
2450: PetscSF indicesSF, matricesSF;
2451: PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices;
2453: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafIndicesSec));
2454: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafMatricesSec));
2455: PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootIndicesSec, &remoteOffsetsIndices, leafIndicesSec));
2456: PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootMatricesSec, &remoteOffsetsMatrices, leafMatricesSec));
2457: PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootIndicesSec, remoteOffsetsIndices, leafIndicesSec, &indicesSF));
2458: PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootMatricesSec, remoteOffsetsMatrices, leafMatricesSec, &matricesSF));
2459: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
2460: PetscCall(PetscFree(remoteOffsetsIndices));
2461: PetscCall(PetscFree(remoteOffsetsMatrices));
2462: PetscCall(PetscSectionGetStorageSize(leafIndicesSec, &numLeafIndices));
2463: PetscCall(PetscSectionGetStorageSize(leafMatricesSec, &numLeafMatrices));
2464: PetscCall(PetscMalloc2(numLeafIndices, &leafIndices, numLeafMatrices, &leafMatrices));
2465: PetscCall(PetscSFBcastBegin(indicesSF, MPIU_INT, rootIndices, leafIndices, MPI_REPLACE));
2466: PetscCall(PetscSFBcastBegin(matricesSF, MPIU_SCALAR, rootMatrices, leafMatrices, MPI_REPLACE));
2467: PetscCall(PetscSFBcastEnd(indicesSF, MPIU_INT, rootIndices, leafIndices, MPI_REPLACE));
2468: PetscCall(PetscSFBcastEnd(matricesSF, MPIU_SCALAR, rootMatrices, leafMatrices, MPI_REPLACE));
2469: PetscCall(PetscSFDestroy(&matricesSF));
2470: PetscCall(PetscSFDestroy(&indicesSF));
2471: PetscCall(PetscFree2(rootIndices, rootMatrices));
2472: PetscCall(PetscSectionDestroy(&rootIndicesSec));
2473: PetscCall(PetscSectionDestroy(&rootMatricesSec));
2474: }
2475: /* count to preallocate */
2476: PetscCall(DMGetLocalSection(fine, &localFine));
2477: {
2478: PetscInt nGlobal;
2479: PetscInt *dnnz, *onnz;
2480: PetscLayout rowMap, colMap;
2481: PetscInt rowStart, rowEnd, colStart, colEnd;
2482: PetscInt maxDof;
2483: PetscInt *rowIndices;
2484: DM refTree;
2485: PetscInt **refPointFieldN;
2486: PetscScalar ***refPointFieldMats;
2487: PetscSection refConSec, refAnSec;
2488: PetscInt pRefStart, pRefEnd, maxConDof, maxColumns, leafStart, leafEnd;
2489: PetscScalar *pointWork;
2491: PetscCall(PetscSectionGetConstrainedStorageSize(globalFine, &nGlobal));
2492: PetscCall(PetscCalloc2(nGlobal, &dnnz, nGlobal, &onnz));
2493: PetscCall(MatGetLayouts(mat, &rowMap, &colMap));
2494: PetscCall(PetscLayoutSetUp(rowMap));
2495: PetscCall(PetscLayoutSetUp(colMap));
2496: PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd));
2497: PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd));
2498: PetscCall(PetscSectionGetMaxDof(localFine, &maxDof));
2499: PetscCall(PetscSectionGetChart(leafIndicesSec, &leafStart, &leafEnd));
2500: PetscCall(DMGetWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
2501: for (p = leafStart; p < leafEnd; p++) {
2502: PetscInt gDof, gcDof, gOff;
2503: PetscInt numColIndices, pIndOff, *pInd;
2504: PetscInt matSize;
2505: PetscInt i;
2507: PetscCall(PetscSectionGetDof(globalFine, p, &gDof));
2508: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &gcDof));
2509: if ((gDof - gcDof) <= 0) continue;
2510: PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
2511: PetscCheck(gOff >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "I though having global dofs meant a non-negative offset");
2512: PetscCheck(gOff >= rowStart && (gOff + gDof - gcDof) <= rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "I thought the row map would constrain the global dofs");
2513: PetscCall(PetscSectionGetDof(leafIndicesSec, p, &numColIndices));
2514: PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &pIndOff));
2515: numColIndices -= 2 * numFields;
2516: PetscCheck(numColIndices > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "global fine dof with no dofs to interpolate from");
2517: pInd = &leafIndices[pIndOff];
2518: offsets[0] = 0;
2519: offsetsCopy[0] = 0;
2520: newOffsets[0] = 0;
2521: newOffsetsCopy[0] = 0;
2522: if (numFields) {
2523: PetscInt f;
2524: for (f = 0; f < numFields; f++) {
2525: PetscInt rowDof;
2527: PetscCall(PetscSectionGetFieldDof(localFine, p, f, &rowDof));
2528: offsets[f + 1] = offsets[f] + rowDof;
2529: offsetsCopy[f + 1] = offsets[f + 1];
2530: newOffsets[f + 1] = pInd[numColIndices + numFields + f];
2531: numD[f] = 0;
2532: numO[f] = 0;
2533: }
2534: PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices));
2535: for (f = 0; f < numFields; f++) {
2536: PetscInt colOffset = newOffsets[f];
2537: PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f];
2539: for (i = 0; i < numFieldCols; i++) {
2540: PetscInt gInd = pInd[i + colOffset];
2542: if (gInd >= colStart && gInd < colEnd) {
2543: numD[f]++;
2544: } else if (gInd >= 0) { /* negative means non-entry */
2545: numO[f]++;
2546: }
2547: }
2548: }
2549: } else {
2550: PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices));
2551: numD[0] = 0;
2552: numO[0] = 0;
2553: for (i = 0; i < numColIndices; i++) {
2554: PetscInt gInd = pInd[i];
2556: if (gInd >= colStart && gInd < colEnd) {
2557: numD[0]++;
2558: } else if (gInd >= 0) { /* negative means non-entry */
2559: numO[0]++;
2560: }
2561: }
2562: }
2563: PetscCall(PetscSectionGetDof(leafMatricesSec, p, &matSize));
2564: if (!matSize) { /* incoming matrix is identity */
2565: PetscInt childId;
2567: childId = childIds[p - pStartF];
2568: if (childId < 0) { /* no child interpolation: one nnz per */
2569: if (numFields) {
2570: PetscInt f;
2571: for (f = 0; f < numFields; f++) {
2572: PetscInt numRows = offsets[f + 1] - offsets[f], row;
2573: for (row = 0; row < numRows; row++) {
2574: PetscInt gIndCoarse = pInd[newOffsets[f] + row];
2575: PetscInt gIndFine = rowIndices[offsets[f] + row];
2576: if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2577: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2578: dnnz[gIndFine - rowStart] = 1;
2579: } else if (gIndCoarse >= 0) { /* remote */
2580: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2581: onnz[gIndFine - rowStart] = 1;
2582: } else { /* constrained */
2583: PetscCheck(gIndFine < 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2584: }
2585: }
2586: }
2587: } else {
2588: PetscInt i;
2589: for (i = 0; i < gDof; i++) {
2590: PetscInt gIndCoarse = pInd[i];
2591: PetscInt gIndFine = rowIndices[i];
2592: if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2593: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2594: dnnz[gIndFine - rowStart] = 1;
2595: } else if (gIndCoarse >= 0) { /* remote */
2596: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2597: onnz[gIndFine - rowStart] = 1;
2598: } else { /* constrained */
2599: PetscCheck(gIndFine < 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2600: }
2601: }
2602: }
2603: } else { /* interpolate from all */
2604: if (numFields) {
2605: PetscInt f;
2606: for (f = 0; f < numFields; f++) {
2607: PetscInt numRows = offsets[f + 1] - offsets[f], row;
2608: for (row = 0; row < numRows; row++) {
2609: PetscInt gIndFine = rowIndices[offsets[f] + row];
2610: if (gIndFine >= 0) {
2611: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2612: dnnz[gIndFine - rowStart] = numD[f];
2613: onnz[gIndFine - rowStart] = numO[f];
2614: }
2615: }
2616: }
2617: } else {
2618: PetscInt i;
2619: for (i = 0; i < gDof; i++) {
2620: PetscInt gIndFine = rowIndices[i];
2621: if (gIndFine >= 0) {
2622: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2623: dnnz[gIndFine - rowStart] = numD[0];
2624: onnz[gIndFine - rowStart] = numO[0];
2625: }
2626: }
2627: }
2628: }
2629: } else { /* interpolate from all */
2630: if (numFields) {
2631: PetscInt f;
2632: for (f = 0; f < numFields; f++) {
2633: PetscInt numRows = offsets[f + 1] - offsets[f], row;
2634: for (row = 0; row < numRows; row++) {
2635: PetscInt gIndFine = rowIndices[offsets[f] + row];
2636: if (gIndFine >= 0) {
2637: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2638: dnnz[gIndFine - rowStart] = numD[f];
2639: onnz[gIndFine - rowStart] = numO[f];
2640: }
2641: }
2642: }
2643: } else { /* every dof get a full row */
2644: PetscInt i;
2645: for (i = 0; i < gDof; i++) {
2646: PetscInt gIndFine = rowIndices[i];
2647: if (gIndFine >= 0) {
2648: PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs");
2649: dnnz[gIndFine - rowStart] = numD[0];
2650: onnz[gIndFine - rowStart] = numO[0];
2651: }
2652: }
2653: }
2654: }
2655: }
2656: PetscCall(MatXAIJSetPreallocation(mat, 1, dnnz, onnz, NULL, NULL));
2657: PetscCall(PetscFree2(dnnz, onnz));
2659: PetscCall(DMPlexGetReferenceTree(fine, &refTree));
2660: PetscCall(DMCopyDisc(fine, refTree));
2661: PetscCall(DMSetLocalSection(refTree, NULL));
2662: PetscCall(DMSetDefaultConstraints(refTree, NULL, NULL, NULL));
2663: PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
2664: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
2665: PetscCall(DMPlexGetAnchors(refTree, &refAnSec, NULL));
2666: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
2667: PetscCall(PetscSectionGetMaxDof(refConSec, &maxConDof));
2668: PetscCall(PetscSectionGetMaxDof(leafIndicesSec, &maxColumns));
2669: PetscCall(PetscMalloc1(maxConDof * maxColumns, &pointWork));
2670: for (p = leafStart; p < leafEnd; p++) {
2671: PetscInt gDof, gcDof, gOff;
2672: PetscInt numColIndices, pIndOff, *pInd;
2673: PetscInt matSize;
2674: PetscInt childId;
2676: PetscCall(PetscSectionGetDof(globalFine, p, &gDof));
2677: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &gcDof));
2678: if ((gDof - gcDof) <= 0) continue;
2679: childId = childIds[p - pStartF];
2680: PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
2681: PetscCall(PetscSectionGetDof(leafIndicesSec, p, &numColIndices));
2682: PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &pIndOff));
2683: numColIndices -= 2 * numFields;
2684: pInd = &leafIndices[pIndOff];
2685: offsets[0] = 0;
2686: offsetsCopy[0] = 0;
2687: newOffsets[0] = 0;
2688: newOffsetsCopy[0] = 0;
2689: rowOffsets[0] = 0;
2690: if (numFields) {
2691: PetscInt f;
2692: for (f = 0; f < numFields; f++) {
2693: PetscInt rowDof;
2695: PetscCall(PetscSectionGetFieldDof(localFine, p, f, &rowDof));
2696: offsets[f + 1] = offsets[f] + rowDof;
2697: offsetsCopy[f + 1] = offsets[f + 1];
2698: rowOffsets[f + 1] = pInd[numColIndices + f];
2699: newOffsets[f + 1] = pInd[numColIndices + numFields + f];
2700: }
2701: PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices));
2702: } else {
2703: PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices));
2704: }
2705: PetscCall(PetscSectionGetDof(leafMatricesSec, p, &matSize));
2706: if (!matSize) { /* incoming matrix is identity */
2707: if (childId < 0) { /* no child interpolation: scatter */
2708: if (numFields) {
2709: PetscInt f;
2710: for (f = 0; f < numFields; f++) {
2711: PetscInt numRows = offsets[f + 1] - offsets[f], row;
2712: for (row = 0; row < numRows; row++) PetscCall(MatSetValue(mat, rowIndices[offsets[f] + row], pInd[newOffsets[f] + row], 1., INSERT_VALUES));
2713: }
2714: } else {
2715: PetscInt numRows = gDof, row;
2716: for (row = 0; row < numRows; row++) PetscCall(MatSetValue(mat, rowIndices[row], pInd[row], 1., INSERT_VALUES));
2717: }
2718: } else { /* interpolate from all */
2719: if (numFields) {
2720: PetscInt f;
2721: for (f = 0; f < numFields; f++) {
2722: PetscInt numRows = offsets[f + 1] - offsets[f];
2723: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2724: PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], refPointFieldMats[childId - pRefStart][f], INSERT_VALUES));
2725: }
2726: } else {
2727: PetscCall(MatSetValues(mat, gDof, rowIndices, numColIndices, pInd, refPointFieldMats[childId - pRefStart][0], INSERT_VALUES));
2728: }
2729: }
2730: } else { /* interpolate from all */
2731: PetscInt pMatOff;
2732: PetscScalar *pMat;
2734: PetscCall(PetscSectionGetOffset(leafMatricesSec, p, &pMatOff));
2735: pMat = &leafMatrices[pMatOff];
2736: if (childId < 0) { /* copy the incoming matrix */
2737: if (numFields) {
2738: PetscInt f, count;
2739: for (f = 0, count = 0; f < numFields; f++) {
2740: PetscInt numRows = offsets[f + 1] - offsets[f];
2741: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2742: PetscInt numInRows = rowOffsets[f + 1] - rowOffsets[f];
2743: PetscScalar *inMat = &pMat[count];
2745: PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], inMat, INSERT_VALUES));
2746: count += numCols * numInRows;
2747: }
2748: } else {
2749: PetscCall(MatSetValues(mat, gDof, rowIndices, numColIndices, pInd, pMat, INSERT_VALUES));
2750: }
2751: } else { /* multiply the incoming matrix by the child interpolation */
2752: if (numFields) {
2753: PetscInt f, count;
2754: for (f = 0, count = 0; f < numFields; f++) {
2755: PetscInt numRows = offsets[f + 1] - offsets[f];
2756: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2757: PetscInt numInRows = rowOffsets[f + 1] - rowOffsets[f];
2758: PetscScalar *inMat = &pMat[count];
2759: PetscInt i, j, k;
2760: PetscCheck(refPointFieldN[childId - pRefStart][f] == numInRows, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point constraint matrix multiply dimension mismatch");
2761: for (i = 0; i < numRows; i++) {
2762: for (j = 0; j < numCols; j++) {
2763: PetscScalar val = 0.;
2764: for (k = 0; k < numInRows; k++) val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j];
2765: pointWork[i * numCols + j] = val;
2766: }
2767: }
2768: PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], pointWork, INSERT_VALUES));
2769: count += numCols * numInRows;
2770: }
2771: } else { /* every dof gets a full row */
2772: PetscInt numRows = gDof;
2773: PetscInt numCols = numColIndices;
2774: PetscInt numInRows = matSize / numColIndices;
2775: PetscInt i, j, k;
2776: PetscCheck(refPointFieldN[childId - pRefStart][0] == numInRows, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point constraint matrix multiply dimension mismatch");
2777: for (i = 0; i < numRows; i++) {
2778: for (j = 0; j < numCols; j++) {
2779: PetscScalar val = 0.;
2780: for (k = 0; k < numInRows; k++) val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j];
2781: pointWork[i * numCols + j] = val;
2782: }
2783: }
2784: PetscCall(MatSetValues(mat, numRows, rowIndices, numCols, pInd, pointWork, INSERT_VALUES));
2785: }
2786: }
2787: }
2788: }
2789: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
2790: PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
2791: PetscCall(PetscFree(pointWork));
2792: }
2793: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
2794: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
2795: PetscCall(PetscSectionDestroy(&leafIndicesSec));
2796: PetscCall(PetscSectionDestroy(&leafMatricesSec));
2797: PetscCall(PetscFree2(leafIndices, leafMatrices));
2798: PetscCall(PetscFree2(*(PetscInt ****)&perms, *(PetscScalar ****)&flips));
2799: PetscCall(PetscFree7(offsets, offsetsCopy, newOffsets, newOffsetsCopy, rowOffsets, numD, numO));
2800: PetscCall(ISRestoreIndices(aIS, &anchors));
2801: PetscFunctionReturn(PETSC_SUCCESS);
2802: }
2804: /*
2805: * Assuming a nodal basis (w.r.t. the dual basis) basis:
2806: *
2807: * for each coarse dof \phi^c_i:
2808: * for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i:
2809: * for each fine dof \phi^f_j;
2810: * a_{i,j} = 0;
2811: * for each fine dof \phi^f_k:
2812: * a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l
2813: * [^^^ this is = \phi^c_i ^^^]
2814: */
2815: PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj)
2816: {
2817: PetscDS ds;
2818: PetscSection section, cSection;
2819: DMLabel canonical, depth;
2820: Mat cMat, mat;
2821: PetscInt *nnz;
2822: PetscInt f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd;
2823: PetscInt m, n;
2824: PetscScalar *pointScalar;
2825: PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent;
2827: PetscFunctionBegin;
2828: PetscCall(DMGetLocalSection(refTree, §ion));
2829: PetscCall(DMGetDimension(refTree, &dim));
2830: PetscCall(PetscMalloc6(dim, &v0, dim, &v0parent, dim, &vtmp, dim * dim, &J, dim * dim, &Jparent, dim * dim, &invJ));
2831: PetscCall(PetscMalloc2(dim, &pointScalar, dim, &pointRef));
2832: PetscCall(DMGetDS(refTree, &ds));
2833: PetscCall(PetscDSGetNumFields(ds, &numFields));
2834: PetscCall(PetscSectionGetNumFields(section, &numSecFields));
2835: PetscCall(DMGetLabel(refTree, "canonical", &canonical));
2836: PetscCall(DMGetLabel(refTree, "depth", &depth));
2837: PetscCall(DMGetDefaultConstraints(refTree, &cSection, &cMat, NULL));
2838: PetscCall(DMPlexGetChart(refTree, &pStart, &pEnd));
2839: PetscCall(DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd));
2840: PetscCall(MatGetSize(cMat, &n, &m)); /* the injector has transpose sizes from the constraint matrix */
2841: /* Step 1: compute non-zero pattern. A proper subset of constraint matrix non-zero */
2842: PetscCall(PetscCalloc1(m, &nnz));
2843: for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */
2844: const PetscInt *children;
2845: PetscInt numChildren;
2846: PetscInt i, numChildDof, numSelfDof;
2848: if (canonical) {
2849: PetscInt pCanonical;
2850: PetscCall(DMLabelGetValue(canonical, p, &pCanonical));
2851: if (p != pCanonical) continue;
2852: }
2853: PetscCall(DMPlexGetTreeChildren(refTree, p, &numChildren, &children));
2854: if (!numChildren) continue;
2855: for (i = 0, numChildDof = 0; i < numChildren; i++) {
2856: PetscInt child = children[i];
2857: PetscInt dof;
2859: PetscCall(PetscSectionGetDof(section, child, &dof));
2860: numChildDof += dof;
2861: }
2862: PetscCall(PetscSectionGetDof(section, p, &numSelfDof));
2863: if (!numChildDof || !numSelfDof) continue;
2864: for (f = 0; f < numFields; f++) {
2865: PetscInt selfOff;
2867: if (numSecFields) { /* count the dofs for just this field */
2868: for (i = 0, numChildDof = 0; i < numChildren; i++) {
2869: PetscInt child = children[i];
2870: PetscInt dof;
2872: PetscCall(PetscSectionGetFieldDof(section, child, f, &dof));
2873: numChildDof += dof;
2874: }
2875: PetscCall(PetscSectionGetFieldDof(section, p, f, &numSelfDof));
2876: PetscCall(PetscSectionGetFieldOffset(section, p, f, &selfOff));
2877: } else {
2878: PetscCall(PetscSectionGetOffset(section, p, &selfOff));
2879: }
2880: for (i = 0; i < numSelfDof; i++) nnz[selfOff + i] = numChildDof;
2881: }
2882: }
2883: PetscCall(MatCreateAIJ(PETSC_COMM_SELF, m, n, m, n, -1, nnz, -1, NULL, &mat));
2884: PetscCall(PetscFree(nnz));
2885: /* Setp 2: compute entries */
2886: for (p = pStart; p < pEnd; p++) {
2887: const PetscInt *children;
2888: PetscInt numChildren;
2889: PetscInt i, numChildDof, numSelfDof;
2891: /* same conditions about when entries occur */
2892: if (canonical) {
2893: PetscInt pCanonical;
2894: PetscCall(DMLabelGetValue(canonical, p, &pCanonical));
2895: if (p != pCanonical) continue;
2896: }
2897: PetscCall(DMPlexGetTreeChildren(refTree, p, &numChildren, &children));
2898: if (!numChildren) continue;
2899: for (i = 0, numChildDof = 0; i < numChildren; i++) {
2900: PetscInt child = children[i];
2901: PetscInt dof;
2903: PetscCall(PetscSectionGetDof(section, child, &dof));
2904: numChildDof += dof;
2905: }
2906: PetscCall(PetscSectionGetDof(section, p, &numSelfDof));
2907: if (!numChildDof || !numSelfDof) continue;
2909: for (f = 0; f < numFields; f++) {
2910: PetscInt pI = -1, cI = -1;
2911: PetscInt selfOff, Nc, parentCell;
2912: PetscInt cellShapeOff;
2913: PetscObject disc;
2914: PetscDualSpace dsp;
2915: PetscClassId classId;
2916: PetscScalar *pointMat;
2917: PetscInt *matRows, *matCols;
2918: PetscInt pO = PETSC_INT_MIN;
2919: const PetscInt *depthNumDof;
2921: if (numSecFields) {
2922: for (i = 0, numChildDof = 0; i < numChildren; i++) {
2923: PetscInt child = children[i];
2924: PetscInt dof;
2926: PetscCall(PetscSectionGetFieldDof(section, child, f, &dof));
2927: numChildDof += dof;
2928: }
2929: PetscCall(PetscSectionGetFieldDof(section, p, f, &numSelfDof));
2930: PetscCall(PetscSectionGetFieldOffset(section, p, f, &selfOff));
2931: } else {
2932: PetscCall(PetscSectionGetOffset(section, p, &selfOff));
2933: }
2935: /* find a cell whose closure contains p */
2936: if (p >= cStart && p < cEnd) {
2937: parentCell = p;
2938: } else {
2939: PetscInt *star = NULL;
2940: PetscInt numStar;
2942: parentCell = -1;
2943: PetscCall(DMPlexGetTransitiveClosure(refTree, p, PETSC_FALSE, &numStar, &star));
2944: for (i = numStar - 1; i >= 0; i--) {
2945: PetscInt c = star[2 * i];
2947: if (c >= cStart && c < cEnd) {
2948: parentCell = c;
2949: break;
2950: }
2951: }
2952: PetscCall(DMPlexRestoreTransitiveClosure(refTree, p, PETSC_FALSE, &numStar, &star));
2953: }
2954: /* determine the offset of p's shape functions within parentCell's shape functions */
2955: PetscCall(PetscDSGetDiscretization(ds, f, &disc));
2956: PetscCall(PetscObjectGetClassId(disc, &classId));
2957: if (classId == PETSCFE_CLASSID) PetscCall(PetscFEGetDualSpace((PetscFE)disc, &dsp));
2958: else {
2959: PetscCheck(classId == PETSCFV_CLASSID, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported discretization object");
2960: PetscCall(PetscFVGetDualSpace((PetscFV)disc, &dsp));
2961: }
2962: PetscCall(PetscDualSpaceGetNumDof(dsp, &depthNumDof));
2963: PetscCall(PetscDualSpaceGetNumComponents(dsp, &Nc));
2964: {
2965: PetscInt *closure = NULL;
2966: PetscInt numClosure;
2968: PetscCall(DMPlexGetTransitiveClosure(refTree, parentCell, PETSC_TRUE, &numClosure, &closure));
2969: for (i = 0, pI = -1, cellShapeOff = 0; i < numClosure; i++) {
2970: PetscInt point = closure[2 * i], pointDepth;
2972: pO = closure[2 * i + 1];
2973: if (point == p) {
2974: pI = i;
2975: break;
2976: }
2977: PetscCall(DMLabelGetValue(depth, point, &pointDepth));
2978: cellShapeOff += depthNumDof[pointDepth];
2979: }
2980: PetscCall(DMPlexRestoreTransitiveClosure(refTree, parentCell, PETSC_TRUE, &numClosure, &closure));
2981: }
2983: PetscCall(DMGetWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR, &pointMat));
2984: PetscCall(DMGetWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT, &matRows));
2985: matCols = matRows + numSelfDof;
2986: for (i = 0; i < numSelfDof; i++) matRows[i] = selfOff + i;
2987: for (i = 0; i < numSelfDof * numChildDof; i++) pointMat[i] = 0.;
2988: {
2989: PetscInt colOff = 0;
2991: for (i = 0; i < numChildren; i++) {
2992: PetscInt child = children[i];
2993: PetscInt dof, off, j;
2995: if (numSecFields) {
2996: PetscCall(PetscSectionGetFieldDof(cSection, child, f, &dof));
2997: PetscCall(PetscSectionGetFieldOffset(cSection, child, f, &off));
2998: } else {
2999: PetscCall(PetscSectionGetDof(cSection, child, &dof));
3000: PetscCall(PetscSectionGetOffset(cSection, child, &off));
3001: }
3003: for (j = 0; j < dof; j++) matCols[colOff++] = off + j;
3004: }
3005: }
3006: if (classId == PETSCFE_CLASSID) {
3007: PetscFE fe = (PetscFE)disc;
3008: PetscInt fSize;
3009: const PetscInt ***perms;
3010: const PetscScalar ***flips;
3011: const PetscInt *pperms;
3013: PetscCall(PetscFEGetDualSpace(fe, &dsp));
3014: PetscCall(PetscDualSpaceGetDimension(dsp, &fSize));
3015: PetscCall(PetscDualSpaceGetSymmetries(dsp, &perms, &flips));
3016: pperms = perms ? perms[pI] ? perms[pI][pO] : NULL : NULL;
3017: for (i = 0; i < numSelfDof; i++) { /* for every shape function */
3018: PetscQuadrature q;
3019: PetscInt dim, thisNc, numPoints, j, k;
3020: const PetscReal *points;
3021: const PetscReal *weights;
3022: PetscInt *closure = NULL;
3023: PetscInt numClosure;
3024: PetscInt iCell = pperms ? pperms[i] : i;
3025: PetscInt parentCellShapeDof = cellShapeOff + iCell;
3026: PetscTabulation Tparent;
3028: PetscCall(PetscDualSpaceGetFunctional(dsp, parentCellShapeDof, &q));
3029: PetscCall(PetscQuadratureGetData(q, &dim, &thisNc, &numPoints, &points, &weights));
3030: PetscCheck(thisNc == Nc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Functional dim %" PetscInt_FMT " does not much basis dim %" PetscInt_FMT, thisNc, Nc);
3031: PetscCall(PetscFECreateTabulation(fe, 1, numPoints, points, 0, &Tparent)); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */
3032: for (j = 0; j < numPoints; j++) {
3033: PetscInt childCell = -1;
3034: PetscReal *parentValAtPoint;
3035: const PetscReal xi0[3] = {-1., -1., -1.};
3036: const PetscReal *pointReal = &points[dim * j];
3037: const PetscScalar *point;
3038: PetscTabulation Tchild;
3039: PetscInt childCellShapeOff, pointMatOff;
3040: #if defined(PETSC_USE_COMPLEX)
3041: PetscInt d;
3043: for (d = 0; d < dim; d++) pointScalar[d] = points[dim * j + d];
3044: point = pointScalar;
3045: #else
3046: point = pointReal;
3047: #endif
3049: parentValAtPoint = &Tparent->T[0][(fSize * j + parentCellShapeDof) * Nc];
3051: for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/
3052: PetscInt child = children[k];
3053: PetscInt *star = NULL;
3054: PetscInt numStar, s;
3056: PetscCall(DMPlexGetTransitiveClosure(refTree, child, PETSC_FALSE, &numStar, &star));
3057: for (s = numStar - 1; s >= 0; s--) {
3058: PetscInt c = star[2 * s];
3060: if (c < cStart || c >= cEnd) continue;
3061: PetscCall(DMPlexLocatePoint_Internal(refTree, dim, point, c, &childCell));
3062: if (childCell >= 0) break;
3063: }
3064: PetscCall(DMPlexRestoreTransitiveClosure(refTree, child, PETSC_FALSE, &numStar, &star));
3065: if (childCell >= 0) break;
3066: }
3067: PetscCheck(childCell >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not locate quadrature point");
3068: PetscCall(DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ));
3069: PetscCall(DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent));
3070: CoordinatesRefToReal(dim, dim, xi0, v0parent, Jparent, pointReal, vtmp);
3071: CoordinatesRealToRef(dim, dim, xi0, v0, invJ, vtmp, pointRef);
3073: PetscCall(PetscFECreateTabulation(fe, 1, 1, pointRef, 0, &Tchild));
3074: PetscCall(DMPlexGetTransitiveClosure(refTree, childCell, PETSC_TRUE, &numClosure, &closure));
3075: for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */
3076: PetscInt child = children[k], childDepth, childDof, childO = PETSC_INT_MIN;
3077: PetscInt l;
3078: const PetscInt *cperms;
3080: PetscCall(DMLabelGetValue(depth, child, &childDepth));
3081: childDof = depthNumDof[childDepth];
3082: for (l = 0, cI = -1, childCellShapeOff = 0; l < numClosure; l++) {
3083: PetscInt point = closure[2 * l];
3084: PetscInt pointDepth;
3086: childO = closure[2 * l + 1];
3087: if (point == child) {
3088: cI = l;
3089: break;
3090: }
3091: PetscCall(DMLabelGetValue(depth, point, &pointDepth));
3092: childCellShapeOff += depthNumDof[pointDepth];
3093: }
3094: if (l == numClosure) {
3095: pointMatOff += childDof;
3096: continue; /* child is not in the closure of the cell: has nothing to contribute to this point */
3097: }
3098: cperms = perms ? perms[cI] ? perms[cI][childO] : NULL : NULL;
3099: for (l = 0; l < childDof; l++) {
3100: PetscInt lCell = cperms ? cperms[l] : l;
3101: PetscInt childCellDof = childCellShapeOff + lCell;
3102: PetscReal *childValAtPoint;
3103: PetscReal val = 0.;
3105: childValAtPoint = &Tchild->T[0][childCellDof * Nc];
3106: for (m = 0; m < Nc; m++) val += weights[j * Nc + m] * parentValAtPoint[m] * childValAtPoint[m];
3108: pointMat[i * numChildDof + pointMatOff + l] += val;
3109: }
3110: pointMatOff += childDof;
3111: }
3112: PetscCall(DMPlexRestoreTransitiveClosure(refTree, childCell, PETSC_TRUE, &numClosure, &closure));
3113: PetscCall(PetscTabulationDestroy(&Tchild));
3114: }
3115: PetscCall(PetscTabulationDestroy(&Tparent));
3116: }
3117: } else { /* just the volume-weighted averages of the children */
3118: PetscReal parentVol;
3119: PetscInt childCell;
3121: PetscCall(DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL));
3122: for (i = 0, childCell = 0; i < numChildren; i++) {
3123: PetscInt child = children[i], j;
3124: PetscReal childVol;
3126: if (child < cStart || child >= cEnd) continue;
3127: PetscCall(DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL));
3128: for (j = 0; j < Nc; j++) pointMat[j * numChildDof + Nc * childCell + j] = childVol / parentVol;
3129: childCell++;
3130: }
3131: }
3132: /* Insert pointMat into mat */
3133: PetscCall(MatSetValues(mat, numSelfDof, matRows, numChildDof, matCols, pointMat, INSERT_VALUES));
3134: PetscCall(DMRestoreWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT, &matRows));
3135: PetscCall(DMRestoreWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR, &pointMat));
3136: }
3137: }
3138: PetscCall(PetscFree6(v0, v0parent, vtmp, J, Jparent, invJ));
3139: PetscCall(PetscFree2(pointScalar, pointRef));
3140: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
3141: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
3142: *inj = mat;
3143: PetscFunctionReturn(PETSC_SUCCESS);
3144: }
3146: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3147: {
3148: PetscDS ds;
3149: PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof;
3150: PetscScalar ***refPointFieldMats;
3151: PetscSection refConSec, refSection;
3153: PetscFunctionBegin;
3154: PetscCall(DMGetDS(refTree, &ds));
3155: PetscCall(PetscDSGetNumFields(ds, &numFields));
3156: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
3157: PetscCall(DMGetLocalSection(refTree, &refSection));
3158: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
3159: PetscCall(PetscMalloc1(pRefEnd - pRefStart, &refPointFieldMats));
3160: PetscCall(PetscSectionGetMaxDof(refConSec, &maxDof));
3161: PetscCall(PetscMalloc1(maxDof, &rows));
3162: PetscCall(PetscMalloc1(maxDof * maxDof, &cols));
3163: for (p = pRefStart; p < pRefEnd; p++) {
3164: PetscInt parent, pDof, parentDof;
3166: PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL));
3167: PetscCall(PetscSectionGetDof(refConSec, p, &pDof));
3168: PetscCall(PetscSectionGetDof(refSection, parent, &parentDof));
3169: if (!pDof || !parentDof || parent == p) continue;
3171: PetscCall(PetscMalloc1(numFields, &refPointFieldMats[p - pRefStart]));
3172: for (f = 0; f < numFields; f++) {
3173: PetscInt cDof, cOff, numCols, r;
3175: if (numFields > 1) {
3176: PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof));
3177: PetscCall(PetscSectionGetFieldOffset(refConSec, p, f, &cOff));
3178: } else {
3179: PetscCall(PetscSectionGetDof(refConSec, p, &cDof));
3180: PetscCall(PetscSectionGetOffset(refConSec, p, &cOff));
3181: }
3183: for (r = 0; r < cDof; r++) rows[r] = cOff + r;
3184: numCols = 0;
3185: {
3186: PetscInt aDof, aOff, j;
3188: if (numFields > 1) {
3189: PetscCall(PetscSectionGetFieldDof(refSection, parent, f, &aDof));
3190: PetscCall(PetscSectionGetFieldOffset(refSection, parent, f, &aOff));
3191: } else {
3192: PetscCall(PetscSectionGetDof(refSection, parent, &aDof));
3193: PetscCall(PetscSectionGetOffset(refSection, parent, &aOff));
3194: }
3196: for (j = 0; j < aDof; j++) cols[numCols++] = aOff + j;
3197: }
3198: PetscCall(PetscMalloc1(cDof * numCols, &refPointFieldMats[p - pRefStart][f]));
3199: /* transpose of constraint matrix */
3200: PetscCall(MatGetValues(inj, numCols, cols, cDof, rows, refPointFieldMats[p - pRefStart][f]));
3201: }
3202: }
3203: *childrenMats = refPointFieldMats;
3204: PetscCall(PetscFree(rows));
3205: PetscCall(PetscFree(cols));
3206: PetscFunctionReturn(PETSC_SUCCESS);
3207: }
3209: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3210: {
3211: PetscDS ds;
3212: PetscScalar ***refPointFieldMats;
3213: PetscInt numFields, pRefStart, pRefEnd, p, f;
3214: PetscSection refConSec, refSection;
3216: PetscFunctionBegin;
3217: refPointFieldMats = *childrenMats;
3218: *childrenMats = NULL;
3219: PetscCall(DMGetDS(refTree, &ds));
3220: PetscCall(DMGetLocalSection(refTree, &refSection));
3221: PetscCall(PetscDSGetNumFields(ds, &numFields));
3222: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
3223: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
3224: for (p = pRefStart; p < pRefEnd; p++) {
3225: PetscInt parent, pDof, parentDof;
3227: PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL));
3228: PetscCall(PetscSectionGetDof(refConSec, p, &pDof));
3229: PetscCall(PetscSectionGetDof(refSection, parent, &parentDof));
3230: if (!pDof || !parentDof || parent == p) continue;
3232: for (f = 0; f < numFields; f++) {
3233: PetscInt cDof;
3235: if (numFields > 1) {
3236: PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof));
3237: } else {
3238: PetscCall(PetscSectionGetDof(refConSec, p, &cDof));
3239: }
3241: PetscCall(PetscFree(refPointFieldMats[p - pRefStart][f]));
3242: }
3243: PetscCall(PetscFree(refPointFieldMats[p - pRefStart]));
3244: }
3245: PetscCall(PetscFree(refPointFieldMats));
3246: PetscFunctionReturn(PETSC_SUCCESS);
3247: }
3249: static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree, Mat *injRef)
3250: {
3251: Mat cMatRef;
3252: PetscObject injRefObj;
3254: PetscFunctionBegin;
3255: PetscCall(DMGetDefaultConstraints(refTree, NULL, &cMatRef, NULL));
3256: PetscCall(PetscObjectQuery((PetscObject)cMatRef, "DMPlexComputeInjectorTree_refTree", &injRefObj));
3257: *injRef = (Mat)injRefObj;
3258: if (!*injRef) {
3259: PetscCall(DMPlexComputeInjectorReferenceTree(refTree, injRef));
3260: PetscCall(PetscObjectCompose((PetscObject)cMatRef, "DMPlexComputeInjectorTree_refTree", (PetscObject)*injRef));
3261: /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */
3262: PetscCall(PetscObjectDereference((PetscObject)*injRef));
3263: }
3264: PetscFunctionReturn(PETSC_SUCCESS);
3265: }
3267: 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)
3268: {
3269: PetscInt pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti;
3270: PetscSection globalCoarse, globalFine;
3271: PetscSection localCoarse, localFine, leafIndicesSec;
3272: PetscSection multiRootSec, rootIndicesSec;
3273: PetscInt *leafInds, *rootInds = NULL;
3274: const PetscInt *rootDegrees;
3275: PetscScalar *leafVals = NULL, *rootVals = NULL;
3276: PetscSF coarseToFineEmbedded;
3278: PetscFunctionBegin;
3279: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
3280: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
3281: PetscCall(DMGetLocalSection(fine, &localFine));
3282: PetscCall(DMGetGlobalSection(fine, &globalFine));
3283: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafIndicesSec));
3284: PetscCall(PetscSectionSetChart(leafIndicesSec, pStartF, pEndF));
3285: PetscCall(PetscSectionGetMaxDof(localFine, &maxDof));
3286: { /* winnow fine points that don't have global dofs out of the sf */
3287: PetscInt l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices;
3288: const PetscInt *leaves;
3290: PetscCall(PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL));
3291: for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
3292: p = leaves ? leaves[l] : l;
3293: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
3294: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
3295: if ((dof - cdof) > 0) {
3296: numPointsWithDofs++;
3298: PetscCall(PetscSectionGetDof(localFine, p, &dof));
3299: PetscCall(PetscSectionSetDof(leafIndicesSec, p, dof + 1));
3300: }
3301: }
3302: PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
3303: PetscCall(PetscSectionSetUp(leafIndicesSec));
3304: PetscCall(PetscSectionGetStorageSize(leafIndicesSec, &numIndices));
3305: PetscCall(PetscMalloc1((gatheredIndices ? numIndices : (maxDof + 1)), &leafInds));
3306: if (gatheredValues) PetscCall(PetscMalloc1(numIndices, &leafVals));
3307: for (l = 0, offset = 0; l < nleaves; l++) {
3308: p = leaves ? leaves[l] : l;
3309: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
3310: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
3311: if ((dof - cdof) > 0) {
3312: PetscInt off, gOff;
3313: PetscInt *pInd;
3314: PetscScalar *pVal = NULL;
3316: pointsWithDofs[offset++] = l;
3318: PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &off));
3320: pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds;
3321: if (gatheredValues) {
3322: PetscInt i;
3324: pVal = &leafVals[off + 1];
3325: for (i = 0; i < dof; i++) pVal[i] = 0.;
3326: }
3327: PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
3329: offsets[0] = 0;
3330: if (numFields) {
3331: PetscInt f;
3333: for (f = 0; f < numFields; f++) {
3334: PetscInt fDof;
3335: PetscCall(PetscSectionGetFieldDof(localFine, p, f, &fDof));
3336: offsets[f + 1] = fDof + offsets[f];
3337: }
3338: PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, -1, NULL, pInd));
3339: } else {
3340: PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, NULL, pInd));
3341: }
3342: if (gatheredValues) PetscCall(VecGetValues(fineVec, dof, pInd, pVal));
3343: }
3344: }
3345: PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded));
3346: PetscCall(PetscFree(pointsWithDofs));
3347: }
3349: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
3350: PetscCall(DMGetLocalSection(coarse, &localCoarse));
3351: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
3353: { /* there may be the case where an sf root has a parent: broadcast parents back to children */
3354: MPI_Datatype threeInt;
3355: PetscMPIInt rank;
3356: PetscInt (*parentNodeAndIdCoarse)[3];
3357: PetscInt (*parentNodeAndIdFine)[3];
3358: PetscInt p, nleaves, nleavesToParents;
3359: PetscSF pointSF, sfToParents;
3360: const PetscInt *ilocal;
3361: const PetscSFNode *iremote;
3362: PetscSFNode *iremoteToParents;
3363: PetscInt *ilocalToParents;
3365: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank));
3366: PetscCallMPI(MPI_Type_contiguous(3, MPIU_INT, &threeInt));
3367: PetscCallMPI(MPI_Type_commit(&threeInt));
3368: PetscCall(PetscMalloc2(pEndC - pStartC, &parentNodeAndIdCoarse, pEndF - pStartF, &parentNodeAndIdFine));
3369: PetscCall(DMGetPointSF(coarse, &pointSF));
3370: PetscCall(PetscSFGetGraph(pointSF, NULL, &nleaves, &ilocal, &iremote));
3371: for (p = pStartC; p < pEndC; p++) {
3372: PetscInt parent, childId;
3373: PetscCall(DMPlexGetTreeParent(coarse, p, &parent, &childId));
3374: parentNodeAndIdCoarse[p - pStartC][0] = rank;
3375: parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC;
3376: parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId;
3377: if (nleaves > 0) {
3378: PetscInt leaf = -1;
3380: if (ilocal) {
3381: PetscCall(PetscFindInt(parent, nleaves, ilocal, &leaf));
3382: } else {
3383: leaf = p - pStartC;
3384: }
3385: if (leaf >= 0) {
3386: parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank;
3387: parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index;
3388: }
3389: }
3390: }
3391: for (p = pStartF; p < pEndF; p++) {
3392: parentNodeAndIdFine[p - pStartF][0] = -1;
3393: parentNodeAndIdFine[p - pStartF][1] = -1;
3394: parentNodeAndIdFine[p - pStartF][2] = -1;
3395: }
3396: PetscCall(PetscSFBcastBegin(coarseToFineEmbedded, threeInt, parentNodeAndIdCoarse, parentNodeAndIdFine, MPI_REPLACE));
3397: PetscCall(PetscSFBcastEnd(coarseToFineEmbedded, threeInt, parentNodeAndIdCoarse, parentNodeAndIdFine, MPI_REPLACE));
3398: for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3399: PetscInt dof;
3401: PetscCall(PetscSectionGetDof(leafIndicesSec, p, &dof));
3402: if (dof) {
3403: PetscInt off;
3405: PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &off));
3406: if (gatheredIndices) {
3407: leafInds[off] = PetscMax(childIds[p - pStartF], parentNodeAndIdFine[p - pStartF][2]);
3408: } else if (gatheredValues) {
3409: leafVals[off] = (PetscScalar)PetscMax(childIds[p - pStartF], parentNodeAndIdFine[p - pStartF][2]);
3410: }
3411: }
3412: if (parentNodeAndIdFine[p - pStartF][0] >= 0) nleavesToParents++;
3413: }
3414: PetscCall(PetscMalloc1(nleavesToParents, &ilocalToParents));
3415: PetscCall(PetscMalloc1(nleavesToParents, &iremoteToParents));
3416: for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3417: if (parentNodeAndIdFine[p - pStartF][0] >= 0) {
3418: ilocalToParents[nleavesToParents] = p - pStartF;
3419: // FIXME PetscCall(PetscMPIIntCast(parentNodeAndIdFine[p - pStartF][0],&iremoteToParents[nleavesToParents].rank));
3420: iremoteToParents[nleavesToParents].rank = parentNodeAndIdFine[p - pStartF][0];
3421: iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p - pStartF][1];
3422: nleavesToParents++;
3423: }
3424: }
3425: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)coarse), &sfToParents));
3426: PetscCall(PetscSFSetGraph(sfToParents, pEndC - pStartC, nleavesToParents, ilocalToParents, PETSC_OWN_POINTER, iremoteToParents, PETSC_OWN_POINTER));
3427: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
3429: coarseToFineEmbedded = sfToParents;
3431: PetscCall(PetscFree2(parentNodeAndIdCoarse, parentNodeAndIdFine));
3432: PetscCallMPI(MPI_Type_free(&threeInt));
3433: }
3435: { /* winnow out coarse points that don't have dofs */
3436: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;
3437: PetscSF sfDofsOnly;
3439: for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) {
3440: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3441: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
3442: if ((dof - cdof) > 0) numPointsWithDofs++;
3443: }
3444: PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
3445: for (p = pStartC, offset = 0; p < pEndC; p++) {
3446: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3447: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
3448: if ((dof - cdof) > 0) pointsWithDofs[offset++] = p - pStartC;
3449: }
3450: PetscCall(PetscSFCreateEmbeddedRootSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly));
3451: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
3452: PetscCall(PetscFree(pointsWithDofs));
3453: coarseToFineEmbedded = sfDofsOnly;
3454: }
3456: /* communicate back to the coarse mesh which coarse points have children (that may require injection) */
3457: PetscCall(PetscSFComputeDegreeBegin(coarseToFineEmbedded, &rootDegrees));
3458: PetscCall(PetscSFComputeDegreeEnd(coarseToFineEmbedded, &rootDegrees));
3459: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &multiRootSec));
3460: PetscCall(PetscSectionSetChart(multiRootSec, pStartC, pEndC));
3461: for (p = pStartC; p < pEndC; p++) PetscCall(PetscSectionSetDof(multiRootSec, p, rootDegrees[p - pStartC]));
3462: PetscCall(PetscSectionSetUp(multiRootSec));
3463: PetscCall(PetscSectionGetStorageSize(multiRootSec, &numMulti));
3464: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootIndicesSec));
3465: { /* distribute the leaf section */
3466: PetscSF multi, multiInv, indicesSF;
3467: PetscInt *remoteOffsets, numRootIndices;
3469: PetscCall(PetscSFGetMultiSF(coarseToFineEmbedded, &multi));
3470: PetscCall(PetscSFCreateInverseSF(multi, &multiInv));
3471: PetscCall(PetscSFDistributeSection(multiInv, leafIndicesSec, &remoteOffsets, rootIndicesSec));
3472: PetscCall(PetscSFCreateSectionSF(multiInv, leafIndicesSec, remoteOffsets, rootIndicesSec, &indicesSF));
3473: PetscCall(PetscFree(remoteOffsets));
3474: PetscCall(PetscSFDestroy(&multiInv));
3475: PetscCall(PetscSectionGetStorageSize(rootIndicesSec, &numRootIndices));
3476: if (gatheredIndices) {
3477: PetscCall(PetscMalloc1(numRootIndices, &rootInds));
3478: PetscCall(PetscSFBcastBegin(indicesSF, MPIU_INT, leafInds, rootInds, MPI_REPLACE));
3479: PetscCall(PetscSFBcastEnd(indicesSF, MPIU_INT, leafInds, rootInds, MPI_REPLACE));
3480: }
3481: if (gatheredValues) {
3482: PetscCall(PetscMalloc1(numRootIndices, &rootVals));
3483: PetscCall(PetscSFBcastBegin(indicesSF, MPIU_SCALAR, leafVals, rootVals, MPI_REPLACE));
3484: PetscCall(PetscSFBcastEnd(indicesSF, MPIU_SCALAR, leafVals, rootVals, MPI_REPLACE));
3485: }
3486: PetscCall(PetscSFDestroy(&indicesSF));
3487: }
3488: PetscCall(PetscSectionDestroy(&leafIndicesSec));
3489: PetscCall(PetscFree(leafInds));
3490: PetscCall(PetscFree(leafVals));
3491: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
3492: *rootMultiSec = multiRootSec;
3493: *multiLeafSec = rootIndicesSec;
3494: if (gatheredIndices) *gatheredIndices = rootInds;
3495: if (gatheredValues) *gatheredValues = rootVals;
3496: PetscFunctionReturn(PETSC_SUCCESS);
3497: }
3499: PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
3500: {
3501: DM refTree;
3502: PetscSection multiRootSec, rootIndicesSec;
3503: PetscSection globalCoarse, globalFine;
3504: PetscSection localCoarse, localFine;
3505: PetscSection cSecRef;
3506: PetscInt *rootIndices = NULL, *parentIndices, pRefStart, pRefEnd;
3507: Mat injRef;
3508: PetscInt numFields, maxDof;
3509: PetscInt pStartC, pEndC, pStartF, pEndF, p;
3510: PetscInt *offsets, *offsetsCopy, *rowOffsets;
3511: PetscLayout rowMap, colMap;
3512: PetscInt rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO;
3513: PetscScalar ***childrenMats = NULL; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */
3515: PetscFunctionBegin;
3516: /* get the templates for the fine-to-coarse injection from the reference tree */
3517: PetscCall(DMPlexGetReferenceTree(coarse, &refTree));
3518: PetscCall(DMCopyDisc(coarse, refTree));
3519: PetscCall(DMSetLocalSection(refTree, NULL));
3520: PetscCall(DMSetDefaultConstraints(refTree, NULL, NULL, NULL));
3521: PetscCall(DMGetDefaultConstraints(refTree, &cSecRef, NULL, NULL));
3522: PetscCall(PetscSectionGetChart(cSecRef, &pRefStart, &pRefEnd));
3523: PetscCall(DMPlexReferenceTreeGetInjector(refTree, &injRef));
3525: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
3526: PetscCall(DMGetLocalSection(fine, &localFine));
3527: PetscCall(DMGetGlobalSection(fine, &globalFine));
3528: PetscCall(PetscSectionGetNumFields(localFine, &numFields));
3529: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
3530: PetscCall(DMGetLocalSection(coarse, &localCoarse));
3531: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
3532: PetscCall(PetscSectionGetMaxDof(localCoarse, &maxDof));
3533: {
3534: PetscInt maxFields = PetscMax(1, numFields) + 1;
3535: PetscCall(PetscMalloc3(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &rowOffsets));
3536: }
3538: PetscCall(DMPlexTransferInjectorTree(coarse, fine, coarseToFine, childIds, NULL, numFields, offsets, &multiRootSec, &rootIndicesSec, &rootIndices, NULL));
3540: PetscCall(PetscMalloc1(maxDof, &parentIndices));
3542: /* count indices */
3543: PetscCall(MatGetLayouts(mat, &rowMap, &colMap));
3544: PetscCall(PetscLayoutSetUp(rowMap));
3545: PetscCall(PetscLayoutSetUp(colMap));
3546: PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd));
3547: PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd));
3548: PetscCall(PetscCalloc2(rowEnd - rowStart, &nnzD, rowEnd - rowStart, &nnzO));
3549: for (p = pStartC; p < pEndC; p++) {
3550: PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
3552: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3553: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
3554: if ((dof - cdof) <= 0) continue;
3555: PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));
3557: rowOffsets[0] = 0;
3558: offsetsCopy[0] = 0;
3559: if (numFields) {
3560: PetscInt f;
3562: for (f = 0; f < numFields; f++) {
3563: PetscInt fDof;
3564: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
3565: rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3566: }
3567: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices));
3568: } else {
3569: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices));
3570: rowOffsets[1] = offsetsCopy[0];
3571: }
3573: PetscCall(PetscSectionGetDof(multiRootSec, p, &numLeaves));
3574: PetscCall(PetscSectionGetOffset(multiRootSec, p, &leafStart));
3575: leafEnd = leafStart + numLeaves;
3576: for (l = leafStart; l < leafEnd; l++) {
3577: PetscInt numIndices, childId, offset;
3578: const PetscInt *childIndices;
3580: PetscCall(PetscSectionGetDof(rootIndicesSec, l, &numIndices));
3581: PetscCall(PetscSectionGetOffset(rootIndicesSec, l, &offset));
3582: childId = rootIndices[offset++];
3583: childIndices = &rootIndices[offset];
3584: numIndices--;
3586: if (childId == -1) { /* equivalent points: scatter */
3587: PetscInt i;
3589: for (i = 0; i < numIndices; i++) {
3590: PetscInt colIndex = childIndices[i];
3591: PetscInt rowIndex = parentIndices[i];
3592: if (rowIndex < 0) continue;
3593: PetscCheck(colIndex >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unconstrained fine and constrained coarse");
3594: if (colIndex >= colStart && colIndex < colEnd) {
3595: nnzD[rowIndex - rowStart] = 1;
3596: } else {
3597: nnzO[rowIndex - rowStart] = 1;
3598: }
3599: }
3600: } else {
3601: PetscInt parentId, f, lim;
3603: PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL));
3605: lim = PetscMax(1, numFields);
3606: offsets[0] = 0;
3607: if (numFields) {
3608: PetscInt f;
3610: for (f = 0; f < numFields; f++) {
3611: PetscInt fDof;
3612: PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof));
3614: offsets[f + 1] = fDof + offsets[f];
3615: }
3616: } else {
3617: PetscInt cDof;
3619: PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof));
3620: offsets[1] = cDof;
3621: }
3622: for (f = 0; f < lim; f++) {
3623: PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1];
3624: PetscInt childStart = offsets[f], childEnd = offsets[f + 1];
3625: PetscInt i, numD = 0, numO = 0;
3627: for (i = childStart; i < childEnd; i++) {
3628: PetscInt colIndex = childIndices[i];
3630: if (colIndex < 0) continue;
3631: if (colIndex >= colStart && colIndex < colEnd) {
3632: numD++;
3633: } else {
3634: numO++;
3635: }
3636: }
3637: for (i = parentStart; i < parentEnd; i++) {
3638: PetscInt rowIndex = parentIndices[i];
3640: if (rowIndex < 0) continue;
3641: nnzD[rowIndex - rowStart] += numD;
3642: nnzO[rowIndex - rowStart] += numO;
3643: }
3644: }
3645: }
3646: }
3647: }
3648: /* preallocate */
3649: PetscCall(MatXAIJSetPreallocation(mat, 1, nnzD, nnzO, NULL, NULL));
3650: PetscCall(PetscFree2(nnzD, nnzO));
3651: /* insert values */
3652: PetscCall(DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree, injRef, &childrenMats));
3653: for (p = pStartC; p < pEndC; p++) {
3654: PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
3656: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3657: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
3658: if ((dof - cdof) <= 0) continue;
3659: PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));
3661: rowOffsets[0] = 0;
3662: offsetsCopy[0] = 0;
3663: if (numFields) {
3664: PetscInt f;
3666: for (f = 0; f < numFields; f++) {
3667: PetscInt fDof;
3668: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
3669: rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3670: }
3671: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices));
3672: } else {
3673: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices));
3674: rowOffsets[1] = offsetsCopy[0];
3675: }
3677: PetscCall(PetscSectionGetDof(multiRootSec, p, &numLeaves));
3678: PetscCall(PetscSectionGetOffset(multiRootSec, p, &leafStart));
3679: leafEnd = leafStart + numLeaves;
3680: for (l = leafStart; l < leafEnd; l++) {
3681: PetscInt numIndices, childId, offset;
3682: const PetscInt *childIndices;
3684: PetscCall(PetscSectionGetDof(rootIndicesSec, l, &numIndices));
3685: PetscCall(PetscSectionGetOffset(rootIndicesSec, l, &offset));
3686: childId = rootIndices[offset++];
3687: childIndices = &rootIndices[offset];
3688: numIndices--;
3690: if (childId == -1) { /* equivalent points: scatter */
3691: PetscInt i;
3693: for (i = 0; i < numIndices; i++) PetscCall(MatSetValue(mat, parentIndices[i], childIndices[i], 1., INSERT_VALUES));
3694: } else {
3695: PetscInt parentId, f, lim;
3697: PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL));
3699: lim = PetscMax(1, numFields);
3700: offsets[0] = 0;
3701: if (numFields) {
3702: PetscInt f;
3704: for (f = 0; f < numFields; f++) {
3705: PetscInt fDof;
3706: PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof));
3708: offsets[f + 1] = fDof + offsets[f];
3709: }
3710: } else {
3711: PetscInt cDof;
3713: PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof));
3714: offsets[1] = cDof;
3715: }
3716: for (f = 0; f < lim; f++) {
3717: PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0];
3718: PetscInt *rowIndices = &parentIndices[rowOffsets[f]];
3719: const PetscInt *colIndices = &childIndices[offsets[f]];
3721: PetscCall(MatSetValues(mat, rowOffsets[f + 1] - rowOffsets[f], rowIndices, offsets[f + 1] - offsets[f], colIndices, childMat, INSERT_VALUES));
3722: }
3723: }
3724: }
3725: }
3726: PetscCall(PetscSectionDestroy(&multiRootSec));
3727: PetscCall(PetscSectionDestroy(&rootIndicesSec));
3728: PetscCall(PetscFree(parentIndices));
3729: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree, injRef, &childrenMats));
3730: PetscCall(PetscFree(rootIndices));
3731: PetscCall(PetscFree3(offsets, offsetsCopy, rowOffsets));
3733: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
3734: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
3735: PetscFunctionReturn(PETSC_SUCCESS);
3736: }
3738: static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom)
3739: {
3740: PetscSF coarseToFineEmbedded;
3741: PetscSection globalCoarse, globalFine;
3742: PetscSection localCoarse, localFine;
3743: PetscSection aSec, cSec;
3744: PetscSection rootValuesSec;
3745: PetscSection leafValuesSec;
3746: PetscScalar *rootValues, *leafValues;
3747: IS aIS;
3748: const PetscInt *anchors;
3749: Mat cMat;
3750: PetscInt numFields;
3751: PetscInt pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd;
3752: PetscInt aStart, aEnd, cStart, cEnd;
3753: PetscInt *maxChildIds;
3754: PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
3755: PetscFV fv = NULL;
3756: PetscInt dim, numFVcomps = -1, fvField = -1;
3757: DM cellDM = NULL, gradDM = NULL;
3758: const PetscScalar *cellGeomArray = NULL;
3759: const PetscScalar *gradArray = NULL;
3761: PetscFunctionBegin;
3762: PetscCall(VecSetOption(vecFine, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE));
3763: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
3764: PetscCall(DMPlexGetSimplexOrBoxCells(coarse, 0, &cellStart, &cellEnd));
3765: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
3766: PetscCall(DMGetGlobalSection(fine, &globalFine));
3767: PetscCall(DMGetCoordinateDim(coarse, &dim));
3768: { /* winnow fine points that don't have global dofs out of the sf */
3769: PetscInt nleaves, l;
3770: const PetscInt *leaves;
3771: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;
3773: PetscCall(PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL));
3775: for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
3776: PetscInt p = leaves ? leaves[l] : l;
3778: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
3779: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
3780: if ((dof - cdof) > 0) numPointsWithDofs++;
3781: }
3782: PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs));
3783: for (l = 0, offset = 0; l < nleaves; l++) {
3784: PetscInt p = leaves ? leaves[l] : l;
3786: PetscCall(PetscSectionGetDof(globalFine, p, &dof));
3787: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof));
3788: if ((dof - cdof) > 0) pointsWithDofs[offset++] = l;
3789: }
3790: PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded));
3791: PetscCall(PetscFree(pointsWithDofs));
3792: }
3793: /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
3794: PetscCall(PetscMalloc1(pEndC - pStartC, &maxChildIds));
3795: for (p = pStartC; p < pEndC; p++) maxChildIds[p - pStartC] = -2;
3796: PetscCall(PetscSFReduceBegin(coarseToFineEmbedded, MPIU_INT, cids, maxChildIds, MPIU_MAX));
3797: PetscCall(PetscSFReduceEnd(coarseToFineEmbedded, MPIU_INT, cids, maxChildIds, MPIU_MAX));
3799: PetscCall(DMGetLocalSection(coarse, &localCoarse));
3800: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
3802: PetscCall(DMPlexGetAnchors(coarse, &aSec, &aIS));
3803: PetscCall(ISGetIndices(aIS, &anchors));
3804: PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd));
3806: PetscCall(DMGetDefaultConstraints(coarse, &cSec, &cMat, NULL));
3807: PetscCall(PetscSectionGetChart(cSec, &cStart, &cEnd));
3809: /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
3810: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootValuesSec));
3811: PetscCall(PetscSectionSetChart(rootValuesSec, pStartC, pEndC));
3812: PetscCall(PetscSectionGetNumFields(localCoarse, &numFields));
3813: {
3814: PetscInt maxFields = PetscMax(1, numFields) + 1;
3815: PetscCall(PetscMalloc7(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &newOffsets, maxFields, &newOffsetsCopy, maxFields, &rowOffsets, maxFields, &numD, maxFields, &numO));
3816: }
3817: if (grad) {
3818: PetscInt i;
3820: PetscCall(VecGetDM(cellGeom, &cellDM));
3821: PetscCall(VecGetArrayRead(cellGeom, &cellGeomArray));
3822: PetscCall(VecGetDM(grad, &gradDM));
3823: PetscCall(VecGetArrayRead(grad, &gradArray));
3824: for (i = 0; i < PetscMax(1, numFields); i++) {
3825: PetscObject obj;
3826: PetscClassId id;
3828: PetscCall(DMGetField(coarse, i, NULL, &obj));
3829: PetscCall(PetscObjectGetClassId(obj, &id));
3830: if (id == PETSCFV_CLASSID) {
3831: fv = (PetscFV)obj;
3832: PetscCall(PetscFVGetNumComponents(fv, &numFVcomps));
3833: fvField = i;
3834: break;
3835: }
3836: }
3837: }
3839: for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
3840: PetscInt dof;
3841: PetscInt maxChildId = maxChildIds[p - pStartC];
3842: PetscInt numValues = 0;
3844: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
3845: if (dof < 0) dof = -(dof + 1);
3846: offsets[0] = 0;
3847: newOffsets[0] = 0;
3848: if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
3849: PetscInt *closure = NULL, closureSize, cl;
3851: PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
3852: for (cl = 0; cl < closureSize; cl++) { /* get the closure */
3853: PetscInt c = closure[2 * cl], clDof;
3855: PetscCall(PetscSectionGetDof(localCoarse, c, &clDof));
3856: numValues += clDof;
3857: }
3858: PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure));
3859: } else if (maxChildId == -1) {
3860: PetscCall(PetscSectionGetDof(localCoarse, p, &numValues));
3861: }
3862: /* we will pack the column indices with the field offsets */
3863: if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) {
3864: /* also send the centroid, and the gradient */
3865: numValues += dim * (1 + numFVcomps);
3866: }
3867: PetscCall(PetscSectionSetDof(rootValuesSec, p, numValues));
3868: }
3869: PetscCall(PetscSectionSetUp(rootValuesSec));
3870: {
3871: PetscInt numRootValues;
3872: const PetscScalar *coarseArray;
3874: PetscCall(PetscSectionGetStorageSize(rootValuesSec, &numRootValues));
3875: PetscCall(PetscMalloc1(numRootValues, &rootValues));
3876: PetscCall(VecGetArrayRead(vecCoarseLocal, &coarseArray));
3877: for (p = pStartC; p < pEndC; p++) {
3878: PetscInt numValues;
3879: PetscInt pValOff;
3880: PetscScalar *pVal;
3881: PetscInt maxChildId = maxChildIds[p - pStartC];
3883: PetscCall(PetscSectionGetDof(rootValuesSec, p, &numValues));
3884: if (!numValues) continue;
3885: PetscCall(PetscSectionGetOffset(rootValuesSec, p, &pValOff));
3886: pVal = &rootValues[pValOff];
3887: if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
3888: PetscInt closureSize = numValues;
3889: PetscCall(DMPlexVecGetClosure(coarse, NULL, vecCoarseLocal, p, &closureSize, &pVal));
3890: if (grad && p >= cellStart && p < cellEnd) {
3891: PetscFVCellGeom *cg;
3892: PetscScalar *gradVals = NULL;
3893: PetscInt i;
3895: pVal += (numValues - dim * (1 + numFVcomps));
3897: PetscCall(DMPlexPointLocalRead(cellDM, p, cellGeomArray, (void *)&cg));
3898: for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i];
3899: pVal += dim;
3900: PetscCall(DMPlexPointGlobalRead(gradDM, p, gradArray, (void *)&gradVals));
3901: for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i];
3902: }
3903: } else if (maxChildId == -1) {
3904: PetscInt lDof, lOff, i;
3906: PetscCall(PetscSectionGetDof(localCoarse, p, &lDof));
3907: PetscCall(PetscSectionGetOffset(localCoarse, p, &lOff));
3908: for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i];
3909: }
3910: }
3911: PetscCall(VecRestoreArrayRead(vecCoarseLocal, &coarseArray));
3912: PetscCall(PetscFree(maxChildIds));
3913: }
3914: {
3915: PetscSF valuesSF;
3916: PetscInt *remoteOffsetsValues, numLeafValues;
3918: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafValuesSec));
3919: PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootValuesSec, &remoteOffsetsValues, leafValuesSec));
3920: PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootValuesSec, remoteOffsetsValues, leafValuesSec, &valuesSF));
3921: PetscCall(PetscSFDestroy(&coarseToFineEmbedded));
3922: PetscCall(PetscFree(remoteOffsetsValues));
3923: PetscCall(PetscSectionGetStorageSize(leafValuesSec, &numLeafValues));
3924: PetscCall(PetscMalloc1(numLeafValues, &leafValues));
3925: PetscCall(PetscSFBcastBegin(valuesSF, MPIU_SCALAR, rootValues, leafValues, MPI_REPLACE));
3926: PetscCall(PetscSFBcastEnd(valuesSF, MPIU_SCALAR, rootValues, leafValues, MPI_REPLACE));
3927: PetscCall(PetscSFDestroy(&valuesSF));
3928: PetscCall(PetscFree(rootValues));
3929: PetscCall(PetscSectionDestroy(&rootValuesSec));
3930: }
3931: PetscCall(DMGetLocalSection(fine, &localFine));
3932: {
3933: PetscInt maxDof;
3934: PetscInt *rowIndices;
3935: DM refTree;
3936: PetscInt **refPointFieldN;
3937: PetscScalar ***refPointFieldMats;
3938: PetscSection refConSec, refAnSec;
3939: PetscInt pRefStart, pRefEnd, leafStart, leafEnd;
3940: PetscScalar *pointWork;
3942: PetscCall(PetscSectionGetMaxDof(localFine, &maxDof));
3943: PetscCall(DMGetWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
3944: PetscCall(DMGetWorkArray(fine, maxDof, MPIU_SCALAR, &pointWork));
3945: PetscCall(DMPlexGetReferenceTree(fine, &refTree));
3946: PetscCall(DMCopyDisc(fine, refTree));
3947: PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
3948: PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL));
3949: PetscCall(DMPlexGetAnchors(refTree, &refAnSec, NULL));
3950: PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd));
3951: PetscCall(PetscSectionGetChart(leafValuesSec, &leafStart, &leafEnd));
3952: PetscCall(DMPlexGetSimplexOrBoxCells(fine, 0, &cellStart, &cellEnd));
3953: for (p = leafStart; p < leafEnd; p++) {
3954: PetscInt gDof, gcDof, gOff, lDof;
3955: PetscInt numValues, pValOff;
3956: PetscInt childId;
3957: const PetscScalar *pVal;
3958: const PetscScalar *fvGradData = NULL;
3960: PetscCall(PetscSectionGetDof(globalFine, p, &gDof));
3961: PetscCall(PetscSectionGetDof(localFine, p, &lDof));
3962: PetscCall(PetscSectionGetConstraintDof(globalFine, p, &gcDof));
3963: if ((gDof - gcDof) <= 0) continue;
3964: PetscCall(PetscSectionGetOffset(globalFine, p, &gOff));
3965: PetscCall(PetscSectionGetDof(leafValuesSec, p, &numValues));
3966: if (!numValues) continue;
3967: PetscCall(PetscSectionGetOffset(leafValuesSec, p, &pValOff));
3968: pVal = &leafValues[pValOff];
3969: offsets[0] = 0;
3970: offsetsCopy[0] = 0;
3971: newOffsets[0] = 0;
3972: newOffsetsCopy[0] = 0;
3973: childId = cids[p - pStartF];
3974: if (numFields) {
3975: PetscInt f;
3976: for (f = 0; f < numFields; f++) {
3977: PetscInt rowDof;
3979: PetscCall(PetscSectionGetFieldDof(localFine, p, f, &rowDof));
3980: offsets[f + 1] = offsets[f] + rowDof;
3981: offsetsCopy[f + 1] = offsets[f + 1];
3982: /* TODO: closure indices */
3983: newOffsets[f + 1] = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]);
3984: }
3985: PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices));
3986: } else {
3987: offsets[0] = 0;
3988: offsets[1] = lDof;
3989: newOffsets[0] = 0;
3990: newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0];
3991: PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices));
3992: }
3993: if (childId == -1) { /* no child interpolation: one nnz per */
3994: PetscCall(VecSetValues(vecFine, numValues, rowIndices, pVal, INSERT_VALUES));
3995: } else {
3996: PetscInt f;
3998: if (grad && p >= cellStart && p < cellEnd) {
3999: numValues -= (dim * (1 + numFVcomps));
4000: fvGradData = &pVal[numValues];
4001: }
4002: for (f = 0; f < PetscMax(1, numFields); f++) {
4003: const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f];
4004: PetscInt numRows = offsets[f + 1] - offsets[f];
4005: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
4006: const PetscScalar *cVal = &pVal[newOffsets[f]];
4007: PetscScalar *rVal = &pointWork[offsets[f]];
4008: PetscInt i, j;
4010: #if 0
4011: 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));
4012: #endif
4013: for (i = 0; i < numRows; i++) {
4014: PetscScalar val = 0.;
4015: for (j = 0; j < numCols; j++) val += childMat[i * numCols + j] * cVal[j];
4016: rVal[i] = val;
4017: }
4018: if (f == fvField && p >= cellStart && p < cellEnd) {
4019: PetscReal centroid[3];
4020: PetscScalar diff[3];
4021: const PetscScalar *parentCentroid = &fvGradData[0];
4022: const PetscScalar *gradient = &fvGradData[dim];
4024: PetscCall(DMPlexComputeCellGeometryFVM(fine, p, NULL, centroid, NULL));
4025: for (i = 0; i < dim; i++) diff[i] = centroid[i] - parentCentroid[i];
4026: for (i = 0; i < numFVcomps; i++) {
4027: PetscScalar val = 0.;
4029: for (j = 0; j < dim; j++) val += gradient[dim * i + j] * diff[j];
4030: rVal[i] += val;
4031: }
4032: }
4033: PetscCall(VecSetValues(vecFine, numRows, &rowIndices[offsets[f]], rVal, INSERT_VALUES));
4034: }
4035: }
4036: }
4037: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN));
4038: PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_SCALAR, &pointWork));
4039: PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_INT, &rowIndices));
4040: }
4041: PetscCall(PetscFree(leafValues));
4042: PetscCall(PetscSectionDestroy(&leafValuesSec));
4043: PetscCall(PetscFree7(offsets, offsetsCopy, newOffsets, newOffsetsCopy, rowOffsets, numD, numO));
4044: PetscCall(ISRestoreIndices(aIS, &anchors));
4045: PetscFunctionReturn(PETSC_SUCCESS);
4046: }
4048: static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids)
4049: {
4050: DM refTree;
4051: PetscSection multiRootSec, rootIndicesSec;
4052: PetscSection globalCoarse, globalFine;
4053: PetscSection localCoarse, localFine;
4054: PetscSection cSecRef;
4055: PetscInt *parentIndices, pRefStart, pRefEnd;
4056: PetscScalar *rootValues, *parentValues;
4057: Mat injRef;
4058: PetscInt numFields, maxDof;
4059: PetscInt pStartC, pEndC, pStartF, pEndF, p;
4060: PetscInt *offsets, *offsetsCopy, *rowOffsets;
4061: PetscLayout rowMap, colMap;
4062: PetscInt rowStart, rowEnd, colStart, colEnd;
4063: PetscScalar ***childrenMats = NULL; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */
4065: PetscFunctionBegin;
4066: /* get the templates for the fine-to-coarse injection from the reference tree */
4067: PetscCall(VecSetOption(vecFine, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE));
4068: PetscCall(VecSetOption(vecCoarse, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE));
4069: PetscCall(DMPlexGetReferenceTree(coarse, &refTree));
4070: PetscCall(DMCopyDisc(coarse, refTree));
4071: PetscCall(DMGetDefaultConstraints(refTree, &cSecRef, NULL, NULL));
4072: PetscCall(PetscSectionGetChart(cSecRef, &pRefStart, &pRefEnd));
4073: PetscCall(DMPlexReferenceTreeGetInjector(refTree, &injRef));
4075: PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF));
4076: PetscCall(DMGetLocalSection(fine, &localFine));
4077: PetscCall(DMGetGlobalSection(fine, &globalFine));
4078: PetscCall(PetscSectionGetNumFields(localFine, &numFields));
4079: PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC));
4080: PetscCall(DMGetLocalSection(coarse, &localCoarse));
4081: PetscCall(DMGetGlobalSection(coarse, &globalCoarse));
4082: PetscCall(PetscSectionGetMaxDof(localCoarse, &maxDof));
4083: {
4084: PetscInt maxFields = PetscMax(1, numFields) + 1;
4085: PetscCall(PetscMalloc3(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &rowOffsets));
4086: }
4088: PetscCall(DMPlexTransferInjectorTree(coarse, fine, coarseToFine, cids, vecFine, numFields, offsets, &multiRootSec, &rootIndicesSec, NULL, &rootValues));
4090: PetscCall(PetscMalloc2(maxDof, &parentIndices, maxDof, &parentValues));
4092: /* count indices */
4093: PetscCall(VecGetLayout(vecFine, &colMap));
4094: PetscCall(VecGetLayout(vecCoarse, &rowMap));
4095: PetscCall(PetscLayoutSetUp(rowMap));
4096: PetscCall(PetscLayoutSetUp(colMap));
4097: PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd));
4098: PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd));
4099: /* insert values */
4100: PetscCall(DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree, injRef, &childrenMats));
4101: for (p = pStartC; p < pEndC; p++) {
4102: PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
4103: PetscBool contribute = PETSC_FALSE;
4105: PetscCall(PetscSectionGetDof(globalCoarse, p, &dof));
4106: PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof));
4107: if ((dof - cdof) <= 0) continue;
4108: PetscCall(PetscSectionGetDof(localCoarse, p, &dof));
4109: PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff));
4111: rowOffsets[0] = 0;
4112: offsetsCopy[0] = 0;
4113: if (numFields) {
4114: PetscInt f;
4116: for (f = 0; f < numFields; f++) {
4117: PetscInt fDof;
4118: PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof));
4119: rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
4120: }
4121: PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices));
4122: } else {
4123: PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices));
4124: rowOffsets[1] = offsetsCopy[0];
4125: }
4127: PetscCall(PetscSectionGetDof(multiRootSec, p, &numLeaves));
4128: PetscCall(PetscSectionGetOffset(multiRootSec, p, &leafStart));
4129: leafEnd = leafStart + numLeaves;
4130: for (l = 0; l < dof; l++) parentValues[l] = 0.;
4131: for (l = leafStart; l < leafEnd; l++) {
4132: PetscInt numIndices, childId, offset;
4133: const PetscScalar *childValues;
4135: PetscCall(PetscSectionGetDof(rootIndicesSec, l, &numIndices));
4136: PetscCall(PetscSectionGetOffset(rootIndicesSec, l, &offset));
4137: childId = (PetscInt)PetscRealPart(rootValues[offset++]);
4138: childValues = &rootValues[offset];
4139: numIndices--;
4141: if (childId == -2) { /* skip */
4142: continue;
4143: } else if (childId == -1) { /* equivalent points: scatter */
4144: PetscInt m;
4146: contribute = PETSC_TRUE;
4147: for (m = 0; m < numIndices; m++) parentValues[m] = childValues[m];
4148: } else { /* contributions from children: sum with injectors from reference tree */
4149: PetscInt parentId, f, lim;
4151: contribute = PETSC_TRUE;
4152: PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL));
4154: lim = PetscMax(1, numFields);
4155: offsets[0] = 0;
4156: if (numFields) {
4157: PetscInt f;
4159: for (f = 0; f < numFields; f++) {
4160: PetscInt fDof;
4161: PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof));
4163: offsets[f + 1] = fDof + offsets[f];
4164: }
4165: } else {
4166: PetscInt cDof;
4168: PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof));
4169: offsets[1] = cDof;
4170: }
4171: for (f = 0; f < lim; f++) {
4172: PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0];
4173: PetscInt n = offsets[f + 1] - offsets[f];
4174: PetscInt m = rowOffsets[f + 1] - rowOffsets[f];
4175: PetscInt i, j;
4176: const PetscScalar *colValues = &childValues[offsets[f]];
4178: for (i = 0; i < m; i++) {
4179: PetscScalar val = 0.;
4180: for (j = 0; j < n; j++) val += childMat[n * i + j] * colValues[j];
4181: parentValues[rowOffsets[f] + i] += val;
4182: }
4183: }
4184: }
4185: }
4186: if (contribute) PetscCall(VecSetValues(vecCoarse, dof, parentIndices, parentValues, INSERT_VALUES));
4187: }
4188: PetscCall(PetscSectionDestroy(&multiRootSec));
4189: PetscCall(PetscSectionDestroy(&rootIndicesSec));
4190: PetscCall(PetscFree2(parentIndices, parentValues));
4191: PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree, injRef, &childrenMats));
4192: PetscCall(PetscFree(rootValues));
4193: PetscCall(PetscFree3(offsets, offsetsCopy, rowOffsets));
4194: PetscFunctionReturn(PETSC_SUCCESS);
4195: }
4197: /*@
4198: DMPlexTransferVecTree - transfer a vector between two meshes that differ from each other by refinement/coarsening
4199: that can be represented by a common reference tree used by both. This routine can be used for a combination of
4200: coarsening and refinement at the same time.
4202: Collective
4204: Input Parameters:
4205: + dmIn - The `DMPLEX` mesh for the input vector
4206: . dmOut - The second `DMPLEX` mesh
4207: . vecIn - The input vector
4208: . sfRefine - A star forest indicating points in the mesh `dmIn` (roots in the star forest) that are parents to points in
4209: the mesh `dmOut` (leaves in the star forest), i.e. where `dmOut` is more refined than `dmIn`
4210: . sfCoarsen - A star forest indicating points in the mesh `dmOut` (roots in the star forest) that are parents to points in
4211: the mesh `dmIn` (leaves in the star forest), i.e. where `dmOut` is more coarsened than `dmIn`
4212: . cidsRefine - The childIds of the points in `dmOut`. These childIds relate back to the reference tree: childid[j] = k implies
4213: that mesh point j of `dmOut` was refined from a point in `dmIn` just as the mesh point k in the reference
4214: tree was refined from its parent. childid[j] = -1 indicates that the point j in `dmOut` is exactly
4215: equivalent to its root in `dmIn`, so no interpolation is necessary. childid[j] = -2 indicates that this
4216: point j in `dmOut` is not a leaf of `sfRefine`.
4217: . cidsCoarsen - The childIds of the points in `dmIn`. These childIds relate back to the reference tree: childid[j] = k implies
4218: that mesh point j of dmIn coarsens to a point in `dmOut` just as the mesh point k in the reference
4219: tree coarsens to its parent. childid[j] = -2 indicates that point j in `dmOut` is not a leaf in `sfCoarsen`.
4220: . useBCs - `PETSC_TRUE` indicates that boundary values should be inserted into `vecIn` before transfer.
4221: - time - Used if boundary values are time dependent.
4223: Output Parameter:
4224: . vecOut - Using interpolation and injection operators calculated on the reference tree, the transferred
4225: projection of `vecIn` from `dmIn` to `dmOut`. Note that any field discretized with a `PetscFV` finite volume
4226: method that uses gradient reconstruction will use reconstructed gradients when interpolating from
4227: coarse points to fine points.
4229: Level: developer
4231: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `Vec`, `PetscFV`, `DMPlexSetReferenceTree()`, `DMPlexGetReferenceTree()`, `PetscFVGetComputeGradients()`
4232: @*/
4233: PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time)
4234: {
4235: PetscFunctionBegin;
4236: PetscCall(VecSet(vecOut, 0.0));
4237: if (sfRefine) {
4238: Vec vecInLocal;
4239: DM dmGrad = NULL;
4240: Vec faceGeom = NULL, cellGeom = NULL, grad = NULL;
4242: PetscCall(DMGetLocalVector(dmIn, &vecInLocal));
4243: PetscCall(VecSet(vecInLocal, 0.0));
4244: {
4245: PetscInt numFields, i;
4247: PetscCall(DMGetNumFields(dmIn, &numFields));
4248: for (i = 0; i < numFields; i++) {
4249: PetscObject obj;
4250: PetscClassId classid;
4252: PetscCall(DMGetField(dmIn, i, NULL, &obj));
4253: PetscCall(PetscObjectGetClassId(obj, &classid));
4254: if (classid == PETSCFV_CLASSID) {
4255: PetscCall(DMPlexGetDataFVM(dmIn, (PetscFV)obj, &cellGeom, &faceGeom, &dmGrad));
4256: break;
4257: }
4258: }
4259: }
4260: if (useBCs) PetscCall(DMPlexInsertBoundaryValues(dmIn, PETSC_TRUE, vecInLocal, time, faceGeom, cellGeom, NULL));
4261: PetscCall(DMGlobalToLocalBegin(dmIn, vecIn, INSERT_VALUES, vecInLocal));
4262: PetscCall(DMGlobalToLocalEnd(dmIn, vecIn, INSERT_VALUES, vecInLocal));
4263: if (dmGrad) {
4264: PetscCall(DMGetGlobalVector(dmGrad, &grad));
4265: PetscCall(DMPlexReconstructGradientsFVM(dmIn, vecInLocal, grad));
4266: }
4267: PetscCall(DMPlexTransferVecTree_Interpolate(dmIn, vecInLocal, dmOut, vecOut, sfRefine, cidsRefine, grad, cellGeom));
4268: PetscCall(DMRestoreLocalVector(dmIn, &vecInLocal));
4269: if (dmGrad) PetscCall(DMRestoreGlobalVector(dmGrad, &grad));
4270: }
4271: if (sfCoarsen) PetscCall(DMPlexTransferVecTree_Inject(dmIn, vecIn, dmOut, vecOut, sfCoarsen, cidsCoarsen));
4272: PetscCall(VecAssemblyBegin(vecOut));
4273: PetscCall(VecAssemblyEnd(vecOut));
4274: PetscFunctionReturn(PETSC_SUCCESS);
4275: }