1 #include "petsc/private/petscimpl.h" 2 #include "petscdmplex.h" 3 #include "petscdmplextransform.h" 4 #include "petscdmplextransformtypes.h" 5 #include "petscerror.h" 6 #include "petscsystypes.h" 7 #include <petsc/private/dmplextransformimpl.h> /*I "petscdmplextransform.h" I*/ 8 9 #include <petsc/private/petscfeimpl.h> /* For PetscFEInterpolate_Static() */ 10 11 PetscClassId DMPLEXTRANSFORM_CLASSID; 12 13 PetscFunctionList DMPlexTransformList = NULL; 14 PetscBool DMPlexTransformRegisterAllCalled = PETSC_FALSE; 15 16 /* Construct cell type order since we must loop over cell types in depth order */ 17 static PetscErrorCode DMPlexCreateCellTypeOrder_Internal(PetscInt dim, PetscInt *ctOrder[], PetscInt *ctOrderInv[]) 18 { 19 PetscInt *ctO, *ctOInv; 20 PetscInt c, d, off = 0; 21 22 PetscFunctionBegin; 23 PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctO, DM_NUM_POLYTOPES + 1, &ctOInv)); 24 for (d = 3; d >= dim; --d) { 25 for (c = 0; c <= DM_NUM_POLYTOPES; ++c) { 26 if (DMPolytopeTypeGetDim((DMPolytopeType)c) != d || c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) continue; 27 ctO[off++] = c; 28 } 29 } 30 if (dim != 0) { 31 for (c = 0; c <= DM_NUM_POLYTOPES; ++c) { 32 if (DMPolytopeTypeGetDim((DMPolytopeType)c) != 0) continue; 33 ctO[off++] = c; 34 } 35 } 36 for (d = dim - 1; d > 0; --d) { 37 for (c = 0; c <= DM_NUM_POLYTOPES; ++c) { 38 if (DMPolytopeTypeGetDim((DMPolytopeType)c) != d || c == DM_POLYTOPE_UNKNOWN_CELL || c == DM_POLYTOPE_UNKNOWN_FACE) continue; 39 ctO[off++] = c; 40 } 41 } 42 for (c = 0; c <= DM_NUM_POLYTOPES; ++c) { 43 if (DMPolytopeTypeGetDim((DMPolytopeType)c) >= 0 && c != DM_POLYTOPE_UNKNOWN_CELL && c != DM_POLYTOPE_UNKNOWN_FACE) continue; 44 ctO[off++] = c; 45 } 46 PetscCheck(off == DM_NUM_POLYTOPES + 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid offset %" PetscInt_FMT " for cell type order", off); 47 for (c = 0; c <= DM_NUM_POLYTOPES; ++c) ctOInv[ctO[c]] = c; 48 *ctOrder = ctO; 49 *ctOrderInv = ctOInv; 50 PetscFunctionReturn(PETSC_SUCCESS); 51 } 52 53 /*@C 54 DMPlexTransformRegister - Adds a new transform component implementation 55 56 Not Collective 57 58 Input Parameters: 59 + name - The name of a new user-defined creation routine 60 - create_func - The creation routine 61 62 Example Usage: 63 .vb 64 DMPlexTransformRegister("my_transform", MyTransformCreate); 65 .ve 66 67 Then, your transform type can be chosen with the procedural interface via 68 .vb 69 DMPlexTransformCreate(MPI_Comm, DMPlexTransform *); 70 DMPlexTransformSetType(DMPlexTransform, "my_transform"); 71 .ve 72 or at runtime via the option 73 .vb 74 -dm_plex_transform_type my_transform 75 .ve 76 77 Level: advanced 78 79 Note: 80 `DMPlexTransformRegister()` may be called multiple times to add several user-defined transforms 81 82 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformRegisterAll()`, `DMPlexTransformRegisterDestroy()` 83 @*/ 84 PetscErrorCode DMPlexTransformRegister(const char name[], PetscErrorCode (*create_func)(DMPlexTransform)) 85 { 86 PetscFunctionBegin; 87 PetscCall(DMInitializePackage()); 88 PetscCall(PetscFunctionListAdd(&DMPlexTransformList, name, create_func)); 89 PetscFunctionReturn(PETSC_SUCCESS); 90 } 91 92 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Filter(DMPlexTransform); 93 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform); 94 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToBox(DMPlexTransform); 95 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform); 96 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform); 97 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform); 98 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_1D(DMPlexTransform); 99 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform); 100 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Cohesive(DMPlexTransform); 101 102 /*@C 103 DMPlexTransformRegisterAll - Registers all of the transform components in the `DM` package. 104 105 Not Collective 106 107 Level: advanced 108 109 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `DMRegisterAll()`, `DMPlexTransformRegisterDestroy()` 110 @*/ 111 PetscErrorCode DMPlexTransformRegisterAll(void) 112 { 113 PetscFunctionBegin; 114 if (DMPlexTransformRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS); 115 DMPlexTransformRegisterAllCalled = PETSC_TRUE; 116 117 PetscCall(DMPlexTransformRegister(DMPLEXTRANSFORMFILTER, DMPlexTransformCreate_Filter)); 118 PetscCall(DMPlexTransformRegister(DMPLEXREFINEREGULAR, DMPlexTransformCreate_Regular)); 119 PetscCall(DMPlexTransformRegister(DMPLEXREFINETOBOX, DMPlexTransformCreate_ToBox)); 120 PetscCall(DMPlexTransformRegister(DMPLEXREFINEALFELD, DMPlexTransformCreate_Alfeld)); 121 PetscCall(DMPlexTransformRegister(DMPLEXREFINEBOUNDARYLAYER, DMPlexTransformCreate_BL)); 122 PetscCall(DMPlexTransformRegister(DMPLEXREFINESBR, DMPlexTransformCreate_SBR)); 123 PetscCall(DMPlexTransformRegister(DMPLEXREFINE1D, DMPlexTransformCreate_1D)); 124 PetscCall(DMPlexTransformRegister(DMPLEXEXTRUDETYPE, DMPlexTransformCreate_Extrude)); 125 PetscCall(DMPlexTransformRegister(DMPLEXCOHESIVEEXTRUDE, DMPlexTransformCreate_Cohesive)); 126 PetscFunctionReturn(PETSC_SUCCESS); 127 } 128 129 /*@C 130 DMPlexTransformRegisterDestroy - This function destroys the registered `DMPlexTransformType`. It is called from `PetscFinalize()`. 131 132 Not collective 133 134 Level: developer 135 136 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRegisterAll()`, `DMPlexTransformType`, `PetscInitialize()` 137 @*/ 138 PetscErrorCode DMPlexTransformRegisterDestroy(void) 139 { 140 PetscFunctionBegin; 141 PetscCall(PetscFunctionListDestroy(&DMPlexTransformList)); 142 DMPlexTransformRegisterAllCalled = PETSC_FALSE; 143 PetscFunctionReturn(PETSC_SUCCESS); 144 } 145 146 /*@ 147 DMPlexTransformCreate - Creates an empty transform object. The type can then be set with `DMPlexTransformSetType()`. 148 149 Collective 150 151 Input Parameter: 152 . comm - The communicator for the transform object 153 154 Output Parameter: 155 . tr - The transform object 156 157 Level: beginner 158 159 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPLEXREFINEREGULAR`, `DMPLEXTRANSFORMFILTER` 160 @*/ 161 PetscErrorCode DMPlexTransformCreate(MPI_Comm comm, DMPlexTransform *tr) 162 { 163 DMPlexTransform t; 164 165 PetscFunctionBegin; 166 PetscAssertPointer(tr, 2); 167 *tr = NULL; 168 PetscCall(DMInitializePackage()); 169 170 PetscCall(PetscHeaderCreate(t, DMPLEXTRANSFORM_CLASSID, "DMPlexTransform", "Mesh Transform", "DMPlexTransform", comm, DMPlexTransformDestroy, DMPlexTransformView)); 171 t->setupcalled = PETSC_FALSE; 172 PetscCall(PetscCalloc2(DM_NUM_POLYTOPES, &t->coordFE, DM_NUM_POLYTOPES, &t->refGeom)); 173 *tr = t; 174 PetscFunctionReturn(PETSC_SUCCESS); 175 } 176 177 /*@ 178 DMPlexTransformSetType - Sets the particular implementation for a transform. 179 180 Collective 181 182 Input Parameters: 183 + tr - The transform 184 - method - The name of the transform type 185 186 Options Database Key: 187 . -dm_plex_transform_type <type> - Sets the transform type; see `DMPlexTransformType` 188 189 Level: intermediate 190 191 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformGetType()`, `DMPlexTransformCreate()` 192 @*/ 193 PetscErrorCode DMPlexTransformSetType(DMPlexTransform tr, DMPlexTransformType method) 194 { 195 PetscErrorCode (*r)(DMPlexTransform); 196 PetscBool match; 197 198 PetscFunctionBegin; 199 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 200 PetscCall(PetscObjectTypeCompare((PetscObject)tr, method, &match)); 201 if (match) PetscFunctionReturn(PETSC_SUCCESS); 202 203 PetscCall(DMPlexTransformRegisterAll()); 204 PetscCall(PetscFunctionListFind(DMPlexTransformList, method, &r)); 205 PetscCheck(r, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMPlexTransform type: %s", method); 206 207 PetscTryTypeMethod(tr, destroy); 208 PetscCall(PetscMemzero(tr->ops, sizeof(*tr->ops))); 209 PetscCall(PetscObjectChangeTypeName((PetscObject)tr, method)); 210 PetscCall((*r)(tr)); 211 PetscFunctionReturn(PETSC_SUCCESS); 212 } 213 214 /*@ 215 DMPlexTransformGetType - Gets the type name (as a string) from the transform. 216 217 Not Collective 218 219 Input Parameter: 220 . tr - The `DMPlexTransform` 221 222 Output Parameter: 223 . type - The `DMPlexTransformType` name 224 225 Level: intermediate 226 227 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPlexTransformCreate()` 228 @*/ 229 PetscErrorCode DMPlexTransformGetType(DMPlexTransform tr, DMPlexTransformType *type) 230 { 231 PetscFunctionBegin; 232 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 233 PetscAssertPointer(type, 2); 234 PetscCall(DMPlexTransformRegisterAll()); 235 *type = ((PetscObject)tr)->type_name; 236 PetscFunctionReturn(PETSC_SUCCESS); 237 } 238 239 static PetscErrorCode DMPlexTransformView_Ascii(DMPlexTransform tr, PetscViewer v) 240 { 241 PetscViewerFormat format; 242 243 PetscFunctionBegin; 244 PetscCall(PetscViewerGetFormat(v, &format)); 245 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 246 const PetscInt *trTypes = NULL; 247 IS trIS; 248 PetscInt cols = 8; 249 PetscInt Nrt = 8, f, g; 250 251 if (tr->trType) PetscCall(DMLabelView(tr->trType, v)); 252 PetscCall(PetscViewerASCIIPrintf(v, "Source Starts\n")); 253 for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g])); 254 PetscCall(PetscViewerASCIIPrintf(v, "\n")); 255 for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStart[f])); 256 PetscCall(PetscViewerASCIIPrintf(v, "\n")); 257 PetscCall(PetscViewerASCIIPrintf(v, "Target Starts\n")); 258 for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g])); 259 PetscCall(PetscViewerASCIIPrintf(v, "\n")); 260 for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStartNew[f])); 261 PetscCall(PetscViewerASCIIPrintf(v, "\n")); 262 263 if (tr->trType) { 264 PetscCall(DMLabelGetNumValues(tr->trType, &Nrt)); 265 PetscCall(DMLabelGetValueIS(tr->trType, &trIS)); 266 PetscCall(ISGetIndices(trIS, &trTypes)); 267 } 268 PetscCall(PetscViewerASCIIPrintf(v, "Offsets\n")); 269 PetscCall(PetscViewerASCIIPrintf(v, " ")); 270 for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g])); 271 PetscCall(PetscViewerASCIIPrintf(v, "\n")); 272 for (f = 0; f < Nrt; ++f) { 273 PetscCall(PetscViewerASCIIPrintf(v, "%2" PetscInt_FMT " |", trTypes ? trTypes[f] : f)); 274 for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->offset[f * DM_NUM_POLYTOPES + g])); 275 PetscCall(PetscViewerASCIIPrintf(v, " |\n")); 276 } 277 if (tr->trType) { 278 PetscCall(ISRestoreIndices(trIS, &trTypes)); 279 PetscCall(ISDestroy(&trIS)); 280 } 281 } 282 PetscFunctionReturn(PETSC_SUCCESS); 283 } 284 285 /*@ 286 DMPlexTransformView - Views a `DMPlexTransform` 287 288 Collective 289 290 Input Parameters: 291 + tr - the `DMPlexTransform` object to view 292 - v - the viewer 293 294 Level: beginner 295 296 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `PetscViewer`, `DMPlexTransformDestroy()`, `DMPlexTransformCreate()` 297 @*/ 298 PetscErrorCode DMPlexTransformView(DMPlexTransform tr, PetscViewer v) 299 { 300 PetscBool isascii; 301 302 PetscFunctionBegin; 303 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 304 if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tr), &v)); 305 PetscValidHeaderSpecific(v, PETSC_VIEWER_CLASSID, 2); 306 PetscCheckSameComm(tr, 1, v, 2); 307 PetscCall(PetscViewerCheckWritable(v)); 308 PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tr, v)); 309 PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii)); 310 if (isascii) PetscCall(DMPlexTransformView_Ascii(tr, v)); 311 PetscTryTypeMethod(tr, view, v); 312 PetscFunctionReturn(PETSC_SUCCESS); 313 } 314 315 /*@ 316 DMPlexTransformSetFromOptions - Sets parameters in a transform from values in the options database 317 318 Collective 319 320 Input Parameter: 321 . tr - the `DMPlexTransform` object to set options for 322 323 Options Database Keys: 324 + -dm_plex_transform_type - Set the transform type, e.g. refine_regular 325 . -dm_plex_transform_label_match_strata - Only label points of the same stratum as the producing point 326 . -dm_plex_transform_label_replica_inc <inc> - Increment for the label value to be multiplied by the replica number, so that the new label value is oldValue + r * inc 327 . -dm_plex_transform_active <name> - Name for active mesh label 328 - -dm_plex_transform_active_values <v0,v1,...> - Values in the active label 329 330 Level: intermediate 331 332 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()` 333 @*/ 334 PetscErrorCode DMPlexTransformSetFromOptions(DMPlexTransform tr) 335 { 336 char typeName[1024], active[PETSC_MAX_PATH_LEN]; 337 const char *defName = DMPLEXREFINEREGULAR; 338 PetscBool flg, match; 339 340 PetscFunctionBegin; 341 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 342 PetscObjectOptionsBegin((PetscObject)tr); 343 PetscCall(PetscOptionsFList("-dm_plex_transform_type", "DMPlexTransform", "DMPlexTransformSetType", DMPlexTransformList, defName, typeName, 1024, &flg)); 344 if (flg) PetscCall(DMPlexTransformSetType(tr, typeName)); 345 else if (!((PetscObject)tr)->type_name) PetscCall(DMPlexTransformSetType(tr, defName)); 346 PetscCall(PetscOptionsBool("-dm_plex_transform_label_match_strata", "Only label points of the same stratum as the producing point", "", tr->labelMatchStrata, &match, &flg)); 347 if (flg) PetscCall(DMPlexTransformSetMatchStrata(tr, match)); 348 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)); 349 PetscCall(PetscOptionsString("-dm_plex_transform_active", "Name for active mesh label", "DMPlexTransformSetActive", active, active, sizeof(active), &flg)); 350 if (flg) { 351 DM dm; 352 DMLabel label; 353 PetscInt values[16]; 354 PetscInt n = 16; 355 356 PetscCall(DMPlexTransformGetDM(tr, &dm)); 357 PetscCall(DMGetLabel(dm, active, &label)); 358 PetscCall(PetscOptionsIntArray("-dm_plex_transform_active_values", "The label values to be active", "DMPlexTransformSetActive", values, &n, &flg)); 359 if (flg && n) { 360 DMLabel newlabel; 361 362 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Active", &newlabel)); 363 for (PetscInt i = 0; i < n; ++i) { 364 IS is; 365 366 PetscCall(DMLabelGetStratumIS(label, values[i], &is)); 367 PetscCall(DMLabelInsertIS(newlabel, is, values[i])); 368 PetscCall(ISDestroy(&is)); 369 } 370 PetscCall(DMPlexTransformSetActive(tr, newlabel)); 371 PetscCall(DMLabelDestroy(&newlabel)); 372 } else { 373 PetscCall(DMPlexTransformSetActive(tr, label)); 374 } 375 } 376 PetscTryTypeMethod(tr, setfromoptions, PetscOptionsObject); 377 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 378 PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)tr, PetscOptionsObject)); 379 PetscOptionsEnd(); 380 PetscFunctionReturn(PETSC_SUCCESS); 381 } 382 383 /*@ 384 DMPlexTransformDestroy - Destroys a `DMPlexTransform` 385 386 Collective 387 388 Input Parameter: 389 . tr - the transform object to destroy 390 391 Level: beginner 392 393 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()` 394 @*/ 395 PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr) 396 { 397 PetscInt c; 398 399 PetscFunctionBegin; 400 if (!*tr) PetscFunctionReturn(PETSC_SUCCESS); 401 PetscValidHeaderSpecific(*tr, DMPLEXTRANSFORM_CLASSID, 1); 402 if (--((PetscObject)*tr)->refct > 0) { 403 *tr = NULL; 404 PetscFunctionReturn(PETSC_SUCCESS); 405 } 406 407 PetscTryTypeMethod(*tr, destroy); 408 PetscCall(DMDestroy(&(*tr)->dm)); 409 PetscCall(DMLabelDestroy(&(*tr)->active)); 410 PetscCall(DMLabelDestroy(&(*tr)->trType)); 411 PetscCall(PetscFree2((*tr)->ctOrderOld, (*tr)->ctOrderInvOld)); 412 PetscCall(PetscFree2((*tr)->ctOrderNew, (*tr)->ctOrderInvNew)); 413 PetscCall(PetscFree2((*tr)->ctStart, (*tr)->ctStartNew)); 414 PetscCall(PetscFree((*tr)->offset)); 415 PetscCall(PetscFree2((*tr)->depthStart, (*tr)->depthEnd)); 416 for (c = 0; c < DM_NUM_POLYTOPES; ++c) { 417 PetscCall(PetscFEDestroy(&(*tr)->coordFE[c])); 418 PetscCall(PetscFEGeomDestroy(&(*tr)->refGeom[c])); 419 } 420 if ((*tr)->trVerts) { 421 for (c = 0; c < DM_NUM_POLYTOPES; ++c) { 422 DMPolytopeType *rct; 423 PetscInt *rsize, *rcone, *rornt, Nct, n, r; 424 425 if (DMPolytopeTypeGetDim((DMPolytopeType)c) > 0 && c != DM_POLYTOPE_UNKNOWN_CELL && c != DM_POLYTOPE_UNKNOWN_FACE) { 426 PetscCall(DMPlexTransformCellTransform(*tr, (DMPolytopeType)c, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 427 for (n = 0; n < Nct; ++n) { 428 if (rct[n] == DM_POLYTOPE_POINT) continue; 429 for (r = 0; r < rsize[n]; ++r) PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]][r])); 430 PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]])); 431 } 432 } 433 PetscCall(PetscFree((*tr)->trSubVerts[c])); 434 PetscCall(PetscFree((*tr)->trVerts[c])); 435 } 436 } 437 PetscCall(PetscFree3((*tr)->trNv, (*tr)->trVerts, (*tr)->trSubVerts)); 438 PetscCall(PetscFree2((*tr)->coordFE, (*tr)->refGeom)); 439 /* We do not destroy (*dm)->data here so that we can reference count backend objects */ 440 PetscCall(PetscHeaderDestroy(tr)); 441 PetscFunctionReturn(PETSC_SUCCESS); 442 } 443 444 static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrderOld[], PetscInt ctStart[], PetscInt **offset) 445 { 446 DMLabel trType = tr->trType; 447 PetscInt c, cN, *off; 448 449 PetscFunctionBegin; 450 if (trType) { 451 DM dm; 452 IS rtIS; 453 const PetscInt *reftypes; 454 PetscInt Nrt, r; 455 456 PetscCall(DMPlexTransformGetDM(tr, &dm)); 457 PetscCall(DMLabelGetNumValues(trType, &Nrt)); 458 PetscCall(DMLabelGetValueIS(trType, &rtIS)); 459 PetscCall(ISGetIndices(rtIS, &reftypes)); 460 PetscCall(PetscCalloc1(Nrt * DM_NUM_POLYTOPES, &off)); 461 for (r = 0; r < Nrt; ++r) { 462 const PetscInt rt = reftypes[r]; 463 IS rtIS; 464 const PetscInt *points; 465 DMPolytopeType ct; 466 PetscInt np, p; 467 468 PetscCall(DMLabelGetStratumIS(trType, rt, &rtIS)); 469 PetscCall(ISGetLocalSize(rtIS, &np)); 470 PetscCall(ISGetIndices(rtIS, &points)); 471 if (!np) continue; 472 p = points[0]; 473 PetscCall(ISRestoreIndices(rtIS, &points)); 474 PetscCall(ISDestroy(&rtIS)); 475 PetscCall(DMPlexGetCellType(dm, p, &ct)); 476 for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) { 477 const DMPolytopeType ctNew = (DMPolytopeType)cN; 478 DMPolytopeType *rct; 479 PetscInt *rsize, *cone, *ornt; 480 PetscInt Nct, n, s; 481 482 if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) { 483 off[r * DM_NUM_POLYTOPES + ctNew] = -1; 484 break; 485 } 486 off[r * DM_NUM_POLYTOPES + ctNew] = 0; 487 for (s = 0; s <= r; ++s) { 488 const PetscInt st = reftypes[s]; 489 DMPolytopeType sct; 490 PetscInt q, qrt; 491 492 PetscCall(DMLabelGetStratumIS(trType, st, &rtIS)); 493 PetscCall(ISGetLocalSize(rtIS, &np)); 494 PetscCall(ISGetIndices(rtIS, &points)); 495 if (!np) continue; 496 q = points[0]; 497 PetscCall(ISRestoreIndices(rtIS, &points)); 498 PetscCall(ISDestroy(&rtIS)); 499 PetscCall(DMPlexGetCellType(dm, q, &sct)); 500 PetscCall(DMPlexTransformCellTransform(tr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt)); 501 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); 502 if (st == rt) { 503 for (n = 0; n < Nct; ++n) 504 if (rct[n] == ctNew) break; 505 if (n == Nct) off[r * DM_NUM_POLYTOPES + ctNew] = -1; 506 break; 507 } 508 for (n = 0; n < Nct; ++n) { 509 if (rct[n] == ctNew) { 510 PetscInt sn; 511 512 PetscCall(DMLabelGetStratumSize(trType, st, &sn)); 513 off[r * DM_NUM_POLYTOPES + ctNew] += sn * rsize[n]; 514 } 515 } 516 } 517 } 518 } 519 PetscCall(ISRestoreIndices(rtIS, &reftypes)); 520 PetscCall(ISDestroy(&rtIS)); 521 } else { 522 PetscCall(PetscCalloc1(DM_NUM_POLYTOPES * DM_NUM_POLYTOPES, &off)); 523 for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) { 524 const DMPolytopeType ct = (DMPolytopeType)c; 525 for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) { 526 const DMPolytopeType ctNew = (DMPolytopeType)cN; 527 DMPolytopeType *rct; 528 PetscInt *rsize, *cone, *ornt; 529 PetscInt Nct, n, i; 530 531 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) { 532 off[ct * DM_NUM_POLYTOPES + ctNew] = -1; 533 continue; 534 } 535 off[ct * DM_NUM_POLYTOPES + ctNew] = 0; 536 for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) { 537 const DMPolytopeType ict = (DMPolytopeType)ctOrderOld[i]; 538 const DMPolytopeType ictn = (DMPolytopeType)ctOrderOld[i + 1]; 539 540 PetscCall(DMPlexTransformCellTransform(tr, ict, PETSC_DETERMINE, NULL, &Nct, &rct, &rsize, &cone, &ornt)); 541 if (ict == ct) { 542 for (n = 0; n < Nct; ++n) 543 if (rct[n] == ctNew) break; 544 if (n == Nct) off[ct * DM_NUM_POLYTOPES + ctNew] = -1; 545 break; 546 } 547 for (n = 0; n < Nct; ++n) 548 if (rct[n] == ctNew) off[ct * DM_NUM_POLYTOPES + ctNew] += (ctStart[ictn] - ctStart[ict]) * rsize[n]; 549 } 550 } 551 } 552 } 553 *offset = off; 554 PetscFunctionReturn(PETSC_SUCCESS); 555 } 556 557 /*@ 558 DMPlexTransformSetUp - Create the tables that drive the transform 559 560 Input Parameter: 561 . tr - The `DMPlexTransform` object 562 563 Level: intermediate 564 565 .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` 566 @*/ 567 PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr) 568 { 569 DM dm; 570 DMPolytopeType ctCell; 571 PetscInt pStart, pEnd, p, c, celldim = 0; 572 573 PetscFunctionBegin; 574 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 575 if (tr->setupcalled) PetscFunctionReturn(PETSC_SUCCESS); 576 PetscTryTypeMethod(tr, setup); 577 PetscCall(DMPlexTransformGetDM(tr, &dm)); 578 PetscCall(DMSetSnapToGeomModel(dm, NULL)); 579 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 580 if (pEnd > pStart) { 581 // Ignore cells hanging off of embedded surfaces 582 PetscInt c = pStart; 583 584 ctCell = DM_POLYTOPE_FV_GHOST; 585 while (DMPolytopeTypeGetDim(ctCell) < 0) PetscCall(DMPlexGetCellType(dm, c++, &ctCell)); 586 } else { 587 PetscInt dim; 588 589 PetscCall(DMGetDimension(dm, &dim)); 590 switch (dim) { 591 case 0: 592 ctCell = DM_POLYTOPE_POINT; 593 break; 594 case 1: 595 ctCell = DM_POLYTOPE_SEGMENT; 596 break; 597 case 2: 598 ctCell = DM_POLYTOPE_TRIANGLE; 599 break; 600 case 3: 601 ctCell = DM_POLYTOPE_TETRAHEDRON; 602 break; 603 default: 604 ctCell = DM_POLYTOPE_UNKNOWN; 605 } 606 } 607 PetscCall(DMPlexCreateCellTypeOrder_Internal(DMPolytopeTypeGetDim(ctCell), &tr->ctOrderOld, &tr->ctOrderInvOld)); 608 for (p = pStart; p < pEnd; ++p) { 609 DMPolytopeType ct; 610 DMPolytopeType *rct; 611 PetscInt *rsize, *cone, *ornt; 612 PetscInt Nct, n; 613 614 PetscCall(DMPlexGetCellType(dm, p, &ct)); 615 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); 616 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt)); 617 for (n = 0; n < Nct; ++n) celldim = PetscMax(celldim, DMPolytopeTypeGetDim(rct[n])); 618 } 619 PetscCall(DMPlexCreateCellTypeOrder_Internal(celldim, &tr->ctOrderNew, &tr->ctOrderInvNew)); 620 /* Construct sizes and offsets for each cell type */ 621 if (!tr->ctStart) { 622 PetscInt *ctS, *ctSN, *ctC, *ctCN; 623 624 PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctS, DM_NUM_POLYTOPES + 1, &ctSN)); 625 PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctC, DM_NUM_POLYTOPES + 1, &ctCN)); 626 for (p = pStart; p < pEnd; ++p) { 627 DMPolytopeType ct; 628 DMPolytopeType *rct; 629 PetscInt *rsize, *cone, *ornt; 630 PetscInt Nct, n; 631 632 PetscCall(DMPlexGetCellType(dm, p, &ct)); 633 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); 634 ++ctC[ct]; 635 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt)); 636 for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n]; 637 } 638 for (c = 0; c < DM_NUM_POLYTOPES; ++c) { 639 const PetscInt cto = tr->ctOrderOld[c]; 640 const PetscInt cton = tr->ctOrderOld[c + 1]; 641 const PetscInt ctn = tr->ctOrderNew[c]; 642 const PetscInt ctnn = tr->ctOrderNew[c + 1]; 643 644 ctS[cton] = ctS[cto] + ctC[cto]; 645 ctSN[ctnn] = ctSN[ctn] + ctCN[ctn]; 646 } 647 PetscCall(PetscFree2(ctC, ctCN)); 648 tr->ctStart = ctS; 649 tr->ctStartNew = ctSN; 650 } 651 PetscCall(DMPlexTransformCreateOffset_Internal(tr, tr->ctOrderOld, tr->ctStart, &tr->offset)); 652 // Compute depth information 653 tr->depth = -1; 654 for (c = 0; c < DM_NUM_POLYTOPES; ++c) 655 if (tr->ctStartNew[tr->ctOrderNew[c + 1]] > tr->ctStartNew[tr->ctOrderNew[c]]) tr->depth = PetscMax(tr->depth, DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c])); 656 PetscCall(PetscMalloc2(tr->depth + 1, &tr->depthStart, tr->depth + 1, &tr->depthEnd)); 657 for (PetscInt d = 0; d <= tr->depth; ++d) { 658 tr->depthStart[d] = PETSC_INT_MAX; 659 tr->depthEnd[d] = -1; 660 } 661 for (c = 0; c < DM_NUM_POLYTOPES; ++c) { 662 const PetscInt dep = DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]); 663 664 if (tr->ctStartNew[tr->ctOrderNew[c + 1]] <= tr->ctStartNew[tr->ctOrderNew[c]]) continue; 665 tr->depthStart[dep] = PetscMin(tr->depthStart[dep], tr->ctStartNew[tr->ctOrderNew[c]]); 666 tr->depthEnd[dep] = PetscMax(tr->depthEnd[dep], tr->ctStartNew[tr->ctOrderNew[c + 1]]); 667 } 668 tr->setupcalled = PETSC_TRUE; 669 PetscFunctionReturn(PETSC_SUCCESS); 670 } 671 672 /*@ 673 DMPlexTransformGetDM - Get the base `DM` for the transform 674 675 Input Parameter: 676 . tr - The `DMPlexTransform` object 677 678 Output Parameter: 679 . dm - The original `DM` which will be transformed 680 681 Level: intermediate 682 683 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` 684 @*/ 685 PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm) 686 { 687 PetscFunctionBegin; 688 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 689 PetscAssertPointer(dm, 2); 690 *dm = tr->dm; 691 PetscFunctionReturn(PETSC_SUCCESS); 692 } 693 694 /*@ 695 DMPlexTransformSetDM - Set the base `DM` for the transform 696 697 Input Parameters: 698 + tr - The `DMPlexTransform` object 699 - dm - The original `DM` which will be transformed 700 701 Level: intermediate 702 703 Note: 704 The user does not typically call this, as it is called by `DMPlexTransformApply()`. 705 706 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetDM()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` 707 @*/ 708 PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm) 709 { 710 PetscFunctionBegin; 711 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 712 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 713 PetscCall(PetscObjectReference((PetscObject)dm)); 714 PetscCall(DMDestroy(&tr->dm)); 715 tr->dm = dm; 716 PetscFunctionReturn(PETSC_SUCCESS); 717 } 718 719 /*@ 720 DMPlexTransformGetActive - Get the `DMLabel` marking the active points for the transform 721 722 Input Parameter: 723 . tr - The `DMPlexTransform` object 724 725 Output Parameter: 726 . active - The `DMLabel` indicating which points will be transformed 727 728 Level: intermediate 729 730 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` 731 @*/ 732 PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active) 733 { 734 PetscFunctionBegin; 735 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 736 PetscAssertPointer(active, 2); 737 *active = tr->active; 738 PetscFunctionReturn(PETSC_SUCCESS); 739 } 740 741 /*@ 742 DMPlexTransformSetActive - Set the `DMLabel` marking the active points for the transform 743 744 Input Parameters: 745 + tr - The `DMPlexTransform` object 746 - active - The `DMLabel` indicating which points will be transformed 747 748 Level: intermediate 749 750 Note: 751 This only applies to transforms listed in [](plex_transform_table) that operate on a subset of the mesh. 752 753 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` 754 @*/ 755 PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active) 756 { 757 PetscFunctionBegin; 758 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 759 if (active) PetscValidHeaderSpecific(active, DMLABEL_CLASSID, 2); 760 PetscCall(PetscObjectReference((PetscObject)active)); 761 PetscCall(DMLabelDestroy(&tr->active)); 762 tr->active = active; 763 PetscFunctionReturn(PETSC_SUCCESS); 764 } 765 766 /*@ 767 DMPlexTransformGetTransformTypes - Get the `DMLabel` marking the transform type of each point for the transform 768 769 Input Parameter: 770 . tr - The `DMPlexTransform` object 771 772 Output Parameter: 773 . trType - The `DMLabel` indicating the transform type for each point 774 775 Level: intermediate 776 777 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexSetTransformType()`, `DMPlexTransformGetActive()`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` 778 @*/ 779 PetscErrorCode DMPlexTransformGetTransformTypes(DMPlexTransform tr, DMLabel *trType) 780 { 781 PetscFunctionBegin; 782 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 783 PetscAssertPointer(trType, 2); 784 *trType = tr->trType; 785 PetscFunctionReturn(PETSC_SUCCESS); 786 } 787 788 /*@ 789 DMPlexTransformSetTransformTypes - Set the `DMLabel` marking the transform type of each point for the transform 790 791 Input Parameters: 792 + tr - The `DMPlexTransform` object 793 - trType - The original `DM` which will be transformed 794 795 Level: intermediate 796 797 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetTransformTypes()`, `DMPlexTransformGetActive())`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` 798 @*/ 799 PetscErrorCode DMPlexTransformSetTransformTypes(DMPlexTransform tr, DMLabel trType) 800 { 801 PetscFunctionBegin; 802 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 803 if (trType) PetscValidHeaderSpecific(trType, DMLABEL_CLASSID, 2); 804 PetscCall(PetscObjectReference((PetscObject)trType)); 805 PetscCall(DMLabelDestroy(&tr->trType)); 806 tr->trType = trType; 807 PetscFunctionReturn(PETSC_SUCCESS); 808 } 809 810 static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe) 811 { 812 PetscFunctionBegin; 813 if (!tr->coordFE[ct]) { 814 PetscInt dim, cdim; 815 816 dim = DMPolytopeTypeGetDim(ct); 817 PetscCall(DMGetCoordinateDim(tr->dm, &cdim)); 818 PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, cdim, ct, 1, PETSC_DETERMINE, &tr->coordFE[ct])); 819 { 820 PetscDualSpace dsp; 821 PetscQuadrature quad; 822 DM K; 823 PetscFEGeom *cg; 824 PetscScalar *Xq; 825 PetscReal *xq, *wq; 826 PetscInt Nq, q; 827 828 PetscCall(DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq)); 829 PetscCall(PetscMalloc1(Nq * cdim, &xq)); 830 for (q = 0; q < Nq * cdim; ++q) xq[q] = PetscRealPart(Xq[q]); 831 PetscCall(PetscMalloc1(Nq, &wq)); 832 for (q = 0; q < Nq; ++q) wq[q] = 1.0; 833 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad)); 834 PetscCall(PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq)); 835 PetscCall(PetscFESetQuadrature(tr->coordFE[ct], quad)); 836 837 PetscCall(PetscFEGetDualSpace(tr->coordFE[ct], &dsp)); 838 PetscCall(PetscDualSpaceGetDM(dsp, &K)); 839 PetscCall(PetscFEGeomCreate(quad, 1, cdim, PETSC_FEGEOM_BASIC, &tr->refGeom[ct])); 840 cg = tr->refGeom[ct]; 841 PetscCall(DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ)); 842 PetscCall(PetscQuadratureDestroy(&quad)); 843 } 844 } 845 *fe = tr->coordFE[ct]; 846 PetscFunctionReturn(PETSC_SUCCESS); 847 } 848 849 PetscErrorCode DMPlexTransformSetDimensions_Internal(DMPlexTransform tr, DM dm, DM tdm) 850 { 851 PetscInt dim, cdim; 852 853 PetscFunctionBegin; 854 PetscCall(DMGetDimension(dm, &dim)); 855 PetscCall(DMSetDimension(tdm, dim)); 856 PetscCall(DMGetCoordinateDim(dm, &cdim)); 857 PetscCall(DMSetCoordinateDim(tdm, cdim)); 858 PetscFunctionReturn(PETSC_SUCCESS); 859 } 860 861 /*@ 862 DMPlexTransformSetDimensions - Set the dimensions for the transformed `DM` 863 864 Input Parameters: 865 + tr - The `DMPlexTransform` object 866 - dm - The original `DM` 867 868 Output Parameter: 869 . tdm - The transformed `DM` 870 871 Level: advanced 872 873 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` 874 @*/ 875 PetscErrorCode DMPlexTransformSetDimensions(DMPlexTransform tr, DM dm, DM tdm) 876 { 877 PetscFunctionBegin; 878 PetscUseTypeMethod(tr, setdimensions, dm, tdm); 879 PetscFunctionReturn(PETSC_SUCCESS); 880 } 881 882 PetscErrorCode DMPlexTransformGetChart(DMPlexTransform tr, PetscInt *pStart, PetscInt *pEnd) 883 { 884 PetscFunctionBegin; 885 if (pStart) *pStart = 0; 886 if (pEnd) *pEnd = tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]]; 887 PetscFunctionReturn(PETSC_SUCCESS); 888 } 889 890 PetscErrorCode DMPlexTransformGetCellType(DMPlexTransform tr, PetscInt cell, DMPolytopeType *celltype) 891 { 892 PetscInt ctNew; 893 894 PetscFunctionBegin; 895 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 896 PetscAssertPointer(celltype, 3); 897 /* TODO Can do bisection since everything is sorted */ 898 for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) { 899 PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]]; 900 901 if (cell >= ctSN && cell < ctEN) break; 902 } 903 PetscCheck(ctNew < DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " cannot be located in the transformed mesh", cell); 904 *celltype = (DMPolytopeType)ctNew; 905 PetscFunctionReturn(PETSC_SUCCESS); 906 } 907 908 PetscErrorCode DMPlexTransformGetCellTypeStratum(DMPlexTransform tr, DMPolytopeType celltype, PetscInt *start, PetscInt *end) 909 { 910 PetscFunctionBegin; 911 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 912 if (start) *start = tr->ctStartNew[celltype]; 913 if (end) *end = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[celltype] + 1]]; 914 PetscFunctionReturn(PETSC_SUCCESS); 915 } 916 917 PetscErrorCode DMPlexTransformGetDepth(DMPlexTransform tr, PetscInt *depth) 918 { 919 PetscFunctionBegin; 920 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 921 *depth = tr->depth; 922 PetscFunctionReturn(PETSC_SUCCESS); 923 } 924 925 PetscErrorCode DMPlexTransformGetDepthStratum(DMPlexTransform tr, PetscInt depth, PetscInt *start, PetscInt *end) 926 { 927 PetscFunctionBegin; 928 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 929 if (start) *start = tr->depthStart[depth]; 930 if (end) *end = tr->depthEnd[depth]; 931 PetscFunctionReturn(PETSC_SUCCESS); 932 } 933 934 /*@ 935 DMPlexTransformGetMatchStrata - Get the flag which determines what points get added to the transformed labels 936 937 Not Collective 938 939 Input Parameter: 940 . tr - The `DMPlexTransform` 941 942 Output Parameter: 943 . match - If `PETSC_TRUE`, only add produced points at the same stratum as the original point to new labels 944 945 Level: intermediate 946 947 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformSetMatchStrata()`, `DMPlexGetPointDepth()` 948 @*/ 949 PetscErrorCode DMPlexTransformGetMatchStrata(DMPlexTransform tr, PetscBool *match) 950 { 951 PetscFunctionBegin; 952 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 953 PetscAssertPointer(match, 2); 954 *match = tr->labelMatchStrata; 955 PetscFunctionReturn(PETSC_SUCCESS); 956 } 957 958 /*@ 959 DMPlexTransformSetMatchStrata - Set the flag which determines what points get added to the transformed labels 960 961 Not Collective 962 963 Input Parameters: 964 + tr - The `DMPlexTransform` 965 - match - If `PETSC_TRUE`, only add produced points at the same stratum as the original point to new labels 966 967 Level: intermediate 968 969 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformGetMatchStrata()`, `DMPlexGetPointDepth()` 970 @*/ 971 PetscErrorCode DMPlexTransformSetMatchStrata(DMPlexTransform tr, PetscBool match) 972 { 973 PetscFunctionBegin; 974 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 975 tr->labelMatchStrata = match; 976 PetscFunctionReturn(PETSC_SUCCESS); 977 } 978 979 /*@ 980 DMPlexTransformGetTargetPoint - Get the number of a point in the transformed mesh based on information from the original mesh. 981 982 Not Collective 983 984 Input Parameters: 985 + tr - The `DMPlexTransform` 986 . ct - The type of the original point which produces the new point 987 . ctNew - The type of the new point 988 . p - The original point which produces the new point 989 - r - The replica number of the new point, meaning it is the rth point of type `ctNew` produced from `p` 990 991 Output Parameter: 992 . pNew - The new point number 993 994 Level: developer 995 996 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSourcePoint()`, `DMPlexTransformCellTransform()` 997 @*/ 998 PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew) 999 { 1000 DMPolytopeType *rct; 1001 PetscInt *rsize, *cone, *ornt; 1002 PetscInt rt, Nct, n, off, rp; 1003 DMLabel trType = tr->trType; 1004 PetscInt ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ct] + 1]]; 1005 PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]]; 1006 PetscInt newp = ctSN, cind; 1007 1008 PetscFunctionBeginHot; 1009 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); 1010 PetscCall(DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt)); 1011 if (trType) { 1012 PetscCall(DMLabelGetValueIndex(trType, rt, &cind)); 1013 PetscCall(DMLabelGetStratumPointIndex(trType, rt, p, &rp)); 1014 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); 1015 } else { 1016 cind = ct; 1017 rp = p - ctS; 1018 } 1019 off = tr->offset[cind * DM_NUM_POLYTOPES + ctNew]; 1020 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); 1021 newp += off; 1022 for (n = 0; n < Nct; ++n) { 1023 if (rct[n] == ctNew) { 1024 if (rsize[n] && r >= rsize[n]) 1025 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]); 1026 newp += rp * rsize[n] + r; 1027 break; 1028 } 1029 } 1030 1031 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); 1032 *pNew = newp; 1033 PetscFunctionReturn(PETSC_SUCCESS); 1034 } 1035 1036 /*@ 1037 DMPlexTransformGetSourcePoint - Get the number of a point in the original mesh based on information from the transformed mesh. 1038 1039 Not Collective 1040 1041 Input Parameters: 1042 + tr - The `DMPlexTransform` 1043 - pNew - The new point number 1044 1045 Output Parameters: 1046 + ct - The type of the original point which produces the new point 1047 . ctNew - The type of the new point 1048 . p - The original point which produces the new point 1049 - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p 1050 1051 Level: developer 1052 1053 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetTargetPoint()`, `DMPlexTransformCellTransform()` 1054 @*/ 1055 PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r) 1056 { 1057 DMLabel trType = tr->trType; 1058 DMPolytopeType *rct, ctN; 1059 PetscInt *rsize, *cone, *ornt; 1060 PetscInt rt = -1, rtTmp, Nct, n, rp = 0, rO = 0, pO; 1061 PetscInt offset = -1, ctS, ctE, ctO = 0, ctTmp, rtS; 1062 1063 PetscFunctionBegin; 1064 PetscCall(DMPlexTransformGetCellType(tr, pNew, &ctN)); 1065 if (trType) { 1066 DM dm; 1067 IS rtIS; 1068 const PetscInt *reftypes; 1069 PetscInt Nrt, r, rtStart; 1070 1071 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1072 PetscCall(DMLabelGetNumValues(trType, &Nrt)); 1073 PetscCall(DMLabelGetValueIS(trType, &rtIS)); 1074 PetscCall(ISGetIndices(rtIS, &reftypes)); 1075 for (r = 0; r < Nrt; ++r) { 1076 const PetscInt off = tr->offset[r * DM_NUM_POLYTOPES + ctN]; 1077 1078 if (tr->ctStartNew[ctN] + off > pNew) continue; 1079 /* Check that any of this refinement type exist */ 1080 /* TODO Actually keep track of the number produced here instead */ 1081 if (off > offset) { 1082 rt = reftypes[r]; 1083 offset = off; 1084 } 1085 } 1086 PetscCall(ISRestoreIndices(rtIS, &reftypes)); 1087 PetscCall(ISDestroy(&rtIS)); 1088 PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew); 1089 /* TODO Map refinement types to cell types */ 1090 PetscCall(DMLabelGetStratumBounds(trType, rt, &rtStart, NULL)); 1091 PetscCheck(rtStart >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Refinement type %" PetscInt_FMT " has no source points", rt); 1092 for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) { 1093 PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]]; 1094 1095 if ((rtStart >= ctS) && (rtStart < ctE)) break; 1096 } 1097 PetscCheck(ctO != DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine a cell type for refinement type %" PetscInt_FMT, rt); 1098 } else { 1099 for (ctTmp = 0; ctTmp < DM_NUM_POLYTOPES; ++ctTmp) { 1100 const PetscInt off = tr->offset[ctTmp * DM_NUM_POLYTOPES + ctN]; 1101 1102 if (tr->ctStartNew[ctN] + off > pNew) continue; 1103 if (tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctTmp] + 1]] <= tr->ctStart[ctTmp]) continue; 1104 /* TODO Actually keep track of the number produced here instead */ 1105 if (off > offset) { 1106 ctO = ctTmp; 1107 offset = off; 1108 } 1109 } 1110 rt = -1; 1111 PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew); 1112 } 1113 ctS = tr->ctStart[ctO]; 1114 ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]]; 1115 if (trType) { 1116 for (rtS = ctS; rtS < ctE; ++rtS) { 1117 PetscInt val; 1118 PetscCall(DMLabelGetValue(trType, rtS, &val)); 1119 if (val == rt) break; 1120 } 1121 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); 1122 } else rtS = ctS; 1123 PetscCall(DMPlexTransformCellTransform(tr, (DMPolytopeType)ctO, rtS, &rtTmp, &Nct, &rct, &rsize, &cone, &ornt)); 1124 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); 1125 for (n = 0; n < Nct; ++n) { 1126 if (rct[n] == ctN) { 1127 PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset, val, c; 1128 1129 if (trType) { 1130 for (c = ctS; c < ctE; ++c) { 1131 PetscCall(DMLabelGetValue(trType, c, &val)); 1132 if (val == rt) { 1133 if (tmp < rsize[n]) break; 1134 tmp -= rsize[n]; 1135 } 1136 } 1137 PetscCheck(c < ctE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Parent point for target point %" PetscInt_FMT " could be not found", pNew); 1138 rp = c - ctS; 1139 rO = tmp; 1140 } else { 1141 // This assumes that all points of type ctO transform the same way 1142 rp = tmp / rsize[n]; 1143 rO = tmp % rsize[n]; 1144 } 1145 break; 1146 } 1147 } 1148 PetscCheck(n != Nct, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number for target point %" PetscInt_FMT " could be not found", pNew); 1149 pO = rp + ctS; 1150 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); 1151 if (ct) *ct = (DMPolytopeType)ctO; 1152 if (ctNew) *ctNew = ctN; 1153 if (p) *p = pO; 1154 if (r) *r = rO; 1155 PetscFunctionReturn(PETSC_SUCCESS); 1156 } 1157 1158 /*@ 1159 DMPlexTransformCellTransform - Describes the transform of a given source cell into a set of other target cells. These produced cells become the new mesh. 1160 1161 Input Parameters: 1162 + tr - The `DMPlexTransform` object 1163 . source - The source cell type 1164 - p - The source point, which can also determine the refine type 1165 1166 Output Parameters: 1167 + rt - The refine type for this point 1168 . Nt - The number of types produced by this point 1169 . target - An array of length `Nt` giving the types produced 1170 . size - An array of length `Nt` giving the number of cells of each type produced 1171 . cone - An array of length `Nt`*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point 1172 - ornt - An array of length `Nt`*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point 1173 1174 Level: advanced 1175 1176 Notes: 1177 The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we 1178 need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical 1179 division (described in these outputs) of the cell in the original mesh. The point identifier is given by 1180 .vb 1181 the number of cones to be taken, or 0 for the current cell 1182 the cell cone point number at each level from which it is subdivided 1183 the replica number r of the subdivision. 1184 .ve 1185 The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is 1186 .vb 1187 Nt = 2 1188 target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT} 1189 size = {1, 2} 1190 cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0} 1191 ornt = { 0, 0, 0, 0} 1192 .ve 1193 1194 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` 1195 @*/ 1196 PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 1197 { 1198 PetscFunctionBegin; 1199 PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt); 1200 PetscFunctionReturn(PETSC_SUCCESS); 1201 } 1202 1203 PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 1204 { 1205 PetscFunctionBegin; 1206 *rnew = r; 1207 *onew = DMPolytopeTypeComposeOrientation(tct, o, so); 1208 PetscFunctionReturn(PETSC_SUCCESS); 1209 } 1210 1211 /* Returns the same thing */ 1212 PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 1213 { 1214 static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT}; 1215 static PetscInt vertexS[] = {1}; 1216 static PetscInt vertexC[] = {0}; 1217 static PetscInt vertexO[] = {0}; 1218 static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT}; 1219 static PetscInt edgeS[] = {1}; 1220 static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 1221 static PetscInt edgeO[] = {0, 0}; 1222 static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR}; 1223 static PetscInt tedgeS[] = {1}; 1224 static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 1225 static PetscInt tedgeO[] = {0, 0}; 1226 static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE}; 1227 static PetscInt triS[] = {1}; 1228 static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0}; 1229 static PetscInt triO[] = {0, 0, 0}; 1230 static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL}; 1231 static PetscInt quadS[] = {1}; 1232 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}; 1233 static PetscInt quadO[] = {0, 0, 0, 0}; 1234 static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR}; 1235 static PetscInt tquadS[] = {1}; 1236 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}; 1237 static PetscInt tquadO[] = {0, 0, 0, 0}; 1238 static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON}; 1239 static PetscInt tetS[] = {1}; 1240 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}; 1241 static PetscInt tetO[] = {0, 0, 0, 0}; 1242 static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON}; 1243 static PetscInt hexS[] = {1}; 1244 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}; 1245 static PetscInt hexO[] = {0, 0, 0, 0, 0, 0}; 1246 static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM}; 1247 static PetscInt tripS[] = {1}; 1248 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}; 1249 static PetscInt tripO[] = {0, 0, 0, 0, 0}; 1250 static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR}; 1251 static PetscInt ttripS[] = {1}; 1252 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}; 1253 static PetscInt ttripO[] = {0, 0, 0, 0, 0}; 1254 static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR}; 1255 static PetscInt tquadpS[] = {1}; 1256 static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, 1257 DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0}; 1258 static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0}; 1259 static DMPolytopeType pyrT[] = {DM_POLYTOPE_PYRAMID}; 1260 static PetscInt pyrS[] = {1}; 1261 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}; 1262 static PetscInt pyrO[] = {0, 0, 0, 0, 0}; 1263 1264 PetscFunctionBegin; 1265 if (rt) *rt = 0; 1266 switch (source) { 1267 case DM_POLYTOPE_POINT: 1268 *Nt = 1; 1269 *target = vertexT; 1270 *size = vertexS; 1271 *cone = vertexC; 1272 *ornt = vertexO; 1273 break; 1274 case DM_POLYTOPE_SEGMENT: 1275 *Nt = 1; 1276 *target = edgeT; 1277 *size = edgeS; 1278 *cone = edgeC; 1279 *ornt = edgeO; 1280 break; 1281 case DM_POLYTOPE_POINT_PRISM_TENSOR: 1282 *Nt = 1; 1283 *target = tedgeT; 1284 *size = tedgeS; 1285 *cone = tedgeC; 1286 *ornt = tedgeO; 1287 break; 1288 case DM_POLYTOPE_TRIANGLE: 1289 *Nt = 1; 1290 *target = triT; 1291 *size = triS; 1292 *cone = triC; 1293 *ornt = triO; 1294 break; 1295 case DM_POLYTOPE_QUADRILATERAL: 1296 *Nt = 1; 1297 *target = quadT; 1298 *size = quadS; 1299 *cone = quadC; 1300 *ornt = quadO; 1301 break; 1302 case DM_POLYTOPE_SEG_PRISM_TENSOR: 1303 *Nt = 1; 1304 *target = tquadT; 1305 *size = tquadS; 1306 *cone = tquadC; 1307 *ornt = tquadO; 1308 break; 1309 case DM_POLYTOPE_TETRAHEDRON: 1310 *Nt = 1; 1311 *target = tetT; 1312 *size = tetS; 1313 *cone = tetC; 1314 *ornt = tetO; 1315 break; 1316 case DM_POLYTOPE_HEXAHEDRON: 1317 *Nt = 1; 1318 *target = hexT; 1319 *size = hexS; 1320 *cone = hexC; 1321 *ornt = hexO; 1322 break; 1323 case DM_POLYTOPE_TRI_PRISM: 1324 *Nt = 1; 1325 *target = tripT; 1326 *size = tripS; 1327 *cone = tripC; 1328 *ornt = tripO; 1329 break; 1330 case DM_POLYTOPE_TRI_PRISM_TENSOR: 1331 *Nt = 1; 1332 *target = ttripT; 1333 *size = ttripS; 1334 *cone = ttripC; 1335 *ornt = ttripO; 1336 break; 1337 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 1338 *Nt = 1; 1339 *target = tquadpT; 1340 *size = tquadpS; 1341 *cone = tquadpC; 1342 *ornt = tquadpO; 1343 break; 1344 case DM_POLYTOPE_PYRAMID: 1345 *Nt = 1; 1346 *target = pyrT; 1347 *size = pyrS; 1348 *cone = pyrC; 1349 *ornt = pyrO; 1350 break; 1351 default: 1352 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 1353 } 1354 PetscFunctionReturn(PETSC_SUCCESS); 1355 } 1356 1357 /*@ 1358 DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point 1359 1360 Not Collective 1361 1362 Input Parameters: 1363 + tr - The `DMPlexTransform` 1364 . sct - The source point cell type, from whom the new cell is being produced 1365 . sp - The source point 1366 . so - The orientation of the source point in its enclosing parent 1367 . tct - The target point cell type 1368 . r - The replica number requested for the produced cell type 1369 - o - The orientation of the replica 1370 1371 Output Parameters: 1372 + rnew - The replica number, given the orientation of the parent 1373 - onew - The replica orientation, given the orientation of the parent 1374 1375 Level: advanced 1376 1377 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformCellTransform()`, `DMPlexTransformApply()` 1378 @*/ 1379 PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 1380 { 1381 PetscFunctionBeginHot; 1382 PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew); 1383 PetscFunctionReturn(PETSC_SUCCESS); 1384 } 1385 1386 static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm) 1387 { 1388 DM dm; 1389 PetscInt pStart, pEnd, p, pNew; 1390 1391 PetscFunctionBegin; 1392 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1393 /* Must create the celltype label here so that we do not automatically try to compute the types */ 1394 PetscCall(DMCreateLabel(rdm, "celltype")); 1395 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1396 for (p = pStart; p < pEnd; ++p) { 1397 DMPolytopeType ct; 1398 DMPolytopeType *rct; 1399 PetscInt *rsize, *rcone, *rornt; 1400 PetscInt Nct, n, r; 1401 1402 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1403 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1404 for (n = 0; n < Nct; ++n) { 1405 for (r = 0; r < rsize[n]; ++r) { 1406 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 1407 PetscCall(DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]))); 1408 PetscCall(DMPlexSetCellType(rdm, pNew, rct[n])); 1409 } 1410 } 1411 } 1412 /* Let the DM know we have set all the cell types */ 1413 { 1414 DMLabel ctLabel; 1415 DM_Plex *plex = (DM_Plex *)rdm->data; 1416 1417 PetscCall(DMPlexGetCellTypeLabel(rdm, &ctLabel)); 1418 PetscCall(PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState)); 1419 } 1420 PetscFunctionReturn(PETSC_SUCCESS); 1421 } 1422 1423 PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize) 1424 { 1425 DMPolytopeType ctNew; 1426 1427 PetscFunctionBegin; 1428 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1429 PetscAssertPointer(coneSize, 3); 1430 PetscCall(DMPlexTransformGetCellType(tr, q, &ctNew)); 1431 *coneSize = DMPolytopeTypeGetConeSize(ctNew); 1432 PetscFunctionReturn(PETSC_SUCCESS); 1433 } 1434 1435 /* The orientation o is for the interior of the cell p */ 1436 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[]) 1437 { 1438 DM dm; 1439 const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); 1440 const PetscInt *cone; 1441 PetscInt c, coff = *coneoff, ooff = *orntoff; 1442 1443 PetscFunctionBegin; 1444 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1445 PetscCall(DMPlexGetOrientedCone(dm, p, &cone, NULL)); 1446 for (c = 0; c < csizeNew; ++c) { 1447 PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */ 1448 PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */ 1449 PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */ 1450 DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */ 1451 const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */ 1452 PetscInt pr = -1; /* Replica number of pp that produces new cone point */ 1453 const DMPolytopeType ft = (DMPolytopeType)rcone[coff++]; /* Cell type for new cone point of pNew */ 1454 const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ 1455 PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */ 1456 PetscInt lc; 1457 1458 /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */ 1459 for (lc = 0; lc < fn; ++lc) { 1460 const PetscInt *parr = DMPolytopeTypeGetArrangement(pct, po); 1461 const PetscInt acp = rcone[coff++]; 1462 const PetscInt pcp = parr[acp * 2]; 1463 const PetscInt pco = parr[acp * 2 + 1]; 1464 const PetscInt *ppornt; 1465 1466 ppp = pp; 1467 pp = pcone[pcp]; 1468 PetscCall(DMPlexGetCellType(dm, pp, &pct)); 1469 // Restore the parent cone from the last iterate 1470 if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, ppp, &pcone, NULL)); 1471 PetscCall(DMPlexGetOrientedCone(dm, pp, &pcone, NULL)); 1472 PetscCall(DMPlexGetOrientedCone(dm, ppp, NULL, &ppornt)); 1473 po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco); 1474 PetscCall(DMPlexRestoreOrientedCone(dm, ppp, NULL, &ppornt)); 1475 } 1476 if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, pp, &pcone, NULL)); 1477 pr = rcone[coff++]; 1478 /* Orientation po of pp maps (pr, fo) -> (pr', fo') */ 1479 PetscCall(DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo)); 1480 PetscCall(DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c])); 1481 orntNew[c] = fo; 1482 } 1483 PetscCall(DMPlexRestoreOrientedCone(dm, p, &cone, NULL)); 1484 *coneoff = coff; 1485 *orntoff = ooff; 1486 PetscFunctionReturn(PETSC_SUCCESS); 1487 } 1488 1489 static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm) 1490 { 1491 DM dm; 1492 DMPolytopeType ct; 1493 PetscInt *coneNew, *orntNew; 1494 PetscInt maxConeSize = 0, pStart, pEnd, p, pNew; 1495 1496 PetscFunctionBegin; 1497 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1498 for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p)); 1499 PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew)); 1500 PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew)); 1501 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1502 for (p = pStart; p < pEnd; ++p) { 1503 PetscInt coff, ooff; 1504 DMPolytopeType *rct; 1505 PetscInt *rsize, *rcone, *rornt; 1506 PetscInt Nct, n, r; 1507 1508 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1509 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1510 for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) { 1511 const DMPolytopeType ctNew = rct[n]; 1512 1513 for (r = 0; r < rsize[n]; ++r) { 1514 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 1515 PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew)); 1516 PetscCall(DMPlexSetCone(rdm, pNew, coneNew)); 1517 PetscCall(DMPlexSetConeOrientation(rdm, pNew, orntNew)); 1518 } 1519 } 1520 } 1521 PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew)); 1522 PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew)); 1523 PetscCall(DMViewFromOptions(rdm, NULL, "-rdm_view")); 1524 PetscCall(DMPlexSymmetrize(rdm)); 1525 PetscCall(DMPlexStratify(rdm)); 1526 PetscFunctionReturn(PETSC_SUCCESS); 1527 } 1528 1529 PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[]) 1530 { 1531 DM dm; 1532 DMPolytopeType ct, qct; 1533 DMPolytopeType *rct; 1534 PetscInt *rsize, *rcone, *rornt, *qcone, *qornt; 1535 PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0; 1536 1537 PetscFunctionBegin; 1538 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1539 PetscAssertPointer(cone, 4); 1540 PetscAssertPointer(ornt, 5); 1541 for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p)); 1542 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1543 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone)); 1544 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt)); 1545 PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r)); 1546 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1547 for (n = 0; n < Nct; ++n) { 1548 const DMPolytopeType ctNew = rct[n]; 1549 const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); 1550 PetscInt Nr = rsize[n], fn, c; 1551 1552 if (ctNew == qct) Nr = r; 1553 for (nr = 0; nr < Nr; ++nr) { 1554 for (c = 0; c < csizeNew; ++c) { 1555 ++coff; /* Cell type of new cone point */ 1556 fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ 1557 coff += fn; 1558 ++coff; /* Replica number of new cone point */ 1559 ++ooff; /* Orientation of new cone point */ 1560 } 1561 } 1562 if (ctNew == qct) break; 1563 } 1564 PetscCall(DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt)); 1565 *cone = qcone; 1566 *ornt = qornt; 1567 PetscFunctionReturn(PETSC_SUCCESS); 1568 } 1569 1570 PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[]) 1571 { 1572 DM dm; 1573 DMPolytopeType ct, qct; 1574 DMPolytopeType *rct; 1575 PetscInt *rsize, *rcone, *rornt, *qcone, *qornt; 1576 PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0; 1577 1578 PetscFunctionBegin; 1579 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1580 if (cone) PetscAssertPointer(cone, 3); 1581 if (ornt) PetscAssertPointer(ornt, 4); 1582 for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p)); 1583 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1584 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone)); 1585 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt)); 1586 PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r)); 1587 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1588 for (n = 0; n < Nct; ++n) { 1589 const DMPolytopeType ctNew = rct[n]; 1590 const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); 1591 PetscInt Nr = rsize[n], fn, c; 1592 1593 if (ctNew == qct) Nr = r; 1594 for (nr = 0; nr < Nr; ++nr) { 1595 for (c = 0; c < csizeNew; ++c) { 1596 ++coff; /* Cell type of new cone point */ 1597 fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ 1598 coff += fn; 1599 ++coff; /* Replica number of new cone point */ 1600 ++ooff; /* Orientation of new cone point */ 1601 } 1602 } 1603 if (ctNew == qct) break; 1604 } 1605 PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt)); 1606 if (cone) *cone = qcone; 1607 else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qcone)); 1608 if (ornt) *ornt = qornt; 1609 else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qornt)); 1610 PetscFunctionReturn(PETSC_SUCCESS); 1611 } 1612 1613 PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[]) 1614 { 1615 DM dm; 1616 1617 PetscFunctionBegin; 1618 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1619 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1620 if (cone) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, cone)); 1621 if (ornt) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, ornt)); 1622 PetscFunctionReturn(PETSC_SUCCESS); 1623 } 1624 1625 static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr) 1626 { 1627 PetscInt ict; 1628 1629 PetscFunctionBegin; 1630 PetscCall(PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts)); 1631 for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) { 1632 const DMPolytopeType ct = (DMPolytopeType)ict; 1633 DMPlexTransform reftr; 1634 DM refdm, trdm; 1635 Vec coordinates; 1636 const PetscScalar *coords; 1637 DMPolytopeType *rct; 1638 PetscInt *rsize, *rcone, *rornt; 1639 PetscInt Nct, n, r, pNew = 0; 1640 PetscInt trdim, vStart, vEnd, Nc; 1641 const PetscInt debug = 0; 1642 const char *typeName; 1643 1644 /* Since points are 0-dimensional, coordinates make no sense */ 1645 if (DMPolytopeTypeGetDim(ct) <= 0 || ct == DM_POLYTOPE_UNKNOWN_CELL || ct == DM_POLYTOPE_UNKNOWN_FACE) continue; 1646 PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm)); 1647 PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &reftr)); 1648 PetscCall(DMPlexTransformSetDM(reftr, refdm)); 1649 PetscCall(DMPlexTransformGetType(tr, &typeName)); 1650 PetscCall(DMPlexTransformSetType(reftr, typeName)); 1651 PetscCall(DMPlexTransformSetUp(reftr)); 1652 PetscCall(DMPlexTransformApply(reftr, refdm, &trdm)); 1653 1654 PetscCall(DMGetDimension(trdm, &trdim)); 1655 PetscCall(DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd)); 1656 tr->trNv[ct] = vEnd - vStart; 1657 PetscCall(DMGetCoordinatesLocal(trdm, &coordinates)); 1658 PetscCall(VecGetLocalSize(coordinates, &Nc)); 1659 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); 1660 PetscCall(PetscCalloc1(Nc, &tr->trVerts[ct])); 1661 PetscCall(VecGetArrayRead(coordinates, &coords)); 1662 PetscCall(PetscArraycpy(tr->trVerts[ct], coords, Nc)); 1663 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 1664 1665 PetscCall(PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct])); 1666 PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1667 for (n = 0; n < Nct; ++n) { 1668 /* Since points are 0-dimensional, coordinates make no sense */ 1669 if (rct[n] == DM_POLYTOPE_POINT) continue; 1670 PetscCall(PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]])); 1671 for (r = 0; r < rsize[n]; ++r) { 1672 PetscInt *closure = NULL; 1673 PetscInt clSize, cl, Nv = 0; 1674 1675 PetscCall(PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r])); 1676 PetscCall(DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew)); 1677 PetscCall(DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure)); 1678 for (cl = 0; cl < clSize * 2; cl += 2) { 1679 const PetscInt sv = closure[cl]; 1680 1681 if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart; 1682 } 1683 PetscCall(DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure)); 1684 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]); 1685 } 1686 } 1687 if (debug) { 1688 DMPolytopeType *rct; 1689 PetscInt *rsize, *rcone, *rornt; 1690 PetscInt v, dE = trdim, d, off = 0; 1691 1692 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct])); 1693 for (v = 0; v < tr->trNv[ct]; ++v) { 1694 PetscCall(PetscPrintf(PETSC_COMM_SELF, " ")); 1695 for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++]))); 1696 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 1697 } 1698 1699 PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1700 for (n = 0; n < Nct; ++n) { 1701 if (rct[n] == DM_POLYTOPE_POINT) continue; 1702 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct])); 1703 for (r = 0; r < rsize[n]; ++r) { 1704 PetscCall(PetscPrintf(PETSC_COMM_SELF, " ")); 1705 for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v])); 1706 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 1707 } 1708 } 1709 } 1710 PetscCall(DMDestroy(&refdm)); 1711 PetscCall(DMDestroy(&trdm)); 1712 PetscCall(DMPlexTransformDestroy(&reftr)); 1713 } 1714 PetscFunctionReturn(PETSC_SUCCESS); 1715 } 1716 1717 /*@C 1718 DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type 1719 1720 Input Parameters: 1721 + tr - The `DMPlexTransform` object 1722 - ct - The cell type 1723 1724 Output Parameters: 1725 + Nv - The number of transformed vertices in the closure of the reference cell of given type 1726 - trVerts - The coordinates of these vertices in the reference cell 1727 1728 Level: developer 1729 1730 .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSubcellVertices()` 1731 @*/ 1732 PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[]) 1733 { 1734 PetscFunctionBegin; 1735 if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr)); 1736 if (Nv) *Nv = tr->trNv[ct]; 1737 if (trVerts) *trVerts = tr->trVerts[ct]; 1738 PetscFunctionReturn(PETSC_SUCCESS); 1739 } 1740 1741 /*@C 1742 DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type 1743 1744 Input Parameters: 1745 + tr - The `DMPlexTransform` object 1746 . ct - The cell type 1747 . rct - The subcell type 1748 - r - The subcell index 1749 1750 Output Parameter: 1751 . subVerts - The indices of these vertices in the set of vertices returned by `DMPlexTransformGetCellVertices()` 1752 1753 Level: developer 1754 1755 .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetCellVertices()` 1756 @*/ 1757 PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[]) 1758 { 1759 PetscFunctionBegin; 1760 if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr)); 1761 PetscCheck(tr->trSubVerts[ct][rct], PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]); 1762 if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r]; 1763 PetscFunctionReturn(PETSC_SUCCESS); 1764 } 1765 1766 /* Computes new vertex as the barycenter, or centroid */ 1767 PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) 1768 { 1769 PetscInt v, d; 1770 1771 PetscFunctionBeginHot; 1772 PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]); 1773 for (d = 0; d < dE; ++d) out[d] = 0.0; 1774 for (v = 0; v < Nv; ++v) 1775 for (d = 0; d < dE; ++d) out[d] += in[v * dE + d]; 1776 for (d = 0; d < dE; ++d) out[d] /= Nv; 1777 PetscFunctionReturn(PETSC_SUCCESS); 1778 } 1779 1780 /*@ 1781 DMPlexTransformMapCoordinates - Calculate new coordinates for produced points 1782 1783 Not collective 1784 1785 Input Parameters: 1786 + tr - The `DMPlexTransform` 1787 . pct - The cell type of the parent, from whom the new cell is being produced 1788 . ct - The type being produced 1789 . p - The original point 1790 . r - The replica number requested for the produced cell type 1791 . Nv - Number of vertices in the closure of the parent cell 1792 . dE - Spatial dimension 1793 - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell 1794 1795 Output Parameter: 1796 . out - The coordinates of the new vertices 1797 1798 Level: intermediate 1799 1800 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()` 1801 @*/ 1802 PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) 1803 { 1804 PetscFunctionBeginHot; 1805 if (Nv) PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out); 1806 PetscFunctionReturn(PETSC_SUCCESS); 1807 } 1808 1809 /* 1810 DMPlexTransformLabelProducedPoint_Private - Label a produced point based on its parent label 1811 1812 Not Collective 1813 1814 Input Parameters: 1815 + tr - The `DMPlexTransform` 1816 . label - The label in the transformed mesh 1817 . pp - The parent point in the original mesh 1818 . pct - The cell type of the parent point 1819 . p - The point in the transformed mesh 1820 . ct - The cell type of the point 1821 . r - The replica number of the point 1822 - val - The label value of the parent point 1823 1824 Level: developer 1825 1826 .seealso: `DMPlexTransformCreateLabels()`, `RefineLabel_Internal()` 1827 */ 1828 static PetscErrorCode DMPlexTransformLabelProducedPoint_Private(DMPlexTransform tr, DMLabel label, PetscInt pp, DMPolytopeType pct, PetscInt p, DMPolytopeType ct, PetscInt r, PetscInt val) 1829 { 1830 PetscFunctionBeginHot; 1831 if (tr->labelMatchStrata && pct != ct) PetscFunctionReturn(PETSC_SUCCESS); 1832 PetscCall(DMLabelSetValue(label, p, val + tr->labelReplicaInc * r)); 1833 PetscFunctionReturn(PETSC_SUCCESS); 1834 } 1835 1836 static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew) 1837 { 1838 DM dm; 1839 IS valueIS; 1840 const PetscInt *values; 1841 PetscInt defVal, Nv, val; 1842 1843 PetscFunctionBegin; 1844 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1845 PetscCall(DMLabelGetDefaultValue(label, &defVal)); 1846 PetscCall(DMLabelSetDefaultValue(labelNew, defVal)); 1847 PetscCall(DMLabelGetValueIS(label, &valueIS)); 1848 PetscCall(ISGetLocalSize(valueIS, &Nv)); 1849 PetscCall(ISGetIndices(valueIS, &values)); 1850 for (val = 0; val < Nv; ++val) { 1851 IS pointIS; 1852 const PetscInt *points; 1853 PetscInt numPoints, p; 1854 1855 /* Ensure refined label is created with same number of strata as 1856 * original (even if no entries here). */ 1857 PetscCall(DMLabelAddStratum(labelNew, values[val])); 1858 PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS)); 1859 PetscCall(ISGetLocalSize(pointIS, &numPoints)); 1860 PetscCall(ISGetIndices(pointIS, &points)); 1861 for (p = 0; p < numPoints; ++p) { 1862 const PetscInt point = points[p]; 1863 DMPolytopeType ct; 1864 DMPolytopeType *rct; 1865 PetscInt *rsize, *rcone, *rornt; 1866 PetscInt Nct, n, r, pNew = 0; 1867 1868 PetscCall(DMPlexGetCellType(dm, point, &ct)); 1869 PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1870 for (n = 0; n < Nct; ++n) { 1871 for (r = 0; r < rsize[n]; ++r) { 1872 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew)); 1873 PetscCall(DMPlexTransformLabelProducedPoint_Private(tr, labelNew, point, ct, pNew, rct[n], r, values[val])); 1874 } 1875 } 1876 } 1877 PetscCall(ISRestoreIndices(pointIS, &points)); 1878 PetscCall(ISDestroy(&pointIS)); 1879 } 1880 PetscCall(ISRestoreIndices(valueIS, &values)); 1881 PetscCall(ISDestroy(&valueIS)); 1882 PetscFunctionReturn(PETSC_SUCCESS); 1883 } 1884 1885 static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm) 1886 { 1887 DM dm; 1888 PetscInt numLabels, l; 1889 1890 PetscFunctionBegin; 1891 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1892 PetscCall(DMGetNumLabels(dm, &numLabels)); 1893 for (l = 0; l < numLabels; ++l) { 1894 DMLabel label, labelNew; 1895 const char *lname; 1896 PetscBool isDepth, isCellType; 1897 1898 PetscCall(DMGetLabelName(dm, l, &lname)); 1899 PetscCall(PetscStrcmp(lname, "depth", &isDepth)); 1900 if (isDepth) continue; 1901 PetscCall(PetscStrcmp(lname, "celltype", &isCellType)); 1902 if (isCellType) continue; 1903 PetscCall(DMCreateLabel(rdm, lname)); 1904 PetscCall(DMGetLabel(dm, lname, &label)); 1905 PetscCall(DMGetLabel(rdm, lname, &labelNew)); 1906 PetscCall(RefineLabel_Internal(tr, label, labelNew)); 1907 } 1908 PetscFunctionReturn(PETSC_SUCCESS); 1909 } 1910 1911 /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */ 1912 PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm) 1913 { 1914 DM dm; 1915 PetscInt Nf, f, Nds, s; 1916 1917 PetscFunctionBegin; 1918 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1919 PetscCall(DMGetNumFields(dm, &Nf)); 1920 for (f = 0; f < Nf; ++f) { 1921 DMLabel label, labelNew; 1922 PetscObject obj; 1923 const char *lname; 1924 1925 PetscCall(DMGetField(rdm, f, &label, &obj)); 1926 if (!label) continue; 1927 PetscCall(PetscObjectGetName((PetscObject)label, &lname)); 1928 PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew)); 1929 PetscCall(RefineLabel_Internal(tr, label, labelNew)); 1930 PetscCall(DMSetField_Internal(rdm, f, labelNew, obj)); 1931 PetscCall(DMLabelDestroy(&labelNew)); 1932 } 1933 PetscCall(DMGetNumDS(dm, &Nds)); 1934 for (s = 0; s < Nds; ++s) { 1935 DMLabel label, labelNew; 1936 const char *lname; 1937 1938 PetscCall(DMGetRegionNumDS(rdm, s, &label, NULL, NULL, NULL)); 1939 if (!label) continue; 1940 PetscCall(PetscObjectGetName((PetscObject)label, &lname)); 1941 PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew)); 1942 PetscCall(RefineLabel_Internal(tr, label, labelNew)); 1943 PetscCall(DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL, NULL)); 1944 PetscCall(DMLabelDestroy(&labelNew)); 1945 } 1946 PetscFunctionReturn(PETSC_SUCCESS); 1947 } 1948 1949 static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm) 1950 { 1951 DM dm; 1952 PetscSF sf, sfNew; 1953 PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m; 1954 const PetscInt *localPoints; 1955 const PetscSFNode *remotePoints; 1956 PetscInt *localPointsNew; 1957 PetscSFNode *remotePointsNew; 1958 PetscInt pStartNew, pEndNew, pNew; 1959 /* Brute force algorithm */ 1960 PetscSF rsf; 1961 PetscSection s; 1962 const PetscInt *rootdegree; 1963 PetscInt *rootPointsNew, *remoteOffsets; 1964 PetscInt numPointsNew, pStart, pEnd, p; 1965 1966 PetscFunctionBegin; 1967 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1968 PetscCall(DMPlexGetChart(rdm, &pStartNew, &pEndNew)); 1969 PetscCall(DMGetPointSF(dm, &sf)); 1970 PetscCall(DMGetPointSF(rdm, &sfNew)); 1971 /* Calculate size of new SF */ 1972 PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints)); 1973 if (numRoots < 0) PetscFunctionReturn(PETSC_SUCCESS); 1974 for (l = 0; l < numLeaves; ++l) { 1975 const PetscInt p = localPoints[l]; 1976 DMPolytopeType ct; 1977 DMPolytopeType *rct; 1978 PetscInt *rsize, *rcone, *rornt; 1979 PetscInt Nct, n; 1980 1981 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1982 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1983 for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n]; 1984 } 1985 /* Send new root point numbers 1986 It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication 1987 */ 1988 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1989 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s)); 1990 PetscCall(PetscSectionSetChart(s, pStart, pEnd)); 1991 for (p = pStart; p < pEnd; ++p) { 1992 DMPolytopeType ct; 1993 DMPolytopeType *rct; 1994 PetscInt *rsize, *rcone, *rornt; 1995 PetscInt Nct, n; 1996 1997 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1998 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1999 for (n = 0; n < Nct; ++n) PetscCall(PetscSectionAddDof(s, p, rsize[n])); 2000 } 2001 PetscCall(PetscSectionSetUp(s)); 2002 PetscCall(PetscSectionGetStorageSize(s, &numPointsNew)); 2003 PetscCall(PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets)); 2004 PetscCall(PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf)); 2005 PetscCall(PetscFree(remoteOffsets)); 2006 PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree)); 2007 PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree)); 2008 PetscCall(PetscMalloc1(numPointsNew, &rootPointsNew)); 2009 for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1; 2010 for (p = pStart; p < pEnd; ++p) { 2011 DMPolytopeType ct; 2012 DMPolytopeType *rct; 2013 PetscInt *rsize, *rcone, *rornt; 2014 PetscInt Nct, n, r, off; 2015 2016 if (!rootdegree[p - pStart]) continue; 2017 PetscCall(PetscSectionGetOffset(s, p, &off)); 2018 PetscCall(DMPlexGetCellType(dm, p, &ct)); 2019 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 2020 for (n = 0, m = 0; n < Nct; ++n) { 2021 for (r = 0; r < rsize[n]; ++r, ++m) { 2022 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 2023 rootPointsNew[off + m] = pNew; 2024 } 2025 } 2026 } 2027 PetscCall(PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE)); 2028 PetscCall(PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE)); 2029 PetscCall(PetscSFDestroy(&rsf)); 2030 PetscCall(PetscMalloc1(numLeavesNew, &localPointsNew)); 2031 PetscCall(PetscMalloc1(numLeavesNew, &remotePointsNew)); 2032 for (l = 0, m = 0; l < numLeaves; ++l) { 2033 const PetscInt p = localPoints[l]; 2034 DMPolytopeType ct; 2035 DMPolytopeType *rct; 2036 PetscInt *rsize, *rcone, *rornt; 2037 PetscInt Nct, n, r, q, off; 2038 2039 PetscCall(PetscSectionGetOffset(s, p, &off)); 2040 PetscCall(DMPlexGetCellType(dm, p, &ct)); 2041 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 2042 for (n = 0, q = 0; n < Nct; ++n) { 2043 for (r = 0; r < rsize[n]; ++r, ++m, ++q) { 2044 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 2045 localPointsNew[m] = pNew; 2046 remotePointsNew[m].index = rootPointsNew[off + q]; 2047 remotePointsNew[m].rank = remotePoints[l].rank; 2048 } 2049 } 2050 } 2051 PetscCall(PetscSectionDestroy(&s)); 2052 PetscCall(PetscFree(rootPointsNew)); 2053 /* SF needs sorted leaves to correctly calculate Gather */ 2054 { 2055 PetscSFNode *rp, *rtmp; 2056 PetscInt *lp, *idx, *ltmp, i; 2057 2058 PetscCall(PetscMalloc1(numLeavesNew, &idx)); 2059 PetscCall(PetscMalloc1(numLeavesNew, &lp)); 2060 PetscCall(PetscMalloc1(numLeavesNew, &rp)); 2061 for (i = 0; i < numLeavesNew; ++i) { 2062 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); 2063 idx[i] = i; 2064 } 2065 PetscCall(PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx)); 2066 for (i = 0; i < numLeavesNew; ++i) { 2067 lp[i] = localPointsNew[idx[i]]; 2068 rp[i] = remotePointsNew[idx[i]]; 2069 } 2070 ltmp = localPointsNew; 2071 localPointsNew = lp; 2072 rtmp = remotePointsNew; 2073 remotePointsNew = rp; 2074 PetscCall(PetscFree(idx)); 2075 PetscCall(PetscFree(ltmp)); 2076 PetscCall(PetscFree(rtmp)); 2077 } 2078 PetscCall(PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER)); 2079 PetscFunctionReturn(PETSC_SUCCESS); 2080 } 2081 2082 /* 2083 DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of `DMPolytopeType` `ct` with localized coordinates `x`, generate localized coordinates `xr` for subcell `r` of type `rct`. 2084 2085 Not Collective 2086 2087 Input Parameters: 2088 + tr - The `DMPlexTransform` 2089 . ct - The type of the parent cell 2090 . rct - The type of the produced cell 2091 . r - The index of the produced cell 2092 - x - The localized coordinates for the parent cell 2093 2094 Output Parameter: 2095 . xr - The localized coordinates for the produced cell 2096 2097 Level: developer 2098 2099 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexCellRefinerSetCoordinates()` 2100 */ 2101 static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[]) 2102 { 2103 PetscFE fe = NULL; 2104 PetscInt cdim, v, *subcellV; 2105 2106 PetscFunctionBegin; 2107 PetscCall(DMPlexTransformGetCoordinateFE(tr, ct, &fe)); 2108 PetscCall(DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV)); 2109 PetscCall(PetscFEGetNumComponents(fe, &cdim)); 2110 for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscCall(PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim])); 2111 PetscFunctionReturn(PETSC_SUCCESS); 2112 } 2113 2114 static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm) 2115 { 2116 DM dm, cdm, cdmCell, cdmNew, cdmCellNew; 2117 PetscSection coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew; 2118 Vec coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew; 2119 const PetscScalar *coords; 2120 PetscScalar *coordsNew; 2121 const PetscReal *maxCell, *Lstart, *L; 2122 PetscBool localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE; 2123 PetscInt dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p; 2124 2125 PetscFunctionBegin; 2126 // Need to clear the DMField for coordinates 2127 PetscCall(DMSetCoordinateField(rdm, NULL)); 2128 PetscCall(DMPlexTransformGetDM(tr, &dm)); 2129 PetscCall(DMGetCoordinateDM(dm, &cdm)); 2130 PetscCall(DMGetCellCoordinateDM(dm, &cdmCell)); 2131 PetscCall(DMGetCoordinatesLocalized(dm, &localized)); 2132 PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L)); 2133 if (localized) { 2134 /* Localize coordinates of new vertices */ 2135 localizeVertices = PETSC_TRUE; 2136 /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */ 2137 if (!maxCell) localizeCells = PETSC_TRUE; 2138 } 2139 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 2140 PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &dEo)); 2141 if (maxCell) { 2142 PetscReal maxCellNew[3]; 2143 2144 for (d = 0; d < dEo; ++d) maxCellNew[d] = maxCell[d] / 2.0; 2145 PetscCall(DMSetPeriodicity(rdm, maxCellNew, Lstart, L)); 2146 } 2147 PetscCall(DMGetCoordinateDim(rdm, &dE)); 2148 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew)); 2149 PetscCall(PetscSectionSetNumFields(coordSectionNew, 1)); 2150 PetscCall(PetscSectionSetFieldComponents(coordSectionNew, 0, dE)); 2151 PetscCall(DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew)); 2152 PetscCall(PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew)); 2153 /* Localization should be inherited */ 2154 /* Stefano calculates parent cells for each new cell for localization */ 2155 /* Localized cells need coordinates of closure */ 2156 for (v = vStartNew; v < vEndNew; ++v) { 2157 PetscCall(PetscSectionSetDof(coordSectionNew, v, dE)); 2158 PetscCall(PetscSectionSetFieldDof(coordSectionNew, v, 0, dE)); 2159 } 2160 PetscCall(PetscSectionSetUp(coordSectionNew)); 2161 PetscCall(DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew)); 2162 2163 if (localizeCells) { 2164 PetscCall(DMGetCoordinateDM(rdm, &cdmNew)); 2165 PetscCall(DMClone(cdmNew, &cdmCellNew)); 2166 PetscCall(DMSetCellCoordinateDM(rdm, cdmCellNew)); 2167 PetscCall(DMDestroy(&cdmCellNew)); 2168 2169 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew)); 2170 PetscCall(PetscSectionSetNumFields(coordSectionCellNew, 1)); 2171 PetscCall(PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE)); 2172 PetscCall(DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew)); 2173 PetscCall(PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew)); 2174 2175 PetscCall(DMGetCellCoordinateSection(dm, &coordSectionCell)); 2176 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 2177 for (c = cStart; c < cEnd; ++c) { 2178 PetscInt dof; 2179 2180 PetscCall(PetscSectionGetDof(coordSectionCell, c, &dof)); 2181 if (dof) { 2182 DMPolytopeType ct; 2183 DMPolytopeType *rct; 2184 PetscInt *rsize, *rcone, *rornt; 2185 PetscInt dim, cNew, Nct, n, r; 2186 2187 PetscCall(DMPlexGetCellType(dm, c, &ct)); 2188 dim = DMPolytopeTypeGetDim(ct); 2189 PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 2190 /* This allows for different cell types */ 2191 for (n = 0; n < Nct; ++n) { 2192 if (dim != DMPolytopeTypeGetDim(rct[n])) continue; 2193 for (r = 0; r < rsize[n]; ++r) { 2194 PetscInt *closure = NULL; 2195 PetscInt clSize, cl, Nv = 0; 2196 2197 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew)); 2198 PetscCall(DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure)); 2199 for (cl = 0; cl < clSize * 2; cl += 2) { 2200 if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv; 2201 } 2202 PetscCall(DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure)); 2203 PetscCall(PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE)); 2204 PetscCall(PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE)); 2205 } 2206 } 2207 } 2208 } 2209 PetscCall(PetscSectionSetUp(coordSectionCellNew)); 2210 PetscCall(DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew)); 2211 } 2212 PetscCall(DMViewFromOptions(dm, NULL, "-coarse_dm_view")); 2213 { 2214 VecType vtype; 2215 PetscInt coordSizeNew, bs; 2216 const char *name; 2217 2218 PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal)); 2219 PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalNew)); 2220 PetscCall(PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew)); 2221 PetscCall(VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE)); 2222 PetscCall(PetscObjectGetName((PetscObject)coordsLocal, &name)); 2223 PetscCall(PetscObjectSetName((PetscObject)coordsLocalNew, name)); 2224 PetscCall(VecGetBlockSize(coordsLocal, &bs)); 2225 PetscCall(VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE)); 2226 PetscCall(VecGetType(coordsLocal, &vtype)); 2227 PetscCall(VecSetType(coordsLocalNew, vtype)); 2228 } 2229 PetscCall(VecGetArrayRead(coordsLocal, &coords)); 2230 PetscCall(VecGetArray(coordsLocalNew, &coordsNew)); 2231 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 2232 /* First set coordinates for vertices */ 2233 for (p = pStart; p < pEnd; ++p) { 2234 DMPolytopeType ct; 2235 DMPolytopeType *rct; 2236 PetscInt *rsize, *rcone, *rornt; 2237 PetscInt Nct, n, r; 2238 PetscBool hasVertex = PETSC_FALSE; 2239 2240 PetscCall(DMPlexGetCellType(dm, p, &ct)); 2241 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 2242 for (n = 0; n < Nct; ++n) { 2243 if (rct[n] == DM_POLYTOPE_POINT) { 2244 hasVertex = PETSC_TRUE; 2245 break; 2246 } 2247 } 2248 if (hasVertex) { 2249 const PetscScalar *icoords = NULL; 2250 const PetscScalar *array = NULL; 2251 PetscScalar *pcoords = NULL; 2252 PetscBool isDG; 2253 PetscInt Nc, Nv, v, d; 2254 2255 PetscCall(DMPlexGetCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords)); 2256 2257 icoords = pcoords; 2258 Nv = Nc / dEo; 2259 if (ct != DM_POLYTOPE_POINT) { 2260 if (localizeVertices && maxCell) { 2261 PetscScalar anchor[3]; 2262 2263 for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d]; 2264 for (v = 0; v < Nv; ++v) PetscCall(DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo])); 2265 } 2266 } 2267 for (n = 0; n < Nct; ++n) { 2268 if (rct[n] != DM_POLYTOPE_POINT) continue; 2269 for (r = 0; r < rsize[n]; ++r) { 2270 PetscScalar vcoords[3]; 2271 PetscInt vNew, off; 2272 2273 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew)); 2274 PetscCall(PetscSectionGetOffset(coordSectionNew, vNew, &off)); 2275 PetscCall(DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords)); 2276 PetscCall(DMSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off])); 2277 } 2278 } 2279 PetscCall(DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords)); 2280 } 2281 } 2282 PetscCall(VecRestoreArrayRead(coordsLocal, &coords)); 2283 PetscCall(VecRestoreArray(coordsLocalNew, &coordsNew)); 2284 PetscCall(DMSetCoordinatesLocal(rdm, coordsLocalNew)); 2285 PetscCall(VecDestroy(&coordsLocalNew)); 2286 PetscCall(PetscSectionDestroy(&coordSectionNew)); 2287 /* Then set coordinates for cells by localizing */ 2288 if (!localizeCells) PetscCall(DMLocalizeCoordinates(rdm)); 2289 else { 2290 VecType vtype; 2291 PetscInt coordSizeNew, bs; 2292 const char *name; 2293 2294 PetscCall(DMGetCellCoordinatesLocal(dm, &coordsLocalCell)); 2295 PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew)); 2296 PetscCall(PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew)); 2297 PetscCall(VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE)); 2298 PetscCall(PetscObjectGetName((PetscObject)coordsLocalCell, &name)); 2299 PetscCall(PetscObjectSetName((PetscObject)coordsLocalCellNew, name)); 2300 PetscCall(VecGetBlockSize(coordsLocalCell, &bs)); 2301 PetscCall(VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE)); 2302 PetscCall(VecGetType(coordsLocalCell, &vtype)); 2303 PetscCall(VecSetType(coordsLocalCellNew, vtype)); 2304 PetscCall(VecGetArrayRead(coordsLocalCell, &coords)); 2305 PetscCall(VecGetArray(coordsLocalCellNew, &coordsNew)); 2306 2307 for (p = pStart; p < pEnd; ++p) { 2308 DMPolytopeType ct; 2309 DMPolytopeType *rct; 2310 PetscInt *rsize, *rcone, *rornt; 2311 PetscInt dof = 0, Nct, n, r; 2312 2313 PetscCall(DMPlexGetCellType(dm, p, &ct)); 2314 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 2315 if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(coordSectionCell, p, &dof)); 2316 if (dof) { 2317 const PetscScalar *pcoords; 2318 2319 PetscCall(DMPlexPointLocalRead(cdmCell, p, coords, &pcoords)); 2320 for (n = 0; n < Nct; ++n) { 2321 const PetscInt Nr = rsize[n]; 2322 2323 if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue; 2324 for (r = 0; r < Nr; ++r) { 2325 PetscInt pNew, offNew; 2326 2327 /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that 2328 DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger 2329 cell to the ones it produces. */ 2330 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 2331 PetscCall(PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew)); 2332 PetscCall(DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew])); 2333 } 2334 } 2335 } 2336 } 2337 PetscCall(VecRestoreArrayRead(coordsLocalCell, &coords)); 2338 PetscCall(VecRestoreArray(coordsLocalCellNew, &coordsNew)); 2339 PetscCall(DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew)); 2340 PetscCall(VecDestroy(&coordsLocalCellNew)); 2341 PetscCall(PetscSectionDestroy(&coordSectionCellNew)); 2342 } 2343 PetscFunctionReturn(PETSC_SUCCESS); 2344 } 2345 2346 /*@ 2347 DMPlexTransformApply - Execute the transformation, producing another `DM` 2348 2349 Collective 2350 2351 Input Parameters: 2352 + tr - The `DMPlexTransform` object 2353 - dm - The original `DM` 2354 2355 Output Parameter: 2356 . tdm - The transformed `DM` 2357 2358 Level: intermediate 2359 2360 Options Database Keys: 2361 + -dm_plex_transform_label_match_strata - Only label points of the same stratum as the producing point 2362 . -dm_plex_transform_label_replica_inc <num> - Increment for the label value to be multiplied by the replica number 2363 - -dm_plex_transform_active <name> - Name for active mesh label 2364 2365 .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformCreate()`, `DMPlexTransformSetDM()` 2366 @*/ 2367 PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm) 2368 { 2369 DM rdm; 2370 DMPlexInterpolatedFlag interp; 2371 PetscInt pStart, pEnd; 2372 2373 PetscFunctionBegin; 2374 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 2375 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 2376 PetscAssertPointer(tdm, 3); 2377 PetscCall(PetscLogEventBegin(DMPLEX_Transform, tr, dm, 0, 0)); 2378 PetscCall(DMPlexTransformSetDM(tr, dm)); 2379 2380 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &rdm)); 2381 PetscCall(DMSetType(rdm, DMPLEX)); 2382 PetscCall(DMPlexTransformSetDimensions(tr, dm, rdm)); 2383 /* Calculate number of new points of each depth */ 2384 PetscCall(DMPlexIsInterpolatedCollective(dm, &interp)); 2385 PetscCheck(interp == DMPLEX_INTERPOLATED_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement"); 2386 /* Step 1: Set chart */ 2387 PetscCall(DMPlexTransformGetChart(tr, &pStart, &pEnd)); 2388 PetscCall(DMPlexSetChart(rdm, pStart, pEnd)); 2389 /* Step 2: Set cone/support sizes (automatically stratifies) */ 2390 PetscCall(DMPlexTransformSetConeSizes(tr, rdm)); 2391 /* Step 3: Setup refined DM */ 2392 PetscCall(DMSetUp(rdm)); 2393 /* Step 4: Set cones and supports (automatically symmetrizes) */ 2394 PetscCall(DMPlexTransformSetCones(tr, rdm)); 2395 /* Step 5: Create pointSF */ 2396 PetscCall(DMPlexTransformCreateSF(tr, rdm)); 2397 /* Step 6: Create labels */ 2398 PetscCall(DMPlexTransformCreateLabels(tr, rdm)); 2399 /* Step 7: Set coordinates */ 2400 PetscCall(DMPlexTransformSetCoordinates(tr, rdm)); 2401 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, rdm)); 2402 // If the original DM was configured from options, the transformed DM should be as well 2403 rdm->setfromoptionscalled = dm->setfromoptionscalled; 2404 PetscCall(PetscLogEventEnd(DMPLEX_Transform, tr, dm, 0, 0)); 2405 *tdm = rdm; 2406 PetscFunctionReturn(PETSC_SUCCESS); 2407 } 2408 2409 PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm) 2410 { 2411 DMPlexTransform tr; 2412 DM cdm, rcdm; 2413 const char *prefix; 2414 2415 PetscFunctionBegin; 2416 PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr)); 2417 PetscCall(PetscObjectSetName((PetscObject)tr, "Adapt Label Transform")); 2418 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix)); 2419 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix)); 2420 PetscCall(DMPlexTransformSetDM(tr, dm)); 2421 PetscCall(DMPlexTransformSetFromOptions(tr)); 2422 PetscCall(DMPlexTransformSetActive(tr, adaptLabel)); 2423 PetscCall(DMPlexTransformSetUp(tr)); 2424 PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view")); 2425 PetscCall(DMPlexTransformApply(tr, dm, rdm)); 2426 PetscCall(DMCopyDisc(dm, *rdm)); 2427 PetscCall(DMGetCoordinateDM(dm, &cdm)); 2428 PetscCall(DMGetCoordinateDM(*rdm, &rcdm)); 2429 PetscCall(DMCopyDisc(cdm, rcdm)); 2430 PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm)); 2431 PetscCall(DMCopyDisc(dm, *rdm)); 2432 PetscCall(DMPlexTransformDestroy(&tr)); 2433 ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation; 2434 PetscFunctionReturn(PETSC_SUCCESS); 2435 } 2436