Getting Started With SFEPY - Finite Element Analysis in Python

by Sukhbinder Singh

This tutorial offers a quick look at using SFEPY to perform finite element analysis in python. This tutorial is meant to give anyone a quick start in sfepy.

The tutorial was written as part of CSRP program of AeSIAA which the author is volunteering as a mentor to few aeronautical students.

Thanks to Robert Cimrman for helping me in this tutorial and for bringing Sfepy to the world

Import the basics

In [2]:
import sfepy

from sfepy.base.base import Struct
from sfepy.discrete import (FieldVariable, Material, Integral, Integrals,
                            Equation, Equations, Problem)
from sfepy.discrete.fem import FEDomain, Field,Mesh
from sfepy.terms import Term
from sfepy.discrete.conditions import Conditions, EssentialBC
from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson,lame_from_youngpoisson
from sfepy.mechanics.tensors import get_von_mises_stress


import numpy as np
from IPython.display import Image
In [3]:
print sfepy.__version__
2014.2

Our Simple Plate Problem

A square plate with
  • bottom edge constrained
  • Unit load applied in Y direction on the top edge
In [42]:
Image('/Users/Sukhbinder/Documents/sfepy-simple_tutorial/simple_plate_geometry.png')
Out[42]:

Read a finite element mesh, that defines the domain $\Omega$

In [5]:
# Load the Mesh
mesh=Mesh.from_file('/Users/Sukhbinder/sqarePlate.vtk')
sfepy: reading mesh [line2, tri3, quad4, tetra4, hexa8] (/Users/Sukhbinder/sqarePlate.vtk)...
sfepy: ...done in 0.01 s
Mesh.from_file(fromfile) Read a mesh from a file.
There is a mesh generator module available in sfepy for simple mesh generation, to explore visit mesh generators

Create a domain. The domain allows defining regions or subdomains.

In [6]:
# Setup the Domain and Regions
domain = FEDomain('domain', mesh)
FEDomain(name, mesh defining the domain) - creates a domain

Define the regions - the whole domain $\Omega$, where the solution is sought, and bottom, top, where the boundary conditions will be applied.

In [7]:
# We call the entire domain omega
omega = domain.create_region('Omega', 'all')

#  We define bottom
bottom = domain.create_region('Bottom','vertices in (y < 0.001)','facet')

# and Top
top = domain.create_region('Top','vertices in (y > 0.99)','facet')
Regions assign names to various parts of the finite element mesh. The region names can later be referred to, for example when specifying portions of the mesh to apply boundary conditions to.

Regions can be specified in a variety of ways, including by element or by node.

domain.create_region(name_of_region,selection_criteria)

selection_criteria can be 'all', or ' vertices in (y < 0.001)' or 'cells in group 1'`

For the materials lets define some values
In [8]:
young=1.0
poisson=0.0
density=1.0
loadval=1.0
The material is linear elastic and its properties are specified as Lame parameters

Define the material

In [9]:
lam,mu = lame_from_youngpoisson(young,poisson,plane='stress')
mtx_d = stiffness_from_youngpoisson(2, young, poisson, plane='stress' )

m = Material('m',lam=lam,mu=mu, D=mtx_d, rho=density)

f = Material('f', val=loadval)
Sfepy solve PDEs if finite elements are used for discretization in space and the problem is expressed as a variational problem

Materials are used to define constitutive parameters (e.g. stiffness, permeability, or viscosity), and other non-field arguments of terms (e.g. known traction or volume forces).

Material(name,Constant_material_values_passed_by_their_names) - Class to hold constitutive and other material parameters

Depending on a particular term, the parameters can be constants, functions defined over FE mesh nodes, functions defined in the elements,

Define the field variables for the variation problem
A field is used to define the approximation on a domain
A displacement field (three DOFs/node) will be computed on a region
In [10]:
field = Field.from_args('fu', np.float64, 'vector',omega,approx_order=1)
u = FieldVariable('u', 'unknown', field)
v = FieldVariable('v', 'test', field, primary_var_name='u')
Field.from_args(field_name,type_of_field,shape_of_field,region_name,fe_approximation_order)
  • shape : int/tuple/str
  • The field shape: 1 or (1,) or 'scalar', space dimension (2, or (2,) or 3 or (3,)) or 'vector', or a tuple. The field shape determines the shape of the FE base functions and is related to the number of components of variables and to the DOF per node count, depending on the field kind.
  • region : Region The region where the field is defined.

  • approx_order : int/str The FE approximation order, e.g. 0, 1, 2, '1B' (1 with bubble).

FieldVariable(name,kind,field) - defines a finite element field variable.

kind can be 'unknown' and 'test'

Define ehe numerical quadrature that will be used to integrate each term over its domain
In [11]:
# Define the integral type Volume/Surface and quadrature rule
integral = Integral('i', order=2)

Integral(name, order)

Define the variational problem
Equations in SfePy are built using terms, which correspond directly to the integral forms of weak formulation of a problem to be solved.
In [12]:
# The weak formulation of the linear elastic problem.
t1 = Term.new('dw_lin_elastic_iso(m.lam, m.mu,v, u)', integral, omega, m=m, v=v, u=u)
t2 = Term.new('dw_surface_ltr(f.val, v)', integral, top, f=f, v=v)
eq1 = Equation('balance_of_forces', t1-t2)

lhs_eqs = Equations([eq1])
Save the problem as plate
In [13]:
pb = Problem('plate', equations=lhs_eqs)

Problem(name,equations) -Problem definition, the top-level class holding all data necessary to solve a problem.

Define Boundary Conditions
In [14]:
fixed_b = EssentialBC('FixedB', bottom, {'u.all' : 0.0})
The essential boundary conditions set values of DOFs in some regions, while the constraints constrain or transform values of DOFs in some regions.

EssentialBC(name,region,dofs) Essential boundary condidion

  • name : str The boundary condition name.
  • region : Region instance The region where the boundary condition is applied.
  • dofs : dict The boundary condition specification defining the constrained DOFs and their values.
Apply the boundary conditions
In [15]:
pb.time_update(ebcs=Conditions([fixed_b])) 
sfepy: updating variables...
sfepy: ...done
sfepy: setting up dof connectivities...
sfepy: ...done in 0.00 s
sfepy: matrix shape: (222, 222)
sfepy: assembling matrix graph...
sfepy: ...done in 0.00 s
sfepy: matrix structural nonzeros: 2796 (5.67e-02% fill)
Solve the problem.
In [16]:
state = pb.solve()
sfepy: updating materials...
sfepy:     m
sfepy:     f
sfepy: ...done in 0.01 s
sfepy: nls: iter: 0, residual: 3.239418e-01 (rel: 1.000000e+00)
sfepy:   rezidual:    0.00 [s]
sfepy:      solve:    0.00 [s]
sfepy:     matrix:    0.00 [s]
sfepy: nls: iter: 1, residual: 3.136819e-15 (rel: 9.683281e-15)
Lets Look at the problem definition
In [17]:
print pb
Problem:plate
  conf:
    Struct
  domain:
    FEDomain:domain
  ebcs:
    Conditions
  epbcs:
    Conditions
  equations:
    Equations
  evaluator:
    BasicEvaluator
  fields:
    dict with keys: ['fu']
  file_per_var:
    False
  float_format:
    None
  functions:
    None
  graph_changed:
    True
  integrals:
    None
  lcbcs:
    Conditions
  linearization:
    Struct
  matrix_hook:
    None
  mtx_a:
    (222, 222) spmatrix of float64, 2796 nonzeros
  name:
    plate
  nls_iter_hook:
    None
  nls_status:
    IndexedStruct
  ofn_trunk:
    domain
  output_dir:
    .
  output_format:
    vtk
  output_modes:
    dict with keys: ['vtk', 'h5']
  solvers:
    Struct:solvers
  ts:
    TimeStepper
Save the result of the solve in a vtk file for visualization
In [18]:
out = state.create_output_dict()
pb.save_state('plate.vtk', out=out)
In [19]:
# Check outout
print out
{'u': Struct:output_data}
In [20]:
# Check solve state
print state
State
  r_vec:
    None
  variables:
    Variables
  vec:
    (242,) array of float64

View the results

In [35]:
from sfepy.postprocess.viewer import Viewer
view = Viewer('plate.vtk')
In [36]:
view()
sfepy: point vectors u at [ 0.  0.  0.]
sfepy: range: -1.50e-15 1.00e+00 l2 norm range: 0.00e+00 1.00e+00
Out[36]:
<sfepy.postprocess.viewer.ViewerGUI at 0x11d9f7ef0>
In [37]:
Image('/Users/Sukhbinder/Documents/sfepy-simple_tutorial/sfepy_simple_plate_simple.png')
Out[37]:
Not as fancy as we would like, so check viewer options
In [39]:
view(vector_mode='warp_norm', rel_scaling=.1,is_scalar_bar=True, is_wireframe=True ,\
     rel_text_width=.1)
sfepy: point vectors u at [ 0.  0.  0.]
sfepy: range: -1.50e-15 1.00e+00 l2 norm range: 0.00e+00 1.00e+00
Out[39]:
<sfepy.postprocess.viewer.ViewerGUI at 0x123628050>
In [40]:
Image('/Users/Sukhbinder/Documents/sfepy-simple_tutorial/sfepy_result_simple_plate.png')
Out[40]:
Save the regions we created and show them with viewer
In [23]:
pb.save_regions_as_groups("plateregions")
view = Viewer('plateregions.vtk')
view(is_scalar_bar=True, is_wireframe=True ,rel_text_width=.2)
sfepy: saving regions as groups...
sfepy:   Omega
sfepy:   Bottom
sfepy:   Top
sfepy: ...done
sfepy: point scalars Bottom at [ 0.  0.  0.]
sfepy: range: 0.00e+00 1.00e+00 l2 norm range: 0.00e+00 1.00e+00
sfepy: point scalars Omega at [ 1.1  0.   0. ]
sfepy: range: 1.00e+00 1.00e+00 l2 norm range: 1.00e+00 1.00e+00
sfepy: point scalars Top at [ 2.2  0.   0. ]
sfepy: range: 0.00e+00 1.00e+00 l2 norm range: 0.00e+00 1.00e+00
Out[23]:
<sfepy.postprocess.viewer.ViewerGUI at 0x1207136b0>
In [41]:
Image('/Users/Sukhbinder/Documents/sfepy-simple_tutorial/sfepy_regions_plot.png')
Out[41]:
Evaluate the strain and stress
In [25]:
strain = pb.evaluate('ev_cauchy_strain.2.Omega(u)', mode='el_avg')
stress = pb.evaluate('ev_cauchy_stress.2.Omega(m.D, u)', mode='el_avg')

# Strain
out['cauchy_strain'] = Struct(name='output_data', mode='cell',data=strain, dofs=None)

# Stress
out['cauchy_stress'] = Struct(name='output_data', mode='cell',data=stress, dofs=None)

# Von mises Stress
vms = get_von_mises_stress(stress.squeeze())
vms.shape = (vms.shape[0], 1, 1, 1)
out['von_mises_stress'] = Struct(name='output_data', mode='cell',data=vms, dofs=None)
sfepy: equation "tmp":
sfepy: ev_cauchy_strain.2.Omega(u)
sfepy: updating materials...
sfepy: ...done in 0.00 s
sfepy: equation "tmp":
sfepy: ev_cauchy_stress.2.Omega(m.D, u)
sfepy: updating materials...
sfepy:     m
sfepy: ...done in 0.00 s
Check what has been calculated
In [26]:
print out 
{'cauchy_stress': Struct:output_data, 'cauchy_strain': Struct:output_data, 'u': Struct:output_data, 'von_mises_stress': Struct:output_data}
Dump solution to a file in VTK format for visualization
In [27]:
pb.save_state('plate_stress_strain.vtk', out=out,extend=True)
View the Stress/Strain/Von Mises Stress
In [28]:
view = Viewer('plate_stress_strain.vtk')
view(vector_mode='warp_norm', rel_scaling=.1,is_scalar_bar=True, is_wireframe=True ,rel_text_width=.1)
sfepy: point vectors u at [ 0.   1.1  0. ]
sfepy: range: -1.50e-15 1.00e+00 l2 norm range: 0.00e+00 1.00e+00
sfepy: cell scalars von_mises_stress at [ 1.1  1.1  0. ]
sfepy: range: 1.00e+00 1.00e+00 l2 norm range: 1.00e+00 1.00e+00
sfepy: cell tensors cauchy_strain at [ 0.  0.  0.]
sfepy: range: -9.67e-15 1.00e+00 l2 norm range: 1.00e+00 1.00e+00
sfepy: cell tensors cauchy_stress at [ 1.1  0.   0. ]
sfepy: range: -4.84e-15 1.00e+00 l2 norm range: 1.00e+00 1.00e+00
Out[28]:
<sfepy.postprocess.viewer.ViewerGUI at 0x1207fb5f0>
In [43]:
Image('/Users/Sukhbinder/Documents/sfepy-simple_tutorial/Sfepy_Stress_Strain.png')
Out[43]:

Question for you?

What is happening here?, Why do the Stress / Strain are all constants?
Hint: check material
Do the same problem with the following materials and check results
  • young=1E5.0
  • poisson=0.3
  • density=1000.0
  • loadval=10.0
Get primary variable for this problem, the displacement vector
In [30]:
u=state()
In [31]:
print u.reshape((121,2))[:100,1]
[ 0.          0.          1.          1.          0.          0.          0.
  0.          0.          0.          0.          0.          0.1111111
  0.2222222   0.3333333   0.4444444   0.5555556   0.6666667   0.7777778
  0.8888889   1.          1.          1.          1.          1.          1.
  1.          1.          0.8888889   0.7777778   0.6666667   0.5555556
  0.4444444   0.3333333   0.2222222   0.1111111   0.5175066   0.2569767
  0.261411    0.8215258   0.6775941   0.8359691   0.9200555   0.1463759
  0.08217954  0.8552634   0.7883398   0.4677138   0.2381719   0.5307371
  0.4299782   0.1149718   0.5198091   0.7735901   0.6304843   0.9165288
  0.9002651   0.3690425   0.1120099   0.09271977  0.2749267   0.1684781
  0.8243221   0.2005828   0.06923705  0.1718829   0.4185083   0.534285
  0.6115431   0.7404318   0.5697951   0.4398279   0.5944333   0.412475
  0.7272115   0.4219179   0.6291189   0.6454831   0.3266441   0.3417535
  0.6957667   0.2244955   0.5342121   0.6763426   0.3669423   0.9015062
  0.1811105   0.4506553   0.7798817   0.8280784   0.355945    0.4048789
  0.3275286   0.4031655   0.8393816   0.6367107   0.8121094   0.4537974
  0.3237081   0.5185756 ]
In [32]:
print u.reshape((121,2))[101:121,1]
[ 0.7197816   0.09507792  0.1491696   0.7097158   0.2676832   0.2902304
  0.4956254   0.9328416   0.08017429  0.5482404   0.09772739  0.9282306
  0.7338181   0.1985919   0.9150268   0.08161845  0.9329709   0.8956666
  0.6212125   0.6199571 ]
In [33]:
print state
State
  r_vec:
    None
  variables:
    Variables
  vec:
    (242,) array of float64

Quick Recap

  • $Mesh.from$$file(from$$file)$ - Read a mesh from a file.
  • $FEDomain(name, mesh.defining.the.domain)$ - Creates a domain
  • $domain.create$_$region(name.of.region,selection.criteria)$ - Creates Region ex: 'all', or ' vertices in (y < 0.001)'
  • $Material(name,Constant.material.values.passed.by.their.names)$ - Class to hold constitutive and other material parameters
  • $Field.from{\_}args(field$_$name,type.of.field,shape.of.field,region.name,fe.approximation.order)$ - define the field variables for the variation problem
  • $FieldVariable(name,kind,field)$ - defines a finite element field variable. Kind can be 'unknown' and 'test'
  • $Integral(name, order)$ - Define the integral type Volume/Surface and quadrature rule
  • $Term.new($'dw_lin_elastic_iso$(m.lam, m.mu,v, u)', integral, omega, m=m, v=v, u=u)$ - Define variation problem equation
  • $Problem(name,equations)$ - Problem definition, the top-level class holding all data necessary to solve a problem.

Further Material

Thank you!!

Ignore the below lines

In [1]:
from IPython.core.display import HTML
def css_styling():
    styles = open("custom.css", "r").read()
    return HTML(styles)
css_styling()
Out[1]:
Below code just to get Mass and Stiffness matrices
In [ ]:
t1 = Term.new('dw_lin_elastic(m.D, v, u)', integral, omega, m=m, v=v, u=u)
t2 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega, m=m, v=v, u=u)
eq1 = Equation('stiffness', t1)
eq2 = Equation('mass', t2)
lhs_eqs = Equations([eq1, eq2])
pbm = Problem('massstiff', equations=lhs_eqs)
pbm.time_update()
pbm.update_materials()


# Assemble stiffness and mass matrices.                                                         
mtx_k = eq1.evaluate(mode='weak', dw_mode='matrix', asm_obj=pbm.mtx_a)
mtx_m = mtx_k.copy()
mtx_m.data[:] = 0.0
mtx_m = eq2.evaluate(mode='weak', dw_mode='matrix', asm_obj=mtx_m)
print mtx_k
print mtx_m
print mesh.coors
print mesh.conns