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 if (tr->ops->destroy) PetscCall((*tr->ops->destroy)(tr)); 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 if (tr->ops->view) PetscCall((*tr->ops->view)(tr, 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 if (tr->ops->setfromoptions) PetscCall((*tr->ops->setfromoptions)(PetscOptionsObject,tr)); 341 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 342 PetscCall(PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) tr)); 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 if (tr->ops->setup) PetscCall((*tr->ops->setup)(tr)); 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) { 679 PetscCall((*tr->ops->setdimensions)(tr, dm, tdm)); 680 } else { 681 PetscInt dim, cdim; 682 683 PetscCall(DMGetDimension(dm, &dim)); 684 PetscCall(DMSetDimension(tdm, dim)); 685 PetscCall(DMGetCoordinateDim(dm, &cdim)); 686 PetscCall(DMSetCoordinateDim(tdm, cdim)); 687 } 688 PetscFunctionReturn(0); 689 } 690 691 /*@ 692 DMPlexTransformGetTargetPoint - Get the number of a point in the transformed mesh based on information from the original mesh. 693 694 Not collective 695 696 Input Parameters: 697 + tr - The DMPlexTransform 698 . ct - The type of the original point which produces the new point 699 . ctNew - The type of the new point 700 . p - The original point which produces the new point 701 - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p 702 703 Output Parameters: 704 . pNew - The new point number 705 706 Level: developer 707 708 .seealso: `DMPlexTransformGetSourcePoint()`, `DMPlexTransformCellTransform()` 709 @*/ 710 PetscErrorCode DMPlexTransformGetTargetPoint(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType ctNew, PetscInt p, PetscInt r, PetscInt *pNew) 711 { 712 DMPolytopeType *rct; 713 PetscInt *rsize, *cone, *ornt; 714 PetscInt rt, Nct, n, off, rp; 715 DMLabel trType = tr->trType; 716 PetscInt ctS = tr->ctStart[ct], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ct]+1]]; 717 PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew]+1]]; 718 PetscInt newp = ctSN, cind; 719 720 PetscFunctionBeginHot; 721 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); 722 PetscCall(DMPlexTransformCellTransform(tr, ct, p, &rt, &Nct, &rct, &rsize, &cone, &ornt)); 723 if (trType) { 724 PetscCall(DMLabelGetValueIndex(trType, rt, &cind)); 725 PetscCall(DMLabelGetStratumPointIndex(trType, rt, p, &rp)); 726 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); 727 } else { 728 cind = ct; 729 rp = p - ctS; 730 } 731 off = tr->offset[cind*DM_NUM_POLYTOPES + ctNew]; 732 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); 733 newp += off; 734 for (n = 0; n < Nct; ++n) { 735 if (rct[n] == ctNew) { 736 if (rsize[n] && r >= rsize[n]) 737 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]); 738 newp += rp * rsize[n] + r; 739 break; 740 } 741 } 742 743 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); 744 *pNew = newp; 745 PetscFunctionReturn(0); 746 } 747 748 /*@ 749 DMPlexTransformGetSourcePoint - Get the number of a point in the original mesh based on information from the transformed mesh. 750 751 Not collective 752 753 Input Parameters: 754 + tr - The DMPlexTransform 755 - pNew - The new point number 756 757 Output Parameters: 758 + ct - The type of the original point which produces the new point 759 . ctNew - The type of the new point 760 . p - The original point which produces the new point 761 - r - The replica number of the new point, meaning it is the rth point of type ctNew produced from p 762 763 Level: developer 764 765 .seealso: `DMPlexTransformGetTargetPoint()`, `DMPlexTransformCellTransform()` 766 @*/ 767 PetscErrorCode DMPlexTransformGetSourcePoint(DMPlexTransform tr, PetscInt pNew, DMPolytopeType *ct, DMPolytopeType *ctNew, PetscInt *p, PetscInt *r) 768 { 769 DMLabel trType = tr->trType; 770 DMPolytopeType *rct; 771 PetscInt *rsize, *cone, *ornt; 772 PetscInt rt, Nct, n, rp = 0, rO = 0, pO; 773 PetscInt offset = -1, ctS, ctE, ctO = 0, ctN, ctTmp; 774 775 PetscFunctionBegin; 776 for (ctN = 0; ctN < DM_NUM_POLYTOPES; ++ctN) { 777 PetscInt ctSN = tr->ctStartNew[ctN], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctN]+1]]; 778 779 if ((pNew >= ctSN) && (pNew < ctEN)) break; 780 } 781 PetscCheck(ctN != DM_NUM_POLYTOPES,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cell type for target point %" PetscInt_FMT " could be not found", pNew); 782 if (trType) { 783 DM dm; 784 IS rtIS; 785 const PetscInt *reftypes; 786 PetscInt Nrt, r, rtStart; 787 788 PetscCall(DMPlexTransformGetDM(tr, &dm)); 789 PetscCall(DMLabelGetNumValues(trType, &Nrt)); 790 PetscCall(DMLabelGetValueIS(trType, &rtIS)); 791 PetscCall(ISGetIndices(rtIS, &reftypes)); 792 for (r = 0; r < Nrt; ++r) { 793 const PetscInt off = tr->offset[r*DM_NUM_POLYTOPES + ctN]; 794 795 if (tr->ctStartNew[ctN] + off > pNew) continue; 796 /* Check that any of this refinement type exist */ 797 /* TODO Actually keep track of the number produced here instead */ 798 if (off > offset) {rt = reftypes[r]; offset = off;} 799 } 800 PetscCall(ISRestoreIndices(rtIS, &reftypes)); 801 PetscCall(ISDestroy(&rtIS)); 802 PetscCheck(offset >= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew); 803 /* TODO Map refinement types to cell types */ 804 PetscCall(DMLabelGetStratumBounds(trType, rt, &rtStart, NULL)); 805 PetscCheck(rtStart >= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Refinement type %" PetscInt_FMT " has no source points", rt); 806 for (ctO = 0; ctO < DM_NUM_POLYTOPES; ++ctO) { 807 PetscInt ctS = tr->ctStart[ctO], ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO]+1]]; 808 809 if ((rtStart >= ctS) && (rtStart < ctE)) break; 810 } 811 PetscCheck(ctO != DM_NUM_POLYTOPES,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine a cell type for refinement type %" PetscInt_FMT, rt); 812 } else { 813 for (ctTmp = 0; ctTmp < DM_NUM_POLYTOPES; ++ctTmp) { 814 const PetscInt off = tr->offset[ctTmp*DM_NUM_POLYTOPES + ctN]; 815 816 if (tr->ctStartNew[ctN] + off > pNew) continue; 817 if (tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctTmp]+1]] <= tr->ctStart[ctTmp]) continue; 818 /* TODO Actually keep track of the number produced here instead */ 819 if (off > offset) {ctO = ctTmp; offset = off;} 820 } 821 PetscCheck(offset >= 0,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Source cell type for target point %" PetscInt_FMT " could be not found", pNew); 822 } 823 ctS = tr->ctStart[ctO]; 824 ctE = tr->ctStart[tr->ctOrderOld[tr->ctOrderInvOld[ctO]+1]]; 825 PetscCall(DMPlexTransformCellTransform(tr, (DMPolytopeType) ctO, ctS, &rt, &Nct, &rct, &rsize, &cone, &ornt)); 826 for (n = 0; n < Nct; ++n) { 827 if ((PetscInt) rct[n] == ctN) { 828 PetscInt tmp = pNew - tr->ctStartNew[ctN] - offset; 829 830 rp = tmp / rsize[n]; 831 rO = tmp % rsize[n]; 832 break; 833 } 834 } 835 PetscCheck(n != Nct,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Replica number for target point %" PetscInt_FMT " could be not found", pNew); 836 pO = rp + ctS; 837 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); 838 if (ct) *ct = (DMPolytopeType) ctO; 839 if (ctNew) *ctNew = (DMPolytopeType) ctN; 840 if (p) *p = pO; 841 if (r) *r = rO; 842 PetscFunctionReturn(0); 843 } 844 845 /*@ 846 DMPlexTransformCellTransform - Describes the transform of a given source cell into a set of other target cells. These produced cells become the new mesh. 847 848 Input Parameters: 849 + tr - The DMPlexTransform object 850 . source - The source cell type 851 - p - The source point, which can also determine the refine type 852 853 Output Parameters: 854 + rt - The refine type for this point 855 . Nt - The number of types produced by this point 856 . target - An array of length Nt giving the types produced 857 . size - An array of length Nt giving the number of cells of each type produced 858 . cone - An array of length Nt*size[t]*coneSize[t] giving the cell type for each point in the cone of each produced point 859 - ornt - An array of length Nt*size[t]*coneSize[t] giving the orientation for each point in the cone of each produced point 860 861 Note: The cone array gives the cone of each subcell listed by the first three outputs. For each cone point, we 862 need the cell type, point identifier, and orientation within the subcell. The orientation is with respect to the canonical 863 division (described in these outputs) of the cell in the original mesh. The point identifier is given by 864 $ the number of cones to be taken, or 0 for the current cell 865 $ the cell cone point number at each level from which it is subdivided 866 $ the replica number r of the subdivision. 867 The orientation is with respect to the canonical cone orientation. For example, the prescription for edge division is 868 $ Nt = 2 869 $ target = {DM_POLYTOPE_POINT, DM_POLYTOPE_SEGMENT} 870 $ size = {1, 2} 871 $ cone = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0} 872 $ ornt = { 0, 0, 0, 0} 873 874 Level: advanced 875 876 .seealso: `DMPlexTransformApply()`, `DMPlexTransformCreate()` 877 @*/ 878 PetscErrorCode DMPlexTransformCellTransform(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 879 { 880 PetscFunctionBegin; 881 PetscCall((*tr->ops->celltransform)(tr, source, p, rt, Nt, target, size, cone, ornt)); 882 PetscFunctionReturn(0); 883 } 884 885 PetscErrorCode DMPlexTransformGetSubcellOrientationIdentity(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 886 { 887 PetscFunctionBegin; 888 *rnew = r; 889 *onew = DMPolytopeTypeComposeOrientation(tct, o, so); 890 PetscFunctionReturn(0); 891 } 892 893 /* Returns the same thing */ 894 PetscErrorCode DMPlexTransformCellTransformIdentity(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 895 { 896 static DMPolytopeType vertexT[] = {DM_POLYTOPE_POINT}; 897 static PetscInt vertexS[] = {1}; 898 static PetscInt vertexC[] = {0}; 899 static PetscInt vertexO[] = {0}; 900 static DMPolytopeType edgeT[] = {DM_POLYTOPE_SEGMENT}; 901 static PetscInt edgeS[] = {1}; 902 static PetscInt edgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 903 static PetscInt edgeO[] = {0, 0}; 904 static DMPolytopeType tedgeT[] = {DM_POLYTOPE_POINT_PRISM_TENSOR}; 905 static PetscInt tedgeS[] = {1}; 906 static PetscInt tedgeC[] = {DM_POLYTOPE_POINT, 1, 0, 0, DM_POLYTOPE_POINT, 1, 1, 0}; 907 static PetscInt tedgeO[] = {0, 0}; 908 static DMPolytopeType triT[] = {DM_POLYTOPE_TRIANGLE}; 909 static PetscInt triS[] = {1}; 910 static PetscInt triC[] = {DM_POLYTOPE_SEGMENT, 1, 0, 0, DM_POLYTOPE_SEGMENT, 1, 1, 0, DM_POLYTOPE_SEGMENT, 1, 2, 0}; 911 static PetscInt triO[] = {0, 0, 0}; 912 static DMPolytopeType quadT[] = {DM_POLYTOPE_QUADRILATERAL}; 913 static PetscInt quadS[] = {1}; 914 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}; 915 static PetscInt quadO[] = {0, 0, 0, 0}; 916 static DMPolytopeType tquadT[] = {DM_POLYTOPE_SEG_PRISM_TENSOR}; 917 static PetscInt tquadS[] = {1}; 918 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}; 919 static PetscInt tquadO[] = {0, 0, 0, 0}; 920 static DMPolytopeType tetT[] = {DM_POLYTOPE_TETRAHEDRON}; 921 static PetscInt tetS[] = {1}; 922 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}; 923 static PetscInt tetO[] = {0, 0, 0, 0}; 924 static DMPolytopeType hexT[] = {DM_POLYTOPE_HEXAHEDRON}; 925 static PetscInt hexS[] = {1}; 926 static PetscInt hexC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, 927 DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0, DM_POLYTOPE_QUADRILATERAL, 1, 5, 0}; 928 static PetscInt hexO[] = {0, 0, 0, 0, 0, 0}; 929 static DMPolytopeType tripT[] = {DM_POLYTOPE_TRI_PRISM}; 930 static PetscInt tripS[] = {1}; 931 static PetscInt tripC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, 932 DM_POLYTOPE_QUADRILATERAL, 1, 2, 0, DM_POLYTOPE_QUADRILATERAL, 1, 3, 0, DM_POLYTOPE_QUADRILATERAL, 1, 4, 0}; 933 static PetscInt tripO[] = {0, 0, 0, 0, 0}; 934 static DMPolytopeType ttripT[] = {DM_POLYTOPE_TRI_PRISM_TENSOR}; 935 static PetscInt ttripS[] = {1}; 936 static PetscInt ttripC[] = {DM_POLYTOPE_TRIANGLE, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, 937 DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 2, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 3, 0, DM_POLYTOPE_SEG_PRISM_TENSOR, 1, 4, 0}; 938 static PetscInt ttripO[] = {0, 0, 0, 0, 0}; 939 static DMPolytopeType tquadpT[] = {DM_POLYTOPE_QUAD_PRISM_TENSOR}; 940 static PetscInt tquadpS[] = {1}; 941 static PetscInt tquadpC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_QUADRILATERAL, 1, 1, 0, 942 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}; 943 static PetscInt tquadpO[] = {0, 0, 0, 0, 0, 0}; 944 static DMPolytopeType pyrT[] = {DM_POLYTOPE_PYRAMID}; 945 static PetscInt pyrS[] = {1}; 946 static PetscInt pyrC[] = {DM_POLYTOPE_QUADRILATERAL, 1, 0, 0, DM_POLYTOPE_TRIANGLE, 1, 1, 0, 947 DM_POLYTOPE_TRIANGLE, 1, 2, 0, DM_POLYTOPE_TRIANGLE, 1, 3, 0, DM_POLYTOPE_TRIANGLE, 1, 4, 0}; 948 static PetscInt pyrO[] = {0, 0, 0, 0, 0}; 949 950 PetscFunctionBegin; 951 if (rt) *rt = 0; 952 switch (source) { 953 case DM_POLYTOPE_POINT: *Nt = 1; *target = vertexT; *size = vertexS; *cone = vertexC; *ornt = vertexO; break; 954 case DM_POLYTOPE_SEGMENT: *Nt = 1; *target = edgeT; *size = edgeS; *cone = edgeC; *ornt = edgeO; break; 955 case DM_POLYTOPE_POINT_PRISM_TENSOR: *Nt = 1; *target = tedgeT; *size = tedgeS; *cone = tedgeC; *ornt = tedgeO; break; 956 case DM_POLYTOPE_TRIANGLE: *Nt = 1; *target = triT; *size = triS; *cone = triC; *ornt = triO; break; 957 case DM_POLYTOPE_QUADRILATERAL: *Nt = 1; *target = quadT; *size = quadS; *cone = quadC; *ornt = quadO; break; 958 case DM_POLYTOPE_SEG_PRISM_TENSOR: *Nt = 1; *target = tquadT; *size = tquadS; *cone = tquadC; *ornt = tquadO; break; 959 case DM_POLYTOPE_TETRAHEDRON: *Nt = 1; *target = tetT; *size = tetS; *cone = tetC; *ornt = tetO; break; 960 case DM_POLYTOPE_HEXAHEDRON: *Nt = 1; *target = hexT; *size = hexS; *cone = hexC; *ornt = hexO; break; 961 case DM_POLYTOPE_TRI_PRISM: *Nt = 1; *target = tripT; *size = tripS; *cone = tripC; *ornt = tripO; break; 962 case DM_POLYTOPE_TRI_PRISM_TENSOR: *Nt = 1; *target = ttripT; *size = ttripS; *cone = ttripC; *ornt = ttripO; break; 963 case DM_POLYTOPE_QUAD_PRISM_TENSOR: *Nt = 1; *target = tquadpT; *size = tquadpS; *cone = tquadpC; *ornt = tquadpO; break; 964 case DM_POLYTOPE_PYRAMID: *Nt = 1; *target = pyrT; *size = pyrS; *cone = pyrC; *ornt = pyrO; break; 965 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No refinement strategy for %s", DMPolytopeTypes[source]); 966 } 967 PetscFunctionReturn(0); 968 } 969 970 /*@ 971 DMPlexTransformGetSubcellOrientation - Transform the replica number and orientation for a target point according to the group action for the source point 972 973 Not collective 974 975 Input Parameters: 976 + tr - The DMPlexTransform 977 . sct - The source point cell type, from whom the new cell is being produced 978 . sp - The source point 979 . so - The orientation of the source point in its enclosing parent 980 . tct - The target point cell type 981 . r - The replica number requested for the produced cell type 982 - o - The orientation of the replica 983 984 Output Parameters: 985 + rnew - The replica number, given the orientation of the parent 986 - onew - The replica orientation, given the orientation of the parent 987 988 Level: advanced 989 990 .seealso: `DMPlexTransformCellTransform()`, `DMPlexTransformApply()` 991 @*/ 992 PetscErrorCode DMPlexTransformGetSubcellOrientation(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 993 { 994 PetscFunctionBeginHot; 995 PetscCall((*tr->ops->getsubcellorientation)(tr, sct, sp, so, tct, r, o, rnew, onew)); 996 PetscFunctionReturn(0); 997 } 998 999 static PetscErrorCode DMPlexTransformSetConeSizes(DMPlexTransform tr, DM rdm) 1000 { 1001 DM dm; 1002 PetscInt pStart, pEnd, p, pNew; 1003 1004 PetscFunctionBegin; 1005 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1006 /* Must create the celltype label here so that we do not automatically try to compute the types */ 1007 PetscCall(DMCreateLabel(rdm, "celltype")); 1008 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1009 for (p = pStart; p < pEnd; ++p) { 1010 DMPolytopeType ct; 1011 DMPolytopeType *rct; 1012 PetscInt *rsize, *rcone, *rornt; 1013 PetscInt Nct, n, r; 1014 1015 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1016 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1017 for (n = 0; n < Nct; ++n) { 1018 for (r = 0; r < rsize[n]; ++r) { 1019 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 1020 PetscCall(DMPlexSetConeSize(rdm, pNew, DMPolytopeTypeGetConeSize(rct[n]))); 1021 PetscCall(DMPlexSetCellType(rdm, pNew, rct[n])); 1022 } 1023 } 1024 } 1025 /* Let the DM know we have set all the cell types */ 1026 { 1027 DMLabel ctLabel; 1028 DM_Plex *plex = (DM_Plex *) rdm->data; 1029 1030 PetscCall(DMPlexGetCellTypeLabel(rdm, &ctLabel)); 1031 PetscCall(PetscObjectStateGet((PetscObject) ctLabel, &plex->celltypeState)); 1032 } 1033 PetscFunctionReturn(0); 1034 } 1035 1036 #if 0 1037 static PetscErrorCode DMPlexTransformGetConeSize(DMPlexTransform tr, PetscInt q, PetscInt *coneSize) 1038 { 1039 PetscInt ctNew; 1040 1041 PetscFunctionBegin; 1042 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1043 PetscValidPointer(coneSize, 3); 1044 /* TODO Can do bisection since everything is sorted */ 1045 for (ctNew = DM_POLYTOPE_POINT; ctNew < DM_NUM_POLYTOPES; ++ctNew) { 1046 PetscInt ctSN = tr->ctStartNew[ctNew], ctEN = tr->ctStartNew[tr->ctOrderNew[tr->ctOrderInvNew[ctNew]+1]]; 1047 1048 if (q >= ctSN && q < ctEN) break; 1049 } 1050 PetscCheck(ctNew < DM_NUM_POLYTOPES,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " cannot be located in the transformed mesh", q); 1051 *coneSize = DMPolytopeTypeGetConeSize((DMPolytopeType) ctNew); 1052 PetscFunctionReturn(0); 1053 } 1054 #endif 1055 1056 /* The orientation o is for the interior of the cell p */ 1057 static PetscErrorCode DMPlexTransformGetCone_Internal(DMPlexTransform tr, PetscInt p, PetscInt o, DMPolytopeType ct, DMPolytopeType ctNew, 1058 const PetscInt rcone[], PetscInt *coneoff, const PetscInt rornt[], PetscInt *orntoff, 1059 PetscInt coneNew[], PetscInt orntNew[]) 1060 { 1061 DM dm; 1062 const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); 1063 const PetscInt *cone; 1064 PetscInt c, coff = *coneoff, ooff = *orntoff; 1065 1066 PetscFunctionBegin; 1067 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1068 PetscCall(DMPlexGetCone(dm, p, &cone)); 1069 for (c = 0; c < csizeNew; ++c) { 1070 PetscInt ppp = -1; /* Parent Parent point: Parent of point pp */ 1071 PetscInt pp = p; /* Parent point: Point in the original mesh producing new cone point */ 1072 PetscInt po = 0; /* Orientation of parent point pp in parent parent point ppp */ 1073 DMPolytopeType pct = ct; /* Parent type: Cell type for parent of new cone point */ 1074 const PetscInt *pcone = cone; /* Parent cone: Cone of parent point pp */ 1075 PetscInt pr = -1; /* Replica number of pp that produces new cone point */ 1076 const DMPolytopeType ft = (DMPolytopeType) rcone[coff++]; /* Cell type for new cone point of pNew */ 1077 const PetscInt fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ 1078 PetscInt fo = rornt[ooff++]; /* Orientation of new cone point in pNew */ 1079 PetscInt lc; 1080 1081 /* Get the type (pct) and point number (pp) of the parent point in the original mesh which produces this cone point */ 1082 for (lc = 0; lc < fn; ++lc) { 1083 const PetscInt *parr = DMPolytopeTypeGetArrangment(pct, po); 1084 const PetscInt acp = rcone[coff++]; 1085 const PetscInt pcp = parr[acp*2]; 1086 const PetscInt pco = parr[acp*2+1]; 1087 const PetscInt *ppornt; 1088 1089 ppp = pp; 1090 pp = pcone[pcp]; 1091 PetscCall(DMPlexGetCellType(dm, pp, &pct)); 1092 PetscCall(DMPlexGetCone(dm, pp, &pcone)); 1093 PetscCall(DMPlexGetConeOrientation(dm, ppp, &ppornt)); 1094 po = DMPolytopeTypeComposeOrientation(pct, ppornt[pcp], pco); 1095 } 1096 pr = rcone[coff++]; 1097 /* Orientation po of pp maps (pr, fo) -> (pr', fo') */ 1098 PetscCall(DMPlexTransformGetSubcellOrientation(tr, pct, pp, fn ? po : o, ft, pr, fo, &pr, &fo)); 1099 PetscCall(DMPlexTransformGetTargetPoint(tr, pct, ft, pp, pr, &coneNew[c])); 1100 orntNew[c] = fo; 1101 } 1102 *coneoff = coff; 1103 *orntoff = ooff; 1104 PetscFunctionReturn(0); 1105 } 1106 1107 static PetscErrorCode DMPlexTransformSetCones(DMPlexTransform tr, DM rdm) 1108 { 1109 DM dm; 1110 DMPolytopeType ct; 1111 PetscInt *coneNew, *orntNew; 1112 PetscInt maxConeSize = 0, pStart, pEnd, p, pNew; 1113 1114 PetscFunctionBegin; 1115 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1116 for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p)); 1117 PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew)); 1118 PetscCall(DMGetWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew)); 1119 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1120 for (p = pStart; p < pEnd; ++p) { 1121 const PetscInt *cone, *ornt; 1122 PetscInt coff, ooff; 1123 DMPolytopeType *rct; 1124 PetscInt *rsize, *rcone, *rornt; 1125 PetscInt Nct, n, r; 1126 1127 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1128 PetscCall(DMPlexGetCone(dm, p, &cone)); 1129 PetscCall(DMPlexGetConeOrientation(dm, p, &ornt)); 1130 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1131 for (n = 0, coff = 0, ooff = 0; n < Nct; ++n) { 1132 const DMPolytopeType ctNew = rct[n]; 1133 1134 for (r = 0; r < rsize[n]; ++r) { 1135 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 1136 PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, ctNew, rcone, &coff, rornt, &ooff, coneNew, orntNew)); 1137 PetscCall(DMPlexSetCone(rdm, pNew, coneNew)); 1138 PetscCall(DMPlexSetConeOrientation(rdm, pNew, orntNew)); 1139 } 1140 } 1141 } 1142 PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &coneNew)); 1143 PetscCall(DMRestoreWorkArray(rdm, maxConeSize, MPIU_INT, &orntNew)); 1144 PetscCall(DMViewFromOptions(rdm, NULL, "-rdm_view")); 1145 PetscCall(DMPlexSymmetrize(rdm)); 1146 PetscCall(DMPlexStratify(rdm)); 1147 PetscFunctionReturn(0); 1148 } 1149 1150 PetscErrorCode DMPlexTransformGetConeOriented(DMPlexTransform tr, PetscInt q, PetscInt po, const PetscInt *cone[], const PetscInt *ornt[]) 1151 { 1152 DM dm; 1153 DMPolytopeType ct, qct; 1154 DMPolytopeType *rct; 1155 PetscInt *rsize, *rcone, *rornt, *qcone, *qornt; 1156 PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0; 1157 1158 PetscFunctionBegin; 1159 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1160 PetscValidPointer(cone, 4); 1161 PetscValidPointer(ornt, 5); 1162 for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p)); 1163 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1164 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone)); 1165 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt)); 1166 PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r)); 1167 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1168 for (n = 0; n < Nct; ++n) { 1169 const DMPolytopeType ctNew = rct[n]; 1170 const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); 1171 PetscInt Nr = rsize[n], fn, c; 1172 1173 if (ctNew == qct) Nr = r; 1174 for (nr = 0; nr < Nr; ++nr) { 1175 for (c = 0; c < csizeNew; ++c) { 1176 ++coff; /* Cell type of new cone point */ 1177 fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ 1178 coff += fn; 1179 ++coff; /* Replica number of new cone point */ 1180 ++ooff; /* Orientation of new cone point */ 1181 } 1182 } 1183 if (ctNew == qct) break; 1184 } 1185 PetscCall(DMPlexTransformGetCone_Internal(tr, p, po, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt)); 1186 *cone = qcone; 1187 *ornt = qornt; 1188 PetscFunctionReturn(0); 1189 } 1190 1191 PetscErrorCode DMPlexTransformGetCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[]) 1192 { 1193 DM dm; 1194 DMPolytopeType ct, qct; 1195 DMPolytopeType *rct; 1196 PetscInt *rsize, *rcone, *rornt, *qcone, *qornt; 1197 PetscInt maxConeSize = 0, Nct, p, r, n, nr, coff = 0, ooff = 0; 1198 1199 PetscFunctionBegin; 1200 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1201 PetscValidPointer(cone, 3); 1202 for (p = 0; p < DM_NUM_POLYTOPES; ++p) maxConeSize = PetscMax(maxConeSize, DMPolytopeTypeGetConeSize((DMPolytopeType) p)); 1203 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1204 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qcone)); 1205 PetscCall(DMGetWorkArray(dm, maxConeSize, MPIU_INT, &qornt)); 1206 PetscCall(DMPlexTransformGetSourcePoint(tr, q, &ct, &qct, &p, &r)); 1207 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1208 for (n = 0; n < Nct; ++n) { 1209 const DMPolytopeType ctNew = rct[n]; 1210 const PetscInt csizeNew = DMPolytopeTypeGetConeSize(ctNew); 1211 PetscInt Nr = rsize[n], fn, c; 1212 1213 if (ctNew == qct) Nr = r; 1214 for (nr = 0; nr < Nr; ++nr) { 1215 for (c = 0; c < csizeNew; ++c) { 1216 ++coff; /* Cell type of new cone point */ 1217 fn = rcone[coff++]; /* Number of cones of p that need to be taken when producing new cone point */ 1218 coff += fn; 1219 ++coff; /* Replica number of new cone point */ 1220 ++ooff; /* Orientation of new cone point */ 1221 } 1222 } 1223 if (ctNew == qct) break; 1224 } 1225 PetscCall(DMPlexTransformGetCone_Internal(tr, p, 0, ct, qct, rcone, &coff, rornt, &ooff, qcone, qornt)); 1226 *cone = qcone; 1227 *ornt = qornt; 1228 PetscFunctionReturn(0); 1229 } 1230 1231 PetscErrorCode DMPlexTransformRestoreCone(DMPlexTransform tr, PetscInt q, const PetscInt *cone[], const PetscInt *ornt[]) 1232 { 1233 DM dm; 1234 1235 PetscFunctionBegin; 1236 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1237 PetscValidPointer(cone, 3); 1238 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1239 PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, cone)); 1240 PetscCall(DMRestoreWorkArray(dm, 0, MPIU_INT, ornt)); 1241 PetscFunctionReturn(0); 1242 } 1243 1244 static PetscErrorCode DMPlexTransformCreateCellVertices_Internal(DMPlexTransform tr) 1245 { 1246 PetscInt ict; 1247 1248 PetscFunctionBegin; 1249 PetscCall(PetscCalloc3(DM_NUM_POLYTOPES, &tr->trNv, DM_NUM_POLYTOPES, &tr->trVerts, DM_NUM_POLYTOPES, &tr->trSubVerts)); 1250 for (ict = DM_POLYTOPE_POINT; ict < DM_NUM_POLYTOPES; ++ict) { 1251 const DMPolytopeType ct = (DMPolytopeType) ict; 1252 DMPlexTransform reftr; 1253 DM refdm, trdm; 1254 Vec coordinates; 1255 const PetscScalar *coords; 1256 DMPolytopeType *rct; 1257 PetscInt *rsize, *rcone, *rornt; 1258 PetscInt Nct, n, r, pNew; 1259 PetscInt vStart, vEnd, Nc; 1260 const PetscInt debug = 0; 1261 const char *typeName; 1262 1263 /* Since points are 0-dimensional, coordinates make no sense */ 1264 if (DMPolytopeTypeGetDim(ct) <= 0) continue; 1265 PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &refdm)); 1266 PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &reftr)); 1267 PetscCall(DMPlexTransformSetDM(reftr, refdm)); 1268 PetscCall(DMPlexTransformGetType(tr, &typeName)); 1269 PetscCall(DMPlexTransformSetType(reftr, typeName)); 1270 PetscCall(DMPlexTransformSetUp(reftr)); 1271 PetscCall(DMPlexTransformApply(reftr, refdm, &trdm)); 1272 1273 PetscCall(DMPlexGetDepthStratum(trdm, 0, &vStart, &vEnd)); 1274 tr->trNv[ct] = vEnd - vStart; 1275 PetscCall(DMGetCoordinatesLocal(trdm, &coordinates)); 1276 PetscCall(VecGetLocalSize(coordinates, &Nc)); 1277 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); 1278 PetscCall(PetscCalloc1(Nc, &tr->trVerts[ct])); 1279 PetscCall(VecGetArrayRead(coordinates, &coords)); 1280 PetscCall(PetscArraycpy(tr->trVerts[ct], coords, Nc)); 1281 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 1282 1283 PetscCall(PetscCalloc1(DM_NUM_POLYTOPES, &tr->trSubVerts[ct])); 1284 PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1285 for (n = 0; n < Nct; ++n) { 1286 1287 /* Since points are 0-dimensional, coordinates make no sense */ 1288 if (rct[n] == DM_POLYTOPE_POINT) continue; 1289 PetscCall(PetscCalloc1(rsize[n], &tr->trSubVerts[ct][rct[n]])); 1290 for (r = 0; r < rsize[n]; ++r) { 1291 PetscInt *closure = NULL; 1292 PetscInt clSize, cl, Nv = 0; 1293 1294 PetscCall(PetscCalloc1(DMPolytopeTypeGetNumVertices(rct[n]), &tr->trSubVerts[ct][rct[n]][r])); 1295 PetscCall(DMPlexTransformGetTargetPoint(reftr, ct, rct[n], 0, r, &pNew)); 1296 PetscCall(DMPlexGetTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure)); 1297 for (cl = 0; cl < clSize*2; cl += 2) { 1298 const PetscInt sv = closure[cl]; 1299 1300 if ((sv >= vStart) && (sv < vEnd)) tr->trSubVerts[ct][rct[n]][r][Nv++] = sv - vStart; 1301 } 1302 PetscCall(DMPlexRestoreTransitiveClosure(trdm, pNew, PETSC_TRUE, &clSize, &closure)); 1303 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]); 1304 } 1305 } 1306 if (debug) { 1307 DMPolytopeType *rct; 1308 PetscInt *rsize, *rcone, *rornt; 1309 PetscInt v, dE = DMPolytopeTypeGetDim(ct), d, off = 0; 1310 1311 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %" PetscInt_FMT " vertices\n", DMPolytopeTypes[ct], tr->trNv[ct])); 1312 for (v = 0; v < tr->trNv[ct]; ++v) { 1313 PetscCall(PetscPrintf(PETSC_COMM_SELF, " ")); 1314 for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g ", (double)PetscRealPart(tr->trVerts[ct][off++]))); 1315 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 1316 } 1317 1318 PetscCall(DMPlexTransformCellTransform(reftr, ct, 0, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1319 for (n = 0; n < Nct; ++n) { 1320 if (rct[n] == DM_POLYTOPE_POINT) continue; 1321 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%s: %s subvertices %" PetscInt_FMT "\n", DMPolytopeTypes[ct], DMPolytopeTypes[rct[n]], tr->trNv[ct])); 1322 for (r = 0; r < rsize[n]; ++r) { 1323 PetscCall(PetscPrintf(PETSC_COMM_SELF, " ")); 1324 for (v = 0; v < DMPolytopeTypeGetNumVertices(rct[n]); ++v) { 1325 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " ", tr->trSubVerts[ct][rct[n]][r][v])); 1326 } 1327 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 1328 } 1329 } 1330 } 1331 PetscCall(DMDestroy(&refdm)); 1332 PetscCall(DMDestroy(&trdm)); 1333 PetscCall(DMPlexTransformDestroy(&reftr)); 1334 } 1335 PetscFunctionReturn(0); 1336 } 1337 1338 /* 1339 DMPlexTransformGetCellVertices - Get the set of transformed vertices lying in the closure of a reference cell of given type 1340 1341 Input Parameters: 1342 + tr - The DMPlexTransform object 1343 - ct - The cell type 1344 1345 Output Parameters: 1346 + Nv - The number of transformed vertices in the closure of the reference cell of given type 1347 - trVerts - The coordinates of these vertices in the reference cell 1348 1349 Level: developer 1350 1351 .seealso: `DMPlexTransformGetSubcellVertices()` 1352 */ 1353 PetscErrorCode DMPlexTransformGetCellVertices(DMPlexTransform tr, DMPolytopeType ct, PetscInt *Nv, PetscScalar *trVerts[]) 1354 { 1355 PetscFunctionBegin; 1356 if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr)); 1357 if (Nv) *Nv = tr->trNv[ct]; 1358 if (trVerts) *trVerts = tr->trVerts[ct]; 1359 PetscFunctionReturn(0); 1360 } 1361 1362 /* 1363 DMPlexTransformGetSubcellVertices - Get the set of transformed vertices defining a subcell in the reference cell of given type 1364 1365 Input Parameters: 1366 + tr - The DMPlexTransform object 1367 . ct - The cell type 1368 . rct - The subcell type 1369 - r - The subcell index 1370 1371 Output Parameter: 1372 . subVerts - The indices of these vertices in the set of vertices returned by DMPlexTransformGetCellVertices() 1373 1374 Level: developer 1375 1376 .seealso: `DMPlexTransformGetCellVertices()` 1377 */ 1378 PetscErrorCode DMPlexTransformGetSubcellVertices(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, PetscInt *subVerts[]) 1379 { 1380 PetscFunctionBegin; 1381 if (!tr->trNv) PetscCall(DMPlexTransformCreateCellVertices_Internal(tr)); 1382 PetscCheck(tr->trSubVerts[ct][rct],PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_WRONG, "Cell type %s does not produce %s", DMPolytopeTypes[ct], DMPolytopeTypes[rct]); 1383 if (subVerts) *subVerts = tr->trSubVerts[ct][rct][r]; 1384 PetscFunctionReturn(0); 1385 } 1386 1387 /* Computes new vertex as the barycenter, or centroid */ 1388 PetscErrorCode DMPlexTransformMapCoordinatesBarycenter_Internal(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) 1389 { 1390 PetscInt v,d; 1391 1392 PetscFunctionBeginHot; 1393 PetscCheck(ct == DM_POLYTOPE_POINT,PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for refined point type %s",DMPolytopeTypes[ct]); 1394 for (d = 0; d < dE; ++d) out[d] = 0.0; 1395 for (v = 0; v < Nv; ++v) for (d = 0; d < dE; ++d) out[d] += in[v*dE+d]; 1396 for (d = 0; d < dE; ++d) out[d] /= Nv; 1397 PetscFunctionReturn(0); 1398 } 1399 1400 /*@ 1401 DMPlexTransformMapCoordinates - Calculate new coordinates for produced points 1402 1403 Not collective 1404 1405 Input Parameters: 1406 + cr - The DMPlexCellRefiner 1407 . pct - The cell type of the parent, from whom the new cell is being produced 1408 . ct - The type being produced 1409 . p - The original point 1410 . r - The replica number requested for the produced cell type 1411 . Nv - Number of vertices in the closure of the parent cell 1412 . dE - Spatial dimension 1413 - in - array of size Nv*dE, holding coordinates of the vertices in the closure of the parent cell 1414 1415 Output Parameter: 1416 . out - The coordinates of the new vertices 1417 1418 Level: intermediate 1419 1420 .seealso: `DMPlexTransformApply()` 1421 @*/ 1422 PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) 1423 { 1424 PetscFunctionBeginHot; 1425 PetscCall((*tr->ops->mapcoordinates)(tr, pct, ct, p, r, Nv, dE, in, out)); 1426 PetscFunctionReturn(0); 1427 } 1428 1429 static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew) 1430 { 1431 DM dm; 1432 IS valueIS; 1433 const PetscInt *values; 1434 PetscInt defVal, Nv, val; 1435 1436 PetscFunctionBegin; 1437 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1438 PetscCall(DMLabelGetDefaultValue(label, &defVal)); 1439 PetscCall(DMLabelSetDefaultValue(labelNew, defVal)); 1440 PetscCall(DMLabelGetValueIS(label, &valueIS)); 1441 PetscCall(ISGetLocalSize(valueIS, &Nv)); 1442 PetscCall(ISGetIndices(valueIS, &values)); 1443 for (val = 0; val < Nv; ++val) { 1444 IS pointIS; 1445 const PetscInt *points; 1446 PetscInt numPoints, p; 1447 1448 /* Ensure refined label is created with same number of strata as 1449 * original (even if no entries here). */ 1450 PetscCall(DMLabelAddStratum(labelNew, values[val])); 1451 PetscCall(DMLabelGetStratumIS(label, values[val], &pointIS)); 1452 PetscCall(ISGetLocalSize(pointIS, &numPoints)); 1453 PetscCall(ISGetIndices(pointIS, &points)); 1454 for (p = 0; p < numPoints; ++p) { 1455 const PetscInt point = points[p]; 1456 DMPolytopeType ct; 1457 DMPolytopeType *rct; 1458 PetscInt *rsize, *rcone, *rornt; 1459 PetscInt Nct, n, r, pNew=0; 1460 1461 PetscCall(DMPlexGetCellType(dm, point, &ct)); 1462 PetscCall(DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1463 for (n = 0; n < Nct; ++n) { 1464 for (r = 0; r < rsize[n]; ++r) { 1465 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew)); 1466 PetscCall(DMLabelSetValue(labelNew, pNew, values[val])); 1467 } 1468 } 1469 } 1470 PetscCall(ISRestoreIndices(pointIS, &points)); 1471 PetscCall(ISDestroy(&pointIS)); 1472 } 1473 PetscCall(ISRestoreIndices(valueIS, &values)); 1474 PetscCall(ISDestroy(&valueIS)); 1475 PetscFunctionReturn(0); 1476 } 1477 1478 static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm) 1479 { 1480 DM dm; 1481 PetscInt numLabels, l; 1482 1483 PetscFunctionBegin; 1484 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1485 PetscCall(DMGetNumLabels(dm, &numLabels)); 1486 for (l = 0; l < numLabels; ++l) { 1487 DMLabel label, labelNew; 1488 const char *lname; 1489 PetscBool isDepth, isCellType; 1490 1491 PetscCall(DMGetLabelName(dm, l, &lname)); 1492 PetscCall(PetscStrcmp(lname, "depth", &isDepth)); 1493 if (isDepth) continue; 1494 PetscCall(PetscStrcmp(lname, "celltype", &isCellType)); 1495 if (isCellType) continue; 1496 PetscCall(DMCreateLabel(rdm, lname)); 1497 PetscCall(DMGetLabel(dm, lname, &label)); 1498 PetscCall(DMGetLabel(rdm, lname, &labelNew)); 1499 PetscCall(RefineLabel_Internal(tr, label, labelNew)); 1500 } 1501 PetscFunctionReturn(0); 1502 } 1503 1504 /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */ 1505 PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm) 1506 { 1507 DM dm; 1508 PetscInt Nf, f, Nds, s; 1509 1510 PetscFunctionBegin; 1511 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1512 PetscCall(DMGetNumFields(dm, &Nf)); 1513 for (f = 0; f < Nf; ++f) { 1514 DMLabel label, labelNew; 1515 PetscObject obj; 1516 const char *lname; 1517 1518 PetscCall(DMGetField(rdm, f, &label, &obj)); 1519 if (!label) continue; 1520 PetscCall(PetscObjectGetName((PetscObject) label, &lname)); 1521 PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew)); 1522 PetscCall(RefineLabel_Internal(tr, label, labelNew)); 1523 PetscCall(DMSetField_Internal(rdm, f, labelNew, obj)); 1524 PetscCall(DMLabelDestroy(&labelNew)); 1525 } 1526 PetscCall(DMGetNumDS(dm, &Nds)); 1527 for (s = 0; s < Nds; ++s) { 1528 DMLabel label, labelNew; 1529 const char *lname; 1530 1531 PetscCall(DMGetRegionNumDS(rdm, s, &label, NULL, NULL)); 1532 if (!label) continue; 1533 PetscCall(PetscObjectGetName((PetscObject) label, &lname)); 1534 PetscCall(DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew)); 1535 PetscCall(RefineLabel_Internal(tr, label, labelNew)); 1536 PetscCall(DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL)); 1537 PetscCall(DMLabelDestroy(&labelNew)); 1538 } 1539 PetscFunctionReturn(0); 1540 } 1541 1542 static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm) 1543 { 1544 DM dm; 1545 PetscSF sf, sfNew; 1546 PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m; 1547 const PetscInt *localPoints; 1548 const PetscSFNode *remotePoints; 1549 PetscInt *localPointsNew; 1550 PetscSFNode *remotePointsNew; 1551 PetscInt pStartNew, pEndNew, pNew; 1552 /* Brute force algorithm */ 1553 PetscSF rsf; 1554 PetscSection s; 1555 const PetscInt *rootdegree; 1556 PetscInt *rootPointsNew, *remoteOffsets; 1557 PetscInt numPointsNew, pStart, pEnd, p; 1558 1559 PetscFunctionBegin; 1560 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1561 PetscCall(DMPlexGetChart(rdm, &pStartNew, &pEndNew)); 1562 PetscCall(DMGetPointSF(dm, &sf)); 1563 PetscCall(DMGetPointSF(rdm, &sfNew)); 1564 /* Calculate size of new SF */ 1565 PetscCall(PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints)); 1566 if (numRoots < 0) PetscFunctionReturn(0); 1567 for (l = 0; l < numLeaves; ++l) { 1568 const PetscInt p = localPoints[l]; 1569 DMPolytopeType ct; 1570 DMPolytopeType *rct; 1571 PetscInt *rsize, *rcone, *rornt; 1572 PetscInt Nct, n; 1573 1574 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1575 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1576 for (n = 0; n < Nct; ++n) { 1577 numLeavesNew += rsize[n]; 1578 } 1579 } 1580 /* Send new root point numbers 1581 It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication 1582 */ 1583 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1584 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject) dm), &s)); 1585 PetscCall(PetscSectionSetChart(s, pStart, pEnd)); 1586 for (p = pStart; p < pEnd; ++p) { 1587 DMPolytopeType ct; 1588 DMPolytopeType *rct; 1589 PetscInt *rsize, *rcone, *rornt; 1590 PetscInt Nct, n; 1591 1592 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1593 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1594 for (n = 0; n < Nct; ++n) { 1595 PetscCall(PetscSectionAddDof(s, p, rsize[n])); 1596 } 1597 } 1598 PetscCall(PetscSectionSetUp(s)); 1599 PetscCall(PetscSectionGetStorageSize(s, &numPointsNew)); 1600 PetscCall(PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets)); 1601 PetscCall(PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf)); 1602 PetscCall(PetscFree(remoteOffsets)); 1603 PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree)); 1604 PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree)); 1605 PetscCall(PetscMalloc1(numPointsNew, &rootPointsNew)); 1606 for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1; 1607 for (p = pStart; p < pEnd; ++p) { 1608 DMPolytopeType ct; 1609 DMPolytopeType *rct; 1610 PetscInt *rsize, *rcone, *rornt; 1611 PetscInt Nct, n, r, off; 1612 1613 if (!rootdegree[p-pStart]) continue; 1614 PetscCall(PetscSectionGetOffset(s, p, &off)); 1615 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1616 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1617 for (n = 0, m = 0; n < Nct; ++n) { 1618 for (r = 0; r < rsize[n]; ++r, ++m) { 1619 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 1620 rootPointsNew[off+m] = pNew; 1621 } 1622 } 1623 } 1624 PetscCall(PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew,MPI_REPLACE)); 1625 PetscCall(PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew,MPI_REPLACE)); 1626 PetscCall(PetscSFDestroy(&rsf)); 1627 PetscCall(PetscMalloc1(numLeavesNew, &localPointsNew)); 1628 PetscCall(PetscMalloc1(numLeavesNew, &remotePointsNew)); 1629 for (l = 0, m = 0; l < numLeaves; ++l) { 1630 const PetscInt p = localPoints[l]; 1631 DMPolytopeType ct; 1632 DMPolytopeType *rct; 1633 PetscInt *rsize, *rcone, *rornt; 1634 PetscInt Nct, n, r, q, off; 1635 1636 PetscCall(PetscSectionGetOffset(s, p, &off)); 1637 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1638 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1639 for (n = 0, q = 0; n < Nct; ++n) { 1640 for (r = 0; r < rsize[n]; ++r, ++m, ++q) { 1641 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 1642 localPointsNew[m] = pNew; 1643 remotePointsNew[m].index = rootPointsNew[off+q]; 1644 remotePointsNew[m].rank = remotePoints[l].rank; 1645 } 1646 } 1647 } 1648 PetscCall(PetscSectionDestroy(&s)); 1649 PetscCall(PetscFree(rootPointsNew)); 1650 /* SF needs sorted leaves to correctly calculate Gather */ 1651 { 1652 PetscSFNode *rp, *rtmp; 1653 PetscInt *lp, *idx, *ltmp, i; 1654 1655 PetscCall(PetscMalloc1(numLeavesNew, &idx)); 1656 PetscCall(PetscMalloc1(numLeavesNew, &lp)); 1657 PetscCall(PetscMalloc1(numLeavesNew, &rp)); 1658 for (i = 0; i < numLeavesNew; ++i) { 1659 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); 1660 idx[i] = i; 1661 } 1662 PetscCall(PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx)); 1663 for (i = 0; i < numLeavesNew; ++i) { 1664 lp[i] = localPointsNew[idx[i]]; 1665 rp[i] = remotePointsNew[idx[i]]; 1666 } 1667 ltmp = localPointsNew; 1668 localPointsNew = lp; 1669 rtmp = remotePointsNew; 1670 remotePointsNew = rp; 1671 PetscCall(PetscFree(idx)); 1672 PetscCall(PetscFree(ltmp)); 1673 PetscCall(PetscFree(rtmp)); 1674 } 1675 PetscCall(PetscSFSetGraph(sfNew, pEndNew-pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER)); 1676 PetscFunctionReturn(0); 1677 } 1678 1679 /* 1680 DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of type ct with localized coordinates x, we generate localized coordinates xr for subcell r of type rct. 1681 1682 Not collective 1683 1684 Input Parameters: 1685 + cr - The DMPlexCellRefiner 1686 . ct - The type of the parent cell 1687 . rct - The type of the produced cell 1688 . r - The index of the produced cell 1689 - x - The localized coordinates for the parent cell 1690 1691 Output Parameter: 1692 . xr - The localized coordinates for the produced cell 1693 1694 Level: developer 1695 1696 .seealso: `DMPlexCellRefinerSetCoordinates()` 1697 */ 1698 static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[]) 1699 { 1700 PetscFE fe = NULL; 1701 PetscInt cdim, v, *subcellV; 1702 1703 PetscFunctionBegin; 1704 PetscCall(DMPlexTransformGetCoordinateFE(tr, ct, &fe)); 1705 PetscCall(DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV)); 1706 PetscCall(PetscFEGetNumComponents(fe, &cdim)); 1707 for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) { 1708 PetscCall(PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v*cdim])); 1709 } 1710 PetscFunctionReturn(0); 1711 } 1712 1713 static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm) 1714 { 1715 DM dm, cdm; 1716 PetscSection coordSection, coordSectionNew; 1717 Vec coordsLocal, coordsLocalNew; 1718 const PetscScalar *coords; 1719 PetscScalar *coordsNew; 1720 const DMBoundaryType *bd; 1721 const PetscReal *maxCell, *L; 1722 PetscBool isperiodic, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE; 1723 PetscInt dE, dEo, d, cStart, cEnd, c, vStartNew, vEndNew, v, pStart, pEnd, p, ocStart, ocEnd; 1724 1725 PetscFunctionBegin; 1726 PetscCall(DMPlexTransformGetDM(tr, &dm)); 1727 PetscCall(DMGetCoordinateDM(dm, &cdm)); 1728 PetscCall(DMGetPeriodicity(dm, &isperiodic, &maxCell, &L, &bd)); 1729 /* Determine if we need to localize coordinates when generating them */ 1730 if (isperiodic) { 1731 localizeVertices = PETSC_TRUE; 1732 if (!maxCell) { 1733 PetscBool localized; 1734 PetscCall(DMGetCoordinatesLocalized(dm, &localized)); 1735 PetscCheck(localized,PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Cannot refine a periodic mesh if coordinates have not been localized"); 1736 localizeCells = PETSC_TRUE; 1737 } 1738 } 1739 1740 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 1741 PetscCall(PetscSectionGetFieldComponents(coordSection, 0, &dEo)); 1742 if (maxCell) { 1743 PetscReal maxCellNew[3]; 1744 1745 for (d = 0; d < dEo; ++d) maxCellNew[d] = maxCell[d]/2.0; 1746 PetscCall(DMSetPeriodicity(rdm, isperiodic, maxCellNew, L, bd)); 1747 } else { 1748 PetscCall(DMSetPeriodicity(rdm, isperiodic, maxCell, L, bd)); 1749 } 1750 PetscCall(DMGetCoordinateDim(rdm, &dE)); 1751 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject) rdm), &coordSectionNew)); 1752 PetscCall(PetscSectionSetNumFields(coordSectionNew, 1)); 1753 PetscCall(PetscSectionSetFieldComponents(coordSectionNew, 0, dE)); 1754 PetscCall(DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew)); 1755 if (localizeCells) PetscCall(PetscSectionSetChart(coordSectionNew, 0, vEndNew)); 1756 else PetscCall(PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew)); 1757 1758 /* Localization should be inherited */ 1759 /* Stefano calculates parent cells for each new cell for localization */ 1760 /* Localized cells need coordinates of closure */ 1761 for (v = vStartNew; v < vEndNew; ++v) { 1762 PetscCall(PetscSectionSetDof(coordSectionNew, v, dE)); 1763 PetscCall(PetscSectionSetFieldDof(coordSectionNew, v, 0, dE)); 1764 } 1765 if (localizeCells) { 1766 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 1767 for (c = cStart; c < cEnd; ++c) { 1768 PetscInt dof; 1769 1770 PetscCall(PetscSectionGetDof(coordSection, c, &dof)); 1771 if (dof) { 1772 DMPolytopeType ct; 1773 DMPolytopeType *rct; 1774 PetscInt *rsize, *rcone, *rornt; 1775 PetscInt dim, cNew, Nct, n, r; 1776 1777 PetscCall(DMPlexGetCellType(dm, c, &ct)); 1778 dim = DMPolytopeTypeGetDim(ct); 1779 PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1780 /* This allows for different cell types */ 1781 for (n = 0; n < Nct; ++n) { 1782 if (dim != DMPolytopeTypeGetDim(rct[n])) continue; 1783 for (r = 0; r < rsize[n]; ++r) { 1784 PetscInt *closure = NULL; 1785 PetscInt clSize, cl, Nv = 0; 1786 1787 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew)); 1788 PetscCall(DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure)); 1789 for (cl = 0; cl < clSize*2; cl += 2) {if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;} 1790 PetscCall(DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure)); 1791 PetscCall(PetscSectionSetDof(coordSectionNew, cNew, Nv * dE)); 1792 PetscCall(PetscSectionSetFieldDof(coordSectionNew, cNew, 0, Nv * dE)); 1793 } 1794 } 1795 } 1796 } 1797 } 1798 PetscCall(PetscSectionSetUp(coordSectionNew)); 1799 PetscCall(DMViewFromOptions(dm, NULL, "-coarse_dm_view")); 1800 PetscCall(DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew)); 1801 { 1802 VecType vtype; 1803 PetscInt coordSizeNew, bs; 1804 const char *name; 1805 1806 PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal)); 1807 PetscCall(VecCreate(PETSC_COMM_SELF, &coordsLocalNew)); 1808 PetscCall(PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew)); 1809 PetscCall(VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE)); 1810 PetscCall(PetscObjectGetName((PetscObject) coordsLocal, &name)); 1811 PetscCall(PetscObjectSetName((PetscObject) coordsLocalNew, name)); 1812 PetscCall(VecGetBlockSize(coordsLocal, &bs)); 1813 PetscCall(VecSetBlockSize(coordsLocalNew, dEo == dE ? bs : dE)); 1814 PetscCall(VecGetType(coordsLocal, &vtype)); 1815 PetscCall(VecSetType(coordsLocalNew, vtype)); 1816 } 1817 PetscCall(VecGetArrayRead(coordsLocal, &coords)); 1818 PetscCall(VecGetArray(coordsLocalNew, &coordsNew)); 1819 PetscCall(PetscSectionGetChart(coordSection, &ocStart, &ocEnd)); 1820 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1821 /* First set coordinates for vertices*/ 1822 for (p = pStart; p < pEnd; ++p) { 1823 DMPolytopeType ct; 1824 DMPolytopeType *rct; 1825 PetscInt *rsize, *rcone, *rornt; 1826 PetscInt Nct, n, r; 1827 PetscBool hasVertex = PETSC_FALSE, isLocalized = PETSC_FALSE; 1828 1829 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1830 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1831 for (n = 0; n < Nct; ++n) { 1832 if (rct[n] == DM_POLYTOPE_POINT) {hasVertex = PETSC_TRUE; break;} 1833 } 1834 if (localizeVertices && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) { 1835 PetscInt dof; 1836 PetscCall(PetscSectionGetDof(coordSection, p, &dof)); 1837 if (dof) isLocalized = PETSC_TRUE; 1838 } 1839 if (hasVertex) { 1840 const PetscScalar *icoords = NULL; 1841 PetscScalar *pcoords = NULL; 1842 PetscInt Nc, Nv, v, d; 1843 1844 PetscCall(DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords)); 1845 1846 icoords = pcoords; 1847 Nv = Nc/dEo; 1848 if (ct != DM_POLYTOPE_POINT) { 1849 if (localizeVertices) { 1850 PetscScalar anchor[3]; 1851 1852 for (d = 0; d < dEo; ++d) anchor[d] = pcoords[d]; 1853 if (!isLocalized) { 1854 for (v = 0; v < Nv; ++v) PetscCall(DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v*dEo], &pcoords[v*dEo])); 1855 } else { 1856 Nv = Nc/(2*dEo); 1857 for (v = Nv; v < Nv*2; ++v) PetscCall(DMLocalizeCoordinate_Internal(dm, dEo, anchor, &pcoords[v*dEo], &pcoords[v*dEo])); 1858 } 1859 } 1860 } 1861 for (n = 0; n < Nct; ++n) { 1862 if (rct[n] != DM_POLYTOPE_POINT) continue; 1863 for (r = 0; r < rsize[n]; ++r) { 1864 PetscScalar vcoords[3]; 1865 PetscInt vNew, off; 1866 1867 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew)); 1868 PetscCall(PetscSectionGetOffset(coordSectionNew, vNew, &off)); 1869 PetscCall(DMPlexTransformMapCoordinates(tr, ct, rct[n], p, r, Nv, dEo, icoords, vcoords)); 1870 PetscCall(DMPlexSnapToGeomModel(dm, p, dE, vcoords, &coordsNew[off])); 1871 } 1872 } 1873 PetscCall(DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords)); 1874 } 1875 } 1876 /* Then set coordinates for cells by localizing */ 1877 for (p = pStart; p < pEnd; ++p) { 1878 DMPolytopeType ct; 1879 DMPolytopeType *rct; 1880 PetscInt *rsize, *rcone, *rornt; 1881 PetscInt Nct, n, r; 1882 PetscBool isLocalized = PETSC_FALSE; 1883 1884 PetscCall(DMPlexGetCellType(dm, p, &ct)); 1885 PetscCall(DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt)); 1886 if (localizeCells && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) { 1887 PetscInt dof; 1888 PetscCall(PetscSectionGetDof(coordSection, p, &dof)); 1889 if (dof) isLocalized = PETSC_TRUE; 1890 } 1891 if (isLocalized) { 1892 const PetscScalar *pcoords; 1893 1894 PetscCall(DMPlexPointLocalRead(cdm, p, coords, &pcoords)); 1895 for (n = 0; n < Nct; ++n) { 1896 const PetscInt Nr = rsize[n]; 1897 1898 if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue; 1899 for (r = 0; r < Nr; ++r) { 1900 PetscInt pNew, offNew; 1901 1902 /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that 1903 DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger 1904 cell to the ones it produces. */ 1905 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew)); 1906 PetscCall(PetscSectionGetOffset(coordSectionNew, pNew, &offNew)); 1907 PetscCall(DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew])); 1908 } 1909 } 1910 } 1911 } 1912 PetscCall(VecRestoreArrayRead(coordsLocal, &coords)); 1913 PetscCall(VecRestoreArray(coordsLocalNew, &coordsNew)); 1914 PetscCall(DMSetCoordinatesLocal(rdm, coordsLocalNew)); 1915 /* TODO Stefano has a final reduction if some hybrid coordinates cannot be found. (needcoords) Should not be needed. */ 1916 PetscCall(VecDestroy(&coordsLocalNew)); 1917 PetscCall(PetscSectionDestroy(&coordSectionNew)); 1918 if (!localizeCells) PetscCall(DMLocalizeCoordinates(rdm)); 1919 PetscFunctionReturn(0); 1920 } 1921 1922 PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm) 1923 { 1924 DM rdm; 1925 DMPlexInterpolatedFlag interp; 1926 1927 PetscFunctionBegin; 1928 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1929 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1930 PetscValidPointer(tdm, 3); 1931 PetscCall(DMPlexTransformSetDM(tr, dm)); 1932 1933 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &rdm)); 1934 PetscCall(DMSetType(rdm, DMPLEX)); 1935 PetscCall(DMPlexTransformSetDimensions(tr, dm, rdm)); 1936 /* Calculate number of new points of each depth */ 1937 PetscCall(DMPlexIsInterpolatedCollective(dm, &interp)); 1938 PetscCheck(interp == DMPLEX_INTERPOLATED_FULL,PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement"); 1939 /* Step 1: Set chart */ 1940 PetscCall(DMPlexSetChart(rdm, 0, tr->ctStartNew[tr->ctOrderNew[DM_NUM_POLYTOPES]])); 1941 /* Step 2: Set cone/support sizes (automatically stratifies) */ 1942 PetscCall(DMPlexTransformSetConeSizes(tr, rdm)); 1943 /* Step 3: Setup refined DM */ 1944 PetscCall(DMSetUp(rdm)); 1945 /* Step 4: Set cones and supports (automatically symmetrizes) */ 1946 PetscCall(DMPlexTransformSetCones(tr, rdm)); 1947 /* Step 5: Create pointSF */ 1948 PetscCall(DMPlexTransformCreateSF(tr, rdm)); 1949 /* Step 6: Create labels */ 1950 PetscCall(DMPlexTransformCreateLabels(tr, rdm)); 1951 /* Step 7: Set coordinates */ 1952 PetscCall(DMPlexTransformSetCoordinates(tr, rdm)); 1953 PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, rdm)); 1954 *tdm = rdm; 1955 PetscFunctionReturn(0); 1956 } 1957 1958 PetscErrorCode DMPlexTransformAdaptLabel(DM dm, PETSC_UNUSED Vec metric, DMLabel adaptLabel, PETSC_UNUSED DMLabel rgLabel, DM *rdm) 1959 { 1960 DMPlexTransform tr; 1961 DM cdm, rcdm; 1962 const char *prefix; 1963 1964 PetscFunctionBegin; 1965 PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject) dm), &tr)); 1966 PetscCall(PetscObjectGetOptionsPrefix((PetscObject) dm, &prefix)); 1967 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) tr, prefix)); 1968 PetscCall(DMPlexTransformSetDM(tr, dm)); 1969 PetscCall(DMPlexTransformSetFromOptions(tr)); 1970 PetscCall(DMPlexTransformSetActive(tr, adaptLabel)); 1971 PetscCall(DMPlexTransformSetUp(tr)); 1972 PetscCall(PetscObjectViewFromOptions((PetscObject) tr, NULL, "-dm_plex_transform_view")); 1973 PetscCall(DMPlexTransformApply(tr, dm, rdm)); 1974 PetscCall(DMCopyDisc(dm, *rdm)); 1975 PetscCall(DMGetCoordinateDM(dm, &cdm)); 1976 PetscCall(DMGetCoordinateDM(*rdm, &rcdm)); 1977 PetscCall(DMCopyDisc(cdm, rcdm)); 1978 PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm)); 1979 PetscCall(DMCopyDisc(dm, *rdm)); 1980 PetscCall(DMPlexTransformDestroy(&tr)); 1981 ((DM_Plex *) (*rdm)->data)->useHashLocation = ((DM_Plex *) dm->data)->useHashLocation; 1982 PetscFunctionReturn(0); 1983 } 1984