Actual source code: plextree.c
petsc-3.7.5 2017-01-01
1: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2: #include <../src/sys/utils/hash.h>
3: #include <petsc/private/isimpl.h>
4: #include <petsc/private/petscfeimpl.h>
5: #include <petscsf.h>
6: #include <petscds.h>
8: /** hierarchy routines */
12: /*@
13: DMPlexSetReferenceTree - set the reference tree for hierarchically non-conforming meshes.
15: Not collective
17: Input Parameters:
18: + dm - The DMPlex object
19: - ref - The reference tree DMPlex object
21: Level: intermediate
23: .seealso: DMPlexGetReferenceTree(), DMPlexCreateDefaultReferenceTree()
24: @*/
25: PetscErrorCode DMPlexSetReferenceTree(DM dm, DM ref)
26: {
27: DM_Plex *mesh = (DM_Plex *)dm->data;
28: PetscErrorCode ierr;
33: PetscObjectReference((PetscObject)ref);
34: DMDestroy(&mesh->referenceTree);
35: mesh->referenceTree = ref;
36: return(0);
37: }
41: /*@
42: DMPlexGetReferenceTree - get the reference tree for hierarchically non-conforming meshes.
44: Not collective
46: Input Parameters:
47: . dm - The DMPlex object
49: Output Parameters
50: . ref - The reference tree DMPlex object
52: Level: intermediate
54: .seealso: DMPlexSetReferenceTree(), DMPlexCreateDefaultReferenceTree()
55: @*/
56: PetscErrorCode DMPlexGetReferenceTree(DM dm, DM *ref)
57: {
58: DM_Plex *mesh = (DM_Plex *)dm->data;
63: *ref = mesh->referenceTree;
64: return(0);
65: }
69: static PetscErrorCode DMPlexReferenceTreeGetChildSymmetry_Default(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
70: {
71: PetscInt coneSize, dStart, dEnd, dim, ABswap, oAvert, oBvert, ABswapVert;
75: if (parentOrientA == parentOrientB) {
76: if (childOrientB) *childOrientB = childOrientA;
77: if (childB) *childB = childA;
78: return(0);
79: }
80: for (dim = 0; dim < 3; dim++) {
81: DMPlexGetDepthStratum(dm,dim,&dStart,&dEnd);
82: if (parent >= dStart && parent <= dEnd) {
83: break;
84: }
85: }
86: if (dim > 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot perform child symmetry for %d-cells",dim);
87: if (!dim) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A vertex has no children");
88: if (childA < dStart || childA >= dEnd) {
89: /* this is a lower-dimensional child: bootstrap */
90: PetscInt size, i, sA = -1, sB, sOrientB, sConeSize;
91: const PetscInt *supp, *coneA, *coneB, *oA, *oB;
93: DMPlexGetSupportSize(dm,childA,&size);
94: DMPlexGetSupport(dm,childA,&supp);
96: /* find a point sA in supp(childA) that has the same parent */
97: for (i = 0; i < size; i++) {
98: PetscInt sParent;
100: sA = supp[i];
101: if (sA == parent) continue;
102: DMPlexGetTreeParent(dm,sA,&sParent,NULL);
103: if (sParent == parent) {
104: break;
105: }
106: }
107: if (i == size) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"could not find support in children");
108: /* find out which point sB is in an equivalent position to sA under
109: * parentOrientB */
110: DMPlexReferenceTreeGetChildSymmetry_Default(dm,parent,parentOrientA,0,sA,parentOrientB,&sOrientB,&sB);
111: DMPlexGetConeSize(dm,sA,&sConeSize);
112: DMPlexGetCone(dm,sA,&coneA);
113: DMPlexGetCone(dm,sB,&coneB);
114: DMPlexGetConeOrientation(dm,sA,&oA);
115: DMPlexGetConeOrientation(dm,sB,&oB);
116: /* step through the cone of sA in natural order */
117: for (i = 0; i < sConeSize; i++) {
118: if (coneA[i] == childA) {
119: /* if childA is at position i in coneA,
120: * then we want the point that is at sOrientB*i in coneB */
121: PetscInt j = (sOrientB >= 0) ? ((sOrientB + i) % sConeSize) : ((sConeSize -(sOrientB+1) - i) % sConeSize);
122: if (childB) *childB = coneB[j];
123: if (childOrientB) {
124: PetscInt oBtrue;
126: DMPlexGetConeSize(dm,childA,&coneSize);
127: /* compose sOrientB and oB[j] */
128: if (coneSize != 0 && coneSize != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expected a vertex or an edge");
129: /* we may have to flip an edge */
130: oBtrue = coneSize ? ((sOrientB >= 0) ? oB[j] : -(oB[j] + 2)) : 0;
131: ABswap = DihedralSwap(coneSize,oA[i],oBtrue);
132: *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap);
133: }
134: break;
135: }
136: }
137: if (i == sConeSize) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"support cone mismatch");
138: return(0);
139: }
140: /* get the cone size and symmetry swap */
141: DMPlexGetConeSize(dm,parent,&coneSize);
142: ABswap = DihedralSwap(coneSize, parentOrientA, parentOrientB);
143: if (dim == 2) {
144: /* orientations refer to cones: we want them to refer to vertices:
145: * if it's a rotation, they are the same, but if the order is reversed, a
146: * permutation that puts side i first does *not* put vertex i first */
147: oAvert = (parentOrientA >= 0) ? parentOrientA : -((-parentOrientA % coneSize) + 1);
148: oBvert = (parentOrientB >= 0) ? parentOrientB : -((-parentOrientB % coneSize) + 1);
149: ABswapVert = DihedralSwap(coneSize, oAvert, oBvert);
150: } else {
151: ABswapVert = ABswap;
152: }
153: if (childB) {
154: /* assume that each child corresponds to a vertex, in the same order */
155: PetscInt p, posA = -1, numChildren, i;
156: const PetscInt *children;
158: /* count which position the child is in */
159: DMPlexGetTreeChildren(dm,parent,&numChildren,&children);
160: for (i = 0; i < numChildren; i++) {
161: p = children[i];
162: if (p == childA) {
163: posA = i;
164: break;
165: }
166: }
167: if (posA >= coneSize) {
168: /* this is the triangle in the middle of a uniformly refined triangle: it
169: * is invariant */
170: if (dim != 2 || posA != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Expected a middle triangle, got something else");
171: *childB = childA;
172: }
173: else {
174: /* figure out position B by applying ABswapVert */
175: PetscInt posB;
177: posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize -(ABswapVert + 1) - posA) % coneSize);
178: if (childB) *childB = children[posB];
179: }
180: }
181: if (childOrientB) *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap);
182: return(0);
183: }
187: /*@
188: DMPlexReferenceTreeGetChildSymmetry - Given a reference tree, transform a childid and orientation from one parent frame to another
190: Input Parameters:
191: + dm - the reference tree DMPlex object
192: . parent - the parent point
193: . parentOrientA - the reference orientation for describing the parent
194: . childOrientA - the reference orientation for describing the child
195: . childA - the reference childID for describing the child
196: - parentOrientB - the new orientation for describing the parent
198: Output Parameters:
199: + childOrientB - if not NULL, set to the new oreintation for describing the child
200: . childB - if not NULL, the new childID for describing the child
202: Level: developer
204: .seealso: DMPlexGetReferenceTree(), DMPlexSetReferenceTree(), DMPlexSetTree()
205: @*/
206: PetscErrorCode DMPlexReferenceTreeGetChildSymmetry(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
207: {
208: DM_Plex *mesh = (DM_Plex *)dm->data;
213: if (!mesh->getchildsymmetry) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"DMPlexReferenceTreeGetChildSymmetry not implemented");
214: mesh->getchildsymmetry(dm,parent,parentOrientA,childOrientA,childA,parentOrientB,childOrientB,childB);
215: return(0);
216: }
218: static PetscErrorCode DMPlexSetTree_Internal(DM,PetscSection,PetscInt*,PetscInt*,PetscBool,PetscBool);
222: PetscErrorCode DMPlexCreateReferenceTree_Union(DM K, DM Kref, const char *labelName, DM *ref)
223: {
224: MPI_Comm comm;
225: PetscInt dim, p, pStart, pEnd, pRefStart, pRefEnd, d, offset, parentSize, *parents, *childIDs;
226: PetscInt *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts;
227: DMLabel identity, identityRef;
228: PetscSection unionSection, unionConeSection, parentSection;
229: PetscScalar *unionCoords;
230: IS perm;
234: comm = PetscObjectComm((PetscObject)K);
235: DMGetDimension(K, &dim);
236: DMPlexGetChart(K, &pStart, &pEnd);
237: DMGetLabel(K, labelName, &identity);
238: DMGetLabel(Kref, labelName, &identityRef);
239: DMPlexGetChart(Kref, &pRefStart, &pRefEnd);
240: PetscSectionCreate(comm, &unionSection);
241: PetscSectionSetChart(unionSection, 0, (pEnd - pStart) + (pRefEnd - pRefStart));
242: /* count points that will go in the union */
243: for (p = pStart; p < pEnd; p++) {
244: PetscSectionSetDof(unionSection, p - pStart, 1);
245: }
246: for (p = pRefStart; p < pRefEnd; p++) {
247: PetscInt q, qSize;
248: DMLabelGetValue(identityRef, p, &q);
249: DMLabelGetStratumSize(identityRef, q, &qSize);
250: if (qSize > 1) {
251: PetscSectionSetDof(unionSection, p - pRefStart + (pEnd - pStart), 1);
252: }
253: }
254: PetscMalloc1(pEnd - pStart + pRefEnd - pRefStart,&permvals);
255: offset = 0;
256: /* stratify points in the union by topological dimension */
257: for (d = 0; d <= dim; d++) {
258: PetscInt cStart, cEnd, c;
260: DMPlexGetHeightStratum(K, d, &cStart, &cEnd);
261: for (c = cStart; c < cEnd; c++) {
262: permvals[offset++] = c;
263: }
265: DMPlexGetHeightStratum(Kref, d, &cStart, &cEnd);
266: for (c = cStart; c < cEnd; c++) {
267: permvals[offset++] = c + (pEnd - pStart);
268: }
269: }
270: ISCreateGeneral(comm, (pEnd - pStart) + (pRefEnd - pRefStart), permvals, PETSC_OWN_POINTER, &perm);
271: PetscSectionSetPermutation(unionSection,perm);
272: PetscSectionSetUp(unionSection);
273: PetscSectionGetStorageSize(unionSection,&numUnionPoints);
274: PetscMalloc2(numUnionPoints,&coneSizes,dim+1,&numDimPoints);
275: /* count dimension points */
276: for (d = 0; d <= dim; d++) {
277: PetscInt cStart, cOff, cOff2;
278: DMPlexGetHeightStratum(K,d,&cStart,NULL);
279: PetscSectionGetOffset(unionSection,cStart-pStart,&cOff);
280: if (d < dim) {
281: DMPlexGetHeightStratum(K,d+1,&cStart,NULL);
282: PetscSectionGetOffset(unionSection,cStart-pStart,&cOff2);
283: }
284: else {
285: cOff2 = numUnionPoints;
286: }
287: numDimPoints[dim - d] = cOff2 - cOff;
288: }
289: PetscSectionCreate(comm, &unionConeSection);
290: PetscSectionSetChart(unionConeSection, 0, numUnionPoints);
291: /* count the cones in the union */
292: for (p = pStart; p < pEnd; p++) {
293: PetscInt dof, uOff;
295: DMPlexGetConeSize(K, p, &dof);
296: PetscSectionGetOffset(unionSection, p - pStart,&uOff);
297: PetscSectionSetDof(unionConeSection, uOff, dof);
298: coneSizes[uOff] = dof;
299: }
300: for (p = pRefStart; p < pRefEnd; p++) {
301: PetscInt dof, uDof, uOff;
303: DMPlexGetConeSize(Kref, p, &dof);
304: PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);
305: PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);
306: if (uDof) {
307: PetscSectionSetDof(unionConeSection, uOff, dof);
308: coneSizes[uOff] = dof;
309: }
310: }
311: PetscSectionSetUp(unionConeSection);
312: PetscSectionGetStorageSize(unionConeSection,&numCones);
313: PetscMalloc2(numCones,&unionCones,numCones,&unionOrientations);
314: /* write the cones in the union */
315: for (p = pStart; p < pEnd; p++) {
316: PetscInt dof, uOff, c, cOff;
317: const PetscInt *cone, *orientation;
319: DMPlexGetConeSize(K, p, &dof);
320: DMPlexGetCone(K, p, &cone);
321: DMPlexGetConeOrientation(K, p, &orientation);
322: PetscSectionGetOffset(unionSection, p - pStart,&uOff);
323: PetscSectionGetOffset(unionConeSection,uOff,&cOff);
324: for (c = 0; c < dof; c++) {
325: PetscInt e, eOff;
326: e = cone[c];
327: PetscSectionGetOffset(unionSection, e - pStart, &eOff);
328: unionCones[cOff + c] = eOff;
329: unionOrientations[cOff + c] = orientation[c];
330: }
331: }
332: for (p = pRefStart; p < pRefEnd; p++) {
333: PetscInt dof, uDof, uOff, c, cOff;
334: const PetscInt *cone, *orientation;
336: DMPlexGetConeSize(Kref, p, &dof);
337: DMPlexGetCone(Kref, p, &cone);
338: DMPlexGetConeOrientation(Kref, p, &orientation);
339: PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);
340: PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);
341: if (uDof) {
342: PetscSectionGetOffset(unionConeSection,uOff,&cOff);
343: for (c = 0; c < dof; c++) {
344: PetscInt e, eOff, eDof;
346: e = cone[c];
347: PetscSectionGetDof(unionSection, e - pRefStart + (pEnd - pStart),&eDof);
348: if (eDof) {
349: PetscSectionGetOffset(unionSection, e - pRefStart + (pEnd - pStart), &eOff);
350: }
351: else {
352: DMLabelGetValue(identityRef, e, &e);
353: PetscSectionGetOffset(unionSection, e - pStart, &eOff);
354: }
355: unionCones[cOff + c] = eOff;
356: unionOrientations[cOff + c] = orientation[c];
357: }
358: }
359: }
360: /* get the coordinates */
361: {
362: PetscInt vStart, vEnd, vRefStart, vRefEnd, v, vDof, vOff;
363: PetscSection KcoordsSec, KrefCoordsSec;
364: Vec KcoordsVec, KrefCoordsVec;
365: PetscScalar *Kcoords;
367: DMGetCoordinateSection(K, &KcoordsSec);
368: DMGetCoordinatesLocal(K, &KcoordsVec);
369: DMGetCoordinateSection(Kref, &KrefCoordsSec);
370: DMGetCoordinatesLocal(Kref, &KrefCoordsVec);
372: numVerts = numDimPoints[0];
373: PetscMalloc1(numVerts * dim,&unionCoords);
374: DMPlexGetDepthStratum(K,0,&vStart,&vEnd);
376: offset = 0;
377: for (v = vStart; v < vEnd; v++) {
378: PetscSectionGetOffset(unionSection,v - pStart,&vOff);
379: VecGetValuesSection(KcoordsVec, KcoordsSec, v, &Kcoords);
380: for (d = 0; d < dim; d++) {
381: unionCoords[offset * dim + d] = Kcoords[d];
382: }
383: offset++;
384: }
385: DMPlexGetDepthStratum(Kref,0,&vRefStart,&vRefEnd);
386: for (v = vRefStart; v < vRefEnd; v++) {
387: PetscSectionGetDof(unionSection,v - pRefStart + (pEnd - pStart),&vDof);
388: PetscSectionGetOffset(unionSection,v - pRefStart + (pEnd - pStart),&vOff);
389: VecGetValuesSection(KrefCoordsVec, KrefCoordsSec, v, &Kcoords);
390: if (vDof) {
391: for (d = 0; d < dim; d++) {
392: unionCoords[offset * dim + d] = Kcoords[d];
393: }
394: offset++;
395: }
396: }
397: }
398: DMCreate(comm,ref);
399: DMSetType(*ref,DMPLEX);
400: DMSetDimension(*ref,dim);
401: DMPlexCreateFromDAG(*ref,dim,numDimPoints,coneSizes,unionCones,unionOrientations,unionCoords);
402: /* set the tree */
403: PetscSectionCreate(comm,&parentSection);
404: PetscSectionSetChart(parentSection,0,numUnionPoints);
405: for (p = pRefStart; p < pRefEnd; p++) {
406: PetscInt uDof, uOff;
408: PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);
409: PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);
410: if (uDof) {
411: PetscSectionSetDof(parentSection,uOff,1);
412: }
413: }
414: PetscSectionSetUp(parentSection);
415: PetscSectionGetStorageSize(parentSection,&parentSize);
416: PetscMalloc2(parentSize,&parents,parentSize,&childIDs);
417: for (p = pRefStart; p < pRefEnd; p++) {
418: PetscInt uDof, uOff;
420: PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);
421: PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);
422: if (uDof) {
423: PetscInt pOff, parent, parentU;
424: PetscSectionGetOffset(parentSection,uOff,&pOff);
425: DMLabelGetValue(identityRef,p,&parent);
426: PetscSectionGetOffset(unionSection, parent - pStart,&parentU);
427: parents[pOff] = parentU;
428: childIDs[pOff] = uOff;
429: }
430: }
431: DMPlexSetTree_Internal(*ref,parentSection,parents,childIDs,PETSC_TRUE,PETSC_FALSE);
432: PetscSectionDestroy(&parentSection);
433: PetscFree2(parents,childIDs);
435: /* clean up */
436: PetscSectionDestroy(&unionSection);
437: PetscSectionDestroy(&unionConeSection);
438: ISDestroy(&perm);
439: PetscFree(unionCoords);
440: PetscFree2(unionCones,unionOrientations);
441: PetscFree2(coneSizes,numDimPoints);
442: return(0);
443: }
447: /*@
448: DMPlexCreateDefaultReferenceTree - create a reference tree for isotropic hierarchical mesh refinement.
450: Collective on comm
452: Input Parameters:
453: + comm - the MPI communicator
454: . dim - the spatial dimension
455: - simplex - Flag for simplex, otherwise use a tensor-product cell
457: Output Parameters:
458: . ref - the reference tree DMPlex object
460: Level: intermediate
462: .keywords: reference cell
463: .seealso: DMPlexSetReferenceTree(), DMPlexGetReferenceTree()
464: @*/
465: PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref)
466: {
467: DM_Plex *mesh;
468: DM K, Kref;
469: PetscInt p, pStart, pEnd;
470: DMLabel identity;
474: #if 1
475: comm = PETSC_COMM_SELF;
476: #endif
477: /* create a reference element */
478: DMPlexCreateReferenceCell(comm, dim, simplex, &K);
479: DMCreateLabel(K, "identity");
480: DMGetLabel(K, "identity", &identity);
481: DMPlexGetChart(K, &pStart, &pEnd);
482: for (p = pStart; p < pEnd; p++) {
483: DMLabelSetValue(identity, p, p);
484: }
485: /* refine it */
486: DMRefine(K,comm,&Kref);
488: /* the reference tree is the union of these two, without duplicating
489: * points that appear in both */
490: DMPlexCreateReferenceTree_Union(K, Kref, "identity", ref);
491: mesh = (DM_Plex *) (*ref)->data;
492: mesh->getchildsymmetry = DMPlexReferenceTreeGetChildSymmetry_Default;
493: DMDestroy(&K);
494: DMDestroy(&Kref);
495: return(0);
496: }
500: static PetscErrorCode DMPlexTreeSymmetrize(DM dm)
501: {
502: DM_Plex *mesh = (DM_Plex *)dm->data;
503: PetscSection childSec, pSec;
504: PetscInt p, pSize, cSize, parMax = PETSC_MIN_INT, parMin = PETSC_MAX_INT;
505: PetscInt *offsets, *children, pStart, pEnd;
510: PetscSectionDestroy(&mesh->childSection);
511: PetscFree(mesh->children);
512: pSec = mesh->parentSection;
513: if (!pSec) return(0);
514: PetscSectionGetStorageSize(pSec,&pSize);
515: for (p = 0; p < pSize; p++) {
516: PetscInt par = mesh->parents[p];
518: parMax = PetscMax(parMax,par+1);
519: parMin = PetscMin(parMin,par);
520: }
521: if (parMin > parMax) {
522: parMin = -1;
523: parMax = -1;
524: }
525: PetscSectionCreate(PetscObjectComm((PetscObject)pSec),&childSec);
526: PetscSectionSetChart(childSec,parMin,parMax);
527: for (p = 0; p < pSize; p++) {
528: PetscInt par = mesh->parents[p];
530: PetscSectionAddDof(childSec,par,1);
531: }
532: PetscSectionSetUp(childSec);
533: PetscSectionGetStorageSize(childSec,&cSize);
534: PetscMalloc1(cSize,&children);
535: PetscCalloc1(parMax-parMin,&offsets);
536: PetscSectionGetChart(pSec,&pStart,&pEnd);
537: for (p = pStart; p < pEnd; p++) {
538: PetscInt dof, off, i;
540: PetscSectionGetDof(pSec,p,&dof);
541: PetscSectionGetOffset(pSec,p,&off);
542: for (i = 0; i < dof; i++) {
543: PetscInt par = mesh->parents[off + i], cOff;
545: PetscSectionGetOffset(childSec,par,&cOff);
546: children[cOff + offsets[par-parMin]++] = p;
547: }
548: }
549: mesh->childSection = childSec;
550: mesh->children = children;
551: PetscFree(offsets);
552: return(0);
553: }
557: static PetscErrorCode AnchorsFlatten (PetscSection section, IS is, PetscSection *sectionNew, IS *isNew)
558: {
559: PetscInt pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL;
560: const PetscInt *vals;
561: PetscSection secNew;
562: PetscBool anyNew, globalAnyNew;
563: PetscBool compress;
567: PetscSectionGetChart(section,&pStart,&pEnd);
568: ISGetLocalSize(is,&size);
569: ISGetIndices(is,&vals);
570: PetscSectionCreate(PetscObjectComm((PetscObject)section),&secNew);
571: PetscSectionSetChart(secNew,pStart,pEnd);
572: for (i = 0; i < size; i++) {
573: PetscInt dof;
575: p = vals[i];
576: if (p < pStart || p >= pEnd) continue;
577: PetscSectionGetDof(section, p, &dof);
578: if (dof) break;
579: }
580: if (i == size) {
581: PetscSectionSetUp(secNew);
582: anyNew = PETSC_FALSE;
583: compress = PETSC_FALSE;
584: sizeNew = 0;
585: }
586: else {
587: anyNew = PETSC_TRUE;
588: for (p = pStart; p < pEnd; p++) {
589: PetscInt dof, off;
591: PetscSectionGetDof(section, p, &dof);
592: PetscSectionGetOffset(section, p, &off);
593: for (i = 0; i < dof; i++) {
594: PetscInt q = vals[off + i], qDof = 0;
596: if (q >= pStart && q < pEnd) {
597: PetscSectionGetDof(section, q, &qDof);
598: }
599: if (qDof) {
600: PetscSectionAddDof(secNew, p, qDof);
601: }
602: else {
603: PetscSectionAddDof(secNew, p, 1);
604: }
605: }
606: }
607: PetscSectionSetUp(secNew);
608: PetscSectionGetStorageSize(secNew,&sizeNew);
609: PetscMalloc1(sizeNew,&valsNew);
610: compress = PETSC_FALSE;
611: for (p = pStart; p < pEnd; p++) {
612: PetscInt dof, off, count, offNew, dofNew;
614: PetscSectionGetDof(section, p, &dof);
615: PetscSectionGetOffset(section, p, &off);
616: PetscSectionGetDof(secNew, p, &dofNew);
617: PetscSectionGetOffset(secNew, p, &offNew);
618: count = 0;
619: for (i = 0; i < dof; i++) {
620: PetscInt q = vals[off + i], qDof = 0, qOff = 0, j;
622: if (q >= pStart && q < pEnd) {
623: PetscSectionGetDof(section, q, &qDof);
624: PetscSectionGetOffset(section, q, &qOff);
625: }
626: if (qDof) {
627: PetscInt oldCount = count;
629: for (j = 0; j < qDof; j++) {
630: PetscInt k, r = vals[qOff + j];
632: for (k = 0; k < oldCount; k++) {
633: if (valsNew[offNew + k] == r) {
634: break;
635: }
636: }
637: if (k == oldCount) {
638: valsNew[offNew + count++] = r;
639: }
640: }
641: }
642: else {
643: PetscInt k, oldCount = count;
645: for (k = 0; k < oldCount; k++) {
646: if (valsNew[offNew + k] == q) {
647: break;
648: }
649: }
650: if (k == oldCount) {
651: valsNew[offNew + count++] = q;
652: }
653: }
654: }
655: if (count < dofNew) {
656: PetscSectionSetDof(secNew, p, count);
657: compress = PETSC_TRUE;
658: }
659: }
660: }
661: ISRestoreIndices(is,&vals);
662: MPIU_Allreduce(&anyNew,&globalAnyNew,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));
663: if (!globalAnyNew) {
664: PetscSectionDestroy(&secNew);
665: *sectionNew = NULL;
666: *isNew = NULL;
667: }
668: else {
669: PetscBool globalCompress;
671: MPIU_Allreduce(&compress,&globalCompress,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));
672: if (compress) {
673: PetscSection secComp;
674: PetscInt *valsComp = NULL;
676: PetscSectionCreate(PetscObjectComm((PetscObject)section),&secComp);
677: PetscSectionSetChart(secComp,pStart,pEnd);
678: for (p = pStart; p < pEnd; p++) {
679: PetscInt dof;
681: PetscSectionGetDof(secNew, p, &dof);
682: PetscSectionSetDof(secComp, p, dof);
683: }
684: PetscSectionSetUp(secComp);
685: PetscSectionGetStorageSize(secComp,&sizeNew);
686: PetscMalloc1(sizeNew,&valsComp);
687: for (p = pStart; p < pEnd; p++) {
688: PetscInt dof, off, offNew, j;
690: PetscSectionGetDof(secNew, p, &dof);
691: PetscSectionGetOffset(secNew, p, &off);
692: PetscSectionGetOffset(secComp, p, &offNew);
693: for (j = 0; j < dof; j++) {
694: valsComp[offNew + j] = valsNew[off + j];
695: }
696: }
697: PetscSectionDestroy(&secNew);
698: secNew = secComp;
699: PetscFree(valsNew);
700: valsNew = valsComp;
701: }
702: ISCreateGeneral(PetscObjectComm((PetscObject)is),sizeNew,valsNew,PETSC_OWN_POINTER,isNew);
703: }
704: return(0);
705: }
709: static PetscErrorCode DMPlexCreateAnchors_Tree(DM dm)
710: {
711: PetscInt p, pStart, pEnd, *anchors, size;
712: PetscInt aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT;
713: PetscSection aSec;
714: DMLabel canonLabel;
715: IS aIS;
720: DMPlexGetChart(dm,&pStart,&pEnd);
721: DMGetLabel(dm,"canonical",&canonLabel);
722: for (p = pStart; p < pEnd; p++) {
723: PetscInt parent;
725: if (canonLabel) {
726: PetscInt canon;
728: DMLabelGetValue(canonLabel,p,&canon);
729: if (p != canon) continue;
730: }
731: DMPlexGetTreeParent(dm,p,&parent,NULL);
732: if (parent != p) {
733: aMin = PetscMin(aMin,p);
734: aMax = PetscMax(aMax,p+1);
735: }
736: }
737: if (aMin > aMax) {
738: aMin = -1;
739: aMax = -1;
740: }
741: PetscSectionCreate(PETSC_COMM_SELF,&aSec);
742: PetscSectionSetChart(aSec,aMin,aMax);
743: for (p = aMin; p < aMax; p++) {
744: PetscInt parent, ancestor = p;
746: if (canonLabel) {
747: PetscInt canon;
749: DMLabelGetValue(canonLabel,p,&canon);
750: if (p != canon) continue;
751: }
752: DMPlexGetTreeParent(dm,p,&parent,NULL);
753: while (parent != ancestor) {
754: ancestor = parent;
755: DMPlexGetTreeParent(dm,ancestor,&parent,NULL);
756: }
757: if (ancestor != p) {
758: PetscInt closureSize, *closure = NULL;
760: DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
761: PetscSectionSetDof(aSec,p,closureSize);
762: DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
763: }
764: }
765: PetscSectionSetUp(aSec);
766: PetscSectionGetStorageSize(aSec,&size);
767: PetscMalloc1(size,&anchors);
768: for (p = aMin; p < aMax; p++) {
769: PetscInt parent, ancestor = p;
771: if (canonLabel) {
772: PetscInt canon;
774: DMLabelGetValue(canonLabel,p,&canon);
775: if (p != canon) continue;
776: }
777: DMPlexGetTreeParent(dm,p,&parent,NULL);
778: while (parent != ancestor) {
779: ancestor = parent;
780: DMPlexGetTreeParent(dm,ancestor,&parent,NULL);
781: }
782: if (ancestor != p) {
783: PetscInt j, closureSize, *closure = NULL, aOff;
785: PetscSectionGetOffset(aSec,p,&aOff);
787: DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
788: for (j = 0; j < closureSize; j++) {
789: anchors[aOff + j] = closure[2*j];
790: }
791: DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
792: }
793: }
794: ISCreateGeneral(PETSC_COMM_SELF,size,anchors,PETSC_OWN_POINTER,&aIS);
795: {
796: PetscSection aSecNew = aSec;
797: IS aISNew = aIS;
799: PetscObjectReference((PetscObject)aSec);
800: PetscObjectReference((PetscObject)aIS);
801: while (aSecNew) {
802: PetscSectionDestroy(&aSec);
803: ISDestroy(&aIS);
804: aSec = aSecNew;
805: aIS = aISNew;
806: aSecNew = NULL;
807: aISNew = NULL;
808: AnchorsFlatten(aSec,aIS,&aSecNew,&aISNew);
809: }
810: }
811: DMPlexSetAnchors(dm,aSec,aIS);
812: PetscSectionDestroy(&aSec);
813: ISDestroy(&aIS);
814: return(0);
815: }
819: static PetscErrorCode DMPlexGetTrueSupportSize(DM dm,PetscInt p,PetscInt *dof,PetscInt *numTrueSupp)
820: {
824: if (numTrueSupp[p] == -1) {
825: PetscInt i, alldof;
826: const PetscInt *supp;
827: PetscInt count = 0;
829: DMPlexGetSupportSize(dm,p,&alldof);
830: DMPlexGetSupport(dm,p,&supp);
831: for (i = 0; i < alldof; i++) {
832: PetscInt q = supp[i], numCones, j;
833: const PetscInt *cone;
835: DMPlexGetConeSize(dm,q,&numCones);
836: DMPlexGetCone(dm,q,&cone);
837: for (j = 0; j < numCones; j++) {
838: if (cone[j] == p) break;
839: }
840: if (j < numCones) count++;
841: }
842: numTrueSupp[p] = count;
843: }
844: *dof = numTrueSupp[p];
845: return(0);
846: }
850: static PetscErrorCode DMPlexTreeExchangeSupports(DM dm)
851: {
852: DM_Plex *mesh = (DM_Plex *)dm->data;
853: PetscSection newSupportSection;
854: PetscInt newSize, *newSupports, pStart, pEnd, p, d, depth;
855: PetscInt *numTrueSupp;
856: PetscInt *offsets;
861: /* symmetrize the hierarchy */
862: DMPlexGetDepth(dm,&depth);
863: PetscSectionCreate(PetscObjectComm((PetscObject)(mesh->supportSection)),&newSupportSection);
864: DMPlexGetChart(dm,&pStart,&pEnd);
865: PetscSectionSetChart(newSupportSection,pStart,pEnd);
866: PetscCalloc1(pEnd,&offsets);
867: PetscMalloc1(pEnd,&numTrueSupp);
868: for (p = 0; p < pEnd; p++) numTrueSupp[p] = -1;
869: /* if a point is in the (true) support of q, it should be in the support of
870: * parent(q) */
871: for (d = 0; d <= depth; d++) {
872: DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);
873: for (p = pStart; p < pEnd; ++p) {
874: PetscInt dof, q, qdof, parent;
876: DMPlexGetTrueSupportSize(dm,p,&dof,numTrueSupp);
877: PetscSectionAddDof(newSupportSection, p, dof);
878: q = p;
879: DMPlexGetTreeParent(dm,q,&parent,NULL);
880: while (parent != q && parent >= pStart && parent < pEnd) {
881: q = parent;
883: DMPlexGetTrueSupportSize(dm,q,&qdof,numTrueSupp);
884: PetscSectionAddDof(newSupportSection,p,qdof);
885: PetscSectionAddDof(newSupportSection,q,dof);
886: DMPlexGetTreeParent(dm,q,&parent,NULL);
887: }
888: }
889: }
890: PetscSectionSetUp(newSupportSection);
891: PetscSectionGetStorageSize(newSupportSection,&newSize);
892: PetscMalloc1(newSize,&newSupports);
893: for (d = 0; d <= depth; d++) {
894: DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);
895: for (p = pStart; p < pEnd; p++) {
896: PetscInt dof, off, q, qdof, qoff, newDof, newOff, newqOff, i, parent;
898: PetscSectionGetDof(mesh->supportSection, p, &dof);
899: PetscSectionGetOffset(mesh->supportSection, p, &off);
900: PetscSectionGetDof(newSupportSection, p, &newDof);
901: PetscSectionGetOffset(newSupportSection, p, &newOff);
902: for (i = 0; i < dof; i++) {
903: PetscInt numCones, j;
904: const PetscInt *cone;
905: PetscInt q = mesh->supports[off + i];
907: DMPlexGetConeSize(dm,q,&numCones);
908: DMPlexGetCone(dm,q,&cone);
909: for (j = 0; j < numCones; j++) {
910: if (cone[j] == p) break;
911: }
912: if (j < numCones) newSupports[newOff+offsets[p]++] = q;
913: }
914: mesh->maxSupportSize = PetscMax(mesh->maxSupportSize,newDof);
916: q = p;
917: DMPlexGetTreeParent(dm,q,&parent,NULL);
918: while (parent != q && parent >= pStart && parent < pEnd) {
919: q = parent;
920: PetscSectionGetDof(mesh->supportSection, q, &qdof);
921: PetscSectionGetOffset(mesh->supportSection, q, &qoff);
922: PetscSectionGetOffset(newSupportSection, q, &newqOff);
923: for (i = 0; i < qdof; i++) {
924: PetscInt numCones, j;
925: const PetscInt *cone;
926: PetscInt r = mesh->supports[qoff + i];
928: DMPlexGetConeSize(dm,r,&numCones);
929: DMPlexGetCone(dm,r,&cone);
930: for (j = 0; j < numCones; j++) {
931: if (cone[j] == q) break;
932: }
933: if (j < numCones) newSupports[newOff+offsets[p]++] = r;
934: }
935: for (i = 0; i < dof; i++) {
936: PetscInt numCones, j;
937: const PetscInt *cone;
938: PetscInt r = mesh->supports[off + i];
940: DMPlexGetConeSize(dm,r,&numCones);
941: DMPlexGetCone(dm,r,&cone);
942: for (j = 0; j < numCones; j++) {
943: if (cone[j] == p) break;
944: }
945: if (j < numCones) newSupports[newqOff+offsets[q]++] = r;
946: }
947: DMPlexGetTreeParent(dm,q,&parent,NULL);
948: }
949: }
950: }
951: PetscSectionDestroy(&mesh->supportSection);
952: mesh->supportSection = newSupportSection;
953: PetscFree(mesh->supports);
954: mesh->supports = newSupports;
955: PetscFree(offsets);
956: PetscFree(numTrueSupp);
958: return(0);
959: }
961: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM,PetscSection,PetscSection,Mat);
962: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM,PetscSection,PetscSection,Mat);
966: static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical, PetscBool exchangeSupports)
967: {
968: DM_Plex *mesh = (DM_Plex *)dm->data;
969: DM refTree;
970: PetscInt size;
976: PetscObjectReference((PetscObject)parentSection);
977: PetscSectionDestroy(&mesh->parentSection);
978: mesh->parentSection = parentSection;
979: PetscSectionGetStorageSize(parentSection,&size);
980: if (parents != mesh->parents) {
981: PetscFree(mesh->parents);
982: PetscMalloc1(size,&mesh->parents);
983: PetscMemcpy(mesh->parents, parents, size * sizeof(*parents));
984: }
985: if (childIDs != mesh->childIDs) {
986: PetscFree(mesh->childIDs);
987: PetscMalloc1(size,&mesh->childIDs);
988: PetscMemcpy(mesh->childIDs, childIDs, size * sizeof(*childIDs));
989: }
990: DMPlexGetReferenceTree(dm,&refTree);
991: if (refTree) {
992: DMLabel canonLabel;
994: DMGetLabel(refTree,"canonical",&canonLabel);
995: if (canonLabel) {
996: PetscInt i;
998: for (i = 0; i < size; i++) {
999: PetscInt canon;
1000: DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon);
1001: if (canon >= 0) {
1002: mesh->childIDs[i] = canon;
1003: }
1004: }
1005: }
1006: mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_FromReference;
1007: }
1008: else {
1009: mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_Direct;
1010: }
1011: DMPlexTreeSymmetrize(dm);
1012: if (computeCanonical) {
1013: PetscInt d, dim;
1015: /* add the canonical label */
1016: DMGetDimension(dm,&dim);
1017: DMCreateLabel(dm,"canonical");
1018: for (d = 0; d <= dim; d++) {
1019: PetscInt p, dStart, dEnd, canon = -1, cNumChildren;
1020: const PetscInt *cChildren;
1022: DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);
1023: for (p = dStart; p < dEnd; p++) {
1024: DMPlexGetTreeChildren(dm,p,&cNumChildren,&cChildren);
1025: if (cNumChildren) {
1026: canon = p;
1027: break;
1028: }
1029: }
1030: if (canon == -1) continue;
1031: for (p = dStart; p < dEnd; p++) {
1032: PetscInt numChildren, i;
1033: const PetscInt *children;
1035: DMPlexGetTreeChildren(dm,p,&numChildren,&children);
1036: if (numChildren) {
1037: 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);
1038: DMSetLabelValue(dm,"canonical",p,canon);
1039: for (i = 0; i < numChildren; i++) {
1040: DMSetLabelValue(dm,"canonical",children[i],cChildren[i]);
1041: }
1042: }
1043: }
1044: }
1045: }
1046: if (exchangeSupports) {
1047: DMPlexTreeExchangeSupports(dm);
1048: }
1049: mesh->createanchors = DMPlexCreateAnchors_Tree;
1050: /* reset anchors */
1051: DMPlexSetAnchors(dm,NULL,NULL);
1052: return(0);
1053: }
1057: /*@
1058: DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points. This routine also creates
1059: the point-to-point constraints determined by the tree: a point is constained to the points in the closure of its
1060: tree root.
1062: Collective on dm
1064: Input Parameters:
1065: + dm - the DMPlex object
1066: . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
1067: offset indexes the parent and childID list; the reference count of parentSection is incremented
1068: . parents - a list of the point parents; copied, can be destroyed
1069: - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
1070: the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed
1072: Level: intermediate
1074: .seealso: DMPlexGetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren()
1075: @*/
1076: PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
1077: {
1081: DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_FALSE,PETSC_TRUE);
1082: return(0);
1083: }
1087: /*@
1088: DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points.
1089: Collective on dm
1091: Input Parameters:
1092: . dm - the DMPlex object
1094: Output Parameters:
1095: + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
1096: offset indexes the parent and childID list
1097: . parents - a list of the point parents
1098: . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
1099: the child corresponds to the point in the reference tree with index childID
1100: . childSection - the inverse of the parent section
1101: - children - a list of the point children
1103: Level: intermediate
1105: .seealso: DMPlexSetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren()
1106: @*/
1107: PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[])
1108: {
1109: DM_Plex *mesh = (DM_Plex *)dm->data;
1113: if (parentSection) *parentSection = mesh->parentSection;
1114: if (parents) *parents = mesh->parents;
1115: if (childIDs) *childIDs = mesh->childIDs;
1116: if (childSection) *childSection = mesh->childSection;
1117: if (children) *children = mesh->children;
1118: return(0);
1119: }
1123: /*@
1124: DMPlexGetTreeParent - get the parent of a point in the tree describing the point hierarchy (not the Sieve DAG)
1126: Input Parameters:
1127: + dm - the DMPlex object
1128: - point - the query point
1130: Output Parameters:
1131: + parent - if not NULL, set to the parent of the point, or the point itself if the point does not have a parent
1132: - childID - if not NULL, set to the child ID of the point with respect to its parent, or 0 if the point
1133: does not have a parent
1135: Level: intermediate
1137: .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeChildren()
1138: @*/
1139: PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID)
1140: {
1141: DM_Plex *mesh = (DM_Plex *)dm->data;
1142: PetscSection pSec;
1147: pSec = mesh->parentSection;
1148: if (pSec && point >= pSec->pStart && point < pSec->pEnd) {
1149: PetscInt dof;
1151: PetscSectionGetDof (pSec, point, &dof);
1152: if (dof) {
1153: PetscInt off;
1155: PetscSectionGetOffset (pSec, point, &off);
1156: if (parent) *parent = mesh->parents[off];
1157: if (childID) *childID = mesh->childIDs[off];
1158: return(0);
1159: }
1160: }
1161: if (parent) {
1162: *parent = point;
1163: }
1164: if (childID) {
1165: *childID = 0;
1166: }
1167: return(0);
1168: }
1172: /*@C
1173: DMPlexGetTreeChildren - get the children of a point in the tree describing the point hierarchy (not the Sieve DAG)
1175: Input Parameters:
1176: + dm - the DMPlex object
1177: - point - the query point
1179: Output Parameters:
1180: + numChildren - if not NULL, set to the number of children
1181: - children - if not NULL, set to a list children, or set to NULL if the point has no children
1183: Level: intermediate
1185: Fortran Notes:
1186: Since it returns an array, this routine is only available in Fortran 90, and you must
1187: include petsc.h90 in your code.
1189: .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeParent()
1190: @*/
1191: PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[])
1192: {
1193: DM_Plex *mesh = (DM_Plex *)dm->data;
1194: PetscSection childSec;
1195: PetscInt dof = 0;
1200: childSec = mesh->childSection;
1201: if (childSec && point >= childSec->pStart && point < childSec->pEnd) {
1202: PetscSectionGetDof (childSec, point, &dof);
1203: }
1204: if (numChildren) *numChildren = dof;
1205: if (children) {
1206: if (dof) {
1207: PetscInt off;
1209: PetscSectionGetOffset (childSec, point, &off);
1210: *children = &mesh->children[off];
1211: }
1212: else {
1213: *children = NULL;
1214: }
1215: }
1216: return(0);
1217: }
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: PetscFE fe = NULL;
1246: PetscFV fv = NULL;
1247: PetscClassId id;
1248: PetscDualSpace space;
1249: PetscInt i, j, k, nPoints, offset;
1250: PetscInt fSize, fComp;
1251: PetscReal *B = NULL;
1252: PetscReal *weights, *pointsRef, *pointsReal;
1253: Mat Amat, Bmat, Xmat;
1255: PetscDSGetDiscretization(ds,f,&disc);
1256: PetscObjectGetClassId(disc,&id);
1257: if (id == PETSCFE_CLASSID) {
1258: fe = (PetscFE) disc;
1259: PetscFEGetDualSpace(fe,&space);
1260: PetscDualSpaceGetDimension(space,&fSize);
1261: PetscFEGetNumComponents(fe,&fComp);
1262: }
1263: else if (id == PETSCFV_CLASSID) {
1264: fv = (PetscFV) disc;
1265: PetscFVGetDualSpace(fv,&space);
1266: PetscDualSpaceGetDimension(space,&fSize);
1267: PetscFVGetNumComponents(fv,&fComp);
1268: }
1269: else SETERRQ1(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id);
1271: MatCreate(PETSC_COMM_SELF,&Amat);
1272: MatSetSizes(Amat,fSize,fSize,fSize,fSize);
1273: MatSetType(Amat,MATSEQDENSE);
1274: MatSetUp(Amat);
1275: MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Bmat);
1276: MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Xmat);
1277: nPoints = 0;
1278: for (i = 0; i < fSize; i++) {
1279: PetscInt qPoints;
1280: PetscQuadrature quad;
1282: PetscDualSpaceGetFunctional(space,i,&quad);
1283: PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);
1284: nPoints += qPoints;
1285: }
1286: PetscMalloc3(nPoints,&weights,spdim*nPoints,&pointsRef,spdim*nPoints,&pointsReal);
1287: offset = 0;
1288: for (i = 0; i < fSize; i++) {
1289: PetscInt qPoints;
1290: const PetscReal *p, *w;
1291: PetscQuadrature quad;
1293: PetscDualSpaceGetFunctional(space,i,&quad);
1294: PetscQuadratureGetData(quad,NULL,&qPoints,&p,&w);
1295: PetscMemcpy(weights+offset,w,qPoints*sizeof(*w));
1296: PetscMemcpy(pointsRef+spdim*offset,p,spdim*qPoints*sizeof(*p));
1297: offset += qPoints;
1298: }
1299: if (id == PETSCFE_CLASSID) {
1300: PetscFEGetTabulation(fe,nPoints,pointsRef,&B,NULL,NULL);
1301: }
1302: else {
1303: PetscFVGetTabulation(fv,nPoints,pointsRef,&B,NULL,NULL);
1304: }
1305: offset = 0;
1306: for (i = 0; i < fSize; i++) {
1307: PetscInt qPoints;
1308: PetscQuadrature quad;
1310: PetscDualSpaceGetFunctional(space,i,&quad);
1311: PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);
1312: for (j = 0; j < fSize; j++) {
1313: PetscScalar val = 0.;
1315: for (k = 0; k < qPoints; k++) {
1316: val += B[((offset + k) * fSize + j) * fComp] * weights[k];
1317: }
1318: MatSetValue(Amat,i,j,val,INSERT_VALUES);
1319: }
1320: offset += qPoints;
1321: }
1322: if (id == PETSCFE_CLASSID) {
1323: PetscFERestoreTabulation(fe,nPoints,pointsRef,&B,NULL,NULL);
1324: }
1325: else {
1326: PetscFVRestoreTabulation(fv,nPoints,pointsRef,&B,NULL,NULL);
1327: }
1328: MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
1329: MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
1330: MatLUFactor(Amat,NULL,NULL,NULL);
1331: for (c = cStart; c < cEnd; c++) {
1332: PetscInt parent;
1333: PetscInt closureSize, closureSizeP, *closure = NULL, *closureP = NULL;
1334: PetscInt *childOffsets, *parentOffsets;
1336: DMPlexGetTreeParent(dm,c,&parent,NULL);
1337: if (parent == c) continue;
1338: DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1339: for (i = 0; i < closureSize; i++) {
1340: PetscInt p = closure[2*i];
1341: PetscInt conDof;
1343: if (p < conStart || p >= conEnd) continue;
1344: if (numFields > 1) {
1345: PetscSectionGetFieldDof(cSec,p,f,&conDof);
1346: }
1347: else {
1348: PetscSectionGetDof(cSec,p,&conDof);
1349: }
1350: if (conDof) break;
1351: }
1352: if (i == closureSize) {
1353: DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1354: continue;
1355: }
1357: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, NULL, &detJ);
1358: DMPlexComputeCellGeometryFEM(dm, parent, NULL, v0parent, Jparent, invJparent, &detJparent);
1359: for (i = 0; i < nPoints; i++) {
1360: CoordinatesRefToReal(spdim, spdim, v0, J, &pointsRef[i*spdim],vtmp);
1361: CoordinatesRealToRef(spdim, spdim, v0parent, invJparent, vtmp, &pointsReal[i*spdim]);
1362: }
1363: if (id == PETSCFE_CLASSID) {
1364: PetscFEGetTabulation(fe,nPoints,pointsReal,&B,NULL,NULL);
1365: }
1366: else {
1367: PetscFVGetTabulation(fv,nPoints,pointsReal,&B,NULL,NULL);
1368: }
1369: offset = 0;
1370: for (i = 0; i < fSize; i++) {
1371: PetscInt qPoints;
1372: PetscQuadrature quad;
1374: PetscDualSpaceGetFunctional(space,i,&quad);
1375: PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);
1376: for (j = 0; j < fSize; j++) {
1377: PetscScalar val = 0.;
1379: for (k = 0; k < qPoints; k++) {
1380: val += B[((offset + k) * fSize + j) * fComp] * weights[k];
1381: }
1382: MatSetValue(Bmat,i,j,val,INSERT_VALUES);
1383: }
1384: offset += qPoints;
1385: }
1386: if (id == PETSCFE_CLASSID) {
1387: PetscFERestoreTabulation(fe,nPoints,pointsReal,&B,NULL,NULL);
1388: }
1389: else {
1390: PetscFVRestoreTabulation(fv,nPoints,pointsReal,&B,NULL,NULL);
1391: }
1392: MatAssemblyBegin(Bmat,MAT_FINAL_ASSEMBLY);
1393: MatAssemblyEnd(Bmat,MAT_FINAL_ASSEMBLY);
1394: MatMatSolve(Amat,Bmat,Xmat);
1395: DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);
1396: PetscMalloc2(closureSize+1,&childOffsets,closureSizeP+1,&parentOffsets);
1397: childOffsets[0] = 0;
1398: for (i = 0; i < closureSize; i++) {
1399: PetscInt p = closure[2*i];
1400: PetscInt dof;
1402: if (numFields > 1) {
1403: PetscSectionGetFieldDof(section,p,f,&dof);
1404: }
1405: else {
1406: PetscSectionGetDof(section,p,&dof);
1407: }
1408: childOffsets[i+1]=childOffsets[i]+dof / fComp;
1409: }
1410: parentOffsets[0] = 0;
1411: for (i = 0; i < closureSizeP; i++) {
1412: PetscInt p = closureP[2*i];
1413: PetscInt dof;
1415: if (numFields > 1) {
1416: PetscSectionGetFieldDof(section,p,f,&dof);
1417: }
1418: else {
1419: PetscSectionGetDof(section,p,&dof);
1420: }
1421: parentOffsets[i+1]=parentOffsets[i]+dof / fComp;
1422: }
1423: for (i = 0; i < closureSize; i++) {
1424: PetscInt conDof, conOff, aDof, aOff;
1425: PetscInt p = closure[2*i];
1426: PetscInt o = closure[2*i+1];
1428: if (p < conStart || p >= conEnd) continue;
1429: if (numFields > 1) {
1430: PetscSectionGetFieldDof(cSec,p,f,&conDof);
1431: PetscSectionGetFieldOffset(cSec,p,f,&conOff);
1432: }
1433: else {
1434: PetscSectionGetDof(cSec,p,&conDof);
1435: PetscSectionGetOffset(cSec,p,&conOff);
1436: }
1437: if (!conDof) continue;
1438: PetscSectionGetDof(aSec,p,&aDof);
1439: PetscSectionGetOffset(aSec,p,&aOff);
1440: for (k = 0; k < aDof; k++) {
1441: PetscInt a = anchors[aOff + k];
1442: PetscInt aSecDof, aSecOff;
1444: if (numFields > 1) {
1445: PetscSectionGetFieldDof(section,a,f,&aSecDof);
1446: PetscSectionGetFieldOffset(section,a,f,&aSecOff);
1447: }
1448: else {
1449: PetscSectionGetDof(section,a,&aSecDof);
1450: PetscSectionGetOffset(section,a,&aSecOff);
1451: }
1452: if (!aSecDof) continue;
1454: for (j = 0; j < closureSizeP; j++) {
1455: PetscInt q = closureP[2*j];
1456: PetscInt oq = closureP[2*j+1];
1458: if (q == a) {
1459: PetscInt r, s, t;
1461: for (r = childOffsets[i]; r < childOffsets[i+1]; r++) {
1462: for (s = parentOffsets[j]; s < parentOffsets[j+1]; s++) {
1463: PetscScalar val;
1464: PetscInt insertCol, insertRow;
1466: MatGetValue(Xmat,r,s,&val);
1467: if (o >= 0) {
1468: insertRow = conOff + fComp * (r - childOffsets[i]);
1469: }
1470: else {
1471: insertRow = conOff + fComp * (childOffsets[i + 1] - 1 - r);
1472: }
1473: if (oq >= 0) {
1474: insertCol = aSecOff + fComp * (s - parentOffsets[j]);
1475: }
1476: else {
1477: insertCol = aSecOff + fComp * (parentOffsets[j + 1] - 1 - s);
1478: }
1479: for (t = 0; t < fComp; t++) {
1480: MatSetValue(cMat,insertRow + t,insertCol + t,val,INSERT_VALUES);
1481: }
1482: }
1483: }
1484: }
1485: }
1486: }
1487: }
1488: PetscFree2(childOffsets,parentOffsets);
1489: DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1490: DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);
1491: }
1492: MatDestroy(&Amat);
1493: MatDestroy(&Bmat);
1494: MatDestroy(&Xmat);
1495: PetscFree3(weights,pointsRef,pointsReal);
1496: }
1497: MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
1498: MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
1499: PetscFree6(v0,v0parent,vtmp,J,Jparent,invJparent);
1500: ISRestoreIndices(aIS,&anchors);
1502: return(0);
1503: }
1507: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1508: {
1509: Mat refCmat;
1510: PetscDS ds;
1511: PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, **refPointFieldN;
1512: PetscScalar ***refPointFieldMats;
1513: PetscSection refConSec, refAnSec, refSection;
1514: IS refAnIS;
1515: const PetscInt *refAnchors;
1519: DMGetDS(refTree,&ds);
1520: PetscDSGetNumFields(ds,&numFields);
1521: DMGetDefaultConstraints(refTree,&refConSec,&refCmat);
1522: DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);
1523: ISGetIndices(refAnIS,&refAnchors);
1524: DMGetDefaultSection(refTree,&refSection);
1525: PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1526: PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);
1527: PetscMalloc1(pRefEnd-pRefStart,&refPointFieldN);
1528: PetscSectionGetMaxDof(refConSec,&maxDof);
1529: PetscSectionGetMaxDof(refAnSec,&maxAnDof);
1530: PetscMalloc1(maxDof,&rows);
1531: PetscMalloc1(maxDof*maxAnDof,&cols);
1532: for (p = pRefStart; p < pRefEnd; p++) {
1533: PetscInt parent, closureSize, *closure = NULL, pDof;
1535: DMPlexGetTreeParent(refTree,p,&parent,NULL);
1536: PetscSectionGetDof(refConSec,p,&pDof);
1537: if (!pDof || parent == p) continue;
1539: PetscMalloc1(numFields,&refPointFieldMats[p-pRefStart]);
1540: PetscCalloc1(numFields,&refPointFieldN[p-pRefStart]);
1541: DMPlexGetTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);
1542: for (f = 0; f < numFields; f++) {
1543: PetscInt cDof, cOff, numCols, r, i, fComp;
1544: PetscObject disc;
1545: PetscClassId id;
1546: PetscFE fe = NULL;
1547: PetscFV fv = NULL;
1549: PetscDSGetDiscretization(ds,f,&disc);
1550: PetscObjectGetClassId(disc,&id);
1551: if (id == PETSCFE_CLASSID) {
1552: fe = (PetscFE) disc;
1553: PetscFEGetNumComponents(fe,&fComp);
1554: }
1555: else if (id == PETSCFV_CLASSID) {
1556: fv = (PetscFV) disc;
1557: PetscFVGetNumComponents(fv,&fComp);
1558: }
1559: else SETERRQ1(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id);
1561: if (numFields > 1) {
1562: PetscSectionGetFieldDof(refConSec,p,f,&cDof);
1563: PetscSectionGetFieldOffset(refConSec,p,f,&cOff);
1564: }
1565: else {
1566: PetscSectionGetDof(refConSec,p,&cDof);
1567: PetscSectionGetOffset(refConSec,p,&cOff);
1568: }
1570: if (!cDof) continue;
1571: for (r = 0; r < cDof; r++) {
1572: rows[r] = cOff + r;
1573: }
1574: numCols = 0;
1575: for (i = 0; i < closureSize; i++) {
1576: PetscInt q = closure[2*i];
1577: PetscInt o = closure[2*i+1];
1578: PetscInt aDof, aOff, j;
1580: if (numFields > 1) {
1581: PetscSectionGetFieldDof(refSection,q,f,&aDof);
1582: PetscSectionGetFieldOffset(refSection,q,f,&aOff);
1583: }
1584: else {
1585: PetscSectionGetDof(refSection,q,&aDof);
1586: PetscSectionGetOffset(refSection,q,&aOff);
1587: }
1589: for (j = 0; j < aDof; j++) {
1590: PetscInt node = (o >= 0) ? (j / fComp) : ((aDof - 1 - j) / fComp);
1591: PetscInt comp = (j % fComp);
1593: cols[numCols++] = aOff + node * fComp + comp;
1594: }
1595: }
1596: refPointFieldN[p-pRefStart][f] = numCols;
1597: PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);
1598: MatGetValues(refCmat,cDof,rows,numCols,cols,refPointFieldMats[p-pRefStart][f]);
1599: }
1600: DMPlexRestoreTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);
1601: }
1602: *childrenMats = refPointFieldMats;
1603: *childrenN = refPointFieldN;
1604: ISRestoreIndices(refAnIS,&refAnchors);
1605: PetscFree(rows);
1606: PetscFree(cols);
1607: return(0);
1608: }
1612: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN)
1613: {
1614: PetscDS ds;
1615: PetscInt **refPointFieldN;
1616: PetscScalar ***refPointFieldMats;
1617: PetscInt numFields, pRefStart, pRefEnd, p, f;
1618: PetscSection refConSec;
1622: refPointFieldN = *childrenN;
1623: *childrenN = NULL;
1624: refPointFieldMats = *childrenMats;
1625: *childrenMats = NULL;
1626: DMGetDS(refTree,&ds);
1627: PetscDSGetNumFields(ds,&numFields);
1628: DMGetDefaultConstraints(refTree,&refConSec,NULL);
1629: PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1630: for (p = pRefStart; p < pRefEnd; p++) {
1631: PetscInt parent, pDof;
1633: DMPlexGetTreeParent(refTree,p,&parent,NULL);
1634: PetscSectionGetDof(refConSec,p,&pDof);
1635: if (!pDof || parent == p) continue;
1637: for (f = 0; f < numFields; f++) {
1638: PetscInt cDof;
1640: if (numFields > 1) {
1641: PetscSectionGetFieldDof(refConSec,p,f,&cDof);
1642: }
1643: else {
1644: PetscSectionGetDof(refConSec,p,&cDof);
1645: }
1647: if (!cDof) continue;
1648: PetscFree(refPointFieldMats[p - pRefStart][f]);
1649: }
1650: PetscFree(refPointFieldMats[p - pRefStart]);
1651: PetscFree(refPointFieldN[p - pRefStart]);
1652: }
1653: PetscFree(refPointFieldMats);
1654: PetscFree(refPointFieldN);
1655: return(0);
1656: }
1660: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm, PetscSection section, PetscSection conSec, Mat cMat)
1661: {
1662: DM refTree;
1663: PetscDS ds;
1664: Mat refCmat;
1665: PetscInt numFields, f, pRefStart, pRefEnd, p, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN;
1666: PetscScalar ***refPointFieldMats, *pointWork;
1667: PetscSection refConSec, refAnSec, anSec;
1668: IS refAnIS, anIS;
1669: const PetscInt *anchors;
1674: DMGetDS(dm,&ds);
1675: PetscDSGetNumFields(ds,&numFields);
1676: DMPlexGetReferenceTree(dm,&refTree);
1677: DMSetDS(refTree,ds);
1678: DMGetDefaultConstraints(refTree,&refConSec,&refCmat);
1679: DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);
1680: DMPlexGetAnchors(dm,&anSec,&anIS);
1681: ISGetIndices(anIS,&anchors);
1682: PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1683: PetscSectionGetChart(conSec,&conStart,&conEnd);
1684: PetscSectionGetMaxDof(refConSec,&maxDof);
1685: PetscSectionGetMaxDof(refAnSec,&maxAnDof);
1686: PetscMalloc1(maxDof*maxDof*maxAnDof,&pointWork);
1688: /* step 1: get submats for every constrained point in the reference tree */
1689: DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
1691: /* step 2: compute the preorder */
1692: DMPlexGetChart(dm,&pStart,&pEnd);
1693: PetscMalloc2(pEnd-pStart,&perm,pEnd-pStart,&iperm);
1694: for (p = pStart; p < pEnd; p++) {
1695: perm[p - pStart] = p;
1696: iperm[p - pStart] = p-pStart;
1697: }
1698: for (p = 0; p < pEnd - pStart;) {
1699: PetscInt point = perm[p];
1700: PetscInt parent;
1702: DMPlexGetTreeParent(dm,point,&parent,NULL);
1703: if (parent == point) {
1704: p++;
1705: }
1706: else {
1707: PetscInt size, closureSize, *closure = NULL, i;
1709: DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1710: for (i = 0; i < closureSize; i++) {
1711: PetscInt q = closure[2*i];
1712: if (iperm[q-pStart] > iperm[point-pStart]) {
1713: /* swap */
1714: perm[p] = q;
1715: perm[iperm[q-pStart]] = point;
1716: iperm[point-pStart] = iperm[q-pStart];
1717: iperm[q-pStart] = p;
1718: break;
1719: }
1720: }
1721: size = closureSize;
1722: DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1723: if (i == size) {
1724: p++;
1725: }
1726: }
1727: }
1729: /* step 3: fill the constraint matrix */
1730: /* we are going to use a preorder progressive fill strategy. Mat doesn't
1731: * allow progressive fill without assembly, so we are going to set up the
1732: * values outside of the Mat first.
1733: */
1734: {
1735: PetscInt nRows, row, nnz;
1736: PetscBool done;
1737: const PetscInt *ia, *ja;
1738: PetscScalar *vals;
1740: MatGetRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);
1741: if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not get RowIJ of constraint matrix");
1742: nnz = ia[nRows];
1743: /* malloc and then zero rows right before we fill them: this way valgrind
1744: * can tell if we are doing progressive fill in the wrong order */
1745: PetscMalloc1(nnz,&vals);
1746: for (p = 0; p < pEnd - pStart; p++) {
1747: PetscInt parent, childid, closureSize, *closure = NULL;
1748: PetscInt point = perm[p], pointDof;
1750: DMPlexGetTreeParent(dm,point,&parent,&childid);
1751: if ((point < conStart) || (point >= conEnd) || (parent == point)) continue;
1752: PetscSectionGetDof(conSec,point,&pointDof);
1753: if (!pointDof) continue;
1754: DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1755: for (f = 0; f < numFields; f++) {
1756: PetscInt cDof, cOff, numCols, numFillCols, i, r, fComp = -1, matOffset, offset;
1757: PetscScalar *pointMat;
1758: PetscObject disc;
1759: PetscClassId id;
1760: PetscFE fe = NULL;
1761: PetscFV fv = NULL;
1763: PetscDSGetDiscretization(ds,f,&disc);
1764: PetscObjectGetClassId(disc,&id);
1765: if (id == PETSCFE_CLASSID) {
1766: fe = (PetscFE) disc;
1767: PetscFEGetNumComponents(fe,&fComp);
1768: }
1769: else if (id == PETSCFV_CLASSID) {
1770: fv = (PetscFV) disc;
1771: PetscFVGetNumComponents(fv,&fComp);
1772: }
1773: else {
1774: SETERRQ(PetscObjectComm((PetscObject)disc),PETSC_ERR_SUP,"Unsupported discretization");
1775: }
1777: if (numFields > 1) {
1778: PetscSectionGetFieldDof(conSec,point,f,&cDof);
1779: PetscSectionGetFieldOffset(conSec,point,f,&cOff);
1780: }
1781: else {
1782: PetscSectionGetDof(conSec,point,&cDof);
1783: PetscSectionGetOffset(conSec,point,&cOff);
1784: }
1785: if (!cDof) continue;
1787: /* make sure that every row for this point is the same size */
1788: #if defined(PETSC_USE_DEBUG)
1789: for (r = 0; r < cDof; r++) {
1790: if (cDof > 1 && r) {
1791: 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]));
1792: }
1793: }
1794: #endif
1795: /* zero rows */
1796: for (i = ia[cOff] ; i< ia[cOff+cDof];i++) {
1797: vals[i] = 0.;
1798: }
1799: matOffset = ia[cOff];
1800: numFillCols = ia[cOff+1] - matOffset;
1801: pointMat = refPointFieldMats[childid-pRefStart][f];
1802: numCols = refPointFieldN[childid-pRefStart][f];
1803: offset = 0;
1804: for (i = 0; i < closureSize; i++) {
1805: PetscInt q = closure[2*i];
1806: PetscInt o = closure[2*i+1];
1807: PetscInt aDof, aOff, j, k, qConDof, qConOff;
1809: qConDof = qConOff = 0;
1810: if (numFields > 1) {
1811: PetscSectionGetFieldDof(section,q,f,&aDof);
1812: PetscSectionGetFieldOffset(section,q,f,&aOff);
1813: if (q >= conStart && q < conEnd) {
1814: PetscSectionGetFieldDof(conSec,q,f,&qConDof);
1815: PetscSectionGetFieldOffset(conSec,q,f,&qConOff);
1816: }
1817: }
1818: else {
1819: PetscSectionGetDof(section,q,&aDof);
1820: PetscSectionGetOffset(section,q,&aOff);
1821: if (q >= conStart && q < conEnd) {
1822: PetscSectionGetDof(conSec,q,&qConDof);
1823: PetscSectionGetOffset(conSec,q,&qConOff);
1824: }
1825: }
1826: if (!aDof) continue;
1827: if (qConDof) {
1828: /* this point has anchors: its rows of the matrix should already
1829: * be filled, thanks to preordering */
1830: /* first multiply into pointWork, then set in matrix */
1831: PetscInt aMatOffset = ia[qConOff];
1832: PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset;
1833: for (r = 0; r < cDof; r++) {
1834: for (j = 0; j < aNumFillCols; j++) {
1835: PetscScalar inVal = 0;
1836: for (k = 0; k < aDof; k++) {
1837: PetscInt node = (o >= 0) ? (k / fComp) : ((aDof - 1 - k) / fComp);
1838: PetscInt comp = (k % fComp);
1839: PetscInt col = node * fComp + comp;
1841: inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j];
1842: }
1843: pointWork[r * aNumFillCols + j] = inVal;
1844: }
1845: }
1846: /* assume that the columns are sorted, spend less time searching */
1847: for (j = 0, k = 0; j < aNumFillCols; j++) {
1848: PetscInt col = ja[aMatOffset + j];
1849: for (;k < numFillCols; k++) {
1850: if (ja[matOffset + k] == col) {
1851: break;
1852: }
1853: }
1854: if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, col);
1855: for (r = 0; r < cDof; r++) {
1856: vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j];
1857: }
1858: }
1859: }
1860: else {
1861: /* find where to put this portion of pointMat into the matrix */
1862: for (k = 0; k < numFillCols; k++) {
1863: if (ja[matOffset + k] == aOff) {
1864: break;
1865: }
1866: }
1867: if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, aOff);
1868: for (j = 0; j < aDof; j++) {
1869: PetscInt node = (o >= 0) ? (j / fComp) : ((aDof - 1 - j) / fComp);
1870: PetscInt comp = (j % fComp);
1871: PetscInt col = node * fComp + comp;
1872: for (r = 0; r < cDof; r++) {
1873: vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + col];
1874: }
1875: }
1876: }
1877: offset += aDof;
1878: }
1879: }
1880: DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1881: }
1882: for (row = 0; row < nRows; row++) {
1883: MatSetValues(cMat,1,&row,ia[row+1]-ia[row],&ja[ia[row]],&vals[ia[row]],INSERT_VALUES);
1884: }
1885: MatRestoreRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);
1886: if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not restore RowIJ of constraint matrix");
1887: MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
1888: MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
1889: PetscFree(vals);
1890: }
1892: /* clean up */
1893: ISRestoreIndices(anIS,&anchors);
1894: PetscFree2(perm,iperm);
1895: PetscFree(pointWork);
1896: DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
1897: return(0);
1898: }
1902: /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of
1903: * a non-conforming mesh. Local refinement comes later */
1904: PetscErrorCode DMPlexTreeRefineCell (DM dm, PetscInt cell, DM *ncdm)
1905: {
1906: DM K;
1907: PetscMPIInt rank;
1908: PetscInt dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd;
1909: PetscInt numNewCones, *newConeSizes, *newCones, *newOrientations;
1910: PetscInt *Kembedding;
1911: PetscInt *cellClosure=NULL, nc;
1912: PetscScalar *newVertexCoords;
1913: PetscInt numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset;
1914: PetscSection parentSection;
1918: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1919: DMGetDimension(dm,&dim);
1920: DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm);
1921: DMSetDimension(*ncdm,dim);
1923: DMPlexGetChart(dm, &pStart, &pEnd);
1924: PetscSectionCreate(PetscObjectComm((PetscObject)dm),&parentSection);
1925: DMPlexGetReferenceTree(dm,&K);
1926: if (!rank) {
1927: /* compute the new charts */
1928: PetscMalloc5(dim+1,&pNewCount,dim+1,&pNewStart,dim+1,&pNewEnd,dim+1,&pOldStart,dim+1,&pOldEnd);
1929: offset = 0;
1930: for (d = 0; d <= dim; d++) {
1931: PetscInt pOldCount, kStart, kEnd, k;
1933: pNewStart[d] = offset;
1934: DMPlexGetHeightStratum(dm,d,&pOldStart[d],&pOldEnd[d]);
1935: DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
1936: pOldCount = pOldEnd[d] - pOldStart[d];
1937: /* adding the new points */
1938: pNewCount[d] = pOldCount + kEnd - kStart;
1939: if (!d) {
1940: /* removing the cell */
1941: pNewCount[d]--;
1942: }
1943: for (k = kStart; k < kEnd; k++) {
1944: PetscInt parent;
1945: DMPlexGetTreeParent(K,k,&parent,NULL);
1946: if (parent == k) {
1947: /* avoid double counting points that won't actually be new */
1948: pNewCount[d]--;
1949: }
1950: }
1951: pNewEnd[d] = pNewStart[d] + pNewCount[d];
1952: offset = pNewEnd[d];
1954: }
1955: 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]);
1956: /* get the current closure of the cell that we are removing */
1957: DMPlexGetTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);
1959: PetscMalloc1(pNewEnd[dim],&newConeSizes);
1960: {
1961: PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j;
1963: DMPlexGetChart(K,&kStart,&kEnd);
1964: PetscMalloc4(kEnd-kStart,&Kembedding,kEnd-kStart,&perm,kEnd-kStart,&iperm,kEnd-kStart,&preOrient);
1966: for (k = kStart; k < kEnd; k++) {
1967: perm[k - kStart] = k;
1968: iperm [k - kStart] = k - kStart;
1969: preOrient[k - kStart] = 0;
1970: }
1972: DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);
1973: for (j = 1; j < closureSizeK; j++) {
1974: PetscInt parentOrientA = closureK[2*j+1];
1975: PetscInt parentOrientB = cellClosure[2*j+1];
1976: PetscInt p, q;
1978: p = closureK[2*j];
1979: q = cellClosure[2*j];
1980: for (d = 0; d <= dim; d++) {
1981: if (q >= pOldStart[d] && q < pOldEnd[d]) {
1982: Kembedding[p] = (q - pOldStart[d]) + pNewStart[d];
1983: }
1984: }
1985: if (parentOrientA != parentOrientB) {
1986: PetscInt numChildren, i;
1987: const PetscInt *children;
1989: DMPlexGetTreeChildren(K,p,&numChildren,&children);
1990: for (i = 0; i < numChildren; i++) {
1991: PetscInt kPerm, oPerm;
1993: k = children[i];
1994: DMPlexReferenceTreeGetChildSymmetry(K,p,parentOrientA,0,k,parentOrientB,&oPerm,&kPerm);
1995: /* perm = what refTree position I'm in */
1996: perm[kPerm-kStart] = k;
1997: /* iperm = who is at this position */
1998: iperm[k-kStart] = kPerm-kStart;
1999: preOrient[kPerm-kStart] = oPerm;
2000: }
2001: }
2002: }
2003: DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);
2004: }
2005: PetscSectionSetChart(parentSection,0,pNewEnd[dim]);
2006: offset = 0;
2007: numNewCones = 0;
2008: for (d = 0; d <= dim; d++) {
2009: PetscInt kStart, kEnd, k;
2010: PetscInt p;
2011: PetscInt size;
2013: for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
2014: /* skip cell 0 */
2015: if (p == cell) continue;
2016: /* old cones to new cones */
2017: DMPlexGetConeSize(dm,p,&size);
2018: newConeSizes[offset++] = size;
2019: numNewCones += size;
2020: }
2022: DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
2023: for (k = kStart; k < kEnd; k++) {
2024: PetscInt kParent;
2026: DMPlexGetTreeParent(K,k,&kParent,NULL);
2027: if (kParent != k) {
2028: Kembedding[k] = offset;
2029: DMPlexGetConeSize(K,k,&size);
2030: newConeSizes[offset++] = size;
2031: numNewCones += size;
2032: if (kParent != 0) {
2033: PetscSectionSetDof(parentSection,Kembedding[k],1);
2034: }
2035: }
2036: }
2037: }
2039: PetscSectionSetUp(parentSection);
2040: PetscSectionGetStorageSize(parentSection,&numPointsWithParents);
2041: PetscMalloc2(numNewCones,&newCones,numNewCones,&newOrientations);
2042: PetscMalloc2(numPointsWithParents,&parents,numPointsWithParents,&childIDs);
2044: /* fill new cones */
2045: offset = 0;
2046: for (d = 0; d <= dim; d++) {
2047: PetscInt kStart, kEnd, k, l;
2048: PetscInt p;
2049: PetscInt size;
2050: const PetscInt *cone, *orientation;
2052: for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
2053: /* skip cell 0 */
2054: if (p == cell) continue;
2055: /* old cones to new cones */
2056: DMPlexGetConeSize(dm,p,&size);
2057: DMPlexGetCone(dm,p,&cone);
2058: DMPlexGetConeOrientation(dm,p,&orientation);
2059: for (l = 0; l < size; l++) {
2060: newCones[offset] = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1];
2061: newOrientations[offset++] = orientation[l];
2062: }
2063: }
2065: DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
2066: for (k = kStart; k < kEnd; k++) {
2067: PetscInt kPerm = perm[k], kParent;
2068: PetscInt preO = preOrient[k];
2070: DMPlexGetTreeParent(K,k,&kParent,NULL);
2071: if (kParent != k) {
2072: /* embed new cones */
2073: DMPlexGetConeSize(K,k,&size);
2074: DMPlexGetCone(K,kPerm,&cone);
2075: DMPlexGetConeOrientation(K,kPerm,&orientation);
2076: for (l = 0; l < size; l++) {
2077: PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size -(preO + 1) - l) % size);
2078: PetscInt newO, lSize, oTrue;
2080: q = iperm[cone[m]];
2081: newCones[offset] = Kembedding[q];
2082: DMPlexGetConeSize(K,q,&lSize);
2083: oTrue = orientation[m];
2084: oTrue = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2);
2085: newO = DihedralCompose(lSize,oTrue,preOrient[q]);
2086: newOrientations[offset++] = newO;
2087: }
2088: if (kParent != 0) {
2089: PetscInt newPoint = Kembedding[kParent];
2090: PetscSectionGetOffset(parentSection,Kembedding[k],&pOffset);
2091: parents[pOffset] = newPoint;
2092: childIDs[pOffset] = k;
2093: }
2094: }
2095: }
2096: }
2098: PetscMalloc1(dim*(pNewEnd[dim]-pNewStart[dim]),&newVertexCoords);
2100: /* fill coordinates */
2101: offset = 0;
2102: {
2103: PetscInt kStart, kEnd, l;
2104: PetscSection vSection;
2105: PetscInt v;
2106: Vec coords;
2107: PetscScalar *coordvals;
2108: PetscInt dof, off;
2109: PetscReal v0[3], J[9], detJ;
2111: #if defined(PETSC_USE_DEBUG)
2112: {
2113: PetscInt k;
2114: DMPlexGetHeightStratum(K,0,&kStart,&kEnd);
2115: for (k = kStart; k < kEnd; k++) {
2116: DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ);
2117: if (detJ <= 0.) SETERRQ1 (PETSC_COMM_SELF,PETSC_ERR_PLIB,"reference tree cell %d has bad determinant",k);
2118: }
2119: }
2120: #endif
2121: DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ);
2122: DMGetCoordinateSection(dm,&vSection);
2123: DMGetCoordinatesLocal(dm,&coords);
2124: VecGetArray(coords,&coordvals);
2125: for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) {
2127: PetscSectionGetDof(vSection,v,&dof);
2128: PetscSectionGetOffset(vSection,v,&off);
2129: for (l = 0; l < dof; l++) {
2130: newVertexCoords[offset++] = coordvals[off + l];
2131: }
2132: }
2133: VecRestoreArray(coords,&coordvals);
2135: DMGetCoordinateSection(K,&vSection);
2136: DMGetCoordinatesLocal(K,&coords);
2137: VecGetArray(coords,&coordvals);
2138: DMPlexGetDepthStratum(K,0,&kStart,&kEnd);
2139: for (v = kStart; v < kEnd; v++) {
2140: PetscReal coord[3], newCoord[3];
2141: PetscInt vPerm = perm[v];
2142: PetscInt kParent;
2144: DMPlexGetTreeParent(K,v,&kParent,NULL);
2145: if (kParent != v) {
2146: /* this is a new vertex */
2147: PetscSectionGetOffset(vSection,vPerm,&off);
2148: for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off+l]);
2149: CoordinatesRefToReal(dim, dim, v0, J, coord, newCoord);
2150: for (l = 0; l < dim; ++l) newVertexCoords[offset+l] = newCoord[l];
2151: offset += dim;
2152: }
2153: }
2154: VecRestoreArray(coords,&coordvals);
2155: }
2157: /* need to reverse the order of pNewCount: vertices first, cells last */
2158: for (d = 0; d < (dim + 1) / 2; d++) {
2159: PetscInt tmp;
2161: tmp = pNewCount[d];
2162: pNewCount[d] = pNewCount[dim - d];
2163: pNewCount[dim - d] = tmp;
2164: }
2166: DMPlexCreateFromDAG(*ncdm,dim,pNewCount,newConeSizes,newCones,newOrientations,newVertexCoords);
2167: DMPlexSetReferenceTree(*ncdm,K);
2168: DMPlexSetTree(*ncdm,parentSection,parents,childIDs);
2170: /* clean up */
2171: DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);
2172: PetscFree5(pNewCount,pNewStart,pNewEnd,pOldStart,pOldEnd);
2173: PetscFree(newConeSizes);
2174: PetscFree2(newCones,newOrientations);
2175: PetscFree(newVertexCoords);
2176: PetscFree2(parents,childIDs);
2177: PetscFree4(Kembedding,perm,iperm,preOrient);
2178: }
2179: else {
2180: PetscInt p, counts[4];
2181: PetscInt *coneSizes, *cones, *orientations;
2182: Vec coordVec;
2183: PetscScalar *coords;
2185: for (d = 0; d <= dim; d++) {
2186: PetscInt dStart, dEnd;
2188: DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);
2189: counts[d] = dEnd - dStart;
2190: }
2191: PetscMalloc1(pEnd-pStart,&coneSizes);
2192: for (p = pStart; p < pEnd; p++) {
2193: DMPlexGetConeSize(dm,p,&coneSizes[p-pStart]);
2194: }
2195: DMPlexGetCones(dm, &cones);
2196: DMPlexGetConeOrientations(dm, &orientations);
2197: DMGetCoordinatesLocal(dm,&coordVec);
2198: VecGetArray(coordVec,&coords);
2200: PetscSectionSetChart(parentSection,pStart,pEnd);
2201: PetscSectionSetUp(parentSection);
2202: DMPlexCreateFromDAG(*ncdm,dim,counts,coneSizes,cones,orientations,NULL);
2203: DMPlexSetReferenceTree(*ncdm,K);
2204: DMPlexSetTree(*ncdm,parentSection,NULL,NULL);
2205: VecRestoreArray(coordVec,&coords);
2206: }
2207: PetscSectionDestroy(&parentSection);
2209: return(0);
2210: }
2212: PETSC_EXTERN PetscErrorCode indicesPointFields_private(PetscSection,PetscInt,PetscInt,PetscInt[],PetscBool,PetscInt,PetscInt[]);
2213: PETSC_EXTERN PetscErrorCode indicesPoint_private(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,PetscInt,PetscInt[]);
2217: PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
2218: {
2219: PetscSF coarseToFineEmbedded;
2220: PetscSection globalCoarse, globalFine;
2221: PetscSection localCoarse, localFine;
2222: PetscSection aSec, cSec;
2223: PetscSection rootIndicesSec, rootMatricesSec;
2224: PetscSection leafIndicesSec, leafMatricesSec;
2225: PetscInt *rootIndices, *leafIndices;
2226: PetscScalar *rootMatrices, *leafMatrices;
2227: IS aIS;
2228: const PetscInt *anchors;
2229: Mat cMat;
2230: PetscInt numFields;
2231: PetscInt pStartC, pEndC, pStartF, pEndF, p;
2232: PetscInt aStart, aEnd, cStart, cEnd;
2233: PetscInt *maxChildIds;
2234: PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO;
2238: DMPlexGetChart(coarse,&pStartC,&pEndC);
2239: DMPlexGetChart(fine,&pStartF,&pEndF);
2240: DMGetDefaultGlobalSection(fine,&globalFine);
2241: { /* winnow fine points that don't have global dofs out of the sf */
2242: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;
2244: for (p = pStartF, numPointsWithDofs = 0; p < pEndF; p++) {
2245: PetscSectionGetDof(globalFine,p,&dof);
2246: PetscSectionGetConstraintDof(globalFine,p,&cdof);
2247: if ((dof - cdof) > 0) {
2248: numPointsWithDofs++;
2249: }
2250: }
2251: PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
2252: for (p = pStartF, offset = 0; p < pEndF; p++) {
2253: PetscSectionGetDof(globalFine,p,&dof);
2254: PetscSectionGetConstraintDof(globalFine,p,&cdof);
2255: if ((dof - cdof) > 0) {
2256: pointsWithDofs[offset++] = p - pStartF;
2257: }
2258: }
2259: PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
2260: PetscFree(pointsWithDofs);
2261: }
2262: /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */
2263: PetscMalloc1(pEndC-pStartC,&maxChildIds);
2264: for (p = pStartC; p < pEndC; p++) {
2265: maxChildIds[p - pStartC] = -2;
2266: }
2267: PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);
2268: PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);
2270: DMGetDefaultSection(coarse,&localCoarse);
2271: DMGetDefaultGlobalSection(coarse,&globalCoarse);
2273: DMPlexGetAnchors(coarse,&aSec,&aIS);
2274: ISGetIndices(aIS,&anchors);
2275: PetscSectionGetChart(aSec,&aStart,&aEnd);
2277: DMGetDefaultConstraints(coarse,&cSec,&cMat);
2278: PetscSectionGetChart(cSec,&cStart,&cEnd);
2280: /* create sections that will send to children the indices and matrices they will need to construct the interpolator */
2281: PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);
2282: PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootMatricesSec);
2283: PetscSectionSetChart(rootIndicesSec,pStartC,pEndC);
2284: PetscSectionSetChart(rootMatricesSec,pStartC,pEndC);
2285: PetscSectionGetNumFields(globalCoarse,&numFields);
2286: {
2287: PetscInt maxFields = PetscMax(1,numFields) + 1;
2288: PetscMalloc7(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&newOffsets,maxFields,&newOffsetsCopy,maxFields,&rowOffsets,maxFields,&numD,maxFields,&numO);
2289: }
2291: for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */
2292: PetscInt dof, matSize = 0;
2293: PetscInt aDof = 0;
2294: PetscInt cDof = 0;
2295: PetscInt maxChildId = maxChildIds[p - pStartC];
2296: PetscInt numRowIndices = 0;
2297: PetscInt numColIndices = 0;
2299: PetscSectionGetDof(globalCoarse,p,&dof);
2300: if (dof < 0) {
2301: dof = -(dof + 1);
2302: }
2303: if (p >= aStart && p < aEnd) {
2304: PetscSectionGetDof(aSec,p,&aDof);
2305: }
2306: if (p >= cStart && p < cEnd) {
2307: PetscSectionGetDof(cSec,p,&cDof);
2308: }
2309: offsets[0] = 0;
2310: newOffsets[0] = 0;
2311: if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */
2312: PetscInt *closure = NULL, closureSize, cl, f;
2314: DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);
2315: for (cl = 0; cl < closureSize; cl++) { /* get the closure */
2316: PetscInt c = closure[2 * cl], clDof;
2318: PetscSectionGetDof(localCoarse,c,&clDof);
2319: numRowIndices += clDof;
2320: for (f = 0; f < numFields; f++) {
2321: PetscSectionGetFieldDof(localCoarse,c,f,&clDof);
2322: offsets[f + 1] += clDof;
2323: }
2324: }
2325: for (f = 0; f < numFields; f++) {
2326: offsets[f + 1] += offsets[f];
2327: newOffsets[f + 1] = offsets[f + 1];
2328: }
2329: /* get the number of indices needed and their field offsets */
2330: DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,NULL,NULL,&numColIndices,NULL,NULL,newOffsets,PETSC_FALSE);
2331: DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);
2332: if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */
2333: numColIndices = numRowIndices;
2334: matSize = 0;
2335: }
2336: else if (numFields) { /* we send one submat for each field: sum their sizes */
2337: matSize = 0;
2338: for (f = 0; f < numFields; f++) {
2339: PetscInt numRow, numCol;
2341: numRow = offsets[f + 1] - offsets[f];
2342: numCol = newOffsets[f + 1] - newOffsets[f + 1];
2343: matSize += numRow * numCol;
2344: }
2345: }
2346: else {
2347: matSize = numRowIndices * numColIndices;
2348: }
2349: }
2350: else if (maxChildId == -1) {
2351: if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */
2352: PetscInt aOff, a, f;
2354: PetscSectionGetOffset(aSec,p,&aOff);
2355: for (f = 0; f < numFields; f++) {
2356: PetscInt fDof;
2358: PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
2359: offsets[f+1] = fDof;
2360: }
2361: for (a = 0; a < aDof; a++) {
2362: PetscInt anchor = anchors[a + aOff], aLocalDof;
2364: PetscSectionGetDof(localCoarse,anchor,&aLocalDof);
2365: numColIndices += aLocalDof;
2366: for (f = 0; f < numFields; f++) {
2367: PetscInt fDof;
2369: PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);
2370: newOffsets[f+1] += fDof;
2371: }
2372: }
2373: if (numFields) {
2374: matSize = 0;
2375: for (f = 0; f < numFields; f++) {
2376: matSize += offsets[f+1] * newOffsets[f+1];
2377: }
2378: }
2379: else {
2380: matSize = numColIndices * dof;
2381: }
2382: }
2383: else { /* no children, and no constraints on dofs: just get the global indices */
2384: numColIndices = dof;
2385: matSize = 0;
2386: }
2387: }
2388: /* we will pack the column indices with the field offsets */
2389: PetscSectionSetDof(rootIndicesSec,p,numColIndices ? numColIndices+2*numFields : 0);
2390: PetscSectionSetDof(rootMatricesSec,p,matSize);
2391: }
2392: PetscSectionSetUp(rootIndicesSec);
2393: PetscSectionSetUp(rootMatricesSec);
2394: {
2395: PetscInt numRootIndices, numRootMatrices;
2397: PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);
2398: PetscSectionGetStorageSize(rootMatricesSec,&numRootMatrices);
2399: PetscMalloc2(numRootIndices,&rootIndices,numRootMatrices,&rootMatrices);
2400: for (p = pStartC; p < pEndC; p++) {
2401: PetscInt numRowIndices, numColIndices, matSize, dof;
2402: PetscInt pIndOff, pMatOff;
2403: PetscInt *pInd;
2404: PetscInt maxChildId = maxChildIds[p - pStartC];
2405: PetscScalar *pMat = NULL;
2407: PetscSectionGetDof(rootIndicesSec,p,&numColIndices);
2408: if (!numColIndices) {
2409: continue;
2410: }
2411: offsets[0] = 0;
2412: newOffsets[0] = 0;
2413: offsetsCopy[0] = 0;
2414: newOffsetsCopy[0] = 0;
2415: numColIndices -= 2 * numFields;
2416: PetscSectionGetOffset(rootIndicesSec,p,&pIndOff);
2417: pInd = &(rootIndices[pIndOff]);
2418: PetscSectionGetDof(rootMatricesSec,p,&matSize);
2419: if (matSize) {
2420: PetscSectionGetOffset(rootMatricesSec,p,&pMatOff);
2421: pMat = &rootMatrices[pMatOff];
2422: }
2423: PetscSectionGetDof(globalCoarse,p,&dof);
2424: if (dof < 0) {
2425: dof = -(dof + 1);
2426: }
2427: if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */
2428: PetscInt i, j;
2429: PetscInt numRowIndices = matSize / numColIndices;
2431: if (!numRowIndices) { /* don't need to calculate the mat, just the indices */
2432: PetscInt numIndices, *indices;
2433: DMPlexGetClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);
2434: if (numIndices != numColIndices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"mismatching constraint indices calculations");
2435: for (i = 0; i < numColIndices; i++) {
2436: pInd[i] = indices[i];
2437: }
2438: for (i = 0; i < numFields; i++) {
2439: pInd[numColIndices + i] = offsets[i+1];
2440: pInd[numColIndices + numFields + i] = offsets[i+1];
2441: }
2442: DMPlexRestoreClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);
2443: }
2444: else {
2445: PetscInt closureSize, *closure = NULL, cl;
2446: PetscScalar *pMatIn, *pMatModified;
2447: PetscInt numPoints,*points;
2449: DMGetWorkArray(coarse,numRowIndices * numRowIndices,PETSC_SCALAR,&pMatIn);
2450: for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */
2451: for (j = 0; j < numRowIndices; j++) {
2452: pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.;
2453: }
2454: }
2455: DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);
2456: if (numFields) {
2457: PetscInt f;
2459: for (cl = 0; cl < closureSize; cl++) {
2460: PetscInt c = closure[2 * cl];
2462: for (f = 0; f < numFields; f++) {
2463: PetscInt fDof;
2465: PetscSectionGetFieldDof(localCoarse,c,f,&fDof);
2466: offsets[f + 1] += fDof;
2467: }
2468: }
2469: for (f = 0; f < numFields; f++) {
2470: offsets[f + 1] += offsets[f];
2471: newOffsets[f + 1] = offsets[f + 1];
2472: }
2473: }
2474: /* apply hanging node constraints on the right, get the new points and the new offsets */
2475: DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,pMatIn,&numPoints,NULL,&points,&pMatModified,newOffsets,PETSC_FALSE);
2476: if (!numFields) {
2477: for (i = 0; i < numRowIndices * numColIndices; i++) {
2478: pMat[i] = pMatModified[i];
2479: }
2480: }
2481: else {
2482: PetscInt f, i, j, count;
2483: for (f = 0, count = 0; f < numFields; f++) {
2484: for (i = offsets[f]; i < offsets[f+1]; i++) {
2485: for (j = newOffsets[f]; j < newOffsets[f+1]; j++, count++) {
2486: pMat[count] = pMatModified[i * numColIndices + j];
2487: }
2488: }
2489: }
2490: }
2491: DMRestoreWorkArray(coarse,numRowIndices * numColIndices,PETSC_SCALAR,&pMatModified);
2492: DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);
2493: DMRestoreWorkArray(coarse,numRowIndices * numColIndices,PETSC_SCALAR,&pMatIn);
2494: if (numFields) {
2495: PetscInt f;
2496: for (f = 0; f < numFields; f++) {
2497: pInd[numColIndices + f] = offsets[f+1];
2498: pInd[numColIndices + numFields + f] = newOffsets[f+1];
2499: }
2500: for (cl = 0; cl < numPoints*2; cl += 2) {
2501: PetscInt o = points[cl+1], globalOff, c = points[cl];
2502: PetscSectionGetOffset(globalCoarse, c, &globalOff);
2503: indicesPointFields_private(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, o, pInd);
2504: }
2505: } else {
2506: for (cl = 0; cl < numPoints*2; cl += 2) {
2507: PetscInt o = points[cl+1], c = points[cl], globalOff;
2508: PetscSectionGetOffset(globalCoarse, c, &globalOff);
2509: indicesPoint_private(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, o, pInd);
2510: }
2511: }
2512: DMRestoreWorkArray(coarse,numPoints,PETSC_SCALAR,&points);
2513: }
2514: }
2515: else if (matSize) {
2516: PetscInt cOff;
2517: PetscInt *rowIndices, *colIndices, a, aDof, aOff;
2519: numRowIndices = matSize / numColIndices;
2520: if (numRowIndices != dof) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Miscounted dofs");
2521: DMGetWorkArray(coarse,numRowIndices,PETSC_INT,&rowIndices);
2522: DMGetWorkArray(coarse,numColIndices,PETSC_INT,&colIndices);
2523: PetscSectionGetOffset(cSec,p,&cOff);
2524: PetscSectionGetDof(aSec,p,&aDof);
2525: PetscSectionGetOffset(aSec,p,&aOff);
2526: if (numFields) {
2527: PetscInt f;
2529: for (f = 0; f < numFields; f++) {
2530: PetscInt fDof;
2531: PetscSectionGetFieldDof(cSec,p,f,&fDof);
2532: offsets[f + 1] = fDof;
2533: for (a = 0; a < aDof; a++) {
2534: PetscInt anchor = anchors[a + aOff];
2535: PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);
2536: newOffsets[f + 1] += fDof;
2537: }
2538: }
2539: for (f = 0; f < numFields; f++) {
2540: offsets[f + 1] += offsets[f];
2541: offsetsCopy[f + 1] = offsets[f + 1];
2542: newOffsets[f + 1] += newOffsets[f];
2543: newOffsetsCopy[f + 1] = newOffsets[f + 1];
2544: }
2545: indicesPointFields_private(cSec,p,cOff,offsetsCopy,PETSC_TRUE,0,rowIndices);
2546: for (a = 0; a < aDof; a++) {
2547: PetscInt anchor = anchors[a + aOff], lOff;
2548: PetscSectionGetOffset(localCoarse,anchor,&lOff);
2549: indicesPointFields_private(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,0,colIndices);
2550: }
2551: }
2552: else {
2553: indicesPoint_private(cSec,p,cOff,offsetsCopy,PETSC_TRUE,0,rowIndices);
2554: for (a = 0; a < aDof; a++) {
2555: PetscInt anchor = anchors[a + aOff], lOff;
2556: PetscSectionGetOffset(localCoarse,anchor,&lOff);
2557: indicesPoint_private(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,0,colIndices);
2558: }
2559: }
2560: if (numFields) {
2561: PetscInt count, f, a;
2562: for (f = 0, count = 0; f < numFields; f++) {
2563: PetscInt iSize = offsets[f + 1] - offsets[f];
2564: PetscInt jSize = newOffsets[f + 1] - newOffsets[f];
2565: MatGetValues(cMat,iSize,&rowIndices[offsets[f]],jSize,&colIndices[newOffsets[f]],&pMat[count]);
2566: count += iSize * jSize;
2567: pInd[numColIndices + f] = offsets[f+1];
2568: pInd[numColIndices + numFields + f] = newOffsets[f+1];
2569: }
2570: for (a = 0; a < aDof; a++) {
2571: PetscInt anchor = anchors[a + aOff];
2572: PetscInt gOff;
2573: PetscSectionGetOffset(globalCoarse,anchor,&gOff);
2574: indicesPointFields_private(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,0,pInd);
2575: }
2576: }
2577: else {
2578: PetscInt a;
2579: MatGetValues(cMat,numRowIndices,rowIndices,numColIndices,colIndices,pMat);
2580: for (a = 0; a < aDof; a++) {
2581: PetscInt anchor = anchors[a + aOff];
2582: PetscInt gOff;
2583: PetscSectionGetOffset(globalCoarse,anchor,&gOff);
2584: indicesPoint_private(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,0,pInd);
2585: }
2586: }
2587: DMRestoreWorkArray(coarse,numColIndices,PETSC_INT,&colIndices);
2588: DMRestoreWorkArray(coarse,numRowIndices,PETSC_INT,&rowIndices);
2589: }
2590: else {
2591: PetscInt gOff;
2593: PetscSectionGetOffset(globalCoarse,p,&gOff);
2594: if (numFields) {
2595: PetscInt f;
2597: for (f = 0; f < numFields; f++) {
2598: PetscInt fDof;
2599: PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
2600: offsets[f + 1] = fDof + offsets[f];
2601: }
2602: for (f = 0; f < numFields; f++) {
2603: pInd[numColIndices + f] = offsets[f+1];
2604: pInd[numColIndices + numFields + f] = offsets[f+1];
2605: }
2606: indicesPointFields_private(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,0,pInd);
2607: }
2608: else {
2609: indicesPoint_private(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,0,pInd);
2610: }
2611: }
2612: }
2613: PetscFree(maxChildIds);
2614: }
2615: {
2616: PetscSF indicesSF, matricesSF;
2617: PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices;
2619: PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);
2620: PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafMatricesSec);
2621: PetscSFDistributeSection(coarseToFineEmbedded,rootIndicesSec,&remoteOffsetsIndices,leafIndicesSec);
2622: PetscSFDistributeSection(coarseToFineEmbedded,rootMatricesSec,&remoteOffsetsMatrices,leafMatricesSec);
2623: PetscSFCreateSectionSF(coarseToFineEmbedded,rootIndicesSec,remoteOffsetsIndices,leafIndicesSec,&indicesSF);
2624: PetscSFCreateSectionSF(coarseToFineEmbedded,rootMatricesSec,remoteOffsetsMatrices,leafMatricesSec,&matricesSF);
2625: PetscSFDestroy(&coarseToFineEmbedded);
2626: PetscFree(remoteOffsetsIndices);
2627: PetscFree(remoteOffsetsMatrices);
2628: PetscSectionGetStorageSize(leafIndicesSec,&numLeafIndices);
2629: PetscSectionGetStorageSize(leafMatricesSec,&numLeafMatrices);
2630: PetscMalloc2(numLeafIndices,&leafIndices,numLeafMatrices,&leafMatrices);
2631: PetscSFBcastBegin(indicesSF,MPIU_INT,rootIndices,leafIndices);
2632: PetscSFBcastBegin(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);
2633: PetscSFBcastEnd(indicesSF,MPIU_INT,rootIndices,leafIndices);
2634: PetscSFBcastEnd(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);
2635: PetscSFDestroy(&matricesSF);
2636: PetscSFDestroy(&indicesSF);
2637: PetscFree2(rootIndices,rootMatrices);
2638: PetscSectionDestroy(&rootIndicesSec);
2639: PetscSectionDestroy(&rootMatricesSec);
2640: }
2641: /* count to preallocate */
2642: DMGetDefaultSection(fine,&localFine);
2643: {
2644: PetscInt nGlobal;
2645: PetscInt *dnnz, *onnz;
2646: PetscLayout rowMap, colMap;
2647: PetscInt rowStart, rowEnd, colStart, colEnd;
2648: PetscInt maxDof;
2649: PetscInt *rowIndices;
2650: DM refTree;
2651: PetscInt **refPointFieldN;
2652: PetscScalar ***refPointFieldMats;
2653: PetscSection refConSec, refAnSec;
2654: PetscInt pRefStart,pRefEnd,maxConDof,maxColumns;
2655: PetscScalar *pointWork;
2657: PetscSectionGetConstrainedStorageSize(globalFine,&nGlobal);
2658: PetscCalloc2(nGlobal,&dnnz,nGlobal,&onnz);
2659: MatGetLayouts(mat,&rowMap,&colMap);
2660: PetscLayoutSetUp(rowMap);
2661: PetscLayoutSetUp(colMap);
2662: PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);
2663: PetscLayoutGetRange(colMap,&colStart,&colEnd);
2664: PetscSectionGetMaxDof(globalFine,&maxDof);
2665: DMGetWorkArray(fine,maxDof,PETSC_INT,&rowIndices);
2666: for (p = pStartF; p < pEndF; p++) {
2667: PetscInt gDof, gcDof, gOff;
2668: PetscInt numColIndices, pIndOff, *pInd;
2669: PetscInt matSize;
2670: PetscInt i;
2672: PetscSectionGetDof(globalFine,p,&gDof);
2673: PetscSectionGetConstraintDof(globalFine,p,&gcDof);
2674: if ((gDof - gcDof) <= 0) {
2675: continue;
2676: }
2677: PetscSectionGetOffset(globalFine,p,&gOff);
2678: if (gOff < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I though having global dofs meant a non-negative offset");
2679: if ((gOff < rowStart) || ((gOff + gDof - gcDof) > rowEnd)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I thought the row map would constrain the global dofs");
2680: PetscSectionGetDof(leafIndicesSec,p,&numColIndices);
2681: PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);
2682: numColIndices -= 2 * numFields;
2683: if (numColIndices <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"global fine dof with no dofs to interpolate from");
2684: pInd = &leafIndices[pIndOff];
2685: offsets[0] = 0;
2686: offsetsCopy[0] = 0;
2687: newOffsets[0] = 0;
2688: newOffsetsCopy[0] = 0;
2689: if (numFields) {
2690: PetscInt f;
2691: for (f = 0; f < numFields; f++) {
2692: PetscInt rowDof;
2694: PetscSectionGetFieldDof(localFine,p,f,&rowDof);
2695: offsets[f + 1] = offsets[f] + rowDof;
2696: offsetsCopy[f + 1] = offsets[f + 1];
2697: newOffsets[f + 1] = pInd[numColIndices + numFields + f];
2698: numD[f] = 0;
2699: numO[f] = 0;
2700: }
2701: indicesPointFields_private(localFine,p,gOff,offsetsCopy,PETSC_FALSE,0,rowIndices);
2702: for (f = 0; f < numFields; f++) {
2703: PetscInt colOffset = newOffsets[f];
2704: PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f];
2706: for (i = 0; i < numFieldCols; i++) {
2707: PetscInt gInd = pInd[i + colOffset];
2709: if (gInd >= colStart && gInd < colEnd) {
2710: numD[f]++;
2711: }
2712: else if (gInd >= 0) { /* negative means non-entry */
2713: numO[f]++;
2714: }
2715: }
2716: }
2717: }
2718: else {
2719: indicesPoint_private(localFine,p,gOff,offsetsCopy,PETSC_FALSE,0,rowIndices);
2720: numD[0] = 0;
2721: numO[0] = 0;
2722: for (i = 0; i < numColIndices; i++) {
2723: PetscInt gInd = pInd[i];
2725: if (gInd >= colStart && gInd < colEnd) {
2726: numD[0]++;
2727: }
2728: else if (gInd >= 0) { /* negative means non-entry */
2729: numO[0]++;
2730: }
2731: }
2732: }
2733: PetscSectionGetDof(leafMatricesSec,p,&matSize);
2734: if (!matSize) { /* incoming matrix is identity */
2735: PetscInt childId;
2737: childId = childIds[p-pStartF];
2738: if (childId < 0) { /* no child interpolation: one nnz per */
2739: if (numFields) {
2740: PetscInt f;
2741: for (f = 0; f < numFields; f++) {
2742: PetscInt numRows = offsets[f+1] - offsets[f], row;
2743: for (row = 0; row < numRows; row++) {
2744: PetscInt gIndCoarse = pInd[newOffsets[f] + row];
2745: PetscInt gIndFine = rowIndices[offsets[f] + row];
2746: if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2747: if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2748: dnnz[gIndFine - rowStart] = 1;
2749: }
2750: else if (gIndCoarse >= 0) { /* remote */
2751: if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2752: onnz[gIndFine - rowStart] = 1;
2753: }
2754: else { /* constrained */
2755: if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2756: }
2757: }
2758: }
2759: }
2760: else {
2761: PetscInt i;
2762: for (i = 0; i < gDof; i++) {
2763: PetscInt gIndCoarse = pInd[i];
2764: PetscInt gIndFine = rowIndices[i];
2765: if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */
2766: if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2767: dnnz[gIndFine - rowStart] = 1;
2768: }
2769: else if (gIndCoarse >= 0) { /* remote */
2770: if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2771: onnz[gIndFine - rowStart] = 1;
2772: }
2773: else { /* constrained */
2774: if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2775: }
2776: }
2777: }
2778: }
2779: else { /* interpolate from all */
2780: if (numFields) {
2781: PetscInt f;
2782: for (f = 0; f < numFields; f++) {
2783: PetscInt numRows = offsets[f+1] - offsets[f], row;
2784: for (row = 0; row < numRows; row++) {
2785: PetscInt gIndFine = rowIndices[offsets[f] + row];
2786: if (gIndFine >= 0) {
2787: if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2788: dnnz[gIndFine - rowStart] = numD[f];
2789: onnz[gIndFine - rowStart] = numO[f];
2790: }
2791: }
2792: }
2793: }
2794: else {
2795: PetscInt i;
2796: for (i = 0; i < gDof; i++) {
2797: PetscInt gIndFine = rowIndices[i];
2798: if (gIndFine >= 0) {
2799: if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2800: dnnz[gIndFine - rowStart] = numD[0];
2801: onnz[gIndFine - rowStart] = numO[0];
2802: }
2803: }
2804: }
2805: }
2806: }
2807: else { /* interpolate from all */
2808: if (numFields) {
2809: PetscInt f;
2810: for (f = 0; f < numFields; f++) {
2811: PetscInt numRows = offsets[f+1] - offsets[f], row;
2812: for (row = 0; row < numRows; row++) {
2813: PetscInt gIndFine = rowIndices[offsets[f] + row];
2814: if (gIndFine >= 0) {
2815: if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2816: dnnz[gIndFine - rowStart] = numD[f];
2817: onnz[gIndFine - rowStart] = numO[f];
2818: }
2819: }
2820: }
2821: }
2822: else { /* every dof get a full row */
2823: PetscInt i;
2824: for (i = 0; i < gDof; i++) {
2825: PetscInt gIndFine = rowIndices[i];
2826: if (gIndFine >= 0) {
2827: if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs");
2828: dnnz[gIndFine - rowStart] = numD[0];
2829: onnz[gIndFine - rowStart] = numO[0];
2830: }
2831: }
2832: }
2833: }
2834: }
2835: MatXAIJSetPreallocation(mat,1,dnnz,onnz,NULL,NULL);
2836: PetscFree2(dnnz,onnz);
2838: DMPlexGetReferenceTree(fine,&refTree);
2839: DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
2840: DMGetDefaultConstraints(refTree,&refConSec,NULL);
2841: DMPlexGetAnchors(refTree,&refAnSec,NULL);
2842: PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
2843: PetscSectionGetMaxDof(refConSec,&maxConDof);
2844: PetscSectionGetMaxDof(leafIndicesSec,&maxColumns);
2845: PetscMalloc1(maxConDof*maxColumns,&pointWork);
2846: for (p = pStartF; p < pEndF; p++) {
2847: PetscInt gDof, gcDof, gOff;
2848: PetscInt numColIndices, pIndOff, *pInd;
2849: PetscInt matSize;
2850: PetscInt childId;
2853: PetscSectionGetDof(globalFine,p,&gDof);
2854: PetscSectionGetConstraintDof(globalFine,p,&gcDof);
2855: if ((gDof - gcDof) <= 0) {
2856: continue;
2857: }
2858: childId = childIds[p-pStartF];
2859: PetscSectionGetOffset(globalFine,p,&gOff);
2860: PetscSectionGetDof(leafIndicesSec,p,&numColIndices);
2861: PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);
2862: numColIndices -= 2 * numFields;
2863: pInd = &leafIndices[pIndOff];
2864: offsets[0] = 0;
2865: offsetsCopy[0] = 0;
2866: newOffsets[0] = 0;
2867: newOffsetsCopy[0] = 0;
2868: rowOffsets[0] = 0;
2869: if (numFields) {
2870: PetscInt f;
2871: for (f = 0; f < numFields; f++) {
2872: PetscInt rowDof;
2874: PetscSectionGetFieldDof(localFine,p,f,&rowDof);
2875: offsets[f + 1] = offsets[f] + rowDof;
2876: offsetsCopy[f + 1] = offsets[f + 1];
2877: rowOffsets[f + 1] = pInd[numColIndices + f];
2878: newOffsets[f + 1] = pInd[numColIndices + numFields + f];
2879: }
2880: indicesPointFields_private(localFine,p,gOff,offsetsCopy,PETSC_FALSE,0,rowIndices);
2881: }
2882: else {
2883: indicesPoint_private(localFine,p,gOff,offsetsCopy,PETSC_FALSE,0,rowIndices);
2884: }
2885: PetscSectionGetDof(leafMatricesSec,p,&matSize);
2886: if (!matSize) { /* incoming matrix is identity */
2887: if (childId < 0) { /* no child interpolation: scatter */
2888: if (numFields) {
2889: PetscInt f;
2890: for (f = 0; f < numFields; f++) {
2891: PetscInt numRows = offsets[f+1] - offsets[f], row;
2892: for (row = 0; row < numRows; row++) {
2893: MatSetValue(mat,rowIndices[offsets[f]+row],pInd[newOffsets[f]+row],1.,INSERT_VALUES);
2894: }
2895: }
2896: }
2897: else {
2898: PetscInt numRows = gDof, row;
2899: for (row = 0; row < numRows; row++) {
2900: MatSetValue(mat,rowIndices[row],pInd[row],1.,INSERT_VALUES);
2901: }
2902: }
2903: }
2904: else { /* interpolate from all */
2905: if (numFields) {
2906: PetscInt f;
2907: for (f = 0; f < numFields; f++) {
2908: PetscInt numRows = offsets[f+1] - offsets[f];
2909: PetscInt numCols = newOffsets[f + 1] - newOffsets[f];
2910: MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],refPointFieldMats[childId - pRefStart][f],INSERT_VALUES);
2911: }
2912: }
2913: else {
2914: MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,refPointFieldMats[childId - pRefStart][0],INSERT_VALUES);
2915: }
2916: }
2917: }
2918: else { /* interpolate from all */
2919: PetscInt pMatOff;
2920: PetscScalar *pMat;
2922: PetscSectionGetOffset(leafMatricesSec,p,&pMatOff);
2923: pMat = &leafMatrices[pMatOff];
2924: if (childId < 0) { /* copy the incoming matrix */
2925: if (numFields) {
2926: PetscInt f, count;
2927: for (f = 0, count = 0; f < numFields; f++) {
2928: PetscInt numRows = offsets[f+1]-offsets[f];
2929: PetscInt numCols = newOffsets[f+1]-newOffsets[f];
2930: PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f];
2931: PetscScalar *inMat = &pMat[count];
2933: MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],inMat,INSERT_VALUES);
2934: count += numCols * numInRows;
2935: }
2936: }
2937: else {
2938: MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,pMat,INSERT_VALUES);
2939: }
2940: }
2941: else { /* multiply the incoming matrix by the child interpolation */
2942: if (numFields) {
2943: PetscInt f, count;
2944: for (f = 0, count = 0; f < numFields; f++) {
2945: PetscInt numRows = offsets[f+1]-offsets[f];
2946: PetscInt numCols = newOffsets[f+1]-newOffsets[f];
2947: PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f];
2948: PetscScalar *inMat = &pMat[count];
2949: PetscInt i, j, k;
2950: if (refPointFieldN[childId - pRefStart][f] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch");
2951: for (i = 0; i < numRows; i++) {
2952: for (j = 0; j < numCols; j++) {
2953: PetscScalar val = 0.;
2954: for (k = 0; k < numInRows; k++) {
2955: val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j];
2956: }
2957: pointWork[i * numCols + j] = val;
2958: }
2959: }
2960: MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],pointWork,INSERT_VALUES);
2961: count += numCols * numInRows;
2962: }
2963: }
2964: else { /* every dof gets a full row */
2965: PetscInt numRows = gDof;
2966: PetscInt numCols = numColIndices;
2967: PetscInt numInRows = matSize / numColIndices;
2968: PetscInt i, j, k;
2969: if (refPointFieldN[childId - pRefStart][0] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch");
2970: for (i = 0; i < numRows; i++) {
2971: for (j = 0; j < numCols; j++) {
2972: PetscScalar val = 0.;
2973: for (k = 0; k < numInRows; k++) {
2974: val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j];
2975: }
2976: pointWork[i * numCols + j] = val;
2977: }
2978: }
2979: MatSetValues(mat,numRows,rowIndices,numCols,pInd,pointWork,INSERT_VALUES);
2980: }
2981: }
2982: }
2983: }
2984: DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);
2985: DMRestoreWorkArray(fine,maxDof,PETSC_INT,&rowIndices);
2986: PetscFree(pointWork);
2987: }
2988: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
2989: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
2990: PetscSectionDestroy(&leafIndicesSec);
2991: PetscSectionDestroy(&leafMatricesSec);
2992: PetscFree2(leafIndices,leafMatrices);
2993: PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);
2994: ISRestoreIndices(aIS,&anchors);
2995: return(0);
2996: }
3000: /*
3001: * Assuming a nodal basis (w.r.t. the dual basis) basis:
3002: *
3003: * for each coarse dof \phi^c_i:
3004: * for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i:
3005: * for each fine dof \phi^f_j;
3006: * a_{i,j} = 0;
3007: * for each fine dof \phi^f_k:
3008: * a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l
3009: * [^^^ this is = \phi^c_i ^^^]
3010: */
3011: PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj)
3012: {
3013: PetscDS ds;
3014: PetscSection section, cSection;
3015: DMLabel canonical, depth;
3016: Mat cMat, mat;
3017: PetscInt *nnz;
3018: PetscInt f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd;
3019: PetscInt m, n;
3020: PetscScalar *pointScalar;
3021: PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent;
3025: DMGetDefaultSection(refTree,§ion);
3026: DMGetDimension(refTree, &dim);
3027: PetscMalloc6(dim,&v0,dim,&v0parent,dim,&vtmp,dim*dim,&J,dim*dim,&Jparent,dim*dim,&invJ);
3028: PetscMalloc2(dim,&pointScalar,dim,&pointRef);
3029: DMGetDS(refTree,&ds);
3030: PetscDSGetNumFields(ds,&numFields);
3031: PetscSectionGetNumFields(section,&numSecFields);
3032: DMGetLabel(refTree,"canonical",&canonical);
3033: DMGetLabel(refTree,"depth",&depth);
3034: DMGetDefaultConstraints(refTree,&cSection,&cMat);
3035: DMPlexGetChart(refTree, &pStart, &pEnd);
3036: DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd);
3037: MatGetSize(cMat,&n,&m); /* the injector has transpose sizes from the constraint matrix */
3038: /* Step 1: compute non-zero pattern. A proper subset of constraint matrix non-zero */
3039: PetscCalloc1(m,&nnz);
3040: for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */
3041: const PetscInt *children;
3042: PetscInt numChildren;
3043: PetscInt i, numChildDof, numSelfDof;
3045: if (canonical) {
3046: PetscInt pCanonical;
3047: DMLabelGetValue(canonical,p,&pCanonical);
3048: if (p != pCanonical) continue;
3049: }
3050: DMPlexGetTreeChildren(refTree,p,&numChildren,&children);
3051: if (!numChildren) continue;
3052: for (i = 0, numChildDof = 0; i < numChildren; i++) {
3053: PetscInt child = children[i];
3054: PetscInt dof;
3056: PetscSectionGetDof(section,child,&dof);
3057: numChildDof += dof;
3058: }
3059: PetscSectionGetDof(section,p,&numSelfDof);
3060: if (!numChildDof || !numSelfDof) continue;
3061: for (f = 0; f < numFields; f++) {
3062: PetscInt selfOff;
3064: if (numSecFields) { /* count the dofs for just this field */
3065: for (i = 0, numChildDof = 0; i < numChildren; i++) {
3066: PetscInt child = children[i];
3067: PetscInt dof;
3069: PetscSectionGetFieldDof(section,child,f,&dof);
3070: numChildDof += dof;
3071: }
3072: PetscSectionGetFieldDof(section,p,f,&numSelfDof);
3073: PetscSectionGetFieldOffset(section,p,f,&selfOff);
3074: }
3075: else {
3076: PetscSectionGetOffset(section,p,&selfOff);
3077: }
3078: for (i = 0; i < numSelfDof; i++) {
3079: nnz[selfOff + i] = numChildDof;
3080: }
3081: }
3082: }
3083: MatCreateAIJ(PETSC_COMM_SELF,m,n,m,n,-1,nnz,-1,NULL,&mat);
3084: PetscFree(nnz);
3085: /* Setp 2: compute entries */
3086: for (p = pStart; p < pEnd; p++) {
3087: const PetscInt *children;
3088: PetscInt numChildren;
3089: PetscInt i, numChildDof, numSelfDof;
3091: /* same conditions about when entries occur */
3092: if (canonical) {
3093: PetscInt pCanonical;
3094: DMLabelGetValue(canonical,p,&pCanonical);
3095: if (p != pCanonical) continue;
3096: }
3097: DMPlexGetTreeChildren(refTree,p,&numChildren,&children);
3098: if (!numChildren) continue;
3099: for (i = 0, numChildDof = 0; i < numChildren; i++) {
3100: PetscInt child = children[i];
3101: PetscInt dof;
3103: PetscSectionGetDof(section,child,&dof);
3104: numChildDof += dof;
3105: }
3106: PetscSectionGetDof(section,p,&numSelfDof);
3107: if (!numChildDof || !numSelfDof) continue;
3109: for (f = 0; f < numFields; f++) {
3110: PetscInt selfOff, fComp, numSelfShapes, numChildShapes, parentCell;
3111: PetscInt cellShapeOff;
3112: PetscObject disc;
3113: PetscDualSpace dsp;
3114: PetscClassId classId;
3115: PetscScalar *pointMat;
3116: PetscInt *matRows, *matCols;
3117: PetscInt pO = PETSC_MIN_INT;
3118: const PetscInt *depthNumDof;
3120: if (numSecFields) {
3121: for (i = 0, numChildDof = 0; i < numChildren; i++) {
3122: PetscInt child = children[i];
3123: PetscInt dof;
3125: PetscSectionGetFieldDof(section,child,f,&dof);
3126: numChildDof += dof;
3127: }
3128: PetscSectionGetFieldDof(section,p,f,&numSelfDof);
3129: PetscSectionGetFieldOffset(section,p,f,&selfOff);
3130: PetscSectionGetFieldComponents(section,f,&fComp);
3131: }
3132: else {
3133: PetscSectionGetOffset(section,p,&selfOff);
3134: fComp = 1;
3135: }
3136: numSelfShapes = numSelfDof / fComp;
3137: numChildShapes = numChildDof / fComp;
3139: /* find a cell whose closure contains p */
3140: if (p >= cStart && p < cEnd) {
3141: parentCell = p;
3142: }
3143: else {
3144: PetscInt *star = NULL;
3145: PetscInt numStar;
3147: parentCell = -1;
3148: DMPlexGetTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);
3149: for (i = numStar - 1; i >= 0; i--) {
3150: PetscInt c = star[2 * i];
3152: if (c >= cStart && c < cEnd) {
3153: parentCell = c;
3154: break;
3155: }
3156: }
3157: DMPlexRestoreTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);
3158: }
3159: /* determine the offset of p's shape functions withing parentCell's shape functions */
3160: PetscDSGetDiscretization(ds,f,&disc);
3161: PetscObjectGetClassId(disc,&classId);
3162: if (classId == PETSCFE_CLASSID) {
3163: PetscFEGetDualSpace((PetscFE)disc,&dsp);
3164: }
3165: else if (classId == PETSCFV_CLASSID) {
3166: PetscFVGetDualSpace((PetscFV)disc,&dsp);
3167: }
3168: else {
3169: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported discretization object");
3170: }
3171: PetscDualSpaceGetNumDof(dsp,&depthNumDof);
3172: {
3173: PetscInt *closure = NULL;
3174: PetscInt numClosure;
3176: DMPlexGetTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);
3177: for (i = 0, cellShapeOff = 0; i < numClosure; i++) {
3178: PetscInt point = closure[2 * i], pointDepth;
3180: pO = closure[2 * i + 1];
3181: if (point == p) break;
3182: DMLabelGetValue(depth,point,&pointDepth);
3183: cellShapeOff += depthNumDof[pointDepth];
3184: }
3185: DMPlexRestoreTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);
3186: }
3188: DMGetWorkArray(refTree, numSelfShapes * numChildShapes, PETSC_SCALAR,&pointMat);
3189: DMGetWorkArray(refTree, numSelfShapes + numChildShapes, PETSC_INT,&matRows);
3190: matCols = matRows + numSelfShapes;
3191: for (i = 0; i < numSelfShapes; i++) {
3192: matRows[i] = selfOff + i * fComp;
3193: }
3194: {
3195: PetscInt colOff = 0;
3197: for (i = 0; i < numChildren; i++) {
3198: PetscInt child = children[i];
3199: PetscInt dof, off, j;
3201: if (numSecFields) {
3202: PetscSectionGetFieldDof(cSection,child,f,&dof);
3203: PetscSectionGetFieldOffset(cSection,child,f,&off);
3204: }
3205: else {
3206: PetscSectionGetDof(cSection,child,&dof);
3207: PetscSectionGetOffset(cSection,child,&off);
3208: }
3210: for (j = 0; j < dof / fComp; j++) {
3211: matCols[colOff++] = off + j * fComp;
3212: }
3213: }
3214: }
3215: if (classId == PETSCFE_CLASSID) {
3216: PetscFE fe = (PetscFE) disc;
3217: PetscInt fSize;
3219: PetscFEGetDualSpace(fe,&dsp);
3220: PetscDualSpaceGetDimension(dsp,&fSize);
3221: for (i = 0; i < numSelfShapes; i++) { /* for every shape function */
3222: PetscQuadrature q;
3223: PetscInt dim, numPoints, j, k;
3224: const PetscReal *points;
3225: const PetscReal *weights;
3226: PetscInt *closure = NULL;
3227: PetscInt numClosure;
3228: PetscInt parentCellShapeDof = cellShapeOff + (pO < 0 ? (numSelfShapes - 1 - i) : i);
3229: PetscReal *Bparent;
3231: PetscDualSpaceGetFunctional(dsp,parentCellShapeDof,&q);
3232: PetscQuadratureGetData(q,&dim,&numPoints,&points,&weights);
3233: PetscFEGetTabulation(fe,numPoints,points,&Bparent,NULL,NULL); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */
3234: for (k = 0; k < numChildShapes; k++) {
3235: pointMat[numChildShapes * i + k] = 0.;
3236: }
3237: for (j = 0; j < numPoints; j++) {
3238: PetscInt childCell = -1;
3239: PetscReal parentValAtPoint;
3240: const PetscReal *pointReal = &points[dim * j];
3241: const PetscScalar *point;
3242: PetscReal *Bchild;
3243: PetscInt childCellShapeOff, pointMatOff;
3244: #if defined(PETSC_USE_COMPLEX)
3245: PetscInt d;
3247: for (d = 0; d < dim; d++) {
3248: pointScalar[d] = points[dim * j + d];
3249: }
3250: point = pointScalar;
3251: #else
3252: point = pointReal;
3253: #endif
3255: parentValAtPoint = Bparent[(fSize * j + parentCellShapeDof) * fComp];
3257: for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/
3258: PetscInt child = children[k];
3259: PetscInt *star = NULL;
3260: PetscInt numStar, s;
3262: DMPlexGetTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);
3263: for (s = numStar - 1; s >= 0; s--) {
3264: PetscInt c = star[2 * s];
3266: if (c < cStart || c >= cEnd) continue;
3267: DMPlexLocatePoint_Internal(refTree,dim,point,c,&childCell);
3268: if (childCell >= 0) break;
3269: }
3270: DMPlexRestoreTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);
3271: if (childCell >= 0) break;
3272: }
3273: if (childCell < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Could not locate quadrature point");
3274: DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ);
3275: DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent);
3276: CoordinatesRefToReal(dim, dim, v0parent, Jparent, pointReal, vtmp);
3277: CoordinatesRealToRef(dim, dim, v0, invJ, vtmp, pointRef);
3279: PetscFEGetTabulation(fe,1,pointRef,&Bchild,NULL,NULL);
3280: DMPlexGetTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);
3281: for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */
3282: PetscInt child = children[k], childDepth, childDof, childO = PETSC_MIN_INT;
3283: PetscInt l;
3285: DMLabelGetValue(depth,child,&childDepth);
3286: childDof = depthNumDof[childDepth];
3287: for (l = 0, childCellShapeOff = 0; l < numClosure; l++) {
3288: PetscInt point = closure[2 * l];
3289: PetscInt pointDepth;
3291: childO = closure[2 * l + 1];
3292: if (point == child) break;
3293: DMLabelGetValue(depth,point,&pointDepth);
3294: childCellShapeOff += depthNumDof[pointDepth];
3295: }
3296: if (l == numClosure) {
3297: pointMatOff += childDof;
3298: continue; /* child is not in the closure of the cell: has nothing to contribute to this point */
3299: }
3300: for (l = 0; l < childDof; l++) {
3301: PetscInt childCellDof = childCellShapeOff + (childO ? (childDof - 1 - l) : l);
3302: PetscReal childValAtPoint = Bchild[childCellDof * fComp];
3304: pointMat[i * numChildShapes + pointMatOff + l] += weights[j] * parentValAtPoint * childValAtPoint;
3305: }
3306: pointMatOff += childDof;
3307: }
3308: DMPlexRestoreTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);
3309: PetscFERestoreTabulation(fe,1,pointRef,&Bchild,NULL,NULL);
3310: }
3311: PetscFERestoreTabulation(fe,numPoints,points,&Bparent,NULL,NULL);
3312: }
3313: }
3314: else { /* just the volume-weighted averages of the children */
3315: PetscInt childShapeOff;
3316: PetscReal parentVol;
3318: DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL);
3319: for (i = 0, childShapeOff = 0; i < numChildren; i++) {
3320: PetscInt child = children[i];
3321: PetscReal childVol;
3323: if (child < cStart || child >= cEnd) continue;
3324: DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL);
3325: pointMat[childShapeOff] = childVol / parentVol;
3326: childShapeOff++;
3327: }
3328: }
3329: /* Insert pointMat into mat */
3330: for (i = 0; i < fComp; i++) {
3331: PetscInt j;
3332: MatSetValues(mat,numSelfShapes,matRows,numChildShapes,matCols,pointMat,INSERT_VALUES);
3334: for (j = 0; j < numSelfShapes; j++) {
3335: matRows[j]++;
3336: }
3337: for (j = 0; j < numChildShapes; j++) {
3338: matCols[j]++;
3339: }
3340: }
3341: DMRestoreWorkArray(refTree, numSelfShapes + numChildShapes, PETSC_INT,&matRows);
3342: DMRestoreWorkArray(refTree, numSelfShapes * numChildShapes, PETSC_SCALAR,&pointMat);
3343: }
3344: }
3345: PetscFree6(v0,v0parent,vtmp,J,Jparent,invJ);
3346: PetscFree2(pointScalar,pointRef);
3347: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
3348: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
3349: *inj = mat;
3350: return(0);
3351: }
3355: static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3356: {
3357: PetscDS ds;
3358: PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof;
3359: PetscScalar ***refPointFieldMats;
3360: PetscSection refConSec, refSection;
3364: DMGetDS(refTree,&ds);
3365: PetscDSGetNumFields(ds,&numFields);
3366: DMGetDefaultConstraints(refTree,&refConSec,NULL);
3367: DMGetDefaultSection(refTree,&refSection);
3368: PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
3369: PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);
3370: PetscSectionGetMaxDof(refConSec,&maxDof);
3371: PetscMalloc1(maxDof,&rows);
3372: PetscMalloc1(maxDof*maxDof,&cols);
3373: for (p = pRefStart; p < pRefEnd; p++) {
3374: PetscInt parent, pDof, parentDof;
3376: DMPlexGetTreeParent(refTree,p,&parent,NULL);
3377: PetscSectionGetDof(refConSec,p,&pDof);
3378: PetscSectionGetDof(refSection,parent,&parentDof);
3379: if (!pDof || !parentDof || parent == p) continue;
3381: PetscMalloc1(numFields,&refPointFieldMats[p-pRefStart]);
3382: for (f = 0; f < numFields; f++) {
3383: PetscInt cDof, cOff, numCols, r, fComp;
3384: PetscObject disc;
3385: PetscClassId id;
3386: PetscFE fe = NULL;
3387: PetscFV fv = NULL;
3389: PetscDSGetDiscretization(ds,f,&disc);
3390: PetscObjectGetClassId(disc,&id);
3391: if (id == PETSCFE_CLASSID) {
3392: fe = (PetscFE) disc;
3393: PetscFEGetNumComponents(fe,&fComp);
3394: }
3395: else if (id == PETSCFV_CLASSID) {
3396: fv = (PetscFV) disc;
3397: PetscFVGetNumComponents(fv,&fComp);
3398: }
3399: else SETERRQ1(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id);
3401: if (numFields > 1) {
3402: PetscSectionGetFieldDof(refConSec,p,f,&cDof);
3403: PetscSectionGetFieldOffset(refConSec,p,f,&cOff);
3404: }
3405: else {
3406: PetscSectionGetDof(refConSec,p,&cDof);
3407: PetscSectionGetOffset(refConSec,p,&cOff);
3408: }
3410: if (!cDof) continue;
3411: for (r = 0; r < cDof; r++) {
3412: rows[r] = cOff + r;
3413: }
3414: numCols = 0;
3415: {
3416: PetscInt aDof, aOff, j;
3418: if (numFields > 1) {
3419: PetscSectionGetFieldDof(refSection,parent,f,&aDof);
3420: PetscSectionGetFieldOffset(refSection,parent,f,&aOff);
3421: }
3422: else {
3423: PetscSectionGetDof(refSection,parent,&aDof);
3424: PetscSectionGetOffset(refSection,parent,&aOff);
3425: }
3427: for (j = 0; j < aDof; j++) {
3428: cols[numCols++] = aOff + j;
3429: }
3430: }
3431: PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);
3432: /* transpose of constraint matrix */
3433: MatGetValues(inj,numCols,cols,cDof,rows,refPointFieldMats[p-pRefStart][f]);
3434: }
3435: }
3436: *childrenMats = refPointFieldMats;
3437: PetscFree(rows);
3438: PetscFree(cols);
3439: return(0);
3440: }
3444: static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats)
3445: {
3446: PetscDS ds;
3447: PetscScalar ***refPointFieldMats;
3448: PetscInt numFields, pRefStart, pRefEnd, p, f;
3449: PetscSection refConSec, refSection;
3453: refPointFieldMats = *childrenMats;
3454: *childrenMats = NULL;
3455: DMGetDS(refTree,&ds);
3456: DMGetDefaultSection(refTree,&refSection);
3457: PetscDSGetNumFields(ds,&numFields);
3458: DMGetDefaultConstraints(refTree,&refConSec,NULL);
3459: PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
3460: for (p = pRefStart; p < pRefEnd; p++) {
3461: PetscInt parent, pDof, parentDof;
3463: DMPlexGetTreeParent(refTree,p,&parent,NULL);
3464: PetscSectionGetDof(refConSec,p,&pDof);
3465: PetscSectionGetDof(refSection,parent,&parentDof);
3466: if (!pDof || !parentDof || parent == p) continue;
3468: for (f = 0; f < numFields; f++) {
3469: PetscInt cDof;
3471: if (numFields > 1) {
3472: PetscSectionGetFieldDof(refConSec,p,f,&cDof);
3473: }
3474: else {
3475: PetscSectionGetDof(refConSec,p,&cDof);
3476: }
3478: if (!cDof) continue;
3479: PetscFree(refPointFieldMats[p - pRefStart][f]);
3480: }
3481: PetscFree(refPointFieldMats[p - pRefStart]);
3482: }
3483: PetscFree(refPointFieldMats);
3484: return(0);
3485: }
3489: PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat)
3490: {
3491: DM refTree;
3492: PetscSF coarseToFineEmbedded;
3493: const PetscInt *rootDegrees;
3494: PetscSection multiRootSec, rootIndicesSec, leafIndicesSec;
3495: PetscSection globalCoarse, globalFine;
3496: PetscSection localCoarse, localFine;
3497: PetscSection cSecRef;
3498: PetscInt *rootIndices, *leafIndices, *parentIndices, pRefStart, pRefEnd;
3499: Mat cMatRef, injRef;
3500: PetscInt numFields, numMulti, maxDof;
3501: PetscInt pStartC, pEndC, pStartF, pEndF, p;
3502: PetscInt *offsets, *offsetsCopy, *rowOffsets;
3503: PetscLayout rowMap, colMap;
3504: PetscInt rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO;
3505: PetscObject injRefObj;
3506: PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */
3511: /* get the templates for the fine-to-coarse injection from the reference tree */
3512: DMPlexGetReferenceTree(coarse,&refTree);
3513: DMGetDefaultConstraints(refTree,&cSecRef,&cMatRef);
3514: PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);
3515: PetscObjectQuery((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",&injRefObj);
3516: injRef = (Mat) injRefObj;
3517: if (!injRef) {
3518: DMPlexComputeInjectorReferenceTree(refTree,&injRef);
3519: PetscObjectCompose((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",(PetscObject)injRef);
3520: /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */
3521: PetscObjectDereference((PetscObject)injRef);
3522: }
3524: DMPlexGetChart(fine,&pStartF,&pEndF);
3525: DMGetDefaultSection(fine,&localFine);
3526: DMGetDefaultGlobalSection(fine,&globalFine);
3528: PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);
3529: PetscSectionSetChart(leafIndicesSec,pStartF, pEndF);
3531: PetscSectionGetNumFields(localFine,&numFields);
3532: {
3533: PetscInt maxFields = PetscMax(1,numFields) + 1;
3534: PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);
3535: }
3537: { /* winnow fine points that don't have global dofs out of the sf */
3538: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices;
3540: for (p = pStartF, numPointsWithDofs = 0; p < pEndF; p++) {
3541: PetscSectionGetDof(globalFine,p,&dof);
3542: PetscSectionGetConstraintDof(globalFine,p,&cdof);
3543: if ((dof - cdof) > 0) {
3544: numPointsWithDofs++;
3546: PetscSectionGetDof(localFine,p,&dof);
3547: PetscSectionSetDof(leafIndicesSec,p,dof + 1);
3548: }
3549: }
3550: PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
3551: PetscSectionSetUp(leafIndicesSec);
3552: PetscSectionGetStorageSize(leafIndicesSec,&numIndices);
3553: PetscMalloc1(numIndices,&leafIndices);
3554: for (p = pStartF, offset = 0; p < pEndF; p++) {
3555: PetscSectionGetDof(globalFine,p,&dof);
3556: PetscSectionGetConstraintDof(globalFine,p,&cdof);
3557: if ((dof - cdof) > 0) {
3558: PetscInt off, gOff;
3559: PetscInt *pInd;
3561: pointsWithDofs[offset++] = p - pStartF;
3563: PetscSectionGetOffset(leafIndicesSec,p,&off);
3565: pInd = &leafIndices[off + 1];
3566: PetscSectionGetOffset(globalFine,p,&gOff);
3568: offsets[0] = 0;
3569: if (numFields) {
3570: PetscInt f;
3572: for (f = 0; f < numFields; f++) {
3573: PetscInt fDof;
3574: PetscSectionGetFieldDof(localFine,p,f,&fDof);
3575: offsets[f + 1] = fDof + offsets[f];
3576: }
3577: indicesPointFields_private(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,0,pInd);
3578: }
3579: else {
3580: indicesPoint_private(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,0,pInd);
3581: }
3582: }
3583: }
3584: PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);
3585: PetscFree(pointsWithDofs);
3586: }
3588: DMPlexGetChart(coarse,&pStartC,&pEndC);
3589: DMGetDefaultSection(coarse,&localCoarse);
3590: DMGetDefaultGlobalSection(coarse,&globalCoarse);
3591: PetscSectionGetMaxDof(localCoarse,&maxDof);
3593: { /* there may be the case where an sf root has a parent: broadcast parents back to children */
3594: MPI_Datatype threeInt;
3595: PetscMPIInt rank;
3596: PetscInt (*parentNodeAndIdCoarse)[3];
3597: PetscInt (*parentNodeAndIdFine)[3];
3598: PetscInt p, nleaves, nleavesToParents;
3599: PetscSF pointSF, sfToParents;
3600: const PetscInt *ilocal;
3601: const PetscSFNode *iremote;
3602: PetscSFNode *iremoteToParents;
3603: PetscInt *ilocalToParents;
3605: MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);
3606: MPI_Type_contiguous(3,MPIU_INT,&threeInt);
3607: MPI_Type_commit(&threeInt);
3608: PetscMalloc2(pEndC-pStartC,&parentNodeAndIdCoarse,pEndF-pStartF,&parentNodeAndIdFine);
3609: DMGetPointSF(coarse,&pointSF);
3610: PetscSFGetGraph(pointSF,NULL,&nleaves,&ilocal,&iremote);
3611: for (p = pStartC; p < pEndC; p++) {
3612: PetscInt parent, childId;
3613: DMPlexGetTreeParent(coarse,p,&parent,&childId);
3614: parentNodeAndIdCoarse[p - pStartC][0] = rank;
3615: parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC;
3616: parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId;
3617: if (nleaves > 0) {
3618: PetscInt leaf = -1;
3620: if (ilocal) {
3621: PetscFindInt(parent,nleaves,ilocal,&leaf);
3622: }
3623: else {
3624: leaf = p - pStartC;
3625: }
3626: if (leaf >= 0) {
3627: parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank;
3628: parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index;
3629: }
3630: }
3631: }
3632: for (p = pStartF; p < pEndF; p++) {
3633: parentNodeAndIdFine[p - pStartF][0] = -1;
3634: parentNodeAndIdFine[p - pStartF][1] = -1;
3635: parentNodeAndIdFine[p - pStartF][2] = -1;
3636: }
3637: PetscSFBcastBegin(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);
3638: PetscSFBcastEnd(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);
3639: for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3640: PetscInt dof;
3642: PetscSectionGetDof(leafIndicesSec,p,&dof);
3643: if (dof) {
3644: PetscInt off;
3646: PetscSectionGetOffset(leafIndicesSec,p,&off);
3647: leafIndices[off] = PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]);
3648: }
3649: if (parentNodeAndIdFine[p-pStartF][0] >= 0) {
3650: nleavesToParents++;
3651: }
3652: }
3653: PetscMalloc1(nleavesToParents,&ilocalToParents);
3654: PetscMalloc1(nleavesToParents,&iremoteToParents);
3655: for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) {
3656: if (parentNodeAndIdFine[p-pStartF][0] >= 0) {
3657: ilocalToParents[nleavesToParents] = p - pStartF;
3658: iremoteToParents[nleavesToParents].rank = parentNodeAndIdFine[p-pStartF][0];
3659: iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p-pStartF][1];
3660: nleavesToParents++;
3661: }
3662: }
3663: PetscSFCreate(PetscObjectComm((PetscObject)coarse),&sfToParents);
3664: PetscSFSetGraph(sfToParents,pEndC-pStartC,nleavesToParents,ilocalToParents,PETSC_OWN_POINTER,iremoteToParents,PETSC_OWN_POINTER);
3665: PetscSFDestroy(&coarseToFineEmbedded);
3667: coarseToFineEmbedded = sfToParents;
3669: PetscFree2(parentNodeAndIdCoarse,parentNodeAndIdFine);
3670: MPI_Type_free(&threeInt);
3671: }
3673: { /* winnow out coarse points that don't have dofs */
3674: PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs;
3675: PetscSF sfDofsOnly;
3677: for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) {
3678: PetscSectionGetDof(globalCoarse,p,&dof);
3679: PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3680: if ((dof - cdof) > 0) {
3681: numPointsWithDofs++;
3682: }
3683: }
3684: PetscMalloc1(numPointsWithDofs,&pointsWithDofs);
3685: for (p = pStartC, offset = 0; p < pEndC; p++) {
3686: PetscSectionGetDof(globalCoarse,p,&dof);
3687: PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3688: if ((dof - cdof) > 0) {
3689: pointsWithDofs[offset++] = p - pStartC;
3690: }
3691: }
3692: PetscSFCreateEmbeddedSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly);
3693: PetscSFDestroy(&coarseToFineEmbedded);
3694: PetscFree(pointsWithDofs);
3695: coarseToFineEmbedded = sfDofsOnly;
3696: }
3698: /* communicate back to the coarse mesh which coarse points have children (that may require injection) */
3699: PetscSFComputeDegreeBegin(coarseToFineEmbedded,&rootDegrees);
3700: PetscSFComputeDegreeEnd(coarseToFineEmbedded,&rootDegrees);
3701: PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&multiRootSec);
3702: PetscSectionSetChart(multiRootSec,pStartC,pEndC);
3703: for (p = pStartC; p < pEndC; p++) {
3704: PetscSectionSetDof(multiRootSec,p,rootDegrees[p-pStartC]);
3705: }
3706: PetscSectionSetUp(multiRootSec);
3707: PetscSectionGetStorageSize(multiRootSec,&numMulti);
3708: PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);
3709: { /* distribute the leaf section */
3710: PetscSF multi, multiInv, indicesSF;
3711: PetscInt *remoteOffsets, numRootIndices;
3713: PetscSFGetMultiSF(coarseToFineEmbedded,&multi);
3714: PetscSFCreateInverseSF(multi,&multiInv);
3715: PetscSFDistributeSection(multiInv,leafIndicesSec,&remoteOffsets,rootIndicesSec);
3716: PetscSFCreateSectionSF(multiInv,leafIndicesSec,remoteOffsets,rootIndicesSec,&indicesSF);
3717: PetscFree(remoteOffsets);
3718: PetscSFDestroy(&multiInv);
3719: PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);
3720: PetscMalloc1(numRootIndices,&rootIndices);
3721: PetscSFBcastBegin(indicesSF,MPIU_INT,leafIndices,rootIndices);
3722: PetscSFBcastEnd(indicesSF,MPIU_INT,leafIndices,rootIndices);
3723: PetscSFDestroy(&indicesSF);
3724: }
3725: PetscSectionDestroy(&leafIndicesSec);
3726: PetscFree(leafIndices);
3727: PetscSFDestroy(&coarseToFineEmbedded);
3729: PetscMalloc1(maxDof,&parentIndices);
3731: /* count indices */
3732: MatGetLayouts(mat,&rowMap,&colMap);
3733: PetscLayoutSetUp(rowMap);
3734: PetscLayoutSetUp(colMap);
3735: PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);
3736: PetscLayoutGetRange(colMap,&colStart,&colEnd);
3737: PetscCalloc2(rowEnd-rowStart,&nnzD,rowEnd-rowStart,&nnzO);
3738: for (p = pStartC; p < pEndC; p++) {
3739: PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
3741: PetscSectionGetDof(globalCoarse,p,&dof);
3742: PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3743: if ((dof - cdof) <= 0) continue;
3744: PetscSectionGetOffset(globalCoarse,p,&gOff);
3746: rowOffsets[0] = 0;
3747: offsetsCopy[0] = 0;
3748: if (numFields) {
3749: PetscInt f;
3751: for (f = 0; f < numFields; f++) {
3752: PetscInt fDof;
3753: PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
3754: rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3755: }
3756: indicesPointFields_private(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,0,parentIndices);
3757: }
3758: else {
3759: indicesPoint_private(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,0,parentIndices);
3760: rowOffsets[1] = offsetsCopy[0];
3761: }
3763: PetscSectionGetDof(multiRootSec,p,&numLeaves);
3764: PetscSectionGetOffset(multiRootSec,p,&leafStart);
3765: leafEnd = leafStart + numLeaves;
3766: for (l = leafStart; l < leafEnd; l++) {
3767: PetscInt numIndices, childId, offset;
3768: const PetscInt *childIndices;
3770: PetscSectionGetDof(rootIndicesSec,l,&numIndices);
3771: PetscSectionGetOffset(rootIndicesSec,l,&offset);
3772: childId = rootIndices[offset++];
3773: childIndices = &rootIndices[offset];
3774: numIndices--;
3776: if (childId == -1) { /* equivalent points: scatter */
3777: PetscInt i;
3779: for (i = 0; i < numIndices; i++) {
3780: PetscInt colIndex = childIndices[i];
3781: PetscInt rowIndex = parentIndices[i];
3782: if (rowIndex < 0) continue;
3783: if (colIndex < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unconstrained fine and constrained coarse");
3784: if (colIndex >= colStart && colIndex < colEnd) {
3785: nnzD[rowIndex - rowStart] = 1;
3786: }
3787: else {
3788: nnzO[rowIndex - rowStart] = 1;
3789: }
3790: }
3791: }
3792: else {
3793: PetscInt parentId, f, lim;
3795: DMPlexGetTreeParent(refTree,childId,&parentId,NULL);
3797: lim = PetscMax(1,numFields);
3798: offsets[0] = 0;
3799: if (numFields) {
3800: PetscInt f;
3802: for (f = 0; f < numFields; f++) {
3803: PetscInt fDof;
3804: PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);
3806: offsets[f + 1] = fDof + offsets[f];
3807: }
3808: }
3809: else {
3810: PetscInt cDof;
3812: PetscSectionGetDof(cSecRef,childId,&cDof);
3813: offsets[1] = cDof;
3814: }
3815: for (f = 0; f < lim; f++) {
3816: PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1];
3817: PetscInt childStart = offsets[f], childEnd = offsets[f + 1];
3818: PetscInt i, numD = 0, numO = 0;
3820: for (i = childStart; i < childEnd; i++) {
3821: PetscInt colIndex = childIndices[i];
3823: if (colIndex < 0) continue;
3824: if (colIndex >= colStart && colIndex < colEnd) {
3825: numD++;
3826: }
3827: else {
3828: numO++;
3829: }
3830: }
3831: for (i = parentStart; i < parentEnd; i++) {
3832: PetscInt rowIndex = parentIndices[i];
3834: if (rowIndex < 0) continue;
3835: nnzD[rowIndex - rowStart] += numD;
3836: nnzO[rowIndex - rowStart] += numO;
3837: }
3838: }
3839: }
3840: }
3841: }
3842: /* preallocate */
3843: MatXAIJSetPreallocation(mat,1,nnzD,nnzO,NULL,NULL);
3844: PetscFree2(nnzD,nnzO);
3845: /* insert values */
3846: DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);
3847: for (p = pStartC; p < pEndC; p++) {
3848: PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff;
3850: PetscSectionGetDof(globalCoarse,p,&dof);
3851: PetscSectionGetConstraintDof(globalCoarse,p,&cdof);
3852: if ((dof - cdof) <= 0) continue;
3853: PetscSectionGetOffset(globalCoarse,p,&gOff);
3855: rowOffsets[0] = 0;
3856: offsetsCopy[0] = 0;
3857: if (numFields) {
3858: PetscInt f;
3860: for (f = 0; f < numFields; f++) {
3861: PetscInt fDof;
3862: PetscSectionGetFieldDof(localCoarse,p,f,&fDof);
3863: rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f];
3864: }
3865: indicesPointFields_private(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,0,parentIndices);
3866: }
3867: else {
3868: indicesPoint_private(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,0,parentIndices);
3869: rowOffsets[1] = offsetsCopy[0];
3870: }
3872: PetscSectionGetDof(multiRootSec,p,&numLeaves);
3873: PetscSectionGetOffset(multiRootSec,p,&leafStart);
3874: leafEnd = leafStart + numLeaves;
3875: for (l = leafStart; l < leafEnd; l++) {
3876: PetscInt numIndices, childId, offset;
3877: const PetscInt *childIndices;
3879: PetscSectionGetDof(rootIndicesSec,l,&numIndices);
3880: PetscSectionGetOffset(rootIndicesSec,l,&offset);
3881: childId = rootIndices[offset++];
3882: childIndices = &rootIndices[offset];
3883: numIndices--;
3885: if (childId == -1) { /* equivalent points: scatter */
3886: PetscInt i;
3888: for (i = 0; i < numIndices; i++) {
3889: MatSetValue(mat,parentIndices[i],childIndices[i],1.,INSERT_VALUES);
3890: }
3891: }
3892: else {
3893: PetscInt parentId, f, lim;
3895: DMPlexGetTreeParent(refTree,childId,&parentId,NULL);
3897: lim = PetscMax(1,numFields);
3898: offsets[0] = 0;
3899: if (numFields) {
3900: PetscInt f;
3902: for (f = 0; f < numFields; f++) {
3903: PetscInt fDof;
3904: PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);
3906: offsets[f + 1] = fDof + offsets[f];
3907: }
3908: }
3909: else {
3910: PetscInt cDof;
3912: PetscSectionGetDof(cSecRef,childId,&cDof);
3913: offsets[1] = cDof;
3914: }
3915: for (f = 0; f < lim; f++) {
3916: PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0];
3917: PetscInt *rowIndices = &parentIndices[rowOffsets[f]];
3918: const PetscInt *colIndices = &childIndices[offsets[f]];
3920: MatSetValues(mat,rowOffsets[f+1]-rowOffsets[f],rowIndices,offsets[f+1]-offsets[f],colIndices,childMat,INSERT_VALUES);
3921: }
3922: }
3923: }
3924: }
3925: PetscSectionDestroy(&multiRootSec);
3926: PetscSectionDestroy(&rootIndicesSec);
3927: PetscFree(parentIndices);
3928: DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);
3929: PetscFree(rootIndices);
3930: PetscFree3(offsets,offsetsCopy,rowOffsets);
3932: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
3933: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
3934: return(0);
3935: }