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/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:

1. modeling the forces at the surface as a material with a function (manual 4.3.9)
2. 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 position-dependent value of the load, you need to use a material parameter given by a function.

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.

• 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 (cell-value) pairs, or something else?

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

You need to use copy_materials=False option, as follows:

 stress = problem.evaluate('ev_cauchy_stress.2.Omega(solid.D, u)',
                           mode='el_avg', copy_materials=False)

Then

pb.save_state('linear_elastic_traction_cylinder.vtk', vec, post_process_hook=post_process)

works. It is also possible to compute directly the stress/strain as in examples/linear_elasticity/modal_analysis.py, around line 291. The post_process_hook way is more useful when using the problem description files - it is not required in an "interactive" script.

r.

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_nod

if __name__ == '__main__': runtime_small_cube, nodes_small_cube = cube_displacement_medium()