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