Hi Robert,
I try to solve a linear elastic problem. I give a example in a python file that is a modified version of the example given in http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_tractio...
The problem I want to solve is the following: I have a object, defined as a 3D mesh of tetrahedrons. In the test example I take the medium size 3D cube. The real 3D object is much bigger and has a more structured surface.
On the top of the object there are certain forces. These forces are defined as 3D vectors, one for each surface triangle. Additionally the object has certain material parameters given by density, poisson ratio and elastic modulus.
Working through the linear elastic examples I figured out two different types of solutions:
- modeling the forces at the surface as a material with a function (manual 4.3.9)
- modeling the forces at the surface as a boundary condition
So the first question is in general: Is this the right way to model my problem? At the moment I think solution type number one is the better approach?
Going through the first approach I have the following more detailed problems:
- Region definition: I need a way to define the surface triangles of the 3D object. Because of the structure of the object, this cannot be done with coordinate ranges.
- Material definition with force vectors: the 3D force vectors are given for each surface triangle (defined at the centroid of each triangle). In the material function examples I found only the definition of forces at quadrature points and only as scalar values.
- After solving the problem I need to extract the node displacement and the stress. For the stress calculation I tried to use a post_process function as visible in the example. But this was not working.
Thank you for your help. I hope the code in my example code makes the problems more clear.
Regards, Kathrin
******************** code example: solution type with force vectors defined in material parameter ******************* # coding: utf-8 from __future__ import print_function from __future__ import absolute_import from argparse import ArgumentParser import numpy as nm import time
import sys sys.path.append('.')
from sfepy.base.base import IndexedStruct from sfepy.discrete import (FieldVariable, Material, Integral, Equation, Equations, Problem, Function) from sfepy.discrete.fem import Mesh, FEDomain, Field from sfepy.terms import Term from sfepy.discrete.conditions import Conditions, EssentialBC from sfepy.solvers.ls import ScipyDirect from sfepy.solvers.nls import Newton from sfepy.postprocess.viewer import Viewer from sfepy.mechanics.matcoefs import stiffness_from_lame
"""u aus: http://sfepy.org/doc-devel/tutorial.html """ help = { 'show' : 'show the results figure', }
r""" Linear elasticity with pressure traction load on a surface and constrained to one-dimensional motion.
Siehe Beispiel als Skript hier!!!!!!!!!1 http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_tractio...
Find :math:\ul{u} such that:
.. math:: \int_{\Omega} D_{ijkl}\ e_{ij}(\ul{v}) e_{kl}(\ul{u}) = - \int_{\Gamma_{right}} \ul{v} \cdot \ull{\sigma} \cdot \ul{n} \;, \quad \forall \ul{v} \;,
where
.. math:: D_{ijkl} = \mu (\delta_{ik} \delta_{jl}+\delta_{il} \delta_{jk}) + \lambda \ \delta_{ij} \delta_{kl} \;.
and :math:\ull{\sigma} \cdot \ul{n} = \bar{p} \ull{I} \cdot \ul{n}
with given traction pressure :math:\bar{p}.
"""
def post_process(out, problem, state, extend=False): """ This will be called after the problem is solved.
Parameters
----------
out : dict
The output dictionary, where this function will store additional
data.
problem : Problem instance
The current Problem instance.
state : State instance
The computed state, containing FE coefficients of all the unknown
variables.
extend : bool
The flag indicating whether to extend the output data to the whole
domain. It can be ignored if the problem is solved on the whole
domain already.
Returns
-------
out : dict
The updated output dictionary.
"""
from sfepy.base.base import Struct
# Cauchy strain averaged in elements.
strain = problem.evaluate('ev_cauchy_strain.2.Omega(u)',
mode='el_avg')
out['cauchy_strain'] = Struct(name='output_data',
mode='cell', data=strain,
dofs=None)
# Cauchy stress averaged in elements.
stress = problem.evaluate('ev_cauchy_stress.2.Omega(solid.D, u)',
mode='el_avg')
out['cauchy_stress'] = Struct(name='output_data',
mode='cell', data=stress,
dofs=None)
return out# # coor: 264 x3 --> qp points in top region! def linear_tension(ts, coor, mode=None, **kwargs): if mode == 'qp': # shape: 128x3 # Array with tension values for each quadrature point (constant -0.8) val = nm.tile(-0.8, (coor.shape[0], 1, 1)) # numpy.tile(1.0, (3, 1, 1)): [[[1]], [[1]], [[1]]]
# ## with vectors it does not work
# val = nm.tile(nm.array([-0.01, 0.05, -0.8]), (coor.shape[0], 1, 1))
return {'val' : val}def cube_displacement_medium(): u""" 3D Box with tetrahedron elements. Problem type: linear elastic problem
linear element approach (later quadratic)
Input:
I have force vectors for the top triangles (top surface of the tetrahedron mesh)
and gravity for each tetrahedron.
Output:
- displacement for each vertex
- stress tensor
"""
from sfepy import data_dir
parser = ArgumentParser()
parser.add_argument('--version', action='version', version='%(prog)s')
parser.add_argument('-s', '--show',
action="store_true", dest='show',
default=True, help=help['show'])
options = parser.parse_args()
# load mesh
mesh = Mesh.from_file(data_dir + '/meshes/3d/cube_medium_tetra.mesh')
print("Mesh number of nodes: {}, number of tetrahedrons: {}".format(mesh.n_nod, mesh.n_el))
# Domain mit Mesh definieren: hier 448 Knoten, 1782 Elemente
domain = FEDomain('domain', mesh)
# why refine?
domain.refine()
omega = domain.create_region('Omega', 'all') # Gesamtbereich
# facets: bounding box here: -0.5, -0.5, -0.5 bis 0.5, 0.5, 0.5
bb = mesh.get_bounding_box()
top = domain.create_region('Top',
'vertices in (z >= 0.48)',
'facet')
bottom = domain.create_region('Bottom',
'vertices in (z <= -0.48)',
'facet')
field = Field.from_args('displacement', nm.float64, 3, omega, approx_order=1)
u = FieldVariable('u', 'unknown', field, order=0)
v = FieldVariable('v', 'test', field, primary_var_name='u')
solid = Material('solid', values={'D': stiffness_from_lame(dim=3, lam=5.769, mu=3.846)})
load = Material('load', function=linear_tension) # material defined by values or function
integral = Integral('i', order=2)
t1 = Term.new('dw_lin_elastic(solid.D, v, u)', integral, omega, solid=solid, v=v, u=u)
t2 = Term.new('dw_surface_ltr(load.val, v)', integral, top, load=load, v=v)
eq = Equation('elasticity', t1 - t2)
eqs = Equations([eq])
fixbottom = EssentialBC('fixbottom', bottom, {'u.all' : 0.0})
ls = ScipyDirect({})
nls_status = IndexedStruct()
nls = Newton({ 'i_max' : 1,
'eps_a' : 1e-10,
'eps_r' : 1.0,
'macheps' : 1e-16,
# Linear system error < (eps_a * lin_red).
'lin_red' : 1e-2,
'ls_red' : 0.1,
'ls_red_warp' : 0.001,
'ls_on' : 1.1,
'ls_min' : 1e-5,
'check' : 0,
'delta' : 1e-6, }, lin_solver=ls, status=nls_status)
pb = Problem('linear_elastic_traction', equations=eqs, nls=nls, ls=ls)
pb.save_regions_as_groups('regions')
# ebcs: elastic boundary conditions
pb.time_update(ebcs=Conditions([fixbottom]))
start_time = time.time()
vec = pb.solve()
run_time = time.time()-start_time
print("runtime medium problem {} {}: {}".format(mesh.n_nod, mesh.n_el, run_time))
print(nls_status)
# von mises stress calculation (not working)
pb.save_state('linear_elastic_simulation2.vtk', vec)
# pb.save_state('linear_elastic_traction_cylinder.vtk', vec, post_process_hook=post_process)
if options.show:
view = Viewer('linear_elastic_simulation2.vtk')
view(vector_mode='warp_norm', rel_scaling=2,
is_scalar_bar=True, is_wireframe=True)
pb.reset()
u.reset()
v.reset()
return run_time, mesh.n_nodif __name__ == '__main__': runtime_small_cube, nodes_small_cube = cube_displacement_medium()