def test_stiffness_tensors(self):
        import numpy as nm

        from sfepy.base.base import assert_
        import sfepy.mechanics.matcoefs as mc

        ok = True

        lam = 1.0
        mu = 4.0

        lam = nm.array([lam] * 3)
        mu = nm.array([mu] * 3)

        d = nm.array([[ 9.,  1.,  1.,  0.,  0.,  0.],
                      [ 1.,  9.,  1.,  0.,  0.,  0.],
                      [ 1.,  1.,  9.,  0.,  0.,  0.],
                      [ 0.,  0.,  0.,  4.,  0.,  0.],
                      [ 0.,  0.,  0.,  0.,  4.,  0.],
                      [ 0.,  0.,  0.,  0.,  0.,  4.]])

        _ds = mc.stiffness_from_lame(3, lam, mu)
        assert_(_ds.shape == (3, 6, 6))

        _ok = True
        for _d in _ds:
            __ok = nm.allclose(_d, d, rtol=0.0, atol=1e-14)
            _ok = _ok and __ok
        self.report('stiffness_from_lame():', _ok)
        ok = ok and _ok

        d = nm.array([[ 6., -2., -2.,  0.,  0.,  0.],
                      [-2.,  6., -2.,  0.,  0.,  0.],
                      [-2., -2.,  6.,  0.,  0.,  0.],
                      [ 0.,  0.,  0.,  4.,  0.,  0.],
                      [ 0.,  0.,  0.,  0.,  4.,  0.],
                      [ 0.,  0.,  0.,  0.,  0.,  4.]])

        _ds = mc.stiffness_from_lame_mixed(3, lam, mu)
        assert_(_ds.shape == (3, 6, 6))

        _ok = True
        for _d in _ds:
            __ok = nm.allclose(_d, d, rtol=0.0, atol=1e-14)
            _ok = _ok and __ok
        self.report('stiffness_from_lame_mixed():', _ok)
        ok = ok and _ok

        blam = 2.0 / 3.0 * (lam - mu)
        _ds = mc.stiffness_from_lame(3, blam, mu)
        assert_(_ds.shape == (3, 6, 6))

        _ok = True
        for _d in _ds:
            __ok = nm.allclose(_d, d, rtol=0.0, atol=1e-14)
            _ok = _ok and __ok
        self.report('stiffness_from_lame() with modified lambda:', _ok)
        ok = ok and _ok

        return ok

def get_inclusion_pars(ts, coor, mode=None, **kwargs):
    """TODO: implement proper 3D -> 2D transformation of constitutive
    matrices."""
    if mode == 'qp':
        n_nod, dim = coor.shape
        sym = (dim + 1) * dim // 2

        dielectric = nm.eye(dim, dtype=nm.float64)
        # !!!
        coupling = nm.ones((dim, sym), dtype=nm.float64)
        #    coupling[0,1] = 0.2

        out = {
            # Lame coefficients in 1e+10 Pa.
            'D' : stiffness_from_lame(dim=2, lam=0.1798, mu=0.148),
            # dielectric tensor
            'dielectric' : dielectric,
            # piezoelectric coupling
            'coupling' : coupling,
            'density' : 0.1142, # in 1e4 kg/m3
        }

        for key, val in six.iteritems(out):
            out[key] = nm.tile(val, (coor.shape[0], 1, 1))
        return out

def get_pars(ts, coors, mode='qp',
             equations=None, term=None, problem=None, **kwargs):
    """
    The material nonlinearity function - the Lamé coefficient `mu`
    depends on the strain.
    """
    if mode != 'qp': return

    val = nm.empty((coors.shape[0], 1, 1), dtype=nm.float64)
    val.fill(1e0)

    order = term.integral.order
    uvar = equations.variables['u']

    strain = problem.evaluate('ev_cauchy_strain.%d.Omega(u)' % order,
                              u=uvar, mode='qp')
    if ts.step > 0:
        strain0 = strains[-1]

    else:
        strain0 = strain

    dstrain = (strain - strain0) / ts.dt
    dstrain.shape = (strain.shape[0] * strain.shape[1], strain.shape[2])

    norm = norm_l2_along_axis(dstrain)

    val += norm[:, None, None]

    # Store history.
    strains[0] = strain

    return {'D': stiffness_from_lame(dim=3, lam=1e1, mu=val), 'mu': val}

    def test_solving(self):
        from sfepy.base.base import IndexedStruct
        from sfepy.discrete import (FieldVariable, Material, Problem, Function,
                                    Equation, Equations, Integral)
        from sfepy.discrete.conditions import Conditions, EssentialBC
        from sfepy.terms import Term
        from sfepy.solvers.ls import ScipyDirect
        from sfepy.solvers.nls import Newton
        from sfepy.mechanics.matcoefs import stiffness_from_lame

        u = FieldVariable('u', 'unknown', self.field)
        v = FieldVariable('v', 'test', self.field, primary_var_name='u')

        m = Material('m', D=stiffness_from_lame(self.dim, 1.0, 1.0))
        f = Material('f', val=[[0.02], [0.01]])

        bc_fun = Function('fix_u_fun', fix_u_fun,
                          extra_args={'extra_arg' : 'hello'})

        fix_u = EssentialBC('fix_u', self.gamma1, {'u.all' : bc_fun})
        shift_u = EssentialBC('shift_u', self.gamma2, {'u.0' : 0.1})

        integral = Integral('i', order=3)

        t1 = Term.new('dw_lin_elastic(m.D, v, u)',
                      integral, self.omega, m=m, v=v, u=u)

        t2 = Term.new('dw_volume_lvf(f.val, v)', integral, self.omega, f=f, v=v)

        eq = Equation('balance', t1 + t2)
        eqs = Equations([eq])

        ls = ScipyDirect({})

        nls_status = IndexedStruct()
        nls = Newton({}, lin_solver=ls, status=nls_status)

        pb = Problem('elasticity', equations=eqs)
        ## pb.save_regions_as_groups('regions')

        pb.set_bcs(ebcs=Conditions([fix_u, shift_u]))
        pb.set_solver(nls)

        state = pb.solve()

        name = op.join(self.options.out_dir, 'test_high_level_solving.vtk')
        pb.save_state(name, state)

        ok = nls_status.condition == 0
        if not ok:
            self.report('solver did not converge!')

        _ok = state.has_ebc()
        if not _ok:
            self.report('EBCs violated!')

        ok = ok and _ok

        return ok

    def get_fargs(self, lam, mu, virtual, state,
                  mode=None, term_mode=None, diff_var=None, **kwargs):
        from sfepy.mechanics.matcoefs import stiffness_from_lame

        mat = stiffness_from_lame(self.region.dim, lam, mu)
        return LinearElasticTerm.get_fargs(self, mat, virtual, state,
                                           mode=mode, term_mode=term_mode,
                                           diff_var=diff_var, **kwargs)

def get_pars(dim, lam, mu):
    c = stiffness_from_lame(3, lam, mu)
    if dim == 2:
        tr = TransformToPlane()
        try:
            c = tr.tensor_plane_stress(c3=c)
        except:
            sym = (dim + 1) * dim / 2
            c = nm.zeros((sym, sym), dtype=nm.float64)

    return c

def get_pars(lam, mu, dim, full=False):
    from sfepy.mechanics.matcoefs import stiffness_from_lame, TransformToPlane

    if full:
        c = stiffness_from_lame(3, lam, mu)
        if dim == 2:
            tr = TransformToPlane()
            try:
                c = tr.tensor_plane_stress(c3=c)
            except:
                sym = (dim + 1) * dim / 2
                c = nm.zeros((sym, sym), dtype=nm.float64)
        return c
    else:
        return lam, mu

        def _material_func_(ts, coors, mode=None, **kwargs):
            if mode == 'qp':
                ijk_out = np.empty_like(coors, dtype=int)
                ijk = np.floor((coors - min_xyz[None]) / self.dx,
                               ijk_out, casting="unsafe")
                ijk_tuple = tuple(ijk.swapaxes(0, 1))
                property_array_qp = property_array[ijk_tuple]
                lam = property_array_qp[..., 0]
                mu = property_array_qp[..., 1]
                lam = np.ascontiguousarray(lam.reshape((lam.shape[0], 1, 1)))
                mu = np.ascontiguousarray(mu.reshape((mu.shape[0], 1, 1)))

                from sfepy.mechanics.matcoefs import stiffness_from_lame
                stiffness = stiffness_from_lame(dims, lam=lam, mu=mu)
                return {'lam': lam, 'mu': mu, 'D': stiffness}
            else:
                return

def get_pars(ts, coor, mode, output_dir=".", **kwargs):
    if mode == "qp":
        n_nod, dim = coor.shape
        sym = (dim + 1) * dim / 2

        out = {}
        out["D"] = nm.tile(stiffness_from_lame(dim, lam=1.7, mu=0.3), (coor.shape[0], 1, 1))

        aa = nm.zeros((sym, 1), dtype=nm.float64)
        aa[:dim] = 0.132
        aa[dim:sym] = 0.092
        out["alpha"] = nm.tile(aa, (coor.shape[0], 1, 1))

        perm = nm.eye(dim, dtype=nm.float64)
        out["K"] = nm.tile(perm, (coor.shape[0], 1, 1))

        return out

def stress_strain(out, pb, state,lam,mu, extend=False):
    """
    Calculate and output strain and stress for given displacements.
    """
    from sfepy.base.base import Struct
    D = stiffness_from_lame(3,lam, mu)    
    #asphalt = Material('Asphalt', D=D)
    mat = Material('Mat', D=D)
    ev = pb.evaluate
    strain = ev('ev_cauchy_strain.2.Omega(u)', mode='el_avg')
    stress = ev('ev_cauchy_stress.2.Omega(Mat.D, u)', mode='el_avg',
                copy_materials=False)

    out['cauchy_strain'] = Struct(name='output_data', mode='cell',
                                  data=strain, dofs=None)
    out['cauchy_stress'] = Struct(name='output_data', mode='cell',
                                  data=stress, dofs=None)

    return out

def get_pars(ts, coor, mode, **kwargs):
    """
    Define the material parameters.
    """
    if mode == 'qp':
        n_nod, dim = coor.shape

        out = {}
        out['D'] = nm.tile(stiffness_from_lame(dim, lam=1.7, mu=0.3),
                           (coor.shape[0], 1, 1))

        alpha = [[0.132, 0.092],
                 [0.052, 0.132]]
        out['alpha'] = nm.tile(alpha, (coor.shape[0], 1, 1))

        perm = nm.eye(dim, dtype=nm.float64)
        out['K'] = nm.tile(perm, (coor.shape[0], 1, 1))

        out['np_eps'] = nm.tile(1e5, (coor.shape[0], 1, 1))

        return out

from sfepy.postprocess.probes_vtk import Probe
from sfepy.mechanics.matcoefs import stiffness_from_lame

# Define options.
options = {
    'output_dir' : '.',
    'post_process_hook' : 'post_process',
}

# Update materials, as ev_cauchy_stress below needs the elastic constants in
# the tensor form.
materials = deepcopy(materials)
solid = materials['solid'][0]
lam, mu = solid['lam'], solid['mu']
solid.update({
    'D' : stiffness_from_lame(3, lam=lam, mu=mu),
})

# The function returning the probe parametrization.
def par_fun(idx):
    return nm.log(idx + 1) / nm.log(20) * 0.04

# Define the function post_process, that will be called after the problem is
# solved.
def post_process(out, problem, state, extend=False):
    """
    This will be called after the problem is solved.

    Parameters
    ----------
    out : dict

import numpy as nm
from sfepy.mechanics.matcoefs import stiffness_from_lame
from sfepy import data_dir

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

regions = {
    'Omega' : 'all',
    'Left' : ('vertices in (x < -4.99)', 'facet'),
    'Omega1' : 'vertices in (x < 0.001)',
    'Omega2' : 'vertices in (x > -0.001)',
}

materials = {
    'solid' : ({
        'D' : stiffness_from_lame(3, lam=1e2, mu=1e1),
        'prestress' : 0.1 * nm.array([[1.0], [1.0], [1.0],
                                      [0.5], [0.5], [0.5]],
                                     dtype=nm.float64),
        'DF' : stiffness_from_lame(3, lam=8e0, mu=8e-1),
        'nu' : nm.array([[-0.5], [0.0], [0.5]], dtype=nm.float64),
    },),
}

fields = {
    'displacement': ('real', 'vector', 'Omega', 1),
}

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

t0 = 0.0
t1 = 20.0
n_step = 21

# Fading memory times.
f_t0 = 0.0
f_t1 = 5.0
f_n_step = 6

decay = 0.8
mode = 'eth'

## Configure above. ##

times = nm.linspace(f_t0, f_t1, f_n_step)
kernel = get_exp_fading_kernel(stiffness_from_lame(3, lam=1.0, mu=1.0),
                               decay, times)

dt = (t1 - t0) / (n_step - 1)
fading_memory_length = min(int((f_t1 - f_t0) / dt) + 1, n_step)
output('fading memory length:', fading_memory_length)

def post_process(out, pb, state, extend=False):
    """
    Calculate and output strain and stress for given displacements.
    """
    from sfepy.base.base import Struct

    ev = pb.evaluate
    strain = ev('ev_cauchy_strain.2.Omega(u)', mode='el_avg')
    out['cauchy_strain'] = Struct(name='output_data', mode='cell',

region_10 = {
    'name' : 'Bottom',
    'select' : 'vertices in (y < %f)' % -wy,
    'kind' : 'facet',
}

region_11 = {
    'name' : 'Top',
    'select' : 'vertices in (y > %f)' % wy,
    'kind' : 'facet',
}

material_1 = {
    'name' : 'solid',
    'values' : {
        'D' : stiffness_from_lame(2, 1e1, 1e0),
        'density' : 1e-1,
    },
}

field_1 = {
    'name' : '2_displacement',
    'dtype' : 'real',
    'shape' : 'vector',
    'region' : 'Y',
    'approx_order' : 2,
}

variable_1 = {
    'name' : 'u',
    'kind' : 'unknown field',

def main():
    from sfepy import data_dir

    parser = ArgumentParser()
    parser.add_argument('--version', action='version', version='%(prog)s')
    parser.add_argument('-s', '--show',
                        action="store_true", dest='show',
                        default=False, help=helps['show'])
    options = parser.parse_args()

    mesh = Mesh.from_file(data_dir + '/meshes/2d/rectangle_tri.mesh')
    domain = FEDomain('domain', mesh)

    min_x, max_x = domain.get_mesh_bounding_box()[:,0]
    eps = 1e-8 * (max_x - min_x)
    omega = domain.create_region('Omega', 'all')
    gamma1 = domain.create_region('Gamma1',
                                  'vertices in x < %.10f' % (min_x + eps),
                                  'facet')
    gamma2 = domain.create_region('Gamma2',
                                  'vertices in x > %.10f' % (max_x - eps),
                                  'facet')

    field = Field.from_args('fu', nm.float64, 'vector', omega,
                            approx_order=2)

    u = FieldVariable('u', 'unknown', field)
    v = FieldVariable('v', 'test', field, primary_var_name='u')

    m = Material('m', D=stiffness_from_lame(dim=2, lam=1.0, mu=1.0))
    f = Material('f', val=[[0.02], [0.01]])

    integral = Integral('i', order=3)

    t1 = Term.new('dw_lin_elastic(m.D, v, u)',
                  integral, omega, m=m, v=v, u=u)
    t2 = Term.new('dw_volume_lvf(f.val, v)', integral, omega, f=f, v=v)
    eq = Equation('balance', t1 + t2)
    eqs = Equations([eq])

    fix_u = EssentialBC('fix_u', gamma1, {'u.all' : 0.0})

    bc_fun = Function('shift_u_fun', shift_u_fun,
                      extra_args={'shift' : 0.01})
    shift_u = EssentialBC('shift_u', gamma2, {'u.0' : bc_fun})

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

    pb = Problem('elasticity', equations=eqs, nls=nls, ls=ls)
    pb.save_regions_as_groups('regions')

    pb.time_update(ebcs=Conditions([fix_u, shift_u]))

    vec = pb.solve()
    print(nls_status)

    pb.save_state('linear_elasticity.vtk', vec)

    if options.show:
        view = Viewer('linear_elasticity.vtk')
        view(vector_mode='warp_norm', rel_scaling=2,
             is_scalar_bar=True, is_wireframe=True)

def create_local_problem(omega_gi, orders):
    """
    Local problem definition using a domain corresponding to the global region
    `omega_gi`.
    """
    order_u, order_p = orders

    mesh = omega_gi.domain.mesh

    # All tasks have the whole mesh.
    bbox = mesh.get_bounding_box()
    min_x, max_x = bbox[:, 0]
    eps_x = 1e-8 * (max_x - min_x)

    min_y, max_y = bbox[:, 1]
    eps_y = 1e-8 * (max_y - min_y)

    mesh_i = Mesh.from_region(omega_gi, mesh, localize=True)
    domain_i = FEDomain('domain_i', mesh_i)
    omega_i = domain_i.create_region('Omega', 'all')

    gamma1_i = domain_i.create_region('Gamma1',
                                      'vertices in (x < %.10f)'
                                      % (min_x + eps_x),
                                      'facet', allow_empty=True)
    gamma2_i = domain_i.create_region('Gamma2',
                                      'vertices in (x > %.10f)'
                                      % (max_x - eps_x),
                                      'facet', allow_empty=True)
    gamma3_i = domain_i.create_region('Gamma3',
                                      'vertices in (y < %.10f)'
                                      % (min_y + eps_y),
                                      'facet', allow_empty=True)

    field1_i = Field.from_args('fu', nm.float64, mesh.dim, omega_i,
                               approx_order=order_u)

    field2_i = Field.from_args('fp', nm.float64, 1, omega_i,
                               approx_order=order_p)

    output('field 1: number of local DOFs:', field1_i.n_nod)
    output('field 2: number of local DOFs:', field2_i.n_nod)

    u_i = FieldVariable('u_i', 'unknown', field1_i, order=0)
    v_i = FieldVariable('v_i', 'test', field1_i, primary_var_name='u_i')
    p_i = FieldVariable('p_i', 'unknown', field2_i, order=1)
    q_i = FieldVariable('q_i', 'test', field2_i, primary_var_name='p_i')

    if mesh.dim == 2:
        alpha = 1e2 * nm.array([[0.132], [0.132], [0.092]])

    else:
        alpha = 1e2 * nm.array([[0.132], [0.132], [0.132],
                                [0.092], [0.092], [0.092]])

    mat = Material('m', D=stiffness_from_lame(mesh.dim, lam=10, mu=5),
                   k=1, alpha=alpha)
    integral = Integral('i', order=2*(max(order_u, order_p)))

    t11 = Term.new('dw_lin_elastic(m.D, v_i, u_i)',
                   integral, omega_i, m=mat, v_i=v_i, u_i=u_i)
    t12 = Term.new('dw_biot(m.alpha, v_i, p_i)',
                   integral, omega_i, m=mat, v_i=v_i, p_i=p_i)
    t21 = Term.new('dw_biot(m.alpha, u_i, q_i)',
                   integral, omega_i, m=mat, u_i=u_i, q_i=q_i)
    t22 = Term.new('dw_laplace(m.k, q_i, p_i)',
                   integral, omega_i, m=mat, q_i=q_i, p_i=p_i)

    eq1 = Equation('eq1', t11 - t12)
    eq2 = Equation('eq1', t21 + t22)
    eqs = Equations([eq1, eq2])

    ebc1 = EssentialBC('ebc1', gamma1_i, {'u_i.all' : 0.0})
    ebc2 = EssentialBC('ebc2', gamma2_i, {'u_i.0' : 0.05})
    def bc_fun(ts, coors, **kwargs):
        val = 0.3 * nm.sin(4 * nm.pi * (coors[:, 0] - min_x) / (max_x - min_x))
        return val

    fun = Function('bc_fun', bc_fun)
    ebc3 = EssentialBC('ebc3', gamma3_i, {'p_i.all' : fun})

    pb = Problem('problem_i', equations=eqs, active_only=False)
    pb.time_update(ebcs=Conditions([ebc1, ebc2, ebc3]))
    pb.update_materials()

    return pb

    'u'       : ('unknown field',   'displacement', 0),
    'v'       : ('test field',      'displacement', 'u'),
    'p'       : ('unknown field',   'pressure', 1),
    'q'       : ('test field',      'pressure', 'p'),
}

ebcs = {
    'fix_u' : ('Bottom', {'u.all' : 0.0}),
    'load_u' : ('Top', {'u.2' : 0.2}),
    'load_p' : ('Left', {'p.all' : 1.0}),
}

material_1 = {
    'name' : 'm',
    'values' : {
        'D': stiffness_from_lame(dim=3, lam=1.7, mu=0.3),
        'alpha' : nm.array( [[0.132], [0.132], [0.132],
                             [0.092], [0.092], [0.092]],
                            dtype = nm.float64 ),
        'K' : nm.array( [[2.0, 0.2, 0.0], [0.2, 1.0, 0.0], [0.0, 0.0, 0.5]],
                        dtype = nm.float64 ),
    }
}

integral_1 = {
    'name' : 'i1',
    'order' : 1,
}

integral_2 = {
    'name' : 'i2',

    return val.T

functions = {
    'get_shift' : (get_shift,),
    'linear_tension' : (linear_tension,),
    'match_y_plane' : (per.match_y_plane,),
    'match_z_plane' : (per.match_z_plane,),
}

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

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

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

regions = {
    'Omega' : 'all',
    'Left' : ('vertices in (x < -4.99)', 'facet'),
    'Right' : ('vertices in (x > 4.99)', 'facet'),
    'Bottom' : ('vertices in (z < -0.99)', 'facet'),
    'Top' : ('vertices in (z > 0.99)', 'facet'),

                              % (vn, vn), verbose=False)
            print 'dw_lin_elastic', vn, val

        except:
            pass

    return out

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

materials = {
    'm' : ({'D' : stiffness_from_lame(3, lam=0.0007, mu=0.0003),
            'one' : 1.0},),
}

regions = {
    'Omega' : ('all', {}),
    'Gamma_Left' : ('nodes in (x < %f)' % xmin, {}),
    'Gamma_Right' : ('nodes in (x > %f)' % xmax, {}),
}

fields = {
    'displacements' : ('real', 'vector', 'Omega', 1),
}

variables = {
    'u' : ('unknown field', 'displacements', 0),