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