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