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 youngs-modulus 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: utf-8 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/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.
 """
 # 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 z-direction. 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 e-modulus 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 = 1e-10
 absolute_tolerance = 0.1
 relative_tolerance = 1.0
 nls = Newton({ 'i_max'      : 1,
                'eps_a'      : absolute_tolerance,
                'eps_r'      : relative_tolerance,
                '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.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()

2017-08-25 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.

2017-08-25 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/doc-devel/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:

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.

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

At the moment my force defintion is exacly as you expected: cell-value 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/