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) {
41: ctOInv[ctO[c]] = c;
42: }
43: *ctOrder = ctO;
44: *ctOrderInv = ctOInv;
45: return 0;
46: }
48: /*@C
49: DMPlexTransformRegister - Adds a new transform component implementation
51: Not Collective
53: Input Parameters:
54: + name - The name of a new user-defined creation routine
55: - create_func - The creation routine itself
57: Notes:
58: DMPlexTransformRegister() may be called multiple times to add several user-defined transforms
60: Sample usage:
61: .vb
62: DMPlexTransformRegister("my_transform", MyTransformCreate);
63: .ve
65: Then, your transform type can be chosen with the procedural interface via
66: .vb
67: DMPlexTransformCreate(MPI_Comm, DMPlexTransform *);
68: DMPlexTransformSetType(DMPlexTransform, "my_transform");
69: .ve
70: or at runtime via the option
71: .vb
72: -dm_plex_transform_type my_transform
73: .ve
75: Level: advanced
77: .seealso: DMPlexTransformRegisterAll(), DMPlexTransformRegisterDestroy()
78: @*/
79: PetscErrorCode DMPlexTransformRegister(const char name[], PetscErrorCode (*create_func)(DMPlexTransform))
80: {
81: DMInitializePackage();
82: PetscFunctionListAdd(&DMPlexTransformList, name, create_func);
83: return 0;
84: }
86: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Filter(DMPlexTransform);
87: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform);
88: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToBox(DMPlexTransform);
89: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform);
90: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform);
91: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform);
92: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_1D(DMPlexTransform);
93: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform);
95: /*@C
96: DMPlexTransformRegisterAll - Registers all of the transform components in the DM package.
98: Not Collective
100: Level: advanced
102: .seealso: DMRegisterAll(), DMPlexTransformRegisterDestroy()
103: @*/
104: PetscErrorCode DMPlexTransformRegisterAll(void)
105: {
106: if (DMPlexTransformRegisterAllCalled) return 0;
107: DMPlexTransformRegisterAllCalled = PETSC_TRUE;
109: DMPlexTransformRegister(DMPLEXTRANSFORMFILTER, DMPlexTransformCreate_Filter);
110: DMPlexTransformRegister(DMPLEXREFINEREGULAR, DMPlexTransformCreate_Regular);
111: DMPlexTransformRegister(DMPLEXREFINETOBOX, DMPlexTransformCreate_ToBox);
112: DMPlexTransformRegister(DMPLEXREFINEALFELD, DMPlexTransformCreate_Alfeld);
113: DMPlexTransformRegister(DMPLEXREFINEBOUNDARYLAYER, DMPlexTransformCreate_BL);
114: DMPlexTransformRegister(DMPLEXREFINESBR, DMPlexTransformCreate_SBR);
115: DMPlexTransformRegister(DMPLEXREFINE1D, DMPlexTransformCreate_1D);
116: DMPlexTransformRegister(DMPLEXEXTRUDE, DMPlexTransformCreate_Extrude);
117: return 0;
118: }
120: /*@C
121: DMPlexTransformRegisterDestroy - This function destroys the . It is called from PetscFinalize().
123: Level: developer
125: .seealso: PetscInitialize()
126: @*/
127: PetscErrorCode DMPlexTransformRegisterDestroy(void)
128: {
129: PetscFunctionListDestroy(&DMPlexTransformList);
130: DMPlexTransformRegisterAllCalled = PETSC_FALSE;
131: return 0;
132: }
134: /*@
135: DMPlexTransformCreate - Creates an empty transform object. The type can then be set with DMPlexTransformSetType().
137: Collective
139: Input Parameter:
140: . comm - The communicator for the transform object
142: Output Parameter:
143: . dm - The transform object
145: Level: beginner
147: .seealso: DMPlexTransformSetType(), DMPLEXREFINEREGULAR, DMPLEXTRANSFORMFILTER
148: @*/
149: PetscErrorCode DMPlexTransformCreate(MPI_Comm comm, DMPlexTransform *tr)
150: {
151: DMPlexTransform t;
154: *tr = NULL;
155: DMInitializePackage();
157: PetscHeaderCreate(t, DMPLEXTRANSFORM_CLASSID, "DMPlexTransform", "Mesh Transform", "DMPlexTransform", comm, DMPlexTransformDestroy, DMPlexTransformView);
158: t->setupcalled = PETSC_FALSE;
159: PetscCalloc2(DM_NUM_POLYTOPES, &t->coordFE, DM_NUM_POLYTOPES, &t->refGeom);
160: *tr = t;
161: return 0;
162: }
164: /*@C
165: DMSetType - Sets the particular implementation for a transform.
167: Collective on tr
169: Input Parameters:
170: + tr - The transform
171: - method - The name of the transform type
173: Options Database Key:
174: . -dm_plex_transform_type <type> - Sets the transform type; use -help for a list of available types
176: Notes:
177: See "petsc/include/petscdmplextransform.h" for available transform types
179: Level: intermediate
181: .seealso: DMPlexTransformGetType(), DMPlexTransformCreate()
182: @*/
183: PetscErrorCode DMPlexTransformSetType(DMPlexTransform tr, DMPlexTransformType method)
184: {
185: PetscErrorCode (*r)(DMPlexTransform);
186: PetscBool match;
189: PetscObjectTypeCompare((PetscObject) tr, method, &match);
190: if (match) return 0;
192: DMPlexTransformRegisterAll();
193: PetscFunctionListFind(DMPlexTransformList, method, &r);
196: if (tr->ops->destroy) (*tr->ops->destroy)(tr);
197: PetscMemzero(tr->ops, sizeof(*tr->ops));
198: PetscObjectChangeTypeName((PetscObject) tr, method);
199: (*r)(tr);
200: return 0;
201: }
203: /*@C
204: DMPlexTransformGetType - Gets the type name (as a string) from the transform.
206: Not Collective
208: Input Parameter:
209: . tr - The DMPlexTransform
211: Output Parameter:
212: . type - The DMPlexTransform type name
214: Level: intermediate
216: .seealso: DMPlexTransformSetType(), DMPlexTransformCreate()
217: @*/
218: PetscErrorCode DMPlexTransformGetType(DMPlexTransform tr, DMPlexTransformType *type)
219: {
222: DMPlexTransformRegisterAll();
223: *type = ((PetscObject) tr)->type_name;
224: return 0;
225: }
227: static PetscErrorCode DMPlexTransformView_Ascii(DMPlexTransform tr, PetscViewer v)
228: {
229: PetscViewerFormat format;
231: PetscViewerGetFormat(v, &format);
232: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
233: const PetscInt *trTypes = NULL;
234: IS trIS;
235: PetscInt cols = 8;
236: PetscInt Nrt = 8, f, g;
238: PetscViewerASCIIPrintf(v, "Source Starts\n");
239: for (g = 0; g <= cols; ++g) PetscViewerASCIIPrintf(v, " % 14s", DMPolytopeTypes[g]);
240: PetscViewerASCIIPrintf(v, "\n");
241: for (f = 0; f <= cols; ++f) PetscViewerASCIIPrintf(v, " % 14d", tr->ctStart[f]);
242: PetscViewerASCIIPrintf(v, "\n");
243: PetscViewerASCIIPrintf(v, "Target Starts\n");
244: for (g = 0; g <= cols; ++g) PetscViewerASCIIPrintf(v, " % 14s", DMPolytopeTypes[g]);
245: PetscViewerASCIIPrintf(v, "\n");
246: for (f = 0; f <= cols; ++f) PetscViewerASCIIPrintf(v, " % 14d", tr->ctStartNew[f]);
247: PetscViewerASCIIPrintf(v, "\n");
249: if (tr->trType) {
250: DMLabelGetNumValues(tr->trType, &Nrt);
251: DMLabelGetValueIS(tr->trType, &trIS);
252: ISGetIndices(trIS, &trTypes);
253: }
254: PetscViewerASCIIPrintf(v, "Offsets\n");
255: PetscViewerASCIIPrintf(v, " ");
256: for (g = 0; g < cols; ++g) {
257: PetscViewerASCIIPrintf(v, " % 14s", DMPolytopeTypes[g]);
258: }
259: PetscViewerASCIIPrintf(v, "\n");
260: for (f = 0; f < Nrt; ++f) {
261: PetscViewerASCIIPrintf(v, "%2d |", trTypes ? trTypes[f] : f);
262: for (g = 0; g < cols; ++g) {
263: PetscViewerASCIIPrintf(v, " % 14D", tr->offset[f*DM_NUM_POLYTOPES+g]);
264: }
265: PetscViewerASCIIPrintf(v, " |\n");
266: }
267: if (trTypes) {
268: ISGetIndices(trIS, &trTypes);
269: ISDestroy(&trIS);
270: }
271: }
272: return 0;
273: }
275: /*@C
276: DMPlexTransformView - Views a DMPlexTransform
278: Collective on tr
280: Input Parameters:
281: + tr - the DMPlexTransform object to view
282: - v - the viewer
284: Level: beginner
286: .seealso DMPlexTransformDestroy(), DMPlexTransformCreate()
287: @*/
288: PetscErrorCode DMPlexTransformView(DMPlexTransform tr, PetscViewer v)
289: {
290: PetscBool isascii;
293: if (!v) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) tr), &v);
296: PetscViewerCheckWritable(v);
297: PetscObjectPrintClassNamePrefixType((PetscObject) tr, v);
298: PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &isascii);
299: if (isascii) DMPlexTransformView_Ascii(tr, v);
300: if (tr->ops->view) (*tr->ops->view)(tr, v);
301: return 0;
302: }
304: /*@
305: DMPlexTransformSetFromOptions - Sets parameters in a transform from the options database
307: Collective on tr
309: Input Parameter:
310: . tr - the DMPlexTransform object to set options for
312: Options Database:
313: . -dm_plex_transform_type - Set the transform type, e.g. refine_regular
315: Level: intermediate
317: .seealso DMPlexTransformView(), DMPlexTransformCreate()
318: @*/
319: PetscErrorCode DMPlexTransformSetFromOptions(DMPlexTransform tr)
320: {
321: char typeName[1024];
322: const char *defName = DMPLEXREFINEREGULAR;
323: PetscBool flg;
327: PetscObjectOptionsBegin((PetscObject)tr);
328: PetscOptionsFList("-dm_plex_transform_type", "DMPlexTransform", "DMPlexTransformSetType", DMPlexTransformList, defName, typeName, 1024, &flg);
329: if (flg) DMPlexTransformSetType(tr, typeName);
330: else if (!((PetscObject) tr)->type_name) DMPlexTransformSetType(tr, defName);
331: if (tr->ops->setfromoptions) (*tr->ops->setfromoptions)(PetscOptionsObject,tr);
332: /* process any options handlers added with PetscObjectAddOptionsHandler() */
333: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) tr);
334: PetscOptionsEnd();
335: return 0;
336: }
338: /*@C
339: DMPlexTransformDestroy - Destroys a transform.
341: Collective on tr
343: Input Parameter:
344: . tr - the transform object to destroy
346: Level: beginner
348: .seealso DMPlexTransformView(), DMPlexTransformCreate()
349: @*/
350: PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr)
351: {
352: PetscInt c;
354: if (!*tr) return 0;
356: if (--((PetscObject) (*tr))->refct > 0) {*tr = NULL; return 0;}
358: if ((*tr)->ops->destroy) {
359: (*(*tr)->ops->destroy)(*tr);
360: }
361: DMDestroy(&(*tr)->dm);
362: DMLabelDestroy(&(*tr)->active);
363: DMLabelDestroy(&(*tr)->trType);
364: PetscFree2((*tr)->ctOrderOld, (*tr)->ctOrderInvOld);
365: PetscFree2((*tr)->ctOrderNew, (*tr)->ctOrderInvNew);
366: PetscFree2((*tr)->ctStart, (*tr)->ctStartNew);
367: PetscFree((*tr)->offset);
368: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
369: PetscFEDestroy(&(*tr)->coordFE[c]);
370: PetscFEGeomDestroy(&(*tr)->refGeom[c]);
371: }
372: if ((*tr)->trVerts) {
373: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
374: DMPolytopeType *rct;
375: PetscInt *rsize, *rcone, *rornt, Nct, n, r;
377: if (DMPolytopeTypeGetDim((DMPolytopeType) c) > 0) {
378: DMPlexTransformCellTransform((*tr), (DMPolytopeType) c, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
379: for (n = 0; n < Nct; ++n) {
381: if (rct[n] == DM_POLYTOPE_POINT) continue;
382: for (r = 0; r < rsize[n]; ++r) PetscFree((*tr)->trSubVerts[c][rct[n]][r]);
383: PetscFree((*tr)->trSubVerts[c][rct[n]]);
384: }
385: }
386: PetscFree((*tr)->trSubVerts[c]);
387: PetscFree((*tr)->trVerts[c]);
388: }
389: }
390: PetscFree3((*tr)->trNv, (*tr)->trVerts, (*tr)->trSubVerts);
391: PetscFree2((*tr)->coordFE, (*tr)->refGeom);
392: /* We do not destroy (*dm)->data here so that we can reference count backend objects */
393: PetscHeaderDestroy(tr);
394: return 0;
395: }
397: static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrderOld[], PetscInt ctStart[], PetscInt **offset)
398: {
399: DMLabel trType = tr->trType;
400: PetscInt c, cN, *off;
402: if (trType) {
403: DM dm;
404: IS rtIS;
405: const PetscInt *reftypes;
406: PetscInt Nrt, r;
408: DMPlexTransformGetDM(tr, &dm);
409: DMLabelGetNumValues(trType, &Nrt);
410: DMLabelGetValueIS(trType, &rtIS);
411: ISGetIndices(rtIS, &reftypes);
412: PetscCalloc1(Nrt*DM_NUM_POLYTOPES, &off);
413: for (r = 0; r < Nrt; ++r) {
414: const PetscInt rt = reftypes[r];
415: IS rtIS;
416: const PetscInt *points;
417: DMPolytopeType ct;
418: PetscInt p;
420: DMLabelGetStratumIS(trType, rt, &rtIS);
421: ISGetIndices(rtIS, &points);
422: p = points[0];
423: ISRestoreIndices(rtIS, &points);
424: ISDestroy(&rtIS);
425: DMPlexGetCellType(dm, p, &ct);
426: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
427: const DMPolytopeType ctNew = (DMPolytopeType) cN;
428: DMPolytopeType *rct;
429: PetscInt *rsize, *cone, *ornt;
430: PetscInt Nct, n, s;
432: if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {off[r*DM_NUM_POLYTOPES+ctNew] = -1; break;}
433: off[r*DM_NUM_POLYTOPES+ctNew] = 0;
434: for (s = 0; s <= r; ++s) {
435: const PetscInt st = reftypes[s];
436: DMPolytopeType sct;
437: PetscInt q, qrt;
439: DMLabelGetStratumIS(trType, st, &rtIS);
440: ISGetIndices(rtIS, &points);
441: q = points[0];
442: ISRestoreIndices(rtIS, &points);
443: ISDestroy(&rtIS);
444: DMPlexGetCellType(dm, q, &sct);
445: DMPlexTransformCellTransform(tr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt);
447: if (st == rt) {
448: for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) break;
449: if (n == Nct) off[r*DM_NUM_POLYTOPES+ctNew] = -1;
450: break;
451: }
452: for (n = 0; n < Nct; ++n) {
453: if (rct[n] == ctNew) {
454: PetscInt sn;
456: DMLabelGetStratumSize(trType, st, &sn);
457: off[r*DM_NUM_POLYTOPES+ctNew] += sn * rsize[n];
458: }
459: }
460: }
461: }
462: }
463: ISRestoreIndices(rtIS, &reftypes);
464: ISDestroy(&rtIS);
465: } else {
466: PetscCalloc1(DM_NUM_POLYTOPES*DM_NUM_POLYTOPES, &off);
467: for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
468: const DMPolytopeType ct = (DMPolytopeType) c;
469: for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
470: const DMPolytopeType ctNew = (DMPolytopeType) cN;
471: DMPolytopeType *rct;
472: PetscInt *rsize, *cone, *ornt;
473: PetscInt Nct, n, i;
475: if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {off[ct*DM_NUM_POLYTOPES+ctNew] = -1; continue;}
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) if (rct[n] == ctNew) break;
484: if (n == Nct) off[ct*DM_NUM_POLYTOPES+ctNew] = -1;
485: break;
486: }
487: for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) off[ct*DM_NUM_POLYTOPES+ctNew] += (ctStart[ictn]-ctStart[ict]) * rsize[n];
488: }
489: }
490: }
491: }
492: *offset = off;
493: return 0;
494: }
496: PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr)
497: {
498: DM dm;
499: DMPolytopeType ctCell;
500: PetscInt pStart, pEnd, p, c, celldim = 0;
503: if (tr->setupcalled) return 0;
504: if (tr->ops->setup) (*tr->ops->setup)(tr);
505: DMPlexTransformGetDM(tr, &dm);
506: DMPlexGetChart(dm, &pStart, &pEnd);
507: if (pEnd > pStart) {
508: DMPlexGetCellType(dm, 0, &ctCell);
509: } else {
510: PetscInt dim;
512: DMGetDimension(dm, &dim);
513: switch (dim) {
514: case 0: ctCell = DM_POLYTOPE_POINT;break;
515: case 1: ctCell = DM_POLYTOPE_SEGMENT;break;
516: case 2: ctCell = DM_POLYTOPE_TRIANGLE;break;
517: case 3: ctCell = DM_POLYTOPE_TETRAHEDRON;break;
518: default: ctCell = DM_POLYTOPE_UNKNOWN;
519: }
520: }
521: DMPlexCreateCellTypeOrder_Internal(DMPolytopeTypeGetDim(ctCell), &tr->ctOrderOld, &tr->ctOrderInvOld);
522: for (p = pStart; p < pEnd; ++p) {
523: DMPolytopeType ct;
524: DMPolytopeType *rct;
525: PetscInt *rsize, *cone, *ornt;
526: PetscInt Nct, n;
528: DMPlexGetCellType(dm, p, &ct);
530: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt);
531: for (n = 0; n < Nct; ++n) celldim = PetscMax(celldim, DMPolytopeTypeGetDim(rct[n]));
532: }
533: DMPlexCreateCellTypeOrder_Internal(celldim, &tr->ctOrderNew, &tr->ctOrderInvNew);
534: /* Construct sizes and offsets for each cell type */
535: if (!tr->ctStart) {
536: PetscInt *ctS, *ctSN, *ctC, *ctCN;
538: PetscCalloc2(DM_NUM_POLYTOPES+1, &ctS, DM_NUM_POLYTOPES+1, &ctSN);
539: PetscCalloc2(DM_NUM_POLYTOPES+1, &ctC, DM_NUM_POLYTOPES+1, &ctCN);
540: for (p = pStart; p < pEnd; ++p) {
541: DMPolytopeType ct;
542: DMPolytopeType *rct;
543: PetscInt *rsize, *cone, *ornt;
544: PetscInt Nct, n;
546: DMPlexGetCellType(dm, p, &ct);
548: ++ctC[ct];
549: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt);
550: for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
551: }
552: for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
553: const PetscInt cto = tr->ctOrderOld[c];
554: const PetscInt cton = tr->ctOrderOld[c+1];
555: const PetscInt ctn = tr->ctOrderNew[c];
556: const PetscInt ctnn = tr->ctOrderNew[c+1];
558: ctS[cton] = ctS[cto] + ctC[cto];
559: ctSN[ctnn] = ctSN[ctn] + ctCN[ctn];
560: }
561: PetscFree2(ctC, ctCN);
562: tr->ctStart = ctS;
563: tr->ctStartNew = ctSN;
564: }
565: DMPlexTransformCreateOffset_Internal(tr, tr->ctOrderOld, tr->ctStart, &tr->offset);
566: tr->setupcalled = PETSC_TRUE;
567: return 0;
568: }
570: PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm)
571: {
574: *dm = tr->dm;
575: return 0;
576: }
578: PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm)
579: {
582: PetscObjectReference((PetscObject) dm);
583: DMDestroy(&tr->dm);
584: tr->dm = dm;
585: return 0;
586: }
588: PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active)
589: {
592: *active = tr->active;
593: return 0;
594: }
596: PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active)
597: {
600: PetscObjectReference((PetscObject) active);
601: DMLabelDestroy(&tr->active);
602: tr->active = active;
603: return 0;
604: }
606: static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe)
607: {
608: if (!tr->coordFE[ct]) {
609: PetscInt dim, cdim;
611: dim = DMPolytopeTypeGetDim(ct);
612: DMGetCoordinateDim(tr->dm, &cdim);
613: PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, cdim, ct, 1, PETSC_DETERMINE, &tr->coordFE[ct]);
614: {
615: PetscDualSpace dsp;
616: PetscQuadrature quad;
617: DM K;
618: PetscFEGeom *cg;
619: PetscScalar *Xq;
620: PetscReal *xq, *wq;
621: PetscInt Nq, q;
623: DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq);
624: PetscMalloc1(Nq*cdim, &xq);
625: for (q = 0; q < Nq*cdim; ++q) xq[q] = PetscRealPart(Xq[q]);
626: PetscMalloc1(Nq, &wq);
627: for (q = 0; q < Nq; ++q) wq[q] = 1.0;
628: PetscQuadratureCreate(PETSC_COMM_SELF, &quad);
629: PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq);
630: PetscFESetQuadrature(tr->coordFE[ct], quad);
632: PetscFEGetDualSpace(tr->coordFE[ct], &dsp);
633: PetscDualSpaceGetDM(dsp, &K);
634: PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &tr->refGeom[ct]);
635: cg = tr->refGeom[ct];
636: DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
637: PetscQuadratureDestroy(&quad);
638: }
639: }
640: *fe = tr->coordFE[ct];
641: return 0;
642: }
644: /*@
645: DMPlexTransformSetDimensions - Set the dimensions for the transformed DM
647: Input Parameters:
648: + tr - The DMPlexTransform object
649: - dm - The original DM
651: Output Parameter:
652: . tdm - The transformed DM
654: Level: advanced
656: .seealso: DMPlexTransformApply(), DMPlexTransformCreate()
657: @*/
658: PetscErrorCode DMPlexTransformSetDimensions(DMPlexTransform tr, DM dm, DM tdm)
659: {
660: if (tr->ops->setdimensions) {
661: (*tr->ops->setdimensions)(tr, dm, tdm);
662: } else {
663: PetscInt dim, cdim;
665: DMGetDimension(dm, &dim);
666: DMSetDimension(tdm, dim);
667: DMGetCoordinateDim(dm, &cdim);
668: DMSetCoordinateDim(tdm, cdim);
669: }
670: return 0;
671: }
673: /*@
674: DMPlexTransformGetTargetPoint - Get the number of a point in the transformed mesh based on information from the original mesh.
676: Not collective
678: Input Parameters:
679: + tr - The DMPlexTransform
680: . ct - The type of the original point which produces the new point
681: . ctNew - The type of the new point
682: . p - The original point which produces the new point
683: - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p
685: Output Parameters:
686: . pNew - The new point number
688: Level: developer
690: .seealso: DMPlexTransformGetSourcePoint(), DMPlexTransformCellTransform()
691: @*/
692: PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
693: {
694: DMPolytopeType *rct;
695: PetscInt *rsize, *cone, *ornt;
696: PetscInt rt, Nct, n, off, rp;
697: DMLabel trType = tr->trType;
698: PetscInt ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ct]+1]];
699: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew]+1]];
700: PetscInt newp = ctSN, cind;
704: DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt);
705: if (trType) {
706: DMLabelGetValueIndex(trType, rt, &cind);
707: DMLabelGetStratumPointIndex(trType, rt, p, &rp);
709: } else {
710: cind = ct;
711: rp = p - ctS;
712: }
713: off = tr->offset[cind*DM_NUM_POLYTOPES + ctNew];
715: newp += off;
716: for (n = 0; n < Nct; ++n) {
717: if (rct[n] == ctNew) {
718: if (rsize[n] && r >= rsize[n])
719: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number %D should be in [0, %D) for subcell type %s in cell type %s", r, rsize[n], DMPolytopeTypes[rct[n]], DMPolytopeTypes[ct]);
720: newp += rp * rsize[n] + r;
721: break;
722: }
723: }
726: *pNew = newp;
727: return 0;
728: }
730: /*@
731: DMPlexTransformGetSourcePoint - Get the number of a point in the original mesh based on information from the transformed mesh.
733: Not collective
735: Input Parameters:
736: + tr - The DMPlexTransform
737: - pNew - The new point number
739: Output Parameters:
740: + ct - The type of the original point which produces the new point
741: . ctNew - The type of the new point
742: . p - The original point which produces the new point
743: - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p
745: Level: developer
747: .seealso: DMPlexTransformGetTargetPoint(), DMPlexTransformCellTransform()
748: @*/
749: PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r)
750: {
751: DMLabel trType = tr->trType;
752: DMPolytopeType *rct;
753: PetscInt *rsize, *cone, *ornt;
754: PetscInt rt, Nct, n, rp = 0, rO = 0, pO;
755: PetscInt offset = -1, ctS, ctE, ctO = 0, ctN, ctTmp;
757: for (ctN = 0; ctN < DM_NUM_POLYTOPES; ++ctN) {
758: PetscInt ctSN = tr->ctStartNew[ctN], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctN]+1]];
760: if ((pNew >= ctSN) && (pNew < ctEN)) break;
761: }
763: if (trType) {
764: DM dm;
765: IS rtIS;
766: const PetscInt *reftypes;
767: PetscInt Nrt, r, rtStart;
769: DMPlexTransformGetDM(tr, &dm);
770: DMLabelGetNumValues(trType, &Nrt);
771: DMLabelGetValueIS(trType, &rtIS);
772: ISGetIndices(rtIS, &reftypes);
773: for (r = 0; r < Nrt; ++r) {
774: const PetscInt off = tr->offset[r*DM_NUM_POLYTOPES + ctN];
776: if (tr->ctStartNew[ctN] + off > pNew) continue;
777: /* Check that any of this refinement type exist */
778: /* TODO Actually keep track of the number produced here instead */
779: if (off > offset) {rt = reftypes[r]; offset = off;}
780: }
781: ISRestoreIndices(rtIS, &reftypes);
782: ISDestroy(&rtIS);
784: /* TODO Map refinement types to cell types */
785: DMLabelGetStratumBounds(trType, rt, &rtStart, NULL);
787: for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) {
788: PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO]+1]];
790: if ((rtStart >= ctS) && (rtStart < ctE)) break;
791: }
793: } else {
794: for (ctTmp = 0; ctTmp < DM_NUM_POLYTOPES; ++ctTmp) {
795: const PetscInt off = tr->offset[ctTmp*DM_NUM_POLYTOPES + ctN];
797: if (tr->ctStartNew[ctN] + off > pNew) continue;
798: if (tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctTmp]+1]] <= tr->ctStart[ctTmp]) continue;
799: /* TODO Actually keep track of the number produced here instead */
800: if (off > offset) {ctO = ctTmp; offset = off;}
801: }
803: }
804: ctS = tr->ctStart[ctO];
805: ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO]+1]];
806: DMPlexTransformCellTransform(tr, (DMPolytopeType) ctO, ctS, &rt, &Nct, &rct, &rsize, &cone, &ornt);
807: for (n = 0; n < Nct; ++n) {
808: if ((PetscInt) rct[n] == ctN) {
809: PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset;
811: rp = tmp / rsize[n];
812: rO = tmp % rsize[n];
813: break;
814: }
815: }
817: pO = rp + ctS;
819: if (ct) *ct = (DMPolytopeType) ctO;
820: if (ctNew) *ctNew = (DMPolytopeType) ctN;
821: if (p) *p = pO;
822: if (r) *r = rO;
823: return 0;
824: }
826: /*@
827: DMPlexTransformCellTransform - Describes the transform of a given source cell into a set of other target cells. These produced cells become the new mesh.
829: Input Parameters:
830: + tr - The DMPlexTransform object
831: . source - The source cell type
832: - p - The source point, which can also determine the refine type
834: Output Parameters:
835: + rt - The refine type for this point
836: . Nt - The number of types produced by this point
837: . target - An array of length Nt giving the types produced
838: . size - An array of length Nt giving the number of cells of each type produced
839: . cone - An array of length Nt*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point
840: - ornt - An array of length Nt*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point
842: Note: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
843: need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
844: division (described in these outputs) of the cell in the original mesh. The point identifier is given by
845: $ the number of cones to be taken, or 0 for the current cell
846: $ the cell cone point number at each level from which it is subdivided
847: $ the replica number r of the subdivision.
848: The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
849: $ Nt = 2
850: $ target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
851: $ size = {1, 2}
852: $ cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
853: $ ornt = { 0, 0, 0, 0}
855: Level: advanced
857: .seealso: DMPlexTransformApply(), DMPlexTransformCreate()
858: @*/
859: PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
860: {
861: (*tr->ops->celltransform)(tr, source, p, rt, Nt, target, size, cone, ornt);
862: return 0;
863: }
865: PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
866: {
867: *rnew = r;
868: *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
869: return 0;
870: }
872: /* Returns the same thing */
873: PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
874: {
875: static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
876: static PetscInt vertexS[] = {1};
877: static PetscInt vertexC[] = {0};
878: static PetscInt vertexO[] = {0};
879: static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT};
880: static PetscInt edgeS[] = {1};
881: static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
882: static PetscInt edgeO[] = {0, 0};
883: static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR};
884: static PetscInt tedgeS[] = {1};
885: static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
886: static PetscInt tedgeO[] = {0, 0};
887: static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE};
888: static PetscInt triS[] = {1};
889: static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
890: static PetscInt triO[] = {0, 0, 0};
891: static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL};
892: static PetscInt quadS[] = {1};
893: 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};
894: static PetscInt quadO[] = {0, 0, 0, 0};
895: static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR};
896: static PetscInt tquadS[] = {1};
897: 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};
898: static PetscInt tquadO[] = {0, 0, 0, 0};
899: static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON};
900: static PetscInt tetS[] = {1};
901: 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};
902: static PetscInt tetO[] = {0, 0, 0, 0};
903: static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON};
904: static PetscInt hexS[] = {1};
905: static PetscInt hexC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0,
906: DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 1, 5, 0};
907: static PetscInt hexO[] = {0, 0, 0, 0, 0, 0};
908: static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM};
909: static PetscInt tripS[] = {1};
910: static PetscInt tripC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0,
911: DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0};
912: static PetscInt tripO[] = {0, 0, 0, 0, 0};
913: static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR};
914: static PetscInt ttripS[] = {1};
915: static PetscInt ttripC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0,
916: DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0};
917: static PetscInt ttripO[] = {0, 0, 0, 0, 0};
918: static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
919: static PetscInt tquadpS[] = {1};
920: static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0,
921: DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0};
922: static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0};
923: static DMPolytopeType pyrT[] = {DM_POLYTOPE_PYRAMID};
924: static PetscInt pyrS[] = {1};
925: static PetscInt pyrC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0,
926: DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 1, 4, 0};
927: static PetscInt pyrO[] = {0, 0, 0, 0, 0};
929: if (rt) *rt = 0;
930: switch (source) {
931: case DM_POLYTOPE_POINT: *Nt = 1; *target = vertexT; *size = vertexS; *cone = vertexC; *ornt = vertexO; break;
932: case DM_POLYTOPE_SEGMENT: *Nt = 1; *target = edgeT; *size = edgeS; *cone = edgeC; *ornt = edgeO; break;
933: case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break;
934: case DM_POLYTOPE_TRIANGLE: *Nt = 1; *target = triT; *size = triS; *cone = triC; *ornt = triO; break;
935: case DM_POLYTOPE_QUADRILATERAL: *Nt = 1; *target = quadT; *size = quadS; *cone = quadC; *ornt = quadO; break;
936: case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 1; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break;
937: case DM_POLYTOPE_TETRAHEDRON: *Nt = 1; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break;
938: case DM_POLYTOPE_HEXAHEDRON: *Nt = 1; *target = hexT; *size = hexS; *cone = hexC; *ornt = hexO; break;
939: case DM_POLYTOPE_TRI_PRISM: *Nt = 1; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break;
940: case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 1; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break;
941: case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 1; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break;
942: case DM_POLYTOPE_PYRAMID: *Nt = 1; *target = pyrT; *size = pyrS; *cone = pyrC; *ornt = pyrO; break;
943: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
944: }
945: return 0;
946: }
948: /*@
949: DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point
951: Not collective
953: Input Parameters:
954: + tr - The DMPlexTransform
955: . sct - The source point cell type, from whom the new cell is being produced
956: . sp - The source point
957: . so - The orientation of the source point in its enclosing parent
958: . tct - The target point cell type
959: . r - The replica number requested for the produced cell type
960: - o - The orientation of the replica
962: Output Parameters:
963: + rnew - The replica number, given the orientation of the parent
964: - onew - The replica orientation, given the orientation of the parent
966: Level: advanced
968: .seealso: DMPlexTransformCellTransform(), DMPlexTransformApply()
969: @*/
970: PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
971: {
973: (*tr->ops->getsubcellorientation)(tr, sct, sp, so, tct, r, o, rnew, onew);
974: return 0;
975: }
977: static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm)
978: {
979: DM dm;
980: PetscInt pStart, pEnd, p, pNew;
982: DMPlexTransformGetDM(tr, &dm);
983: /* Must create the celltype label here so that we do not automatically try to compute the types */
984: DMCreateLabel(rdm, "celltype");
985: DMPlexGetChart(dm, &pStart, &pEnd);
986: for (p = pStart; p < pEnd; ++p) {
987: DMPolytopeType ct;
988: DMPolytopeType *rct;
989: PetscInt *rsize, *rcone, *rornt;
990: PetscInt Nct, n, r;
992: DMPlexGetCellType(dm, p, &ct);
993: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
994: for (n = 0; n < Nct; ++n) {
995: for (r = 0; r < rsize[n]; ++r) {
996: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
997: DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]));
998: DMPlexSetCellType(rdm, pNew, rct[n]);
999: }
1000: }
1001: }
1002: /* Let the DM know we have set all the cell types */
1003: {
1004: DMLabel ctLabel;
1005: DM_Plex *plex = (DM_Plex *) rdm->data;
1007: DMPlexGetCellTypeLabel(rdm, &ctLabel);
1008: PetscObjectStateGet((PetscObject) ctLabel, &plex->celltypeState);
1009: }
1010: return 0;
1011: }
1013: #if 0
1014: static PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize)
1015: {
1016: PetscInt ctNew;
1020: /* TODO Can do bisection since everything is sorted */
1021: for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) {
1022: PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew]+1]];
1024: if (q >= ctSN && q < ctEN) break;
1025: }
1027: *coneSize = DMPolytopeTypeGetConeSize((DMPolytopeType) ctNew);
1028: return 0;
1029: }
1030: #endif
1032: /* The orientation o is for the interior of the cell p */
1033: static PetscErrorCode DMPlexTransformGetCone_Internal(DMPlexTransform tr, PetscInt p, PetscInt o, DMPolytopeType ct, DMPolytopeType ctNew,
1034: const PetscInt rcone[], PetscInt *coneoff, const PetscInt rornt[], PetscInt *orntoff,
1035: PetscInt coneNew[], PetscInt orntNew[])
1036: {
1037: DM dm;
1038: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1039: const PetscInt *cone;
1040: PetscInt c, coff = *coneoff, ooff = *orntoff;
1042: DMPlexTransformGetDM(tr, &dm);
1043: DMPlexGetCone(dm, p, &cone);
1044: for (c = 0; c < csizeNew; ++c) {
1045: PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */
1046: PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */
1047: PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */
1048: DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */
1049: const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */
1050: PetscInt pr = -1; /* Replica number of pp that produces new cone point */
1051: const DMPolytopeType ft = (DMPolytopeType) rcone[coff++]; /* Cell type for new cone point of pNew */
1052: const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1053: PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */
1054: PetscInt lc;
1056: /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
1057: for (lc = 0; lc < fn; ++lc) {
1058: const PetscInt *parr = DMPolytopeTypeGetArrangment(pct, po);
1059: const PetscInt acp = rcone[coff++];
1060: const PetscInt pcp = parr[acp*2];
1061: const PetscInt pco = parr[acp*2+1];
1062: const PetscInt *ppornt;
1064: ppp = pp;
1065: pp = pcone[pcp];
1066: DMPlexGetCellType(dm, pp, &pct);
1067: DMPlexGetCone(dm, pp, &pcone);
1068: DMPlexGetConeOrientation(dm, ppp, &ppornt);
1069: po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco);
1070: }
1071: pr = rcone[coff++];
1072: /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
1073: DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo);
1074: DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c]);
1075: orntNew[c] = fo;
1076: }
1077: *coneoff = coff;
1078: *orntoff = ooff;
1079: return 0;
1080: }
1082: static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm)
1083: {
1084: DM dm;
1085: DMPolytopeType ct;
1086: PetscInt *coneNew, *orntNew;
1087: PetscInt maxConeSize = 0, pStart, pEnd, p, pNew;
1089: DMPlexTransformGetDM(tr, &dm);
1090: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p));
1091: DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew);
1092: DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew);
1093: DMPlexGetChart(dm, &pStart, &pEnd);
1094: for (p = pStart; p < pEnd; ++p) {
1095: const PetscInt *cone, *ornt;
1096: PetscInt coff, ooff;
1097: DMPolytopeType *rct;
1098: PetscInt *rsize, *rcone, *rornt;
1099: PetscInt Nct, n, r;
1101: DMPlexGetCellType(dm, p, &ct);
1102: DMPlexGetCone(dm, p, &cone);
1103: DMPlexGetConeOrientation(dm, p, &ornt);
1104: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1105: for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) {
1106: const DMPolytopeType ctNew = rct[n];
1108: for (r = 0; r < rsize[n]; ++r) {
1109: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1110: DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew);
1111: DMPlexSetCone(rdm, pNew, coneNew);
1112: DMPlexSetConeOrientation(rdm, pNew, orntNew);
1113: }
1114: }
1115: }
1116: DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew);
1117: DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew);
1118: DMViewFromOptions(rdm, NULL, "-rdm_view");
1119: DMPlexSymmetrize(rdm);
1120: DMPlexStratify(rdm);
1121: return 0;
1122: }
1124: PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[])
1125: {
1126: DM dm;
1127: DMPolytopeType ct, qct;
1128: DMPolytopeType *rct;
1129: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1130: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1135: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p));
1136: DMPlexTransformGetDM(tr, &dm);
1137: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone);
1138: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt);
1139: DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r);
1140: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1141: for (n = 0; n < Nct; ++n) {
1142: const DMPolytopeType ctNew = rct[n];
1143: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1144: PetscInt Nr = rsize[n], fn, c;
1146: if (ctNew == qct) Nr = r;
1147: for (nr = 0; nr < Nr; ++nr) {
1148: for (c = 0; c < csizeNew; ++c) {
1149: ++coff; /* Cell type of new cone point */
1150: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1151: coff += fn;
1152: ++coff; /* Replica number of new cone point */
1153: ++ooff; /* Orientation of new cone point */
1154: }
1155: }
1156: if (ctNew == qct) break;
1157: }
1158: DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt);
1159: *cone = qcone;
1160: *ornt = qornt;
1161: return 0;
1162: }
1164: PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1165: {
1166: DM dm;
1167: DMPolytopeType ct, qct;
1168: DMPolytopeType *rct;
1169: PetscInt *rsize, *rcone, *rornt, *qcone, *qornt;
1170: PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1174: for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p));
1175: DMPlexTransformGetDM(tr, &dm);
1176: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone);
1177: DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt);
1178: DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r);
1179: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1180: for (n = 0; n < Nct; ++n) {
1181: const DMPolytopeType ctNew = rct[n];
1182: const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1183: PetscInt Nr = rsize[n], fn, c;
1185: if (ctNew == qct) Nr = r;
1186: for (nr = 0; nr < Nr; ++nr) {
1187: for (c = 0; c < csizeNew; ++c) {
1188: ++coff; /* Cell type of new cone point */
1189: fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1190: coff += fn;
1191: ++coff; /* Replica number of new cone point */
1192: ++ooff; /* Orientation of new cone point */
1193: }
1194: }
1195: if (ctNew == qct) break;
1196: }
1197: DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt);
1198: *cone = qcone;
1199: *ornt = qornt;
1200: return 0;
1201: }
1203: PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1204: {
1205: DM dm;
1209: DMPlexTransformGetDM(tr, &dm);
1210: DMRestoreWorkArray(dm, 0, MPIU_INT, cone);
1211: DMRestoreWorkArray(dm, 0, MPIU_INT, ornt);
1212: return 0;
1213: }
1215: static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr)
1216: {
1217: PetscInt ict;
1219: PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts);
1220: for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) {
1221: const DMPolytopeType ct = (DMPolytopeType) ict;
1222: DMPlexTransform reftr;
1223: DM refdm, trdm;
1224: Vec coordinates;
1225: const PetscScalar *coords;
1226: DMPolytopeType *rct;
1227: PetscInt *rsize, *rcone, *rornt;
1228: PetscInt Nct, n, r, pNew;
1229: PetscInt vStart, vEnd, Nc;
1230: const PetscInt debug = 0;
1231: const char *typeName;
1233: /* Since points are 0-dimensional, coordinates make no sense */
1234: if (DMPolytopeTypeGetDim(ct) <= 0) continue;
1235: DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm);
1236: DMPlexTransformCreate(PETSC_COMM_SELF, &reftr);
1237: DMPlexTransformSetDM(reftr, refdm);
1238: DMPlexTransformGetType(tr, &typeName);
1239: DMPlexTransformSetType(reftr, typeName);
1240: DMPlexTransformSetUp(reftr);
1241: DMPlexTransformApply(reftr, refdm, &trdm);
1243: DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd);
1244: tr->trNv[ct] = vEnd - vStart;
1245: DMGetCoordinatesLocal(trdm, &coordinates);
1246: VecGetLocalSize(coordinates, &Nc);
1248: PetscCalloc1(Nc, &tr->trVerts[ct]);
1249: VecGetArrayRead(coordinates, &coords);
1250: PetscArraycpy(tr->trVerts[ct], coords, Nc);
1251: VecRestoreArrayRead(coordinates, &coords);
1253: PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct]);
1254: DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1255: for (n = 0; n < Nct; ++n) {
1257: /* Since points are 0-dimensional, coordinates make no sense */
1258: if (rct[n] == DM_POLYTOPE_POINT) continue;
1259: PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]]);
1260: for (r = 0; r < rsize[n]; ++r) {
1261: PetscInt *closure = NULL;
1262: PetscInt clSize, cl, Nv = 0;
1264: PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r]);
1265: DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew);
1266: DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure);
1267: for (cl = 0; cl < clSize*2; cl += 2) {
1268: const PetscInt sv = closure[cl];
1270: if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart;
1271: }
1272: DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure);
1274: }
1275: }
1276: if (debug) {
1277: DMPolytopeType *rct;
1278: PetscInt *rsize, *rcone, *rornt;
1279: PetscInt v, dE = DMPolytopeTypeGetDim(ct), d, off = 0;
1281: PetscPrintf(PETSC_COMM_SELF, "%s: %D vertices\n", DMPolytopeTypes[ct], tr->trNv[ct]);
1282: for (v = 0; v < tr->trNv[ct]; ++v) {
1283: PetscPrintf(PETSC_COMM_SELF, " ");
1284: for (d = 0; d < dE; ++d) PetscPrintf(PETSC_COMM_SELF, "%g ", tr->trVerts[ct][off++]);
1285: PetscPrintf(PETSC_COMM_SELF, "\n");
1286: }
1288: DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1289: for (n = 0; n < Nct; ++n) {
1290: if (rct[n] == DM_POLYTOPE_POINT) continue;
1291: PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct]);
1292: for (r = 0; r < rsize[n]; ++r) {
1293: PetscPrintf(PETSC_COMM_SELF, " ");
1294: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) {
1295: PetscPrintf(PETSC_COMM_SELF, "%D ", tr->trSubVerts[ct][rct[n]][r][v]);
1296: }
1297: PetscPrintf(PETSC_COMM_SELF, "\n");
1298: }
1299: }
1300: }
1301: DMDestroy(&refdm);
1302: DMDestroy(&trdm);
1303: DMPlexTransformDestroy(&reftr);
1304: }
1305: return 0;
1306: }
1308: /*
1309: DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type
1311: Input Parameters:
1312: + tr - The DMPlexTransform object
1313: - ct - The cell type
1315: Output Parameters:
1316: + Nv - The number of transformed vertices in the closure of the reference cell of given type
1317: - trVerts - The coordinates of these vertices in the reference cell
1319: Level: developer
1321: .seealso: DMPlexTransformGetSubcellVertices()
1322: */
1323: PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[])
1324: {
1325: if (!tr->trNv) DMPlexTransformCreateCellVertices_Internal(tr);
1326: if (Nv) *Nv = tr->trNv[ct];
1327: if (trVerts) *trVerts = tr->trVerts[ct];
1328: return 0;
1329: }
1331: /*
1332: DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type
1334: Input Parameters:
1335: + tr - The DMPlexTransform object
1336: . ct - The cell type
1337: . rct - The subcell type
1338: - r - The subcell index
1340: Output Parameter:
1341: . subVerts - The indices of these vertices in the set of vertices returned by DMPlexTransformGetCellVertices()
1343: Level: developer
1345: .seealso: DMPlexTransformGetCellVertices()
1346: */
1347: PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[])
1348: {
1349: if (!tr->trNv) DMPlexTransformCreateCellVertices_Internal(tr);
1351: if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r];
1352: return 0;
1353: }
1355: /* Computes new vertex as the barycenter, or centroid */
1356: PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1357: {
1358: PetscInt v,d;
1362: for (d = 0; d < dE; ++d) out[d] = 0.0;
1363: for (v = 0; v < Nv; ++v) for (d = 0; d < dE; ++d) out[d] += in[v*dE+d];
1364: for (d = 0; d < dE; ++d) out[d] /= Nv;
1365: return 0;
1366: }
1368: /*@
1369: DMPlexTransformMapCoordinates - Calculate new coordinates for produced points
1371: Not collective
1373: Input Parameters:
1374: + cr - The DMPlexCellRefiner
1375: . pct - The cell type of the parent, from whom the new cell is being produced
1376: . ct - The type being produced
1377: . p - The original point
1378: . r - The replica number requested for the produced cell type
1379: . Nv - Number of vertices in the closure of the parent cell
1380: . dE - Spatial dimension
1381: - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell
1383: Output Parameter:
1384: . out - The coordinates of the new vertices
1386: Level: intermediate
1388: .seealso: DMPlexTransformApply()
1389: @*/
1390: PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1391: {
1393: (*tr->ops->mapcoordinates)(tr, pct, ct, p, r, Nv, dE, in, out);
1394: return 0;
1395: }
1397: static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew)
1398: {
1399: DM dm;
1400: IS valueIS;
1401: const PetscInt *values;
1402: PetscInt defVal, Nv, val;
1404: DMPlexTransformGetDM(tr, &dm);
1405: DMLabelGetDefaultValue(label, &defVal);
1406: DMLabelSetDefaultValue(labelNew, defVal);
1407: DMLabelGetValueIS(label, &valueIS);
1408: ISGetLocalSize(valueIS, &Nv);
1409: ISGetIndices(valueIS, &values);
1410: for (val = 0; val < Nv; ++val) {
1411: IS pointIS;
1412: const PetscInt *points;
1413: PetscInt numPoints, p;
1415: /* Ensure refined label is created with same number of strata as
1416: * original (even if no entries here). */
1417: DMLabelAddStratum(labelNew, values[val]);
1418: DMLabelGetStratumIS(label, values[val], &pointIS);
1419: ISGetLocalSize(pointIS, &numPoints);
1420: ISGetIndices(pointIS, &points);
1421: for (p = 0; p < numPoints; ++p) {
1422: const PetscInt point = points[p];
1423: DMPolytopeType ct;
1424: DMPolytopeType *rct;
1425: PetscInt *rsize, *rcone, *rornt;
1426: PetscInt Nct, n, r, pNew=0;
1428: DMPlexGetCellType(dm, point, &ct);
1429: DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1430: for (n = 0; n < Nct; ++n) {
1431: for (r = 0; r < rsize[n]; ++r) {
1432: DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew);
1433: DMLabelSetValue(labelNew, pNew, values[val]);
1434: }
1435: }
1436: }
1437: ISRestoreIndices(pointIS, &points);
1438: ISDestroy(&pointIS);
1439: }
1440: ISRestoreIndices(valueIS, &values);
1441: ISDestroy(&valueIS);
1442: return 0;
1443: }
1445: static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm)
1446: {
1447: DM dm;
1448: PetscInt numLabels, l;
1450: DMPlexTransformGetDM(tr, &dm);
1451: DMGetNumLabels(dm, &numLabels);
1452: for (l = 0; l < numLabels; ++l) {
1453: DMLabel label, labelNew;
1454: const char *lname;
1455: PetscBool isDepth, isCellType;
1457: DMGetLabelName(dm, l, &lname);
1458: PetscStrcmp(lname, "depth", &isDepth);
1459: if (isDepth) continue;
1460: PetscStrcmp(lname, "celltype", &isCellType);
1461: if (isCellType) continue;
1462: DMCreateLabel(rdm, lname);
1463: DMGetLabel(dm, lname, &label);
1464: DMGetLabel(rdm, lname, &labelNew);
1465: RefineLabel_Internal(tr, label, labelNew);
1466: }
1467: return 0;
1468: }
1470: /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */
1471: PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm)
1472: {
1473: DM dm;
1474: PetscInt Nf, f, Nds, s;
1476: DMPlexTransformGetDM(tr, &dm);
1477: DMGetNumFields(dm, &Nf);
1478: for (f = 0; f < Nf; ++f) {
1479: DMLabel label, labelNew;
1480: PetscObject obj;
1481: const char *lname;
1483: DMGetField(rdm, f, &label, &obj);
1484: if (!label) continue;
1485: PetscObjectGetName((PetscObject) label, &lname);
1486: DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
1487: RefineLabel_Internal(tr, label, labelNew);
1488: DMSetField_Internal(rdm, f, labelNew, obj);
1489: DMLabelDestroy(&labelNew);
1490: }
1491: DMGetNumDS(dm, &Nds);
1492: for (s = 0; s < Nds; ++s) {
1493: DMLabel label, labelNew;
1494: const char *lname;
1496: DMGetRegionNumDS(rdm, s, &label, NULL, NULL);
1497: if (!label) continue;
1498: PetscObjectGetName((PetscObject) label, &lname);
1499: DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
1500: RefineLabel_Internal(tr, label, labelNew);
1501: DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL);
1502: DMLabelDestroy(&labelNew);
1503: }
1504: return 0;
1505: }
1507: static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm)
1508: {
1509: DM dm;
1510: PetscSF sf, sfNew;
1511: PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m;
1512: const PetscInt *localPoints;
1513: const PetscSFNode *remotePoints;
1514: PetscInt *localPointsNew;
1515: PetscSFNode *remotePointsNew;
1516: PetscInt pStartNew, pEndNew, pNew;
1517: /* Brute force algorithm */
1518: PetscSF rsf;
1519: PetscSection s;
1520: const PetscInt *rootdegree;
1521: PetscInt *rootPointsNew, *remoteOffsets;
1522: PetscInt numPointsNew, pStart, pEnd, p;
1524: DMPlexTransformGetDM(tr, &dm);
1525: DMPlexGetChart(rdm, &pStartNew, &pEndNew);
1526: DMGetPointSF(dm, &sf);
1527: DMGetPointSF(rdm, &sfNew);
1528: /* Calculate size of new SF */
1529: PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);
1530: if (numRoots < 0) return 0;
1531: for (l = 0; l < numLeaves; ++l) {
1532: const PetscInt p = localPoints[l];
1533: DMPolytopeType ct;
1534: DMPolytopeType *rct;
1535: PetscInt *rsize, *rcone, *rornt;
1536: PetscInt Nct, n;
1538: DMPlexGetCellType(dm, p, &ct);
1539: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1540: for (n = 0; n < Nct; ++n) {
1541: numLeavesNew += rsize[n];
1542: }
1543: }
1544: /* Send new root point numbers
1545: It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication
1546: */
1547: DMPlexGetChart(dm, &pStart, &pEnd);
1548: PetscSectionCreate(PetscObjectComm((PetscObject) dm), &s);
1549: PetscSectionSetChart(s, pStart, pEnd);
1550: for (p = pStart; p < pEnd; ++p) {
1551: DMPolytopeType ct;
1552: DMPolytopeType *rct;
1553: PetscInt *rsize, *rcone, *rornt;
1554: PetscInt Nct, n;
1556: DMPlexGetCellType(dm, p, &ct);
1557: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1558: for (n = 0; n < Nct; ++n) {
1559: PetscSectionAddDof(s, p, rsize[n]);
1560: }
1561: }
1562: PetscSectionSetUp(s);
1563: PetscSectionGetStorageSize(s, &numPointsNew);
1564: PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets);
1565: PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf);
1566: PetscFree(remoteOffsets);
1567: PetscSFComputeDegreeBegin(sf, &rootdegree);
1568: PetscSFComputeDegreeEnd(sf, &rootdegree);
1569: PetscMalloc1(numPointsNew, &rootPointsNew);
1570: for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1;
1571: for (p = pStart; p < pEnd; ++p) {
1572: DMPolytopeType ct;
1573: DMPolytopeType *rct;
1574: PetscInt *rsize, *rcone, *rornt;
1575: PetscInt Nct, n, r, off;
1577: if (!rootdegree[p-pStart]) continue;
1578: PetscSectionGetOffset(s, p, &off);
1579: DMPlexGetCellType(dm, p, &ct);
1580: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1581: for (n = 0, m = 0; n < Nct; ++n) {
1582: for (r = 0; r < rsize[n]; ++r, ++m) {
1583: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1584: rootPointsNew[off+m] = pNew;
1585: }
1586: }
1587: }
1588: PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew,MPI_REPLACE);
1589: PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew,MPI_REPLACE);
1590: PetscSFDestroy(&rsf);
1591: PetscMalloc1(numLeavesNew, &localPointsNew);
1592: PetscMalloc1(numLeavesNew, &remotePointsNew);
1593: for (l = 0, m = 0; l < numLeaves; ++l) {
1594: const PetscInt p = localPoints[l];
1595: DMPolytopeType ct;
1596: DMPolytopeType *rct;
1597: PetscInt *rsize, *rcone, *rornt;
1598: PetscInt Nct, n, r, q, off;
1600: PetscSectionGetOffset(s, p, &off);
1601: DMPlexGetCellType(dm, p, &ct);
1602: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1603: for (n = 0, q = 0; n < Nct; ++n) {
1604: for (r = 0; r < rsize[n]; ++r, ++m, ++q) {
1605: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1606: localPointsNew[m] = pNew;
1607: remotePointsNew[m].index = rootPointsNew[off+q];
1608: remotePointsNew[m].rank = remotePoints[l].rank;
1609: }
1610: }
1611: }
1612: PetscSectionDestroy(&s);
1613: PetscFree(rootPointsNew);
1614: /* SF needs sorted leaves to correctly calculate Gather */
1615: {
1616: PetscSFNode *rp, *rtmp;
1617: PetscInt *lp, *idx, *ltmp, i;
1619: PetscMalloc1(numLeavesNew, &idx);
1620: PetscMalloc1(numLeavesNew, &lp);
1621: PetscMalloc1(numLeavesNew, &rp);
1622: for (i = 0; i < numLeavesNew; ++i) {
1624: idx[i] = i;
1625: }
1626: PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx);
1627: for (i = 0; i < numLeavesNew; ++i) {
1628: lp[i] = localPointsNew[idx[i]];
1629: rp[i] = remotePointsNew[idx[i]];
1630: }
1631: ltmp = localPointsNew;
1632: localPointsNew = lp;
1633: rtmp = remotePointsNew;
1634: remotePointsNew = rp;
1635: PetscFree(idx);
1636: PetscFree(ltmp);
1637: PetscFree(rtmp);
1638: }
1639: PetscSFSetGraph(sfNew, pEndNew-pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
1640: return 0;
1641: }
1643: /*
1644: DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of type ct with localized coordinates x, we generate localized coordinates xr for subcell r of type rct.
1646: Not collective
1648: Input Parameters:
1649: + cr - The DMPlexCellRefiner
1650: . ct - The type of the parent cell
1651: . rct - The type of the produced cell
1652: . r - The index of the produced cell
1653: - x - The localized coordinates for the parent cell
1655: Output Parameter:
1656: . xr - The localized coordinates for the produced cell
1658: Level: developer
1660: .seealso: DMPlexCellRefinerSetCoordinates()
1661: */
1662: static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
1663: {
1664: PetscFE fe = NULL;
1665: PetscInt cdim, v, *subcellV;
1667: DMPlexTransformGetCoordinateFE(tr, ct, &fe);
1668: DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV);
1669: PetscFEGetNumComponents(fe, &cdim);
1670: for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) {
1671: PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v*cdim]);
1672: }
1673: return 0;
1674: }
1676: static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm)
1677: {
1678: DM dm, cdm;
1679: PetscSection coordSection, coordSectionNew;
1680: Vec coordsLocal, coordsLocalNew;
1681: const PetscScalar *coords;
1682: PetscScalar *coordsNew;
1683: const DMBoundaryType *bd;
1684: const PetscReal *maxCell, *L;
1685: PetscBool isperiodic, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE;
1686: PetscInt dE, dEo, d, cStart, cEnd, c, vStartNew, vEndNew, v, pStart, pEnd, p, ocStart, ocEnd;
1688: DMPlexTransformGetDM(tr, &dm);
1689: DMGetCoordinateDM(dm, &cdm);
1690: DMGetPeriodicity(dm, &isperiodic, &maxCell, &L, &bd);
1691: /* Determine if we need to localize coordinates when generating them */
1692: if (isperiodic) {
1693: localizeVertices = PETSC_TRUE;
1694: if (!maxCell) {
1695: PetscBool localized;
1696: DMGetCoordinatesLocalized(dm, &localized);
1698: localizeCells = PETSC_TRUE;
1699: }
1700: }
1702: DMGetCoordinateSection(dm, &coordSection);
1703: PetscSectionGetFieldComponents(coordSection, 0, &dEo);
1704: if (maxCell) {
1705: PetscReal maxCellNew[3];
1707: for (d = 0; d < dEo; ++d) maxCellNew[d] = maxCell[d]/2.0;
1708: DMSetPeriodicity(rdm, isperiodic, maxCellNew, L, bd);
1709: } else {
1710: DMSetPeriodicity(rdm, isperiodic, maxCell, L, bd);
1711: }
1712: DMGetCoordinateDim(rdm, &dE);
1713: PetscSectionCreate(PetscObjectComm((PetscObject) rdm), &coordSectionNew);
1714: PetscSectionSetNumFields(coordSectionNew, 1);
1715: PetscSectionSetFieldComponents(coordSectionNew, 0, dE);
1716: DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew);
1717: if (localizeCells) PetscSectionSetChart(coordSectionNew, 0, vEndNew);
1718: else PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew);
1720: /* Localization should be inherited */
1721: /* Stefano calculates parent cells for each new cell for localization */
1722: /* Localized cells need coordinates of closure */
1723: for (v = vStartNew; v < vEndNew; ++v) {
1724: PetscSectionSetDof(coordSectionNew, v, dE);
1725: PetscSectionSetFieldDof(coordSectionNew, v, 0, dE);
1726: }
1727: if (localizeCells) {
1728: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1729: for (c = cStart; c < cEnd; ++c) {
1730: PetscInt dof;
1732: PetscSectionGetDof(coordSection, c, &dof);
1733: if (dof) {
1734: DMPolytopeType ct;
1735: DMPolytopeType *rct;
1736: PetscInt *rsize, *rcone, *rornt;
1737: PetscInt dim, cNew, Nct, n, r;
1739: DMPlexGetCellType(dm, c, &ct);
1740: dim = DMPolytopeTypeGetDim(ct);
1741: DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1742: /* This allows for different cell types */
1743: for (n = 0; n < Nct; ++n) {
1744: if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
1745: for (r = 0; r < rsize[n]; ++r) {
1746: PetscInt *closure = NULL;
1747: PetscInt clSize, cl, Nv = 0;
1749: DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew);
1750: DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
1751: for (cl = 0; cl < clSize*2; cl += 2) {if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;}
1752: DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
1753: PetscSectionSetDof(coordSectionNew, cNew, Nv * dE);
1754: PetscSectionSetFieldDof(coordSectionNew, cNew, 0, Nv * dE);
1755: }
1756: }
1757: }
1758: }
1759: }
1760: PetscSectionSetUp(coordSectionNew);
1761: DMViewFromOptions(dm, NULL, "-coarse_dm_view");
1762: DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew);
1763: {
1764: VecType vtype;
1765: PetscInt coordSizeNew, bs;
1766: const char *name;
1768: DMGetCoordinatesLocal(dm, &coordsLocal);
1769: VecCreate(PETSC_COMM_SELF, &coordsLocalNew);
1770: PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);
1771: VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE);
1772: PetscObjectGetName((PetscObject) coordsLocal, &name);
1773: PetscObjectSetName((PetscObject) coordsLocalNew, name);
1774: VecGetBlockSize(coordsLocal, &bs);
1775: VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE);
1776: VecGetType(coordsLocal, &vtype);
1777: VecSetType(coordsLocalNew, vtype);
1778: }
1779: VecGetArrayRead(coordsLocal, &coords);
1780: VecGetArray(coordsLocalNew, &coordsNew);
1781: PetscSectionGetChart(coordSection, &ocStart, &ocEnd);
1782: DMPlexGetChart(dm, &pStart, &pEnd);
1783: /* First set coordinates for vertices*/
1784: for (p = pStart; p < pEnd; ++p) {
1785: DMPolytopeType ct;
1786: DMPolytopeType *rct;
1787: PetscInt *rsize, *rcone, *rornt;
1788: PetscInt Nct, n, r;
1789: PetscBool hasVertex = PETSC_FALSE, isLocalized = PETSC_FALSE;
1791: DMPlexGetCellType(dm, p, &ct);
1792: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1793: for (n = 0; n < Nct; ++n) {
1794: if (rct[n] == DM_POLYTOPE_POINT) {hasVertex = PETSC_TRUE; break;}
1795: }
1796: if (localizeVertices && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) {
1797: PetscInt dof;
1798: PetscSectionGetDof(coordSection, p, &dof);
1799: if (dof) isLocalized = PETSC_TRUE;
1800: }
1801: if (hasVertex) {
1802: const PetscScalar *icoords = NULL;
1803: PetscScalar *pcoords = NULL;
1804: PetscInt Nc, Nv, v, d;
1806: DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);
1808: icoords = pcoords;
1809: Nv = Nc/dEo;
1810: if (ct != DM_POLYTOPE_POINT) {
1811: if (localizeVertices) {
1812: PetscScalar anchor[3];
1814: for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d];
1815: if (!isLocalized) {
1816: for (v = 0; v < Nv; ++v) DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v*dEo], &pcoords[v*dEo]);
1817: } else {
1818: Nv = Nc/(2*dEo);
1819: for (v = Nv; v < Nv*2; ++v) DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v*dEo], &pcoords[v*dEo]);
1820: }
1821: }
1822: }
1823: for (n = 0; n < Nct; ++n) {
1824: if (rct[n] != DM_POLYTOPE_POINT) continue;
1825: for (r = 0; r < rsize[n]; ++r) {
1826: PetscScalar vcoords[3];
1827: PetscInt vNew, off;
1829: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew);
1830: PetscSectionGetOffset(coordSectionNew, vNew, &off);
1831: DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords);
1832: DMPlexSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off]);
1833: }
1834: }
1835: DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);
1836: }
1837: }
1838: /* Then set coordinates for cells by localizing */
1839: for (p = pStart; p < pEnd; ++p) {
1840: DMPolytopeType ct;
1841: DMPolytopeType *rct;
1842: PetscInt *rsize, *rcone, *rornt;
1843: PetscInt Nct, n, r;
1844: PetscBool isLocalized = PETSC_FALSE;
1846: DMPlexGetCellType(dm, p, &ct);
1847: DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1848: if (localizeCells && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) {
1849: PetscInt dof;
1850: PetscSectionGetDof(coordSection, p, &dof);
1851: if (dof) isLocalized = PETSC_TRUE;
1852: }
1853: if (isLocalized) {
1854: const PetscScalar *pcoords;
1856: DMPlexPointLocalRead(cdm, p, coords, &pcoords);
1857: for (n = 0; n < Nct; ++n) {
1858: const PetscInt Nr = rsize[n];
1860: if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
1861: for (r = 0; r < Nr; ++r) {
1862: PetscInt pNew, offNew;
1864: /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
1865: DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
1866: cell to the ones it produces. */
1867: DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1868: PetscSectionGetOffset(coordSectionNew, pNew, &offNew);
1869: DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew]);
1870: }
1871: }
1872: }
1873: }
1874: VecRestoreArrayRead(coordsLocal, &coords);
1875: VecRestoreArray(coordsLocalNew, &coordsNew);
1876: DMSetCoordinatesLocal(rdm, coordsLocalNew);
1877: /* TODO Stefano has a final reduction if some hybrid coordinates cannot be found. (needcoords) Should not be needed. */
1878: VecDestroy(&coordsLocalNew);
1879: PetscSectionDestroy(&coordSectionNew);
1880: if (!localizeCells) DMLocalizeCoordinates(rdm);
1881: return 0;
1882: }
1884: PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm)
1885: {
1886: DM rdm;
1887: DMPlexInterpolatedFlag interp;
1892: DMPlexTransformSetDM(tr, dm);
1894: DMCreate(PetscObjectComm((PetscObject)dm), &rdm);
1895: DMSetType(rdm, DMPLEX);
1896: DMPlexTransformSetDimensions(tr, dm, rdm);
1897: /* Calculate number of new points of each depth */
1898: DMPlexIsInterpolated(dm, &interp);
1900: /* Step 1: Set chart */
1901: DMPlexSetChart(rdm, 0, tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]]);
1902: /* Step 2: Set cone/support sizes (automatically stratifies) */
1903: DMPlexTransformSetConeSizes(tr, rdm);
1904: /* Step 3: Setup refined DM */
1905: DMSetUp(rdm);
1906: /* Step 4: Set cones and supports (automatically symmetrizes) */
1907: DMPlexTransformSetCones(tr, rdm);
1908: /* Step 5: Create pointSF */
1909: DMPlexTransformCreateSF(tr, rdm);
1910: /* Step 6: Create labels */
1911: DMPlexTransformCreateLabels(tr, rdm);
1912: /* Step 7: Set coordinates */
1913: DMPlexTransformSetCoordinates(tr, rdm);
1914: DMPlexCopy_Internal(dm, PETSC_TRUE, rdm);
1915: *tdm = rdm;
1916: return 0;
1917: }
1919: PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm)
1920: {
1921: DMPlexTransform tr;
1922: DM cdm, rcdm;
1923: const char *prefix;
1925: DMPlexTransformCreate(PetscObjectComm((PetscObject) dm), &tr);
1926: PetscObjectGetOptionsPrefix((PetscObject) dm, &prefix);
1927: PetscObjectSetOptionsPrefix((PetscObject) tr, prefix);
1928: DMPlexTransformSetDM(tr, dm);
1929: DMPlexTransformSetFromOptions(tr);
1930: DMPlexTransformSetActive(tr, adaptLabel);
1931: DMPlexTransformSetUp(tr);
1932: PetscObjectViewFromOptions((PetscObject) tr, NULL, "-dm_plex_transform_view");
1933: DMPlexTransformApply(tr, dm, rdm);
1934: DMCopyDisc(dm, *rdm);
1935: DMGetCoordinateDM(dm, &cdm);
1936: DMGetCoordinateDM(*rdm, &rcdm);
1937: DMCopyDisc(cdm, rcdm);
1938: DMPlexTransformCreateDiscLabels(tr, *rdm);
1939: DMCopyDisc(dm, *rdm);
1940: DMPlexTransformDestroy(&tr);
1941: ((DM_Plex *) (*rdm)->data)->useHashLocation = ((DM_Plex *) dm->data)->useHashLocation;
1942: return 0;
1943: }