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