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: Sample usage:
56: .vb
57: DMPlexTransformRegister("my_transform", MyTransformCreate);
58: .ve
60: Then, your transform type can be chosen with the procedural interface via
61: .vb
62: DMPlexTransformCreate(MPI_Comm, DMPlexTransform *);
63: DMPlexTransformSetType(DMPlexTransform, "my_transform");
64: .ve
65: or at runtime via the option
66: .vb
67: -dm_plex_transform_type my_transform
68: .ve
70: Level: advanced
72: Note:
73: `DMPlexTransformRegister()` may be called multiple times to add several user-defined transforms
75: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `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: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `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 registered `DMPlexTransformType`. It is called from `PetscFinalize()`.
121: Level: developer
123: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMRegisterAll()`, `DMPlexTransformType`, `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: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `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: Level: intermediate
176: Notes:
177: See "petsc/include/petscdmplextransform.h" for available transform types
179: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `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 `DMPlexTransformType` name
212: Level: intermediate
214: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `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: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `PetscViewer`, `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 Key:
307: . -dm_plex_transform_type - Set the transform type, e.g. refine_regular
309: Level: intermediate
311: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `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 `DMPlexTransform`
334: Collective on tr
336: Input Parameter:
337: . tr - the transform object to destroy
339: Level: beginner
341: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `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: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `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: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `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: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `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: Level: advanced
860: Notes:
861: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
862: need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
863: division (described in these outputs) of the cell in the original mesh. The point identifier is given by
864: .vb
865: the number of cones to be taken, or 0 for the current cell
866: the cell cone point number at each level from which it is subdivided
867: the replica number r of the subdivision.
868: .ve
869: The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
870: .vb
871: Nt = 2
872: target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
873: size = {1, 2}
874: cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
875: ornt = { 0, 0, 0, 0}
876: .ve
878: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`, `DMPlexTransformCreate()`
879: @*/
880: PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
881: {
882: PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt);
883: return 0;
884: }
886: PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
887: {
888: *rnew = r;
889: *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
890: return 0;
891: }
893: /* Returns the same thing */
894: PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
895: {
896: static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
897: static PetscInt vertexS[] = {1};
898: static PetscInt vertexC[] = {0};
899: static PetscInt vertexO[] = {0};
900: static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT};
901: static PetscInt edgeS[] = {1};
902: static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
903: static PetscInt edgeO[] = {0, 0};
904: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
905: static PetscInt tedgeS[] = {1};
906: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
907: static PetscInt tedgeO[] = {0, 0};
908: static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE};
909: static PetscInt triS[] = {1};
910: static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
911: static PetscInt triO[] = {0, 0, 0};
912: static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL};
913: static PetscInt quadS[] = {1};
914: 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};
915: static PetscInt quadO[] = {0, 0, 0, 0};
916: static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR};
917: static PetscInt tquadS[] = {1};
918: 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};
919: static PetscInt tquadO[] = {0, 0, 0, 0};
920: static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON};
921: static PetscInt tetS[] = {1};
922: 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};
923: static PetscInt tetO[] = {0, 0, 0, 0};
924: static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON};
925: static PetscInt hexS[] = {1};
926: 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};
927: static PetscInt hexO[] = {0, 0, 0, 0, 0, 0};
928: static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM};
929: static PetscInt tripS[] = {1};
930: 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};
931: static PetscInt tripO[] = {0, 0, 0, 0, 0};
932: static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR};
933: static PetscInt ttripS[] = {1};
934: 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};
935: static PetscInt ttripO[] = {0, 0, 0, 0, 0};
936: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
937: static PetscInt tquadpS[] = {1};
938: static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0,
939: DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
940: static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0};
941: static DMPolytopeType pyrT[] = {DM_POLYTOPE_PYRAMID};
942: static PetscInt pyrS[] = {1};
943: 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};
944: static PetscInt pyrO[] = {0, 0, 0, 0, 0};
946: if (rt) *rt = 0;
947: switch (source) {
948: case DM_POLYTOPE_POINT:
949: *Nt = 1;
950: *target = vertexT;
951: *size = vertexS;
952: *cone = vertexC;
953: *ornt = vertexO;
954: break;
955: case DM_POLYTOPE_SEGMENT:
956: *Nt = 1;
957: *target = edgeT;
958: *size = edgeS;
959: *cone = edgeC;
960: *ornt = edgeO;
961: break;
962: case DM_POLYTOPE_POINT_PRISM_TENSOR:
963: *Nt = 1;
964: *target = tedgeT;
965: *size = tedgeS;
966: *cone = tedgeC;
967: *ornt = tedgeO;
968: break;
969: case DM_POLYTOPE_TRIANGLE:
970: *Nt = 1;
971: *target = triT;
972: *size = triS;
973: *cone = triC;
974: *ornt = triO;
975: break;
976: case DM_POLYTOPE_QUADRILATERAL:
977: *Nt = 1;
978: *target = quadT;
979: *size = quadS;
980: *cone = quadC;
981: *ornt = quadO;
982: break;
983: case DM_POLYTOPE_SEG_PRISM_TENSOR:
984: *Nt = 1;
985: *target = tquadT;
986: *size = tquadS;
987: *cone = tquadC;
988: *ornt = tquadO;
989: break;
990: case DM_POLYTOPE_TETRAHEDRON:
991: *Nt = 1;
992: *target = tetT;
993: *size = tetS;
994: *cone = tetC;
995: *ornt = tetO;
996: break;
997: case DM_POLYTOPE_HEXAHEDRON:
998: *Nt = 1;
999: *target = hexT;
1000: *size = hexS;
1001: *cone = hexC;
1002: *ornt = hexO;
1003: break;
1004: case DM_POLYTOPE_TRI_PRISM:
1005: *Nt = 1;
1006: *target = tripT;
1007: *size = tripS;
1008: *cone = tripC;
1009: *ornt = tripO;
1010: break;
1011: case DM_POLYTOPE_TRI_PRISM_TENSOR:
1012: *Nt = 1;
1013: *target = ttripT;
1014: *size = ttripS;
1015: *cone = ttripC;
1016: *ornt = ttripO;
1017: break;
1018: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1019: *Nt = 1;
1020: *target = tquadpT;
1021: *size = tquadpS;
1022: *cone = tquadpC;
1023: *ornt = tquadpO;
1024: break;
1025: case DM_POLYTOPE_PYRAMID:
1026: *Nt = 1;
1027: *target = pyrT;
1028: *size = pyrS;
1029: *cone = pyrC;
1030: *ornt = pyrO;
1031: break;
1032: default:
1033: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
1034: }
1035: return 0;
1036: }
1038: /*@
1039: DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point
1041: Not Collective
1043: Input Parameters:
1044: + tr - The `DMPlexTransform`
1045: . sct - The source point cell type, from whom the new cell is being produced
1046: . sp - The source point
1047: . so - The orientation of the source point in its enclosing parent
1048: . tct - The target point cell type
1049: . r - The replica number requested for the produced cell type
1050: - o - The orientation of the replica
1052: Output Parameters:
1053: + rnew - The replica number, given the orientation of the parent
1054: - onew - The replica orientation, given the orientation of the parent
1056: Level: advanced
1058: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformCellTransform()`, `DMPlexTransformApply()`
1059: @*/
1060: PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
1061: {
1063: PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew);
1064: return 0;
1065: }
1067: static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm)
1068: {
1069: DM dm;
1070: PetscInt pStart, pEnd, p, pNew;
1072: DMPlexTransformGetDM(tr, &dm);
1073: /* Must create the celltype label here so that we do not automatically try to compute the types */
1074: DMCreateLabel(rdm, "celltype");
1075: DMPlexGetChart(dm, &pStart, &pEnd);
1076: for (p = pStart; p < pEnd; ++p) {
1077: DMPolytopeType ct;
1078: DMPolytopeType *rct;
1079: PetscInt *rsize, *rcone, *rornt;
1080: PetscInt Nct, n, r;
1082: DMPlexGetCellType(dm, p, &ct);
1083: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1084: for (n = 0; n < Nct; ++n) {
1085: for (r = 0; r < rsize[n]; ++r) {
1086: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1087: DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]));
1088: DMPlexSetCellType(rdm, pNew, rct[n]);
1089: }
1090: }
1091: }
1092: /* Let the DM know we have set all the cell types */
1093: {
1094: DMLabel ctLabel;
1095: DM_Plex *plex = (DM_Plex *)rdm->data;
1097: DMPlexGetCellTypeLabel(rdm, &ctLabel);
1098: PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState);
1099: }
1100: return 0;
1101: }
1103: #if 0
1104: static PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize)
1105: {
1106: PetscInt ctNew;
1110: /* TODO Can do bisection since everything is sorted */
1111: for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) {
1112: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew]+1]];
1114: if (q >= ctSN && q < ctEN) break;
1115: }
1117: *coneSize = DMPolytopeTypeGetConeSize((DMPolytopeType) ctNew);
1118: return 0;
1119: }
1120: #endif
1122: /* The orientation o is for the interior of the cell p */
1123: 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[])
1124: {
1125: DM dm;
1126: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1127: const PetscInt *cone;
1128: PetscInt c, coff = *coneoff, ooff = *orntoff;
1130: DMPlexTransformGetDM(tr, &dm);
1131: DMPlexGetCone(dm, p, &cone);
1132: for (c = 0; c < csizeNew; ++c) {
1133: PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */
1134: PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */
1135: PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */
1136: DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */
1137: const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */
1138: PetscInt pr = -1; /* Replica number of pp that produces new cone point */
1139: const DMPolytopeType ft = (DMPolytopeType)rcone[coff++]; /* Cell type for new cone point of pNew */
1140: const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1141: PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */
1142: PetscInt lc;
1144: /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
1145: for (lc = 0; lc < fn; ++lc) {
1146: const PetscInt *parr = DMPolytopeTypeGetArrangment(pct, po);
1147: const PetscInt acp = rcone[coff++];
1148: const PetscInt pcp = parr[acp * 2];
1149: const PetscInt pco = parr[acp * 2 + 1];
1150: const PetscInt *ppornt;
1152: ppp = pp;
1153: pp = pcone[pcp];
1154: DMPlexGetCellType(dm, pp, &pct);
1155: DMPlexGetCone(dm, pp, &pcone);
1156: DMPlexGetConeOrientation(dm, ppp, &ppornt);
1157: po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco);
1158: }
1159: pr = rcone[coff++];
1160: /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
1161: DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo);
1162: DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c]);
1163: orntNew[c] = fo;
1164: }
1165: *coneoff = coff;
1166: *orntoff = ooff;
1167: return 0;
1168: }
1170: static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm)
1171: {
1172: DM dm;
1173: DMPolytopeType ct;
1174: PetscInt *coneNew, *orntNew;
1175: PetscInt maxConeSize = 0, pStart, pEnd, p, pNew;
1177: DMPlexTransformGetDM(tr, &dm);
1178: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1179: DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew);
1180: DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew);
1181: DMPlexGetChart(dm, &pStart, &pEnd);
1182: for (p = pStart; p < pEnd; ++p) {
1183: const PetscInt *cone, *ornt;
1184: PetscInt coff, ooff;
1185: DMPolytopeType *rct;
1186: PetscInt *rsize, *rcone, *rornt;
1187: PetscInt Nct, n, r;
1189: DMPlexGetCellType(dm, p, &ct);
1190: DMPlexGetCone(dm, p, &cone);
1191: DMPlexGetConeOrientation(dm, p, &ornt);
1192: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1193: for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
1194: const DMPolytopeType ctNew = rct[n];
1196: for (r = 0; r < rsize[n]; ++r) {
1197: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1198: DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew);
1199: DMPlexSetCone(rdm, pNew, coneNew);
1200: DMPlexSetConeOrientation(rdm, pNew, orntNew);
1201: }
1202: }
1203: }
1204: DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew);
1205: DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew);
1206: DMViewFromOptions(rdm, NULL, "-rdm_view");
1207: DMPlexSymmetrize(rdm);
1208: DMPlexStratify(rdm);
1209: return 0;
1210: }
1212: PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[])
1213: {
1214: DM dm;
1215: DMPolytopeType ct, qct;
1216: DMPolytopeType *rct;
1217: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1218: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1223: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1224: DMPlexTransformGetDM(tr, &dm);
1225: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone);
1226: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt);
1227: DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r);
1228: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1229: for (n = 0; n < Nct; ++n) {
1230: const DMPolytopeType ctNew = rct[n];
1231: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1232: PetscInt Nr = rsize[n], fn, c;
1234: if (ctNew == qct) Nr = r;
1235: for (nr = 0; nr < Nr; ++nr) {
1236: for (c = 0; c < csizeNew; ++c) {
1237: ++coff; /* Cell type of new cone point */
1238: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1239: coff += fn;
1240: ++coff; /* Replica number of new cone point */
1241: ++ooff; /* Orientation of new cone point */
1242: }
1243: }
1244: if (ctNew == qct) break;
1245: }
1246: DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt);
1247: *cone = qcone;
1248: *ornt = qornt;
1249: return 0;
1250: }
1252: PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1253: {
1254: DM dm;
1255: DMPolytopeType ct, qct;
1256: DMPolytopeType *rct;
1257: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1258: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1262: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p));
1263: DMPlexTransformGetDM(tr, &dm);
1264: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone);
1265: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt);
1266: DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r);
1267: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1268: for (n = 0; n < Nct; ++n) {
1269: const DMPolytopeType ctNew = rct[n];
1270: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1271: PetscInt Nr = rsize[n], fn, c;
1273: if (ctNew == qct) Nr = r;
1274: for (nr = 0; nr < Nr; ++nr) {
1275: for (c = 0; c < csizeNew; ++c) {
1276: ++coff; /* Cell type of new cone point */
1277: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1278: coff += fn;
1279: ++coff; /* Replica number of new cone point */
1280: ++ooff; /* Orientation of new cone point */
1281: }
1282: }
1283: if (ctNew == qct) break;
1284: }
1285: DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt);
1286: *cone = qcone;
1287: *ornt = qornt;
1288: return 0;
1289: }
1291: PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1292: {
1293: DM dm;
1297: DMPlexTransformGetDM(tr, &dm);
1298: DMRestoreWorkArray(dm, 0, MPIU_INT, cone);
1299: DMRestoreWorkArray(dm, 0, MPIU_INT, ornt);
1300: return 0;
1301: }
1303: static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr)
1304: {
1305: PetscInt ict;
1307: PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts);
1308: for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) {
1309: const DMPolytopeType ct = (DMPolytopeType)ict;
1310: DMPlexTransform reftr;
1311: DM refdm, trdm;
1312: Vec coordinates;
1313: const PetscScalar *coords;
1314: DMPolytopeType *rct;
1315: PetscInt *rsize, *rcone, *rornt;
1316: PetscInt Nct, n, r, pNew;
1317: PetscInt trdim, vStart, vEnd, Nc;
1318: const PetscInt debug = 0;
1319: const char *typeName;
1321: /* Since points are 0-dimensional, coordinates make no sense */
1322: if (DMPolytopeTypeGetDim(ct) <= 0) continue;
1323: DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm);
1324: DMPlexTransformCreate(PETSC_COMM_SELF, &reftr);
1325: DMPlexTransformSetDM(reftr, refdm);
1326: DMPlexTransformGetType(tr, &typeName);
1327: DMPlexTransformSetType(reftr, typeName);
1328: DMPlexTransformSetUp(reftr);
1329: DMPlexTransformApply(reftr, refdm, &trdm);
1331: DMGetDimension(trdm, &trdim);
1332: DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd);
1333: tr->trNv[ct] = vEnd - vStart;
1334: DMGetCoordinatesLocal(trdm, &coordinates);
1335: VecGetLocalSize(coordinates, &Nc);
1337: PetscCalloc1(Nc, &tr->trVerts[ct]);
1338: VecGetArrayRead(coordinates, &coords);
1339: PetscArraycpy(tr->trVerts[ct], coords, Nc);
1340: VecRestoreArrayRead(coordinates, &coords);
1342: PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct]);
1343: DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1344: for (n = 0; n < Nct; ++n) {
1345: /* Since points are 0-dimensional, coordinates make no sense */
1346: if (rct[n] == DM_POLYTOPE_POINT) continue;
1347: PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]]);
1348: for (r = 0; r < rsize[n]; ++r) {
1349: PetscInt *closure = NULL;
1350: PetscInt clSize, cl, Nv = 0;
1352: PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r]);
1353: DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew);
1354: DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure);
1355: for (cl = 0; cl < clSize * 2; cl += 2) {
1356: const PetscInt sv = closure[cl];
1358: if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart;
1359: }
1360: DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure);
1362: }
1363: }
1364: if (debug) {
1365: DMPolytopeType *rct;
1366: PetscInt *rsize, *rcone, *rornt;
1367: PetscInt v, dE = trdim, d, off = 0;
1369: PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct]);
1370: for (v = 0; v < tr->trNv[ct]; ++v) {
1371: PetscPrintf(PETSC_COMM_SELF, " ");
1372: for (d = 0; d < dE; ++d) PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++]));
1373: PetscPrintf(PETSC_COMM_SELF, "\n");
1374: }
1376: DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1377: for (n = 0; n < Nct; ++n) {
1378: if (rct[n] == DM_POLYTOPE_POINT) continue;
1379: PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct]);
1380: for (r = 0; r < rsize[n]; ++r) {
1381: PetscPrintf(PETSC_COMM_SELF, " ");
1382: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v]);
1383: PetscPrintf(PETSC_COMM_SELF, "\n");
1384: }
1385: }
1386: }
1387: DMDestroy(&refdm);
1388: DMDestroy(&trdm);
1389: DMPlexTransformDestroy(&reftr);
1390: }
1391: return 0;
1392: }
1394: /*
1395: DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type
1397: Input Parameters:
1398: + tr - The DMPlexTransform object
1399: - ct - The cell type
1401: Output Parameters:
1402: + Nv - The number of transformed vertices in the closure of the reference cell of given type
1403: - trVerts - The coordinates of these vertices in the reference cell
1405: Level: developer
1407: .seealso: `DMPlexTransformGetSubcellVertices()`
1408: */
1409: PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[])
1410: {
1411: if (!tr->trNv) DMPlexTransformCreateCellVertices_Internal(tr);
1412: if (Nv) *Nv = tr->trNv[ct];
1413: if (trVerts) *trVerts = tr->trVerts[ct];
1414: return 0;
1415: }
1417: /*
1418: DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type
1420: Input Parameters:
1421: + tr - The DMPlexTransform object
1422: . ct - The cell type
1423: . rct - The subcell type
1424: - r - The subcell index
1426: Output Parameter:
1427: . subVerts - The indices of these vertices in the set of vertices returned by DMPlexTransformGetCellVertices()
1429: Level: developer
1431: .seealso: `DMPlexTransformGetCellVertices()`
1432: */
1433: PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[])
1434: {
1435: if (!tr->trNv) DMPlexTransformCreateCellVertices_Internal(tr);
1437: if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r];
1438: return 0;
1439: }
1441: /* Computes new vertex as the barycenter, or centroid */
1442: PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1443: {
1444: PetscInt v, d;
1448: for (d = 0; d < dE; ++d) out[d] = 0.0;
1449: for (v = 0; v < Nv; ++v)
1450: for (d = 0; d < dE; ++d) out[d] += in[v * dE + d];
1451: for (d = 0; d < dE; ++d) out[d] /= Nv;
1452: return 0;
1453: }
1455: /*@
1456: DMPlexTransformMapCoordinates - Calculate new coordinates for produced points
1458: Not collective
1460: Input Parameters:
1461: + tr - The `DMPlexTransform`
1462: . pct - The cell type of the parent, from whom the new cell is being produced
1463: . ct - The type being produced
1464: . p - The original point
1465: . r - The replica number requested for the produced cell type
1466: . Nv - Number of vertices in the closure of the parent cell
1467: . dE - Spatial dimension
1468: - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell
1470: Output Parameter:
1471: . out - The coordinates of the new vertices
1473: Level: intermediate
1475: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransform`, `DMPlexTransformApply()`
1476: @*/
1477: PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1478: {
1480: PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out);
1481: return 0;
1482: }
1484: static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew)
1485: {
1486: DM dm;
1487: IS valueIS;
1488: const PetscInt *values;
1489: PetscInt defVal, Nv, val;
1491: DMPlexTransformGetDM(tr, &dm);
1492: DMLabelGetDefaultValue(label, &defVal);
1493: DMLabelSetDefaultValue(labelNew, defVal);
1494: DMLabelGetValueIS(label, &valueIS);
1495: ISGetLocalSize(valueIS, &Nv);
1496: ISGetIndices(valueIS, &values);
1497: for (val = 0; val < Nv; ++val) {
1498: IS pointIS;
1499: const PetscInt *points;
1500: PetscInt numPoints, p;
1502: /* Ensure refined label is created with same number of strata as
1503: * original (even if no entries here). */
1504: DMLabelAddStratum(labelNew, values[val]);
1505: DMLabelGetStratumIS(label, values[val], &pointIS);
1506: ISGetLocalSize(pointIS, &numPoints);
1507: ISGetIndices(pointIS, &points);
1508: for (p = 0; p < numPoints; ++p) {
1509: const PetscInt point = points[p];
1510: DMPolytopeType ct;
1511: DMPolytopeType *rct;
1512: PetscInt *rsize, *rcone, *rornt;
1513: PetscInt Nct, n, r, pNew = 0;
1515: DMPlexGetCellType(dm, point, &ct);
1516: DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1517: for (n = 0; n < Nct; ++n) {
1518: for (r = 0; r < rsize[n]; ++r) {
1519: DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew);
1520: DMLabelSetValue(labelNew, pNew, values[val]);
1521: }
1522: }
1523: }
1524: ISRestoreIndices(pointIS, &points);
1525: ISDestroy(&pointIS);
1526: }
1527: ISRestoreIndices(valueIS, &values);
1528: ISDestroy(&valueIS);
1529: return 0;
1530: }
1532: static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm)
1533: {
1534: DM dm;
1535: PetscInt numLabels, l;
1537: DMPlexTransformGetDM(tr, &dm);
1538: DMGetNumLabels(dm, &numLabels);
1539: for (l = 0; l < numLabels; ++l) {
1540: DMLabel label, labelNew;
1541: const char *lname;
1542: PetscBool isDepth, isCellType;
1544: DMGetLabelName(dm, l, &lname);
1545: PetscStrcmp(lname, "depth", &isDepth);
1546: if (isDepth) continue;
1547: PetscStrcmp(lname, "celltype", &isCellType);
1548: if (isCellType) continue;
1549: DMCreateLabel(rdm, lname);
1550: DMGetLabel(dm, lname, &label);
1551: DMGetLabel(rdm, lname, &labelNew);
1552: RefineLabel_Internal(tr, label, labelNew);
1553: }
1554: return 0;
1555: }
1557: /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */
1558: PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm)
1559: {
1560: DM dm;
1561: PetscInt Nf, f, Nds, s;
1563: DMPlexTransformGetDM(tr, &dm);
1564: DMGetNumFields(dm, &Nf);
1565: for (f = 0; f < Nf; ++f) {
1566: DMLabel label, labelNew;
1567: PetscObject obj;
1568: const char *lname;
1570: DMGetField(rdm, f, &label, &obj);
1571: if (!label) continue;
1572: PetscObjectGetName((PetscObject)label, &lname);
1573: DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
1574: RefineLabel_Internal(tr, label, labelNew);
1575: DMSetField_Internal(rdm, f, labelNew, obj);
1576: DMLabelDestroy(&labelNew);
1577: }
1578: DMGetNumDS(dm, &Nds);
1579: for (s = 0; s < Nds; ++s) {
1580: DMLabel label, labelNew;
1581: const char *lname;
1583: DMGetRegionNumDS(rdm, s, &label, NULL, NULL);
1584: if (!label) continue;
1585: PetscObjectGetName((PetscObject)label, &lname);
1586: DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
1587: RefineLabel_Internal(tr, label, labelNew);
1588: DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL);
1589: DMLabelDestroy(&labelNew);
1590: }
1591: return 0;
1592: }
1594: static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm)
1595: {
1596: DM dm;
1597: PetscSF sf, sfNew;
1598: PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m;
1599: const PetscInt *localPoints;
1600: const PetscSFNode *remotePoints;
1601: PetscInt *localPointsNew;
1602: PetscSFNode *remotePointsNew;
1603: PetscInt pStartNew, pEndNew, pNew;
1604: /* Brute force algorithm */
1605: PetscSF rsf;
1606: PetscSection s;
1607: const PetscInt *rootdegree;
1608: PetscInt *rootPointsNew, *remoteOffsets;
1609: PetscInt numPointsNew, pStart, pEnd, p;
1611: DMPlexTransformGetDM(tr, &dm);
1612: DMPlexGetChart(rdm, &pStartNew, &pEndNew);
1613: DMGetPointSF(dm, &sf);
1614: DMGetPointSF(rdm, &sfNew);
1615: /* Calculate size of new SF */
1616: PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);
1617: if (numRoots < 0) return 0;
1618: for (l = 0; l < numLeaves; ++l) {
1619: const PetscInt p = localPoints[l];
1620: DMPolytopeType ct;
1621: DMPolytopeType *rct;
1622: PetscInt *rsize, *rcone, *rornt;
1623: PetscInt Nct, n;
1625: DMPlexGetCellType(dm, p, &ct);
1626: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1627: for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n];
1628: }
1629: /* Send new root point numbers
1630: It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication
1631: */
1632: DMPlexGetChart(dm, &pStart, &pEnd);
1633: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s);
1634: PetscSectionSetChart(s, pStart, pEnd);
1635: for (p = pStart; p < pEnd; ++p) {
1636: DMPolytopeType ct;
1637: DMPolytopeType *rct;
1638: PetscInt *rsize, *rcone, *rornt;
1639: PetscInt Nct, n;
1641: DMPlexGetCellType(dm, p, &ct);
1642: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1643: for (n = 0; n < Nct; ++n) PetscSectionAddDof(s, p, rsize[n]);
1644: }
1645: PetscSectionSetUp(s);
1646: PetscSectionGetStorageSize(s, &numPointsNew);
1647: PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets);
1648: PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf);
1649: PetscFree(remoteOffsets);
1650: PetscSFComputeDegreeBegin(sf, &rootdegree);
1651: PetscSFComputeDegreeEnd(sf, &rootdegree);
1652: PetscMalloc1(numPointsNew, &rootPointsNew);
1653: for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1;
1654: for (p = pStart; p < pEnd; ++p) {
1655: DMPolytopeType ct;
1656: DMPolytopeType *rct;
1657: PetscInt *rsize, *rcone, *rornt;
1658: PetscInt Nct, n, r, off;
1660: if (!rootdegree[p - pStart]) continue;
1661: PetscSectionGetOffset(s, p, &off);
1662: DMPlexGetCellType(dm, p, &ct);
1663: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1664: for (n = 0, m = 0; n < Nct; ++n) {
1665: for (r = 0; r < rsize[n]; ++r, ++m) {
1666: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1667: rootPointsNew[off + m] = pNew;
1668: }
1669: }
1670: }
1671: PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE);
1672: PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE);
1673: PetscSFDestroy(&rsf);
1674: PetscMalloc1(numLeavesNew, &localPointsNew);
1675: PetscMalloc1(numLeavesNew, &remotePointsNew);
1676: for (l = 0, m = 0; l < numLeaves; ++l) {
1677: const PetscInt p = localPoints[l];
1678: DMPolytopeType ct;
1679: DMPolytopeType *rct;
1680: PetscInt *rsize, *rcone, *rornt;
1681: PetscInt Nct, n, r, q, off;
1683: PetscSectionGetOffset(s, p, &off);
1684: DMPlexGetCellType(dm, p, &ct);
1685: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1686: for (n = 0, q = 0; n < Nct; ++n) {
1687: for (r = 0; r < rsize[n]; ++r, ++m, ++q) {
1688: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1689: localPointsNew[m] = pNew;
1690: remotePointsNew[m].index = rootPointsNew[off + q];
1691: remotePointsNew[m].rank = remotePoints[l].rank;
1692: }
1693: }
1694: }
1695: PetscSectionDestroy(&s);
1696: PetscFree(rootPointsNew);
1697: /* SF needs sorted leaves to correctly calculate Gather */
1698: {
1699: PetscSFNode *rp, *rtmp;
1700: PetscInt *lp, *idx, *ltmp, i;
1702: PetscMalloc1(numLeavesNew, &idx);
1703: PetscMalloc1(numLeavesNew, &lp);
1704: PetscMalloc1(numLeavesNew, &rp);
1705: for (i = 0; i < numLeavesNew; ++i) {
1707: idx[i] = i;
1708: }
1709: PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx);
1710: for (i = 0; i < numLeavesNew; ++i) {
1711: lp[i] = localPointsNew[idx[i]];
1712: rp[i] = remotePointsNew[idx[i]];
1713: }
1714: ltmp = localPointsNew;
1715: localPointsNew = lp;
1716: rtmp = remotePointsNew;
1717: remotePointsNew = rp;
1718: PetscFree(idx);
1719: PetscFree(ltmp);
1720: PetscFree(rtmp);
1721: }
1722: PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
1723: return 0;
1724: }
1726: /*@C
1727: DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of type ct with localized coordinates x, we generate localized coordinates xr for subcell r of type rct.
1729: Not collective
1731: Input Parameters:
1732: + cr - The DMPlexCellRefiner
1733: . ct - The type of the parent cell
1734: . rct - The type of the produced cell
1735: . r - The index of the produced cell
1736: - x - The localized coordinates for the parent cell
1738: Output Parameter:
1739: . xr - The localized coordinates for the produced cell
1741: Level: developer
1743: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexCellRefinerSetCoordinates()`
1744: @*/
1745: static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
1746: {
1747: PetscFE fe = NULL;
1748: PetscInt cdim, v, *subcellV;
1750: DMPlexTransformGetCoordinateFE(tr, ct, &fe);
1751: DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV);
1752: PetscFEGetNumComponents(fe, &cdim);
1753: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim]);
1754: return 0;
1755: }
1757: static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm)
1758: {
1759: DM dm, cdm, cdmCell, cdmNew, cdmCellNew;
1760: PetscSection coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew;
1761: Vec coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew;
1762: const PetscScalar *coords;
1763: PetscScalar *coordsNew;
1764: const PetscReal *maxCell, *Lstart, *L;
1765: PetscBool localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE;
1766: PetscInt dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p;
1768: DMPlexTransformGetDM(tr, &dm);
1769: DMGetCoordinateDM(dm, &cdm);
1770: DMGetCellCoordinateDM(dm, &cdmCell);
1771: DMGetCoordinatesLocalized(dm, &localized);
1772: DMGetPeriodicity(dm, &maxCell, &Lstart, &L);
1773: if (localized) {
1774: /* Localize coordinates of new vertices */
1775: localizeVertices = PETSC_TRUE;
1776: /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */
1777: if (!maxCell) localizeCells = PETSC_TRUE;
1778: }
1779: DMGetCoordinateSection(dm, &coordSection);
1780: PetscSectionGetFieldComponents(coordSection, 0, &dEo);
1781: if (maxCell) {
1782: PetscReal maxCellNew[3];
1784: for (d = 0; d < dEo; ++d) maxCellNew[d] = maxCell[d] / 2.0;
1785: DMSetPeriodicity(rdm, maxCellNew, Lstart, L);
1786: }
1787: DMGetCoordinateDim(rdm, &dE);
1788: PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew);
1789: PetscSectionSetNumFields(coordSectionNew, 1);
1790: PetscSectionSetFieldComponents(coordSectionNew, 0, dE);
1791: DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew);
1792: PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew);
1793: /* Localization should be inherited */
1794: /* Stefano calculates parent cells for each new cell for localization */
1795: /* Localized cells need coordinates of closure */
1796: for (v = vStartNew; v < vEndNew; ++v) {
1797: PetscSectionSetDof(coordSectionNew, v, dE);
1798: PetscSectionSetFieldDof(coordSectionNew, v, 0, dE);
1799: }
1800: PetscSectionSetUp(coordSectionNew);
1801: DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew);
1803: if (localizeCells) {
1804: DMGetCoordinateDM(rdm, &cdmNew);
1805: DMClone(cdmNew, &cdmCellNew);
1806: DMSetCellCoordinateDM(rdm, cdmCellNew);
1807: DMDestroy(&cdmCellNew);
1809: PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew);
1810: PetscSectionSetNumFields(coordSectionCellNew, 1);
1811: PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE);
1812: DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew);
1813: PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew);
1815: DMGetCellCoordinateSection(dm, &coordSectionCell);
1816: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1817: for (c = cStart; c < cEnd; ++c) {
1818: PetscInt dof;
1820: PetscSectionGetDof(coordSectionCell, c, &dof);
1821: if (dof) {
1822: DMPolytopeType ct;
1823: DMPolytopeType *rct;
1824: PetscInt *rsize, *rcone, *rornt;
1825: PetscInt dim, cNew, Nct, n, r;
1827: DMPlexGetCellType(dm, c, &ct);
1828: dim = DMPolytopeTypeGetDim(ct);
1829: DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1830: /* This allows for different cell types */
1831: for (n = 0; n < Nct; ++n) {
1832: if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
1833: for (r = 0; r < rsize[n]; ++r) {
1834: PetscInt *closure = NULL;
1835: PetscInt clSize, cl, Nv = 0;
1837: DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew);
1838: DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
1839: for (cl = 0; cl < clSize * 2; cl += 2) {
1840: if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;
1841: }
1842: DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
1843: PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE);
1844: PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE);
1845: }
1846: }
1847: }
1848: }
1849: PetscSectionSetUp(coordSectionCellNew);
1850: DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew);
1851: }
1852: DMViewFromOptions(dm, NULL, "-coarse_dm_view");
1853: {
1854: VecType vtype;
1855: PetscInt coordSizeNew, bs;
1856: const char *name;
1858: DMGetCoordinatesLocal(dm, &coordsLocal);
1859: VecCreate(PETSC_COMM_SELF, &coordsLocalNew);
1860: PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);
1861: VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE);
1862: PetscObjectGetName((PetscObject)coordsLocal, &name);
1863: PetscObjectSetName((PetscObject)coordsLocalNew, name);
1864: VecGetBlockSize(coordsLocal, &bs);
1865: VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE);
1866: VecGetType(coordsLocal, &vtype);
1867: VecSetType(coordsLocalNew, vtype);
1868: }
1869: VecGetArrayRead(coordsLocal, &coords);
1870: VecGetArray(coordsLocalNew, &coordsNew);
1871: DMPlexGetChart(dm, &pStart, &pEnd);
1872: /* First set coordinates for vertices */
1873: for (p = pStart; p < pEnd; ++p) {
1874: DMPolytopeType ct;
1875: DMPolytopeType *rct;
1876: PetscInt *rsize, *rcone, *rornt;
1877: PetscInt Nct, n, r;
1878: PetscBool hasVertex = PETSC_FALSE;
1880: DMPlexGetCellType(dm, p, &ct);
1881: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1882: for (n = 0; n < Nct; ++n) {
1883: if (rct[n] == DM_POLYTOPE_POINT) {
1884: hasVertex = PETSC_TRUE;
1885: break;
1886: }
1887: }
1888: if (hasVertex) {
1889: const PetscScalar *icoords = NULL;
1890: const PetscScalar *array = NULL;
1891: PetscScalar *pcoords = NULL;
1892: PetscBool isDG;
1893: PetscInt Nc, Nv, v, d;
1895: DMPlexGetCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords);
1897: icoords = pcoords;
1898: Nv = Nc / dEo;
1899: if (ct != DM_POLYTOPE_POINT) {
1900: if (localizeVertices && maxCell) {
1901: PetscScalar anchor[3];
1903: for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d];
1904: for (v = 0; v < Nv; ++v) DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo]);
1905: }
1906: }
1907: for (n = 0; n < Nct; ++n) {
1908: if (rct[n] != DM_POLYTOPE_POINT) continue;
1909: for (r = 0; r < rsize[n]; ++r) {
1910: PetscScalar vcoords[3];
1911: PetscInt vNew, off;
1913: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew);
1914: PetscSectionGetOffset(coordSectionNew, vNew, &off);
1915: DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords);
1916: DMPlexSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off]);
1917: }
1918: }
1919: DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords);
1920: }
1921: }
1922: VecRestoreArrayRead(coordsLocal, &coords);
1923: VecRestoreArray(coordsLocalNew, &coordsNew);
1924: DMSetCoordinatesLocal(rdm, coordsLocalNew);
1925: VecDestroy(&coordsLocalNew);
1926: PetscSectionDestroy(&coordSectionNew);
1927: /* Then set coordinates for cells by localizing */
1928: if (!localizeCells) DMLocalizeCoordinates(rdm);
1929: else {
1930: VecType vtype;
1931: PetscInt coordSizeNew, bs;
1932: const char *name;
1934: DMGetCellCoordinatesLocal(dm, &coordsLocalCell);
1935: VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew);
1936: PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew);
1937: VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE);
1938: PetscObjectGetName((PetscObject)coordsLocalCell, &name);
1939: PetscObjectSetName((PetscObject)coordsLocalCellNew, name);
1940: VecGetBlockSize(coordsLocalCell, &bs);
1941: VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE);
1942: VecGetType(coordsLocalCell, &vtype);
1943: VecSetType(coordsLocalCellNew, vtype);
1944: VecGetArrayRead(coordsLocalCell, &coords);
1945: VecGetArray(coordsLocalCellNew, &coordsNew);
1947: for (p = pStart; p < pEnd; ++p) {
1948: DMPolytopeType ct;
1949: DMPolytopeType *rct;
1950: PetscInt *rsize, *rcone, *rornt;
1951: PetscInt dof = 0, Nct, n, r;
1953: DMPlexGetCellType(dm, p, &ct);
1954: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1955: if (p >= cStart && p < cEnd) PetscSectionGetDof(coordSectionCell, p, &dof);
1956: if (dof) {
1957: const PetscScalar *pcoords;
1959: DMPlexPointLocalRead(cdmCell, p, coords, &pcoords);
1960: for (n = 0; n < Nct; ++n) {
1961: const PetscInt Nr = rsize[n];
1963: if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
1964: for (r = 0; r < Nr; ++r) {
1965: PetscInt pNew, offNew;
1967: /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
1968: DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
1969: cell to the ones it produces. */
1970: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1971: PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew);
1972: DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew]);
1973: }
1974: }
1975: }
1976: }
1977: VecRestoreArrayRead(coordsLocalCell, &coords);
1978: VecRestoreArray(coordsLocalCellNew, &coordsNew);
1979: DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew);
1980: VecDestroy(&coordsLocalCellNew);
1981: PetscSectionDestroy(&coordSectionCellNew);
1982: }
1983: return 0;
1984: }
1986: PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm)
1987: {
1988: DM rdm;
1989: DMPlexInterpolatedFlag interp;
1994: DMPlexTransformSetDM(tr, dm);
1996: DMCreate(PetscObjectComm((PetscObject)dm), &rdm);
1997: DMSetType(rdm, DMPLEX);
1998: DMPlexTransformSetDimensions(tr, dm, rdm);
1999: /* Calculate number of new points of each depth */
2000: DMPlexIsInterpolatedCollective(dm, &interp);
2002: /* Step 1: Set chart */
2003: DMPlexSetChart(rdm, 0, tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]]);
2004: /* Step 2: Set cone/support sizes (automatically stratifies) */
2005: DMPlexTransformSetConeSizes(tr, rdm);
2006: /* Step 3: Setup refined DM */
2007: DMSetUp(rdm);
2008: /* Step 4: Set cones and supports (automatically symmetrizes) */
2009: DMPlexTransformSetCones(tr, rdm);
2010: /* Step 5: Create pointSF */
2011: DMPlexTransformCreateSF(tr, rdm);
2012: /* Step 6: Create labels */
2013: DMPlexTransformCreateLabels(tr, rdm);
2014: /* Step 7: Set coordinates */
2015: DMPlexTransformSetCoordinates(tr, rdm);
2016: DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, rdm);
2017: *tdm = rdm;
2018: return 0;
2019: }
2021: PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm)
2022: {
2023: DMPlexTransform tr;
2024: DM cdm, rcdm;
2025: const char *prefix;
2027: DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr);
2028: PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
2029: PetscObjectSetOptionsPrefix((PetscObject)tr, prefix);
2030: DMPlexTransformSetDM(tr, dm);
2031: DMPlexTransformSetFromOptions(tr);
2032: DMPlexTransformSetActive(tr, adaptLabel);
2033: DMPlexTransformSetUp(tr);
2034: PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view");
2035: DMPlexTransformApply(tr, dm, rdm);
2036: DMCopyDisc(dm, *rdm);
2037: DMGetCoordinateDM(dm, &cdm);
2038: DMGetCoordinateDM(*rdm, &rcdm);
2039: DMCopyDisc(cdm, rcdm);
2040: DMPlexTransformCreateDiscLabels(tr, *rdm);
2041: DMCopyDisc(dm, *rdm);
2042: DMPlexTransformDestroy(&tr);
2043: ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation;
2044: return 0;
2045: }