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