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