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