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 Notes: 57 DMPlexTransformRegister() may be called multiple times to add several user-defined transforms 58 59 Sample usage: 60 .vb 61 DMPlexTransformRegister("my_transform", MyTransformCreate); 62 .ve 63 64 Then, your transform type can be chosen with the procedural interface via 65 .vb 66 DMPlexTransformCreate(MPI_Comm, DMPlexTransform *); 67 DMPlexTransformSetType(DMPlexTransform, "my_transform"); 68 .ve 69 or at runtime via the option 70 .vb 71 -dm_plex_transform_type my_transform 72 .ve 73 74 Level: advanced 75 76 .seealso: `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: `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 . It is called from PetscFinalize(). 123 124 Level: developer 125 126 .seealso: `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: `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 DMSetType - 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 Notes: 180 See "petsc/include/petscdmplextransform.h" for available transform types 181 182 Level: intermediate 183 184 .seealso: `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 DMPlexTransform type name 217 218 Level: intermediate 219 220 .seealso: `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 `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: 317 . -dm_plex_transform_type - Set the transform type, e.g. refine_regular 318 319 Level: intermediate 320 321 .seealso `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 transform. 344 345 Collective on tr 346 347 Input Parameter: 348 . tr - the transform object to destroy 349 350 Level: beginner 351 352 .seealso `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: `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: `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: `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 Note: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we 967 need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical 968 division (described in these outputs) of the cell in the original mesh. The point identifier is given by 969 $ the number of cones to be taken, or 0 for the current cell 970 $ the cell cone point number at each level from which it is subdivided 971 $ the replica number r of the subdivision. 972 The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is 973 $ Nt = 2 974 $ target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT} 975 $ size = {1, 2} 976 $ cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0} 977 $ ornt = { 0, 0, 0, 0} 978 979 Level: advanced 980 981 .seealso: `DMPlexTransformApply()`, `DMPlexTransformCreate()` 982 @*/ 983 PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 984 { 985 PetscFunctionBegin; 986 PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt); 987 PetscFunctionReturn(0); 988 } 989 990 PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 991 { 992 PetscFunctionBegin; 993 *rnew = r; 994 *onew = DMPolytopeTypeComposeOrientation(tct, o, so); 995 PetscFunctionReturn(0); 996 } 997 998 /* Returns the same thing */ 999 PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 1000 { 1001 static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT}; 1002 static PetscInt vertexS[] = {1}; 1003 static PetscInt vertexC[] = {0}; 1004 static PetscInt vertexO[] = {0}; 1005 static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT}; 1006 static PetscInt edgeS[] = {1}; 1007 static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 1008 static PetscInt edgeO[] = {0, 0}; 1009 static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR}; 1010 static PetscInt tedgeS[] = {1}; 1011 static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 1012 static PetscInt tedgeO[] = {0, 0}; 1013 static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE}; 1014 static PetscInt triS[] = {1}; 1015 static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0}; 1016 static PetscInt triO[] = {0, 0, 0}; 1017 static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL}; 1018 static PetscInt quadS[] = {1}; 1019 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}; 1020 static PetscInt quadO[] = {0, 0, 0, 0}; 1021 static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR}; 1022 static PetscInt tquadS[] = {1}; 1023 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}; 1024 static PetscInt tquadO[] = {0, 0, 0, 0}; 1025 static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON}; 1026 static PetscInt tetS[] = {1}; 1027 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}; 1028 static PetscInt tetO[] = {0, 0, 0, 0}; 1029 static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON}; 1030 static PetscInt hexS[] = {1}; 1031 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}; 1032 static PetscInt hexO[] = {0, 0, 0, 0, 0, 0}; 1033 static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM}; 1034 static PetscInt tripS[] = {1}; 1035 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}; 1036 static PetscInt tripO[] = {0, 0, 0, 0, 0}; 1037 static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR}; 1038 static PetscInt ttripS[] = {1}; 1039 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}; 1040 static PetscInt ttripO[] = {0, 0, 0, 0, 0}; 1041 static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR}; 1042 static PetscInt tquadpS[] = {1}; 1043 static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, 1044 DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0}; 1045 static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0}; 1046 static DMPolytopeType pyrT[] = {DM_POLYTOPE_PYRAMID}; 1047 static PetscInt pyrS[] = {1}; 1048 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}; 1049 static PetscInt pyrO[] = {0, 0, 0, 0, 0}; 1050 1051 PetscFunctionBegin; 1052 if (rt) *rt = 0; 1053 switch (source) { 1054 case DM_POLYTOPE_POINT: 1055 *Nt = 1; 1056 *target = vertexT; 1057 *size = vertexS; 1058 *cone = vertexC; 1059 *ornt = vertexO; 1060 break; 1061 case DM_POLYTOPE_SEGMENT: 1062 *Nt = 1; 1063 *target = edgeT; 1064 *size = edgeS; 1065 *cone = edgeC; 1066 *ornt = edgeO; 1067 break; 1068 case DM_POLYTOPE_POINT_PRISM_TENSOR: 1069 *Nt = 1; 1070 *target = tedgeT; 1071 *size = tedgeS; 1072 *cone = tedgeC; 1073 *ornt = tedgeO; 1074 break; 1075 case DM_POLYTOPE_TRIANGLE: 1076 *Nt = 1; 1077 *target = triT; 1078 *size = triS; 1079 *cone = triC; 1080 *ornt = triO; 1081 break; 1082 case DM_POLYTOPE_QUADRILATERAL: 1083 *Nt = 1; 1084 *target = quadT; 1085 *size = quadS; 1086 *cone = quadC; 1087 *ornt = quadO; 1088 break; 1089 case DM_POLYTOPE_SEG_PRISM_TENSOR: 1090 *Nt = 1; 1091 *target = tquadT; 1092 *size = tquadS; 1093 *cone = tquadC; 1094 *ornt = tquadO; 1095 break; 1096 case DM_POLYTOPE_TETRAHEDRON: 1097 *Nt = 1; 1098 *target = tetT; 1099 *size = tetS; 1100 *cone = tetC; 1101 *ornt = tetO; 1102 break; 1103 case DM_POLYTOPE_HEXAHEDRON: 1104 *Nt = 1; 1105 *target = hexT; 1106 *size = hexS; 1107 *cone = hexC; 1108 *ornt = hexO; 1109 break; 1110 case DM_POLYTOPE_TRI_PRISM: 1111 *Nt = 1; 1112 *target = tripT; 1113 *size = tripS; 1114 *cone = tripC; 1115 *ornt = tripO; 1116 break; 1117 case DM_POLYTOPE_TRI_PRISM_TENSOR: 1118 *Nt = 1; 1119 *target = ttripT; 1120 *size = ttripS; 1121 *cone = ttripC; 1122 *ornt = ttripO; 1123 break; 1124 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 1125 *Nt = 1; 1126 *target = tquadpT; 1127 *size = tquadpS; 1128 *cone = tquadpC; 1129 *ornt = tquadpO; 1130 break; 1131 case DM_POLYTOPE_PYRAMID: 1132 *Nt = 1; 1133 *target = pyrT; 1134 *size = pyrS; 1135 *cone = pyrC; 1136 *ornt = pyrO; 1137 break; 1138 default: 1139 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 1140 } 1141 PetscFunctionReturn(0); 1142 } 1143 1144 /*@ 1145 DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point 1146 1147 Not collective 1148 1149 Input Parameters: 1150 + tr - The DMPlexTransform 1151 . sct - The source point cell type, from whom the new cell is being produced 1152 . sp - The source point 1153 . so - The orientation of the source point in its enclosing parent 1154 . tct - The target point cell type 1155 . r - The replica number requested for the produced cell type 1156 - o - The orientation of the replica 1157 1158 Output Parameters: 1159 + rnew - The replica number, given the orientation of the parent 1160 - onew - The replica orientation, given the orientation of the parent 1161 1162 Level: advanced 1163 1164 .seealso: `DMPlexTransformCellTransform()`, `DMPlexTransformApply()` 1165 @*/ 1166 PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 1167 { 1168 PetscFunctionBeginHot; 1169 PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew); 1170 PetscFunctionReturn(0); 1171 } 1172 1173 static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm) 1174 { 1175 DM dm; 1176 PetscInt pStart, pEnd, p, pNew; 1177 1178 PetscFunctionBegin; 1179 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1180 /* Must create the celltype label here so that we do not automatically try to compute the types */ 1181 PetscCall(DMCreateLabel(rdm, "celltype")); 1182 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1183 for (p = pStart; p < pEnd; ++p) { 1184 DMPolytopeType ct; 1185 DMPolytopeType *rct; 1186 PetscInt *rsize, *rcone, *rornt; 1187 PetscInt Nct, n, r; 1188 1189 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1190 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1191 for (n = 0; n < Nct; ++n) { 1192 for (r = 0; r < rsize[n]; ++r) { 1193 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 1194 PetscCall(DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]))); 1195 PetscCall(DMPlexSetCellType(rdm, pNew, rct[n])); 1196 } 1197 } 1198 } 1199 /* Let the DM know we have set all the cell types */ 1200 { 1201 DMLabel ctLabel; 1202 DM_Plex *plex = (DM_Plex *)rdm->data; 1203 1204 PetscCall(DMPlexGetCellTypeLabel(rdm, &ctLabel)); 1205 PetscCall(PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState)); 1206 } 1207 PetscFunctionReturn(0); 1208 } 1209 1210 PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize) 1211 { 1212 DMPolytopeType ctNew; 1213 1214 PetscFunctionBegin; 1215 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1216 PetscValidIntPointer(coneSize, 3); 1217 PetscCall(DMPlexTransformGetCellType(tr, q, &ctNew)); 1218 *coneSize = DMPolytopeTypeGetConeSize((DMPolytopeType)ctNew); 1219 PetscFunctionReturn(0); 1220 } 1221 1222 /* The orientation o is for the interior of the cell p */ 1223 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[]) 1224 { 1225 DM dm; 1226 const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); 1227 const PetscInt *cone; 1228 PetscInt c, coff = *coneoff, ooff = *orntoff; 1229 1230 PetscFunctionBegin; 1231 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1232 PetscCall(DMPlexGetOrientedCone(dm, p, &cone, NULL)); 1233 for (c = 0; c < csizeNew; ++c) { 1234 PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */ 1235 PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */ 1236 PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */ 1237 DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */ 1238 const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */ 1239 PetscInt pr = -1; /* Replica number of pp that produces new cone point */ 1240 const DMPolytopeType ft = (DMPolytopeType)rcone[coff++]; /* Cell type for new cone point of pNew */ 1241 const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ 1242 PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */ 1243 PetscInt lc; 1244 1245 /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */ 1246 for (lc = 0; lc < fn; ++lc) { 1247 const PetscInt *parr = DMPolytopeTypeGetArrangment(pct, po); 1248 const PetscInt acp = rcone[coff++]; 1249 const PetscInt pcp = parr[acp * 2]; 1250 const PetscInt pco = parr[acp * 2 + 1]; 1251 const PetscInt *ppornt; 1252 1253 ppp = pp; 1254 pp = pcone[pcp]; 1255 PetscCall(DMPlexGetCellType(dm, pp, &pct)); 1256 // Restore the parent cone from the last iterate 1257 if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, ppp, &pcone, NULL)); 1258 PetscCall(DMPlexGetOrientedCone(dm, pp, &pcone, NULL)); 1259 PetscCall(DMPlexGetOrientedCone(dm, ppp, NULL, &ppornt)); 1260 po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco); 1261 PetscCall(DMPlexRestoreOrientedCone(dm, ppp, NULL, &ppornt)); 1262 } 1263 if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, pp, &pcone, NULL)); 1264 pr = rcone[coff++]; 1265 /* Orientation po of pp maps (pr, fo) -> (pr', fo') */ 1266 PetscCall(DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo)); 1267 PetscCall(DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c])); 1268 orntNew[c] = fo; 1269 } 1270 PetscCall(DMPlexRestoreOrientedCone(dm, p, &cone, NULL)); 1271 *coneoff = coff; 1272 *orntoff = ooff; 1273 PetscFunctionReturn(0); 1274 } 1275 1276 static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm) 1277 { 1278 DM dm; 1279 DMPolytopeType ct; 1280 PetscInt *coneNew, *orntNew; 1281 PetscInt maxConeSize = 0, pStart, pEnd, p, pNew; 1282 1283 PetscFunctionBegin; 1284 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1285 for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p)); 1286 PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew)); 1287 PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew)); 1288 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1289 for (p = pStart; p < pEnd; ++p) { 1290 PetscInt coff, ooff; 1291 DMPolytopeType *rct; 1292 PetscInt *rsize, *rcone, *rornt; 1293 PetscInt Nct, n, r; 1294 1295 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1296 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1297 for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) { 1298 const DMPolytopeType ctNew = rct[n]; 1299 1300 for (r = 0; r < rsize[n]; ++r) { 1301 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 1302 PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew)); 1303 PetscCall(DMPlexSetCone(rdm, pNew, coneNew)); 1304 PetscCall(DMPlexSetConeOrientation(rdm, pNew, orntNew)); 1305 } 1306 } 1307 } 1308 PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew)); 1309 PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew)); 1310 PetscCall(DMViewFromOptions(rdm, NULL, "-rdm_view")); 1311 PetscCall(DMPlexSymmetrize(rdm)); 1312 PetscCall(DMPlexStratify(rdm)); 1313 PetscFunctionReturn(0); 1314 } 1315 1316 PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[]) 1317 { 1318 DM dm; 1319 DMPolytopeType ct, qct; 1320 DMPolytopeType *rct; 1321 PetscInt *rsize, *rcone, *rornt, *qcone, *qornt; 1322 PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0; 1323 1324 PetscFunctionBegin; 1325 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1326 PetscValidPointer(cone, 4); 1327 PetscValidPointer(ornt, 5); 1328 for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p)); 1329 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1330 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone)); 1331 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt)); 1332 PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r)); 1333 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1334 for (n = 0; n < Nct; ++n) { 1335 const DMPolytopeType ctNew = rct[n]; 1336 const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); 1337 PetscInt Nr = rsize[n], fn, c; 1338 1339 if (ctNew == qct) Nr = r; 1340 for (nr = 0; nr < Nr; ++nr) { 1341 for (c = 0; c < csizeNew; ++c) { 1342 ++coff; /* Cell type of new cone point */ 1343 fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ 1344 coff += fn; 1345 ++coff; /* Replica number of new cone point */ 1346 ++ooff; /* Orientation of new cone point */ 1347 } 1348 } 1349 if (ctNew == qct) break; 1350 } 1351 PetscCall(DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt)); 1352 *cone = qcone; 1353 *ornt = qornt; 1354 PetscFunctionReturn(0); 1355 } 1356 1357 PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[]) 1358 { 1359 DM dm; 1360 DMPolytopeType ct, qct; 1361 DMPolytopeType *rct; 1362 PetscInt *rsize, *rcone, *rornt, *qcone, *qornt; 1363 PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0; 1364 1365 PetscFunctionBegin; 1366 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1367 if (cone) PetscValidPointer(cone, 3); 1368 if (ornt) PetscValidPointer(ornt, 4); 1369 for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p)); 1370 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1371 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone)); 1372 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt)); 1373 PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r)); 1374 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1375 for (n = 0; n < Nct; ++n) { 1376 const DMPolytopeType ctNew = rct[n]; 1377 const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); 1378 PetscInt Nr = rsize[n], fn, c; 1379 1380 if (ctNew == qct) Nr = r; 1381 for (nr = 0; nr < Nr; ++nr) { 1382 for (c = 0; c < csizeNew; ++c) { 1383 ++coff; /* Cell type of new cone point */ 1384 fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ 1385 coff += fn; 1386 ++coff; /* Replica number of new cone point */ 1387 ++ooff; /* Orientation of new cone point */ 1388 } 1389 } 1390 if (ctNew == qct) break; 1391 } 1392 PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt)); 1393 if (cone) *cone = qcone; 1394 else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qcone)); 1395 if (ornt) *ornt = qornt; 1396 else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qornt)); 1397 PetscFunctionReturn(0); 1398 } 1399 1400 PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[]) 1401 { 1402 DM dm; 1403 1404 PetscFunctionBegin; 1405 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1406 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1407 if (cone) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, cone)); 1408 if (ornt) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, ornt)); 1409 PetscFunctionReturn(0); 1410 } 1411 1412 static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr) 1413 { 1414 PetscInt ict; 1415 1416 PetscFunctionBegin; 1417 PetscCall(PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts)); 1418 for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) { 1419 const DMPolytopeType ct = (DMPolytopeType)ict; 1420 DMPlexTransform reftr; 1421 DM refdm, trdm; 1422 Vec coordinates; 1423 const PetscScalar *coords; 1424 DMPolytopeType *rct; 1425 PetscInt *rsize, *rcone, *rornt; 1426 PetscInt Nct, n, r, pNew; 1427 PetscInt trdim, vStart, vEnd, Nc; 1428 const PetscInt debug = 0; 1429 const char *typeName; 1430 1431 /* Since points are 0-dimensional, coordinates make no sense */ 1432 if (DMPolytopeTypeGetDim(ct) <= 0) continue; 1433 PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm)); 1434 PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &reftr)); 1435 PetscCall(DMPlexTransformSetDM(reftr, refdm)); 1436 PetscCall(DMPlexTransformGetType(tr, &typeName)); 1437 PetscCall(DMPlexTransformSetType(reftr, typeName)); 1438 PetscCall(DMPlexTransformSetUp(reftr)); 1439 PetscCall(DMPlexTransformApply(reftr, refdm, &trdm)); 1440 1441 PetscCall(DMGetDimension(trdm, &trdim)); 1442 PetscCall(DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd)); 1443 tr->trNv[ct] = vEnd - vStart; 1444 PetscCall(DMGetCoordinatesLocal(trdm, &coordinates)); 1445 PetscCall(VecGetLocalSize(coordinates, &Nc)); 1446 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); 1447 PetscCall(PetscCalloc1(Nc, &tr->trVerts[ct])); 1448 PetscCall(VecGetArrayRead(coordinates, &coords)); 1449 PetscCall(PetscArraycpy(tr->trVerts[ct], coords, Nc)); 1450 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 1451 1452 PetscCall(PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct])); 1453 PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1454 for (n = 0; n < Nct; ++n) { 1455 /* Since points are 0-dimensional, coordinates make no sense */ 1456 if (rct[n] == DM_POLYTOPE_POINT) continue; 1457 PetscCall(PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]])); 1458 for (r = 0; r < rsize[n]; ++r) { 1459 PetscInt *closure = NULL; 1460 PetscInt clSize, cl, Nv = 0; 1461 1462 PetscCall(PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r])); 1463 PetscCall(DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew)); 1464 PetscCall(DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure)); 1465 for (cl = 0; cl < clSize * 2; cl += 2) { 1466 const PetscInt sv = closure[cl]; 1467 1468 if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart; 1469 } 1470 PetscCall(DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure)); 1471 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]); 1472 } 1473 } 1474 if (debug) { 1475 DMPolytopeType *rct; 1476 PetscInt *rsize, *rcone, *rornt; 1477 PetscInt v, dE = trdim, d, off = 0; 1478 1479 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct])); 1480 for (v = 0; v < tr->trNv[ct]; ++v) { 1481 PetscCall(PetscPrintf(PETSC_COMM_SELF, " ")); 1482 for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++]))); 1483 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 1484 } 1485 1486 PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1487 for (n = 0; n < Nct; ++n) { 1488 if (rct[n] == DM_POLYTOPE_POINT) continue; 1489 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct])); 1490 for (r = 0; r < rsize[n]; ++r) { 1491 PetscCall(PetscPrintf(PETSC_COMM_SELF, " ")); 1492 for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v])); 1493 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 1494 } 1495 } 1496 } 1497 PetscCall(DMDestroy(&refdm)); 1498 PetscCall(DMDestroy(&trdm)); 1499 PetscCall(DMPlexTransformDestroy(&reftr)); 1500 } 1501 PetscFunctionReturn(0); 1502 } 1503 1504 /* 1505 DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type 1506 1507 Input Parameters: 1508 + tr - The DMPlexTransform object 1509 - ct - The cell type 1510 1511 Output Parameters: 1512 + Nv - The number of transformed vertices in the closure of the reference cell of given type 1513 - trVerts - The coordinates of these vertices in the reference cell 1514 1515 Level: developer 1516 1517 .seealso: `DMPlexTransformGetSubcellVertices()` 1518 */ 1519 PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[]) 1520 { 1521 PetscFunctionBegin; 1522 if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr)); 1523 if (Nv) *Nv = tr->trNv[ct]; 1524 if (trVerts) *trVerts = tr->trVerts[ct]; 1525 PetscFunctionReturn(0); 1526 } 1527 1528 /* 1529 DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type 1530 1531 Input Parameters: 1532 + tr - The DMPlexTransform object 1533 . ct - The cell type 1534 . rct - The subcell type 1535 - r - The subcell index 1536 1537 Output Parameter: 1538 . subVerts - The indices of these vertices in the set of vertices returned by DMPlexTransformGetCellVertices() 1539 1540 Level: developer 1541 1542 .seealso: `DMPlexTransformGetCellVertices()` 1543 */ 1544 PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[]) 1545 { 1546 PetscFunctionBegin; 1547 if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr)); 1548 PetscCheck(tr->trSubVerts[ct][rct], PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]); 1549 if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r]; 1550 PetscFunctionReturn(0); 1551 } 1552 1553 /* Computes new vertex as the barycenter, or centroid */ 1554 PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) 1555 { 1556 PetscInt v, d; 1557 1558 PetscFunctionBeginHot; 1559 PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]); 1560 for (d = 0; d < dE; ++d) out[d] = 0.0; 1561 for (v = 0; v < Nv; ++v) 1562 for (d = 0; d < dE; ++d) out[d] += in[v * dE + d]; 1563 for (d = 0; d < dE; ++d) out[d] /= Nv; 1564 PetscFunctionReturn(0); 1565 } 1566 1567 /*@ 1568 DMPlexTransformMapCoordinates - Calculate new coordinates for produced points 1569 1570 Not collective 1571 1572 Input Parameters: 1573 + cr - The DMPlexCellRefiner 1574 . pct - The cell type of the parent, from whom the new cell is being produced 1575 . ct - The type being produced 1576 . p - The original point 1577 . r - The replica number requested for the produced cell type 1578 . Nv - Number of vertices in the closure of the parent cell 1579 . dE - Spatial dimension 1580 - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell 1581 1582 Output Parameter: 1583 . out - The coordinates of the new vertices 1584 1585 Level: intermediate 1586 1587 .seealso: `DMPlexTransformApply()` 1588 @*/ 1589 PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) 1590 { 1591 PetscFunctionBeginHot; 1592 PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out); 1593 PetscFunctionReturn(0); 1594 } 1595 1596 static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew) 1597 { 1598 DM dm; 1599 IS valueIS; 1600 const PetscInt *values; 1601 PetscInt defVal, Nv, val; 1602 1603 PetscFunctionBegin; 1604 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1605 PetscCall(DMLabelGetDefaultValue(label, &defVal)); 1606 PetscCall(DMLabelSetDefaultValue(labelNew, defVal)); 1607 PetscCall(DMLabelGetValueIS(label, &valueIS)); 1608 PetscCall(ISGetLocalSize(valueIS, &Nv)); 1609 PetscCall(ISGetIndices(valueIS, &values)); 1610 for (val = 0; val < Nv; ++val) { 1611 IS pointIS; 1612 const PetscInt *points; 1613 PetscInt numPoints, p; 1614 1615 /* Ensure refined label is created with same number of strata as 1616 * original (even if no entries here). */ 1617 PetscCall(DMLabelAddStratum(labelNew, values[val])); 1618 PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS)); 1619 PetscCall(ISGetLocalSize(pointIS, &numPoints)); 1620 PetscCall(ISGetIndices(pointIS, &points)); 1621 for (p = 0; p < numPoints; ++p) { 1622 const PetscInt point = points[p]; 1623 DMPolytopeType ct; 1624 DMPolytopeType *rct; 1625 PetscInt *rsize, *rcone, *rornt; 1626 PetscInt Nct, n, r, pNew = 0; 1627 1628 PetscCall(DMPlexGetCellType(dm, point, &ct)); 1629 PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1630 for (n = 0; n < Nct; ++n) { 1631 for (r = 0; r < rsize[n]; ++r) { 1632 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew)); 1633 PetscCall(DMLabelSetValue(labelNew, pNew, values[val])); 1634 } 1635 } 1636 } 1637 PetscCall(ISRestoreIndices(pointIS, &points)); 1638 PetscCall(ISDestroy(&pointIS)); 1639 } 1640 PetscCall(ISRestoreIndices(valueIS, &values)); 1641 PetscCall(ISDestroy(&valueIS)); 1642 PetscFunctionReturn(0); 1643 } 1644 1645 static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm) 1646 { 1647 DM dm; 1648 PetscInt numLabels, l; 1649 1650 PetscFunctionBegin; 1651 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1652 // TODO Figure out labels for ephemeral meshes 1653 if (((DM_Plex *)dm->data)->tr) PetscFunctionReturn(0); 1654 PetscCall(DMGetNumLabels(dm, &numLabels)); 1655 for (l = 0; l < numLabels; ++l) { 1656 DMLabel label, labelNew; 1657 const char *lname; 1658 PetscBool isDepth, isCellType; 1659 1660 PetscCall(DMGetLabelName(dm, l, &lname)); 1661 PetscCall(PetscStrcmp(lname, "depth", &isDepth)); 1662 if (isDepth) continue; 1663 PetscCall(PetscStrcmp(lname, "celltype", &isCellType)); 1664 if (isCellType) continue; 1665 PetscCall(DMCreateLabel(rdm, lname)); 1666 PetscCall(DMGetLabel(dm, lname, &label)); 1667 PetscCall(DMGetLabel(rdm, lname, &labelNew)); 1668 PetscCall(RefineLabel_Internal(tr, label, labelNew)); 1669 } 1670 PetscFunctionReturn(0); 1671 } 1672 1673 /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */ 1674 PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm) 1675 { 1676 DM dm; 1677 PetscInt Nf, f, Nds, s; 1678 1679 PetscFunctionBegin; 1680 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1681 PetscCall(DMGetNumFields(dm, &Nf)); 1682 for (f = 0; f < Nf; ++f) { 1683 DMLabel label, labelNew; 1684 PetscObject obj; 1685 const char *lname; 1686 1687 PetscCall(DMGetField(rdm, f, &label, &obj)); 1688 if (!label) continue; 1689 PetscCall(PetscObjectGetName((PetscObject)label, &lname)); 1690 PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew)); 1691 PetscCall(RefineLabel_Internal(tr, label, labelNew)); 1692 PetscCall(DMSetField_Internal(rdm, f, labelNew, obj)); 1693 PetscCall(DMLabelDestroy(&labelNew)); 1694 } 1695 PetscCall(DMGetNumDS(dm, &Nds)); 1696 for (s = 0; s < Nds; ++s) { 1697 DMLabel label, labelNew; 1698 const char *lname; 1699 1700 PetscCall(DMGetRegionNumDS(rdm, s, &label, NULL, NULL)); 1701 if (!label) continue; 1702 PetscCall(PetscObjectGetName((PetscObject)label, &lname)); 1703 PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew)); 1704 PetscCall(RefineLabel_Internal(tr, label, labelNew)); 1705 PetscCall(DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL)); 1706 PetscCall(DMLabelDestroy(&labelNew)); 1707 } 1708 PetscFunctionReturn(0); 1709 } 1710 1711 static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm) 1712 { 1713 DM dm; 1714 PetscSF sf, sfNew; 1715 PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m; 1716 const PetscInt *localPoints; 1717 const PetscSFNode *remotePoints; 1718 PetscInt *localPointsNew; 1719 PetscSFNode *remotePointsNew; 1720 PetscInt pStartNew, pEndNew, pNew; 1721 /* Brute force algorithm */ 1722 PetscSF rsf; 1723 PetscSection s; 1724 const PetscInt *rootdegree; 1725 PetscInt *rootPointsNew, *remoteOffsets; 1726 PetscInt numPointsNew, pStart, pEnd, p; 1727 1728 PetscFunctionBegin; 1729 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1730 PetscCall(DMPlexGetChart(rdm, &pStartNew, &pEndNew)); 1731 PetscCall(DMGetPointSF(dm, &sf)); 1732 PetscCall(DMGetPointSF(rdm, &sfNew)); 1733 /* Calculate size of new SF */ 1734 PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints)); 1735 if (numRoots < 0) PetscFunctionReturn(0); 1736 for (l = 0; l < numLeaves; ++l) { 1737 const PetscInt p = localPoints[l]; 1738 DMPolytopeType ct; 1739 DMPolytopeType *rct; 1740 PetscInt *rsize, *rcone, *rornt; 1741 PetscInt Nct, n; 1742 1743 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1744 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1745 for (n = 0; n < Nct; ++n) numLeavesNew += rsize[n]; 1746 } 1747 /* Send new root point numbers 1748 It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication 1749 */ 1750 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1751 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s)); 1752 PetscCall(PetscSectionSetChart(s, pStart, pEnd)); 1753 for (p = pStart; p < pEnd; ++p) { 1754 DMPolytopeType ct; 1755 DMPolytopeType *rct; 1756 PetscInt *rsize, *rcone, *rornt; 1757 PetscInt Nct, n; 1758 1759 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1760 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1761 for (n = 0; n < Nct; ++n) PetscCall(PetscSectionAddDof(s, p, rsize[n])); 1762 } 1763 PetscCall(PetscSectionSetUp(s)); 1764 PetscCall(PetscSectionGetStorageSize(s, &numPointsNew)); 1765 PetscCall(PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets)); 1766 PetscCall(PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf)); 1767 PetscCall(PetscFree(remoteOffsets)); 1768 PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree)); 1769 PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree)); 1770 PetscCall(PetscMalloc1(numPointsNew, &rootPointsNew)); 1771 for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1; 1772 for (p = pStart; p < pEnd; ++p) { 1773 DMPolytopeType ct; 1774 DMPolytopeType *rct; 1775 PetscInt *rsize, *rcone, *rornt; 1776 PetscInt Nct, n, r, off; 1777 1778 if (!rootdegree[p - pStart]) continue; 1779 PetscCall(PetscSectionGetOffset(s, p, &off)); 1780 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1781 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1782 for (n = 0, m = 0; n < Nct; ++n) { 1783 for (r = 0; r < rsize[n]; ++r, ++m) { 1784 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 1785 rootPointsNew[off + m] = pNew; 1786 } 1787 } 1788 } 1789 PetscCall(PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE)); 1790 PetscCall(PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE)); 1791 PetscCall(PetscSFDestroy(&rsf)); 1792 PetscCall(PetscMalloc1(numLeavesNew, &localPointsNew)); 1793 PetscCall(PetscMalloc1(numLeavesNew, &remotePointsNew)); 1794 for (l = 0, m = 0; l < numLeaves; ++l) { 1795 const PetscInt p = localPoints[l]; 1796 DMPolytopeType ct; 1797 DMPolytopeType *rct; 1798 PetscInt *rsize, *rcone, *rornt; 1799 PetscInt Nct, n, r, q, off; 1800 1801 PetscCall(PetscSectionGetOffset(s, p, &off)); 1802 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1803 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1804 for (n = 0, q = 0; n < Nct; ++n) { 1805 for (r = 0; r < rsize[n]; ++r, ++m, ++q) { 1806 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 1807 localPointsNew[m] = pNew; 1808 remotePointsNew[m].index = rootPointsNew[off + q]; 1809 remotePointsNew[m].rank = remotePoints[l].rank; 1810 } 1811 } 1812 } 1813 PetscCall(PetscSectionDestroy(&s)); 1814 PetscCall(PetscFree(rootPointsNew)); 1815 /* SF needs sorted leaves to correctly calculate Gather */ 1816 { 1817 PetscSFNode *rp, *rtmp; 1818 PetscInt *lp, *idx, *ltmp, i; 1819 1820 PetscCall(PetscMalloc1(numLeavesNew, &idx)); 1821 PetscCall(PetscMalloc1(numLeavesNew, &lp)); 1822 PetscCall(PetscMalloc1(numLeavesNew, &rp)); 1823 for (i = 0; i < numLeavesNew; ++i) { 1824 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); 1825 idx[i] = i; 1826 } 1827 PetscCall(PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx)); 1828 for (i = 0; i < numLeavesNew; ++i) { 1829 lp[i] = localPointsNew[idx[i]]; 1830 rp[i] = remotePointsNew[idx[i]]; 1831 } 1832 ltmp = localPointsNew; 1833 localPointsNew = lp; 1834 rtmp = remotePointsNew; 1835 remotePointsNew = rp; 1836 PetscCall(PetscFree(idx)); 1837 PetscCall(PetscFree(ltmp)); 1838 PetscCall(PetscFree(rtmp)); 1839 } 1840 PetscCall(PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER)); 1841 PetscFunctionReturn(0); 1842 } 1843 1844 /* 1845 DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of type ct with localized coordinates x, we generate localized coordinates xr for subcell r of type rct. 1846 1847 Not collective 1848 1849 Input Parameters: 1850 + cr - The DMPlexCellRefiner 1851 . ct - The type of the parent cell 1852 . rct - The type of the produced cell 1853 . r - The index of the produced cell 1854 - x - The localized coordinates for the parent cell 1855 1856 Output Parameter: 1857 . xr - The localized coordinates for the produced cell 1858 1859 Level: developer 1860 1861 .seealso: `DMPlexCellRefinerSetCoordinates()` 1862 */ 1863 static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[]) 1864 { 1865 PetscFE fe = NULL; 1866 PetscInt cdim, v, *subcellV; 1867 1868 PetscFunctionBegin; 1869 PetscCall(DMPlexTransformGetCoordinateFE(tr, ct, &fe)); 1870 PetscCall(DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV)); 1871 PetscCall(PetscFEGetNumComponents(fe, &cdim)); 1872 for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) PetscCall(PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim])); 1873 PetscFunctionReturn(0); 1874 } 1875 1876 static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm) 1877 { 1878 DM dm, cdm, cdmCell, cdmNew, cdmCellNew; 1879 PetscSection coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew; 1880 Vec coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew; 1881 const PetscScalar *coords; 1882 PetscScalar *coordsNew; 1883 const PetscReal *maxCell, *Lstart, *L; 1884 PetscBool localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE; 1885 PetscInt dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p; 1886 1887 PetscFunctionBegin; 1888 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1889 PetscCall(DMGetCoordinateDM(dm, &cdm)); 1890 PetscCall(DMGetCellCoordinateDM(dm, &cdmCell)); 1891 PetscCall(DMGetCoordinatesLocalized(dm, &localized)); 1892 PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L)); 1893 if (localized) { 1894 /* Localize coordinates of new vertices */ 1895 localizeVertices = PETSC_TRUE; 1896 /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */ 1897 if (!maxCell) localizeCells = PETSC_TRUE; 1898 } 1899 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 1900 PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &dEo)); 1901 if (maxCell) { 1902 PetscReal maxCellNew[3]; 1903 1904 for (d = 0; d < dEo; ++d) maxCellNew[d] = maxCell[d] / 2.0; 1905 PetscCall(DMSetPeriodicity(rdm, maxCellNew, Lstart, L)); 1906 } 1907 PetscCall(DMGetCoordinateDim(rdm, &dE)); 1908 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew)); 1909 PetscCall(PetscSectionSetNumFields(coordSectionNew, 1)); 1910 PetscCall(PetscSectionSetFieldComponents(coordSectionNew, 0, dE)); 1911 PetscCall(DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew)); 1912 PetscCall(PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew)); 1913 /* Localization should be inherited */ 1914 /* Stefano calculates parent cells for each new cell for localization */ 1915 /* Localized cells need coordinates of closure */ 1916 for (v = vStartNew; v < vEndNew; ++v) { 1917 PetscCall(PetscSectionSetDof(coordSectionNew, v, dE)); 1918 PetscCall(PetscSectionSetFieldDof(coordSectionNew, v, 0, dE)); 1919 } 1920 PetscCall(PetscSectionSetUp(coordSectionNew)); 1921 PetscCall(DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew)); 1922 1923 if (localizeCells) { 1924 PetscCall(DMGetCoordinateDM(rdm, &cdmNew)); 1925 PetscCall(DMClone(cdmNew, &cdmCellNew)); 1926 PetscCall(DMSetCellCoordinateDM(rdm, cdmCellNew)); 1927 PetscCall(DMDestroy(&cdmCellNew)); 1928 1929 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew)); 1930 PetscCall(PetscSectionSetNumFields(coordSectionCellNew, 1)); 1931 PetscCall(PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE)); 1932 PetscCall(DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew)); 1933 PetscCall(PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew)); 1934 1935 PetscCall(DMGetCellCoordinateSection(dm, &coordSectionCell)); 1936 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 1937 for (c = cStart; c < cEnd; ++c) { 1938 PetscInt dof; 1939 1940 PetscCall(PetscSectionGetDof(coordSectionCell, c, &dof)); 1941 if (dof) { 1942 DMPolytopeType ct; 1943 DMPolytopeType *rct; 1944 PetscInt *rsize, *rcone, *rornt; 1945 PetscInt dim, cNew, Nct, n, r; 1946 1947 PetscCall(DMPlexGetCellType(dm, c, &ct)); 1948 dim = DMPolytopeTypeGetDim(ct); 1949 PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1950 /* This allows for different cell types */ 1951 for (n = 0; n < Nct; ++n) { 1952 if (dim != DMPolytopeTypeGetDim(rct[n])) continue; 1953 for (r = 0; r < rsize[n]; ++r) { 1954 PetscInt *closure = NULL; 1955 PetscInt clSize, cl, Nv = 0; 1956 1957 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew)); 1958 PetscCall(DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure)); 1959 for (cl = 0; cl < clSize * 2; cl += 2) { 1960 if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv; 1961 } 1962 PetscCall(DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure)); 1963 PetscCall(PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE)); 1964 PetscCall(PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE)); 1965 } 1966 } 1967 } 1968 } 1969 PetscCall(PetscSectionSetUp(coordSectionCellNew)); 1970 PetscCall(DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew)); 1971 } 1972 PetscCall(DMViewFromOptions(dm, NULL, "-coarse_dm_view")); 1973 { 1974 VecType vtype; 1975 PetscInt coordSizeNew, bs; 1976 const char *name; 1977 1978 PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal)); 1979 PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalNew)); 1980 PetscCall(PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew)); 1981 PetscCall(VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE)); 1982 PetscCall(PetscObjectGetName((PetscObject)coordsLocal, &name)); 1983 PetscCall(PetscObjectSetName((PetscObject)coordsLocalNew, name)); 1984 PetscCall(VecGetBlockSize(coordsLocal, &bs)); 1985 PetscCall(VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE)); 1986 PetscCall(VecGetType(coordsLocal, &vtype)); 1987 PetscCall(VecSetType(coordsLocalNew, vtype)); 1988 } 1989 PetscCall(VecGetArrayRead(coordsLocal, &coords)); 1990 PetscCall(VecGetArray(coordsLocalNew, &coordsNew)); 1991 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1992 /* First set coordinates for vertices */ 1993 for (p = pStart; p < pEnd; ++p) { 1994 DMPolytopeType ct; 1995 DMPolytopeType *rct; 1996 PetscInt *rsize, *rcone, *rornt; 1997 PetscInt Nct, n, r; 1998 PetscBool hasVertex = PETSC_FALSE; 1999 2000 PetscCall(DMPlexGetCellType(dm, p, &ct)); 2001 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 2002 for (n = 0; n < Nct; ++n) { 2003 if (rct[n] == DM_POLYTOPE_POINT) { 2004 hasVertex = PETSC_TRUE; 2005 break; 2006 } 2007 } 2008 if (hasVertex) { 2009 const PetscScalar *icoords = NULL; 2010 const PetscScalar *array = NULL; 2011 PetscScalar *pcoords = NULL; 2012 PetscBool isDG; 2013 PetscInt Nc, Nv, v, d; 2014 2015 PetscCall(DMPlexGetCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords)); 2016 2017 icoords = pcoords; 2018 Nv = Nc / dEo; 2019 if (ct != DM_POLYTOPE_POINT) { 2020 if (localizeVertices && maxCell) { 2021 PetscScalar anchor[3]; 2022 2023 for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d]; 2024 for (v = 0; v < Nv; ++v) PetscCall(DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo])); 2025 } 2026 } 2027 for (n = 0; n < Nct; ++n) { 2028 if (rct[n] != DM_POLYTOPE_POINT) continue; 2029 for (r = 0; r < rsize[n]; ++r) { 2030 PetscScalar vcoords[3]; 2031 PetscInt vNew, off; 2032 2033 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew)); 2034 PetscCall(PetscSectionGetOffset(coordSectionNew, vNew, &off)); 2035 PetscCall(DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords)); 2036 PetscCall(DMPlexSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off])); 2037 } 2038 } 2039 PetscCall(DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords)); 2040 } 2041 } 2042 PetscCall(VecRestoreArrayRead(coordsLocal, &coords)); 2043 PetscCall(VecRestoreArray(coordsLocalNew, &coordsNew)); 2044 PetscCall(DMSetCoordinatesLocal(rdm, coordsLocalNew)); 2045 PetscCall(VecDestroy(&coordsLocalNew)); 2046 PetscCall(PetscSectionDestroy(&coordSectionNew)); 2047 /* Then set coordinates for cells by localizing */ 2048 if (!localizeCells) PetscCall(DMLocalizeCoordinates(rdm)); 2049 else { 2050 VecType vtype; 2051 PetscInt coordSizeNew, bs; 2052 const char *name; 2053 2054 PetscCall(DMGetCellCoordinatesLocal(dm, &coordsLocalCell)); 2055 PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew)); 2056 PetscCall(PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew)); 2057 PetscCall(VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE)); 2058 PetscCall(PetscObjectGetName((PetscObject)coordsLocalCell, &name)); 2059 PetscCall(PetscObjectSetName((PetscObject)coordsLocalCellNew, name)); 2060 PetscCall(VecGetBlockSize(coordsLocalCell, &bs)); 2061 PetscCall(VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE)); 2062 PetscCall(VecGetType(coordsLocalCell, &vtype)); 2063 PetscCall(VecSetType(coordsLocalCellNew, vtype)); 2064 PetscCall(VecGetArrayRead(coordsLocalCell, &coords)); 2065 PetscCall(VecGetArray(coordsLocalCellNew, &coordsNew)); 2066 2067 for (p = pStart; p < pEnd; ++p) { 2068 DMPolytopeType ct; 2069 DMPolytopeType *rct; 2070 PetscInt *rsize, *rcone, *rornt; 2071 PetscInt dof = 0, Nct, n, r; 2072 2073 PetscCall(DMPlexGetCellType(dm, p, &ct)); 2074 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 2075 if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(coordSectionCell, p, &dof)); 2076 if (dof) { 2077 const PetscScalar *pcoords; 2078 2079 PetscCall(DMPlexPointLocalRead(cdmCell, p, coords, &pcoords)); 2080 for (n = 0; n < Nct; ++n) { 2081 const PetscInt Nr = rsize[n]; 2082 2083 if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue; 2084 for (r = 0; r < Nr; ++r) { 2085 PetscInt pNew, offNew; 2086 2087 /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that 2088 DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger 2089 cell to the ones it produces. */ 2090 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 2091 PetscCall(PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew)); 2092 PetscCall(DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew])); 2093 } 2094 } 2095 } 2096 } 2097 PetscCall(VecRestoreArrayRead(coordsLocalCell, &coords)); 2098 PetscCall(VecRestoreArray(coordsLocalCellNew, &coordsNew)); 2099 PetscCall(DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew)); 2100 PetscCall(VecDestroy(&coordsLocalCellNew)); 2101 PetscCall(PetscSectionDestroy(&coordSectionCellNew)); 2102 } 2103 PetscFunctionReturn(0); 2104 } 2105 2106 PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm) 2107 { 2108 DM rdm; 2109 DMPlexInterpolatedFlag interp; 2110 PetscInt pStart, pEnd; 2111 2112 PetscFunctionBegin; 2113 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 2114 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 2115 PetscValidPointer(tdm, 3); 2116 PetscCall(DMPlexTransformSetDM(tr, dm)); 2117 2118 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &rdm)); 2119 PetscCall(DMSetType(rdm, DMPLEX)); 2120 PetscCall(DMPlexTransformSetDimensions(tr, dm, rdm)); 2121 /* Calculate number of new points of each depth */ 2122 PetscCall(DMPlexIsInterpolatedCollective(dm, &interp)); 2123 PetscCheck(interp == DMPLEX_INTERPOLATED_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement"); 2124 /* Step 1: Set chart */ 2125 PetscCall(DMPlexTransformGetChart(tr, &pStart, &pEnd)); 2126 PetscCall(DMPlexSetChart(rdm, pStart, pEnd)); 2127 /* Step 2: Set cone/support sizes (automatically stratifies) */ 2128 PetscCall(DMPlexTransformSetConeSizes(tr, rdm)); 2129 /* Step 3: Setup refined DM */ 2130 PetscCall(DMSetUp(rdm)); 2131 /* Step 4: Set cones and supports (automatically symmetrizes) */ 2132 PetscCall(DMPlexTransformSetCones(tr, rdm)); 2133 /* Step 5: Create pointSF */ 2134 PetscCall(DMPlexTransformCreateSF(tr, rdm)); 2135 /* Step 6: Create labels */ 2136 PetscCall(DMPlexTransformCreateLabels(tr, rdm)); 2137 /* Step 7: Set coordinates */ 2138 PetscCall(DMPlexTransformSetCoordinates(tr, rdm)); 2139 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, rdm)); 2140 *tdm = rdm; 2141 PetscFunctionReturn(0); 2142 } 2143 2144 PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm) 2145 { 2146 DMPlexTransform tr; 2147 DM cdm, rcdm; 2148 const char *prefix; 2149 2150 PetscFunctionBegin; 2151 PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr)); 2152 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix)); 2153 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix)); 2154 PetscCall(DMPlexTransformSetDM(tr, dm)); 2155 PetscCall(DMPlexTransformSetFromOptions(tr)); 2156 PetscCall(DMPlexTransformSetActive(tr, adaptLabel)); 2157 PetscCall(DMPlexTransformSetUp(tr)); 2158 PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view")); 2159 PetscCall(DMPlexTransformApply(tr, dm, rdm)); 2160 PetscCall(DMCopyDisc(dm, *rdm)); 2161 PetscCall(DMGetCoordinateDM(dm, &cdm)); 2162 PetscCall(DMGetCoordinateDM(*rdm, &rcdm)); 2163 PetscCall(DMCopyDisc(cdm, rcdm)); 2164 PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm)); 2165 PetscCall(DMCopyDisc(dm, *rdm)); 2166 PetscCall(DMPlexTransformDestroy(&tr)); 2167 ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation; 2168 PetscFunctionReturn(0); 2169 } 2170