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