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