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