#! Linear Elasticity
#! =================
#$ \centerline{Example input file, \today}

#! This file models a cylinder that is fixed at one end while the
#! second end has a specified displacement of 0.01 in the x direction
#! (this boundary condition is named Displaced). There is also a specified
#! displacement of 0.005 in the z direction for points in
#! the region labeled SomewhereTop. This boundary condition is named
#! PerturbedSurface.  The region SomewhereTop is specified as those nodes for
#! which
#!              (z > 0.017) & (x > 0.03) & (x <  0.07).
#! The output is the displacement for each node, saved by default to
#! simple_out.vtk. The material is linear elastic and its properties are
#! specified as Lame parameters (see
#! http://en.wikipedia.org/wiki/Lam%C3%A9_parameters)
#!

#! Mesh
#! ----
from sfepy import data_dir

filename_mesh = data_dir + '/meshes/3d/cylinder.mesh'
#! Regions
#! -------
#! Whole domain 'Omega', left and right ends.
regions = {
    'Omega' : ('all', {}),
    'Left' : ('nodes in (x < 0.001)', {}),
    'Right' : ('nodes in (x > 0.099)', {}),
    'SomewhereTop' : ('nodes in (z > 0.017) & (x > 0.03) & (x < 0.07)', {}),
}
#! Materials
#! ---------
#! The linear elastic material model is used. Properties are
#! specified as Lame parameters.
materials = {
    'solid' : ({'lam' : 1e1, 'mu' : 1e0},),
}
#! Fields
#! ------
#! A field is used to define the approximation on a (sub)domain
#! A displacement field (three DOFs/node) will be computed on a region
#! called 'Omega' using P1 (four-node tetrahedral) finite elements.
fields = {
    'displacement': ('real', 'vector', 'Omega', 1),
}
#! Integrals
#! ---------
#! Define the integral type Volume/Surface and quadrature rule
#! (here: dim=3, order=1).
integrals = {
    'i1' : ('v', 'gauss_o1_d3'),
}
#! Variables
#! ---------
#! One field is used for unknown variable (generate discrete degrees
#! of freedom) and the seccond field for the corresponding test variable of
#! the weak formulation.
variables = {
    'u' : ('unknown field', 'displacement', 0),
    'v' : ('test field', 'displacement', 'u'),
}
#! Boundary Conditions
#! -------------------
#! The left end of the cilinder is fixed (all DOFs are zero) and
#! the 'right' end has non-zero displacements only in the x direction.
ebcs = {
    'Fixed' : ('Left', {'u.all' : 0.0}),
    'Displaced' : ('Right', {'u.0' : 0.01, 'u.[1,2]' : 0.0}),
    'PerturbedSurface' : ('SomewhereTop', {'u.2' : 0.005}),
}
#! Equations
#! ---------
#! The weak formulation of the linear elastic problem.
equations = {
    'balance_of_forces' :
    """dw_lin_elastic_iso.i1.Omega( solid.lam, solid.mu, v, u ) = 0""",
}
#! Solvers
#! -------
#! Define linear and nonlinear solver.
#! Even linear problems are solved by a nonlinear solver (KISS rule) - only one
#! iteration is needed and the final rezidual is obtained for free.
solvers = {
    'ls' : ('ls.scipy_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,
                  'is_plot'    : False,
                  # 'nonlinear' or 'linear' (ignore i_max)
                  'problem'   : 'nonlinear'}),
}

"""
Linear elasticity with pressure traction load on a surface and
constrained to one-dimensional motion.
"""
import numpy as nm

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

        return {'val' : val}

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

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

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

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

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

    materials = {
        'solid' : ({
            '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' : ('nodes in (x < -4.99)', {}),
        'Right' : ('nodes in (x > 4.99)', {}),
    }

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

    ##
    # Balance of forces.
    equations = {
        'elasticity' :
        """dw_lin_elastic_iso.2.Omega( solid.lam, solid.mu, v, u )
         = - dw_surface_ltr.2.Right( load.val, v )""",
    }

    ##
    # Solvers etc.
    solvers = {
        'ls' : ('ls.scipy_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,
                      'is_plot'    : False,
                      # 'nonlinear' or 'linear' (ignore i_max)
                      'problem'   : 'nonlinear'}),
    }

    ##
    # FE assembling parameters.
    fe = {
        'chunk_size' : 1000,
        'cache_override' : False,
    }

    return locals()

# 30.04.2009
#!
#! Linear Elasticity
#! =================
#$ \centerline{Example input file, \today}

#! This file models a cylinder that is fixed at one end while the
#! second end has a specified displacement of 0.02 in the x direction
#! (this boundary condition is named PerturbedSurface).
#! The output is the displacement for each node, saved by default to
#! simple_out.vtk. The material is linear elastic.
from sfepy import data_dir

from sfepy.mechanics.matcoefs import stiffness_tensor_youngpoisson_mixed, bulk_modulus_youngpoisson

#! Mesh
#! ----

dim = 3
approx_u = '3_4_P1'
approx_p = '3_4_P0'
order = 2
filename_mesh = data_dir + '/meshes/3d/cylinder.mesh'
#! Regions
#! -------
#! Whole domain 'Omega', left and right ends.
regions = {
    'Omega' : ('all', {}),
    'Left' : ('nodes in (x < 0.001)', {}),
    'Right' : ('nodes in (x > 0.099)', {}),
}
#! Materials
#! ---------
#! The linear elastic material model is used.
materials = {
    'solid' : ({'D' : stiffness_tensor_youngpoisson_mixed( dim, 0.7e9, 0.4 ),
                'gamma' : 1.0/bulk_modulus_youngpoisson( 0.7e9, 0.4 )},),
}
#! Fields
#! ------
#! A field is used to define the approximation on a (sub)domain
fields = {
    'displacement': ('real', 'vector', 'Omega', 1),
    'pressure' : ('real', 'scalar', 'Omega', 0),
}
#! Integrals
#! ---------
#! Define the integral type Volume/Surface and quadrature rule.
integrals = {
    'i1' : ('v', 'gauss_o%d_d%d' % (order, dim) ),
}
#! Variables
#! ---------
#! Define displacement and pressure fields and corresponding fields
#! for test variables.
variables = {
    'u' : ('unknown field', 'displacement'),
    'v' : ('test field', 'displacement', 'u'),
    'p' : ('unknown field', 'pressure'),
    'q' : ('test field', 'pressure', 'p'),
}
#! Boundary Conditions
#! -------------------
#! The left end of the cylinder is fixed (all DOFs are zero) and
#! the 'right' end has non-zero displacements only in the x direction.
ebcs = {
    'Fixed' : ('Left', {'u.all' : 0.0}),
    'PerturbedSurface' : ('Right', {'u.0' : 0.02, 'u.1' : 0.0}),
}
#! Equations
#! ---------
#! The weak formulation of the linear elastic problem.
equations = {
    'balance_of_forces' :
    """  dw_lin_elastic.i1.Omega( solid.D, v, u )
       - dw_stokes.i1.Omega( v, p )
       = 0 """,
    'pressure constraint' :
    """- dw_stokes.i1.Omega( u, q )
       - dw_mass_scalar_w.i1.Omega( solid.gamma, q, p )
       = 0""",
}
#! Solvers
#! -------
#! Define linear and nonlinear solver.
#! Even linear problems are solved by a nonlinear solver - only one
#! iteration is needed and the final rezidual is obtained for free.
solvers = {
    'ls' : ('ls.scipy_direct', {}),
    'newton' : ('nls.newton',
                { 'i_max'      : 1,
                  'eps_a'      : 1e-2,
                  'eps_r'      : 1e-10,
                  'problem'   : 'nonlinear'}),
}
#! FE assembling parameters
#! ------------------------
#! 'chunk_size' determines maximum number of elements to assemble in one C
#! function call. Higher values mean faster assembling, but also more memory
#! usage.
fe = {
    'chunk_size' : 1000
}
#! Options
#! -------
#! Various problem-specific options.
options = {
    'output_dir' : './output',
    'absolute_mesh_path' : True,
}