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