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 PetscErrorCode DMPlexTransformMapCoordinates(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) 1398 { 1399 PetscErrorCode ierr; 1400 1401 PetscFunctionBeginHot; 1402 ierr = (*tr->ops->mapcoordinates)(tr, pct, ct, r, Nv, dE, in, out);CHKERRQ(ierr); 1403 PetscFunctionReturn(0); 1404 } 1405 1406 static PetscErrorCode RefineLabel_Internal(DMPlexTransform tr, DMLabel label, DMLabel labelNew) 1407 { 1408 DM dm; 1409 IS valueIS; 1410 const PetscInt *values; 1411 PetscInt defVal, Nv, val; 1412 PetscErrorCode ierr; 1413 1414 PetscFunctionBegin; 1415 ierr = DMPlexTransformGetDM(tr, &dm);CHKERRQ(ierr); 1416 ierr = DMLabelGetDefaultValue(label, &defVal);CHKERRQ(ierr); 1417 ierr = DMLabelSetDefaultValue(labelNew, defVal);CHKERRQ(ierr); 1418 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 1419 ierr = ISGetLocalSize(valueIS, &Nv);CHKERRQ(ierr); 1420 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 1421 for (val = 0; val < Nv; ++val) { 1422 IS pointIS; 1423 const PetscInt *points; 1424 PetscInt numPoints, p; 1425 1426 /* Ensure refined label is created with same number of strata as 1427 * original (even if no entries here). */ 1428 ierr = DMLabelAddStratum(labelNew, values[val]);CHKERRQ(ierr); 1429 ierr = DMLabelGetStratumIS(label, values[val], &pointIS);CHKERRQ(ierr); 1430 ierr = ISGetLocalSize(pointIS, &numPoints);CHKERRQ(ierr); 1431 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 1432 for (p = 0; p < numPoints; ++p) { 1433 const PetscInt point = points[p]; 1434 DMPolytopeType ct; 1435 DMPolytopeType *rct; 1436 PetscInt *rsize, *rcone, *rornt; 1437 PetscInt Nct, n, r, pNew=0; 1438 1439 ierr = DMPlexGetCellType(dm, point, &ct);CHKERRQ(ierr); 1440 ierr = DMPlexTransformCellTransform(tr, ct, point, NULL, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 1441 for (n = 0; n < Nct; ++n) { 1442 for (r = 0; r < rsize[n]; ++r) { 1443 ierr = DMPlexTransformGetTargetPoint(tr, ct, rct[n], point, r, &pNew);CHKERRQ(ierr); 1444 ierr = DMLabelSetValue(labelNew, pNew, values[val]);CHKERRQ(ierr); 1445 } 1446 } 1447 } 1448 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 1449 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 1450 } 1451 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 1452 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 1453 PetscFunctionReturn(0); 1454 } 1455 1456 static PetscErrorCode DMPlexTransformCreateLabels(DMPlexTransform tr, DM rdm) 1457 { 1458 DM dm; 1459 PetscInt numLabels, l; 1460 PetscErrorCode ierr; 1461 1462 PetscFunctionBegin; 1463 ierr = DMPlexTransformGetDM(tr, &dm);CHKERRQ(ierr); 1464 ierr = DMGetNumLabels(dm, &numLabels);CHKERRQ(ierr); 1465 for (l = 0; l < numLabels; ++l) { 1466 DMLabel label, labelNew; 1467 const char *lname; 1468 PetscBool isDepth, isCellType; 1469 1470 ierr = DMGetLabelName(dm, l, &lname);CHKERRQ(ierr); 1471 ierr = PetscStrcmp(lname, "depth", &isDepth);CHKERRQ(ierr); 1472 if (isDepth) continue; 1473 ierr = PetscStrcmp(lname, "celltype", &isCellType);CHKERRQ(ierr); 1474 if (isCellType) continue; 1475 ierr = DMCreateLabel(rdm, lname);CHKERRQ(ierr); 1476 ierr = DMGetLabel(dm, lname, &label);CHKERRQ(ierr); 1477 ierr = DMGetLabel(rdm, lname, &labelNew);CHKERRQ(ierr); 1478 ierr = RefineLabel_Internal(tr, label, labelNew);CHKERRQ(ierr); 1479 } 1480 PetscFunctionReturn(0); 1481 } 1482 1483 /* This refines the labels which define regions for fields and DSes since they are not in the list of labels for the DM */ 1484 PetscErrorCode DMPlexTransformCreateDiscLabels(DMPlexTransform tr, DM rdm) 1485 { 1486 DM dm; 1487 PetscInt Nf, f, Nds, s; 1488 PetscErrorCode ierr; 1489 1490 PetscFunctionBegin; 1491 ierr = DMPlexTransformGetDM(tr, &dm);CHKERRQ(ierr); 1492 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 1493 for (f = 0; f < Nf; ++f) { 1494 DMLabel label, labelNew; 1495 PetscObject obj; 1496 const char *lname; 1497 1498 ierr = DMGetField(rdm, f, &label, &obj);CHKERRQ(ierr); 1499 if (!label) continue; 1500 ierr = PetscObjectGetName((PetscObject) label, &lname);CHKERRQ(ierr); 1501 ierr = DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);CHKERRQ(ierr); 1502 ierr = RefineLabel_Internal(tr, label, labelNew);CHKERRQ(ierr); 1503 ierr = DMSetField_Internal(rdm, f, labelNew, obj);CHKERRQ(ierr); 1504 ierr = DMLabelDestroy(&labelNew);CHKERRQ(ierr); 1505 } 1506 ierr = DMGetNumDS(dm, &Nds);CHKERRQ(ierr); 1507 for (s = 0; s < Nds; ++s) { 1508 DMLabel label, labelNew; 1509 const char *lname; 1510 1511 ierr = DMGetRegionNumDS(rdm, s, &label, NULL, NULL);CHKERRQ(ierr); 1512 if (!label) continue; 1513 ierr = PetscObjectGetName((PetscObject) label, &lname);CHKERRQ(ierr); 1514 ierr = DMLabelCreate(PETSC_COMM_SELF, lname, &labelNew);CHKERRQ(ierr); 1515 ierr = RefineLabel_Internal(tr, label, labelNew);CHKERRQ(ierr); 1516 ierr = DMSetRegionNumDS(rdm, s, labelNew, NULL, NULL);CHKERRQ(ierr); 1517 ierr = DMLabelDestroy(&labelNew);CHKERRQ(ierr); 1518 } 1519 PetscFunctionReturn(0); 1520 } 1521 1522 static PetscErrorCode DMPlexTransformCreateSF(DMPlexTransform tr, DM rdm) 1523 { 1524 DM dm; 1525 PetscSF sf, sfNew; 1526 PetscInt numRoots, numLeaves, numLeavesNew = 0, l, m; 1527 const PetscInt *localPoints; 1528 const PetscSFNode *remotePoints; 1529 PetscInt *localPointsNew; 1530 PetscSFNode *remotePointsNew; 1531 PetscInt pStartNew, pEndNew, pNew; 1532 /* Brute force algorithm */ 1533 PetscSF rsf; 1534 PetscSection s; 1535 const PetscInt *rootdegree; 1536 PetscInt *rootPointsNew, *remoteOffsets; 1537 PetscInt numPointsNew, pStart, pEnd, p; 1538 PetscErrorCode ierr; 1539 1540 PetscFunctionBegin; 1541 ierr = DMPlexTransformGetDM(tr, &dm);CHKERRQ(ierr); 1542 ierr = DMPlexGetChart(rdm, &pStartNew, &pEndNew);CHKERRQ(ierr); 1543 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 1544 ierr = DMGetPointSF(rdm, &sfNew);CHKERRQ(ierr); 1545 /* Calculate size of new SF */ 1546 ierr = PetscSFGetGraph(sf, &numRoots, &numLeaves, &localPoints, &remotePoints);CHKERRQ(ierr); 1547 if (numRoots < 0) PetscFunctionReturn(0); 1548 for (l = 0; l < numLeaves; ++l) { 1549 const PetscInt p = localPoints[l]; 1550 DMPolytopeType ct; 1551 DMPolytopeType *rct; 1552 PetscInt *rsize, *rcone, *rornt; 1553 PetscInt Nct, n; 1554 1555 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 1556 ierr = DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 1557 for (n = 0; n < Nct; ++n) { 1558 numLeavesNew += rsize[n]; 1559 } 1560 } 1561 /* Send new root point numbers 1562 It is possible to optimize for regular transforms by sending only the cell type offsets, but it seems a needless complication 1563 */ 1564 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1565 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &s);CHKERRQ(ierr); 1566 ierr = PetscSectionSetChart(s, pStart, pEnd);CHKERRQ(ierr); 1567 for (p = pStart; p < pEnd; ++p) { 1568 DMPolytopeType ct; 1569 DMPolytopeType *rct; 1570 PetscInt *rsize, *rcone, *rornt; 1571 PetscInt Nct, n; 1572 1573 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 1574 ierr = DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 1575 for (n = 0; n < Nct; ++n) { 1576 ierr = PetscSectionAddDof(s, p, rsize[n]);CHKERRQ(ierr); 1577 } 1578 } 1579 ierr = PetscSectionSetUp(s);CHKERRQ(ierr); 1580 ierr = PetscSectionGetStorageSize(s, &numPointsNew);CHKERRQ(ierr); 1581 ierr = PetscSFCreateRemoteOffsets(sf, s, s, &remoteOffsets);CHKERRQ(ierr); 1582 ierr = PetscSFCreateSectionSF(sf, s, remoteOffsets, s, &rsf);CHKERRQ(ierr); 1583 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 1584 ierr = PetscSFComputeDegreeBegin(sf, &rootdegree);CHKERRQ(ierr); 1585 ierr = PetscSFComputeDegreeEnd(sf, &rootdegree);CHKERRQ(ierr); 1586 ierr = PetscMalloc1(numPointsNew, &rootPointsNew);CHKERRQ(ierr); 1587 for (p = 0; p < numPointsNew; ++p) rootPointsNew[p] = -1; 1588 for (p = pStart; p < pEnd; ++p) { 1589 DMPolytopeType ct; 1590 DMPolytopeType *rct; 1591 PetscInt *rsize, *rcone, *rornt; 1592 PetscInt Nct, n, r, off; 1593 1594 if (!rootdegree[p-pStart]) continue; 1595 ierr = PetscSectionGetOffset(s, p, &off);CHKERRQ(ierr); 1596 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 1597 ierr = DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 1598 for (n = 0, m = 0; n < Nct; ++n) { 1599 for (r = 0; r < rsize[n]; ++r, ++m) { 1600 ierr = DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);CHKERRQ(ierr); 1601 rootPointsNew[off+m] = pNew; 1602 } 1603 } 1604 } 1605 ierr = PetscSFBcastBegin(rsf, MPIU_INT, rootPointsNew, rootPointsNew,MPI_REPLACE);CHKERRQ(ierr); 1606 ierr = PetscSFBcastEnd(rsf, MPIU_INT, rootPointsNew, rootPointsNew,MPI_REPLACE);CHKERRQ(ierr); 1607 ierr = PetscSFDestroy(&rsf);CHKERRQ(ierr); 1608 ierr = PetscMalloc1(numLeavesNew, &localPointsNew);CHKERRQ(ierr); 1609 ierr = PetscMalloc1(numLeavesNew, &remotePointsNew);CHKERRQ(ierr); 1610 for (l = 0, m = 0; l < numLeaves; ++l) { 1611 const PetscInt p = localPoints[l]; 1612 DMPolytopeType ct; 1613 DMPolytopeType *rct; 1614 PetscInt *rsize, *rcone, *rornt; 1615 PetscInt Nct, n, r, q, off; 1616 1617 ierr = PetscSectionGetOffset(s, p, &off);CHKERRQ(ierr); 1618 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 1619 ierr = DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 1620 for (n = 0, q = 0; n < Nct; ++n) { 1621 for (r = 0; r < rsize[n]; ++r, ++m, ++q) { 1622 ierr = DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);CHKERRQ(ierr); 1623 localPointsNew[m] = pNew; 1624 remotePointsNew[m].index = rootPointsNew[off+q]; 1625 remotePointsNew[m].rank = remotePoints[l].rank; 1626 } 1627 } 1628 } 1629 ierr = PetscSectionDestroy(&s);CHKERRQ(ierr); 1630 ierr = PetscFree(rootPointsNew);CHKERRQ(ierr); 1631 /* SF needs sorted leaves to correctly calculate Gather */ 1632 { 1633 PetscSFNode *rp, *rtmp; 1634 PetscInt *lp, *idx, *ltmp, i; 1635 1636 ierr = PetscMalloc1(numLeavesNew, &idx);CHKERRQ(ierr); 1637 ierr = PetscMalloc1(numLeavesNew, &lp);CHKERRQ(ierr); 1638 ierr = PetscMalloc1(numLeavesNew, &rp);CHKERRQ(ierr); 1639 for (i = 0; i < numLeavesNew; ++i) { 1640 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); 1641 idx[i] = i; 1642 } 1643 ierr = PetscSortIntWithPermutation(numLeavesNew, localPointsNew, idx);CHKERRQ(ierr); 1644 for (i = 0; i < numLeavesNew; ++i) { 1645 lp[i] = localPointsNew[idx[i]]; 1646 rp[i] = remotePointsNew[idx[i]]; 1647 } 1648 ltmp = localPointsNew; 1649 localPointsNew = lp; 1650 rtmp = remotePointsNew; 1651 remotePointsNew = rp; 1652 ierr = PetscFree(idx);CHKERRQ(ierr); 1653 ierr = PetscFree(ltmp);CHKERRQ(ierr); 1654 ierr = PetscFree(rtmp);CHKERRQ(ierr); 1655 } 1656 ierr = PetscSFSetGraph(sfNew, pEndNew-pStartNew, numLeavesNew, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER);CHKERRQ(ierr); 1657 PetscFunctionReturn(0); 1658 } 1659 1660 /* 1661 DMPlexCellRefinerMapLocalizedCoordinates - Given a cell of type ct with localized coordinates x, we generate localized coordinates xr for subcell r of type rct. 1662 1663 Not collective 1664 1665 Input Parameters: 1666 + tr - The DMPlexTransform 1667 . ct - The type of the parent cell 1668 . rct - The type of the produced cell 1669 . r - The index of the produced cell 1670 - x - The localized coordinates for the parent cell 1671 1672 Output Parameter: 1673 . xr - The localized coordinates for the produced cell 1674 1675 Level: developer 1676 1677 .seealso: DMPlexCellRefinerSetCoordinates() 1678 */ 1679 static PetscErrorCode DMPlexTransformMapLocalizedCoordinates(DMPlexTransform tr, DMPolytopeType ct, DMPolytopeType rct, PetscInt r, const PetscScalar x[], PetscScalar xr[]) 1680 { 1681 PetscFE fe = NULL; 1682 PetscInt cdim, v, *subcellV; 1683 PetscErrorCode ierr; 1684 1685 PetscFunctionBegin; 1686 ierr = DMPlexTransformGetCoordinateFE(tr, ct, &fe);CHKERRQ(ierr); 1687 ierr = DMPlexTransformGetSubcellVertices(tr, ct, rct, r, &subcellV);CHKERRQ(ierr); 1688 ierr = PetscFEGetNumComponents(fe, &cdim);CHKERRQ(ierr); 1689 for (v = 0; v < DMPolytopeTypeGetNumVertices(rct); ++v) { 1690 ierr = PetscFEInterpolate_Static(fe, x, tr->refGeom[ct], subcellV[v], &xr[v*cdim]);CHKERRQ(ierr); 1691 } 1692 PetscFunctionReturn(0); 1693 } 1694 1695 static PetscErrorCode DMPlexTransformSetCoordinates(DMPlexTransform tr, DM rdm) 1696 { 1697 DM dm, cdm; 1698 PetscSection coordSection, coordSectionNew; 1699 Vec coordsLocal, coordsLocalNew; 1700 const PetscScalar *coords; 1701 PetscScalar *coordsNew; 1702 const DMBoundaryType *bd; 1703 const PetscReal *maxCell, *L; 1704 PetscBool isperiodic, localizeVertices = PETSC_FALSE, localizeCells = PETSC_FALSE; 1705 PetscInt dE, d, cStart, cEnd, c, vStartNew, vEndNew, v, pStart, pEnd, p, ocStart, ocEnd; 1706 PetscErrorCode ierr; 1707 1708 PetscFunctionBegin; 1709 ierr = DMPlexTransformGetDM(tr, &dm);CHKERRQ(ierr); 1710 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 1711 ierr = DMGetPeriodicity(dm, &isperiodic, &maxCell, &L, &bd);CHKERRQ(ierr); 1712 /* Determine if we need to localize coordinates when generating them */ 1713 if (isperiodic) { 1714 localizeVertices = PETSC_TRUE; 1715 if (!maxCell) { 1716 PetscBool localized; 1717 ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr); 1718 if (!localized) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Cannot refine a periodic mesh if coordinates have not been localized"); 1719 localizeCells = PETSC_TRUE; 1720 } 1721 } 1722 1723 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 1724 ierr = PetscSectionGetFieldComponents(coordSection, 0, &dE);CHKERRQ(ierr); 1725 if (maxCell) { 1726 PetscReal maxCellNew[3]; 1727 1728 for (d = 0; d < dE; ++d) maxCellNew[d] = maxCell[d]/2.0; 1729 ierr = DMSetPeriodicity(rdm, isperiodic, maxCellNew, L, bd);CHKERRQ(ierr); 1730 } else { 1731 ierr = DMSetPeriodicity(rdm, isperiodic, maxCell, L, bd);CHKERRQ(ierr); 1732 } 1733 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &coordSectionNew);CHKERRQ(ierr); 1734 ierr = PetscSectionSetNumFields(coordSectionNew, 1);CHKERRQ(ierr); 1735 ierr = PetscSectionSetFieldComponents(coordSectionNew, 0, dE);CHKERRQ(ierr); 1736 ierr = DMPlexGetDepthStratum(rdm, 0, &vStartNew, &vEndNew);CHKERRQ(ierr); 1737 if (localizeCells) {ierr = PetscSectionSetChart(coordSectionNew, 0, vEndNew);CHKERRQ(ierr);} 1738 else {ierr = PetscSectionSetChart(coordSectionNew, vStartNew, vEndNew);CHKERRQ(ierr);} 1739 1740 /* Localization should be inherited */ 1741 /* Stefano calculates parent cells for each new cell for localization */ 1742 /* Localized cells need coordinates of closure */ 1743 for (v = vStartNew; v < vEndNew; ++v) { 1744 ierr = PetscSectionSetDof(coordSectionNew, v, dE);CHKERRQ(ierr); 1745 ierr = PetscSectionSetFieldDof(coordSectionNew, v, 0, dE);CHKERRQ(ierr); 1746 } 1747 if (localizeCells) { 1748 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1749 for (c = cStart; c < cEnd; ++c) { 1750 PetscInt dof; 1751 1752 ierr = PetscSectionGetDof(coordSection, c, &dof);CHKERRQ(ierr); 1753 if (dof) { 1754 DMPolytopeType ct; 1755 DMPolytopeType *rct; 1756 PetscInt *rsize, *rcone, *rornt; 1757 PetscInt dim, cNew, Nct, n, r; 1758 1759 ierr = DMPlexGetCellType(dm, c, &ct);CHKERRQ(ierr); 1760 dim = DMPolytopeTypeGetDim(ct); 1761 ierr = DMPlexTransformCellTransform(tr, ct, c, NULL, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 1762 /* This allows for different cell types */ 1763 for (n = 0; n < Nct; ++n) { 1764 if (dim != DMPolytopeTypeGetDim(rct[n])) continue; 1765 for (r = 0; r < rsize[n]; ++r) { 1766 PetscInt *closure = NULL; 1767 PetscInt clSize, cl, Nv = 0; 1768 1769 ierr = DMPlexTransformGetTargetPoint(tr, ct, rct[n], c, r, &cNew);CHKERRQ(ierr); 1770 ierr = DMPlexGetTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 1771 for (cl = 0; cl < clSize*2; cl += 2) {if ((closure[cl] >= vStartNew) && (closure[cl] < vEndNew)) ++Nv;} 1772 ierr = DMPlexRestoreTransitiveClosure(rdm, cNew, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 1773 ierr = PetscSectionSetDof(coordSectionNew, cNew, Nv * dE);CHKERRQ(ierr); 1774 ierr = PetscSectionSetFieldDof(coordSectionNew, cNew, 0, Nv * dE);CHKERRQ(ierr); 1775 } 1776 } 1777 } 1778 } 1779 } 1780 ierr = PetscSectionSetUp(coordSectionNew);CHKERRQ(ierr); 1781 ierr = DMViewFromOptions(dm, NULL, "-coarse_dm_view");CHKERRQ(ierr); 1782 ierr = DMSetCoordinateSection(rdm, PETSC_DETERMINE, coordSectionNew);CHKERRQ(ierr); 1783 { 1784 VecType vtype; 1785 PetscInt coordSizeNew, bs; 1786 const char *name; 1787 1788 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr); 1789 ierr = VecCreate(PETSC_COMM_SELF, &coordsLocalNew);CHKERRQ(ierr); 1790 ierr = PetscSectionGetStorageSize(coordSectionNew, &coordSizeNew);CHKERRQ(ierr); 1791 ierr = VecSetSizes(coordsLocalNew, coordSizeNew, PETSC_DETERMINE);CHKERRQ(ierr); 1792 ierr = PetscObjectGetName((PetscObject) coordsLocal, &name);CHKERRQ(ierr); 1793 ierr = PetscObjectSetName((PetscObject) coordsLocalNew, name);CHKERRQ(ierr); 1794 ierr = VecGetBlockSize(coordsLocal, &bs);CHKERRQ(ierr); 1795 ierr = VecSetBlockSize(coordsLocalNew, bs);CHKERRQ(ierr); 1796 ierr = VecGetType(coordsLocal, &vtype);CHKERRQ(ierr); 1797 ierr = VecSetType(coordsLocalNew, vtype);CHKERRQ(ierr); 1798 } 1799 ierr = VecGetArrayRead(coordsLocal, &coords);CHKERRQ(ierr); 1800 ierr = VecGetArray(coordsLocalNew, &coordsNew);CHKERRQ(ierr); 1801 ierr = PetscSectionGetChart(coordSection, &ocStart, &ocEnd);CHKERRQ(ierr); 1802 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1803 /* First set coordinates for vertices*/ 1804 for (p = pStart; p < pEnd; ++p) { 1805 DMPolytopeType ct; 1806 DMPolytopeType *rct; 1807 PetscInt *rsize, *rcone, *rornt; 1808 PetscInt Nct, n, r; 1809 PetscBool hasVertex = PETSC_FALSE, isLocalized = PETSC_FALSE; 1810 1811 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 1812 ierr = DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 1813 for (n = 0; n < Nct; ++n) { 1814 if (rct[n] == DM_POLYTOPE_POINT) {hasVertex = PETSC_TRUE; break;} 1815 } 1816 if (localizeVertices && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) { 1817 PetscInt dof; 1818 ierr = PetscSectionGetDof(coordSection, p, &dof);CHKERRQ(ierr); 1819 if (dof) isLocalized = PETSC_TRUE; 1820 } 1821 if (hasVertex) { 1822 const PetscScalar *icoords = NULL; 1823 PetscScalar *pcoords = NULL; 1824 PetscInt Nc, Nv, v, d; 1825 1826 ierr = DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);CHKERRQ(ierr); 1827 1828 icoords = pcoords; 1829 Nv = Nc/dE; 1830 if (ct != DM_POLYTOPE_POINT) { 1831 if (localizeVertices) { 1832 PetscScalar anchor[3]; 1833 1834 for (d = 0; d < dE; ++d) anchor[d] = pcoords[d]; 1835 if (!isLocalized) { 1836 for (v = 0; v < Nv; ++v) {ierr = DMLocalizeCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], &pcoords[v*dE]);CHKERRQ(ierr);} 1837 } else { 1838 Nv = Nc/(2*dE); 1839 for (v = Nv; v < Nv*2; ++v) {ierr = DMLocalizeCoordinate_Internal(dm, dE, anchor, &pcoords[v*dE], &pcoords[v*dE]);CHKERRQ(ierr);} 1840 } 1841 } 1842 } 1843 for (n = 0; n < Nct; ++n) { 1844 if (rct[n] != DM_POLYTOPE_POINT) continue; 1845 for (r = 0; r < rsize[n]; ++r) { 1846 PetscScalar vcoords[3]; 1847 PetscInt vNew, off; 1848 1849 ierr = DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &vNew);CHKERRQ(ierr); 1850 ierr = PetscSectionGetOffset(coordSectionNew, vNew, &off);CHKERRQ(ierr); 1851 ierr = DMPlexTransformMapCoordinates(tr, ct, rct[n], r, Nv, dE, icoords, vcoords);CHKERRQ(ierr); 1852 ierr = DMPlexSnapToGeomModel(dm, p, vcoords, &coordsNew[off]);CHKERRQ(ierr); 1853 } 1854 } 1855 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &Nc, &pcoords);CHKERRQ(ierr); 1856 } 1857 } 1858 /* Then set coordinates for cells by localizing */ 1859 for (p = pStart; p < pEnd; ++p) { 1860 DMPolytopeType ct; 1861 DMPolytopeType *rct; 1862 PetscInt *rsize, *rcone, *rornt; 1863 PetscInt Nct, n, r; 1864 PetscBool isLocalized = PETSC_FALSE; 1865 1866 ierr = DMPlexGetCellType(dm, p, &ct);CHKERRQ(ierr); 1867 ierr = DMPlexTransformCellTransform(tr, ct, p, NULL, &Nct, &rct, &rsize, &rcone, &rornt);CHKERRQ(ierr); 1868 if (localizeCells && ct != DM_POLYTOPE_POINT && (p >= ocStart) && (p < ocEnd)) { 1869 PetscInt dof; 1870 ierr = PetscSectionGetDof(coordSection, p, &dof);CHKERRQ(ierr); 1871 if (dof) isLocalized = PETSC_TRUE; 1872 } 1873 if (isLocalized) { 1874 const PetscScalar *pcoords; 1875 1876 ierr = DMPlexPointLocalRead(cdm, p, coords, &pcoords);CHKERRQ(ierr); 1877 for (n = 0; n < Nct; ++n) { 1878 const PetscInt Nr = rsize[n]; 1879 1880 if (DMPolytopeTypeGetDim(ct) != DMPolytopeTypeGetDim(rct[n])) continue; 1881 for (r = 0; r < Nr; ++r) { 1882 PetscInt pNew, offNew; 1883 1884 /* It looks like Stefano and Lisandro are allowing localized coordinates without defining the periodic boundary, which means that 1885 DMLocalizeCoordinate_Internal() will not work. Localized coordinates will have to have obtained by the affine map of the larger 1886 cell to the ones it produces. */ 1887 ierr = DMPlexTransformGetTargetPoint(tr, ct, rct[n], p, r, &pNew);CHKERRQ(ierr); 1888 ierr = PetscSectionGetOffset(coordSectionNew, pNew, &offNew);CHKERRQ(ierr); 1889 ierr = DMPlexTransformMapLocalizedCoordinates(tr, ct, rct[n], r, pcoords, &coordsNew[offNew]);CHKERRQ(ierr); 1890 } 1891 } 1892 } 1893 } 1894 ierr = VecRestoreArrayRead(coordsLocal, &coords);CHKERRQ(ierr); 1895 ierr = VecRestoreArray(coordsLocalNew, &coordsNew);CHKERRQ(ierr); 1896 ierr = DMSetCoordinatesLocal(rdm, coordsLocalNew);CHKERRQ(ierr); 1897 /* TODO Stefano has a final reduction if some hybrid coordinates cannot be found. (needcoords) Should not be needed. */ 1898 ierr = VecDestroy(&coordsLocalNew);CHKERRQ(ierr); 1899 ierr = PetscSectionDestroy(&coordSectionNew);CHKERRQ(ierr); 1900 if (!localizeCells) {ierr = DMLocalizeCoordinates(rdm);CHKERRQ(ierr);} 1901 PetscFunctionReturn(0); 1902 } 1903 1904 PetscErrorCode DMPlexTransformApply(DMPlexTransform tr, DM dm, DM *tdm) 1905 { 1906 DM rdm; 1907 DMPlexInterpolatedFlag interp; 1908 PetscInt dim, embedDim; 1909 PetscErrorCode ierr; 1910 1911 PetscFunctionBegin; 1912 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1913 PetscValidHeaderSpecific(dm, DM_CLASSID, 2); 1914 PetscValidPointer(tdm, 3); 1915 ierr = DMPlexTransformSetDM(tr, dm);CHKERRQ(ierr); 1916 1917 ierr = DMCreate(PetscObjectComm((PetscObject)dm), &rdm);CHKERRQ(ierr); 1918 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 1919 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1920 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 1921 ierr = DMGetCoordinateDim(dm, &embedDim);CHKERRQ(ierr); 1922 ierr = DMSetCoordinateDim(rdm, embedDim);CHKERRQ(ierr); 1923 /* Calculate number of new points of each depth */ 1924 ierr = DMPlexIsInterpolated(dm, &interp);CHKERRQ(ierr); 1925 if (interp != DMPLEX_INTERPOLATED_FULL) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Mesh must be fully interpolated for regular refinement"); 1926 /* Step 1: Set chart */ 1927 ierr = DMPlexSetChart(rdm, 0, tr->ctStartNew[tr->ctOrder[DM_NUM_POLYTOPES]]);CHKERRQ(ierr); 1928 /* Step 2: Set cone/support sizes (automatically stratifies) */ 1929 ierr = DMPlexTransformSetConeSizes(tr, rdm);CHKERRQ(ierr); 1930 /* Step 3: Setup refined DM */ 1931 ierr = DMSetUp(rdm);CHKERRQ(ierr); 1932 /* Step 4: Set cones and supports (automatically symmetrizes) */ 1933 ierr = DMPlexTransformSetCones(tr, rdm);CHKERRQ(ierr); 1934 /* Step 5: Create pointSF */ 1935 ierr = DMPlexTransformCreateSF(tr, rdm);CHKERRQ(ierr); 1936 /* Step 6: Create labels */ 1937 ierr = DMPlexTransformCreateLabels(tr, rdm);CHKERRQ(ierr); 1938 /* Step 7: Set coordinates */ 1939 ierr = DMPlexTransformSetCoordinates(tr, rdm);CHKERRQ(ierr); 1940 *tdm = rdm; 1941 PetscFunctionReturn(0); 1942 } 1943 1944 PetscErrorCode DMPlexTransformAdaptLabel(DM dm, DMLabel adaptLabel, DM *rdm) 1945 { 1946 DMPlexTransform tr; 1947 DM cdm, rcdm; 1948 const char *prefix; 1949 PetscErrorCode ierr; 1950 1951 PetscFunctionBegin; 1952 ierr = DMPlexTransformCreate(PetscObjectComm((PetscObject) dm), &tr);CHKERRQ(ierr); 1953 ierr = PetscObjectGetOptionsPrefix((PetscObject) dm, &prefix);CHKERRQ(ierr); 1954 ierr = PetscObjectSetOptionsPrefix((PetscObject) tr, prefix);CHKERRQ(ierr); 1955 ierr = DMPlexTransformSetDM(tr, dm);CHKERRQ(ierr); 1956 ierr = DMPlexTransformSetFromOptions(tr);CHKERRQ(ierr); 1957 ierr = DMPlexTransformSetActive(tr, adaptLabel);CHKERRQ(ierr); 1958 ierr = DMPlexTransformSetUp(tr);CHKERRQ(ierr); 1959 ierr = PetscObjectViewFromOptions((PetscObject) tr, NULL, "-dm_plex_transform_view");CHKERRQ(ierr); 1960 ierr = DMPlexTransformApply(tr, dm, rdm);CHKERRQ(ierr); 1961 ierr = DMCopyDisc(dm, *rdm);CHKERRQ(ierr); 1962 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 1963 ierr = DMGetCoordinateDM(*rdm, &rcdm);CHKERRQ(ierr); 1964 ierr = DMCopyDisc(cdm, rcdm);CHKERRQ(ierr); 1965 ierr = DMPlexTransformCreateDiscLabels(tr, *rdm);CHKERRQ(ierr); 1966 ierr = DMCopyDisc(dm, *rdm);CHKERRQ(ierr); 1967 ierr = DMPlexTransformDestroy(&tr);CHKERRQ(ierr); 1968 PetscFunctionReturn(0); 1969 } 1970