 Revision 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (3 years, 8 months ago) by jfenwick
File MIME type: application/x-tex
File size: 33469 byte(s)
Make everyone sad by touching all the files


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