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

filename_mesh ='silo.mesh'

def get_traction(ts, coors, mode=None, **kwargs):
    #print args
    if mode == 'qp':
        val = [1.0,1.0,1.0] 

        shape = (coors.shape[0], 1, 1)
                        
        out = {
            'val' : np.tile(val, shape),
        }

        return out
	

regions = {
	'Omega' : ('all', {}),
	'Bottom' : ('nodes in (y < 0.01)', {}),
    'Onur' : ('nodes in (y > 0.01)', {}),
}	

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

materials = {
	'solid' : ({
		'lam' : 20.769,
		'mu' : 3.846,
	},),
	'load' : (None,'get_traction'),
}

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

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

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

integrals = {
    'i1' : ('v', 'gauss_o1_d3'),
    'i2' : ('s3', 'gauss_o2_d2'),
}


equations = {
	'elasticity' : 
	"""dw_lin_elastic_iso.i1.Omega( solid.lam, solid.mu, v, u )
	 = - dw_surface_ltr.i1.Onur( load.val, v, u )""",
}

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