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