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