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