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