#!/usr/bin/env python
"""
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."""
    from sfepy.discrete.fem.meshio import UserMeshIO
    from sfepy.mesh.mesh_generators import gen_block_mesh
    from sfepy.mechanics.matcoefs import stiffness_from_lame

    def mesh_hook(mesh, mode):
        """
        Generate the block mesh.
        """
        if mode == 'read':
            mesh = gen_block_mesh([2, 2, 3], [2, 2, 4], [0, 0, 1.5], name='el3',
                                  verbose=False)
            return mesh

        elif mode == 'write':
            pass

    filename_mesh = UserMeshIO(mesh_hook)

    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,),
    }

    fields = {
        'displacement' : ('real', 3, 'Omega', 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
    materials = {
        'solid' : ({
            '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
            # elasticity for LE term
            'D' : 
stiffness_from_lame(dim=3, lam=5.769, mu=3.846)