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(PETSC_SUCCESS); 45 } 46 47 /*@C 48 DMPlexTransformRegister - Adds a new transform component implementation 49 50 Not Collective 51 52 Input Parameters: 53 + name - The name of a new user-defined creation routine 54 - create_func - The creation routine 55 56 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(PETSC_SUCCESS); 84 } 85 86 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Filter(DMPlexTransform); 87 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Regular(DMPlexTransform); 88 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_ToBox(DMPlexTransform); 89 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Alfeld(DMPlexTransform); 90 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_SBR(DMPlexTransform); 91 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_BL(DMPlexTransform); 92 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_1D(DMPlexTransform); 93 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform); 94 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 165 } 166 167 /*@C 168 DMPlexTransformSetType - Sets the particular implementation for a transform. 169 170 Collective 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; see `DMPlexTransformType` 178 179 Level: intermediate 180 181 .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformGetType()`, `DMPlexTransformCreate()` 182 @*/ 183 PetscErrorCode DMPlexTransformSetType(DMPlexTransform tr, DMPlexTransformType method) 184 { 185 PetscErrorCode (*r)(DMPlexTransform); 186 PetscBool match; 187 188 PetscFunctionBegin; 189 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 190 PetscCall(PetscObjectTypeCompare((PetscObject)tr, method, &match)); 191 if (match) PetscFunctionReturn(PETSC_SUCCESS); 192 193 PetscCall(DMPlexTransformRegisterAll()); 194 PetscCall(PetscFunctionListFind(DMPlexTransformList, method, &r)); 195 PetscCheck(r, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMPlexTransform type: %s", method); 196 197 PetscTryTypeMethod(tr, destroy); 198 PetscCall(PetscMemzero(tr->ops, sizeof(*tr->ops))); 199 PetscCall(PetscObjectChangeTypeName((PetscObject)tr, method)); 200 PetscCall((*r)(tr)); 201 PetscFunctionReturn(PETSC_SUCCESS); 202 } 203 204 /*@C 205 DMPlexTransformGetType - Gets the type name (as a string) from the transform. 206 207 Not Collective 208 209 Input Parameter: 210 . tr - The `DMPlexTransform` 211 212 Output Parameter: 213 . type - The `DMPlexTransformType` name 214 215 Level: intermediate 216 217 .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `DMPlexTransformSetType()`, `DMPlexTransformCreate()` 218 @*/ 219 PetscErrorCode DMPlexTransformGetType(DMPlexTransform tr, DMPlexTransformType *type) 220 { 221 PetscFunctionBegin; 222 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 223 PetscValidPointer(type, 2); 224 PetscCall(DMPlexTransformRegisterAll()); 225 *type = ((PetscObject)tr)->type_name; 226 PetscFunctionReturn(PETSC_SUCCESS); 227 } 228 229 static PetscErrorCode DMPlexTransformView_Ascii(DMPlexTransform tr, PetscViewer v) 230 { 231 PetscViewerFormat format; 232 233 PetscFunctionBegin; 234 PetscCall(PetscViewerGetFormat(v, &format)); 235 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 236 const PetscInt *trTypes = NULL; 237 IS trIS; 238 PetscInt cols = 8; 239 PetscInt Nrt = 8, f, g; 240 241 if (tr->trType) PetscCall(DMLabelView(tr->trType, v)); 242 PetscCall(PetscViewerASCIIPrintf(v, "Source Starts\n")); 243 for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g])); 244 PetscCall(PetscViewerASCIIPrintf(v, "\n")); 245 for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStart[f])); 246 PetscCall(PetscViewerASCIIPrintf(v, "\n")); 247 PetscCall(PetscViewerASCIIPrintf(v, "Target Starts\n")); 248 for (g = 0; g <= cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g])); 249 PetscCall(PetscViewerASCIIPrintf(v, "\n")); 250 for (f = 0; f <= cols; ++f) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->ctStartNew[f])); 251 PetscCall(PetscViewerASCIIPrintf(v, "\n")); 252 253 if (tr->trType) { 254 PetscCall(DMLabelGetNumValues(tr->trType, &Nrt)); 255 PetscCall(DMLabelGetValueIS(tr->trType, &trIS)); 256 PetscCall(ISGetIndices(trIS, &trTypes)); 257 } 258 PetscCall(PetscViewerASCIIPrintf(v, "Offsets\n")); 259 PetscCall(PetscViewerASCIIPrintf(v, " ")); 260 for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14s", DMPolytopeTypes[g])); 261 PetscCall(PetscViewerASCIIPrintf(v, "\n")); 262 for (f = 0; f < Nrt; ++f) { 263 PetscCall(PetscViewerASCIIPrintf(v, "%2" PetscInt_FMT " |", trTypes ? trTypes[f] : f)); 264 for (g = 0; g < cols; ++g) PetscCall(PetscViewerASCIIPrintf(v, " %14" PetscInt_FMT, tr->offset[f * DM_NUM_POLYTOPES + g])); 265 PetscCall(PetscViewerASCIIPrintf(v, " |\n")); 266 } 267 if (tr->trType) { 268 PetscCall(ISRestoreIndices(trIS, &trTypes)); 269 PetscCall(ISDestroy(&trIS)); 270 } 271 } 272 PetscFunctionReturn(PETSC_SUCCESS); 273 } 274 275 /*@C 276 DMPlexTransformView - Views a `DMPlexTransform` 277 278 Collective 279 280 Input Parameters: 281 + tr - the `DMPlexTransform` object to view 282 - v - the viewer 283 284 Level: beginner 285 286 .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformType`, `PetscViewer`, `DMPlexTransformDestroy()`, `DMPlexTransformCreate()` 287 @*/ 288 PetscErrorCode DMPlexTransformView(DMPlexTransform tr, PetscViewer v) 289 { 290 PetscBool isascii; 291 292 PetscFunctionBegin; 293 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 294 if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tr), &v)); 295 PetscValidHeaderSpecific(v, PETSC_VIEWER_CLASSID, 2); 296 PetscCheckSameComm(tr, 1, v, 2); 297 PetscCall(PetscViewerCheckWritable(v)); 298 PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tr, v)); 299 PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii)); 300 if (isascii) PetscCall(DMPlexTransformView_Ascii(tr, v)); 301 PetscTryTypeMethod(tr, view, v); 302 PetscFunctionReturn(PETSC_SUCCESS); 303 } 304 305 /*@ 306 DMPlexTransformSetFromOptions - Sets parameters in a transform from values in the options database 307 308 Collective 309 310 Input Parameter: 311 . tr - the `DMPlexTransform` object to set options for 312 313 Options Database Key: 314 . -dm_plex_transform_type - Set the transform type, e.g. refine_regular 315 316 Level: intermediate 317 318 .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()` 319 @*/ 320 PetscErrorCode DMPlexTransformSetFromOptions(DMPlexTransform tr) 321 { 322 char typeName[1024]; 323 const char *defName = DMPLEXREFINEREGULAR; 324 PetscBool flg; 325 326 PetscFunctionBegin; 327 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 328 PetscObjectOptionsBegin((PetscObject)tr); 329 PetscCall(PetscOptionsFList("-dm_plex_transform_type", "DMPlexTransform", "DMPlexTransformSetType", DMPlexTransformList, defName, typeName, 1024, &flg)); 330 if (flg) PetscCall(DMPlexTransformSetType(tr, typeName)); 331 else if (!((PetscObject)tr)->type_name) PetscCall(DMPlexTransformSetType(tr, defName)); 332 PetscTryTypeMethod(tr, setfromoptions, PetscOptionsObject); 333 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 334 PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)tr, PetscOptionsObject)); 335 PetscOptionsEnd(); 336 PetscFunctionReturn(PETSC_SUCCESS); 337 } 338 339 /*@C 340 DMPlexTransformDestroy - Destroys a `DMPlexTransform` 341 342 Collective 343 344 Input Parameter: 345 . tr - the transform object to destroy 346 347 Level: beginner 348 349 .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformView()`, `DMPlexTransformCreate()` 350 @*/ 351 PetscErrorCode DMPlexTransformDestroy(DMPlexTransform *tr) 352 { 353 PetscInt c; 354 355 PetscFunctionBegin; 356 if (!*tr) PetscFunctionReturn(PETSC_SUCCESS); 357 PetscValidHeaderSpecific((*tr), DMPLEXTRANSFORM_CLASSID, 1); 358 if (--((PetscObject)(*tr))->refct > 0) { 359 *tr = NULL; 360 PetscFunctionReturn(PETSC_SUCCESS); 361 } 362 363 if ((*tr)->ops->destroy) PetscCall((*(*tr)->ops->destroy)(*tr)); 364 PetscCall(DMDestroy(&(*tr)->dm)); 365 PetscCall(DMLabelDestroy(&(*tr)->active)); 366 PetscCall(DMLabelDestroy(&(*tr)->trType)); 367 PetscCall(PetscFree2((*tr)->ctOrderOld, (*tr)->ctOrderInvOld)); 368 PetscCall(PetscFree2((*tr)->ctOrderNew, (*tr)->ctOrderInvNew)); 369 PetscCall(PetscFree2((*tr)->ctStart, (*tr)->ctStartNew)); 370 PetscCall(PetscFree((*tr)->offset)); 371 PetscCall(PetscFree2((*tr)->depthStart, (*tr)->depthEnd)); 372 for (c = 0; c < DM_NUM_POLYTOPES; ++c) { 373 PetscCall(PetscFEDestroy(&(*tr)->coordFE[c])); 374 PetscCall(PetscFEGeomDestroy(&(*tr)->refGeom[c])); 375 } 376 if ((*tr)->trVerts) { 377 for (c = 0; c < DM_NUM_POLYTOPES; ++c) { 378 DMPolytopeType *rct; 379 PetscInt *rsize, *rcone, *rornt, Nct, n, r; 380 381 if (DMPolytopeTypeGetDim((DMPolytopeType)c) > 0) { 382 PetscCall(DMPlexTransformCellTransform((*tr), (DMPolytopeType)c, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 383 for (n = 0; n < Nct; ++n) { 384 if (rct[n] == DM_POLYTOPE_POINT) continue; 385 for (r = 0; r < rsize[n]; ++r) PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]][r])); 386 PetscCall(PetscFree((*tr)->trSubVerts[c][rct[n]])); 387 } 388 } 389 PetscCall(PetscFree((*tr)->trSubVerts[c])); 390 PetscCall(PetscFree((*tr)->trVerts[c])); 391 } 392 } 393 PetscCall(PetscFree3((*tr)->trNv, (*tr)->trVerts, (*tr)->trSubVerts)); 394 PetscCall(PetscFree2((*tr)->coordFE, (*tr)->refGeom)); 395 /* We do not destroy (*dm)->data here so that we can reference count backend objects */ 396 PetscCall(PetscHeaderDestroy(tr)); 397 PetscFunctionReturn(PETSC_SUCCESS); 398 } 399 400 static PetscErrorCode DMPlexTransformCreateOffset_Internal(DMPlexTransform tr, PetscInt ctOrderOld[], PetscInt ctStart[], PetscInt **offset) 401 { 402 DMLabel trType = tr->trType; 403 PetscInt c, cN, *off; 404 405 PetscFunctionBegin; 406 if (trType) { 407 DM dm; 408 IS rtIS; 409 const PetscInt *reftypes; 410 PetscInt Nrt, r; 411 412 PetscCall(DMPlexTransformGetDM(tr, &dm)); 413 PetscCall(DMLabelGetNumValues(trType, &Nrt)); 414 PetscCall(DMLabelGetValueIS(trType, &rtIS)); 415 PetscCall(ISGetIndices(rtIS, &reftypes)); 416 PetscCall(PetscCalloc1(Nrt * DM_NUM_POLYTOPES, &off)); 417 for (r = 0; r < Nrt; ++r) { 418 const PetscInt rt = reftypes[r]; 419 IS rtIS; 420 const PetscInt *points; 421 DMPolytopeType ct; 422 PetscInt np, p; 423 424 PetscCall(DMLabelGetStratumIS(trType, rt, &rtIS)); 425 PetscCall(ISGetLocalSize(rtIS, &np)); 426 PetscCall(ISGetIndices(rtIS, &points)); 427 if (!np) continue; 428 p = points[0]; 429 PetscCall(ISRestoreIndices(rtIS, &points)); 430 PetscCall(ISDestroy(&rtIS)); 431 PetscCall(DMPlexGetCellType(dm, p, &ct)); 432 for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) { 433 const DMPolytopeType ctNew = (DMPolytopeType)cN; 434 DMPolytopeType *rct; 435 PetscInt *rsize, *cone, *ornt; 436 PetscInt Nct, n, s; 437 438 if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) { 439 off[r * DM_NUM_POLYTOPES + ctNew] = -1; 440 break; 441 } 442 off[r * DM_NUM_POLYTOPES + ctNew] = 0; 443 for (s = 0; s <= r; ++s) { 444 const PetscInt st = reftypes[s]; 445 DMPolytopeType sct; 446 PetscInt q, qrt; 447 448 PetscCall(DMLabelGetStratumIS(trType, st, &rtIS)); 449 PetscCall(ISGetLocalSize(rtIS, &np)); 450 PetscCall(ISGetIndices(rtIS, &points)); 451 if (!np) continue; 452 q = points[0]; 453 PetscCall(ISRestoreIndices(rtIS, &points)); 454 PetscCall(ISDestroy(&rtIS)); 455 PetscCall(DMPlexGetCellType(dm, q, &sct)); 456 PetscCall(DMPlexTransformCellTransform(tr, sct, q, &qrt, &Nct, &rct, &rsize, &cone, &ornt)); 457 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); 458 if (st == rt) { 459 for (n = 0; n < Nct; ++n) 460 if (rct[n] == ctNew) break; 461 if (n == Nct) off[r * DM_NUM_POLYTOPES + ctNew] = -1; 462 break; 463 } 464 for (n = 0; n < Nct; ++n) { 465 if (rct[n] == ctNew) { 466 PetscInt sn; 467 468 PetscCall(DMLabelGetStratumSize(trType, st, &sn)); 469 off[r * DM_NUM_POLYTOPES + ctNew] += sn * rsize[n]; 470 } 471 } 472 } 473 } 474 } 475 PetscCall(ISRestoreIndices(rtIS, &reftypes)); 476 PetscCall(ISDestroy(&rtIS)); 477 } else { 478 PetscCall(PetscCalloc1(DM_NUM_POLYTOPES * DM_NUM_POLYTOPES, &off)); 479 for (c = DM_POLYTOPE_POINT; c < DM_NUM_POLYTOPES; ++c) { 480 const DMPolytopeType ct = (DMPolytopeType)c; 481 for (cN = DM_POLYTOPE_POINT; cN < DM_NUM_POLYTOPES; ++cN) { 482 const DMPolytopeType ctNew = (DMPolytopeType)cN; 483 DMPolytopeType *rct; 484 PetscInt *rsize, *cone, *ornt; 485 PetscInt Nct, n, i; 486 487 if (DMPolytopeTypeGetDim(ct) < 0 || DMPolytopeTypeGetDim(ctNew) < 0) { 488 off[ct * DM_NUM_POLYTOPES + ctNew] = -1; 489 continue; 490 } 491 off[ct * DM_NUM_POLYTOPES + ctNew] = 0; 492 for (i = DM_POLYTOPE_POINT; i < DM_NUM_POLYTOPES; ++i) { 493 const DMPolytopeType ict = (DMPolytopeType)ctOrderOld[i]; 494 const DMPolytopeType ictn = (DMPolytopeType)ctOrderOld[i + 1]; 495 496 PetscCall(DMPlexTransformCellTransform(tr, ict, PETSC_DETERMINE, NULL, &Nct, &rct, &rsize, &cone, &ornt)); 497 if (ict == ct) { 498 for (n = 0; n < Nct; ++n) 499 if (rct[n] == ctNew) break; 500 if (n == Nct) off[ct * DM_NUM_POLYTOPES + ctNew] = -1; 501 break; 502 } 503 for (n = 0; n < Nct; ++n) 504 if (rct[n] == ctNew) off[ct * DM_NUM_POLYTOPES + ctNew] += (ctStart[ictn] - ctStart[ict]) * rsize[n]; 505 } 506 } 507 } 508 } 509 *offset = off; 510 PetscFunctionReturn(PETSC_SUCCESS); 511 } 512 513 PetscErrorCode DMPlexTransformSetUp(DMPlexTransform tr) 514 { 515 DM dm; 516 DMPolytopeType ctCell; 517 PetscInt pStart, pEnd, p, c, celldim = 0; 518 519 PetscFunctionBegin; 520 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 521 if (tr->setupcalled) PetscFunctionReturn(PETSC_SUCCESS); 522 PetscTryTypeMethod(tr, setup); 523 PetscCall(DMPlexTransformGetDM(tr, &dm)); 524 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 525 if (pEnd > pStart) { 526 PetscCall(DMPlexGetCellType(dm, 0, &ctCell)); 527 } else { 528 PetscInt dim; 529 530 PetscCall(DMGetDimension(dm, &dim)); 531 switch (dim) { 532 case 0: 533 ctCell = DM_POLYTOPE_POINT; 534 break; 535 case 1: 536 ctCell = DM_POLYTOPE_SEGMENT; 537 break; 538 case 2: 539 ctCell = DM_POLYTOPE_TRIANGLE; 540 break; 541 case 3: 542 ctCell = DM_POLYTOPE_TETRAHEDRON; 543 break; 544 default: 545 ctCell = DM_POLYTOPE_UNKNOWN; 546 } 547 } 548 PetscCall(DMPlexCreateCellTypeOrder_Internal(DMPolytopeTypeGetDim(ctCell), &tr->ctOrderOld, &tr->ctOrderInvOld)); 549 for (p = pStart; p < pEnd; ++p) { 550 DMPolytopeType ct; 551 DMPolytopeType *rct; 552 PetscInt *rsize, *cone, *ornt; 553 PetscInt Nct, n; 554 555 PetscCall(DMPlexGetCellType(dm, p, &ct)); 556 PetscCheck(ct != DM_POLYTOPE_UNKNOWN, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %" PetscInt_FMT, p); 557 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt)); 558 for (n = 0; n < Nct; ++n) celldim = PetscMax(celldim, DMPolytopeTypeGetDim(rct[n])); 559 } 560 PetscCall(DMPlexCreateCellTypeOrder_Internal(celldim, &tr->ctOrderNew, &tr->ctOrderInvNew)); 561 /* Construct sizes and offsets for each cell type */ 562 if (!tr->ctStart) { 563 PetscInt *ctS, *ctSN, *ctC, *ctCN; 564 565 PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctS, DM_NUM_POLYTOPES + 1, &ctSN)); 566 PetscCall(PetscCalloc2(DM_NUM_POLYTOPES + 1, &ctC, DM_NUM_POLYTOPES + 1, &ctCN)); 567 for (p = pStart; p < pEnd; ++p) { 568 DMPolytopeType ct; 569 DMPolytopeType *rct; 570 PetscInt *rsize, *cone, *ornt; 571 PetscInt Nct, n; 572 573 PetscCall(DMPlexGetCellType(dm, p, &ct)); 574 PetscCheck(ct != DM_POLYTOPE_UNKNOWN, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No cell type for point %" PetscInt_FMT, p); 575 ++ctC[ct]; 576 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &cone, &ornt)); 577 for (n = 0; n < Nct; ++n) ctCN[rct[n]] += rsize[n]; 578 } 579 for (c = 0; c < DM_NUM_POLYTOPES; ++c) { 580 const PetscInt cto = tr->ctOrderOld[c]; 581 const PetscInt cton = tr->ctOrderOld[c + 1]; 582 const PetscInt ctn = tr->ctOrderNew[c]; 583 const PetscInt ctnn = tr->ctOrderNew[c + 1]; 584 585 ctS[cton] = ctS[cto] + ctC[cto]; 586 ctSN[ctnn] = ctSN[ctn] + ctCN[ctn]; 587 } 588 PetscCall(PetscFree2(ctC, ctCN)); 589 tr->ctStart = ctS; 590 tr->ctStartNew = ctSN; 591 } 592 PetscCall(DMPlexTransformCreateOffset_Internal(tr, tr->ctOrderOld, tr->ctStart, &tr->offset)); 593 // Compute depth information 594 tr->depth = -1; 595 for (c = 0; c < DM_NUM_POLYTOPES; ++c) 596 if (tr->ctStartNew[tr->ctOrderNew[c + 1]] > tr->ctStartNew[tr->ctOrderNew[c]]) tr->depth = PetscMax(tr->depth, DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c])); 597 PetscCall(PetscMalloc2(tr->depth + 1, &tr->depthStart, tr->depth + 1, &tr->depthEnd)); 598 for (PetscInt d = 0; d <= tr->depth; ++d) { 599 tr->depthStart[d] = PETSC_MAX_INT; 600 tr->depthEnd[d] = -1; 601 } 602 for (c = 0; c < DM_NUM_POLYTOPES; ++c) { 603 const PetscInt dep = DMPolytopeTypeGetDim((DMPolytopeType)tr->ctOrderNew[c]); 604 605 if (tr->ctStartNew[tr->ctOrderNew[c + 1]] <= tr->ctStartNew[tr->ctOrderNew[c]]) continue; 606 tr->depthStart[dep] = PetscMin(tr->depthStart[dep], tr->ctStartNew[tr->ctOrderNew[c]]); 607 tr->depthEnd[dep] = PetscMax(tr->depthEnd[dep], tr->ctStartNew[tr->ctOrderNew[c + 1]]); 608 } 609 tr->setupcalled = PETSC_TRUE; 610 PetscFunctionReturn(PETSC_SUCCESS); 611 } 612 613 PetscErrorCode DMPlexTransformGetDM(DMPlexTransform tr, DM *dm) 614 { 615 PetscFunctionBegin; 616 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 617 PetscValidPointer(dm, 2); 618 *dm = tr->dm; 619 PetscFunctionReturn(PETSC_SUCCESS); 620 } 621 622 PetscErrorCode DMPlexTransformSetDM(DMPlexTransform tr, DM dm) 623 { 624 PetscFunctionBegin; 625 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 626 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 627 PetscCall(PetscObjectReference((PetscObject)dm)); 628 PetscCall(DMDestroy(&tr->dm)); 629 tr->dm = dm; 630 PetscFunctionReturn(PETSC_SUCCESS); 631 } 632 633 PetscErrorCode DMPlexTransformGetActive(DMPlexTransform tr, DMLabel *active) 634 { 635 PetscFunctionBegin; 636 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 637 PetscValidPointer(active, 2); 638 *active = tr->active; 639 PetscFunctionReturn(PETSC_SUCCESS); 640 } 641 642 PetscErrorCode DMPlexTransformSetActive(DMPlexTransform tr, DMLabel active) 643 { 644 PetscFunctionBegin; 645 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 646 if (active) PetscValidHeaderSpecific(active, DMLABEL_CLASSID, 2); 647 PetscCall(PetscObjectReference((PetscObject)active)); 648 PetscCall(DMLabelDestroy(&tr->active)); 649 tr->active = active; 650 PetscFunctionReturn(PETSC_SUCCESS); 651 } 652 653 static PetscErrorCode DMPlexTransformGetCoordinateFE(DMPlexTransform tr, DMPolytopeType ct, PetscFE *fe) 654 { 655 PetscFunctionBegin; 656 if (!tr->coordFE[ct]) { 657 PetscInt dim, cdim; 658 659 dim = DMPolytopeTypeGetDim(ct); 660 PetscCall(DMGetCoordinateDim(tr->dm, &cdim)); 661 PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, cdim, ct, 1, PETSC_DETERMINE, &tr->coordFE[ct])); 662 { 663 PetscDualSpace dsp; 664 PetscQuadrature quad; 665 DM K; 666 PetscFEGeom *cg; 667 PetscScalar *Xq; 668 PetscReal *xq, *wq; 669 PetscInt Nq, q; 670 671 PetscCall(DMPlexTransformGetCellVertices(tr, ct, &Nq, &Xq)); 672 PetscCall(PetscMalloc1(Nq * cdim, &xq)); 673 for (q = 0; q < Nq * cdim; ++q) xq[q] = PetscRealPart(Xq[q]); 674 PetscCall(PetscMalloc1(Nq, &wq)); 675 for (q = 0; q < Nq; ++q) wq[q] = 1.0; 676 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &quad)); 677 PetscCall(PetscQuadratureSetData(quad, dim, 1, Nq, xq, wq)); 678 PetscCall(PetscFESetQuadrature(tr->coordFE[ct], quad)); 679 680 PetscCall(PetscFEGetDualSpace(tr->coordFE[ct], &dsp)); 681 PetscCall(PetscDualSpaceGetDM(dsp, &K)); 682 PetscCall(PetscFEGeomCreate(quad, 1, cdim, PETSC_FALSE, &tr->refGeom[ct])); 683 cg = tr->refGeom[ct]; 684 PetscCall(DMPlexComputeCellGeometryFEM(K, 0, NULL, cg->v, cg->J, cg->invJ, cg->detJ)); 685 PetscCall(PetscQuadratureDestroy(&quad)); 686 } 687 } 688 *fe = tr->coordFE[ct]; 689 PetscFunctionReturn(PETSC_SUCCESS); 690 } 691 692 PetscErrorCode DMPlexTransformSetDimensions_Internal(DMPlexTransform tr, DM dm, DM tdm) 693 { 694 PetscInt dim, cdim; 695 696 PetscFunctionBegin; 697 PetscCall(DMGetDimension(dm, &dim)); 698 PetscCall(DMSetDimension(tdm, dim)); 699 PetscCall(DMGetCoordinateDim(dm, &cdim)); 700 PetscCall(DMSetCoordinateDim(tdm, cdim)); 701 PetscFunctionReturn(PETSC_SUCCESS); 702 } 703 704 /*@ 705 DMPlexTransformSetDimensions - Set the dimensions for the transformed `DM` 706 707 Input Parameters: 708 + tr - The `DMPlexTransform` object 709 - dm - The original `DM` 710 711 Output Parameter: 712 . tdm - The transformed `DM` 713 714 Level: advanced 715 716 .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` 717 @*/ 718 PetscErrorCode DMPlexTransformSetDimensions(DMPlexTransform tr, DM dm, DM tdm) 719 { 720 PetscFunctionBegin; 721 PetscUseTypeMethod(tr, setdimensions, dm, tdm); 722 PetscFunctionReturn(PETSC_SUCCESS); 723 } 724 725 PetscErrorCode DMPlexTransformGetChart(DMPlexTransform tr, PetscInt *pStart, PetscInt *pEnd) 726 { 727 PetscFunctionBegin; 728 if (pStart) *pStart = 0; 729 if (pEnd) *pEnd = tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]]; 730 PetscFunctionReturn(PETSC_SUCCESS); 731 } 732 733 PetscErrorCode DMPlexTransformGetCellType(DMPlexTransform tr, PetscInt cell, DMPolytopeType *celltype) 734 { 735 PetscInt ctNew; 736 737 PetscFunctionBegin; 738 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 739 PetscValidPointer(celltype, 3); 740 /* TODO Can do bisection since everything is sorted */ 741 for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) { 742 PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]]; 743 744 if (cell >= ctSN && cell < ctEN) break; 745 } 746 PetscCheck(ctNew < DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " cannot be located in the transformed mesh", cell); 747 *celltype = (DMPolytopeType)ctNew; 748 PetscFunctionReturn(PETSC_SUCCESS); 749 } 750 751 PetscErrorCode DMPlexTransformGetDepth(DMPlexTransform tr, PetscInt *depth) 752 { 753 PetscFunctionBegin; 754 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 755 *depth = tr->depth; 756 PetscFunctionReturn(PETSC_SUCCESS); 757 } 758 759 PetscErrorCode DMPlexTransformGetDepthStratum(DMPlexTransform tr, PetscInt depth, PetscInt *start, PetscInt *end) 760 { 761 PetscFunctionBegin; 762 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 763 if (start) *start = tr->depthStart[depth]; 764 if (end) *end = tr->depthEnd[depth]; 765 PetscFunctionReturn(PETSC_SUCCESS); 766 } 767 768 /*@ 769 DMPlexTransformGetTargetPoint - Get the number of a point in the transformed mesh based on information from the original mesh. 770 771 Not Collective 772 773 Input Parameters: 774 + tr - The `DMPlexTransform` 775 . ct - The type of the original point which produces the new point 776 . ctNew - The type of the new point 777 . p - The original point which produces the new point 778 - r - The replica number of the new point, meaning it is the rth point of type `ctNew` produced from `p` 779 780 Output Parameter: 781 . pNew - The new point number 782 783 Level: developer 784 785 .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSourcePoint()`, `DMPlexTransformCellTransform()` 786 @*/ 787 PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew) 788 { 789 DMPolytopeType *rct; 790 PetscInt *rsize, *cone, *ornt; 791 PetscInt rt, Nct, n, off, rp; 792 DMLabel trType = tr->trType; 793 PetscInt ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ct] + 1]]; 794 PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew] + 1]]; 795 PetscInt newp = ctSN, cind; 796 797 PetscFunctionBeginHot; 798 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); 799 PetscCall(DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt)); 800 if (trType) { 801 PetscCall(DMLabelGetValueIndex(trType, rt, &cind)); 802 PetscCall(DMLabelGetStratumPointIndex(trType, rt, p, &rp)); 803 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); 804 } else { 805 cind = ct; 806 rp = p - ctS; 807 } 808 off = tr->offset[cind * DM_NUM_POLYTOPES + ctNew]; 809 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); 810 newp += off; 811 for (n = 0; n < Nct; ++n) { 812 if (rct[n] == ctNew) { 813 if (rsize[n] && r >= rsize[n]) 814 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]); 815 newp += rp * rsize[n] + r; 816 break; 817 } 818 } 819 820 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); 821 *pNew = newp; 822 PetscFunctionReturn(PETSC_SUCCESS); 823 } 824 825 /*@ 826 DMPlexTransformGetSourcePoint - Get the number of a point in the original mesh based on information from the transformed mesh. 827 828 Not Collective 829 830 Input Parameters: 831 + tr - The `DMPlexTransform` 832 - pNew - The new point number 833 834 Output Parameters: 835 + ct - The type of the original point which produces the new point 836 . ctNew - The type of the new point 837 . p - The original point which produces the new point 838 - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p 839 840 Level: developer 841 842 .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetTargetPoint()`, `DMPlexTransformCellTransform()` 843 @*/ 844 PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r) 845 { 846 DMLabel trType = tr->trType; 847 DMPolytopeType *rct, ctN; 848 PetscInt *rsize, *cone, *ornt; 849 PetscInt rt = -1, rtTmp, Nct, n, rp = 0, rO = 0, pO; 850 PetscInt offset = -1, ctS, ctE, ctO = 0, ctTmp, rtS; 851 852 PetscFunctionBegin; 853 PetscCall(DMPlexTransformGetCellType(tr, pNew, &ctN)); 854 if (trType) { 855 DM dm; 856 IS rtIS; 857 const PetscInt *reftypes; 858 PetscInt Nrt, r, rtStart; 859 860 PetscCall(DMPlexTransformGetDM(tr, &dm)); 861 PetscCall(DMLabelGetNumValues(trType, &Nrt)); 862 PetscCall(DMLabelGetValueIS(trType, &rtIS)); 863 PetscCall(ISGetIndices(rtIS, &reftypes)); 864 for (r = 0; r < Nrt; ++r) { 865 const PetscInt off = tr->offset[r * DM_NUM_POLYTOPES + ctN]; 866 867 if (tr->ctStartNew[ctN] + off > pNew) continue; 868 /* Check that any of this refinement type exist */ 869 /* TODO Actually keep track of the number produced here instead */ 870 if (off > offset) { 871 rt = reftypes[r]; 872 offset = off; 873 } 874 } 875 PetscCall(ISRestoreIndices(rtIS, &reftypes)); 876 PetscCall(ISDestroy(&rtIS)); 877 PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew); 878 /* TODO Map refinement types to cell types */ 879 PetscCall(DMLabelGetStratumBounds(trType, rt, &rtStart, NULL)); 880 PetscCheck(rtStart >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Refinement type %" PetscInt_FMT " has no source points", rt); 881 for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) { 882 PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]]; 883 884 if ((rtStart >= ctS) && (rtStart < ctE)) break; 885 } 886 PetscCheck(ctO != DM_NUM_POLYTOPES, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine a cell type for refinement type %" PetscInt_FMT, rt); 887 } else { 888 for (ctTmp = 0; ctTmp < DM_NUM_POLYTOPES; ++ctTmp) { 889 const PetscInt off = tr->offset[ctTmp * DM_NUM_POLYTOPES + ctN]; 890 891 if (tr->ctStartNew[ctN] + off > pNew) continue; 892 if (tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctTmp] + 1]] <= tr->ctStart[ctTmp]) continue; 893 /* TODO Actually keep track of the number produced here instead */ 894 if (off > offset) { 895 ctO = ctTmp; 896 offset = off; 897 } 898 } 899 rt = -1; 900 PetscCheck(offset >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew); 901 } 902 ctS = tr->ctStart[ctO]; 903 ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO] + 1]]; 904 if (trType) { 905 for (rtS = ctS; rtS < ctE; ++rtS) { 906 PetscInt val; 907 PetscCall(DMLabelGetValue(trType, rtS, &val)); 908 if (val == rt) break; 909 } 910 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); 911 } else rtS = ctS; 912 PetscCall(DMPlexTransformCellTransform(tr, (DMPolytopeType)ctO, rtS, &rtTmp, &Nct, &rct, &rsize, &cone, &ornt)); 913 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); 914 for (n = 0; n < Nct; ++n) { 915 if (rct[n] == ctN) { 916 PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset, val, c; 917 918 if (trType) { 919 for (c = ctS; c < ctE; ++c) { 920 PetscCall(DMLabelGetValue(trType, c, &val)); 921 if (val == rt) { 922 if (tmp < rsize[n]) break; 923 tmp -= rsize[n]; 924 } 925 } 926 PetscCheck(c < ctE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Parent point for target point %" PetscInt_FMT " could be not found", pNew); 927 rp = c - ctS; 928 rO = tmp; 929 } else { 930 // This assumes that all points of type ctO transform the same way 931 rp = tmp / rsize[n]; 932 rO = tmp % rsize[n]; 933 } 934 break; 935 } 936 } 937 PetscCheck(n != Nct, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number for target point %" PetscInt_FMT " could be not found", pNew); 938 pO = rp + ctS; 939 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); 940 if (ct) *ct = (DMPolytopeType)ctO; 941 if (ctNew) *ctNew = (DMPolytopeType)ctN; 942 if (p) *p = pO; 943 if (r) *r = rO; 944 PetscFunctionReturn(PETSC_SUCCESS); 945 } 946 947 /*@ 948 DMPlexTransformCellTransform - Describes the transform of a given source cell into a set of other target cells. These produced cells become the new mesh. 949 950 Input Parameters: 951 + tr - The `DMPlexTransform` object 952 . source - The source cell type 953 - p - The source point, which can also determine the refine type 954 955 Output Parameters: 956 + rt - The refine type for this point 957 . Nt - The number of types produced by this point 958 . target - An array of length `Nt` giving the types produced 959 . size - An array of length `Nt` giving the number of cells of each type produced 960 . cone - An array of length `Nt`*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point 961 - ornt - An array of length `Nt`*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point 962 963 Level: advanced 964 965 Notes: 966 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 .vb 970 the number of cones to be taken, or 0 for the current cell 971 the cell cone point number at each level from which it is subdivided 972 the replica number r of the subdivision. 973 .ve 974 The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is 975 .vb 976 Nt = 2 977 target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT} 978 size = {1, 2} 979 cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0} 980 ornt = { 0, 0, 0, 0} 981 .ve 982 983 .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformApply()`, `DMPlexTransformCreate()` 984 @*/ 985 PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 986 { 987 PetscFunctionBegin; 988 PetscUseTypeMethod(tr, celltransform, source, p, rt, Nt, target, size, cone, ornt); 989 PetscFunctionReturn(PETSC_SUCCESS); 990 } 991 992 PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 993 { 994 PetscFunctionBegin; 995 *rnew = r; 996 *onew = DMPolytopeTypeComposeOrientation(tct, o, so); 997 PetscFunctionReturn(PETSC_SUCCESS); 998 } 999 1000 /* Returns the same thing */ 1001 PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 1002 { 1003 static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT}; 1004 static PetscInt vertexS[] = {1}; 1005 static PetscInt vertexC[] = {0}; 1006 static PetscInt vertexO[] = {0}; 1007 static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT}; 1008 static PetscInt edgeS[] = {1}; 1009 static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 1010 static PetscInt edgeO[] = {0, 0}; 1011 static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR}; 1012 static PetscInt tedgeS[] = {1}; 1013 static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 1014 static PetscInt tedgeO[] = {0, 0}; 1015 static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE}; 1016 static PetscInt triS[] = {1}; 1017 static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0}; 1018 static PetscInt triO[] = {0, 0, 0}; 1019 static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL}; 1020 static PetscInt quadS[] = {1}; 1021 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}; 1022 static PetscInt quadO[] = {0, 0, 0, 0}; 1023 static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR}; 1024 static PetscInt tquadS[] = {1}; 1025 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}; 1026 static PetscInt tquadO[] = {0, 0, 0, 0}; 1027 static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON}; 1028 static PetscInt tetS[] = {1}; 1029 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}; 1030 static PetscInt tetO[] = {0, 0, 0, 0}; 1031 static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON}; 1032 static PetscInt hexS[] = {1}; 1033 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}; 1034 static PetscInt hexO[] = {0, 0, 0, 0, 0, 0}; 1035 static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM}; 1036 static PetscInt tripS[] = {1}; 1037 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}; 1038 static PetscInt tripO[] = {0, 0, 0, 0, 0}; 1039 static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR}; 1040 static PetscInt ttripS[] = {1}; 1041 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}; 1042 static PetscInt ttripO[] = {0, 0, 0, 0, 0}; 1043 static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR}; 1044 static PetscInt tquadpS[] = {1}; 1045 static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, 1046 DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 5, 0}; 1047 static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0}; 1048 static DMPolytopeType pyrT[] = {DM_POLYTOPE_PYRAMID}; 1049 static PetscInt pyrS[] = {1}; 1050 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}; 1051 static PetscInt pyrO[] = {0, 0, 0, 0, 0}; 1052 1053 PetscFunctionBegin; 1054 if (rt) *rt = 0; 1055 switch (source) { 1056 case DM_POLYTOPE_POINT: 1057 *Nt = 1; 1058 *target = vertexT; 1059 *size = vertexS; 1060 *cone = vertexC; 1061 *ornt = vertexO; 1062 break; 1063 case DM_POLYTOPE_SEGMENT: 1064 *Nt = 1; 1065 *target = edgeT; 1066 *size = edgeS; 1067 *cone = edgeC; 1068 *ornt = edgeO; 1069 break; 1070 case DM_POLYTOPE_POINT_PRISM_TENSOR: 1071 *Nt = 1; 1072 *target = tedgeT; 1073 *size = tedgeS; 1074 *cone = tedgeC; 1075 *ornt = tedgeO; 1076 break; 1077 case DM_POLYTOPE_TRIANGLE: 1078 *Nt = 1; 1079 *target = triT; 1080 *size = triS; 1081 *cone = triC; 1082 *ornt = triO; 1083 break; 1084 case DM_POLYTOPE_QUADRILATERAL: 1085 *Nt = 1; 1086 *target = quadT; 1087 *size = quadS; 1088 *cone = quadC; 1089 *ornt = quadO; 1090 break; 1091 case DM_POLYTOPE_SEG_PRISM_TENSOR: 1092 *Nt = 1; 1093 *target = tquadT; 1094 *size = tquadS; 1095 *cone = tquadC; 1096 *ornt = tquadO; 1097 break; 1098 case DM_POLYTOPE_TETRAHEDRON: 1099 *Nt = 1; 1100 *target = tetT; 1101 *size = tetS; 1102 *cone = tetC; 1103 *ornt = tetO; 1104 break; 1105 case DM_POLYTOPE_HEXAHEDRON: 1106 *Nt = 1; 1107 *target = hexT; 1108 *size = hexS; 1109 *cone = hexC; 1110 *ornt = hexO; 1111 break; 1112 case DM_POLYTOPE_TRI_PRISM: 1113 *Nt = 1; 1114 *target = tripT; 1115 *size = tripS; 1116 *cone = tripC; 1117 *ornt = tripO; 1118 break; 1119 case DM_POLYTOPE_TRI_PRISM_TENSOR: 1120 *Nt = 1; 1121 *target = ttripT; 1122 *size = ttripS; 1123 *cone = ttripC; 1124 *ornt = ttripO; 1125 break; 1126 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 1127 *Nt = 1; 1128 *target = tquadpT; 1129 *size = tquadpS; 1130 *cone = tquadpC; 1131 *ornt = tquadpO; 1132 break; 1133 case DM_POLYTOPE_PYRAMID: 1134 *Nt = 1; 1135 *target = pyrT; 1136 *size = pyrS; 1137 *cone = pyrC; 1138 *ornt = pyrO; 1139 break; 1140 default: 1141 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 1142 } 1143 PetscFunctionReturn(PETSC_SUCCESS); 1144 } 1145 1146 /*@ 1147 DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point 1148 1149 Not Collective 1150 1151 Input Parameters: 1152 + tr - The `DMPlexTransform` 1153 . sct - The source point cell type, from whom the new cell is being produced 1154 . sp - The source point 1155 . so - The orientation of the source point in its enclosing parent 1156 . tct - The target point cell type 1157 . r - The replica number requested for the produced cell type 1158 - o - The orientation of the replica 1159 1160 Output Parameters: 1161 + rnew - The replica number, given the orientation of the parent 1162 - onew - The replica orientation, given the orientation of the parent 1163 1164 Level: advanced 1165 1166 .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformCellTransform()`, `DMPlexTransformApply()` 1167 @*/ 1168 PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 1169 { 1170 PetscFunctionBeginHot; 1171 PetscUseTypeMethod(tr, getsubcellorientation, sct, sp, so, tct, r, o, rnew, onew); 1172 PetscFunctionReturn(PETSC_SUCCESS); 1173 } 1174 1175 static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm) 1176 { 1177 DM dm; 1178 PetscInt pStart, pEnd, p, pNew; 1179 1180 PetscFunctionBegin; 1181 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1182 /* Must create the celltype label here so that we do not automatically try to compute the types */ 1183 PetscCall(DMCreateLabel(rdm, "celltype")); 1184 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1185 for (p = pStart; p < pEnd; ++p) { 1186 DMPolytopeType ct; 1187 DMPolytopeType *rct; 1188 PetscInt *rsize, *rcone, *rornt; 1189 PetscInt Nct, n, r; 1190 1191 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1192 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1193 for (n = 0; n < Nct; ++n) { 1194 for (r = 0; r < rsize[n]; ++r) { 1195 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 1196 PetscCall(DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]))); 1197 PetscCall(DMPlexSetCellType(rdm, pNew, rct[n])); 1198 } 1199 } 1200 } 1201 /* Let the DM know we have set all the cell types */ 1202 { 1203 DMLabel ctLabel; 1204 DM_Plex *plex = (DM_Plex *)rdm->data; 1205 1206 PetscCall(DMPlexGetCellTypeLabel(rdm, &ctLabel)); 1207 PetscCall(PetscObjectStateGet((PetscObject)ctLabel, &plex->celltypeState)); 1208 } 1209 PetscFunctionReturn(PETSC_SUCCESS); 1210 } 1211 1212 PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize) 1213 { 1214 DMPolytopeType ctNew; 1215 1216 PetscFunctionBegin; 1217 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1218 PetscValidIntPointer(coneSize, 3); 1219 PetscCall(DMPlexTransformGetCellType(tr, q, &ctNew)); 1220 *coneSize = DMPolytopeTypeGetConeSize((DMPolytopeType)ctNew); 1221 PetscFunctionReturn(PETSC_SUCCESS); 1222 } 1223 1224 /* The orientation o is for the interior of the cell p */ 1225 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[]) 1226 { 1227 DM dm; 1228 const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); 1229 const PetscInt *cone; 1230 PetscInt c, coff = *coneoff, ooff = *orntoff; 1231 1232 PetscFunctionBegin; 1233 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1234 PetscCall(DMPlexGetOrientedCone(dm, p, &cone, NULL)); 1235 for (c = 0; c < csizeNew; ++c) { 1236 PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */ 1237 PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */ 1238 PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */ 1239 DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */ 1240 const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */ 1241 PetscInt pr = -1; /* Replica number of pp that produces new cone point */ 1242 const DMPolytopeType ft = (DMPolytopeType)rcone[coff++]; /* Cell type for new cone point of pNew */ 1243 const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ 1244 PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */ 1245 PetscInt lc; 1246 1247 /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */ 1248 for (lc = 0; lc < fn; ++lc) { 1249 const PetscInt *parr = DMPolytopeTypeGetArrangment(pct, po); 1250 const PetscInt acp = rcone[coff++]; 1251 const PetscInt pcp = parr[acp * 2]; 1252 const PetscInt pco = parr[acp * 2 + 1]; 1253 const PetscInt *ppornt; 1254 1255 ppp = pp; 1256 pp = pcone[pcp]; 1257 PetscCall(DMPlexGetCellType(dm, pp, &pct)); 1258 // Restore the parent cone from the last iterate 1259 if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, ppp, &pcone, NULL)); 1260 PetscCall(DMPlexGetOrientedCone(dm, pp, &pcone, NULL)); 1261 PetscCall(DMPlexGetOrientedCone(dm, ppp, NULL, &ppornt)); 1262 po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco); 1263 PetscCall(DMPlexRestoreOrientedCone(dm, ppp, NULL, &ppornt)); 1264 } 1265 if (lc) PetscCall(DMPlexRestoreOrientedCone(dm, pp, &pcone, NULL)); 1266 pr = rcone[coff++]; 1267 /* Orientation po of pp maps (pr, fo) -> (pr', fo') */ 1268 PetscCall(DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo)); 1269 PetscCall(DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c])); 1270 orntNew[c] = fo; 1271 } 1272 PetscCall(DMPlexRestoreOrientedCone(dm, p, &cone, NULL)); 1273 *coneoff = coff; 1274 *orntoff = ooff; 1275 PetscFunctionReturn(PETSC_SUCCESS); 1276 } 1277 1278 static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm) 1279 { 1280 DM dm; 1281 DMPolytopeType ct; 1282 PetscInt *coneNew, *orntNew; 1283 PetscInt maxConeSize = 0, pStart, pEnd, p, pNew; 1284 1285 PetscFunctionBegin; 1286 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1287 for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p)); 1288 PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew)); 1289 PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew)); 1290 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1291 for (p = pStart; p < pEnd; ++p) { 1292 PetscInt coff, ooff; 1293 DMPolytopeType *rct; 1294 PetscInt *rsize, *rcone, *rornt; 1295 PetscInt Nct, n, r; 1296 1297 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1298 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1299 for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) { 1300 const DMPolytopeType ctNew = rct[n]; 1301 1302 for (r = 0; r < rsize[n]; ++r) { 1303 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 1304 PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew)); 1305 PetscCall(DMPlexSetCone(rdm, pNew, coneNew)); 1306 PetscCall(DMPlexSetConeOrientation(rdm, pNew, orntNew)); 1307 } 1308 } 1309 } 1310 PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew)); 1311 PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew)); 1312 PetscCall(DMViewFromOptions(rdm, NULL, "-rdm_view")); 1313 PetscCall(DMPlexSymmetrize(rdm)); 1314 PetscCall(DMPlexStratify(rdm)); 1315 PetscFunctionReturn(PETSC_SUCCESS); 1316 } 1317 1318 PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[]) 1319 { 1320 DM dm; 1321 DMPolytopeType ct, qct; 1322 DMPolytopeType *rct; 1323 PetscInt *rsize, *rcone, *rornt, *qcone, *qornt; 1324 PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0; 1325 1326 PetscFunctionBegin; 1327 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1328 PetscValidPointer(cone, 4); 1329 PetscValidPointer(ornt, 5); 1330 for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p)); 1331 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1332 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone)); 1333 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt)); 1334 PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r)); 1335 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1336 for (n = 0; n < Nct; ++n) { 1337 const DMPolytopeType ctNew = rct[n]; 1338 const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); 1339 PetscInt Nr = rsize[n], fn, c; 1340 1341 if (ctNew == qct) Nr = r; 1342 for (nr = 0; nr < Nr; ++nr) { 1343 for (c = 0; c < csizeNew; ++c) { 1344 ++coff; /* Cell type of new cone point */ 1345 fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ 1346 coff += fn; 1347 ++coff; /* Replica number of new cone point */ 1348 ++ooff; /* Orientation of new cone point */ 1349 } 1350 } 1351 if (ctNew == qct) break; 1352 } 1353 PetscCall(DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt)); 1354 *cone = qcone; 1355 *ornt = qornt; 1356 PetscFunctionReturn(PETSC_SUCCESS); 1357 } 1358 1359 PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[]) 1360 { 1361 DM dm; 1362 DMPolytopeType ct, qct; 1363 DMPolytopeType *rct; 1364 PetscInt *rsize, *rcone, *rornt, *qcone, *qornt; 1365 PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0; 1366 1367 PetscFunctionBegin; 1368 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1369 if (cone) PetscValidPointer(cone, 3); 1370 if (ornt) PetscValidPointer(ornt, 4); 1371 for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType)p)); 1372 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1373 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone)); 1374 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt)); 1375 PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r)); 1376 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1377 for (n = 0; n < Nct; ++n) { 1378 const DMPolytopeType ctNew = rct[n]; 1379 const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); 1380 PetscInt Nr = rsize[n], fn, c; 1381 1382 if (ctNew == qct) Nr = r; 1383 for (nr = 0; nr < Nr; ++nr) { 1384 for (c = 0; c < csizeNew; ++c) { 1385 ++coff; /* Cell type of new cone point */ 1386 fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ 1387 coff += fn; 1388 ++coff; /* Replica number of new cone point */ 1389 ++ooff; /* Orientation of new cone point */ 1390 } 1391 } 1392 if (ctNew == qct) break; 1393 } 1394 PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt)); 1395 if (cone) *cone = qcone; 1396 else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qcone)); 1397 if (ornt) *ornt = qornt; 1398 else PetscCall(DMRestoreWorkArray(dm, maxConeSize, MPIU_INT, &qornt)); 1399 PetscFunctionReturn(PETSC_SUCCESS); 1400 } 1401 1402 PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[]) 1403 { 1404 DM dm; 1405 1406 PetscFunctionBegin; 1407 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1408 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1409 if (cone) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, cone)); 1410 if (ornt) PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, ornt)); 1411 PetscFunctionReturn(PETSC_SUCCESS); 1412 } 1413 1414 static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr) 1415 { 1416 PetscInt ict; 1417 1418 PetscFunctionBegin; 1419 PetscCall(PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts)); 1420 for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) { 1421 const DMPolytopeType ct = (DMPolytopeType)ict; 1422 DMPlexTransform reftr; 1423 DM refdm, trdm; 1424 Vec coordinates; 1425 const PetscScalar *coords; 1426 DMPolytopeType *rct; 1427 PetscInt *rsize, *rcone, *rornt; 1428 PetscInt Nct, n, r, pNew; 1429 PetscInt trdim, vStart, vEnd, Nc; 1430 const PetscInt debug = 0; 1431 const char *typeName; 1432 1433 /* Since points are 0-dimensional, coordinates make no sense */ 1434 if (DMPolytopeTypeGetDim(ct) <= 0) continue; 1435 PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm)); 1436 PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &reftr)); 1437 PetscCall(DMPlexTransformSetDM(reftr, refdm)); 1438 PetscCall(DMPlexTransformGetType(tr, &typeName)); 1439 PetscCall(DMPlexTransformSetType(reftr, typeName)); 1440 PetscCall(DMPlexTransformSetUp(reftr)); 1441 PetscCall(DMPlexTransformApply(reftr, refdm, &trdm)); 1442 1443 PetscCall(DMGetDimension(trdm, &trdim)); 1444 PetscCall(DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd)); 1445 tr->trNv[ct] = vEnd - vStart; 1446 PetscCall(DMGetCoordinatesLocal(trdm, &coordinates)); 1447 PetscCall(VecGetLocalSize(coordinates, &Nc)); 1448 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); 1449 PetscCall(PetscCalloc1(Nc, &tr->trVerts[ct])); 1450 PetscCall(VecGetArrayRead(coordinates, &coords)); 1451 PetscCall(PetscArraycpy(tr->trVerts[ct], coords, Nc)); 1452 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 1453 1454 PetscCall(PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct])); 1455 PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1456 for (n = 0; n < Nct; ++n) { 1457 /* Since points are 0-dimensional, coordinates make no sense */ 1458 if (rct[n] == DM_POLYTOPE_POINT) continue; 1459 PetscCall(PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]])); 1460 for (r = 0; r < rsize[n]; ++r) { 1461 PetscInt *closure = NULL; 1462 PetscInt clSize, cl, Nv = 0; 1463 1464 PetscCall(PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r])); 1465 PetscCall(DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew)); 1466 PetscCall(DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure)); 1467 for (cl = 0; cl < clSize * 2; cl += 2) { 1468 const PetscInt sv = closure[cl]; 1469 1470 if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart; 1471 } 1472 PetscCall(DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure)); 1473 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]); 1474 } 1475 } 1476 if (debug) { 1477 DMPolytopeType *rct; 1478 PetscInt *rsize, *rcone, *rornt; 1479 PetscInt v, dE = trdim, d, off = 0; 1480 1481 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct])); 1482 for (v = 0; v < tr->trNv[ct]; ++v) { 1483 PetscCall(PetscPrintf(PETSC_COMM_SELF, " ")); 1484 for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++]))); 1485 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 1486 } 1487 1488 PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1489 for (n = 0; n < Nct; ++n) { 1490 if (rct[n] == DM_POLYTOPE_POINT) continue; 1491 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct])); 1492 for (r = 0; r < rsize[n]; ++r) { 1493 PetscCall(PetscPrintf(PETSC_COMM_SELF, " ")); 1494 for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v])); 1495 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 1496 } 1497 } 1498 } 1499 PetscCall(DMDestroy(&refdm)); 1500 PetscCall(DMDestroy(&trdm)); 1501 PetscCall(DMPlexTransformDestroy(&reftr)); 1502 } 1503 PetscFunctionReturn(PETSC_SUCCESS); 1504 } 1505 1506 /*@C 1507 DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type 1508 1509 Input Parameters: 1510 + tr - The `DMPlexTransform` object 1511 - ct - The cell type 1512 1513 Output Parameters: 1514 + Nv - The number of transformed vertices in the closure of the reference cell of given type 1515 - trVerts - The coordinates of these vertices in the reference cell 1516 1517 Level: developer 1518 1519 .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetSubcellVertices()` 1520 @*/ 1521 PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[]) 1522 { 1523 PetscFunctionBegin; 1524 if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr)); 1525 if (Nv) *Nv = tr->trNv[ct]; 1526 if (trVerts) *trVerts = tr->trVerts[ct]; 1527 PetscFunctionReturn(PETSC_SUCCESS); 1528 } 1529 1530 /*@C 1531 DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type 1532 1533 Input Parameters: 1534 + tr - The `DMPlexTransform` object 1535 . ct - The cell type 1536 . rct - The subcell type 1537 - r - The subcell index 1538 1539 Output Parameter: 1540 . subVerts - The indices of these vertices in the set of vertices returned by `DMPlexTransformGetCellVertices()` 1541 1542 Level: developer 1543 1544 .seealso: `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransformGetCellVertices()` 1545 @*/ 1546 PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[]) 1547 { 1548 PetscFunctionBegin; 1549 if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr)); 1550 PetscCheck(tr->trSubVerts[ct][rct], PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]); 1551 if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r]; 1552 PetscFunctionReturn(PETSC_SUCCESS); 1553 } 1554 1555 /* Computes new vertex as the barycenter, or centroid */ 1556 PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) 1557 { 1558 PetscInt v, d; 1559 1560 PetscFunctionBeginHot; 1561 PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]); 1562 for (d = 0; d < dE; ++d) out[d] = 0.0; 1563 for (v = 0; v < Nv; ++v) 1564 for (d = 0; d < dE; ++d) out[d] += in[v * dE + d]; 1565 for (d = 0; d < dE; ++d) out[d] /= Nv; 1566 PetscFunctionReturn(PETSC_SUCCESS); 1567 } 1568 1569 /*@ 1570 DMPlexTransformMapCoordinates - Calculate new coordinates for produced points 1571 1572 Not collective 1573 1574 Input Parameters: 1575 + tr - The `DMPlexTransform` 1576 . pct - The cell type of the parent, from whom the new cell is being produced 1577 . ct - The type being produced 1578 . p - The original point 1579 . r - The replica number requested for the produced cell type 1580 . Nv - Number of vertices in the closure of the parent cell 1581 . dE - Spatial dimension 1582 - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell 1583 1584 Output Parameter: 1585 . out - The coordinates of the new vertices 1586 1587 Level: intermediate 1588 1589 .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `DMPlexTransform`, `DMPlexTransformApply()` 1590 @*/ 1591 PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) 1592 { 1593 PetscFunctionBeginHot; 1594 PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out); 1595 PetscFunctionReturn(PETSC_SUCCESS); 1596 } 1597 1598 static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew) 1599 { 1600 DM dm; 1601 IS valueIS; 1602 const PetscInt *values; 1603 PetscInt defVal, Nv, val; 1604 1605 PetscFunctionBegin; 1606 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1607 PetscCall(DMLabelGetDefaultValue(label, &defVal)); 1608 PetscCall(DMLabelSetDefaultValue(labelNew, defVal)); 1609 PetscCall(DMLabelGetValueIS(label, &valueIS)); 1610 PetscCall(ISGetLocalSize(valueIS, &Nv)); 1611 PetscCall(ISGetIndices(valueIS, &values)); 1612 for (val = 0; val < Nv; ++val) { 1613 IS pointIS; 1614 const PetscInt *points; 1615 PetscInt numPoints, p; 1616 1617 /* Ensure refined label is created with same number of strata as 1618 * original (even if no entries here). */ 1619 PetscCall(DMLabelAddStratum(labelNew, values[val])); 1620 PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS)); 1621 PetscCall(ISGetLocalSize(pointIS, &numPoints)); 1622 PetscCall(ISGetIndices(pointIS, &points)); 1623 for (p = 0; p < numPoints; ++p) { 1624 const PetscInt point = points[p]; 1625 DMPolytopeType ct; 1626 DMPolytopeType *rct; 1627 PetscInt *rsize, *rcone, *rornt; 1628 PetscInt Nct, n, r, pNew = 0; 1629 1630 PetscCall(DMPlexGetCellType(dm, point, &ct)); 1631 PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1632 for (n = 0; n < Nct; ++n) { 1633 for (r = 0; r < rsize[n]; ++r) { 1634 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew)); 1635 PetscCall(DMLabelSetValue(labelNew, pNew, values[val])); 1636 } 1637 } 1638 } 1639 PetscCall(ISRestoreIndices(pointIS, &points)); 1640 PetscCall(ISDestroy(&pointIS)); 1641 } 1642 PetscCall(ISRestoreIndices(valueIS, &values)); 1643 PetscCall(ISDestroy(&valueIS)); 1644 PetscFunctionReturn(PETSC_SUCCESS); 1645 } 1646 1647 static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm) 1648 { 1649 DM dm; 1650 PetscInt numLabels, l; 1651 1652 PetscFunctionBegin; 1653 PetscCall(DMPlexTransformGetDM(tr, &dm)); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 1842 } 1843 1844 /*@C 1845 DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of `DMPolytopeType` ct with localized coordinates x, generate localized coordinates xr for subcell r of type rct. 1846 1847 Not Collective 1848 1849 Input Parameters: 1850 + tr - The `DMPlexTransform` 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: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexTransform`, `DMPolytopeType`, `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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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 // If the original DM was configured from options, the transformed DM should be as well 2141 rdm->setfromoptionscalled = dm->setfromoptionscalled; 2142 2143 *tdm = rdm; 2144 PetscFunctionReturn(PETSC_SUCCESS); 2145 } 2146 2147 PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm) 2148 { 2149 DMPlexTransform tr; 2150 DM cdm, rcdm; 2151 const char *prefix; 2152 2153 PetscFunctionBegin; 2154 PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr)); 2155 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix)); 2156 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix)); 2157 PetscCall(DMPlexTransformSetDM(tr, dm)); 2158 PetscCall(DMPlexTransformSetFromOptions(tr)); 2159 PetscCall(DMPlexTransformSetActive(tr, adaptLabel)); 2160 PetscCall(DMPlexTransformSetUp(tr)); 2161 PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view")); 2162 PetscCall(DMPlexTransformApply(tr, dm, rdm)); 2163 PetscCall(DMCopyDisc(dm, *rdm)); 2164 PetscCall(DMGetCoordinateDM(dm, &cdm)); 2165 PetscCall(DMGetCoordinateDM(*rdm, &rcdm)); 2166 PetscCall(DMCopyDisc(cdm, rcdm)); 2167 PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm)); 2168 PetscCall(DMCopyDisc(dm, *rdm)); 2169 PetscCall(DMPlexTransformDestroy(&tr)); 2170 ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation; 2171 PetscFunctionReturn(PETSC_SUCCESS); 2172 } 2173