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