pycad allows noe to set the number of elements per line and the generation of quadrilateral elements over serfaces.

 1 ksteube 1811 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % 4 % Copyright (c) 2003-2008 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 15 ksteube 1318 \chapter{The Module \pycad} \label{PYCAD CHAP} 16 gross 993 17 gross 999 18 \section{Introduction} 19 20 ksteube 1066 \pycad provides a simple way to build a mesh for your finite element 21 simulation. You begin by building what we call a {\it Design} using 22 primitive geometric objects, and then to go on to build a mesh from 23 the {\it Design}. The final step of generating the mesh from a {\it 24 Design} uses freely available mesh generation software, such as \gmshextern. 25 gross 999 26 ksteube 1066 A {\it 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 what we 30 call a {\it Design}, you can build meshes for arbitrarily complex 2-D 31 and 3-D structures. 32 33 The example included below shows how to use {\it pycad} to create a 2-D mesh 34 in the shape of a trapezoid with a cutout area. 35 36 \begin{python} 37 from esys.pycad import * 38 from esys.pycad.gmsh import Design 39 from esys.finley import MakeDomain 40 41 # A trapezoid 42 p0=Point(0.0, 0.0, 0.0) 43 p1=Point(1.0, 0.0, 0.0) 44 p2=Point(1.0, 0.5, 0.0) 45 p3=Point(0.0, 1.0, 0.0) 46 l01=Line(p0, p1) 47 l12=Line(p1, p2) 48 l23=Line(p2, p3) 49 l30=Line(p3, p0) 50 c=CurveLoop(l01, l12, l23, l30) 51 52 # A small triangular cutout 53 x0=Point(0.1, 0.1, 0.0) 54 x1=Point(0.5, 0.1, 0.0) 55 x2=Point(0.5, 0.2, 0.0) 56 x01=Line(x0, x1) 57 x12=Line(x1, x2) 58 x20=Line(x2, x0) 59 cutout=CurveLoop(x01, x12, x20) 60 61 # Create the surface with cutout 62 s=PlaneSurface(c, holes=[cutout]) 63 64 # Create a Design which can make the mesh 65 d=Design(dim=2, element_size=0.05) 66 67 # Add the trapezoid with cutout 68 d.addItems(s) 69 70 # Create the geometry, mesh and Escript domain 71 d.setScriptFileName("trapezoid.geo") 72 d.setMeshFileName("trapezoid.msh") 73 domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1, optimizeLabeling=True) 74 75 # Create a file that can be read back in to python with mesh=ReadMesh(fileName) 76 domain.write("trapezoid.fly") 77 \end{python} 78 79 This example is included with the software in 80 \code{pycad/examples/trapezoid.py}. If you have gmsh installed you can 81 run the example and view the geometry and mesh with: 82 83 \begin{python} 84 python trapezoid.py 85 gmsh trapezoid.geo 86 gmsh trapezoid.msh 87 \end{python} 88 89 A \code{CurveLoop} is used to connect several lines into a single curve. 90 It is used in the example above to create the trapezoidal outline for the grid 91 and also for the triangular cutout area. 92 You can use any number of lines when creating a \code{CurveLoop}, but 93 the end of one line must be identical to the start of the next. 94 95 Sometimes you might see us write \code{-c} where \code{c} is a 96 \code{CurveLoop}. This is the reverse curve of the curve \code{c}. 97 It is identical to the original except that its points are traversed 98 in the opposite order. This may make it easier to connect two curves 99 in a \code{CurveLoop}. 100 101 The example python script above calls both 102 \code{d.setScriptFileName()} and \code{d.setMeshFileName()}. You need 103 only call these if you wish to save the gmsh geometry and mesh files. 104 105 Note that the underlying mesh generation software will not accept all 106 the geometries you can create with {\it pycad}. For example, {\it 107 pycad} will happily allow you to create a 2-D {\it Design} that is a 108 closed loop with some additional points or lines lying outside of the 109 enclosed area, but gmsh will fail to create a mesh for it. 110 111 112 113 114 115 116 gross 999 \section{\pycad Classes} 117 \declaremodule{extension}{esys.pycad} 118 gross 993 \modulesynopsis{Python geometry description and meshing interface} 119 120 gross 999 \subsection{Primitives} 121 gross 993 122 ksteube 1066 Some of the most commonly-used objects in {\it pycad} are listed here. For a more complete 123 list see the full API documentation. 124 gross 999 125 gross 2429 \begin{classdesc}{Point}{x=0.,y=0.,z=0.\optional{,local_scale=1.}} 126 Create a point with from coordinates with local characteristic length \var{local_scale} 127 gross 999 \end{classdesc} 128 129 ksteube 1066 \begin{classdesc}{Line}{point1, point2} 130 Create a line with between starting and ending points. 131 \end{classdesc} 132 gross 2429 \begin{methoddesc}[Line]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}} 133 Defines the number of elements on the line. If set it overwrites the local length setting which would be applied. The progression factor \var{progression} defines the change of element size between naighboured elements. If \var{createBump} is set 134 progression is applied towards the center of the line. 135 \end{methoddesc} 136 \begin{methoddesc}[Line]{resetElementDistribution}{} 137 removes a previously set element distribution from the line. 138 \end{methoddesc} 139 \begin{methoddesc}[Line]{getElemenofDistribution}{} 140 Returns the element distribution as tuple of 141 number of elements, progression factor and bump flag. If 142 no element distribution is set None is returned. 143 \end{methoddesc} 144 gross 999 145 gross 2429 146 \begin{classdesc}{Spline}{point0, point1, ...} 147 A spline curve defined by a list of points \var{point0}, \var{point1},.... 148 gross 999 \end{classdesc} 149 gross 2429 \begin{methoddesc}[Spline]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}} 150 Defines the number of elements on the line. If set it overwrites the local length setting which would be applied. The progression factor \var{progression} defines the change of element size between naighboured elements. If \var{createBump} is set 151 progression is applied towards the center of the line. 152 \end{methoddesc} 153 \begin{methoddesc}[Spline]{resetElementDistribution}{} 154 removes a previously set element distribution from the line. 155 \end{methoddesc} 156 \begin{methoddesc}[Spline]{getElemenofDistribution}{} 157 Returns the element distribution as tuple of 158 number of elements, progression factor and bump flag. If 159 no element distribution is set None is returned. 160 \end{methoddesc} 161 gross 999 162 gross 2429 \begin{classdesc}{BSpline}{point0, point1, ...} 163 A B-spline curve defined by a list of points \var{point0}, \var{point1},.... 164 ksteube 1066 \end{classdesc} 165 gross 2429 \begin{methoddesc}[BSpline]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}} 166 Defines the number of elements on the line. If set it overwrites the local length setting which would be applied. The progression factor \var{progression} defines the change of element size between naighboured elements. If \var{createBump} is set 167 progression is applied towards the center of the line. 168 \end{methoddesc} 169 \begin{methoddesc}[BSpline]{resetElementDistribution}{} 170 removes a previously set element distribution from the line. 171 \end{methoddesc} 172 \begin{methoddesc}[BSpline]{getElemenofDistribution}{} 173 Returns the element distribution as tuple of 174 number of elements, progression factor and bump flag. If 175 no element distribution is set None is returned. 176 \end{methoddesc} 177 gross 999 178 gross 2429 \begin{classdesc}{BezierCurve}{point0, point1, ...} 179 A Brezier spline curve defined by a list of points \var{point0}, \var{point1},.... 180 gross 999 \end{classdesc} 181 gross 2429 \begin{methoddesc}[BezierCurve]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}} 182 Defines the number of elements on the line. If set it overwrites the local length setting which would be applied. The progression factor \var{progression} defines the change of element size between naighboured elements. If \var{createBump} is set 183 progression is applied towards the center of the line. 184 \end{methoddesc} 185 \begin{methoddesc}[BezierCurve]{resetElementDistribution}{} 186 removes a previously set element distribution from the line. 187 \end{methoddesc} 188 \begin{methoddesc}[BezierCurve]{getElemenofDistribution}{} 189 Returns the element distribution as tuple of 190 number of elements, progression factor and bump flag. If 191 no element distribution is set None is returned. 192 \end{methoddesc} 193 gross 999 194 gross 2429 \begin{classdesc}{Arc}{center_point, start_point, end_point} 195 Create an arc by specifying a center for a circle and start and end points. An arc may subtend an angle of at most $\pi$ radians. 196 ksteube 1066 \end{classdesc} 197 gross 2429 \begin{methoddesc}[Arc]{setElementDistribution}{n\optional{,progression=1\optional{,createBump=\False}}} 198 Defines the number of elements on the line. If set it overwrites the local length setting which would be applied. The progression factor \var{progression} defines the change of element size between naighboured elements. If \var{createBump} is set 199 progression is applied towards the center of the line. 200 \end{methoddesc} 201 \begin{methoddesc}[Arc]{resetElementDistribution}{} 202 removes a previously set element distribution from the line. 203 \end{methoddesc} 204 \begin{methoddesc}[Arc]{getElemenofDistribution}{} 205 Returns the element distribution as tuple of 206 number of elements, progression factor and bump flag. If 207 no element distribution is set None is returned. 208 \end{methoddesc} 209 gross 999 210 ksteube 1066 \begin{classdesc}{CurveLoop}{list} 211 gross 2429 Create a closed curve from the \code{list}. of 212 \class{Line}, \class{Arc}, \class{Spline}, \class{BSpline}, 213 \class{BrezierSpline}. 214 gross 999 \end{classdesc} 215 216 ksteube 1066 \begin{classdesc}{PlaneSurface}{loop, \optional{holes=[list]}} 217 gross 2429 Create a plane surface from a \class{CurveLoop}, which may have one or more holes 218 described by \var{list} of \class{CurveLoop}. 219 ksteube 1066 \end{classdesc} 220 gross 2429 \begin{methoddesc}[PlaneSurface]{setRecombination}{max_deviation} 221 the mesh generator will try to recombine triangular elements 222 into quadrilateral elements. \var{max_deviation} (in radians) defines the 223 maximum deviation of any angle in the quadrilaterals from the right angle. 224 Set \var{max_deviation}=\var{None} to remove recombination. 225 \end{methoddesc} 226 \begin{methoddesc}[PlaneSurface]{setTransfiniteMeshing}{\optional{orientation="Left"}} 227 applies 2D transfinite meshing to the surface. 228 \var{orientation} defines the orientation of triangles. Allowed values 229 are \var{Left''}, \var{Right''} or \var{Alternate''}. The 230 boundary of the surface muist be defined by three or four lines where an 231 element distribution must be defined on all faces where opposite 232 faces uses the same element distribution. No holes must be present. 233 \end{methoddesc} 234 ksteube 1066 235 gross 2429 236 237 ksteube 1066 \begin{classdesc}{RuledSurface}{list} 238 Create a surface that can be interpolated using transfinite interpolation. 239 gross 2429 \var{list} gives a list of three or four lines defining the boundary of the 240 surface. 241 ksteube 1066 \end{classdesc} 242 gross 2429 \begin{methoddesc}[RuledSurface]{setRecombination}{max_deviation} 243 the mesh generator will try to recombine triangular elements 244 into quadrilateral elements. \var{max_deviation} (in radians) defines the 245 maximum deviation of any angle in the quadrilaterals from the right angle. 246 Set \var{max_deviation}=\var{None} to remove recombination. 247 \end{methoddesc} 248 \begin{methoddesc}[RuledSurface]{setTransfiniteMeshing}{\optional{orientation="Left"}} 249 applies 2D transfinite meshing to the surface. 250 \var{orientation} defines the orientation of triangles. Allowed values 251 are \var{Left''}, \var{Right''} or \var{Alternate''}. The 252 boundary of the surface muist be defined by three or four lines where an 253 element distribution must be defined on all faces where opposite 254 faces uses the same element distribution. No holes must be present. 255 \end{methoddesc} 256 ksteube 1066 257 gross 2429 258 ksteube 1066 \begin{classdesc}{SurfaceLoop}{list} 259 gross 2429 Create a loop of \class{PlaneSurface} or \class{RuledSurface}, which defines the shell of a volume. 260 ksteube 1066 \end{classdesc} 261 262 \begin{classdesc}{Volume}{loop, \optional{holes=[list]}} 263 gross 2429 Create a volume given a \class{SurfaceLoop}, which may have one or more holes 264 define by the list of \class{SurfaceLoop}. 265 ksteube 1066 \end{classdesc} 266 267 \begin{classdesc}{PropertySet}{list} 268 Create a PropertySet given a list of 1-D, 2-D or 3-D items. See the section on Properties below for more information. 269 \end{classdesc} 270 271 gross 999 %============================================================================================================ 272 \subsection{Transformations} 273 274 ksteube 1066 Sometimes it's convenient to create an object and then make copies at 275 different orientations and in different sizes. Transformations are 276 used to move geometrical objects in the 3-dimensional space and to 277 resize them. 278 gross 999 279 \begin{classdesc}{Translation}{\optional{b=[0,0,0]}} 280 defines a translation $x \to x+b$. \var{b} can be any object that can be converted 281 into a \numarray object of shape $(3,)$. 282 \end{classdesc} 283 284 \begin{classdesc}{Rotatation}{\optional{axis=[1,1,1], \optional{ point = [0,0,0], \optional{angle=0*RAD} } } } 285 defines a rotation by \var{angle} around axis through point \var{point} and direction \var{axis}. 286 \var{axis} and \var{point} can be any object that can be converted 287 into a \numarray object of shape $(3,)$. 288 \var{axis} does not have to be normalized but must have positive length. The right hand rule~\cite{RIGHTHANDRULE} 289 applies. 290 \end{classdesc} 291 292 293 \begin{classdesc}{Dilation}{\optional{factor=1., \optional{center=[0,0,0]}}} 294 defines a dilation by the expansion/contraction \var{factor} with 295 \var{center} as the dilation center. 296 \var{center} can be any object that can be converted 297 into a \numarray object of shape $(3,)$. 298 \end{classdesc} 299 300 \begin{classdesc}{Reflection}{\optional{normal=[1,1,1], \optional{offset=0}}} 301 defines a reflection on a plane defined in normal form $n^t x = d$ 302 where $n$ is the surface normal \var{normal} and $d$ is the plane \var{offset}. 303 \var{normal} can be any object that can be converted 304 into a \numarray object of shape $(3,)$. 305 \var{normal} does not have to be normalized but must have positive length. 306 \end{classdesc} 307 308 ksteube 1066 \begin{datadesc}{DEG} 309 A constant to convert from degrees to an internal angle representation in radians. For instance use \code{90*DEG} for $90$ degrees. 310 \end{datadesc} 311 312 gross 999 \subsection{Properties} 313 314 ksteube 1066 If you are building a larger geometry you may find it convenient to 315 create it in smaller pieces and then assemble them into the whole. 316 Property sets make this easy, and they allow you to name the smaller 317 pieces for convenience. 318 gross 999 319 ksteube 1066 Property sets are used to bundle a set of geometrical objects in a 320 group. The group is identified by a name. Typically a property set 321 is used to mark subregions with share the same material properties or 322 to mark portions of the boundary. For efficiency, the \Design class 323 object assigns a integer to each of its property sets, a so-called tag 324 \index{tag}. The appropriate tag is attached to the elements at 325 generation time. 326 gross 1044 327 ksteube 1066 See the file \code{pycad/examples/quad.py} for an example using a {\it PropertySet}. 328 329 330 gross 999 \begin{classdesc}{PropertySet}{name,*items} 331 defines a group geometrical objects which can be accessed through a \var{name} 332 The objects in the tuple \var{items} mast all be \ManifoldOneD, \ManifoldTwoD or \ManifoldThreeD objects. 333 \end{classdesc} 334 335 336 \begin{methoddesc}[PropertySet]{getManifoldClass}{} 337 returns the manifold class \ManifoldOneD, \ManifoldTwoD or \ManifoldThreeD expected from the items 338 in the property set. 339 \end{methoddesc} 340 341 \begin{methoddesc}[PropertySet]{getDim}{} 342 returns the spatial dimension of the items 343 in the property set. 344 \end{methoddesc} 345 346 \begin{methoddesc}[PropertySet]{getName}{} 347 returns the name of the set 348 \end{methoddesc} 349 350 \begin{methoddesc}[PropertySet]{setName}{name} 351 sets the name. This name should be unique within a \Design. 352 \end{methoddesc} 353 354 \begin{methoddesc}[PropertySet]{addItem}{*items} 355 adds a tuple of items. They need to be objects of class \ManifoldOneD, \ManifoldTwoD or \ManifoldThreeD. 356 \end{methoddesc} 357 358 \begin{methoddesc}[PropertySet]{getItems}{} 359 returns the list of items 360 \end{methoddesc} 361 362 \begin{methoddesc}[PropertySet]{clearItems}{} 363 clears the list of items 364 \end{methoddesc} 365 366 \begin{methoddesc}[PropertySet]{getTag}{} 367 returns the tag used for this property set 368 \end{methoddesc} 369 370 ksteube 1066 \section{Interface to the mesh generation software} 371 gross 999 \declaremodule{extension}{esys.pycad.gmsh} 372 \modulesynopsis{Python geometry description and meshing interface} 373 374 ksteube 1066 The class and methods described here provide an interface to the mesh 375 generation software, which is currently gmsh. This interface could be 376 adopted to triangle or another mesh generation package if this is 377 lgraham 1700 deemed to be desirable in the future. 378 ksteube 1066 379 gross 999 \begin{classdesc}{Design}{ 380 \optional{dim=3, \optional{element_size=1., \optional{order=1, \optional{keep_files=False}}}}} 381 The \class{Design} describes the geometry defined by primitives to be meshed. 382 The \var{dim} specifies the spatial dimension. The argument \var{element_size} defines the global 383 element size which is multiplied by the local scale to set the element size at each \Point. 384 The argument \var{order} defines the element order to be used. If \var{keep_files} is set to 385 \True temporary files a kept otherwise they are removed when the instance of the class is deleted. 386 \end{classdesc} 387 388 389 \begin{methoddesc}[Design]{setDim}{\optional{dim=3}} 390 sets the spatial dimension which needs to be $1$, $2$ or $3$. 391 \end{methoddesc} 392 393 \begin{methoddesc}[Design]{getDim}{} 394 returns the spatial dimension. 395 \end{methoddesc} 396 397 \begin{methoddesc}[Design]{setElementOrder}{\optional{order=1}} 398 sets the element order which needs to be $1$ or $2$. 399 \end{methoddesc} 400 401 \begin{methoddesc}[Design]{getElementOrder}{} 402 returns the element order. 403 \end{methoddesc} 404 405 406 \begin{methoddesc}[Design]{setElementSize}{\optional{element_size=1}} 407 set the global element size. The local element size at a point is defined as 408 lgraham 1700 the global element size multiplied by the local scale. The element size must be positive. 409 gross 999 \end{methoddesc} 410 411 412 \begin{methoddesc}[Design]{getElementSize}{} 413 returns the global element size. 414 \end{methoddesc} 415 416 \begin{memberdesc}[Design]{DELAUNAY} 417 ksteube 1066 the gmsh Delauny triangulator. 418 gross 999 \end{memberdesc} 419 420 \begin{memberdesc}[Design]{TETGEN} 421 the TetGen~\cite{TETGEN} triangulator. 422 \end{memberdesc} 423 424 gross 2371 \begin{memberdesc}[Design]{NETGEN} 425 gross 999 the NETGEN~\cite{NETGEN} triangulator. 426 \end{memberdesc} 427 428 \begin{methoddesc}[Design]{setKeepFilesOn}{} 429 work files are kept at the end of the generation. 430 \end{methoddesc} 431 432 \begin{methoddesc}[Design]{setKeepFilesOff}{} 433 work files are deleted at the end of the generation. 434 \end{methoddesc} 435 436 \begin{methoddesc}[Design]{keepFiles}{} 437 returns \True if work files are kept. Otherwise \False is returned. 438 \end{methoddesc} 439 440 \begin{methoddesc}[Design]{setScriptFileName}{\optional{name=None}} 441 ksteube 1066 set the filename for the gmsh input script. if no name is given a name with extension "geo" is generated. 442 gross 999 \end{methoddesc} 443 444 \begin{methoddesc}[Design]{getScriptFileName}{} 445 ksteube 1066 returns the name of the file for the gmsh script. 446 gross 999 \end{methoddesc} 447 448 449 \begin{methoddesc}[Design]{setMeshFileName}{\optional{name=None}} 450 ksteube 1066 sets the name for the gmsh mesh file. if no name is given a name with extension "msh" is generated. 451 gross 999 \end{methoddesc} 452 453 \begin{methoddesc}[Design]{getMeshFileName}{} 454 returns the name of the file for the gmsh msh 455 \end{methoddesc} 456 457 458 \begin{methoddesc}[Design]{addItems}{*items} 459 adds the tuple of var{items}. An item can be any primitive or a \class{PropertySet}. 460 \warning{If a \PropertySet is added as an item added object that are not 461 part of a \PropertySet are not considered in the messing. 462 } 463 464 \end{methoddesc} 465 466 \begin{methoddesc}[Design]{getItems}{} 467 returns a list of the items 468 \end{methoddesc} 469 470 \begin{methoddesc}[Design]{clearItems}{} 471 resets the items in design 472 \end{methoddesc} 473 474 \begin{methoddesc}[Design]{getMeshHandler}{} 475 returns a handle to the mesh. The call of this method generates the mesh from the geometry and 476 lgraham 1700 returns a mechanism to access the mesh data. In the current implementation this 477 ksteube 1066 method returns a file name for a gmsh file containing the mesh data. 478 gross 999 \end{methoddesc} 479 480 \begin{methoddesc}[Design]{getScriptString}{} 481 ksteube 1066 returns the gmsh script to generate the mesh as a string. 482 gross 999 \end{methoddesc} 483 484 \begin{methoddesc}[Design]{getCommandString}{} 485 ksteube 1066 returns the gmsh command used to generate the mesh as string. 486 gross 999 \end{methoddesc} 487 488 \begin{methoddesc}[Design]{setOptions}{\optional{algorithm=None, \optional{ optimize_quality=True,\optional{ smoothing=1}}}} 489 sets options for the mesh generator. \var{algorithm} sets the algorithm to be used. 490 The algorithm needs to be \var{Design.DELAUNAY} 491 \var{Design.TETGEN} 492 or \var{Design.NETGEN}. By default \var{Design.DELAUNAY} is used. \var{optimize_quality}=\True invokes an optimization of the mesh quality. \var{smoothing} sets the number of smoothing steps to be applied to the mesh. 493 \end{methoddesc} 494 495 \begin{methoddesc}[Design]{getTagMap}{} 496 ksteube 1066 returns a \class{TagMap} to map the name \class{PropertySet} in the class to tag numbers generated by gmsh. 497 gross 999 \end{methoddesc}