# 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.

## 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


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

In [6]:
# Setup the Domain and Regions
domain = FEDomain('domain', mesh)


## 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 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


## 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)


##### 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,

##### 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)

##### 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?

##### Do the same problem with the following materials and check results
• young=1E5.0
• poisson=0.3
• density=1000.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.

# Thank you!!

Ignore the below lines

In [1]:
from IPython.core.display import HTML
def css_styling():
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
`