Fwd: Re: Re: Several questions for solving a linear elastic problem
Forwarding to the list...
 Forwarded Message  Subject: Re: [sfepy] Re: Several questions for solving a linear elastic problem Date: Tue, 29 Aug 2017 10:59:10 +0200 From: Kathrin Sobe <kathrinsbriefkasten@googlemail.com> To: Robert Cimrman <cimrman3@ntc.zcu.cz>
Hi Robert,
just to let you know: This is the result of the example code after I changed the things according to your suggestions. Thank you.
I think the questions are answered and I also get a nice result picture. In the next phase I will try the same with the real mesh object.
My modifications are as follows: per triangle (facet). I hope this is the right way to distributet the
 the surface region is extracted directly with: surface = domain.create_region('Toppart', 'vertices of surface', 'facet') # 528 facets
 the material function returns a force vector for each quadrature point: I just splitted the triangle force vector by the number of quadrature points
forces.
 result extraction with post processing: this works after adding copy_materials = False: The second approach from examples/linear_elasticity/modal_analysis.py, around line 291 did not work in my example. The code is there in comment lines. The code was working, but the result was zero for all elements? However with post_processing I am able to access the results.
 Furthermore, I changed the definition of the linear elastic material to solid = Material('solid', values={'D': stiffness_from_youngpoisson(dim=3, young=200, poisson=0.3)}). In my case the material is defined by the youngsmodulus and the poissions ratio.
 the gravity is changed to: t3 = Term.new('dw_volume_lvf(f.val, v)', integral, omega, f=f, v=v) The values for the material are also defined by a function.
Here is the complete modified code example:
# coding: utf8 from __future__ import print_function from __future__ import absolute_import from argparse import ArgumentParser import dill import functools import numpy import time
import sys from sfepy.discrete.common.mappings import get_physical_qps sys.path.append('.')
from sfepy import data_dir from sfepy.base.base import IndexedStruct, Struct from sfepy.discrete import (FieldVariable, Material, Integral, Integrals, 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_youngpoisson
"""u aus: http://sfepy.org/docdevel/tutorial.html """ help = { 'show' : 'show the results figure', }
r""" Linear elasticity with pressure traction load on a surface and constrained to onedimensional motion.
Siehe Beispiel als Skript hier!!!!!!!!!1 http://sfepy.org/docdevel/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.
"""
# 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. problem.equations.materials[0]
stress = problem.evaluate('ev_cauchy_stress.2.Omega(solid.D, u)',
mode='el_avg', copy_materials=False)
out['cauchy_stress'] = Struct(name='output_data',
mode='cell', data=stress,
dofs=None) # # size 1782 mal 6 Werte
Pro Element eine Zeile (=1782) return out
def _get_loadcase(): u""" Returns a list of facets and a list of force vectors. The order defines the matching between facet and force vector. """ facets = dill.load(open(fname, "r")) # #force vectors from another example: temporary solution: take the necessary number of vectors force_vectors_example = dill.load(open(fname2, "r")) force_vectors = force_vectors_example[:len(facets)] return facets, force_vectors
def _set_not_top_region_zero(facets, force_vectors, facet_indices_top): u""" Surface facets that are not in the top region get the zero force vector. """ zero_vector = numpy.array([0.0, 0.0, 0.0]) for i, f in enumerate(facets): if f not in facet_indices_top: force_vectors[i] = zero_vector else: pass return facets, force_vectors
def _get_qp_force_vectors_from_loadcase(force_vectors, facets, region, qps, number_qps_facet): u""" Ensure order of region_facets and divide force vectors for quadrature points. Number of quadrature points per facets is always the same? Even if the triangles have a different size? """ force_vectors_qps = numpy.repeat(force_vectors, 3, axis=0) force_vectors_qps = force_vectors_qps / float(number_qps_facet) return force_vectors_qps
def linear_tension(ts, coors, mode=None, force_vectors_qp=None, **kwargs): u""" mapping between force vectors (one for every triangle) and quadrature points.
coors are quadrature points: here are three quadrature points per
tetrahedron
that means: the force vector is divided by three and set to every qp
what happens when the number of quadrature points is not even for all
triangles? """ if mode == 'qp': # shape: 128x3 # Array with tension values for each quadrature point (constant 0.8) # val = numpy.tile(10.0, (coor.shape[0], 1, 1)) # numpy.tile(1.0, (3, 1, 1)): [[[1]], [[1]], [[1]]] # val = numpy.tile([[0.001], [0.01], [1.8]], (coors.shape[0], 1, 1)) # 1584 mal [[0.5 0.45284157 0.4803068], ...] val = numpy.reshape(force_vectors_qp, (coors.shape[0], 3, 1)) return {'val' : val}
def get_gravity(ts, coors, mode=None, **kwargs): u""" Gravity in zdirection. Gravity is a volume load: here constant for every quadrature point. Mapping also necessary (see above) """ if mode == 'qp': val = numpy.tile([[0.0], [0.0], [10.5]], (coors.shape[0], 1, 1)) return {'val' : val}
def cube_displacement_medium(): u""" 3D Box with tetrahedron elements. Problem type: linear elastic problem
Input:
force vectors for the top triangles (top surface of the tetrahedron
mesh) and gravity for each tetrahedron.
Output:
 displacement for each vertex
 stress tensor
"""
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 = 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') # 88 facets
# gives you vertices of each face in the top region, that can be used
to index mesh.coors (mesh.coors[vof]). vertices_indexes_facets_top = mesh.cmesh.get_incident(0, top.faces, 2).reshape((1, 3)) # vertices_facets_top = mesh.coors[vertices_indexes_facets_top] # facet_indices_top = top.get_facet_indices() surface = domain.create_region('Toppart', 'vertices of surface', 'facet') # 528 facets # vertices of facets # defines the order of facets: order has to be applied on facets load vectors vertices_indexes_facets_surface = mesh.cmesh.get_incident(0, surface.faces, 2).reshape((1, 3)) # gives you vertices of each face in the top region, that can be used to index mesh.coors (mesh.coors[vof]).
z_min = bb[0][2]
z_max = bb[1][2]
z_range = z_max  z_min
z_tol = 0.01 * z_range
z_min_tol = z_min + z_tol
bottom = domain.create_region('Bottom',
'vertices in (z <= {})'.format(z_min_tol),
'vertex')
field = Field.from_args('displacement', numpy.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)}) # changed to the given matarial parameter emodulus and poisson solid = Material('solid', values={'D': stiffness_from_youngpoisson(dim=3, young=200, poisson=0.3)}) # young=29 oder 200
# # loadcase with forces on the surface
# # currently: 528 elements and 528 force vectors
facets, force_vectors = _get_loadcase()
facets, force_vectors = _set_not_top_region_zero(facets, force_vectors,
vertices_indexes_facets_top)
integral = Integral('i', order=2)
# # quadrature points in the order of surface facets
# # number of quadrature points for each facet defines the splitting
factor? # # e.g. three qps per facet: divide force vector by three # # returns a force vector array with one force vector for each qp qps1 = get_physical_qps(surface, integral) qps_shape = qps1.shape # (528, 3, 3) number_qps_facet = qps_shape[1] force_vectors_qp = _get_qp_force_vectors_from_loadcase(force_vectors, facets, surface, qps1, number_qps_facet) force_vectors_qp = force_vectors_qp * numpy.array([0, 0, 1]) linear_tension_function = functools.partial(linear_tension, force_vectors_qp=force_vectors_qp) load = Material('load', function=linear_tension_function) # material defined by values or function # volume forces f (gravity modelling)
f = Material('f', function=Function('get_gravity', get_gravity))
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, surface, load=load, v=v) # old: # t3 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega, m=solid, v=v, u=u) t3 = Term.new('dw_volume_lvf(f.val, v)', integral, omega, f=f, v=v)
# ???? is this right??
eq1 = Equation('balance_of_forces', t1 + t2 + t3)
eqs = Equations([eq1])
fixbottom = EssentialBC('fixbottom', bottom, {'u.all' : 0.0})
ls = ScipyDirect({})
nls_status = IndexedStruct()
# configurable residuum
# absolute_tolerance = 1e10
absolute_tolerance = 0.1
relative_tolerance = 1.0
nls = Newton({ 'i_max' : 1,
'eps_a' : absolute_tolerance,
'eps_r' : relative_tolerance,
'macheps' : 1e16,
# Linear system error < (eps_a * lin_red).
'lin_red' : 1e2,
'ls_red' : 0.1,
'ls_red_warp' : 0.001,
'ls_on' : 1.1,
'ls_min' : 1e5,
'check' : 0,
'delta' : 1e6, }, lin_solver=ls, status=nls_status)
pb = Problem('linear_elastic_traction', equations=eqs, nls=nls, ls=ls)
# pb.set_materials([load, solid])
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 berechnen: is working now!!! copy_materials = False
in function! # pb.save_state('linear_elastic_traction_cylinder.vtk', vec, post_process_hook=post_process)
name = 'linear_elastic_traction_medium_box.vtk'
pb.save_state(name, vec, post_process_hook=post_process)
######### stress and strain results are zero ####### # # alternative for calculating the stress tensor # state = pb.create_state() # aux = state.create_output_dict() # strain = pb.evaluate('ev_cauchy_strain.2.Omega(u)', # mode='el_avg', verbose=False) # # # aux['cauchy_strain'] = Struct(name='output_data', # mode='cell', data=strain, # dofs=None) # # # stress = pb.evaluate('ev_cauchy_stress.2.Omega(solid.D, u)', # mode='el_avg', verbose=False, copy_materials=False) # # aux['cauchy_stress'] = Struct(name='output_data', # mode='cell', data=stress, # dofs=None) # # size 1782 mal 6 Werte Pro Element eine Zeile (=1782)
if options.show:
view = Viewer(name)
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_nod
if __name__ == '__main__': runtime_small_cube, nodes_small_cube = cube_displacement_medium()
20170825 17:03 GMT+02:00 Robert Cimrman <cimrman3@ntc.zcu.cz>:
Hi Kathrin,
On 08/25/2017 04:42 PM, Kathrin Sobe wrote:
Hi Robert,
thanks a lot for your help.
20170825 14:37 GMT+02:00 Robert Cimrman <cimrman3@ntc.zcu.cz>:
Hi Kathrin,
On 08/25/2017 10:37 AM, Kathrin Sobe wrote:
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/docdevel/examples/linear_elasticity/linear _elastic_tractions.html
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?
Yes. The dw_surface_ltr term you use corresponds to the Neumann boundary conditions = applied forces/tractions/pressure in the weak form of the elasticity PDE. Because you need a positiondependent value of the load, you need to use a material parameter given by a function.
ok, than I will concentrate on that approach. Additionally, I want to model the gravity. For that I added the rho parameter to the material defintion and defined a mass equation like this: # # rho ist die Dichte: bei geringerer Dichte, wird der Block mehr verformt: ist klar solid = Material('solid', values={'D': stiffness_from_lame(dim=3, lam=5.769, mu=3.846), 'rho': 0.501}) 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, surface, load=load, v=v) t3 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega, m=solid, v=v, u=u) eq1 = Equation('elasticity', t1  t2) eq2 = Equation('mass', t3) eqs = Equations([eq1, eq2])
I am not sure if this is the right way to add the gravity to the problem.
Use the dw_volume_lvf term for volume forces just as you use dw_surface_ltr for surface forces. The material parameter here is a force per unit volume. The mass matrix would be useful in a dynamic simulation, not here.
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.
Is it going to be applied on the whole surface? Anyway, you can select the whole surface of a mesh with 'vertices of surface', for example:
surf = domain.create_region('surf', 'vertices of surface', 'facet')
Check the manual section 4.3.3 (on regions) for additional possibilities of defining regions and combining existing regions using various operations.
Or, if your mesh comes from some mesh generator/cad software, you can mark various vertex / cell groups there and use those for region definitions. This depends on the mesher/format you use  let me know if you need help with that.
The surface extraction is exactly what I need, cool. I found that out as well. Thanks for clarifying the indexing. I will try that.
Just as an information: Not all surface triangles have a force vector, only the triangles in the upper region of the geometry.
Then you may need to remove the extra faces from the surface region (or simply define the force zero there).
 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.
You can use vectors there too, for example:
val = nm.tile([[0.8], [0.1], [0.2]], (coor.shape[0], 1, 1))
To have a constant value per element, you need to know in which element or face each quadrature point is. This information is available in the material function (for example, kwargs['term'].region.faces contains the list of faces of the region, etc.), but the actual code depends on how your material function is represented  is it a table of (cellvalue) pairs, or something else?
At the moment my force defintion is exacly as you expected: cellvalue pairs. The cell is the surface triangle and the value is the force vector.
Is it that easy: Find the cell for each quadrature point and return the cell force vector as its load? Do I need to split the force vector of the cell between its quadrature points? Is there always exacly one cell for each quadrature point?
You need to provide a force per unit area vector (so yes, we may call it splitting) in each quadrature point. Each quadrature point is in exactly one face, yes. The number of quadrature point per face can be found as follows:
qps = t2.get_physical_qps() print qps.shape (88, 3, 3)
= 88 faces, 3 QPs per face, 3 coordinates (in 3D)
So this tells you, that the coor argument of linear_tension has 3 rows for each face (in the order of top.faces (= t2.region.faces).
r.
SfePy mailing list sfepy@python.org https://mail.python.org/mm3/mailman3/lists/sfepy.python.org/
Hi Kathrin,
On 08/29/2017 11:09 AM, Robert Cimrman wrote:
Forwarding to the list...
 Forwarded Message  Subject: Re: [sfepy] Re: Several questions for solving a linear elastic problem Date: Tue, 29 Aug 2017 10:59:10 +0200 From: Kathrin Sobe <kathrinsbriefkasten@googlemail.com> To: Robert Cimrman <cimrman3@ntc.zcu.cz>
Hi Robert,
just to let you know: This is the result of the example code after I changed the things according to your suggestions. Thank you.
Glad to hear that!
I think the questions are answered and I also get a nice result picture. In the next phase I will try the same with the real mesh object.
My modifications are as follows: per triangle (facet). I hope this is the right way to distributet the
 the surface region is extracted directly with: surface = domain.create_region('Toppart', 'vertices of surface', 'facet') # 528 facets
 the material function returns a force vector for each quadrature point: I just splitted the triangle force vector by the number of quadrature points
forces.
IMHO you should divide the force vector (magnitude) by the triangle area, and put that value to all the quadrature points in that triangle.
BTW. the current git master version supports getting volumes (areas for faces) of all mesh entities by simply calling: mesh.cmesh.get_volumes(dim), where dim=2 for faces.
 result extraction with post processing: this works after adding copy_materials = False: The second approach from examples/linear_elasticity/modal_analysis.py, around line 291 did not work in my example. The code is there in comment lines. The code was working, but the result was zero for all elements? However with post_processing I am able to access the results.
Note that you need to call
pb.save_state(name, out=out)
where out is the dict of the output data (aux in your code), and not
pb.save_state(name, vec, post_process_hook=post_process)
 Furthermore, I changed the definition of the linear elastic material to solid = Material('solid', values={'D': stiffness_from_youngpoisson(dim=3, young=200, poisson=0.3)}). In my case the material is defined by the youngsmodulus and the poissions ratio.
 the gravity is changed to: t3 = Term.new('dw_volume_lvf(f.val, v)', integral, omega, f=f, v=v) The values for the material are also defined by a function.
Looks OK.
Here is the complete modified code example:
I cannot run the example, as it uses dill to load some additional files. Also, please, try sending the code as an attachment, pasting it into an email introduced line breaks at wrong places. (Coding style note: try following PEP 8
 lines max. 79 characters long :))
Let us know how it goes with the real mesh. r.
participants (1)

Robert Cimrman