Mesh generation with CALFEM

Included in the Python version of CALFEM is a mesh generation library based on GMSH. The library encapsulates the usage of GMSH transparently for the user. It will also parse the output from GMSH and create the necessary data structures required by CALFEM for solving finite element problems.

Mesh generation in CALFEM is divided in three steps:

  1. Defining the geometry of the model.

  2. Creating a finite element mesh with the desired elements and properties

  3. Extracting data structures that can be used by CALFEM.

The following sections describe these steps.

Required modules for geometry and mesh generation

To use the CALFEM geometry and mesh generation routines, we use the following import directives:

import calfem.geometry as cfg
import calfem.mesh as cfm
import calfem.vis_mpl as cfv

Defining geometry

Geometry in CALFEM is described using the Geometry class. A Geometry-object will hold all points, lines and surfaces describing the geometry. This object is also what is passed to the mesh generation routines in the following sections.

A Geometry-object, g, is created with the following code:

g = cfg.Geometry()

This creates a Geometry object which will be used to described our geometry. Next we define the points that will be used to define lines, splines or ellipses. In this example we define a simple triangle:

g.point([0.0, 0.0]) # point 0
g.point([5.0, 0.0]) # point 1
g.point([2.5, 4.0]) # point 2

The points are connected together with spline-elements. These can have 2 or three nodes. In this case we only use 2 node splines (lines):

g.spline([0, 1]) # line 0
g.spline([1, 2]) # line 1
g.spline([2, 0]) # line 2

Finally we create the surface by defining what lines make up the surface:

g.surface([0, 1, 2])

To display our geometry, we use the calfem.vis module:

cfv.drawGeometry(g)
cfv.showAndWait()

Running this example will show the following window with a simple rectangle:

_images/mesh1.png

Creating a mesh

To create a mesh we need to create a GmshMesh object and initialize this with our geometry:

mesh = cfm.GmshMesh(g)

Next, we need to set some desired properties on our mesh:

mesh.el_type = 3          # Element type is quadrangle
mesh.dofs_per_node = 1     # Degrees of freedom per node
mesh.el_size_factor = 0.15 # Element size Factor

The el_type property determines the element used for mesh generation. Elements that can be generated are:

  • 2 - 3 node triangle element

  • 3 - 4 node quadrangle element

  • 5 - 8 node hexahedron

  • 16 - 8 node second order quadrangle

The dofs_per_node defines the number of degrees of freedom for each node. el_size_factor determines the coarseness of the mesh.

To generate the mesh and at the same time get the needed data structures for use with CALFEM we call the .create() method of the mesh object:

coords, edof, dofs, bdofs, elementmarkers = mesh.create()

The returned data structures are:

  • coords - Element coordinates

  • edof - Element topology

  • dofs - Degrees of freedom per node

  • bdofs - Boundary degrees of freedom. Dictionary containing the dofs for each boundary marker. More on markers in the next section.

  • elementmarkers - List of integer markers. Row i contains the marker of element i. Can be used to determine what region an element is in.

To display the generated mesh we can use the drawMesh() function of the calfem.vis module:

cfv.figure()

# Draw the mesh.

cfv.drawMesh(
    coords=coords,
    edof=edof,
    dofs_per_node=mesh.dofsPerNode,
    el_type=mesh.elType,
    filled=True,
    title="Example 01"
        )

Running the example will produce the following mesh with quad elements:

_images/mesh2.png

Changing the elType property to 2 (mesh.elType = 2) will produce a mesh with triangle elements instead:

_images/mesh3.png

Specifying boundary markers

To assist in assigning boundary conditions, markers can be defined on the geometry, which can be used to identify which dofs are assigned to nodes, lines and surfaces.

In this example we add a marker, 10, to line number 2. Markers are added as a parameter to the .spline() method of the Geometry object as shown in the following code:

g.spline([0, 1]) # line 0
g.spline([1, 2]) # line 1
g.spline([2, 0], marker=10) # line 2 with marker 10

It is also possible to assign markers to points. The marker parameter is added to the .point() method of the Geometry object.

g.point([0.0, 0.0])             # point 0
g.point([5.0, 0.0], marker=20)  # point 1
g.point([2.5, 4.0])             # point 2

In the same way markers can be added to surfaces as well.

Extracting dofs defined by markers

To extract the dofs defined by the marker we use the bdofs dictionary returned when the mesh was created by the .create() method. If we print the bdofs dictionary we get the following output:

{20: [2], 0: [1, 2, ... , 67], 10: [1, 3, 68, ... , 98]}

If we examine the output we see that there is a key, 10, containing the dofs of the number 2 line. We also have the key 20 with a single dof 2 in this case. If the dofsPerNode property in the mesh generator was set to 2 the marker 20 would have contained 2 integers.

Complete example with a solver

To illustrate the workflow of the mesh generation modules we implement a complete 2D plane stress solver.

Updated module imports

We need to add some additional import directives, such as the core calfem module as well as the calfem.utils module. We will also need NumPy as well as the standard math routines:

import calfem.core as cfc
import calfem.geometry as cfg
import calfem.mesh as cfm
import calfem.vis_mpl as cfv
import calfem.utils as cfu

import numpy as np

Problem variables and constants

To make it easier to update our example we define a number of variables describing our problem. First some geometric parameters describing our module, in this case a simple rectangular beam:

l = 5.0     # length
h = 1.0     # height
t = 0.2     # thickness

Next, we define our material properties we will need later in the code:

v = 0.35            # Poisson
E = 2.1e9           # Young modulus
ptype = 1           # plane stress
ep = [ptype,t]
D = cfc.hooke(ptype, E, v)

To make it easier to read our code we define 3 constants, which we will use instead of numerical markers.

left_support = 10
right_support = 20
top_line = 30

Creating a Geometry object

We are now ready to create a Geometry object describing our geometry:

g = cfg.Geometry()

g.point([0.0, 0.0], marker = left_support)  # point 0
g.point([l, 0.0], marker = right_support)   # point 1
g.point([l, h])                                     # point 2
g.point([0.0, h])                                   # point 3

g.spline([0, 1])                                    # line 0
g.spline([1, 2])                                    # line 1
g.spline([2, 3], marker = top_line)                 # line 3
g.spline([3, 0])                                    # line 3

g.surface([0, 1, 2, 3])

To show the geometry you can use:

cfv.draw_geometry(g)
cfv.showAndWait()

The finished geometry is shown in below:

_images/tut2-1.png

Creating a quad mesh

A quadrilateral mesh is now created with the following code. Please not that we use the dofsPerNode property to specify 2 dofs per node as this is a mechanical example.

mesh = cfm.GmshMesh(g)

mesh.elType = 3             # Type of mesh
mesh.dofsPerNode = 2        # Factor that changes element sizes
mesh.elSizeFactor = 0.10    # Factor that changes element sizes

coords, edof, dofs, bdofs, elementmarkers = mesh.create()

As previous, to show the finished mesh you can use:

cfv.figure()
cfv.draw_mesh(
    coords=coords,
    edof=edof,
    dofs_per_node=mesh.dofsPerNode,
    el_type=mesh.elType,
    filled=True)
cfv.showAndWait()

The finished mesh is shown belo:

_images/exm0_2.png

Implementing a CALFEM solver

We now have the necessary input to implement a simple CALFEM solver. First, we create some convenience variables, nDofs (total number of dofs), ex and ey (element x and y coordinates).

nDofs = np.size(dofs)
ex, ey = cfc.coordxtr(edof, coords, dofs)

The global stiffness matrix can now be allocated:

K = np.zeros([nDofs,nDofs])

For larger problems please consider using sparse matrices instead.

To make the assembly loop less cluttered we use the zip() method to extract rows from edof, ex and ey to eltopo, elx and ely. The loop then becomes:

for eltopo, elx, ely in zip(edof, ex, ey):
    Ke = cfc.planqe(elx, ely, ep, D)
    cfc.assem(eltopo, K, Ke)

Please not the difference from standard MATLAB CALFEM that the assem routine does not require returning the K matrix on the left side as the assembly is done in place.

Next, we need to setup our boundary conditions. Two empty arrays are created, bc for storing the dof to prescribe and and a second bcVal for storing the prescribed values.

bc = np.array([],'i')
bcVal = np.array([],'f')

To prescribe a boundary condition the utility function applybc() is used. This function takes the boundary dictionary as input and applies the boundary condition to the correct dofs. Here we prescribe the left support as fixed and the right support fixed in y-direction.

bc, bcVal = cfu.applybc(bdofs, bc, bcVal, left_support, 0.0, 0)
bc, bcVal = cfu.applybc(bdofs, bc, bcVal, right_support, 0.0, 2)

Now, we have some data structures that can be used by CALFEM. If there is no force and displacement, nothing happens. If we applied force f and calculate the displacement, then we can solve the equation to obtain more data to plot. These next steps are left in the next examples.