from sfepy.mechanics.matcoefs import lame_from_youngpoisson

from sfepy.discrete.fem.utils import refine_mesh

from sfepy import data_dir

# Tell SfePy to use our .mesh file

filename_mesh = data_dir + '/meshes/3d/block.mesh'


# Tell SfePy where u want to save the output

output_dir = '.'


# Material parameters.

young = 200000.0 # Young's modulus [MPa]

poisson = 0.26 # Poisson's ratio


options = {

    'output_dir' : output_dir,

}

#Regions are used to define the boundary conditions, the domains of terms
#and materials etc.

regions = {

    'Omega' : 'all',

    'Bottom' : ('vertices in (y<=0.0001)', 'facet'),

    'Top':('vertices in (y>=0.9999)', 'facet'),

}

#Define constitutive parameters (e.g. stiffness, permeability, or
#viscosity),

#Also other non-field arguments of terms (e.g. known traction or volume
#forces).

materials = {

    'Steel' : ({

            'lam' : lame_from_youngpoisson(young, poisson)[0],

            'mu' : lame_from_youngpoisson(young, poisson)[1],

            },),

    'Load' : ({'val' : [[0.0], [-100.0], [0.0]]},),

}


#FE field.Here 'real' is datatype, 3 is dof per node, field is defined
#over omega & '1' is approximation order

fields = {

    'displacement': ('real', 3, 'Omega', 1),

}


#Here we use linear elastic spring equation

equations = {
    'balance_of_forces' :

    """dw_lin_elastic_iso.2.Omega(Steel.lam, Steel.mu, v, u )
    = dw_surface_ltr.0.Top(Load.val, v)""",

}


# Specify the variables that use the FE approximation given by the
#specified field

variables = {

    'u' : ('unknown field', 'displacement', 0),

    'v' : ('test field', 'displacement', 'u'),

}


#Since the bottom is fixed corresponding nodal displacement are zero

ebcs = {

    'Fixed' : ('Bottom', {'u.all' : 0.0}),

}


#In SfePy, a non-linear solver has to be specified even when solving a
#linear problem.

#The linear problem is/should be then solved in one iteration of the
#nonlinear solver

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,
                     
    }),
}


