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