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 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(DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd)); 1316 tr->trNv[ct] = vEnd - vStart; 1317 PetscCall(DMGetCoordinatesLocal(trdm, &coordinates)); 1318 PetscCall(VecGetLocalSize(coordinates, &Nc)); 1319 PetscCheck(tr->trNv[ct] * DMPolytopeTypeGetDim(ct) == 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] * DMPolytopeTypeGetDim(ct), Nc); 1320 PetscCall(PetscCalloc1(Nc, &tr->trVerts[ct])); 1321 PetscCall(VecGetArrayRead(coordinates, &coords)); 1322 PetscCall(PetscArraycpy(tr->trVerts[ct], coords, Nc)); 1323 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 1324 1325 PetscCall(PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct])); 1326 PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1327 for (n = 0; n < Nct; ++n) { 1328 /* Since points are 0-dimensional, coordinates make no sense */ 1329 if (rct[n] == DM_POLYTOPE_POINT) continue; 1330 PetscCall(PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]])); 1331 for (r = 0; r < rsize[n]; ++r) { 1332 PetscInt *closure = NULL; 1333 PetscInt clSize, cl, Nv = 0; 1334 1335 PetscCall(PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r])); 1336 PetscCall(DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew)); 1337 PetscCall(DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure)); 1338 for (cl = 0; cl < clSize * 2; cl += 2) { 1339 const PetscInt sv = closure[cl]; 1340 1341 if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart; 1342 } 1343 PetscCall(DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure)); 1344 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]); 1345 } 1346 } 1347 if (debug) { 1348 DMPolytopeType *rct; 1349 PetscInt *rsize, *rcone, *rornt; 1350 PetscInt v, dE = DMPolytopeTypeGetDim(ct), d, off = 0; 1351 1352 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct])); 1353 for (v = 0; v < tr->trNv[ct]; ++v) { 1354 PetscCall(PetscPrintf(PETSC_COMM_SELF, " ")); 1355 for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++]))); 1356 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 1357 } 1358 1359 PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1360 for (n = 0; n < Nct; ++n) { 1361 if (rct[n] == DM_POLYTOPE_POINT) continue; 1362 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct])); 1363 for (r = 0; r < rsize[n]; ++r) { 1364 PetscCall(PetscPrintf(PETSC_COMM_SELF, " ")); 1365 for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) { PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v])); } 1366 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 1367 } 1368 } 1369 } 1370 PetscCall(DMDestroy(&refdm)); 1371 PetscCall(DMDestroy(&trdm)); 1372 PetscCall(DMPlexTransformDestroy(&reftr)); 1373 } 1374 PetscFunctionReturn(0); 1375 } 1376 1377 /* 1378 DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type 1379 1380 Input Parameters: 1381 + tr - The DMPlexTransform object 1382 - ct - The cell type 1383 1384 Output Parameters: 1385 + Nv - The number of transformed vertices in the closure of the reference cell of given type 1386 - trVerts - The coordinates of these vertices in the reference cell 1387 1388 Level: developer 1389 1390 .seealso: `DMPlexTransformGetSubcellVertices()` 1391 */ 1392 PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[]) { 1393 PetscFunctionBegin; 1394 if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr)); 1395 if (Nv) *Nv = tr->trNv[ct]; 1396 if (trVerts) *trVerts = tr->trVerts[ct]; 1397 PetscFunctionReturn(0); 1398 } 1399 1400 /* 1401 DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type 1402 1403 Input Parameters: 1404 + tr - The DMPlexTransform object 1405 . ct - The cell type 1406 . rct - The subcell type 1407 - r - The subcell index 1408 1409 Output Parameter: 1410 . subVerts - The indices of these vertices in the set of vertices returned by DMPlexTransformGetCellVertices() 1411 1412 Level: developer 1413 1414 .seealso: `DMPlexTransformGetCellVertices()` 1415 */ 1416 PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[]) { 1417 PetscFunctionBegin; 1418 if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr)); 1419 PetscCheck(tr->trSubVerts[ct][rct], PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]); 1420 if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r]; 1421 PetscFunctionReturn(0); 1422 } 1423 1424 /* Computes new vertex as the barycenter, or centroid */ 1425 PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) { 1426 PetscInt v, d; 1427 1428 PetscFunctionBeginHot; 1429 PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]); 1430 for (d = 0; d < dE; ++d) out[d] = 0.0; 1431 for (v = 0; v < Nv; ++v) 1432 for (d = 0; d < dE; ++d) out[d] += in[v * dE + d]; 1433 for (d = 0; d < dE; ++d) out[d] /= Nv; 1434 PetscFunctionReturn(0); 1435 } 1436 1437 /*@ 1438 DMPlexTransformMapCoordinates - Calculate new coordinates for produced points 1439 1440 Not collective 1441 1442 Input Parameters: 1443 + cr - The DMPlexCellRefiner 1444 . pct - The cell type of the parent, from whom the new cell is being produced 1445 . ct - The type being produced 1446 . p - The original point 1447 . r - The replica number requested for the produced cell type 1448 . Nv - Number of vertices in the closure of the parent cell 1449 . dE - Spatial dimension 1450 - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell 1451 1452 Output Parameter: 1453 . out - The coordinates of the new vertices 1454 1455 Level: intermediate 1456 1457 .seealso: `DMPlexTransformApply()` 1458 @*/ 1459 PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) { 1460 PetscFunctionBeginHot; 1461 PetscUseTypeMethod(tr, mapcoordinates, pct, ct, p, r, Nv, dE, in, out); 1462 PetscFunctionReturn(0); 1463 } 1464 1465 static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew) { 1466 DM dm; 1467 IS valueIS; 1468 const PetscInt *values; 1469 PetscInt defVal, Nv, val; 1470 1471 PetscFunctionBegin; 1472 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1473 PetscCall(DMLabelGetDefaultValue(label, &defVal)); 1474 PetscCall(DMLabelSetDefaultValue(labelNew, defVal)); 1475 PetscCall(DMLabelGetValueIS(label, &valueIS)); 1476 PetscCall(ISGetLocalSize(valueIS, &Nv)); 1477 PetscCall(ISGetIndices(valueIS, &values)); 1478 for (val = 0; val < Nv; ++val) { 1479 IS pointIS; 1480 const PetscInt *points; 1481 PetscInt numPoints, p; 1482 1483 /* Ensure refined label is created with same number of strata as 1484 * original (even if no entries here). */ 1485 PetscCall(DMLabelAddStratum(labelNew, values[val])); 1486 PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS)); 1487 PetscCall(ISGetLocalSize(pointIS, &numPoints)); 1488 PetscCall(ISGetIndices(pointIS, &points)); 1489 for (p = 0; p < numPoints; ++p) { 1490 const PetscInt point = points[p]; 1491 DMPolytopeType ct; 1492 DMPolytopeType *rct; 1493 PetscInt *rsize, *rcone, *rornt; 1494 PetscInt Nct, n, r, pNew = 0; 1495 1496 PetscCall(DMPlexGetCellType(dm, point, &ct)); 1497 PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1498 for (n = 0; n < Nct; ++n) { 1499 for (r = 0; r < rsize[n]; ++r) { 1500 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew)); 1501 PetscCall(DMLabelSetValue(labelNew, pNew, values[val])); 1502 } 1503 } 1504 } 1505 PetscCall(ISRestoreIndices(pointIS, &points)); 1506 PetscCall(ISDestroy(&pointIS)); 1507 } 1508 PetscCall(ISRestoreIndices(valueIS, &values)); 1509 PetscCall(ISDestroy(&valueIS)); 1510 PetscFunctionReturn(0); 1511 } 1512 1513 static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm) { 1514 DM dm; 1515 PetscInt numLabels, l; 1516 1517 PetscFunctionBegin; 1518 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1519 PetscCall(DMGetNumLabels(dm, &numLabels)); 1520 for (l = 0; l < numLabels; ++l) { 1521 DMLabel label, labelNew; 1522 const char *lname; 1523 PetscBool isDepth, isCellType; 1524 1525 PetscCall(DMGetLabelName(dm, l, &lname)); 1526 PetscCall(PetscStrcmp(lname, "depth", &isDepth)); 1527 if (isDepth) continue; 1528 PetscCall(PetscStrcmp(lname, "celltype", &isCellType)); 1529 if (isCellType) continue; 1530 PetscCall(DMCreateLabel(rdm, lname)); 1531 PetscCall(DMGetLabel(dm, lname, &label)); 1532 PetscCall(DMGetLabel(rdm, lname, &labelNew)); 1533 PetscCall(RefineLabel_Internal(tr, label, labelNew)); 1534 } 1535 PetscFunctionReturn(0); 1536 } 1537 1538 /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */ 1539 PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm) { 1540 DM dm; 1541 PetscInt Nf, f, Nds, s; 1542 1543 PetscFunctionBegin; 1544 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1545 PetscCall(DMGetNumFields(dm, &Nf)); 1546 for (f = 0; f < Nf; ++f) { 1547 DMLabel label, labelNew; 1548 PetscObject obj; 1549 const char *lname; 1550 1551 PetscCall(DMGetField(rdm, f, &label, &obj)); 1552 if (!label) continue; 1553 PetscCall(PetscObjectGetName((PetscObject)label, &lname)); 1554 PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew)); 1555 PetscCall(RefineLabel_Internal(tr, label, labelNew)); 1556 PetscCall(DMSetField_Internal(rdm, f, labelNew, obj)); 1557 PetscCall(DMLabelDestroy(&labelNew)); 1558 } 1559 PetscCall(DMGetNumDS(dm, &Nds)); 1560 for (s = 0; s < Nds; ++s) { 1561 DMLabel label, labelNew; 1562 const char *lname; 1563 1564 PetscCall(DMGetRegionNumDS(rdm, s, &label, NULL, NULL)); 1565 if (!label) continue; 1566 PetscCall(PetscObjectGetName((PetscObject)label, &lname)); 1567 PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew)); 1568 PetscCall(RefineLabel_Internal(tr, label, labelNew)); 1569 PetscCall(DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL)); 1570 PetscCall(DMLabelDestroy(&labelNew)); 1571 } 1572 PetscFunctionReturn(0); 1573 } 1574 1575 static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm) { 1576 DM dm; 1577 PetscSF sf, sfNew; 1578 PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m; 1579 const PetscInt *localPoints; 1580 const PetscSFNode *remotePoints; 1581 PetscInt *localPointsNew; 1582 PetscSFNode *remotePointsNew; 1583 PetscInt pStartNew, pEndNew, pNew; 1584 /* Brute force algorithm */ 1585 PetscSF rsf; 1586 PetscSection s; 1587 const PetscInt *rootdegree; 1588 PetscInt *rootPointsNew, *remoteOffsets; 1589 PetscInt numPointsNew, pStart, pEnd, p; 1590 1591 PetscFunctionBegin; 1592 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1593 PetscCall(DMPlexGetChart(rdm, &pStartNew, &pEndNew)); 1594 PetscCall(DMGetPointSF(dm, &sf)); 1595 PetscCall(DMGetPointSF(rdm, &sfNew)); 1596 /* Calculate size of new SF */ 1597 PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints)); 1598 if (numRoots < 0) PetscFunctionReturn(0); 1599 for (l = 0; l < numLeaves; ++l) { 1600 const PetscInt p = localPoints[l]; 1601 DMPolytopeType ct; 1602 DMPolytopeType *rct; 1603 PetscInt *rsize, *rcone, *rornt; 1604 PetscInt Nct, n; 1605 1606 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1607 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1608 for (n = 0; n < Nct; ++n) { numLeavesNew += rsize[n]; } 1609 } 1610 /* Send new root point numbers 1611 It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication 1612 */ 1613 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1614 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &s)); 1615 PetscCall(PetscSectionSetChart(s, pStart, pEnd)); 1616 for (p = pStart; p < pEnd; ++p) { 1617 DMPolytopeType ct; 1618 DMPolytopeType *rct; 1619 PetscInt *rsize, *rcone, *rornt; 1620 PetscInt Nct, n; 1621 1622 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1623 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1624 for (n = 0; n < Nct; ++n) { PetscCall(PetscSectionAddDof(s, p, rsize[n])); } 1625 } 1626 PetscCall(PetscSectionSetUp(s)); 1627 PetscCall(PetscSectionGetStorageSize(s, &numPointsNew)); 1628 PetscCall(PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets)); 1629 PetscCall(PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf)); 1630 PetscCall(PetscFree(remoteOffsets)); 1631 PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree)); 1632 PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree)); 1633 PetscCall(PetscMalloc1(numPointsNew, &rootPointsNew)); 1634 for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1; 1635 for (p = pStart; p < pEnd; ++p) { 1636 DMPolytopeType ct; 1637 DMPolytopeType *rct; 1638 PetscInt *rsize, *rcone, *rornt; 1639 PetscInt Nct, n, r, off; 1640 1641 if (!rootdegree[p - pStart]) continue; 1642 PetscCall(PetscSectionGetOffset(s, p, &off)); 1643 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1644 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1645 for (n = 0, m = 0; n < Nct; ++n) { 1646 for (r = 0; r < rsize[n]; ++r, ++m) { 1647 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 1648 rootPointsNew[off + m] = pNew; 1649 } 1650 } 1651 } 1652 PetscCall(PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE)); 1653 PetscCall(PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew, MPI_REPLACE)); 1654 PetscCall(PetscSFDestroy(&rsf)); 1655 PetscCall(PetscMalloc1(numLeavesNew, &localPointsNew)); 1656 PetscCall(PetscMalloc1(numLeavesNew, &remotePointsNew)); 1657 for (l = 0, m = 0; l < numLeaves; ++l) { 1658 const PetscInt p = localPoints[l]; 1659 DMPolytopeType ct; 1660 DMPolytopeType *rct; 1661 PetscInt *rsize, *rcone, *rornt; 1662 PetscInt Nct, n, r, q, off; 1663 1664 PetscCall(PetscSectionGetOffset(s, p, &off)); 1665 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1666 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1667 for (n = 0, q = 0; n < Nct; ++n) { 1668 for (r = 0; r < rsize[n]; ++r, ++m, ++q) { 1669 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 1670 localPointsNew[m] = pNew; 1671 remotePointsNew[m].index = rootPointsNew[off + q]; 1672 remotePointsNew[m].rank = remotePoints[l].rank; 1673 } 1674 } 1675 } 1676 PetscCall(PetscSectionDestroy(&s)); 1677 PetscCall(PetscFree(rootPointsNew)); 1678 /* SF needs sorted leaves to correctly calculate Gather */ 1679 { 1680 PetscSFNode *rp, *rtmp; 1681 PetscInt *lp, *idx, *ltmp, i; 1682 1683 PetscCall(PetscMalloc1(numLeavesNew, &idx)); 1684 PetscCall(PetscMalloc1(numLeavesNew, &lp)); 1685 PetscCall(PetscMalloc1(numLeavesNew, &rp)); 1686 for (i = 0; i < numLeavesNew; ++i) { 1687 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); 1688 idx[i] = i; 1689 } 1690 PetscCall(PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx)); 1691 for (i = 0; i < numLeavesNew; ++i) { 1692 lp[i] = localPointsNew[idx[i]]; 1693 rp[i] = remotePointsNew[idx[i]]; 1694 } 1695 ltmp = localPointsNew; 1696 localPointsNew = lp; 1697 rtmp = remotePointsNew; 1698 remotePointsNew = rp; 1699 PetscCall(PetscFree(idx)); 1700 PetscCall(PetscFree(ltmp)); 1701 PetscCall(PetscFree(rtmp)); 1702 } 1703 PetscCall(PetscSFSetGraph(sfNew, pEndNew - pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER)); 1704 PetscFunctionReturn(0); 1705 } 1706 1707 /* 1708 DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of type ct with localized coordinates x, we generate localized coordinates xr for subcell r of type rct. 1709 1710 Not collective 1711 1712 Input Parameters: 1713 + cr - The DMPlexCellRefiner 1714 . ct - The type of the parent cell 1715 . rct - The type of the produced cell 1716 . r - The index of the produced cell 1717 - x - The localized coordinates for the parent cell 1718 1719 Output Parameter: 1720 . xr - The localized coordinates for the produced cell 1721 1722 Level: developer 1723 1724 .seealso: `DMPlexCellRefinerSetCoordinates()` 1725 */ 1726 static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[]) { 1727 PetscFE fe = NULL; 1728 PetscInt cdim, v, *subcellV; 1729 1730 PetscFunctionBegin; 1731 PetscCall(DMPlexTransformGetCoordinateFE(tr, ct, &fe)); 1732 PetscCall(DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV)); 1733 PetscCall(PetscFEGetNumComponents(fe, &cdim)); 1734 for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) { PetscCall(PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v * cdim])); } 1735 PetscFunctionReturn(0); 1736 } 1737 1738 static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm) { 1739 DM dm, cdm, cdmCell, cdmNew, cdmCellNew; 1740 PetscSection coordSection, coordSectionNew, coordSectionCell, coordSectionCellNew; 1741 Vec coordsLocal, coordsLocalNew, coordsLocalCell = NULL, coordsLocalCellNew; 1742 const PetscScalar *coords; 1743 PetscScalar *coordsNew; 1744 const PetscReal *maxCell, *Lstart, *L; 1745 PetscBool localized, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE; 1746 PetscInt dE, dEo, d, cStart, cEnd, c, cStartNew, cEndNew, vStartNew, vEndNew, v, pStart, pEnd, p; 1747 1748 PetscFunctionBegin; 1749 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1750 PetscCall(DMGetCoordinateDM(dm, &cdm)); 1751 PetscCall(DMGetCellCoordinateDM(dm, &cdmCell)); 1752 PetscCall(DMGetCoordinatesLocalized(dm, &localized)); 1753 PetscCall(DMGetPeriodicity(dm, &maxCell, &Lstart, &L)); 1754 if (localized) { 1755 /* Localize coordinates of new vertices */ 1756 localizeVertices = PETSC_TRUE; 1757 /* If we do not have a mechanism for automatically localizing cell coordinates, we need to compute them explicitly for every divided cell */ 1758 if (!maxCell) localizeCells = PETSC_TRUE; 1759 } 1760 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 1761 PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &dEo)); 1762 if (maxCell) { 1763 PetscReal maxCellNew[3]; 1764 1765 for (d = 0; d < dEo; ++d) maxCellNew[d] = maxCell[d] / 2.0; 1766 PetscCall(DMSetPeriodicity(rdm, maxCellNew, Lstart, L)); 1767 } 1768 PetscCall(DMGetCoordinateDim(rdm, &dE)); 1769 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionNew)); 1770 PetscCall(PetscSectionSetNumFields(coordSectionNew, 1)); 1771 PetscCall(PetscSectionSetFieldComponents(coordSectionNew, 0, dE)); 1772 PetscCall(DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew)); 1773 PetscCall(PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew)); 1774 /* Localization should be inherited */ 1775 /* Stefano calculates parent cells for each new cell for localization */ 1776 /* Localized cells need coordinates of closure */ 1777 for (v = vStartNew; v < vEndNew; ++v) { 1778 PetscCall(PetscSectionSetDof(coordSectionNew, v, dE)); 1779 PetscCall(PetscSectionSetFieldDof(coordSectionNew, v, 0, dE)); 1780 } 1781 PetscCall(PetscSectionSetUp(coordSectionNew)); 1782 PetscCall(DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew)); 1783 1784 if (localizeCells) { 1785 PetscCall(DMGetCoordinateDM(rdm, &cdmNew)); 1786 PetscCall(DMClone(cdmNew, &cdmCellNew)); 1787 PetscCall(DMSetCellCoordinateDM(rdm, cdmCellNew)); 1788 PetscCall(DMDestroy(&cdmCellNew)); 1789 1790 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)rdm), &coordSectionCellNew)); 1791 PetscCall(PetscSectionSetNumFields(coordSectionCellNew, 1)); 1792 PetscCall(PetscSectionSetFieldComponents(coordSectionCellNew, 0, dE)); 1793 PetscCall(DMPlexGetHeightStratum(rdm, 0, &cStartNew, &cEndNew)); 1794 PetscCall(PetscSectionSetChart(coordSectionCellNew, cStartNew, cEndNew)); 1795 1796 PetscCall(DMGetCellCoordinateSection(dm, &coordSectionCell)); 1797 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 1798 for (c = cStart; c < cEnd; ++c) { 1799 PetscInt dof; 1800 1801 PetscCall(PetscSectionGetDof(coordSectionCell, c, &dof)); 1802 if (dof) { 1803 DMPolytopeType ct; 1804 DMPolytopeType *rct; 1805 PetscInt *rsize, *rcone, *rornt; 1806 PetscInt dim, cNew, Nct, n, r; 1807 1808 PetscCall(DMPlexGetCellType(dm, c, &ct)); 1809 dim = DMPolytopeTypeGetDim(ct); 1810 PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1811 /* This allows for different cell types */ 1812 for (n = 0; n < Nct; ++n) { 1813 if (dim != DMPolytopeTypeGetDim(rct[n])) continue; 1814 for (r = 0; r < rsize[n]; ++r) { 1815 PetscInt *closure = NULL; 1816 PetscInt clSize, cl, Nv = 0; 1817 1818 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew)); 1819 PetscCall(DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure)); 1820 for (cl = 0; cl < clSize * 2; cl += 2) { 1821 if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv; 1822 } 1823 PetscCall(DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure)); 1824 PetscCall(PetscSectionSetDof(coordSectionCellNew, cNew, Nv * dE)); 1825 PetscCall(PetscSectionSetFieldDof(coordSectionCellNew, cNew, 0, Nv * dE)); 1826 } 1827 } 1828 } 1829 } 1830 PetscCall(PetscSectionSetUp(coordSectionCellNew)); 1831 PetscCall(DMSetCellCoordinateSection(rdm, PETSC_DETERMINE, coordSectionCellNew)); 1832 } 1833 PetscCall(DMViewFromOptions(dm, NULL, "-coarse_dm_view")); 1834 { 1835 VecType vtype; 1836 PetscInt coordSizeNew, bs; 1837 const char *name; 1838 1839 PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal)); 1840 PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalNew)); 1841 PetscCall(PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew)); 1842 PetscCall(VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE)); 1843 PetscCall(PetscObjectGetName((PetscObject)coordsLocal, &name)); 1844 PetscCall(PetscObjectSetName((PetscObject)coordsLocalNew, name)); 1845 PetscCall(VecGetBlockSize(coordsLocal, &bs)); 1846 PetscCall(VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE)); 1847 PetscCall(VecGetType(coordsLocal, &vtype)); 1848 PetscCall(VecSetType(coordsLocalNew, vtype)); 1849 } 1850 PetscCall(VecGetArrayRead(coordsLocal, &coords)); 1851 PetscCall(VecGetArray(coordsLocalNew, &coordsNew)); 1852 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1853 /* First set coordinates for vertices */ 1854 for (p = pStart; p < pEnd; ++p) { 1855 DMPolytopeType ct; 1856 DMPolytopeType *rct; 1857 PetscInt *rsize, *rcone, *rornt; 1858 PetscInt Nct, n, r; 1859 PetscBool hasVertex = PETSC_FALSE; 1860 1861 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1862 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1863 for (n = 0; n < Nct; ++n) { 1864 if (rct[n] == DM_POLYTOPE_POINT) { 1865 hasVertex = PETSC_TRUE; 1866 break; 1867 } 1868 } 1869 if (hasVertex) { 1870 const PetscScalar *icoords = NULL; 1871 const PetscScalar *array = NULL; 1872 PetscScalar *pcoords = NULL; 1873 PetscBool isDG; 1874 PetscInt Nc, Nv, v, d; 1875 1876 PetscCall(DMPlexGetCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords)); 1877 1878 icoords = pcoords; 1879 Nv = Nc / dEo; 1880 if (ct != DM_POLYTOPE_POINT) { 1881 if (localizeVertices && maxCell) { 1882 PetscScalar anchor[3]; 1883 1884 for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d]; 1885 for (v = 0; v < Nv; ++v) PetscCall(DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v * dEo], &pcoords[v * dEo])); 1886 } 1887 } 1888 for (n = 0; n < Nct; ++n) { 1889 if (rct[n] != DM_POLYTOPE_POINT) continue; 1890 for (r = 0; r < rsize[n]; ++r) { 1891 PetscScalar vcoords[3]; 1892 PetscInt vNew, off; 1893 1894 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew)); 1895 PetscCall(PetscSectionGetOffset(coordSectionNew, vNew, &off)); 1896 PetscCall(DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords)); 1897 PetscCall(DMPlexSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off])); 1898 } 1899 } 1900 PetscCall(DMPlexRestoreCellCoordinates(dm, p, &isDG, &Nc, &array, &pcoords)); 1901 } 1902 } 1903 PetscCall(VecRestoreArrayRead(coordsLocal, &coords)); 1904 PetscCall(VecRestoreArray(coordsLocalNew, &coordsNew)); 1905 PetscCall(DMSetCoordinatesLocal(rdm, coordsLocalNew)); 1906 PetscCall(VecDestroy(&coordsLocalNew)); 1907 PetscCall(PetscSectionDestroy(&coordSectionNew)); 1908 /* Then set coordinates for cells by localizing */ 1909 if (!localizeCells) PetscCall(DMLocalizeCoordinates(rdm)); 1910 else { 1911 VecType vtype; 1912 PetscInt coordSizeNew, bs; 1913 const char *name; 1914 1915 PetscCall(DMGetCellCoordinatesLocal(dm, &coordsLocalCell)); 1916 PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalCellNew)); 1917 PetscCall(PetscSectionGetStorageSize(coordSectionCellNew, &coordSizeNew)); 1918 PetscCall(VecSetSizes(coordsLocalCellNew, coordSizeNew, PETSC_DETERMINE)); 1919 PetscCall(PetscObjectGetName((PetscObject)coordsLocalCell, &name)); 1920 PetscCall(PetscObjectSetName((PetscObject)coordsLocalCellNew, name)); 1921 PetscCall(VecGetBlockSize(coordsLocalCell, &bs)); 1922 PetscCall(VecSetBlockSize(coordsLocalCellNew, dEo == dE ? bs : dE)); 1923 PetscCall(VecGetType(coordsLocalCell, &vtype)); 1924 PetscCall(VecSetType(coordsLocalCellNew, vtype)); 1925 PetscCall(VecGetArrayRead(coordsLocalCell, &coords)); 1926 PetscCall(VecGetArray(coordsLocalCellNew, &coordsNew)); 1927 1928 for (p = pStart; p < pEnd; ++p) { 1929 DMPolytopeType ct; 1930 DMPolytopeType *rct; 1931 PetscInt *rsize, *rcone, *rornt; 1932 PetscInt dof = 0, Nct, n, r; 1933 1934 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1935 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1936 if (p >= cStart && p < cEnd) PetscCall(PetscSectionGetDof(coordSectionCell, p, &dof)); 1937 if (dof) { 1938 const PetscScalar *pcoords; 1939 1940 PetscCall(DMPlexPointLocalRead(cdmCell, p, coords, &pcoords)); 1941 for (n = 0; n < Nct; ++n) { 1942 const PetscInt Nr = rsize[n]; 1943 1944 if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue; 1945 for (r = 0; r < Nr; ++r) { 1946 PetscInt pNew, offNew; 1947 1948 /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that 1949 DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger 1950 cell to the ones it produces. */ 1951 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 1952 PetscCall(PetscSectionGetOffset(coordSectionCellNew, pNew, &offNew)); 1953 PetscCall(DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew])); 1954 } 1955 } 1956 } 1957 } 1958 PetscCall(VecRestoreArrayRead(coordsLocalCell, &coords)); 1959 PetscCall(VecRestoreArray(coordsLocalCellNew, &coordsNew)); 1960 PetscCall(DMSetCellCoordinatesLocal(rdm, coordsLocalCellNew)); 1961 PetscCall(VecDestroy(&coordsLocalCellNew)); 1962 PetscCall(PetscSectionDestroy(&coordSectionCellNew)); 1963 } 1964 PetscFunctionReturn(0); 1965 } 1966 1967 PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm) { 1968 DM rdm; 1969 DMPlexInterpolatedFlag interp; 1970 1971 PetscFunctionBegin; 1972 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1973 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1974 PetscValidPointer(tdm, 3); 1975 PetscCall(DMPlexTransformSetDM(tr, dm)); 1976 1977 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &rdm)); 1978 PetscCall(DMSetType(rdm, DMPLEX)); 1979 PetscCall(DMPlexTransformSetDimensions(tr, dm, rdm)); 1980 /* Calculate number of new points of each depth */ 1981 PetscCall(DMPlexIsInterpolatedCollective(dm, &interp)); 1982 PetscCheck(interp == DMPLEX_INTERPOLATED_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement"); 1983 /* Step 1: Set chart */ 1984 PetscCall(DMPlexSetChart(rdm, 0, tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]])); 1985 /* Step 2: Set cone/support sizes (automatically stratifies) */ 1986 PetscCall(DMPlexTransformSetConeSizes(tr, rdm)); 1987 /* Step 3: Setup refined DM */ 1988 PetscCall(DMSetUp(rdm)); 1989 /* Step 4: Set cones and supports (automatically symmetrizes) */ 1990 PetscCall(DMPlexTransformSetCones(tr, rdm)); 1991 /* Step 5: Create pointSF */ 1992 PetscCall(DMPlexTransformCreateSF(tr, rdm)); 1993 /* Step 6: Create labels */ 1994 PetscCall(DMPlexTransformCreateLabels(tr, rdm)); 1995 /* Step 7: Set coordinates */ 1996 PetscCall(DMPlexTransformSetCoordinates(tr, rdm)); 1997 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, rdm)); 1998 *tdm = rdm; 1999 PetscFunctionReturn(0); 2000 } 2001 2002 PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm) { 2003 DMPlexTransform tr; 2004 DM cdm, rcdm; 2005 const char *prefix; 2006 2007 PetscFunctionBegin; 2008 PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr)); 2009 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix)); 2010 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix)); 2011 PetscCall(DMPlexTransformSetDM(tr, dm)); 2012 PetscCall(DMPlexTransformSetFromOptions(tr)); 2013 PetscCall(DMPlexTransformSetActive(tr, adaptLabel)); 2014 PetscCall(DMPlexTransformSetUp(tr)); 2015 PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view")); 2016 PetscCall(DMPlexTransformApply(tr, dm, rdm)); 2017 PetscCall(DMCopyDisc(dm, *rdm)); 2018 PetscCall(DMGetCoordinateDM(dm, &cdm)); 2019 PetscCall(DMGetCoordinateDM(*rdm, &rcdm)); 2020 PetscCall(DMCopyDisc(cdm, rcdm)); 2021 PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm)); 2022 PetscCall(DMCopyDisc(dm, *rdm)); 2023 PetscCall(DMPlexTransformDestroy(&tr)); 2024 ((DM_Plex *)(*rdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation; 2025 PetscFunctionReturn(0); 2026 } 2027