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