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