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