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