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