from sfepy.mechanics.matcoefs import ElasticConstants, stiffness_tensor_youngpoisson
filename_mesh = 'cyl.inp'

output_dir = '.'

young = 2000.0
poisson = 0.49

ec = ElasticConstants(young=young,poisson=poisson)


def zdisp(ts, coor, **kwargs):
    ret = -ts.dt*ts.step
    print '='*51
    print 'Z-Displacement: %s mm' % str(ret)
    return ret

def stress_strain(out, pb, state, extend=False):
    from sfepy.base.base import Struct

    ev = pb.evaluate
    strain = ev('de_cauchy_strain.2.Omega(u)', mode='el_avg')
    stress = ev('de_cauchy_stress.2.Omega(Asphalt.D, u)', mode='el_avg')

    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
    
options = {
    'ts' : 'ts',
    'save_steps' : -1,
    'output_dir' : output_dir,
    #'post_process_hook' : 'stress_strain',
}

regions = {
    'Omega' : ('all', {}),
    'X-Face' : ('nodes in (x < 0.001)', {}),
    'Y-Face' : ('nodes in (y < 0.001)', {}),
    'Z-Face' : ('nodes in (z < 0.001)', {}),
    'Top' : ('nodes in (z > 74.999)', {}),
}

materials = {
    'Asphalt' : ({
        'lam' : ec.lam,
        'mu' : ec.mu,
        'D' : stiffness_tensor_youngpoisson(3, young, poisson),
    },),
}

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

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

equations = {
   'balance_of_forces' :
   #"""dw_lin_elastic.2.Omega(Asphalt.D,v,u) + dw_lin_elastic_th.2.Omega(ts,Asphalt.D,v,du/dt) = 0"""
   """dw_lin_elastic_th.2.Omega(ts,Asphalt.D,v,du/dt) = 0"""
}

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

ebcs = {
    'X-Sym' : ('X-Face', {'u.0' : 0.0}),
    'Y-Sym' : ('Y-Face', {'u.1' : 0.0}),
    'Z-Sym' : ('Z-Face', {'u.2' : 0.0}),
    'U-22' : ('Top', {'u.2' : 'zdisp'}),
}

solvers = {
    'ls' : ('ls.scipy_direct', {'presolve' : True}),
    'newton' : ('nls.newton', {
        'i_max' : 1,
        'eps_a' : 1e-6,
        'problem' : 'nonlinear'
    }),
    'ts' : ('ts.simple', {
        't0' : 0.0,
        't1' : 1.0,
        'dt' : None,
        'n_step' : 5,
    }),
}