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