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 PetscFunctionReturn(PETSC_SUCCESS); 1521 } 1522 1523 PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[]) 1524 { 1525 DM dm; 1526 DMPolytopeType ct, qct; 1527 DMPolytopeType *rct; 1528 PetscInt *rsize, *rcone, *rornt, *qcone, *qornt; 1529 PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0; 1530 1531 PetscFunctionBegin; 1532 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1533 PetscAssertPointer(cone, 4); 1534 PetscAssertPointer(ornt, 5); 1535 for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p)); 1536 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1537 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone)); 1538 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt)); 1539 PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r)); 1540 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1541 for (n = 0; n < Nct; ++n) { 1542 const DMPolytopeType ctNew = rct[n]; 1543 const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); 1544 PetscInt Nr = rsize[n], fn, c; 1545 1546 if (ctNew == qct) Nr = r; 1547 for (nr = 0; nr < Nr; ++nr) { 1548 for (c = 0; c < csizeNew; ++c) { 1549 ++coff; /* Cell type of new cone point */ 1550 fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ 1551 coff += fn; 1552 ++coff; /* Replica number of new cone point */ 1553 ++ooff; /* Orientation of new cone point */ 1554 } 1555 } 1556 if (ctNew == qct) break; 1557 } 1558 PetscCall(DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt)); 1559 *cone = qcone; 1560 *ornt = qornt; 1561 PetscFunctionReturn(PETSC_SUCCESS); 1562 } 1563 1564 PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[]) 1565 { 1566 DM dm; 1567 DMPolytopeType ct, qct; 1568 DMPolytopeType *rct; 1569 PetscInt *rsize, *rcone, *rornt, *qcone, *qornt; 1570 PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0; 1571 1572 PetscFunctionBegin; 1573 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1574 if (cone) PetscAssertPointer(cone, 3); 1575 if (ornt) PetscAssertPointer(ornt, 4); 1576 for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p)); 1577 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1578 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone)); 1579 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt)); 1580 PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r)); 1581 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1582 for (n = 0; n < Nct; ++n) { 1583 const DMPolytopeType ctNew = rct[n]; 1584 const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); 1585 PetscInt Nr = rsize[n], fn, c; 1586 1587 if (ctNew == qct) Nr = r; 1588 for (nr = 0; nr < Nr; ++nr) { 1589 for (c = 0; c < csizeNew; ++c) { 1590 ++coff; /* Cell type of new cone point */ 1591 fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ 1592 coff += fn; 1593 ++coff; /* Replica number of new cone point */ 1594 ++ooff; /* Orientation of new cone point */ 1595 } 1596 } 1597 if (ctNew == qct) break; 1598 } 1599 PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt)); 1600 if (cone) *cone = qcone; 1601 else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qcone)); 1602 if (ornt) *ornt = qornt; 1603 else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qornt)); 1604 PetscFunctionReturn(PETSC_SUCCESS); 1605 } 1606 1607 PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[]) 1608 { 1609 DM dm; 1610 1611 PetscFunctionBegin; 1612 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1613 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1614 if (cone) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, cone)); 1615 if (ornt) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, ornt)); 1616 PetscFunctionReturn(PETSC_SUCCESS); 1617 } 1618 1619 static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr) 1620 { 1621 PetscInt ict; 1622 1623 PetscFunctionBegin; 1624 PetscCall(PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts)); 1625 for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) { 1626 const DMPolytopeType ct = (DMPolytopeType)ict; 1627 DMPlexTransform reftr; 1628 DM refdm, trdm; 1629 Vec coordinates; 1630 const PetscScalar *coords; 1631 DMPolytopeType *rct; 1632 PetscInt *rsize, *rcone, *rornt; 1633 PetscInt Nct, n, r, pNew = 0; 1634 PetscInt trdim, vStart, vEnd, Nc; 1635 const PetscInt debug = 0; 1636 const char *typeName; 1637 1638 /* Since points are 0-dimensional, coordinates make no sense */ 1639 if (DMPolytopeTypeGetDim(ct) <= 0 || ct == DM_POLYTOPE_UNKNOWN_CELL || ct == DM_POLYTOPE_UNKNOWN_FACE) continue; 1640 PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm)); 1641 PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &reftr)); 1642 PetscCall(DMPlexTransformSetDM(reftr, refdm)); 1643 PetscCall(DMPlexTransformGetType(tr, &typeName)); 1644 PetscCall(DMPlexTransformSetType(reftr, typeName)); 1645 PetscCall(DMPlexTransformSetUp(reftr)); 1646 PetscCall(DMPlexTransformApply(reftr, refdm, &trdm)); 1647 1648 PetscCall(DMGetDimension(trdm, &trdim)); 1649 PetscCall(DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd)); 1650 tr->trNv[ct] = vEnd - vStart; 1651 PetscCall(DMGetCoordinatesLocal(trdm, &coordinates)); 1652 PetscCall(VecGetLocalSize(coordinates, &Nc)); 1653 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); 1654 PetscCall(PetscCalloc1(Nc, &tr->trVerts[ct])); 1655 PetscCall(VecGetArrayRead(coordinates, &coords)); 1656 PetscCall(PetscArraycpy(tr->trVerts[ct], coords, Nc)); 1657 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 1658 1659 PetscCall(PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct])); 1660 PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1661 for (n = 0; n < Nct; ++n) { 1662 /* Since points are 0-dimensional, coordinates make no sense */ 1663 if (rct[n] == DM_POLYTOPE_POINT) continue; 1664 PetscCall(PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]])); 1665 for (r = 0; r < rsize[n]; ++r) { 1666 PetscInt *closure = NULL; 1667 PetscInt clSize, cl, Nv = 0; 1668 1669 PetscCall(PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r])); 1670 PetscCall(DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew)); 1671 PetscCall(DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure)); 1672 for (cl = 0; cl < clSize * 2; cl += 2) { 1673 const PetscInt sv = closure[cl]; 1674 1675 if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart; 1676 } 1677 PetscCall(DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure)); 1678 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]); 1679 } 1680 } 1681 if (debug) { 1682 DMPolytopeType *rct; 1683 PetscInt *rsize, *rcone, *rornt; 1684 PetscInt v, dE = trdim, d, off = 0; 1685 1686 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct])); 1687 for (v = 0; v < tr->trNv[ct]; ++v) { 1688 PetscCall(PetscPrintf(PETSC_COMM_SELF, " ")); 1689 for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++]))); 1690 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 1691 } 1692 1693 PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1694 for (n = 0; n < Nct; ++n) { 1695 if (rct[n] == DM_POLYTOPE_POINT) continue; 1696 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct])); 1697 for (r = 0; r < rsize[n]; ++r) { 1698 PetscCall(PetscPrintf(PETSC_COMM_SELF, " ")); 1699 for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v])); 1700 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 1701 } 1702 } 1703 } 1704 PetscCall(DMDestroy(&refdm)); 1705 PetscCall(DMDestroy(&trdm)); 1706 PetscCall(DMPlexTransformDestroy(&reftr)); 1707 } 1708 PetscFunctionReturn(PETSC_SUCCESS); 1709 } 1710 1711 /*@C 1712 DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type 1713 1714 Input Parameters: 1715 + tr - The `DMPlexTransform` object 1716 - ct - The cell type 1717 1718 Output Parameters: 1719 + Nv - The number of transformed vertices in the closure of the reference cell of given type 1720 - trVerts - The coordinates of these vertices in the reference cell 1721 1722 Level: developer 1723 1724 .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSubcellVertices()` 1725 @*/ 1726 PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[]) 1727 { 1728 PetscFunctionBegin; 1729 if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr)); 1730 if (Nv) *Nv = tr->trNv[ct]; 1731 if (trVerts) *trVerts = tr->trVerts[ct]; 1732 PetscFunctionReturn(PETSC_SUCCESS); 1733 } 1734 1735 /*@C 1736 DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type 1737 1738 Input Parameters: 1739 + tr - The `DMPlexTransform` object 1740 . ct - The cell type 1741 . rct - The subcell type 1742 - r - The subcell index 1743 1744 Output Parameter: 1745 . subVerts - The indices of these vertices in the set of vertices returned by `DMPlexTransformGetCellVertices()` 1746 1747 Level: developer 1748 1749 .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetCellVertices()` 1750 @*/ 1751 PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[]) 1752 { 1753 PetscFunctionBegin; 1754 if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr)); 1755 PetscCheck(tr->trSubVerts[ct][rct], PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]); 1756 if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r]; 1757 PetscFunctionReturn(PETSC_SUCCESS); 1758 } 1759 1760 /* Computes new vertex as the barycenter, or centroid */ 1761 PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) 1762 { 1763 PetscInt v, d; 1764 1765 PetscFunctionBeginHot; 1766 PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]); 1767 for (d = 0; d < dE; ++d) out[d] = 0.0; 1768 for (v = 0; v < Nv; ++v) 1769 for (d = 0; d < dE; ++d) out[d] += in[v * dE + d]; 1770 for (d = 0; d < dE; ++d) out[d] /= Nv; 1771 PetscFunctionReturn(PETSC_SUCCESS); 1772 } 1773 1774 /*@ 1775 DMPlexTransformMapCoordinates - Calculate new coordinates for produced points 1776 1777 Not collective 1778 1779 Input Parameters: 1780 + tr - The `DMPlexTransform` 1781 . pct - The cell type of the parent, from whom the new cell is being produced 1782 . ct - The type being produced 1783 . p - The original point 1784 . r - The replica number requested for the produced cell type 1785 . Nv - Number of vertices in the closure of the parent cell 1786 . dE - Spatial dimension 1787 - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell 1788 1789 Output Parameter: 1790 . out - The coordinates of the new vertices 1791 1792 Level: intermediate 1793 1794 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()` 1795 @*/ 1796 PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) 1797 { 1798 PetscFunctionBeginHot; 1799 if (Nv) PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out); 1800 PetscFunctionReturn(PETSC_SUCCESS); 1801 } 1802 1803 /* 1804 DMPlexTransformLabelProducedPoint_Private - Label a produced point based on its parent label 1805 1806 Not Collective 1807 1808 Input Parameters: 1809 + tr - The `DMPlexTransform` 1810 . label - The label in the transformed mesh 1811 . pp - The parent point in the original mesh 1812 . pct - The cell type of the parent point 1813 . p - The point in the transformed mesh 1814 . ct - The cell type of the point 1815 . r - The replica number of the point 1816 - val - The label value of the parent point 1817 1818 Level: developer 1819 1820 .seealso: `DMPlexTransformCreateLabels()`, `RefineLabel_Internal()` 1821 */ 1822 static PetscErrorCode DMPlexTransformLabelProducedPoint_Private(DMPlexTransform tr, DMLabel label, PetscInt pp, DMPolytopeType pct, PetscInt p, DMPolytopeType ct, PetscInt r, PetscInt val) 1823 { 1824 PetscFunctionBeginHot; 1825 if (tr->labelMatchStrata && pct != ct) PetscFunctionReturn(PETSC_SUCCESS); 1826 PetscCall(DMLabelSetValue(label, p, val + tr->labelReplicaInc * r)); 1827 PetscFunctionReturn(PETSC_SUCCESS); 1828 } 1829 1830 static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew) 1831 { 1832 DM dm; 1833 IS valueIS; 1834 const PetscInt *values; 1835 PetscInt defVal, Nv, val; 1836 1837 PetscFunctionBegin; 1838 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1839 PetscCall(DMLabelGetDefaultValue(label, &defVal)); 1840 PetscCall(DMLabelSetDefaultValue(labelNew, defVal)); 1841 PetscCall(DMLabelGetValueIS(label, &valueIS)); 1842 PetscCall(ISGetLocalSize(valueIS, &Nv)); 1843 PetscCall(ISGetIndices(valueIS, &values)); 1844 for (val = 0; val < Nv; ++val) { 1845 IS pointIS; 1846 const PetscInt *points; 1847 PetscInt numPoints, p; 1848 1849 /* Ensure refined label is created with same number of strata as 1850 * original (even if no entries here). */ 1851 PetscCall(DMLabelAddStratum(labelNew, values[val])); 1852 PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS)); 1853 PetscCall(ISGetLocalSize(pointIS, &numPoints)); 1854 PetscCall(ISGetIndices(pointIS, &points)); 1855 for (p = 0; p < numPoints; ++p) { 1856 const PetscInt point = points[p]; 1857 DMPolytopeType ct; 1858 DMPolytopeType *rct; 1859 PetscInt *rsize, *rcone, *rornt; 1860 PetscInt Nct, n, r, pNew = 0; 1861 1862 PetscCall(DMPlexGetCellType(dm, point, &ct)); 1863 PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1864 for (n = 0; n < Nct; ++n) { 1865 for (r = 0; r < rsize[n]; ++r) { 1866 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew)); 1867 PetscCall(DMPlexTransformLabelProducedPoint_Private(tr, labelNew, point, ct, pNew, rct[n], r, values[val])); 1868 } 1869 } 1870 } 1871 PetscCall(ISRestoreIndices(pointIS, &points)); 1872 PetscCall(ISDestroy(&pointIS)); 1873 } 1874 PetscCall(ISRestoreIndices(valueIS, &values)); 1875 PetscCall(ISDestroy(&valueIS)); 1876 PetscFunctionReturn(PETSC_SUCCESS); 1877 } 1878 1879 static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm) 1880 { 1881 DM dm; 1882 PetscInt numLabels, l; 1883 1884 PetscFunctionBegin; 1885 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1886 PetscCall(DMGetNumLabels(dm, &numLabels)); 1887 for (l = 0; l < numLabels; ++l) { 1888 DMLabel label, labelNew; 1889 const char *lname; 1890 PetscBool isDepth, isCellType; 1891 1892 PetscCall(DMGetLabelName(dm, l, &lname)); 1893 PetscCall(PetscStrcmp(lname, "depth", &isDepth)); 1894 if (isDepth) continue; 1895 PetscCall(PetscStrcmp(lname, "celltype", &isCellType)); 1896 if (isCellType) continue; 1897 PetscCall(DMCreateLabel(rdm, lname)); 1898 PetscCall(DMGetLabel(dm, lname, &label)); 1899 PetscCall(DMGetLabel(rdm, lname, &labelNew)); 1900 PetscCall(RefineLabel_Internal(tr, label, labelNew)); 1901 } 1902 PetscFunctionReturn(PETSC_SUCCESS); 1903 } 1904 1905 /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */ 1906 PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm) 1907 { 1908 DM dm; 1909 PetscInt Nf, f, Nds, s; 1910 1911 PetscFunctionBegin; 1912 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1913 PetscCall(DMGetNumFields(dm, &Nf)); 1914 for (f = 0; f < Nf; ++f) { 1915 DMLabel label, labelNew; 1916 PetscObject obj; 1917 const char *lname; 1918 1919 PetscCall(DMGetField(rdm, f, &label, &obj)); 1920 if (!label) continue; 1921 PetscCall(PetscObjectGetName((PetscObject)label, &lname)); 1922 PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew)); 1923 PetscCall(RefineLabel_Internal(tr, label, labelNew)); 1924 PetscCall(DMSetField_Internal(rdm, f, labelNew, obj)); 1925 PetscCall(DMLabelDestroy(&labelNew)); 1926 } 1927 PetscCall(DMGetNumDS(dm, &Nds)); 1928 for (s = 0; s < Nds; ++s) { 1929 DMLabel label, labelNew; 1930 const char *lname; 1931 1932 PetscCall(DMGetRegionNumDS(rdm, s, &label, NULL, NULL, NULL)); 1933 if (!label) continue; 1934 PetscCall(PetscObjectGetName((PetscObject)label, &lname)); 1935 PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew)); 1936 PetscCall(RefineLabel_Internal(tr, label, labelNew)); 1937 PetscCall(DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL, NULL)); 1938 PetscCall(DMLabelDestroy(&labelNew)); 1939 } 1940 PetscFunctionReturn(PETSC_SUCCESS); 1941 } 1942 1943 static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm) 1944 { 1945 DM dm; 1946 PetscSF sf, sfNew; 1947 PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m; 1948 const PetscInt *localPoints; 1949 const PetscSFNode *remotePoints; 1950 PetscInt *localPointsNew; 1951 PetscSFNode *remotePointsNew; 1952 PetscInt pStartNew, pEndNew, pNew; 1953 /* Brute force algorithm */ 1954 PetscSF rsf; 1955 PetscSection s; 1956 const PetscInt *rootdegree; 1957 PetscInt *rootPointsNew, *remoteOffsets; 1958 PetscInt numPointsNew, pStart, pEnd, p; 1959 1960 PetscFunctionBegin; 1961 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1962 PetscCall(DMPlexGetChart(rdm, &pStartNew, &pEndNew)); 1963 PetscCall(DMGetPointSF(dm, &sf)); 1964 PetscCall(DMGetPointSF(rdm, &sfNew)); 1965 /* Calculate size of new SF */ 1966 PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints)); 1967 if (numRoots < 0) PetscFunctionReturn(PETSC_SUCCESS); 1968 for (l = 0; l < numLeaves; ++l) { 1969 const PetscInt p = localPoints[l]; 1970 DMPolytopeType ct; 1971 DMPolytopeType *rct; 1972 PetscInt *rsize, *rcone, *rornt; 1973 PetscInt Nct, n; 1974 1975 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1976 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1977 for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n]; 1978 } 1979 /* Send new root point numbers 1980 It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication 1981 */ 1982 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1983 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s)); 1984 PetscCall(PetscSectionSetChart(s, pStart, pEnd)); 1985 for (p = pStart; p < pEnd; ++p) { 1986 DMPolytopeType ct; 1987 DMPolytopeType *rct; 1988 PetscInt *rsize, *rcone, *rornt; 1989 PetscInt Nct, n; 1990 1991 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1992 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1993 for (n = 0; n < Nct; ++n) PetscCall(PetscSectionAddDof(s, p, rsize[n])); 1994 } 1995 PetscCall(PetscSectionSetUp(s)); 1996 PetscCall(PetscSectionGetStorageSize(s, &numPointsNew)); 1997 PetscCall(PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets)); 1998 PetscCall(PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf)); 1999 PetscCall(PetscFree(remoteOffsets)); 2000 PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree)); 2001 PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree)); 2002 PetscCall(PetscMalloc1(numPointsNew, &rootPointsNew)); 2003 for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1; 2004 for (p = pStart; p < pEnd; ++p) { 2005 DMPolytopeType ct; 2006 DMPolytopeType *rct; 2007 PetscInt *rsize, *rcone, *rornt; 2008 PetscInt Nct, n, r, off; 2009 2010 if (!rootdegree[p - pStart]) continue; 2011 PetscCall(PetscSectionGetOffset(s, p, &off)); 2012 PetscCall(DMPlexGetCellType(dm, p, &ct)); 2013 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 2014 for (n = 0, m = 0; n < Nct; ++n) { 2015 for (r = 0; r < rsize[n]; ++r, ++m) { 2016 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 2017 rootPointsNew[off + m] = pNew; 2018 } 2019 } 2020 } 2021 PetscCall(PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE)); 2022 PetscCall(PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE)); 2023 PetscCall(PetscSFDestroy(&rsf)); 2024 PetscCall(PetscMalloc1(numLeavesNew, &localPointsNew)); 2025 PetscCall(PetscMalloc1(numLeavesNew, &remotePointsNew)); 2026 for (l = 0, m = 0; l < numLeaves; ++l) { 2027 const PetscInt p = localPoints[l]; 2028 DMPolytopeType ct; 2029 DMPolytopeType *rct; 2030 PetscInt *rsize, *rcone, *rornt; 2031 PetscInt Nct, n, r, q, off; 2032 2033 PetscCall(PetscSectionGetOffset(s, p, &off)); 2034 PetscCall(DMPlexGetCellType(dm, p, &ct)); 2035 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 2036 for (n = 0, q = 0; n < Nct; ++n) { 2037 for (r = 0; r < rsize[n]; ++r, ++m, ++q) { 2038 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 2039 localPointsNew[m] = pNew; 2040 remotePointsNew[m].index = rootPointsNew[off + q]; 2041 remotePointsNew[m].rank = remotePoints[l].rank; 2042 } 2043 } 2044 } 2045 PetscCall(PetscSectionDestroy(&s)); 2046 PetscCall(PetscFree(rootPointsNew)); 2047 /* SF needs sorted leaves to correctly calculate Gather */ 2048 { 2049 PetscSFNode *rp, *rtmp; 2050 PetscInt *lp, *idx, *ltmp, i; 2051 2052 PetscCall(PetscMalloc1(numLeavesNew, &idx)); 2053 PetscCall(PetscMalloc1(numLeavesNew, &lp)); 2054 PetscCall(PetscMalloc1(numLeavesNew, &rp)); 2055 for (i = 0; i < numLeavesNew; ++i) { 2056 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); 2057 idx[i] = i; 2058 } 2059 PetscCall(PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx)); 2060 for (i = 0; i < numLeavesNew; ++i) { 2061 lp[i] = localPointsNew[idx[i]]; 2062 rp[i] = remotePointsNew[idx[i]]; 2063 } 2064 ltmp = localPointsNew; 2065 localPointsNew = lp; 2066 rtmp = remotePointsNew; 2067 remotePointsNew = rp; 2068 PetscCall(PetscFree(idx)); 2069 PetscCall(PetscFree(ltmp)); 2070 PetscCall(PetscFree(rtmp)); 2071 } 2072 PetscCall(PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER)); 2073 PetscFunctionReturn(PETSC_SUCCESS); 2074 } 2075 2076 /* 2077 DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of `DMPolytopeType` `ct` with localized coordinates `x`, generate localized coordinates `xr` for subcell `r` of type `rct`. 2078 2079 Not Collective 2080 2081 Input Parameters: 2082 + tr - The `DMPlexTransform` 2083 . ct - The type of the parent cell 2084 . rct - The type of the produced cell 2085 . r - The index of the produced cell 2086 - x - The localized coordinates for the parent cell 2087 2088 Output Parameter: 2089 . xr - The localized coordinates for the produced cell 2090 2091 Level: developer 2092 2093 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexCellRefinerSetCoordinates()` 2094 */ 2095 static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[]) 2096 { 2097 PetscFE fe = NULL; 2098 PetscInt cdim, v, *subcellV; 2099 2100 PetscFunctionBegin; 2101 PetscCall(DMPlexTransformGetCoordinateFE(tr, ct, &fe)); 2102 PetscCall(DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV)); 2103 PetscCall(PetscFEGetNumComponents(fe, &cdim)); 2104 for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscCall(PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim])); 2105 PetscFunctionReturn(PETSC_SUCCESS); 2106 } 2107 2108 static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm) 2109 { 2110 DM dm, cdm, cdmCell, cdmNew, cdmCellNew; 2111 PetscSection coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew; 2112 Vec coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew; 2113 const PetscScalar *coords; 2114 PetscScalar *coordsNew; 2115 const PetscReal *maxCell, *Lstart, *L; 2116 PetscBool localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE; 2117 PetscInt dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p; 2118 2119 PetscFunctionBegin; 2120 // Need to clear the DMField for coordinates 2121 PetscCall(DMSetCoordinateField(rdm, NULL)); 2122 PetscCall(DMPlexTransformGetDM(tr, &dm)); 2123 PetscCall(DMGetCoordinateDM(dm, &cdm)); 2124 PetscCall(DMGetCellCoordinateDM(dm, &cdmCell)); 2125 PetscCall(DMGetCoordinatesLocalized(dm, &localized)); 2126 PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L)); 2127 if (localized) { 2128 /* Localize coordinates of new vertices */ 2129 localizeVertices = PETSC_TRUE; 2130 /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */ 2131 if (!maxCell) localizeCells = PETSC_TRUE; 2132 } 2133 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 2134 PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &dEo)); 2135 if (maxCell) { 2136 PetscReal maxCellNew[3]; 2137 2138 for (d = 0; d < dEo; ++d) maxCellNew[d] = maxCell[d] / 2.0; 2139 PetscCall(DMSetPeriodicity(rdm, maxCellNew, Lstart, L)); 2140 } 2141 PetscCall(DMGetCoordinateDim(rdm, &dE)); 2142 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew)); 2143 PetscCall(PetscSectionSetNumFields(coordSectionNew, 1)); 2144 PetscCall(PetscSectionSetFieldComponents(coordSectionNew, 0, dE)); 2145 PetscCall(DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew)); 2146 PetscCall(PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew)); 2147 /* Localization should be inherited */ 2148 /* Stefano calculates parent cells for each new cell for localization */ 2149 /* Localized cells need coordinates of closure */ 2150 for (v = vStartNew; v < vEndNew; ++v) { 2151 PetscCall(PetscSectionSetDof(coordSectionNew, v, dE)); 2152 PetscCall(PetscSectionSetFieldDof(coordSectionNew, v, 0, dE)); 2153 } 2154 PetscCall(PetscSectionSetUp(coordSectionNew)); 2155 PetscCall(DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew)); 2156 2157 if (localizeCells) { 2158 PetscCall(DMGetCoordinateDM(rdm, &cdmNew)); 2159 PetscCall(DMClone(cdmNew, &cdmCellNew)); 2160 PetscCall(DMSetCellCoordinateDM(rdm, cdmCellNew)); 2161 PetscCall(DMDestroy(&cdmCellNew)); 2162 2163 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew)); 2164 PetscCall(PetscSectionSetNumFields(coordSectionCellNew, 1)); 2165 PetscCall(PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE)); 2166 PetscCall(DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew)); 2167 PetscCall(PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew)); 2168 2169 PetscCall(DMGetCellCoordinateSection(dm, &coordSectionCell)); 2170 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 2171 for (c = cStart; c < cEnd; ++c) { 2172 PetscInt dof; 2173 2174 PetscCall(PetscSectionGetDof(coordSectionCell, c, &dof)); 2175 if (dof) { 2176 DMPolytopeType ct; 2177 DMPolytopeType *rct; 2178 PetscInt *rsize, *rcone, *rornt; 2179 PetscInt dim, cNew, Nct, n, r; 2180 2181 PetscCall(DMPlexGetCellType(dm, c, &ct)); 2182 dim = DMPolytopeTypeGetDim(ct); 2183 PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 2184 /* This allows for different cell types */ 2185 for (n = 0; n < Nct; ++n) { 2186 if (dim != DMPolytopeTypeGetDim(rct[n])) continue; 2187 for (r = 0; r < rsize[n]; ++r) { 2188 PetscInt *closure = NULL; 2189 PetscInt clSize, cl, Nv = 0; 2190 2191 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew)); 2192 PetscCall(DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure)); 2193 for (cl = 0; cl < clSize * 2; cl += 2) { 2194 if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv; 2195 } 2196 PetscCall(DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure)); 2197 PetscCall(PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE)); 2198 PetscCall(PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE)); 2199 } 2200 } 2201 } 2202 } 2203 PetscCall(PetscSectionSetUp(coordSectionCellNew)); 2204 PetscCall(DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew)); 2205 } 2206 PetscCall(DMViewFromOptions(dm, NULL, "-coarse_dm_view")); 2207 { 2208 VecType vtype; 2209 PetscInt coordSizeNew, bs; 2210 const char *name; 2211 2212 PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal)); 2213 PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalNew)); 2214 PetscCall(PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew)); 2215 PetscCall(VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE)); 2216 PetscCall(PetscObjectGetName((PetscObject)coordsLocal, &name)); 2217 PetscCall(PetscObjectSetName((PetscObject)coordsLocalNew, name)); 2218 PetscCall(VecGetBlockSize(coordsLocal, &bs)); 2219 PetscCall(VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE)); 2220 PetscCall(VecGetType(coordsLocal, &vtype)); 2221 PetscCall(VecSetType(coordsLocalNew, vtype)); 2222 } 2223 PetscCall(VecGetArrayRead(coordsLocal, &coords)); 2224 PetscCall(VecGetArray(coordsLocalNew, &coordsNew)); 2225 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 2226 /* First set coordinates for vertices */ 2227 for (p = pStart; p < pEnd; ++p) { 2228 DMPolytopeType ct; 2229 DMPolytopeType *rct; 2230 PetscInt *rsize, *rcone, *rornt; 2231 PetscInt Nct, n, r; 2232 PetscBool hasVertex = PETSC_FALSE; 2233 2234 PetscCall(DMPlexGetCellType(dm, p, &ct)); 2235 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 2236 for (n = 0; n < Nct; ++n) { 2237 if (rct[n] == DM_POLYTOPE_POINT) { 2238 hasVertex = PETSC_TRUE; 2239 break; 2240 } 2241 } 2242 if (hasVertex) { 2243 const PetscScalar *icoords = NULL; 2244 const PetscScalar *array = NULL; 2245 PetscScalar *pcoords = NULL; 2246 PetscBool isDG; 2247 PetscInt Nc, Nv, v, d; 2248 2249 PetscCall(DMPlexGetCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords)); 2250 2251 icoords = pcoords; 2252 Nv = Nc / dEo; 2253 if (ct != DM_POLYTOPE_POINT) { 2254 if (localizeVertices && maxCell) { 2255 PetscScalar anchor[3]; 2256 2257 for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d]; 2258 for (v = 0; v < Nv; ++v) PetscCall(DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo])); 2259 } 2260 } 2261 for (n = 0; n < Nct; ++n) { 2262 if (rct[n] != DM_POLYTOPE_POINT) continue; 2263 for (r = 0; r < rsize[n]; ++r) { 2264 PetscScalar vcoords[3]; 2265 PetscInt vNew, off; 2266 2267 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew)); 2268 PetscCall(PetscSectionGetOffset(coordSectionNew, vNew, &off)); 2269 PetscCall(DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords)); 2270 PetscCall(DMSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off])); 2271 } 2272 } 2273 PetscCall(DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords)); 2274 } 2275 } 2276 PetscCall(VecRestoreArrayRead(coordsLocal, &coords)); 2277 PetscCall(VecRestoreArray(coordsLocalNew, &coordsNew)); 2278 PetscCall(DMSetCoordinatesLocal(rdm, coordsLocalNew)); 2279 PetscCall(VecDestroy(&coordsLocalNew)); 2280 PetscCall(PetscSectionDestroy(&coordSectionNew)); 2281 /* Then set coordinates for cells by localizing */ 2282 if (!localizeCells) PetscCall(DMLocalizeCoordinates(rdm)); 2283 else { 2284 VecType vtype; 2285 PetscInt coordSizeNew, bs; 2286 const char *name; 2287 2288 PetscCall(DMGetCellCoordinatesLocal(dm, &coordsLocalCell)); 2289 PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew)); 2290 PetscCall(PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew)); 2291 PetscCall(VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE)); 2292 PetscCall(PetscObjectGetName((PetscObject)coordsLocalCell, &name)); 2293 PetscCall(PetscObjectSetName((PetscObject)coordsLocalCellNew, name)); 2294 PetscCall(VecGetBlockSize(coordsLocalCell, &bs)); 2295 PetscCall(VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE)); 2296 PetscCall(VecGetType(coordsLocalCell, &vtype)); 2297 PetscCall(VecSetType(coordsLocalCellNew, vtype)); 2298 PetscCall(VecGetArrayRead(coordsLocalCell, &coords)); 2299 PetscCall(VecGetArray(coordsLocalCellNew, &coordsNew)); 2300 2301 for (p = pStart; p < pEnd; ++p) { 2302 DMPolytopeType ct; 2303 DMPolytopeType *rct; 2304 PetscInt *rsize, *rcone, *rornt; 2305 PetscInt dof = 0, Nct, n, r; 2306 2307 PetscCall(DMPlexGetCellType(dm, p, &ct)); 2308 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 2309 if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(coordSectionCell, p, &dof)); 2310 if (dof) { 2311 const PetscScalar *pcoords; 2312 2313 PetscCall(DMPlexPointLocalRead(cdmCell, p, coords, &pcoords)); 2314 for (n = 0; n < Nct; ++n) { 2315 const PetscInt Nr = rsize[n]; 2316 2317 if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue; 2318 for (r = 0; r < Nr; ++r) { 2319 PetscInt pNew, offNew; 2320 2321 /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that 2322 DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger 2323 cell to the ones it produces. */ 2324 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 2325 PetscCall(PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew)); 2326 PetscCall(DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew])); 2327 } 2328 } 2329 } 2330 } 2331 PetscCall(VecRestoreArrayRead(coordsLocalCell, &coords)); 2332 PetscCall(VecRestoreArray(coordsLocalCellNew, &coordsNew)); 2333 PetscCall(DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew)); 2334 PetscCall(VecDestroy(&coordsLocalCellNew)); 2335 PetscCall(PetscSectionDestroy(&coordSectionCellNew)); 2336 } 2337 PetscFunctionReturn(PETSC_SUCCESS); 2338 } 2339 2340 /*@ 2341 DMPlexTransformApply - Execute the transformation, producing another `DM` 2342 2343 Collective 2344 2345 Input Parameters: 2346 + tr - The `DMPlexTransform` object 2347 - dm - The original `DM` 2348 2349 Output Parameter: 2350 . tdm - The transformed `DM` 2351 2352 Level: intermediate 2353 2354 Options Database Keys: 2355 + -dm_plex_transform_label_match_strata - Only label points of the same stratum as the producing point 2356 . -dm_plex_transform_label_replica_inc <num> - Increment for the label value to be multiplied by the replica number 2357 - -dm_plex_transform_active <name> - Name for active mesh label 2358 2359 .seealso: [](plex_transform_table), [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformCreate()`, `DMPlexTransformSetDM()` 2360 @*/ 2361 PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm) 2362 { 2363 DM rdm; 2364 DMPlexInterpolatedFlag interp; 2365 PetscInt pStart, pEnd; 2366 2367 PetscFunctionBegin; 2368 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 2369 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 2370 PetscAssertPointer(tdm, 3); 2371 PetscCall(PetscLogEventBegin(DMPLEX_Transform, tr, dm, 0, 0)); 2372 PetscCall(DMPlexTransformSetDM(tr, dm)); 2373 2374 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &rdm)); 2375 PetscCall(DMSetType(rdm, DMPLEX)); 2376 PetscCall(DMPlexTransformSetDimensions(tr, dm, rdm)); 2377 /* Calculate number of new points of each depth */ 2378 PetscCall(DMPlexIsInterpolatedCollective(dm, &interp)); 2379 PetscCheck(interp == DMPLEX_INTERPOLATED_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement"); 2380 /* Step 1: Set chart */ 2381 PetscCall(DMPlexTransformGetChart(tr, &pStart, &pEnd)); 2382 PetscCall(DMPlexSetChart(rdm, pStart, pEnd)); 2383 /* Step 2: Set cone/support sizes (automatically stratifies) */ 2384 PetscCall(DMPlexTransformSetConeSizes(tr, rdm)); 2385 /* Step 3: Setup refined DM */ 2386 PetscCall(DMSetUp(rdm)); 2387 /* Step 4: Set cones and supports (automatically symmetrizes) */ 2388 PetscCall(DMPlexTransformSetCones(tr, rdm)); 2389 /* Step 5: Create pointSF */ 2390 PetscCall(DMPlexTransformCreateSF(tr, rdm)); 2391 /* Step 6: Create labels */ 2392 PetscCall(DMPlexTransformCreateLabels(tr, rdm)); 2393 /* Step 7: Set coordinates */ 2394 PetscCall(DMPlexTransformSetCoordinates(tr, rdm)); 2395 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, rdm)); 2396 // If the original DM was configured from options, the transformed DM should be as well 2397 rdm->setfromoptionscalled = dm->setfromoptionscalled; 2398 PetscCall(PetscLogEventEnd(DMPLEX_Transform, tr, dm, 0, 0)); 2399 *tdm = rdm; 2400 PetscFunctionReturn(PETSC_SUCCESS); 2401 } 2402 2403 PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm) 2404 { 2405 DMPlexTransform tr; 2406 DM cdm, rcdm; 2407 const char *prefix; 2408 PetscBool save; 2409 2410 PetscFunctionBegin; 2411 PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr)); 2412 PetscCall(PetscObjectSetName((PetscObject)tr, "Adapt Label Transform")); 2413 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix)); 2414 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix)); 2415 PetscCall(DMPlexTransformSetDM(tr, dm)); 2416 PetscCall(DMPlexTransformSetFromOptions(tr)); 2417 PetscCall(DMPlexTransformSetActive(tr, adaptLabel)); 2418 PetscCall(DMPlexTransformSetUp(tr)); 2419 PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view")); 2420 PetscCall(DMPlexTransformApply(tr, dm, rdm)); 2421 PetscCall(DMCopyDisc(dm, *rdm)); 2422 PetscCall(DMGetCoordinateDM(dm, &cdm)); 2423 PetscCall(DMGetCoordinateDM(*rdm, &rcdm)); 2424 PetscCall(DMCopyDisc(cdm, rcdm)); 2425 PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm)); 2426 PetscCall(DMCopyDisc(dm, *rdm)); 2427 PetscCall(DMPlexGetSaveTransform(dm, &save)); 2428 if (save) PetscCall(DMPlexSetTransform(*rdm, tr)); 2429 PetscCall(DMPlexTransformDestroy(&tr)); 2430 ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation; 2431 PetscFunctionReturn(PETSC_SUCCESS); 2432 } 2433