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 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
print sfepy.__version__
2014.2
Image('/Users/Sukhbinder/Documents/sfepy-simple_tutorial/simple_plate_geometry.png')
# 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
# Setup the Domain and Regions
domain = FEDomain('domain', mesh)
# 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')
domain.create_region(name_of_region,selection_criteria)
selection_criteria can be 'all', or ' vertices in (y < 0.001)' or 'cells in group 1'`
young=1.0
poisson=0.0
density=1.0
loadval=1.0
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)
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).
Depending on a particular term, the parameters can be constants, functions defined over FE mesh nodes, functions defined in the elements,
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')
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).
kind can be 'unknown'* and 'test' *
# Define the integral type Volume/Surface and quadrature rule
integral = Integral('i', order=2)
Integral(name, order)
# 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])
pb = Problem('plate', equations=lhs_eqs)
Problem(name,equations) -Problem definition, the top-level class holding all data necessary to solve a problem.
fixed_b = EssentialBC('FixedB', bottom, {'u.all' : 0.0})
EssentialBC(name,region,dofs) Essential boundary condidion
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)
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)
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
out = state.create_output_dict()
pb.save_state('plate.vtk', out=out)
# Check outout
print out
{'u': Struct:output_data}
# Check solve state
print state
State r_vec: None variables: Variables vec: (242,) array of float64
from sfepy.postprocess.viewer import Viewer
view = Viewer('plate.vtk')
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
<sfepy.postprocess.viewer.ViewerGUI at 0x11d9f7ef0>
Image('/Users/Sukhbinder/Documents/sfepy-simple_tutorial/sfepy_simple_plate_simple.png')
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
<sfepy.postprocess.viewer.ViewerGUI at 0x123628050>
Image('/Users/Sukhbinder/Documents/sfepy-simple_tutorial/sfepy_result_simple_plate.png')
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
<sfepy.postprocess.viewer.ViewerGUI at 0x1207136b0>
Image('/Users/Sukhbinder/Documents/sfepy-simple_tutorial/sfepy_regions_plot.png')
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
print out
{'cauchy_stress': Struct:output_data, 'cauchy_strain': Struct:output_data, 'u': Struct:output_data, 'von_mises_stress': Struct:output_data}
pb.save_state('plate_stress_strain.vtk', out=out,extend=True)
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
<sfepy.postprocess.viewer.ViewerGUI at 0x1207fb5f0>
Image('/Users/Sukhbinder/Documents/sfepy-simple_tutorial/Sfepy_Stress_Strain.png')
u=state()
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 ]
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 ]
print state
State r_vec: None variables: Variables vec: (242,) array of float64
Ignore the below lines
from IPython.core.display import HTML
def css_styling():
styles = open("custom.css", "r").read()
return HTML(styles)
css_styling()
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