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: static PetscErrorCode DMPlexCreateCellTypeOrder_Internal(PetscInt dim, PetscInt *ctOrder[], PetscInt *ctOrderInv[])
12: {
13: PetscInt *ctO, *ctOInv;
14: PetscInt c, d, off = 0;
16: PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctO, DM_NUM_POLYTOPES + 1, &ctOInv);
17: for (d = 3; d >= dim; --d) {
18: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
19: if (DMPolytopeTypeGetDim((DMPolytopeType)c) != d) continue;
20: ctO[off++] = c;
21: }
22: }
23: if (dim != 0) {
24: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
25: if (DMPolytopeTypeGetDim((DMPolytopeType)c) != 0) continue;
26: ctO[off++] = c;
27: }
28: }
29: for (d = dim - 1; d > 0; --d) {
30: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
31: if (DMPolytopeTypeGetDim((DMPolytopeType)c) != d) continue;
32: ctO[off++] = c;
33: }
34: }
35: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
36: if (DMPolytopeTypeGetDim((DMPolytopeType)c) >= 0) continue;
37: ctO[off++] = c;
38: }
40: for (c = 0; c <= DM_NUM_POLYTOPES; ++c) ctOInv[ctO[c]] = c;
41: *ctOrder = ctO;
42: *ctOrderInv = ctOInv;
43: return 0;
44: }
46: /*@C
47: DMPlexTransformRegister - Adds a new transform component implementation
49: Not Collective
51: Input Parameters:
52: + name - The name of a new user-defined creation routine
53: - create_func - The creation routine itself
55: Notes:
56: DMPlexTransformRegister() may be called multiple times to add several user-defined transforms
58: Sample usage:
59: .vb
60: DMPlexTransformRegister("my_transform", MyTransformCreate);
61: .ve
63: Then, your transform type can be chosen with the procedural interface via
64: .vb
65: DMPlexTransformCreate(MPI_Comm, DMPlexTransform *);
66: DMPlexTransformSetType(DMPlexTransform, "my_transform");
67: .ve
68: or at runtime via the option
69: .vb
70: -dm_plex_transform_type my_transform
71: .ve
73: Level: advanced
75: .seealso: `DMPlexTransformRegisterAll()`, `DMPlexTransformRegisterDestroy()`
76: @*/
77: PetscErrorCode DMPlexTransformRegister(const char name[], PetscErrorCode (*create_func)(DMPlexTransform))
78: {
79: DMInitializePackage();
80: PetscFunctionListAdd(&DMPlexTransformList, name, create_func);
81: return 0;
82: }
84: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Filter(DMPlexTransform);
85: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform);
86: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToBox(DMPlexTransform);
87: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform);
88: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform);
89: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform);
90: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_1D(DMPlexTransform);
91: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform);
93: /*@C
94: DMPlexTransformRegisterAll - Registers all of the transform components in the DM package.
96: Not Collective
98: Level: advanced
100: .seealso: `DMRegisterAll()`, `DMPlexTransformRegisterDestroy()`
101: @*/
102: PetscErrorCode DMPlexTransformRegisterAll(void)
103: {
104: if (DMPlexTransformRegisterAllCalled) return 0;
105: DMPlexTransformRegisterAllCalled = PETSC_TRUE;
107: DMPlexTransformRegister(DMPLEXTRANSFORMFILTER, DMPlexTransformCreate_Filter);
108: DMPlexTransformRegister(DMPLEXREFINEREGULAR, DMPlexTransformCreate_Regular);
109: DMPlexTransformRegister(DMPLEXREFINETOBOX, DMPlexTransformCreate_ToBox);
110: DMPlexTransformRegister(DMPLEXREFINEALFELD, DMPlexTransformCreate_Alfeld);
111: DMPlexTransformRegister(DMPLEXREFINEBOUNDARYLAYER, DMPlexTransformCreate_BL);
112: DMPlexTransformRegister(DMPLEXREFINESBR, DMPlexTransformCreate_SBR);
113: DMPlexTransformRegister(DMPLEXREFINE1D, DMPlexTransformCreate_1D);
114: DMPlexTransformRegister(DMPLEXEXTRUDE, DMPlexTransformCreate_Extrude);
115: return 0;
116: }
118: /*@C
119: DMPlexTransformRegisterDestroy - This function destroys the . It is called from PetscFinalize().
121: Level: developer
123: .seealso: `PetscInitialize()`
124: @*/
125: PetscErrorCode DMPlexTransformRegisterDestroy(void)
126: {
127: PetscFunctionListDestroy(&DMPlexTransformList);
128: DMPlexTransformRegisterAllCalled = PETSC_FALSE;
129: return 0;
130: }
132: /*@
133: DMPlexTransformCreate - Creates an empty transform object. The type can then be set with DMPlexTransformSetType().
135: Collective
137: Input Parameter:
138: . comm - The communicator for the transform object
140: Output Parameter:
141: . dm - The transform object
143: Level: beginner
145: .seealso: `DMPlexTransformSetType()`, `DMPLEXREFINEREGULAR`, `DMPLEXTRANSFORMFILTER`
146: @*/
147: PetscErrorCode DMPlexTransformCreate(MPI_Comm comm, DMPlexTransform *tr)
148: {
149: DMPlexTransform t;
152: *tr = NULL;
153: DMInitializePackage();
155: PetscHeaderCreate(t, DMPLEXTRANSFORM_CLASSID, "DMPlexTransform", "Mesh Transform", "DMPlexTransform", comm, DMPlexTransformDestroy, DMPlexTransformView);
156: t->setupcalled = PETSC_FALSE;
157: PetscCalloc2(DM_NUM_POLYTOPES, &t->coordFE, DM_NUM_POLYTOPES, &t->refGeom);
158: *tr = t;
159: return 0;
160: }
162: /*@C
163: DMSetType - Sets the particular implementation for a transform.
165: Collective on tr
167: Input Parameters:
168: + tr - The transform
169: - method - The name of the transform type
171: Options Database Key:
172: . -dm_plex_transform_type <type> - Sets the transform type; use -help for a list of available types
174: Notes:
175: See "petsc/include/petscdmplextransform.h" for available transform types
177: Level: intermediate
179: .seealso: `DMPlexTransformGetType()`, `DMPlexTransformCreate()`
180: @*/
181: PetscErrorCode DMPlexTransformSetType(DMPlexTransform tr, DMPlexTransformType method)
182: {
183: PetscErrorCode (*r)(DMPlexTransform);
184: PetscBool match;
187: PetscObjectTypeCompare((PetscObject)tr, method, &match);
188: if (match) return 0;
190: DMPlexTransformRegisterAll();
191: PetscFunctionListFind(DMPlexTransformList, method, &r);
194: PetscTryTypeMethod(tr, destroy);
195: PetscMemzero(tr->ops, sizeof(*tr->ops));
196: PetscObjectChangeTypeName((PetscObject)tr, method);
197: (*r)(tr);
198: return 0;
199: }
201: /*@C
202: DMPlexTransformGetType - Gets the type name (as a string) from the transform.
204: Not Collective
206: Input Parameter:
207: . tr - The DMPlexTransform
209: Output Parameter:
210: . type - The DMPlexTransform type name
212: Level: intermediate
214: .seealso: `DMPlexTransformSetType()`, `DMPlexTransformCreate()`
215: @*/
216: PetscErrorCode DMPlexTransformGetType(DMPlexTransform tr, DMPlexTransformType *type)
217: {
220: DMPlexTransformRegisterAll();
221: *type = ((PetscObject)tr)->type_name;
222: return 0;
223: }
225: static PetscErrorCode DMPlexTransformView_Ascii(DMPlexTransform tr, PetscViewer v)
226: {
227: PetscViewerFormat format;
229: PetscViewerGetFormat(v, &format);
230: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
231: const PetscInt *trTypes = NULL;
232: IS trIS;
233: PetscInt cols = 8;
234: PetscInt Nrt = 8, f, g;
236: PetscViewerASCIIPrintf(v, "Source Starts\n");
237: for (g = 0; g <= cols; ++g) PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]);
238: PetscViewerASCIIPrintf(v, "\n");
239: for (f = 0; f <= cols; ++f) PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStart[f]);
240: PetscViewerASCIIPrintf(v, "\n");
241: PetscViewerASCIIPrintf(v, "Target Starts\n");
242: for (g = 0; g <= cols; ++g) PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]);
243: PetscViewerASCIIPrintf(v, "\n");
244: for (f = 0; f <= cols; ++f) PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStartNew[f]);
245: PetscViewerASCIIPrintf(v, "\n");
247: if (tr->trType) {
248: DMLabelGetNumValues(tr->trType, &Nrt);
249: DMLabelGetValueIS(tr->trType, &trIS);
250: ISGetIndices(trIS, &trTypes);
251: }
252: PetscViewerASCIIPrintf(v, "Offsets\n");
253: PetscViewerASCIIPrintf(v, " ");
254: for (g = 0; g < cols; ++g) PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g]);
255: PetscViewerASCIIPrintf(v, "\n");
256: for (f = 0; f < Nrt; ++f) {
257: PetscViewerASCIIPrintf(v, "%2" PetscInt_FMT " |", trTypes ? trTypes[f] : f);
258: for (g = 0; g < cols; ++g) PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->offset[f * DM_NUM_POLYTOPES + g]);
259: PetscViewerASCIIPrintf(v, " |\n");
260: }
261: if (trTypes) {
262: ISGetIndices(trIS, &trTypes);
263: ISDestroy(&trIS);
264: }
265: }
266: return 0;
267: }
269: /*@C
270: DMPlexTransformView - Views a DMPlexTransform
272: Collective on tr
274: Input Parameters:
275: + tr - the DMPlexTransform object to view
276: - v - the viewer
278: Level: beginner
280: .seealso `DMPlexTransformDestroy()`, `DMPlexTransformCreate()`
281: @*/
282: PetscErrorCode DMPlexTransformView(DMPlexTransform tr, PetscViewer v)
283: {
284: PetscBool isascii;
287: if (!v) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tr), &v);
290: PetscViewerCheckWritable(v);
291: PetscObjectPrintClassNamePrefixType((PetscObject)tr, v);
292: PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii);
293: if (isascii) DMPlexTransformView_Ascii(tr, v);
294: PetscTryTypeMethod(tr, view, v);
295: return 0;
296: }
298: /*@
299: DMPlexTransformSetFromOptions - Sets parameters in a transform from the options database
301: Collective on tr
303: Input Parameter:
304: . tr - the DMPlexTransform object to set options for
306: Options Database:
307: . -dm_plex_transform_type - Set the transform type, e.g. refine_regular
309: Level: intermediate
311: .seealso `DMPlexTransformView()`, `DMPlexTransformCreate()`
312: @*/
313: PetscErrorCode DMPlexTransformSetFromOptions(DMPlexTransform tr)
314: {
315: char typeName[1024];
316: const char *defName = DMPLEXREFINEREGULAR;
317: PetscBool flg;
320: PetscObjectOptionsBegin((PetscObject)tr);
321: PetscOptionsFList("-dm_plex_transform_type", "DMPlexTransform", "DMPlexTransformSetType", DMPlexTransformList, defName, typeName, 1024, &flg);
322: if (flg) DMPlexTransformSetType(tr, typeName);
323: else if (!((PetscObject)tr)->type_name) DMPlexTransformSetType(tr, defName);
324: PetscTryTypeMethod(tr, setfromoptions, PetscOptionsObject);
325: /* process any options handlers added with PetscObjectAddOptionsHandler() */
326: PetscObjectProcessOptionsHandlers((PetscObject)tr, PetscOptionsObject);
327: PetscOptionsEnd();
328: return 0;
329: }
331: /*@C
332: DMPlexTransformDestroy - Destroys a transform.
334: Collective on tr
336: Input Parameter:
337: . tr - the transform object to destroy
339: Level: beginner
341: .seealso `DMPlexTransformView()`, `DMPlexTransformCreate()`
342: @*/
343: PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr)
344: {
345: PetscInt c;
347: if (!*tr) return 0;
349: if (--((PetscObject)(*tr))->refct > 0) {
350: *tr = NULL;
351: return 0;
352: }
354: if ((*tr)->ops->destroy) (*(*tr)->ops->destroy)(*tr);
355: DMDestroy(&(*tr)->dm);
356: DMLabelDestroy(&(*tr)->active);
357: DMLabelDestroy(&(*tr)->trType);
358: PetscFree2((*tr)->ctOrderOld, (*tr)->ctOrderInvOld);
359: PetscFree2((*tr)->ctOrderNew, (*tr)->ctOrderInvNew);
360: PetscFree2((*tr)->ctStart, (*tr)->ctStartNew);
361: PetscFree((*tr)->offset);
362: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
363: PetscFEDestroy(&(*tr)->coordFE[c]);
364: PetscFEGeomDestroy(&(*tr)->refGeom[c]);
365: }
366: if ((*tr)->trVerts) {
367: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
368: DMPolytopeType *rct;
369: PetscInt *rsize, *rcone, *rornt, Nct, n, r;
371: if (DMPolytopeTypeGetDim((DMPolytopeType)c) > 0) {
372: DMPlexTransformCellTransform((*tr), (DMPolytopeType)c, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
373: for (n = 0; n < Nct; ++n) {
374: if (rct[n] == DM_POLYTOPE_POINT) continue;
375: for (r = 0; r < rsize[n]; ++r) PetscFree((*tr)->trSubVerts[c][rct[n]][r]);
376: PetscFree((*tr)->trSubVerts[c][rct[n]]);
377: }
378: }
379: PetscFree((*tr)->trSubVerts[c]);
380: PetscFree((*tr)->trVerts[c]);
381: }
382: }
383: PetscFree3((*tr)->trNv, (*tr)->trVerts, (*tr)->trSubVerts);
384: PetscFree2((*tr)->coordFE, (*tr)->refGeom);
385: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
386: PetscHeaderDestroy(tr);
387: return 0;
388: }
390: static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrderOld[], PetscInt ctStart[], PetscInt **offset)
391: {
392: DMLabel trType = tr->trType;
393: PetscInt c, cN, *off;
395: if (trType) {
396: DM dm;
397: IS rtIS;
398: const PetscInt *reftypes;
399: PetscInt Nrt, r;
401: DMPlexTransformGetDM(tr, &dm);
402: DMLabelGetNumValues(trType, &Nrt);
403: DMLabelGetValueIS(trType, &rtIS);
404: ISGetIndices(rtIS, &reftypes);
405: PetscCalloc1(Nrt * DM_NUM_POLYTOPES, &off);
406: for (r = 0; r < Nrt; ++r) {
407: const PetscInt rt = reftypes[r];
408: IS rtIS;
409: const PetscInt *points;
410: DMPolytopeType ct;
411: PetscInt p;
413: DMLabelGetStratumIS(trType, rt, &rtIS);
414: ISGetIndices(rtIS, &points);
415: p = points[0];
416: ISRestoreIndices(rtIS, &points);
417: ISDestroy(&rtIS);
418: DMPlexGetCellType(dm, p, &ct);
419: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
420: const DMPolytopeType ctNew = (DMPolytopeType)cN;
421: DMPolytopeType *rct;
422: PetscInt *rsize, *cone, *ornt;
423: PetscInt Nct, n, s;
425: if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {
426: off[r * DM_NUM_POLYTOPES + ctNew] = -1;
427: break;
428: }
429: off[r * DM_NUM_POLYTOPES + ctNew] = 0;
430: for (s = 0; s <= r; ++s) {
431: const PetscInt st = reftypes[s];
432: DMPolytopeType sct;
433: PetscInt q, qrt;
435: DMLabelGetStratumIS(trType, st, &rtIS);
436: ISGetIndices(rtIS, &points);
437: q = points[0];
438: ISRestoreIndices(rtIS, &points);
439: ISDestroy(&rtIS);
440: DMPlexGetCellType(dm, q, &sct);
441: DMPlexTransformCellTransform(tr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt);
443: if (st == rt) {
444: for (n = 0; n < Nct; ++n)
445: if (rct[n] == ctNew) break;
446: if (n == Nct) off[r * DM_NUM_POLYTOPES + ctNew] = -1;
447: break;
448: }
449: for (n = 0; n < Nct; ++n) {
450: if (rct[n] == ctNew) {
451: PetscInt sn;
453: DMLabelGetStratumSize(trType, st, &sn);
454: off[r * DM_NUM_POLYTOPES + ctNew] += sn * rsize[n];
455: }
456: }
457: }
458: }
459: }
460: ISRestoreIndices(rtIS, &reftypes);
461: ISDestroy(&rtIS);
462: } else {
463: PetscCalloc1(DM_NUM_POLYTOPES * DM_NUM_POLYTOPES, &off);
464: for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
465: const DMPolytopeType ct = (DMPolytopeType)c;
466: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
467: const DMPolytopeType ctNew = (DMPolytopeType)cN;
468: DMPolytopeType *rct;
469: PetscInt *rsize, *cone, *ornt;
470: PetscInt Nct, n, i;
472: if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {
473: off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
474: continue;
475: }
476: off[ct * DM_NUM_POLYTOPES + ctNew] = 0;
477: for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
478: const DMPolytopeType ict = (DMPolytopeType)ctOrderOld[i];
479: const DMPolytopeType ictn = (DMPolytopeType)ctOrderOld[i + 1];
481: DMPlexTransformCellTransform(tr, ict, PETSC_DETERMINE, NULL, &Nct, &rct, &rsize, &cone, &ornt);
482: if (ict == ct) {
483: for (n = 0; n < Nct; ++n)
484: if (rct[n] == ctNew) break;
485: if (n == Nct) off[ct * DM_NUM_POLYTOPES + ctNew] = -1;
486: break;
487: }
488: for (n = 0; n < Nct; ++n)
489: if (rct[n] == ctNew) off[ct * DM_NUM_POLYTOPES + ctNew] += (ctStart[ictn] - ctStart[ict]) * rsize[n];
490: }
491: }
492: }
493: }
494: *offset = off;
495: return 0;
496: }
498: PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr)
499: {
500: DM dm;
501: DMPolytopeType ctCell;
502: PetscInt pStart, pEnd, p, c, celldim = 0;
505: if (tr->setupcalled) return 0;
506: PetscTryTypeMethod(tr, setup);
507: DMPlexTransformGetDM(tr, &dm);
508: DMPlexGetChart(dm, &pStart, &pEnd);
509: if (pEnd > pStart) {
510: DMPlexGetCellType(dm, 0, &ctCell);
511: } else {
512: PetscInt dim;
514: DMGetDimension(dm, &dim);
515: switch (dim) {
516: case 0:
517: ctCell = DM_POLYTOPE_POINT;
518: break;
519: case 1:
520: ctCell = DM_POLYTOPE_SEGMENT;
521: break;
522: case 2:
523: ctCell = DM_POLYTOPE_TRIANGLE;
524: break;
525: case 3:
526: ctCell = DM_POLYTOPE_TETRAHEDRON;
527: break;
528: default:
529: ctCell = DM_POLYTOPE_UNKNOWN;
530: }
531: }
532: DMPlexCreateCellTypeOrder_Internal(DMPolytopeTypeGetDim(ctCell), &tr->ctOrderOld, &tr->ctOrderInvOld);
533: for (p = pStart; p < pEnd; ++p) {
534: DMPolytopeType ct;
535: DMPolytopeType *rct;
536: PetscInt *rsize, *cone, *ornt;
537: PetscInt Nct, n;
539: DMPlexGetCellType(dm, p, &ct);
541: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt);
542: for (n = 0; n < Nct; ++n) celldim = PetscMax(celldim, DMPolytopeTypeGetDim(rct[n]));
543: }
544: DMPlexCreateCellTypeOrder_Internal(celldim, &tr->ctOrderNew, &tr->ctOrderInvNew);
545: /* Construct sizes and offsets for each cell type */
546: if (!tr->ctStart) {
547: PetscInt *ctS, *ctSN, *ctC, *ctCN;
549: PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctS, DM_NUM_POLYTOPES + 1, &ctSN);
550: PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctC, DM_NUM_POLYTOPES + 1, &ctCN);
551: for (p = pStart; p < pEnd; ++p) {
552: DMPolytopeType ct;
553: DMPolytopeType *rct;
554: PetscInt *rsize, *cone, *ornt;
555: PetscInt Nct, n;
557: DMPlexGetCellType(dm, p, &ct);
559: ++ctC[ct];
560: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt);
561: for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
562: }
563: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
564: const PetscInt cto = tr->ctOrderOld[c];
565: const PetscInt cton = tr->ctOrderOld[c + 1];
566: const PetscInt ctn = tr->ctOrderNew[c];
567: const PetscInt ctnn = tr->ctOrderNew[c + 1];
569: ctS[cton] = ctS[cto] + ctC[cto];
570: ctSN[ctnn] = ctSN[ctn] + ctCN[ctn];
571: }
572: PetscFree2(ctC, ctCN);
573: tr->ctStart = ctS;
574: tr->ctStartNew = ctSN;
575: }
576: DMPlexTransformCreateOffset_Internal(tr, tr->ctOrderOld, tr->ctStart, &tr->offset);
577: tr->setupcalled = PETSC_TRUE;
578: return 0;
579: }
581: PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm)
582: {
585: *dm = tr->dm;
586: return 0;
587: }
589: PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm)
590: {
593: PetscObjectReference((PetscObject)dm);
594: DMDestroy(&tr->dm);
595: tr->dm = dm;
596: return 0;
597: }
599: PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active)
600: {
603: *active = tr->active;
604: return 0;
605: }
607: PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active)
608: {
611: PetscObjectReference((PetscObject)active);
612: DMLabelDestroy(&tr->active);
613: tr->active = active;
614: return 0;
615: }
617: static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe)
618: {
619: if (!tr->coordFE[ct]) {
620: PetscInt dim, cdim;
622: dim = DMPolytopeTypeGetDim(ct);
623: DMGetCoordinateDim(tr->dm, &cdim);
624: PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, cdim, ct, 1, PETSC_DETERMINE, &tr->coordFE[ct]);
625: {
626: PetscDualSpace dsp;
627: PetscQuadrature quad;
628: DM K;
629: PetscFEGeom *cg;
630: PetscScalar *Xq;
631: PetscReal *xq, *wq;
632: PetscInt Nq, q;
634: DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq);
635: PetscMalloc1(Nq * cdim, &xq);
636: for (q = 0; q < Nq * cdim; ++q) xq[q] = PetscRealPart(Xq[q]);
637: PetscMalloc1(Nq, &wq);
638: for (q = 0; q < Nq; ++q) wq[q] = 1.0;
639: PetscQuadratureCreate(PETSC_COMM_SELF, &quad);
640: PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq);
641: PetscFESetQuadrature(tr->coordFE[ct], quad);
643: PetscFEGetDualSpace(tr->coordFE[ct], &dsp);
644: PetscDualSpaceGetDM(dsp, &K);
645: PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &tr->refGeom[ct]);
646: cg = tr->refGeom[ct];
647: DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
648: PetscQuadratureDestroy(&quad);
649: }
650: }
651: *fe = tr->coordFE[ct];
652: return 0;
653: }
655: /*@
656: DMPlexTransformSetDimensions - Set the dimensions for the transformed DM
658: Input Parameters:
659: + tr - The DMPlexTransform object
660: - dm - The original DM
662: Output Parameter:
663: . tdm - The transformed DM
665: Level: advanced
667: .seealso: `DMPlexTransformApply()`, `DMPlexTransformCreate()`
668: @*/
669: PetscErrorCode DMPlexTransformSetDimensions(DMPlexTransform tr, DM dm, DM tdm)
670: {
671: if (tr->ops->setdimensions) PetscUseTypeMethod(tr, setdimensions, dm, tdm);
672: else {
673: PetscInt dim, cdim;
675: DMGetDimension(dm, &dim);
676: DMSetDimension(tdm, dim);
677: DMGetCoordinateDim(dm, &cdim);
678: DMSetCoordinateDim(tdm, cdim);
679: }
680: return 0;
681: }
683: /*@
684: DMPlexTransformGetTargetPoint - Get the number of a point in the transformed mesh based on information from the original mesh.
686: Not collective
688: Input Parameters:
689: + tr - The DMPlexTransform
690: . ct - The type of the original point which produces the new point
691: . ctNew - The type of the new point
692: . p - The original point which produces the new point
693: - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p
695: Output Parameters:
696: . pNew - The new point number
698: Level: developer
700: .seealso: `DMPlexTransformGetSourcePoint()`, `DMPlexTransformCellTransform()`
701: @*/
702: PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
703: {
704: DMPolytopeType *rct;
705: PetscInt *rsize, *cone, *ornt;
706: PetscInt rt, Nct, n, off, rp;
707: DMLabel trType = tr->trType;
708: PetscInt ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ct] + 1]];
709: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]];
710: PetscInt newp = ctSN, cind;
714: DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt);
715: if (trType) {
716: DMLabelGetValueIndex(trType, rt, &cind);
717: DMLabelGetStratumPointIndex(trType, rt, p, &rp);
719: } else {
720: cind = ct;
721: rp = p - ctS;
722: }
723: off = tr->offset[cind * DM_NUM_POLYTOPES + ctNew];
725: newp += off;
726: for (n = 0; n < Nct; ++n) {
727: if (rct[n] == ctNew) {
728: if (rsize[n] && r >= rsize[n])
729: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ") for subcell type %s in cell type %s", r, rsize[n], DMPolytopeTypes[rct[n]], DMPolytopeTypes[ct]);
730: newp += rp * rsize[n] + r;
731: break;
732: }
733: }
736: *pNew = newp;
737: return 0;
738: }
740: /*@
741: DMPlexTransformGetSourcePoint - Get the number of a point in the original mesh based on information from the transformed mesh.
743: Not collective
745: Input Parameters:
746: + tr - The DMPlexTransform
747: - pNew - The new point number
749: Output Parameters:
750: + ct - The type of the original point which produces the new point
751: . ctNew - The type of the new point
752: . p - The original point which produces the new point
753: - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p
755: Level: developer
757: .seealso: `DMPlexTransformGetTargetPoint()`, `DMPlexTransformCellTransform()`
758: @*/
759: PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r)
760: {
761: DMLabel trType = tr->trType;
762: DMPolytopeType *rct;
763: PetscInt *rsize, *cone, *ornt;
764: PetscInt rt, Nct, n, rp = 0, rO = 0, pO;
765: PetscInt offset = -1, ctS, ctE, ctO = 0, ctN, ctTmp;
767: for (ctN = 0; ctN < DM_NUM_POLYTOPES; ++ctN) {
768: PetscInt ctSN = tr->ctStartNew[ctN], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctN] + 1]];
770: if ((pNew >= ctSN) && (pNew < ctEN)) break;
771: }
773: if (trType) {
774: DM dm;
775: IS rtIS;
776: const PetscInt *reftypes;
777: PetscInt Nrt, r, rtStart;
779: DMPlexTransformGetDM(tr, &dm);
780: DMLabelGetNumValues(trType, &Nrt);
781: DMLabelGetValueIS(trType, &rtIS);
782: ISGetIndices(rtIS, &reftypes);
783: for (r = 0; r < Nrt; ++r) {
784: const PetscInt off = tr->offset[r * DM_NUM_POLYTOPES + ctN];
786: if (tr->ctStartNew[ctN] + off > pNew) continue;
787: /* Check that any of this refinement type exist */
788: /* TODO Actually keep track of the number produced here instead */
789: if (off > offset) {
790: rt = reftypes[r];
791: offset = off;
792: }
793: }
794: ISRestoreIndices(rtIS, &reftypes);
795: ISDestroy(&rtIS);
797: /* TODO Map refinement types to cell types */
798: DMLabelGetStratumBounds(trType, rt, &rtStart, NULL);
800: for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) {
801: PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
803: if ((rtStart >= ctS) && (rtStart < ctE)) break;
804: }
806: } else {
807: for (ctTmp = 0; ctTmp < DM_NUM_POLYTOPES; ++ctTmp) {
808: const PetscInt off = tr->offset[ctTmp * DM_NUM_POLYTOPES + ctN];
810: if (tr->ctStartNew[ctN] + off > pNew) continue;
811: if (tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctTmp] + 1]] <= tr->ctStart[ctTmp]) continue;
812: /* TODO Actually keep track of the number produced here instead */
813: if (off > offset) {
814: ctO = ctTmp;
815: offset = off;
816: }
817: }
819: }
820: ctS = tr->ctStart[ctO];
821: ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]];
822: DMPlexTransformCellTransform(tr, (DMPolytopeType)ctO, ctS, &rt, &Nct, &rct, &rsize, &cone, &ornt);
823: for (n = 0; n < Nct; ++n) {
824: if ((PetscInt)rct[n] == ctN) {
825: PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset;
827: rp = tmp / rsize[n];
828: rO = tmp % rsize[n];
829: break;
830: }
831: }
833: pO = rp + ctS;
835: if (ct) *ct = (DMPolytopeType)ctO;
836: if (ctNew) *ctNew = (DMPolytopeType)ctN;
837: if (p) *p = pO;
838: if (r) *r = rO;
839: return 0;
840: }
842: /*@
843: DMPlexTransformCellTransform - Describes the transform of a given source cell into a set of other target cells. These produced cells become the new mesh.
845: Input Parameters:
846: + tr - The DMPlexTransform object
847: . source - The source cell type
848: - p - The source point, which can also determine the refine type
850: Output Parameters:
851: + rt - The refine type for this point
852: . Nt - The number of types produced by this point
853: . target - An array of length Nt giving the types produced
854: . size - An array of length Nt giving the number of cells of each type produced
855: . cone - An array of length Nt*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point
856: - ornt - An array of length Nt*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point
858: Note: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
859: need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
860: division (described in these outputs) of the cell in the original mesh. The point identifier is given by
861: $ the number of cones to be taken, or 0 for the current cell
862: $ the cell cone point number at each level from which it is subdivided
863: $ the replica number r of the subdivision.
864: The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
865: $ Nt = 2
866: $ target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
867: $ size = {1, 2}
868: $ cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
869: $ ornt = { 0, 0, 0, 0}
871: Level: advanced
873: .seealso: `DMPlexTransformApply()`, `DMPlexTransformCreate()`
874: @*/
875: PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
876: {
877: PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt);
878: return 0;
879: }
881: PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
882: {
883: *rnew = r;
884: *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
885: return 0;
886: }
888: /* Returns the same thing */
889: PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
890: {
891: static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
892: static PetscInt vertexS[] = {1};
893: static PetscInt vertexC[] = {0};
894: static PetscInt vertexO[] = {0};
895: static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT};
896: static PetscInt edgeS[] = {1};
897: static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
898: static PetscInt edgeO[] = {0, 0};
899: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
900: static PetscInt tedgeS[] = {1};
901: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
902: static PetscInt tedgeO[] = {0, 0};
903: static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE};
904: static PetscInt triS[] = {1};
905: static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
906: static PetscInt triO[] = {0, 0, 0};
907: static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL};
908: static PetscInt quadS[] = {1};
909: 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};
910: static PetscInt quadO[] = {0, 0, 0, 0};
911: static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR};
912: static PetscInt tquadS[] = {1};
913: 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};
914: static PetscInt tquadO[] = {0, 0, 0, 0};
915: static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON};
916: static PetscInt tetS[] = {1};
917: 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};
918: static PetscInt tetO[] = {0, 0, 0, 0};
919: static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON};
920: static PetscInt hexS[] = {1};
921: static PetscInt hexC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 1, 5, 0};
922: static PetscInt hexO[] = {0, 0, 0, 0, 0, 0};
923: static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM};
924: static PetscInt tripS[] = {1};
925: static PetscInt tripC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0};
926: static PetscInt tripO[] = {0, 0, 0, 0, 0};
927: static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR};
928: static PetscInt ttripS[] = {1};
929: static PetscInt ttripC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0};
930: static PetscInt ttripO[] = {0, 0, 0, 0, 0};
931: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
932: static PetscInt tquadpS[] = {1};
933: static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0,
934: DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
935: static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0};
936: static DMPolytopeType pyrT[] = {DM_POLYTOPE_PYRAMID};
937: static PetscInt pyrS[] = {1};
938: static PetscInt pyrC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 1, 4, 0};
939: static PetscInt pyrO[] = {0, 0, 0, 0, 0};
941: if (rt) *rt = 0;
942: switch (source) {
943: case DM_POLYTOPE_POINT:
944: *Nt = 1;
945: *target = vertexT;
946: *size = vertexS;
947: *cone = vertexC;
948: *ornt = vertexO;
949: break;
950: case DM_POLYTOPE_SEGMENT:
951: *Nt = 1;
952: *target = edgeT;
953: *size = edgeS;
954: *cone = edgeC;
955: *ornt = edgeO;
956: break;
957: case DM_POLYTOPE_POINT_PRISM_TENSOR:
958: *Nt = 1;
959: *target = tedgeT;
960: *size = tedgeS;
961: *cone = tedgeC;
962: *ornt = tedgeO;
963: break;
964: case DM_POLYTOPE_TRIANGLE:
965: *Nt = 1;
966: *target = triT;
967: *size = triS;
968: *cone = triC;
969: *ornt = triO;
970: break;
971: case DM_POLYTOPE_QUADRILATERAL:
972: *Nt = 1;
973: *target = quadT;
974: *size = quadS;
975: *cone = quadC;
976: *ornt = quadO;
977: break;
978: case DM_POLYTOPE_SEG_PRISM_TENSOR:
979: *Nt = 1;
980: *target = tquadT;
981: *size = tquadS;
982: *cone = tquadC;
983: *ornt = tquadO;
984: break;
985: case DM_POLYTOPE_TETRAHEDRON:
986: *Nt = 1;
987: *target = tetT;
988: *size = tetS;
989: *cone = tetC;
990: *ornt = tetO;
991: break;
992: case DM_POLYTOPE_HEXAHEDRON:
993: *Nt = 1;
994: *target = hexT;
995: *size = hexS;
996: *cone = hexC;
997: *ornt = hexO;
998: break;
999: case DM_POLYTOPE_TRI_PRISM:
1000: *Nt = 1;
1001: *target = tripT;
1002: *size = tripS;
1003: *cone = tripC;
1004: *ornt = tripO;
1005: break;
1006: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1007: *Nt = 1;
1008: *target = ttripT;
1009: *size = ttripS;
1010: *cone = ttripC;
1011: *ornt = ttripO;
1012: break;
1013: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1014: *Nt = 1;
1015: *target = tquadpT;
1016: *size = tquadpS;
1017: *cone = tquadpC;
1018: *ornt = tquadpO;
1019: break;
1020: case DM_POLYTOPE_PYRAMID:
1021: *Nt = 1;
1022: *target = pyrT;
1023: *size = pyrS;
1024: *cone = pyrC;
1025: *ornt = pyrO;
1026: break;
1027: default:
1028: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1029: }
1030: return 0;
1031: }
1033: /*@
1034: DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point
1036: Not collective
1038: Input Parameters:
1039: + tr - The DMPlexTransform
1040: . sct - The source point cell type, from whom the new cell is being produced
1041: . sp - The source point
1042: . so - The orientation of the source point in its enclosing parent
1043: . tct - The target point cell type
1044: . r - The replica number requested for the produced cell type
1045: - o - The orientation of the replica
1047: Output Parameters:
1048: + rnew - The replica number, given the orientation of the parent
1049: - onew - The replica orientation, given the orientation of the parent
1051: Level: advanced
1053: .seealso: `DMPlexTransformCellTransform()`, `DMPlexTransformApply()`
1054: @*/
1055: PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1056: {
1058: PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew);
1059: return 0;
1060: }
1062: static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm)
1063: {
1064: DM dm;
1065: PetscInt pStart, pEnd, p, pNew;
1067: DMPlexTransformGetDM(tr, &dm);
1068: /* Must create the celltype label here so that we do not automatically try to compute the types */
1069: DMCreateLabel(rdm, "celltype");
1070: DMPlexGetChart(dm, &pStart, &pEnd);
1071: for (p = pStart; p < pEnd; ++p) {
1072: DMPolytopeType ct;
1073: DMPolytopeType *rct;
1074: PetscInt *rsize, *rcone, *rornt;
1075: PetscInt Nct, n, r;
1077: DMPlexGetCellType(dm, p, &ct);
1078: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1079: for (n = 0; n < Nct; ++n) {
1080: for (r = 0; r < rsize[n]; ++r) {
1081: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1082: DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]));
1083: DMPlexSetCellType(rdm, pNew, rct[n]);
1084: }
1085: }
1086: }
1087: /* Let the DM know we have set all the cell types */
1088: {
1089: DMLabel ctLabel;
1090: DM_Plex *plex = (DM_Plex *)rdm->data;
1092: DMPlexGetCellTypeLabel(rdm, &ctLabel);
1093: PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState);
1094: }
1095: return 0;
1096: }
1098: #if 0
1099: static PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize)
1100: {
1101: PetscInt ctNew;
1105: /* TODO Can do bisection since everything is sorted */
1106: for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) {
1107: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew]+1]];
1109: if (q >= ctSN && q < ctEN) break;
1110: }
1112: *coneSize = DMPolytopeTypeGetConeSize((DMPolytopeType) ctNew);
1113: return 0;
1114: }
1115: #endif
1117: /* The orientation o is for the interior of the cell p */
1118: static PetscErrorCode DMPlexTransformGetCone_Internal(DMPlexTransform tr, PetscInt p, PetscInt o, DMPolytopeType ct, DMPolytopeType ctNew, const PetscInt rcone[], PetscInt *coneoff, const PetscInt rornt[], PetscInt *orntoff, PetscInt coneNew[], PetscInt orntNew[])
1119: {
1120: DM dm;
1121: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1122: const PetscInt *cone;
1123: PetscInt c, coff = *coneoff, ooff = *orntoff;
1125: DMPlexTransformGetDM(tr, &dm);
1126: DMPlexGetCone(dm, p, &cone);
1127: for (c = 0; c < csizeNew; ++c) {
1128: PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */
1129: PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */
1130: PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */
1131: DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */
1132: const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */
1133: PetscInt pr = -1; /* Replica number of pp that produces new cone point */
1134: const DMPolytopeType ft = (DMPolytopeType)rcone[coff++]; /* Cell type for new cone point of pNew */
1135: const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1136: PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */
1137: PetscInt lc;
1139: /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
1140: for (lc = 0; lc < fn; ++lc) {
1141: const PetscInt *parr = DMPolytopeTypeGetArrangment(pct, po);
1142: const PetscInt acp = rcone[coff++];
1143: const PetscInt pcp = parr[acp * 2];
1144: const PetscInt pco = parr[acp * 2 + 1];
1145: const PetscInt *ppornt;
1147: ppp = pp;
1148: pp = pcone[pcp];
1149: DMPlexGetCellType(dm, pp, &pct);
1150: DMPlexGetCone(dm, pp, &pcone);
1151: DMPlexGetConeOrientation(dm, ppp, &ppornt);
1152: po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco);
1153: }
1154: pr = rcone[coff++];
1155: /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
1156: DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo);
1157: DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c]);
1158: orntNew[c] = fo;
1159: }
1160: *coneoff = coff;
1161: *orntoff = ooff;
1162: return 0;
1163: }
1165: static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm)
1166: {
1167: DM dm;
1168: DMPolytopeType ct;
1169: PetscInt *coneNew, *orntNew;
1170: PetscInt maxConeSize = 0, pStart, pEnd, p, pNew;
1172: DMPlexTransformGetDM(tr, &dm);
1173: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1174: DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew);
1175: DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew);
1176: DMPlexGetChart(dm, &pStart, &pEnd);
1177: for (p = pStart; p < pEnd; ++p) {
1178: const PetscInt *cone, *ornt;
1179: PetscInt coff, ooff;
1180: DMPolytopeType *rct;
1181: PetscInt *rsize, *rcone, *rornt;
1182: PetscInt Nct, n, r;
1184: DMPlexGetCellType(dm, p, &ct);
1185: DMPlexGetCone(dm, p, &cone);
1186: DMPlexGetConeOrientation(dm, p, &ornt);
1187: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1188: for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
1189: const DMPolytopeType ctNew = rct[n];
1191: for (r = 0; r < rsize[n]; ++r) {
1192: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1193: DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew);
1194: DMPlexSetCone(rdm, pNew, coneNew);
1195: DMPlexSetConeOrientation(rdm, pNew, orntNew);
1196: }
1197: }
1198: }
1199: DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew);
1200: DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew);
1201: DMViewFromOptions(rdm, NULL, "-rdm_view");
1202: DMPlexSymmetrize(rdm);
1203: DMPlexStratify(rdm);
1204: return 0;
1205: }
1207: PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[])
1208: {
1209: DM dm;
1210: DMPolytopeType ct, qct;
1211: DMPolytopeType *rct;
1212: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1213: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1218: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1219: DMPlexTransformGetDM(tr, &dm);
1220: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone);
1221: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt);
1222: DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r);
1223: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1224: for (n = 0; n < Nct; ++n) {
1225: const DMPolytopeType ctNew = rct[n];
1226: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1227: PetscInt Nr = rsize[n], fn, c;
1229: if (ctNew == qct) Nr = r;
1230: for (nr = 0; nr < Nr; ++nr) {
1231: for (c = 0; c < csizeNew; ++c) {
1232: ++coff; /* Cell type of new cone point */
1233: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1234: coff += fn;
1235: ++coff; /* Replica number of new cone point */
1236: ++ooff; /* Orientation of new cone point */
1237: }
1238: }
1239: if (ctNew == qct) break;
1240: }
1241: DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt);
1242: *cone = qcone;
1243: *ornt = qornt;
1244: return 0;
1245: }
1247: PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1248: {
1249: DM dm;
1250: DMPolytopeType ct, qct;
1251: DMPolytopeType *rct;
1252: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1253: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1257: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1258: DMPlexTransformGetDM(tr, &dm);
1259: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone);
1260: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt);
1261: DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r);
1262: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1263: for (n = 0; n < Nct; ++n) {
1264: const DMPolytopeType ctNew = rct[n];
1265: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1266: PetscInt Nr = rsize[n], fn, c;
1268: if (ctNew == qct) Nr = r;
1269: for (nr = 0; nr < Nr; ++nr) {
1270: for (c = 0; c < csizeNew; ++c) {
1271: ++coff; /* Cell type of new cone point */
1272: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1273: coff += fn;
1274: ++coff; /* Replica number of new cone point */
1275: ++ooff; /* Orientation of new cone point */
1276: }
1277: }
1278: if (ctNew == qct) break;
1279: }
1280: DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt);
1281: *cone = qcone;
1282: *ornt = qornt;
1283: return 0;
1284: }
1286: PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1287: {
1288: DM dm;
1292: DMPlexTransformGetDM(tr, &dm);
1293: DMRestoreWorkArray(dm, 0, MPIU_INT, cone);
1294: DMRestoreWorkArray(dm, 0, MPIU_INT, ornt);
1295: return 0;
1296: }
1298: static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr)
1299: {
1300: PetscInt ict;
1302: PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts);
1303: for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) {
1304: const DMPolytopeType ct = (DMPolytopeType)ict;
1305: DMPlexTransform reftr;
1306: DM refdm, trdm;
1307: Vec coordinates;
1308: const PetscScalar *coords;
1309: DMPolytopeType *rct;
1310: PetscInt *rsize, *rcone, *rornt;
1311: PetscInt Nct, n, r, pNew;
1312: PetscInt trdim, vStart, vEnd, Nc;
1313: const PetscInt debug = 0;
1314: const char *typeName;
1316: /* Since points are 0-dimensional, coordinates make no sense */
1317: if (DMPolytopeTypeGetDim(ct) <= 0) continue;
1318: DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm);
1319: DMPlexTransformCreate(PETSC_COMM_SELF, &reftr);
1320: DMPlexTransformSetDM(reftr, refdm);
1321: DMPlexTransformGetType(tr, &typeName);
1322: DMPlexTransformSetType(reftr, typeName);
1323: DMPlexTransformSetUp(reftr);
1324: DMPlexTransformApply(reftr, refdm, &trdm);
1326: DMGetDimension(trdm, &trdim);
1327: DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd);
1328: tr->trNv[ct] = vEnd - vStart;
1329: DMGetCoordinatesLocal(trdm, &coordinates);
1330: VecGetLocalSize(coordinates, &Nc);
1332: PetscCalloc1(Nc, &tr->trVerts[ct]);
1333: VecGetArrayRead(coordinates, &coords);
1334: PetscArraycpy(tr->trVerts[ct], coords, Nc);
1335: VecRestoreArrayRead(coordinates, &coords);
1337: PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct]);
1338: DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1339: for (n = 0; n < Nct; ++n) {
1340: /* Since points are 0-dimensional, coordinates make no sense */
1341: if (rct[n] == DM_POLYTOPE_POINT) continue;
1342: PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]]);
1343: for (r = 0; r < rsize[n]; ++r) {
1344: PetscInt *closure = NULL;
1345: PetscInt clSize, cl, Nv = 0;
1347: PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r]);
1348: DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew);
1349: DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure);
1350: for (cl = 0; cl < clSize * 2; cl += 2) {
1351: const PetscInt sv = closure[cl];
1353: if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart;
1354: }
1355: DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure);
1357: }
1358: }
1359: if (debug) {
1360: DMPolytopeType *rct;
1361: PetscInt *rsize, *rcone, *rornt;
1362: PetscInt v, dE = trdim, d, off = 0;
1364: PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct]);
1365: for (v = 0; v < tr->trNv[ct]; ++v) {
1366: PetscPrintf(PETSC_COMM_SELF, " ");
1367: for (d = 0; d < dE; ++d) PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++]));
1368: PetscPrintf(PETSC_COMM_SELF, "\n");
1369: }
1371: DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1372: for (n = 0; n < Nct; ++n) {
1373: if (rct[n] == DM_POLYTOPE_POINT) continue;
1374: PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct]);
1375: for (r = 0; r < rsize[n]; ++r) {
1376: PetscPrintf(PETSC_COMM_SELF, " ");
1377: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v]);
1378: PetscPrintf(PETSC_COMM_SELF, "\n");
1379: }
1380: }
1381: }
1382: DMDestroy(&refdm);
1383: DMDestroy(&trdm);
1384: DMPlexTransformDestroy(&reftr);
1385: }
1386: return 0;
1387: }
1389: /*
1390: DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type
1392: Input Parameters:
1393: + tr - The DMPlexTransform object
1394: - ct - The cell type
1396: Output Parameters:
1397: + Nv - The number of transformed vertices in the closure of the reference cell of given type
1398: - trVerts - The coordinates of these vertices in the reference cell
1400: Level: developer
1402: .seealso: `DMPlexTransformGetSubcellVertices()`
1403: */
1404: PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[])
1405: {
1406: if (!tr->trNv) DMPlexTransformCreateCellVertices_Internal(tr);
1407: if (Nv) *Nv = tr->trNv[ct];
1408: if (trVerts) *trVerts = tr->trVerts[ct];
1409: return 0;
1410: }
1412: /*
1413: DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type
1415: Input Parameters:
1416: + tr - The DMPlexTransform object
1417: . ct - The cell type
1418: . rct - The subcell type
1419: - r - The subcell index
1421: Output Parameter:
1422: . subVerts - The indices of these vertices in the set of vertices returned by DMPlexTransformGetCellVertices()
1424: Level: developer
1426: .seealso: `DMPlexTransformGetCellVertices()`
1427: */
1428: PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[])
1429: {
1430: if (!tr->trNv) DMPlexTransformCreateCellVertices_Internal(tr);
1432: if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r];
1433: return 0;
1434: }
1436: /* Computes new vertex as the barycenter, or centroid */
1437: PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1438: {
1439: PetscInt v, d;
1443: for (d = 0; d < dE; ++d) out[d] = 0.0;
1444: for (v = 0; v < Nv; ++v)
1445: for (d = 0; d < dE; ++d) out[d] += in[v * dE + d];
1446: for (d = 0; d < dE; ++d) out[d] /= Nv;
1447: return 0;
1448: }
1450: /*@
1451: DMPlexTransformMapCoordinates - Calculate new coordinates for produced points
1453: Not collective
1455: Input Parameters:
1456: + cr - The DMPlexCellRefiner
1457: . pct - The cell type of the parent, from whom the new cell is being produced
1458: . ct - The type being produced
1459: . p - The original point
1460: . r - The replica number requested for the produced cell type
1461: . Nv - Number of vertices in the closure of the parent cell
1462: . dE - Spatial dimension
1463: - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell
1465: Output Parameter:
1466: . out - The coordinates of the new vertices
1468: Level: intermediate
1470: .seealso: `DMPlexTransformApply()`
1471: @*/
1472: PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1473: {
1475: PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out);
1476: return 0;
1477: }
1479: static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew)
1480: {
1481: DM dm;
1482: IS valueIS;
1483: const PetscInt *values;
1484: PetscInt defVal, Nv, val;
1486: DMPlexTransformGetDM(tr, &dm);
1487: DMLabelGetDefaultValue(label, &defVal);
1488: DMLabelSetDefaultValue(labelNew, defVal);
1489: DMLabelGetValueIS(label, &valueIS);
1490: ISGetLocalSize(valueIS, &Nv);
1491: ISGetIndices(valueIS, &values);
1492: for (val = 0; val < Nv; ++val) {
1493: IS pointIS;
1494: const PetscInt *points;
1495: PetscInt numPoints, p;
1497: /* Ensure refined label is created with same number of strata as
1498: * original (even if no entries here). */
1499: DMLabelAddStratum(labelNew, values[val]);
1500: DMLabelGetStratumIS(label, values[val], &pointIS);
1501: ISGetLocalSize(pointIS, &numPoints);
1502: ISGetIndices(pointIS, &points);
1503: for (p = 0; p < numPoints; ++p) {
1504: const PetscInt point = points[p];
1505: DMPolytopeType ct;
1506: DMPolytopeType *rct;
1507: PetscInt *rsize, *rcone, *rornt;
1508: PetscInt Nct, n, r, pNew = 0;
1510: DMPlexGetCellType(dm, point, &ct);
1511: DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1512: for (n = 0; n < Nct; ++n) {
1513: for (r = 0; r < rsize[n]; ++r) {
1514: DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew);
1515: DMLabelSetValue(labelNew, pNew, values[val]);
1516: }
1517: }
1518: }
1519: ISRestoreIndices(pointIS, &points);
1520: ISDestroy(&pointIS);
1521: }
1522: ISRestoreIndices(valueIS, &values);
1523: ISDestroy(&valueIS);
1524: return 0;
1525: }
1527: static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm)
1528: {
1529: DM dm;
1530: PetscInt numLabels, l;
1532: DMPlexTransformGetDM(tr, &dm);
1533: DMGetNumLabels(dm, &numLabels);
1534: for (l = 0; l < numLabels; ++l) {
1535: DMLabel label, labelNew;
1536: const char *lname;
1537: PetscBool isDepth, isCellType;
1539: DMGetLabelName(dm, l, &lname);
1540: PetscStrcmp(lname, "depth", &isDepth);
1541: if (isDepth) continue;
1542: PetscStrcmp(lname, "celltype", &isCellType);
1543: if (isCellType) continue;
1544: DMCreateLabel(rdm, lname);
1545: DMGetLabel(dm, lname, &label);
1546: DMGetLabel(rdm, lname, &labelNew);
1547: RefineLabel_Internal(tr, label, labelNew);
1548: }
1549: return 0;
1550: }
1552: /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */
1553: PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm)
1554: {
1555: DM dm;
1556: PetscInt Nf, f, Nds, s;
1558: DMPlexTransformGetDM(tr, &dm);
1559: DMGetNumFields(dm, &Nf);
1560: for (f = 0; f < Nf; ++f) {
1561: DMLabel label, labelNew;
1562: PetscObject obj;
1563: const char *lname;
1565: DMGetField(rdm, f, &label, &obj);
1566: if (!label) continue;
1567: PetscObjectGetName((PetscObject)label, &lname);
1568: DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
1569: RefineLabel_Internal(tr, label, labelNew);
1570: DMSetField_Internal(rdm, f, labelNew, obj);
1571: DMLabelDestroy(&labelNew);
1572: }
1573: DMGetNumDS(dm, &Nds);
1574: for (s = 0; s < Nds; ++s) {
1575: DMLabel label, labelNew;
1576: const char *lname;
1578: DMGetRegionNumDS(rdm, s, &label, NULL, NULL);
1579: if (!label) continue;
1580: PetscObjectGetName((PetscObject)label, &lname);
1581: DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
1582: RefineLabel_Internal(tr, label, labelNew);
1583: DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL);
1584: DMLabelDestroy(&labelNew);
1585: }
1586: return 0;
1587: }
1589: static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm)
1590: {
1591: DM dm;
1592: PetscSF sf, sfNew;
1593: PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m;
1594: const PetscInt *localPoints;
1595: const PetscSFNode *remotePoints;
1596: PetscInt *localPointsNew;
1597: PetscSFNode *remotePointsNew;
1598: PetscInt pStartNew, pEndNew, pNew;
1599: /* Brute force algorithm */
1600: PetscSF rsf;
1601: PetscSection s;
1602: const PetscInt *rootdegree;
1603: PetscInt *rootPointsNew, *remoteOffsets;
1604: PetscInt numPointsNew, pStart, pEnd, p;
1606: DMPlexTransformGetDM(tr, &dm);
1607: DMPlexGetChart(rdm, &pStartNew, &pEndNew);
1608: DMGetPointSF(dm, &sf);
1609: DMGetPointSF(rdm, &sfNew);
1610: /* Calculate size of new SF */
1611: PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);
1612: if (numRoots < 0) return 0;
1613: for (l = 0; l < numLeaves; ++l) {
1614: const PetscInt p = localPoints[l];
1615: DMPolytopeType ct;
1616: DMPolytopeType *rct;
1617: PetscInt *rsize, *rcone, *rornt;
1618: PetscInt Nct, n;
1620: DMPlexGetCellType(dm, p, &ct);
1621: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1622: for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n];
1623: }
1624: /* Send new root point numbers
1625: It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication
1626: */
1627: DMPlexGetChart(dm, &pStart, &pEnd);
1628: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s);
1629: PetscSectionSetChart(s, pStart, pEnd);
1630: for (p = pStart; p < pEnd; ++p) {
1631: DMPolytopeType ct;
1632: DMPolytopeType *rct;
1633: PetscInt *rsize, *rcone, *rornt;
1634: PetscInt Nct, n;
1636: DMPlexGetCellType(dm, p, &ct);
1637: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1638: for (n = 0; n < Nct; ++n) PetscSectionAddDof(s, p, rsize[n]);
1639: }
1640: PetscSectionSetUp(s);
1641: PetscSectionGetStorageSize(s, &numPointsNew);
1642: PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets);
1643: PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf);
1644: PetscFree(remoteOffsets);
1645: PetscSFComputeDegreeBegin(sf, &rootdegree);
1646: PetscSFComputeDegreeEnd(sf, &rootdegree);
1647: PetscMalloc1(numPointsNew, &rootPointsNew);
1648: for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1;
1649: for (p = pStart; p < pEnd; ++p) {
1650: DMPolytopeType ct;
1651: DMPolytopeType *rct;
1652: PetscInt *rsize, *rcone, *rornt;
1653: PetscInt Nct, n, r, off;
1655: if (!rootdegree[p - pStart]) continue;
1656: PetscSectionGetOffset(s, p, &off);
1657: DMPlexGetCellType(dm, p, &ct);
1658: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1659: for (n = 0, m = 0; n < Nct; ++n) {
1660: for (r = 0; r < rsize[n]; ++r, ++m) {
1661: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1662: rootPointsNew[off + m] = pNew;
1663: }
1664: }
1665: }
1666: PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE);
1667: PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE);
1668: PetscSFDestroy(&rsf);
1669: PetscMalloc1(numLeavesNew, &localPointsNew);
1670: PetscMalloc1(numLeavesNew, &remotePointsNew);
1671: for (l = 0, m = 0; l < numLeaves; ++l) {
1672: const PetscInt p = localPoints[l];
1673: DMPolytopeType ct;
1674: DMPolytopeType *rct;
1675: PetscInt *rsize, *rcone, *rornt;
1676: PetscInt Nct, n, r, q, off;
1678: PetscSectionGetOffset(s, p, &off);
1679: DMPlexGetCellType(dm, p, &ct);
1680: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1681: for (n = 0, q = 0; n < Nct; ++n) {
1682: for (r = 0; r < rsize[n]; ++r, ++m, ++q) {
1683: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1684: localPointsNew[m] = pNew;
1685: remotePointsNew[m].index = rootPointsNew[off + q];
1686: remotePointsNew[m].rank = remotePoints[l].rank;
1687: }
1688: }
1689: }
1690: PetscSectionDestroy(&s);
1691: PetscFree(rootPointsNew);
1692: /* SF needs sorted leaves to correctly calculate Gather */
1693: {
1694: PetscSFNode *rp, *rtmp;
1695: PetscInt *lp, *idx, *ltmp, i;
1697: PetscMalloc1(numLeavesNew, &idx);
1698: PetscMalloc1(numLeavesNew, &lp);
1699: PetscMalloc1(numLeavesNew, &rp);
1700: for (i = 0; i < numLeavesNew; ++i) {
1702: idx[i] = i;
1703: }
1704: PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx);
1705: for (i = 0; i < numLeavesNew; ++i) {
1706: lp[i] = localPointsNew[idx[i]];
1707: rp[i] = remotePointsNew[idx[i]];
1708: }
1709: ltmp = localPointsNew;
1710: localPointsNew = lp;
1711: rtmp = remotePointsNew;
1712: remotePointsNew = rp;
1713: PetscFree(idx);
1714: PetscFree(ltmp);
1715: PetscFree(rtmp);
1716: }
1717: PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
1718: return 0;
1719: }
1721: /*
1722: DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of type ct with localized coordinates x, we generate localized coordinates xr for subcell r of type rct.
1724: Not collective
1726: Input Parameters:
1727: + cr - The DMPlexCellRefiner
1728: . ct - The type of the parent cell
1729: . rct - The type of the produced cell
1730: . r - The index of the produced cell
1731: - x - The localized coordinates for the parent cell
1733: Output Parameter:
1734: . xr - The localized coordinates for the produced cell
1736: Level: developer
1738: .seealso: `DMPlexCellRefinerSetCoordinates()`
1739: */
1740: static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
1741: {
1742: PetscFE fe = NULL;
1743: PetscInt cdim, v, *subcellV;
1745: DMPlexTransformGetCoordinateFE(tr, ct, &fe);
1746: DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV);
1747: PetscFEGetNumComponents(fe, &cdim);
1748: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim]);
1749: return 0;
1750: }
1752: static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm)
1753: {
1754: DM dm, cdm, cdmCell, cdmNew, cdmCellNew;
1755: PetscSection coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew;
1756: Vec coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew;
1757: const PetscScalar *coords;
1758: PetscScalar *coordsNew;
1759: const PetscReal *maxCell, *Lstart, *L;
1760: PetscBool localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE;
1761: PetscInt dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p;
1763: DMPlexTransformGetDM(tr, &dm);
1764: DMGetCoordinateDM(dm, &cdm);
1765: DMGetCellCoordinateDM(dm, &cdmCell);
1766: DMGetCoordinatesLocalized(dm, &localized);
1767: DMGetPeriodicity(dm, &maxCell, &Lstart, &L);
1768: if (localized) {
1769: /* Localize coordinates of new vertices */
1770: localizeVertices = PETSC_TRUE;
1771: /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */
1772: if (!maxCell) localizeCells = PETSC_TRUE;
1773: }
1774: DMGetCoordinateSection(dm, &coordSection);
1775: PetscSectionGetFieldComponents(coordSection, 0, &dEo);
1776: if (maxCell) {
1777: PetscReal maxCellNew[3];
1779: for (d = 0; d < dEo; ++d) maxCellNew[d] = maxCell[d] / 2.0;
1780: DMSetPeriodicity(rdm, maxCellNew, Lstart, L);
1781: }
1782: DMGetCoordinateDim(rdm, &dE);
1783: PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew);
1784: PetscSectionSetNumFields(coordSectionNew, 1);
1785: PetscSectionSetFieldComponents(coordSectionNew, 0, dE);
1786: DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew);
1787: PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew);
1788: /* Localization should be inherited */
1789: /* Stefano calculates parent cells for each new cell for localization */
1790: /* Localized cells need coordinates of closure */
1791: for (v = vStartNew; v < vEndNew; ++v) {
1792: PetscSectionSetDof(coordSectionNew, v, dE);
1793: PetscSectionSetFieldDof(coordSectionNew, v, 0, dE);
1794: }
1795: PetscSectionSetUp(coordSectionNew);
1796: DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew);
1798: if (localizeCells) {
1799: DMGetCoordinateDM(rdm, &cdmNew);
1800: DMClone(cdmNew, &cdmCellNew);
1801: DMSetCellCoordinateDM(rdm, cdmCellNew);
1802: DMDestroy(&cdmCellNew);
1804: PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew);
1805: PetscSectionSetNumFields(coordSectionCellNew, 1);
1806: PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE);
1807: DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew);
1808: PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew);
1810: DMGetCellCoordinateSection(dm, &coordSectionCell);
1811: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1812: for (c = cStart; c < cEnd; ++c) {
1813: PetscInt dof;
1815: PetscSectionGetDof(coordSectionCell, c, &dof);
1816: if (dof) {
1817: DMPolytopeType ct;
1818: DMPolytopeType *rct;
1819: PetscInt *rsize, *rcone, *rornt;
1820: PetscInt dim, cNew, Nct, n, r;
1822: DMPlexGetCellType(dm, c, &ct);
1823: dim = DMPolytopeTypeGetDim(ct);
1824: DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1825: /* This allows for different cell types */
1826: for (n = 0; n < Nct; ++n) {
1827: if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
1828: for (r = 0; r < rsize[n]; ++r) {
1829: PetscInt *closure = NULL;
1830: PetscInt clSize, cl, Nv = 0;
1832: DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew);
1833: DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
1834: for (cl = 0; cl < clSize * 2; cl += 2) {
1835: if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;
1836: }
1837: DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
1838: PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE);
1839: PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE);
1840: }
1841: }
1842: }
1843: }
1844: PetscSectionSetUp(coordSectionCellNew);
1845: DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew);
1846: }
1847: DMViewFromOptions(dm, NULL, "-coarse_dm_view");
1848: {
1849: VecType vtype;
1850: PetscInt coordSizeNew, bs;
1851: const char *name;
1853: DMGetCoordinatesLocal(dm, &coordsLocal);
1854: VecCreate(PETSC_COMM_SELF, &coordsLocalNew);
1855: PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);
1856: VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE);
1857: PetscObjectGetName((PetscObject)coordsLocal, &name);
1858: PetscObjectSetName((PetscObject)coordsLocalNew, name);
1859: VecGetBlockSize(coordsLocal, &bs);
1860: VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE);
1861: VecGetType(coordsLocal, &vtype);
1862: VecSetType(coordsLocalNew, vtype);
1863: }
1864: VecGetArrayRead(coordsLocal, &coords);
1865: VecGetArray(coordsLocalNew, &coordsNew);
1866: DMPlexGetChart(dm, &pStart, &pEnd);
1867: /* First set coordinates for vertices */
1868: for (p = pStart; p < pEnd; ++p) {
1869: DMPolytopeType ct;
1870: DMPolytopeType *rct;
1871: PetscInt *rsize, *rcone, *rornt;
1872: PetscInt Nct, n, r;
1873: PetscBool hasVertex = PETSC_FALSE;
1875: DMPlexGetCellType(dm, p, &ct);
1876: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1877: for (n = 0; n < Nct; ++n) {
1878: if (rct[n] == DM_POLYTOPE_POINT) {
1879: hasVertex = PETSC_TRUE;
1880: break;
1881: }
1882: }
1883: if (hasVertex) {
1884: const PetscScalar *icoords = NULL;
1885: const PetscScalar *array = NULL;
1886: PetscScalar *pcoords = NULL;
1887: PetscBool isDG;
1888: PetscInt Nc, Nv, v, d;
1890: DMPlexGetCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords);
1892: icoords = pcoords;
1893: Nv = Nc / dEo;
1894: if (ct != DM_POLYTOPE_POINT) {
1895: if (localizeVertices && maxCell) {
1896: PetscScalar anchor[3];
1898: for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d];
1899: for (v = 0; v < Nv; ++v) DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo]);
1900: }
1901: }
1902: for (n = 0; n < Nct; ++n) {
1903: if (rct[n] != DM_POLYTOPE_POINT) continue;
1904: for (r = 0; r < rsize[n]; ++r) {
1905: PetscScalar vcoords[3];
1906: PetscInt vNew, off;
1908: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew);
1909: PetscSectionGetOffset(coordSectionNew, vNew, &off);
1910: DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords);
1911: DMPlexSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off]);
1912: }
1913: }
1914: DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords);
1915: }
1916: }
1917: VecRestoreArrayRead(coordsLocal, &coords);
1918: VecRestoreArray(coordsLocalNew, &coordsNew);
1919: DMSetCoordinatesLocal(rdm, coordsLocalNew);
1920: VecDestroy(&coordsLocalNew);
1921: PetscSectionDestroy(&coordSectionNew);
1922: /* Then set coordinates for cells by localizing */
1923: if (!localizeCells) DMLocalizeCoordinates(rdm);
1924: else {
1925: VecType vtype;
1926: PetscInt coordSizeNew, bs;
1927: const char *name;
1929: DMGetCellCoordinatesLocal(dm, &coordsLocalCell);
1930: VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew);
1931: PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew);
1932: VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE);
1933: PetscObjectGetName((PetscObject)coordsLocalCell, &name);
1934: PetscObjectSetName((PetscObject)coordsLocalCellNew, name);
1935: VecGetBlockSize(coordsLocalCell, &bs);
1936: VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE);
1937: VecGetType(coordsLocalCell, &vtype);
1938: VecSetType(coordsLocalCellNew, vtype);
1939: VecGetArrayRead(coordsLocalCell, &coords);
1940: VecGetArray(coordsLocalCellNew, &coordsNew);
1942: for (p = pStart; p < pEnd; ++p) {
1943: DMPolytopeType ct;
1944: DMPolytopeType *rct;
1945: PetscInt *rsize, *rcone, *rornt;
1946: PetscInt dof = 0, Nct, n, r;
1948: DMPlexGetCellType(dm, p, &ct);
1949: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1950: if (p >= cStart && p < cEnd) PetscSectionGetDof(coordSectionCell, p, &dof);
1951: if (dof) {
1952: const PetscScalar *pcoords;
1954: DMPlexPointLocalRead(cdmCell, p, coords, &pcoords);
1955: for (n = 0; n < Nct; ++n) {
1956: const PetscInt Nr = rsize[n];
1958: if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
1959: for (r = 0; r < Nr; ++r) {
1960: PetscInt pNew, offNew;
1962: /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
1963: DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
1964: cell to the ones it produces. */
1965: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1966: PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew);
1967: DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew]);
1968: }
1969: }
1970: }
1971: }
1972: VecRestoreArrayRead(coordsLocalCell, &coords);
1973: VecRestoreArray(coordsLocalCellNew, &coordsNew);
1974: DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew);
1975: VecDestroy(&coordsLocalCellNew);
1976: PetscSectionDestroy(&coordSectionCellNew);
1977: }
1978: return 0;
1979: }
1981: PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm)
1982: {
1983: DM rdm;
1984: DMPlexInterpolatedFlag interp;
1989: DMPlexTransformSetDM(tr, dm);
1991: DMCreate(PetscObjectComm((PetscObject)dm), &rdm);
1992: DMSetType(rdm, DMPLEX);
1993: DMPlexTransformSetDimensions(tr, dm, rdm);
1994: /* Calculate number of new points of each depth */
1995: DMPlexIsInterpolatedCollective(dm, &interp);
1997: /* Step 1: Set chart */
1998: DMPlexSetChart(rdm, 0, tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]]);
1999: /* Step 2: Set cone/support sizes (automatically stratifies) */
2000: DMPlexTransformSetConeSizes(tr, rdm);
2001: /* Step 3: Setup refined DM */
2002: DMSetUp(rdm);
2003: /* Step 4: Set cones and supports (automatically symmetrizes) */
2004: DMPlexTransformSetCones(tr, rdm);
2005: /* Step 5: Create pointSF */
2006: DMPlexTransformCreateSF(tr, rdm);
2007: /* Step 6: Create labels */
2008: DMPlexTransformCreateLabels(tr, rdm);
2009: /* Step 7: Set coordinates */
2010: DMPlexTransformSetCoordinates(tr, rdm);
2011: DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, rdm);
2012: *tdm = rdm;
2013: return 0;
2014: }
2016: PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm)
2017: {
2018: DMPlexTransform tr;
2019: DM cdm, rcdm;
2020: const char *prefix;
2022: DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr);
2023: PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
2024: PetscObjectSetOptionsPrefix((PetscObject)tr, prefix);
2025: DMPlexTransformSetDM(tr, dm);
2026: DMPlexTransformSetFromOptions(tr);
2027: DMPlexTransformSetActive(tr, adaptLabel);
2028: DMPlexTransformSetUp(tr);
2029: PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view");
2030: DMPlexTransformApply(tr, dm, rdm);
2031: DMCopyDisc(dm, *rdm);
2032: DMGetCoordinateDM(dm, &cdm);
2033: DMGetCoordinateDM(*rdm, &rcdm);
2034: DMCopyDisc(cdm, rcdm);
2035: DMPlexTransformCreateDiscLabels(tr, *rdm);
2036: DMCopyDisc(dm, *rdm);
2037: DMPlexTransformDestroy(&tr);
2038: ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation;
2039: return 0;
2040: }