#include /*I "petscdmplextransform.h" I*/ #include /* For PetscFEInterpolate_Static() */ PetscClassId DMPLEXTRANSFORM_CLASSID; PetscFunctionList DMPlexTransformList = NULL; PetscBool DMPlexTransformRegisterAllCalled = PETSC_FALSE; PetscLogEvent DMPLEXTRANSFORM_SetUp, DMPLEXTRANSFORM_Apply, DMPLEXTRANSFORM_SetConeSizes, DMPLEXTRANSFORM_SetCones, DMPLEXTRANSFORM_CreateSF, DMPLEXTRANSFORM_CreateLabels, DMPLEXTRANSFORM_SetCoordinates; /* Construct cell type order since we must loop over cell types in the same dimensional order they are stored in the plex if dm != NULL OR in standard plex ordering if dm == NULL */ static PetscErrorCode DMPlexCreateCellTypeOrder_Internal(DM dm, PetscInt dim, PetscInt *ctOrder[], PetscInt *ctOrderInv[]) { PetscInt *ctO, *ctOInv; PetscInt d, c, off = 0; PetscInt dimOrder[5] = {3, 2, 1, 0, -1}; PetscFunctionBegin; PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctO, DM_NUM_POLYTOPES + 1, &ctOInv)); if (dm) { // Order the dimensions by their starting location PetscInt hStart[4] = {-1, -1, -1, -1}; for (d = 0; d <= dim; ++d) PetscCall(DMPlexGetDepthStratum(dm, dim - d, &hStart[d], NULL)); PetscCall(PetscSortIntWithArray(dim + 1, hStart, &dimOrder[3 - dim])); } else if (dim > 1) { // Standard plex ordering. dimOrder is in correct order if dim > 1 off = 4 - dim; dimOrder[off++] = 0; for (d = dim - 1; d > 0; --d) dimOrder[off++] = d; } off = 0; for (d = 0; d < 5; ++d) { for (c = 0; c < DM_NUM_POLYTOPES; ++c) { if (c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) continue; if (DMPolytopeTypeGetDim((DMPolytopeType)c) == dimOrder[d]) ctO[off++] = c; } } for (c = 0; c < DM_NUM_POLYTOPES; ++c) { if (c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) ctO[off++] = c; } ctO[off++] = DM_NUM_POLYTOPES; PetscCheck(off == DM_NUM_POLYTOPES + 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid offset %" PetscInt_FMT " for cell type order", off); for (c = 0; c <= DM_NUM_POLYTOPES; ++c) ctOInv[ctO[c]] = c; *ctOrder = ctO; *ctOrderInv = ctOInv; PetscFunctionReturn(PETSC_SUCCESS); } /*@C DMPlexTransformRegister - Adds a new transform component implementation Not Collective Input Parameters: + name - The name of a new user-defined creation routine - create_func - The creation routine Example Usage: .vb DMPlexTransformRegister("my_transform", MyTransformCreate); .ve Then, your transform type can be chosen with the procedural interface via .vb DMPlexTransformCreate(MPI_Comm, DMPlexTransform *); DMPlexTransformSetType(DMPlexTransform, "my_transform"); .ve or at runtime via the option .vb -dm_plex_transform_type my_transform .ve Level: advanced Note: `DMPlexTransformRegister()` may be called multiple times to add several user-defined transforms .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformRegisterAll()`, `DMPlexTransformRegisterDestroy()` @*/ PetscErrorCode DMPlexTransformRegister(const char name[], PetscErrorCode (*create_func)(DMPlexTransform)) { PetscFunctionBegin; PetscCall(DMInitializePackage()); PetscCall(PetscFunctionListAdd(&DMPlexTransformList, name, create_func)); PetscFunctionReturn(PETSC_SUCCESS); } PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Filter(DMPlexTransform); PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform); PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToBox(DMPlexTransform); PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform); PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform); PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform); PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_1D(DMPlexTransform); PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform); PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Cohesive(DMPlexTransform); /*@C DMPlexTransformRegisterAll - Registers all of the transform components in the `DM` package. Not Collective Level: advanced .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `DMRegisterAll()`, `DMPlexTransformRegisterDestroy()` @*/ PetscErrorCode DMPlexTransformRegisterAll(void) { PetscFunctionBegin; if (DMPlexTransformRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS); DMPlexTransformRegisterAllCalled = PETSC_TRUE; PetscCall(DMPlexTransformRegister(DMPLEXTRANSFORMFILTER, DMPlexTransformCreate_Filter)); PetscCall(DMPlexTransformRegister(DMPLEXREFINEREGULAR, DMPlexTransformCreate_Regular)); PetscCall(DMPlexTransformRegister(DMPLEXREFINETOBOX, DMPlexTransformCreate_ToBox)); PetscCall(DMPlexTransformRegister(DMPLEXREFINEALFELD, DMPlexTransformCreate_Alfeld)); PetscCall(DMPlexTransformRegister(DMPLEXREFINEBOUNDARYLAYER, DMPlexTransformCreate_BL)); PetscCall(DMPlexTransformRegister(DMPLEXREFINESBR, DMPlexTransformCreate_SBR)); PetscCall(DMPlexTransformRegister(DMPLEXREFINE1D, DMPlexTransformCreate_1D)); PetscCall(DMPlexTransformRegister(DMPLEXEXTRUDETYPE, DMPlexTransformCreate_Extrude)); PetscCall(DMPlexTransformRegister(DMPLEXCOHESIVEEXTRUDE, DMPlexTransformCreate_Cohesive)); PetscFunctionReturn(PETSC_SUCCESS); } /*@C DMPlexTransformRegisterDestroy - This function destroys the registered `DMPlexTransformType`. It is called from `PetscFinalize()`. Not collective Level: developer .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRegisterAll()`, `DMPlexTransformType`, `PetscInitialize()` @*/ PetscErrorCode DMPlexTransformRegisterDestroy(void) { PetscFunctionBegin; PetscCall(PetscFunctionListDestroy(&DMPlexTransformList)); DMPlexTransformRegisterAllCalled = PETSC_FALSE; PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformCreate - Creates an empty transform object. The type can then be set with `DMPlexTransformSetType()`. Collective Input Parameter: . comm - The communicator for the transform object Output Parameter: . tr - The transform object Level: beginner .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPLEXREFINEREGULAR`, `DMPLEXTRANSFORMFILTER` @*/ PetscErrorCode DMPlexTransformCreate(MPI_Comm comm, DMPlexTransform *tr) { DMPlexTransform t; PetscFunctionBegin; PetscAssertPointer(tr, 2); *tr = NULL; PetscCall(DMInitializePackage()); PetscCall(PetscHeaderCreate(t, DMPLEXTRANSFORM_CLASSID, "DMPlexTransform", "Mesh Transform", "DMPlexTransform", comm, DMPlexTransformDestroy, DMPlexTransformView)); t->setupcalled = PETSC_FALSE; PetscCall(PetscCalloc2(DM_NUM_POLYTOPES, &t->coordFE, DM_NUM_POLYTOPES, &t->refGeom)); *tr = t; PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformSetType - Sets the particular implementation for a transform. Collective Input Parameters: + tr - The transform - method - The name of the transform type Options Database Key: . -dm_plex_transform_type - Sets the transform type; see `DMPlexTransformType` Level: intermediate .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformGetType()`, `DMPlexTransformCreate()` @*/ PetscErrorCode DMPlexTransformSetType(DMPlexTransform tr, DMPlexTransformType method) { PetscErrorCode (*r)(DMPlexTransform); PetscBool match; PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); PetscCall(PetscObjectTypeCompare((PetscObject)tr, method, &match)); if (match) PetscFunctionReturn(PETSC_SUCCESS); PetscCall(DMPlexTransformRegisterAll()); PetscCall(PetscFunctionListFind(DMPlexTransformList, method, &r)); PetscCheck(r, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMPlexTransform type: %s", method); PetscTryTypeMethod(tr, destroy); PetscCall(PetscMemzero(tr->ops, sizeof(*tr->ops))); PetscCall(PetscObjectChangeTypeName((PetscObject)tr, method)); PetscCall((*r)(tr)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformGetType - Gets the type name (as a string) from the transform. Not Collective Input Parameter: . tr - The `DMPlexTransform` Output Parameter: . type - The `DMPlexTransformType` name Level: intermediate .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPlexTransformCreate()` @*/ PetscErrorCode DMPlexTransformGetType(DMPlexTransform tr, DMPlexTransformType *type) { PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); PetscAssertPointer(type, 2); PetscCall(DMPlexTransformRegisterAll()); *type = ((PetscObject)tr)->type_name; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMPlexTransformView_Ascii(DMPlexTransform tr, PetscViewer v) { PetscViewerFormat format; PetscFunctionBegin; PetscCall(PetscViewerGetFormat(v, &format)); if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { const PetscInt *trTypes = NULL; IS trIS; PetscInt cols = 8; PetscInt Nrt = 8, f, g; if (tr->trType) PetscCall(DMLabelView(tr->trType, v)); PetscCall(PetscViewerASCIIPrintf(v, "Source Starts\n")); for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g])); PetscCall(PetscViewerASCIIPrintf(v, "\n")); for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStart[f])); PetscCall(PetscViewerASCIIPrintf(v, "\n")); PetscCall(PetscViewerASCIIPrintf(v, "Target Starts\n")); for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g])); PetscCall(PetscViewerASCIIPrintf(v, "\n")); for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStartNew[f])); PetscCall(PetscViewerASCIIPrintf(v, "\n")); if (tr->trType) { PetscCall(DMLabelGetNumValues(tr->trType, &Nrt)); PetscCall(DMLabelGetValueIS(tr->trType, &trIS)); PetscCall(ISGetIndices(trIS, &trTypes)); } PetscCall(PetscViewerASCIIPrintf(v, "Offsets\n")); PetscCall(PetscViewerASCIIPrintf(v, " ")); for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g])); PetscCall(PetscViewerASCIIPrintf(v, "\n")); for (f = 0; f < Nrt; ++f) { PetscCall(PetscViewerASCIIPrintf(v, "%2" PetscInt_FMT " |", trTypes ? trTypes[f] : f)); for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->offset[f * DM_NUM_POLYTOPES + g])); PetscCall(PetscViewerASCIIPrintf(v, " |\n")); } if (tr->trType) { PetscCall(ISRestoreIndices(trIS, &trTypes)); PetscCall(ISDestroy(&trIS)); } } PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformView - Views a `DMPlexTransform` Collective Input Parameters: + tr - the `DMPlexTransform` object to view - v - the viewer Level: beginner .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `PetscViewer`, `DMPlexTransformDestroy()`, `DMPlexTransformCreate()` @*/ PetscErrorCode DMPlexTransformView(DMPlexTransform tr, PetscViewer v) { PetscBool isascii; PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tr), &v)); PetscValidHeaderSpecific(v, PETSC_VIEWER_CLASSID, 2); PetscCheckSameComm(tr, 1, v, 2); PetscCall(PetscViewerCheckWritable(v)); PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tr, v)); PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii)); if (isascii) PetscCall(DMPlexTransformView_Ascii(tr, v)); PetscTryTypeMethod(tr, view, v); PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformSetFromOptions - Sets parameters in a transform from values in the options database Collective Input Parameter: . tr - the `DMPlexTransform` object to set options for Options Database Keys: + -dm_plex_transform_type - Set the transform type, e.g. refine_regular . -dm_plex_transform_label_match_strata - Only label points of the same stratum as the producing point . -dm_plex_transform_label_replica_inc - Increment for the label value to be multiplied by the replica number, so that the new label value is oldValue + r * inc . -dm_plex_transform_active - Name for active mesh label - -dm_plex_transform_active_values - Values in the active label Level: intermediate .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()` @*/ PetscErrorCode DMPlexTransformSetFromOptions(DMPlexTransform tr) { char typeName[1024], active[PETSC_MAX_PATH_LEN]; const char *defName = DMPLEXREFINEREGULAR; PetscBool flg, match; PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); PetscObjectOptionsBegin((PetscObject)tr); PetscCall(PetscOptionsFList("-dm_plex_transform_type", "DMPlexTransform", "DMPlexTransformSetType", DMPlexTransformList, defName, typeName, 1024, &flg)); if (flg) PetscCall(DMPlexTransformSetType(tr, typeName)); else if (!((PetscObject)tr)->type_name) PetscCall(DMPlexTransformSetType(tr, defName)); PetscCall(PetscOptionsBool("-dm_plex_transform_label_match_strata", "Only label points of the same stratum as the producing point", "", tr->labelMatchStrata, &match, &flg)); if (flg) PetscCall(DMPlexTransformSetMatchStrata(tr, match)); PetscCall(PetscOptionsInt("-dm_plex_transform_label_replica_inc", "Increment for the label value to be multiplied by the replica number", "", tr->labelReplicaInc, &tr->labelReplicaInc, NULL)); PetscCall(PetscOptionsString("-dm_plex_transform_active", "Name for active mesh label", "DMPlexTransformSetActive", active, active, sizeof(active), &flg)); if (flg) { DM dm; DMLabel label; PetscInt values[16]; PetscInt n = 16; PetscCall(DMPlexTransformGetDM(tr, &dm)); PetscCall(DMGetLabel(dm, active, &label)); PetscCall(PetscOptionsIntArray("-dm_plex_transform_active_values", "The label values to be active", "DMPlexTransformSetActive", values, &n, &flg)); if (flg && n) { DMLabel newlabel; PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Active", &newlabel)); for (PetscInt i = 0; i < n; ++i) { IS is; PetscCall(DMLabelGetStratumIS(label, values[i], &is)); PetscCall(DMLabelInsertIS(newlabel, is, values[i])); PetscCall(ISDestroy(&is)); } PetscCall(DMPlexTransformSetActive(tr, newlabel)); PetscCall(DMLabelDestroy(&newlabel)); } else { PetscCall(DMPlexTransformSetActive(tr, label)); } } PetscTryTypeMethod(tr, setfromoptions, PetscOptionsObject); /* process any options handlers added with PetscObjectAddOptionsHandler() */ PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)tr, PetscOptionsObject)); PetscOptionsEnd(); PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformDestroy - Destroys a `DMPlexTransform` Collective Input Parameter: . tr - the transform object to destroy Level: beginner .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()` @*/ PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr) { PetscInt c; PetscFunctionBegin; if (!*tr) PetscFunctionReturn(PETSC_SUCCESS); PetscValidHeaderSpecific(*tr, DMPLEXTRANSFORM_CLASSID, 1); if (--((PetscObject)*tr)->refct > 0) { *tr = NULL; PetscFunctionReturn(PETSC_SUCCESS); } PetscTryTypeMethod(*tr, destroy); PetscCall(DMDestroy(&(*tr)->dm)); PetscCall(DMLabelDestroy(&(*tr)->active)); PetscCall(DMLabelDestroy(&(*tr)->trType)); PetscCall(PetscFree2((*tr)->ctOrderOld, (*tr)->ctOrderInvOld)); PetscCall(PetscFree2((*tr)->ctOrderNew, (*tr)->ctOrderInvNew)); PetscCall(PetscFree2((*tr)->ctStart, (*tr)->ctStartNew)); PetscCall(PetscFree((*tr)->offset)); PetscCall(PetscFree2((*tr)->depthStart, (*tr)->depthEnd)); for (c = 0; c < DM_NUM_POLYTOPES; ++c) { PetscCall(PetscFEDestroy(&(*tr)->coordFE[c])); PetscCall(PetscFEGeomDestroy(&(*tr)->refGeom[c])); } if ((*tr)->trVerts) { for (c = 0; c < DM_NUM_POLYTOPES; ++c) { DMPolytopeType *rct; PetscInt *rsize, *rcone, *rornt, Nct, n, r; if (DMPolytopeTypeGetDim((DMPolytopeType)c) > 0 && c != DM_POLYTOPE_UNKNOWN_CELL && c != DM_POLYTOPE_UNKNOWN_FACE) { PetscCall(DMPlexTransformCellTransform(*tr, (DMPolytopeType)c, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); for (n = 0; n < Nct; ++n) { if (rct[n] == DM_POLYTOPE_POINT) continue; for (r = 0; r < rsize[n]; ++r) PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]][r])); PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]])); } } PetscCall(PetscFree((*tr)->trSubVerts[c])); PetscCall(PetscFree((*tr)->trVerts[c])); } } PetscCall(PetscFree3((*tr)->trNv, (*tr)->trVerts, (*tr)->trSubVerts)); PetscCall(PetscFree2((*tr)->coordFE, (*tr)->refGeom)); /* We do not destroy (*dm)->data here so that we can reference count backend objects */ PetscCall(PetscHeaderDestroy(tr)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrderOld[], PetscInt ctStart[], PetscInt **offset) { DMLabel trType = tr->trType; PetscInt c, cN, *off; PetscFunctionBegin; if (trType) { DM dm; IS rtIS; const PetscInt *reftypes; PetscInt Nrt, r; PetscCall(DMPlexTransformGetDM(tr, &dm)); PetscCall(DMLabelGetNumValues(trType, &Nrt)); PetscCall(DMLabelGetValueIS(trType, &rtIS)); PetscCall(ISGetIndices(rtIS, &reftypes)); PetscCall(PetscCalloc1(Nrt * DM_NUM_POLYTOPES, &off)); for (r = 0; r < Nrt; ++r) { const PetscInt rt = reftypes[r]; IS rtIS; const PetscInt *points; DMPolytopeType ct; PetscInt np, p; PetscCall(DMLabelGetStratumIS(trType, rt, &rtIS)); PetscCall(ISGetLocalSize(rtIS, &np)); PetscCall(ISGetIndices(rtIS, &points)); if (!np) continue; p = points[0]; PetscCall(ISRestoreIndices(rtIS, &points)); PetscCall(ISDestroy(&rtIS)); PetscCall(DMPlexGetCellType(dm, p, &ct)); for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) { const DMPolytopeType ctNew = (DMPolytopeType)cN; DMPolytopeType *rct; PetscInt *rsize, *cone, *ornt; PetscInt Nct, n, s; if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) { off[r * DM_NUM_POLYTOPES + ctNew] = -1; break; } off[r * DM_NUM_POLYTOPES + ctNew] = 0; for (s = 0; s <= r; ++s) { const PetscInt st = reftypes[s]; DMPolytopeType sct; PetscInt q, qrt; PetscCall(DMLabelGetStratumIS(trType, st, &rtIS)); PetscCall(ISGetLocalSize(rtIS, &np)); PetscCall(ISGetIndices(rtIS, &points)); if (!np) continue; q = points[0]; PetscCall(ISRestoreIndices(rtIS, &points)); PetscCall(ISDestroy(&rtIS)); PetscCall(DMPlexGetCellType(dm, q, &sct)); PetscCall(DMPlexTransformCellTransform(tr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt)); PetscCheck(st == qrt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Refine type %" PetscInt_FMT " of point %" PetscInt_FMT " does not match predicted type %" PetscInt_FMT, qrt, q, st); if (st == rt) { for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) break; if (n == Nct) off[r * DM_NUM_POLYTOPES + ctNew] = -1; break; } for (n = 0; n < Nct; ++n) { if (rct[n] == ctNew) { PetscInt sn; PetscCall(DMLabelGetStratumSize(trType, st, &sn)); off[r * DM_NUM_POLYTOPES + ctNew] += sn * rsize[n]; } } } } } PetscCall(ISRestoreIndices(rtIS, &reftypes)); PetscCall(ISDestroy(&rtIS)); } else { PetscCall(PetscCalloc1(DM_NUM_POLYTOPES * DM_NUM_POLYTOPES, &off)); for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) { const DMPolytopeType ct = (DMPolytopeType)c; for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) { const DMPolytopeType ctNew = (DMPolytopeType)cN; DMPolytopeType *rct; PetscInt *rsize, *cone, *ornt; PetscInt Nct, n, i; if (DMPolytopeTypeGetDim(ct) < 0 || ct == DM_POLYTOPE_UNKNOWN_CELL || ct == DM_POLYTOPE_UNKNOWN_FACE || DMPolytopeTypeGetDim(ctNew) < 0 || ctNew == DM_POLYTOPE_UNKNOWN_CELL || ctNew == DM_POLYTOPE_UNKNOWN_FACE) { off[ct * DM_NUM_POLYTOPES + ctNew] = -1; continue; } off[ct * DM_NUM_POLYTOPES + ctNew] = 0; for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) { const DMPolytopeType ict = (DMPolytopeType)ctOrderOld[i]; const DMPolytopeType ictn = (DMPolytopeType)ctOrderOld[i + 1]; PetscCall(DMPlexTransformCellTransform(tr, ict, PETSC_DETERMINE, NULL, &Nct, &rct, &rsize, &cone, &ornt)); if (ict == ct) { for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) break; if (n == Nct) off[ct * DM_NUM_POLYTOPES + ctNew] = -1; break; } for (n = 0; n < Nct; ++n) if (rct[n] == ctNew) off[ct * DM_NUM_POLYTOPES + ctNew] += (ctStart[ictn] - ctStart[ict]) * rsize[n]; } } } } *offset = off; PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformSetUp - Create the tables that drive the transform Input Parameter: . tr - The `DMPlexTransform` object Level: intermediate .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` @*/ PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr) { DMPolytopeType ctCell; DM dm; PetscInt pStart, pEnd, p, c, celldim = 0; PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); if (tr->setupcalled) PetscFunctionReturn(PETSC_SUCCESS); PetscCall(DMPlexTransformGetDM(tr, &dm)); PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetUp, tr, dm, 0, 0)); PetscTryTypeMethod(tr, setup); PetscCall(DMSetSnapToGeomModel(dm, NULL)); PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); if (pEnd > pStart) { // Ignore cells hanging off of embedded surfaces PetscInt c = pStart; ctCell = DM_POLYTOPE_FV_GHOST; while (DMPolytopeTypeGetDim(ctCell) < 0) PetscCall(DMPlexGetCellType(dm, c++, &ctCell)); } else { PetscInt dim; PetscCall(DMGetDimension(dm, &dim)); switch (dim) { case 0: ctCell = DM_POLYTOPE_POINT; break; case 1: ctCell = DM_POLYTOPE_SEGMENT; break; case 2: ctCell = DM_POLYTOPE_TRIANGLE; break; case 3: ctCell = DM_POLYTOPE_TETRAHEDRON; break; default: ctCell = DM_POLYTOPE_UNKNOWN; } } PetscCall(DMPlexCreateCellTypeOrder_Internal(dm, DMPolytopeTypeGetDim(ctCell), &tr->ctOrderOld, &tr->ctOrderInvOld)); for (p = pStart; p < pEnd; ++p) { DMPolytopeType ct; DMPolytopeType *rct; PetscInt *rsize, *cone, *ornt; PetscInt Nct, n; PetscCall(DMPlexGetCellType(dm, p, &ct)); PetscCheck(ct != DM_POLYTOPE_UNKNOWN && ct != DM_POLYTOPE_UNKNOWN_CELL && ct != DM_POLYTOPE_UNKNOWN_FACE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %" PetscInt_FMT, p); PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt)); for (n = 0; n < Nct; ++n) celldim = PetscMax(celldim, DMPolytopeTypeGetDim(rct[n])); } PetscCall(DMPlexCreateCellTypeOrder_Internal(NULL, celldim, &tr->ctOrderNew, &tr->ctOrderInvNew)); /* Construct sizes and offsets for each cell type */ if (!tr->ctStart) { PetscInt *ctS, *ctSN, *ctC, *ctCN; PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctS, DM_NUM_POLYTOPES + 1, &ctSN)); PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctC, DM_NUM_POLYTOPES + 1, &ctCN)); for (p = pStart; p < pEnd; ++p) { DMPolytopeType ct; DMPolytopeType *rct; PetscInt *rsize, *cone, *ornt; PetscInt Nct, n; PetscCall(DMPlexGetCellType(dm, p, &ct)); PetscCheck(ct != DM_POLYTOPE_UNKNOWN && ct != DM_POLYTOPE_UNKNOWN_CELL && ct != DM_POLYTOPE_UNKNOWN_FACE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %" PetscInt_FMT, p); ++ctC[ct]; PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt)); for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n]; } for (c = 0; c < DM_NUM_POLYTOPES; ++c) { const PetscInt cto = tr->ctOrderOld[c]; const PetscInt cton = tr->ctOrderOld[c + 1]; const PetscInt ctn = tr->ctOrderNew[c]; const PetscInt ctnn = tr->ctOrderNew[c + 1]; ctS[cton] = ctS[cto] + ctC[cto]; ctSN[ctnn] = ctSN[ctn] + ctCN[ctn]; } PetscCall(PetscFree2(ctC, ctCN)); tr->ctStart = ctS; tr->ctStartNew = ctSN; } PetscCall(DMPlexTransformCreateOffset_Internal(tr, tr->ctOrderOld, tr->ctStart, &tr->offset)); // Compute depth information tr->depth = -1; for (c = 0; c < DM_NUM_POLYTOPES; ++c) if (tr->ctStartNew[tr->ctOrderNew[c + 1]] > tr->ctStartNew[tr->ctOrderNew[c]]) tr->depth = PetscMax(tr->depth, DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c])); PetscCall(PetscMalloc2(tr->depth + 1, &tr->depthStart, tr->depth + 1, &tr->depthEnd)); for (PetscInt d = 0; d <= tr->depth; ++d) { tr->depthStart[d] = PETSC_INT_MAX; tr->depthEnd[d] = -1; } for (c = 0; c < DM_NUM_POLYTOPES; ++c) { const PetscInt dep = DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]); if (tr->ctStartNew[tr->ctOrderNew[c + 1]] <= tr->ctStartNew[tr->ctOrderNew[c]]) continue; tr->depthStart[dep] = PetscMin(tr->depthStart[dep], tr->ctStartNew[tr->ctOrderNew[c]]); tr->depthEnd[dep] = PetscMax(tr->depthEnd[dep], tr->ctStartNew[tr->ctOrderNew[c + 1]]); } tr->setupcalled = PETSC_TRUE; PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetUp, tr, dm, 0, 0)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformGetDM - Get the base `DM` for the transform Input Parameter: . tr - The `DMPlexTransform` object Output Parameter: . dm - The original `DM` which will be transformed Level: intermediate .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` @*/ PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm) { PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); PetscAssertPointer(dm, 2); *dm = tr->dm; PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformSetDM - Set the base `DM` for the transform Input Parameters: + tr - The `DMPlexTransform` object - dm - The original `DM` which will be transformed Level: intermediate Note: The user does not typically call this, as it is called by `DMPlexTransformApply()`. .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` @*/ PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm) { PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); PetscValidHeaderSpecific(dm, DM_CLASSID, 2); PetscCall(PetscObjectReference((PetscObject)dm)); PetscCall(DMDestroy(&tr->dm)); tr->dm = dm; PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformGetActive - Get the `DMLabel` marking the active points for the transform Input Parameter: . tr - The `DMPlexTransform` object Output Parameter: . active - The `DMLabel` indicating which points will be transformed Level: intermediate .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` @*/ PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active) { PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); PetscAssertPointer(active, 2); *active = tr->active; PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformSetActive - Set the `DMLabel` marking the active points for the transform Input Parameters: + tr - The `DMPlexTransform` object - active - The `DMLabel` indicating which points will be transformed Level: intermediate Note: This only applies to transforms listed in [](plex_transform_table) that operate on a subset of the mesh. .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` @*/ PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active) { PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); if (active) PetscValidHeaderSpecific(active, DMLABEL_CLASSID, 2); PetscCall(PetscObjectReference((PetscObject)active)); PetscCall(DMLabelDestroy(&tr->active)); tr->active = active; PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformGetTransformTypes - Get the `DMLabel` marking the transform type of each point for the transform Input Parameter: . tr - The `DMPlexTransform` object Output Parameter: . trType - The `DMLabel` indicating the transform type for each point Level: intermediate .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexSetTransformType()`, `DMPlexTransformGetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` @*/ PetscErrorCode DMPlexTransformGetTransformTypes(DMPlexTransform tr, DMLabel *trType) { PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); PetscAssertPointer(trType, 2); *trType = tr->trType; PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformSetTransformTypes - Set the `DMLabel` marking the transform type of each point for the transform Input Parameters: + tr - The `DMPlexTransform` object - trType - The original `DM` which will be transformed Level: intermediate .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetTransformTypes()`, `DMPlexTransformGetActive())`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` @*/ PetscErrorCode DMPlexTransformSetTransformTypes(DMPlexTransform tr, DMLabel trType) { PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); if (trType) PetscValidHeaderSpecific(trType, DMLABEL_CLASSID, 2); PetscCall(PetscObjectReference((PetscObject)trType)); PetscCall(DMLabelDestroy(&tr->trType)); tr->trType = trType; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe) { PetscFunctionBegin; if (!tr->coordFE[ct]) { PetscInt dim, cdim; dim = DMPolytopeTypeGetDim(ct); PetscCall(DMGetCoordinateDim(tr->dm, &cdim)); PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, cdim, ct, 1, PETSC_DETERMINE, &tr->coordFE[ct])); { PetscDualSpace dsp; PetscQuadrature quad; DM K; PetscFEGeom *cg; PetscScalar *Xq; PetscReal *xq, *wq; PetscInt Nq, q; PetscCall(DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq)); PetscCall(PetscMalloc1(Nq * cdim, &xq)); for (q = 0; q < Nq * cdim; ++q) xq[q] = PetscRealPart(Xq[q]); PetscCall(PetscMalloc1(Nq, &wq)); for (q = 0; q < Nq; ++q) wq[q] = 1.0; PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad)); PetscCall(PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq)); PetscCall(PetscFESetQuadrature(tr->coordFE[ct], quad)); PetscCall(PetscFEGetDualSpace(tr->coordFE[ct], &dsp)); PetscCall(PetscDualSpaceGetDM(dsp, &K)); PetscCall(PetscFEGeomCreate(quad, 1, cdim, PETSC_FEGEOM_BASIC, &tr->refGeom[ct])); cg = tr->refGeom[ct]; PetscCall(DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ)); PetscCall(PetscQuadratureDestroy(&quad)); } } *fe = tr->coordFE[ct]; PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode DMPlexTransformSetDimensions_Internal(DMPlexTransform tr, DM dm, DM tdm) { PetscInt dim, cdim; PetscFunctionBegin; PetscCall(DMGetDimension(dm, &dim)); PetscCall(DMSetDimension(tdm, dim)); PetscCall(DMGetCoordinateDim(dm, &cdim)); PetscCall(DMSetCoordinateDim(tdm, cdim)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformSetDimensions - Set the dimensions for the transformed `DM` Input Parameters: + tr - The `DMPlexTransform` object - dm - The original `DM` Output Parameter: . trdm - The transformed `DM` Level: advanced .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` @*/ PetscErrorCode DMPlexTransformSetDimensions(DMPlexTransform tr, DM dm, DM trdm) { PetscFunctionBegin; PetscUseTypeMethod(tr, setdimensions, dm, trdm); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode DMPlexTransformGetChart(DMPlexTransform tr, PetscInt *pStart, PetscInt *pEnd) { PetscFunctionBegin; if (pStart) *pStart = 0; if (pEnd) *pEnd = tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]]; PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode DMPlexTransformGetCellType(DMPlexTransform tr, PetscInt cell, DMPolytopeType *celltype) { PetscInt ctNew; PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); PetscAssertPointer(celltype, 3); /* TODO Can do bisection since everything is sorted */ for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) { PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]]; if (cell >= ctSN && cell < ctEN) break; } PetscCheck(ctNew < DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " cannot be located in the transformed mesh", cell); *celltype = (DMPolytopeType)ctNew; PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode DMPlexTransformGetCellTypeStratum(DMPlexTransform tr, DMPolytopeType celltype, PetscInt *start, PetscInt *end) { PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); if (start) *start = tr->ctStartNew[celltype]; if (end) *end = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[celltype] + 1]]; PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode DMPlexTransformGetDepth(DMPlexTransform tr, PetscInt *depth) { PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); *depth = tr->depth; PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode DMPlexTransformGetDepthStratum(DMPlexTransform tr, PetscInt depth, PetscInt *start, PetscInt *end) { PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); if (start) *start = tr->depthStart[depth]; if (end) *end = tr->depthEnd[depth]; PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformGetMatchStrata - Get the flag which determines what points get added to the transformed labels Not Collective Input Parameter: . tr - The `DMPlexTransform` Output Parameter: . match - If `PETSC_TRUE`, only add produced points at the same stratum as the original point to new labels Level: intermediate .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetMatchStrata()`, `DMPlexGetPointDepth()` @*/ PetscErrorCode DMPlexTransformGetMatchStrata(DMPlexTransform tr, PetscBool *match) { PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); PetscAssertPointer(match, 2); *match = tr->labelMatchStrata; PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformSetMatchStrata - Set the flag which determines what points get added to the transformed labels Not Collective Input Parameters: + tr - The `DMPlexTransform` - match - If `PETSC_TRUE`, only add produced points at the same stratum as the original point to new labels Level: intermediate .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetMatchStrata()`, `DMPlexGetPointDepth()` @*/ PetscErrorCode DMPlexTransformSetMatchStrata(DMPlexTransform tr, PetscBool match) { PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); tr->labelMatchStrata = match; PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformGetTargetPoint - Get the number of a point in the transformed mesh based on information from the original mesh. Not Collective Input Parameters: + tr - The `DMPlexTransform` . ct - The type of the original point which produces the new point . ctNew - The type of the new point . p - The original point which produces the new point - r - The replica number of the new point, meaning it is the rth point of type `ctNew` produced from `p` Output Parameter: . pNew - The new point number Level: developer .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSourcePoint()`, `DMPlexTransformCellTransform()` @*/ PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew) { DMPolytopeType *rct; PetscInt *rsize, *cone, *ornt; PetscInt rt, Nct, n, off, rp; DMLabel trType = tr->trType; PetscInt ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ct] + 1]]; PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]]; PetscInt newp = ctSN, cind; PetscFunctionBeginHot; PetscCheck(p >= ctS && p < ctE, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %" PetscInt_FMT " is not a %s [%" PetscInt_FMT ", %" PetscInt_FMT ")", p, DMPolytopeTypes[ct], ctS, ctE); PetscCall(DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt)); if (trType) { PetscCall(DMLabelGetValueIndex(trType, rt, &cind)); PetscCall(DMLabelGetStratumPointIndex(trType, rt, p, &rp)); PetscCheck(rp >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s point %" PetscInt_FMT " does not have refine type %" PetscInt_FMT, DMPolytopeTypes[ct], p, rt); } else { cind = ct; rp = p - ctS; } off = tr->offset[cind * DM_NUM_POLYTOPES + ctNew]; PetscCheck(off >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s (%" PetscInt_FMT ") of point %" PetscInt_FMT " does not produce type %s for transform %s", DMPolytopeTypes[ct], rt, p, DMPolytopeTypes[ctNew], tr->hdr.type_name); newp += off; for (n = 0; n < Nct; ++n) { if (rct[n] == ctNew) { if (rsize[n] && r >= rsize[n]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ") for subcell type %s in cell type %s", r, rsize[n], DMPolytopeTypes[rct[n]], DMPolytopeTypes[ct]); newp += rp * rsize[n] + r; break; } } PetscCheck(!(newp < ctSN) && !(newp >= ctEN), PETSC_COMM_SELF, PETSC_ERR_PLIB, "New point %" PetscInt_FMT " is not a %s [%" PetscInt_FMT ", %" PetscInt_FMT ")", newp, DMPolytopeTypes[ctNew], ctSN, ctEN); *pNew = newp; PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformGetSourcePoint - Get the number of a point in the original mesh based on information from the transformed mesh. Not Collective Input Parameters: + tr - The `DMPlexTransform` - pNew - The new point number Output Parameters: + ct - The type of the original point which produces the new point . ctNew - The type of the new point . p - The original point which produces the new point - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p Level: developer .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetTargetPoint()`, `DMPlexTransformCellTransform()` @*/ PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r) { DMLabel trType = tr->trType; DMPolytopeType *rct, ctN; PetscInt *rsize, *cone, *ornt; PetscInt rt = -1, rtTmp, Nct, n, rp = 0, rO = 0, pO; PetscInt offset = -1, ctS, ctE, ctO = 0, ctTmp, rtS; PetscFunctionBegin; PetscCall(DMPlexTransformGetCellType(tr, pNew, &ctN)); if (trType) { DM dm; IS rtIS; const PetscInt *reftypes; PetscInt Nrt, r, rtStart; PetscCall(DMPlexTransformGetDM(tr, &dm)); PetscCall(DMLabelGetNumValues(trType, &Nrt)); PetscCall(DMLabelGetValueIS(trType, &rtIS)); PetscCall(ISGetIndices(rtIS, &reftypes)); for (r = 0; r < Nrt; ++r) { const PetscInt off = tr->offset[r * DM_NUM_POLYTOPES + ctN]; if (tr->ctStartNew[ctN] + off > pNew) continue; /* Check that any of this refinement type exist */ /* TODO Actually keep track of the number produced here instead */ if (off > offset) { rt = reftypes[r]; offset = off; } } PetscCall(ISRestoreIndices(rtIS, &reftypes)); PetscCall(ISDestroy(&rtIS)); PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew); /* TODO Map refinement types to cell types */ PetscCall(DMLabelGetStratumBounds(trType, rt, &rtStart, NULL)); PetscCheck(rtStart >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Refinement type %" PetscInt_FMT " has no source points", rt); for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) { PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]]; if ((rtStart >= ctS) && (rtStart < ctE)) break; } PetscCheck(ctO != DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine a cell type for refinement type %" PetscInt_FMT, rt); } else { for (ctTmp = 0; ctTmp < DM_NUM_POLYTOPES; ++ctTmp) { const PetscInt off = tr->offset[ctTmp * DM_NUM_POLYTOPES + ctN]; if (tr->ctStartNew[ctN] + off > pNew) continue; if (tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctTmp] + 1]] <= tr->ctStart[ctTmp]) continue; /* TODO Actually keep track of the number produced here instead */ if (off > offset) { ctO = ctTmp; offset = off; } } rt = -1; PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew); } ctS = tr->ctStart[ctO]; ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]]; if (trType) { for (rtS = ctS; rtS < ctE; ++rtS) { PetscInt val; PetscCall(DMLabelGetValue(trType, rtS, &val)); if (val == rt) break; } PetscCheck(rtS < ctE, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find point of type %s with refine type %" PetscInt_FMT, DMPolytopeTypes[ctO], rt); } else rtS = ctS; PetscCall(DMPlexTransformCellTransform(tr, (DMPolytopeType)ctO, rtS, &rtTmp, &Nct, &rct, &rsize, &cone, &ornt)); PetscCheck(!trType || rt == rtTmp, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " has refine type %" PetscInt_FMT " != %" PetscInt_FMT " refine type which produced point %" PetscInt_FMT, rtS, rtTmp, rt, pNew); for (n = 0; n < Nct; ++n) { if (rct[n] == ctN) { PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset, val, c; if (trType) { for (c = ctS; c < ctE; ++c) { PetscCall(DMLabelGetValue(trType, c, &val)); if (val == rt) { if (tmp < rsize[n]) break; tmp -= rsize[n]; } } PetscCheck(c < ctE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Parent point for target point %" PetscInt_FMT " could be not found", pNew); rp = c - ctS; rO = tmp; } else { // This assumes that all points of type ctO transform the same way rp = tmp / rsize[n]; rO = tmp % rsize[n]; } break; } } PetscCheck(n != Nct, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number for target point %" PetscInt_FMT " could be not found", pNew); pO = rp + ctS; PetscCheck(!(pO < ctS) && !(pO >= ctE), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Source point %" PetscInt_FMT " is not a %s [%" PetscInt_FMT ", %" PetscInt_FMT ")", pO, DMPolytopeTypes[ctO], ctS, ctE); if (ct) *ct = (DMPolytopeType)ctO; if (ctNew) *ctNew = ctN; if (p) *p = pO; if (r) *r = rO; PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformCellTransform - Describes the transform of a given source cell into a set of other target cells. These produced cells become the new mesh. Input Parameters: + tr - The `DMPlexTransform` object . source - The source cell type - p - The source point, which can also determine the refine type Output Parameters: + rt - The refine type for this point . Nt - The number of types produced by this point . target - An array of length `Nt` giving the types produced . size - An array of length `Nt` giving the number of cells of each type produced . cone - An array of length `Nt`*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point - ornt - An array of length `Nt`*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point Level: advanced Notes: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical division (described in these outputs) of the cell in the original mesh. The point identifier is given by .vb the number of cones to be taken, or 0 for the current cell the cell cone point number at each level from which it is subdivided the replica number r of the subdivision. .ve The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is .vb Nt = 2 target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT} size = {1, 2} cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0} ornt = { 0, 0, 0, 0} .ve .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` @*/ PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) { PetscFunctionBegin; PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) { PetscFunctionBegin; *rnew = r; *onew = DMPolytopeTypeComposeOrientation(tct, o, so); PetscFunctionReturn(PETSC_SUCCESS); } /* Returns the same thing */ PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) { static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT}; static PetscInt vertexS[] = {1}; static PetscInt vertexC[] = {0}; static PetscInt vertexO[] = {0}; static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT}; static PetscInt edgeS[] = {1}; static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; static PetscInt edgeO[] = {0, 0}; static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR}; static PetscInt tedgeS[] = {1}; static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; static PetscInt tedgeO[] = {0, 0}; static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE}; static PetscInt triS[] = {1}; static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0}; static PetscInt triO[] = {0, 0, 0}; static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL}; static PetscInt quadS[] = {1}; 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}; static PetscInt quadO[] = {0, 0, 0, 0}; static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR}; static PetscInt tquadS[] = {1}; 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}; static PetscInt tquadO[] = {0, 0, 0, 0}; static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON}; static PetscInt tetS[] = {1}; 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}; static PetscInt tetO[] = {0, 0, 0, 0}; static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON}; static PetscInt hexS[] = {1}; static PetscInt hexC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 1, 5, 0}; static PetscInt hexO[] = {0, 0, 0, 0, 0, 0}; static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM}; static PetscInt tripS[] = {1}; static PetscInt tripC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0}; static PetscInt tripO[] = {0, 0, 0, 0, 0}; static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR}; static PetscInt ttripS[] = {1}; static PetscInt ttripC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0}; static PetscInt ttripO[] = {0, 0, 0, 0, 0}; static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR}; static PetscInt tquadpS[] = {1}; static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0}; static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0}; static DMPolytopeType pyrT[] = {DM_POLYTOPE_PYRAMID}; static PetscInt pyrS[] = {1}; static PetscInt pyrC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 1, 4, 0}; static PetscInt pyrO[] = {0, 0, 0, 0, 0}; PetscFunctionBegin; if (rt) *rt = 0; switch (source) { case DM_POLYTOPE_POINT: *Nt = 1; *target = vertexT; *size = vertexS; *cone = vertexC; *ornt = vertexO; break; case DM_POLYTOPE_SEGMENT: *Nt = 1; *target = edgeT; *size = edgeS; *cone = edgeC; *ornt = edgeO; break; case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break; case DM_POLYTOPE_TRIANGLE: *Nt = 1; *target = triT; *size = triS; *cone = triC; *ornt = triO; break; case DM_POLYTOPE_QUADRILATERAL: *Nt = 1; *target = quadT; *size = quadS; *cone = quadC; *ornt = quadO; break; case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 1; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break; case DM_POLYTOPE_TETRAHEDRON: *Nt = 1; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break; case DM_POLYTOPE_HEXAHEDRON: *Nt = 1; *target = hexT; *size = hexS; *cone = hexC; *ornt = hexO; break; case DM_POLYTOPE_TRI_PRISM: *Nt = 1; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break; case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 1; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break; case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 1; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break; case DM_POLYTOPE_PYRAMID: *Nt = 1; *target = pyrT; *size = pyrS; *cone = pyrC; *ornt = pyrO; break; default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); } PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point Not Collective Input Parameters: + tr - The `DMPlexTransform` . sct - The source point cell type, from whom the new cell is being produced . sp - The source point . so - The orientation of the source point in its enclosing parent . tct - The target point cell type . r - The replica number requested for the produced cell type - o - The orientation of the replica Output Parameters: + rnew - The replica number, given the orientation of the parent - onew - The replica orientation, given the orientation of the parent Level: advanced .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformCellTransform()`, `DMPlexTransformApply()` @*/ PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) { PetscFunctionBeginHot; PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm) { DM dm; PetscInt pStart, pEnd, pNew; PetscFunctionBegin; PetscCall(DMPlexTransformGetDM(tr, &dm)); PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetConeSizes, tr, dm, 0, 0)); /* Must create the celltype label here so that we do not automatically try to compute the types */ PetscCall(DMCreateLabel(rdm, "celltype")); PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); for (PetscInt p = pStart; p < pEnd; ++p) { DMPolytopeType ct; DMPolytopeType *rct; PetscInt *rsize, *rcone, *rornt; PetscInt Nct, n, r; PetscCall(DMPlexGetCellType(dm, p, &ct)); PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); for (n = 0; n < Nct; ++n) { for (r = 0; r < rsize[n]; ++r) { PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); PetscCall(DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]))); PetscCall(DMPlexSetCellType(rdm, pNew, rct[n])); } } } /* Let the DM know we have set all the cell types */ { DMLabel ctLabel; DM_Plex *plex = (DM_Plex *)rdm->data; PetscCall(DMPlexGetCellTypeLabel(rdm, &ctLabel)); PetscCall(PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState)); } PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetConeSizes, tr, dm, 0, 0)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize) { DMPolytopeType ctNew; PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); PetscAssertPointer(coneSize, 3); PetscCall(DMPlexTransformGetCellType(tr, q, &ctNew)); *coneSize = DMPolytopeTypeGetConeSize(ctNew); PetscFunctionReturn(PETSC_SUCCESS); } /* The orientation o is for the interior of the cell p */ static PetscErrorCode DMPlexTransformGetCone_Internal(DMPlexTransform tr, PetscInt p, PetscInt o, DMPolytopeType ct, DMPolytopeType ctNew, const PetscInt rcone[], PetscInt *coneoff, const PetscInt rornt[], PetscInt *orntoff, PetscInt coneNew[], PetscInt orntNew[]) { DM dm; const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); const PetscInt *cone; DMPolytopeType *newft = NULL; PetscInt c, coff = *coneoff, ooff = *orntoff; PetscInt dim, cr = 0, co = 0, nr, no; PetscFunctionBegin; PetscCall(DMPlexTransformGetDM(tr, &dm)); PetscCall(DMPlexGetOrientedCone(dm, p, &cone, NULL)); // Check if we have to permute this cell PetscCall(DMGetDimension(dm, &dim)); if (DMPolytopeTypeGetDim(ctNew) == dim && DMPolytopeTypeGetDim(ct) == dim - 1) { PetscCall(DMPlexTransformGetSubcellOrientation(tr, ct, p, o, ctNew, cr, co, &nr, &no)); if (cr != nr || co != no) PetscCall(DMGetWorkArray(dm, csizeNew, MPIU_INT, &newft)); } for (c = 0; c < csizeNew; ++c) { PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */ PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */ PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */ DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */ const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */ PetscInt pr = -1; /* Replica number of pp that produces new cone point */ const DMPolytopeType ft = (DMPolytopeType)rcone[coff++]; /* Cell type for new cone point of pNew */ const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */ PetscInt lc; /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */ for (lc = 0; lc < fn; ++lc) { const PetscInt *parr = DMPolytopeTypeGetArrangement(pct, po); const PetscInt acp = rcone[coff++]; const PetscInt pcp = parr[acp * 2]; const PetscInt pco = parr[acp * 2 + 1]; const PetscInt *ppornt; ppp = pp; pp = pcone[pcp]; PetscCall(DMPlexGetCellType(dm, pp, &pct)); // Restore the parent cone from the last iterate if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, ppp, &pcone, NULL)); PetscCall(DMPlexGetOrientedCone(dm, pp, &pcone, NULL)); PetscCall(DMPlexGetOrientedCone(dm, ppp, NULL, &ppornt)); po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco); PetscCall(DMPlexRestoreOrientedCone(dm, ppp, NULL, &ppornt)); } if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, pp, &pcone, NULL)); pr = rcone[coff++]; /* Orientation po of pp maps (pr, fo) -> (pr', fo') */ PetscCall(DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo)); PetscCall(DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c])); orntNew[c] = fo; if (newft) newft[c] = ft; } PetscCall(DMPlexRestoreOrientedCone(dm, p, &cone, NULL)); if (newft) { const PetscInt *arr; PetscInt *newcone, *newornt; arr = DMPolytopeTypeGetArrangement(ctNew, no); PetscCall(DMGetWorkArray(dm, csizeNew, MPIU_INT, &newcone)); PetscCall(DMGetWorkArray(dm, csizeNew, MPIU_INT, &newornt)); for (PetscInt c = 0; c < csizeNew; ++c) { DMPolytopeType ft = newft[c]; PetscInt nO; nO = DMPolytopeTypeGetNumArrangements(ft) / 2; newcone[c] = coneNew[arr[c * 2 + 0]]; newornt[c] = DMPolytopeTypeComposeOrientation(ft, arr[c * 2 + 1], orntNew[arr[c * 2 + 0]]); PetscCheck(!newornt[c] || !(newornt[c] >= nO || newornt[c] < -nO), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid orientation %" PetscInt_FMT " not in [%" PetscInt_FMT ",%" PetscInt_FMT ") for %s %" PetscInt_FMT, newornt[c], -nO, nO, DMPolytopeTypes[ft], coneNew[c]); } for (PetscInt c = 0; c < csizeNew; ++c) { coneNew[c] = newcone[c]; orntNew[c] = newornt[c]; } PetscCall(DMRestoreWorkArray(dm, csizeNew, MPIU_INT, &newcone)); PetscCall(DMRestoreWorkArray(dm, csizeNew, MPIU_INT, &newornt)); PetscCall(DMRestoreWorkArray(dm, csizeNew, MPIU_INT, &newft)); } *coneoff = coff; *orntoff = ooff; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm) { DM dm; DMPolytopeType ct; PetscInt *coneNew, *orntNew; PetscInt maxConeSize = 0, pStart, pEnd, p, pNew; PetscFunctionBegin; PetscCall(DMPlexTransformGetDM(tr, &dm)); PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetCones, tr, dm, 0, 0)); for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p)); PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew)); PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew)); PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); for (p = pStart; p < pEnd; ++p) { PetscInt coff, ooff; DMPolytopeType *rct; PetscInt *rsize, *rcone, *rornt; PetscInt Nct, n, r; PetscCall(DMPlexGetCellType(dm, p, &ct)); PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) { const DMPolytopeType ctNew = rct[n]; for (r = 0; r < rsize[n]; ++r) { PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew)); PetscCall(DMPlexSetCone(rdm, pNew, coneNew)); PetscCall(DMPlexSetConeOrientation(rdm, pNew, orntNew)); } } } PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew)); PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew)); PetscCall(DMViewFromOptions(rdm, NULL, "-rdm_view")); PetscCall(DMPlexSymmetrize(rdm)); PetscCall(DMPlexStratify(rdm)); PetscTryTypeMethod(tr, ordersupports, dm, rdm); PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetCones, tr, dm, 0, 0)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[]) { DM dm; DMPolytopeType ct, qct; DMPolytopeType *rct; PetscInt *rsize, *rcone, *rornt, *qcone, *qornt; PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0; PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); PetscAssertPointer(cone, 4); PetscAssertPointer(ornt, 5); for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p)); PetscCall(DMPlexTransformGetDM(tr, &dm)); PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone)); PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt)); PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r)); PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); for (n = 0; n < Nct; ++n) { const DMPolytopeType ctNew = rct[n]; const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); PetscInt Nr = rsize[n], fn, c; if (ctNew == qct) Nr = r; for (nr = 0; nr < Nr; ++nr) { for (c = 0; c < csizeNew; ++c) { ++coff; /* Cell type of new cone point */ fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ coff += fn; ++coff; /* Replica number of new cone point */ ++ooff; /* Orientation of new cone point */ } } if (ctNew == qct) break; } PetscCall(DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt)); *cone = qcone; *ornt = qornt; PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[]) { DM dm; DMPolytopeType ct, qct; DMPolytopeType *rct; PetscInt *rsize, *rcone, *rornt, *qcone, *qornt; PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0; PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); if (cone) PetscAssertPointer(cone, 3); if (ornt) PetscAssertPointer(ornt, 4); for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p)); PetscCall(DMPlexTransformGetDM(tr, &dm)); PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone)); PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt)); PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r)); PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); for (n = 0; n < Nct; ++n) { const DMPolytopeType ctNew = rct[n]; const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); PetscInt Nr = rsize[n], fn, c; if (ctNew == qct) Nr = r; for (nr = 0; nr < Nr; ++nr) { for (c = 0; c < csizeNew; ++c) { ++coff; /* Cell type of new cone point */ fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ coff += fn; ++coff; /* Replica number of new cone point */ ++ooff; /* Orientation of new cone point */ } } if (ctNew == qct) break; } PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt)); if (cone) *cone = qcone; else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qcone)); if (ornt) *ornt = qornt; else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qornt)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[]) { DM dm; PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); PetscCall(DMPlexTransformGetDM(tr, &dm)); if (cone) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, cone)); if (ornt) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, ornt)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr) { PetscInt ict; PetscFunctionBegin; PetscCall(PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts)); for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) { const DMPolytopeType ct = (DMPolytopeType)ict; DMPlexTransform reftr; DM refdm, trdm; Vec coordinates; const PetscScalar *coords; DMPolytopeType *rct; PetscInt *rsize, *rcone, *rornt; PetscInt Nct, n, r, pNew = 0; PetscInt trdim, vStart, vEnd, Nc; const PetscInt debug = 0; const char *typeName; /* Since points are 0-dimensional, coordinates make no sense */ if (DMPolytopeTypeGetDim(ct) <= 0 || ct == DM_POLYTOPE_UNKNOWN_CELL || ct == DM_POLYTOPE_UNKNOWN_FACE) continue; PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm)); PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &reftr)); PetscCall(DMPlexTransformSetDM(reftr, refdm)); PetscCall(DMPlexTransformGetType(tr, &typeName)); PetscCall(DMPlexTransformSetType(reftr, typeName)); PetscCall(DMPlexTransformSetUp(reftr)); PetscCall(DMPlexTransformApply(reftr, refdm, &trdm)); PetscCall(DMGetDimension(trdm, &trdim)); PetscCall(DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd)); tr->trNv[ct] = vEnd - vStart; PetscCall(DMGetCoordinatesLocal(trdm, &coordinates)); PetscCall(VecGetLocalSize(coordinates, &Nc)); PetscCheck(tr->trNv[ct] * trdim == Nc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell type %s, transformed coordinate size %" PetscInt_FMT " != %" PetscInt_FMT " size of coordinate storage", DMPolytopeTypes[ct], tr->trNv[ct] * trdim, Nc); PetscCall(PetscCalloc1(Nc, &tr->trVerts[ct])); PetscCall(VecGetArrayRead(coordinates, &coords)); PetscCall(PetscArraycpy(tr->trVerts[ct], coords, Nc)); PetscCall(VecRestoreArrayRead(coordinates, &coords)); PetscCall(PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct])); PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); for (n = 0; n < Nct; ++n) { /* Since points are 0-dimensional, coordinates make no sense */ if (rct[n] == DM_POLYTOPE_POINT) continue; PetscCall(PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]])); for (r = 0; r < rsize[n]; ++r) { PetscInt *closure = NULL; PetscInt clSize, cl, Nv = 0; PetscCall(PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r])); PetscCall(DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew)); PetscCall(DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure)); for (cl = 0; cl < clSize * 2; cl += 2) { const PetscInt sv = closure[cl]; if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart; } PetscCall(DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure)); PetscCheck(Nv == DMPolytopeTypeGetNumVertices(rct[n]), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of vertices %" PetscInt_FMT " != %" PetscInt_FMT " for %s subcell %" PetscInt_FMT " from cell %s", Nv, DMPolytopeTypeGetNumVertices(rct[n]), DMPolytopeTypes[rct[n]], r, DMPolytopeTypes[ct]); } } if (debug) { DMPolytopeType *rct; PetscInt *rsize, *rcone, *rornt; PetscInt v, dE = trdim, d, off = 0; PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct])); for (v = 0; v < tr->trNv[ct]; ++v) { PetscCall(PetscPrintf(PETSC_COMM_SELF, " ")); for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++]))); PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); } PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); for (n = 0; n < Nct; ++n) { if (rct[n] == DM_POLYTOPE_POINT) continue; PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct])); for (r = 0; r < rsize[n]; ++r) { PetscCall(PetscPrintf(PETSC_COMM_SELF, " ")); for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v])); PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); } } } PetscCall(DMDestroy(&refdm)); PetscCall(DMDestroy(&trdm)); PetscCall(DMPlexTransformDestroy(&reftr)); } PetscFunctionReturn(PETSC_SUCCESS); } /*@C DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type Input Parameters: + tr - The `DMPlexTransform` object - ct - The cell type Output Parameters: + Nv - The number of transformed vertices in the closure of the reference cell of given type - trVerts - The coordinates of these vertices in the reference cell Level: developer .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSubcellVertices()` @*/ PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[]) { PetscFunctionBegin; if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr)); if (Nv) *Nv = tr->trNv[ct]; if (trVerts) *trVerts = tr->trVerts[ct]; PetscFunctionReturn(PETSC_SUCCESS); } /*@C DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type Input Parameters: + tr - The `DMPlexTransform` object . ct - The cell type . rct - The subcell type - r - The subcell index Output Parameter: . subVerts - The indices of these vertices in the set of vertices returned by `DMPlexTransformGetCellVertices()` Level: developer .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetCellVertices()` @*/ PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[]) { PetscFunctionBegin; if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr)); PetscCheck(tr->trSubVerts[ct][rct], PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]); if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r]; PetscFunctionReturn(PETSC_SUCCESS); } /* Computes new vertex as the barycenter, or centroid */ PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) { PetscInt v, d; PetscFunctionBeginHot; PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]); for (d = 0; d < dE; ++d) out[d] = 0.0; for (v = 0; v < Nv; ++v) for (d = 0; d < dE; ++d) out[d] += in[v * dE + d]; for (d = 0; d < dE; ++d) out[d] /= Nv; PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformMapCoordinates - Calculate new coordinates for produced points Not collective Input Parameters: + tr - The `DMPlexTransform` . pct - The cell type of the parent, from whom the new cell is being produced . ct - The type being produced . p - The original point . r - The replica number requested for the produced cell type . Nv - Number of vertices in the closure of the parent cell . dE - Spatial dimension - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell Output Parameter: . out - The coordinates of the new vertices Level: intermediate .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()` @*/ PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) { PetscFunctionBeginHot; if (Nv) PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out); PetscFunctionReturn(PETSC_SUCCESS); } /* DMPlexTransformLabelProducedPoint_Private - Label a produced point based on its parent label Not Collective Input Parameters: + tr - The `DMPlexTransform` . label - The label in the transformed mesh . pp - The parent point in the original mesh . pct - The cell type of the parent point . p - The point in the transformed mesh . ct - The cell type of the point . r - The replica number of the point - val - The label value of the parent point Level: developer .seealso: `DMPlexTransformCreateLabels()`, `RefineLabel_Internal()` */ static PetscErrorCode DMPlexTransformLabelProducedPoint_Private(DMPlexTransform tr, DMLabel label, PetscInt pp, DMPolytopeType pct, PetscInt p, DMPolytopeType ct, PetscInt r, PetscInt val) { PetscFunctionBeginHot; if (tr->labelMatchStrata && pct != ct) PetscFunctionReturn(PETSC_SUCCESS); PetscCall(DMLabelSetValue(label, p, val + tr->labelReplicaInc * r)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew) { DM dm; IS valueIS; const PetscInt *values; PetscInt defVal, Nv, val; PetscFunctionBegin; PetscCall(DMPlexTransformGetDM(tr, &dm)); PetscCall(DMLabelGetDefaultValue(label, &defVal)); PetscCall(DMLabelSetDefaultValue(labelNew, defVal)); PetscCall(DMLabelGetValueIS(label, &valueIS)); PetscCall(ISGetLocalSize(valueIS, &Nv)); PetscCall(ISGetIndices(valueIS, &values)); for (val = 0; val < Nv; ++val) { IS pointIS; const PetscInt *points; PetscInt numPoints, p; /* Ensure refined label is created with same number of strata as * original (even if no entries here). */ PetscCall(DMLabelAddStratum(labelNew, values[val])); PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS)); PetscCall(ISGetLocalSize(pointIS, &numPoints)); PetscCall(ISGetIndices(pointIS, &points)); for (p = 0; p < numPoints; ++p) { const PetscInt point = points[p]; DMPolytopeType ct; DMPolytopeType *rct; PetscInt *rsize, *rcone, *rornt; PetscInt Nct, n, r, pNew = 0; PetscCall(DMPlexGetCellType(dm, point, &ct)); PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); for (n = 0; n < Nct; ++n) { for (r = 0; r < rsize[n]; ++r) { PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew)); PetscCall(DMPlexTransformLabelProducedPoint_Private(tr, labelNew, point, ct, pNew, rct[n], r, values[val])); } } } PetscCall(ISRestoreIndices(pointIS, &points)); PetscCall(ISDestroy(&pointIS)); } PetscCall(ISRestoreIndices(valueIS, &values)); PetscCall(ISDestroy(&valueIS)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm) { DM dm; PetscInt numLabels, l; PetscFunctionBegin; PetscCall(DMPlexTransformGetDM(tr, &dm)); PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_CreateLabels, tr, dm, 0, 0)); PetscCall(DMGetNumLabels(dm, &numLabels)); for (l = 0; l < numLabels; ++l) { DMLabel label, labelNew; const char *lname; PetscBool isDepth, isCellType; PetscCall(DMGetLabelName(dm, l, &lname)); PetscCall(PetscStrcmp(lname, "depth", &isDepth)); if (isDepth) continue; PetscCall(PetscStrcmp(lname, "celltype", &isCellType)); if (isCellType) continue; PetscCall(DMCreateLabel(rdm, lname)); PetscCall(DMGetLabel(dm, lname, &label)); PetscCall(DMGetLabel(rdm, lname, &labelNew)); PetscCall(RefineLabel_Internal(tr, label, labelNew)); } PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_CreateLabels, tr, dm, 0, 0)); PetscFunctionReturn(PETSC_SUCCESS); } /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */ PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm) { DM dm; PetscInt Nf, f, Nds, s; PetscFunctionBegin; PetscCall(DMPlexTransformGetDM(tr, &dm)); PetscCall(DMGetNumFields(dm, &Nf)); for (f = 0; f < Nf; ++f) { DMLabel label, labelNew; PetscObject obj; const char *lname; PetscCall(DMGetField(rdm, f, &label, &obj)); if (!label) continue; PetscCall(PetscObjectGetName((PetscObject)label, &lname)); PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew)); PetscCall(RefineLabel_Internal(tr, label, labelNew)); PetscCall(DMSetField_Internal(rdm, f, labelNew, obj)); PetscCall(DMLabelDestroy(&labelNew)); } PetscCall(DMGetNumDS(dm, &Nds)); for (s = 0; s < Nds; ++s) { DMLabel label, labelNew; const char *lname; PetscCall(DMGetRegionNumDS(rdm, s, &label, NULL, NULL, NULL)); if (!label) continue; PetscCall(PetscObjectGetName((PetscObject)label, &lname)); PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew)); PetscCall(RefineLabel_Internal(tr, label, labelNew)); PetscCall(DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL, NULL)); PetscCall(DMLabelDestroy(&labelNew)); } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm) { DM dm; PetscSF sf, sfNew; PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m; const PetscInt *localPoints; const PetscSFNode *remotePoints; PetscInt *localPointsNew; PetscSFNode *remotePointsNew; PetscInt pStartNew, pEndNew, pNew; /* Brute force algorithm */ PetscSF rsf; PetscSection s; const PetscInt *rootdegree; PetscInt *rootPointsNew, *remoteOffsets; PetscInt numPointsNew, pStart, pEnd, p; PetscFunctionBegin; PetscCall(DMPlexTransformGetDM(tr, &dm)); PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_CreateSF, tr, dm, 0, 0)); PetscCall(DMPlexGetChart(rdm, &pStartNew, &pEndNew)); PetscCall(DMGetPointSF(dm, &sf)); PetscCall(DMGetPointSF(rdm, &sfNew)); /* Calculate size of new SF */ PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints)); if (numRoots < 0) { PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_CreateSF, tr, dm, 0, 0)); PetscFunctionReturn(PETSC_SUCCESS); } for (l = 0; l < numLeaves; ++l) { const PetscInt p = localPoints[l]; DMPolytopeType ct; DMPolytopeType *rct; PetscInt *rsize, *rcone, *rornt; PetscInt Nct, n; PetscCall(DMPlexGetCellType(dm, p, &ct)); PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n]; } /* Send new root point numbers It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication */ PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s)); PetscCall(PetscSectionSetChart(s, pStart, pEnd)); for (p = pStart; p < pEnd; ++p) { DMPolytopeType ct; DMPolytopeType *rct; PetscInt *rsize, *rcone, *rornt; PetscInt Nct, n; PetscCall(DMPlexGetCellType(dm, p, &ct)); PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); for (n = 0; n < Nct; ++n) PetscCall(PetscSectionAddDof(s, p, rsize[n])); } PetscCall(PetscSectionSetUp(s)); PetscCall(PetscSectionGetStorageSize(s, &numPointsNew)); PetscCall(PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets)); PetscCall(PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf)); PetscCall(PetscFree(remoteOffsets)); PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree)); PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree)); PetscCall(PetscMalloc1(numPointsNew, &rootPointsNew)); for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1; for (p = pStart; p < pEnd; ++p) { DMPolytopeType ct; DMPolytopeType *rct; PetscInt *rsize, *rcone, *rornt; PetscInt Nct, n, r, off; if (!rootdegree[p - pStart]) continue; PetscCall(PetscSectionGetOffset(s, p, &off)); PetscCall(DMPlexGetCellType(dm, p, &ct)); PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); for (n = 0, m = 0; n < Nct; ++n) { for (r = 0; r < rsize[n]; ++r, ++m) { PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); rootPointsNew[off + m] = pNew; } } } PetscCall(PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE)); PetscCall(PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE)); PetscCall(PetscSFDestroy(&rsf)); PetscCall(PetscMalloc1(numLeavesNew, &localPointsNew)); PetscCall(PetscMalloc1(numLeavesNew, &remotePointsNew)); for (l = 0, m = 0; l < numLeaves; ++l) { const PetscInt p = localPoints[l]; DMPolytopeType ct; DMPolytopeType *rct; PetscInt *rsize, *rcone, *rornt; PetscInt Nct, n, r, q, off; PetscCall(PetscSectionGetOffset(s, p, &off)); PetscCall(DMPlexGetCellType(dm, p, &ct)); PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); for (n = 0, q = 0; n < Nct; ++n) { for (r = 0; r < rsize[n]; ++r, ++m, ++q) { PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); localPointsNew[m] = pNew; remotePointsNew[m].index = rootPointsNew[off + q]; remotePointsNew[m].rank = remotePoints[l].rank; } } } PetscCall(PetscSectionDestroy(&s)); PetscCall(PetscFree(rootPointsNew)); /* SF needs sorted leaves to correctly calculate Gather */ { PetscSFNode *rp, *rtmp; PetscInt *lp, *idx, *ltmp, i; PetscCall(PetscMalloc1(numLeavesNew, &idx)); PetscCall(PetscMalloc1(numLeavesNew, &lp)); PetscCall(PetscMalloc1(numLeavesNew, &rp)); for (i = 0; i < numLeavesNew; ++i) { PetscCheck(!(localPointsNew[i] < pStartNew) && !(localPointsNew[i] >= pEndNew), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local SF point %" PetscInt_FMT " (%" PetscInt_FMT ") not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", localPointsNew[i], i, pStartNew, pEndNew); idx[i] = i; } PetscCall(PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx)); for (i = 0; i < numLeavesNew; ++i) { lp[i] = localPointsNew[idx[i]]; rp[i] = remotePointsNew[idx[i]]; } ltmp = localPointsNew; localPointsNew = lp; rtmp = remotePointsNew; remotePointsNew = rp; PetscCall(PetscFree(idx)); PetscCall(PetscFree(ltmp)); PetscCall(PetscFree(rtmp)); } PetscCall(PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER)); PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_CreateSF, tr, dm, 0, 0)); PetscFunctionReturn(PETSC_SUCCESS); } /* DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of `DMPolytopeType` `ct` with localized coordinates `x`, generate localized coordinates `xr` for subcell `r` of type `rct`. Not Collective Input Parameters: + tr - The `DMPlexTransform` . ct - The type of the parent cell . rct - The type of the produced cell . r - The index of the produced cell - x - The localized coordinates for the parent cell Output Parameter: . xr - The localized coordinates for the produced cell Level: developer .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexCellRefinerSetCoordinates()` */ static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[]) { PetscFE fe = NULL; PetscInt cdim, v, *subcellV; PetscFunctionBegin; PetscCall(DMPlexTransformGetCoordinateFE(tr, ct, &fe)); PetscCall(DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV)); PetscCall(PetscFEGetNumComponents(fe, &cdim)); for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscCall(PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim])); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm) { DM dm, cdm, cdmCell, cdmNew, cdmCellNew; PetscSection coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew; Vec coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew; const PetscScalar *coords; PetscScalar *coordsNew; const PetscReal *maxCell, *Lstart, *L; PetscBool localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE; PetscInt dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p; PetscFunctionBegin; // Need to clear the DMField for coordinates PetscCall(DMSetCoordinateField(rdm, NULL)); PetscCall(DMPlexTransformGetDM(tr, &dm)); PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_SetCoordinates, tr, dm, 0, 0)); PetscCall(DMGetCoordinateDM(dm, &cdm)); PetscCall(DMGetCellCoordinateDM(dm, &cdmCell)); PetscCall(DMGetCoordinatesLocalized(dm, &localized)); PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L)); if (localized) { /* Localize coordinates of new vertices */ localizeVertices = PETSC_TRUE; /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */ if (!maxCell) localizeCells = PETSC_TRUE; } PetscCall(DMGetCoordinateSection(dm, &coordSection)); PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &dEo)); if (maxCell) { PetscReal maxCellNew[3]; for (d = 0; d < dEo; ++d) maxCellNew[d] = maxCell[d] / 2.0; PetscCall(DMSetPeriodicity(rdm, maxCellNew, Lstart, L)); } PetscCall(DMGetCoordinateDim(rdm, &dE)); PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew)); PetscCall(PetscSectionSetNumFields(coordSectionNew, 1)); PetscCall(PetscSectionSetFieldComponents(coordSectionNew, 0, dE)); PetscCall(DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew)); PetscCall(PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew)); /* Localization should be inherited */ /* Stefano calculates parent cells for each new cell for localization */ /* Localized cells need coordinates of closure */ for (v = vStartNew; v < vEndNew; ++v) { PetscCall(PetscSectionSetDof(coordSectionNew, v, dE)); PetscCall(PetscSectionSetFieldDof(coordSectionNew, v, 0, dE)); } PetscCall(PetscSectionSetUp(coordSectionNew)); PetscCall(DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew)); if (localizeCells) { PetscCall(DMGetCoordinateDM(rdm, &cdmNew)); PetscCall(DMClone(cdmNew, &cdmCellNew)); PetscCall(DMSetCellCoordinateDM(rdm, cdmCellNew)); PetscCall(DMDestroy(&cdmCellNew)); PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew)); PetscCall(PetscSectionSetNumFields(coordSectionCellNew, 1)); PetscCall(PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE)); PetscCall(DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew)); PetscCall(PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew)); PetscCall(DMGetCellCoordinateSection(dm, &coordSectionCell)); PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); for (c = cStart; c < cEnd; ++c) { PetscInt dof; PetscCall(PetscSectionGetDof(coordSectionCell, c, &dof)); if (dof) { DMPolytopeType ct; DMPolytopeType *rct; PetscInt *rsize, *rcone, *rornt; PetscInt dim, cNew, Nct, n, r; PetscCall(DMPlexGetCellType(dm, c, &ct)); dim = DMPolytopeTypeGetDim(ct); PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); /* This allows for different cell types */ for (n = 0; n < Nct; ++n) { if (dim != DMPolytopeTypeGetDim(rct[n])) continue; for (r = 0; r < rsize[n]; ++r) { PetscInt *closure = NULL; PetscInt clSize, cl, Nv = 0; PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew)); PetscCall(DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure)); for (cl = 0; cl < clSize * 2; cl += 2) { if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv; } PetscCall(DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure)); PetscCall(PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE)); PetscCall(PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE)); } } } } PetscCall(PetscSectionSetUp(coordSectionCellNew)); PetscCall(DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew)); } PetscCall(DMViewFromOptions(dm, NULL, "-coarse_dm_view")); { VecType vtype; PetscInt coordSizeNew, bs; const char *name; PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal)); PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalNew)); PetscCall(PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew)); PetscCall(VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE)); PetscCall(PetscObjectGetName((PetscObject)coordsLocal, &name)); PetscCall(PetscObjectSetName((PetscObject)coordsLocalNew, name)); PetscCall(VecGetBlockSize(coordsLocal, &bs)); PetscCall(VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE)); PetscCall(VecGetType(coordsLocal, &vtype)); PetscCall(VecSetType(coordsLocalNew, vtype)); } PetscCall(VecGetArrayRead(coordsLocal, &coords)); PetscCall(VecGetArray(coordsLocalNew, &coordsNew)); PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); /* First set coordinates for vertices */ for (p = pStart; p < pEnd; ++p) { DMPolytopeType ct; DMPolytopeType *rct; PetscInt *rsize, *rcone, *rornt; PetscInt Nct, n, r; PetscBool hasVertex = PETSC_FALSE; PetscCall(DMPlexGetCellType(dm, p, &ct)); PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); for (n = 0; n < Nct; ++n) { if (rct[n] == DM_POLYTOPE_POINT) { hasVertex = PETSC_TRUE; break; } } if (hasVertex) { const PetscScalar *icoords = NULL; const PetscScalar *array = NULL; PetscScalar *pcoords = NULL; PetscBool isDG; PetscInt Nc, Nv, v, d; PetscCall(DMPlexGetCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords)); icoords = pcoords; Nv = Nc / dEo; if (ct != DM_POLYTOPE_POINT) { if (localizeVertices && maxCell) { PetscScalar anchor[3]; for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d]; for (v = 0; v < Nv; ++v) PetscCall(DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo])); } } for (n = 0; n < Nct; ++n) { if (rct[n] != DM_POLYTOPE_POINT) continue; for (r = 0; r < rsize[n]; ++r) { PetscScalar vcoords[3]; PetscInt vNew, off; PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew)); PetscCall(PetscSectionGetOffset(coordSectionNew, vNew, &off)); PetscCall(DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords)); PetscCall(DMSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off])); } } PetscCall(DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords)); } } PetscCall(VecRestoreArrayRead(coordsLocal, &coords)); PetscCall(VecRestoreArray(coordsLocalNew, &coordsNew)); PetscCall(DMSetCoordinatesLocal(rdm, coordsLocalNew)); PetscCall(VecDestroy(&coordsLocalNew)); PetscCall(PetscSectionDestroy(&coordSectionNew)); /* Then set coordinates for cells by localizing */ if (!localizeCells) PetscCall(DMLocalizeCoordinates(rdm)); else { VecType vtype; PetscInt coordSizeNew, bs; const char *name; PetscCall(DMGetCellCoordinatesLocal(dm, &coordsLocalCell)); PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew)); PetscCall(PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew)); PetscCall(VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE)); PetscCall(PetscObjectGetName((PetscObject)coordsLocalCell, &name)); PetscCall(PetscObjectSetName((PetscObject)coordsLocalCellNew, name)); PetscCall(VecGetBlockSize(coordsLocalCell, &bs)); PetscCall(VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE)); PetscCall(VecGetType(coordsLocalCell, &vtype)); PetscCall(VecSetType(coordsLocalCellNew, vtype)); PetscCall(VecGetArrayRead(coordsLocalCell, &coords)); PetscCall(VecGetArray(coordsLocalCellNew, &coordsNew)); for (p = pStart; p < pEnd; ++p) { DMPolytopeType ct; DMPolytopeType *rct; PetscInt *rsize, *rcone, *rornt; PetscInt dof = 0, Nct, n, r; PetscCall(DMPlexGetCellType(dm, p, &ct)); PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(coordSectionCell, p, &dof)); if (dof) { const PetscScalar *pcoords; PetscCall(DMPlexPointLocalRead(cdmCell, p, coords, &pcoords)); for (n = 0; n < Nct; ++n) { const PetscInt Nr = rsize[n]; if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue; for (r = 0; r < Nr; ++r) { PetscInt pNew, offNew; /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger cell to the ones it produces. */ PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); PetscCall(PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew)); PetscCall(DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew])); } } } } PetscCall(VecRestoreArrayRead(coordsLocalCell, &coords)); PetscCall(VecRestoreArray(coordsLocalCellNew, &coordsNew)); PetscCall(DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew)); PetscCall(VecDestroy(&coordsLocalCellNew)); PetscCall(PetscSectionDestroy(&coordSectionCellNew)); } PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_SetCoordinates, tr, dm, 0, 0)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexTransformApply - Execute the transformation, producing another `DM` Collective Input Parameters: + tr - The `DMPlexTransform` object - dm - The original `DM` Output Parameter: . trdm - The transformed `DM` Level: intermediate Options Database Keys: + -dm_plex_transform_label_match_strata - Only label points of the same stratum as the producing point . -dm_plex_transform_label_replica_inc - Increment for the label value to be multiplied by the replica number - -dm_plex_transform_active - Name for active mesh label .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformCreate()`, `DMPlexTransformSetDM()` @*/ PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *trdm) { DM rdm; DMPlexInterpolatedFlag interp; PetscInt pStart, pEnd; PetscFunctionBegin; PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); PetscValidHeaderSpecific(dm, DM_CLASSID, 2); PetscAssertPointer(trdm, 3); PetscCall(PetscLogEventBegin(DMPLEXTRANSFORM_Apply, tr, dm, 0, 0)); PetscCall(DMPlexTransformSetDM(tr, dm)); PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &rdm)); PetscCall(DMSetType(rdm, DMPLEX)); PetscCall(DMPlexTransformSetDimensions(tr, dm, rdm)); /* Calculate number of new points of each depth */ PetscCall(DMPlexIsInterpolatedCollective(dm, &interp)); PetscCheck(interp == DMPLEX_INTERPOLATED_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement"); /* Step 1: Set chart */ PetscCall(DMPlexTransformGetChart(tr, &pStart, &pEnd)); PetscCall(DMPlexSetChart(rdm, pStart, pEnd)); /* Step 2: Set cone/support sizes (automatically stratifies) */ PetscCall(DMPlexTransformSetConeSizes(tr, rdm)); /* Step 3: Setup refined DM */ PetscCall(DMSetUp(rdm)); /* Step 4: Set cones and supports (automatically symmetrizes) */ PetscCall(DMPlexTransformSetCones(tr, rdm)); /* Step 5: Create pointSF */ PetscCall(DMPlexTransformCreateSF(tr, rdm)); /* Step 6: Create labels */ PetscCall(DMPlexTransformCreateLabels(tr, rdm)); /* Step 7: Set coordinates */ PetscCall(DMPlexTransformSetCoordinates(tr, rdm)); PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, rdm)); // If the original DM was configured from options, the transformed DM should be as well rdm->setfromoptionscalled = dm->setfromoptionscalled; PetscCall(PetscLogEventEnd(DMPLEXTRANSFORM_Apply, tr, dm, 0, 0)); *trdm = rdm; PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm) { DMPlexTransform tr; DM cdm, rcdm; const char *prefix; PetscBool save; PetscFunctionBegin; PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr)); PetscCall(PetscObjectSetName((PetscObject)tr, "Adapt Label Transform")); PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix)); PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix)); PetscCall(DMPlexTransformSetDM(tr, dm)); PetscCall(DMPlexTransformSetFromOptions(tr)); if (adaptLabel) PetscCall(DMPlexTransformSetActive(tr, adaptLabel)); PetscCall(DMPlexTransformSetUp(tr)); PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view")); PetscCall(DMPlexTransformApply(tr, dm, rdm)); PetscCall(DMCopyDisc(dm, *rdm)); PetscCall(DMGetCoordinateDM(dm, &cdm)); PetscCall(DMGetCoordinateDM(*rdm, &rcdm)); PetscCall(DMCopyDisc(cdm, rcdm)); PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm)); PetscCall(DMCopyDisc(dm, *rdm)); PetscCall(DMPlexGetSaveTransform(dm, &save)); if (save) PetscCall(DMPlexSetTransform(*rdm, tr)); PetscCall(DMPlexTransformDestroy(&tr)); ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation; PetscFunctionReturn(PETSC_SUCCESS); }