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