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