1 #include <petsc/private/dmplextransformimpl.h> /*I "petscdmplextransform.h" I*/ 2 3 static PetscErrorCode DMPlexTransformView_Extrude(DMPlexTransform tr, PetscViewer viewer) 4 { 5 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 6 PetscBool isascii; 7 8 PetscFunctionBegin; 9 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 10 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 11 PetscCall(PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii)); 12 if (isascii) { 13 const char *name; 14 15 PetscCall(PetscObjectGetName((PetscObject) tr, &name)); 16 PetscCall(PetscViewerASCIIPrintf(viewer, "Extrusion transformation %s\n", name ? name : "")); 17 PetscCall(PetscViewerASCIIPrintf(viewer, " number of layers: %" PetscInt_FMT "\n", ex->layers)); 18 PetscCall(PetscViewerASCIIPrintf(viewer, " create tensor cells: %s\n", ex->useTensor ? "YES" : "NO")); 19 } else { 20 SETERRQ(PetscObjectComm((PetscObject) tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject) viewer)->type_name); 21 } 22 PetscFunctionReturn(0); 23 } 24 25 static PetscErrorCode DMPlexTransformSetFromOptions_Extrude(PetscOptionItems *PetscOptionsObject, DMPlexTransform tr) 26 { 27 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 28 PetscReal th, normal[3], *thicknesses; 29 PetscInt nl, Nc; 30 PetscBool tensor, sym, flg; 31 char funcname[PETSC_MAX_PATH_LEN]; 32 33 PetscFunctionBegin; 34 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 2); 35 PetscOptionsHeadBegin(PetscOptionsObject, "DMPlexTransform Extrusion Options"); 36 PetscCall(PetscOptionsBoundedInt("-dm_plex_transform_extrude_layers", "Number of layers to extrude", "", ex->layers, &nl, &flg, 1)); 37 if (flg) PetscCall(DMPlexTransformExtrudeSetLayers(tr, nl)); 38 PetscCall(PetscOptionsReal("-dm_plex_transform_extrude_thickness", "Total thickness of extruded layers", "", ex->thickness, &th, &flg)); 39 if (flg) PetscCall(DMPlexTransformExtrudeSetThickness(tr, th)); 40 PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_use_tensor", "Create tensor cells", "", ex->useTensor, &tensor, &flg)); 41 if (flg) PetscCall(DMPlexTransformExtrudeSetTensor(tr, tensor)); 42 PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_symmetric", "Extrude layers symmetrically about the surface", "", ex->symmetric, &sym, &flg)); 43 if (flg) PetscCall(DMPlexTransformExtrudeSetSymmetric(tr, sym)); 44 Nc = 3; 45 PetscCall(PetscOptionsRealArray("-dm_plex_transform_extrude_normal", "Input normal vector for extrusion", "DMPlexTransformExtrudeSetNormal", normal, &Nc, &flg)); 46 if (flg) { 47 // Extrusion dimension might not yet be determined 48 PetscCheck(!ex->cdimEx || Nc == ex->cdimEx,PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_SIZ, "Input normal has size %" PetscInt_FMT " != %" PetscInt_FMT " extruded coordinate dimension", Nc, ex->cdimEx); 49 PetscCall(DMPlexTransformExtrudeSetNormal(tr, normal)); 50 } 51 PetscCall(PetscOptionsString("-dm_plex_transform_extrude_normal_function", "Function to determine normal vector", "DMPlexTransformExtrudeSetNormalFunction", funcname, funcname, sizeof(funcname), &flg)); 52 if (flg) { 53 PetscSimplePointFunc normalFunc; 54 55 PetscCall(PetscDLSym(NULL, funcname, (void **) &normalFunc)); 56 PetscCall(DMPlexTransformExtrudeSetNormalFunction(tr, normalFunc)); 57 } 58 nl = ex->layers; 59 PetscCall(PetscMalloc1(nl, &thicknesses)); 60 PetscCall(PetscOptionsRealArray("-dm_plex_transform_extrude_thicknesses", "Thickness of each individual extruded layer", "", thicknesses, &nl, &flg)); 61 if (flg) { 62 PetscCheck(nl,PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_OUTOFRANGE, "Must give at least one thickness for -dm_plex_transform_extrude_thicknesses"); 63 PetscCall(DMPlexTransformExtrudeSetThicknesses(tr, nl, thicknesses)); 64 } 65 PetscCall(PetscFree(thicknesses)); 66 PetscOptionsHeadEnd(); 67 PetscFunctionReturn(0); 68 } 69 70 /* Determine the implicit dimension pre-extrusion (either the implicit dimension of the DM or of a point in the active set for the transform). 71 If that dimension is the same as the current coordinate dimension (ex->dim), the extruded mesh will have a coordinate dimension one greater; 72 Otherwise the coordinate dimension will be kept. */ 73 static PetscErrorCode DMPlexTransformExtrudeComputeExtrusionDim(DMPlexTransform tr) 74 { 75 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 76 DM dm; 77 DMLabel active; 78 PetscInt dim, dimOrig; 79 80 PetscFunctionBegin; 81 PetscCall(DMPlexTransformGetDM(tr, &dm)); 82 PetscCall(DMGetDimension(dm, &dim)); 83 dimOrig = dim; 84 PetscCall(DMPlexTransformGetActive(tr, &active)); 85 if (active) { 86 PetscInt pStart, pEnd, p; 87 88 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 89 for (p = pStart; p < pEnd; ++p) { 90 DMPolytopeType ct; 91 PetscInt val; 92 93 PetscCall(DMLabelGetValue(active, p, &val)); 94 if (val < 0) continue; 95 PetscCall(DMPlexGetCellType(dm, p, &ct)); 96 dim = DMPolytopeTypeGetDim(ct); 97 break; 98 } 99 } 100 ex->dimEx = PetscMax(dimOrig, dim+1); 101 ex->cdimEx = ex->cdim == dim ? ex->cdim+1 : ex->cdim; 102 PetscFunctionReturn(0); 103 } 104 105 static PetscErrorCode DMPlexTransformSetDimensions_Extrude(DMPlexTransform tr, DM dm, DM tdm) 106 { 107 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 108 PetscInt dim; 109 110 PetscFunctionBegin; 111 PetscCall(DMGetDimension(dm, &dim)); 112 PetscCall(DMSetDimension(tdm, ex->dimEx)); 113 PetscCall(DMSetCoordinateDim(tdm, ex->cdimEx)); 114 PetscFunctionReturn(0); 115 } 116 117 /* 118 The refine types for extrusion are: 119 120 ct: For any normally extruded point 121 ct + 100: For any point which should just return itself 122 */ 123 static PetscErrorCode DMPlexTransformSetUp_Extrude(DMPlexTransform tr) 124 { 125 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 126 DM dm; 127 DMLabel active; 128 DMPolytopeType ct; 129 PetscInt Nl = ex->layers, l, i, ict, Nc, No, coff, ooff; 130 131 PetscFunctionBegin; 132 PetscCall(DMPlexTransformExtrudeComputeExtrusionDim(tr)); 133 PetscCall(DMPlexTransformGetDM(tr, &dm)); 134 PetscCall(DMPlexTransformGetActive(tr, &active)); 135 if (active) { 136 DMLabel celltype; 137 PetscInt pStart, pEnd, p; 138 139 PetscCall(DMPlexGetCellTypeLabel(dm, &celltype)); 140 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType)); 141 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 142 for (p = pStart; p < pEnd; ++p) { 143 PetscInt ct, val; 144 145 PetscCall(DMLabelGetValue(celltype, p, &ct)); 146 PetscCall(DMLabelGetValue(active, p, &val)); 147 if (val < 0) {PetscCall(DMLabelSetValue(tr->trType, p, ct + 100));} 148 else {PetscCall(DMLabelSetValue(tr->trType, p, ct));} 149 } 150 } 151 PetscCall(PetscMalloc5(DM_NUM_POLYTOPES, &ex->Nt, DM_NUM_POLYTOPES, &ex->target, DM_NUM_POLYTOPES, &ex->size, DM_NUM_POLYTOPES, &ex->cone, DM_NUM_POLYTOPES, &ex->ornt)); 152 for (ict = 0; ict < DM_NUM_POLYTOPES; ++ict) { 153 ex->Nt[ict] = -1; 154 ex->target[ict] = NULL; 155 ex->size[ict] = NULL; 156 ex->cone[ict] = NULL; 157 ex->ornt[ict] = NULL; 158 } 159 /* DM_POLYTOPE_POINT produces Nl+1 points and Nl segments/tensor segments */ 160 ct = DM_POLYTOPE_POINT; 161 ex->Nt[ct] = 2; 162 Nc = 6*Nl; 163 No = 2*Nl; 164 PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct])); 165 ex->target[ct][0] = DM_POLYTOPE_POINT; 166 ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_POINT_PRISM_TENSOR : DM_POLYTOPE_SEGMENT; 167 ex->size[ct][0] = Nl+1; 168 ex->size[ct][1] = Nl; 169 /* cones for segments/tensor segments */ 170 for (i = 0; i < Nl; ++i) { 171 ex->cone[ct][6*i+0] = DM_POLYTOPE_POINT; 172 ex->cone[ct][6*i+1] = 0; 173 ex->cone[ct][6*i+2] = i; 174 ex->cone[ct][6*i+3] = DM_POLYTOPE_POINT; 175 ex->cone[ct][6*i+4] = 0; 176 ex->cone[ct][6*i+5] = i+1; 177 } 178 for (i = 0; i < No; ++i) ex->ornt[ct][i] = 0; 179 /* DM_POLYTOPE_SEGMENT produces Nl+1 segments and Nl quads/tensor quads */ 180 ct = DM_POLYTOPE_SEGMENT; 181 ex->Nt[ct] = 2; 182 Nc = 8*(Nl+1) + 14*Nl; 183 No = 2*(Nl+1) + 4*Nl; 184 PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct])); 185 ex->target[ct][0] = DM_POLYTOPE_SEGMENT; 186 ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_SEG_PRISM_TENSOR : DM_POLYTOPE_QUADRILATERAL; 187 ex->size[ct][0] = Nl+1; 188 ex->size[ct][1] = Nl; 189 /* cones for segments */ 190 for (i = 0; i < Nl+1; ++i) { 191 ex->cone[ct][8*i+0] = DM_POLYTOPE_POINT; 192 ex->cone[ct][8*i+1] = 1; 193 ex->cone[ct][8*i+2] = 0; 194 ex->cone[ct][8*i+3] = i; 195 ex->cone[ct][8*i+4] = DM_POLYTOPE_POINT; 196 ex->cone[ct][8*i+5] = 1; 197 ex->cone[ct][8*i+6] = 1; 198 ex->cone[ct][8*i+7] = i; 199 } 200 for (i = 0; i < 2*(Nl+1); ++i) ex->ornt[ct][i] = 0; 201 /* cones for quads/tensor quads */ 202 coff = 8*(Nl+1); 203 ooff = 2*(Nl+1); 204 for (i = 0; i < Nl; ++i) { 205 if (ex->useTensor) { 206 ex->cone[ct][coff+14*i+0] = DM_POLYTOPE_SEGMENT; 207 ex->cone[ct][coff+14*i+1] = 0; 208 ex->cone[ct][coff+14*i+2] = i; 209 ex->cone[ct][coff+14*i+3] = DM_POLYTOPE_SEGMENT; 210 ex->cone[ct][coff+14*i+4] = 0; 211 ex->cone[ct][coff+14*i+5] = i+1; 212 ex->cone[ct][coff+14*i+6] = DM_POLYTOPE_POINT_PRISM_TENSOR; 213 ex->cone[ct][coff+14*i+7] = 1; 214 ex->cone[ct][coff+14*i+8] = 0; 215 ex->cone[ct][coff+14*i+9] = i; 216 ex->cone[ct][coff+14*i+10] = DM_POLYTOPE_POINT_PRISM_TENSOR; 217 ex->cone[ct][coff+14*i+11] = 1; 218 ex->cone[ct][coff+14*i+12] = 1; 219 ex->cone[ct][coff+14*i+13] = i; 220 ex->ornt[ct][ooff+4*i+0] = 0; 221 ex->ornt[ct][ooff+4*i+1] = 0; 222 ex->ornt[ct][ooff+4*i+2] = 0; 223 ex->ornt[ct][ooff+4*i+3] = 0; 224 } else { 225 ex->cone[ct][coff+14*i+0] = DM_POLYTOPE_SEGMENT; 226 ex->cone[ct][coff+14*i+1] = 0; 227 ex->cone[ct][coff+14*i+2] = i; 228 ex->cone[ct][coff+14*i+3] = DM_POLYTOPE_SEGMENT; 229 ex->cone[ct][coff+14*i+4] = 1; 230 ex->cone[ct][coff+14*i+5] = 1; 231 ex->cone[ct][coff+14*i+6] = i; 232 ex->cone[ct][coff+14*i+7] = DM_POLYTOPE_SEGMENT; 233 ex->cone[ct][coff+14*i+8] = 0; 234 ex->cone[ct][coff+14*i+9] = i+1; 235 ex->cone[ct][coff+14*i+10] = DM_POLYTOPE_SEGMENT; 236 ex->cone[ct][coff+14*i+11] = 1; 237 ex->cone[ct][coff+14*i+12] = 0; 238 ex->cone[ct][coff+14*i+13] = i; 239 ex->ornt[ct][ooff+4*i+0] = 0; 240 ex->ornt[ct][ooff+4*i+1] = 0; 241 ex->ornt[ct][ooff+4*i+2] = -1; 242 ex->ornt[ct][ooff+4*i+3] = -1; 243 } 244 } 245 /* DM_POLYTOPE_TRIANGLE produces Nl+1 triangles and Nl triangular prisms/tensor triangular prisms */ 246 ct = DM_POLYTOPE_TRIANGLE; 247 ex->Nt[ct] = 2; 248 Nc = 12*(Nl+1) + 18*Nl; 249 No = 3*(Nl+1) + 5*Nl; 250 PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct])); 251 ex->target[ct][0] = DM_POLYTOPE_TRIANGLE; 252 ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_TRI_PRISM_TENSOR : DM_POLYTOPE_TRI_PRISM; 253 ex->size[ct][0] = Nl+1; 254 ex->size[ct][1] = Nl; 255 /* cones for triangles */ 256 for (i = 0; i < Nl+1; ++i) { 257 ex->cone[ct][12*i+0] = DM_POLYTOPE_SEGMENT; 258 ex->cone[ct][12*i+1] = 1; 259 ex->cone[ct][12*i+2] = 0; 260 ex->cone[ct][12*i+3] = i; 261 ex->cone[ct][12*i+4] = DM_POLYTOPE_SEGMENT; 262 ex->cone[ct][12*i+5] = 1; 263 ex->cone[ct][12*i+6] = 1; 264 ex->cone[ct][12*i+7] = i; 265 ex->cone[ct][12*i+8] = DM_POLYTOPE_SEGMENT; 266 ex->cone[ct][12*i+9] = 1; 267 ex->cone[ct][12*i+10] = 2; 268 ex->cone[ct][12*i+11] = i; 269 } 270 for (i = 0; i < 3*(Nl+1); ++i) ex->ornt[ct][i] = 0; 271 /* cones for triangular prisms/tensor triangular prisms */ 272 coff = 12*(Nl+1); 273 ooff = 3*(Nl+1); 274 for (i = 0; i < Nl; ++i) { 275 if (ex->useTensor) { 276 ex->cone[ct][coff+18*i+0] = DM_POLYTOPE_TRIANGLE; 277 ex->cone[ct][coff+18*i+1] = 0; 278 ex->cone[ct][coff+18*i+2] = i; 279 ex->cone[ct][coff+18*i+3] = DM_POLYTOPE_TRIANGLE; 280 ex->cone[ct][coff+18*i+4] = 0; 281 ex->cone[ct][coff+18*i+5] = i+1; 282 ex->cone[ct][coff+18*i+6] = DM_POLYTOPE_SEG_PRISM_TENSOR; 283 ex->cone[ct][coff+18*i+7] = 1; 284 ex->cone[ct][coff+18*i+8] = 0; 285 ex->cone[ct][coff+18*i+9] = i; 286 ex->cone[ct][coff+18*i+10] = DM_POLYTOPE_SEG_PRISM_TENSOR; 287 ex->cone[ct][coff+18*i+11] = 1; 288 ex->cone[ct][coff+18*i+12] = 1; 289 ex->cone[ct][coff+18*i+13] = i; 290 ex->cone[ct][coff+18*i+14] = DM_POLYTOPE_SEG_PRISM_TENSOR; 291 ex->cone[ct][coff+18*i+15] = 1; 292 ex->cone[ct][coff+18*i+16] = 2; 293 ex->cone[ct][coff+18*i+17] = i; 294 ex->ornt[ct][ooff+5*i+0] = 0; 295 ex->ornt[ct][ooff+5*i+1] = 0; 296 ex->ornt[ct][ooff+5*i+2] = 0; 297 ex->ornt[ct][ooff+5*i+3] = 0; 298 ex->ornt[ct][ooff+5*i+4] = 0; 299 } else { 300 ex->cone[ct][coff+18*i+0] = DM_POLYTOPE_TRIANGLE; 301 ex->cone[ct][coff+18*i+1] = 0; 302 ex->cone[ct][coff+18*i+2] = i; 303 ex->cone[ct][coff+18*i+3] = DM_POLYTOPE_TRIANGLE; 304 ex->cone[ct][coff+18*i+4] = 0; 305 ex->cone[ct][coff+18*i+5] = i+1; 306 ex->cone[ct][coff+18*i+6] = DM_POLYTOPE_QUADRILATERAL; 307 ex->cone[ct][coff+18*i+7] = 1; 308 ex->cone[ct][coff+18*i+8] = 0; 309 ex->cone[ct][coff+18*i+9] = i; 310 ex->cone[ct][coff+18*i+10] = DM_POLYTOPE_QUADRILATERAL; 311 ex->cone[ct][coff+18*i+11] = 1; 312 ex->cone[ct][coff+18*i+12] = 1; 313 ex->cone[ct][coff+18*i+13] = i; 314 ex->cone[ct][coff+18*i+14] = DM_POLYTOPE_QUADRILATERAL; 315 ex->cone[ct][coff+18*i+15] = 1; 316 ex->cone[ct][coff+18*i+16] = 2; 317 ex->cone[ct][coff+18*i+17] = i; 318 ex->ornt[ct][ooff+5*i+0] = -2; 319 ex->ornt[ct][ooff+5*i+1] = 0; 320 ex->ornt[ct][ooff+5*i+2] = 0; 321 ex->ornt[ct][ooff+5*i+3] = 0; 322 ex->ornt[ct][ooff+5*i+4] = 0; 323 } 324 } 325 /* DM_POLYTOPE_QUADRILATERAL produces Nl+1 quads and Nl hexes/tensor hexes */ 326 ct = DM_POLYTOPE_QUADRILATERAL; 327 ex->Nt[ct] = 2; 328 Nc = 16*(Nl+1) + 22*Nl; 329 No = 4*(Nl+1) + 6*Nl; 330 PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct])); 331 ex->target[ct][0] = DM_POLYTOPE_QUADRILATERAL; 332 ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_QUAD_PRISM_TENSOR : DM_POLYTOPE_HEXAHEDRON; 333 ex->size[ct][0] = Nl+1; 334 ex->size[ct][1] = Nl; 335 /* cones for quads */ 336 for (i = 0; i < Nl+1; ++i) { 337 ex->cone[ct][16*i+0] = DM_POLYTOPE_SEGMENT; 338 ex->cone[ct][16*i+1] = 1; 339 ex->cone[ct][16*i+2] = 0; 340 ex->cone[ct][16*i+3] = i; 341 ex->cone[ct][16*i+4] = DM_POLYTOPE_SEGMENT; 342 ex->cone[ct][16*i+5] = 1; 343 ex->cone[ct][16*i+6] = 1; 344 ex->cone[ct][16*i+7] = i; 345 ex->cone[ct][16*i+8] = DM_POLYTOPE_SEGMENT; 346 ex->cone[ct][16*i+9] = 1; 347 ex->cone[ct][16*i+10] = 2; 348 ex->cone[ct][16*i+11] = i; 349 ex->cone[ct][16*i+12] = DM_POLYTOPE_SEGMENT; 350 ex->cone[ct][16*i+13] = 1; 351 ex->cone[ct][16*i+14] = 3; 352 ex->cone[ct][16*i+15] = i; 353 } 354 for (i = 0; i < 4*(Nl+1); ++i) ex->ornt[ct][i] = 0; 355 /* cones for hexes/tensor hexes */ 356 coff = 16*(Nl+1); 357 ooff = 4*(Nl+1); 358 for (i = 0; i < Nl; ++i) { 359 if (ex->useTensor) { 360 ex->cone[ct][coff+22*i+0] = DM_POLYTOPE_QUADRILATERAL; 361 ex->cone[ct][coff+22*i+1] = 0; 362 ex->cone[ct][coff+22*i+2] = i; 363 ex->cone[ct][coff+22*i+3] = DM_POLYTOPE_QUADRILATERAL; 364 ex->cone[ct][coff+22*i+4] = 0; 365 ex->cone[ct][coff+22*i+5] = i+1; 366 ex->cone[ct][coff+22*i+6] = DM_POLYTOPE_SEG_PRISM_TENSOR; 367 ex->cone[ct][coff+22*i+7] = 1; 368 ex->cone[ct][coff+22*i+8] = 0; 369 ex->cone[ct][coff+22*i+9] = i; 370 ex->cone[ct][coff+22*i+10] = DM_POLYTOPE_SEG_PRISM_TENSOR; 371 ex->cone[ct][coff+22*i+11] = 1; 372 ex->cone[ct][coff+22*i+12] = 1; 373 ex->cone[ct][coff+22*i+13] = i; 374 ex->cone[ct][coff+22*i+14] = DM_POLYTOPE_SEG_PRISM_TENSOR; 375 ex->cone[ct][coff+22*i+15] = 1; 376 ex->cone[ct][coff+22*i+16] = 2; 377 ex->cone[ct][coff+22*i+17] = i; 378 ex->cone[ct][coff+22*i+18] = DM_POLYTOPE_SEG_PRISM_TENSOR; 379 ex->cone[ct][coff+22*i+19] = 1; 380 ex->cone[ct][coff+22*i+20] = 3; 381 ex->cone[ct][coff+22*i+21] = i; 382 ex->ornt[ct][ooff+6*i+0] = 0; 383 ex->ornt[ct][ooff+6*i+1] = 0; 384 ex->ornt[ct][ooff+6*i+2] = 0; 385 ex->ornt[ct][ooff+6*i+3] = 0; 386 ex->ornt[ct][ooff+6*i+4] = 0; 387 ex->ornt[ct][ooff+6*i+5] = 0; 388 } else { 389 ex->cone[ct][coff+22*i+0] = DM_POLYTOPE_QUADRILATERAL; 390 ex->cone[ct][coff+22*i+1] = 0; 391 ex->cone[ct][coff+22*i+2] = i; 392 ex->cone[ct][coff+22*i+3] = DM_POLYTOPE_QUADRILATERAL; 393 ex->cone[ct][coff+22*i+4] = 0; 394 ex->cone[ct][coff+22*i+5] = i+1; 395 ex->cone[ct][coff+22*i+6] = DM_POLYTOPE_QUADRILATERAL; 396 ex->cone[ct][coff+22*i+7] = 1; 397 ex->cone[ct][coff+22*i+8] = 0; 398 ex->cone[ct][coff+22*i+9] = i; 399 ex->cone[ct][coff+22*i+10] = DM_POLYTOPE_QUADRILATERAL; 400 ex->cone[ct][coff+22*i+11] = 1; 401 ex->cone[ct][coff+22*i+12] = 2; 402 ex->cone[ct][coff+22*i+13] = i; 403 ex->cone[ct][coff+22*i+14] = DM_POLYTOPE_QUADRILATERAL; 404 ex->cone[ct][coff+22*i+15] = 1; 405 ex->cone[ct][coff+22*i+16] = 1; 406 ex->cone[ct][coff+22*i+17] = i; 407 ex->cone[ct][coff+22*i+18] = DM_POLYTOPE_QUADRILATERAL; 408 ex->cone[ct][coff+22*i+19] = 1; 409 ex->cone[ct][coff+22*i+20] = 3; 410 ex->cone[ct][coff+22*i+21] = i; 411 ex->ornt[ct][ooff+6*i+0] = -2; 412 ex->ornt[ct][ooff+6*i+1] = 0; 413 ex->ornt[ct][ooff+6*i+2] = 0; 414 ex->ornt[ct][ooff+6*i+3] = 0; 415 ex->ornt[ct][ooff+6*i+4] = 0; 416 ex->ornt[ct][ooff+6*i+5] = 1; 417 } 418 } 419 /* Layers positions */ 420 if (!ex->Nth) { 421 if (ex->symmetric) for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness*l)/ex->layers - ex->thickness/2; 422 else for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness*l)/ex->layers; 423 } else { 424 ex->thickness = 0.; 425 ex->layerPos[0] = 0.; 426 for (l = 0; l < ex->layers; ++l) { 427 const PetscReal t = ex->thicknesses[PetscMin(l, ex->Nth-1)]; 428 ex->layerPos[l+1] = ex->layerPos[l] + t; 429 ex->thickness += t; 430 } 431 if (ex->symmetric) for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] -= ex->thickness/2.; 432 } 433 PetscFunctionReturn(0); 434 } 435 436 static PetscErrorCode DMPlexTransformDestroy_Extrude(DMPlexTransform tr) 437 { 438 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 439 PetscInt ct; 440 441 PetscFunctionBegin; 442 for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) { 443 PetscCall(PetscFree4(ex->target[ct], ex->size[ct], ex->cone[ct], ex->ornt[ct])); 444 } 445 PetscCall(PetscFree5(ex->Nt, ex->target, ex->size, ex->cone, ex->ornt)); 446 PetscCall(PetscFree(ex->layerPos)); 447 PetscCall(PetscFree(ex)); 448 PetscFunctionReturn(0); 449 } 450 451 static PetscErrorCode DMPlexTransformGetSubcellOrientation_Extrude(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 452 { 453 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 454 DMLabel trType = tr->trType; 455 PetscInt rt; 456 457 PetscFunctionBeginHot; 458 *rnew = r; 459 *onew = DMPolytopeTypeComposeOrientation(tct, o, so); 460 if (!so) PetscFunctionReturn(0); 461 if (trType) { 462 PetscCall(DMLabelGetValue(tr->trType, sp, &rt)); 463 if (rt >= 100) PetscFunctionReturn(0); 464 } 465 if (ex->useTensor) { 466 switch (sct) { 467 case DM_POLYTOPE_POINT: break; 468 case DM_POLYTOPE_SEGMENT: 469 switch (tct) { 470 case DM_POLYTOPE_SEGMENT: break; 471 case DM_POLYTOPE_SEG_PRISM_TENSOR: 472 *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0); 473 break; 474 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 475 } 476 break; 477 // We need to handle identity extrusions from volumes (TET, HEX, etc) when boundary faces are being extruded 478 case DM_POLYTOPE_TRIANGLE: break; 479 case DM_POLYTOPE_QUADRILATERAL: break; 480 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]); 481 } 482 } else { 483 switch (sct) { 484 case DM_POLYTOPE_POINT: break; 485 case DM_POLYTOPE_SEGMENT: 486 switch (tct) { 487 case DM_POLYTOPE_SEGMENT: break; 488 case DM_POLYTOPE_QUADRILATERAL: 489 *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -3 : 0); 490 break; 491 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 492 } 493 break; 494 // We need to handle identity extrusions from volumes (TET, HEX, etc) when boundary faces are being extruded 495 case DM_POLYTOPE_TRIANGLE: break; 496 case DM_POLYTOPE_QUADRILATERAL: break; 497 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]); 498 } 499 } 500 PetscFunctionReturn(0); 501 } 502 503 static PetscErrorCode DMPlexTransformCellTransform_Extrude(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 504 { 505 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 506 DMLabel trType = tr->trType; 507 PetscBool ignore = PETSC_FALSE, identity = PETSC_FALSE; 508 PetscInt val = 0; 509 510 PetscFunctionBegin; 511 if (trType) { 512 PetscCall(DMLabelGetValue(trType, p, &val)); 513 identity = val >= 100 ? PETSC_TRUE : PETSC_FALSE; 514 } else { 515 ignore = ex->Nt[source] < 0 ? PETSC_TRUE : PETSC_FALSE; 516 } 517 if (rt) *rt = val; 518 if (ignore) { 519 /* Ignore cells that cannot be extruded */ 520 *Nt = 0; 521 *target = NULL; 522 *size = NULL; 523 *cone = NULL; 524 *ornt = NULL; 525 } else if (identity) { 526 PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt)); 527 } else { 528 *Nt = ex->Nt[source]; 529 *target = ex->target[source]; 530 *size = ex->size[source]; 531 *cone = ex->cone[source]; 532 *ornt = ex->ornt[source]; 533 } 534 PetscFunctionReturn(0); 535 } 536 537 /* Computes new vertex along normal */ 538 static PetscErrorCode DMPlexTransformMapCoordinates_Extrude(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) 539 { 540 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 541 DM dm; 542 DMLabel active; 543 PetscReal ones2[2] = {0., 1.}, ones3[3] = { 0., 0., 1.}; 544 PetscReal normal[3] = {0., 0., 0.}, norm; 545 PetscBool computeNormal; 546 PetscInt dim, dEx = ex->cdimEx, cStart, cEnd, d; 547 548 PetscFunctionBeginHot; 549 PetscCheck(pct == DM_POLYTOPE_POINT,PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for parent point type %s",DMPolytopeTypes[pct]); 550 PetscCheck(ct == DM_POLYTOPE_POINT,PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for refined point type %s",DMPolytopeTypes[ct]); 551 PetscCheck(Nv == 1,PETSC_COMM_SELF,PETSC_ERR_SUP,"Vertices should be produced from a single vertex, not %" PetscInt_FMT,Nv); 552 PetscCheck(dE == ex->cdim,PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Coordinate dim %" PetscInt_FMT " != %" PetscInt_FMT " original dimension", dE, ex->cdim); 553 554 PetscCall(DMPlexTransformGetDM(tr, &dm)); 555 PetscCall(DMGetDimension(dm, &dim)); 556 PetscCall(DMPlexTransformGetActive(tr, &active)); 557 computeNormal = dim != ex->cdim && !ex->useNormal ? PETSC_TRUE : PETSC_FALSE; 558 if (computeNormal) { 559 PetscInt *closure = NULL; 560 PetscInt closureSize, cl; 561 562 PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 563 PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure)); 564 for (cl = 0; cl < closureSize*2; cl += 2) { 565 if ((closure[cl] >= cStart) && (closure[cl] < cEnd)) { 566 PetscReal cnormal[3] = {0, 0, 0}; 567 568 PetscCall(DMPlexComputeCellGeometryFVM(dm, closure[cl], NULL, NULL, cnormal)); 569 for (d = 0; d < dEx; ++d) normal[d] += cnormal[d]; 570 } 571 } 572 PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure)); 573 } else if (ex->useNormal) { 574 for (d = 0; d < dEx; ++d) normal[d] = ex->normal[d]; 575 } else if (active) { // find an active point in the closure of p and use its coordinate normal as the extrusion direction 576 PetscInt *closure = NULL; 577 PetscInt closureSize, cl, pStart, pEnd; 578 579 PetscCall(DMPlexGetDepthStratum(dm, ex->cdimEx-1, &pStart, &pEnd)); 580 PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure)); 581 for (cl = 0; cl < closureSize*2; cl += 2) { 582 if ((closure[cl] >= pStart) && (closure[cl] < pEnd)) { 583 PetscReal cnormal[3] = {0, 0, 0}; 584 const PetscInt *supp; 585 PetscInt suppSize; 586 587 PetscCall(DMPlexComputeCellGeometryFVM(dm, closure[cl], NULL, NULL, cnormal)); 588 PetscCall(DMPlexGetSupportSize(dm, closure[cl], &suppSize)); 589 PetscCall(DMPlexGetSupport(dm, closure[cl], &supp)); 590 // Only use external faces, so I can get the orientation from any cell 591 if (suppSize == 1) { 592 const PetscInt *cone, *ornt; 593 PetscInt coneSize, c; 594 595 PetscCall(DMPlexGetConeSize(dm, supp[0], &coneSize)); 596 PetscCall(DMPlexGetCone(dm, supp[0], &cone)); 597 PetscCall(DMPlexGetConeOrientation(dm, supp[0], &ornt)); 598 for (c = 0; c < coneSize; ++c) if (cone[c] == closure[cl]) break; 599 PetscCheck(c < coneSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Asymmetry in cone/support"); 600 if (ornt[c] < 0) for (d = 0; d < dEx; ++d) cnormal[d] *= -1.; 601 for (d = 0; d < dEx; ++d) normal[d] += cnormal[d]; 602 } 603 } 604 } 605 PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure)); 606 } else if (ex->cdimEx == 2) { 607 for (d = 0; d < dEx; ++d) normal[d] = ones2[d]; 608 } else if (ex->cdimEx == 3) { 609 for (d = 0; d < dEx; ++d) normal[d] = ones3[d]; 610 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to determine normal for extrusion"); 611 if (ex->normalFunc) { 612 PetscScalar n[3]; 613 PetscReal x[3], dot; 614 615 for (d = 0; d < ex->cdim; ++d) x[d] = PetscRealPart(in[d]); 616 PetscCall((*ex->normalFunc)(ex->cdim, 0., x, r, n, NULL)); 617 for (dot=0,d=0; d < dEx; d++) dot += PetscRealPart(n[d]) * normal[d]; 618 for (d = 0; d < dEx; ++d) normal[d] = PetscSign(dot) * PetscRealPart(n[d]); 619 } 620 621 for (d = 0, norm = 0.0; d < dEx; ++d) norm += PetscSqr(normal[d]); 622 for (d = 0; d < dEx; ++d) normal[d] *= norm == 0.0 ? 1.0 : 1./PetscSqrtReal(norm); 623 for (d = 0; d < dEx; ++d) out[d] = normal[d]*ex->layerPos[r]; 624 for (d = 0; d < ex->cdim; ++d) out[d] += in[d]; 625 PetscFunctionReturn(0); 626 } 627 628 static PetscErrorCode DMPlexTransformInitialize_Extrude(DMPlexTransform tr) 629 { 630 PetscFunctionBegin; 631 tr->ops->view = DMPlexTransformView_Extrude; 632 tr->ops->setfromoptions = DMPlexTransformSetFromOptions_Extrude; 633 tr->ops->setup = DMPlexTransformSetUp_Extrude; 634 tr->ops->destroy = DMPlexTransformDestroy_Extrude; 635 tr->ops->setdimensions = DMPlexTransformSetDimensions_Extrude; 636 tr->ops->celltransform = DMPlexTransformCellTransform_Extrude; 637 tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Extrude; 638 tr->ops->mapcoordinates = DMPlexTransformMapCoordinates_Extrude; 639 PetscFunctionReturn(0); 640 } 641 642 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform tr) 643 { 644 DMPlexTransform_Extrude *ex; 645 DM dm; 646 PetscInt dim; 647 648 PetscFunctionBegin; 649 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 650 PetscCall(PetscNewLog(tr, &ex)); 651 tr->data = ex; 652 ex->thickness = 1.; 653 ex->useTensor = PETSC_TRUE; 654 ex->symmetric = PETSC_FALSE; 655 ex->useNormal = PETSC_FALSE; 656 ex->layerPos = NULL; 657 PetscCall(DMPlexTransformGetDM(tr, &dm)); 658 PetscCall(DMGetDimension(dm, &dim)); 659 PetscCall(DMGetCoordinateDim(dm, &ex->cdim)); 660 PetscCall(DMPlexTransformExtrudeSetLayers(tr, 1)); 661 PetscCall(DMPlexTransformInitialize_Extrude(tr)); 662 PetscFunctionReturn(0); 663 } 664 665 /*@ 666 DMPlexTransformExtrudeGetLayers - Get the number of extruded layers. 667 668 Not collective 669 670 Input Parameter: 671 . tr - The DMPlexTransform 672 673 Output Parameter: 674 . layers - The number of layers 675 676 Level: intermediate 677 678 .seealso: `DMPlexTransformExtrudeSetLayers()` 679 @*/ 680 PetscErrorCode DMPlexTransformExtrudeGetLayers(DMPlexTransform tr, PetscInt *layers) 681 { 682 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 683 684 PetscFunctionBegin; 685 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 686 PetscValidIntPointer(layers, 2); 687 *layers = ex->layers; 688 PetscFunctionReturn(0); 689 } 690 691 /*@ 692 DMPlexTransformExtrudeSetLayers - Set the number of extruded layers. 693 694 Not collective 695 696 Input Parameters: 697 + tr - The DMPlexTransform 698 - layers - The number of layers 699 700 Level: intermediate 701 702 .seealso: `DMPlexTransformExtrudeGetLayers()` 703 @*/ 704 PetscErrorCode DMPlexTransformExtrudeSetLayers(DMPlexTransform tr, PetscInt layers) 705 { 706 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 707 708 PetscFunctionBegin; 709 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 710 ex->layers = layers; 711 PetscCall(PetscFree(ex->layerPos)); 712 PetscCall(PetscCalloc1(ex->layers+1, &ex->layerPos)); 713 PetscFunctionReturn(0); 714 } 715 716 /*@ 717 DMPlexTransformExtrudeGetThickness - Get the total thickness of the layers 718 719 Not collective 720 721 Input Parameter: 722 . tr - The DMPlexTransform 723 724 Output Parameter: 725 . thickness - The total thickness of the layers 726 727 Level: intermediate 728 729 .seealso: `DMPlexTransformExtrudeSetThickness()` 730 @*/ 731 PetscErrorCode DMPlexTransformExtrudeGetThickness(DMPlexTransform tr, PetscReal *thickness) 732 { 733 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 734 735 PetscFunctionBegin; 736 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 737 PetscValidRealPointer(thickness, 2); 738 *thickness = ex->thickness; 739 PetscFunctionReturn(0); 740 } 741 742 /*@ 743 DMPlexTransformExtrudeSetThickness - Set the total thickness of the layers 744 745 Not collective 746 747 Input Parameters: 748 + tr - The DMPlexTransform 749 - thickness - The total thickness of the layers 750 751 Level: intermediate 752 753 .seealso: `DMPlexTransformExtrudeGetThickness()` 754 @*/ 755 PetscErrorCode DMPlexTransformExtrudeSetThickness(DMPlexTransform tr, PetscReal thickness) 756 { 757 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 758 759 PetscFunctionBegin; 760 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 761 PetscCheck(thickness > 0.,PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_OUTOFRANGE, "Height of layers %g must be positive", (double) thickness); 762 ex->thickness = thickness; 763 PetscFunctionReturn(0); 764 } 765 766 /*@ 767 DMPlexTransformExtrudeGetTensor - Get the flag to use tensor cells 768 769 Not collective 770 771 Input Parameter: 772 . tr - The DMPlexTransform 773 774 Output Parameter: 775 . useTensor - The flag to use tensor cells 776 777 Note: This flag determines the orientation behavior of the created points. For example, if tensor is PETSC_TRUE, then DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT, DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON. 778 779 Level: intermediate 780 781 .seealso: `DMPlexTransformExtrudeSetTensor()` 782 @*/ 783 PetscErrorCode DMPlexTransformExtrudeGetTensor(DMPlexTransform tr, PetscBool *useTensor) 784 { 785 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 786 787 PetscFunctionBegin; 788 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 789 PetscValidBoolPointer(useTensor, 2); 790 *useTensor = ex->useTensor; 791 PetscFunctionReturn(0); 792 } 793 794 /*@ 795 DMPlexTransformExtrudeSetTensor - Set the flag to use tensor cells 796 797 Not collective 798 799 Input Parameters: 800 + tr - The DMPlexTransform 801 - useTensor - The flag for tensor cells 802 803 Note: This flag determines the orientation behavior of the created points. For example, if tensor is PETSC_TRUE, then DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT, DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL, DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON. 804 805 Level: intermediate 806 807 .seealso: `DMPlexTransformExtrudeGetTensor()` 808 @*/ 809 PetscErrorCode DMPlexTransformExtrudeSetTensor(DMPlexTransform tr, PetscBool useTensor) 810 { 811 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 812 813 PetscFunctionBegin; 814 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 815 ex->useTensor = useTensor; 816 PetscFunctionReturn(0); 817 } 818 819 /*@ 820 DMPlexTransformExtrudeGetSymmetric - Get the flag to extrude symmetrically from the initial surface 821 822 Not collective 823 824 Input Parameter: 825 . tr - The DMPlexTransform 826 827 Output Parameter: 828 . symmetric - The flag to extrude symmetrically 829 830 Level: intermediate 831 832 .seealso: `DMPlexTransformExtrudeSetSymmetric()` 833 @*/ 834 PetscErrorCode DMPlexTransformExtrudeGetSymmetric(DMPlexTransform tr, PetscBool *symmetric) 835 { 836 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 837 838 PetscFunctionBegin; 839 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 840 PetscValidBoolPointer(symmetric, 2); 841 *symmetric = ex->symmetric; 842 PetscFunctionReturn(0); 843 } 844 845 /*@ 846 DMPlexTransformExtrudeSetSymmetric - Set the flag to extrude symmetrically from the initial surface 847 848 Not collective 849 850 Input Parameters: 851 + tr - The DMPlexTransform 852 - symmetric - The flag to extrude symmetrically 853 854 Level: intermediate 855 856 .seealso: `DMPlexTransformExtrudeGetSymmetric()` 857 @*/ 858 PetscErrorCode DMPlexTransformExtrudeSetSymmetric(DMPlexTransform tr, PetscBool symmetric) 859 { 860 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 861 862 PetscFunctionBegin; 863 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 864 ex->symmetric = symmetric; 865 PetscFunctionReturn(0); 866 } 867 868 /*@ 869 DMPlexTransformExtrudeGetNormal - Get the extrusion normal vector 870 871 Not collective 872 873 Input Parameter: 874 . tr - The DMPlexTransform 875 876 Output Parameter: 877 . normal - The extrusion direction 878 879 Note: The user passes in an array, which is filled by the library. 880 881 Level: intermediate 882 883 .seealso: `DMPlexTransformExtrudeSetNormal()` 884 @*/ 885 PetscErrorCode DMPlexTransformExtrudeGetNormal(DMPlexTransform tr, PetscReal normal[]) 886 { 887 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 888 PetscInt d; 889 890 PetscFunctionBegin; 891 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 892 if (ex->useNormal) {for (d = 0; d < ex->cdimEx; ++d) normal[d] = ex->normal[d];} 893 else {for (d = 0; d < ex->cdimEx; ++d) normal[d] = 0.;} 894 PetscFunctionReturn(0); 895 } 896 897 /*@ 898 DMPlexTransformExtrudeSetNormal - Set the extrusion normal 899 900 Not collective 901 902 Input Parameters: 903 + tr - The DMPlexTransform 904 - normal - The extrusion direction 905 906 Level: intermediate 907 908 .seealso: `DMPlexTransformExtrudeGetNormal()` 909 @*/ 910 PetscErrorCode DMPlexTransformExtrudeSetNormal(DMPlexTransform tr, const PetscReal normal[]) 911 { 912 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 913 PetscInt d; 914 915 PetscFunctionBegin; 916 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 917 ex->useNormal = PETSC_TRUE; 918 for (d = 0; d < ex->cdimEx; ++d) ex->normal[d] = normal[d]; 919 PetscFunctionReturn(0); 920 } 921 922 /*@C 923 DMPlexTransformExtrudeSetNormalFunction - Set a function to determine the extrusion normal 924 925 Not collective 926 927 Input Parameters: 928 + tr - The DMPlexTransform 929 - normalFunc - A function determining the extrusion direction 930 931 Note: The calling sequence for the function is normalFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt r, PetscScalar u[], void *ctx) 932 $ dim - The coordinate dimension of the original mesh (usually a surface) 933 $ time - The current time, or 0. 934 $ x - The location of the current normal, in the coordinate space of the original mesh 935 $ r - The extrusion replica number (layer number) of this point 936 $ u - The user provides the computed normal on output; the sign and magnitude is not significant 937 $ ctx - An optional user context 938 939 Level: intermediate 940 941 .seealso: `DMPlexTransformExtrudeGetNormal()` 942 @*/ 943 PetscErrorCode DMPlexTransformExtrudeSetNormalFunction(DMPlexTransform tr, PetscSimplePointFunc normalFunc) 944 { 945 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 946 947 PetscFunctionBegin; 948 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 949 ex->normalFunc = normalFunc; 950 PetscFunctionReturn(0); 951 } 952 953 /*@ 954 DMPlexTransformExtrudeSetThicknesses - Set the thickness of each layer 955 956 Not collective 957 958 Input Parameters: 959 + tr - The DMPlexTransform 960 . Nth - The number of thicknesses 961 - thickness - The array of thicknesses 962 963 Level: intermediate 964 965 .seealso: `DMPlexTransformExtrudeSetThickness()`, `DMPlexTransformExtrudeGetThickness()` 966 @*/ 967 PetscErrorCode DMPlexTransformExtrudeSetThicknesses(DMPlexTransform tr, PetscInt Nth, const PetscReal thicknesses[]) 968 { 969 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 970 PetscInt t; 971 972 PetscFunctionBegin; 973 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 974 PetscCheck(Nth > 0,PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_OUTOFRANGE, "Number of thicknesses %" PetscInt_FMT " must be positive", Nth); 975 ex->Nth = PetscMin(Nth, ex->layers); 976 PetscCall(PetscFree(ex->thicknesses)); 977 PetscCall(PetscMalloc1(ex->Nth, &ex->thicknesses)); 978 for (t = 0; t < ex->Nth; ++t) { 979 PetscCheck(thicknesses[t] > 0.,PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_OUTOFRANGE, "Thickness %g of layer %" PetscInt_FMT " must be positive", (double) thicknesses[t], t); 980 ex->thicknesses[t] = thicknesses[t]; 981 } 982 PetscFunctionReturn(0); 983 } 984