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