Actual source code: plextree.c
petsc-3.10.5 2019-03-28
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;
25: PetscErrorCode ierr;
30: PetscObjectReference((PetscObject)ref);
31: DMDestroy(&mesh->referenceTree);
32: mesh->referenceTree = ref;
33: return(0);
34: }
36: /*@
37: DMPlexGetReferenceTree - get the reference tree for hierarchically non-conforming meshes.
39: Not collective
41: Input Parameters:
42: . dm - The DMPlex object
44: Output Parameters
45: . ref - The reference tree DMPlex object
47: Level: intermediate
49: .seealso: DMPlexSetReferenceTree(), DMPlexCreateDefaultReferenceTree()
50: @*/
51: PetscErrorCode DMPlexGetReferenceTree(DM dm, DM *ref)
52: {
53: DM_Plex *mesh = (DM_Plex *)dm->data;
58: *ref = mesh->referenceTree;
59: return(0);
60: }
62: static PetscErrorCode DMPlexReferenceTreeGetChildSymmetry_Default(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
63: {
64: PetscInt coneSize, dStart, dEnd, dim, ABswap, oAvert, oBvert, ABswapVert;
68: if (parentOrientA == parentOrientB) {
69: if (childOrientB) *childOrientB = childOrientA;
70: if (childB) *childB = childA;
71: return(0);
72: }
73: for (dim = 0; dim < 3; dim++) {
74: DMPlexGetDepthStratum(dm,dim,&dStart,&dEnd);
75: if (parent >= dStart && parent <= dEnd) {
76: break;
77: }
78: }
79: if (dim > 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot perform child symmetry for %d-cells",dim);
80: if (!dim) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A vertex has no children");
81: if (childA < dStart || childA >= dEnd) {
82: /* this is a lower-dimensional child: bootstrap */
83: PetscInt size, i, sA = -1, sB, sOrientB, sConeSize;
84: const PetscInt *supp, *coneA, *coneB, *oA, *oB;
86: DMPlexGetSupportSize(dm,childA,&size);
87: DMPlexGetSupport(dm,childA,&supp);
89: /* find a point sA in supp(childA) that has the same parent */
90: for (i = 0; i < size; i++) {
91: PetscInt sParent;
93: sA = supp[i];
94: if (sA == parent) continue;
95: DMPlexGetTreeParent(dm,sA,&sParent,NULL);
96: if (sParent == parent) {
97: break;
98: }
99: }
100: if (i == size) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"could not find support in children");
101: /* find out which point sB is in an equivalent position to sA under
102: * parentOrientB */
103: DMPlexReferenceTreeGetChildSymmetry_Default(dm,parent,parentOrientA,0,sA,parentOrientB,&sOrientB,&sB);
104: DMPlexGetConeSize(dm,sA,&sConeSize);
105: DMPlexGetCone(dm,sA,&coneA);
106: DMPlexGetCone(dm,sB,&coneB);
107: DMPlexGetConeOrientation(dm,sA,&oA);
108: DMPlexGetConeOrientation(dm,sB,&oB);
109: /* step through the cone of sA in natural order */
110: for (i = 0; i < sConeSize; i++) {
111: if (coneA[i] == childA) {
112: /* if childA is at position i in coneA,
113: * then we want the point that is at sOrientB*i in coneB */
114: PetscInt j = (sOrientB >= 0) ? ((sOrientB + i) % sConeSize) : ((sConeSize -(sOrientB+1) - i) % sConeSize);
115: if (childB) *childB = coneB[j];
116: if (childOrientB) {
117: PetscInt oBtrue;
119: DMPlexGetConeSize(dm,childA,&coneSize);
120: /* compose sOrientB and oB[j] */
121: if (coneSize != 0 && coneSize != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expected a vertex or an edge");
122: /* we may have to flip an edge */
123: oBtrue = coneSize ? ((sOrientB >= 0) ? oB[j] : -(oB[j] + 2)) : 0;
124: ABswap = DihedralSwap(coneSize,oA[i],oBtrue);
125: *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap);
126: }
127: break;
128: }
129: }
130: if (i == sConeSize) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"support cone mismatch");
131: return(0);
132: }
133: /* get the cone size and symmetry swap */
134: DMPlexGetConeSize(dm,parent,&coneSize);
135: ABswap = DihedralSwap(coneSize, parentOrientA, parentOrientB);
136: if (dim == 2) {
137: /* orientations refer to cones: we want them to refer to vertices:
138: * if it's a rotation, they are the same, but if the order is reversed, a
139: * permutation that puts side i first does *not* put vertex i first */
140: oAvert = (parentOrientA >= 0) ? parentOrientA : -((-parentOrientA % coneSize) + 1);
141: oBvert = (parentOrientB >= 0) ? parentOrientB : -((-parentOrientB % coneSize) + 1);
142: ABswapVert = DihedralSwap(coneSize, oAvert, oBvert);
143: } else {
144: ABswapVert = ABswap;
145: }
146: if (childB) {
147: /* assume that each child corresponds to a vertex, in the same order */
148: PetscInt p, posA = -1, numChildren, i;
149: const PetscInt *children;
151: /* count which position the child is in */
152: DMPlexGetTreeChildren(dm,parent,&numChildren,&children);
153: for (i = 0; i < numChildren; i++) {
154: p = children[i];
155: if (p == childA) {
156: posA = i;
157: break;
158: }
159: }
160: if (posA >= coneSize) {
161: /* this is the triangle in the middle of a uniformly refined triangle: it
162: * is invariant */
163: if (dim != 2 || posA != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Expected a middle triangle, got something else");
164: *childB = childA;
165: }
166: else {
167: /* figure out position B by applying ABswapVert */
168: PetscInt posB;
170: posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize -(ABswapVert + 1) - posA) % coneSize);
171: if (childB) *childB = children[posB];
172: }
173: }
174: if (childOrientB) *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap);
175: return(0);
176: }
178: /*@
179: DMPlexReferenceTreeGetChildSymmetry - Given a reference tree, transform a childid and orientation from one parent frame to another
181: Input Parameters:
182: + dm - the reference tree DMPlex object
183: . parent - the parent point
184: . parentOrientA - the reference orientation for describing the parent
185: . childOrientA - the reference orientation for describing the child
186: . childA - the reference childID for describing the child
187: - parentOrientB - the new orientation for describing the parent
189: Output Parameters:
190: + childOrientB - if not NULL, set to the new oreintation for describing the child
191: - childB - if not NULL, the new childID for describing the child
193: Level: developer
195: .seealso: DMPlexGetReferenceTree(), DMPlexSetReferenceTree(), DMPlexSetTree()
196: @*/
197: PetscErrorCode DMPlexReferenceTreeGetChildSymmetry(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
198: {
199: DM_Plex *mesh = (DM_Plex *)dm->data;
204: if (!mesh->getchildsymmetry) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"DMPlexReferenceTreeGetChildSymmetry not implemented");
205: mesh->getchildsymmetry(dm,parent,parentOrientA,childOrientA,childA,parentOrientB,childOrientB,childB);
206: return(0);
207: }
209: static PetscErrorCode DMPlexSetTree_Internal(DM,PetscSection,PetscInt*,PetscInt*,PetscBool,PetscBool);
211: PetscErrorCode DMPlexCreateReferenceTree_Union(DM K, DM Kref, const char *labelName, DM *ref)
212: {
213: MPI_Comm comm;
214: PetscInt dim, p, pStart, pEnd, pRefStart, pRefEnd, d, offset, parentSize, *parents, *childIDs;
215: PetscInt *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts;
216: DMLabel identity, identityRef;
217: PetscSection unionSection, unionConeSection, parentSection;
218: PetscScalar *unionCoords;
219: 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: DMPlexSetTree_Internal(*ref,parentSection,parents,childIDs,PETSC_TRUE,PETSC_FALSE);
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 on comm
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: .keywords: reference cell
450: .seealso: DMPlexSetReferenceTree(), DMPlexGetReferenceTree()
451: @*/
452: PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref)
453: {
454: DM_Plex *mesh;
455: DM K, Kref;
456: PetscInt p, pStart, pEnd;
457: DMLabel identity;
461: #if 1
462: comm = PETSC_COMM_SELF;
463: #endif
464: /* create a reference element */
465: DMPlexCreateReferenceCell(comm, dim, simplex, &K);
466: DMCreateLabel(K, "identity");
467: DMGetLabel(K, "identity", &identity);
468: DMPlexGetChart(K, &pStart, &pEnd);
469: for (p = pStart; p < pEnd; p++) {
470: DMLabelSetValue(identity, p, p);
471: }
472: /* refine it */
473: DMRefine(K,comm,&Kref);
475: /* the reference tree is the union of these two, without duplicating
476: * points that appear in both */
477: DMPlexCreateReferenceTree_Union(K, Kref, "identity", ref);
478: mesh = (DM_Plex *) (*ref)->data;
479: mesh->getchildsymmetry = DMPlexReferenceTreeGetChildSymmetry_Default;
480: DMDestroy(&K);
481: DMDestroy(&Kref);
482: return(0);
483: }
485: static PetscErrorCode DMPlexTreeSymmetrize(DM dm)
486: {
487: DM_Plex *mesh = (DM_Plex *)dm->data;
488: PetscSection childSec, pSec;
489: PetscInt p, pSize, cSize, parMax = PETSC_MIN_INT, parMin = PETSC_MAX_INT;
490: PetscInt *offsets, *children, pStart, pEnd;
495: PetscSectionDestroy(&mesh->childSection);
496: PetscFree(mesh->children);
497: pSec = mesh->parentSection;
498: if (!pSec) return(0);
499: PetscSectionGetStorageSize(pSec,&pSize);
500: for (p = 0; p < pSize; p++) {
501: PetscInt par = mesh->parents[p];
503: parMax = PetscMax(parMax,par+1);
504: parMin = PetscMin(parMin,par);
505: }
506: if (parMin > parMax) {
507: parMin = -1;
508: parMax = -1;
509: }
510: PetscSectionCreate(PetscObjectComm((PetscObject)pSec),&childSec);
511: PetscSectionSetChart(childSec,parMin,parMax);
512: for (p = 0; p < pSize; p++) {
513: PetscInt par = mesh->parents[p];
515: PetscSectionAddDof(childSec,par,1);
516: }
517: PetscSectionSetUp(childSec);
518: PetscSectionGetStorageSize(childSec,&cSize);
519: PetscMalloc1(cSize,&children);
520: PetscCalloc1(parMax-parMin,&offsets);
521: PetscSectionGetChart(pSec,&pStart,&pEnd);
522: for (p = pStart; p < pEnd; p++) {
523: PetscInt dof, off, i;
525: PetscSectionGetDof(pSec,p,&dof);
526: PetscSectionGetOffset(pSec,p,&off);
527: for (i = 0; i < dof; i++) {
528: PetscInt par = mesh->parents[off + i], cOff;
530: PetscSectionGetOffset(childSec,par,&cOff);
531: children[cOff + offsets[par-parMin]++] = p;
532: }
533: }
534: mesh->childSection = childSec;
535: mesh->children = children;
536: PetscFree(offsets);
537: return(0);
538: }
540: static PetscErrorCode AnchorsFlatten (PetscSection section, IS is, PetscSection *sectionNew, IS *isNew)
541: {
542: PetscInt pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL;
543: const PetscInt *vals;
544: PetscSection secNew;
545: PetscBool anyNew, globalAnyNew;
546: PetscBool compress;
550: PetscSectionGetChart(section,&pStart,&pEnd);
551: ISGetLocalSize(is,&size);
552: ISGetIndices(is,&vals);
553: PetscSectionCreate(PetscObjectComm((PetscObject)section),&secNew);
554: PetscSectionSetChart(secNew,pStart,pEnd);
555: for (i = 0; i < size; i++) {
556: PetscInt dof;
558: p = vals[i];
559: if (p < pStart || p >= pEnd) continue;
560: PetscSectionGetDof(section, p, &dof);
561: if (dof) break;
562: }
563: if (i == size) {
564: PetscSectionSetUp(secNew);
565: anyNew = PETSC_FALSE;
566: compress = PETSC_FALSE;
567: sizeNew = 0;
568: }
569: else {
570: anyNew = PETSC_TRUE;
571: for (p = pStart; p < pEnd; p++) {
572: PetscInt dof, off;
574: PetscSectionGetDof(section, p, &dof);
575: PetscSectionGetOffset(section, p, &off);
576: for (i = 0; i < dof; i++) {
577: PetscInt q = vals[off + i], qDof = 0;
579: if (q >= pStart && q < pEnd) {
580: PetscSectionGetDof(section, q, &qDof);
581: }
582: if (qDof) {
583: PetscSectionAddDof(secNew, p, qDof);
584: }
585: else {
586: PetscSectionAddDof(secNew, p, 1);
587: }
588: }
589: }
590: PetscSectionSetUp(secNew);
591: PetscSectionGetStorageSize(secNew,&sizeNew);
592: PetscMalloc1(sizeNew,&valsNew);
593: compress = PETSC_FALSE;
594: for (p = pStart; p < pEnd; p++) {
595: PetscInt dof, off, count, offNew, dofNew;
597: PetscSectionGetDof(section, p, &dof);
598: PetscSectionGetOffset(section, p, &off);
599: PetscSectionGetDof(secNew, p, &dofNew);
600: PetscSectionGetOffset(secNew, p, &offNew);
601: count = 0;
602: for (i = 0; i < dof; i++) {
603: PetscInt q = vals[off + i], qDof = 0, qOff = 0, j;
605: if (q >= pStart && q < pEnd) {
606: PetscSectionGetDof(section, q, &qDof);
607: PetscSectionGetOffset(section, q, &qOff);
608: }
609: if (qDof) {
610: PetscInt oldCount = count;
612: for (j = 0; j < qDof; j++) {
613: PetscInt k, r = vals[qOff + j];
615: for (k = 0; k < oldCount; k++) {
616: if (valsNew[offNew + k] == r) {
617: break;
618: }
619: }
620: if (k == oldCount) {
621: valsNew[offNew + count++] = r;
622: }
623: }
624: }
625: else {
626: PetscInt k, oldCount = count;
628: for (k = 0; k < oldCount; k++) {
629: if (valsNew[offNew + k] == q) {
630: break;
631: }
632: }
633: if (k == oldCount) {
634: valsNew[offNew + count++] = q;
635: }
636: }
637: }
638: if (count < dofNew) {
639: PetscSectionSetDof(secNew, p, count);
640: compress = PETSC_TRUE;
641: }
642: }
643: }
644: ISRestoreIndices(is,&vals);
645: MPIU_Allreduce(&anyNew,&globalAnyNew,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));
646: if (!globalAnyNew) {
647: PetscSectionDestroy(&secNew);
648: *sectionNew = NULL;
649: *isNew = NULL;
650: }
651: else {
652: PetscBool globalCompress;
654: MPIU_Allreduce(&compress,&globalCompress,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));
655: if (compress) {
656: PetscSection secComp;
657: PetscInt *valsComp = NULL;
659: PetscSectionCreate(PetscObjectComm((PetscObject)section),&secComp);
660: PetscSectionSetChart(secComp,pStart,pEnd);
661: for (p = pStart; p < pEnd; p++) {
662: PetscInt dof;
664: PetscSectionGetDof(secNew, p, &dof);
665: PetscSectionSetDof(secComp, p, dof);
666: }
667: PetscSectionSetUp(secComp);
668: PetscSectionGetStorageSize(secComp,&sizeNew);
669: PetscMalloc1(sizeNew,&valsComp);
670: for (p = pStart; p < pEnd; p++) {
671: PetscInt dof, off, offNew, j;
673: PetscSectionGetDof(secNew, p, &dof);
674: PetscSectionGetOffset(secNew, p, &off);
675: PetscSectionGetOffset(secComp, p, &offNew);
676: for (j = 0; j < dof; j++) {
677: valsComp[offNew + j] = valsNew[off + j];
678: }
679: }
680: PetscSectionDestroy(&secNew);
681: secNew = secComp;
682: PetscFree(valsNew);
683: valsNew = valsComp;
684: }
685: ISCreateGeneral(PetscObjectComm((PetscObject)is),sizeNew,valsNew,PETSC_OWN_POINTER,isNew);
686: }
687: return(0);
688: }
690: static PetscErrorCode DMPlexCreateAnchors_Tree(DM dm)
691: {
692: PetscInt p, pStart, pEnd, *anchors, size;
693: PetscInt aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT;
694: PetscSection aSec;
695: DMLabel canonLabel;
696: IS aIS;
701: DMPlexGetChart(dm,&pStart,&pEnd);
702: DMGetLabel(dm,"canonical",&canonLabel);
703: for (p = pStart; p < pEnd; p++) {
704: PetscInt parent;
706: if (canonLabel) {
707: PetscInt canon;
709: DMLabelGetValue(canonLabel,p,&canon);
710: if (p != canon) continue;
711: }
712: DMPlexGetTreeParent(dm,p,&parent,NULL);
713: if (parent != p) {
714: aMin = PetscMin(aMin,p);
715: aMax = PetscMax(aMax,p+1);
716: }
717: }
718: if (aMin > aMax) {
719: aMin = -1;
720: aMax = -1;
721: }
722: PetscSectionCreate(PETSC_COMM_SELF,&aSec);
723: PetscSectionSetChart(aSec,aMin,aMax);
724: for (p = aMin; p < aMax; p++) {
725: PetscInt parent, ancestor = p;
727: if (canonLabel) {
728: PetscInt canon;
730: DMLabelGetValue(canonLabel,p,&canon);
731: if (p != canon) continue;
732: }
733: DMPlexGetTreeParent(dm,p,&parent,NULL);
734: while (parent != ancestor) {
735: ancestor = parent;
736: DMPlexGetTreeParent(dm,ancestor,&parent,NULL);
737: }
738: if (ancestor != p) {
739: PetscInt closureSize, *closure = NULL;
741: DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
742: PetscSectionSetDof(aSec,p,closureSize);
743: DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
744: }
745: }
746: PetscSectionSetUp(aSec);
747: PetscSectionGetStorageSize(aSec,&size);
748: PetscMalloc1(size,&anchors);
749: for (p = aMin; p < aMax; p++) {
750: PetscInt parent, ancestor = p;
752: if (canonLabel) {
753: PetscInt canon;
755: DMLabelGetValue(canonLabel,p,&canon);
756: if (p != canon) continue;
757: }
758: DMPlexGetTreeParent(dm,p,&parent,NULL);
759: while (parent != ancestor) {
760: ancestor = parent;
761: DMPlexGetTreeParent(dm,ancestor,&parent,NULL);
762: }
763: if (ancestor != p) {
764: PetscInt j, closureSize, *closure = NULL, aOff;
766: PetscSectionGetOffset(aSec,p,&aOff);
768: DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
769: for (j = 0; j < closureSize; j++) {
770: anchors[aOff + j] = closure[2*j];
771: }
772: DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
773: }
774: }
775: ISCreateGeneral(PETSC_COMM_SELF,size,anchors,PETSC_OWN_POINTER,&aIS);
776: {
777: PetscSection aSecNew = aSec;
778: IS aISNew = aIS;
780: PetscObjectReference((PetscObject)aSec);
781: PetscObjectReference((PetscObject)aIS);
782: while (aSecNew) {
783: PetscSectionDestroy(&aSec);
784: ISDestroy(&aIS);
785: aSec = aSecNew;
786: aIS = aISNew;
787: aSecNew = NULL;
788: aISNew = NULL;
789: AnchorsFlatten(aSec,aIS,&aSecNew,&aISNew);
790: }
791: }
792: DMPlexSetAnchors(dm,aSec,aIS);
793: PetscSectionDestroy(&aSec);
794: ISDestroy(&aIS);
795: return(0);
796: }
798: static PetscErrorCode DMPlexGetTrueSupportSize(DM dm,PetscInt p,PetscInt *dof,PetscInt *numTrueSupp)
799: {
803: if (numTrueSupp[p] == -1) {
804: PetscInt i, alldof;
805: const PetscInt *supp;
806: PetscInt count = 0;
808: DMPlexGetSupportSize(dm,p,&alldof);
809: DMPlexGetSupport(dm,p,&supp);
810: for (i = 0; i < alldof; i++) {
811: PetscInt q = supp[i], numCones, j;
812: const PetscInt *cone;
814: DMPlexGetConeSize(dm,q,&numCones);
815: DMPlexGetCone(dm,q,&cone);
816: for (j = 0; j < numCones; j++) {
817: if (cone[j] == p) break;
818: }
819: if (j < numCones) count++;
820: }
821: numTrueSupp[p] = count;
822: }
823: *dof = numTrueSupp[p];
824: return(0);
825: }
827: static PetscErrorCode DMPlexTreeExchangeSupports(DM dm)
828: {
829: DM_Plex *mesh = (DM_Plex *)dm->data;
830: PetscSection newSupportSection;
831: PetscInt newSize, *newSupports, pStart, pEnd, p, d, depth;
832: PetscInt *numTrueSupp;
833: PetscInt *offsets;
838: /* symmetrize the hierarchy */
839: DMPlexGetDepth(dm,&depth);
840: PetscSectionCreate(PetscObjectComm((PetscObject)(mesh->supportSection)),&newSupportSection);
841: DMPlexGetChart(dm,&pStart,&pEnd);
842: PetscSectionSetChart(newSupportSection,pStart,pEnd);
843: PetscCalloc1(pEnd,&offsets);
844: PetscMalloc1(pEnd,&numTrueSupp);
845: for (p = 0; p < pEnd; p++) numTrueSupp[p] = -1;
846: /* if a point is in the (true) support of q, it should be in the support of
847: * parent(q) */
848: for (d = 0; d <= depth; d++) {
849: DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);
850: for (p = pStart; p < pEnd; ++p) {
851: PetscInt dof, q, qdof, parent;
853: DMPlexGetTrueSupportSize(dm,p,&dof,numTrueSupp);
854: PetscSectionAddDof(newSupportSection, p, dof);
855: q = p;
856: DMPlexGetTreeParent(dm,q,&parent,NULL);
857: while (parent != q && parent >= pStart && parent < pEnd) {
858: q = parent;
860: DMPlexGetTrueSupportSize(dm,q,&qdof,numTrueSupp);
861: PetscSectionAddDof(newSupportSection,p,qdof);
862: PetscSectionAddDof(newSupportSection,q,dof);
863: DMPlexGetTreeParent(dm,q,&parent,NULL);
864: }
865: }
866: }
867: PetscSectionSetUp(newSupportSection);
868: PetscSectionGetStorageSize(newSupportSection,&newSize);
869: PetscMalloc1(newSize,&newSupports);
870: for (d = 0; d <= depth; d++) {
871: DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);
872: for (p = pStart; p < pEnd; p++) {
873: PetscInt dof, off, q, qdof, qoff, newDof, newOff, newqOff, i, parent;
875: PetscSectionGetDof(mesh->supportSection, p, &dof);
876: PetscSectionGetOffset(mesh->supportSection, p, &off);
877: PetscSectionGetDof(newSupportSection, p, &newDof);
878: PetscSectionGetOffset(newSupportSection, p, &newOff);
879: for (i = 0; i < dof; i++) {
880: PetscInt numCones, j;
881: const PetscInt *cone;
882: PetscInt q = mesh->supports[off + i];
884: DMPlexGetConeSize(dm,q,&numCones);
885: DMPlexGetCone(dm,q,&cone);
886: for (j = 0; j < numCones; j++) {
887: if (cone[j] == p) break;
888: }
889: if (j < numCones) newSupports[newOff+offsets[p]++] = q;
890: }
891: mesh->maxSupportSize = PetscMax(mesh->maxSupportSize,newDof);
893: q = p;
894: DMPlexGetTreeParent(dm,q,&parent,NULL);
895: while (parent != q && parent >= pStart && parent < pEnd) {
896: q = parent;
897: PetscSectionGetDof(mesh->supportSection, q, &qdof);
898: PetscSectionGetOffset(mesh->supportSection, q, &qoff);
899: PetscSectionGetOffset(newSupportSection, q, &newqOff);
900: for (i = 0; i < qdof; i++) {
901: PetscInt numCones, j;
902: const PetscInt *cone;
903: PetscInt r = mesh->supports[qoff + i];
905: DMPlexGetConeSize(dm,r,&numCones);
906: DMPlexGetCone(dm,r,&cone);
907: for (j = 0; j < numCones; j++) {
908: if (cone[j] == q) break;
909: }
910: if (j < numCones) newSupports[newOff+offsets[p]++] = r;
911: }
912: for (i = 0; i < dof; i++) {
913: PetscInt numCones, j;
914: const PetscInt *cone;
915: PetscInt r = mesh->supports[off + i];
917: DMPlexGetConeSize(dm,r,&numCones);
918: DMPlexGetCone(dm,r,&cone);
919: for (j = 0; j < numCones; j++) {
920: if (cone[j] == p) break;
921: }
922: if (j < numCones) newSupports[newqOff+offsets[q]++] = r;
923: }
924: DMPlexGetTreeParent(dm,q,&parent,NULL);
925: }
926: }
927: }
928: PetscSectionDestroy(&mesh->supportSection);
929: mesh->supportSection = newSupportSection;
930: PetscFree(mesh->supports);
931: mesh->supports = newSupports;
932: PetscFree(offsets);
933: PetscFree(numTrueSupp);
935: return(0);
936: }
938: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM,PetscSection,PetscSection,Mat);
939: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM,PetscSection,PetscSection,Mat);
941: static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical, PetscBool exchangeSupports)
942: {
943: DM_Plex *mesh = (DM_Plex *)dm->data;
944: DM refTree;
945: PetscInt size;
951: PetscObjectReference((PetscObject)parentSection);
952: PetscSectionDestroy(&mesh->parentSection);
953: mesh->parentSection = parentSection;
954: PetscSectionGetStorageSize(parentSection,&size);
955: if (parents != mesh->parents) {
956: PetscFree(mesh->parents);
957: PetscMalloc1(size,&mesh->parents);
958: PetscMemcpy(mesh->parents, parents, size * sizeof(*parents));
959: }
960: if (childIDs != mesh->childIDs) {
961: PetscFree(mesh->childIDs);
962: PetscMalloc1(size,&mesh->childIDs);
963: PetscMemcpy(mesh->childIDs, childIDs, size * sizeof(*childIDs));
964: }
965: DMPlexGetReferenceTree(dm,&refTree);
966: if (refTree) {
967: DMLabel canonLabel;
969: DMGetLabel(refTree,"canonical",&canonLabel);
970: if (canonLabel) {
971: PetscInt i;
973: for (i = 0; i < size; i++) {
974: PetscInt canon;
975: DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon);
976: if (canon >= 0) {
977: mesh->childIDs[i] = canon;
978: }
979: }
980: }
981: mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_FromReference;
982: }
983: else {
984: mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_Direct;
985: }
986: DMPlexTreeSymmetrize(dm);
987: if (computeCanonical) {
988: PetscInt d, dim;
990: /* add the canonical label */
991: DMGetDimension(dm,&dim);
992: DMCreateLabel(dm,"canonical");
993: for (d = 0; d <= dim; d++) {
994: PetscInt p, dStart, dEnd, canon = -1, cNumChildren;
995: const PetscInt *cChildren;
997: DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);
998: for (p = dStart; p < dEnd; p++) {
999: DMPlexGetTreeChildren(dm,p,&cNumChildren,&cChildren);
1000: if (cNumChildren) {
1001: canon = p;
1002: break;
1003: }
1004: }
1005: if (canon == -1) continue;
1006: for (p = dStart; p < dEnd; p++) {
1007: PetscInt numChildren, i;
1008: const PetscInt *children;
1010: DMPlexGetTreeChildren(dm,p,&numChildren,&children);
1011: if (numChildren) {
1012: if (numChildren != cNumChildren) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"All parent points in a stratum should have the same number of children: %d != %d", numChildren, cNumChildren);
1013: DMSetLabelValue(dm,"canonical",p,canon);
1014: for (i = 0; i < numChildren; i++) {
1015: DMSetLabelValue(dm,"canonical",children[i],cChildren[i]);
1016: }
1017: }
1018: }
1019: }
1020: }
1021: if (exchangeSupports) {
1022: DMPlexTreeExchangeSupports(dm);
1023: }
1024: mesh->createanchors = DMPlexCreateAnchors_Tree;
1025: /* reset anchors */
1026: DMPlexSetAnchors(dm,NULL,NULL);
1027: return(0);
1028: }
1030: /*@
1031: DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points. This routine also creates
1032: the point-to-point constraints determined by the tree: a point is constained to the points in the closure of its
1033: tree root.
1035: Collective on dm
1037: Input Parameters:
1038: + dm - the DMPlex object
1039: . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
1040: offset indexes the parent and childID list; the reference count of parentSection is incremented
1041: . parents - a list of the point parents; copied, can be destroyed
1042: - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
1043: the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed
1045: Level: intermediate
1047: .seealso: DMPlexGetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren()
1048: @*/
1049: PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
1050: {
1054: DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_FALSE,PETSC_TRUE);
1055: return(0);
1056: }
1058: /*@
1059: DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points.
1060: Collective on dm
1062: Input Parameters:
1063: . dm - the DMPlex object
1065: Output Parameters:
1066: + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
1067: offset indexes the parent and childID list
1068: . parents - a list of the point parents
1069: . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
1070: the child corresponds to the point in the reference tree with index childID
1071: . childSection - the inverse of the parent section
1072: - children - a list of the point children
1074: Level: intermediate
1076: .seealso: DMPlexSetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren()
1077: @*/
1078: PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[])
1079: {
1080: DM_Plex *mesh = (DM_Plex *)dm->data;
1084: if (parentSection) *parentSection = mesh->parentSection;
1085: if (parents) *parents = mesh->parents;
1086: if (childIDs) *childIDs = mesh->childIDs;
1087: if (childSection) *childSection = mesh->childSection;
1088: if (children) *children = mesh->children;
1089: return(0);
1090: }
1092: /*@
1093: DMPlexGetTreeParent - get the parent of a point in the tree describing the point hierarchy (not the DAG)
1095: Input Parameters:
1096: + dm - the DMPlex object
1097: - point - the query point
1099: Output Parameters:
1100: + parent - if not NULL, set to the parent of the point, or the point itself if the point does not have a parent
1101: - childID - if not NULL, set to the child ID of the point with respect to its parent, or 0 if the point
1102: does not have a parent
1104: Level: intermediate
1106: .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeChildren()
1107: @*/
1108: PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID)
1109: {
1110: DM_Plex *mesh = (DM_Plex *)dm->data;
1111: PetscSection pSec;
1116: pSec = mesh->parentSection;
1117: if (pSec && point >= pSec->pStart && point < pSec->pEnd) {
1118: PetscInt dof;
1120: PetscSectionGetDof (pSec, point, &dof);
1121: if (dof) {
1122: PetscInt off;
1124: PetscSectionGetOffset (pSec, point, &off);
1125: if (parent) *parent = mesh->parents[off];
1126: if (childID) *childID = mesh->childIDs[off];
1127: return(0);
1128: }
1129: }
1130: if (parent) {
1131: *parent = point;
1132: }
1133: if (childID) {
1134: *childID = 0;
1135: }
1136: return(0);
1137: }
1139: /*@C
1140: DMPlexGetTreeChildren - get the children of a point in the tree describing the point hierarchy (not the DAG)
1142: Input Parameters:
1143: + dm - the DMPlex object
1144: - point - the query point
1146: Output Parameters:
1147: + numChildren - if not NULL, set to the number of children
1148: - children - if not NULL, set to a list children, or set to NULL if the point has no children
1150: Level: intermediate
1152: Fortran Notes:
1153: Since it returns an array, this routine is only available in Fortran 90, and you must
1154: include petsc.h90 in your code.
1156: .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeParent()
1157: @*/
1158: PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[])
1159: {
1160: DM_Plex *mesh = (DM_Plex *)dm->data;
1161: PetscSection childSec;
1162: PetscInt dof = 0;
1167: childSec = mesh->childSection;
1168: if (childSec && point >= childSec->pStart && point < childSec->pEnd) {
1169: PetscSectionGetDof (childSec, point, &dof);
1170: }
1171: if (numChildren) *numChildren = dof;
1172: if (children) {
1173: if (dof) {
1174: PetscInt off;
1176: PetscSectionGetOffset (childSec, point, &off);
1177: *children = &mesh->children[off];
1178: }
1179: else {
1180: *children = NULL;
1181: }
1182: }
1183: return(0);
1184: }
1186: 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)
1187: {
1188: PetscInt f, b, p, c, offset, qPoints;
1192: PetscSpaceEvaluate(space,nPoints,points,work,NULL,NULL);
1193: for (f = 0, offset = 0; f < nFunctionals; f++) {
1194: qPoints = pointsPerFn[f];
1195: for (b = 0; b < nBasis; b++) {
1196: PetscScalar val = 0.;
1198: for (p = 0; p < qPoints; p++) {
1199: for (c = 0; c < nComps; c++) {
1200: val += work[((offset + p) * nBasis + b) * nComps + c] * weights[(offset + p) * nComps + c];
1201: }
1202: }
1203: MatSetValue(basisAtPoints,b,f,val,INSERT_VALUES);
1204: }
1205: offset += qPoints;
1206: }
1207: MatAssemblyBegin(basisAtPoints,MAT_FINAL_ASSEMBLY);
1208: MatAssemblyEnd(basisAtPoints,MAT_FINAL_ASSEMBLY);
1209: return(0);
1210: }
1212: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM dm, PetscSection section, PetscSection cSec, Mat cMat)
1213: {
1214: PetscDS ds;
1215: PetscInt spdim;
1216: PetscInt numFields, f, c, cStart, cEnd, pStart, pEnd, conStart, conEnd;
1217: const PetscInt *anchors;
1218: PetscSection aSec;
1219: PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJparent, detJ, detJparent;
1220: IS aIS;
1224: DMPlexGetChart(dm,&pStart,&pEnd);
1225: DMGetDS(dm,&ds);
1226: PetscDSGetNumFields(ds,&numFields);
1227: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
1228: DMPlexGetAnchors(dm,&aSec,&aIS);
1229: ISGetIndices(aIS,&anchors);
1230: PetscSectionGetChart(cSec,&conStart,&conEnd);
1231: DMGetDimension(dm,&spdim);
1232: PetscMalloc6(spdim,&v0,spdim,&v0parent,spdim,&vtmp,spdim*spdim,&J,spdim*spdim,&Jparent,spdim*spdim,&invJparent);
1234: for (f = 0; f < numFields; f++) {
1235: PetscObject disc;
1236: PetscClassId id;
1237: PetscSpace bspace;
1238: PetscDualSpace dspace;
1239: PetscInt i, j, k, nPoints, Nc, offset;
1240: PetscInt fSize, maxDof;
1241: PetscReal *weights, *pointsRef, *pointsReal, *work;
1242: PetscScalar *scwork, *X;
1243: PetscInt *sizes, *workIndRow, *workIndCol;
1244: Mat Amat, Bmat, Xmat;
1245: const PetscInt *numDof = NULL;
1246: const PetscInt ***perms = NULL;
1247: const PetscScalar ***flips = NULL;
1249: PetscDSGetDiscretization(ds,f,&disc);
1250: PetscObjectGetClassId(disc,&id);
1251: if (id == PETSCFE_CLASSID) {
1252: PetscFE fe = (PetscFE) disc;
1254: PetscFEGetBasisSpace(fe,&bspace);
1255: PetscFEGetDualSpace(fe,&dspace);
1256: PetscDualSpaceGetDimension(dspace,&fSize);
1257: PetscFEGetNumComponents(fe,&Nc);
1258: }
1259: else if (id == PETSCFV_CLASSID) {
1260: PetscFV fv = (PetscFV) disc;
1262: PetscFVGetNumComponents(fv,&Nc);
1263: PetscSpaceCreate(PetscObjectComm((PetscObject)fv),&bspace);
1264: PetscSpaceSetType(bspace,PETSCSPACEPOLYNOMIAL);
1265: PetscSpaceSetDegree(bspace,0,PETSC_DETERMINE);
1266: PetscSpaceSetNumComponents(bspace,Nc);
1267: PetscSpaceSetNumVariables(bspace,spdim);
1268: PetscSpaceSetUp(bspace);
1269: PetscFVGetDualSpace(fv,&dspace);
1270: PetscDualSpaceGetDimension(dspace,&fSize);
1271: }
1272: else SETERRQ1(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id);
1273: PetscDualSpaceGetNumDof(dspace,&numDof);
1274: for (i = 0, maxDof = 0; i <= spdim; i++) {maxDof = PetscMax(maxDof,numDof[i]);}
1275: PetscDualSpaceGetSymmetries(dspace,&perms,&flips);
1277: MatCreate(PETSC_COMM_SELF,&Amat);
1278: MatSetSizes(Amat,fSize,fSize,fSize,fSize);
1279: MatSetType(Amat,MATSEQDENSE);
1280: MatSetUp(Amat);
1281: MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Bmat);
1282: MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Xmat);
1283: nPoints = 0;
1284: for (i = 0; i < fSize; i++) {
1285: PetscInt qPoints, thisNc;
1286: PetscQuadrature quad;
1288: PetscDualSpaceGetFunctional(dspace,i,&quad);
1289: PetscQuadratureGetData(quad,NULL,&thisNc,&qPoints,NULL,NULL);
1290: if (thisNc != Nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Functional dim %D does not much basis dim %D\n",thisNc,Nc);
1291: nPoints += qPoints;
1292: }
1293: PetscMalloc7(fSize,&sizes,nPoints*Nc,&weights,spdim*nPoints,&pointsRef,spdim*nPoints,&pointsReal,nPoints*fSize*Nc,&work,maxDof,&workIndRow,maxDof,&workIndCol);
1294: PetscMalloc1(maxDof * maxDof,&scwork);
1295: offset = 0;
1296: for (i = 0; i < fSize; i++) {
1297: PetscInt qPoints;
1298: const PetscReal *p, *w;
1299: PetscQuadrature quad;
1301: PetscDualSpaceGetFunctional(dspace,i,&quad);
1302: PetscQuadratureGetData(quad,NULL,NULL,&qPoints,&p,&w);
1303: PetscMemcpy(weights+Nc*offset,w,Nc*qPoints*sizeof(*w));
1304: PetscMemcpy(pointsRef+spdim*offset,p,spdim*qPoints*sizeof(*p));
1305: sizes[i] = qPoints;
1306: offset += qPoints;
1307: }
1308: EvaluateBasis(bspace,fSize,fSize,Nc,nPoints,sizes,pointsRef,weights,work,Amat);
1309: MatLUFactor(Amat,NULL,NULL,NULL);
1310: for (c = cStart; c < cEnd; c++) {
1311: PetscInt parent;
1312: PetscInt closureSize, closureSizeP, *closure = NULL, *closureP = NULL;
1313: PetscInt *childOffsets, *parentOffsets;
1315: DMPlexGetTreeParent(dm,c,&parent,NULL);
1316: if (parent == c) continue;
1317: DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1318: for (i = 0; i < closureSize; i++) {
1319: PetscInt p = closure[2*i];
1320: PetscInt conDof;
1322: if (p < conStart || p >= conEnd) continue;
1323: if (numFields) {
1324: PetscSectionGetFieldDof(cSec,p,f,&conDof);
1325: }
1326: else {
1327: PetscSectionGetDof(cSec,p,&conDof);
1328: }
1329: if (conDof) break;
1330: }
1331: if (i == closureSize) {
1332: DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1333: continue;
1334: }
1336: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, NULL, &detJ);
1337: DMPlexComputeCellGeometryFEM(dm, parent, NULL, v0parent, Jparent, invJparent, &detJparent);
1338: for (i = 0; i < nPoints; i++) {
1339: const PetscReal xi0[3] = {-1.,-1.,-1.};
1341: CoordinatesRefToReal(spdim, spdim, xi0, v0, J, &pointsRef[i*spdim],vtmp);
1342: CoordinatesRealToRef(spdim, spdim, xi0, v0parent, invJparent, vtmp, &pointsReal[i*spdim]);
1343: }
1344: EvaluateBasis(bspace,fSize,fSize,Nc,nPoints,sizes,pointsReal,weights,work,Bmat);
1345: MatMatSolve(Amat,Bmat,Xmat);
1346: MatDenseGetArray(Xmat,&X);
1347: DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);
1348: PetscMalloc2(closureSize+1,&childOffsets,closureSizeP+1,&parentOffsets);
1349: childOffsets[0] = 0;
1350: for (i = 0; i < closureSize; i++) {
1351: PetscInt p = closure[2*i];
1352: PetscInt dof;
1354: if (numFields) {
1355: PetscSectionGetFieldDof(section,p,f,&dof);
1356: }
1357: else {
1358: PetscSectionGetDof(section,p,&dof);
1359: }
1360: childOffsets[i+1]=childOffsets[i]+dof;
1361: }
1362: parentOffsets[0] = 0;
1363: for (i = 0; i < closureSizeP; i++) {
1364: PetscInt p = closureP[2*i];
1365: PetscInt dof;
1367: if (numFields) {
1368: PetscSectionGetFieldDof(section,p,f,&dof);
1369: }
1370: else {
1371: PetscSectionGetDof(section,p,&dof);
1372: }
1373: parentOffsets[i+1]=parentOffsets[i]+dof;
1374: }
1375: for (i = 0; i < closureSize; i++) {
1376: PetscInt conDof, conOff, aDof, aOff, nWork;
1377: PetscInt p = closure[2*i];
1378: PetscInt o = closure[2*i+1];
1379: const PetscInt *perm;
1380: const PetscScalar *flip;
1382: if (p < conStart || p >= conEnd) continue;
1383: if (numFields) {
1384: PetscSectionGetFieldDof(cSec,p,f,&conDof);
1385: PetscSectionGetFieldOffset(cSec,p,f,&conOff);
1386: }
1387: else {
1388: PetscSectionGetDof(cSec,p,&conDof);
1389: PetscSectionGetOffset(cSec,p,&conOff);
1390: }
1391: if (!conDof) continue;
1392: perm = (perms && perms[i]) ? perms[i][o] : NULL;
1393: flip = (flips && flips[i]) ? flips[i][o] : NULL;
1394: PetscSectionGetDof(aSec,p,&aDof);
1395: PetscSectionGetOffset(aSec,p,&aOff);
1396: nWork = childOffsets[i+1]-childOffsets[i];
1397: for (k = 0; k < aDof; k++) {
1398: PetscInt a = anchors[aOff + k];
1399: PetscInt aSecDof, aSecOff;
1401: if (numFields) {
1402: PetscSectionGetFieldDof(section,a,f,&aSecDof);
1403: PetscSectionGetFieldOffset(section,a,f,&aSecOff);
1404: }
1405: else {
1406: PetscSectionGetDof(section,a,&aSecDof);
1407: PetscSectionGetOffset(section,a,&aSecOff);
1408: }
1409: if (!aSecDof) continue;
1411: for (j = 0; j < closureSizeP; j++) {
1412: PetscInt q = closureP[2*j];
1413: PetscInt oq = closureP[2*j+1];
1415: if (q == a) {
1416: PetscInt r, s, nWorkP;
1417: const PetscInt *permP;
1418: const PetscScalar *flipP;
1420: permP = (perms && perms[j]) ? perms[j][oq] : NULL;
1421: flipP = (flips && flips[j]) ? flips[j][oq] : NULL;
1422: nWorkP = parentOffsets[j+1]-parentOffsets[j];
1423: /* get a copy of the child-to-anchor portion of the matrix, and transpose so that rows correspond to the
1424: * child and columns correspond to the anchor: BUT the maxrix returned by MatDenseGetArray is
1425: * column-major, so transpose-transpose = do nothing */
1426: for (r = 0; r < nWork; r++) {
1427: for (s = 0; s < nWorkP; s++) {
1428: scwork[r * nWorkP + s] = X[fSize * (r + childOffsets[i]) + (s + parentOffsets[j])];
1429: }
1430: }
1431: for (r = 0; r < nWork; r++) {workIndRow[perm ? perm[r] : r] = conOff + r;}
1432: for (s = 0; s < nWorkP; s++) {workIndCol[permP ? permP[s] : s] = aSecOff + s;}
1433: if (flip) {
1434: for (r = 0; r < nWork; r++) {
1435: for (s = 0; s < nWorkP; s++) {
1436: scwork[r * nWorkP + s] *= flip[r];
1437: }
1438: }
1439: }
1440: if (flipP) {
1441: for (r = 0; r < nWork; r++) {
1442: for (s = 0; s < nWorkP; s++) {
1443: scwork[r * nWorkP + s] *= flipP[s];
1444: }
1445: }
1446: }
1447: MatSetValues(cMat,nWork,workIndRow,nWorkP,workIndCol,scwork,INSERT_VALUES);
1448: break;
1449: }
1450: }
1451: }
1452: }
1453: MatDenseRestoreArray(Xmat,&X);
1454: PetscFree2(childOffsets,parentOffsets);
1455: DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1456: DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);
1457: }
1458: MatDestroy(&Amat);
1459: MatDestroy(&Bmat);
1460: MatDestroy(&Xmat);
1461: PetscFree(scwork);
1462: PetscFree7(sizes,weights,pointsRef,pointsReal,work,workIndRow,workIndCol);
1463: if (id == PETSCFV_CLASSID) {
1464: PetscSpaceDestroy(&bspace);
1465: }
1466: }
1467: MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
1468: MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
1469: PetscFree6(v0,v0parent,vtmp,J,Jparent,invJparent);
1470: ISRestoreIndices(aIS,&anchors);
1472: return(0);
1473: }
1475: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1476: {
1477: Mat refCmat;
1478: PetscDS ds;
1479: PetscInt numFields, maxFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, **refPointFieldN;
1480: PetscScalar ***refPointFieldMats;
1481: PetscSection refConSec, refAnSec, refSection;
1482: IS refAnIS;
1483: const PetscInt *refAnchors;
1484: const PetscInt **perms;
1485: const PetscScalar **flips;
1486: PetscErrorCode ierr;
1489: DMGetDS(refTree,&ds);
1490: PetscDSGetNumFields(ds,&numFields);
1491: maxFields = PetscMax(1,numFields);
1492: DMGetDefaultConstraints(refTree,&refConSec,&refCmat);
1493: DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);
1494: ISGetIndices(refAnIS,&refAnchors);
1495: DMGetSection(refTree,&refSection);
1496: PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1497: PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);
1498: PetscMalloc1(pRefEnd-pRefStart,&refPointFieldN);
1499: PetscSectionGetMaxDof(refConSec,&maxDof);
1500: PetscSectionGetMaxDof(refAnSec,&maxAnDof);
1501: PetscMalloc1(maxDof,&rows);
1502: PetscMalloc1(maxDof*maxAnDof,&cols);
1503: for (p = pRefStart; p < pRefEnd; p++) {
1504: PetscInt parent, closureSize, *closure = NULL, pDof;
1506: DMPlexGetTreeParent(refTree,p,&parent,NULL);
1507: PetscSectionGetDof(refConSec,p,&pDof);
1508: if (!pDof || parent == p) continue;
1510: PetscMalloc1(maxFields,&refPointFieldMats[p-pRefStart]);
1511: PetscCalloc1(maxFields,&refPointFieldN[p-pRefStart]);
1512: DMPlexGetTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);
1513: for (f = 0; f < maxFields; f++) {
1514: PetscInt cDof, cOff, numCols, r, i;
1516: if (f < numFields) {
1517: PetscSectionGetFieldDof(refConSec,p,f,&cDof);
1518: PetscSectionGetFieldOffset(refConSec,p,f,&cOff);
1519: PetscSectionGetFieldPointSyms(refSection,f,closureSize,closure,&perms,&flips);
1520: } else {
1521: PetscSectionGetDof(refConSec,p,&cDof);
1522: PetscSectionGetOffset(refConSec,p,&cOff);
1523: PetscSectionGetPointSyms(refSection,closureSize,closure,&perms,&flips);
1524: }
1526: for (r = 0; r < cDof; r++) {
1527: rows[r] = cOff + r;
1528: }
1529: numCols = 0;
1530: for (i = 0; i < closureSize; i++) {
1531: PetscInt q = closure[2*i];
1532: PetscInt aDof, aOff, j;
1533: const PetscInt *perm = perms ? perms[i] : NULL;
1535: if (numFields) {
1536: PetscSectionGetFieldDof(refSection,q,f,&aDof);
1537: PetscSectionGetFieldOffset(refSection,q,f,&aOff);
1538: }
1539: else {
1540: PetscSectionGetDof(refSection,q,&aDof);
1541: PetscSectionGetOffset(refSection,q,&aOff);
1542: }
1544: for (j = 0; j < aDof; j++) {
1545: cols[numCols++] = aOff + (perm ? perm[j] : j);
1546: }
1547: }
1548: refPointFieldN[p-pRefStart][f] = numCols;
1549: PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);
1550: MatGetValues(refCmat,cDof,rows,numCols,cols,refPointFieldMats[p-pRefStart][f]);
1551: if (flips) {
1552: PetscInt colOff = 0;
1554: for (i = 0; i < closureSize; i++) {
1555: PetscInt q = closure[2*i];
1556: PetscInt aDof, aOff, j;
1557: const PetscScalar *flip = flips ? flips[i] : NULL;
1559: if (numFields) {
1560: PetscSectionGetFieldDof(refSection,q,f,&aDof);
1561: PetscSectionGetFieldOffset(refSection,q,f,&aOff);
1562: }
1563: else {
1564: PetscSectionGetDof(refSection,q,&aDof);
1565: PetscSectionGetOffset(refSection,q,&aOff);
1566: }
1567: if (flip) {
1568: PetscInt k;
1569: for (k = 0; k < cDof; k++) {
1570: for (j = 0; j < aDof; j++) {
1571: refPointFieldMats[p-pRefStart][f][k * numCols + colOff + j] *= flip[j];
1572: }
1573: }
1574: }
1575: colOff += aDof;
1576: }
1577: }
1578: if (numFields) {
1579: PetscSectionRestoreFieldPointSyms(refSection,f,closureSize,closure,&perms,&flips);
1580: } else {
1581: PetscSectionRestorePointSyms(refSection,closureSize,closure,&perms,&flips);
1582: }
1583: }
1584: DMPlexRestoreTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);
1585: }
1586: *childrenMats = refPointFieldMats;
1587: *childrenN = refPointFieldN;
1588: ISRestoreIndices(refAnIS,&refAnchors);
1589: PetscFree(rows);
1590: PetscFree(cols);
1591: return(0);
1592: }
1594: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1595: {
1596: PetscDS ds;
1597: PetscInt **refPointFieldN;
1598: PetscScalar ***refPointFieldMats;
1599: PetscInt numFields, maxFields, pRefStart, pRefEnd, p, f;
1600: PetscSection refConSec;
1604: refPointFieldN = *childrenN;
1605: *childrenN = NULL;
1606: refPointFieldMats = *childrenMats;
1607: *childrenMats = NULL;
1608: DMGetDS(refTree,&ds);
1609: PetscDSGetNumFields(ds,&numFields);
1610: maxFields = PetscMax(1,numFields);
1611: DMGetDefaultConstraints(refTree,&refConSec,NULL);
1612: PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1613: for (p = pRefStart; p < pRefEnd; p++) {
1614: PetscInt parent, pDof;
1616: DMPlexGetTreeParent(refTree,p,&parent,NULL);
1617: PetscSectionGetDof(refConSec,p,&pDof);
1618: if (!pDof || parent == p) continue;
1620: for (f = 0; f < maxFields; f++) {
1621: PetscInt cDof;
1623: if (numFields) {
1624: PetscSectionGetFieldDof(refConSec,p,f,&cDof);
1625: }
1626: else {
1627: PetscSectionGetDof(refConSec,p,&cDof);
1628: }
1630: PetscFree(refPointFieldMats[p - pRefStart][f]);
1631: }
1632: PetscFree(refPointFieldMats[p - pRefStart]);
1633: PetscFree(refPointFieldN[p - pRefStart]);
1634: }
1635: PetscFree(refPointFieldMats);
1636: PetscFree(refPointFieldN);
1637: return(0);
1638: }
1640: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm, PetscSection section, PetscSection conSec, Mat cMat)
1641: {
1642: DM refTree;
1643: PetscDS ds;
1644: Mat refCmat;
1645: PetscInt numFields, maxFields, f, pRefStart, pRefEnd, p, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN;
1646: PetscScalar ***refPointFieldMats, *pointWork;
1647: PetscSection refConSec, refAnSec, anSec;
1648: IS refAnIS, anIS;
1649: const PetscInt *anchors;
1654: DMGetDS(dm,&ds);
1655: PetscDSGetNumFields(ds,&numFields);
1656: maxFields = PetscMax(1,numFields);
1657: DMPlexGetReferenceTree(dm,&refTree);
1658: DMSetDS(refTree,ds);
1659: DMGetDefaultConstraints(refTree,&refConSec,&refCmat);
1660: DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);
1661: DMPlexGetAnchors(dm,&anSec,&anIS);
1662: ISGetIndices(anIS,&anchors);
1663: PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1664: PetscSectionGetChart(conSec,&conStart,&conEnd);
1665: PetscSectionGetMaxDof(refConSec,&maxDof);
1666: PetscSectionGetMaxDof(refAnSec,&maxAnDof);
1667: PetscMalloc1(maxDof*maxDof*maxAnDof,&pointWork);
1669: /* step 1: get submats for every constrained point in the reference tree */
1670: DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
1672: /* step 2: compute the preorder */
1673: DMPlexGetChart(dm,&pStart,&pEnd);
1674: PetscMalloc2(pEnd-pStart,&perm,pEnd-pStart,&iperm);
1675: for (p = pStart; p < pEnd; p++) {
1676: perm[p - pStart] = p;
1677: iperm[p - pStart] = p-pStart;
1678: }
1679: for (p = 0; p < pEnd - pStart;) {
1680: PetscInt point = perm[p];
1681: PetscInt parent;
1683: DMPlexGetTreeParent(dm,point,&parent,NULL);
1684: if (parent == point) {
1685: p++;
1686: }
1687: else {
1688: PetscInt size, closureSize, *closure = NULL, i;
1690: DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1691: for (i = 0; i < closureSize; i++) {
1692: PetscInt q = closure[2*i];
1693: if (iperm[q-pStart] > iperm[point-pStart]) {
1694: /* swap */
1695: perm[p] = q;
1696: perm[iperm[q-pStart]] = point;
1697: iperm[point-pStart] = iperm[q-pStart];
1698: iperm[q-pStart] = p;
1699: break;
1700: }
1701: }
1702: size = closureSize;
1703: DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1704: if (i == size) {
1705: p++;
1706: }
1707: }
1708: }
1710: /* step 3: fill the constraint matrix */
1711: /* we are going to use a preorder progressive fill strategy. Mat doesn't
1712: * allow progressive fill without assembly, so we are going to set up the
1713: * values outside of the Mat first.
1714: */
1715: {
1716: PetscInt nRows, row, nnz;
1717: PetscBool done;
1718: const PetscInt *ia, *ja;
1719: PetscScalar *vals;
1721: MatGetRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);
1722: if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not get RowIJ of constraint matrix");
1723: nnz = ia[nRows];
1724: /* malloc and then zero rows right before we fill them: this way valgrind
1725: * can tell if we are doing progressive fill in the wrong order */
1726: PetscMalloc1(nnz,&vals);
1727: for (p = 0; p < pEnd - pStart; p++) {
1728: PetscInt parent, childid, closureSize, *closure = NULL;
1729: PetscInt point = perm[p], pointDof;
1731: DMPlexGetTreeParent(dm,point,&parent,&childid);
1732: if ((point < conStart) || (point >= conEnd) || (parent == point)) continue;
1733: PetscSectionGetDof(conSec,point,&pointDof);
1734: if (!pointDof) continue;
1735: DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1736: for (f = 0; f < maxFields; f++) {
1737: PetscInt cDof, cOff, numCols, numFillCols, i, r, matOffset, offset;
1738: PetscScalar *pointMat;
1739: const PetscInt **perms;
1740: const PetscScalar **flips;
1742: if (numFields) {
1743: PetscSectionGetFieldDof(conSec,point,f,&cDof);
1744: PetscSectionGetFieldOffset(conSec,point,f,&cOff);
1745: }
1746: else {
1747: PetscSectionGetDof(conSec,point,&cDof);
1748: PetscSectionGetOffset(conSec,point,&cOff);
1749: }
1750: if (!cDof) continue;
1751: if (numFields) {PetscSectionGetFieldPointSyms(section,f,closureSize,closure,&perms,&flips);}
1752: else {PetscSectionGetPointSyms(section,closureSize,closure,&perms,&flips);}
1754: /* make sure that every row for this point is the same size */
1755: #if defined(PETSC_USE_DEBUG)
1756: for (r = 0; r < cDof; r++) {
1757: if (cDof > 1 && r) {
1758: if ((ia[cOff+r+1]-ia[cOff+r]) != (ia[cOff+r]-ia[cOff+r-1])) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Two point rows have different nnz: %D vs. %D", (ia[cOff+r+1]-ia[cOff+r]), (ia[cOff+r]-ia[cOff+r-1]));
1759: }
1760: }
1761: #endif
1762: /* zero rows */
1763: for (i = ia[cOff] ; i< ia[cOff+cDof];i++) {
1764: vals[i] = 0.;
1765: }
1766: matOffset = ia[cOff];
1767: numFillCols = ia[cOff+1] - matOffset;
1768: pointMat = refPointFieldMats[childid-pRefStart][f];
1769: numCols = refPointFieldN[childid-pRefStart][f];
1770: offset = 0;
1771: for (i = 0; i < closureSize; i++) {
1772: PetscInt q = closure[2*i];
1773: PetscInt aDof, aOff, j, k, qConDof, qConOff;
1774: const PetscInt *perm = perms ? perms[i] : NULL;
1775: const PetscScalar *flip = flips ? flips[i] : NULL;
1777: qConDof = qConOff = 0;
1778: if (numFields) {
1779: PetscSectionGetFieldDof(section,q,f,&aDof);
1780: PetscSectionGetFieldOffset(section,q,f,&aOff);
1781: if (q >= conStart && q < conEnd) {
1782: PetscSectionGetFieldDof(conSec,q,f,&qConDof);
1783: PetscSectionGetFieldOffset(conSec,q,f,&qConOff);
1784: }
1785: }
1786: else {
1787: PetscSectionGetDof(section,q,&aDof);
1788: PetscSectionGetOffset(section,q,&aOff);
1789: if (q >= conStart && q < conEnd) {
1790: PetscSectionGetDof(conSec,q,&qConDof);
1791: PetscSectionGetOffset(conSec,q,&qConOff);
1792: }
1793: }
1794: if (!aDof) continue;
1795: if (qConDof) {
1796: /* this point has anchors: its rows of the matrix should already
1797: * be filled, thanks to preordering */
1798: /* first multiply into pointWork, then set in matrix */
1799: PetscInt aMatOffset = ia[qConOff];
1800: PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset;
1801: for (r = 0; r < cDof; r++) {
1802: for (j = 0; j < aNumFillCols; j++) {
1803: PetscScalar inVal = 0;
1804: for (k = 0; k < aDof; k++) {
1805: PetscInt col = perm ? perm[k] : k;
1807: inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j] * (flip ? flip[col] : 1.);
1808: }
1809: pointWork[r * aNumFillCols + j] = inVal;
1810: }
1811: }
1812: /* assume that the columns are sorted, spend less time searching */
1813: for (j = 0, k = 0; j < aNumFillCols; j++) {
1814: PetscInt col = ja[aMatOffset + j];
1815: for (;k < numFillCols; k++) {
1816: if (ja[matOffset + k] == col) {
1817: break;
1818: }
1819: }
1820: if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, col);
1821: for (r = 0; r < cDof; r++) {
1822: vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j];
1823: }
1824: }
1825: }
1826: else {
1827: /* find where to put this portion of pointMat into the matrix */
1828: for (k = 0; k < numFillCols; k++) {
1829: if (ja[matOffset + k] == aOff) {
1830: break;
1831: }
1832: }
1833: if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, aOff);
1834: for (r = 0; r < cDof; r++) {
1835: for (j = 0; j < aDof; j++) {
1836: PetscInt col = perm ? perm[j] : j;
1838: vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + j] * (flip ? flip[col] : 1.);
1839: }
1840: }
1841: }
1842: offset += aDof;
1843: }
1844: if (numFields) {
1845: PetscSectionRestoreFieldPointSyms(section,f,closureSize,closure,&perms,&flips);
1846: } else {
1847: PetscSectionRestorePointSyms(section,closureSize,closure,&perms,&flips);
1848: }
1849: }
1850: DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1851: }
1852: for (row = 0; row < nRows; row++) {
1853: MatSetValues(cMat,1,&row,ia[row+1]-ia[row],&ja[ia[row]],&vals[ia[row]],INSERT_VALUES);
1854: }
1855: MatRestoreRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);
1856: if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not restore RowIJ of constraint matrix");
1857: MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
1858: MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
1859: PetscFree(vals);
1860: }
1862: /* clean up */
1863: ISRestoreIndices(anIS,&anchors);
1864: PetscFree2(perm,iperm);
1865: PetscFree(pointWork);
1866: DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
1867: return(0);
1868: }
1870: /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of
1871: * a non-conforming mesh. Local refinement comes later */
1872: PetscErrorCode DMPlexTreeRefineCell (DM dm, PetscInt cell, DM *ncdm)
1873: {
1874: DM K;
1875: PetscMPIInt rank;
1876: PetscInt dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd;
1877: PetscInt numNewCones, *newConeSizes, *newCones, *newOrientations;
1878: PetscInt *Kembedding;
1879: PetscInt *cellClosure=NULL, nc;
1880: PetscScalar *newVertexCoords;
1881: PetscInt numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset;
1882: PetscSection parentSection;
1886: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1887: DMGetDimension(dm,&dim);
1888: DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm);
1889: DMSetDimension(*ncdm,dim);
1891: DMPlexGetChart(dm, &pStart, &pEnd);
1892: PetscSectionCreate(PetscObjectComm((PetscObject)dm),&parentSection);
1893: DMPlexGetReferenceTree(dm,&K);
1894: if (!rank) {
1895: /* compute the new charts */
1896: PetscMalloc5(dim+1,&pNewCount,dim+1,&pNewStart,dim+1,&pNewEnd,dim+1,&pOldStart,dim+1,&pOldEnd);
1897: offset = 0;
1898: for (d = 0; d <= dim; d++) {
1899: PetscInt pOldCount, kStart, kEnd, k;
1901: pNewStart[d] = offset;
1902: DMPlexGetHeightStratum(dm,d,&pOldStart[d],&pOldEnd[d]);
1903: DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
1904: pOldCount = pOldEnd[d] - pOldStart[d];
1905: /* adding the new points */
1906: pNewCount[d] = pOldCount + kEnd - kStart;
1907: if (!d) {
1908: /* removing the cell */
1909: pNewCount[d]--;
1910: }
1911: for (k = kStart; k < kEnd; k++) {
1912: PetscInt parent;
1913: DMPlexGetTreeParent(K,k,&parent,NULL);
1914: if (parent == k) {
1915: /* avoid double counting points that won't actually be new */
1916: pNewCount[d]--;
1917: }
1918: }
1919: pNewEnd[d] = pNewStart[d] + pNewCount[d];
1920: offset = pNewEnd[d];
1922: }
1923: if (cell < pOldStart[0] || cell >= pOldEnd[0]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"%d not in cell range [%d, %d)", cell, pOldStart[0], pOldEnd[0]);
1924: /* get the current closure of the cell that we are removing */
1925: DMPlexGetTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);
1927: PetscMalloc1(pNewEnd[dim],&newConeSizes);
1928: {
1929: PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j;
1931: DMPlexGetChart(K,&kStart,&kEnd);
1932: PetscMalloc4(kEnd-kStart,&Kembedding,kEnd-kStart,&perm,kEnd-kStart,&iperm,kEnd-kStart,&preOrient);
1934: for (k = kStart; k < kEnd; k++) {
1935: perm[k - kStart] = k;
1936: iperm [k - kStart] = k - kStart;
1937: preOrient[k - kStart] = 0;
1938: }
1940: DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);
1941: for (j = 1; j < closureSizeK; j++) {
1942: PetscInt parentOrientA = closureK[2*j+1];
1943: PetscInt parentOrientB = cellClosure[2*j+1];
1944: PetscInt p, q;
1946: p = closureK[2*j];
1947: q = cellClosure[2*j];
1948: for (d = 0; d <= dim; d++) {
1949: if (q >= pOldStart[d] && q < pOldEnd[d]) {
1950: Kembedding[p] = (q - pOldStart[d]) + pNewStart[d];
1951: }
1952: }
1953: if (parentOrientA != parentOrientB) {
1954: PetscInt numChildren, i;
1955: const PetscInt *children;
1957: DMPlexGetTreeChildren(K,p,&numChildren,&children);
1958: for (i = 0; i < numChildren; i++) {
1959: PetscInt kPerm, oPerm;
1961: k = children[i];
1962: DMPlexReferenceTreeGetChildSymmetry(K,p,parentOrientA,0,k,parentOrientB,&oPerm,&kPerm);
1963: /* perm = what refTree position I'm in */
1964: perm[kPerm-kStart] = k;
1965: /* iperm = who is at this position */
1966: iperm[k-kStart] = kPerm-kStart;
1967: preOrient[kPerm-kStart] = oPerm;
1968: }
1969: }
1970: }
1971: DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);
1972: }
1973: PetscSectionSetChart(parentSection,0,pNewEnd[dim]);
1974: offset = 0;
1975: numNewCones = 0;
1976: for (d = 0; d <= dim; d++) {
1977: PetscInt kStart, kEnd, k;
1978: PetscInt p;
1979: PetscInt size;
1981: for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1982: /* skip cell 0 */
1983: if (p == cell) continue;
1984: /* old cones to new cones */
1985: DMPlexGetConeSize(dm,p,&size);
1986: newConeSizes[offset++] = size;
1987: numNewCones += size;
1988: }
1990: DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
1991: for (k = kStart; k < kEnd; k++) {
1992: PetscInt kParent;
1994: DMPlexGetTreeParent(K,k,&kParent,NULL);
1995: if (kParent != k) {
1996: Kembedding[k] = offset;
1997: DMPlexGetConeSize(K,k,&size);
1998: newConeSizes[offset++] = size;
1999: numNewCones += size;
2000: if (kParent != 0) {
2001: PetscSectionSetDof(parentSection,Kembedding[k],1);
2002: }
2003: }
2004: }
2005: }
2007: PetscSectionSetUp(parentSection);
2008: PetscSectionGetStorageSize(parentSection,&numPointsWithParents);
2009: PetscMalloc2(numNewCones,&newCones,numNewCones,&newOrientations);
2010: PetscMalloc2(numPointsWithParents,&parents,numPointsWithParents,&childIDs);
2012: /* fill new cones */
2013: offset = 0;
2014: for (d = 0; d <= dim; d++) {
2015: PetscInt kStart, kEnd, k, l;
2016: PetscInt p;
2017: PetscInt size;
2018: const PetscInt *cone, *orientation;
2020: for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
2021: /* skip cell 0 */
2022: if (p == cell) continue;
2023: /* old cones to new cones */
2024: DMPlexGetConeSize(dm,p,&size);
2025: DMPlexGetCone(dm,p,&cone);
2026: DMPlexGetConeOrientation(dm,p,&orientation);
2027: for (l = 0; l < size; l++) {
2028: newCones[offset] = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1];
2029: newOrientations[offset++] = orientation[l];
2030: }
2031: }
2033: DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
2034: for (k = kStart; k < kEnd; k++) {
2035: PetscInt kPerm = perm[k], kParent;
2036: PetscInt preO = preOrient[k];
2038: DMPlexGetTreeParent(K,k,&kParent,NULL);
2039: if (kParent != k) {
2040: /* embed new cones */
2041: DMPlexGetConeSize(K,k,&size);
2042: DMPlexGetCone(K,kPerm,&cone);
2043: DMPlexGetConeOrientation(K,kPerm,&orientation);
2044: for (l = 0; l < size; l++) {
2045: PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size -(preO + 1) - l) % size);
2046: PetscInt newO, lSize, oTrue;
2048: q = iperm[cone[m]];
2049: newCones[offset] = Kembedding[q];
2050: DMPlexGetConeSize(K,q,&lSize);
2051: oTrue = orientation[m];
2052: oTrue = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2);
2053: newO = DihedralCompose(lSize,oTrue,preOrient[q]);
2054: newOrientations[offset++] = newO;
2055: }
2056: if (kParent != 0) {
2057: PetscInt newPoint = Kembedding[kParent];
2058: PetscSectionGetOffset(parentSection,Kembedding[k],&pOffset);
2059: parents[pOffset] = newPoint;
2060: childIDs[pOffset] = k;
2061: }
2062: }
2063: }
2064: }
2066: PetscMalloc1(dim*(pNewEnd[dim]-pNewStart[dim]),&newVertexCoords);
2068: /* fill coordinates */
2069: offset = 0;
2070: {
2071: PetscInt kStart, kEnd, l;
2072: PetscSection vSection;
2073: PetscInt v;
2074: Vec coords;
2075: PetscScalar *coordvals;
2076: PetscInt dof, off;
2077: PetscReal v0[3], J[9], detJ;
2079: #if defined(PETSC_USE_DEBUG)
2080: {
2081: PetscInt k;
2082: DMPlexGetHeightStratum(K,0,&kStart,&kEnd);
2083: for (k = kStart; k < kEnd; k++) {
2084: DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ);
2085: if (detJ <= 0.) SETERRQ1 (PETSC_COMM_SELF,PETSC_ERR_PLIB,"reference tree cell %d has bad determinant",k);
2086: }
2087: }
2088: #endif
2089: DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ);
2090: DMGetCoordinateSection(dm,&vSection);
2091: DMGetCoordinatesLocal(dm,&coords);
2092: VecGetArray(coords,&coordvals);
2093: for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) {
2095: PetscSectionGetDof(vSection,v,&dof);
2096: PetscSectionGetOffset(vSection,v,&off);
2097: for (l = 0; l < dof; l++) {
2098: newVertexCoords[offset++] = coordvals[off + l];
2099: }
2100: }
2101: VecRestoreArray(coords,&coordvals);
2103: DMGetCoordinateSection(K,&vSection);
2104: DMGetCoordinatesLocal(K,&coords);
2105: VecGetArray(coords,&coordvals);
2106: DMPlexGetDepthStratum(K,0,&kStart,&kEnd);
2107: for (v = kStart; v < kEnd; v++) {
2108: PetscReal coord[3], newCoord[3];
2109: PetscInt vPerm = perm[v];
2110: PetscInt kParent;
2111: const PetscReal xi0[3] = {-1.,-1.,-1.};
2113: DMPlexGetTreeParent(K,v,&kParent,NULL);
2114: if (kParent != v) {
2115: /* this is a new vertex */
2116: PetscSectionGetOffset(vSection,vPerm,&off);
2117: for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off+l]);
2118: CoordinatesRefToReal(dim, dim, xi0, v0, J, coord, newCoord);
2119: for (l = 0; l < dim; ++l) newVertexCoords[offset+l] = newCoord[l];
2120: offset += dim;
2121: }
2122: }
2123: VecRestoreArray(coords,&coordvals);
2124: }
2126: /* need to reverse the order of pNewCount: vertices first, cells last */
2127: for (d = 0; d < (dim + 1) / 2; d++) {
2128: PetscInt tmp;
2130: tmp = pNewCount[d];
2131: pNewCount[d] = pNewCount[dim - d];
2132: pNewCount[dim - d] = tmp;
2133: }
2135: DMPlexCreateFromDAG(*ncdm,dim,pNewCount,newConeSizes,newCones,newOrientations,newVertexCoords);
2136: DMPlexSetReferenceTree(*ncdm,K);
2137: DMPlexSetTree(*ncdm,parentSection,parents,childIDs);
2139: /* clean up */
2140: DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);
2141: PetscFree5(pNewCount,pNewStart,pNewEnd,pOldStart,pOldEnd);
2142: PetscFree(newConeSizes);
2143: PetscFree2(newCones,newOrientations);
2144: PetscFree(newVertexCoords);
2145: PetscFree2(parents,childIDs);
2146: PetscFree4(Kembedding,perm,iperm,preOrient);
2147: }
2148: else {
2149: PetscInt p, counts[4];
2150: PetscInt *coneSizes, *cones, *orientations;
2151: Vec coordVec;
2152: PetscScalar *coords;
2154: for (d = 0; d <= dim; d++) {
2155: PetscInt dStart, dEnd;
2157: DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);
2158: counts[d] = dEnd - dStart;
2159: }
2160: PetscMalloc1(pEnd-pStart,&coneSizes);
2161: for (p = pStart; p < pEnd; p++) {
2162: DMPlexGetConeSize(dm,p,&coneSizes[p-pStart]);
2163: }
2164: DMPlexGetCones(dm, &cones);
2165: DMPlexGetConeOrientations(dm, &orientations);
2166: DMGetCoordinatesLocal(dm,&coordVec);
2167: VecGetArray(coordVec,&coords);
2169: PetscSectionSetChart(parentSection,pStart,pEnd);
2170: PetscSectionSetUp(parentSection);
2171: DMPlexCreateFromDAG(*ncdm,dim,counts,coneSizes,cones,orientations,NULL);
2172: DMPlexSetReferenceTree(*ncdm,K);
2173: DMPlexSetTree(*ncdm,parentSection,NULL,NULL);
2174: VecRestoreArray(coordVec,&coords);
2175: }
2176: PetscSectionDestroy(&parentSection);
2178: return(0);
2179: }
2181: PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
2182: {
2183: PetscSF coarseToFineEmbedded;
2184: PetscSection globalCoarse, globalFine;
2185: PetscSection localCoarse, localFine;
2186: PetscSection aSec, cSec;
2187: PetscSection rootIndicesSec, rootMatricesSec;
2188: PetscSection leafIndicesSec, leafMatricesSec;
2189: PetscInt *rootIndices, *leafIndices;
2190: PetscScalar *rootMatrices, *leafMatrices;
2191: IS aIS;
2192: const PetscInt *anchors;
2193: Mat cMat;
2194: PetscInt numFields, maxFields;
2195: PetscInt pStartC, pEndC, pStartF, pEndF, p;
2196: PetscInt aStart, aEnd, cStart, cEnd;
2197: PetscInt *maxChildIds;
2198: PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
2199: const PetscInt ***perms;
2200: const PetscScalar ***flips;
2201: PetscErrorCode ierr;
2204: DMPlexGetChart(coarse,&pStartC,&pEndC);
2205: DMPlexGetChart(fine,&pStartF,&pEndF);
2206: DMGetGlobalSection(fine,&globalFine);
2207: { /* winnow fine points that don't have global dofs out of the sf */
2208: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, nleaves, l;
2209: const PetscInt *leaves;
2211: PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);
2212: for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
2213: p = leaves ? leaves[l] : l;
2214: PetscSectionGetDof(globalFine,p,&dof);
2215: PetscSectionGetConstraintDof(globalFine,p,&cdof);
2216: if ((dof - cdof) > 0) {
2217: numPointsWithDofs++;
2218: }
2219: }
2220: PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
2221: for (l = 0, offset = 0; l < nleaves; l++) {
2222: p = leaves ? leaves[l] : l;
2223: PetscSectionGetDof(globalFine,p,&dof);
2224: PetscSectionGetConstraintDof(globalFine,p,&cdof);
2225: if ((dof - cdof) > 0) {
2226: pointsWithDofs[offset++] = l;
2227: }
2228: }
2229: PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
2230: PetscFree(pointsWithDofs);
2231: }
2232: /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
2233: PetscMalloc1(pEndC-pStartC,&maxChildIds);
2234: for (p = pStartC; p < pEndC; p++) {
2235: maxChildIds[p - pStartC] = -2;
2236: }
2237: PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);
2238: PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);
2240: DMGetSection(coarse,&localCoarse);
2241: DMGetGlobalSection(coarse,&globalCoarse);
2243: DMPlexGetAnchors(coarse,&aSec,&aIS);
2244: ISGetIndices(aIS,&anchors);
2245: PetscSectionGetChart(aSec,&aStart,&aEnd);
2247: DMGetDefaultConstraints(coarse,&cSec,&cMat);
2248: PetscSectionGetChart(cSec,&cStart,&cEnd);
2250: /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
2251: PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);
2252: PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootMatricesSec);
2253: PetscSectionSetChart(rootIndicesSec,pStartC,pEndC);
2254: PetscSectionSetChart(rootMatricesSec,pStartC,pEndC);
2255: PetscSectionGetNumFields(localCoarse,&numFields);
2256: maxFields = PetscMax(1,numFields);
2257: PetscMalloc7(maxFields+1,&offsets,maxFields+1,&offsetsCopy,maxFields+1,&newOffsets,maxFields+1,&newOffsetsCopy,maxFields+1,&rowOffsets,maxFields+1,&numD,maxFields+1,&numO);
2258: PetscMalloc2(maxFields+1,(PetscInt****)&perms,maxFields+1,(PetscScalar****)&flips);
2259: PetscMemzero((void *) perms, (maxFields+1) * sizeof(const PetscInt **));
2260: PetscMemzero((void *) flips, (maxFields+1) * sizeof(const PetscScalar **));
2262: for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
2263: PetscInt dof, matSize = 0;
2264: PetscInt aDof = 0;
2265: PetscInt cDof = 0;
2266: PetscInt maxChildId = maxChildIds[p - pStartC];
2267: PetscInt numRowIndices = 0;
2268: PetscInt numColIndices = 0;
2269: PetscInt f;
2271: PetscSectionGetDof(globalCoarse,p,&dof);
2272: if (dof < 0) {
2273: dof = -(dof + 1);
2274: }
2275: if (p >= aStart && p < aEnd) {
2276: PetscSectionGetDof(aSec,p,&aDof);
2277: }
2278: if (p >= cStart && p < cEnd) {
2279: PetscSectionGetDof(cSec,p,&cDof);
2280: }
2281: for (f = 0; f <= numFields; f++) offsets[f] = 0;
2282: for (f = 0; f <= numFields; f++) newOffsets[f] = 0;
2283: if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
2284: PetscInt *closure = NULL, closureSize, cl;
2286: DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);
2287: for (cl = 0; cl < closureSize; cl++) { /* get the closure */
2288: PetscInt c = closure[2 * cl], clDof;
2290: PetscSectionGetDof(localCoarse,c,&clDof);
2291: numRowIndices += clDof;
2292: for (f = 0; f < numFields; f++) {
2293: PetscSectionGetFieldDof(localCoarse,c,f,&clDof);
2294: offsets[f + 1] += clDof;
2295: }
2296: }
2297: for (f = 0; f < numFields; f++) {
2298: offsets[f + 1] += offsets[f];
2299: newOffsets[f + 1] = offsets[f + 1];
2300: }
2301: /* get the number of indices needed and their field offsets */
2302: DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,NULL,NULL,NULL,&numColIndices,NULL,NULL,newOffsets,PETSC_FALSE);
2303: DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);
2304: if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */
2305: numColIndices = numRowIndices;
2306: matSize = 0;
2307: }
2308: else if (numFields) { /* we send one submat for each field: sum their sizes */
2309: matSize = 0;
2310: for (f = 0; f < numFields; f++) {
2311: PetscInt numRow, numCol;
2313: numRow = offsets[f + 1] - offsets[f];
2314: numCol = newOffsets[f + 1] - newOffsets[f];
2315: matSize += numRow * numCol;
2316: }
2317: }
2318: else {
2319: matSize = numRowIndices * numColIndices;
2320: }
2321: } else if (maxChildId == -1) {
2322: if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */
2323: PetscInt aOff, a;
2325: PetscSectionGetOffset(aSec,p,&aOff);
2326: for (f = 0; f < numFields; f++) {
2327: PetscInt fDof;
2329: PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
2330: offsets[f+1] = fDof;
2331: }
2332: for (a = 0; a < aDof; a++) {
2333: PetscInt anchor = anchors[a + aOff], aLocalDof;
2335: PetscSectionGetDof(localCoarse,anchor,&aLocalDof);
2336: numColIndices += aLocalDof;
2337: for (f = 0; f < numFields; f++) {
2338: PetscInt fDof;
2340: PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);
2341: newOffsets[f+1] += fDof;
2342: }
2343: }
2344: if (numFields) {
2345: matSize = 0;
2346: for (f = 0; f < numFields; f++) {
2347: matSize += offsets[f+1] * newOffsets[f+1];
2348: }
2349: }
2350: else {
2351: matSize = numColIndices * dof;
2352: }
2353: }
2354: else { /* no children, and no constraints on dofs: just get the global indices */
2355: numColIndices = dof;
2356: matSize = 0;
2357: }
2358: }
2359: /* we will pack the column indices with the field offsets */
2360: PetscSectionSetDof(rootIndicesSec,p,numColIndices ? numColIndices+2*numFields : 0);
2361: PetscSectionSetDof(rootMatricesSec,p,matSize);
2362: }
2363: PetscSectionSetUp(rootIndicesSec);
2364: PetscSectionSetUp(rootMatricesSec);
2365: {
2366: PetscInt numRootIndices, numRootMatrices;
2368: PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);
2369: PetscSectionGetStorageSize(rootMatricesSec,&numRootMatrices);
2370: PetscMalloc2(numRootIndices,&rootIndices,numRootMatrices,&rootMatrices);
2371: for (p = pStartC; p < pEndC; p++) {
2372: PetscInt numRowIndices, numColIndices, matSize, dof;
2373: PetscInt pIndOff, pMatOff, f;
2374: PetscInt *pInd;
2375: PetscInt maxChildId = maxChildIds[p - pStartC];
2376: PetscScalar *pMat = NULL;
2378: PetscSectionGetDof(rootIndicesSec,p,&numColIndices);
2379: if (!numColIndices) {
2380: continue;
2381: }
2382: for (f = 0; f <= numFields; f++) {
2383: offsets[f] = 0;
2384: newOffsets[f] = 0;
2385: offsetsCopy[f] = 0;
2386: newOffsetsCopy[f] = 0;
2387: }
2388: numColIndices -= 2 * numFields;
2389: PetscSectionGetOffset(rootIndicesSec,p,&pIndOff);
2390: pInd = &(rootIndices[pIndOff]);
2391: PetscSectionGetDof(rootMatricesSec,p,&matSize);
2392: if (matSize) {
2393: PetscSectionGetOffset(rootMatricesSec,p,&pMatOff);
2394: pMat = &rootMatrices[pMatOff];
2395: }
2396: PetscSectionGetDof(globalCoarse,p,&dof);
2397: if (dof < 0) {
2398: dof = -(dof + 1);
2399: }
2400: if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
2401: PetscInt i, j;
2402: PetscInt numRowIndices = matSize / numColIndices;
2404: if (!numRowIndices) { /* don't need to calculate the mat, just the indices */
2405: PetscInt numIndices, *indices;
2406: DMPlexGetClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);
2407: if (numIndices != numColIndices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"mismatching constraint indices calculations");
2408: for (i = 0; i < numColIndices; i++) {
2409: pInd[i] = indices[i];
2410: }
2411: for (i = 0; i < numFields; i++) {
2412: pInd[numColIndices + i] = offsets[i+1];
2413: pInd[numColIndices + numFields + i] = offsets[i+1];
2414: }
2415: DMPlexRestoreClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);
2416: }
2417: else {
2418: PetscInt closureSize, *closure = NULL, cl;
2419: PetscScalar *pMatIn, *pMatModified;
2420: PetscInt numPoints,*points;
2422: DMGetWorkArray(coarse,numRowIndices * numRowIndices,MPIU_SCALAR,&pMatIn);
2423: for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */
2424: for (j = 0; j < numRowIndices; j++) {
2425: pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.;
2426: }
2427: }
2428: DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);
2429: for (f = 0; f < maxFields; f++) {
2430: if (numFields) {PetscSectionGetFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f]);}
2431: else {PetscSectionGetPointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f]);}
2432: }
2433: if (numFields) {
2434: for (cl = 0; cl < closureSize; cl++) {
2435: PetscInt c = closure[2 * cl];
2437: for (f = 0; f < numFields; f++) {
2438: PetscInt fDof;
2440: PetscSectionGetFieldDof(localCoarse,c,f,&fDof);
2441: offsets[f + 1] += fDof;
2442: }
2443: }
2444: for (f = 0; f < numFields; f++) {
2445: offsets[f + 1] += offsets[f];
2446: newOffsets[f + 1] = offsets[f + 1];
2447: }
2448: }
2449: /* TODO : flips here ? */
2450: /* apply hanging node constraints on the right, get the new points and the new offsets */
2451: DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,perms,pMatIn,&numPoints,NULL,&points,&pMatModified,newOffsets,PETSC_FALSE);
2452: for (f = 0; f < maxFields; f++) {
2453: if (numFields) {PetscSectionRestoreFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f]);}
2454: else {PetscSectionRestorePointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f]);}
2455: }
2456: for (f = 0; f < maxFields; f++) {
2457: if (numFields) {PetscSectionGetFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f]);}
2458: else {PetscSectionGetPointSyms(localCoarse,numPoints,points,&perms[f],&flips[f]);}
2459: }
2460: if (!numFields) {
2461: for (i = 0; i < numRowIndices * numColIndices; i++) {
2462: pMat[i] = pMatModified[i];
2463: }
2464: }
2465: else {
2466: PetscInt i, j, count;
2467: for (f = 0, count = 0; f < numFields; f++) {
2468: for (i = offsets[f]; i < offsets[f+1]; i++) {
2469: for (j = newOffsets[f]; j < newOffsets[f+1]; j++, count++) {
2470: pMat[count] = pMatModified[i * numColIndices + j];
2471: }
2472: }
2473: }
2474: }
2475: DMRestoreWorkArray(coarse,numRowIndices * numColIndices,MPIU_SCALAR,&pMatModified);
2476: DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);
2477: DMRestoreWorkArray(coarse,numRowIndices * numColIndices,MPIU_SCALAR,&pMatIn);
2478: if (numFields) {
2479: for (f = 0; f < numFields; f++) {
2480: pInd[numColIndices + f] = offsets[f+1];
2481: pInd[numColIndices + numFields + f] = newOffsets[f+1];
2482: }
2483: for (cl = 0; cl < numPoints; cl++) {
2484: PetscInt globalOff, c = points[2*cl];
2485: PetscSectionGetOffset(globalCoarse, c, &globalOff);
2486: DMPlexGetIndicesPointFields_Internal(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perms, cl, pInd);
2487: }
2488: } else {
2489: for (cl = 0; cl < numPoints; cl++) {
2490: PetscInt c = points[2*cl], globalOff;
2491: const PetscInt *perm = perms[0] ? perms[0][cl] : NULL;
2493: PetscSectionGetOffset(globalCoarse, c, &globalOff);
2494: DMPlexGetIndicesPoint_Internal(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perm, pInd);
2495: }
2496: }
2497: for (f = 0; f < maxFields; f++) {
2498: if (numFields) {PetscSectionRestoreFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f]);}
2499: else {PetscSectionRestorePointSyms(localCoarse,numPoints,points,&perms[f],&flips[f]);}
2500: }
2501: DMRestoreWorkArray(coarse,numPoints,MPIU_SCALAR,&points);
2502: }
2503: }
2504: else if (matSize) {
2505: PetscInt cOff;
2506: PetscInt *rowIndices, *colIndices, a, aDof, aOff;
2508: numRowIndices = matSize / numColIndices;
2509: if (numRowIndices != dof) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Miscounted dofs");
2510: DMGetWorkArray(coarse,numRowIndices,MPIU_INT,&rowIndices);
2511: DMGetWorkArray(coarse,numColIndices,MPIU_INT,&colIndices);
2512: PetscSectionGetOffset(cSec,p,&cOff);
2513: PetscSectionGetDof(aSec,p,&aDof);
2514: PetscSectionGetOffset(aSec,p,&aOff);
2515: if (numFields) {
2516: for (f = 0; f < numFields; f++) {
2517: PetscInt fDof;
2519: PetscSectionGetFieldDof(cSec,p,f,&fDof);
2520: offsets[f + 1] = fDof;
2521: for (a = 0; a < aDof; a++) {
2522: PetscInt anchor = anchors[a + aOff];
2523: PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);
2524: newOffsets[f + 1] += fDof;
2525: }
2526: }
2527: for (f = 0; f < numFields; f++) {
2528: offsets[f + 1] += offsets[f];
2529: offsetsCopy[f + 1] = offsets[f + 1];
2530: newOffsets[f + 1] += newOffsets[f];
2531: newOffsetsCopy[f + 1] = newOffsets[f + 1];
2532: }
2533: DMPlexGetIndicesPointFields_Internal(cSec,p,cOff,offsetsCopy,PETSC_TRUE,NULL,-1,rowIndices);
2534: for (a = 0; a < aDof; a++) {
2535: PetscInt anchor = anchors[a + aOff], lOff;
2536: PetscSectionGetOffset(localCoarse,anchor,&lOff);
2537: DMPlexGetIndicesPointFields_Internal(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,-1,colIndices);
2538: }
2539: }
2540: else {
2541: DMPlexGetIndicesPoint_Internal(cSec,p,cOff,offsetsCopy,PETSC_TRUE,NULL,rowIndices);
2542: for (a = 0; a < aDof; a++) {
2543: PetscInt anchor = anchors[a + aOff], lOff;
2544: PetscSectionGetOffset(localCoarse,anchor,&lOff);
2545: DMPlexGetIndicesPoint_Internal(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,colIndices);
2546: }
2547: }
2548: if (numFields) {
2549: PetscInt count, a;
2551: for (f = 0, count = 0; f < numFields; f++) {
2552: PetscInt iSize = offsets[f + 1] - offsets[f];
2553: PetscInt jSize = newOffsets[f + 1] - newOffsets[f];
2554: MatGetValues(cMat,iSize,&rowIndices[offsets[f]],jSize,&colIndices[newOffsets[f]],&pMat[count]);
2555: count += iSize * jSize;
2556: pInd[numColIndices + f] = offsets[f+1];
2557: pInd[numColIndices + numFields + f] = newOffsets[f+1];
2558: }
2559: for (a = 0; a < aDof; a++) {
2560: PetscInt anchor = anchors[a + aOff];
2561: PetscInt gOff;
2562: PetscSectionGetOffset(globalCoarse,anchor,&gOff);
2563: DMPlexGetIndicesPointFields_Internal(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,-1,pInd);
2564: }
2565: }
2566: else {
2567: PetscInt a;
2568: MatGetValues(cMat,numRowIndices,rowIndices,numColIndices,colIndices,pMat);
2569: for (a = 0; a < aDof; a++) {
2570: PetscInt anchor = anchors[a + aOff];
2571: PetscInt gOff;
2572: PetscSectionGetOffset(globalCoarse,anchor,&gOff);
2573: DMPlexGetIndicesPoint_Internal(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,pInd);
2574: }
2575: }
2576: DMRestoreWorkArray(coarse,numColIndices,MPIU_INT,&colIndices);
2577: DMRestoreWorkArray(coarse,numRowIndices,MPIU_INT,&rowIndices);
2578: }
2579: else {
2580: PetscInt gOff;
2582: PetscSectionGetOffset(globalCoarse,p,&gOff);
2583: if (numFields) {
2584: for (f = 0; f < numFields; f++) {
2585: PetscInt fDof;
2586: PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
2587: offsets[f + 1] = fDof + offsets[f];
2588: }
2589: for (f = 0; f < numFields; f++) {
2590: pInd[numColIndices + f] = offsets[f+1];
2591: pInd[numColIndices + numFields + f] = offsets[f+1];
2592: }
2593: DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1,pInd);
2594: }
2595: else {
2596: DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,pInd);
2597: }
2598: }
2599: }
2600: PetscFree(maxChildIds);
2601: }
2602: {
2603: PetscSF indicesSF, matricesSF;
2604: PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices;
2606: PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);
2607: PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafMatricesSec);
2608: PetscSFDistributeSection(coarseToFineEmbedded,rootIndicesSec,&remoteOffsetsIndices,leafIndicesSec);
2609: PetscSFDistributeSection(coarseToFineEmbedded,rootMatricesSec,&remoteOffsetsMatrices,leafMatricesSec);
2610: PetscSFCreateSectionSF(coarseToFineEmbedded,rootIndicesSec,remoteOffsetsIndices,leafIndicesSec,&indicesSF);
2611: PetscSFCreateSectionSF(coarseToFineEmbedded,rootMatricesSec,remoteOffsetsMatrices,leafMatricesSec,&matricesSF);
2612: PetscSFDestroy(&coarseToFineEmbedded);
2613: PetscFree(remoteOffsetsIndices);
2614: PetscFree(remoteOffsetsMatrices);
2615: PetscSectionGetStorageSize(leafIndicesSec,&numLeafIndices);
2616: PetscSectionGetStorageSize(leafMatricesSec,&numLeafMatrices);
2617: PetscMalloc2(numLeafIndices,&leafIndices,numLeafMatrices,&leafMatrices);
2618: PetscSFBcastBegin(indicesSF,MPIU_INT,rootIndices,leafIndices);
2619: PetscSFBcastBegin(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);
2620: PetscSFBcastEnd(indicesSF,MPIU_INT,rootIndices,leafIndices);
2621: PetscSFBcastEnd(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);
2622: PetscSFDestroy(&matricesSF);
2623: PetscSFDestroy(&indicesSF);
2624: PetscFree2(rootIndices,rootMatrices);
2625: PetscSectionDestroy(&rootIndicesSec);
2626: PetscSectionDestroy(&rootMatricesSec);
2627: }
2628: /* count to preallocate */
2629: DMGetSection(fine,&localFine);
2630: {
2631: PetscInt nGlobal;
2632: PetscInt *dnnz, *onnz;
2633: PetscLayout rowMap, colMap;
2634: PetscInt rowStart, rowEnd, colStart, colEnd;
2635: PetscInt maxDof;
2636: PetscInt *rowIndices;
2637: DM refTree;
2638: PetscInt **refPointFieldN;
2639: PetscScalar ***refPointFieldMats;
2640: PetscSection refConSec, refAnSec;
2641: PetscInt pRefStart,pRefEnd,maxConDof,maxColumns,leafStart,leafEnd;
2642: PetscScalar *pointWork;
2644: PetscSectionGetConstrainedStorageSize(globalFine,&nGlobal);
2645: PetscCalloc2(nGlobal,&dnnz,nGlobal,&onnz);
2646: MatGetLayouts(mat,&rowMap,&colMap);
2647: PetscLayoutSetUp(rowMap);
2648: PetscLayoutSetUp(colMap);
2649: PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);
2650: PetscLayoutGetRange(colMap,&colStart,&colEnd);
2651: PetscSectionGetMaxDof(globalFine,&maxDof);
2652: PetscSectionGetChart(leafIndicesSec,&leafStart,&leafEnd);
2653: DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices);
2654: for (p = leafStart; p < leafEnd; p++) {
2655: PetscInt gDof, gcDof, gOff;
2656: PetscInt numColIndices, pIndOff, *pInd;
2657: PetscInt matSize;
2658: PetscInt i;
2660: PetscSectionGetDof(globalFine,p,&gDof);
2661: PetscSectionGetConstraintDof(globalFine,p,&gcDof);
2662: if ((gDof - gcDof) <= 0) {
2663: continue;
2664: }
2665: PetscSectionGetOffset(globalFine,p,&gOff);
2666: if (gOff < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I though having global dofs meant a non-negative offset");
2667: if ((gOff < rowStart) || ((gOff + gDof - gcDof) > rowEnd)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I thought the row map would constrain the global dofs");
2668: PetscSectionGetDof(leafIndicesSec,p,&numColIndices);
2669: PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);
2670: numColIndices -= 2 * numFields;
2671: if (numColIndices <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"global fine dof with no dofs to interpolate from");
2672: pInd = &leafIndices[pIndOff];
2673: offsets[0] = 0;
2674: offsetsCopy[0] = 0;
2675: newOffsets[0] = 0;
2676: newOffsetsCopy[0] = 0;
2677: if (numFields) {
2678: PetscInt f;
2679: for (f = 0; f < numFields; f++) {
2680: PetscInt rowDof;
2682: PetscSectionGetFieldDof(localFine,p,f,&rowDof);
2683: offsets[f + 1] = offsets[f] + rowDof;
2684: offsetsCopy[f + 1] = offsets[f + 1];
2685: newOffsets[f + 1] = pInd[numColIndices + numFields + f];
2686: numD[f] = 0;
2687: numO[f] = 0;
2688: }
2689: DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);
2690: for (f = 0; f < numFields; f++) {
2691: PetscInt colOffset = newOffsets[f];
2692: PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f];
2694: for (i = 0; i < numFieldCols; i++) {
2695: PetscInt gInd = pInd[i + colOffset];
2697: if (gInd >= colStart && gInd < colEnd) {
2698: numD[f]++;
2699: }
2700: else if (gInd >= 0) { /* negative means non-entry */
2701: numO[f]++;
2702: }
2703: }
2704: }
2705: }
2706: else {
2707: DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);
2708: numD[0] = 0;
2709: numO[0] = 0;
2710: for (i = 0; i < numColIndices; i++) {
2711: PetscInt gInd = pInd[i];
2713: if (gInd >= colStart && gInd < colEnd) {
2714: numD[0]++;
2715: }
2716: else if (gInd >= 0) { /* negative means non-entry */
2717: numO[0]++;
2718: }
2719: }
2720: }
2721: PetscSectionGetDof(leafMatricesSec,p,&matSize);
2722: if (!matSize) { /* incoming matrix is identity */
2723: PetscInt childId;
2725: childId = childIds[p-pStartF];
2726: if (childId < 0) { /* no child interpolation: one nnz per */
2727: if (numFields) {
2728: PetscInt f;
2729: for (f = 0; f < numFields; f++) {
2730: PetscInt numRows = offsets[f+1] - offsets[f], row;
2731: for (row = 0; row < numRows; row++) {
2732: PetscInt gIndCoarse = pInd[newOffsets[f] + row];
2733: PetscInt gIndFine = rowIndices[offsets[f] + row];
2734: if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2735: if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2736: dnnz[gIndFine - rowStart] = 1;
2737: }
2738: else if (gIndCoarse >= 0) { /* remote */
2739: if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2740: onnz[gIndFine - rowStart] = 1;
2741: }
2742: else { /* constrained */
2743: if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2744: }
2745: }
2746: }
2747: }
2748: else {
2749: PetscInt i;
2750: for (i = 0; i < gDof; i++) {
2751: PetscInt gIndCoarse = pInd[i];
2752: PetscInt gIndFine = rowIndices[i];
2753: if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2754: if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2755: dnnz[gIndFine - rowStart] = 1;
2756: }
2757: else if (gIndCoarse >= 0) { /* remote */
2758: if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2759: onnz[gIndFine - rowStart] = 1;
2760: }
2761: else { /* constrained */
2762: if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2763: }
2764: }
2765: }
2766: }
2767: else { /* interpolate from all */
2768: if (numFields) {
2769: PetscInt f;
2770: for (f = 0; f < numFields; f++) {
2771: PetscInt numRows = offsets[f+1] - offsets[f], row;
2772: for (row = 0; row < numRows; row++) {
2773: PetscInt gIndFine = rowIndices[offsets[f] + row];
2774: if (gIndFine >= 0) {
2775: if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2776: dnnz[gIndFine - rowStart] = numD[f];
2777: onnz[gIndFine - rowStart] = numO[f];
2778: }
2779: }
2780: }
2781: }
2782: else {
2783: PetscInt i;
2784: for (i = 0; i < gDof; i++) {
2785: PetscInt gIndFine = rowIndices[i];
2786: if (gIndFine >= 0) {
2787: if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2788: dnnz[gIndFine - rowStart] = numD[0];
2789: onnz[gIndFine - rowStart] = numO[0];
2790: }
2791: }
2792: }
2793: }
2794: }
2795: else { /* interpolate from all */
2796: if (numFields) {
2797: PetscInt f;
2798: for (f = 0; f < numFields; f++) {
2799: PetscInt numRows = offsets[f+1] - offsets[f], row;
2800: for (row = 0; row < numRows; row++) {
2801: PetscInt gIndFine = rowIndices[offsets[f] + row];
2802: if (gIndFine >= 0) {
2803: if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2804: dnnz[gIndFine - rowStart] = numD[f];
2805: onnz[gIndFine - rowStart] = numO[f];
2806: }
2807: }
2808: }
2809: }
2810: else { /* every dof get a full row */
2811: PetscInt i;
2812: for (i = 0; i < gDof; i++) {
2813: PetscInt gIndFine = rowIndices[i];
2814: if (gIndFine >= 0) {
2815: if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2816: dnnz[gIndFine - rowStart] = numD[0];
2817: onnz[gIndFine - rowStart] = numO[0];
2818: }
2819: }
2820: }
2821: }
2822: }
2823: MatXAIJSetPreallocation(mat,1,dnnz,onnz,NULL,NULL);
2824: PetscFree2(dnnz,onnz);
2826: DMPlexGetReferenceTree(fine,&refTree);
2827: DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
2828: DMGetDefaultConstraints(refTree,&refConSec,NULL);
2829: DMPlexGetAnchors(refTree,&refAnSec,NULL);
2830: PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
2831: PetscSectionGetMaxDof(refConSec,&maxConDof);
2832: PetscSectionGetMaxDof(leafIndicesSec,&maxColumns);
2833: PetscMalloc1(maxConDof*maxColumns,&pointWork);
2834: for (p = leafStart; p < leafEnd; p++) {
2835: PetscInt gDof, gcDof, gOff;
2836: PetscInt numColIndices, pIndOff, *pInd;
2837: PetscInt matSize;
2838: PetscInt childId;
2841: PetscSectionGetDof(globalFine,p,&gDof);
2842: PetscSectionGetConstraintDof(globalFine,p,&gcDof);
2843: if ((gDof - gcDof) <= 0) {
2844: continue;
2845: }
2846: childId = childIds[p-pStartF];
2847: PetscSectionGetOffset(globalFine,p,&gOff);
2848: PetscSectionGetDof(leafIndicesSec,p,&numColIndices);
2849: PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);
2850: numColIndices -= 2 * numFields;
2851: pInd = &leafIndices[pIndOff];
2852: offsets[0] = 0;
2853: offsetsCopy[0] = 0;
2854: newOffsets[0] = 0;
2855: newOffsetsCopy[0] = 0;
2856: rowOffsets[0] = 0;
2857: if (numFields) {
2858: PetscInt f;
2859: for (f = 0; f < numFields; f++) {
2860: PetscInt rowDof;
2862: PetscSectionGetFieldDof(localFine,p,f,&rowDof);
2863: offsets[f + 1] = offsets[f] + rowDof;
2864: offsetsCopy[f + 1] = offsets[f + 1];
2865: rowOffsets[f + 1] = pInd[numColIndices + f];
2866: newOffsets[f + 1] = pInd[numColIndices + numFields + f];
2867: }
2868: DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);
2869: }
2870: else {
2871: DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);
2872: }
2873: PetscSectionGetDof(leafMatricesSec,p,&matSize);
2874: if (!matSize) { /* incoming matrix is identity */
2875: if (childId < 0) { /* no child interpolation: scatter */
2876: if (numFields) {
2877: PetscInt f;
2878: for (f = 0; f < numFields; f++) {
2879: PetscInt numRows = offsets[f+1] - offsets[f], row;
2880: for (row = 0; row < numRows; row++) {
2881: MatSetValue(mat,rowIndices[offsets[f]+row],pInd[newOffsets[f]+row],1.,INSERT_VALUES);
2882: }
2883: }
2884: }
2885: else {
2886: PetscInt numRows = gDof, row;
2887: for (row = 0; row < numRows; row++) {
2888: MatSetValue(mat,rowIndices[row],pInd[row],1.,INSERT_VALUES);
2889: }
2890: }
2891: }
2892: else { /* interpolate from all */
2893: if (numFields) {
2894: PetscInt f;
2895: for (f = 0; f < numFields; f++) {
2896: PetscInt numRows = offsets[f+1] - offsets[f];
2897: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2898: MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],refPointFieldMats[childId - pRefStart][f],INSERT_VALUES);
2899: }
2900: }
2901: else {
2902: MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,refPointFieldMats[childId - pRefStart][0],INSERT_VALUES);
2903: }
2904: }
2905: }
2906: else { /* interpolate from all */
2907: PetscInt pMatOff;
2908: PetscScalar *pMat;
2910: PetscSectionGetOffset(leafMatricesSec,p,&pMatOff);
2911: pMat = &leafMatrices[pMatOff];
2912: if (childId < 0) { /* copy the incoming matrix */
2913: if (numFields) {
2914: PetscInt f, count;
2915: for (f = 0, count = 0; f < numFields; f++) {
2916: PetscInt numRows = offsets[f+1]-offsets[f];
2917: PetscInt numCols = newOffsets[f+1]-newOffsets[f];
2918: PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f];
2919: PetscScalar *inMat = &pMat[count];
2921: MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],inMat,INSERT_VALUES);
2922: count += numCols * numInRows;
2923: }
2924: }
2925: else {
2926: MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,pMat,INSERT_VALUES);
2927: }
2928: }
2929: else { /* multiply the incoming matrix by the child interpolation */
2930: if (numFields) {
2931: PetscInt f, count;
2932: for (f = 0, count = 0; f < numFields; f++) {
2933: PetscInt numRows = offsets[f+1]-offsets[f];
2934: PetscInt numCols = newOffsets[f+1]-newOffsets[f];
2935: PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f];
2936: PetscScalar *inMat = &pMat[count];
2937: PetscInt i, j, k;
2938: if (refPointFieldN[childId - pRefStart][f] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch");
2939: for (i = 0; i < numRows; i++) {
2940: for (j = 0; j < numCols; j++) {
2941: PetscScalar val = 0.;
2942: for (k = 0; k < numInRows; k++) {
2943: val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j];
2944: }
2945: pointWork[i * numCols + j] = val;
2946: }
2947: }
2948: MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],pointWork,INSERT_VALUES);
2949: count += numCols * numInRows;
2950: }
2951: }
2952: else { /* every dof gets a full row */
2953: PetscInt numRows = gDof;
2954: PetscInt numCols = numColIndices;
2955: PetscInt numInRows = matSize / numColIndices;
2956: PetscInt i, j, k;
2957: if (refPointFieldN[childId - pRefStart][0] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch");
2958: for (i = 0; i < numRows; i++) {
2959: for (j = 0; j < numCols; j++) {
2960: PetscScalar val = 0.;
2961: for (k = 0; k < numInRows; k++) {
2962: val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j];
2963: }
2964: pointWork[i * numCols + j] = val;
2965: }
2966: }
2967: MatSetValues(mat,numRows,rowIndices,numCols,pInd,pointWork,INSERT_VALUES);
2968: }
2969: }
2970: }
2971: }
2972: DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
2973: DMRestoreWorkArray(fine,maxDof,MPIU_INT,&rowIndices);
2974: PetscFree(pointWork);
2975: }
2976: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
2977: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
2978: PetscSectionDestroy(&leafIndicesSec);
2979: PetscSectionDestroy(&leafMatricesSec);
2980: PetscFree2(leafIndices,leafMatrices);
2981: PetscFree2(*(PetscInt****)&perms,*(PetscScalar****)&flips);
2982: PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);
2983: ISRestoreIndices(aIS,&anchors);
2984: return(0);
2985: }
2987: /*
2988: * Assuming a nodal basis (w.r.t. the dual basis) basis:
2989: *
2990: * for each coarse dof \phi^c_i:
2991: * for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i:
2992: * for each fine dof \phi^f_j;
2993: * a_{i,j} = 0;
2994: * for each fine dof \phi^f_k:
2995: * a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l
2996: * [^^^ this is = \phi^c_i ^^^]
2997: */
2998: PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj)
2999: {
3000: PetscDS ds;
3001: PetscSection section, cSection;
3002: DMLabel canonical, depth;
3003: Mat cMat, mat;
3004: PetscInt *nnz;
3005: PetscInt f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd;
3006: PetscInt m, n;
3007: PetscScalar *pointScalar;
3008: PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent;
3012: DMGetSection(refTree,§ion);
3013: DMGetDimension(refTree, &dim);
3014: PetscMalloc6(dim,&v0,dim,&v0parent,dim,&vtmp,dim*dim,&J,dim*dim,&Jparent,dim*dim,&invJ);
3015: PetscMalloc2(dim,&pointScalar,dim,&pointRef);
3016: DMGetDS(refTree,&ds);
3017: PetscDSGetNumFields(ds,&numFields);
3018: PetscSectionGetNumFields(section,&numSecFields);
3019: DMGetLabel(refTree,"canonical",&canonical);
3020: DMGetLabel(refTree,"depth",&depth);
3021: DMGetDefaultConstraints(refTree,&cSection,&cMat);
3022: DMPlexGetChart(refTree, &pStart, &pEnd);
3023: DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd);
3024: MatGetSize(cMat,&n,&m); /* the injector has transpose sizes from the constraint matrix */
3025: /* Step 1: compute non-zero pattern. A proper subset of constraint matrix non-zero */
3026: PetscCalloc1(m,&nnz);
3027: for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */
3028: const PetscInt *children;
3029: PetscInt numChildren;
3030: PetscInt i, numChildDof, numSelfDof;
3032: if (canonical) {
3033: PetscInt pCanonical;
3034: DMLabelGetValue(canonical,p,&pCanonical);
3035: if (p != pCanonical) continue;
3036: }
3037: DMPlexGetTreeChildren(refTree,p,&numChildren,&children);
3038: if (!numChildren) continue;
3039: for (i = 0, numChildDof = 0; i < numChildren; i++) {
3040: PetscInt child = children[i];
3041: PetscInt dof;
3043: PetscSectionGetDof(section,child,&dof);
3044: numChildDof += dof;
3045: }
3046: PetscSectionGetDof(section,p,&numSelfDof);
3047: if (!numChildDof || !numSelfDof) continue;
3048: for (f = 0; f < numFields; f++) {
3049: PetscInt selfOff;
3051: if (numSecFields) { /* count the dofs for just this field */
3052: for (i = 0, numChildDof = 0; i < numChildren; i++) {
3053: PetscInt child = children[i];
3054: PetscInt dof;
3056: PetscSectionGetFieldDof(section,child,f,&dof);
3057: numChildDof += dof;
3058: }
3059: PetscSectionGetFieldDof(section,p,f,&numSelfDof);
3060: PetscSectionGetFieldOffset(section,p,f,&selfOff);
3061: }
3062: else {
3063: PetscSectionGetOffset(section,p,&selfOff);
3064: }
3065: for (i = 0; i < numSelfDof; i++) {
3066: nnz[selfOff + i] = numChildDof;
3067: }
3068: }
3069: }
3070: MatCreateAIJ(PETSC_COMM_SELF,m,n,m,n,-1,nnz,-1,NULL,&mat);
3071: PetscFree(nnz);
3072: /* Setp 2: compute entries */
3073: for (p = pStart; p < pEnd; p++) {
3074: const PetscInt *children;
3075: PetscInt numChildren;
3076: PetscInt i, numChildDof, numSelfDof;
3078: /* same conditions about when entries occur */
3079: if (canonical) {
3080: PetscInt pCanonical;
3081: DMLabelGetValue(canonical,p,&pCanonical);
3082: if (p != pCanonical) continue;
3083: }
3084: DMPlexGetTreeChildren(refTree,p,&numChildren,&children);
3085: if (!numChildren) continue;
3086: for (i = 0, numChildDof = 0; i < numChildren; i++) {
3087: PetscInt child = children[i];
3088: PetscInt dof;
3090: PetscSectionGetDof(section,child,&dof);
3091: numChildDof += dof;
3092: }
3093: PetscSectionGetDof(section,p,&numSelfDof);
3094: if (!numChildDof || !numSelfDof) continue;
3096: for (f = 0; f < numFields; f++) {
3097: PetscInt selfOff, Nc, parentCell;
3098: PetscInt cellShapeOff;
3099: PetscObject disc;
3100: PetscDualSpace dsp;
3101: PetscClassId classId;
3102: PetscScalar *pointMat;
3103: PetscInt *matRows, *matCols;
3104: PetscInt pO = PETSC_MIN_INT;
3105: const PetscInt *depthNumDof;
3107: if (numSecFields) {
3108: for (i = 0, numChildDof = 0; i < numChildren; i++) {
3109: PetscInt child = children[i];
3110: PetscInt dof;
3112: PetscSectionGetFieldDof(section,child,f,&dof);
3113: numChildDof += dof;
3114: }
3115: PetscSectionGetFieldDof(section,p,f,&numSelfDof);
3116: PetscSectionGetFieldOffset(section,p,f,&selfOff);
3117: }
3118: else {
3119: PetscSectionGetOffset(section,p,&selfOff);
3120: }
3122: /* find a cell whose closure contains p */
3123: if (p >= cStart && p < cEnd) {
3124: parentCell = p;
3125: }
3126: else {
3127: PetscInt *star = NULL;
3128: PetscInt numStar;
3130: parentCell = -1;
3131: DMPlexGetTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);
3132: for (i = numStar - 1; i >= 0; i--) {
3133: PetscInt c = star[2 * i];
3135: if (c >= cStart && c < cEnd) {
3136: parentCell = c;
3137: break;
3138: }
3139: }
3140: DMPlexRestoreTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);
3141: }
3142: /* determine the offset of p's shape functions withing parentCell's shape functions */
3143: PetscDSGetDiscretization(ds,f,&disc);
3144: PetscObjectGetClassId(disc,&classId);
3145: if (classId == PETSCFE_CLASSID) {
3146: PetscFEGetDualSpace((PetscFE)disc,&dsp);
3147: }
3148: else if (classId == PETSCFV_CLASSID) {
3149: PetscFVGetDualSpace((PetscFV)disc,&dsp);
3150: }
3151: else {
3152: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported discretization object");
3153: }
3154: PetscDualSpaceGetNumDof(dsp,&depthNumDof);
3155: PetscDualSpaceGetNumComponents(dsp,&Nc);
3156: {
3157: PetscInt *closure = NULL;
3158: PetscInt numClosure;
3160: DMPlexGetTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);
3161: for (i = 0, cellShapeOff = 0; i < numClosure; i++) {
3162: PetscInt point = closure[2 * i], pointDepth;
3164: pO = closure[2 * i + 1];
3165: if (point == p) break;
3166: DMLabelGetValue(depth,point,&pointDepth);
3167: cellShapeOff += depthNumDof[pointDepth];
3168: }
3169: DMPlexRestoreTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);
3170: }
3172: DMGetWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat);
3173: DMGetWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows);
3174: matCols = matRows + numSelfDof;
3175: for (i = 0; i < numSelfDof; i++) {
3176: matRows[i] = selfOff + i;
3177: }
3178: for (i = 0; i < numSelfDof * numChildDof; i++) pointMat[i] = 0.;
3179: {
3180: PetscInt colOff = 0;
3182: for (i = 0; i < numChildren; i++) {
3183: PetscInt child = children[i];
3184: PetscInt dof, off, j;
3186: if (numSecFields) {
3187: PetscSectionGetFieldDof(cSection,child,f,&dof);
3188: PetscSectionGetFieldOffset(cSection,child,f,&off);
3189: }
3190: else {
3191: PetscSectionGetDof(cSection,child,&dof);
3192: PetscSectionGetOffset(cSection,child,&off);
3193: }
3195: for (j = 0; j < dof; j++) {
3196: matCols[colOff++] = off + j;
3197: }
3198: }
3199: }
3200: if (classId == PETSCFE_CLASSID) {
3201: PetscFE fe = (PetscFE) disc;
3202: PetscInt fSize;
3204: PetscFEGetDualSpace(fe,&dsp);
3205: PetscDualSpaceGetDimension(dsp,&fSize);
3206: for (i = 0; i < numSelfDof; i++) { /* for every shape function */
3207: PetscQuadrature q;
3208: PetscInt dim, thisNc, numPoints, j, k;
3209: const PetscReal *points;
3210: const PetscReal *weights;
3211: PetscInt *closure = NULL;
3212: PetscInt numClosure;
3213: PetscInt parentCellShapeDof = cellShapeOff + (pO < 0 ? (numSelfDof - 1 - i) : i);
3214: PetscReal *Bparent;
3216: PetscDualSpaceGetFunctional(dsp,parentCellShapeDof,&q);
3217: PetscQuadratureGetData(q,&dim,&thisNc,&numPoints,&points,&weights);
3218: if (thisNc != Nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Functional dim %D does not much basis dim %D\n",thisNc,Nc);
3219: PetscFEGetTabulation(fe,numPoints,points,&Bparent,NULL,NULL); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */
3220: for (j = 0; j < numPoints; j++) {
3221: PetscInt childCell = -1;
3222: PetscReal *parentValAtPoint;
3223: const PetscReal xi0[3] = {-1.,-1.,-1.};
3224: const PetscReal *pointReal = &points[dim * j];
3225: const PetscScalar *point;
3226: PetscReal *Bchild;
3227: PetscInt childCellShapeOff, pointMatOff;
3228: #if defined(PETSC_USE_COMPLEX)
3229: PetscInt d;
3231: for (d = 0; d < dim; d++) {
3232: pointScalar[d] = points[dim * j + d];
3233: }
3234: point = pointScalar;
3235: #else
3236: point = pointReal;
3237: #endif
3239: parentValAtPoint = &Bparent[(fSize * j + parentCellShapeDof) * Nc];
3241: for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/
3242: PetscInt child = children[k];
3243: PetscInt *star = NULL;
3244: PetscInt numStar, s;
3246: DMPlexGetTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);
3247: for (s = numStar - 1; s >= 0; s--) {
3248: PetscInt c = star[2 * s];
3250: if (c < cStart || c >= cEnd) continue;
3251: DMPlexLocatePoint_Internal(refTree,dim,point,c,&childCell);
3252: if (childCell >= 0) break;
3253: }
3254: DMPlexRestoreTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);
3255: if (childCell >= 0) break;
3256: }
3257: if (childCell < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Could not locate quadrature point");
3258: DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ);
3259: DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent);
3260: CoordinatesRefToReal(dim, dim, xi0, v0parent, Jparent, pointReal, vtmp);
3261: CoordinatesRealToRef(dim, dim, xi0, v0, invJ, vtmp, pointRef);
3263: PetscFEGetTabulation(fe,1,pointRef,&Bchild,NULL,NULL);
3264: DMPlexGetTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);
3265: for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */
3266: PetscInt child = children[k], childDepth, childDof, childO = PETSC_MIN_INT;
3267: PetscInt l;
3269: DMLabelGetValue(depth,child,&childDepth);
3270: childDof = depthNumDof[childDepth];
3271: for (l = 0, childCellShapeOff = 0; l < numClosure; l++) {
3272: PetscInt point = closure[2 * l];
3273: PetscInt pointDepth;
3275: childO = closure[2 * l + 1];
3276: if (point == child) break;
3277: DMLabelGetValue(depth,point,&pointDepth);
3278: childCellShapeOff += depthNumDof[pointDepth];
3279: }
3280: if (l == numClosure) {
3281: pointMatOff += childDof;
3282: continue; /* child is not in the closure of the cell: has nothing to contribute to this point */
3283: }
3284: for (l = 0; l < childDof; l++) {
3285: PetscInt childCellDof = childCellShapeOff + (childO ? (childDof - 1 - l) : l), m;
3286: PetscReal *childValAtPoint;
3287: PetscReal val = 0.;
3289: childValAtPoint = &Bchild[childCellDof * Nc];
3290: for (m = 0; m < Nc; m++) {
3291: val += weights[j * Nc + m] * parentValAtPoint[m] * childValAtPoint[m];
3292: }
3294: pointMat[i * numChildDof + pointMatOff + l] += val;
3295: }
3296: pointMatOff += childDof;
3297: }
3298: DMPlexRestoreTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);
3299: PetscFERestoreTabulation(fe,1,pointRef,&Bchild,NULL,NULL);
3300: }
3301: PetscFERestoreTabulation(fe,numPoints,points,&Bparent,NULL,NULL);
3302: }
3303: }
3304: else { /* just the volume-weighted averages of the children */
3305: PetscReal parentVol;
3306: PetscInt childCell;
3308: DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL);
3309: for (i = 0, childCell = 0; i < numChildren; i++) {
3310: PetscInt child = children[i], j;
3311: PetscReal childVol;
3313: if (child < cStart || child >= cEnd) continue;
3314: DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL);
3315: for (j = 0; j < Nc; j++) {
3316: pointMat[j * numChildDof + Nc * childCell + j] = childVol / parentVol;
3317: }
3318: childCell++;
3319: }
3320: }
3321: /* Insert pointMat into mat */
3322: MatSetValues(mat,numSelfDof,matRows,numChildDof,matCols,pointMat,INSERT_VALUES);
3323: DMRestoreWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows);
3324: DMRestoreWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat);
3325: }
3326: }
3327: PetscFree6(v0,v0parent,vtmp,J,Jparent,invJ);
3328: PetscFree2(pointScalar,pointRef);
3329: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
3330: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
3331: *inj = mat;
3332: return(0);
3333: }
3335: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3336: {
3337: PetscDS ds;
3338: PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof;
3339: PetscScalar ***refPointFieldMats;
3340: PetscSection refConSec, refSection;
3344: DMGetDS(refTree,&ds);
3345: PetscDSGetNumFields(ds,&numFields);
3346: DMGetDefaultConstraints(refTree,&refConSec,NULL);
3347: DMGetSection(refTree,&refSection);
3348: PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
3349: PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);
3350: PetscSectionGetMaxDof(refConSec,&maxDof);
3351: PetscMalloc1(maxDof,&rows);
3352: PetscMalloc1(maxDof*maxDof,&cols);
3353: for (p = pRefStart; p < pRefEnd; p++) {
3354: PetscInt parent, pDof, parentDof;
3356: DMPlexGetTreeParent(refTree,p,&parent,NULL);
3357: PetscSectionGetDof(refConSec,p,&pDof);
3358: PetscSectionGetDof(refSection,parent,&parentDof);
3359: if (!pDof || !parentDof || parent == p) continue;
3361: PetscMalloc1(numFields,&refPointFieldMats[p-pRefStart]);
3362: for (f = 0; f < numFields; f++) {
3363: PetscInt cDof, cOff, numCols, r;
3365: if (numFields > 1) {
3366: PetscSectionGetFieldDof(refConSec,p,f,&cDof);
3367: PetscSectionGetFieldOffset(refConSec,p,f,&cOff);
3368: }
3369: else {
3370: PetscSectionGetDof(refConSec,p,&cDof);
3371: PetscSectionGetOffset(refConSec,p,&cOff);
3372: }
3374: for (r = 0; r < cDof; r++) {
3375: rows[r] = cOff + r;
3376: }
3377: numCols = 0;
3378: {
3379: PetscInt aDof, aOff, j;
3381: if (numFields > 1) {
3382: PetscSectionGetFieldDof(refSection,parent,f,&aDof);
3383: PetscSectionGetFieldOffset(refSection,parent,f,&aOff);
3384: }
3385: else {
3386: PetscSectionGetDof(refSection,parent,&aDof);
3387: PetscSectionGetOffset(refSection,parent,&aOff);
3388: }
3390: for (j = 0; j < aDof; j++) {
3391: cols[numCols++] = aOff + j;
3392: }
3393: }
3394: PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);
3395: /* transpose of constraint matrix */
3396: MatGetValues(inj,numCols,cols,cDof,rows,refPointFieldMats[p-pRefStart][f]);
3397: }
3398: }
3399: *childrenMats = refPointFieldMats;
3400: PetscFree(rows);
3401: PetscFree(cols);
3402: return(0);
3403: }
3405: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3406: {
3407: PetscDS ds;
3408: PetscScalar ***refPointFieldMats;
3409: PetscInt numFields, pRefStart, pRefEnd, p, f;
3410: PetscSection refConSec, refSection;
3414: refPointFieldMats = *childrenMats;
3415: *childrenMats = NULL;
3416: DMGetDS(refTree,&ds);
3417: DMGetSection(refTree,&refSection);
3418: PetscDSGetNumFields(ds,&numFields);
3419: DMGetDefaultConstraints(refTree,&refConSec,NULL);
3420: PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
3421: for (p = pRefStart; p < pRefEnd; p++) {
3422: PetscInt parent, pDof, parentDof;
3424: DMPlexGetTreeParent(refTree,p,&parent,NULL);
3425: PetscSectionGetDof(refConSec,p,&pDof);
3426: PetscSectionGetDof(refSection,parent,&parentDof);
3427: if (!pDof || !parentDof || parent == p) continue;
3429: for (f = 0; f < numFields; f++) {
3430: PetscInt cDof;
3432: if (numFields > 1) {
3433: PetscSectionGetFieldDof(refConSec,p,f,&cDof);
3434: }
3435: else {
3436: PetscSectionGetDof(refConSec,p,&cDof);
3437: }
3439: PetscFree(refPointFieldMats[p - pRefStart][f]);
3440: }
3441: PetscFree(refPointFieldMats[p - pRefStart]);
3442: }
3443: PetscFree(refPointFieldMats);
3444: return(0);
3445: }
3447: static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree,Mat *injRef)
3448: {
3449: Mat cMatRef;
3450: PetscObject injRefObj;
3454: DMGetDefaultConstraints(refTree,NULL,&cMatRef);
3455: PetscObjectQuery((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",&injRefObj);
3456: *injRef = (Mat) injRefObj;
3457: if (!*injRef) {
3458: DMPlexComputeInjectorReferenceTree(refTree,injRef);
3459: PetscObjectCompose((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",(PetscObject)*injRef);
3460: /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */
3461: PetscObjectDereference((PetscObject)*injRef);
3462: }
3463: return(0);
3464: }
3466: 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)
3467: {
3468: PetscInt pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti;
3469: PetscSection globalCoarse, globalFine;
3470: PetscSection localCoarse, localFine, leafIndicesSec;
3471: PetscSection multiRootSec, rootIndicesSec;
3472: PetscInt *leafInds, *rootInds = NULL;
3473: const PetscInt *rootDegrees;
3474: PetscScalar *leafVals = NULL, *rootVals = NULL;
3475: PetscSF coarseToFineEmbedded;
3479: DMPlexGetChart(coarse,&pStartC,&pEndC);
3480: DMPlexGetChart(fine,&pStartF,&pEndF);
3481: DMGetSection(fine,&localFine);
3482: DMGetGlobalSection(fine,&globalFine);
3483: PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);
3484: PetscSectionSetChart(leafIndicesSec,pStartF, pEndF);
3485: PetscSectionGetMaxDof(localFine,&maxDof);
3486: { /* winnow fine points that don't have global dofs out of the sf */
3487: PetscInt l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices;
3488: const PetscInt *leaves;
3490: PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);
3491: for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
3492: p = leaves ? leaves[l] : l;
3493: PetscSectionGetDof(globalFine,p,&dof);
3494: PetscSectionGetConstraintDof(globalFine,p,&cdof);
3495: if ((dof - cdof) > 0) {
3496: numPointsWithDofs++;
3498: PetscSectionGetDof(localFine,p,&dof);
3499: PetscSectionSetDof(leafIndicesSec,p,dof + 1);
3500: }
3501: }
3502: PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
3503: PetscSectionSetUp(leafIndicesSec);
3504: PetscSectionGetStorageSize(leafIndicesSec,&numIndices);
3505: PetscMalloc1(gatheredIndices ? numIndices : (maxDof + 1),&leafInds);
3506: if (gatheredValues) {PetscMalloc1(numIndices,&leafVals);}
3507: for (l = 0, offset = 0; l < nleaves; l++) {
3508: p = leaves ? leaves[l] : l;
3509: PetscSectionGetDof(globalFine,p,&dof);
3510: PetscSectionGetConstraintDof(globalFine,p,&cdof);
3511: if ((dof - cdof) > 0) {
3512: PetscInt off, gOff;
3513: PetscInt *pInd;
3514: PetscScalar *pVal = NULL;
3516: pointsWithDofs[offset++] = l;
3518: PetscSectionGetOffset(leafIndicesSec,p,&off);
3520: pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds;
3521: if (gatheredValues) {
3522: PetscInt i;
3524: pVal = &leafVals[off + 1];
3525: for (i = 0; i < dof; i++) pVal[i] = 0.;
3526: }
3527: PetscSectionGetOffset(globalFine,p,&gOff);
3529: offsets[0] = 0;
3530: if (numFields) {
3531: PetscInt f;
3533: for (f = 0; f < numFields; f++) {
3534: PetscInt fDof;
3535: PetscSectionGetFieldDof(localFine,p,f,&fDof);
3536: offsets[f + 1] = fDof + offsets[f];
3537: }
3538: DMPlexGetIndicesPointFields_Internal(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1,pInd);
3539: }
3540: else {
3541: DMPlexGetIndicesPoint_Internal(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,pInd);
3542: }
3543: if (gatheredValues) {VecGetValues(fineVec,dof,pInd,pVal);}
3544: }
3545: }
3546: PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
3547: PetscFree(pointsWithDofs);
3548: }
3550: DMPlexGetChart(coarse,&pStartC,&pEndC);
3551: DMGetSection(coarse,&localCoarse);
3552: DMGetGlobalSection(coarse,&globalCoarse);
3554: { /* there may be the case where an sf root has a parent: broadcast parents back to children */
3555: MPI_Datatype threeInt;
3556: PetscMPIInt rank;
3557: PetscInt (*parentNodeAndIdCoarse)[3];
3558: PetscInt (*parentNodeAndIdFine)[3];
3559: PetscInt p, nleaves, nleavesToParents;
3560: PetscSF pointSF, sfToParents;
3561: const PetscInt *ilocal;
3562: const PetscSFNode *iremote;
3563: PetscSFNode *iremoteToParents;
3564: PetscInt *ilocalToParents;
3566: MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);
3567: MPI_Type_contiguous(3,MPIU_INT,&threeInt);
3568: MPI_Type_commit(&threeInt);
3569: PetscMalloc2(pEndC-pStartC,&parentNodeAndIdCoarse,pEndF-pStartF,&parentNodeAndIdFine);
3570: DMGetPointSF(coarse,&pointSF);
3571: PetscSFGetGraph(pointSF,NULL,&nleaves,&ilocal,&iremote);
3572: for (p = pStartC; p < pEndC; p++) {
3573: PetscInt parent, childId;
3574: DMPlexGetTreeParent(coarse,p,&parent,&childId);
3575: parentNodeAndIdCoarse[p - pStartC][0] = rank;
3576: parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC;
3577: parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId;
3578: if (nleaves > 0) {
3579: PetscInt leaf = -1;
3581: if (ilocal) {
3582: PetscFindInt(parent,nleaves,ilocal,&leaf);
3583: }
3584: else {
3585: leaf = p - pStartC;
3586: }
3587: if (leaf >= 0) {
3588: parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank;
3589: parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index;
3590: }
3591: }
3592: }
3593: for (p = pStartF; p < pEndF; p++) {
3594: parentNodeAndIdFine[p - pStartF][0] = -1;
3595: parentNodeAndIdFine[p - pStartF][1] = -1;
3596: parentNodeAndIdFine[p - pStartF][2] = -1;
3597: }
3598: PetscSFBcastBegin(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);
3599: PetscSFBcastEnd(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);
3600: for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3601: PetscInt dof;
3603: PetscSectionGetDof(leafIndicesSec,p,&dof);
3604: if (dof) {
3605: PetscInt off;
3607: PetscSectionGetOffset(leafIndicesSec,p,&off);
3608: if (gatheredIndices) {
3609: leafInds[off] = PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]);
3610: } else if (gatheredValues) {
3611: leafVals[off] = (PetscScalar) PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]);
3612: }
3613: }
3614: if (parentNodeAndIdFine[p-pStartF][0] >= 0) {
3615: nleavesToParents++;
3616: }
3617: }
3618: PetscMalloc1(nleavesToParents,&ilocalToParents);
3619: PetscMalloc1(nleavesToParents,&iremoteToParents);
3620: for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3621: if (parentNodeAndIdFine[p-pStartF][0] >= 0) {
3622: ilocalToParents[nleavesToParents] = p - pStartF;
3623: iremoteToParents[nleavesToParents].rank = parentNodeAndIdFine[p-pStartF][0];
3624: iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p-pStartF][1];
3625: nleavesToParents++;
3626: }
3627: }
3628: PetscSFCreate(PetscObjectComm((PetscObject)coarse),&sfToParents);
3629: PetscSFSetGraph(sfToParents,pEndC-pStartC,nleavesToParents,ilocalToParents,PETSC_OWN_POINTER,iremoteToParents,PETSC_OWN_POINTER);
3630: PetscSFDestroy(&coarseToFineEmbedded);
3632: coarseToFineEmbedded = sfToParents;
3634: PetscFree2(parentNodeAndIdCoarse,parentNodeAndIdFine);
3635: MPI_Type_free(&threeInt);
3636: }
3638: { /* winnow out coarse points that don't have dofs */
3639: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;
3640: PetscSF sfDofsOnly;
3642: for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) {
3643: PetscSectionGetDof(globalCoarse,p,&dof);
3644: PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3645: if ((dof - cdof) > 0) {
3646: numPointsWithDofs++;
3647: }
3648: }
3649: PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
3650: for (p = pStartC, offset = 0; p < pEndC; p++) {
3651: PetscSectionGetDof(globalCoarse,p,&dof);
3652: PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3653: if ((dof - cdof) > 0) {
3654: pointsWithDofs[offset++] = p - pStartC;
3655: }
3656: }
3657: PetscSFCreateEmbeddedSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly);
3658: PetscSFDestroy(&coarseToFineEmbedded);
3659: PetscFree(pointsWithDofs);
3660: coarseToFineEmbedded = sfDofsOnly;
3661: }
3663: /* communicate back to the coarse mesh which coarse points have children (that may require injection) */
3664: PetscSFComputeDegreeBegin(coarseToFineEmbedded,&rootDegrees);
3665: PetscSFComputeDegreeEnd(coarseToFineEmbedded,&rootDegrees);
3666: PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&multiRootSec);
3667: PetscSectionSetChart(multiRootSec,pStartC,pEndC);
3668: for (p = pStartC; p < pEndC; p++) {
3669: PetscSectionSetDof(multiRootSec,p,rootDegrees[p-pStartC]);
3670: }
3671: PetscSectionSetUp(multiRootSec);
3672: PetscSectionGetStorageSize(multiRootSec,&numMulti);
3673: PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);
3674: { /* distribute the leaf section */
3675: PetscSF multi, multiInv, indicesSF;
3676: PetscInt *remoteOffsets, numRootIndices;
3678: PetscSFGetMultiSF(coarseToFineEmbedded,&multi);
3679: PetscSFCreateInverseSF(multi,&multiInv);
3680: PetscSFDistributeSection(multiInv,leafIndicesSec,&remoteOffsets,rootIndicesSec);
3681: PetscSFCreateSectionSF(multiInv,leafIndicesSec,remoteOffsets,rootIndicesSec,&indicesSF);
3682: PetscFree(remoteOffsets);
3683: PetscSFDestroy(&multiInv);
3684: PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);
3685: if (gatheredIndices) {
3686: PetscMalloc1(numRootIndices,&rootInds);
3687: PetscSFBcastBegin(indicesSF,MPIU_INT,leafInds,rootInds);
3688: PetscSFBcastEnd(indicesSF,MPIU_INT,leafInds,rootInds);
3689: }
3690: if (gatheredValues) {
3691: PetscMalloc1(numRootIndices,&rootVals);
3692: PetscSFBcastBegin(indicesSF,MPIU_SCALAR,leafVals,rootVals);
3693: PetscSFBcastEnd(indicesSF,MPIU_SCALAR,leafVals,rootVals);
3694: }
3695: PetscSFDestroy(&indicesSF);
3696: }
3697: PetscSectionDestroy(&leafIndicesSec);
3698: PetscFree(leafInds);
3699: PetscFree(leafVals);
3700: PetscSFDestroy(&coarseToFineEmbedded);
3701: *rootMultiSec = multiRootSec;
3702: *multiLeafSec = rootIndicesSec;
3703: if (gatheredIndices) *gatheredIndices = rootInds;
3704: if (gatheredValues) *gatheredValues = rootVals;
3705: return(0);
3706: }
3708: PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
3709: {
3710: DM refTree;
3711: PetscSection multiRootSec, rootIndicesSec;
3712: PetscSection globalCoarse, globalFine;
3713: PetscSection localCoarse, localFine;
3714: PetscSection cSecRef;
3715: PetscInt *rootIndices, *parentIndices, pRefStart, pRefEnd;
3716: Mat injRef;
3717: PetscInt numFields, maxDof;
3718: PetscInt pStartC, pEndC, pStartF, pEndF, p;
3719: PetscInt *offsets, *offsetsCopy, *rowOffsets;
3720: PetscLayout rowMap, colMap;
3721: PetscInt rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO;
3722: PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */
3727: /* get the templates for the fine-to-coarse injection from the reference tree */
3728: DMPlexGetReferenceTree(coarse,&refTree);
3729: DMGetDefaultConstraints(refTree,&cSecRef,NULL);
3730: PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);
3731: DMPlexReferenceTreeGetInjector(refTree,&injRef);
3733: DMPlexGetChart(fine,&pStartF,&pEndF);
3734: DMGetSection(fine,&localFine);
3735: DMGetGlobalSection(fine,&globalFine);
3736: PetscSectionGetNumFields(localFine,&numFields);
3737: DMPlexGetChart(coarse,&pStartC,&pEndC);
3738: DMGetSection(coarse,&localCoarse);
3739: DMGetGlobalSection(coarse,&globalCoarse);
3740: PetscSectionGetMaxDof(localCoarse,&maxDof);
3741: {
3742: PetscInt maxFields = PetscMax(1,numFields) + 1;
3743: PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);
3744: }
3746: DMPlexTransferInjectorTree(coarse,fine,coarseToFine,childIds,NULL,numFields,offsets,&multiRootSec,&rootIndicesSec,&rootIndices,NULL);
3748: PetscMalloc1(maxDof,&parentIndices);
3750: /* count indices */
3751: MatGetLayouts(mat,&rowMap,&colMap);
3752: PetscLayoutSetUp(rowMap);
3753: PetscLayoutSetUp(colMap);
3754: PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);
3755: PetscLayoutGetRange(colMap,&colStart,&colEnd);
3756: PetscCalloc2(rowEnd-rowStart,&nnzD,rowEnd-rowStart,&nnzO);
3757: for (p = pStartC; p < pEndC; p++) {
3758: PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
3760: PetscSectionGetDof(globalCoarse,p,&dof);
3761: PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3762: if ((dof - cdof) <= 0) continue;
3763: PetscSectionGetOffset(globalCoarse,p,&gOff);
3765: rowOffsets[0] = 0;
3766: offsetsCopy[0] = 0;
3767: if (numFields) {
3768: PetscInt f;
3770: for (f = 0; f < numFields; f++) {
3771: PetscInt fDof;
3772: PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
3773: rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3774: }
3775: DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);
3776: }
3777: else {
3778: DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);
3779: rowOffsets[1] = offsetsCopy[0];
3780: }
3782: PetscSectionGetDof(multiRootSec,p,&numLeaves);
3783: PetscSectionGetOffset(multiRootSec,p,&leafStart);
3784: leafEnd = leafStart + numLeaves;
3785: for (l = leafStart; l < leafEnd; l++) {
3786: PetscInt numIndices, childId, offset;
3787: const PetscInt *childIndices;
3789: PetscSectionGetDof(rootIndicesSec,l,&numIndices);
3790: PetscSectionGetOffset(rootIndicesSec,l,&offset);
3791: childId = rootIndices[offset++];
3792: childIndices = &rootIndices[offset];
3793: numIndices--;
3795: if (childId == -1) { /* equivalent points: scatter */
3796: PetscInt i;
3798: for (i = 0; i < numIndices; i++) {
3799: PetscInt colIndex = childIndices[i];
3800: PetscInt rowIndex = parentIndices[i];
3801: if (rowIndex < 0) continue;
3802: if (colIndex < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unconstrained fine and constrained coarse");
3803: if (colIndex >= colStart && colIndex < colEnd) {
3804: nnzD[rowIndex - rowStart] = 1;
3805: }
3806: else {
3807: nnzO[rowIndex - rowStart] = 1;
3808: }
3809: }
3810: }
3811: else {
3812: PetscInt parentId, f, lim;
3814: DMPlexGetTreeParent(refTree,childId,&parentId,NULL);
3816: lim = PetscMax(1,numFields);
3817: offsets[0] = 0;
3818: if (numFields) {
3819: PetscInt f;
3821: for (f = 0; f < numFields; f++) {
3822: PetscInt fDof;
3823: PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);
3825: offsets[f + 1] = fDof + offsets[f];
3826: }
3827: }
3828: else {
3829: PetscInt cDof;
3831: PetscSectionGetDof(cSecRef,childId,&cDof);
3832: offsets[1] = cDof;
3833: }
3834: for (f = 0; f < lim; f++) {
3835: PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1];
3836: PetscInt childStart = offsets[f], childEnd = offsets[f + 1];
3837: PetscInt i, numD = 0, numO = 0;
3839: for (i = childStart; i < childEnd; i++) {
3840: PetscInt colIndex = childIndices[i];
3842: if (colIndex < 0) continue;
3843: if (colIndex >= colStart && colIndex < colEnd) {
3844: numD++;
3845: }
3846: else {
3847: numO++;
3848: }
3849: }
3850: for (i = parentStart; i < parentEnd; i++) {
3851: PetscInt rowIndex = parentIndices[i];
3853: if (rowIndex < 0) continue;
3854: nnzD[rowIndex - rowStart] += numD;
3855: nnzO[rowIndex - rowStart] += numO;
3856: }
3857: }
3858: }
3859: }
3860: }
3861: /* preallocate */
3862: MatXAIJSetPreallocation(mat,1,nnzD,nnzO,NULL,NULL);
3863: PetscFree2(nnzD,nnzO);
3864: /* insert values */
3865: DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);
3866: for (p = pStartC; p < pEndC; p++) {
3867: PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
3869: PetscSectionGetDof(globalCoarse,p,&dof);
3870: PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3871: if ((dof - cdof) <= 0) continue;
3872: PetscSectionGetOffset(globalCoarse,p,&gOff);
3874: rowOffsets[0] = 0;
3875: offsetsCopy[0] = 0;
3876: if (numFields) {
3877: PetscInt f;
3879: for (f = 0; f < numFields; f++) {
3880: PetscInt fDof;
3881: PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
3882: rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3883: }
3884: DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);
3885: }
3886: else {
3887: DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);
3888: rowOffsets[1] = offsetsCopy[0];
3889: }
3891: PetscSectionGetDof(multiRootSec,p,&numLeaves);
3892: PetscSectionGetOffset(multiRootSec,p,&leafStart);
3893: leafEnd = leafStart + numLeaves;
3894: for (l = leafStart; l < leafEnd; l++) {
3895: PetscInt numIndices, childId, offset;
3896: const PetscInt *childIndices;
3898: PetscSectionGetDof(rootIndicesSec,l,&numIndices);
3899: PetscSectionGetOffset(rootIndicesSec,l,&offset);
3900: childId = rootIndices[offset++];
3901: childIndices = &rootIndices[offset];
3902: numIndices--;
3904: if (childId == -1) { /* equivalent points: scatter */
3905: PetscInt i;
3907: for (i = 0; i < numIndices; i++) {
3908: MatSetValue(mat,parentIndices[i],childIndices[i],1.,INSERT_VALUES);
3909: }
3910: }
3911: else {
3912: PetscInt parentId, f, lim;
3914: DMPlexGetTreeParent(refTree,childId,&parentId,NULL);
3916: lim = PetscMax(1,numFields);
3917: offsets[0] = 0;
3918: if (numFields) {
3919: PetscInt f;
3921: for (f = 0; f < numFields; f++) {
3922: PetscInt fDof;
3923: PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);
3925: offsets[f + 1] = fDof + offsets[f];
3926: }
3927: }
3928: else {
3929: PetscInt cDof;
3931: PetscSectionGetDof(cSecRef,childId,&cDof);
3932: offsets[1] = cDof;
3933: }
3934: for (f = 0; f < lim; f++) {
3935: PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0];
3936: PetscInt *rowIndices = &parentIndices[rowOffsets[f]];
3937: const PetscInt *colIndices = &childIndices[offsets[f]];
3939: MatSetValues(mat,rowOffsets[f+1]-rowOffsets[f],rowIndices,offsets[f+1]-offsets[f],colIndices,childMat,INSERT_VALUES);
3940: }
3941: }
3942: }
3943: }
3944: PetscSectionDestroy(&multiRootSec);
3945: PetscSectionDestroy(&rootIndicesSec);
3946: PetscFree(parentIndices);
3947: DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);
3948: PetscFree(rootIndices);
3949: PetscFree3(offsets,offsetsCopy,rowOffsets);
3951: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
3952: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
3953: return(0);
3954: }
3956: static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom)
3957: {
3958: PetscDS ds;
3959: PetscSF coarseToFineEmbedded;
3960: PetscSection globalCoarse, globalFine;
3961: PetscSection localCoarse, localFine;
3962: PetscSection aSec, cSec;
3963: PetscSection rootValuesSec;
3964: PetscSection leafValuesSec;
3965: PetscScalar *rootValues, *leafValues;
3966: IS aIS;
3967: const PetscInt *anchors;
3968: Mat cMat;
3969: PetscInt numFields;
3970: PetscInt pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd, cellEndInterior;
3971: PetscInt aStart, aEnd, cStart, cEnd;
3972: PetscInt *maxChildIds;
3973: PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
3974: PetscFV fv = NULL;
3975: PetscInt dim, numFVcomps = -1, fvField = -1;
3976: DM cellDM = NULL, gradDM = NULL;
3977: const PetscScalar *cellGeomArray = NULL;
3978: const PetscScalar *gradArray = NULL;
3979: PetscErrorCode ierr;
3982: VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);
3983: DMPlexGetChart(coarse,&pStartC,&pEndC);
3984: DMPlexGetHeightStratum(coarse,0,&cellStart,&cellEnd);
3985: DMPlexGetHybridBounds(coarse,&cellEndInterior,NULL,NULL,NULL);
3986: cellEnd = (cellEndInterior < 0) ? cellEnd : cellEndInterior;
3987: DMPlexGetChart(fine,&pStartF,&pEndF);
3988: DMGetGlobalSection(fine,&globalFine);
3989: DMGetCoordinateDim(coarse,&dim);
3990: { /* winnow fine points that don't have global dofs out of the sf */
3991: PetscInt nleaves, l;
3992: const PetscInt *leaves;
3993: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;
3995: PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);
3997: for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) {
3998: PetscInt p = leaves ? leaves[l] : l;
4000: PetscSectionGetDof(globalFine,p,&dof);
4001: PetscSectionGetConstraintDof(globalFine,p,&cdof);
4002: if ((dof - cdof) > 0) {
4003: numPointsWithDofs++;
4004: }
4005: }
4006: PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
4007: for (l = 0, offset = 0; l < nleaves; l++) {
4008: PetscInt p = leaves ? leaves[l] : l;
4010: PetscSectionGetDof(globalFine,p,&dof);
4011: PetscSectionGetConstraintDof(globalFine,p,&cdof);
4012: if ((dof - cdof) > 0) {
4013: pointsWithDofs[offset++] = l;
4014: }
4015: }
4016: PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
4017: PetscFree(pointsWithDofs);
4018: }
4019: /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
4020: PetscMalloc1(pEndC-pStartC,&maxChildIds);
4021: for (p = pStartC; p < pEndC; p++) {
4022: maxChildIds[p - pStartC] = -2;
4023: }
4024: PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);
4025: PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);
4027: DMGetSection(coarse,&localCoarse);
4028: DMGetGlobalSection(coarse,&globalCoarse);
4030: DMPlexGetAnchors(coarse,&aSec,&aIS);
4031: ISGetIndices(aIS,&anchors);
4032: PetscSectionGetChart(aSec,&aStart,&aEnd);
4034: DMGetDefaultConstraints(coarse,&cSec,&cMat);
4035: PetscSectionGetChart(cSec,&cStart,&cEnd);
4037: /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
4038: PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootValuesSec);
4039: PetscSectionSetChart(rootValuesSec,pStartC,pEndC);
4040: PetscSectionGetNumFields(localCoarse,&numFields);
4041: {
4042: PetscInt maxFields = PetscMax(1,numFields) + 1;
4043: PetscMalloc7(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&newOffsets,maxFields,&newOffsetsCopy,maxFields,&rowOffsets,maxFields,&numD,maxFields,&numO);
4044: }
4045: if (grad) {
4046: PetscInt i;
4048: VecGetDM(cellGeom,&cellDM);
4049: VecGetArrayRead(cellGeom,&cellGeomArray);
4050: VecGetDM(grad,&gradDM);
4051: VecGetArrayRead(grad,&gradArray);
4052: for (i = 0; i < PetscMax(1,numFields); i++) {
4053: PetscObject obj;
4054: PetscClassId id;
4056: DMGetField(coarse, i, &obj);
4057: PetscObjectGetClassId(obj,&id);
4058: if (id == PETSCFV_CLASSID) {
4059: fv = (PetscFV) obj;
4060: PetscFVGetNumComponents(fv,&numFVcomps);
4061: fvField = i;
4062: break;
4063: }
4064: }
4065: }
4067: for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
4068: PetscInt dof;
4069: PetscInt maxChildId = maxChildIds[p - pStartC];
4070: PetscInt numValues = 0;
4072: PetscSectionGetDof(globalCoarse,p,&dof);
4073: if (dof < 0) {
4074: dof = -(dof + 1);
4075: }
4076: offsets[0] = 0;
4077: newOffsets[0] = 0;
4078: if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
4079: PetscInt *closure = NULL, closureSize, cl;
4081: DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);
4082: for (cl = 0; cl < closureSize; cl++) { /* get the closure */
4083: PetscInt c = closure[2 * cl], clDof;
4085: PetscSectionGetDof(localCoarse,c,&clDof);
4086: numValues += clDof;
4087: }
4088: DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);
4089: }
4090: else if (maxChildId == -1) {
4091: PetscSectionGetDof(localCoarse,p,&numValues);
4092: }
4093: /* we will pack the column indices with the field offsets */
4094: if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) {
4095: /* also send the centroid, and the gradient */
4096: numValues += dim * (1 + numFVcomps);
4097: }
4098: PetscSectionSetDof(rootValuesSec,p,numValues);
4099: }
4100: PetscSectionSetUp(rootValuesSec);
4101: {
4102: PetscInt numRootValues;
4103: const PetscScalar *coarseArray;
4105: PetscSectionGetStorageSize(rootValuesSec,&numRootValues);
4106: PetscMalloc1(numRootValues,&rootValues);
4107: VecGetArrayRead(vecCoarseLocal,&coarseArray);
4108: for (p = pStartC; p < pEndC; p++) {
4109: PetscInt numValues;
4110: PetscInt pValOff;
4111: PetscScalar *pVal;
4112: PetscInt maxChildId = maxChildIds[p - pStartC];
4114: PetscSectionGetDof(rootValuesSec,p,&numValues);
4115: if (!numValues) {
4116: continue;
4117: }
4118: PetscSectionGetOffset(rootValuesSec,p,&pValOff);
4119: pVal = &(rootValues[pValOff]);
4120: if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
4121: PetscInt closureSize = numValues;
4122: DMPlexVecGetClosure(coarse,NULL,vecCoarseLocal,p,&closureSize,&pVal);
4123: if (grad && p >= cellStart && p < cellEnd) {
4124: PetscFVCellGeom *cg;
4125: PetscScalar *gradVals = NULL;
4126: PetscInt i;
4128: pVal += (numValues - dim * (1 + numFVcomps));
4130: DMPlexPointLocalRead(cellDM,p,cellGeomArray,(void *) &cg);
4131: for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i];
4132: pVal += dim;
4133: DMPlexPointGlobalRead(gradDM,p,gradArray,(void *) &gradVals);
4134: for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i];
4135: }
4136: }
4137: else if (maxChildId == -1) {
4138: PetscInt lDof, lOff, i;
4140: PetscSectionGetDof(localCoarse,p,&lDof);
4141: PetscSectionGetOffset(localCoarse,p,&lOff);
4142: for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i];
4143: }
4144: }
4145: VecRestoreArrayRead(vecCoarseLocal,&coarseArray);
4146: PetscFree(maxChildIds);
4147: }
4148: {
4149: PetscSF valuesSF;
4150: PetscInt *remoteOffsetsValues, numLeafValues;
4152: PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafValuesSec);
4153: PetscSFDistributeSection(coarseToFineEmbedded,rootValuesSec,&remoteOffsetsValues,leafValuesSec);
4154: PetscSFCreateSectionSF(coarseToFineEmbedded,rootValuesSec,remoteOffsetsValues,leafValuesSec,&valuesSF);
4155: PetscSFDestroy(&coarseToFineEmbedded);
4156: PetscFree(remoteOffsetsValues);
4157: PetscSectionGetStorageSize(leafValuesSec,&numLeafValues);
4158: PetscMalloc1(numLeafValues,&leafValues);
4159: PetscSFBcastBegin(valuesSF,MPIU_SCALAR,rootValues,leafValues);
4160: PetscSFBcastEnd(valuesSF,MPIU_SCALAR,rootValues,leafValues);
4161: PetscSFDestroy(&valuesSF);
4162: PetscFree(rootValues);
4163: PetscSectionDestroy(&rootValuesSec);
4164: }
4165: DMGetSection(fine,&localFine);
4166: {
4167: PetscInt maxDof;
4168: PetscInt *rowIndices;
4169: DM refTree;
4170: PetscInt **refPointFieldN;
4171: PetscScalar ***refPointFieldMats;
4172: PetscSection refConSec, refAnSec;
4173: PetscInt pRefStart,pRefEnd,leafStart,leafEnd;
4174: PetscScalar *pointWork;
4176: PetscSectionGetMaxDof(globalFine,&maxDof);
4177: DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices);
4178: DMGetWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork);
4179: DMPlexGetReferenceTree(fine,&refTree);
4180: DMGetDS(fine,&ds);
4181: DMSetDS(refTree,ds);
4182: DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
4183: DMGetDefaultConstraints(refTree,&refConSec,NULL);
4184: DMPlexGetAnchors(refTree,&refAnSec,NULL);
4185: PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
4186: PetscSectionGetChart(leafValuesSec,&leafStart,&leafEnd);
4187: DMPlexGetHeightStratum(fine,0,&cellStart,&cellEnd);
4188: DMPlexGetHybridBounds(fine,&cellEndInterior,NULL,NULL,NULL);
4189: cellEnd = (cellEndInterior < 0) ? cellEnd : cellEndInterior;
4190: for (p = leafStart; p < leafEnd; p++) {
4191: PetscInt gDof, gcDof, gOff, lDof;
4192: PetscInt numValues, pValOff;
4193: PetscInt childId;
4194: const PetscScalar *pVal;
4195: const PetscScalar *fvGradData = NULL;
4197: PetscSectionGetDof(globalFine,p,&gDof);
4198: PetscSectionGetDof(localFine,p,&lDof);
4199: PetscSectionGetConstraintDof(globalFine,p,&gcDof);
4200: if ((gDof - gcDof) <= 0) {
4201: continue;
4202: }
4203: PetscSectionGetOffset(globalFine,p,&gOff);
4204: PetscSectionGetDof(leafValuesSec,p,&numValues);
4205: if (!numValues) continue;
4206: PetscSectionGetOffset(leafValuesSec,p,&pValOff);
4207: pVal = &leafValues[pValOff];
4208: offsets[0] = 0;
4209: offsetsCopy[0] = 0;
4210: newOffsets[0] = 0;
4211: newOffsetsCopy[0] = 0;
4212: childId = cids[p - pStartF];
4213: if (numFields) {
4214: PetscInt f;
4215: for (f = 0; f < numFields; f++) {
4216: PetscInt rowDof;
4218: PetscSectionGetFieldDof(localFine,p,f,&rowDof);
4219: offsets[f + 1] = offsets[f] + rowDof;
4220: offsetsCopy[f + 1] = offsets[f + 1];
4221: /* TODO: closure indices */
4222: newOffsets[f + 1] = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]);
4223: }
4224: DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);
4225: }
4226: else {
4227: offsets[0] = 0;
4228: offsets[1] = lDof;
4229: newOffsets[0] = 0;
4230: newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0];
4231: DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);
4232: }
4233: if (childId == -1) { /* no child interpolation: one nnz per */
4234: VecSetValues(vecFine,numValues,rowIndices,pVal,INSERT_VALUES);
4235: } else {
4236: PetscInt f;
4238: if (grad && p >= cellStart && p < cellEnd) {
4239: numValues -= (dim * (1 + numFVcomps));
4240: fvGradData = &pVal[numValues];
4241: }
4242: for (f = 0; f < PetscMax(1,numFields); f++) {
4243: const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f];
4244: PetscInt numRows = offsets[f+1] - offsets[f];
4245: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
4246: const PetscScalar *cVal = &pVal[newOffsets[f]];
4247: PetscScalar *rVal = &pointWork[offsets[f]];
4248: PetscInt i, j;
4250: #if 0
4251: PetscInfo6(coarse,"field %D, numFields %D, childId %D, numRows %D, numCols %D, refPointFieldN %D\n",f,numFields,childId,numRows,numCols,refPointFieldN[childId - pRefStart][f]);
4252: #endif
4253: for (i = 0; i < numRows; i++) {
4254: PetscScalar val = 0.;
4255: for (j = 0; j < numCols; j++) {
4256: val += childMat[i * numCols + j] * cVal[j];
4257: }
4258: rVal[i] = val;
4259: }
4260: if (f == fvField && p >= cellStart && p < cellEnd) {
4261: PetscReal centroid[3];
4262: PetscScalar diff[3];
4263: const PetscScalar *parentCentroid = &fvGradData[0];
4264: const PetscScalar *gradient = &fvGradData[dim];
4266: DMPlexComputeCellGeometryFVM(fine,p,NULL,centroid,NULL);
4267: for (i = 0; i < dim; i++) {
4268: diff[i] = centroid[i] - parentCentroid[i];
4269: }
4270: for (i = 0; i < numFVcomps; i++) {
4271: PetscScalar val = 0.;
4273: for (j = 0; j < dim; j++) {
4274: val += gradient[dim * i + j] * diff[j];
4275: }
4276: rVal[i] += val;
4277: }
4278: }
4279: VecSetValues(vecFine,numRows,&rowIndices[offsets[f]],rVal,INSERT_VALUES);
4280: }
4281: }
4282: }
4283: DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
4284: DMRestoreWorkArray(fine,maxDof,MPIU_INT,&rowIndices);
4285: DMRestoreWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork);
4286: }
4287: PetscFree(leafValues);
4288: PetscSectionDestroy(&leafValuesSec);
4289: PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);
4290: ISRestoreIndices(aIS,&anchors);
4291: return(0);
4292: }
4294: static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids)
4295: {
4296: PetscDS ds;
4297: DM refTree;
4298: PetscSection multiRootSec, rootIndicesSec;
4299: PetscSection globalCoarse, globalFine;
4300: PetscSection localCoarse, localFine;
4301: PetscSection cSecRef;
4302: PetscInt *parentIndices, pRefStart, pRefEnd;
4303: PetscScalar *rootValues, *parentValues;
4304: Mat injRef;
4305: PetscInt numFields, maxDof;
4306: PetscInt pStartC, pEndC, pStartF, pEndF, p;
4307: PetscInt *offsets, *offsetsCopy, *rowOffsets;
4308: PetscLayout rowMap, colMap;
4309: PetscInt rowStart, rowEnd, colStart, colEnd;
4310: PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */
4315: /* get the templates for the fine-to-coarse injection from the reference tree */
4316: VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);
4317: VecSetOption(vecCoarse,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);
4318: DMPlexGetReferenceTree(coarse,&refTree);
4319: DMGetDS(coarse,&ds);
4320: DMSetDS(refTree,ds);
4321: DMGetDefaultConstraints(refTree,&cSecRef,NULL);
4322: PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);
4323: DMPlexReferenceTreeGetInjector(refTree,&injRef);
4325: DMPlexGetChart(fine,&pStartF,&pEndF);
4326: DMGetSection(fine,&localFine);
4327: DMGetGlobalSection(fine,&globalFine);
4328: PetscSectionGetNumFields(localFine,&numFields);
4329: DMPlexGetChart(coarse,&pStartC,&pEndC);
4330: DMGetSection(coarse,&localCoarse);
4331: DMGetGlobalSection(coarse,&globalCoarse);
4332: PetscSectionGetMaxDof(localCoarse,&maxDof);
4333: {
4334: PetscInt maxFields = PetscMax(1,numFields) + 1;
4335: PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);
4336: }
4338: DMPlexTransferInjectorTree(coarse,fine,coarseToFine,cids,vecFine,numFields,offsets,&multiRootSec,&rootIndicesSec,NULL,&rootValues);
4340: PetscMalloc2(maxDof,&parentIndices,maxDof,&parentValues);
4342: /* count indices */
4343: VecGetLayout(vecFine,&colMap);
4344: VecGetLayout(vecCoarse,&rowMap);
4345: PetscLayoutSetUp(rowMap);
4346: PetscLayoutSetUp(colMap);
4347: PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);
4348: PetscLayoutGetRange(colMap,&colStart,&colEnd);
4349: /* insert values */
4350: DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);
4351: for (p = pStartC; p < pEndC; p++) {
4352: PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
4353: PetscBool contribute = PETSC_FALSE;
4355: PetscSectionGetDof(globalCoarse,p,&dof);
4356: PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
4357: if ((dof - cdof) <= 0) continue;
4358: PetscSectionGetDof(localCoarse,p,&dof);
4359: PetscSectionGetOffset(globalCoarse,p,&gOff);
4361: rowOffsets[0] = 0;
4362: offsetsCopy[0] = 0;
4363: if (numFields) {
4364: PetscInt f;
4366: for (f = 0; f < numFields; f++) {
4367: PetscInt fDof;
4368: PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
4369: rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
4370: }
4371: DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);
4372: }
4373: else {
4374: DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);
4375: rowOffsets[1] = offsetsCopy[0];
4376: }
4378: PetscSectionGetDof(multiRootSec,p,&numLeaves);
4379: PetscSectionGetOffset(multiRootSec,p,&leafStart);
4380: leafEnd = leafStart + numLeaves;
4381: for (l = 0; l < dof; l++) parentValues[l] = 0.;
4382: for (l = leafStart; l < leafEnd; l++) {
4383: PetscInt numIndices, childId, offset;
4384: const PetscScalar *childValues;
4386: PetscSectionGetDof(rootIndicesSec,l,&numIndices);
4387: PetscSectionGetOffset(rootIndicesSec,l,&offset);
4388: childId = (PetscInt) PetscRealPart(rootValues[offset++]);
4389: childValues = &rootValues[offset];
4390: numIndices--;
4392: if (childId == -2) { /* skip */
4393: continue;
4394: } else if (childId == -1) { /* equivalent points: scatter */
4395: PetscInt m;
4397: contribute = PETSC_TRUE;
4398: for (m = 0; m < numIndices; m++) parentValues[m] = childValues[m];
4399: } else { /* contributions from children: sum with injectors from reference tree */
4400: PetscInt parentId, f, lim;
4402: contribute = PETSC_TRUE;
4403: DMPlexGetTreeParent(refTree,childId,&parentId,NULL);
4405: lim = PetscMax(1,numFields);
4406: offsets[0] = 0;
4407: if (numFields) {
4408: PetscInt f;
4410: for (f = 0; f < numFields; f++) {
4411: PetscInt fDof;
4412: PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);
4414: offsets[f + 1] = fDof + offsets[f];
4415: }
4416: }
4417: else {
4418: PetscInt cDof;
4420: PetscSectionGetDof(cSecRef,childId,&cDof);
4421: offsets[1] = cDof;
4422: }
4423: for (f = 0; f < lim; f++) {
4424: PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0];
4425: PetscInt n = offsets[f+1]-offsets[f];
4426: PetscInt m = rowOffsets[f+1]-rowOffsets[f];
4427: PetscInt i, j;
4428: const PetscScalar *colValues = &childValues[offsets[f]];
4430: for (i = 0; i < m; i++) {
4431: PetscScalar val = 0.;
4432: for (j = 0; j < n; j++) {
4433: val += childMat[n * i + j] * colValues[j];
4434: }
4435: parentValues[rowOffsets[f] + i] += val;
4436: }
4437: }
4438: }
4439: }
4440: if (contribute) {VecSetValues(vecCoarse,dof,parentIndices,parentValues,INSERT_VALUES);}
4441: }
4442: PetscSectionDestroy(&multiRootSec);
4443: PetscSectionDestroy(&rootIndicesSec);
4444: PetscFree2(parentIndices,parentValues);
4445: DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);
4446: PetscFree(rootValues);
4447: PetscFree3(offsets,offsetsCopy,rowOffsets);
4448: return(0);
4449: }
4451: /*@
4452: DMPlexTransferVecTree - transfer a vector between two meshes that differ from each other by refinement/coarsening
4453: that can be represented by a common reference tree used by both. This routine can be used for a combination of
4454: coarsening and refinement at the same time.
4456: collective
4458: Input Parameters:
4459: + dmIn - The DMPlex mesh for the input vector
4460: . vecIn - The input vector
4461: . sfRefine - A star forest indicating points in the mesh dmIn (roots in the star forest) that are parents to points in
4462: the mesh dmOut (leaves in the star forest), i.e. where dmOut is more refined than dmIn
4463: . sfCoarsen - A star forest indicating points in the mesh dmOut (roots in the star forest) that are parents to points in
4464: the mesh dmIn (leaves in the star forest), i.e. where dmOut is more coarsened than dmIn
4465: . cidsRefine - The childIds of the points in dmOut. These childIds relate back to the reference tree: childid[j] = k implies
4466: that mesh point j of dmOut was refined from a point in dmIn just as the mesh point k in the reference
4467: tree was refined from its parent. childid[j] = -1 indicates that the point j in dmOut is exactly
4468: equivalent to its root in dmIn, so no interpolation is necessary. childid[j] = -2 indicates that this
4469: point j in dmOut is not a leaf of sfRefine.
4470: . cidsCoarsen - The childIds of the points in dmIn. These childIds relate back to the reference tree: childid[j] = k implies
4471: that mesh point j of dmIn coarsens to a point in dmOut just as the mesh point k in the reference
4472: tree coarsens to its parent. childid[j] = -2 indicates that point j in dmOut is not a leaf in sfCoarsen.
4473: . useBCs - PETSC_TRUE indicates that boundary values should be inserted into vecIn before transfer.
4474: - time - Used if boundary values are time dependent.
4476: Output Parameters:
4477: . vecOut - Using interpolation and injection operators calculated on the reference tree, the transfered
4478: projection of vecIn from dmIn to dmOut. Note that any field discretized with a PetscFV finite volume
4479: method that uses gradient reconstruction will use reconstructed gradients when interpolating from
4480: coarse points to fine points.
4482: Level: developer
4484: .seealso(): DMPlexSetReferenceTree(), DMPlexGetReferenceTree(), PetscFVGetComputeGradients()
4485: @*/
4486: PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time)
4487: {
4491: VecSet(vecOut,0.0);
4492: if (sfRefine) {
4493: Vec vecInLocal;
4494: DM dmGrad = NULL;
4495: Vec faceGeom = NULL, cellGeom = NULL, grad = NULL;
4497: DMGetLocalVector(dmIn,&vecInLocal);
4498: VecSet(vecInLocal,0.0);
4499: {
4500: PetscInt numFields, i;
4502: DMGetNumFields(dmIn, &numFields);
4503: for (i = 0; i < numFields; i++) {
4504: PetscObject obj;
4505: PetscClassId classid;
4507: DMGetField(dmIn, i, &obj);
4508: PetscObjectGetClassId(obj, &classid);
4509: if (classid == PETSCFV_CLASSID) {
4510: DMPlexGetDataFVM(dmIn,(PetscFV)obj,&cellGeom,&faceGeom,&dmGrad);
4511: break;
4512: }
4513: }
4514: }
4515: if (useBCs) {
4516: DMPlexInsertBoundaryValues(dmIn,PETSC_TRUE,vecInLocal,time,faceGeom,cellGeom,NULL);
4517: }
4518: DMGlobalToLocalBegin(dmIn,vecIn,INSERT_VALUES,vecInLocal);
4519: DMGlobalToLocalEnd(dmIn,vecIn,INSERT_VALUES,vecInLocal);
4520: if (dmGrad) {
4521: DMGetGlobalVector(dmGrad,&grad);
4522: DMPlexReconstructGradientsFVM(dmIn,vecInLocal,grad);
4523: }
4524: DMPlexTransferVecTree_Interpolate(dmIn,vecInLocal,dmOut,vecOut,sfRefine,cidsRefine,grad,cellGeom);
4525: DMRestoreLocalVector(dmIn,&vecInLocal);
4526: if (dmGrad) {
4527: DMRestoreGlobalVector(dmGrad,&grad);
4528: }
4529: }
4530: if (sfCoarsen) {
4531: DMPlexTransferVecTree_Inject(dmIn,vecIn,dmOut,vecOut,sfCoarsen,cidsCoarsen);
4532: }
4533: VecAssemblyBegin(vecOut);
4534: VecAssemblyEnd(vecOut);
4535: return(0);
4536: }