standalone/elastic_materials/compare_elastic_materials.pyΒΆ
Description
Compare various elastic materials w.r.t. uniaxial tension/compression test.
Requires Matplotlib.
"""
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 = {
'name' : 'load',
'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 )
= dw_surface_ltr.isurf.Top( load.val, v )""",
'neoHookean' : """dw_tl_he_neohook.i1.Omega( solid.mu_nh, v, u )
+ dw_tl_bulk_penalty.i1.Omega( solid.K, v, u )
= dw_surface_ltr.isurf.Top( load.val, v )""",
'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 )
= dw_surface_ltr.isurf.Top( load.val, v )""",
}
##
# 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} )
load = problem.get_materials()['load']
load.set_function(branch_function)
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')
parser.add_option('-n', '--no-plot',
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()
load_t = nm.array([linear_tension(ts, nm.array([[0.0]]), 'qp')['val']
for aux in ts.iter_from( 0 )],
dtype=nm.float64).squeeze()
load_c = nm.array([linear_compression(ts, nm.array([[0.0]]), 'qp')['val']
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]]
load = nm.r_[load_c[::-1], load_t]
if plt is None:
print 'matplotlib cannot be imported, printing raw data!'
print displacements
print load
else:
legend = []
for key, val in displacements.iteritems():
plt.plot( load, val )
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()