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