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