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