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 PetscErrorCode ierr; 8 9 PetscFunctionBegin; 10 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 11 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 12 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr); 13 if (isascii) { 14 const char *name; 15 16 ierr = PetscObjectGetName((PetscObject) tr, &name);CHKERRQ(ierr); 17 ierr = PetscViewerASCIIPrintf(viewer, "Extrusion transformation %s\n", name ? name : "");CHKERRQ(ierr); 18 ierr = PetscViewerASCIIPrintf(viewer, " number of layers: %D\n", ex->layers);CHKERRQ(ierr); 19 ierr = PetscViewerASCIIPrintf(viewer, " create tensor cells: %s\n", ex->useTensor ? "YES" : "NO");CHKERRQ(ierr); 20 } else { 21 SETERRQ(PetscObjectComm((PetscObject) tr), PETSC_ERR_SUP, "Viewer type %s not yet supported for DMPlexTransform writing", ((PetscObject) viewer)->type_name); 22 } 23 PetscFunctionReturn(0); 24 } 25 26 static PetscErrorCode DMPlexTransformSetFromOptions_Extrude(PetscOptionItems *PetscOptionsObject, DMPlexTransform tr) 27 { 28 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 29 PetscReal th, normal[3], *thicknesses; 30 PetscInt nl, Nc; 31 PetscBool tensor, sym, flg; 32 PetscErrorCode ierr; 33 34 PetscFunctionBegin; 35 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 2); 36 ierr = PetscOptionsHead(PetscOptionsObject, "DMPlexTransform Extrusion Options");CHKERRQ(ierr); 37 ierr = PetscOptionsBoundedInt("-dm_plex_transform_extrude_layers", "Number of layers to extrude", "", ex->layers, &nl, &flg, 1);CHKERRQ(ierr); 38 if (flg) {ierr = DMPlexTransformExtrudeSetLayers(tr, nl);CHKERRQ(ierr);} 39 ierr = PetscOptionsReal("-dm_plex_transform_extrude_thickness", "Total thickness of extruded layers", "", ex->thickness, &th, &flg);CHKERRQ(ierr); 40 if (flg) {ierr = DMPlexTransformExtrudeSetThickness(tr, th);CHKERRQ(ierr);} 41 ierr = PetscOptionsBool("-dm_plex_transform_extrude_use_tensor", "Create tensor cells", "", ex->useTensor, &tensor, &flg);CHKERRQ(ierr); 42 if (flg) {ierr = DMPlexTransformExtrudeSetTensor(tr, tensor);CHKERRQ(ierr);} 43 ierr = PetscOptionsBool("-dm_plex_transform_extrude_symmetric", "Extrude layers symmetrically about the surface", "", ex->symmetric, &sym, &flg);CHKERRQ(ierr); 44 if (flg) {ierr = DMPlexTransformExtrudeSetSymmetric(tr, sym);CHKERRQ(ierr);} 45 Nc = 3; 46 ierr = PetscOptionsRealArray("-dm_plex_transform_extrude_normal", "Input normal vector for extrusion", "", normal, &Nc, &flg);CHKERRQ(ierr); 47 if (flg) { 48 PetscCheckFalse(Nc != ex->cdimEx,PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_SIZ, "Input normal has size %D != %D extruded coordinate dimension", Nc, ex->cdimEx); 49 ierr = DMPlexTransformExtrudeSetNormal(tr, normal);CHKERRQ(ierr); 50 } 51 nl = ex->layers; 52 ierr = PetscMalloc1(nl, &thicknesses);CHKERRQ(ierr); 53 ierr = PetscOptionsRealArray("-dm_plex_transform_extrude_thicknesses", "Thickness of each individual extruded layer", "", thicknesses, &nl, &flg);CHKERRQ(ierr); 54 if (flg) { 55 PetscCheckFalse(!nl,PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_OUTOFRANGE, "Must give at least one thickness for -dm_plex_transform_extrude_thicknesses"); 56 ierr = DMPlexTransformExtrudeSetThicknesses(tr, nl, thicknesses);CHKERRQ(ierr); 57 } 58 ierr = PetscFree(thicknesses);CHKERRQ(ierr); 59 ierr = PetscOptionsTail();CHKERRQ(ierr); 60 PetscFunctionReturn(0); 61 } 62 63 static PetscErrorCode DMPlexTransformSetDimensions_Extrude(DMPlexTransform tr, DM dm, DM tdm) 64 { 65 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 66 PetscInt dim; 67 PetscErrorCode ierr; 68 69 PetscFunctionBegin; 70 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 71 ierr = DMSetDimension(tdm, dim+1);CHKERRQ(ierr); 72 ierr = DMSetCoordinateDim(tdm, ex->cdimEx);CHKERRQ(ierr); 73 PetscFunctionReturn(0); 74 } 75 76 static PetscErrorCode DMPlexTransformSetUp_Extrude(DMPlexTransform tr) 77 { 78 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 79 DM dm; 80 DMPolytopeType ct; 81 PetscInt Nl = ex->layers, l, i, ict, Nc, No, coff, ooff; 82 PetscErrorCode ierr; 83 84 PetscFunctionBegin; 85 ierr = DMPlexTransformGetDM(tr, &dm);CHKERRQ(ierr); 86 ierr = 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);CHKERRQ(ierr); 87 for (ict = 0; ict < DM_NUM_POLYTOPES; ++ict) { 88 ex->Nt[ict] = -1; 89 ex->target[ict] = NULL; 90 ex->size[ict] = NULL; 91 ex->cone[ict] = NULL; 92 ex->ornt[ict] = NULL; 93 } 94 /* DM_POLYTOPE_POINT produces Nl+1 points and Nl segments/tensor segments */ 95 ct = DM_POLYTOPE_POINT; 96 ex->Nt[ct] = 2; 97 Nc = 6*Nl; 98 No = 2*Nl; 99 ierr = PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]);CHKERRQ(ierr); 100 ex->target[ct][0] = DM_POLYTOPE_POINT; 101 ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_POINT_PRISM_TENSOR : DM_POLYTOPE_SEGMENT; 102 ex->size[ct][0] = Nl+1; 103 ex->size[ct][1] = Nl; 104 /* cones for segments/tensor segments */ 105 for (i = 0; i < Nl; ++i) { 106 ex->cone[ct][6*i+0] = DM_POLYTOPE_POINT; 107 ex->cone[ct][6*i+1] = 0; 108 ex->cone[ct][6*i+2] = i; 109 ex->cone[ct][6*i+3] = DM_POLYTOPE_POINT; 110 ex->cone[ct][6*i+4] = 0; 111 ex->cone[ct][6*i+5] = i+1; 112 } 113 for (i = 0; i < No; ++i) ex->ornt[ct][i] = 0; 114 /* DM_POLYTOPE_SEGMENT produces Nl+1 segments and Nl quads/tensor quads */ 115 ct = DM_POLYTOPE_SEGMENT; 116 ex->Nt[ct] = 2; 117 Nc = 8*(Nl+1) + 14*Nl; 118 No = 2*(Nl+1) + 4*Nl; 119 ierr = PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]);CHKERRQ(ierr); 120 ex->target[ct][0] = DM_POLYTOPE_SEGMENT; 121 ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_SEG_PRISM_TENSOR : DM_POLYTOPE_QUADRILATERAL; 122 ex->size[ct][0] = Nl+1; 123 ex->size[ct][1] = Nl; 124 /* cones for segments */ 125 for (i = 0; i < Nl+1; ++i) { 126 ex->cone[ct][8*i+0] = DM_POLYTOPE_POINT; 127 ex->cone[ct][8*i+1] = 1; 128 ex->cone[ct][8*i+2] = 0; 129 ex->cone[ct][8*i+3] = i; 130 ex->cone[ct][8*i+4] = DM_POLYTOPE_POINT; 131 ex->cone[ct][8*i+5] = 1; 132 ex->cone[ct][8*i+6] = 1; 133 ex->cone[ct][8*i+7] = i; 134 } 135 for (i = 0; i < 2*(Nl+1); ++i) ex->ornt[ct][i] = 0; 136 /* cones for quads/tensor quads */ 137 coff = 8*(Nl+1); 138 ooff = 2*(Nl+1); 139 for (i = 0; i < Nl; ++i) { 140 if (ex->useTensor) { 141 ex->cone[ct][coff+14*i+0] = DM_POLYTOPE_SEGMENT; 142 ex->cone[ct][coff+14*i+1] = 0; 143 ex->cone[ct][coff+14*i+2] = i; 144 ex->cone[ct][coff+14*i+3] = DM_POLYTOPE_SEGMENT; 145 ex->cone[ct][coff+14*i+4] = 0; 146 ex->cone[ct][coff+14*i+5] = i+1; 147 ex->cone[ct][coff+14*i+6] = DM_POLYTOPE_POINT_PRISM_TENSOR; 148 ex->cone[ct][coff+14*i+7] = 1; 149 ex->cone[ct][coff+14*i+8] = 0; 150 ex->cone[ct][coff+14*i+9] = i; 151 ex->cone[ct][coff+14*i+10] = DM_POLYTOPE_POINT_PRISM_TENSOR; 152 ex->cone[ct][coff+14*i+11] = 1; 153 ex->cone[ct][coff+14*i+12] = 1; 154 ex->cone[ct][coff+14*i+13] = i; 155 ex->ornt[ct][ooff+4*i+0] = 0; 156 ex->ornt[ct][ooff+4*i+1] = 0; 157 ex->ornt[ct][ooff+4*i+2] = 0; 158 ex->ornt[ct][ooff+4*i+3] = 0; 159 } else { 160 ex->cone[ct][coff+14*i+0] = DM_POLYTOPE_SEGMENT; 161 ex->cone[ct][coff+14*i+1] = 0; 162 ex->cone[ct][coff+14*i+2] = i; 163 ex->cone[ct][coff+14*i+3] = DM_POLYTOPE_SEGMENT; 164 ex->cone[ct][coff+14*i+4] = 1; 165 ex->cone[ct][coff+14*i+5] = 1; 166 ex->cone[ct][coff+14*i+6] = i; 167 ex->cone[ct][coff+14*i+7] = DM_POLYTOPE_SEGMENT; 168 ex->cone[ct][coff+14*i+8] = 0; 169 ex->cone[ct][coff+14*i+9] = i+1; 170 ex->cone[ct][coff+14*i+10] = DM_POLYTOPE_SEGMENT; 171 ex->cone[ct][coff+14*i+11] = 1; 172 ex->cone[ct][coff+14*i+12] = 0; 173 ex->cone[ct][coff+14*i+13] = i; 174 ex->ornt[ct][ooff+4*i+0] = 0; 175 ex->ornt[ct][ooff+4*i+1] = 0; 176 ex->ornt[ct][ooff+4*i+2] = -1; 177 ex->ornt[ct][ooff+4*i+3] = -1; 178 } 179 } 180 /* DM_POLYTOPE_TRIANGLE produces Nl+1 triangles and Nl triangular prisms/tensor triangular prisms */ 181 ct = DM_POLYTOPE_TRIANGLE; 182 ex->Nt[ct] = 2; 183 Nc = 12*(Nl+1) + 18*Nl; 184 No = 3*(Nl+1) + 5*Nl; 185 ierr = PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]);CHKERRQ(ierr); 186 ex->target[ct][0] = DM_POLYTOPE_TRIANGLE; 187 ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_TRI_PRISM_TENSOR : DM_POLYTOPE_TRI_PRISM; 188 ex->size[ct][0] = Nl+1; 189 ex->size[ct][1] = Nl; 190 /* cones for triangles */ 191 for (i = 0; i < Nl+1; ++i) { 192 ex->cone[ct][12*i+0] = DM_POLYTOPE_SEGMENT; 193 ex->cone[ct][12*i+1] = 1; 194 ex->cone[ct][12*i+2] = 0; 195 ex->cone[ct][12*i+3] = i; 196 ex->cone[ct][12*i+4] = DM_POLYTOPE_SEGMENT; 197 ex->cone[ct][12*i+5] = 1; 198 ex->cone[ct][12*i+6] = 1; 199 ex->cone[ct][12*i+7] = i; 200 ex->cone[ct][12*i+8] = DM_POLYTOPE_SEGMENT; 201 ex->cone[ct][12*i+9] = 1; 202 ex->cone[ct][12*i+10] = 2; 203 ex->cone[ct][12*i+11] = i; 204 } 205 for (i = 0; i < 3*(Nl+1); ++i) ex->ornt[ct][i] = 0; 206 /* cones for triangular prisms/tensor triangular prisms */ 207 coff = 12*(Nl+1); 208 ooff = 3*(Nl+1); 209 for (i = 0; i < Nl; ++i) { 210 if (ex->useTensor) { 211 ex->cone[ct][coff+18*i+0] = DM_POLYTOPE_TRIANGLE; 212 ex->cone[ct][coff+18*i+1] = 0; 213 ex->cone[ct][coff+18*i+2] = i; 214 ex->cone[ct][coff+18*i+3] = DM_POLYTOPE_TRIANGLE; 215 ex->cone[ct][coff+18*i+4] = 0; 216 ex->cone[ct][coff+18*i+5] = i+1; 217 ex->cone[ct][coff+18*i+6] = DM_POLYTOPE_SEG_PRISM_TENSOR; 218 ex->cone[ct][coff+18*i+7] = 1; 219 ex->cone[ct][coff+18*i+8] = 0; 220 ex->cone[ct][coff+18*i+9] = i; 221 ex->cone[ct][coff+18*i+10] = DM_POLYTOPE_SEG_PRISM_TENSOR; 222 ex->cone[ct][coff+18*i+11] = 1; 223 ex->cone[ct][coff+18*i+12] = 1; 224 ex->cone[ct][coff+18*i+13] = i; 225 ex->cone[ct][coff+18*i+14] = DM_POLYTOPE_SEG_PRISM_TENSOR; 226 ex->cone[ct][coff+18*i+15] = 1; 227 ex->cone[ct][coff+18*i+16] = 2; 228 ex->cone[ct][coff+18*i+17] = i; 229 ex->ornt[ct][ooff+5*i+0] = 0; 230 ex->ornt[ct][ooff+5*i+1] = 0; 231 ex->ornt[ct][ooff+5*i+2] = 0; 232 ex->ornt[ct][ooff+5*i+3] = 0; 233 ex->ornt[ct][ooff+5*i+4] = 0; 234 } else { 235 ex->cone[ct][coff+18*i+0] = DM_POLYTOPE_TRIANGLE; 236 ex->cone[ct][coff+18*i+1] = 0; 237 ex->cone[ct][coff+18*i+2] = i; 238 ex->cone[ct][coff+18*i+3] = DM_POLYTOPE_TRIANGLE; 239 ex->cone[ct][coff+18*i+4] = 0; 240 ex->cone[ct][coff+18*i+5] = i+1; 241 ex->cone[ct][coff+18*i+6] = DM_POLYTOPE_QUADRILATERAL; 242 ex->cone[ct][coff+18*i+7] = 1; 243 ex->cone[ct][coff+18*i+8] = 0; 244 ex->cone[ct][coff+18*i+9] = i; 245 ex->cone[ct][coff+18*i+10] = DM_POLYTOPE_QUADRILATERAL; 246 ex->cone[ct][coff+18*i+11] = 1; 247 ex->cone[ct][coff+18*i+12] = 1; 248 ex->cone[ct][coff+18*i+13] = i; 249 ex->cone[ct][coff+18*i+14] = DM_POLYTOPE_QUADRILATERAL; 250 ex->cone[ct][coff+18*i+15] = 1; 251 ex->cone[ct][coff+18*i+16] = 2; 252 ex->cone[ct][coff+18*i+17] = i; 253 ex->ornt[ct][ooff+5*i+0] = -2; 254 ex->ornt[ct][ooff+5*i+1] = 0; 255 ex->ornt[ct][ooff+5*i+2] = 0; 256 ex->ornt[ct][ooff+5*i+3] = 0; 257 ex->ornt[ct][ooff+5*i+4] = 0; 258 } 259 } 260 /* DM_POLYTOPE_QUADRILATERAL produces Nl+1 quads and Nl hexes/tensor hexes */ 261 ct = DM_POLYTOPE_QUADRILATERAL; 262 ex->Nt[ct] = 2; 263 Nc = 16*(Nl+1) + 22*Nl; 264 No = 4*(Nl+1) + 6*Nl; 265 ierr = PetscMalloc4(ex->Nt[ct], &ex->target[ct], ex->Nt[ct], &ex->size[ct], Nc, &ex->cone[ct], No, &ex->ornt[ct]);CHKERRQ(ierr); 266 ex->target[ct][0] = DM_POLYTOPE_QUADRILATERAL; 267 ex->target[ct][1] = ex->useTensor ? DM_POLYTOPE_QUAD_PRISM_TENSOR : DM_POLYTOPE_HEXAHEDRON; 268 ex->size[ct][0] = Nl+1; 269 ex->size[ct][1] = Nl; 270 /* cones for quads */ 271 for (i = 0; i < Nl+1; ++i) { 272 ex->cone[ct][16*i+0] = DM_POLYTOPE_SEGMENT; 273 ex->cone[ct][16*i+1] = 1; 274 ex->cone[ct][16*i+2] = 0; 275 ex->cone[ct][16*i+3] = i; 276 ex->cone[ct][16*i+4] = DM_POLYTOPE_SEGMENT; 277 ex->cone[ct][16*i+5] = 1; 278 ex->cone[ct][16*i+6] = 1; 279 ex->cone[ct][16*i+7] = i; 280 ex->cone[ct][16*i+8] = DM_POLYTOPE_SEGMENT; 281 ex->cone[ct][16*i+9] = 1; 282 ex->cone[ct][16*i+10] = 2; 283 ex->cone[ct][16*i+11] = i; 284 ex->cone[ct][16*i+12] = DM_POLYTOPE_SEGMENT; 285 ex->cone[ct][16*i+13] = 1; 286 ex->cone[ct][16*i+14] = 3; 287 ex->cone[ct][16*i+15] = i; 288 } 289 for (i = 0; i < 4*(Nl+1); ++i) ex->ornt[ct][i] = 0; 290 /* cones for hexes/tensor hexes */ 291 coff = 16*(Nl+1); 292 ooff = 4*(Nl+1); 293 for (i = 0; i < Nl; ++i) { 294 if (ex->useTensor) { 295 ex->cone[ct][coff+22*i+0] = DM_POLYTOPE_QUADRILATERAL; 296 ex->cone[ct][coff+22*i+1] = 0; 297 ex->cone[ct][coff+22*i+2] = i; 298 ex->cone[ct][coff+22*i+3] = DM_POLYTOPE_QUADRILATERAL; 299 ex->cone[ct][coff+22*i+4] = 0; 300 ex->cone[ct][coff+22*i+5] = i+1; 301 ex->cone[ct][coff+22*i+6] = DM_POLYTOPE_SEG_PRISM_TENSOR; 302 ex->cone[ct][coff+22*i+7] = 1; 303 ex->cone[ct][coff+22*i+8] = 0; 304 ex->cone[ct][coff+22*i+9] = i; 305 ex->cone[ct][coff+22*i+10] = DM_POLYTOPE_SEG_PRISM_TENSOR; 306 ex->cone[ct][coff+22*i+11] = 1; 307 ex->cone[ct][coff+22*i+12] = 1; 308 ex->cone[ct][coff+22*i+13] = i; 309 ex->cone[ct][coff+22*i+14] = DM_POLYTOPE_SEG_PRISM_TENSOR; 310 ex->cone[ct][coff+22*i+15] = 1; 311 ex->cone[ct][coff+22*i+16] = 2; 312 ex->cone[ct][coff+22*i+17] = i; 313 ex->cone[ct][coff+22*i+18] = DM_POLYTOPE_SEG_PRISM_TENSOR; 314 ex->cone[ct][coff+22*i+19] = 1; 315 ex->cone[ct][coff+22*i+20] = 3; 316 ex->cone[ct][coff+22*i+21] = i; 317 ex->ornt[ct][ooff+6*i+0] = 0; 318 ex->ornt[ct][ooff+6*i+1] = 0; 319 ex->ornt[ct][ooff+6*i+2] = 0; 320 ex->ornt[ct][ooff+6*i+3] = 0; 321 ex->ornt[ct][ooff+6*i+4] = 0; 322 ex->ornt[ct][ooff+6*i+5] = 0; 323 } else { 324 ex->cone[ct][coff+22*i+0] = DM_POLYTOPE_QUADRILATERAL; 325 ex->cone[ct][coff+22*i+1] = 0; 326 ex->cone[ct][coff+22*i+2] = i; 327 ex->cone[ct][coff+22*i+3] = DM_POLYTOPE_QUADRILATERAL; 328 ex->cone[ct][coff+22*i+4] = 0; 329 ex->cone[ct][coff+22*i+5] = i+1; 330 ex->cone[ct][coff+22*i+6] = DM_POLYTOPE_QUADRILATERAL; 331 ex->cone[ct][coff+22*i+7] = 1; 332 ex->cone[ct][coff+22*i+8] = 0; 333 ex->cone[ct][coff+22*i+9] = i; 334 ex->cone[ct][coff+22*i+10] = DM_POLYTOPE_QUADRILATERAL; 335 ex->cone[ct][coff+22*i+11] = 1; 336 ex->cone[ct][coff+22*i+12] = 2; 337 ex->cone[ct][coff+22*i+13] = i; 338 ex->cone[ct][coff+22*i+14] = DM_POLYTOPE_QUADRILATERAL; 339 ex->cone[ct][coff+22*i+15] = 1; 340 ex->cone[ct][coff+22*i+16] = 1; 341 ex->cone[ct][coff+22*i+17] = i; 342 ex->cone[ct][coff+22*i+18] = DM_POLYTOPE_QUADRILATERAL; 343 ex->cone[ct][coff+22*i+19] = 1; 344 ex->cone[ct][coff+22*i+20] = 3; 345 ex->cone[ct][coff+22*i+21] = i; 346 ex->ornt[ct][ooff+6*i+0] = -2; 347 ex->ornt[ct][ooff+6*i+1] = 0; 348 ex->ornt[ct][ooff+6*i+2] = 0; 349 ex->ornt[ct][ooff+6*i+3] = 0; 350 ex->ornt[ct][ooff+6*i+4] = 0; 351 ex->ornt[ct][ooff+6*i+5] = 1; 352 } 353 } 354 /* Layers positions */ 355 if (!ex->Nth) { 356 if (ex->symmetric) for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness*l)/ex->layers - ex->thickness/2; 357 else for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] = (ex->thickness*l)/ex->layers; 358 } else { 359 ex->thickness = 0.; 360 ex->layerPos[0] = 0.; 361 for (l = 0; l < ex->layers; ++l) { 362 const PetscReal t = ex->thicknesses[PetscMin(l, ex->Nth-1)]; 363 ex->layerPos[l+1] = ex->layerPos[l] + t; 364 ex->thickness += t; 365 } 366 if (ex->symmetric) for (l = 0; l <= ex->layers; ++l) ex->layerPos[l] -= ex->thickness/2.; 367 } 368 PetscFunctionReturn(0); 369 } 370 371 static PetscErrorCode DMPlexTransformDestroy_Extrude(DMPlexTransform tr) 372 { 373 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 374 PetscInt ct; 375 PetscErrorCode ierr; 376 377 PetscFunctionBegin; 378 for (ct = 0; ct < DM_NUM_POLYTOPES; ++ct) { 379 ierr = PetscFree4(ex->target[ct], ex->size[ct], ex->cone[ct], ex->ornt[ct]);CHKERRQ(ierr); 380 } 381 ierr = PetscFree5(ex->Nt, ex->target, ex->size, ex->cone, ex->ornt);CHKERRQ(ierr); 382 ierr = PetscFree(ex->layerPos);CHKERRQ(ierr); 383 ierr = PetscFree(ex);CHKERRQ(ierr); 384 PetscFunctionReturn(0); 385 } 386 387 static PetscErrorCode DMPlexTransformGetSubcellOrientation_Extrude(DMPlexTransform tr, DMPolytopeType sct, PetscInt sp, PetscInt so, DMPolytopeType tct, PetscInt r, PetscInt o, PetscInt *rnew, PetscInt *onew) 388 { 389 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 390 391 PetscFunctionBeginHot; 392 *rnew = r; 393 *onew = DMPolytopeTypeComposeOrientation(tct, o, so); 394 if (!so) PetscFunctionReturn(0); 395 if (ex->useTensor) { 396 switch (sct) { 397 case DM_POLYTOPE_POINT: break; 398 case DM_POLYTOPE_SEGMENT: 399 switch (tct) { 400 case DM_POLYTOPE_SEGMENT: break; 401 case DM_POLYTOPE_SEG_PRISM_TENSOR: 402 *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -1 : 0); 403 break; 404 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 405 } 406 break; 407 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]); 408 } 409 } else { 410 switch (sct) { 411 case DM_POLYTOPE_POINT: break; 412 case DM_POLYTOPE_SEGMENT: 413 switch (tct) { 414 case DM_POLYTOPE_SEGMENT: break; 415 case DM_POLYTOPE_QUADRILATERAL: 416 *onew = DMPolytopeTypeComposeOrientation(tct, o, so ? -3 : 0); 417 break; 418 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type %s is not produced by %s", DMPolytopeTypes[tct], DMPolytopeTypes[sct]); 419 } 420 break; 421 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell type %s", DMPolytopeTypes[sct]); 422 } 423 } 424 PetscFunctionReturn(0); 425 } 426 427 static PetscErrorCode DMPlexTransformCellTransform_Extrude(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 428 { 429 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 430 431 PetscFunctionBegin; 432 if (rt) *rt = 0; 433 if (ex->Nt[source] < 0) { 434 /* Ignore cells that cannot be extruded */ 435 *Nt = 0; 436 *target = NULL; 437 *size = NULL; 438 *cone = NULL; 439 *ornt = NULL; 440 } else { 441 *Nt = ex->Nt[source]; 442 *target = ex->target[source]; 443 *size = ex->size[source]; 444 *cone = ex->cone[source]; 445 *ornt = ex->ornt[source]; 446 } 447 PetscFunctionReturn(0); 448 } 449 450 /* Computes new vertex along normal */ 451 static PetscErrorCode DMPlexTransformMapCoordinates_Extrude(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) 452 { 453 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 454 DM dm; 455 PetscReal ones2[2] = {0., 1.}, ones3[3] = { 0., 0., 1.}; 456 PetscReal normal[3] = {0., 0., 0.}, norm; 457 PetscBool computeNormal; 458 PetscInt dim, dEx = ex->cdimEx, cStart, cEnd, d; 459 PetscErrorCode ierr; 460 461 PetscFunctionBeginHot; 462 PetscCheckFalse(pct != DM_POLYTOPE_POINT,PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for parent point type %s",DMPolytopeTypes[pct]); 463 PetscCheckFalse(ct != DM_POLYTOPE_POINT,PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for refined point type %s",DMPolytopeTypes[ct]); 464 PetscCheckFalse(Nv != 1,PETSC_COMM_SELF,PETSC_ERR_SUP,"Vertices should be produced from a single vertex, not %D",Nv); 465 PetscCheckFalse(dE != ex->cdim,PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Coordinate dim %D != %D original dimension", dE, ex->cdim); 466 467 ierr = DMPlexTransformGetDM(tr, &dm);CHKERRQ(ierr); 468 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 469 computeNormal = dim != ex->cdim && !ex->useNormal ? PETSC_TRUE : PETSC_FALSE; 470 if (computeNormal) { 471 PetscInt *closure = NULL; 472 PetscInt closureSize, cl; 473 474 ierr = DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 475 ierr = DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure);CHKERRQ(ierr); 476 for (cl = 0; cl < closureSize*2; cl += 2) { 477 if ((closure[cl] >= cStart) && (closure[cl] < cEnd)) { 478 PetscReal cnormal[3] = {0, 0, 0}; 479 480 ierr = DMPlexComputeCellGeometryFVM(dm, closure[cl], NULL, NULL, cnormal);CHKERRQ(ierr); 481 for (d = 0; d < dEx; ++d) normal[d] += cnormal[d]; 482 } 483 } 484 ierr = DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure);CHKERRQ(ierr); 485 } else if (ex->useNormal) { 486 for (d = 0; d < dEx; ++d) normal[d] = ex->normal[d]; 487 } else if (ex->cdimEx == 2) { 488 for (d = 0; d < dEx; ++d) normal[d] = ones2[d]; 489 } else if (ex->cdimEx == 3) { 490 for (d = 0; d < dEx; ++d) normal[d] = ones3[d]; 491 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to determine normal for extrusion"); 492 493 for (d = 0, norm = 0.0; d < dEx; ++d) norm += PetscSqr(normal[d]); 494 for (d = 0; d < dEx; ++d) normal[d] *= 1./PetscSqrtReal(norm); 495 for (d = 0; d < dEx; ++d) out[d] = normal[d]*ex->layerPos[r]; 496 for (d = 0; d < ex->cdim; ++d) out[d] += in[d]; 497 PetscFunctionReturn(0); 498 } 499 500 static PetscErrorCode DMPlexTransformInitialize_Extrude(DMPlexTransform tr) 501 { 502 PetscFunctionBegin; 503 tr->ops->view = DMPlexTransformView_Extrude; 504 tr->ops->setfromoptions = DMPlexTransformSetFromOptions_Extrude; 505 tr->ops->setup = DMPlexTransformSetUp_Extrude; 506 tr->ops->destroy = DMPlexTransformDestroy_Extrude; 507 tr->ops->setdimensions = DMPlexTransformSetDimensions_Extrude; 508 tr->ops->celltransform = DMPlexTransformCellTransform_Extrude; 509 tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientation_Extrude; 510 tr->ops->mapcoordinates = DMPlexTransformMapCoordinates_Extrude; 511 PetscFunctionReturn(0); 512 } 513 514 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform tr) 515 { 516 DMPlexTransform_Extrude *ex; 517 DM dm; 518 PetscInt dim; 519 PetscErrorCode ierr; 520 521 PetscFunctionBegin; 522 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 523 ierr = PetscNewLog(tr, &ex);CHKERRQ(ierr); 524 tr->data = ex; 525 ex->thickness = 1.; 526 ex->useTensor = PETSC_TRUE; 527 ex->symmetric = PETSC_FALSE; 528 ex->useNormal = PETSC_FALSE; 529 ex->layerPos = NULL; 530 ierr = DMPlexTransformGetDM(tr, &dm);CHKERRQ(ierr); 531 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 532 ierr = DMGetCoordinateDim(dm, &ex->cdim);CHKERRQ(ierr); 533 ex->cdimEx = ex->cdim == dim ? ex->cdim+1 : ex->cdim; 534 ierr = DMPlexTransformExtrudeSetLayers(tr, 1);CHKERRQ(ierr); 535 ierr = DMPlexTransformInitialize_Extrude(tr);CHKERRQ(ierr); 536 PetscFunctionReturn(0); 537 } 538 539 /*@ 540 DMPlexTransformExtrudeGetLayers - Get the number of extruded layers. 541 542 Not collective 543 544 Input Parameter: 545 . tr - The DMPlexTransform 546 547 Output Parameter: 548 . layers - The number of layers 549 550 Level: intermediate 551 552 .seealso: DMPlexTransformExtrudeSetLayers() 553 @*/ 554 PetscErrorCode DMPlexTransformExtrudeGetLayers(DMPlexTransform tr, PetscInt *layers) 555 { 556 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 557 558 PetscFunctionBegin; 559 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 560 PetscValidIntPointer(layers, 2); 561 *layers = ex->layers; 562 PetscFunctionReturn(0); 563 } 564 565 /*@ 566 DMPlexTransformExtrudeSetLayers - Set the number of extruded layers. 567 568 Not collective 569 570 Input Parameters: 571 + tr - The DMPlexTransform 572 - layers - The number of layers 573 574 Level: intermediate 575 576 .seealso: DMPlexTransformExtrudeGetLayers() 577 @*/ 578 PetscErrorCode DMPlexTransformExtrudeSetLayers(DMPlexTransform tr, PetscInt layers) 579 { 580 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 581 PetscErrorCode ierr; 582 583 PetscFunctionBegin; 584 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 585 ex->layers = layers; 586 ierr = PetscFree(ex->layerPos);CHKERRQ(ierr); 587 ierr = PetscCalloc1(ex->layers+1, &ex->layerPos);CHKERRQ(ierr); 588 PetscFunctionReturn(0); 589 } 590 591 /*@ 592 DMPlexTransformExtrudeGetThickness - Get the total thickness of the layers 593 594 Not collective 595 596 Input Parameter: 597 . tr - The DMPlexTransform 598 599 Output Parameter: 600 . thickness - The total thickness of the layers 601 602 Level: intermediate 603 604 .seealso: DMPlexTransformExtrudeSetThickness() 605 @*/ 606 PetscErrorCode DMPlexTransformExtrudeGetThickness(DMPlexTransform tr, PetscReal *thickness) 607 { 608 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 609 610 PetscFunctionBegin; 611 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 612 PetscValidRealPointer(thickness, 2); 613 *thickness = ex->thickness; 614 PetscFunctionReturn(0); 615 } 616 617 /*@ 618 DMPlexTransformExtrudeSetThickness - Set the total thickness of the layers 619 620 Not collective 621 622 Input Parameters: 623 + tr - The DMPlexTransform 624 - thickness - The total thickness of the layers 625 626 Level: intermediate 627 628 .seealso: DMPlexTransformExtrudeGetThickness() 629 @*/ 630 PetscErrorCode DMPlexTransformExtrudeSetThickness(DMPlexTransform tr, PetscReal thickness) 631 { 632 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 633 634 PetscFunctionBegin; 635 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 636 PetscCheckFalse(thickness <= 0.,PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_OUTOFRANGE, "Height of layers %g must be positive", (double) thickness); 637 ex->thickness = thickness; 638 PetscFunctionReturn(0); 639 } 640 641 PetscErrorCode DMPlexTransformExtrudeGetTensor(DMPlexTransform tr, PetscBool *useTensor) 642 { 643 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 644 645 PetscFunctionBegin; 646 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 647 PetscValidBoolPointer(useTensor, 2); 648 *useTensor = ex->useTensor; 649 PetscFunctionReturn(0); 650 } 651 652 PetscErrorCode DMPlexTransformExtrudeSetTensor(DMPlexTransform tr, PetscBool useTensor) 653 { 654 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 655 656 PetscFunctionBegin; 657 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 658 ex->useTensor = useTensor; 659 PetscFunctionReturn(0); 660 } 661 662 PetscErrorCode DMPlexTransformExtrudeGetSymmetric(DMPlexTransform tr, PetscBool *symmetric) 663 { 664 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 665 666 PetscFunctionBegin; 667 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 668 PetscValidBoolPointer(symmetric, 2); 669 *symmetric = ex->symmetric; 670 PetscFunctionReturn(0); 671 } 672 673 PetscErrorCode DMPlexTransformExtrudeSetSymmetric(DMPlexTransform tr, PetscBool symmetric) 674 { 675 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 676 677 PetscFunctionBegin; 678 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 679 ex->symmetric = symmetric; 680 PetscFunctionReturn(0); 681 } 682 683 PetscErrorCode DMPlexTransformExtrudeSetNormal(DMPlexTransform tr, const PetscReal normal[]) 684 { 685 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 686 PetscInt d; 687 688 PetscFunctionBegin; 689 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 690 ex->useNormal = PETSC_TRUE; 691 for (d = 0; d < ex->cdimEx; ++d) ex->normal[d] = normal[d]; 692 PetscFunctionReturn(0); 693 } 694 695 PetscErrorCode DMPlexTransformExtrudeSetThicknesses(DMPlexTransform tr, PetscInt Nth, const PetscReal thicknesses[]) 696 { 697 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 698 PetscInt t; 699 PetscErrorCode ierr; 700 701 PetscFunctionBegin; 702 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 703 PetscCheckFalse(Nth <= 0,PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_OUTOFRANGE, "Number of thicknesses %D must be positive", Nth); 704 ex->Nth = PetscMin(Nth, ex->layers); 705 ierr = PetscFree(ex->thicknesses);CHKERRQ(ierr); 706 ierr = PetscMalloc1(ex->Nth, &ex->thicknesses);CHKERRQ(ierr); 707 for (t = 0; t < ex->Nth; ++t) { 708 PetscCheckFalse(thicknesses[t] <= 0.,PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_OUTOFRANGE, "Thickness %g of layer %D must be positive", (double) thicknesses[t], t); 709 ex->thicknesses[t] = thicknesses[t]; 710 } 711 PetscFunctionReturn(0); 712 } 713