Actual source code: plextree.c
petsc-3.6.1 2015-08-06
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: }
151: else {
152: oAvert = parentOrientA;
153: oBvert = parentOrientB;
154: ABswapVert = ABswap;
155: }
156: if (childB) {
157: /* assume that each child corresponds to a vertex, in the same order */
158: PetscInt p, posA = -1, numChildren, i;
159: const PetscInt *children;
161: /* count which position the child is in */
162: DMPlexGetTreeChildren(dm,parent,&numChildren,&children);
163: for (i = 0; i < numChildren; i++) {
164: p = children[i];
165: if (p == childA) {
166: posA = i;
167: break;
168: }
169: }
170: if (posA >= coneSize) {
171: /* this is the triangle in the middle of a uniformly refined triangle: it
172: * is invariant */
173: if (dim != 2 || posA != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Expected a middle triangle, got something else");
174: *childB = childA;
175: }
176: else {
177: /* figure out position B by applying ABswapVert */
178: PetscInt posB;
180: posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize -(ABswapVert + 1) - posA) % coneSize);
181: if (childB) *childB = children[posB];
182: }
183: }
184: if (childOrientB) *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap);
185: return(0);
186: }
190: /*@
191: DMPlexReferenceTreeGetChildSymmetry - Given a reference tree, transform a childid and orientation from one parent frame to another
193: Input Parameters:
194: + dm - the reference tree DMPlex object
195: . parent - the parent point
196: . parentOrientA - the reference orientation for describing the parent
197: . childOrientA - the reference orientation for describing the child
198: . childA - the reference childID for describing the child
199: - parentOrientB - the new orientation for describing the parent
201: Output Parameters:
202: + childOrientB - if not NULL, set to the new oreintation for describing the child
203: . childB - if not NULL, the new childID for describing the child
205: Level: developer
207: .seealso: DMPlexGetReferenceTree(), DMPlexSetReferenceTree(), DMPlexSetTree()
208: @*/
209: PetscErrorCode DMPlexReferenceTreeGetChildSymmetry(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB)
210: {
211: DM_Plex *mesh = (DM_Plex *)dm->data;
216: if (!mesh->getchildsymmetry) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"DMPlexReferenceTreeGetChildSymmetry not implemented");
217: mesh->getchildsymmetry(dm,parent,parentOrientA,childOrientA,childA,parentOrientB,childOrientB,childB);
218: return(0);
219: }
221: static PetscErrorCode DMPlexSetTree_Internal(DM,PetscSection,PetscInt*,PetscInt*,PetscBool,PetscBool);
225: /*@
226: DMPlexCreateDefaultReferenceTree - create a reference tree for isotropic hierarchical mesh refinement.
228: Collective on comm
230: Input Parameters:
231: + comm - the MPI communicator
232: . dim - the spatial dimension
233: - simplex - Flag for simplex, otherwise use a tensor-product cell
235: Output Parameters:
236: . ref - the reference tree DMPlex object
238: Level: intermediate
240: .keywords: reference cell
241: .seealso: DMPlexSetReferenceTree(), DMPlexGetReferenceTree()
242: @*/
243: PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref)
244: {
245: DM_Plex *mesh;
246: DM K, Kref;
247: PetscInt p, pStart, pEnd, pRefStart, pRefEnd, d, offset, parentSize, *parents, *childIDs;
248: PetscInt *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts;
249: DMLabel identity, identityRef;
250: PetscSection unionSection, unionConeSection, parentSection;
251: PetscScalar *unionCoords;
252: IS perm;
256: #if 1
257: comm = PETSC_COMM_SELF;
258: #endif
259: /* create a reference element */
260: DMPlexCreateReferenceCell(comm, dim, simplex, &K);
261: DMPlexCreateLabel(K, "identity");
262: DMPlexGetLabel(K, "identity", &identity);
263: DMPlexGetChart(K, &pStart, &pEnd);
264: for (p = pStart; p < pEnd; p++) {
265: DMLabelSetValue(identity, p, p);
266: }
267: /* refine it */
268: DMRefine(K,comm,&Kref);
270: /* the reference tree is the union of these two, without duplicating
271: * points that appear in both */
272: DMPlexGetLabel(Kref, "identity", &identityRef);
273: DMPlexGetChart(Kref, &pRefStart, &pRefEnd);
274: PetscSectionCreate(comm, &unionSection);
275: PetscSectionSetChart(unionSection, 0, (pEnd - pStart) + (pRefEnd - pRefStart));
276: /* count points that will go in the union */
277: for (p = pStart; p < pEnd; p++) {
278: PetscSectionSetDof(unionSection, p - pStart, 1);
279: }
280: for (p = pRefStart; p < pRefEnd; p++) {
281: PetscInt q, qSize;
282: DMLabelGetValue(identityRef, p, &q);
283: DMLabelGetStratumSize(identityRef, q, &qSize);
284: if (qSize > 1) {
285: PetscSectionSetDof(unionSection, p - pRefStart + (pEnd - pStart), 1);
286: }
287: }
288: PetscMalloc1(pEnd - pStart + pRefEnd - pRefStart,&permvals);
289: offset = 0;
290: /* stratify points in the union by topological dimension */
291: for (d = 0; d <= dim; d++) {
292: PetscInt cStart, cEnd, c;
294: DMPlexGetHeightStratum(K, d, &cStart, &cEnd);
295: for (c = cStart; c < cEnd; c++) {
296: permvals[offset++] = c;
297: }
299: DMPlexGetHeightStratum(Kref, d, &cStart, &cEnd);
300: for (c = cStart; c < cEnd; c++) {
301: permvals[offset++] = c + (pEnd - pStart);
302: }
303: }
304: ISCreateGeneral(comm, (pEnd - pStart) + (pRefEnd - pRefStart), permvals, PETSC_OWN_POINTER, &perm);
305: PetscSectionSetPermutation(unionSection,perm);
306: PetscSectionSetUp(unionSection);
307: PetscSectionGetStorageSize(unionSection,&numUnionPoints);
308: PetscMalloc2(numUnionPoints,&coneSizes,dim+1,&numDimPoints);
309: /* count dimension points */
310: for (d = 0; d <= dim; d++) {
311: PetscInt cStart, cOff, cOff2;
312: DMPlexGetHeightStratum(K,d,&cStart,NULL);
313: PetscSectionGetOffset(unionSection,cStart-pStart,&cOff);
314: if (d < dim) {
315: DMPlexGetHeightStratum(K,d+1,&cStart,NULL);
316: PetscSectionGetOffset(unionSection,cStart-pStart,&cOff2);
317: }
318: else {
319: cOff2 = numUnionPoints;
320: }
321: numDimPoints[dim - d] = cOff2 - cOff;
322: }
323: PetscSectionCreate(comm, &unionConeSection);
324: PetscSectionSetChart(unionConeSection, 0, numUnionPoints);
325: /* count the cones in the union */
326: for (p = pStart; p < pEnd; p++) {
327: PetscInt dof, uOff;
329: DMPlexGetConeSize(K, p, &dof);
330: PetscSectionGetOffset(unionSection, p - pStart,&uOff);
331: PetscSectionSetDof(unionConeSection, uOff, dof);
332: coneSizes[uOff] = dof;
333: }
334: for (p = pRefStart; p < pRefEnd; p++) {
335: PetscInt dof, uDof, uOff;
337: DMPlexGetConeSize(Kref, p, &dof);
338: PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);
339: PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);
340: if (uDof) {
341: PetscSectionSetDof(unionConeSection, uOff, dof);
342: coneSizes[uOff] = dof;
343: }
344: }
345: PetscSectionSetUp(unionConeSection);
346: PetscSectionGetStorageSize(unionConeSection,&numCones);
347: PetscMalloc2(numCones,&unionCones,numCones,&unionOrientations);
348: /* write the cones in the union */
349: for (p = pStart; p < pEnd; p++) {
350: PetscInt dof, uOff, c, cOff;
351: const PetscInt *cone, *orientation;
353: DMPlexGetConeSize(K, p, &dof);
354: DMPlexGetCone(K, p, &cone);
355: DMPlexGetConeOrientation(K, p, &orientation);
356: PetscSectionGetOffset(unionSection, p - pStart,&uOff);
357: PetscSectionGetOffset(unionConeSection,uOff,&cOff);
358: for (c = 0; c < dof; c++) {
359: PetscInt e, eOff;
360: e = cone[c];
361: PetscSectionGetOffset(unionSection, e - pStart, &eOff);
362: unionCones[cOff + c] = eOff;
363: unionOrientations[cOff + c] = orientation[c];
364: }
365: }
366: for (p = pRefStart; p < pRefEnd; p++) {
367: PetscInt dof, uDof, uOff, c, cOff;
368: const PetscInt *cone, *orientation;
370: DMPlexGetConeSize(Kref, p, &dof);
371: DMPlexGetCone(Kref, p, &cone);
372: DMPlexGetConeOrientation(Kref, p, &orientation);
373: PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);
374: PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);
375: if (uDof) {
376: PetscSectionGetOffset(unionConeSection,uOff,&cOff);
377: for (c = 0; c < dof; c++) {
378: PetscInt e, eOff, eDof;
380: e = cone[c];
381: PetscSectionGetDof(unionSection, e - pRefStart + (pEnd - pStart),&eDof);
382: if (eDof) {
383: PetscSectionGetOffset(unionSection, e - pRefStart + (pEnd - pStart), &eOff);
384: }
385: else {
386: DMLabelGetValue(identityRef, e, &e);
387: PetscSectionGetOffset(unionSection, e - pStart, &eOff);
388: }
389: unionCones[cOff + c] = eOff;
390: unionOrientations[cOff + c] = orientation[c];
391: }
392: }
393: }
394: /* get the coordinates */
395: {
396: PetscInt vStart, vEnd, vRefStart, vRefEnd, v, vDof, vOff;
397: PetscSection KcoordsSec, KrefCoordsSec;
398: Vec KcoordsVec, KrefCoordsVec;
399: PetscScalar *Kcoords;
401: DMGetCoordinateSection(K, &KcoordsSec);
402: DMGetCoordinatesLocal(K, &KcoordsVec);
403: DMGetCoordinateSection(Kref, &KrefCoordsSec);
404: DMGetCoordinatesLocal(Kref, &KrefCoordsVec);
406: numVerts = numDimPoints[0];
407: PetscMalloc1(numVerts * dim,&unionCoords);
408: DMPlexGetDepthStratum(K,0,&vStart,&vEnd);
410: offset = 0;
411: for (v = vStart; v < vEnd; v++) {
412: PetscSectionGetOffset(unionSection,v - pStart,&vOff);
413: VecGetValuesSection(KcoordsVec, KcoordsSec, v, &Kcoords);
414: for (d = 0; d < dim; d++) {
415: unionCoords[offset * dim + d] = Kcoords[d];
416: }
417: offset++;
418: }
419: DMPlexGetDepthStratum(Kref,0,&vRefStart,&vRefEnd);
420: for (v = vRefStart; v < vRefEnd; v++) {
421: PetscSectionGetDof(unionSection,v - pRefStart + (pEnd - pStart),&vDof);
422: PetscSectionGetOffset(unionSection,v - pRefStart + (pEnd - pStart),&vOff);
423: VecGetValuesSection(KrefCoordsVec, KrefCoordsSec, v, &Kcoords);
424: if (vDof) {
425: for (d = 0; d < dim; d++) {
426: unionCoords[offset * dim + d] = Kcoords[d];
427: }
428: offset++;
429: }
430: }
431: }
432: DMCreate(comm,ref);
433: DMSetType(*ref,DMPLEX);
434: DMSetDimension(*ref,dim);
435: DMPlexCreateFromDAG(*ref,dim,numDimPoints,coneSizes,unionCones,unionOrientations,unionCoords);
436: /* set the tree */
437: PetscSectionCreate(comm,&parentSection);
438: PetscSectionSetChart(parentSection,0,numUnionPoints);
439: for (p = pRefStart; p < pRefEnd; p++) {
440: PetscInt uDof, uOff;
442: PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);
443: PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);
444: if (uDof) {
445: PetscSectionSetDof(parentSection,uOff,1);
446: }
447: }
448: PetscSectionSetUp(parentSection);
449: PetscSectionGetStorageSize(parentSection,&parentSize);
450: PetscMalloc2(parentSize,&parents,parentSize,&childIDs);
451: for (p = pRefStart; p < pRefEnd; p++) {
452: PetscInt uDof, uOff;
454: PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);
455: PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);
456: if (uDof) {
457: PetscInt pOff, parent, parentU;
458: PetscSectionGetOffset(parentSection,uOff,&pOff);
459: DMLabelGetValue(identityRef,p,&parent);
460: PetscSectionGetOffset(unionSection, parent - pStart,&parentU);
461: parents[pOff] = parentU;
462: childIDs[pOff] = uOff;
463: }
464: }
465: DMPlexSetTree_Internal(*ref,parentSection,parents,childIDs,PETSC_TRUE,PETSC_FALSE);
466: mesh = (DM_Plex *) (*ref)->data;
467: mesh->getchildsymmetry = DMPlexReferenceTreeGetChildSymmetry_Default;
468: PetscSectionDestroy(&parentSection);
469: PetscFree2(parents,childIDs);
471: /* clean up */
472: PetscSectionDestroy(&unionSection);
473: PetscSectionDestroy(&unionConeSection);
474: ISDestroy(&perm);
475: PetscFree(unionCoords);
476: PetscFree2(unionCones,unionOrientations);
477: PetscFree2(coneSizes,numDimPoints);
478: DMDestroy(&K);
479: DMDestroy(&Kref);
480: return(0);
481: }
485: static PetscErrorCode DMPlexTreeSymmetrize(DM dm)
486: {
487: DM_Plex *mesh = (DM_Plex *)dm->data;
488: PetscSection childSec, pSec;
489: PetscInt p, pSize, cSize, parMax = PETSC_MIN_INT, parMin = PETSC_MAX_INT;
490: PetscInt *offsets, *children, pStart, pEnd;
495: PetscSectionDestroy(&mesh->childSection);
496: PetscFree(mesh->children);
497: pSec = mesh->parentSection;
498: if (!pSec) return(0);
499: PetscSectionGetStorageSize(pSec,&pSize);
500: for (p = 0; p < pSize; p++) {
501: PetscInt par = mesh->parents[p];
503: parMax = PetscMax(parMax,par+1);
504: parMin = PetscMin(parMin,par);
505: }
506: if (parMin > parMax) {
507: parMin = -1;
508: parMax = -1;
509: }
510: PetscSectionCreate(PetscObjectComm((PetscObject)pSec),&childSec);
511: PetscSectionSetChart(childSec,parMin,parMax);
512: for (p = 0; p < pSize; p++) {
513: PetscInt par = mesh->parents[p];
515: PetscSectionAddDof(childSec,par,1);
516: }
517: PetscSectionSetUp(childSec);
518: PetscSectionGetStorageSize(childSec,&cSize);
519: PetscMalloc1(cSize,&children);
520: PetscCalloc1(parMax-parMin,&offsets);
521: PetscSectionGetChart(pSec,&pStart,&pEnd);
522: for (p = pStart; p < pEnd; p++) {
523: PetscInt dof, off, i;
525: PetscSectionGetDof(pSec,p,&dof);
526: PetscSectionGetOffset(pSec,p,&off);
527: for (i = 0; i < dof; i++) {
528: PetscInt par = mesh->parents[off + i], cOff;
530: PetscSectionGetOffset(childSec,par,&cOff);
531: children[cOff + offsets[par-parMin]++] = p;
532: }
533: }
534: mesh->childSection = childSec;
535: mesh->children = children;
536: PetscFree(offsets);
537: return(0);
538: }
542: static PetscErrorCode AnchorsFlatten (PetscSection section, IS is, PetscSection *sectionNew, IS *isNew)
543: {
544: PetscInt pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL;
545: const PetscInt *vals;
546: PetscSection secNew;
547: PetscBool anyNew, globalAnyNew;
548: PetscBool compress;
552: PetscSectionGetChart(section,&pStart,&pEnd);
553: ISGetLocalSize(is,&size);
554: ISGetIndices(is,&vals);
555: PetscSectionCreate(PetscObjectComm((PetscObject)section),&secNew);
556: PetscSectionSetChart(secNew,pStart,pEnd);
557: for (i = 0; i < size; i++) {
558: PetscInt dof;
560: p = vals[i];
561: if (p < pStart || p >= pEnd) continue;
562: PetscSectionGetDof(section, p, &dof);
563: if (dof) break;
564: }
565: if (i == size) {
566: PetscSectionSetUp(secNew);
567: anyNew = PETSC_FALSE;
568: compress = PETSC_FALSE;
569: sizeNew = 0;
570: }
571: else {
572: anyNew = PETSC_TRUE;
573: for (p = pStart; p < pEnd; p++) {
574: PetscInt dof, off;
576: PetscSectionGetDof(section, p, &dof);
577: PetscSectionGetOffset(section, p, &off);
578: for (i = 0; i < dof; i++) {
579: PetscInt q = vals[off + i], qDof = 0;
581: if (q >= pStart && q < pEnd) {
582: PetscSectionGetDof(section, q, &qDof);
583: }
584: if (qDof) {
585: PetscSectionAddDof(secNew, p, qDof);
586: }
587: else {
588: PetscSectionAddDof(secNew, p, 1);
589: }
590: }
591: }
592: PetscSectionSetUp(secNew);
593: PetscSectionGetStorageSize(secNew,&sizeNew);
594: PetscMalloc1(sizeNew,&valsNew);
595: compress = PETSC_FALSE;
596: for (p = pStart; p < pEnd; p++) {
597: PetscInt dof, off, count, offNew, dofNew;
599: PetscSectionGetDof(section, p, &dof);
600: PetscSectionGetOffset(section, p, &off);
601: PetscSectionGetDof(secNew, p, &dofNew);
602: PetscSectionGetOffset(secNew, p, &offNew);
603: count = 0;
604: for (i = 0; i < dof; i++) {
605: PetscInt q = vals[off + i], qDof = 0, qOff = 0, j;
607: if (q >= pStart && q < pEnd) {
608: PetscSectionGetDof(section, q, &qDof);
609: PetscSectionGetOffset(section, q, &qOff);
610: }
611: if (qDof) {
612: PetscInt oldCount = count;
614: for (j = 0; j < qDof; j++) {
615: PetscInt k, r = vals[qOff + j];
617: for (k = 0; k < oldCount; k++) {
618: if (valsNew[offNew + k] == r) {
619: break;
620: }
621: }
622: if (k == oldCount) {
623: valsNew[offNew + count++] = r;
624: }
625: }
626: }
627: else {
628: PetscInt k, oldCount = count;
630: for (k = 0; k < oldCount; k++) {
631: if (valsNew[offNew + k] == q) {
632: break;
633: }
634: }
635: if (k == oldCount) {
636: valsNew[offNew + count++] = q;
637: }
638: }
639: }
640: if (count < dofNew) {
641: PetscSectionSetDof(secNew, p, count);
642: compress = PETSC_TRUE;
643: }
644: }
645: }
646: ISRestoreIndices(is,&vals);
647: MPI_Allreduce(&anyNew,&globalAnyNew,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));
648: if (!globalAnyNew) {
649: PetscSectionDestroy(&secNew);
650: *sectionNew = NULL;
651: *isNew = NULL;
652: }
653: else {
654: PetscBool globalCompress;
656: MPI_Allreduce(&compress,&globalCompress,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));
657: if (compress) {
658: PetscSection secComp;
659: PetscInt *valsComp = NULL;
661: PetscSectionCreate(PetscObjectComm((PetscObject)section),&secComp);
662: PetscSectionSetChart(secComp,pStart,pEnd);
663: for (p = pStart; p < pEnd; p++) {
664: PetscInt dof;
666: PetscSectionGetDof(secNew, p, &dof);
667: PetscSectionSetDof(secComp, p, dof);
668: }
669: PetscSectionSetUp(secComp);
670: PetscSectionGetStorageSize(secComp,&sizeNew);
671: PetscMalloc1(sizeNew,&valsComp);
672: for (p = pStart; p < pEnd; p++) {
673: PetscInt dof, off, offNew, j;
675: PetscSectionGetDof(secNew, p, &dof);
676: PetscSectionGetOffset(secNew, p, &off);
677: PetscSectionGetOffset(secComp, p, &offNew);
678: for (j = 0; j < dof; j++) {
679: valsComp[offNew + j] = valsNew[off + j];
680: }
681: }
682: PetscSectionDestroy(&secNew);
683: secNew = secComp;
684: PetscFree(valsNew);
685: valsNew = valsComp;
686: }
687: ISCreateGeneral(PetscObjectComm((PetscObject)is),sizeNew,valsNew,PETSC_OWN_POINTER,isNew);
688: }
689: return(0);
690: }
694: static PetscErrorCode DMPlexCreateAnchors_Tree(DM dm)
695: {
696: PetscInt p, pStart, pEnd, *anchors, size;
697: PetscInt aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT;
698: PetscSection aSec;
699: DMLabel canonLabel;
700: IS aIS;
705: DMPlexGetChart(dm,&pStart,&pEnd);
706: DMPlexGetLabel(dm,"canonical",&canonLabel);
707: for (p = pStart; p < pEnd; p++) {
708: PetscInt parent;
710: if (canonLabel) {
711: PetscInt canon;
713: DMLabelGetValue(canonLabel,p,&canon);
714: if (p != canon) continue;
715: }
716: DMPlexGetTreeParent(dm,p,&parent,NULL);
717: if (parent != p) {
718: aMin = PetscMin(aMin,p);
719: aMax = PetscMax(aMax,p+1);
720: }
721: }
722: if (aMin > aMax) {
723: aMin = -1;
724: aMax = -1;
725: }
726: PetscSectionCreate(PETSC_COMM_SELF,&aSec);
727: PetscSectionSetChart(aSec,aMin,aMax);
728: for (p = aMin; p < aMax; p++) {
729: PetscInt parent, ancestor = p;
731: if (canonLabel) {
732: PetscInt canon;
734: DMLabelGetValue(canonLabel,p,&canon);
735: if (p != canon) continue;
736: }
737: DMPlexGetTreeParent(dm,p,&parent,NULL);
738: while (parent != ancestor) {
739: ancestor = parent;
740: DMPlexGetTreeParent(dm,ancestor,&parent,NULL);
741: }
742: if (ancestor != p) {
743: PetscInt closureSize, *closure = NULL;
745: DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
746: PetscSectionSetDof(aSec,p,closureSize);
747: DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
748: }
749: }
750: PetscSectionSetUp(aSec);
751: PetscSectionGetStorageSize(aSec,&size);
752: PetscMalloc1(size,&anchors);
753: for (p = aMin; p < aMax; p++) {
754: PetscInt parent, ancestor = p;
756: if (canonLabel) {
757: PetscInt canon;
759: DMLabelGetValue(canonLabel,p,&canon);
760: if (p != canon) continue;
761: }
762: DMPlexGetTreeParent(dm,p,&parent,NULL);
763: while (parent != ancestor) {
764: ancestor = parent;
765: DMPlexGetTreeParent(dm,ancestor,&parent,NULL);
766: }
767: if (ancestor != p) {
768: PetscInt j, closureSize, *closure = NULL, aOff;
770: PetscSectionGetOffset(aSec,p,&aOff);
772: DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
773: for (j = 0; j < closureSize; j++) {
774: anchors[aOff + j] = closure[2*j];
775: }
776: DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);
777: }
778: }
779: ISCreateGeneral(PETSC_COMM_SELF,size,anchors,PETSC_OWN_POINTER,&aIS);
780: {
781: PetscSection aSecNew = aSec;
782: IS aISNew = aIS;
784: PetscObjectReference((PetscObject)aSec);
785: PetscObjectReference((PetscObject)aIS);
786: while (aSecNew) {
787: PetscSectionDestroy(&aSec);
788: ISDestroy(&aIS);
789: aSec = aSecNew;
790: aIS = aISNew;
791: aSecNew = NULL;
792: aISNew = NULL;
793: AnchorsFlatten(aSec,aIS,&aSecNew,&aISNew);
794: }
795: }
796: DMPlexSetAnchors(dm,aSec,aIS);
797: PetscSectionDestroy(&aSec);
798: ISDestroy(&aIS);
799: return(0);
800: }
804: static PetscErrorCode DMPlexTreeExchangeSupports(DM dm)
805: {
806: DM_Plex *mesh = (DM_Plex *)dm->data;
807: PetscSection newSupportSection;
808: PetscInt newSize, *newSupports, pStart, pEnd, p, d, depth;
809: PetscInt *offsets;
814: /* symmetrize the hierarchy */
815: DMPlexGetDepth(dm,&depth);
816: PetscSectionCreate(PetscObjectComm((PetscObject)(mesh->supportSection)),&newSupportSection);
817: DMPlexGetChart(dm,&pStart,&pEnd);
818: PetscSectionSetChart(newSupportSection,pStart,pEnd);
819: PetscCalloc1(pEnd,&offsets);
820: /* if a point is in the support of q, it should be in the support of
821: * parent(q) */
822: for (d = 0; d <= depth; d++) {
823: DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);
824: for (p = pStart; p < pEnd; ++p) {
825: PetscInt dof, q, qdof, parent;
827: PetscSectionGetDof(mesh->supportSection, p, &dof);
828: PetscSectionAddDof(newSupportSection, p, dof);
829: q = p;
830: DMPlexGetTreeParent(dm,q,&parent,NULL);
831: while (parent != q && parent >= pStart && parent < pEnd) {
832: q = parent;
834: PetscSectionGetDof(mesh->supportSection, q, &qdof);
835: PetscSectionAddDof(newSupportSection,p,qdof);
836: PetscSectionAddDof(newSupportSection,q,dof);
837: DMPlexGetTreeParent(dm,q,&parent,NULL);
838: }
839: }
840: }
841: PetscSectionSetUp(newSupportSection);
842: PetscSectionGetStorageSize(newSupportSection,&newSize);
843: PetscMalloc1(newSize,&newSupports);
844: for (d = 0; d <= depth; d++) {
845: DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);
846: for (p = pStart; p < pEnd; p++) {
847: PetscInt dof, off, q, qdof, qoff, newDof, newOff, newqOff, i, parent;
849: PetscSectionGetDof(mesh->supportSection, p, &dof);
850: PetscSectionGetOffset(mesh->supportSection, p, &off);
851: PetscSectionGetDof(newSupportSection, p, &newDof);
852: PetscSectionGetOffset(newSupportSection, p, &newOff);
853: for (i = 0; i < dof; i++) {
854: newSupports[newOff+offsets[p]++] = mesh->supports[off + i];
855: }
856: mesh->maxSupportSize = PetscMax(mesh->maxSupportSize,newDof);
858: q = p;
859: DMPlexGetTreeParent(dm,q,&parent,NULL);
860: while (parent != q && parent >= pStart && parent < pEnd) {
861: q = parent;
862: PetscSectionGetDof(mesh->supportSection, q, &qdof);
863: PetscSectionGetOffset(mesh->supportSection, q, &qoff);
864: PetscSectionGetOffset(newSupportSection, q, &newqOff);
865: for (i = 0; i < qdof; i++) {
866: newSupports[newOff+offsets[p]++] = mesh->supports[qoff + i];
867: }
868: for (i = 0; i < dof; i++) {
869: newSupports[newqOff+offsets[q]++] = mesh->supports[off + i];
870: }
871: DMPlexGetTreeParent(dm,q,&parent,NULL);
872: }
873: }
874: }
875: PetscSectionDestroy(&mesh->supportSection);
876: mesh->supportSection = newSupportSection;
877: PetscFree(mesh->supports);
878: mesh->supports = newSupports;
879: PetscFree(offsets);
881: return(0);
882: }
884: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM,PetscSection,PetscSection,Mat);
885: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM,PetscSection,PetscSection,Mat);
889: static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical, PetscBool exchangeSupports)
890: {
891: DM_Plex *mesh = (DM_Plex *)dm->data;
892: DM refTree;
893: PetscInt size;
899: PetscObjectReference((PetscObject)parentSection);
900: PetscSectionDestroy(&mesh->parentSection);
901: mesh->parentSection = parentSection;
902: PetscSectionGetStorageSize(parentSection,&size);
903: if (parents != mesh->parents) {
904: PetscFree(mesh->parents);
905: PetscMalloc1(size,&mesh->parents);
906: PetscMemcpy(mesh->parents, parents, size * sizeof(*parents));
907: }
908: if (childIDs != mesh->childIDs) {
909: PetscFree(mesh->childIDs);
910: PetscMalloc1(size,&mesh->childIDs);
911: PetscMemcpy(mesh->childIDs, childIDs, size * sizeof(*childIDs));
912: }
913: DMPlexGetReferenceTree(dm,&refTree);
914: if (refTree) {
915: DMLabel canonLabel;
917: DMPlexGetLabel(refTree,"canonical",&canonLabel);
918: if (canonLabel) {
919: PetscInt i;
921: for (i = 0; i < size; i++) {
922: PetscInt canon;
923: DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon);
924: if (canon >= 0) {
925: mesh->childIDs[i] = canon;
926: }
927: }
928: }
929: mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_FromReference;
930: }
931: else {
932: mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_Direct;
933: }
934: DMPlexTreeSymmetrize(dm);
935: if (computeCanonical) {
936: PetscInt d, dim;
938: /* add the canonical label */
939: DMGetDimension(dm,&dim);
940: DMPlexCreateLabel(dm,"canonical");
941: for (d = 0; d <= dim; d++) {
942: PetscInt p, dStart, dEnd, canon = -1, cNumChildren;
943: const PetscInt *cChildren;
945: DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);
946: for (p = dStart; p < dEnd; p++) {
947: DMPlexGetTreeChildren(dm,p,&cNumChildren,&cChildren);
948: if (cNumChildren) {
949: canon = p;
950: break;
951: }
952: }
953: if (canon == -1) continue;
954: for (p = dStart; p < dEnd; p++) {
955: PetscInt numChildren, i;
956: const PetscInt *children;
958: DMPlexGetTreeChildren(dm,p,&numChildren,&children);
959: if (numChildren) {
960: 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);
961: DMPlexSetLabelValue(dm,"canonical",p,canon);
962: for (i = 0; i < numChildren; i++) {
963: DMPlexSetLabelValue(dm,"canonical",children[i],cChildren[i]);
964: }
965: }
966: }
967: }
968: }
969: if (exchangeSupports) {
970: DMPlexTreeExchangeSupports(dm);
971: }
972: mesh->createanchors = DMPlexCreateAnchors_Tree;
973: /* reset anchors */
974: DMPlexSetAnchors(dm,NULL,NULL);
975: return(0);
976: }
980: /*@
981: DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points. This routine also creates
982: the point-to-point constraints determined by the tree: a point is constained to the points in the closure of its
983: tree root.
985: Collective on dm
987: Input Parameters:
988: + dm - the DMPlex object
989: . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
990: offset indexes the parent and childID list; the reference count of parentSection is incremented
991: . parents - a list of the point parents; copied, can be destroyed
992: - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
993: the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed
995: Level: intermediate
997: .seealso: DMPlexGetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren()
998: @*/
999: PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[])
1000: {
1004: DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_FALSE,PETSC_TRUE);
1005: return(0);
1006: }
1010: /*@
1011: DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points.
1012: Collective on dm
1014: Input Parameters:
1015: . dm - the DMPlex object
1017: Output Parameters:
1018: + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section
1019: offset indexes the parent and childID list
1020: . parents - a list of the point parents
1021: . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then
1022: the child corresponds to the point in the reference tree with index childID
1023: . childSection - the inverse of the parent section
1024: - children - a list of the point children
1026: Level: intermediate
1028: .seealso: DMPlexSetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren()
1029: @*/
1030: PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[])
1031: {
1032: DM_Plex *mesh = (DM_Plex *)dm->data;
1036: if (parentSection) *parentSection = mesh->parentSection;
1037: if (parents) *parents = mesh->parents;
1038: if (childIDs) *childIDs = mesh->childIDs;
1039: if (childSection) *childSection = mesh->childSection;
1040: if (children) *children = mesh->children;
1041: return(0);
1042: }
1046: /*@
1047: DMPlexGetTreeParent - get the parent of a point in the tree describing the point hierarchy (not the Sieve DAG)
1049: Input Parameters:
1050: + dm - the DMPlex object
1051: - point - the query point
1053: Output Parameters:
1054: + parent - if not NULL, set to the parent of the point, or the point itself if the point does not have a parent
1055: - childID - if not NULL, set to the child ID of the point with respect to its parent, or 0 if the point
1056: does not have a parent
1058: Level: intermediate
1060: .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeChildren()
1061: @*/
1062: PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID)
1063: {
1064: DM_Plex *mesh = (DM_Plex *)dm->data;
1065: PetscSection pSec;
1070: pSec = mesh->parentSection;
1071: if (pSec && point >= pSec->pStart && point < pSec->pEnd) {
1072: PetscInt dof;
1074: PetscSectionGetDof (pSec, point, &dof);
1075: if (dof) {
1076: PetscInt off;
1078: PetscSectionGetOffset (pSec, point, &off);
1079: if (parent) *parent = mesh->parents[off];
1080: if (childID) *childID = mesh->childIDs[off];
1081: return(0);
1082: }
1083: }
1084: if (parent) {
1085: *parent = point;
1086: }
1087: if (childID) {
1088: *childID = 0;
1089: }
1090: return(0);
1091: }
1095: /*@C
1096: DMPlexGetTreeChildren - get the children of a point in the tree describing the point hierarchy (not the Sieve DAG)
1098: Input Parameters:
1099: + dm - the DMPlex object
1100: - point - the query point
1102: Output Parameters:
1103: + numChildren - if not NULL, set to the number of children
1104: - children - if not NULL, set to a list children, or set to NULL if the point has no children
1106: Level: intermediate
1108: Fortran Notes:
1109: Since it returns an array, this routine is only available in Fortran 90, and you must
1110: include petsc.h90 in your code.
1112: .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeParent()
1113: @*/
1114: PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[])
1115: {
1116: DM_Plex *mesh = (DM_Plex *)dm->data;
1117: PetscSection childSec;
1118: PetscInt dof = 0;
1123: childSec = mesh->childSection;
1124: if (childSec && point >= childSec->pStart && point < childSec->pEnd) {
1125: PetscSectionGetDof (childSec, point, &dof);
1126: }
1127: if (numChildren) *numChildren = dof;
1128: if (children) {
1129: if (dof) {
1130: PetscInt off;
1132: PetscSectionGetOffset (childSec, point, &off);
1133: *children = &mesh->children[off];
1134: }
1135: else {
1136: *children = NULL;
1137: }
1138: }
1139: return(0);
1140: }
1144: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM dm, PetscSection section, PetscSection cSec, Mat cMat)
1145: {
1146: PetscDS ds;
1147: PetscInt spdim;
1148: PetscInt numFields, f, c, cStart, cEnd, pStart, pEnd, conStart, conEnd;
1149: const PetscInt *anchors;
1150: PetscSection aSec;
1151: PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJparent, detJ, detJparent;
1152: IS aIS;
1156: DMPlexGetChart(dm,&pStart,&pEnd);
1157: DMGetDS(dm,&ds);
1158: PetscDSGetNumFields(ds,&numFields);
1159: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
1160: DMPlexGetAnchors(dm,&aSec,&aIS);
1161: ISGetIndices(aIS,&anchors);
1162: PetscSectionGetChart(cSec,&conStart,&conEnd);
1163: DMGetDimension(dm,&spdim);
1164: PetscMalloc6(spdim,&v0,spdim,&v0parent,spdim,&vtmp,spdim*spdim,&J,spdim*spdim,&Jparent,spdim*spdim,&invJparent);
1166: for (f = 0; f < numFields; f++) {
1167: PetscFE fe;
1168: PetscDualSpace space;
1169: PetscInt i, j, k, nPoints, offset;
1170: PetscInt fSize, fComp;
1171: PetscReal *B = NULL;
1172: PetscReal *weights, *pointsRef, *pointsReal;
1173: Mat Amat, Bmat, Xmat;
1175: PetscDSGetDiscretization(ds,f,(PetscObject *)(&fe));
1176: PetscFEGetDualSpace(fe,&space);
1177: PetscDualSpaceGetDimension(space,&fSize);
1178: PetscFEGetNumComponents(fe,&fComp);
1179: MatCreate(PETSC_COMM_SELF,&Amat);
1180: MatSetSizes(Amat,fSize,fSize,fSize,fSize);
1181: MatSetType(Amat,MATSEQDENSE);
1182: MatSetUp(Amat);
1183: MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Bmat);
1184: MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Xmat);
1185: nPoints = 0;
1186: for (i = 0; i < fSize; i++) {
1187: PetscInt qPoints;
1188: PetscQuadrature quad;
1190: PetscDualSpaceGetFunctional(space,i,&quad);
1191: PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);
1192: nPoints += qPoints;
1193: }
1194: PetscMalloc3(nPoints,&weights,spdim*nPoints,&pointsRef,spdim*nPoints,&pointsReal);
1195: offset = 0;
1196: for (i = 0; i < fSize; i++) {
1197: PetscInt qPoints;
1198: const PetscReal *p, *w;
1199: PetscQuadrature quad;
1201: PetscDualSpaceGetFunctional(space,i,&quad);
1202: PetscQuadratureGetData(quad,NULL,&qPoints,&p,&w);
1203: PetscMemcpy(weights+offset,w,qPoints*sizeof(*w));
1204: PetscMemcpy(pointsRef+spdim*offset,p,spdim*qPoints*sizeof(*p));
1205: offset += qPoints;
1206: }
1207: PetscFEGetTabulation(fe,nPoints,pointsRef,&B,NULL,NULL);
1208: offset = 0;
1209: for (i = 0; i < fSize; i++) {
1210: PetscInt qPoints;
1211: PetscQuadrature quad;
1213: PetscDualSpaceGetFunctional(space,i,&quad);
1214: PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);
1215: for (j = 0; j < fSize; j++) {
1216: PetscScalar val = 0.;
1218: for (k = 0; k < qPoints; k++) {
1219: val += B[((offset + k) * fSize + j) * fComp] * weights[k];
1220: }
1221: MatSetValue(Amat,i,j,val,INSERT_VALUES);
1222: }
1223: offset += qPoints;
1224: }
1225: PetscFERestoreTabulation(fe,nPoints,pointsRef,&B,NULL,NULL);
1226: MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
1227: MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);
1228: MatLUFactor(Amat,NULL,NULL,NULL);
1229: for (c = cStart; c < cEnd; c++) {
1230: PetscInt parent;
1231: PetscInt closureSize, closureSizeP, *closure = NULL, *closureP = NULL;
1232: PetscInt *childOffsets, *parentOffsets;
1234: DMPlexGetTreeParent(dm,c,&parent,NULL);
1235: if (parent == c) continue;
1236: DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1237: for (i = 0; i < closureSize; i++) {
1238: PetscInt p = closure[2*i];
1239: PetscInt conDof;
1241: if (p < conStart || p >= conEnd) continue;
1242: if (numFields > 1) {
1243: PetscSectionGetFieldDof(cSec,p,f,&conDof);
1244: }
1245: else {
1246: PetscSectionGetDof(cSec,p,&conDof);
1247: }
1248: if (conDof) break;
1249: }
1250: if (i == closureSize) {
1251: DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1252: continue;
1253: }
1255: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, NULL, &detJ);
1256: DMPlexComputeCellGeometryFEM(dm, parent, NULL, v0parent, Jparent, invJparent, &detJparent);
1257: for (i = 0; i < nPoints; i++) {
1258: CoordinatesRefToReal(spdim, spdim, v0, J, &pointsRef[i*spdim],vtmp);
1259: CoordinatesRealToRef(spdim, spdim, v0parent, invJparent, vtmp, &pointsReal[i*spdim]);
1260: }
1261: PetscFEGetTabulation(fe,nPoints,pointsReal,&B,NULL,NULL);
1262: offset = 0;
1263: for (i = 0; i < fSize; i++) {
1264: PetscInt qPoints;
1265: PetscQuadrature quad;
1267: PetscDualSpaceGetFunctional(space,i,&quad);
1268: PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);
1269: for (j = 0; j < fSize; j++) {
1270: PetscScalar val = 0.;
1272: for (k = 0; k < qPoints; k++) {
1273: val += B[((offset + k) * fSize + j) * fComp] * weights[k];
1274: }
1275: MatSetValue(Bmat,i,j,val,INSERT_VALUES);
1276: }
1277: offset += qPoints;
1278: }
1279: PetscFERestoreTabulation(fe,nPoints,pointsReal,&B,NULL,NULL);
1280: MatAssemblyBegin(Bmat,MAT_FINAL_ASSEMBLY);
1281: MatAssemblyEnd(Bmat,MAT_FINAL_ASSEMBLY);
1282: MatMatSolve(Amat,Bmat,Xmat);
1283: DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);
1284: PetscMalloc2(closureSize+1,&childOffsets,closureSizeP+1,&parentOffsets);
1285: childOffsets[0] = 0;
1286: for (i = 0; i < closureSize; i++) {
1287: PetscInt p = closure[2*i];
1288: PetscInt dof;
1290: if (numFields > 1) {
1291: PetscSectionGetFieldDof(section,p,f,&dof);
1292: }
1293: else {
1294: PetscSectionGetDof(section,p,&dof);
1295: }
1296: childOffsets[i+1]=childOffsets[i]+dof / fComp;
1297: }
1298: parentOffsets[0] = 0;
1299: for (i = 0; i < closureSizeP; i++) {
1300: PetscInt p = closureP[2*i];
1301: PetscInt dof;
1303: if (numFields > 1) {
1304: PetscSectionGetFieldDof(section,p,f,&dof);
1305: }
1306: else {
1307: PetscSectionGetDof(section,p,&dof);
1308: }
1309: parentOffsets[i+1]=parentOffsets[i]+dof / fComp;
1310: }
1311: for (i = 0; i < closureSize; i++) {
1312: PetscInt conDof, conOff, aDof, aOff;
1313: PetscInt p = closure[2*i];
1314: PetscInt o = closure[2*i+1];
1316: if (p < conStart || p >= conEnd) continue;
1317: if (numFields > 1) {
1318: PetscSectionGetFieldDof(cSec,p,f,&conDof);
1319: PetscSectionGetFieldOffset(cSec,p,f,&conOff);
1320: }
1321: else {
1322: PetscSectionGetDof(cSec,p,&conDof);
1323: PetscSectionGetOffset(cSec,p,&conOff);
1324: }
1325: if (!conDof) continue;
1326: PetscSectionGetDof(aSec,p,&aDof);
1327: PetscSectionGetOffset(aSec,p,&aOff);
1328: for (k = 0; k < aDof; k++) {
1329: PetscInt a = anchors[aOff + k];
1330: PetscInt aSecDof, aSecOff;
1332: if (numFields > 1) {
1333: PetscSectionGetFieldDof(section,a,f,&aSecDof);
1334: PetscSectionGetFieldOffset(section,a,f,&aSecOff);
1335: }
1336: else {
1337: PetscSectionGetDof(section,a,&aSecDof);
1338: PetscSectionGetOffset(section,a,&aSecOff);
1339: }
1340: if (!aSecDof) continue;
1342: for (j = 0; j < closureSizeP; j++) {
1343: PetscInt q = closureP[2*j];
1344: PetscInt oq = closureP[2*j+1];
1346: if (q == a) {
1347: PetscInt r, s, t;
1349: for (r = childOffsets[i]; r < childOffsets[i+1]; r++) {
1350: for (s = parentOffsets[j]; s < parentOffsets[j+1]; s++) {
1351: PetscScalar val;
1352: PetscInt insertCol, insertRow;
1354: MatGetValue(Xmat,r,s,&val);
1355: if (o >= 0) {
1356: insertRow = conOff + fComp * (r - childOffsets[i]);
1357: }
1358: else {
1359: insertRow = conOff + fComp * (childOffsets[i + 1] - 1 - r);
1360: }
1361: if (oq >= 0) {
1362: insertCol = aSecOff + fComp * (s - parentOffsets[j]);
1363: }
1364: else {
1365: insertCol = aSecOff + fComp * (parentOffsets[j + 1] - 1 - s);
1366: }
1367: for (t = 0; t < fComp; t++) {
1368: MatSetValue(cMat,insertRow + t,insertCol + t,val,INSERT_VALUES);
1369: }
1370: }
1371: }
1372: }
1373: }
1374: }
1375: }
1376: PetscFree2(childOffsets,parentOffsets);
1377: DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);
1378: DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);
1379: }
1380: MatDestroy(&Amat);
1381: MatDestroy(&Bmat);
1382: MatDestroy(&Xmat);
1383: PetscFree3(weights,pointsRef,pointsReal);
1384: }
1385: MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
1386: MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
1387: PetscFree6(v0,v0parent,vtmp,J,Jparent,invJparent);
1388: ISRestoreIndices(aIS,&anchors);
1390: return(0);
1391: }
1395: static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm, PetscSection section, PetscSection conSec, Mat cMat)
1396: {
1397: DM refTree;
1398: PetscDS ds;
1399: Mat refCmat;
1400: PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN;
1401: PetscScalar ***refPointFieldMats, *pointWork;
1402: PetscSection refConSec, refAnSec, anSec, refSection;
1403: IS refAnIS, anIS;
1404: const PetscInt *refAnchors, *anchors;
1409: DMGetDS(dm,&ds);
1410: PetscDSGetNumFields(ds,&numFields);
1411: DMPlexGetReferenceTree(dm,&refTree);
1412: DMSetDS(refTree,ds);
1413: DMGetDefaultConstraints(refTree,&refConSec,&refCmat);
1414: DMPlexGetChart(refTree,&pRefStart,&pRefEnd);
1415: DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);
1416: DMPlexGetAnchors(dm,&anSec,&anIS);
1417: ISGetIndices(refAnIS,&refAnchors);
1418: ISGetIndices(anIS,&anchors);
1419: DMGetDefaultSection(refTree,&refSection);
1420: PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);
1421: PetscSectionGetChart(conSec,&conStart,&conEnd);
1422: PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);
1423: PetscMalloc1(pRefEnd-pRefStart,&refPointFieldN);
1424: PetscSectionGetMaxDof(refConSec,&maxDof);
1425: PetscSectionGetMaxDof(refAnSec,&maxAnDof);
1426: PetscMalloc1(maxDof,&rows);
1427: PetscMalloc1(maxDof*maxAnDof,&cols);
1428: PetscMalloc1(maxDof*maxDof*maxAnDof,&pointWork);
1430: /* step 1: get submats for every constrained point in the reference tree */
1431: for (p = pRefStart; p < pRefEnd; p++) {
1432: PetscInt parent, closureSize, *closure = NULL, pDof;
1434: DMPlexGetTreeParent(refTree,p,&parent,NULL);
1435: PetscSectionGetDof(refConSec,p,&pDof);
1436: if (!pDof || parent == p) continue;
1438: PetscMalloc1(numFields,&refPointFieldMats[p-pRefStart]);
1439: PetscCalloc1(numFields,&refPointFieldN[p-pRefStart]);
1440: DMPlexGetTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);
1441: for (f = 0; f < numFields; f++) {
1442: PetscInt cDof, cOff, numCols, r, i, fComp;
1443: PetscFE fe;
1445: PetscDSGetDiscretization(ds,f,(PetscObject *) &fe);
1446: PetscFEGetNumComponents(fe,&fComp);
1448: if (numFields > 1) {
1449: PetscSectionGetFieldDof(refConSec,p,f,&cDof);
1450: PetscSectionGetFieldOffset(refConSec,p,f,&cOff);
1451: }
1452: else {
1453: PetscSectionGetDof(refConSec,p,&cDof);
1454: PetscSectionGetOffset(refConSec,p,&cOff);
1455: }
1457: if (!cDof) continue;
1458: for (r = 0; r < cDof; r++) {
1459: rows[r] = cOff + r;
1460: }
1461: numCols = 0;
1462: for (i = 0; i < closureSize; i++) {
1463: PetscInt q = closure[2*i];
1464: PetscInt o = closure[2*i+1];
1465: PetscInt aDof, aOff, j;
1467: if (numFields > 1) {
1468: PetscSectionGetFieldDof(refSection,q,f,&aDof);
1469: PetscSectionGetFieldOffset(refSection,q,f,&aOff);
1470: }
1471: else {
1472: PetscSectionGetDof(refSection,q,&aDof);
1473: PetscSectionGetOffset(refSection,q,&aOff);
1474: }
1476: for (j = 0; j < aDof; j++) {
1477: PetscInt node = (o >= 0) ? (j / fComp) : ((aDof - 1 - j) / fComp);
1478: PetscInt comp = (j % fComp);
1480: cols[numCols++] = aOff + node * fComp + comp;
1481: }
1482: }
1483: refPointFieldN[p-pRefStart][f] = numCols;
1484: PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);
1485: MatGetValues(refCmat,cDof,rows,numCols,cols,refPointFieldMats[p-pRefStart][f]);
1486: }
1487: DMPlexRestoreTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);
1488: }
1490: /* step 2: compute the preorder */
1491: DMPlexGetChart(dm,&pStart,&pEnd);
1492: PetscMalloc2(pEnd-pStart,&perm,pEnd-pStart,&iperm);
1493: for (p = pStart; p < pEnd; p++) {
1494: perm[p - pStart] = p;
1495: iperm[p - pStart] = p-pStart;
1496: }
1497: for (p = 0; p < pEnd - pStart;) {
1498: PetscInt point = perm[p];
1499: PetscInt parent;
1501: DMPlexGetTreeParent(dm,point,&parent,NULL);
1502: if (parent == point) {
1503: p++;
1504: }
1505: else {
1506: PetscInt size, closureSize, *closure = NULL, i;
1508: DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1509: for (i = 0; i < closureSize; i++) {
1510: PetscInt q = closure[2*i];
1511: if (iperm[q-pStart] > iperm[point-pStart]) {
1512: /* swap */
1513: perm[p] = q;
1514: perm[iperm[q-pStart]] = point;
1515: iperm[point-pStart] = iperm[q-pStart];
1516: iperm[q-pStart] = p;
1517: break;
1518: }
1519: }
1520: size = closureSize;
1521: DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1522: if (i == size) {
1523: p++;
1524: }
1525: }
1526: }
1528: /* step 3: fill the constraint matrix */
1529: /* we are going to use a preorder progressive fill strategy. Mat doesn't
1530: * allow progressive fill without assembly, so we are going to set up the
1531: * values outside of the Mat first.
1532: */
1533: {
1534: PetscInt nRows, row, nnz;
1535: PetscBool done;
1536: const PetscInt *ia, *ja;
1537: PetscScalar *vals;
1539: MatGetRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);
1540: if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not get RowIJ of constraint matrix");
1541: nnz = ia[nRows];
1542: /* malloc and then zero rows right before we fill them: this way valgrind
1543: * can tell if we are doing progressive fill in the wrong order */
1544: PetscMalloc1(nnz,&vals);
1545: for (p = 0; p < pEnd - pStart; p++) {
1546: PetscInt parent, childid, closureSize, *closure = NULL;
1547: PetscInt point = perm[p], pointDof;
1549: DMPlexGetTreeParent(dm,point,&parent,&childid);
1550: if ((point < conStart) || (point >= conEnd) || (parent == point)) continue;
1551: PetscSectionGetDof(conSec,point,&pointDof);
1552: if (!pointDof) continue;
1553: DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1554: for (f = 0; f < numFields; f++) {
1555: PetscInt cDof, cOff, numCols, numFillCols, i, r, fComp, matOffset, offset;
1556: PetscScalar *pointMat;
1557: PetscFE fe;
1559: PetscDSGetDiscretization(ds,f,(PetscObject *) &fe);
1560: PetscFEGetNumComponents(fe,&fComp);
1562: if (numFields > 1) {
1563: PetscSectionGetFieldDof(conSec,point,f,&cDof);
1564: PetscSectionGetFieldOffset(conSec,point,f,&cOff);
1565: }
1566: else {
1567: PetscSectionGetDof(conSec,point,&cDof);
1568: PetscSectionGetOffset(conSec,point,&cOff);
1569: }
1570: if (!cDof) continue;
1572: /* make sure that every row for this point is the same size */
1573: #if defined(PETSC_USE_DEBUG)
1574: for (r = 0; r < cDof; r++) {
1575: if (cDof > 1 && r) {
1576: if ((ia[cOff+r+1]-ia[cOff+r]) != (ia[cOff+r]-ia[cOff+r-1])) {
1577: SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Two point rows have different nnz: %d vs. %d", (ia[rows[r]+1]-ia[rows[r]]), (ia[rows[r]]-ia[rows[r]-1]));
1578: }
1579: }
1580: }
1581: #endif
1582: /* zero rows */
1583: for (i = ia[cOff] ; i< ia[cOff+cDof];i++) {
1584: vals[i] = 0.;
1585: }
1586: matOffset = ia[cOff];
1587: numFillCols = ia[cOff+1] - matOffset;
1588: pointMat = refPointFieldMats[childid-pRefStart][f];
1589: numCols = refPointFieldN[childid-pRefStart][f];
1590: offset = 0;
1591: for (i = 0; i < closureSize; i++) {
1592: PetscInt q = closure[2*i];
1593: PetscInt o = closure[2*i+1];
1594: PetscInt aDof, aOff, j, k, qConDof, qConOff;
1596: qConDof = qConOff = 0;
1597: if (numFields > 1) {
1598: PetscSectionGetFieldDof(section,q,f,&aDof);
1599: PetscSectionGetFieldOffset(section,q,f,&aOff);
1600: if (q >= conStart && q < conEnd) {
1601: PetscSectionGetFieldDof(conSec,q,f,&qConDof);
1602: PetscSectionGetFieldOffset(conSec,q,f,&qConOff);
1603: }
1604: }
1605: else {
1606: PetscSectionGetDof(section,q,&aDof);
1607: PetscSectionGetOffset(section,q,&aOff);
1608: if (q >= conStart && q < conEnd) {
1609: PetscSectionGetDof(conSec,q,&qConDof);
1610: PetscSectionGetOffset(conSec,q,&qConOff);
1611: }
1612: }
1613: if (!aDof) continue;
1614: if (qConDof) {
1615: /* this point has anchors: its rows of the matrix should already
1616: * be filled, thanks to preordering */
1617: /* first multiply into pointWork, then set in matrix */
1618: PetscInt aMatOffset = ia[qConOff];
1619: PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset;
1620: for (r = 0; r < cDof; r++) {
1621: for (j = 0; j < aNumFillCols; j++) {
1622: PetscScalar inVal = 0;
1623: for (k = 0; k < aDof; k++) {
1624: PetscInt node = (o >= 0) ? (k / fComp) : ((aDof - 1 - k) / fComp);
1625: PetscInt comp = (k % fComp);
1626: PetscInt col = node * fComp + comp;
1628: inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j];
1629: }
1630: pointWork[r * aNumFillCols + j] = inVal;
1631: }
1632: }
1633: /* assume that the columns are sorted, spend less time searching */
1634: for (j = 0, k = 0; j < aNumFillCols; j++) {
1635: PetscInt col = ja[aMatOffset + j];
1636: for (;k < numFillCols; k++) {
1637: if (ja[matOffset + k] == col) {
1638: break;
1639: }
1640: }
1641: if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, col);
1642: for (r = 0; r < cDof; r++) {
1643: vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j];
1644: }
1645: }
1646: }
1647: else {
1648: /* find where to put this portion of pointMat into the matrix */
1649: for (k = 0; k < numFillCols; k++) {
1650: if (ja[matOffset + k] == aOff) {
1651: break;
1652: }
1653: }
1654: if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, aOff);
1655: for (j = 0; j < aDof; j++) {
1656: PetscInt node = (o >= 0) ? (j / fComp) : ((aDof - 1 - j) / fComp);
1657: PetscInt comp = (j % fComp);
1658: PetscInt col = node * fComp + comp;
1659: for (r = 0; r < cDof; r++) {
1660: vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + col];
1661: }
1662: }
1663: }
1664: offset += aDof;
1665: }
1666: }
1667: DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
1668: }
1669: for (row = 0; row < nRows; row++) {
1670: MatSetValues(cMat,1,&row,ia[row+1]-ia[row],&ja[ia[row]],&vals[ia[row]],INSERT_VALUES);
1671: }
1672: MatRestoreRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);
1673: if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not restore RowIJ of constraint matrix");
1674: MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
1675: MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
1676: PetscFree(vals);
1677: }
1679: /* clean up */
1680: ISRestoreIndices(refAnIS,&refAnchors);
1681: ISRestoreIndices(anIS,&anchors);
1682: PetscFree2(perm,iperm);
1683: PetscFree(rows);
1684: PetscFree(cols);
1685: PetscFree(pointWork);
1686: for (p = pRefStart; p < pRefEnd; p++) {
1687: PetscInt parent, pDof;
1689: DMPlexGetTreeParent(refTree,p,&parent,NULL);
1690: PetscSectionGetDof(refConSec,p,&pDof);
1691: if (!pDof || parent == p) continue;
1693: for (f = 0; f < numFields; f++) {
1694: PetscInt cDof;
1696: if (numFields > 1) {
1697: PetscSectionGetFieldDof(refConSec,p,f,&cDof);
1698: }
1699: else {
1700: PetscSectionGetDof(refConSec,p,&cDof);
1701: }
1703: if (!cDof) continue;
1704: PetscFree(refPointFieldMats[p - pRefStart][f]);
1705: }
1706: PetscFree(refPointFieldMats[p - pRefStart]);
1707: PetscFree(refPointFieldN[p - pRefStart]);
1708: }
1709: PetscFree(refPointFieldMats);
1710: PetscFree(refPointFieldN);
1711: return(0);
1712: }
1716: /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of
1717: * a non-conforming mesh. Local refinement comes later */
1718: PetscErrorCode DMPlexTreeRefineCell (DM dm, PetscInt cell, DM *ncdm)
1719: {
1720: DM K;
1721: PetscMPIInt rank;
1722: PetscInt dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd;
1723: PetscInt numNewCones, *newConeSizes, *newCones, *newOrientations;
1724: PetscInt *Kembedding;
1725: PetscInt *cellClosure=NULL, nc;
1726: PetscScalar *newVertexCoords;
1727: PetscInt numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset;
1728: PetscSection parentSection;
1732: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1733: DMGetDimension(dm,&dim);
1734: DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm);
1735: DMSetDimension(*ncdm,dim);
1737: DMPlexGetChart(dm, &pStart, &pEnd);
1738: PetscSectionCreate(PetscObjectComm((PetscObject)dm),&parentSection);
1739: DMPlexGetReferenceTree(dm,&K);
1740: if (!rank) {
1741: /* compute the new charts */
1742: PetscMalloc5(dim+1,&pNewCount,dim+1,&pNewStart,dim+1,&pNewEnd,dim+1,&pOldStart,dim+1,&pOldEnd);
1743: offset = 0;
1744: for (d = 0; d <= dim; d++) {
1745: PetscInt pOldCount, kStart, kEnd, k;
1747: pNewStart[d] = offset;
1748: DMPlexGetHeightStratum(dm,d,&pOldStart[d],&pOldEnd[d]);
1749: DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
1750: pOldCount = pOldEnd[d] - pOldStart[d];
1751: /* adding the new points */
1752: pNewCount[d] = pOldCount + kEnd - kStart;
1753: if (!d) {
1754: /* removing the cell */
1755: pNewCount[d]--;
1756: }
1757: for (k = kStart; k < kEnd; k++) {
1758: PetscInt parent;
1759: DMPlexGetTreeParent(K,k,&parent,NULL);
1760: if (parent == k) {
1761: /* avoid double counting points that won't actually be new */
1762: pNewCount[d]--;
1763: }
1764: }
1765: pNewEnd[d] = pNewStart[d] + pNewCount[d];
1766: offset = pNewEnd[d];
1768: }
1769: 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]);
1770: /* get the current closure of the cell that we are removing */
1771: DMPlexGetTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);
1773: PetscMalloc1(pNewEnd[dim],&newConeSizes);
1774: {
1775: PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j;
1777: DMPlexGetChart(K,&kStart,&kEnd);
1778: PetscMalloc4(kEnd-kStart,&Kembedding,kEnd-kStart,&perm,kEnd-kStart,&iperm,kEnd-kStart,&preOrient);
1780: for (k = kStart; k < kEnd; k++) {
1781: perm[k - kStart] = k;
1782: iperm [k - kStart] = k - kStart;
1783: preOrient[k - kStart] = 0;
1784: }
1786: DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);
1787: for (j = 1; j < closureSizeK; j++) {
1788: PetscInt parentOrientA = closureK[2*j+1];
1789: PetscInt parentOrientB = cellClosure[2*j+1];
1790: PetscInt p, q;
1792: p = closureK[2*j];
1793: q = cellClosure[2*j];
1794: for (d = 0; d <= dim; d++) {
1795: if (q >= pOldStart[d] && q < pOldEnd[d]) {
1796: Kembedding[p] = (q - pOldStart[d]) + pNewStart[d];
1797: }
1798: }
1799: if (parentOrientA != parentOrientB) {
1800: PetscInt numChildren, i;
1801: const PetscInt *children;
1803: DMPlexGetTreeChildren(K,p,&numChildren,&children);
1804: for (i = 0; i < numChildren; i++) {
1805: PetscInt kPerm, oPerm;
1807: k = children[i];
1808: DMPlexReferenceTreeGetChildSymmetry(K,p,parentOrientA,0,k,parentOrientB,&oPerm,&kPerm);
1809: /* perm = what refTree position I'm in */
1810: perm[kPerm-kStart] = k;
1811: /* iperm = who is at this position */
1812: iperm[k-kStart] = kPerm-kStart;
1813: preOrient[kPerm-kStart] = oPerm;
1814: }
1815: }
1816: }
1817: DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);
1818: }
1819: PetscSectionSetChart(parentSection,0,pNewEnd[dim]);
1820: offset = 0;
1821: numNewCones = 0;
1822: for (d = 0; d <= dim; d++) {
1823: PetscInt kStart, kEnd, k;
1824: PetscInt p;
1825: PetscInt size;
1827: for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1828: /* skip cell 0 */
1829: if (p == cell) continue;
1830: /* old cones to new cones */
1831: DMPlexGetConeSize(dm,p,&size);
1832: newConeSizes[offset++] = size;
1833: numNewCones += size;
1834: }
1836: DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
1837: for (k = kStart; k < kEnd; k++) {
1838: PetscInt kParent;
1840: DMPlexGetTreeParent(K,k,&kParent,NULL);
1841: if (kParent != k) {
1842: Kembedding[k] = offset;
1843: DMPlexGetConeSize(K,k,&size);
1844: newConeSizes[offset++] = size;
1845: numNewCones += size;
1846: if (kParent != 0) {
1847: PetscSectionSetDof(parentSection,Kembedding[k],1);
1848: }
1849: }
1850: }
1851: }
1853: PetscSectionSetUp(parentSection);
1854: PetscSectionGetStorageSize(parentSection,&numPointsWithParents);
1855: PetscMalloc2(numNewCones,&newCones,numNewCones,&newOrientations);
1856: PetscMalloc2(numPointsWithParents,&parents,numPointsWithParents,&childIDs);
1858: /* fill new cones */
1859: offset = 0;
1860: for (d = 0; d <= dim; d++) {
1861: PetscInt kStart, kEnd, k, l;
1862: PetscInt p;
1863: PetscInt size;
1864: const PetscInt *cone, *orientation;
1866: for (p = pOldStart[d]; p < pOldEnd[d]; p++) {
1867: /* skip cell 0 */
1868: if (p == cell) continue;
1869: /* old cones to new cones */
1870: DMPlexGetConeSize(dm,p,&size);
1871: DMPlexGetCone(dm,p,&cone);
1872: DMPlexGetConeOrientation(dm,p,&orientation);
1873: for (l = 0; l < size; l++) {
1874: newCones[offset] = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1];
1875: newOrientations[offset++] = orientation[l];
1876: }
1877: }
1879: DMPlexGetHeightStratum(K,d,&kStart,&kEnd);
1880: for (k = kStart; k < kEnd; k++) {
1881: PetscInt kPerm = perm[k], kParent;
1882: PetscInt preO = preOrient[k];
1884: DMPlexGetTreeParent(K,k,&kParent,NULL);
1885: if (kParent != k) {
1886: /* embed new cones */
1887: DMPlexGetConeSize(K,k,&size);
1888: DMPlexGetCone(K,kPerm,&cone);
1889: DMPlexGetConeOrientation(K,kPerm,&orientation);
1890: for (l = 0; l < size; l++) {
1891: PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size -(preO + 1) - l) % size);
1892: PetscInt newO, lSize, oTrue;
1894: q = iperm[cone[m]];
1895: newCones[offset] = Kembedding[q];
1896: DMPlexGetConeSize(K,q,&lSize);
1897: oTrue = orientation[m];
1898: oTrue = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2);
1899: newO = DihedralCompose(lSize,oTrue,preOrient[q]);
1900: newOrientations[offset++] = newO;
1901: }
1902: if (kParent != 0) {
1903: PetscInt newPoint = Kembedding[kParent];
1904: PetscSectionGetOffset(parentSection,Kembedding[k],&pOffset);
1905: parents[pOffset] = newPoint;
1906: childIDs[pOffset] = k;
1907: }
1908: }
1909: }
1910: }
1912: PetscMalloc1(dim*(pNewEnd[dim]-pNewStart[dim]),&newVertexCoords);
1914: /* fill coordinates */
1915: offset = 0;
1916: {
1917: PetscInt kStart, kEnd, l;
1918: PetscSection vSection;
1919: PetscInt v;
1920: Vec coords;
1921: PetscScalar *coordvals;
1922: PetscInt dof, off;
1923: PetscReal v0[3], J[9], detJ;
1925: #if defined(PETSC_USE_DEBUG)
1926: {
1927: PetscInt k;
1928: DMPlexGetHeightStratum(K,0,&kStart,&kEnd);
1929: for (k = kStart; k < kEnd; k++) {
1930: DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ);
1931: if (detJ <= 0.) SETERRQ1 (PETSC_COMM_SELF,PETSC_ERR_PLIB,"reference tree cell %d has bad determinant",k);
1932: }
1933: }
1934: #endif
1935: DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ);
1936: DMGetCoordinateSection(dm,&vSection);
1937: DMGetCoordinatesLocal(dm,&coords);
1938: VecGetArray(coords,&coordvals);
1939: for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) {
1941: PetscSectionGetDof(vSection,v,&dof);
1942: PetscSectionGetOffset(vSection,v,&off);
1943: for (l = 0; l < dof; l++) {
1944: newVertexCoords[offset++] = coordvals[off + l];
1945: }
1946: }
1947: VecRestoreArray(coords,&coordvals);
1949: DMGetCoordinateSection(K,&vSection);
1950: DMGetCoordinatesLocal(K,&coords);
1951: VecGetArray(coords,&coordvals);
1952: DMPlexGetDepthStratum(K,0,&kStart,&kEnd);
1953: for (v = kStart; v < kEnd; v++) {
1954: PetscReal coord[3], newCoord[3];
1955: PetscInt vPerm = perm[v];
1956: PetscInt kParent;
1958: DMPlexGetTreeParent(K,v,&kParent,NULL);
1959: if (kParent != v) {
1960: /* this is a new vertex */
1961: PetscSectionGetOffset(vSection,vPerm,&off);
1962: for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off+l]);
1963: CoordinatesRefToReal(dim, dim, v0, J, coord, newCoord);
1964: for (l = 0; l < dim; ++l) newVertexCoords[offset+l] = newCoord[l];
1965: offset += dim;
1966: }
1967: }
1968: VecRestoreArray(coords,&coordvals);
1969: }
1971: /* need to reverse the order of pNewCount: vertices first, cells last */
1972: for (d = 0; d < (dim + 1) / 2; d++) {
1973: PetscInt tmp;
1975: tmp = pNewCount[d];
1976: pNewCount[d] = pNewCount[dim - d];
1977: pNewCount[dim - d] = tmp;
1978: }
1980: DMPlexCreateFromDAG(*ncdm,dim,pNewCount,newConeSizes,newCones,newOrientations,newVertexCoords);
1981: DMPlexSetReferenceTree(*ncdm,K);
1982: DMPlexSetTree(*ncdm,parentSection,parents,childIDs);
1984: /* clean up */
1985: DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);
1986: PetscFree5(pNewCount,pNewStart,pNewEnd,pOldStart,pOldEnd);
1987: PetscFree(newConeSizes);
1988: PetscFree2(newCones,newOrientations);
1989: PetscFree(newVertexCoords);
1990: PetscFree2(parents,childIDs);
1991: PetscFree4(Kembedding,perm,iperm,preOrient);
1992: }
1993: else {
1994: PetscInt p, counts[4];
1995: PetscInt *coneSizes, *cones, *orientations;
1996: Vec coordVec;
1997: PetscScalar *coords;
1999: for (d = 0; d <= dim; d++) {
2000: PetscInt dStart, dEnd;
2002: DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);
2003: counts[d] = dEnd - dStart;
2004: }
2005: PetscMalloc1(pEnd-pStart,&coneSizes);
2006: for (p = pStart; p < pEnd; p++) {
2007: DMPlexGetConeSize(dm,p,&coneSizes[p-pStart]);
2008: }
2009: DMPlexGetCones(dm, &cones);
2010: DMPlexGetConeOrientations(dm, &orientations);
2011: DMGetCoordinatesLocal(dm,&coordVec);
2012: VecGetArray(coordVec,&coords);
2014: PetscSectionSetChart(parentSection,pStart,pEnd);
2015: PetscSectionSetUp(parentSection);
2016: DMPlexCreateFromDAG(*ncdm,dim,counts,coneSizes,cones,orientations,NULL);
2017: DMPlexSetReferenceTree(*ncdm,K);
2018: DMPlexSetTree(*ncdm,parentSection,NULL,NULL);
2019: VecRestoreArray(coordVec,&coords);
2020: }
2021: PetscSectionDestroy(&parentSection);
2023: return(0);
2024: }