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