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