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