Actual source code: plextransform.c
1: #include <petsc/private/dmplextransformimpl.h>
3: #include <petsc/private/petscfeimpl.h>
5: PetscClassId DMPLEXTRANSFORM_CLASSID;
7: PetscFunctionList DMPlexTransformList = NULL;
8: PetscBool DMPlexTransformRegisterAllCalled = PETSC_FALSE;
10: /* Construct cell type order since we must loop over cell types in depth order */
11: PetscErrorCode DMPlexCreateCellTypeOrder_Internal(DMPolytopeType ctCell, PetscInt *ctOrder[], PetscInt *ctOrderInv[])
12: {
13: PetscInt *ctO, *ctOInv;
14: PetscInt c, d, off = 0;
18: PetscCalloc2(DM_NUM_POLYTOPES+1, &ctO, DM_NUM_POLYTOPES+1, &ctOInv);
19: for (d = 3; d >= DMPolytopeTypeGetDim(ctCell); --d) {
20: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
21: if (DMPolytopeTypeGetDim((DMPolytopeType) c) != d) continue;
22: ctO[off++] = c;
23: }
24: }
25: if (DMPolytopeTypeGetDim(ctCell) != 0) {
26: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
27: if (DMPolytopeTypeGetDim((DMPolytopeType) c) != 0) continue;
28: ctO[off++] = c;
29: }
30: }
31: for (d = DMPolytopeTypeGetDim(ctCell)-1; d > 0; --d) {
32: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
33: if (DMPolytopeTypeGetDim((DMPolytopeType) c) != d) continue;
34: ctO[off++] = c;
35: }
36: }
37: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
38: if (DMPolytopeTypeGetDim((DMPolytopeType) c) >= 0) continue;
39: ctO[off++] = c;
40: }
41: if (off != DM_NUM_POLYTOPES+1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid offset %D for cell type order", off);
42: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
43: ctOInv[ctO[c]] = c;
44: }
45: *ctOrder = ctO;
46: *ctOrderInv = ctOInv;
47: return(0);
48: }
50: /*@C
51: DMPlexTransformRegister - Adds a new transform component implementation
53: Not Collective
55: Input Parameters:
56: + name - The name of a new user-defined creation routine
57: - create_func - The creation routine itself
59: Notes:
60: DMPlexTransformRegister() may be called multiple times to add several user-defined transforms
62: Sample usage:
63: .vb
64: DMPlexTransformRegister("my_transform", MyTransformCreate);
65: .ve
67: Then, your transform type can be chosen with the procedural interface via
68: .vb
69: DMPlexTransformCreate(MPI_Comm, DMPlexTransform *);
70: DMPlexTransformSetType(DMPlexTransform, "my_transform");
71: .ve
72: or at runtime via the option
73: .vb
74: -dm_plex_transform_type my_transform
75: .ve
77: Level: advanced
79: .seealso: DMPlexTransformRegisterAll(), DMPlexTransformRegisterDestroy()
80: @*/
81: PetscErrorCode DMPlexTransformRegister(const char name[], PetscErrorCode (*create_func)(DMPlexTransform))
82: {
86: DMInitializePackage();
87: PetscFunctionListAdd(&DMPlexTransformList, name, create_func);
88: return(0);
89: }
91: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Filter(DMPlexTransform);
92: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform);
93: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToBox(DMPlexTransform);
94: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform);
95: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform);
96: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform);
98: /*@C
99: DMPlexTransformRegisterAll - Registers all of the transform components in the DM package.
101: Not Collective
103: Level: advanced
105: .seealso: DMRegisterAll(), DMPlexTransformRegisterDestroy()
106: @*/
107: PetscErrorCode DMPlexTransformRegisterAll(void)
108: {
112: if (DMPlexTransformRegisterAllCalled) return(0);
113: DMPlexTransformRegisterAllCalled = PETSC_TRUE;
115: DMPlexTransformRegister(DMPLEXTRANSFORMFILTER, DMPlexTransformCreate_Filter);
116: DMPlexTransformRegister(DMPLEXREFINEREGULAR, DMPlexTransformCreate_Regular);
117: DMPlexTransformRegister(DMPLEXREFINETOBOX, DMPlexTransformCreate_ToBox);
118: DMPlexTransformRegister(DMPLEXREFINEALFELD, DMPlexTransformCreate_Alfeld);
119: DMPlexTransformRegister(DMPLEXREFINEBOUNDARYLAYER, DMPlexTransformCreate_BL);
120: DMPlexTransformRegister(DMPLEXREFINESBR, DMPlexTransformCreate_SBR);
121: return(0);
122: }
124: /*@C
125: DMPlexTransformRegisterDestroy - This function destroys the . It is called from PetscFinalize().
127: Level: developer
129: .seealso: PetscInitialize()
130: @*/
131: PetscErrorCode DMPlexTransformRegisterDestroy(void)
132: {
136: PetscFunctionListDestroy(&DMPlexTransformList);
137: DMPlexTransformRegisterAllCalled = PETSC_FALSE;
138: return(0);
139: }
141: /*@
142: DMPlexTransformCreate - Creates an empty transform object. The type can then be set with DMPlexTransformSetType().
144: Collective
146: Input Parameter:
147: . comm - The communicator for the transform object
149: Output Parameter:
150: . dm - The transform object
152: Level: beginner
154: .seealso: DMPlexTransformSetType(), DMPLEXREFINEREGULAR, DMPLEXTRANSFORMFILTER
155: @*/
156: PetscErrorCode DMPlexTransformCreate(MPI_Comm comm, DMPlexTransform *tr)
157: {
158: DMPlexTransform t;
159: PetscErrorCode ierr;
163: *tr = NULL;
164: DMInitializePackage();
166: PetscHeaderCreate(t, DMPLEXTRANSFORM_CLASSID, "DMPlexTransform", "Mesh Transform", "DMPlexTransform", comm, DMPlexTransformDestroy, DMPlexTransformView);
167: t->setupcalled = PETSC_FALSE;
168: PetscCalloc2(DM_NUM_POLYTOPES, &t->coordFE, DM_NUM_POLYTOPES, &t->refGeom);
169: *tr = t;
170: return(0);
171: }
173: /*@C
174: DMSetType - Sets the particular implementation for a transform.
176: Collective on tr
178: Input Parameters:
179: + tr - The transform
180: - method - The name of the transform type
182: Options Database Key:
183: . -dm_plex_transform_type <type> - Sets the transform type; use -help for a list of available types
185: Notes:
186: See "petsc/include/petscdmplextransform.h" for available transform types
188: Level: intermediate
190: .seealso: DMPlexTransformGetType(), DMPlexTransformCreate()
191: @*/
192: PetscErrorCode DMPlexTransformSetType(DMPlexTransform tr, DMPlexTransformType method)
193: {
194: PetscErrorCode (*r)(DMPlexTransform);
195: PetscBool match;
200: PetscObjectTypeCompare((PetscObject) tr, method, &match);
201: if (match) return(0);
203: DMPlexTransformRegisterAll();
204: PetscFunctionListFind(DMPlexTransformList, method, &r);
205: if (!r) SETERRQ1(PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMPlexTransform type: %s", method);
207: if (tr->ops->destroy) {(*tr->ops->destroy)(tr);}
208: PetscMemzero(tr->ops, sizeof(*tr->ops));
209: PetscObjectChangeTypeName((PetscObject) tr, method);
210: (*r)(tr);
211: return(0);
212: }
214: /*@C
215: DMPlexTransformGetType - Gets the type name (as a string) from the transform.
217: Not Collective
219: Input Parameter:
220: . tr - The DMPlexTransform
222: Output Parameter:
223: . type - The DMPlexTransform type name
225: Level: intermediate
227: .seealso: DMPlexTransformSetType(), DMPlexTransformCreate()
228: @*/
229: PetscErrorCode DMPlexTransformGetType(DMPlexTransform tr, DMPlexTransformType *type)
230: {
236: DMPlexTransformRegisterAll();
237: *type = ((PetscObject) tr)->type_name;
238: return(0);
239: }
241: static PetscErrorCode DMPlexTransformView_Ascii(DMPlexTransform tr, PetscViewer v)
242: {
243: PetscViewerFormat format;
244: PetscErrorCode ierr;
247: PetscViewerGetFormat(v, &format);
248: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
249: const PetscInt *trTypes = NULL;
250: IS trIS;
251: PetscInt cols = 8;
252: PetscInt Nrt = 8, f, g;
254: PetscViewerASCIIPrintf(v, "Source Starts\n");
255: for (g = 0; g <= cols; ++g) {PetscViewerASCIIPrintf(v, " % 14s", DMPolytopeTypes[g]);}
256: PetscViewerASCIIPrintf(v, "\n");
257: for (f = 0; f <= cols; ++f) {PetscViewerASCIIPrintf(v, " % 14d", tr->ctStart[f]);}
258: PetscViewerASCIIPrintf(v, "\n");
259: PetscViewerASCIIPrintf(v, "Target Starts\n");
260: for (g = 0; g <= cols; ++g) {PetscViewerASCIIPrintf(v, " % 14s", DMPolytopeTypes[g]);}
261: PetscViewerASCIIPrintf(v, "\n");
262: for (f = 0; f <= cols; ++f) {PetscViewerASCIIPrintf(v, " % 14d", tr->ctStartNew[f]);}
263: PetscViewerASCIIPrintf(v, "\n");
265: if (tr->trType) {
266: DMLabelGetNumValues(tr->trType, &Nrt);
267: DMLabelGetValueIS(tr->trType, &trIS);
268: ISGetIndices(trIS, &trTypes);
269: }
270: PetscViewerASCIIPrintf(v, "Offsets\n");
271: PetscViewerASCIIPrintf(v, " ");
272: for (g = 0; g < cols; ++g) {
273: PetscViewerASCIIPrintf(v, " % 14s", DMPolytopeTypes[g]);
274: }
275: PetscViewerASCIIPrintf(v, "\n");
276: for (f = 0; f < Nrt; ++f) {
277: PetscViewerASCIIPrintf(v, "%2d |", trTypes ? trTypes[f] : f);
278: for (g = 0; g < cols; ++g) {
279: PetscViewerASCIIPrintf(v, " % 14D", tr->offset[f*DM_NUM_POLYTOPES+g]);
280: }
281: PetscViewerASCIIPrintf(v, " |\n");
282: }
283: if (trTypes) {
284: ISGetIndices(trIS, &trTypes);
285: ISDestroy(&trIS);
286: }
287: }
288: return(0);
289: }
291: /*@C
292: DMPlexTransformView - Views a DMPlexTransform
294: Collective on tr
296: Input Parameters:
297: + tr - the DMPlexTransform object to view
298: - v - the viewer
300: Level: beginner
302: .seealso DMPlexTransformDestroy(), DMPlexTransformCreate()
303: @*/
304: PetscErrorCode DMPlexTransformView(DMPlexTransform tr, PetscViewer v)
305: {
306: PetscBool isascii;
311: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) tr), &v);}
314: PetscViewerCheckWritable(v);
315: PetscObjectPrintClassNamePrefixType((PetscObject) tr, v);
316: PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &isascii);
317: if (isascii) {DMPlexTransformView_Ascii(tr, v);}
318: if (tr->ops->view) {(*tr->ops->view)(tr, v);}
319: return(0);
320: }
322: /*@
323: DMPlexTransformSetFromOptions - Sets parameters in a transform from the options database
325: Collective on tr
327: Input Parameter:
328: . tr - the DMPlexTransform object to set options for
330: Options Database:
331: . -dm_plex_transform_type - Set the transform type, e.g. refine_regular
333: Level: intermediate
335: .seealso DMPlexTransformView(), DMPlexTransformCreate()
336: @*/
337: PetscErrorCode DMPlexTransformSetFromOptions(DMPlexTransform tr)
338: {
339: char typeName[1024];
340: const char *defName;
341: PetscBool flg;
346: DMPlexTransformGetType(tr, &defName);
347: if (!defName) defName = DMPLEXREFINEREGULAR;
348: PetscObjectOptionsBegin((PetscObject)tr);
349: PetscOptionsFList("-dm_plex_transform_type", "DMPlexTransform", "DMPlexTransformSetType", DMPlexTransformList, defName, typeName, 1024, &flg);
350: if (flg) {DMPlexTransformSetType(tr, typeName);}
351: else if (!((PetscObject) tr)->type_name) {DMPlexTransformSetType(tr, defName);}
352: if (tr->ops->setfromoptions) {(*tr->ops->setfromoptions)(PetscOptionsObject,tr);}
353: /* process any options handlers added with PetscObjectAddOptionsHandler() */
354: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) tr);
355: PetscOptionsEnd();
356: return(0);
357: }
359: /*@C
360: DMPlexTransformDestroy - Destroys a transform.
362: Collective on tr
364: Input Parameter:
365: . tr - the transform object to destroy
367: Level: beginner
369: .seealso DMPlexTransformView(), DMPlexTransformCreate()
370: @*/
371: PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr)
372: {
373: PetscInt c;
377: if (!*tr) return(0);
379: if (--((PetscObject) (*tr))->refct > 0) {*tr = NULL; return(0);}
381: if ((*tr)->ops->destroy) {
382: (*(*tr)->ops->destroy)(*tr);
383: }
384: DMDestroy(&(*tr)->dm);
385: DMLabelDestroy(&(*tr)->active);
386: DMLabelDestroy(&(*tr)->trType);
387: PetscFree2((*tr)->ctOrder, (*tr)->ctOrderInv);
388: PetscFree2((*tr)->ctStart, (*tr)->ctStartNew);
389: PetscFree((*tr)->offset);
390: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
391: PetscFEDestroy(&(*tr)->coordFE[c]);
392: PetscFEGeomDestroy(&(*tr)->refGeom[c]);
393: }
394: if ((*tr)->trVerts) {
395: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
396: DMPolytopeType *rct;
397: PetscInt *rsize, *rcone, *rornt, Nct, n, r;
399: if (DMPolytopeTypeGetDim((DMPolytopeType) c) > 0) {
400: DMPlexTransformCellTransform((*tr), (DMPolytopeType) c, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
401: for (n = 0; n < Nct; ++n) {
403: if (rct[n] == DM_POLYTOPE_POINT) continue;
404: for (r = 0; r < rsize[n]; ++r) {PetscFree((*tr)->trSubVerts[c][rct[n]][r]);}
405: PetscFree((*tr)->trSubVerts[c][rct[n]]);
406: }
407: }
408: PetscFree((*tr)->trSubVerts[c]);
409: PetscFree((*tr)->trVerts[c]);
410: }
411: }
412: PetscFree3((*tr)->trNv, (*tr)->trVerts, (*tr)->trSubVerts);
413: PetscFree2((*tr)->coordFE, (*tr)->refGeom);
414: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
415: PetscHeaderDestroy(tr);
416: return(0);
417: }
419: static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrder[], PetscInt ctStart[], PetscInt **offset)
420: {
421: DMLabel trType = tr->trType;
422: PetscInt c, cN, *off;
426: if (trType) {
427: DM dm;
428: IS rtIS;
429: const PetscInt *reftypes;
430: PetscInt Nrt, r;
432: DMPlexTransformGetDM(tr, &dm);
433: DMLabelGetNumValues(trType, &Nrt);
434: DMLabelGetValueIS(trType, &rtIS);
435: ISGetIndices(rtIS, &reftypes);
436: PetscCalloc1(Nrt*DM_NUM_POLYTOPES, &off);
437: for (r = 0; r < Nrt; ++r) {
438: const PetscInt rt = reftypes[r];
439: IS rtIS;
440: const PetscInt *points;
441: DMPolytopeType ct;
442: PetscInt p;
444: DMLabelGetStratumIS(trType, rt, &rtIS);
445: ISGetIndices(rtIS, &points);
446: p = points[0];
447: ISRestoreIndices(rtIS, &points);
448: ISDestroy(&rtIS);
449: DMPlexGetCellType(dm, p, &ct);
450: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
451: const DMPolytopeType ctNew = (DMPolytopeType) cN;
452: DMPolytopeType *rct;
453: PetscInt *rsize, *cone, *ornt;
454: PetscInt Nct, n, s;
456: if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {off[r*DM_NUM_POLYTOPES+ctNew] = -1; break;}
457: off[r*DM_NUM_POLYTOPES+ctNew] = 0;
458: for (s = 0; s <= r; ++s) {
459: const PetscInt st = reftypes[s];
460: DMPolytopeType sct;
461: PetscInt q, qrt;
463: DMLabelGetStratumIS(trType, st, &rtIS);
464: ISGetIndices(rtIS, &points);
465: q = points[0];
466: ISRestoreIndices(rtIS, &points);
467: ISDestroy(&rtIS);
468: DMPlexGetCellType(dm, q, &sct);
469: DMPlexTransformCellTransform(tr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt);
470: if (st != qrt) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Refine type %D of point %D does not match predicted type %D", qrt, q, st);
471: if (st == rt) {
472: for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) break;
473: if (n == Nct) off[r*DM_NUM_POLYTOPES+ctNew] = -1;
474: break;
475: }
476: for (n = 0; n < Nct; ++n) {
477: if (rct[n] == ctNew) {
478: PetscInt sn;
480: DMLabelGetStratumSize(trType, st, &sn);
481: off[r*DM_NUM_POLYTOPES+ctNew] += sn * rsize[n];
482: }
483: }
484: }
485: }
486: }
487: ISRestoreIndices(rtIS, &reftypes);
488: ISDestroy(&rtIS);
489: } else {
490: PetscCalloc1(DM_NUM_POLYTOPES*DM_NUM_POLYTOPES, &off);
491: for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
492: const DMPolytopeType ct = (DMPolytopeType) c;
493: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
494: const DMPolytopeType ctNew = (DMPolytopeType) cN;
495: DMPolytopeType *rct;
496: PetscInt *rsize, *cone, *ornt;
497: PetscInt Nct, n, i;
499: if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {off[ct*DM_NUM_POLYTOPES+ctNew] = -1; continue;}
500: off[ct*DM_NUM_POLYTOPES+ctNew] = 0;
501: for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
502: const DMPolytopeType ict = (DMPolytopeType) ctOrder[i];
503: const DMPolytopeType ictn = (DMPolytopeType) ctOrder[i+1];
505: DMPlexTransformCellTransform(tr, ict, PETSC_DETERMINE, NULL, &Nct, &rct, &rsize, &cone, &ornt);
506: if (ict == ct) {
507: for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) break;
508: if (n == Nct) off[ct*DM_NUM_POLYTOPES+ctNew] = -1;
509: break;
510: }
511: for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) off[ct*DM_NUM_POLYTOPES+ctNew] += (ctStart[ictn]-ctStart[ict]) * rsize[n];
512: }
513: }
514: }
515: }
516: *offset = off;
517: return(0);
518: }
520: PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr)
521: {
522: DM dm;
523: DMPolytopeType ctCell;
524: PetscInt pStart, pEnd, p, c;
529: if (tr->setupcalled) return(0);
530: if (tr->ops->setup) {(*tr->ops->setup)(tr);}
531: DMPlexTransformGetDM(tr, &dm);
532: DMPlexGetChart(dm, &pStart, &pEnd);
533: if (pEnd > pStart) {
534: DMPlexGetCellType(dm, 0, &ctCell);
535: } else {
536: PetscInt dim;
538: DMGetDimension(dm, &dim);
539: switch (dim) {
540: case 0: ctCell = DM_POLYTOPE_POINT;break;
541: case 1: ctCell = DM_POLYTOPE_SEGMENT;break;
542: case 2: ctCell = DM_POLYTOPE_TRIANGLE;break;
543: case 3: ctCell = DM_POLYTOPE_TETRAHEDRON;break;
544: default: ctCell = DM_POLYTOPE_UNKNOWN;
545: }
546: }
547: DMPlexCreateCellTypeOrder_Internal(ctCell, &tr->ctOrder, &tr->ctOrderInv);
548: /* Construct sizes and offsets for each cell type */
549: if (!tr->ctStart) {
550: PetscInt *ctS, *ctSN, *ctC, *ctCN;
552: PetscCalloc2(DM_NUM_POLYTOPES+1, &ctS, DM_NUM_POLYTOPES+1, &ctSN);
553: PetscCalloc2(DM_NUM_POLYTOPES+1, &ctC, DM_NUM_POLYTOPES+1, &ctCN);
554: for (p = pStart; p < pEnd; ++p) {
555: DMPolytopeType ct;
556: DMPolytopeType *rct;
557: PetscInt *rsize, *cone, *ornt;
558: PetscInt Nct, n;
560: DMPlexGetCellType(dm, p, &ct);
561: if (ct == DM_POLYTOPE_UNKNOWN) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %D", p);
562: ++ctC[ct];
563: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt);
564: for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
565: }
566: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
567: const PetscInt ct = tr->ctOrder[c];
568: const PetscInt ctn = tr->ctOrder[c+1];
570: ctS[ctn] = ctS[ct] + ctC[ct];
571: ctSN[ctn] = ctSN[ct] + ctCN[ct];
572: }
573: PetscFree2(ctC, ctCN);
574: tr->ctStart = ctS;
575: tr->ctStartNew = ctSN;
576: }
577: DMPlexTransformCreateOffset_Internal(tr, tr->ctOrder, tr->ctStart, &tr->offset);
578: tr->setupcalled = PETSC_TRUE;
579: return(0);
580: }
582: PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm)
583: {
587: *dm = tr->dm;
588: return(0);
589: }
591: PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm)
592: {
598: PetscObjectReference((PetscObject) dm);
599: DMDestroy(&tr->dm);
600: tr->dm = dm;
601: return(0);
602: }
604: PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active)
605: {
609: *active = tr->active;
610: return(0);
611: }
613: PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active)
614: {
620: PetscObjectReference((PetscObject) active);
621: DMLabelDestroy(&tr->active);
622: tr->active = active;
623: return(0);
624: }
626: static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe)
627: {
631: if (!tr->coordFE[ct]) {
632: PetscInt dim, cdim;
633: PetscBool isSimplex;
635: switch (ct) {
636: case DM_POLYTOPE_SEGMENT: dim = 1; isSimplex = PETSC_TRUE; break;
637: case DM_POLYTOPE_TRIANGLE: dim = 2; isSimplex = PETSC_TRUE; break;
638: case DM_POLYTOPE_QUADRILATERAL: dim = 2; isSimplex = PETSC_FALSE; break;
639: case DM_POLYTOPE_TETRAHEDRON: dim = 3; isSimplex = PETSC_TRUE; break;
640: case DM_POLYTOPE_HEXAHEDRON: dim = 3; isSimplex = PETSC_FALSE; break;
641: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No coordinate FE for cell type %s", DMPolytopeTypes[ct]);
642: }
643: DMGetCoordinateDim(tr->dm, &cdim);
644: PetscFECreateLagrange(PETSC_COMM_SELF, dim, cdim, isSimplex, 1, PETSC_DETERMINE, &tr->coordFE[ct]);
645: {
646: PetscDualSpace dsp;
647: PetscQuadrature quad;
648: DM K;
649: PetscFEGeom *cg;
650: PetscScalar *Xq;
651: PetscReal *xq, *wq;
652: PetscInt Nq, q;
654: DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq);
655: PetscMalloc1(Nq*cdim, &xq);
656: for (q = 0; q < Nq*cdim; ++q) xq[q] = PetscRealPart(Xq[q]);
657: PetscMalloc1(Nq, &wq);
658: for (q = 0; q < Nq; ++q) wq[q] = 1.0;
659: PetscQuadratureCreate(PETSC_COMM_SELF, &quad);
660: PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq);
661: PetscFESetQuadrature(tr->coordFE[ct], quad);
663: PetscFEGetDualSpace(tr->coordFE[ct], &dsp);
664: PetscDualSpaceGetDM(dsp, &K);
665: PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &tr->refGeom[ct]);
666: cg = tr->refGeom[ct];
667: DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
668: PetscQuadratureDestroy(&quad);
669: }
670: }
671: *fe = tr->coordFE[ct];
672: return(0);
673: }
675: /*@
676: DMPlexTransformGetTargetPoint - Get the number of a point in the transformed mesh based on information from the original mesh.
678: Not collective
680: Input Parameters:
681: + tr - The DMPlexTransform
682: . ct - The type of the original point which produces the new point
683: . ctNew - The type of the new point
684: . p - The original point which produces the new point
685: - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p
687: Output Parameters:
688: . pNew - The new point number
690: Level: developer
692: .seealso: DMPlexTransformGetSourcePoint(), DMPlexTransformCellTransform()
693: @*/
694: PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
695: {
696: DMPolytopeType *rct;
697: PetscInt *rsize, *cone, *ornt;
698: PetscInt rt, Nct, n, off, rp;
699: DMLabel trType = tr->trType;
700: PetscInt ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrder[tr->ctOrderInv[ct]+1]];
701: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrder[tr->ctOrderInv[ctNew]+1]];
702: PetscInt newp = ctSN, cind;
706: if ((p < ctS) || (p >= ctE)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D is not a %s [%D, %D)", p, DMPolytopeTypes[ct], ctS, ctE);
707: DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt);
708: if (trType) {
709: DMLabelGetValueIndex(trType, rt, &cind);
710: DMLabelGetStratumPointIndex(trType, rt, p, &rp);
711: if (rp < 0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s point %D does not have refine type %D", DMPolytopeTypes[ct], p, rt);
712: } else {
713: cind = ct;
714: rp = p - ctS;
715: }
716: off = tr->offset[cind*DM_NUM_POLYTOPES + ctNew];
717: if (off < 0) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s (%D) of point %D does not produce type %s for transform %s", DMPolytopeTypes[ct], rt, p, DMPolytopeTypes[ctNew], tr->hdr.type_name);
718: newp += off;
719: for (n = 0; n < Nct; ++n) {
720: if (rct[n] == ctNew) {
721: if (rsize[n] && r >= rsize[n])
722: SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number %D should be in [0, %D) for subcell type %s in cell type %s", r, rsize[n], DMPolytopeTypes[rct[n]], DMPolytopeTypes[ct]);
723: newp += rp * rsize[n] + r;
724: break;
725: }
726: }
728: if ((newp < ctSN) || (newp >= ctEN)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "New point %D is not a %s [%D, %D)", newp, DMPolytopeTypes[ctNew], ctSN, ctEN);
729: *pNew = newp;
730: return(0);
731: }
733: /*@
734: DMPlexTransformGetSourcePoint - Get the number of a point in the original mesh based on information from the transformed mesh.
736: Not collective
738: Input Parameters:
739: + tr - The DMPlexTransform
740: - pNew - The new point number
742: Output Parameters:
743: + ct - The type of the original point which produces the new point
744: . ctNew - The type of the new point
745: . p - The original point which produces the new point
746: - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p
748: Level: developer
750: .seealso: DMPlexTransformGetTargetPoint(), DMPlexTransformCellTransform()
751: @*/
752: PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r)
753: {
754: DMLabel trType = tr->trType;
755: DMPolytopeType *rct;
756: PetscInt *rsize, *cone, *ornt;
757: PetscInt rt, Nct, n, rp = 0, rO = 0, pO;
758: PetscInt offset = -1, ctS, ctE, ctO, ctN, ctTmp;
759: PetscErrorCode ierr;
762: for (ctN = 0; ctN < DM_NUM_POLYTOPES; ++ctN) {
763: PetscInt ctSN = tr->ctStartNew[ctN], ctEN = tr->ctStartNew[tr->ctOrder[tr->ctOrderInv[ctN]+1]];
765: if ((pNew >= ctSN) && (pNew < ctEN)) break;
766: }
767: if (ctN == DM_NUM_POLYTOPES) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cell type for target point %D could be not found", pNew);
768: if (trType) {
769: DM dm;
770: IS rtIS;
771: const PetscInt *reftypes;
772: PetscInt Nrt, r, rtStart;
774: DMPlexTransformGetDM(tr, &dm);
775: DMLabelGetNumValues(trType, &Nrt);
776: DMLabelGetValueIS(trType, &rtIS);
777: ISGetIndices(rtIS, &reftypes);
778: for (r = 0; r < Nrt; ++r) {
779: const PetscInt off = tr->offset[r*DM_NUM_POLYTOPES + ctN];
781: if (tr->ctStartNew[ctN] + off > pNew) continue;
782: /* Check that any of this refinement type exist */
783: /* TODO Actually keep track of the number produced here instead */
784: if (off > offset) {rt = reftypes[r]; offset = off;}
785: }
786: ISRestoreIndices(rtIS, &reftypes);
787: ISDestroy(&rtIS);
788: if (offset < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %D could be not found", pNew);
789: /* TODO Map refinement types to cell types */
790: DMLabelGetStratumBounds(trType, rt, &rtStart, NULL);
791: if (rtStart < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Refinement type %D has no source points", rt);
792: for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) {
793: PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrder[tr->ctOrderInv[ctO]+1]];
795: if ((rtStart >= ctS) && (rtStart < ctE)) break;
796: }
797: if (ctO == DM_NUM_POLYTOPES) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine a cell type for refinement type %D", rt);
798: } else {
799: for (ctTmp = 0; ctTmp < DM_NUM_POLYTOPES; ++ctTmp) {
800: const PetscInt off = tr->offset[ctTmp*DM_NUM_POLYTOPES + ctN];
802: if (tr->ctStartNew[ctN] + off > pNew) continue;
803: if (tr->ctStart[tr->ctOrder[tr->ctOrderInv[ctTmp]+1]] <= tr->ctStart[ctTmp]) continue;
804: /* TODO Actually keep track of the number produced here instead */
805: if (off > offset) {ctO = ctTmp; offset = off;}
806: }
807: if (offset < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %D could be not found", pNew);
808: }
809: ctS = tr->ctStart[ctO];
810: ctE = tr->ctStart[tr->ctOrder[tr->ctOrderInv[ctO]+1]];
811: DMPlexTransformCellTransform(tr, (DMPolytopeType) ctO, ctS, &rt, &Nct, &rct, &rsize, &cone, &ornt);
812: for (n = 0; n < Nct; ++n) {
813: if (rct[n] == ctN) {
814: PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset;
816: rp = tmp / rsize[n];
817: rO = tmp % rsize[n];
818: break;
819: }
820: }
821: if (n == Nct) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number for target point %D could be not found", pNew);
822: pO = rp + ctS;
823: if ((pO < ctS) || (pO >= ctE)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Source point %D is not a %s [%D, %D)", pO, DMPolytopeTypes[ctO], ctS, ctE);
824: if (ct) *ct = (DMPolytopeType) ctO;
825: if (ctNew) *ctNew = (DMPolytopeType) ctN;
826: if (p) *p = pO;
827: if (r) *r = rO;
828: return(0);
829: }
831: /*@
832: DMPlexTransformCellTransform - Describes the transform of a given source cell into a set of other target cells. These produced cells become the new mesh.
834: Input Parameters:
835: + tr - The DMPlexTransform object
836: . source - The source cell type
837: - p - The source point, which can also determine the refine type
839: Output Parameters:
840: + rt - The refine type for this point
841: . Nt - The number of types produced by this point
842: . target - An array of length Nt giving the types produced
843: . size - An array of length Nt giving the number of cells of each type produced
844: . cone - An array of length Nt*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point
845: - ornt - An array of length Nt*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point
847: Note: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
848: need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
849: division (described in these outputs) of the cell in the original mesh. The point identifier is given by
850: $ the number of cones to be taken, or 0 for the current cell
851: $ the cell cone point number at each level from which it is subdivided
852: $ the replica number r of the subdivision.
853: The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
854: $ Nt = 2
855: $ target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
856: $ size = {1, 2}
857: $ cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
858: $ ornt = { 0, 0, 0, 0}
860: Level: advanced
862: .seealso: DMPlexTransformApply(), DMPlexTransformCreate()
863: @*/
864: PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
865: {
869: (*tr->ops->celltransform)(tr, source, p, rt, Nt, target, size, cone, ornt);
870: return(0);
871: }
873: PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
874: {
876: *rnew = r;
877: *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
878: return(0);
879: }
881: /* Returns the same thing */
882: PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
883: {
884: static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
885: static PetscInt vertexS[] = {1};
886: static PetscInt vertexC[] = {0};
887: static PetscInt vertexO[] = {0};
888: static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT};
889: static PetscInt edgeS[] = {1};
890: static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
891: static PetscInt edgeO[] = {0, 0};
892: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
893: static PetscInt tedgeS[] = {1};
894: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
895: static PetscInt tedgeO[] = {0, 0};
896: static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE};
897: static PetscInt triS[] = {1};
898: static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
899: static PetscInt triO[] = {0, 0, 0};
900: static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL};
901: static PetscInt quadS[] = {1};
902: static PetscInt quadC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0, DM_POLYTOPE_SEGMENT, 1, 3, 0};
903: static PetscInt quadO[] = {0, 0, 0, 0};
904: static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR};
905: static PetscInt tquadS[] = {1};
906: static PetscInt tquadC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_POINT_PRISM_TENSOR, 1, 3, 0};
907: static PetscInt tquadO[] = {0, 0, 0, 0};
908: static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON};
909: static PetscInt tetS[] = {1};
910: static PetscInt tetC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 3, 0};
911: static PetscInt tetO[] = {0, 0, 0, 0};
912: static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON};
913: static PetscInt hexS[] = {1};
914: static PetscInt hexC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0,
915: DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 1, 5, 0};
916: static PetscInt hexO[] = {0, 0, 0, 0, 0, 0};
917: static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM};
918: static PetscInt tripS[] = {1};
919: static PetscInt tripC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0,
920: DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0};
921: static PetscInt tripO[] = {0, 0, 0, 0, 0};
922: static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR};
923: static PetscInt ttripS[] = {1};
924: static PetscInt ttripC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0,
925: DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0};
926: static PetscInt ttripO[] = {0, 0, 0, 0, 0};
927: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
928: static PetscInt tquadpS[] = {1};
929: static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0,
930: DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
931: static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0};
932: static DMPolytopeType pyrT[] = {DM_POLYTOPE_PYRAMID};
933: static PetscInt pyrS[] = {1};
934: static PetscInt pyrC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0,
935: DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 1, 4, 0};
936: static PetscInt pyrO[] = {0, 0, 0, 0, 0};
939: if (rt) *rt = 0;
940: switch (source) {
941: case DM_POLYTOPE_POINT: *Nt = 1; *target = vertexT; *size = vertexS; *cone = vertexC; *ornt = vertexO; break;
942: case DM_POLYTOPE_SEGMENT: *Nt = 1; *target = edgeT; *size = edgeS; *cone = edgeC; *ornt = edgeO; break;
943: case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break;
944: case DM_POLYTOPE_TRIANGLE: *Nt = 1; *target = triT; *size = triS; *cone = triC; *ornt = triO; break;
945: case DM_POLYTOPE_QUADRILATERAL: *Nt = 1; *target = quadT; *size = quadS; *cone = quadC; *ornt = quadO; break;
946: case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 1; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break;
947: case DM_POLYTOPE_TETRAHEDRON: *Nt = 1; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break;
948: case DM_POLYTOPE_HEXAHEDRON: *Nt = 1; *target = hexT; *size = hexS; *cone = hexC; *ornt = hexO; break;
949: case DM_POLYTOPE_TRI_PRISM: *Nt = 1; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break;
950: case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 1; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break;
951: case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 1; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break;
952: case DM_POLYTOPE_PYRAMID: *Nt = 1; *target = pyrT; *size = pyrS; *cone = pyrC; *ornt = pyrO; break;
953: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
954: }
955: return(0);
956: }
958: /*@
959: DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point
961: Not collective
963: Input Parameters:
964: + tr - The DMPlexTransform
965: . sct - The source point cell type, from whom the new cell is being produced
966: . sp - The source point
967: . so - The orientation of the source point in its enclosing parent
968: . tct - The target point cell type
969: . r - The replica number requested for the produced cell type
970: - o - The orientation of the replica
972: Output Parameters:
973: + rnew - The replica number, given the orientation of the parent
974: - onew - The replica orientation, given the orientation of the parent
976: Level: advanced
978: .seealso: DMPlexTransformCellTransform(), DMPlexTransformApply()
979: @*/
980: PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
981: {
985: (*tr->ops->getsubcellorientation)(tr, sct, sp, so, tct, r, o, rnew, onew);
986: return(0);
987: }
989: static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm)
990: {
991: DM dm;
992: PetscInt pStart, pEnd, p, pNew;
993: PetscErrorCode ierr;
996: DMPlexTransformGetDM(tr, &dm);
997: /* Must create the celltype label here so that we do not automatically try to compute the types */
998: DMCreateLabel(rdm, "celltype");
999: DMPlexGetChart(dm, &pStart, &pEnd);
1000: for (p = pStart; p < pEnd; ++p) {
1001: DMPolytopeType ct;
1002: DMPolytopeType *rct;
1003: PetscInt *rsize, *rcone, *rornt;
1004: PetscInt Nct, n, r;
1006: DMPlexGetCellType(dm, p, &ct);
1007: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1008: for (n = 0; n < Nct; ++n) {
1009: for (r = 0; r < rsize[n]; ++r) {
1010: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1011: DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]));
1012: DMPlexSetCellType(rdm, pNew, rct[n]);
1013: }
1014: }
1015: }
1016: /* Let the DM know we have set all the cell types */
1017: {
1018: DMLabel ctLabel;
1019: DM_Plex *plex = (DM_Plex *) rdm->data;
1021: DMPlexGetCellTypeLabel(rdm, &ctLabel);
1022: PetscObjectStateGet((PetscObject) ctLabel, &plex->celltypeState);
1023: }
1024: return(0);
1025: }
1027: #if 0
1028: static PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize)
1029: {
1030: PetscInt ctNew;
1035: /* TODO Can do bisection since everything is sorted */
1036: for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) {
1037: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrder[tr->ctOrderInv[ctNew]+1]];
1039: if (q >= ctSN && q < ctEN) break;
1040: }
1041: if (ctNew >= DM_NUM_POLYTOPES) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %D cannot be located in the transformed mesh", q);
1042: *coneSize = DMPolytopeTypeGetConeSize((DMPolytopeType) ctNew);
1043: return(0);
1044: }
1045: #endif
1047: /* The orientation o is for the interior of the cell p */
1048: static PetscErrorCode DMPlexTransformGetCone_Internal(DMPlexTransform tr, PetscInt p, PetscInt o, DMPolytopeType ct, DMPolytopeType ctNew,
1049: const PetscInt rcone[], PetscInt *coneoff, const PetscInt rornt[], PetscInt *orntoff,
1050: PetscInt coneNew[], PetscInt orntNew[])
1051: {
1052: DM dm;
1053: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1054: const PetscInt *cone;
1055: PetscInt c, coff = *coneoff, ooff = *orntoff;
1056: PetscErrorCode ierr;
1059: DMPlexTransformGetDM(tr, &dm);
1060: DMPlexGetCone(dm, p, &cone);
1061: for (c = 0; c < csizeNew; ++c) {
1062: PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */
1063: PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */
1064: PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */
1065: DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */
1066: const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */
1067: PetscInt pr = -1; /* Replica number of pp that produces new cone point */
1068: const DMPolytopeType ft = (DMPolytopeType) rcone[coff++]; /* Cell type for new cone point of pNew */
1069: const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1070: PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */
1071: PetscInt lc;
1073: /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
1074: for (lc = 0; lc < fn; ++lc) {
1075: const PetscInt *parr = DMPolytopeTypeGetArrangment(pct, po);
1076: const PetscInt acp = rcone[coff++];
1077: const PetscInt pcp = parr[acp*2];
1078: const PetscInt pco = parr[acp*2+1];
1079: const PetscInt *ppornt;
1081: ppp = pp;
1082: pp = pcone[pcp];
1083: DMPlexGetCellType(dm, pp, &pct);
1084: DMPlexGetCone(dm, pp, &pcone);
1085: DMPlexGetConeOrientation(dm, ppp, &ppornt);
1086: po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco);
1087: }
1088: pr = rcone[coff++];
1089: /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
1090: DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo);
1091: DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c]);
1092: orntNew[c] = fo;
1093: }
1094: *coneoff = coff;
1095: *orntoff = ooff;
1096: return(0);
1097: }
1099: static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm)
1100: {
1101: DM dm;
1102: DMPolytopeType ct;
1103: PetscInt *coneNew, *orntNew;
1104: PetscInt maxConeSize = 0, pStart, pEnd, p, pNew;
1108: DMPlexTransformGetDM(tr, &dm);
1109: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p));
1110: DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew);
1111: DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew);
1112: DMPlexGetChart(dm, &pStart, &pEnd);
1113: for (p = pStart; p < pEnd; ++p) {
1114: const PetscInt *cone, *ornt;
1115: PetscInt coff, ooff;
1116: DMPolytopeType *rct;
1117: PetscInt *rsize, *rcone, *rornt;
1118: PetscInt Nct, n, r;
1120: DMPlexGetCellType(dm, p, &ct);
1121: DMPlexGetCone(dm, p, &cone);
1122: DMPlexGetConeOrientation(dm, p, &ornt);
1123: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1124: for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
1125: const DMPolytopeType ctNew = rct[n];
1127: for (r = 0; r < rsize[n]; ++r) {
1128: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1129: DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew);
1130: DMPlexSetCone(rdm, pNew, coneNew);
1131: DMPlexSetConeOrientation(rdm, pNew, orntNew);
1132: }
1133: }
1134: }
1135: DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew);
1136: DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew);
1137: DMViewFromOptions(rdm, NULL, "-rdm_view");
1138: DMPlexSymmetrize(rdm);
1139: DMPlexStratify(rdm);
1140: return(0);
1141: }
1143: PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[])
1144: {
1145: DM dm;
1146: DMPolytopeType ct, qct;
1147: DMPolytopeType *rct;
1148: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1149: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1150: PetscErrorCode ierr;
1156: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p));
1157: DMPlexTransformGetDM(tr, &dm);
1158: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone);
1159: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt);
1160: DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r);
1161: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1162: for (n = 0; n < Nct; ++n) {
1163: const DMPolytopeType ctNew = rct[n];
1164: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1165: PetscInt Nr = rsize[n], fn, c;
1167: if (ctNew == qct) Nr = r;
1168: for (nr = 0; nr < Nr; ++nr) {
1169: for (c = 0; c < csizeNew; ++c) {
1170: ++coff; /* Cell type of new cone point */
1171: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1172: coff += fn;
1173: ++coff; /* Replica number of new cone point */
1174: ++ooff; /* Orientation of new cone point */
1175: }
1176: }
1177: if (ctNew == qct) break;
1178: }
1179: DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt);
1180: *cone = qcone;
1181: *ornt = qornt;
1182: return(0);
1183: }
1185: PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1186: {
1187: DM dm;
1188: DMPolytopeType ct, qct;
1189: DMPolytopeType *rct;
1190: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1191: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1192: PetscErrorCode ierr;
1197: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p));
1198: DMPlexTransformGetDM(tr, &dm);
1199: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone);
1200: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt);
1201: DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r);
1202: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1203: for (n = 0; n < Nct; ++n) {
1204: const DMPolytopeType ctNew = rct[n];
1205: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1206: PetscInt Nr = rsize[n], fn, c;
1208: if (ctNew == qct) Nr = r;
1209: for (nr = 0; nr < Nr; ++nr) {
1210: for (c = 0; c < csizeNew; ++c) {
1211: ++coff; /* Cell type of new cone point */
1212: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1213: coff += fn;
1214: ++coff; /* Replica number of new cone point */
1215: ++ooff; /* Orientation of new cone point */
1216: }
1217: }
1218: if (ctNew == qct) break;
1219: }
1220: DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt);
1221: *cone = qcone;
1222: *ornt = qornt;
1223: return(0);
1224: }
1226: PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1227: {
1228: DM dm;
1234: DMPlexTransformGetDM(tr, &dm);
1235: DMRestoreWorkArray(dm, 0, MPIU_INT, cone);
1236: DMRestoreWorkArray(dm, 0, MPIU_INT, ornt);
1237: return(0);
1238: }
1240: static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr)
1241: {
1242: PetscInt ict;
1246: PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts);
1247: for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) {
1248: const DMPolytopeType ct = (DMPolytopeType) ict;
1249: DMPlexTransform reftr;
1250: DM refdm, trdm;
1251: Vec coordinates;
1252: const PetscScalar *coords;
1253: DMPolytopeType *rct;
1254: PetscInt *rsize, *rcone, *rornt;
1255: PetscInt Nct, n, r, pNew;
1256: PetscInt vStart, vEnd, Nc;
1257: const PetscInt debug = 0;
1258: const char *typeName;
1260: /* Since points are 0-dimensional, coordinates make no sense */
1261: if (DMPolytopeTypeGetDim(ct) <= 0) continue;
1262: DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm);
1263: DMPlexTransformCreate(PETSC_COMM_SELF, &reftr);
1264: DMPlexTransformSetDM(reftr, refdm);
1265: DMPlexTransformGetType(tr, &typeName);
1266: DMPlexTransformSetType(reftr, typeName);
1267: DMPlexTransformSetUp(reftr);
1268: DMPlexTransformApply(reftr, refdm, &trdm);
1270: DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd);
1271: tr->trNv[ct] = vEnd - vStart;
1272: DMGetCoordinatesLocal(trdm, &coordinates);
1273: VecGetLocalSize(coordinates, &Nc);
1274: if (tr->trNv[ct] * DMPolytopeTypeGetDim(ct) != Nc) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell type %s, transformed coordinate size %D != %D size of coordinate storage", DMPolytopeTypes[ct], tr->trNv[ct] * DMPolytopeTypeGetDim(ct), Nc);
1275: PetscCalloc1(Nc, &tr->trVerts[ct]);
1276: VecGetArrayRead(coordinates, &coords);
1277: PetscArraycpy(tr->trVerts[ct], coords, Nc);
1278: VecRestoreArrayRead(coordinates, &coords);
1280: PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct]);
1281: DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1282: for (n = 0; n < Nct; ++n) {
1284: /* Since points are 0-dimensional, coordinates make no sense */
1285: if (rct[n] == DM_POLYTOPE_POINT) continue;
1286: PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]]);
1287: for (r = 0; r < rsize[n]; ++r) {
1288: PetscInt *closure = NULL;
1289: PetscInt clSize, cl, Nv = 0;
1291: PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r]);
1292: DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew);
1293: DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure);
1294: for (cl = 0; cl < clSize*2; cl += 2) {
1295: const PetscInt sv = closure[cl];
1297: if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart;
1298: }
1299: DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure);
1300: if (Nv != DMPolytopeTypeGetNumVertices(rct[n])) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of vertices %D != %D for %s subcell %D from cell %s", Nv, DMPolytopeTypeGetNumVertices(rct[n]), DMPolytopeTypes[rct[n]], r, DMPolytopeTypes[ct]);
1301: }
1302: }
1303: if (debug) {
1304: DMPolytopeType *rct;
1305: PetscInt *rsize, *rcone, *rornt;
1306: PetscInt v, dE = DMPolytopeTypeGetDim(ct), d, off = 0;
1308: PetscPrintf(PETSC_COMM_SELF, "%s: %D vertices\n", DMPolytopeTypes[ct], tr->trNv[ct]);
1309: for (v = 0; v < tr->trNv[ct]; ++v) {
1310: PetscPrintf(PETSC_COMM_SELF, " ");
1311: for (d = 0; d < dE; ++d) {PetscPrintf(PETSC_COMM_SELF, "%g ", tr->trVerts[ct][off++]);}
1312: PetscPrintf(PETSC_COMM_SELF, "\n");
1313: }
1315: DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1316: for (n = 0; n < Nct; ++n) {
1317: if (rct[n] == DM_POLYTOPE_POINT) continue;
1318: PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct]);
1319: for (r = 0; r < rsize[n]; ++r) {
1320: PetscPrintf(PETSC_COMM_SELF, " ");
1321: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) {
1322: PetscPrintf(PETSC_COMM_SELF, "%D ", tr->trSubVerts[ct][rct[n]][r][v]);
1323: }
1324: PetscPrintf(PETSC_COMM_SELF, "\n");
1325: }
1326: }
1327: }
1328: DMDestroy(&refdm);
1329: DMDestroy(&trdm);
1330: DMPlexTransformDestroy(&reftr);
1331: }
1332: return(0);
1333: }
1335: /*
1336: DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type
1338: Input Parameters:
1339: + tr - The DMPlexTransform object
1340: - ct - The cell type
1342: Output Parameters:
1343: + Nv - The number of transformed vertices in the closure of the reference cell of given type
1344: - trVerts - The coordinates of these vertices in the reference cell
1346: Level: developer
1348: .seealso: DMPlexTransformGetSubcellVertices()
1349: */
1350: PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[])
1351: {
1355: if (!tr->trNv) {DMPlexTransformCreateCellVertices_Internal(tr);}
1356: if (Nv) *Nv = tr->trNv[ct];
1357: if (trVerts) *trVerts = tr->trVerts[ct];
1358: return(0);
1359: }
1361: /*
1362: DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type
1364: Input Parameters:
1365: + tr - The DMPlexTransform object
1366: . ct - The cell type
1367: . rct - The subcell type
1368: - r - The subcell index
1370: Output Parameter:
1371: . subVerts - The indices of these vertices in the set of vertices returned by DMPlexTransformGetCellVertices()
1373: Level: developer
1375: .seealso: DMPlexTransformGetCellVertices()
1376: */
1377: PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[])
1378: {
1382: if (!tr->trNv) {DMPlexTransformCreateCellVertices_Internal(tr);}
1383: if (!tr->trSubVerts[ct][rct]) SETERRQ2(PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
1384: if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r];
1385: return(0);
1386: }
1388: /* Computes new vertex as the barycenter, or centroid */
1389: PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1390: {
1391: PetscInt v,d;
1394: if (ct != DM_POLYTOPE_POINT) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for refined point type %s",DMPolytopeTypes[ct]);
1395: for (d = 0; d < dE; ++d) out[d] = 0.0;
1396: for (v = 0; v < Nv; ++v) for (d = 0; d < dE; ++d) out[d] += in[v*dE+d];
1397: for (d = 0; d < dE; ++d) out[d] /= Nv;
1398: return(0);
1399: }
1401: /*@
1402: DMPlexTransformMapCoordinates -
1404: Input Parameters:
1405: + tr - The DMPlexTransform
1406: . pct - The cell type of the parent, from whom the new cell is being produced
1407: . ct - The type being produced
1408: . r - The replica number requested for the produced cell type
1409: . Nv - Number of vertices in the closure of the parent cell
1410: . dE - Spatial dimension
1411: - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell
1413: Output Parameter:
1414: . out - The coordinates of the new vertices
1416: Level: intermediate
1417: @*/
1418: PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1419: {
1423: (*tr->ops->mapcoordinates)(tr, pct, ct, r, Nv, dE, in, out);
1424: return(0);
1425: }
1427: static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew)
1428: {
1429: DM dm;
1430: IS valueIS;
1431: const PetscInt *values;
1432: PetscInt defVal, Nv, val;
1433: PetscErrorCode ierr;
1436: DMPlexTransformGetDM(tr, &dm);
1437: DMLabelGetDefaultValue(label, &defVal);
1438: DMLabelSetDefaultValue(labelNew, defVal);
1439: DMLabelGetValueIS(label, &valueIS);
1440: ISGetLocalSize(valueIS, &Nv);
1441: ISGetIndices(valueIS, &values);
1442: for (val = 0; val < Nv; ++val) {
1443: IS pointIS;
1444: const PetscInt *points;
1445: PetscInt numPoints, p;
1447: /* Ensure refined label is created with same number of strata as
1448: * original (even if no entries here). */
1449: DMLabelAddStratum(labelNew, values[val]);
1450: DMLabelGetStratumIS(label, values[val], &pointIS);
1451: ISGetLocalSize(pointIS, &numPoints);
1452: ISGetIndices(pointIS, &points);
1453: for (p = 0; p < numPoints; ++p) {
1454: const PetscInt point = points[p];
1455: DMPolytopeType ct;
1456: DMPolytopeType *rct;
1457: PetscInt *rsize, *rcone, *rornt;
1458: PetscInt Nct, n, r, pNew=0;
1460: DMPlexGetCellType(dm, point, &ct);
1461: DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1462: for (n = 0; n < Nct; ++n) {
1463: for (r = 0; r < rsize[n]; ++r) {
1464: DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew);
1465: DMLabelSetValue(labelNew, pNew, values[val]);
1466: }
1467: }
1468: }
1469: ISRestoreIndices(pointIS, &points);
1470: ISDestroy(&pointIS);
1471: }
1472: ISRestoreIndices(valueIS, &values);
1473: ISDestroy(&valueIS);
1474: return(0);
1475: }
1477: static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm)
1478: {
1479: DM dm;
1480: PetscInt numLabels, l;
1484: DMPlexTransformGetDM(tr, &dm);
1485: DMGetNumLabels(dm, &numLabels);
1486: for (l = 0; l < numLabels; ++l) {
1487: DMLabel label, labelNew;
1488: const char *lname;
1489: PetscBool isDepth, isCellType;
1491: DMGetLabelName(dm, l, &lname);
1492: PetscStrcmp(lname, "depth", &isDepth);
1493: if (isDepth) continue;
1494: PetscStrcmp(lname, "celltype", &isCellType);
1495: if (isCellType) continue;
1496: DMCreateLabel(rdm, lname);
1497: DMGetLabel(dm, lname, &label);
1498: DMGetLabel(rdm, lname, &labelNew);
1499: RefineLabel_Internal(tr, label, labelNew);
1500: }
1501: return(0);
1502: }
1504: /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */
1505: PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm)
1506: {
1507: DM dm;
1508: PetscInt Nf, f, Nds, s;
1512: DMPlexTransformGetDM(tr, &dm);
1513: DMGetNumFields(dm, &Nf);
1514: for (f = 0; f < Nf; ++f) {
1515: DMLabel label, labelNew;
1516: PetscObject obj;
1517: const char *lname;
1519: DMGetField(rdm, f, &label, &obj);
1520: if (!label) continue;
1521: PetscObjectGetName((PetscObject) label, &lname);
1522: DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
1523: RefineLabel_Internal(tr, label, labelNew);
1524: DMSetField_Internal(rdm, f, labelNew, obj);
1525: DMLabelDestroy(&labelNew);
1526: }
1527: DMGetNumDS(dm, &Nds);
1528: for (s = 0; s < Nds; ++s) {
1529: DMLabel label, labelNew;
1530: const char *lname;
1532: DMGetRegionNumDS(rdm, s, &label, NULL, NULL);
1533: if (!label) continue;
1534: PetscObjectGetName((PetscObject) label, &lname);
1535: DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
1536: RefineLabel_Internal(tr, label, labelNew);
1537: DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL);
1538: DMLabelDestroy(&labelNew);
1539: }
1540: return(0);
1541: }
1543: static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm)
1544: {
1545: DM dm;
1546: PetscSF sf, sfNew;
1547: PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m;
1548: const PetscInt *localPoints;
1549: const PetscSFNode *remotePoints;
1550: PetscInt *localPointsNew;
1551: PetscSFNode *remotePointsNew;
1552: PetscInt pStartNew, pEndNew, pNew;
1553: /* Brute force algorithm */
1554: PetscSF rsf;
1555: PetscSection s;
1556: const PetscInt *rootdegree;
1557: PetscInt *rootPointsNew, *remoteOffsets;
1558: PetscInt numPointsNew, pStart, pEnd, p;
1559: PetscErrorCode ierr;
1562: DMPlexTransformGetDM(tr, &dm);
1563: DMPlexGetChart(rdm, &pStartNew, &pEndNew);
1564: DMGetPointSF(dm, &sf);
1565: DMGetPointSF(rdm, &sfNew);
1566: /* Calculate size of new SF */
1567: PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);
1568: if (numRoots < 0) return(0);
1569: for (l = 0; l < numLeaves; ++l) {
1570: const PetscInt p = localPoints[l];
1571: DMPolytopeType ct;
1572: DMPolytopeType *rct;
1573: PetscInt *rsize, *rcone, *rornt;
1574: PetscInt Nct, n;
1576: DMPlexGetCellType(dm, p, &ct);
1577: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1578: for (n = 0; n < Nct; ++n) {
1579: numLeavesNew += rsize[n];
1580: }
1581: }
1582: /* Send new root point numbers
1583: It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication
1584: */
1585: DMPlexGetChart(dm, &pStart, &pEnd);
1586: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &s);
1587: PetscSectionSetChart(s, pStart, pEnd);
1588: for (p = pStart; p < pEnd; ++p) {
1589: DMPolytopeType ct;
1590: DMPolytopeType *rct;
1591: PetscInt *rsize, *rcone, *rornt;
1592: PetscInt Nct, n;
1594: DMPlexGetCellType(dm, p, &ct);
1595: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1596: for (n = 0; n < Nct; ++n) {
1597: PetscSectionAddDof(s, p, rsize[n]);
1598: }
1599: }
1600: PetscSectionSetUp(s);
1601: PetscSectionGetStorageSize(s, &numPointsNew);
1602: PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets);
1603: PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf);
1604: PetscFree(remoteOffsets);
1605: PetscSFComputeDegreeBegin(sf, &rootdegree);
1606: PetscSFComputeDegreeEnd(sf, &rootdegree);
1607: PetscMalloc1(numPointsNew, &rootPointsNew);
1608: for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1;
1609: for (p = pStart; p < pEnd; ++p) {
1610: DMPolytopeType ct;
1611: DMPolytopeType *rct;
1612: PetscInt *rsize, *rcone, *rornt;
1613: PetscInt Nct, n, r, off;
1615: if (!rootdegree[p-pStart]) continue;
1616: PetscSectionGetOffset(s, p, &off);
1617: DMPlexGetCellType(dm, p, &ct);
1618: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1619: for (n = 0, m = 0; n < Nct; ++n) {
1620: for (r = 0; r < rsize[n]; ++r, ++m) {
1621: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1622: rootPointsNew[off+m] = pNew;
1623: }
1624: }
1625: }
1626: PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew,MPI_REPLACE);
1627: PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew,MPI_REPLACE);
1628: PetscSFDestroy(&rsf);
1629: PetscMalloc1(numLeavesNew, &localPointsNew);
1630: PetscMalloc1(numLeavesNew, &remotePointsNew);
1631: for (l = 0, m = 0; l < numLeaves; ++l) {
1632: const PetscInt p = localPoints[l];
1633: DMPolytopeType ct;
1634: DMPolytopeType *rct;
1635: PetscInt *rsize, *rcone, *rornt;
1636: PetscInt Nct, n, r, q, off;
1638: PetscSectionGetOffset(s, p, &off);
1639: DMPlexGetCellType(dm, p, &ct);
1640: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1641: for (n = 0, q = 0; n < Nct; ++n) {
1642: for (r = 0; r < rsize[n]; ++r, ++m, ++q) {
1643: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1644: localPointsNew[m] = pNew;
1645: remotePointsNew[m].index = rootPointsNew[off+q];
1646: remotePointsNew[m].rank = remotePoints[l].rank;
1647: }
1648: }
1649: }
1650: PetscSectionDestroy(&s);
1651: PetscFree(rootPointsNew);
1652: /* SF needs sorted leaves to correctly calculate Gather */
1653: {
1654: PetscSFNode *rp, *rtmp;
1655: PetscInt *lp, *idx, *ltmp, i;
1657: PetscMalloc1(numLeavesNew, &idx);
1658: PetscMalloc1(numLeavesNew, &lp);
1659: PetscMalloc1(numLeavesNew, &rp);
1660: for (i = 0; i < numLeavesNew; ++i) {
1661: if ((localPointsNew[i] < pStartNew) || (localPointsNew[i] >= pEndNew)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local SF point %D (%D) not in [%D, %D)", localPointsNew[i], i, pStartNew, pEndNew);
1662: idx[i] = i;
1663: }
1664: PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx);
1665: for (i = 0; i < numLeavesNew; ++i) {
1666: lp[i] = localPointsNew[idx[i]];
1667: rp[i] = remotePointsNew[idx[i]];
1668: }
1669: ltmp = localPointsNew;
1670: localPointsNew = lp;
1671: rtmp = remotePointsNew;
1672: remotePointsNew = rp;
1673: PetscFree(idx);
1674: PetscFree(ltmp);
1675: PetscFree(rtmp);
1676: }
1677: PetscSFSetGraph(sfNew, pEndNew-pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
1678: return(0);
1679: }
1681: /*
1682: DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of type ct with localized coordinates x, we generate localized coordinates xr for subcell r of type rct.
1684: Not collective
1686: Input Parameters:
1687: + tr - The DMPlexTransform
1688: . ct - The type of the parent cell
1689: . rct - The type of the produced cell
1690: . r - The index of the produced cell
1691: - x - The localized coordinates for the parent cell
1693: Output Parameter:
1694: . xr - The localized coordinates for the produced cell
1696: Level: developer
1698: .seealso: DMPlexCellRefinerSetCoordinates()
1699: */
1700: static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
1701: {
1702: PetscFE fe = NULL;
1703: PetscInt cdim, v, *subcellV;
1707: DMPlexTransformGetCoordinateFE(tr, ct, &fe);
1708: DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV);
1709: PetscFEGetNumComponents(fe, &cdim);
1710: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) {
1711: PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v*cdim]);
1712: }
1713: return(0);
1714: }
1716: static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm)
1717: {
1718: DM dm, cdm;
1719: PetscSection coordSection, coordSectionNew;
1720: Vec coordsLocal, coordsLocalNew;
1721: const PetscScalar *coords;
1722: PetscScalar *coordsNew;
1723: const DMBoundaryType *bd;
1724: const PetscReal *maxCell, *L;
1725: PetscBool isperiodic, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE;
1726: PetscInt dE, d, cStart, cEnd, c, vStartNew, vEndNew, v, pStart, pEnd, p, ocStart, ocEnd;
1727: PetscErrorCode ierr;
1730: DMPlexTransformGetDM(tr, &dm);
1731: DMGetCoordinateDM(dm, &cdm);
1732: DMGetPeriodicity(dm, &isperiodic, &maxCell, &L, &bd);
1733: /* Determine if we need to localize coordinates when generating them */
1734: if (isperiodic) {
1735: localizeVertices = PETSC_TRUE;
1736: if (!maxCell) {
1737: PetscBool localized;
1738: DMGetCoordinatesLocalized(dm, &localized);
1739: if (!localized) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Cannot refine a periodic mesh if coordinates have not been localized");
1740: localizeCells = PETSC_TRUE;
1741: }
1742: }
1744: DMGetCoordinateSection(dm, &coordSection);
1745: PetscSectionGetFieldComponents(coordSection, 0, &dE);
1746: if (maxCell) {
1747: PetscReal maxCellNew[3];
1749: for (d = 0; d < dE; ++d) maxCellNew[d] = maxCell[d]/2.0;
1750: DMSetPeriodicity(rdm, isperiodic, maxCellNew, L, bd);
1751: } else {
1752: DMSetPeriodicity(rdm, isperiodic, maxCell, L, bd);
1753: }
1754: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &coordSectionNew);
1755: PetscSectionSetNumFields(coordSectionNew, 1);
1756: PetscSectionSetFieldComponents(coordSectionNew, 0, dE);
1757: DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew);
1758: if (localizeCells) {PetscSectionSetChart(coordSectionNew, 0, vEndNew);}
1759: else {PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew);}
1761: /* Localization should be inherited */
1762: /* Stefano calculates parent cells for each new cell for localization */
1763: /* Localized cells need coordinates of closure */
1764: for (v = vStartNew; v < vEndNew; ++v) {
1765: PetscSectionSetDof(coordSectionNew, v, dE);
1766: PetscSectionSetFieldDof(coordSectionNew, v, 0, dE);
1767: }
1768: if (localizeCells) {
1769: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1770: for (c = cStart; c < cEnd; ++c) {
1771: PetscInt dof;
1773: PetscSectionGetDof(coordSection, c, &dof);
1774: if (dof) {
1775: DMPolytopeType ct;
1776: DMPolytopeType *rct;
1777: PetscInt *rsize, *rcone, *rornt;
1778: PetscInt dim, cNew, Nct, n, r;
1780: DMPlexGetCellType(dm, c, &ct);
1781: dim = DMPolytopeTypeGetDim(ct);
1782: DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1783: /* This allows for different cell types */
1784: for (n = 0; n < Nct; ++n) {
1785: if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
1786: for (r = 0; r < rsize[n]; ++r) {
1787: PetscInt *closure = NULL;
1788: PetscInt clSize, cl, Nv = 0;
1790: DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew);
1791: DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
1792: for (cl = 0; cl < clSize*2; cl += 2) {if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;}
1793: DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
1794: PetscSectionSetDof(coordSectionNew, cNew, Nv * dE);
1795: PetscSectionSetFieldDof(coordSectionNew, cNew, 0, Nv * dE);
1796: }
1797: }
1798: }
1799: }
1800: }
1801: PetscSectionSetUp(coordSectionNew);
1802: DMViewFromOptions(dm, NULL, "-coarse_dm_view");
1803: DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew);
1804: {
1805: VecType vtype;
1806: PetscInt coordSizeNew, bs;
1807: const char *name;
1809: DMGetCoordinatesLocal(dm, &coordsLocal);
1810: VecCreate(PETSC_COMM_SELF, &coordsLocalNew);
1811: PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);
1812: VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE);
1813: PetscObjectGetName((PetscObject) coordsLocal, &name);
1814: PetscObjectSetName((PetscObject) coordsLocalNew, name);
1815: VecGetBlockSize(coordsLocal, &bs);
1816: VecSetBlockSize(coordsLocalNew, bs);
1817: VecGetType(coordsLocal, &vtype);
1818: VecSetType(coordsLocalNew, vtype);
1819: }
1820: VecGetArrayRead(coordsLocal, &coords);
1821: VecGetArray(coordsLocalNew, &coordsNew);
1822: PetscSectionGetChart(coordSection, &ocStart, &ocEnd);
1823: DMPlexGetChart(dm, &pStart, &pEnd);
1824: /* First set coordinates for vertices*/
1825: for (p = pStart; p < pEnd; ++p) {
1826: DMPolytopeType ct;
1827: DMPolytopeType *rct;
1828: PetscInt *rsize, *rcone, *rornt;
1829: PetscInt Nct, n, r;
1830: PetscBool hasVertex = PETSC_FALSE, isLocalized = PETSC_FALSE;
1832: DMPlexGetCellType(dm, p, &ct);
1833: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1834: for (n = 0; n < Nct; ++n) {
1835: if (rct[n] == DM_POLYTOPE_POINT) {hasVertex = PETSC_TRUE; break;}
1836: }
1837: if (localizeVertices && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) {
1838: PetscInt dof;
1839: PetscSectionGetDof(coordSection, p, &dof);
1840: if (dof) isLocalized = PETSC_TRUE;
1841: }
1842: if (hasVertex) {
1843: const PetscScalar *icoords = NULL;
1844: PetscScalar *pcoords = NULL;
1845: PetscInt Nc, Nv, v, d;
1847: DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);
1849: icoords = pcoords;
1850: Nv = Nc/dE;
1851: if (ct != DM_POLYTOPE_POINT) {
1852: if (localizeVertices) {
1853: PetscScalar anchor[3];
1855: for (d = 0; d < dE; ++d) anchor[d] = pcoords[d];
1856: if (!isLocalized) {
1857: for (v = 0; v < Nv; ++v) {DMLocalizeCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], &pcoords[v*dE]);}
1858: } else {
1859: Nv = Nc/(2*dE);
1860: for (v = Nv; v < Nv*2; ++v) {DMLocalizeCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], &pcoords[v*dE]);}
1861: }
1862: }
1863: }
1864: for (n = 0; n < Nct; ++n) {
1865: if (rct[n] != DM_POLYTOPE_POINT) continue;
1866: for (r = 0; r < rsize[n]; ++r) {
1867: PetscScalar vcoords[3];
1868: PetscInt vNew, off;
1870: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew);
1871: PetscSectionGetOffset(coordSectionNew, vNew, &off);
1872: DMPlexTransformMapCoordinates(tr, ct, rct[n], r, Nv, dE, icoords, vcoords);
1873: DMPlexSnapToGeomModel(dm, p, vcoords, &coordsNew[off]);
1874: }
1875: }
1876: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);
1877: }
1878: }
1879: /* Then set coordinates for cells by localizing */
1880: for (p = pStart; p < pEnd; ++p) {
1881: DMPolytopeType ct;
1882: DMPolytopeType *rct;
1883: PetscInt *rsize, *rcone, *rornt;
1884: PetscInt Nct, n, r;
1885: PetscBool isLocalized = PETSC_FALSE;
1887: DMPlexGetCellType(dm, p, &ct);
1888: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1889: if (localizeCells && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) {
1890: PetscInt dof;
1891: PetscSectionGetDof(coordSection, p, &dof);
1892: if (dof) isLocalized = PETSC_TRUE;
1893: }
1894: if (isLocalized) {
1895: const PetscScalar *pcoords;
1897: DMPlexPointLocalRead(cdm, p, coords, &pcoords);
1898: for (n = 0; n < Nct; ++n) {
1899: const PetscInt Nr = rsize[n];
1901: if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
1902: for (r = 0; r < Nr; ++r) {
1903: PetscInt pNew, offNew;
1905: /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
1906: DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
1907: cell to the ones it produces. */
1908: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1909: PetscSectionGetOffset(coordSectionNew, pNew, &offNew);
1910: DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew]);
1911: }
1912: }
1913: }
1914: }
1915: VecRestoreArrayRead(coordsLocal, &coords);
1916: VecRestoreArray(coordsLocalNew, &coordsNew);
1917: DMSetCoordinatesLocal(rdm, coordsLocalNew);
1918: /* TODO Stefano has a final reduction if some hybrid coordinates cannot be found. (needcoords) Should not be needed. */
1919: VecDestroy(&coordsLocalNew);
1920: PetscSectionDestroy(&coordSectionNew);
1921: if (!localizeCells) {DMLocalizeCoordinates(rdm);}
1922: return(0);
1923: }
1925: PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm)
1926: {
1927: DM rdm;
1928: DMPlexInterpolatedFlag interp;
1929: PetscInt dim, embedDim;
1930: PetscErrorCode ierr;
1936: DMPlexTransformSetDM(tr, dm);
1938: DMCreate(PetscObjectComm((PetscObject)dm), &rdm);
1939: DMSetType(rdm, DMPLEX);
1940: DMGetDimension(dm, &dim);
1941: DMSetDimension(rdm, dim);
1942: DMGetCoordinateDim(dm, &embedDim);
1943: DMSetCoordinateDim(rdm, embedDim);
1944: /* Calculate number of new points of each depth */
1945: DMPlexIsInterpolated(dm, &interp);
1946: if (interp != DMPLEX_INTERPOLATED_FULL) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement");
1947: /* Step 1: Set chart */
1948: DMPlexSetChart(rdm, 0, tr->ctStartNew[tr->ctOrder[DM_NUM_POLYTOPES]]);
1949: /* Step 2: Set cone/support sizes (automatically stratifies) */
1950: DMPlexTransformSetConeSizes(tr, rdm);
1951: /* Step 3: Setup refined DM */
1952: DMSetUp(rdm);
1953: /* Step 4: Set cones and supports (automatically symmetrizes) */
1954: DMPlexTransformSetCones(tr, rdm);
1955: /* Step 5: Create pointSF */
1956: DMPlexTransformCreateSF(tr, rdm);
1957: /* Step 6: Create labels */
1958: DMPlexTransformCreateLabels(tr, rdm);
1959: /* Step 7: Set coordinates */
1960: DMPlexTransformSetCoordinates(tr, rdm);
1961: *tdm = rdm;
1962: return(0);
1963: }
1965: PetscErrorCode DMPlexTransformAdaptLabel(DM dm, DMLabel adaptLabel, DM *rdm)
1966: {
1967: DMPlexTransform tr;
1968: DM cdm, rcdm;
1969: const char *prefix;
1970: PetscErrorCode ierr;
1973: DMPlexTransformCreate(PetscObjectComm((PetscObject) dm), &tr);
1974: PetscObjectGetOptionsPrefix((PetscObject) dm, &prefix);
1975: PetscObjectSetOptionsPrefix((PetscObject) tr, prefix);
1976: DMPlexTransformSetDM(tr, dm);
1977: DMPlexTransformSetFromOptions(tr);
1978: DMPlexTransformSetActive(tr, adaptLabel);
1979: DMPlexTransformSetUp(tr);
1980: PetscObjectViewFromOptions((PetscObject) tr, NULL, "-dm_plex_transform_view");
1981: DMPlexTransformApply(tr, dm, rdm);
1982: DMCopyDisc(dm, *rdm);
1983: DMGetCoordinateDM(dm, &cdm);
1984: DMGetCoordinateDM(*rdm, &rcdm);
1985: DMCopyDisc(cdm, rcdm);
1986: DMPlexTransformCreateDiscLabels(tr, *rdm);
1987: DMCopyDisc(dm, *rdm);
1988: DMPlexTransformDestroy(&tr);
1989: return(0);
1990: }