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