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