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