Faultsystems2d. not finished yet.

 1 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % 4 % Copyright (c) 2003-2010 by University of Queensland 5 % Earth Systems Science Computational Center (ESSCC) 6 7 % 8 % Primary Business: Queensland, Australia 9 % Licensed under the Open Software License version 3.0 10 11 % 12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 13 14 \chapter{The \pycad Module}\label{PYCAD CHAP} 15 16 \section{Introduction} 17 18 \pycad provides a simple way to build a mesh for your finite element 19 simulation. You begin by building what we call a \class{Design} using 20 primitive geometric objects, and then go on to build a mesh from this. 21 The final step of generating the mesh from a \class{Design} uses freely 22 available mesh generation software, such as \gmshextern. 23 24 A \class{Design} is built by defining points, which are used to specify 25 the corners of geometric objects and the vertices of curves. Using 26 points you construct more interesting objects such as lines, 27 rectangles, and arcs. By adding many of these objects into a \class{Design}, 28 you can build meshes for arbitrarily complex 2-D and 3-D structures. 29 30 \section{The Unit Square} 31 The simplest geometry is the unit square. First we generate the corner points 32 \begin{python} 33 from esys.pycad import * 34 p0=Point(0.,0.,0.) 35 p1=Point(1.,0.,0.) 36 p2=Point(1.,1.,0.) 37 p3=Point(0.,1.,0.) 38 \end{python} 39 which are then linked to define the edges of the square 40 \begin{python} 41 l01=Line(p0,p1) 42 l12=Line(p1,p2) 43 l23=Line(p2,p3) 44 l30=Line(p3,p0) 45 \end{python} 46 The lines are put together to form a loop 47 \begin{python} 48 c=CurveLoop(l01,l12,l23,l30) 49 \end{python} 50 The orientation of the line defining the \class{CurveLoop} is important. 51 It is assumed that the surrounded area is to the left when moving along the 52 lines from their starting points towards the end points. 53 Moreover, the line needs to form a closed loop. 54 We now use the \class{CurveLoop} to define a surface 55 \begin{python} 56 s=PlaneSurface(c) 57 \end{python} 58 Note that there is a difference between the \class{CurveLoop}, which defines 59 the boundary of the surface, and the actual surface \class{PlaneSurface}. 60 This difference becomes clearer in the next example with a hole. 61 Now we are ready to define the geometry which is described by an instance of 62 the \class{Design} class: 63 \begin{python} 64 d=Design(dim=2,element_size=0.05) 65 \end{python} 66 Here we use the two-dimensional domain with a local element size in the finite 67 element mesh of $0.05$. 68 We then add the surface \code{s} to the geometry 69 \begin{python} 70 d.addItems(s) 71 \end{python} 72 This will automatically import all items used to construct \code{s} into the 73 \class{Design} \code{d}. 74 Now we are ready to construct a \finley FEM mesh and then write it to the file 75 \file{quad.fly}: 76 \begin{python} 77 from esys.finley import MakeDomain 78 dom=MakeDomain(d) 79 dom.write("quad.fly") 80 \end{python} 81 In some cases it is useful to access the script used to generate the geometry. 82 You can specify a specific name for the script file. In our case we use 83 \begin{python} 84 d.setScriptFileName("quad.geo") 85 \end{python} 86 It is also useful to check error messages generated during the mesh generation 87 process. \gmshextern writes messages to the file \file{.gmsh-errors} in your 88 home directory. 89 Putting everything together we get the script 90 \begin{python} 91 from esys.pycad import * 92 from esys.pycad.gmsh import Design 93 from esys.finley import MakeDomain 94 p0=Point(0.,0.,0.) 95 p1=Point(1.,0.,0.) 96 p2=Point(1.,1.,0.) 97 p3=Point(0.,1.,0.) 98 l01=Line(p0,p1) 99 l12=Line(p1,p2) 100 l23=Line(p2,p3) 101 l30=Line(p3,p0) 102 c=CurveLoop(l01,l12,l23,l30) 103 s=PlaneSurface(c) 104 d=Design(dim=2,element_size=0.05) 105 d.setScriptFileName("quad.geo") 106 d.addItems(s) 107 pl1=PropertySet("sides",l01,l23) 108 pl2=PropertySet("top_and_bottom",l12,l30) 109 d.addItems(pl1, pl2) 110 dom=MakeDomain(d) 111 dom.write("quad.fly") 112 \end{python} 113 This example is included with the software in \file{quad.py} in the 114 \ExampleDirectory. 115 116 There are three extra statements which we have not discussed yet. 117 By default the mesh used to subdivide the boundary is not written into the 118 mesh file mainly to reduce the size of the data file. 119 One needs to explicitly add the lines to the \Design which should be present 120 in the mesh data. Here we additionally labeled the lines on the top and the 121 bottom with the name top_and_bottom and the lines on the left and right 122 hand side with the name sides using \class{PropertySet} objects. 123 The labeling is convenient when using tagging\index{tagging}, see \Chap{ESCRIPT CHAP}. 124 125 \begin{figure} 126 \centerline{\includegraphics{quad}} 127 \caption{Quadrilateral with mesh of triangles} 128 \label{fig:PYCAD 0} 129 \end{figure} 130 131 If you have \gmshextern installed you can run the example and view the 132 geometry and mesh with: 133 \begin{verbatim} 134 run-escript quad.py 135 gmsh quad.geo 136 gmsh quad.msh 137 \end{verbatim} 138 See Figure~\ref{fig:PYCAD 0} for a result. 139 140 In most cases it is best practice to generate the mesh and solve the 141 mathematical model in two separate scripts. In our example you can read the 142 \finley mesh into your simulation code\footnote{\gmshextern files can be 143 directly read using \function{ReadGmsh}, see \Chap{CHAPTER ON FINLEY}} using 144 \begin{python} 145 from finley import ReadMesh 146 mesh=ReadMesh("quad.fly") 147 \end{python} 148 Note that the underlying mesh generation software will not accept all 149 the geometries you can create. For example, \pycad will happily allow 150 you to create a 2-D \class{Design} that is a closed loop with some 151 additional points or lines lying outside of the enclosed area, but 152 \gmshextern will fail to create a mesh for it. 153 154 \section{Holes} 155 The example included below shows how to use \pycad to create a 2-D mesh 156 in the shape of a trapezoid with a cut-out area as in Figure~\ref{fig:PYCAD 1}. 157 158 \begin{python} 159 from esys.pycad import * 160 from esys.pycad.gmsh import Design 161 from esys.finley import MakeDomain 162 163 # A trapezoid 164 p0=Point(0.0, 0.0, 0.0) 165 p1=Point(1.0, 0.0, 0.0) 166 p2=Point(1.0, 0.5, 0.0) 167 p3=Point(0.0, 1.0, 0.0) 168 l01=Line(p0, p1) 169 l12=Line(p1, p2) 170 l23=Line(p2, p3) 171 l30=Line(p3, p0) 172 c=CurveLoop(l01, l12, l23, l30) 173 174 # A small triangular cutout 175 x0=Point(0.1, 0.1, 0.0) 176 x1=Point(0.5, 0.1, 0.0) 177 x2=Point(0.5, 0.2, 0.0) 178 x01=Line(x0, x1) 179 x12=Line(x1, x2) 180 x20=Line(x2, x0) 181 cutout=CurveLoop(x01, x12, x20) 182 183 # Create the surface with cutout 184 s=PlaneSurface(c, holes=[cutout]) 185 186 # Create a Design which can make the mesh 187 d=Design(dim=2, element_size=0.05) 188 189 # Add the trapezoid with cutout 190 d.addItems(s) 191 192 # Create the geometry, mesh and Escript domain 193 d.setScriptFileName("trapezoid.geo") 194 d.setMeshFileName("trapezoid.msh") 195 domain=MakeDomain(d) 196 # write mesh to a finley file: 197 domain.write("trapezoid.fly") 198 \end{python} 199 This example is included with the software in \file{trapezoid.py} in the 200 \ExampleDirectory. 201 202 \begin{figure} 203 \centerline{\includegraphics{trapezoid}} 204 \caption{Trapezoid with triangular hole} 205 \label{fig:PYCAD 1} 206 \end{figure} 207 208 A \code{CurveLoop} is used to connect several lines into a single curve. 209 It is used in the example above to create the trapezoidal outline for the grid 210 and also for the triangular cutout area. 211 You can define any number of lines when creating a \class{CurveLoop}, but 212 the end of one line must be identical to the start of the next. 213 214 \section{A 3D example} 215 216 In this section we discuss the definition of 3D geometries. As an example 217 consider the unit cube as shown in Figure~\ref{fig:PYCAD 2}. 218 First we generate the vertices of the cube: 219 \begin{python} 220 from esys.pycad import * 221 p0=Point(0.,0.,0.) 222 p1=Point(1.,0.,0.) 223 p2=Point(0.,1.,0.) 224 p3=Point(1.,1.,0.) 225 p4=Point(0.,0.,1.) 226 p5=Point(1.,0.,1.) 227 p6=Point(0.,1.,1.) 228 p7=Point(1.,1.,1.) 229 \end{python} 230 231 We connect the points to form the bottom and top surfaces of the cube: 232 \begin{python} 233 l01=Line(p0,p1) 234 l13=Line(p1,p3) 235 l32=Line(p3,p2) 236 l20=Line(p2,p0) 237 bottom=PlaneSurface(-CurveLoop(l01,l13,l32,l20)) 238 \end{python} 239 240 \begin{figure} 241 \centerline{\includegraphics{brick}} 242 \caption{Three dimensional block} 243 \label{fig:PYCAD 2} 244 \end{figure} 245 246 Similar to the definition of a \code{CurvedLoop} the orientation of the 247 surfaces in a \code{SurfaceLoop} is relevant. 248 In fact, the surface normal direction defined by the right-hand rule needs to 249 point outwards as indicated by the surface normals in Figure~\ref{fig:PYCAD 2}. 250 As the \code{bottom} face is directed upwards it is inserted with the minus 251 sign into the \code{SurfaceLoop} in order to adjust the orientation of the 252 surface. Similarly we set 253 \begin{python} 254 l45=Line(p4,p5) 255 l57=Line(p5,p7) 256 l76=Line(p7,p6) 257 l64=Line(p6,p4) 258 top=PlaneSurface(CurveLoop(l45,l57,l76,l64)) 259 \end{python} 260 To form the front face we introduce the two additional lines connecting the 261 left and right front points of the \code{top} and \code{bottom} face: 262 \begin{python} 263 l15=Line(p1,p5) 264 l40=Line(p4,p0) 265 \end{python} 266 To form the front face we encounter the problem as the line \code{l45} used 267 to define the \code{top} face is pointing the wrong direction. 268 In \pycad you can reversing direction of an object by changing its sign. 269 So we write \code{-l45} to indicate that the direction is to be reversed. 270 With this notation we can write 271 \begin{python} 272 front=PlaneSurface(CurveLoop(l01,l15,-l45,l40)) 273 \end{python} 274 Keep in mind that if you use \code{Line(p4,p5)} instead of \code{-l45} both 275 objects are treated as different although connecting the same points with a 276 straight line in the same direction. The resulting geometry would include an 277 opening along the \code{p4}--\code{p5} connection. 278 This will lead to an inconsistent mesh and may result in a failure of the 279 volumetric mesh generator. Similarly we can define the other sides of the cube: 280 \begin{python} 281 l37=Line(p3,p7) 282 l62=Line(p6,p2) 283 back=PlaneSurface(CurveLoop(l32,-l62,-l76,-l37)) 284 left=PlaneSurface(CurveLoop(-l40,-l64,l62,l20)) 285 right=PlaneSurface(CurveLoop(-l15,l13,l37,-l57)) 286 \end{python} 287 We can now put the six surfaces together to form a \class{SurfaceLoop} 288 defining the boundary of the volume of the cube: 289 \begin{python} 290 sl=SurfaceLoop(top,bottom,front,back,left,right) 291 v=Volume(sl) 292 \end{python} 293 As in the 2D case, the \class{Design} class is used to define the geometry: 294 \begin{python} 295 from esys.pycad.gmsh import Design 296 from esys.finley import MakeDomain 297 298 des=Design(dim=3, element_size = 0.1, keep_files=True) 299 des.setScriptFileName("brick.geo") 300 des.addItems(v, top, bottom, back, front, left, right) 301 302 dom=MakeDomain(des) 303 dom.write("brick.fly") 304 \end{python} 305 Note that the \finley mesh file \file{brick.fly} will contain the 306 triangles used to define the surfaces as they are added to the \class{Design}. 307 The example script of the cube is included with the software in 308 \file{brick.py} in the \ExampleDirectory. 309 310 \section{Alternative File Formats} 311 \pycad supports other file formats including I-DEAS universal file, VRML, 312 Nastran and STL. The following example shows how to generate the STL file 313 \file{brick.stl}: 314 \begin{python} 315 from esys.pycad.gmsh import Design 316 317 des=Design(dim=3, element_size = 0.1, keep_files=True) 318 des.addItems(v, top, bottom, back, front, left , right) 319 320 des.setFileFormat(des.STL) 321 des.setMeshFileName("brick.stl") 322 des.generate() 323 \end{python} 324 The example script of the cube is included with the software in 325 \file{brick_stl.py} in the \ExampleDirectory. 326 327 \section{Element Sizes} 328 The element size used globally is defined by the \code{element_size} argument 329 of the \class{Design}. The mesh generator will try to use this mesh size 330 everywhere in the geometry. In some cases it can be desirable to use a finer 331 mesh locally. A local refinement can be defined at each \class{Point}: 332 \begin{python} 333 p0=Point(0., 0., 0., local_scale=0.01) 334 \end{python} 335 Here the mesh generator will create a mesh with an element size which is by 336 the factor \code{0.01} times smaller than the global mesh size 337 \code{element_size=0.3}, see Figure~\ref{fig:PYCAD 5}. 338 The point where a refinement is defined must be a point on a curve used to 339 define the geometry. 340 341 \begin{figure} 342 \centerline{\includegraphics{refine}} 343 \caption{Local refinement at the origin by \var{local_scale=0.01} 344 with \var{element_size=0.3} and number of elements on the top set to 10} 345 \label{fig:PYCAD 5} 346 \end{figure} 347 348 Alternatively, one can define a mesh size along a curve by defining the number 349 of elements to be used to subdivide the curve. For instance, to use $20$ 350 elements on line \code{l23}: 351 \begin{python} 352 l23=Line(p2, p3) 353 l23.setElementDistribution(20) 354 \end{python} 355 Setting the number of elements on a curve overwrites the global mesh size 356 \code{element_size}. The result is shown in Figure~\ref{fig:PYCAD 5}. 357 358 \section{\pycad Classes} 359 %\declaremodule{extension}{esys.pycad} 360 %\modulesynopsis{Python geometry description and meshing interface} 361 362 \subsection{Primitives} 363 Some of the most commonly-used objects in \pycad are listed here. 364 For a more complete list see the full API documentation. 365 366 \begin{classdesc}{Point}{x=0.,y=0.,z=0.\optional{,local_scale=1.}} 367 creates a point at the given coordinates with local characteristic length \var{local_scale} 368 \end{classdesc} 369 370 \begin{classdesc}{CurveLoop}{list} 371 creates a closed curve from a \code{list} of 372 \class{Line}, \class{Arc}, \class{Spline}, \class{BSpline}, 373 \class{BezierSpline} objects. 374 \end{classdesc} 375 376 \begin{classdesc}{SurfaceLoop}{list} 377 creates a loop of \class{PlaneSurface} or \class{RuledSurface}, which defines 378 the shell of a volume. 379 \end{classdesc} 380 381 \subsubsection{Lines} 382 \begin{classdesc}{Line}{point1, point2} 383 creates a line between two points. 384 \end{classdesc} 385 386 \begin{methoddesc}[Line]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}} 387 defines the number of elements on the line. If set, it overwrites the local 388 length setting which would be applied. 389 The progression factor \var{progression} defines the change of element size 390 between neighboured elements. If \var{createBump} is set progression is 391 applied towards the centre of the line. 392 \end{methoddesc} 393 394 \begin{methoddesc}[Line]{resetElementDistribution}{} 395 removes a previously set element distribution from the line. 396 \end{methoddesc} 397 \begin{methoddesc}[Line]{getElementDistribution}{} 398 returns the element distribution as a tuple of number of elements, progression 399 factor and bump flag. If no element distribution is set None is returned. 400 \end{methoddesc} 401 402 \subsubsection{Splines} 403 \begin{classdesc}{Spline}{point0, point1, ...} 404 A spline curve defined by a list of points \var{point0}, \var{point1}, \ldots 405 \end{classdesc} 406 407 \begin{methoddesc}[Spline]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}} 408 defines the number of elements on the spline. If set, it overwrites the local 409 length setting which would be applied. 410 The progression factor \var{progression} defines the change of element size 411 between neighboured elements. If \var{createBump} is set progression is 412 applied towards the centre of the spline. 413 \end{methoddesc} 414 415 \begin{methoddesc}[Spline]{resetElementDistribution}{} 416 removes a previously set element distribution from the spline. 417 \end{methoddesc} 418 419 \begin{methoddesc}[Spline]{getElementDistribution}{} 420 returns the element distribution as a tuple of number of elements, progression 421 factor and bump flag. If no element distribution is set None is returned. 422 \end{methoddesc} 423 424 \subsubsection{BSplines} 425 \begin{classdesc}{BSpline}{point0, point1, ...} 426 A B-spline curve defined by a list of points \var{point0}, \var{point1}, \ldots 427 \end{classdesc} 428 429 \begin{methoddesc}[BSpline]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}} 430 defines the number of elements on the curve. If set, it overwrites the local 431 length setting which would be applied. 432 The progression factor \var{progression} defines the change of element size 433 between neighboured elements. If \var{createBump} is set progression is 434 applied towards the centre of the curve. 435 \end{methoddesc} 436 437 \begin{methoddesc}[BSpline]{resetElementDistribution}{} 438 removes a previously set element distribution from the curve. 439 \end{methoddesc} 440 441 \begin{methoddesc}[BSpline]{getElementDistribution}{} 442 returns the element distribution as a tuple of number of elements, progression 443 factor and bump flag. If no element distribution is set None is returned. 444 \end{methoddesc} 445 446 \subsubsection{Bezier Curves} 447 \begin{classdesc}{BezierCurve}{point0, point1, ...} 448 A Bezier spline curve defined by a list of points \var{point0}, \var{point1}, \ldots 449 \end{classdesc} 450 451 \begin{methoddesc}[BezierCurve]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}} 452 defines the number of elements on the curve. If set, it overwrites the local 453 length setting which would be applied. 454 The progression factor \var{progression} defines the change of element size 455 between neighboured elements. If \var{createBump} is set progression is 456 applied towards the centre of the curve. 457 \end{methoddesc} 458 459 \begin{methoddesc}[BezierCurve]{resetElementDistribution}{} 460 removes a previously set element distribution from the curve. 461 \end{methoddesc} 462 463 \begin{methoddesc}[BezierCurve]{getElementDistribution}{} 464 returns the element distribution as a tuple of number of elements, progression 465 factor and bump flag. If no element distribution is set None is returned. 466 \end{methoddesc} 467 468 \subsubsection{Arcs} 469 \begin{classdesc}{Arc}{centre_point, start_point, end_point} 470 creates an arc by specifying a centre for a circle and start and end points. 471 An arc may subtend an angle of at most $\pi$ radians. 472 \end{classdesc} 473 474 \begin{methoddesc}[Arc]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}} 475 defines the number of elements on the arc. If set, it overwrites the local 476 length setting which would be applied. 477 The progression factor \var{progression} defines the change of element size 478 between neighboured elements. If \var{createBump} is set progression is 479 applied towards the centre of the arc. 480 \end{methoddesc} 481 482 \begin{methoddesc}[Arc]{resetElementDistribution}{} 483 removes a previously set element distribution from the arc. 484 \end{methoddesc} 485 486 \begin{methoddesc}[Arc]{getElementDistribution}{} 487 returns the element distribution as a tuple of number of elements, progression 488 factor and bump flag. If no element distribution is set None is returned. 489 \end{methoddesc} 490 491 \subsubsection{Plane surfaces} 492 \begin{classdesc}{PlaneSurface}{loop, \optional{holes=[list]}} 493 creates a plane surface from a \class{CurveLoop}, which may have one or more 494 holes described by a \var{list} of \class{CurveLoop} objects. 495 \end{classdesc} 496 497 \begin{methoddesc}[PlaneSurface]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}} 498 defines the number of elements on all lines. 499 \end{methoddesc} 500 501 \begin{methoddesc}[PlaneSurface]{setRecombination}{\optional{max_deviation=45 * \var{DEG} }} 502 the mesh generator will try to recombine triangular elements into 503 quadrilateral elements. \var{max_deviation} (in radians) defines the 504 maximum deviation of any angle in the quadrilaterals from the right angle. 505 Set \var{max_deviation}=\var{None} to remove recombination. 506 \end{methoddesc} 507 508 \begin{methoddesc}[PlaneSurface]{setTransfiniteMeshing}{\optional{orientation="Left"}} 509 applies 2D transfinite meshing to the surface. 510 \var{orientation} defines the orientation of triangles. Allowed values 511 are \var{Left''}, \var{Right''} and \var{Alternate''}. 512 The boundary of the surface must be defined by three or four lines and an 513 element distribution must be defined on all faces where opposite 514 faces use the same element distribution. No holes must be present. 515 \end{methoddesc} 516 517 \subsubsection{Ruled Surfaces} 518 \begin{classdesc}{RuledSurface}{list} 519 creates a surface that can be interpolated using transfinite interpolation. 520 \var{list} gives a list of three or four lines defining the boundary of the 521 surface. 522 \end{classdesc} 523 524 \begin{methoddesc}[RuledSurface]{setRecombination}{\optional{max_deviation=45 * \var{DEG} }} 525 the mesh generator will try to recombine triangular elements into 526 quadrilateral elements. \var{max_deviation} (in radians) defines the 527 maximum deviation of any angle in the quadrilaterals from the right angle. 528 Set \var{max_deviation}=\var{None} to remove recombination. 529 \end{methoddesc} 530 531 \begin{methoddesc}[RuledSurface]{setTransfiniteMeshing}{\optional{orientation="Left"}} 532 applies 2D transfinite meshing to the surface. 533 \var{orientation} defines the orientation of triangles. Allowed values 534 are \var{Left''}, \var{Right''} and \var{Alternate''}. 535 The boundary of the surface must be defined by three or four lines and an 536 element distribution must be defined on all faces where opposite 537 faces use the same element distribution. No holes must be present. 538 \end{methoddesc} 539 540 \begin{methoddesc}[RuledSurface]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}} 541 defines the number of elements on all lines. 542 \end{methoddesc} 543 544 \subsubsection{Volumes} 545 \begin{classdesc}{Volume}{loop, \optional{holes=[list]}} 546 creates a volume given a \class{SurfaceLoop}, which may have one or more holes 547 define by the list of \class{SurfaceLoop}. 548 \end{classdesc} 549 550 \begin{methoddesc}[Volume]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}} 551 defines the number of elements on all lines. 552 \end{methoddesc} 553 554 \begin{methoddesc}[Volume]{setRecombination}{\optional{max_deviation=45 * \var{DEG} }} 555 the mesh generator will try to recombine triangular elements into 556 quadrilateral elements. These meshes are then used to generate the volume mesh 557 if possible. 558 Together with transfinite meshing one can construct rectangular meshes. 559 \var{max_deviation} (in radians) defines the maximum deviation of any angle in 560 the quadrilaterals from the right angle. 561 Set \var{max_deviation}=\var{None} to remove recombination. 562 \end{methoddesc} 563 564 \begin{methoddesc}[Volume]{setTransfiniteMeshing}{\optional{orientation="Left"}} 565 applies transfinite meshing to the volume and all surfaces (if 566 \var{orientation} is not equal to \var{None}). 567 \var{orientation} defines the orientation of triangles. Allowed values 568 are \var{Left''}, \var{Right''} and \var{Alternate''}. 569 The boundary of the surface must be defined by three or four lines and an 570 element distribution must be defined on all faces where opposite 571 faces use the same element distribution. 572 If \var{orientation} is equal to \var{None} transfinite meshing is not 573 switched on for the surfaces but needs to be set by the user. 574 No holes must be present. 575 \textbf{Warning: The functionality of transfinite meshing without 576 recombination is not entirely clear in \gmshextern. 577 So please apply this method with care.} 578 \end{methoddesc} 579 580 %============================================================================ 581 \subsection{Transformations} 582 583 Sometimes it is convenient to create an object and then make copies at 584 different orientations or in different sizes. This can be achieved by 585 applying transformations which are used to move geometrical objects in the 586 3-dimensional space and to resize them. 587 588 \begin{classdesc}{Translation}{\optional{b=[0,0,0]}} 589 defines a translation $x \to x+b$. \var{b} can be any object that can be 590 converted into a \numpy object of shape $(3,)$. 591 \end{classdesc} 592 593 \begin{classdesc}{Rotation}{\optional{axis=[1,1,1], \optional{ point = [0,0,0], \optional{angle=0*RAD} } } } 594 defines a rotation by \var{angle} around axis through point \var{point} and direction \var{axis}. 595 \var{axis} and \var{point} can be any object that can be converted into a 596 \numpy object of shape $(3,)$. 597 \var{axis} does not have to be normalised but must have positive length. 598 The right-hand rule~\cite{RIGHTHANDRULE} applies. 599 \end{classdesc} 600 601 \begin{classdesc}{Dilation}{\optional{factor=1., \optional{centre=[0,0,0]}}} 602 defines a dilation by the expansion/contraction \var{factor} with 603 \var{centre} as the dilation centre. 604 \var{centre} can be any object that can be converted into a \numpy object of 605 shape $(3,)$. 606 \end{classdesc} 607 608 \begin{classdesc}{Reflection}{\optional{normal=[1,1,1], \optional{offset=0}}} 609 defines a reflection on a plane defined in normal form $n^t x = d$ 610 where $n$ is the surface normal \var{normal} and $d$ is the plane \var{offset}. 611 \var{normal} can be any object that can be converted into a \numpy object of 612 shape $(3,)$. 613 \var{normal} does not have to be normalised but must have positive length. 614 \end{classdesc} 615 616 \begin{datadesc}{DEG} 617 a constant to convert from degrees to an internal angle representation in 618 radians. For instance use \code{90*DEG} for $90$ degrees. 619 \end{datadesc} 620 621 \subsection{Properties} 622 If you are building a larger geometry you may find it convenient to 623 create it in smaller pieces and then assemble them. 624 Property Sets make this easy, and they allow you to name the smaller 625 pieces for convenience. 626 627 Property Sets are used to bundle a set of geometrical objects in a 628 group. The group is identified by a name. Typically a Property Set 629 is used to mark subregions which share the same material properties or 630 to mark portions of the boundary. For efficiency, the \Design class 631 assigns an integer to each of its Property Sets, a so-called tag\index{tag}. 632 The appropriate tag is attached to the elements at generation time. 633 634 See the file \code{pycad/examples/quad.py} for an example using a {\it PropertySet}. 635 636 \begin{classdesc}{PropertySet}{name,*items} 637 defines a group geometrical objects which can be accessed through a \var{name} 638 The objects in the tuple \var{items} mast all be \ManifoldOneD, \ManifoldTwoD or \ManifoldThreeD objects. 639 \end{classdesc} 640 641 \begin{methoddesc}[PropertySet]{getManifoldClass}{} 642 returns the manifold class \ManifoldOneD, \ManifoldTwoD or \ManifoldThreeD 643 expected from the items in the property set. 644 \end{methoddesc} 645 646 \begin{methoddesc}[PropertySet]{getDim}{} 647 returns the spatial dimension of the items in the property set. 648 \end{methoddesc} 649 650 \begin{methoddesc}[PropertySet]{getName}{} 651 returns the name of the set 652 \end{methoddesc} 653 654 \begin{methoddesc}[PropertySet]{setName}{name} 655 sets the name. This name should be unique within a \Design. 656 \end{methoddesc} 657 658 \begin{methoddesc}[PropertySet]{addItem}{*items} 659 adds a tuple of items. They need to be objects of class \ManifoldOneD, 660 \ManifoldTwoD or \ManifoldThreeD. 661 \end{methoddesc} 662 663 \begin{methoddesc}[PropertySet]{getItems}{} 664 returns the list of items 665 \end{methoddesc} 666 667 \begin{methoddesc}[PropertySet]{clearItems}{} 668 clears the list of items 669 \end{methoddesc} 670 671 \begin{methoddesc}[PropertySet]{getTag}{} 672 returns the tag used for this property set 673 \end{methoddesc} 674 675 \section{Interface to the mesh generation software} 676 %\declaremodule{extension}{esys.pycad.gmsh} 677 %\modulesynopsis{Python geometry description and meshing interface} 678 679 The class and methods described here provide an interface to the mesh 680 generation software, which is currently \gmshextern. This interface could be 681 adopted to \emph{triangle} or another mesh generation package if this is 682 deemed to be desirable in the future. 683 684 \begin{classdesc}{Design}{ 685 \optional{dim=3, \optional{element_size=1., \optional{order=1, \optional{keep_files=False}}}}} 686 describes the geometry defined by primitives to be meshed. 687 \var{dim} specifies the spatial dimension, while \var{element_size} defines 688 the global element size which is multiplied by the local scale to set the 689 element size at each \Point. 690 The argument \var{order} defines the element order to be used. 691 If \var{keep_files} is set to \True temporary files are kept, otherwise they 692 are removed when the instance of the class is deleted. 693 \end{classdesc} 694 695 \begin{memberdesc}[Design]{GMSH} 696 gmsh file format~\cite{GMSH} 697 \end{memberdesc} 698 699 \begin{memberdesc}[Design]{IDEAS} 700 I-DEAS universal file format~\cite{IDEAS} 701 \end{memberdesc} 702 703 \begin{memberdesc}[Design]{VRML} 704 VRML file format, \cite{VRML} 705 \end{memberdesc} 706 707 \begin{memberdesc}[Design]{STL} 708 STL file format~\cite{STL} 709 \end{memberdesc} 710 711 \begin{memberdesc}[Design]{NASTRAN} 712 NASTRAN bulk data format~\cite{NASTRAN} 713 \end{memberdesc} 714 715 \begin{memberdesc}[Design]{MEDIT} 716 Medit file format~\cite{MEDIT} 717 \end{memberdesc} 718 719 \begin{memberdesc}[Design]{CGNS} 720 CGNS file format~\cite{CGNS} 721 \end{memberdesc} 722 723 \begin{memberdesc}[Design]{PLOT3D} 724 Plot3D file format~\cite{PLOT3D} 725 \end{memberdesc} 726 727 \begin{memberdesc}[Design]{DIFFPACK} 728 Diffpack 3D file format~\cite{DIFFPACK} 729 \end{memberdesc} 730 731 \begin{memberdesc}[Design]{DELAUNAY} 732 the Delaunay triangulator, see \gmshextern and \cite{TETGEN} 733 \end{memberdesc} 734 735 \begin{memberdesc}[Design]{MESHADAPT} 736 the gmsh triangulator, see \gmshextern 737 \end{memberdesc} 738 739 \begin{memberdesc}[Design]{FRONTAL} 740 the NETGEN~\cite{NETGEN} triangulator 741 \end{memberdesc} 742 743 \begin{methoddesc}[Design]{generate}{} 744 generates the mesh file. The data are written to the file \var{Design.getMeshFileName}. 745 \end{methoddesc} 746 747 \begin{methoddesc}[Design]{setDim}{\optional{dim=3}} 748 sets the spatial dimension which needs to be $1$, $2$ or $3$. 749 \end{methoddesc} 750 751 \begin{methoddesc}[Design]{getDim}{} 752 returns the spatial dimension. 753 \end{methoddesc} 754 755 \begin{methoddesc}[Design]{setElementOrder}{\optional{order=1}} 756 sets the element order which needs to be $1$ or $2$. 757 \end{methoddesc} 758 759 \begin{methoddesc}[Design]{getElementOrder}{} 760 returns the element order. 761 \end{methoddesc} 762 763 \begin{methoddesc}[Design]{setElementSize}{\optional{element_size=1}} 764 sets the global element size. The local element size at a point is defined as 765 the global element size multiplied by the local scale. 766 The element size must be positive. 767 \end{methoddesc} 768 769 \begin{methoddesc}[Design]{getElementSize}{} 770 returns the global element size. 771 \end{methoddesc} 772 773 \begin{methoddesc}[Design]{setKeepFilesOn}{} 774 work files are kept at the end of the generation. 775 \end{methoddesc} 776 777 \begin{methoddesc}[Design]{setKeepFilesOff}{} 778 work files are deleted at the end of the generation. 779 \end{methoddesc} 780 781 \begin{methoddesc}[Design]{keepFiles}{} 782 returns \True if work files are kept, \False otherwise. 783 \end{methoddesc} 784 785 \begin{methoddesc}[Design]{setScriptFileName}{\optional{name=None}} 786 sets the file name for the gmsh input script. 787 If no name is given a name with extension "geo" is generated. 788 \end{methoddesc} 789 790 \begin{methoddesc}[Design]{getScriptFileName}{} 791 returns the name of the file for the gmsh script. 792 \end{methoddesc} 793 794 \begin{methoddesc}[Design]{setMeshFileName}{\optional{name=None}} 795 sets the name for the mesh file. If no name is given a name is generated. 796 The format is set by \\\var{Design.setFileFormat}. 797 \end{methoddesc} 798 799 \begin{methoddesc}[Design]{getMeshFileName}{} 800 returns the name of the mesh file. 801 \end{methoddesc} 802 803 \begin{methoddesc}[Design]{addItems}{*items} 804 adds the tuple of \var{items}. An item can be any primitive or a 805 \class{PropertySet}. 806 \warning{If a \PropertySet is added which includes an object that is not 807 part of a \PropertySet, it may not be considered in the meshing.} 808 \end{methoddesc} 809 810 \begin{methoddesc}[Design]{getItems}{} 811 returns a list of the items. 812 \end{methoddesc} 813 814 \begin{methoddesc}[Design]{clearItems}{} 815 resets the items in this design. 816 \end{methoddesc} 817 818 \begin{methoddesc}[Design]{getMeshHandler}{} 819 returns a handle to the mesh. 820 Calling this method generates the mesh from the geometry and returns a 821 mechanism to access the mesh data. In the current implementation this 822 method returns a file name for a file containing the mesh data. 823 \end{methoddesc} 824 825 \begin{methoddesc}[Design]{getScriptString}{} 826 returns the gmsh script to generate the mesh as a string. 827 \end{methoddesc} 828 829 \begin{methoddesc}[Design]{getCommandString}{} 830 returns the gmsh command used to generate the mesh as a string. 831 \end{methoddesc} 832 833 \begin{methoddesc}[Design]{setOptions}{ 834 \optional{optimize_quality=\True 835 \optional{, smoothing=1 836 \optional{, curvature_based_element_size=\False\\ 837 \optional{, algorithm2D=\var{Design.MESHADAPT} 838 \optional{, algorithm3D=\var{Design.FRONTAL}\\ 839 \optional{, generate_hexahedra=False}}}}}}} 840 sets options for the mesh generator. 841 \var{algorithm2D} sets the 2D meshing algorithm to be used. 842 The algorithm needs to be \var{Design.DELAUNAY}, \var{Design.FRONTAL}, 843 or \var{Design.MESHADAPT}. By default \var{Design.MESHADAPT} is used. 844 \var{algorithm3D} sets the 3D meshing algorithm to be used. 845 The algorithm needs to be \var{Design.DELAUNAY} or \var{Design.FRONTAL}. 846 By default \var{Design.FRONTAL} is used. 847 \var{optimize_quality}=\True invokes an optimization of the mesh quality. 848 \var{smoothing} sets the number of smoothing steps to be applied to the mesh. 849 \var{curvature_based_element_size}=\True switches on curvature based 850 definition of element size. 851 \var{generate_hexahedra}=\True switches on the usage of quadrilateral or 852 hexahedral elements. 853 \end{methoddesc} 854 855 \begin{methoddesc}[Design]{getTagMap}{} 856 returns a \class{TagMap} to map the \class{PropertySet} names to tag numbers 857 generated by gmsh. 858 \end{methoddesc} 859 860 \begin{methoddesc}[Design]{setFileFormat}{\optional{format=\var{Design.GMSH}}} 861 sets the file format. \var{format} must be one of the values:\\ 862 \var{Design.GMSH}\\ 863 \var{Design.IDEAS}\\ 864 \var{Design.VRML}\\ 865 \var{Design.STL}\\ 866 \var{Design.NASTRAN}\\ 867 \var{Design.MEDIT}\\ 868 \var{Design.CGNS}\\ 869 \var{Design.PLOT3D}\\ 870 \var{Design.DIFFPACK}. 871 \end{methoddesc} 872 873 \begin{methoddesc}[Design]{getFileFormat}{} 874 returns the file format. 875 \end{methoddesc} 876