# linear_elasticity/linear_elastic_tractions.py¶

Description

Linear elasticity with pressure traction load on a surface and constrained to one-dimensional motion.

Find such that:

where

and with given traction pressure .

The function verify_tractions() is called after the solution to verify that the inner surface tractions correspond to the load applied to the external surface. Try running the example with different approximation orders and/or uniform refinement levels:

• the default options:

python simple.py examples/linear_elasticity/linear_elastic_tractions.py -O refinement_level=0 -d approx_order=1

• refine once:

python simple.py examples/linear_elasticity/linear_elastic_tractions.py -O refinement_level=1 -d approx_order=1

• use the tri-quadratic approximation (Q2):

python simple.py examples/linear_elasticity/linear_elastic_tractions.py -O refinement_level=0 -d approx_order=2

source code

r"""
Linear elasticity with pressure traction load on a surface and constrained to
one-dimensional motion.

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}

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

The function :func:verify_tractions() is called after the solution to verify
that the inner surface tractions correspond to the load applied to the external
surface. Try running the example with different approximation orders and/or uniform refinement levels:

- the default options::

python simple.py examples/linear_elasticity/linear_elastic_tractions.py -O refinement_level=0 -d approx_order=1

- refine once::

python simple.py examples/linear_elasticity/linear_elastic_tractions.py -O refinement_level=1 -d approx_order=1

- use the tri-quadratic approximation (Q2)::

python simple.py examples/linear_elasticity/linear_elastic_tractions.py -O refinement_level=0 -d approx_order=2
"""
from __future__ import absolute_import
import numpy as nm
from sfepy.base.base import output
from sfepy.mechanics.matcoefs import stiffness_from_lame

def linear_tension(ts, coor, mode=None, **kwargs):
if mode == 'qp':
val = nm.tile(1.0, (coor.shape[0], 1, 1))

return {'val' : val}

def verify_tractions(out, problem, state, extend=False):
"""
Verify that the inner surface tractions correspond to the load applied
to the external surface.
"""
from sfepy.mechanics.tensors import get_full_indices
from sfepy.discrete import Material, Function

)

def eval_force(region_name):
strain = problem.evaluate(
'ev_cauchy_strain_s.i.%s(u)' % region_name, mode='qp',
verbose=False,
)
D = problem.evaluate(
'ev_surface_integrate_mat.i.%s(solid.D, u)' % region_name,
mode='qp',
verbose=False,
)

normal = nm.array([1, 0, 0], dtype=nm.float64)

s2f = get_full_indices(len(normal))
stress = nm.einsum('cqij,cqjk->cqik', D, strain)
# Full (matrix) form of stress.
mstress = stress[..., s2f, 0]

# Force in normal direction.
force = nm.einsum('cqij,i,j->cq', mstress, normal, normal)

def get_force(ts, coors, mode=None, **kwargs):
if mode == 'qp':
return {'force' : force.reshape(coors.shape[0], 1, 1)}
aux = Material('aux', function=Function('get_force', get_force))

middle_force = - problem.evaluate(
'ev_surface_integrate_mat.i.%s(aux.force, u)' % region_name,
aux=aux,
verbose=False,
)
output('%s section axial force:' % region_name, middle_force)

eval_force('Left')
eval_force('Middle')
eval_force('Right')

return out

def define(approx_order=1):
"""Define the problem to solve."""
from sfepy import data_dir

filename_mesh = data_dir + '/meshes/3d/block.mesh'

options = {
'nls' : 'newton',
'ls' : 'ls',
'post_process_hook' : 'verify_tractions',
}

functions = {
'linear_tension' : (linear_tension,),
}

fields = {
'displacement': ('real', 3, 'Omega', approx_order),
}

materials = {
'solid' : ({'D': stiffness_from_lame(3, lam=5.769, mu=3.846)},),
}

variables = {
'u' : ('unknown field', 'displacement', 0),
'v' : ('test field', 'displacement', 'u'),
}

regions = {
'Omega' : 'all',
'Left' : ('vertices in (x < -4.99)', 'facet'),
# Use a parent region to select only facets belonging to cells in the
# parent region. Otherwise, each facet is in the region two times, with
# opposite normals.
'Middle' : ('vertices in (x > -1e-10) & (x < 1e-10)', 'facet', 'Rhalf'),
'Rhalf' : 'vertices in x > -1e-10',
'Right' : ('vertices in (x > 4.99)', 'facet'),
}

ebcs = {
'fixb' : ('Left', {'u.all' : 0.0}),
'fixt' : ('Right', {'u.[1,2]' : 0.0}),
}

integrals = {
'i' : 2 * approx_order,
}

##
# Balance of forces.
equations = {
'elasticity' :
"""dw_lin_elastic.i.Omega( solid.D, v, u )
= - dw_surface_ltr.i.Right( load.val, v )""",
}

##
# Solvers etc.
solvers = {
'ls' : ('ls.auto_direct', {}),
'newton' : ('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,
})
}

return locals()