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: }