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