#### Previous topic

standalone/elastic_materials examples

#### Next topic

standalone/homogenized_elasticity examples

# standalone/elastic_materials/compare_elastic_materials.pyΒΆ

Description

Compare various elastic materials w.r.t. uniaxial tension/compression test.

Requires Matplotlib.

source code

```"""
Compare various elastic materials w.r.t. uniaxial tension/compression test.

Requires Matplotlib.
"""
from optparse import OptionParser
import sys
sys.path.append( '.' )

import numpy as nm

def define():
"""Define the problem to solve."""

filename_mesh = 'el3.mesh'

options = {
'nls' : 'newton',
'ls' : 'ls',
'ts' : 'ts',
'save_steps' : -1,
}

functions = {
'linear_tension' : (linear_tension,),
'linear_compression' : (linear_compression,),
'empty' : (lambda ts, coor, mode, region, ig: None,),
}

field_1 = {
'name' : 'displacement',
'dtype' : nm.float64,
'shape' : (3,),
'region' : 'Omega',
'approx_order' : 1,
}

# Coefficients are chosen so that the tangent stiffness is the same for all
# material for zero strains.
# Young modulus = 10 kPa, Poisson's ratio = 0.3
material_1 = {
'name' : 'solid',

'values' : {
'K'  : 8.333, # bulk modulus
'mu_nh' : 3.846, # shear modulus of neoHookean term
'mu_mr' : 1.923, # shear modulus of Mooney-Rivlin term
'kappa' : 1.923, # second modulus of Mooney-Rivlin term
'lam' : 5.769, # Lame coefficients for LE term
'mu_le' : 3.846,
}
}

material_2 = {
'function' : 'empty'
}

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

regions = {
'Omega' : ('all', {}),
'Bottom' : ('nodes in (z < 0.1)', {}),
'Top' : ('nodes in (z > 2.9)', {}),
}

ebcs = {
'fixb' : ('Bottom', {'u.all' : 0.0}),
'fixt' : ('Top', {'u.[0,1]' : 0.0}),
}

##
# Balance of forces.
integral_1 = {
'name' : 'i1',
'kind' : 'v',
'order' : 1,
}
integral_3 = {
'name' : 'isurf',
'kind' : 's',
'order' : 2,
}
equations = {
'linear' : """dw_lin_elastic_iso.i1.Omega( solid.lam, solid.mu_le, v, u )
'neoHookean' : """dw_tl_he_neohook.i1.Omega( solid.mu_nh, v, u )
+ dw_tl_bulk_penalty.i1.Omega( solid.K, v, u )
'Mooney-Rivlin' : """dw_tl_he_neohook.i1.Omega( solid.mu_mr, v, u )
+ dw_tl_he_mooney_rivlin.i1.Omega( solid.kappa, v, u )
+ dw_tl_bulk_penalty.i1.Omega( solid.K, v, u )
}

##
# Solvers etc.
solver_0 = {
'name' : 'ls',
'kind' : 'ls.scipy_direct',
}

solver_1 = {
'name' : 'newton',
'kind' : 'nls.newton',

'i_max'      : 5,
'eps_a'      : 1e-10,
'eps_r'      : 1.0,
'macheps'    : 1e-16,
'lin_red'    : 1e-2, # Linear system error < (eps_a * lin_red).
'ls_red'     : 0.1,
'ls_red_warp': 0.001,
'ls_on'      : 1.1,
'ls_min'     : 1e-5,
'check'      : 0,
'delta'      : 1e-6,
'is_plot'    : False,
'problem'    : 'nonlinear', # 'nonlinear' or 'linear' (ignore i_max)
}

solver_2 = {
'name' : 'ts',
'kind' : 'ts.simple',

't0'    : 0,
't1'    : 1,
'dt'    : None,
'n_step' : 101, # has precedence over dt!
}

return locals()

##
# Pressure tractions.
def linear_tension(ts, coor, mode=None, **kwargs):
if mode == 'qp':
val = nm.tile(0.1 * ts.step, (coor.shape[0], 1, 1))
return {'val' : val}

def linear_compression(ts, coor, mode=None, **kwargs):
if mode == 'qp':
val = nm.tile(-0.1 * ts.step, (coor.shape[0], 1, 1))
return {'val' : val}

def store_top_u( displacements ):
"""Function _store() will be called at the end of each loading step. Top
displacements will be stored into `displacements`."""
def _store( problem, ts, state ):

top = problem.domain.regions['Top']
top_u = problem.get_variables()['u'].get_state_in_region( top )
displacements.append( nm.mean( top_u[:,-1] ) )

return _store

def solve_branch(problem, branch_function):
displacements = {}
for key, eq in problem.conf.equations.iteritems():
problem.set_equations( {key : eq} )

time_solver = problem.get_time_solver()

out = []
time_solver(save_results=False, step_hook=store_top_u(out))
displacements[key] = nm.array(out, dtype=nm.float64)

return displacements

usage = """%prog [options]"""
helps = {
'no_plot' : 'do not show plot window',
}

def main():
from sfepy.base.conf import ProblemConf, get_standard_keywords
from sfepy.fem import ProblemDefinition
from sfepy.base.plotutils import plt

parser = OptionParser(usage=usage, version='%prog')
action="store_true", dest='no_plot',
default=False, help=helps['no_plot'])
options, args = parser.parse_args()

required, other = get_standard_keywords()
# Use this file as the input file.
conf = ProblemConf.from_file( __file__, required, other )

# Create problem instance, but do not set equations.
problem = ProblemDefinition.from_conf( conf,
init_equations = False )

# Solve the problem. Output is ignored, results stored by using the
# step_hook.
u_t = solve_branch(problem, linear_tension)
u_c = solve_branch(problem, linear_compression)

# Get pressure load by calling linear_*() for each time step.
ts = problem.get_timestepper()
for aux in ts.iter_from( 0 )],
dtype=nm.float64).squeeze()
for aux in ts.iter_from( 0 )],
dtype=nm.float64).squeeze()

# Join the branches.
displacements = {}
for key in u_t.keys():
displacements[key] = nm.r_[u_c[key][::-1], u_t[key]]

if plt is None:
print 'matplotlib cannot be imported, printing raw data!'
print displacements
else:
legend = []
for key, val in displacements.iteritems():
legend.append( key )

plt.legend( legend, loc = 2 )
plt.xlabel( 'tension [kPa]' )
plt.ylabel( 'displacement [mm]' )
plt.grid( True )

plt.gcf().savefig( 'pressure_displacement.png' )

if not options.no_plot:
plt.show()

if __name__ == '__main__':
main()
```