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