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 SETERRQ1(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 if (Nc != ex->cdimEx) SETERRQ2(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 if (!nl) SETERRQ(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] = 1; 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] = 2; 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 DMPlexTransformCellTransform_Extrude(DMPlexTransform tr, DMPolytopeType source, PetscInt p, PetscInt *rt, PetscInt *Nt, DMPolytopeType *target[], PetscInt *size[], PetscInt *cone[], PetscInt *ornt[]) 388 { 389 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 390 391 PetscFunctionBegin; 392 if (rt) *rt = 0; 393 if (ex->Nt[source] < 0) { 394 /* Ignore cells that cannot be extruded */ 395 *Nt = 0; 396 *target = NULL; 397 *size = NULL; 398 *cone = NULL; 399 *ornt = NULL; 400 } else { 401 *Nt = ex->Nt[source]; 402 *target = ex->target[source]; 403 *size = ex->size[source]; 404 *cone = ex->cone[source]; 405 *ornt = ex->ornt[source]; 406 } 407 PetscFunctionReturn(0); 408 } 409 410 /* Computes new vertex along normal */ 411 static PetscErrorCode DMPlexTransformMapCoordinates_Extrude(DMPlexTransform tr, DMPolytopeType pct, DMPolytopeType ct, PetscInt p, PetscInt r, PetscInt Nv, PetscInt dE, const PetscScalar in[], PetscScalar out[]) 412 { 413 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 414 DM dm; 415 PetscReal ones2[2] = {0., 1.}, ones3[3] = { 0., 0., 1.}; 416 PetscReal normal[3] = {0., 0., 0.}, norm; 417 PetscBool computeNormal; 418 PetscInt dim, dEx = ex->cdimEx, cStart, cEnd, d; 419 PetscErrorCode ierr; 420 421 PetscFunctionBeginHot; 422 if (pct != DM_POLYTOPE_POINT) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for parent point type %s",DMPolytopeTypes[pct]); 423 if (ct != DM_POLYTOPE_POINT) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for refined point type %s",DMPolytopeTypes[ct]); 424 if (Nv != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vertices should be produced from a single vertex, not %D",Nv); 425 if (dE != ex->cdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Coordinate dim %D != %D original dimension", dE, ex->cdim); 426 427 ierr = DMPlexTransformGetDM(tr, &dm);CHKERRQ(ierr); 428 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 429 computeNormal = dim != ex->cdim && !ex->useNormal ? PETSC_TRUE : PETSC_FALSE; 430 if (computeNormal) { 431 PetscInt *closure = NULL; 432 PetscInt closureSize, cl; 433 434 ierr = DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 435 ierr = DMPlexGetTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure);CHKERRQ(ierr); 436 for (cl = 0; cl < closureSize*2; cl += 2) { 437 if ((closure[cl] >= cStart) && (closure[cl] < cEnd)) { 438 PetscReal cnormal[3] = {0, 0, 0}; 439 440 ierr = DMPlexComputeCellGeometryFVM(dm, closure[cl], NULL, NULL, cnormal);CHKERRQ(ierr); 441 for (d = 0; d < dEx; ++d) normal[d] += cnormal[d]; 442 } 443 } 444 ierr = DMPlexRestoreTransitiveClosure(dm, p, PETSC_FALSE, &closureSize, &closure);CHKERRQ(ierr); 445 } else if (ex->useNormal) { 446 for (d = 0; d < dEx; ++d) normal[d] = ex->normal[d]; 447 } else if (ex->cdimEx == 2) { 448 for (d = 0; d < dEx; ++d) normal[d] = ones2[d]; 449 } else if (ex->cdimEx == 3) { 450 for (d = 0; d < dEx; ++d) normal[d] = ones3[d]; 451 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to determine normal for extrusion"); 452 453 for (d = 0, norm = 0.0; d < dEx; ++d) norm += PetscSqr(normal[d]); 454 for (d = 0; d < dEx; ++d) normal[d] *= 1./PetscSqrtReal(norm); 455 for (d = 0; d < dEx; ++d) out[d] = normal[d]*ex->layerPos[r]; 456 for (d = 0; d < ex->cdim; ++d) out[d] += in[d]; 457 PetscFunctionReturn(0); 458 } 459 460 static PetscErrorCode DMPlexTransformInitialize_Extrude(DMPlexTransform tr) 461 { 462 PetscFunctionBegin; 463 tr->ops->view = DMPlexTransformView_Extrude; 464 tr->ops->setfromoptions = DMPlexTransformSetFromOptions_Extrude; 465 tr->ops->setup = DMPlexTransformSetUp_Extrude; 466 tr->ops->destroy = DMPlexTransformDestroy_Extrude; 467 tr->ops->setdimensions = DMPlexTransformSetDimensions_Extrude; 468 tr->ops->celltransform = DMPlexTransformCellTransform_Extrude; 469 tr->ops->getsubcellorientation = DMPlexTransformGetSubcellOrientationIdentity; 470 tr->ops->mapcoordinates = DMPlexTransformMapCoordinates_Extrude; 471 PetscFunctionReturn(0); 472 } 473 474 PETSC_EXTERN PetscErrorCode DMPlexTransformCreate_Extrude(DMPlexTransform tr) 475 { 476 DMPlexTransform_Extrude *ex; 477 DM dm; 478 PetscInt dim; 479 PetscErrorCode ierr; 480 481 PetscFunctionBegin; 482 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 483 ierr = PetscNewLog(tr, &ex);CHKERRQ(ierr); 484 tr->data = ex; 485 ex->thickness = 1.; 486 ex->useTensor = PETSC_TRUE; 487 ex->symmetric = PETSC_FALSE; 488 ex->useNormal = PETSC_FALSE; 489 ex->layerPos = NULL; 490 ierr = DMPlexTransformGetDM(tr, &dm);CHKERRQ(ierr); 491 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 492 ierr = DMGetCoordinateDim(dm, &ex->cdim);CHKERRQ(ierr); 493 ex->cdimEx = ex->cdim == dim ? ex->cdim+1 : ex->cdim; 494 ierr = DMPlexTransformExtrudeSetLayers(tr, 1);CHKERRQ(ierr); 495 ierr = DMPlexTransformInitialize_Extrude(tr);CHKERRQ(ierr); 496 PetscFunctionReturn(0); 497 } 498 499 /*@ 500 DMPlexTransformExtrudeGetLayers - Get the number of extruded layers. 501 502 Not collective 503 504 Input Parameter: 505 . tr - The DMPlexTransform 506 507 Output Parameter: 508 . layers - The number of layers 509 510 Level: intermediate 511 512 .seealso: DMPlexTransformExtrudeSetLayers() 513 @*/ 514 PetscErrorCode DMPlexTransformExtrudeGetLayers(DMPlexTransform tr, PetscInt *layers) 515 { 516 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 517 518 PetscFunctionBegin; 519 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 520 PetscValidIntPointer(layers, 2); 521 *layers = ex->layers; 522 PetscFunctionReturn(0); 523 } 524 525 /*@ 526 DMPlexTransformExtrudeSetLayers - Set the number of extruded layers. 527 528 Not collective 529 530 Input Parameters: 531 + tr - The DMPlexTransform 532 - layers - The number of layers 533 534 Level: intermediate 535 536 .seealso: DMPlexTransformExtrudeGetLayers() 537 @*/ 538 PetscErrorCode DMPlexTransformExtrudeSetLayers(DMPlexTransform tr, PetscInt layers) 539 { 540 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 541 PetscErrorCode ierr; 542 543 PetscFunctionBegin; 544 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 545 ex->layers = layers; 546 ierr = PetscFree(ex->layerPos);CHKERRQ(ierr); 547 ierr = PetscCalloc1(ex->layers+1, &ex->layerPos);CHKERRQ(ierr); 548 PetscFunctionReturn(0); 549 } 550 551 /*@ 552 DMPlexTransformExtrudeGetThickness - Get the total thickness of the layers 553 554 Not collective 555 556 Input Parameter: 557 . tr - The DMPlexTransform 558 559 Output Parameter: 560 . thickness - The total thickness of the layers 561 562 Level: intermediate 563 564 .seealso: DMPlexTransformExtrudeSetThickness() 565 @*/ 566 PetscErrorCode DMPlexTransformExtrudeGetThickness(DMPlexTransform tr, PetscReal *thickness) 567 { 568 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 569 570 PetscFunctionBegin; 571 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 572 PetscValidRealPointer(thickness, 2); 573 *thickness = ex->thickness; 574 PetscFunctionReturn(0); 575 } 576 577 /*@ 578 DMPlexTransformExtrudeSetThickness - Set the total thickness of the layers 579 580 Not collective 581 582 Input Parameters: 583 + tr - The DMPlexTransform 584 - thickness - The total thickness of the layers 585 586 Level: intermediate 587 588 .seealso: DMPlexTransformExtrudeGetThickness() 589 @*/ 590 PetscErrorCode DMPlexTransformExtrudeSetThickness(DMPlexTransform tr, PetscReal thickness) 591 { 592 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 593 594 PetscFunctionBegin; 595 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 596 if (thickness <= 0.) SETERRQ1(PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_OUTOFRANGE, "Height of layers %g must be positive", (double) thickness); 597 ex->thickness = thickness; 598 PetscFunctionReturn(0); 599 } 600 601 PetscErrorCode DMPlexTransformExtrudeGetTensor(DMPlexTransform tr, PetscBool *useTensor) 602 { 603 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 604 605 PetscFunctionBegin; 606 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 607 PetscValidBoolPointer(useTensor, 2); 608 *useTensor = ex->useTensor; 609 PetscFunctionReturn(0); 610 } 611 612 PetscErrorCode DMPlexTransformExtrudeSetTensor(DMPlexTransform tr, PetscBool useTensor) 613 { 614 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 615 616 PetscFunctionBegin; 617 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 618 ex->useTensor = useTensor; 619 PetscFunctionReturn(0); 620 } 621 622 PetscErrorCode DMPlexTransformExtrudeGetSymmetric(DMPlexTransform tr, PetscBool *symmetric) 623 { 624 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 625 626 PetscFunctionBegin; 627 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 628 PetscValidBoolPointer(symmetric, 2); 629 *symmetric = ex->symmetric; 630 PetscFunctionReturn(0); 631 } 632 633 PetscErrorCode DMPlexTransformExtrudeSetSymmetric(DMPlexTransform tr, PetscBool symmetric) 634 { 635 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 636 637 PetscFunctionBegin; 638 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 639 ex->symmetric = symmetric; 640 PetscFunctionReturn(0); 641 } 642 643 PetscErrorCode DMPlexTransformExtrudeSetNormal(DMPlexTransform tr, const PetscReal normal[]) 644 { 645 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 646 PetscInt d; 647 648 PetscFunctionBegin; 649 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 650 ex->useNormal = PETSC_TRUE; 651 for (d = 0; d < ex->cdimEx; ++d) ex->normal[d] = normal[d]; 652 PetscFunctionReturn(0); 653 } 654 655 PetscErrorCode DMPlexTransformExtrudeSetThicknesses(DMPlexTransform tr, PetscInt Nth, const PetscReal thicknesses[]) 656 { 657 DMPlexTransform_Extrude *ex = (DMPlexTransform_Extrude *) tr->data; 658 PetscInt t; 659 PetscErrorCode ierr; 660 661 PetscFunctionBegin; 662 PetscValidHeaderSpecific(tr, DMPLEXTRANSFORM_CLASSID, 1); 663 if (Nth <= 0) SETERRQ1(PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_OUTOFRANGE, "Number of thicknesses %D must be positive", Nth); 664 ex->Nth = PetscMin(Nth, ex->layers); 665 ierr = PetscFree(ex->thicknesses);CHKERRQ(ierr); 666 ierr = PetscMalloc1(ex->Nth, &ex->thicknesses);CHKERRQ(ierr); 667 for (t = 0; t < ex->Nth; ++t) { 668 if (thicknesses[t] <= 0.) SETERRQ2(PetscObjectComm((PetscObject) tr), PETSC_ERR_ARG_OUTOFRANGE, "Thickness %g of layer %D must be positive", (double) thicknesses[t], t); 669 ex->thicknesses[t] = thicknesses[t]; 670 } 671 PetscFunctionReturn(0); 672 } 673