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