.. highlight:: python :linenothreshold: 3 .. include:: links.inc .. _sec-preprocessing: Preprocessing: *FreeCAD*/*OpenSCAD* + *Gmsh* ============================================ .. only:: html .. contents:: Table of Contents :local: :backlinks: top Introduction ------------ There are several open source tools for preparing 2D and 3D finite element meshes like `Salome`_, `FreeCAD`_, `Gmsh`_, `Netgen`_, etc. Most of them are GUI based geometrical modeling and meshing environments/tools but they also usually allow using their libraries in user scripts. Some of the above mentioned tools are handy for solid modeling, some of them are great for meshing. This tutorial shows how to combine solid geometry modeling functions provided by *FreeCAD* or `OpenSCAD`_ with meshing functions of *Gmsh*. The collaboration of modeling, meshing and conversion tools and the workflow are illustrated in the following scheme. .. image:: images/preprocessing/workflow.png :width: 90 % :align: center Creating geometry using *FreeCAD* --------------------------------- Functionalities of *FreeCAD* are accessible to Python and can be used to define geometrical models in simple Python scripts. There is a tutorial related to `Python scripting in FreeCAD`_. The first step in creating a Python script is to set up a path to the *FreeCAD* libraries and import all required modules:: import sys FREECADPATH = '/usr/lib/freecad/lib/' sys.path.append(FREECADPATH) from FreeCAD import Base, newDocument import Part import Draft import ProfileLib.RegularPolygon as Poly Now, a new empty *FreeCAD* document can be defined as:: doc = newDocument() All new objects describing the geometry will be added to this document. In the following lines a geometrical model of a screwdriver handle will be created. Let's start by defining a sphere and a cylinder and join these objects into the one called ``uni``:: radius = 0.01 height = 0.1 cyl = doc.addObject("Part::Cylinder", "cyl") cyl.Radius = radius cyl.Height = height sph = doc.addObject("Part::Sphere", "sph") sph.Radius = radius uni = doc.addObject("Part::MultiFuse", "uni") uni.Shapes = [cyl, sph] Create a polygon, revolve it around the *z*-axis to create a solid and use the result as the cutting tool applied to ``uni`` object:: ske = doc.addObject('Sketcher::SketchObject', 'Sketch') ske.Placement = Base.Placement(Base.Vector(0, 0, 0), Base.Rotation(-0.707107, 0, 0, -0.707107)) Poly.makeRegularPolygon('Sketch', 5, Base.Vector(-1.2 * radius, 0.9 * height, 0), Base.Vector(-0.8 * radius, 0.9 * height, 0)) cut = doc.addObject("PartDesign::Revolution", "Revolution") cut.Sketch = ske cut.ReferenceAxis = (ske, ['V_Axis']) cut.Angle = 360.0 dif = doc.addObject("Part::Cut", "dif") dif.Base = uni dif.Tool = cut Create a cylinder, make a polar array of the cylinder objects and subtract it from the previous result:: cyl1 = doc.addObject("Part::Cylinder", "cyl1") cyl1.Radius = 0.2 * radius cyl1.Height = 1.1 * height cyl1.Placement = Base.Placement(Base.Vector(-1.1 * radius, 0, -0.2 * height), Base.Rotation(0, 0, 0, 1)) arr = Draft.makeArray(cyl1, Base.Vector(1, 0, 0), Base.Vector(0, 1, 0), 2, 2) arr.ArrayType = "polar" arr.NumberPolar = 6 dif2 = doc.addObject("Part::Cut", "dif2") dif2.Base = dif dif2.Tool = arr Create a middle hole for the screwdriver metal part:: cyl2 = doc.addObject("Part::Cylinder", "cyl2") cyl2.Radius = 0.3 * radius cyl2.Height = height dif3 = doc.addObject("Part::Cut", "dif3") dif3.Base = dif2 dif3.Tool = cyl2 Finally, recompute the geometry, export the part to the *STEP* file and save the document in *FreeCAD* format (not really needed for subsequent mesh generation, but may be useful for visualization and geometry check):: doc.recompute() Part.export([dif3], 'screwdriver_handle.step') doc.saveAs('screwdriver_handle.FCStd') A finite element mesh can be generated directly in *FreeCAD* using ``MeshPart`` module:: import MeshPart mesh = doc.addObject("Mesh::Feature", "Mesh") mesh.Mesh = MeshPart.meshFromShape(Shape=dif3.Shape, MaxLength=0.002) mesh.Mesh.write("./screwdriver_handle.bdf", "NAS", "mesh") The meshing function of ``MeshPart`` module is limited to triangular grids so it is better to use `Gmsh`_ mesh generator which can provide triangular and quadrilateral meshes in 2D or tetrahedral and hexahedral meshes in 3D. *Gmsh* allows to control the meshing process through a wide range of parameters. Meshing by *Gmsh* will be described in section :ref:`preprocessing_gmsh`. .. image:: images/preprocessing/sdh_freecad.png :width: 80 % :align: center The example of screwdriver handle: :download:`screwdriver_handle.py `. There are two simple ways how to discover Python calls of *FreeCAD* functions. You can enable "show script commands in python console" in ``Edit->Preferences->General->Macro`` and the Python console by selecting ``View->Views->Python Console`` and all subsequent operations will be printed in the console as the Python code. The second way is to switch on the macro recording function (``Macro->Macro recording ...``) which generates a Python script (*FCMacro* file) containing all the code related to actions in the *FreeCAD* graphical interface. Creating geometry using *OpenSCAD* ---------------------------------- The alternative tool for solid geometrical modeling is *OpenSCAD* - "The Programmers Solid 3D CAD Modeller". It has its own description language based on functional programming that is used to construct solid models using geometrical primitives similar to *FreeCAD*. Solid geometries can be exported to several file formats including *STL* and *CSG*. *OpenSCAD* allows solid modeling based on Constructive Solid Geometry (CSG) principles and extrusion of 2D objects into 3D. The model of a screwdriver handle presented in the previous section can be defined in *OpenSCAD* by the following code (:download:`screwdriver_handle.scad `):: radius = 0.01; height = 0.1; $fn = 50; difference() { difference() { difference() { union() { cylinder(center=false, h=height, r=radius); sphere(radius); }; translate([0, 0, 0.9*height]) rotate_extrude() polygon([[0.8*radius, 0], [1.8*radius, -0.577*radius], [1.8*radius, 0.577*radius]]); } cylinder(center=false, h=1.1*height, r=0.3*radius); } for (i = [1:6]) { rotate([0, 0, 360/6*i]) translate([-1.1*radius, 0.0, -0.2*height]) cylinder(center=false, h=1.1*height, r=0.2*radius); } } .. image:: images/preprocessing/sdh_openscad.png :width: 80 % :align: center To generate a finite element mesh of the solid geometry the model must be exported to a suitable file format. *OpenSCAD* has limited export options, but by using *FreeCAD* import/export functions, it is possible to find a workaround. The *OpenSCAD* model can be exported to the *CSG* file format and *FreeCAD* can be used as a mesh converter to the *STEP* format:: import sys sys.path.append('/usr/lib/freecad/lib/') sys.path.append('/usr/lib/freecad/Mod/OpenSCAD/') import FreeCAD import Part import importCSG importCSG.open('screwdriver_handle.csg') Part.export([FreeCAD.ActiveDocument.Objects[-1]], 'screwdriver_handle.step') .. _preprocessing_gmsh: Gmsh - generating finite element mesh ------------------------------------- *Gmsh* can create finite element meshes using geometrical models imported from *STEP*, *IGES* and *BRep* files (has to be compiled with *OpenCASCADE* support). The following *GEO* file imports ``screwdriver_handle.step`` file and defines a field controlling the mesh size (:download:`screwdriver_handle.geo `):: Merge "screwdriver_handle.step"; Field[1] = MathEval; Field[1].F = "0.002"; Background Field = 1; Now, run *Gmsh* generator and export the mesh into the *MSH* format in which all surface and volumetric elements are stored:: gmsh -3 -format msh -o screwdriver_handle.msh screwdriver_handle.geo By converting the *MSH* file into the *VTK* format using ``sfepy-convert``:: sfepy-convert -d 3 screwdriver_handle.msh screwdriver_handle.vtk the surface elements are discarded and only the volumetric mesh is preserved. .. image:: images/preprocessing/sdh_mesh.png :width: 40 % :align: center Note: planar 2D meshes ^^^^^^^^^^^^^^^^^^^^^^ To create a planar 2D mesh, such as .. image:: images/preprocessing/circle_in_square.png :width: 40 % :align: center that can be described by :download:`this ` *Gmsh* code, the mesh generator can be called as follows:: gmsh -2 -format msh -o circle_in_square.msh circle_in_square.geo This, however is not enough to create a truly 2D mesh - the created mesh vertices still have the third, :math:`z`, component which is equal to zero. In order to remove the third component, use:: sfepy-convert --2d circle_in_square.msh circle_in_square.h5 Now, in the resulting ``circle_in_square.h5``, each vertex has only two coordinates. Another way of generating the 2D mesh is to use the legacy VTK format as follows:: gmsh -2 -format vtk -o circle_in_square.vtk circle_in_square.geo sfepy-convert circle_in_square.vtk circle_in_square.h5 This is due to the fact that the legacy VTK does not support 2D vertices and so the :class:`VTKMeshIO ` reader tries to detect the planar geometry by comparing the :math:`z` components to zero - the ``--2d`` option of ``sfepy-convert`` is not needed in this case. Multipart models ---------------- Meshing models composed of parts with different material groups is a little bit tricky task. But there are some more or less general ways of doing that. Here, the method using functions of *Gmsh* for periodic meshes will be shown. The screwdriver handle example is extended by adding a screwdriver shank. The new part is composed of a cylinder trimmed at one end:: cyl3 = doc.addObject("Part::Cylinder", "cyl3") cyl3.Radius = 0.3 * radius cyl3.Height = height cyl3.Placement = Base.Placement(Base.Vector(0, 0, height), Base.Rotation(0, 0, 0, 1)) tip1 = doc.addObject("Part::Box", "tip1") tip1.Length = radius tip1.Width = 2 * radius tip1.Height = 3 * radius tip1.Placement = Base.Placement(Base.Vector(0, -radius, 1.71 * height), Base.Rotation(Base.Vector(0, 1, 0), -10), Base.Vector(0, 0, 3 * radius)) tip2 = doc.addObject("Part::Mirroring", "tip2") tip2.Source = tip1 tip2.Normal = (1, 0, 0) tip3 = doc.addObject("Part::MultiFuse", "tip3") tip3.Shapes = [tip1, tip2] dif4 = doc.addObject("Part::Cut", "dif4") dif4.Base = cyl3 dif4.Tool = tip3 uni2 = doc.addObject("Part::MultiFuse", "uni2") uni2.Shapes = [cyl2, dif4] The handle and shank are exported to the *STEP* file as two separated parts:: doc.recompute() Part.export([dif3, uni2], 'screwdriver_full.step') doc.saveAs('screwdriver_full.FCStd') The full screwdriver example (handle + shank): :download:`screwdriver_full.py `. To create a coincidence mesh on the handle and shank interface, it is necessary to identify the interface surfaces and declare them to be periodic in the *GEO* file. The identification has to be done manually in the *Gmsh* graphical interface. .. image:: images/preprocessing/sdf_gmsh1.png :width: 60% :align: center .. image:: images/preprocessing/sdf_gmsh2.png :width: 60% :align: center The input file for *Gmsh* is than as follows (:download:`screwdriver_full.geo `):: Merge "screwdriver_full.step"; Periodic Surface 5 {7} = 26 {67}; Periodic Surface 3 {6, 2, -6, 7} = 27 {68, 69, -68, 67}; Physical Volume(1) = {1}; Physical Volume(2) = {2}; Field[1] = MathEval; Field[1].F = "0.0015"; Background Field = 1; where the first pair of periodic surfaces corresponds to the common circle faces (bottom of the shank) and the second pair to the common cylindrical surfaces. See `Gmsh Reference manual`_ for details on periodic meshing. Using the above stated *GEO* file, *Gmsh* creates a mesh containing duplicate vertices on the handle/shank interface. These duplicate vertices can be removed during the conversion to the *VTK* format by giving ``--merge`` (or just ``-m``) argument to `convert_mesh.py` script:: sfepy-convert -m screwdriver_full.msh screwdriver_full.vtk In order to extract the cells by the physical groups use the conversion script with ``--save-per-mat`` argument:: sfepy-convert --save-per-mat screwdriver_full.vtk screwdriver.vtk It produces `screwdriver.vtk` contaning the original mesh and `screwdriver_matid_1.vtk`, `screwdriver_matid_2.vtk` files containing only the cells of a given physical group and all vertices of the original mesh. .. image:: images/preprocessing/sdf_mesh.png :width: 60 % :align: center When using *OpenSCAD*, define the full screwdriver geometry as (:download:`screwdriver_full.scad `):: radius = 0.01; height = 0.1; $fn = 50; module tip() { rotate([0, -10, 0]) translate([0, -radius, -3*radius]) cube([radius, 2*radius, 3*radius], center=false); } difference() { difference() { difference() { union() { cylinder(center=false, h=height, r=radius); sphere(radius); }; translate([0, 0, 0.9*height]) rotate_extrude() polygon([[0.8*radius, 0], [1.8*radius, -0.577*radius], [1.8*radius, 0.577*radius]]); } cylinder(center=false, h=height, r=0.3*radius); } for (i = [1:6]) { rotate([0, 0, 360/6*i]) translate([-1.1*radius, 0.0, -0.2*height]) cylinder(center=false, h=1.1*height, r=0.2*radius); } } union() { difference() { translate([0, 0, height]) cylinder(center=false, h=height, r=0.3*radius); translate([0, 0, 1.71*height + 3*radius]) union() { tip(); mirror ([1, 0, 0]) tip(); } } cylinder(center=false, h=height, r=0.3*radius); } and convert the *CSG* file to the *STEP* file by:: importCSG.open('screwdriver_full.csg') top_group = FreeCAD.ActiveDocument.Objects[-1] Part.export(top_group.OutList, 'screwdriver_full.step') Since the different tools for geometry definition have been used, the numbering of geometric objects may differ and the surface and edge numbers have to be changed in the *GEO* file:: Periodic Surface 5 {6} = 26 {66}; Periodic Surface 3 {5, 2, -5, 6} = 27 {67, 68, -67, 66}; Note: The numbering of objects may vary between *FreeCAD*, *OpenSCAD* and *Gmsh* versions.