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