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