1 #include <petsc/private/dmplextransformimpl.h> /*I "petscdmplextransform.h" I*/ 2 3 const char *const PlexNormalAlgs[] = {"default", "input", "compute", "compute_bd"}; 4 5 static PetscErrorCode DMPlexTransformView_Extrude(DMPlexTransform tr, PetscViewer viewer) 6 { 7 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 8 PetscBool isascii; 9 10 PetscFunctionBegin; 11 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 12 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 13 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 14 if (isascii) { 15 DM dm; 16 DMLabel active; 17 PetscInt dim; 18 const char *name; 19 20 PetscCall(PetscObjectGetName((PetscObject)tr, &name)); 21 PetscCall(DMPlexTransformGetDM(tr, &dm)); 22 PetscCall(DMGetDimension(dm, &dim)); 23 PetscCall(DMPlexTransformGetActive(tr, &active)); 24 25 PetscCall(PetscViewerASCIIPrintf(viewer, "Extrusion transformation %s\n", name ? name : "")); 26 PetscCall(PetscViewerASCIIPrintf(viewer, " number of layers: %" PetscInt_FMT "\n", ex->layers)); 27 PetscCall(PetscViewerASCIIPrintf(viewer, " create tensor cells: %s\n", ex->useTensor ? "YES" : "NO")); 28 if (ex->periodic) PetscCall(PetscViewerASCIIPrintf(viewer, " periodic\n")); 29 PetscCall(PetscViewerASCIIPrintf(viewer, " normal algorithm: %s\n", PlexNormalAlgs[ex->normalAlg])); 30 if (ex->normalFunc) PetscCall(PetscViewerASCIIPrintf(viewer, " normal modified by user function\n")); 31 } else { 32 SETERRQ(PetscObjectComm((PetscObject)tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject)viewer)->type_name); 33 } 34 PetscFunctionReturn(PETSC_SUCCESS); 35 } 36 37 static PetscErrorCode DMPlexTransformSetFromOptions_Extrude(DMPlexTransform tr, PetscOptionItems PetscOptionsObject) 38 { 39 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 40 PetscReal th, normal[3], *thicknesses; 41 PetscInt nl, Nc; 42 PetscBool tensor, sym, per, flg; 43 char funcname[PETSC_MAX_PATH_LEN]; 44 45 PetscFunctionBegin; 46 PetscOptionsHeadBegin(PetscOptionsObject, "DMPlexTransform Extrusion Options"); 47 PetscCall(PetscOptionsBoundedInt("-dm_plex_transform_extrude_layers", "Number of layers to extrude", "", ex->layers, &nl, &flg, 1)); 48 if (flg) PetscCall(DMPlexTransformExtrudeSetLayers(tr, nl)); 49 PetscCall(PetscOptionsReal("-dm_plex_transform_extrude_thickness", "Total thickness of extruded layers", "", ex->thickness, &th, &flg)); 50 if (flg) PetscCall(DMPlexTransformExtrudeSetThickness(tr, th)); 51 PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_use_tensor", "Create tensor cells", "", ex->useTensor, &tensor, &flg)); 52 if (flg) PetscCall(DMPlexTransformExtrudeSetTensor(tr, tensor)); 53 PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_symmetric", "Extrude layers symmetrically about the surface", "", ex->symmetric, &sym, &flg)); 54 if (flg) PetscCall(DMPlexTransformExtrudeSetSymmetric(tr, sym)); 55 PetscCall(PetscOptionsBool("-dm_plex_transform_extrude_periodic", "Extrude layers periodically about the surface", "", ex->periodic, &per, &flg)); 56 if (flg) PetscCall(DMPlexTransformExtrudeSetPeriodic(tr, per)); 57 Nc = 3; 58 PetscCall(PetscOptionsRealArray("-dm_plex_transform_extrude_normal", "Input normal vector for extrusion", "DMPlexTransformExtrudeSetNormal", normal, &Nc, &flg)); 59 if (flg) { 60 // Extrusion dimension might not yet be determined 61 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); 62 PetscCall(DMPlexTransformExtrudeSetNormal(tr, normal)); 63 } 64 PetscCall(PetscOptionsString("-dm_plex_transform_extrude_normal_function", "Function to determine normal vector", "DMPlexTransformExtrudeSetNormalFunction", funcname, funcname, sizeof(funcname), &flg)); 65 if (flg) { 66 PetscSimplePointFn *normalFunc; 67 68 PetscCall(PetscDLSym(NULL, funcname, (void **)&normalFunc)); 69 PetscCall(DMPlexTransformExtrudeSetNormalFunction(tr, normalFunc)); 70 } 71 nl = ex->layers; 72 PetscCall(PetscMalloc1(nl, &thicknesses)); 73 PetscCall(PetscOptionsRealArray("-dm_plex_transform_extrude_thicknesses", "Thickness of each individual extruded layer", "", thicknesses, &nl, &flg)); 74 if (flg) { 75 PetscCheck(nl, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Must give at least one thickness for -dm_plex_transform_extrude_thicknesses"); 76 PetscCall(DMPlexTransformExtrudeSetThicknesses(tr, nl, thicknesses)); 77 } 78 PetscCall(PetscFree(thicknesses)); 79 PetscOptionsHeadEnd(); 80 PetscFunctionReturn(PETSC_SUCCESS); 81 } 82 83 /* Determine the implicit dimension pre-extrusion (either the implicit dimension of the DM or of a point in the active set for the transform). 84 If that dimension is the same as the current coordinate dimension (ex->dim), the extruded mesh will have a coordinate dimension one greater; 85 Otherwise the coordinate dimension will be kept. */ 86 static PetscErrorCode DMPlexTransformExtrudeComputeExtrusionDim(DMPlexTransform tr) 87 { 88 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 89 DM dm; 90 DMLabel active; 91 PetscInt dim, dimExtPoint, dimExtPointG; 92 93 PetscFunctionBegin; 94 PetscCall(DMPlexTransformGetDM(tr, &dm)); 95 PetscCall(DMGetDimension(dm, &dim)); 96 PetscCall(DMPlexTransformGetActive(tr, &active)); 97 if (active) { 98 IS valueIS, pointIS; 99 const PetscInt *values, *points; 100 DMPolytopeType ct; 101 PetscInt Nv, Np; 102 103 dimExtPoint = 0; 104 PetscCall(DMLabelGetValueIS(active, &valueIS)); 105 PetscCall(ISGetLocalSize(valueIS, &Nv)); 106 PetscCall(ISGetIndices(valueIS, &values)); 107 for (PetscInt v = 0; v < Nv; ++v) { 108 PetscCall(DMLabelGetStratumIS(active, values[v], &pointIS)); 109 PetscCall(ISGetLocalSize(pointIS, &Np)); 110 PetscCall(ISGetIndices(pointIS, &points)); 111 for (PetscInt p = 0; p < Np; ++p) { 112 PetscCall(DMPlexGetCellType(dm, points[p], &ct)); 113 dimExtPoint = PetscMax(dimExtPoint, DMPolytopeTypeGetDim(ct)); 114 } 115 PetscCall(ISRestoreIndices(pointIS, &points)); 116 PetscCall(ISDestroy(&pointIS)); 117 } 118 PetscCall(ISRestoreIndices(valueIS, &values)); 119 PetscCall(ISDestroy(&valueIS)); 120 } else dimExtPoint = dim; 121 PetscCallMPI(MPIU_Allreduce(&dimExtPoint, &dimExtPointG, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)tr))); 122 ex->dimEx = PetscMax(dim, dimExtPointG + 1); 123 ex->cdimEx = ex->cdim == dimExtPointG ? ex->cdim + 1 : ex->cdim; 124 PetscCheck(ex->dimEx <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Topological dimension for extruded mesh %" PetscInt_FMT " must not exceed 3", ex->dimEx); 125 PetscCheck(ex->cdimEx <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension for extruded mesh %" PetscInt_FMT " must not exceed 3", ex->cdimEx); 126 PetscFunctionReturn(PETSC_SUCCESS); 127 } 128 129 static PetscErrorCode DMPlexTransformSetDimensions_Extrude(DMPlexTransform tr, DM dm, DM tdm) 130 { 131 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 132 PetscInt dim; 133 134 PetscFunctionBegin; 135 PetscCall(DMGetDimension(dm, &dim)); 136 PetscCall(DMSetDimension(tdm, ex->dimEx)); 137 PetscCall(DMSetCoordinateDim(tdm, ex->cdimEx)); 138 PetscFunctionReturn(PETSC_SUCCESS); 139 } 140 141 /* DM_POLYTOPE_POINT produces 142 Nl+1 points, or Nl points periodically and 143 Nl segments, or tensor segments 144 */ 145 static PetscErrorCode DMPlexTransformExtrudeSetUp_Point(DMPlexTransform_Extrude *ex, PetscInt Nl) 146 { 147 const DMPolytopeType ct = DM_POLYTOPE_POINT; 148 const PetscInt Np = ex->periodic ? Nl : Nl + 1; 149 PetscInt Nc, No; 150 151 PetscFunctionBegin; 152 ex->Nt[ct] = 2; 153 Nc = 6 * Nl; 154 No = 2 * Nl; 155 PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct])); 156 ex->target[ct][0] = DM_POLYTOPE_POINT; 157 ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_POINT_PRISM_TENSOR : DM_POLYTOPE_SEGMENT; 158 ex->size[ct][0] = Np; 159 ex->size[ct][1] = Nl; 160 /* cones for segments/tensor segments */ 161 for (PetscInt i = 0; i < Nl; ++i) { 162 ex->cone[ct][6 * i + 0] = DM_POLYTOPE_POINT; 163 ex->cone[ct][6 * i + 1] = 0; 164 ex->cone[ct][6 * i + 2] = i; 165 ex->cone[ct][6 * i + 3] = DM_POLYTOPE_POINT; 166 ex->cone[ct][6 * i + 4] = 0; 167 ex->cone[ct][6 * i + 5] = (i + 1) % Np; 168 } 169 for (PetscInt i = 0; i < No; ++i) ex->ornt[ct][i] = 0; 170 PetscFunctionReturn(PETSC_SUCCESS); 171 } 172 173 /* DM_POLYTOPE_SEGMENT produces 174 Nl+1 segments, or Nl segments periodically and 175 Nl quads, or tensor quads 176 */ 177 static PetscErrorCode DMPlexTransformExtrudeSetUp_Segment(DMPlexTransform_Extrude *ex, PetscInt Nl) 178 { 179 const DMPolytopeType ct = DM_POLYTOPE_SEGMENT; 180 const PetscInt Np = ex->periodic ? Nl : Nl + 1; 181 PetscInt Nc, No, coff, ooff; 182 183 PetscFunctionBegin; 184 ex->Nt[ct] = 2; 185 Nc = 8 * Np + 14 * Nl; 186 No = 2 * Np + 4 * Nl; 187 PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct])); 188 ex->target[ct][0] = DM_POLYTOPE_SEGMENT; 189 ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_SEG_PRISM_TENSOR : DM_POLYTOPE_QUADRILATERAL; 190 ex->size[ct][0] = Np; 191 ex->size[ct][1] = Nl; 192 /* cones for segments */ 193 for (PetscInt i = 0; i < Np; ++i) { 194 ex->cone[ct][8 * i + 0] = DM_POLYTOPE_POINT; 195 ex->cone[ct][8 * i + 1] = 1; 196 ex->cone[ct][8 * i + 2] = 0; 197 ex->cone[ct][8 * i + 3] = i; 198 ex->cone[ct][8 * i + 4] = DM_POLYTOPE_POINT; 199 ex->cone[ct][8 * i + 5] = 1; 200 ex->cone[ct][8 * i + 6] = 1; 201 ex->cone[ct][8 * i + 7] = i; 202 } 203 for (PetscInt i = 0; i < 2 * Np; ++i) ex->ornt[ct][i] = 0; 204 /* cones for quads/tensor quads */ 205 coff = 8 * Np; 206 ooff = 2 * Np; 207 for (PetscInt i = 0; i < Nl; ++i) { 208 if (ex->useTensor) { 209 ex->cone[ct][coff + 14 * i + 0] = DM_POLYTOPE_SEGMENT; 210 ex->cone[ct][coff + 14 * i + 1] = 0; 211 ex->cone[ct][coff + 14 * i + 2] = i; 212 ex->cone[ct][coff + 14 * i + 3] = DM_POLYTOPE_SEGMENT; 213 ex->cone[ct][coff + 14 * i + 4] = 0; 214 ex->cone[ct][coff + 14 * i + 5] = (i + 1) % Np; 215 ex->cone[ct][coff + 14 * i + 6] = DM_POLYTOPE_POINT_PRISM_TENSOR; 216 ex->cone[ct][coff + 14 * i + 7] = 1; 217 ex->cone[ct][coff + 14 * i + 8] = 0; 218 ex->cone[ct][coff + 14 * i + 9] = i; 219 ex->cone[ct][coff + 14 * i + 10] = DM_POLYTOPE_POINT_PRISM_TENSOR; 220 ex->cone[ct][coff + 14 * i + 11] = 1; 221 ex->cone[ct][coff + 14 * i + 12] = 1; 222 ex->cone[ct][coff + 14 * i + 13] = i; 223 ex->ornt[ct][ooff + 4 * i + 0] = 0; 224 ex->ornt[ct][ooff + 4 * i + 1] = 0; 225 ex->ornt[ct][ooff + 4 * i + 2] = 0; 226 ex->ornt[ct][ooff + 4 * i + 3] = 0; 227 } else { 228 ex->cone[ct][coff + 14 * i + 0] = DM_POLYTOPE_SEGMENT; 229 ex->cone[ct][coff + 14 * i + 1] = 0; 230 ex->cone[ct][coff + 14 * i + 2] = i; 231 ex->cone[ct][coff + 14 * i + 3] = DM_POLYTOPE_SEGMENT; 232 ex->cone[ct][coff + 14 * i + 4] = 1; 233 ex->cone[ct][coff + 14 * i + 5] = 1; 234 ex->cone[ct][coff + 14 * i + 6] = i; 235 ex->cone[ct][coff + 14 * i + 7] = DM_POLYTOPE_SEGMENT; 236 ex->cone[ct][coff + 14 * i + 8] = 0; 237 ex->cone[ct][coff + 14 * i + 9] = (i + 1) % Np; 238 ex->cone[ct][coff + 14 * i + 10] = DM_POLYTOPE_SEGMENT; 239 ex->cone[ct][coff + 14 * i + 11] = 1; 240 ex->cone[ct][coff + 14 * i + 12] = 0; 241 ex->cone[ct][coff + 14 * i + 13] = i; 242 ex->ornt[ct][ooff + 4 * i + 0] = 0; 243 ex->ornt[ct][ooff + 4 * i + 1] = 0; 244 ex->ornt[ct][ooff + 4 * i + 2] = -1; 245 ex->ornt[ct][ooff + 4 * i + 3] = -1; 246 } 247 } 248 PetscFunctionReturn(PETSC_SUCCESS); 249 } 250 251 /* DM_POLYTOPE_TRIANGLE produces 252 Nl+1 triangles, or Nl triangles periodically and 253 Nl triangular prisms/tensor triangular prisms 254 */ 255 static PetscErrorCode DMPlexTransformExtrudeSetUp_Triangle(DMPlexTransform_Extrude *ex, PetscInt Nl) 256 { 257 const DMPolytopeType ct = DM_POLYTOPE_TRIANGLE; 258 const PetscInt Np = ex->periodic ? Nl : Nl + 1; 259 PetscInt Nc, No, coff, ooff; 260 261 PetscFunctionBegin; 262 ex->Nt[ct] = 2; 263 Nc = 12 * Np + 18 * Nl; 264 No = 3 * Np + 5 * Nl; 265 PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct])); 266 ex->target[ct][0] = DM_POLYTOPE_TRIANGLE; 267 ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_TRI_PRISM_TENSOR : DM_POLYTOPE_TRI_PRISM; 268 ex->size[ct][0] = Np; 269 ex->size[ct][1] = Nl; 270 /* cones for triangles */ 271 for (PetscInt i = 0; i < Np; ++i) { 272 ex->cone[ct][12 * i + 0] = DM_POLYTOPE_SEGMENT; 273 ex->cone[ct][12 * i + 1] = 1; 274 ex->cone[ct][12 * i + 2] = 0; 275 ex->cone[ct][12 * i + 3] = i; 276 ex->cone[ct][12 * i + 4] = DM_POLYTOPE_SEGMENT; 277 ex->cone[ct][12 * i + 5] = 1; 278 ex->cone[ct][12 * i + 6] = 1; 279 ex->cone[ct][12 * i + 7] = i; 280 ex->cone[ct][12 * i + 8] = DM_POLYTOPE_SEGMENT; 281 ex->cone[ct][12 * i + 9] = 1; 282 ex->cone[ct][12 * i + 10] = 2; 283 ex->cone[ct][12 * i + 11] = i; 284 } 285 for (PetscInt i = 0; i < 3 * Np; ++i) ex->ornt[ct][i] = 0; 286 /* cones for triangular prisms/tensor triangular prisms */ 287 coff = 12 * Np; 288 ooff = 3 * Np; 289 for (PetscInt i = 0; i < Nl; ++i) { 290 if (ex->useTensor) { 291 ex->cone[ct][coff + 18 * i + 0] = DM_POLYTOPE_TRIANGLE; 292 ex->cone[ct][coff + 18 * i + 1] = 0; 293 ex->cone[ct][coff + 18 * i + 2] = i; 294 ex->cone[ct][coff + 18 * i + 3] = DM_POLYTOPE_TRIANGLE; 295 ex->cone[ct][coff + 18 * i + 4] = 0; 296 ex->cone[ct][coff + 18 * i + 5] = (i + 1) % Np; 297 ex->cone[ct][coff + 18 * i + 6] = DM_POLYTOPE_SEG_PRISM_TENSOR; 298 ex->cone[ct][coff + 18 * i + 7] = 1; 299 ex->cone[ct][coff + 18 * i + 8] = 0; 300 ex->cone[ct][coff + 18 * i + 9] = i; 301 ex->cone[ct][coff + 18 * i + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR; 302 ex->cone[ct][coff + 18 * i + 11] = 1; 303 ex->cone[ct][coff + 18 * i + 12] = 1; 304 ex->cone[ct][coff + 18 * i + 13] = i; 305 ex->cone[ct][coff + 18 * i + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR; 306 ex->cone[ct][coff + 18 * i + 15] = 1; 307 ex->cone[ct][coff + 18 * i + 16] = 2; 308 ex->cone[ct][coff + 18 * i + 17] = i; 309 ex->ornt[ct][ooff + 5 * i + 0] = 0; 310 ex->ornt[ct][ooff + 5 * i + 1] = 0; 311 ex->ornt[ct][ooff + 5 * i + 2] = 0; 312 ex->ornt[ct][ooff + 5 * i + 3] = 0; 313 ex->ornt[ct][ooff + 5 * i + 4] = 0; 314 } else { 315 ex->cone[ct][coff + 18 * i + 0] = DM_POLYTOPE_TRIANGLE; 316 ex->cone[ct][coff + 18 * i + 1] = 0; 317 ex->cone[ct][coff + 18 * i + 2] = i; 318 ex->cone[ct][coff + 18 * i + 3] = DM_POLYTOPE_TRIANGLE; 319 ex->cone[ct][coff + 18 * i + 4] = 0; 320 ex->cone[ct][coff + 18 * i + 5] = (i + 1) % Np; 321 ex->cone[ct][coff + 18 * i + 6] = DM_POLYTOPE_QUADRILATERAL; 322 ex->cone[ct][coff + 18 * i + 7] = 1; 323 ex->cone[ct][coff + 18 * i + 8] = 0; 324 ex->cone[ct][coff + 18 * i + 9] = i; 325 ex->cone[ct][coff + 18 * i + 10] = DM_POLYTOPE_QUADRILATERAL; 326 ex->cone[ct][coff + 18 * i + 11] = 1; 327 ex->cone[ct][coff + 18 * i + 12] = 1; 328 ex->cone[ct][coff + 18 * i + 13] = i; 329 ex->cone[ct][coff + 18 * i + 14] = DM_POLYTOPE_QUADRILATERAL; 330 ex->cone[ct][coff + 18 * i + 15] = 1; 331 ex->cone[ct][coff + 18 * i + 16] = 2; 332 ex->cone[ct][coff + 18 * i + 17] = i; 333 ex->ornt[ct][ooff + 5 * i + 0] = -2; 334 ex->ornt[ct][ooff + 5 * i + 1] = 0; 335 ex->ornt[ct][ooff + 5 * i + 2] = 0; 336 ex->ornt[ct][ooff + 5 * i + 3] = 0; 337 ex->ornt[ct][ooff + 5 * i + 4] = 0; 338 } 339 } 340 PetscFunctionReturn(PETSC_SUCCESS); 341 } 342 343 /* DM_POLYTOPE_QUADRILATERAL produces 344 Nl+1 quads, or Nl quads periodically and 345 Nl hexes/tensor hexes 346 */ 347 static PetscErrorCode DMPlexTransformExtrudeSetUp_Quadrilateral(DMPlexTransform_Extrude *ex, PetscInt Nl) 348 { 349 const DMPolytopeType ct = DM_POLYTOPE_QUADRILATERAL; 350 const PetscInt Np = ex->periodic ? Nl : Nl + 1; 351 PetscInt Nc, No, coff, ooff; 352 353 PetscFunctionBegin; 354 ex->Nt[ct] = 2; 355 Nc = 16 * Np + 22 * Nl; 356 No = 4 * Np + 6 * Nl; 357 PetscCall(PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct])); 358 ex->target[ct][0] = DM_POLYTOPE_QUADRILATERAL; 359 ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_QUAD_PRISM_TENSOR : DM_POLYTOPE_HEXAHEDRON; 360 ex->size[ct][0] = Np; 361 ex->size[ct][1] = Nl; 362 /* cones for quads */ 363 for (PetscInt i = 0; i < Np; ++i) { 364 ex->cone[ct][16 * i + 0] = DM_POLYTOPE_SEGMENT; 365 ex->cone[ct][16 * i + 1] = 1; 366 ex->cone[ct][16 * i + 2] = 0; 367 ex->cone[ct][16 * i + 3] = i; 368 ex->cone[ct][16 * i + 4] = DM_POLYTOPE_SEGMENT; 369 ex->cone[ct][16 * i + 5] = 1; 370 ex->cone[ct][16 * i + 6] = 1; 371 ex->cone[ct][16 * i + 7] = i; 372 ex->cone[ct][16 * i + 8] = DM_POLYTOPE_SEGMENT; 373 ex->cone[ct][16 * i + 9] = 1; 374 ex->cone[ct][16 * i + 10] = 2; 375 ex->cone[ct][16 * i + 11] = i; 376 ex->cone[ct][16 * i + 12] = DM_POLYTOPE_SEGMENT; 377 ex->cone[ct][16 * i + 13] = 1; 378 ex->cone[ct][16 * i + 14] = 3; 379 ex->cone[ct][16 * i + 15] = i; 380 } 381 for (PetscInt i = 0; i < 4 * Np; ++i) ex->ornt[ct][i] = 0; 382 /* cones for hexes/tensor hexes */ 383 coff = 16 * Np; 384 ooff = 4 * Np; 385 for (PetscInt i = 0; i < Nl; ++i) { 386 if (ex->useTensor) { 387 ex->cone[ct][coff + 22 * i + 0] = DM_POLYTOPE_QUADRILATERAL; 388 ex->cone[ct][coff + 22 * i + 1] = 0; 389 ex->cone[ct][coff + 22 * i + 2] = i; 390 ex->cone[ct][coff + 22 * i + 3] = DM_POLYTOPE_QUADRILATERAL; 391 ex->cone[ct][coff + 22 * i + 4] = 0; 392 ex->cone[ct][coff + 22 * i + 5] = (i + 1) % Np; 393 ex->cone[ct][coff + 22 * i + 6] = DM_POLYTOPE_SEG_PRISM_TENSOR; 394 ex->cone[ct][coff + 22 * i + 7] = 1; 395 ex->cone[ct][coff + 22 * i + 8] = 0; 396 ex->cone[ct][coff + 22 * i + 9] = i; 397 ex->cone[ct][coff + 22 * i + 10] = DM_POLYTOPE_SEG_PRISM_TENSOR; 398 ex->cone[ct][coff + 22 * i + 11] = 1; 399 ex->cone[ct][coff + 22 * i + 12] = 1; 400 ex->cone[ct][coff + 22 * i + 13] = i; 401 ex->cone[ct][coff + 22 * i + 14] = DM_POLYTOPE_SEG_PRISM_TENSOR; 402 ex->cone[ct][coff + 22 * i + 15] = 1; 403 ex->cone[ct][coff + 22 * i + 16] = 2; 404 ex->cone[ct][coff + 22 * i + 17] = i; 405 ex->cone[ct][coff + 22 * i + 18] = DM_POLYTOPE_SEG_PRISM_TENSOR; 406 ex->cone[ct][coff + 22 * i + 19] = 1; 407 ex->cone[ct][coff + 22 * i + 20] = 3; 408 ex->cone[ct][coff + 22 * i + 21] = i; 409 ex->ornt[ct][ooff + 6 * i + 0] = 0; 410 ex->ornt[ct][ooff + 6 * i + 1] = 0; 411 ex->ornt[ct][ooff + 6 * i + 2] = 0; 412 ex->ornt[ct][ooff + 6 * i + 3] = 0; 413 ex->ornt[ct][ooff + 6 * i + 4] = 0; 414 ex->ornt[ct][ooff + 6 * i + 5] = 0; 415 } else { 416 ex->cone[ct][coff + 22 * i + 0] = DM_POLYTOPE_QUADRILATERAL; 417 ex->cone[ct][coff + 22 * i + 1] = 0; 418 ex->cone[ct][coff + 22 * i + 2] = i; 419 ex->cone[ct][coff + 22 * i + 3] = DM_POLYTOPE_QUADRILATERAL; 420 ex->cone[ct][coff + 22 * i + 4] = 0; 421 ex->cone[ct][coff + 22 * i + 5] = (i + 1) % Np; 422 ex->cone[ct][coff + 22 * i + 6] = DM_POLYTOPE_QUADRILATERAL; 423 ex->cone[ct][coff + 22 * i + 7] = 1; 424 ex->cone[ct][coff + 22 * i + 8] = 0; 425 ex->cone[ct][coff + 22 * i + 9] = i; 426 ex->cone[ct][coff + 22 * i + 10] = DM_POLYTOPE_QUADRILATERAL; 427 ex->cone[ct][coff + 22 * i + 11] = 1; 428 ex->cone[ct][coff + 22 * i + 12] = 2; 429 ex->cone[ct][coff + 22 * i + 13] = i; 430 ex->cone[ct][coff + 22 * i + 14] = DM_POLYTOPE_QUADRILATERAL; 431 ex->cone[ct][coff + 22 * i + 15] = 1; 432 ex->cone[ct][coff + 22 * i + 16] = 1; 433 ex->cone[ct][coff + 22 * i + 17] = i; 434 ex->cone[ct][coff + 22 * i + 18] = DM_POLYTOPE_QUADRILATERAL; 435 ex->cone[ct][coff + 22 * i + 19] = 1; 436 ex->cone[ct][coff + 22 * i + 20] = 3; 437 ex->cone[ct][coff + 22 * i + 21] = i; 438 ex->ornt[ct][ooff + 6 * i + 0] = -2; 439 ex->ornt[ct][ooff + 6 * i + 1] = 0; 440 ex->ornt[ct][ooff + 6 * i + 2] = 0; 441 ex->ornt[ct][ooff + 6 * i + 3] = 0; 442 ex->ornt[ct][ooff + 6 * i + 4] = 0; 443 ex->ornt[ct][ooff + 6 * i + 5] = 1; 444 } 445 } 446 PetscFunctionReturn(PETSC_SUCCESS); 447 } 448 449 /* 450 The refine types for extrusion are: 451 452 ct: For any normally extruded point 453 ct + 100: For any point which should just return itself 454 */ 455 static PetscErrorCode DMPlexTransformSetUp_Extrude(DMPlexTransform tr) 456 { 457 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 458 DM dm; 459 DMLabel active; 460 PetscInt Nl = ex->layers, l, ict, dim; 461 462 PetscFunctionBegin; 463 PetscCall(DMPlexTransformExtrudeComputeExtrusionDim(tr)); 464 PetscCall(DMPlexTransformGetDM(tr, &dm)); 465 PetscCall(DMGetDimension(dm, &dim)); 466 PetscCall(DMPlexTransformGetActive(tr, &active)); 467 if (active) { 468 DMLabel celltype; 469 PetscInt pStart, pEnd, p; 470 471 PetscCall(DMPlexGetCellTypeLabel(dm, &celltype)); 472 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Refine Type", &tr->trType)); 473 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 474 for (p = pStart; p < pEnd; ++p) { 475 PetscInt ct, val; 476 477 PetscCall(DMLabelGetValue(celltype, p, &ct)); 478 PetscCall(DMLabelGetValue(active, p, &val)); 479 if (val < 0) { 480 PetscCall(DMLabelSetValue(tr->trType, p, ct + 100)); 481 } else { 482 PetscCall(DMLabelSetValue(tr->trType, p, ct)); 483 } 484 } 485 } 486 if (ex->normalAlg != NORMAL_INPUT) { 487 if (dim != ex->cdim) ex->normalAlg = NORMAL_COMPUTE; 488 else if (active) ex->normalAlg = NORMAL_COMPUTE_BD; 489 } 490 // Need this to determine face sharing 491 PetscMPIInt size; 492 493 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size)); 494 if (size > 1) { 495 PetscSF sf; 496 PetscInt Nr; 497 498 PetscCall(DMGetPointSF(dm, &sf)); 499 PetscCall(PetscSFGetGraph(sf, &Nr, NULL, NULL, NULL)); 500 if (Nr >= 0) { 501 PetscCall(PetscSFComputeDegreeBegin(sf, &ex->degree)); 502 PetscCall(PetscSFComputeDegreeEnd(sf, &ex->degree)); 503 } 504 } 505 // Create normal field 506 if (ex->normalAlg == NORMAL_COMPUTE_BD) { 507 PetscSection s; 508 PetscInt vStart, vEnd; 509 510 PetscMPIInt rank; 511 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 512 513 PetscCall(DMClone(dm, &ex->dmNormal)); 514 PetscCall(DMGetLocalSection(ex->dmNormal, &s)); 515 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 516 PetscCall(PetscSectionSetNumFields(s, 1)); 517 PetscCall(PetscSectionSetChart(s, vStart, vEnd)); 518 for (PetscInt v = vStart; v < vEnd; ++v) PetscCall(PetscSectionSetDof(s, v, ex->cdimEx)); 519 PetscCall(PetscSectionSetUp(s)); 520 PetscCall(DMCreateLocalVector(ex->dmNormal, &ex->vecNormal)); 521 PetscCall(PetscObjectSetName((PetscObject)ex->vecNormal, "Normal Field")); 522 523 // find an active point in the closure of v and use its coordinate normal as the extrusion direction 524 PetscSF sf; 525 const PetscInt *leaves; 526 PetscScalar *a, *normal; 527 PetscInt Nl; 528 529 PetscCall(DMGetPointSF(ex->dmNormal, &sf)); 530 PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL)); 531 PetscCall(VecGetArrayWrite(ex->vecNormal, &a)); 532 for (PetscInt v = vStart; v < vEnd; ++v) { 533 PetscInt *star = NULL; 534 PetscInt starSize, pStart, pEnd; 535 536 PetscCall(DMPlexGetDepthStratum(ex->dmNormal, ex->cdimEx - 1, &pStart, &pEnd)); 537 PetscCall(DMPlexGetTransitiveClosure(ex->dmNormal, v, PETSC_FALSE, &starSize, &star)); 538 PetscCall(DMPlexPointLocalRef(ex->dmNormal, v, a, &normal)); 539 for (PetscInt st = 0; st < starSize * 2; st += 2) { 540 const PetscInt face = star[st]; 541 if ((face >= pStart) && (face < pEnd)) { 542 PetscReal cnormal[3] = {0, 0, 0}; 543 const PetscInt *supp; 544 PetscInt suppSize, floc = -1; 545 PetscBool shared; 546 547 PetscCall(DMPlexComputeCellGeometryFVM(ex->dmNormal, face, NULL, NULL, cnormal)); 548 PetscCall(DMPlexGetSupportSize(ex->dmNormal, face, &suppSize)); 549 PetscCall(DMPlexGetSupport(ex->dmNormal, face, &supp)); 550 // Only use external faces, so I can get the orientation from any cell 551 if (leaves) PetscCall(PetscFindInt(face, Nl, leaves, &floc)); 552 shared = floc >= 0 || (ex->degree && ex->degree[face]) ? PETSC_TRUE : PETSC_FALSE; 553 if (suppSize == 1 && !shared) { 554 const PetscInt *cone, *ornt; 555 PetscInt coneSize, c; 556 557 PetscCall(DMPlexGetConeSize(ex->dmNormal, supp[0], &coneSize)); 558 PetscCall(DMPlexGetCone(ex->dmNormal, supp[0], &cone)); 559 PetscCall(DMPlexGetConeOrientation(ex->dmNormal, supp[0], &ornt)); 560 for (c = 0; c < coneSize; ++c) 561 if (cone[c] == face) break; 562 PetscCheck(c < coneSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Asymmetry in cone/support"); 563 if (ornt[c] < 0) 564 for (PetscInt d = 0; d < ex->cdimEx; ++d) cnormal[d] *= -1.; 565 for (PetscInt d = 0; d < ex->cdimEx; ++d) normal[d] += cnormal[d]; 566 } 567 } 568 } 569 PetscCall(DMPlexRestoreTransitiveClosure(ex->dmNormal, v, PETSC_FALSE, &starSize, &star)); 570 } 571 PetscCall(VecRestoreArrayWrite(ex->vecNormal, &a)); 572 573 Vec g; 574 575 PetscCall(DMGetGlobalVector(ex->dmNormal, &g)); 576 PetscCall(VecSet(g, 0.)); 577 PetscCall(DMLocalToGlobal(ex->dmNormal, ex->vecNormal, ADD_VALUES, g)); 578 PetscCall(DMGlobalToLocal(ex->dmNormal, g, INSERT_VALUES, ex->vecNormal)); 579 PetscCall(DMRestoreGlobalVector(ex->dmNormal, &g)); 580 } 581 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)); 582 for (ict = 0; ict < DM_NUM_POLYTOPES; ++ict) { 583 ex->Nt[ict] = -1; 584 ex->target[ict] = NULL; 585 ex->size[ict] = NULL; 586 ex->cone[ict] = NULL; 587 ex->ornt[ict] = NULL; 588 } 589 PetscCall(DMPlexTransformExtrudeSetUp_Point(ex, Nl)); 590 PetscCall(DMPlexTransformExtrudeSetUp_Segment(ex, Nl)); 591 PetscCall(DMPlexTransformExtrudeSetUp_Triangle(ex, Nl)); 592 PetscCall(DMPlexTransformExtrudeSetUp_Quadrilateral(ex, Nl)); 593 /* Layers positions */ 594 if (!ex->Nth) { 595 if (ex->symmetric) 596 for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness * l) / ex->layers - ex->thickness / 2; 597 else 598 for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness * l) / ex->layers; 599 } else { 600 ex->thickness = 0.; 601 ex->layerPos[0] = 0.; 602 for (l = 0; l < ex->layers; ++l) { 603 const PetscReal t = ex->thicknesses[PetscMin(l, ex->Nth - 1)]; 604 ex->layerPos[l + 1] = ex->layerPos[l] + t; 605 ex->thickness += t; 606 } 607 if (ex->symmetric) 608 for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] -= ex->thickness / 2.; 609 } 610 PetscFunctionReturn(PETSC_SUCCESS); 611 } 612 613 static PetscErrorCode DMPlexTransformDestroy_Extrude(DMPlexTransform tr) 614 { 615 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 616 PetscInt ct; 617 618 PetscFunctionBegin; 619 if (ex->target) { 620 for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) PetscCall(PetscFree4(ex->target[ct], ex->size[ct], ex->cone[ct], ex->ornt[ct])); 621 } 622 PetscCall(PetscFree5(ex->Nt, ex->target, ex->size, ex->cone, ex->ornt)); 623 PetscCall(PetscFree(ex->layerPos)); 624 PetscCall(DMDestroy(&ex->dmNormal)); 625 PetscCall(VecDestroy(&ex->vecNormal)); 626 PetscCall(PetscFree(ex)); 627 PetscFunctionReturn(PETSC_SUCCESS); 628 } 629 630 static PetscErrorCode DMPlexTransformGetSubcellOrientation_Extrude(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 631 { 632 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 633 DMLabel trType = tr->trType, active; 634 PetscBool onBd = PETSC_FALSE; 635 PetscInt rt; 636 637 PetscFunctionBeginHot; 638 *rnew = r; 639 *onew = DMPolytopeTypeComposeOrientation(tct, o, so); 640 PetscCall(DMPlexTransformGetActive(tr, &active)); 641 if (!so && !active) PetscFunctionReturn(PETSC_SUCCESS); 642 if (trType) { 643 PetscCall(DMLabelGetValue(tr->trType, sp, &rt)); 644 if (rt >= 100) PetscFunctionReturn(PETSC_SUCCESS); 645 } 646 if (active) { 647 // Get orientation of boundary face in cell 648 if (DMPolytopeTypeGetDim(sct) == ex->dimEx - 1) { 649 DM dm; 650 const PetscInt *supp, *cone, *ornt; 651 PetscInt suppSize, coneSize, c; 652 653 PetscCall(DMPlexTransformGetDM(tr, &dm)); 654 PetscCall(DMPlexGetSupportSize(dm, sp, &suppSize)); 655 PetscCall(DMPlexGetSupport(dm, sp, &supp)); 656 PetscCheck(suppSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Source point %" PetscInt_FMT " is not a boundary face", sp); 657 PetscCall(DMPlexGetConeSize(dm, supp[0], &coneSize)); 658 PetscCall(DMPlexGetOrientedCone(dm, supp[0], &cone, &ornt)); 659 for (c = 0; c < coneSize; ++c) 660 if (cone[c] == sp) break; 661 PetscCheck(c < coneSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Source point %" PetscInt_FMT " not found in cone of support %" PetscInt_FMT, sp, supp[0]); 662 o = ornt[c]; 663 onBd = PETSC_TRUE; 664 PetscCall(DMPlexRestoreOrientedCone(dm, supp[0], &cone, &ornt)); 665 } 666 } 667 if (ex->useTensor) { 668 switch (sct) { 669 case DM_POLYTOPE_POINT: 670 break; 671 case DM_POLYTOPE_SEGMENT: 672 switch (tct) { 673 case DM_POLYTOPE_SEGMENT: 674 break; 675 case DM_POLYTOPE_SEG_PRISM_TENSOR: 676 *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -1) : DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0); 677 break; 678 default: 679 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 680 } 681 break; 682 // We need to handle identity extrusions from volumes (TET, HEX, etc) when boundary faces are being extruded 683 case DM_POLYTOPE_TRIANGLE: 684 switch (tct) { 685 case DM_POLYTOPE_TRIANGLE: 686 break; 687 case DM_POLYTOPE_TRI_PRISM_TENSOR: 688 *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -1) : DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0); 689 break; 690 default: 691 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 692 } 693 break; 694 case DM_POLYTOPE_QUADRILATERAL: 695 switch (tct) { 696 case DM_POLYTOPE_QUADRILATERAL: 697 break; 698 case DM_POLYTOPE_QUAD_PRISM_TENSOR: 699 *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -1) : DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0); 700 break; 701 default: 702 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 703 } 704 break; 705 default: 706 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]); 707 } 708 } else { 709 switch (sct) { 710 case DM_POLYTOPE_POINT: 711 break; 712 case DM_POLYTOPE_SEGMENT: 713 switch (tct) { 714 case DM_POLYTOPE_SEGMENT: 715 break; 716 case DM_POLYTOPE_QUADRILATERAL: 717 *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -3) : DMPolytopeTypeComposeOrientation(tct, o, so ? -3 : 0); 718 break; 719 default: 720 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 721 } 722 break; 723 case DM_POLYTOPE_TRIANGLE: 724 switch (tct) { 725 case DM_POLYTOPE_TRIANGLE: 726 break; 727 case DM_POLYTOPE_TRI_PRISM: 728 *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -1) : DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0); 729 break; 730 default: 731 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 732 } 733 break; 734 case DM_POLYTOPE_QUADRILATERAL: 735 switch (tct) { 736 case DM_POLYTOPE_QUADRILATERAL: 737 break; 738 case DM_POLYTOPE_HEXAHEDRON: 739 *onew = onBd ? DMPolytopeTypeComposeOrientation(tct, o, so ? 0 : -1) : DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0); 740 break; 741 default: 742 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 743 } 744 break; 745 default: 746 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]); 747 } 748 } 749 PetscFunctionReturn(PETSC_SUCCESS); 750 } 751 752 static PetscErrorCode DMPlexTransformCellTransform_Extrude(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 753 { 754 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 755 DMLabel trType = tr->trType; 756 PetscBool ignore = PETSC_FALSE, identity = PETSC_FALSE; 757 PetscInt val = 0; 758 759 PetscFunctionBegin; 760 if (trType) { 761 PetscCall(DMLabelGetValue(trType, p, &val)); 762 identity = val >= 100 ? PETSC_TRUE : PETSC_FALSE; 763 } else { 764 ignore = ex->Nt[source] < 0 ? PETSC_TRUE : PETSC_FALSE; 765 } 766 if (rt) *rt = val; 767 if (ignore) { 768 /* Ignore cells that cannot be extruded */ 769 *Nt = 0; 770 *target = NULL; 771 *size = NULL; 772 *cone = NULL; 773 *ornt = NULL; 774 } else if (identity) { 775 PetscCall(DMPlexTransformCellTransformIdentity(tr, source, p, NULL, Nt, target, size, cone, ornt)); 776 } else { 777 *Nt = ex->Nt[source]; 778 *target = ex->target[source]; 779 *size = ex->size[source]; 780 *cone = ex->cone[source]; 781 *ornt = ex->ornt[source]; 782 } 783 PetscFunctionReturn(PETSC_SUCCESS); 784 } 785 786 /* Computes new vertex along normal */ 787 static PetscErrorCode DMPlexTransformMapCoordinates_Extrude(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) 788 { 789 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 790 DM dm; 791 DMLabel active; 792 PetscReal ones2[2] = {0., 1.}; 793 PetscReal ones3[3] = {0., 0., 1.}; 794 PetscReal normal[3] = {0., 0., 0.}; 795 PetscReal norm = 0.; 796 PetscInt dEx = ex->cdimEx; 797 PetscInt dim, cStart, cEnd; 798 799 PetscFunctionBeginHot; 800 PetscCheck(pct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for parent point type %s", DMPolytopeTypes[pct]); 801 PetscCheck(ct == DM_POLYTOPE_POINT, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for refined point type %s", DMPolytopeTypes[ct]); 802 PetscCheck(Nv == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Vertices should be produced from a single vertex, not %" PetscInt_FMT, Nv); 803 PetscCheck(dE == ex->cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Coordinate dim %" PetscInt_FMT " != %" PetscInt_FMT " original dimension", dE, ex->cdim); 804 PetscCheck(dEx <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension for extruded mesh %" PetscInt_FMT " must not exceed 3", dEx); 805 806 PetscCall(DMPlexTransformGetDM(tr, &dm)); 807 PetscCall(DMGetDimension(dm, &dim)); 808 PetscCall(DMPlexTransformGetActive(tr, &active)); 809 switch (ex->normalAlg) { 810 case NORMAL_DEFAULT: 811 switch (ex->cdimEx) { 812 case 2: 813 for (PetscInt d = 0; d < dEx; ++d) normal[d] = ones2[d]; 814 break; 815 case 3: 816 for (PetscInt d = 0; d < dEx; ++d) normal[d] = ones3[d]; 817 break; 818 default: 819 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No default normal for dimension %" PetscInt_FMT, ex->cdimEx); 820 } 821 break; 822 case NORMAL_INPUT: 823 for (PetscInt d = 0; d < dEx; ++d) normal[d] = ex->normal[d]; 824 break; 825 case NORMAL_COMPUTE: { 826 PetscInt *star = NULL; 827 PetscInt starSize; 828 829 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 830 PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &starSize, &star)); 831 for (PetscInt st = 0; st < starSize * 2; st += 2) { 832 if ((star[st] >= cStart) && (star[st] < cEnd)) { 833 PetscReal cnormal[3] = {0, 0, 0}; 834 835 PetscCall(DMPlexComputeCellGeometryFVM(dm, star[st], NULL, NULL, cnormal)); 836 for (PetscInt d = 0; d < dEx; ++d) normal[d] += cnormal[d]; 837 } 838 } 839 PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &starSize, &star)); 840 } break; 841 case NORMAL_COMPUTE_BD: { 842 const PetscScalar *a; 843 PetscScalar *vnormal; 844 845 PetscCall(VecGetArrayRead(ex->vecNormal, &a)); 846 PetscCall(DMPlexPointLocalRead(ex->dmNormal, p, a, (void *)&vnormal)); 847 for (PetscInt d = 0; d < dEx; ++d) normal[d] = PetscRealPart(vnormal[d]); 848 PetscCall(VecRestoreArrayRead(ex->vecNormal, &a)); 849 } break; 850 default: 851 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to determine normal for extrusion"); 852 } 853 if (ex->normalFunc) { 854 PetscScalar n[3]; 855 PetscReal x[3], dot = 0.; 856 857 for (PetscInt d = 0; d < ex->cdim; ++d) x[d] = PetscRealPart(in[d]); 858 PetscCall((*ex->normalFunc)(ex->cdim, 0., x, r, n, NULL)); 859 for (PetscInt d = 0; d < dEx; ++d) dot += PetscRealPart(n[d]) * normal[d]; 860 for (PetscInt d = 0; d < dEx; ++d) normal[d] = PetscSign(dot) * PetscRealPart(n[d]); 861 } 862 for (PetscInt d = 0; d < dEx; ++d) norm += PetscSqr(normal[d]); 863 for (PetscInt d = 0; d < dEx; ++d) normal[d] *= norm == 0.0 ? 1.0 : 1. / PetscSqrtReal(norm); 864 for (PetscInt d = 0; d < dEx; ++d) out[d] = normal[d] * ex->layerPos[r]; 865 for (PetscInt d = 0; d < ex->cdim; ++d) out[d] += in[d]; 866 PetscFunctionReturn(PETSC_SUCCESS); 867 } 868 869 static PetscErrorCode DMPlexTransformInitialize_Extrude(DMPlexTransform tr) 870 { 871 PetscFunctionBegin; 872 tr->ops->view = DMPlexTransformView_Extrude; 873 tr->ops->setfromoptions = DMPlexTransformSetFromOptions_Extrude; 874 tr->ops->setup = DMPlexTransformSetUp_Extrude; 875 tr->ops->destroy = DMPlexTransformDestroy_Extrude; 876 tr->ops->setdimensions = DMPlexTransformSetDimensions_Extrude; 877 tr->ops->celltransform = DMPlexTransformCellTransform_Extrude; 878 tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Extrude; 879 tr->ops->mapcoordinates = DMPlexTransformMapCoordinates_Extrude; 880 PetscFunctionReturn(PETSC_SUCCESS); 881 } 882 883 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform tr) 884 { 885 DMPlexTransform_Extrude *ex; 886 DM dm; 887 PetscInt dim; 888 889 PetscFunctionBegin; 890 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 891 PetscCall(PetscNew(&ex)); 892 tr->data = ex; 893 ex->thickness = 1.; 894 ex->useTensor = PETSC_TRUE; 895 ex->symmetric = PETSC_FALSE; 896 ex->periodic = PETSC_FALSE; 897 ex->normalAlg = NORMAL_DEFAULT; 898 ex->layerPos = NULL; 899 PetscCall(DMPlexTransformGetDM(tr, &dm)); 900 PetscCall(DMGetDimension(dm, &dim)); 901 PetscCall(DMGetCoordinateDim(dm, &ex->cdim)); 902 PetscCall(DMPlexTransformExtrudeSetLayers(tr, 1)); 903 PetscCall(DMPlexTransformInitialize_Extrude(tr)); 904 PetscFunctionReturn(PETSC_SUCCESS); 905 } 906 907 /*@ 908 DMPlexTransformExtrudeGetLayers - Get the number of extruded layers. 909 910 Not Collective 911 912 Input Parameter: 913 . tr - The `DMPlexTransform` 914 915 Output Parameter: 916 . layers - The number of layers 917 918 Level: intermediate 919 920 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetLayers()` 921 @*/ 922 PetscErrorCode DMPlexTransformExtrudeGetLayers(DMPlexTransform tr, PetscInt *layers) 923 { 924 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 925 926 PetscFunctionBegin; 927 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 928 PetscAssertPointer(layers, 2); 929 *layers = ex->layers; 930 PetscFunctionReturn(PETSC_SUCCESS); 931 } 932 933 /*@ 934 DMPlexTransformExtrudeSetLayers - Set the number of extruded layers. 935 936 Not Collective 937 938 Input Parameters: 939 + tr - The `DMPlexTransform` 940 - layers - The number of layers 941 942 Level: intermediate 943 944 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetLayers()` 945 @*/ 946 PetscErrorCode DMPlexTransformExtrudeSetLayers(DMPlexTransform tr, PetscInt layers) 947 { 948 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 949 950 PetscFunctionBegin; 951 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 952 ex->layers = layers; 953 PetscCall(PetscFree(ex->layerPos)); 954 PetscCall(PetscCalloc1(ex->layers + 1, &ex->layerPos)); 955 PetscFunctionReturn(PETSC_SUCCESS); 956 } 957 958 /*@ 959 DMPlexTransformExtrudeGetThickness - Get the total thickness of the layers 960 961 Not Collective 962 963 Input Parameter: 964 . tr - The `DMPlexTransform` 965 966 Output Parameter: 967 . thickness - The total thickness of the layers 968 969 Level: intermediate 970 971 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetThickness()` 972 @*/ 973 PetscErrorCode DMPlexTransformExtrudeGetThickness(DMPlexTransform tr, PetscReal *thickness) 974 { 975 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 976 977 PetscFunctionBegin; 978 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 979 PetscAssertPointer(thickness, 2); 980 *thickness = ex->thickness; 981 PetscFunctionReturn(PETSC_SUCCESS); 982 } 983 984 /*@ 985 DMPlexTransformExtrudeSetThickness - Set the total thickness of the layers 986 987 Not Collective 988 989 Input Parameters: 990 + tr - The `DMPlexTransform` 991 - thickness - The total thickness of the layers 992 993 Level: intermediate 994 995 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetThickness()` 996 @*/ 997 PetscErrorCode DMPlexTransformExtrudeSetThickness(DMPlexTransform tr, PetscReal thickness) 998 { 999 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 1000 1001 PetscFunctionBegin; 1002 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1003 PetscCheck(thickness > 0., PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Height of layers %g must be positive", (double)thickness); 1004 ex->thickness = thickness; 1005 PetscFunctionReturn(PETSC_SUCCESS); 1006 } 1007 1008 /*@ 1009 DMPlexTransformExtrudeGetTensor - Get the flag to use tensor cells 1010 1011 Not Collective 1012 1013 Input Parameter: 1014 . tr - The `DMPlexTransform` 1015 1016 Output Parameter: 1017 . useTensor - The flag to use tensor cells 1018 1019 Note: 1020 This flag determines the orientation behavior of the created points. 1021 1022 For example, if tensor is `PETSC_TRUE`, then 1023 .vb 1024 DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT, 1025 DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL, 1026 DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and 1027 DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON. 1028 .ve 1029 1030 Level: intermediate 1031 1032 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetTensor()` 1033 @*/ 1034 PetscErrorCode DMPlexTransformExtrudeGetTensor(DMPlexTransform tr, PetscBool *useTensor) 1035 { 1036 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 1037 1038 PetscFunctionBegin; 1039 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1040 PetscAssertPointer(useTensor, 2); 1041 *useTensor = ex->useTensor; 1042 PetscFunctionReturn(PETSC_SUCCESS); 1043 } 1044 1045 /*@ 1046 DMPlexTransformExtrudeSetTensor - Set the flag to use tensor cells 1047 1048 Not Collective 1049 1050 Input Parameters: 1051 + tr - The `DMPlexTransform` 1052 - useTensor - The flag for tensor cells 1053 1054 Note: 1055 This flag determines the orientation behavior of the created points 1056 For example, if tensor is `PETSC_TRUE`, then 1057 .vb 1058 DM_POLYTOPE_POINT_PRISM_TENSOR is made instead of DM_POLYTOPE_SEGMENT, 1059 DM_POLYTOPE_SEG_PRISM_TENSOR instead of DM_POLYTOPE_QUADRILATERAL, 1060 DM_POLYTOPE_TRI_PRISM_TENSOR instead of DM_POLYTOPE_TRI_PRISM, and 1061 DM_POLYTOPE_QUAD_PRISM_TENSOR instead of DM_POLYTOPE_HEXAHEDRON. 1062 .ve 1063 1064 Level: intermediate 1065 1066 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetTensor()` 1067 @*/ 1068 PetscErrorCode DMPlexTransformExtrudeSetTensor(DMPlexTransform tr, PetscBool useTensor) 1069 { 1070 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 1071 1072 PetscFunctionBegin; 1073 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1074 ex->useTensor = useTensor; 1075 PetscFunctionReturn(PETSC_SUCCESS); 1076 } 1077 1078 /*@ 1079 DMPlexTransformExtrudeGetSymmetric - Get the flag to extrude symmetrically from the initial surface 1080 1081 Not Collective 1082 1083 Input Parameter: 1084 . tr - The `DMPlexTransform` 1085 1086 Output Parameter: 1087 . symmetric - The flag to extrude symmetrically 1088 1089 Level: intermediate 1090 1091 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetSymmetric()`, `DMPlexTransformExtrudeGetPeriodic()` 1092 @*/ 1093 PetscErrorCode DMPlexTransformExtrudeGetSymmetric(DMPlexTransform tr, PetscBool *symmetric) 1094 { 1095 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 1096 1097 PetscFunctionBegin; 1098 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1099 PetscAssertPointer(symmetric, 2); 1100 *symmetric = ex->symmetric; 1101 PetscFunctionReturn(PETSC_SUCCESS); 1102 } 1103 1104 /*@ 1105 DMPlexTransformExtrudeSetSymmetric - Set the flag to extrude symmetrically from the initial surface 1106 1107 Not Collective 1108 1109 Input Parameters: 1110 + tr - The `DMPlexTransform` 1111 - symmetric - The flag to extrude symmetrically 1112 1113 Level: intermediate 1114 1115 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetSymmetric()`, `DMPlexTransformExtrudeSetPeriodic()` 1116 @*/ 1117 PetscErrorCode DMPlexTransformExtrudeSetSymmetric(DMPlexTransform tr, PetscBool symmetric) 1118 { 1119 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 1120 1121 PetscFunctionBegin; 1122 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1123 ex->symmetric = symmetric; 1124 PetscFunctionReturn(PETSC_SUCCESS); 1125 } 1126 1127 /*@ 1128 DMPlexTransformExtrudeGetPeriodic - Get the flag to extrude periodically from the initial surface 1129 1130 Not Collective 1131 1132 Input Parameter: 1133 . tr - The `DMPlexTransform` 1134 1135 Output Parameter: 1136 . periodic - The flag to extrude periodically 1137 1138 Level: intermediate 1139 1140 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetPeriodic()`, `DMPlexTransformExtrudeGetSymmetric()` 1141 @*/ 1142 PetscErrorCode DMPlexTransformExtrudeGetPeriodic(DMPlexTransform tr, PetscBool *periodic) 1143 { 1144 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 1145 1146 PetscFunctionBegin; 1147 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1148 PetscAssertPointer(periodic, 2); 1149 *periodic = ex->periodic; 1150 PetscFunctionReturn(PETSC_SUCCESS); 1151 } 1152 1153 /*@ 1154 DMPlexTransformExtrudeSetPeriodic - Set the flag to extrude periodically from the initial surface 1155 1156 Not Collective 1157 1158 Input Parameters: 1159 + tr - The `DMPlexTransform` 1160 - periodic - The flag to extrude periodically 1161 1162 Level: intermediate 1163 1164 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetPeriodic()`, `DMPlexTransformExtrudeSetSymmetric()` 1165 @*/ 1166 PetscErrorCode DMPlexTransformExtrudeSetPeriodic(DMPlexTransform tr, PetscBool periodic) 1167 { 1168 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 1169 1170 PetscFunctionBegin; 1171 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1172 ex->periodic = periodic; 1173 PetscFunctionReturn(PETSC_SUCCESS); 1174 } 1175 1176 /*@ 1177 DMPlexTransformExtrudeGetNormal - Get the extrusion normal vector 1178 1179 Not Collective 1180 1181 Input Parameter: 1182 . tr - The `DMPlexTransform` 1183 1184 Output Parameter: 1185 . normal - The extrusion direction 1186 1187 Note: 1188 The user passes in an array, which is filled by the library. 1189 1190 Level: intermediate 1191 1192 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetNormal()` 1193 @*/ 1194 PetscErrorCode DMPlexTransformExtrudeGetNormal(DMPlexTransform tr, PetscReal normal[]) 1195 { 1196 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 1197 PetscInt d; 1198 1199 PetscFunctionBegin; 1200 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1201 if (ex->normalAlg == NORMAL_INPUT) { 1202 for (d = 0; d < ex->cdimEx; ++d) normal[d] = ex->normal[d]; 1203 } else { 1204 for (d = 0; d < ex->cdimEx; ++d) normal[d] = 0.; 1205 } 1206 PetscFunctionReturn(PETSC_SUCCESS); 1207 } 1208 1209 /*@ 1210 DMPlexTransformExtrudeSetNormal - Set the extrusion normal 1211 1212 Not Collective 1213 1214 Input Parameters: 1215 + tr - The `DMPlexTransform` 1216 - normal - The extrusion direction 1217 1218 Level: intermediate 1219 1220 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetNormal()` 1221 @*/ 1222 PetscErrorCode DMPlexTransformExtrudeSetNormal(DMPlexTransform tr, const PetscReal normal[]) 1223 { 1224 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 1225 PetscInt d; 1226 1227 PetscFunctionBegin; 1228 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1229 ex->normalAlg = NORMAL_INPUT; 1230 if (!ex->cdimEx) PetscCall(DMPlexTransformExtrudeComputeExtrusionDim(tr)); 1231 for (d = 0; d < ex->cdimEx; ++d) ex->normal[d] = normal[d]; 1232 PetscFunctionReturn(PETSC_SUCCESS); 1233 } 1234 1235 /*@C 1236 DMPlexTransformExtrudeSetNormalFunction - Set a function to determine the extrusion normal 1237 1238 Not Collective 1239 1240 Input Parameters: 1241 + tr - The `DMPlexTransform` 1242 - normalFunc - A function determining the extrusion direction, see `PetscSimplePointFn` for the calling sequence 1243 1244 Level: intermediate 1245 1246 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeGetNormal()`, `PetscSimplePointFn` 1247 @*/ 1248 PetscErrorCode DMPlexTransformExtrudeSetNormalFunction(DMPlexTransform tr, PetscSimplePointFn *normalFunc) 1249 { 1250 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 1251 1252 PetscFunctionBegin; 1253 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1254 ex->normalFunc = normalFunc; 1255 PetscFunctionReturn(PETSC_SUCCESS); 1256 } 1257 1258 /*@ 1259 DMPlexTransformExtrudeSetThicknesses - Set the thickness of each layer 1260 1261 Not Collective 1262 1263 Input Parameters: 1264 + tr - The `DMPlexTransform` 1265 . Nth - The number of thicknesses 1266 - thicknesses - The array of thicknesses 1267 1268 Level: intermediate 1269 1270 .seealso: `DMPlexTransform`, `DMPlexTransformExtrudeSetThickness()`, `DMPlexTransformExtrudeGetThickness()` 1271 @*/ 1272 PetscErrorCode DMPlexTransformExtrudeSetThicknesses(DMPlexTransform tr, PetscInt Nth, const PetscReal thicknesses[]) 1273 { 1274 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *)tr->data; 1275 PetscInt t; 1276 1277 PetscFunctionBegin; 1278 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 1279 PetscCheck(Nth > 0, PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Number of thicknesses %" PetscInt_FMT " must be positive", Nth); 1280 ex->Nth = PetscMin(Nth, ex->layers); 1281 PetscCall(PetscFree(ex->thicknesses)); 1282 PetscCall(PetscMalloc1(ex->Nth, &ex->thicknesses)); 1283 for (t = 0; t < ex->Nth; ++t) { 1284 PetscCheck(thicknesses[t] > 0., PetscObjectComm((PetscObject)tr), PETSC_ERR_ARG_OUTOFRANGE, "Thickness %g of layer %" PetscInt_FMT " must be positive", (double)thicknesses[t], t); 1285 ex->thicknesses[t] = thicknesses[t]; 1286 } 1287 PetscFunctionReturn(PETSC_SUCCESS); 1288 } 1289