Actual source code: plextransform.c

  1: #include <petsc/private/dmplextransformimpl.h>

  3: #include <petsc/private/petscfeimpl.h>

  5: PetscClassId DMPLEXTRANSFORM_CLASSID;

  7: PetscFunctionList DMPlexTransformList = NULL;
  8: PetscBool         DMPlexTransformRegisterAllCalled = PETSC_FALSE;

 10: /* Construct cell type order since we must loop over cell types in depth order */
 11: PetscErrorCode DMPlexCreateCellTypeOrder_Internal(DMPolytopeType ctCell, PetscInt *ctOrder[], PetscInt *ctOrderInv[])
 12: {
 13:   PetscInt      *ctO, *ctOInv;
 14:   PetscInt       c, d, off = 0;

 18:   PetscCalloc2(DM_NUM_POLYTOPES+1, &ctO, DM_NUM_POLYTOPES+1, &ctOInv);
 19:   for (d = 3; d >= DMPolytopeTypeGetDim(ctCell); --d) {
 20:     for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
 21:       if (DMPolytopeTypeGetDim((DMPolytopeType) c) != d) continue;
 22:       ctO[off++] = c;
 23:     }
 24:   }
 25:   if (DMPolytopeTypeGetDim(ctCell) != 0) {
 26:     for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
 27:       if (DMPolytopeTypeGetDim((DMPolytopeType) c) != 0) continue;
 28:       ctO[off++] = c;
 29:     }
 30:   }
 31:   for (d = DMPolytopeTypeGetDim(ctCell)-1; d > 0; --d) {
 32:     for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
 33:       if (DMPolytopeTypeGetDim((DMPolytopeType) c) != d) continue;
 34:       ctO[off++] = c;
 35:     }
 36:   }
 37:   for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
 38:     if (DMPolytopeTypeGetDim((DMPolytopeType) c) >= 0) continue;
 39:     ctO[off++] = c;
 40:   }
 41:   if (off != DM_NUM_POLYTOPES+1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid offset %D for cell type order", off);
 42:   for (c = 0; c <= DM_NUM_POLYTOPES; ++c) {
 43:     ctOInv[ctO[c]] = c;
 44:   }
 45:   *ctOrder    = ctO;
 46:   *ctOrderInv = ctOInv;
 47:   return(0);
 48: }

 50: /*@C
 51:   DMPlexTransformRegister - Adds a new transform component implementation

 53:   Not Collective

 55:   Input Parameters:
 56: + name        - The name of a new user-defined creation routine
 57: - create_func - The creation routine itself

 59:   Notes:
 60:   DMPlexTransformRegister() may be called multiple times to add several user-defined transforms

 62:   Sample usage:
 63: .vb
 64:   DMPlexTransformRegister("my_transform", MyTransformCreate);
 65: .ve

 67:   Then, your transform type can be chosen with the procedural interface via
 68: .vb
 69:   DMPlexTransformCreate(MPI_Comm, DMPlexTransform *);
 70:   DMPlexTransformSetType(DMPlexTransform, "my_transform");
 71: .ve
 72:   or at runtime via the option
 73: .vb
 74:   -dm_plex_transform_type my_transform
 75: .ve

 77:   Level: advanced

 79: .seealso: DMPlexTransformRegisterAll(), DMPlexTransformRegisterDestroy()
 80: @*/
 81: PetscErrorCode DMPlexTransformRegister(const char name[], PetscErrorCode (*create_func)(DMPlexTransform))
 82: {

 86:   DMInitializePackage();
 87:   PetscFunctionListAdd(&DMPlexTransformList, name, create_func);
 88:   return(0);
 89: }

 91: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Filter(DMPlexTransform);
 92: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform);
 93: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToBox(DMPlexTransform);
 94: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform);
 95: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform);
 96: PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform);

 98: /*@C
 99:   DMPlexTransformRegisterAll - Registers all of the transform components in the DM package.

101:   Not Collective

103:   Level: advanced

105: .seealso: DMRegisterAll(), DMPlexTransformRegisterDestroy()
106: @*/
107: PetscErrorCode DMPlexTransformRegisterAll(void)
108: {

112:   if (DMPlexTransformRegisterAllCalled) return(0);
113:   DMPlexTransformRegisterAllCalled = PETSC_TRUE;

115:   DMPlexTransformRegister(DMPLEXTRANSFORMFILTER,     DMPlexTransformCreate_Filter);
116:   DMPlexTransformRegister(DMPLEXREFINEREGULAR,       DMPlexTransformCreate_Regular);
117:   DMPlexTransformRegister(DMPLEXREFINETOBOX,         DMPlexTransformCreate_ToBox);
118:   DMPlexTransformRegister(DMPLEXREFINEALFELD,        DMPlexTransformCreate_Alfeld);
119:   DMPlexTransformRegister(DMPLEXREFINEBOUNDARYLAYER, DMPlexTransformCreate_BL);
120:   DMPlexTransformRegister(DMPLEXREFINESBR,           DMPlexTransformCreate_SBR);
121:   return(0);
122: }

124: /*@C
125:   DMPlexTransformRegisterDestroy - This function destroys the . It is called from PetscFinalize().

127:   Level: developer

129: .seealso: PetscInitialize()
130: @*/
131: PetscErrorCode DMPlexTransformRegisterDestroy(void)
132: {

136:   PetscFunctionListDestroy(&DMPlexTransformList);
137:   DMPlexTransformRegisterAllCalled = PETSC_FALSE;
138:   return(0);
139: }

141: /*@
142:   DMPlexTransformCreate - Creates an empty transform object. The type can then be set with DMPlexTransformSetType().

144:   Collective

146:   Input Parameter:
147: . comm - The communicator for the transform object

149:   Output Parameter:
150: . dm - The transform object

152:   Level: beginner

154: .seealso: DMPlexTransformSetType(), DMPLEXREFINEREGULAR, DMPLEXTRANSFORMFILTER
155: @*/
156: PetscErrorCode DMPlexTransformCreate(MPI_Comm comm, DMPlexTransform *tr)
157: {
158:   DMPlexTransform t;
159:   PetscErrorCode  ierr;

163:   *tr = NULL;
164:   DMInitializePackage();

166:   PetscHeaderCreate(t, DMPLEXTRANSFORM_CLASSID, "DMPlexTransform", "Mesh Transform", "DMPlexTransform", comm, DMPlexTransformDestroy, DMPlexTransformView);
167:   t->setupcalled = PETSC_FALSE;
168:   PetscCalloc2(DM_NUM_POLYTOPES, &t->coordFE, DM_NUM_POLYTOPES, &t->refGeom);
169:   *tr = t;
170:   return(0);
171: }

173: /*@C
174:   DMSetType - Sets the particular implementation for a transform.

176:   Collective on tr

178:   Input Parameters:
179: + tr     - The transform
180: - method - The name of the transform type

182:   Options Database Key:
183: . -dm_plex_transform_type <type> - Sets the transform type; use -help for a list of available types

185:   Notes:
186:   See "petsc/include/petscdmplextransform.h" for available transform types

188:   Level: intermediate

190: .seealso: DMPlexTransformGetType(), DMPlexTransformCreate()
191: @*/
192: PetscErrorCode DMPlexTransformSetType(DMPlexTransform tr, DMPlexTransformType method)
193: {
194:   PetscErrorCode (*r)(DMPlexTransform);
195:   PetscBool      match;

200:   PetscObjectTypeCompare((PetscObject) tr, method, &match);
201:   if (match) return(0);

203:   DMPlexTransformRegisterAll();
204:   PetscFunctionListFind(DMPlexTransformList, method, &r);
205:   if (!r) SETERRQ1(PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMPlexTransform type: %s", method);

207:   if (tr->ops->destroy) {(*tr->ops->destroy)(tr);}
208:   PetscMemzero(tr->ops, sizeof(*tr->ops));
209:   PetscObjectChangeTypeName((PetscObject) tr, method);
210:   (*r)(tr);
211:   return(0);
212: }

214: /*@C
215:   DMPlexTransformGetType - Gets the type name (as a string) from the transform.

217:   Not Collective

219:   Input Parameter:
220: . tr  - The DMPlexTransform

222:   Output Parameter:
223: . type - The DMPlexTransform type name

225:   Level: intermediate

227: .seealso: DMPlexTransformSetType(), DMPlexTransformCreate()
228: @*/
229: PetscErrorCode DMPlexTransformGetType(DMPlexTransform tr, DMPlexTransformType *type)
230: {

236:   DMPlexTransformRegisterAll();
237:   *type = ((PetscObject) tr)->type_name;
238:   return(0);
239: }

241: static PetscErrorCode DMPlexTransformView_Ascii(DMPlexTransform tr, PetscViewer v)
242: {
243:   PetscViewerFormat format;
244:   PetscErrorCode    ierr;

247:   PetscViewerGetFormat(v, &format);
248:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
249:     const PetscInt *trTypes = NULL;
250:     IS              trIS;
251:     PetscInt        cols = 8;
252:     PetscInt        Nrt = 8, f, g;

254:     PetscViewerASCIIPrintf(v, "Source Starts\n");
255:     for (g = 0; g <= cols; ++g) {PetscViewerASCIIPrintf(v, " % 14s", DMPolytopeTypes[g]);}
256:     PetscViewerASCIIPrintf(v, "\n");
257:     for (f = 0; f <= cols; ++f) {PetscViewerASCIIPrintf(v, " % 14d", tr->ctStart[f]);}
258:     PetscViewerASCIIPrintf(v, "\n");
259:     PetscViewerASCIIPrintf(v, "Target Starts\n");
260:     for (g = 0; g <= cols; ++g) {PetscViewerASCIIPrintf(v, " % 14s", DMPolytopeTypes[g]);}
261:     PetscViewerASCIIPrintf(v, "\n");
262:     for (f = 0; f <= cols; ++f) {PetscViewerASCIIPrintf(v, " % 14d", tr->ctStartNew[f]);}
263:     PetscViewerASCIIPrintf(v, "\n");

265:     if (tr->trType) {
266:       DMLabelGetNumValues(tr->trType, &Nrt);
267:       DMLabelGetValueIS(tr->trType, &trIS);
268:       ISGetIndices(trIS, &trTypes);
269:     }
270:     PetscViewerASCIIPrintf(v, "Offsets\n");
271:     PetscViewerASCIIPrintf(v, "     ");
272:     for (g = 0; g < cols; ++g) {
273:       PetscViewerASCIIPrintf(v, " % 14s", DMPolytopeTypes[g]);
274:     }
275:     PetscViewerASCIIPrintf(v, "\n");
276:     for (f = 0; f < Nrt; ++f) {
277:       PetscViewerASCIIPrintf(v, "%2d  |", trTypes ? trTypes[f] : f);
278:       for (g = 0; g < cols; ++g) {
279:         PetscViewerASCIIPrintf(v, " % 14D", tr->offset[f*DM_NUM_POLYTOPES+g]);
280:       }
281:       PetscViewerASCIIPrintf(v, " |\n");
282:     }
283:     if (trTypes) {
284:       ISGetIndices(trIS, &trTypes);
285:       ISDestroy(&trIS);
286:     }
287:   }
288:   return(0);
289: }

291: /*@C
292:   DMPlexTransformView - Views a DMPlexTransform

294:   Collective on tr

296:   Input Parameters:
297: + tr - the DMPlexTransform object to view
298: - v  - the viewer

300:   Level: beginner

302: .seealso DMPlexTransformDestroy(), DMPlexTransformCreate()
303: @*/
304: PetscErrorCode DMPlexTransformView(DMPlexTransform tr, PetscViewer v)
305: {
306:   PetscBool      isascii;

311:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) tr), &v);}
314:   PetscViewerCheckWritable(v);
315:   PetscObjectPrintClassNamePrefixType((PetscObject) tr, v);
316:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &isascii);
317:   if (isascii) {DMPlexTransformView_Ascii(tr, v);}
318:   if (tr->ops->view) {(*tr->ops->view)(tr, v);}
319:   return(0);
320: }

322: /*@
323:   DMPlexTransformSetFromOptions - Sets parameters in a transform from the options database

325:   Collective on tr

327:   Input Parameter:
328: . tr - the DMPlexTransform object to set options for

330:   Options Database:
331: . -dm_plex_transform_type - Set the transform type, e.g. refine_regular

333:   Level: intermediate

335: .seealso DMPlexTransformView(), DMPlexTransformCreate()
336: @*/
337: PetscErrorCode DMPlexTransformSetFromOptions(DMPlexTransform tr)
338: {
339:   char           typeName[1024];
340:   const char    *defName;
341:   PetscBool      flg;

346:   DMPlexTransformGetType(tr, &defName);
347:   if (!defName) defName = DMPLEXREFINEREGULAR;
348:   PetscObjectOptionsBegin((PetscObject)tr);
349:   PetscOptionsFList("-dm_plex_transform_type", "DMPlexTransform", "DMPlexTransformSetType", DMPlexTransformList, defName, typeName, 1024, &flg);
350:   if (flg) {DMPlexTransformSetType(tr, typeName);}
351:   else if (!((PetscObject) tr)->type_name) {DMPlexTransformSetType(tr, defName);}
352:   if (tr->ops->setfromoptions) {(*tr->ops->setfromoptions)(PetscOptionsObject,tr);}
353:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
354:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) tr);
355:   PetscOptionsEnd();
356:   return(0);
357: }

359: /*@C
360:   DMPlexTransformDestroy - Destroys a transform.

362:   Collective on tr

364:   Input Parameter:
365: . tr - the transform object to destroy

367:   Level: beginner

369: .seealso DMPlexTransformView(), DMPlexTransformCreate()
370: @*/
371: PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr)
372: {
373:   PetscInt       c;

377:   if (!*tr) return(0);
379:   if (--((PetscObject) (*tr))->refct > 0) {*tr = NULL; return(0);}

381:   if ((*tr)->ops->destroy) {
382:     (*(*tr)->ops->destroy)(*tr);
383:   }
384:   DMDestroy(&(*tr)->dm);
385:   DMLabelDestroy(&(*tr)->active);
386:   DMLabelDestroy(&(*tr)->trType);
387:   PetscFree2((*tr)->ctOrder, (*tr)->ctOrderInv);
388:   PetscFree2((*tr)->ctStart, (*tr)->ctStartNew);
389:   PetscFree((*tr)->offset);
390:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
391:     PetscFEDestroy(&(*tr)->coordFE[c]);
392:     PetscFEGeomDestroy(&(*tr)->refGeom[c]);
393:   }
394:   if ((*tr)->trVerts) {
395:     for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
396:       DMPolytopeType *rct;
397:       PetscInt       *rsize, *rcone, *rornt, Nct, n, r;

399:       if (DMPolytopeTypeGetDim((DMPolytopeType) c) > 0) {
400:         DMPlexTransformCellTransform((*tr), (DMPolytopeType) c, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
401:         for (n = 0; n < Nct; ++n) {

403:           if (rct[n] == DM_POLYTOPE_POINT) continue;
404:           for (r = 0; r < rsize[n]; ++r) {PetscFree((*tr)->trSubVerts[c][rct[n]][r]);}
405:           PetscFree((*tr)->trSubVerts[c][rct[n]]);
406:         }
407:       }
408:       PetscFree((*tr)->trSubVerts[c]);
409:       PetscFree((*tr)->trVerts[c]);
410:     }
411:   }
412:   PetscFree3((*tr)->trNv, (*tr)->trVerts, (*tr)->trSubVerts);
413:   PetscFree2((*tr)->coordFE, (*tr)->refGeom);
414:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
415:   PetscHeaderDestroy(tr);
416:   return(0);
417: }

419: static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrder[], PetscInt ctStart[], PetscInt **offset)
420: {
421:   DMLabel        trType = tr->trType;
422:   PetscInt       c, cN, *off;

426:   if (trType) {
427:     DM              dm;
428:     IS              rtIS;
429:     const PetscInt *reftypes;
430:     PetscInt        Nrt, r;

432:     DMPlexTransformGetDM(tr, &dm);
433:     DMLabelGetNumValues(trType, &Nrt);
434:     DMLabelGetValueIS(trType, &rtIS);
435:     ISGetIndices(rtIS, &reftypes);
436:     PetscCalloc1(Nrt*DM_NUM_POLYTOPES, &off);
437:     for (r = 0; r < Nrt; ++r) {
438:       const PetscInt  rt = reftypes[r];
439:       IS              rtIS;
440:       const PetscInt *points;
441:       DMPolytopeType  ct;
442:       PetscInt        p;

444:       DMLabelGetStratumIS(trType, rt, &rtIS);
445:       ISGetIndices(rtIS, &points);
446:       p    = points[0];
447:       ISRestoreIndices(rtIS, &points);
448:       ISDestroy(&rtIS);
449:       DMPlexGetCellType(dm, p, &ct);
450:       for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
451:         const DMPolytopeType ctNew = (DMPolytopeType) cN;
452:         DMPolytopeType      *rct;
453:         PetscInt            *rsize, *cone, *ornt;
454:         PetscInt             Nct, n, s;

456:         if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {off[r*DM_NUM_POLYTOPES+ctNew] = -1; break;}
457:         off[r*DM_NUM_POLYTOPES+ctNew] = 0;
458:         for (s = 0; s <= r; ++s) {
459:           const PetscInt st = reftypes[s];
460:           DMPolytopeType sct;
461:           PetscInt       q, qrt;

463:           DMLabelGetStratumIS(trType, st, &rtIS);
464:           ISGetIndices(rtIS, &points);
465:           q    = points[0];
466:           ISRestoreIndices(rtIS, &points);
467:           ISDestroy(&rtIS);
468:           DMPlexGetCellType(dm, q, &sct);
469:           DMPlexTransformCellTransform(tr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt);
470:           if (st != qrt) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Refine type %D of point %D does not match predicted type %D", qrt, q, st);
471:           if (st == rt) {
472:             for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) break;
473:             if (n == Nct) off[r*DM_NUM_POLYTOPES+ctNew] = -1;
474:             break;
475:           }
476:           for (n = 0; n < Nct; ++n) {
477:             if (rct[n] == ctNew) {
478:               PetscInt sn;

480:               DMLabelGetStratumSize(trType, st, &sn);
481:               off[r*DM_NUM_POLYTOPES+ctNew] += sn * rsize[n];
482:             }
483:           }
484:         }
485:       }
486:     }
487:     ISRestoreIndices(rtIS, &reftypes);
488:     ISDestroy(&rtIS);
489:   } else {
490:     PetscCalloc1(DM_NUM_POLYTOPES*DM_NUM_POLYTOPES, &off);
491:     for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) {
492:       const DMPolytopeType ct = (DMPolytopeType) c;
493:       for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) {
494:         const DMPolytopeType ctNew = (DMPolytopeType) cN;
495:         DMPolytopeType      *rct;
496:         PetscInt            *rsize, *cone, *ornt;
497:         PetscInt             Nct, n, i;

499:         if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) {off[ct*DM_NUM_POLYTOPES+ctNew] = -1; continue;}
500:         off[ct*DM_NUM_POLYTOPES+ctNew] = 0;
501:         for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) {
502:           const DMPolytopeType ict  = (DMPolytopeType) ctOrder[i];
503:           const DMPolytopeType ictn = (DMPolytopeType) ctOrder[i+1];

505:           DMPlexTransformCellTransform(tr, ict, PETSC_DETERMINE, NULL, &Nct, &rct, &rsize, &cone, &ornt);
506:           if (ict == ct) {
507:             for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) break;
508:             if (n == Nct) off[ct*DM_NUM_POLYTOPES+ctNew] = -1;
509:             break;
510:           }
511:           for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) off[ct*DM_NUM_POLYTOPES+ctNew] += (ctStart[ictn]-ctStart[ict]) * rsize[n];
512:         }
513:       }
514:     }
515:   }
516:   *offset = off;
517:   return(0);
518: }

520: PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr)
521: {
522:   DM             dm;
523:   DMPolytopeType ctCell;
524:   PetscInt       pStart, pEnd, p, c;

529:   if (tr->setupcalled) return(0);
530:   if (tr->ops->setup) {(*tr->ops->setup)(tr);}
531:   DMPlexTransformGetDM(tr, &dm);
532:   DMPlexGetChart(dm, &pStart, &pEnd);
533:   if (pEnd > pStart) {
534:     DMPlexGetCellType(dm, 0, &ctCell);
535:   } else {
536:     PetscInt dim;

538:     DMGetDimension(dm, &dim);
539:     switch (dim) {
540:       case 0: ctCell = DM_POLYTOPE_POINT;break;
541:       case 1: ctCell = DM_POLYTOPE_SEGMENT;break;
542:       case 2: ctCell = DM_POLYTOPE_TRIANGLE;break;
543:       case 3: ctCell = DM_POLYTOPE_TETRAHEDRON;break;
544:       default: ctCell = DM_POLYTOPE_UNKNOWN;
545:     }
546:   }
547:   DMPlexCreateCellTypeOrder_Internal(ctCell, &tr->ctOrder, &tr->ctOrderInv);
548:   /* Construct sizes and offsets for each cell type */
549:   if (!tr->ctStart) {
550:     PetscInt *ctS, *ctSN, *ctC, *ctCN;

552:     PetscCalloc2(DM_NUM_POLYTOPES+1, &ctS, DM_NUM_POLYTOPES+1, &ctSN);
553:     PetscCalloc2(DM_NUM_POLYTOPES+1, &ctC, DM_NUM_POLYTOPES+1, &ctCN);
554:     for (p = pStart; p < pEnd; ++p) {
555:       DMPolytopeType  ct;
556:       DMPolytopeType *rct;
557:       PetscInt       *rsize, *cone, *ornt;
558:       PetscInt        Nct, n;

560:       DMPlexGetCellType(dm, p, &ct);
561:       if (ct == DM_POLYTOPE_UNKNOWN) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %D", p);
562:       ++ctC[ct];
563:       DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt);
564:       for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n];
565:     }
566:     for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
567:       const PetscInt ct  = tr->ctOrder[c];
568:       const PetscInt ctn = tr->ctOrder[c+1];

570:       ctS[ctn]  = ctS[ct]  + ctC[ct];
571:       ctSN[ctn] = ctSN[ct] + ctCN[ct];
572:     }
573:     PetscFree2(ctC, ctCN);
574:     tr->ctStart    = ctS;
575:     tr->ctStartNew = ctSN;
576:   }
577:   DMPlexTransformCreateOffset_Internal(tr, tr->ctOrder, tr->ctStart, &tr->offset);
578:   tr->setupcalled = PETSC_TRUE;
579:   return(0);
580: }

582: PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm)
583: {
587:   *dm = tr->dm;
588:   return(0);
589: }

591: PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm)
592: {

598:   PetscObjectReference((PetscObject) dm);
599:   DMDestroy(&tr->dm);
600:   tr->dm = dm;
601:   return(0);
602: }

604: PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active)
605: {
609:   *active = tr->active;
610:   return(0);
611: }

613: PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active)
614: {

620:   PetscObjectReference((PetscObject) active);
621:   DMLabelDestroy(&tr->active);
622:   tr->active = active;
623:   return(0);
624: }

626: static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe)
627: {

631:   if (!tr->coordFE[ct]) {
632:     PetscInt  dim, cdim;
633:     PetscBool isSimplex;

635:     switch (ct) {
636:       case DM_POLYTOPE_SEGMENT:       dim = 1; isSimplex = PETSC_TRUE;  break;
637:       case DM_POLYTOPE_TRIANGLE:      dim = 2; isSimplex = PETSC_TRUE;  break;
638:       case DM_POLYTOPE_QUADRILATERAL: dim = 2; isSimplex = PETSC_FALSE; break;
639:       case DM_POLYTOPE_TETRAHEDRON:   dim = 3; isSimplex = PETSC_TRUE;  break;
640:       case DM_POLYTOPE_HEXAHEDRON:    dim = 3; isSimplex = PETSC_FALSE; break;
641:       default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No coordinate FE for cell type %s", DMPolytopeTypes[ct]);
642:     }
643:     DMGetCoordinateDim(tr->dm, &cdim);
644:     PetscFECreateLagrange(PETSC_COMM_SELF, dim, cdim, isSimplex, 1, PETSC_DETERMINE, &tr->coordFE[ct]);
645:     {
646:       PetscDualSpace  dsp;
647:       PetscQuadrature quad;
648:       DM              K;
649:       PetscFEGeom    *cg;
650:       PetscScalar    *Xq;
651:       PetscReal      *xq, *wq;
652:       PetscInt        Nq, q;

654:       DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq);
655:       PetscMalloc1(Nq*cdim, &xq);
656:       for (q = 0; q < Nq*cdim; ++q) xq[q] = PetscRealPart(Xq[q]);
657:       PetscMalloc1(Nq, &wq);
658:       for (q = 0; q < Nq; ++q) wq[q] = 1.0;
659:       PetscQuadratureCreate(PETSC_COMM_SELF, &quad);
660:       PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq);
661:       PetscFESetQuadrature(tr->coordFE[ct], quad);

663:       PetscFEGetDualSpace(tr->coordFE[ct], &dsp);
664:       PetscDualSpaceGetDM(dsp, &K);
665:       PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &tr->refGeom[ct]);
666:       cg   = tr->refGeom[ct];
667:       DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
668:       PetscQuadratureDestroy(&quad);
669:     }
670:   }
671:   *fe = tr->coordFE[ct];
672:   return(0);
673: }

675: /*@
676:   DMPlexTransformGetTargetPoint - Get the number of a point in the transformed mesh based on information from the original mesh.

678:   Not collective

680:   Input Parameters:
681: + tr    - The DMPlexTransform
682: . ct    - The type of the original point which produces the new point
683: . ctNew - The type of the new point
684: . p     - The original point which produces the new point
685: - r     - The replica number of the new point, meaning it is the rth point of type ctNew produced from p

687:   Output Parameters:
688: . pNew  - The new point number

690:   Level: developer

692: .seealso: DMPlexTransformGetSourcePoint(), DMPlexTransformCellTransform()
693: @*/
694: PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew)
695: {
696:   DMPolytopeType *rct;
697:   PetscInt       *rsize, *cone, *ornt;
698:   PetscInt       rt, Nct, n, off, rp;
699:   DMLabel        trType = tr->trType;
700:   PetscInt       ctS  = tr->ctStart[ct],       ctE  = tr->ctStart[tr->ctOrder[tr->ctOrderInv[ct]+1]];
701:   PetscInt       ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrder[tr->ctOrderInv[ctNew]+1]];
702:   PetscInt       newp = ctSN, cind;

706:   if ((p < ctS) || (p >= ctE)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %D is not a %s [%D, %D)", p, DMPolytopeTypes[ct], ctS, ctE);
707:   DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt);
708:   if (trType) {
709:     DMLabelGetValueIndex(trType, rt, &cind);
710:     DMLabelGetStratumPointIndex(trType, rt, p, &rp);
711:     if (rp < 0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s point %D does not have refine type %D", DMPolytopeTypes[ct], p, rt);
712:   } else {
713:     cind = ct;
714:     rp   = p - ctS;
715:   }
716:   off = tr->offset[cind*DM_NUM_POLYTOPES + ctNew];
717:   if (off < 0) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s (%D) of point %D does not produce type %s for transform %s", DMPolytopeTypes[ct], rt, p, DMPolytopeTypes[ctNew], tr->hdr.type_name);
718:   newp += off;
719:   for (n = 0; n < Nct; ++n) {
720:     if (rct[n] == ctNew) {
721:       if (rsize[n] && r >= rsize[n])
722:         SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number %D should be in [0, %D) for subcell type %s in cell type %s", r, rsize[n], DMPolytopeTypes[rct[n]], DMPolytopeTypes[ct]);
723:       newp += rp * rsize[n] + r;
724:       break;
725:     }
726:   }

728:   if ((newp < ctSN) || (newp >= ctEN)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "New point %D is not a %s [%D, %D)", newp, DMPolytopeTypes[ctNew], ctSN, ctEN);
729:   *pNew = newp;
730:   return(0);
731: }

733: PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r)
734: {
735:   DMLabel         trType = tr->trType;
736:   DMPolytopeType *rct;
737:   PetscInt       *rsize, *cone, *ornt;
738:   PetscInt        rt, Nct, n, rp = 0, rO = 0, pO;
739:   PetscInt        offset = -1, ctS, ctE, ctO, ctN, ctTmp;
740:   PetscErrorCode  ierr;

743:   for (ctN = 0; ctN < DM_NUM_POLYTOPES; ++ctN) {
744:     PetscInt ctSN = tr->ctStartNew[ctN], ctEN = tr->ctStartNew[tr->ctOrder[tr->ctOrderInv[ctN]+1]];

746:     if ((pNew >= ctSN) && (pNew < ctEN)) break;
747:   }
748:   if (ctN == DM_NUM_POLYTOPES) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cell type for target point %D could be not found", pNew);
749:   if (trType) {
750:     DM              dm;
751:     IS              rtIS;
752:     const PetscInt *reftypes;
753:     PetscInt        Nrt, r, rtStart;

755:     DMPlexTransformGetDM(tr, &dm);
756:     DMLabelGetNumValues(trType, &Nrt);
757:     DMLabelGetValueIS(trType, &rtIS);
758:     ISGetIndices(rtIS, &reftypes);
759:     for (r = 0; r < Nrt; ++r) {
760:       const PetscInt off = tr->offset[r*DM_NUM_POLYTOPES + ctN];

762:       if (tr->ctStartNew[ctN] + off > pNew) continue;
763:       /* Check that any of this refinement type exist */
764:       /* TODO Actually keep track of the number produced here instead */
765:       if (off > offset) {rt = reftypes[r]; offset = off;}
766:     }
767:     ISRestoreIndices(rtIS, &reftypes);
768:     ISDestroy(&rtIS);
769:     if (offset < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %D could be not found", pNew);
770:     /* TODO Map refinement types to cell types */
771:     DMLabelGetStratumBounds(trType, rt, &rtStart, NULL);
772:     if (rtStart < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Refinement type %D has no source points", rt);
773:     for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) {
774:       PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrder[tr->ctOrderInv[ctO]+1]];

776:       if ((rtStart >= ctS) && (rtStart < ctE)) break;
777:     }
778:     if (ctO == DM_NUM_POLYTOPES) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine a cell type for refinement type %D", rt);
779:   } else {
780:     for (ctTmp = 0; ctTmp < DM_NUM_POLYTOPES; ++ctTmp) {
781:       const PetscInt off = tr->offset[ctTmp*DM_NUM_POLYTOPES + ctN];

783:       if (tr->ctStartNew[ctN] + off > pNew) continue;
784:       if (tr->ctStart[tr->ctOrder[tr->ctOrderInv[ctTmp]+1]] <= tr->ctStart[ctTmp]) continue;
785:       /* TODO Actually keep track of the number produced here instead */
786:       if (off > offset) {ctO = ctTmp; offset = off;}
787:     }
788:     if (offset < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %D could be not found", pNew);
789:   }
790:   ctS = tr->ctStart[ctO];
791:   ctE = tr->ctStart[tr->ctOrder[tr->ctOrderInv[ctO]+1]];
792:   DMPlexTransformCellTransform(tr, (DMPolytopeType) ctO, ctS, &rt, &Nct, &rct, &rsize, &cone, &ornt);
793:   for (n = 0; n < Nct; ++n) {
794:     if (rct[n] == ctN) {
795:       PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset;

797:       rp = tmp / rsize[n];
798:       rO = tmp % rsize[n];
799:       break;
800:     }
801:   }
802:   if (n == Nct) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number for target point %D could be not found", pNew);
803:   pO = rp + ctS;
804:   if ((pO < ctS) || (pO >= ctE)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Source point %D is not a %s [%D, %D)", pO, DMPolytopeTypes[ctO], ctS, ctE);
805:   if (ct)    *ct    = (DMPolytopeType) ctO;
806:   if (ctNew) *ctNew = (DMPolytopeType) ctN;
807:   if (p)     *p     = pO;
808:   if (r)     *r     = rO;
809:   return(0);
810: }

812: /*@
813:   DMPlexTransformCellTransform - Describes the transform of a given source cell into a set of other target cells. These produced cells become the new mesh.

815:   Input Parameters:
816: + tr     - The DMPlexTransform object
817: . source - The source cell type
818: - p      - The source point, which can also determine the refine type

820:   Output Parameters:
821: + rt     - The refine type for this point
822: . Nt     - The number of types produced by this point
823: . target - An array of length Nt giving the types produced
824: . size   - An array of length Nt giving the number of cells of each type produced
825: . cone   - An array of length Nt*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point
826: - ornt   - An array of length Nt*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point

828:   Note: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we
829:   need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical
830:   division (described in these outputs) of the cell in the original mesh. The point identifier is given by
831: $   the number of cones to be taken, or 0 for the current cell
832: $   the cell cone point number at each level from which it is subdivided
833: $   the replica number r of the subdivision.
834: The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is
835: $   Nt     = 2
836: $   target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT}
837: $   size   = {1, 2}
838: $   cone   = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0,  DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}
839: $   ornt   = {                         0,                       0,                        0,                          0}

841:   Level: advanced

843: .seealso: DMPlexTransformApply(), DMPlexTransformCreate()
844: @*/
845: PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
846: {

850:   (*tr->ops->celltransform)(tr, source, p, rt, Nt, target, size, cone, ornt);
851:   return(0);
852: }

854: PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
855: {
857:   *rnew = r;
858:   *onew = DMPolytopeTypeComposeOrientation(tct, o, so);
859:   return(0);
860: }

862: /* Returns the same thing */
863: PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[])
864: {
865:   static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT};
866:   static PetscInt       vertexS[] = {1};
867:   static PetscInt       vertexC[] = {0};
868:   static PetscInt       vertexO[] = {0};
869:   static DMPolytopeType edgeT[]   = {DM_POLYTOPE_SEGMENT};
870:   static PetscInt       edgeS[]   = {1};
871:   static PetscInt       edgeC[]   = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
872:   static PetscInt       edgeO[]   = {0, 0};
873:   static DMPolytopeType tedgeT[]  = {DM_POLYTOPE_POINT_PRISM_TENSOR};
874:   static PetscInt       tedgeS[]  = {1};
875:   static PetscInt       tedgeC[]  = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0};
876:   static PetscInt       tedgeO[]  = {0, 0};
877:   static DMPolytopeType triT[]    = {DM_POLYTOPE_TRIANGLE};
878:   static PetscInt       triS[]    = {1};
879:   static PetscInt       triC[]    = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0};
880:   static PetscInt       triO[]    = {0, 0, 0};
881:   static DMPolytopeType quadT[]   = {DM_POLYTOPE_QUADRILATERAL};
882:   static PetscInt       quadS[]   = {1};
883:   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};
884:   static PetscInt       quadO[]   = {0, 0, 0, 0};
885:   static DMPolytopeType tquadT[]  = {DM_POLYTOPE_SEG_PRISM_TENSOR};
886:   static PetscInt       tquadS[]  = {1};
887:   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};
888:   static PetscInt       tquadO[]  = {0, 0, 0, 0};
889:   static DMPolytopeType tetT[]    = {DM_POLYTOPE_TETRAHEDRON};
890:   static PetscInt       tetS[]    = {1};
891:   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};
892:   static PetscInt       tetO[]    = {0, 0, 0, 0};
893:   static DMPolytopeType hexT[]    = {DM_POLYTOPE_HEXAHEDRON};
894:   static PetscInt       hexS[]    = {1};
895:   static PetscInt       hexC[]    = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0,
896:                                      DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 1, 5, 0};
897:   static PetscInt       hexO[]    = {0, 0, 0, 0, 0, 0};
898:   static DMPolytopeType tripT[]   = {DM_POLYTOPE_TRI_PRISM};
899:   static PetscInt       tripS[]   = {1};
900:   static PetscInt       tripC[]   = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0,
901:                                      DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0};
902:   static PetscInt       tripO[]   = {0, 0, 0, 0, 0};
903:   static DMPolytopeType ttripT[]  = {DM_POLYTOPE_TRI_PRISM_TENSOR};
904:   static PetscInt       ttripS[]  = {1};
905:   static PetscInt       ttripC[]  = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0,
906:                                      DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0};
907:   static PetscInt       ttripO[]  = {0, 0, 0, 0, 0};
908:   static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR};
909:   static PetscInt       tquadpS[] = {1};
910:   static PetscInt       tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0,
911:                                      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};
912:   static PetscInt       tquadpO[] = {0, 0, 0, 0, 0, 0};
913:   static DMPolytopeType pyrT[]    = {DM_POLYTOPE_PYRAMID};
914:   static PetscInt       pyrS[]    = {1};
915:   static PetscInt       pyrC[]    = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0,
916:                                      DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 1, 4, 0};
917:   static PetscInt       pyrO[]    = {0, 0, 0, 0, 0};

920:   if (rt) *rt = 0;
921:   switch (source) {
922:     case DM_POLYTOPE_POINT:              *Nt = 1; *target = vertexT; *size = vertexS; *cone = vertexC; *ornt = vertexO; break;
923:     case DM_POLYTOPE_SEGMENT:            *Nt = 1; *target = edgeT;   *size = edgeS;   *cone = edgeC;   *ornt = edgeO;   break;
924:     case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT;  *size = tedgeS;  *cone = tedgeC;  *ornt = tedgeO;  break;
925:     case DM_POLYTOPE_TRIANGLE:           *Nt = 1; *target = triT;    *size = triS;    *cone = triC;    *ornt = triO;    break;
926:     case DM_POLYTOPE_QUADRILATERAL:      *Nt = 1; *target = quadT;   *size = quadS;   *cone = quadC;   *ornt = quadO;   break;
927:     case DM_POLYTOPE_SEG_PRISM_TENSOR:   *Nt = 1; *target = tquadT;  *size = tquadS;  *cone = tquadC;  *ornt = tquadO;  break;
928:     case DM_POLYTOPE_TETRAHEDRON:        *Nt = 1; *target = tetT;    *size = tetS;    *cone = tetC;    *ornt = tetO;    break;
929:     case DM_POLYTOPE_HEXAHEDRON:         *Nt = 1; *target = hexT;    *size = hexS;    *cone = hexC;    *ornt = hexO;    break;
930:     case DM_POLYTOPE_TRI_PRISM:          *Nt = 1; *target = tripT;   *size = tripS;   *cone = tripC;   *ornt = tripO;   break;
931:     case DM_POLYTOPE_TRI_PRISM_TENSOR:   *Nt = 1; *target = ttripT;  *size = ttripS;  *cone = ttripC;  *ornt = ttripO;  break;
932:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:  *Nt = 1; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break;
933:     case DM_POLYTOPE_PYRAMID:            *Nt = 1; *target = pyrT;    *size = pyrS;    *cone = pyrC;    *ornt = pyrO;    break;
934:     default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]);
935:   }
936:   return(0);
937: }

939: /*@
940:   DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point

942:   Not collective

944:   Input Parameters:
945: + tr  - The DMPlexTransform
946: . sct - The source point cell type, from whom the new cell is being produced
947: . sp  - The source point
948: . so  - The orientation of the source point in its enclosing parent
949: . tct - The target point cell type
950: . r   - The replica number requested for the produced cell type
951: - o   - The orientation of the replica

953:   Output Parameters:
954: + rnew - The replica number, given the orientation of the parent
955: - onew - The replica orientation, given the orientation of the parent

957:   Level: advanced

959: .seealso: DMPlexTransformCellTransform(), DMPlexTransformApply()
960: @*/
961: PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew)
962: {

966:   (*tr->ops->getsubcellorientation)(tr, sct, sp, so, tct, r, o, rnew, onew);
967:   return(0);
968: }

970: static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm)
971: {
972:   DM              dm;
973:   PetscInt        pStart, pEnd, p, pNew;
974:   PetscErrorCode  ierr;

977:   DMPlexTransformGetDM(tr, &dm);
978:   /* Must create the celltype label here so that we do not automatically try to compute the types */
979:   DMCreateLabel(rdm, "celltype");
980:   DMPlexGetChart(dm, &pStart, &pEnd);
981:   for (p = pStart; p < pEnd; ++p) {
982:     DMPolytopeType  ct;
983:     DMPolytopeType *rct;
984:     PetscInt       *rsize, *rcone, *rornt;
985:     PetscInt        Nct, n, r;

987:     DMPlexGetCellType(dm, p, &ct);
988:     DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
989:     for (n = 0; n < Nct; ++n) {
990:       for (r = 0; r < rsize[n]; ++r) {
991:         DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
992:         DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]));
993:         DMPlexSetCellType(rdm, pNew, rct[n]);
994:       }
995:     }
996:   }
997:   /* Let the DM know we have set all the cell types */
998:   {
999:     DMLabel  ctLabel;
1000:     DM_Plex *plex = (DM_Plex *) rdm->data;

1002:     DMPlexGetCellTypeLabel(rdm, &ctLabel);
1003:     PetscObjectStateGet((PetscObject) ctLabel, &plex->celltypeState);
1004:   }
1005:   return(0);
1006: }

1008: #if 0
1009: static PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize)
1010: {
1011:   PetscInt ctNew;

1016:   /* TODO Can do bisection since everything is sorted */
1017:   for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) {
1018:     PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrder[tr->ctOrderInv[ctNew]+1]];

1020:     if (q >= ctSN && q < ctEN) break;
1021:   }
1022:   if (ctNew >= DM_NUM_POLYTOPES) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %D cannot be located in the transformed mesh", q);
1023:   *coneSize = DMPolytopeTypeGetConeSize((DMPolytopeType) ctNew);
1024:   return(0);
1025: }
1026: #endif

1028: /* The orientation o is for the interior of the cell p */
1029: static PetscErrorCode DMPlexTransformGetCone_Internal(DMPlexTransform tr, PetscInt p, PetscInt o, DMPolytopeType ct, DMPolytopeType ctNew,
1030:                                                       const PetscInt rcone[], PetscInt *coneoff, const PetscInt rornt[], PetscInt *orntoff,
1031:                                                       PetscInt coneNew[], PetscInt orntNew[])
1032: {
1033:   DM              dm;
1034:   const PetscInt  csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1035:   const PetscInt *cone;
1036:   PetscInt        c, coff = *coneoff, ooff = *orntoff;
1037:   PetscErrorCode  ierr;

1040:   DMPlexTransformGetDM(tr, &dm);
1041:   DMPlexGetCone(dm, p, &cone);
1042:   for (c = 0; c < csizeNew; ++c) {
1043:     PetscInt             ppp   = -1;                             /* Parent Parent point: Parent of point pp */
1044:     PetscInt             pp    = p;                              /* Parent point: Point in the original mesh producing new cone point */
1045:     PetscInt             po    = 0;                              /* Orientation of parent point pp in parent parent point ppp */
1046:     DMPolytopeType       pct   = ct;                             /* Parent type: Cell type for parent of new cone point */
1047:     const PetscInt      *pcone = cone;                           /* Parent cone: Cone of parent point pp */
1048:     PetscInt             pr    = -1;                             /* Replica number of pp that produces new cone point  */
1049:     const DMPolytopeType ft    = (DMPolytopeType) rcone[coff++]; /* Cell type for new cone point of pNew */
1050:     const PetscInt       fn    = rcone[coff++];                  /* Number of cones of p that need to be taken when producing new cone point */
1051:     PetscInt             fo    = rornt[ooff++];                  /* Orientation of new cone point in pNew */
1052:     PetscInt             lc;

1054:     /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */
1055:     for (lc = 0; lc < fn; ++lc) {
1056:       const PetscInt *parr = DMPolytopeTypeGetArrangment(pct, po);
1057:       const PetscInt  acp  = rcone[coff++];
1058:       const PetscInt  pcp  = parr[acp*2];
1059:       const PetscInt  pco  = parr[acp*2+1];
1060:       const PetscInt *ppornt;

1062:       ppp  = pp;
1063:       pp   = pcone[pcp];
1064:       DMPlexGetCellType(dm, pp, &pct);
1065:       DMPlexGetCone(dm, pp, &pcone);
1066:       DMPlexGetConeOrientation(dm, ppp, &ppornt);
1067:       po   = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco);
1068:     }
1069:     pr = rcone[coff++];
1070:     /* Orientation po of pp maps (pr, fo) -> (pr', fo') */
1071:     DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo);
1072:     DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c]);
1073:     orntNew[c] = fo;
1074:   }
1075:   *coneoff = coff;
1076:   *orntoff = ooff;
1077:   return(0);
1078: }

1080: static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm)
1081: {
1082:   DM             dm;
1083:   DMPolytopeType ct;
1084:   PetscInt      *coneNew, *orntNew;
1085:   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;
1131:   PetscErrorCode  ierr;

1137:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p));
1138:   DMPlexTransformGetDM(tr, &dm);
1139:   DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone);
1140:   DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt);
1141:   DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r);
1142:   DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1143:   for (n = 0; n < Nct; ++n) {
1144:     const DMPolytopeType ctNew    = rct[n];
1145:     const PetscInt       csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1146:     PetscInt             Nr = rsize[n], fn, c;

1148:     if (ctNew == qct) Nr = r;
1149:     for (nr = 0; nr < Nr; ++nr) {
1150:       for (c = 0; c < csizeNew; ++c) {
1151:         ++coff;             /* Cell type of new cone point */
1152:         fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1153:         coff += fn;
1154:         ++coff;             /* Replica number of new cone point */
1155:         ++ooff;             /* Orientation of new cone point */
1156:       }
1157:     }
1158:     if (ctNew == qct) break;
1159:   }
1160:   DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt);
1161:   *cone = qcone;
1162:   *ornt = qornt;
1163:   return(0);
1164: }

1166: PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1167: {
1168:   DM              dm;
1169:   DMPolytopeType  ct, qct;
1170:   DMPolytopeType *rct;
1171:   PetscInt       *rsize, *rcone, *rornt, *qcone, *qornt;
1172:   PetscInt        maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0;
1173:   PetscErrorCode  ierr;

1178:   for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p));
1179:   DMPlexTransformGetDM(tr, &dm);
1180:   DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone);
1181:   DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt);
1182:   DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r);
1183:   DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1184:   for (n = 0; n < Nct; ++n) {
1185:     const DMPolytopeType ctNew    = rct[n];
1186:     const PetscInt       csizeNew = DMPolytopeTypeGetConeSize(ctNew);
1187:     PetscInt             Nr = rsize[n], fn, c;

1189:     if (ctNew == qct) Nr = r;
1190:     for (nr = 0; nr < Nr; ++nr) {
1191:       for (c = 0; c < csizeNew; ++c) {
1192:         ++coff;             /* Cell type of new cone point */
1193:         fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */
1194:         coff += fn;
1195:         ++coff;             /* Replica number of new cone point */
1196:         ++ooff;             /* Orientation of new cone point */
1197:       }
1198:     }
1199:     if (ctNew == qct) break;
1200:   }
1201:   DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt);
1202:   *cone = qcone;
1203:   *ornt = qornt;
1204:   return(0);
1205: }

1207: PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[])
1208: {
1209:   DM             dm;

1215:   DMPlexTransformGetDM(tr, &dm);
1216:   DMRestoreWorkArray(dm, 0, MPIU_INT, cone);
1217:   DMRestoreWorkArray(dm, 0, MPIU_INT, ornt);
1218:   return(0);
1219: }

1221: static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr)
1222: {
1223:   PetscInt       ict;

1227:   PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts);
1228:   for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) {
1229:     const DMPolytopeType ct = (DMPolytopeType) ict;
1230:     DMPlexTransform    reftr;
1231:     DM                 refdm, trdm;
1232:     Vec                coordinates;
1233:     const PetscScalar *coords;
1234:     DMPolytopeType    *rct;
1235:     PetscInt          *rsize, *rcone, *rornt;
1236:     PetscInt           Nct, n, r, pNew;
1237:     PetscInt           vStart, vEnd, Nc;
1238:     const PetscInt     debug = 0;
1239:     const char        *typeName;

1241:     /* Since points are 0-dimensional, coordinates make no sense */
1242:     if (DMPolytopeTypeGetDim(ct) <= 0) continue;
1243:     DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm);
1244:     DMPlexTransformCreate(PETSC_COMM_SELF, &reftr);
1245:     DMPlexTransformSetDM(reftr, refdm);
1246:     DMPlexTransformGetType(tr, &typeName);
1247:     DMPlexTransformSetType(reftr, typeName);
1248:     DMPlexTransformSetUp(reftr);
1249:     DMPlexTransformApply(reftr, refdm, &trdm);

1251:     DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd);
1252:     tr->trNv[ct] = vEnd - vStart;
1253:     DMGetCoordinatesLocal(trdm, &coordinates);
1254:     VecGetLocalSize(coordinates, &Nc);
1255:     if (tr->trNv[ct] * DMPolytopeTypeGetDim(ct) != Nc) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell type %s, transformed coordinate size %D != %D size of coordinate storage", DMPolytopeTypes[ct], tr->trNv[ct] * DMPolytopeTypeGetDim(ct), Nc);
1256:     PetscCalloc1(Nc, &tr->trVerts[ct]);
1257:     VecGetArrayRead(coordinates, &coords);
1258:     PetscArraycpy(tr->trVerts[ct], coords, Nc);
1259:     VecRestoreArrayRead(coordinates, &coords);

1261:     PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct]);
1262:     DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1263:     for (n = 0; n < Nct; ++n) {

1265:       /* Since points are 0-dimensional, coordinates make no sense */
1266:       if (rct[n] == DM_POLYTOPE_POINT) continue;
1267:       PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]]);
1268:       for (r = 0; r < rsize[n]; ++r) {
1269:         PetscInt *closure = NULL;
1270:         PetscInt  clSize, cl, Nv = 0;

1272:         PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r]);
1273:         DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew);
1274:         DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure);
1275:         for (cl = 0; cl < clSize*2; cl += 2) {
1276:           const PetscInt sv = closure[cl];

1278:           if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart;
1279:         }
1280:         DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure);
1281:         if (Nv != DMPolytopeTypeGetNumVertices(rct[n])) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of vertices %D != %D for %s subcell %D from cell %s", Nv, DMPolytopeTypeGetNumVertices(rct[n]), DMPolytopeTypes[rct[n]], r, DMPolytopeTypes[ct]);
1282:       }
1283:     }
1284:     if (debug) {
1285:       DMPolytopeType *rct;
1286:       PetscInt       *rsize, *rcone, *rornt;
1287:       PetscInt        v, dE = DMPolytopeTypeGetDim(ct), d, off = 0;

1289:       PetscPrintf(PETSC_COMM_SELF, "%s: %D vertices\n", DMPolytopeTypes[ct], tr->trNv[ct]);
1290:       for (v = 0; v < tr->trNv[ct]; ++v) {
1291:         PetscPrintf(PETSC_COMM_SELF, "  ");
1292:         for (d = 0; d < dE; ++d) {PetscPrintf(PETSC_COMM_SELF, "%g ", tr->trVerts[ct][off++]);}
1293:         PetscPrintf(PETSC_COMM_SELF, "\n");
1294:       }

1296:       DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1297:       for (n = 0; n < Nct; ++n) {
1298:         if (rct[n] == DM_POLYTOPE_POINT) continue;
1299:         PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct]);
1300:         for (r = 0; r < rsize[n]; ++r) {
1301:           PetscPrintf(PETSC_COMM_SELF, "  ");
1302:           for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) {
1303:             PetscPrintf(PETSC_COMM_SELF, "%D ", tr->trSubVerts[ct][rct[n]][r][v]);
1304:           }
1305:           PetscPrintf(PETSC_COMM_SELF, "\n");
1306:         }
1307:       }
1308:     }
1309:     DMDestroy(&refdm);
1310:     DMDestroy(&trdm);
1311:     DMPlexTransformDestroy(&reftr);
1312:   }
1313:   return(0);
1314: }

1316: /*
1317:   DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type

1319:   Input Parameters:
1320: + tr - The DMPlexTransform object
1321: - ct - The cell type

1323:   Output Parameters:
1324: + Nv      - The number of transformed vertices in the closure of the reference cell of given type
1325: - trVerts - The coordinates of these vertices in the reference cell

1327:   Level: developer

1329: .seealso: DMPlexTransformGetSubcellVertices()
1330: */
1331: PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[])
1332: {

1336:   if (!tr->trNv) {DMPlexTransformCreateCellVertices_Internal(tr);}
1337:   if (Nv)      *Nv      = tr->trNv[ct];
1338:   if (trVerts) *trVerts = tr->trVerts[ct];
1339:   return(0);
1340: }

1342: /*
1343:   DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type

1345:   Input Parameters:
1346: + tr  - The DMPlexTransform object
1347: . ct  - The cell type
1348: . rct - The subcell type
1349: - r   - The subcell index

1351:   Output Parameter:
1352: . subVerts - The indices of these vertices in the set of vertices returned by DMPlexTransformGetCellVertices()

1354:   Level: developer

1356: .seealso: DMPlexTransformGetCellVertices()
1357: */
1358: PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[])
1359: {

1363:   if (!tr->trNv) {DMPlexTransformCreateCellVertices_Internal(tr);}
1364:   if (!tr->trSubVerts[ct][rct]) SETERRQ2(PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]);
1365:   if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r];
1366:   return(0);
1367: }

1369: /* Computes new vertex as the barycenter, or centroid */
1370: PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1371: {
1372:   PetscInt v,d;

1375:   if (ct != DM_POLYTOPE_POINT) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for refined point type %s",DMPolytopeTypes[ct]);
1376:   for (d = 0; d < dE; ++d) out[d] = 0.0;
1377:   for (v = 0; v < Nv; ++v) for (d = 0; d < dE; ++d) out[d] += in[v*dE+d];
1378:   for (d = 0; d < dE; ++d) out[d] /= Nv;
1379:   return(0);
1380: }

1382: /*@
1383:   DMPlexTransformMapCoordinates -

1385:   Input Parameters:
1386: + tr   - The DMPlexTransform
1387: . pct  - The cell type of the parent, from whom the new cell is being produced
1388: . ct   - The type being produced
1389: . r    - The replica number requested for the produced cell type
1390: . Nv   - Number of vertices in the closure of the parent cell
1391: . dE   - Spatial dimension
1392: - in   - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell

1394:   Output Parameter:
1395: . out - The coordinates of the new vertices
1396: @*/
1397: PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[])
1398: {

1402:   (*tr->ops->mapcoordinates)(tr, pct, ct, r, Nv, dE, in, out);
1403:   return(0);
1404: }

1406: static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew)
1407: {
1408:   DM              dm;
1409:   IS              valueIS;
1410:   const PetscInt *values;
1411:   PetscInt        defVal, Nv, val;
1412:   PetscErrorCode  ierr;

1415:   DMPlexTransformGetDM(tr, &dm);
1416:   DMLabelGetDefaultValue(label, &defVal);
1417:   DMLabelSetDefaultValue(labelNew, defVal);
1418:   DMLabelGetValueIS(label, &valueIS);
1419:   ISGetLocalSize(valueIS, &Nv);
1420:   ISGetIndices(valueIS, &values);
1421:   for (val = 0; val < Nv; ++val) {
1422:     IS              pointIS;
1423:     const PetscInt *points;
1424:     PetscInt        numPoints, p;

1426:     /* Ensure refined label is created with same number of strata as
1427:      * original (even if no entries here). */
1428:     DMLabelAddStratum(labelNew, values[val]);
1429:     DMLabelGetStratumIS(label, values[val], &pointIS);
1430:     ISGetLocalSize(pointIS, &numPoints);
1431:     ISGetIndices(pointIS, &points);
1432:     for (p = 0; p < numPoints; ++p) {
1433:       const PetscInt  point = points[p];
1434:       DMPolytopeType  ct;
1435:       DMPolytopeType *rct;
1436:       PetscInt       *rsize, *rcone, *rornt;
1437:       PetscInt        Nct, n, r, pNew=0;

1439:       DMPlexGetCellType(dm, point, &ct);
1440:       DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1441:       for (n = 0; n < Nct; ++n) {
1442:         for (r = 0; r < rsize[n]; ++r) {
1443:           DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew);
1444:           DMLabelSetValue(labelNew, pNew, values[val]);
1445:         }
1446:       }
1447:     }
1448:     ISRestoreIndices(pointIS, &points);
1449:     ISDestroy(&pointIS);
1450:   }
1451:   ISRestoreIndices(valueIS, &values);
1452:   ISDestroy(&valueIS);
1453:   return(0);
1454: }

1456: static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm)
1457: {
1458:   DM             dm;
1459:   PetscInt       numLabels, l;

1463:   DMPlexTransformGetDM(tr, &dm);
1464:   DMGetNumLabels(dm, &numLabels);
1465:   for (l = 0; l < numLabels; ++l) {
1466:     DMLabel         label, labelNew;
1467:     const char     *lname;
1468:     PetscBool       isDepth, isCellType;

1470:     DMGetLabelName(dm, l, &lname);
1471:     PetscStrcmp(lname, "depth", &isDepth);
1472:     if (isDepth) continue;
1473:     PetscStrcmp(lname, "celltype", &isCellType);
1474:     if (isCellType) continue;
1475:     DMCreateLabel(rdm, lname);
1476:     DMGetLabel(dm, lname, &label);
1477:     DMGetLabel(rdm, lname, &labelNew);
1478:     RefineLabel_Internal(tr, label, labelNew);
1479:   }
1480:   return(0);
1481: }

1483: /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */
1484: PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm)
1485: {
1486:   DM             dm;
1487:   PetscInt       Nf, f, Nds, s;

1491:   DMPlexTransformGetDM(tr, &dm);
1492:   DMGetNumFields(dm, &Nf);
1493:   for (f = 0; f < Nf; ++f) {
1494:     DMLabel     label, labelNew;
1495:     PetscObject obj;
1496:     const char *lname;

1498:     DMGetField(rdm, f, &label, &obj);
1499:     if (!label) continue;
1500:     PetscObjectGetName((PetscObject) label, &lname);
1501:     DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
1502:     RefineLabel_Internal(tr, label, labelNew);
1503:     DMSetField_Internal(rdm, f, labelNew, obj);
1504:     DMLabelDestroy(&labelNew);
1505:   }
1506:   DMGetNumDS(dm, &Nds);
1507:   for (s = 0; s < Nds; ++s) {
1508:     DMLabel     label, labelNew;
1509:     const char *lname;

1511:     DMGetRegionNumDS(rdm, s, &label, NULL, NULL);
1512:     if (!label) continue;
1513:     PetscObjectGetName((PetscObject) label, &lname);
1514:     DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);
1515:     RefineLabel_Internal(tr, label, labelNew);
1516:     DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL);
1517:     DMLabelDestroy(&labelNew);
1518:   }
1519:   return(0);
1520: }

1522: static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm)
1523: {
1524:   DM                 dm;
1525:   PetscSF            sf, sfNew;
1526:   PetscInt           numRoots, numLeaves, numLeavesNew = 0, l, m;
1527:   const PetscInt    *localPoints;
1528:   const PetscSFNode *remotePoints;
1529:   PetscInt          *localPointsNew;
1530:   PetscSFNode       *remotePointsNew;
1531:   PetscInt           pStartNew, pEndNew, pNew;
1532:   /* Brute force algorithm */
1533:   PetscSF            rsf;
1534:   PetscSection       s;
1535:   const PetscInt    *rootdegree;
1536:   PetscInt          *rootPointsNew, *remoteOffsets;
1537:   PetscInt           numPointsNew, pStart, pEnd, p;
1538:   PetscErrorCode     ierr;

1541:   DMPlexTransformGetDM(tr, &dm);
1542:   DMPlexGetChart(rdm, &pStartNew, &pEndNew);
1543:   DMGetPointSF(dm, &sf);
1544:   DMGetPointSF(rdm, &sfNew);
1545:   /* Calculate size of new SF */
1546:   PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);
1547:   if (numRoots < 0) return(0);
1548:   for (l = 0; l < numLeaves; ++l) {
1549:     const PetscInt  p = localPoints[l];
1550:     DMPolytopeType  ct;
1551:     DMPolytopeType *rct;
1552:     PetscInt       *rsize, *rcone, *rornt;
1553:     PetscInt        Nct, n;

1555:     DMPlexGetCellType(dm, p, &ct);
1556:     DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1557:     for (n = 0; n < Nct; ++n) {
1558:       numLeavesNew += rsize[n];
1559:     }
1560:   }
1561:   /* Send new root point numbers
1562:        It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication
1563:   */
1564:   DMPlexGetChart(dm, &pStart, &pEnd);
1565:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &s);
1566:   PetscSectionSetChart(s, pStart, pEnd);
1567:   for (p = pStart; p < pEnd; ++p) {
1568:     DMPolytopeType  ct;
1569:     DMPolytopeType *rct;
1570:     PetscInt       *rsize, *rcone, *rornt;
1571:     PetscInt        Nct, n;

1573:     DMPlexGetCellType(dm, p, &ct);
1574:     DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1575:     for (n = 0; n < Nct; ++n) {
1576:       PetscSectionAddDof(s, p, rsize[n]);
1577:     }
1578:   }
1579:   PetscSectionSetUp(s);
1580:   PetscSectionGetStorageSize(s, &numPointsNew);
1581:   PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets);
1582:   PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf);
1583:   PetscFree(remoteOffsets);
1584:   PetscSFComputeDegreeBegin(sf, &rootdegree);
1585:   PetscSFComputeDegreeEnd(sf, &rootdegree);
1586:   PetscMalloc1(numPointsNew, &rootPointsNew);
1587:   for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1;
1588:   for (p = pStart; p < pEnd; ++p) {
1589:     DMPolytopeType  ct;
1590:     DMPolytopeType *rct;
1591:     PetscInt       *rsize, *rcone, *rornt;
1592:     PetscInt        Nct, n, r, off;

1594:     if (!rootdegree[p-pStart]) continue;
1595:     PetscSectionGetOffset(s, p, &off);
1596:     DMPlexGetCellType(dm, p, &ct);
1597:     DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1598:     for (n = 0, m = 0; n < Nct; ++n) {
1599:       for (r = 0; r < rsize[n]; ++r, ++m) {
1600:         DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1601:         rootPointsNew[off+m] = pNew;
1602:       }
1603:     }
1604:   }
1605:   PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew,MPI_REPLACE);
1606:   PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew,MPI_REPLACE);
1607:   PetscSFDestroy(&rsf);
1608:   PetscMalloc1(numLeavesNew, &localPointsNew);
1609:   PetscMalloc1(numLeavesNew, &remotePointsNew);
1610:   for (l = 0, m = 0; l < numLeaves; ++l) {
1611:     const PetscInt  p = localPoints[l];
1612:     DMPolytopeType  ct;
1613:     DMPolytopeType *rct;
1614:     PetscInt       *rsize, *rcone, *rornt;
1615:     PetscInt        Nct, n, r, q, off;

1617:     PetscSectionGetOffset(s, p, &off);
1618:     DMPlexGetCellType(dm, p, &ct);
1619:     DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1620:     for (n = 0, q = 0; n < Nct; ++n) {
1621:       for (r = 0; r < rsize[n]; ++r, ++m, ++q) {
1622:         DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1623:         localPointsNew[m]        = pNew;
1624:         remotePointsNew[m].index = rootPointsNew[off+q];
1625:         remotePointsNew[m].rank  = remotePoints[l].rank;
1626:       }
1627:     }
1628:   }
1629:   PetscSectionDestroy(&s);
1630:   PetscFree(rootPointsNew);
1631:   /* SF needs sorted leaves to correctly calculate Gather */
1632:   {
1633:     PetscSFNode *rp, *rtmp;
1634:     PetscInt    *lp, *idx, *ltmp, i;

1636:     PetscMalloc1(numLeavesNew, &idx);
1637:     PetscMalloc1(numLeavesNew, &lp);
1638:     PetscMalloc1(numLeavesNew, &rp);
1639:     for (i = 0; i < numLeavesNew; ++i) {
1640:       if ((localPointsNew[i] < pStartNew) || (localPointsNew[i] >= pEndNew)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local SF point %D (%D) not in [%D, %D)", localPointsNew[i], i, pStartNew, pEndNew);
1641:       idx[i] = i;
1642:     }
1643:     PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx);
1644:     for (i = 0; i < numLeavesNew; ++i) {
1645:       lp[i] = localPointsNew[idx[i]];
1646:       rp[i] = remotePointsNew[idx[i]];
1647:     }
1648:     ltmp            = localPointsNew;
1649:     localPointsNew  = lp;
1650:     rtmp            = remotePointsNew;
1651:     remotePointsNew = rp;
1652:     PetscFree(idx);
1653:     PetscFree(ltmp);
1654:     PetscFree(rtmp);
1655:   }
1656:   PetscSFSetGraph(sfNew, pEndNew-pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);
1657:   return(0);
1658: }

1660: /*
1661:   DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of type ct with localized coordinates x, we generate localized coordinates xr for subcell r of type rct.

1663:   Not collective

1665:   Input Parameters:
1666: + tr  - The DMPlexTransform
1667: . ct  - The type of the parent cell
1668: . rct - The type of the produced cell
1669: . r   - The index of the produced cell
1670: - x   - The localized coordinates for the parent cell

1672:   Output Parameter:
1673: . xr  - The localized coordinates for the produced cell

1675:   Level: developer

1677: .seealso: DMPlexCellRefinerSetCoordinates()
1678: */
1679: static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[])
1680: {
1681:   PetscFE        fe = NULL;
1682:   PetscInt       cdim, v, *subcellV;

1686:   DMPlexTransformGetCoordinateFE(tr, ct, &fe);
1687:   DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV);
1688:   PetscFEGetNumComponents(fe, &cdim);
1689:   for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) {
1690:     PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v*cdim]);
1691:   }
1692:   return(0);
1693: }

1695: static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm)
1696: {
1697:   DM                    dm, cdm;
1698:   PetscSection          coordSection, coordSectionNew;
1699:   Vec                   coordsLocal, coordsLocalNew;
1700:   const PetscScalar    *coords;
1701:   PetscScalar          *coordsNew;
1702:   const DMBoundaryType *bd;
1703:   const PetscReal      *maxCell, *L;
1704:   PetscBool             isperiodic, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE;
1705:   PetscInt              dE, d, cStart, cEnd, c, vStartNew, vEndNew, v, pStart, pEnd, p, ocStart, ocEnd;
1706:   PetscErrorCode        ierr;

1709:   DMPlexTransformGetDM(tr, &dm);
1710:   DMGetCoordinateDM(dm, &cdm);
1711:   DMGetPeriodicity(dm, &isperiodic, &maxCell, &L, &bd);
1712:   /* Determine if we need to localize coordinates when generating them */
1713:   if (isperiodic) {
1714:     localizeVertices = PETSC_TRUE;
1715:     if (!maxCell) {
1716:       PetscBool localized;
1717:       DMGetCoordinatesLocalized(dm, &localized);
1718:       if (!localized) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Cannot refine a periodic mesh if coordinates have not been localized");
1719:       localizeCells = PETSC_TRUE;
1720:     }
1721:   }

1723:   DMGetCoordinateSection(dm, &coordSection);
1724:   PetscSectionGetFieldComponents(coordSection, 0, &dE);
1725:   if (maxCell) {
1726:     PetscReal maxCellNew[3];

1728:     for (d = 0; d < dE; ++d) maxCellNew[d] = maxCell[d]/2.0;
1729:     DMSetPeriodicity(rdm, isperiodic, maxCellNew, L, bd);
1730:   } else {
1731:     DMSetPeriodicity(rdm, isperiodic, maxCell, L, bd);
1732:   }
1733:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &coordSectionNew);
1734:   PetscSectionSetNumFields(coordSectionNew, 1);
1735:   PetscSectionSetFieldComponents(coordSectionNew, 0, dE);
1736:   DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew);
1737:   if (localizeCells) {PetscSectionSetChart(coordSectionNew, 0,         vEndNew);}
1738:   else               {PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew);}

1740:   /* Localization should be inherited */
1741:   /*   Stefano calculates parent cells for each new cell for localization */
1742:   /*   Localized cells need coordinates of closure */
1743:   for (v = vStartNew; v < vEndNew; ++v) {
1744:     PetscSectionSetDof(coordSectionNew, v, dE);
1745:     PetscSectionSetFieldDof(coordSectionNew, v, 0, dE);
1746:   }
1747:   if (localizeCells) {
1748:     DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1749:     for (c = cStart; c < cEnd; ++c) {
1750:       PetscInt dof;

1752:       PetscSectionGetDof(coordSection, c, &dof);
1753:       if (dof) {
1754:         DMPolytopeType  ct;
1755:         DMPolytopeType *rct;
1756:         PetscInt       *rsize, *rcone, *rornt;
1757:         PetscInt        dim, cNew, Nct, n, r;

1759:         DMPlexGetCellType(dm, c, &ct);
1760:         dim  = DMPolytopeTypeGetDim(ct);
1761:         DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1762:         /* This allows for different cell types */
1763:         for (n = 0; n < Nct; ++n) {
1764:           if (dim != DMPolytopeTypeGetDim(rct[n])) continue;
1765:           for (r = 0; r < rsize[n]; ++r) {
1766:             PetscInt *closure = NULL;
1767:             PetscInt  clSize, cl, Nv = 0;

1769:             DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew);
1770:             DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
1771:             for (cl = 0; cl < clSize*2; cl += 2) {if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;}
1772:             DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);
1773:             PetscSectionSetDof(coordSectionNew, cNew, Nv * dE);
1774:             PetscSectionSetFieldDof(coordSectionNew, cNew, 0, Nv * dE);
1775:           }
1776:         }
1777:       }
1778:     }
1779:   }
1780:   PetscSectionSetUp(coordSectionNew);
1781:   DMViewFromOptions(dm, NULL, "-coarse_dm_view");
1782:   DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew);
1783:   {
1784:     VecType     vtype;
1785:     PetscInt    coordSizeNew, bs;
1786:     const char *name;

1788:     DMGetCoordinatesLocal(dm, &coordsLocal);
1789:     VecCreate(PETSC_COMM_SELF, &coordsLocalNew);
1790:     PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);
1791:     VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE);
1792:     PetscObjectGetName((PetscObject) coordsLocal, &name);
1793:     PetscObjectSetName((PetscObject) coordsLocalNew, name);
1794:     VecGetBlockSize(coordsLocal, &bs);
1795:     VecSetBlockSize(coordsLocalNew, bs);
1796:     VecGetType(coordsLocal, &vtype);
1797:     VecSetType(coordsLocalNew, vtype);
1798:   }
1799:   VecGetArrayRead(coordsLocal, &coords);
1800:   VecGetArray(coordsLocalNew, &coordsNew);
1801:   PetscSectionGetChart(coordSection, &ocStart, &ocEnd);
1802:   DMPlexGetChart(dm, &pStart, &pEnd);
1803:   /* First set coordinates for vertices*/
1804:   for (p = pStart; p < pEnd; ++p) {
1805:     DMPolytopeType  ct;
1806:     DMPolytopeType *rct;
1807:     PetscInt       *rsize, *rcone, *rornt;
1808:     PetscInt        Nct, n, r;
1809:     PetscBool       hasVertex = PETSC_FALSE, isLocalized = PETSC_FALSE;

1811:     DMPlexGetCellType(dm, p, &ct);
1812:     DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1813:     for (n = 0; n < Nct; ++n) {
1814:       if (rct[n] == DM_POLYTOPE_POINT) {hasVertex = PETSC_TRUE; break;}
1815:     }
1816:     if (localizeVertices && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) {
1817:       PetscInt dof;
1818:       PetscSectionGetDof(coordSection, p, &dof);
1819:       if (dof) isLocalized = PETSC_TRUE;
1820:     }
1821:     if (hasVertex) {
1822:       const PetscScalar *icoords = NULL;
1823:       PetscScalar       *pcoords = NULL;
1824:       PetscInt          Nc, Nv, v, d;

1826:       DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);

1828:       icoords = pcoords;
1829:       Nv      = Nc/dE;
1830:       if (ct != DM_POLYTOPE_POINT) {
1831:         if (localizeVertices) {
1832:           PetscScalar anchor[3];

1834:           for (d = 0; d < dE; ++d) anchor[d] = pcoords[d];
1835:           if (!isLocalized) {
1836:             for (v = 0; v < Nv; ++v) {DMLocalizeCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], &pcoords[v*dE]);}
1837:           } else {
1838:             Nv = Nc/(2*dE);
1839:             for (v = Nv; v < Nv*2; ++v) {DMLocalizeCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], &pcoords[v*dE]);}
1840:           }
1841:         }
1842:       }
1843:       for (n = 0; n < Nct; ++n) {
1844:         if (rct[n] != DM_POLYTOPE_POINT) continue;
1845:         for (r = 0; r < rsize[n]; ++r) {
1846:           PetscScalar vcoords[3];
1847:           PetscInt    vNew, off;

1849:           DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew);
1850:           PetscSectionGetOffset(coordSectionNew, vNew, &off);
1851:           DMPlexTransformMapCoordinates(tr, ct, rct[n], r, Nv, dE, icoords, vcoords);
1852:           DMPlexSnapToGeomModel(dm, p, vcoords, &coordsNew[off]);
1853:         }
1854:       }
1855:       DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);
1856:     }
1857:   }
1858:   /* Then set coordinates for cells by localizing */
1859:   for (p = pStart; p < pEnd; ++p) {
1860:     DMPolytopeType  ct;
1861:     DMPolytopeType *rct;
1862:     PetscInt       *rsize, *rcone, *rornt;
1863:     PetscInt        Nct, n, r;
1864:     PetscBool       isLocalized = PETSC_FALSE;

1866:     DMPlexGetCellType(dm, p, &ct);
1867:     DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);
1868:     if (localizeCells && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) {
1869:       PetscInt dof;
1870:       PetscSectionGetDof(coordSection, p, &dof);
1871:       if (dof) isLocalized = PETSC_TRUE;
1872:     }
1873:     if (isLocalized) {
1874:       const PetscScalar *pcoords;

1876:       DMPlexPointLocalRead(cdm, p, coords, &pcoords);
1877:       for (n = 0; n < Nct; ++n) {
1878:         const PetscInt Nr = rsize[n];

1880:         if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue;
1881:         for (r = 0; r < Nr; ++r) {
1882:           PetscInt pNew, offNew;

1884:           /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that
1885:              DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger
1886:              cell to the ones it produces. */
1887:           DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);
1888:           PetscSectionGetOffset(coordSectionNew, pNew, &offNew);
1889:           DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew]);
1890:         }
1891:       }
1892:     }
1893:   }
1894:   VecRestoreArrayRead(coordsLocal, &coords);
1895:   VecRestoreArray(coordsLocalNew, &coordsNew);
1896:   DMSetCoordinatesLocal(rdm, coordsLocalNew);
1897:   /* TODO Stefano has a final reduction if some hybrid coordinates cannot be found. (needcoords) Should not be needed. */
1898:   VecDestroy(&coordsLocalNew);
1899:   PetscSectionDestroy(&coordSectionNew);
1900:   if (!localizeCells) {DMLocalizeCoordinates(rdm);}
1901:   return(0);
1902: }

1904: PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm)
1905: {
1906:   DM                     rdm;
1907:   DMPlexInterpolatedFlag interp;
1908:   PetscInt               dim, embedDim;
1909:   PetscErrorCode         ierr;

1915:   DMPlexTransformSetDM(tr, dm);

1917:   DMCreate(PetscObjectComm((PetscObject)dm), &rdm);
1918:   DMSetType(rdm, DMPLEX);
1919:   DMGetDimension(dm, &dim);
1920:   DMSetDimension(rdm, dim);
1921:   DMGetCoordinateDim(dm, &embedDim);
1922:   DMSetCoordinateDim(rdm, embedDim);
1923:   /* Calculate number of new points of each depth */
1924:   DMPlexIsInterpolated(dm, &interp);
1925:   if (interp != DMPLEX_INTERPOLATED_FULL) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement");
1926:   /* Step 1: Set chart */
1927:   DMPlexSetChart(rdm, 0, tr->ctStartNew[tr->ctOrder[DM_NUM_POLYTOPES]]);
1928:   /* Step 2: Set cone/support sizes (automatically stratifies) */
1929:   DMPlexTransformSetConeSizes(tr, rdm);
1930:   /* Step 3: Setup refined DM */
1931:   DMSetUp(rdm);
1932:   /* Step 4: Set cones and supports (automatically symmetrizes) */
1933:   DMPlexTransformSetCones(tr, rdm);
1934:   /* Step 5: Create pointSF */
1935:   DMPlexTransformCreateSF(tr, rdm);
1936:   /* Step 6: Create labels */
1937:   DMPlexTransformCreateLabels(tr, rdm);
1938:   /* Step 7: Set coordinates */
1939:   DMPlexTransformSetCoordinates(tr, rdm);
1940:   *tdm = rdm;
1941:   return(0);
1942: }

1944: PetscErrorCode DMPlexTransformAdaptLabel(DM dm, DMLabel adaptLabel, DM *rdm)
1945: {
1946:   DMPlexTransform tr;
1947:   DM              cdm, rcdm;
1948:   const char     *prefix;
1949:   PetscErrorCode  ierr;

1952:   DMPlexTransformCreate(PetscObjectComm((PetscObject) dm), &tr);
1953:   PetscObjectGetOptionsPrefix((PetscObject) dm, &prefix);
1954:   PetscObjectSetOptionsPrefix((PetscObject) tr,  prefix);
1955:   DMPlexTransformSetDM(tr, dm);
1956:   DMPlexTransformSetFromOptions(tr);
1957:   DMPlexTransformSetActive(tr, adaptLabel);
1958:   DMPlexTransformSetUp(tr);
1959:   PetscObjectViewFromOptions((PetscObject) tr, NULL, "-dm_plex_transform_view");
1960:   DMPlexTransformApply(tr, dm, rdm);
1961:   DMCopyDisc(dm, *rdm);
1962:   DMGetCoordinateDM(dm, &cdm);
1963:   DMGetCoordinateDM(*rdm, &rcdm);
1964:   DMCopyDisc(cdm, rcdm);
1965:   DMPlexTransformCreateDiscLabels(tr, *rdm);
1966:   DMCopyDisc(dm, *rdm);
1967:   DMPlexTransformDestroy(&tr);
1968:   return(0);
1969: }