# usefull packages
import numpy as np
import getfem as gf
import time

gf.util_trace_level(3)

#########################
# Constants Declaration #
#########################
# Constants for the Steel
E_S = 2e11                                    				# Young's modulus
nu_S = 0.33                                   				# Poisson's ratio
clambda_S = E_S*nu_S/((1+nu_S)*(1-2*nu_S))           		# First Lame coefficient (N/cm^2)
cmu_S = E_S/(2*(1+nu_S))                          			# Second Lame coefficient (N/cm^2)
clambdastar_S = 2*clambda_S*cmu_S/(clambda_S+2*cmu_S);		# Used to compute von Mises stress from displacement
# Constants for the Steel
E_Al = 5.7e10                                    			# Young's modulus
nu_Al = 0.35                                   				# Poisson's ratio
clambda_Al = E_Al*nu_Al/((1+nu_Al)*(1-2*nu_Al))           	# First Lame coefficient (N/cm^2)
cmu_Al = E_Al/(2*(1+nu_Al))                          		# Second Lame coefficient (N/cm^2)
clambdastar_Al = 2*clambda_Al*cmu_Al/(clambda_Al+2*cmu_Al);	# Used to compute von Mises stress from displacement
# Constants for the simulation
poly_order = 2                                          	# Order of polynomial basis function on elements
integ_order = 5                                         	# Order of integration
max_res = 1e-9												# Maximal resolution
max_iter = 30												# Maximal number of iterration to convergence

########################################
# Managing the geometry of the problem #
########################################
# importing the mesh
Mesh_file_name = 'problem_step2.msh'
mesh = gf.Mesh('import','gmsh',Mesh_file_name)
elements_degree = mesh.dim()                                # Degree of the finite element methods
# Recovering specific regions for boundary conditions
# see problem.geo file
BOTTOM_BOUND=10
RIGHT_BOUND=11
TOP_BOUND=12
LEFT_BOUND=13
STEEL_CANT=14
ALUMINIUM_CANT=15

##################################################
# Creating the FEM for the considered quantities #
# and managing the associated methods            #
##################################################
# create a MeshFem objects for the quantities involvedin the problem 
displacement = gf.MeshFem(mesh,elements_degree)  			# displacement
von_Mises = gf.MeshFem(mesh,1)								# von Mises Stress
# assigning basis functions for the quantities on defined elements
method = 'FEM_PK(%d,%d)' %(elements_degree, poly_order)
displacement.set_fem(gf.Fem(method))
von_Mises.set_classical_discontinuous_fem(elements_degree) 
# definition of the integration method
int_meth = 'IM_TRIANGLE(%d)'%integ_order
mim = gf.MeshIm(mesh, gf.Integ(int_meth))

######################
# Creating the model #
######################
model = gf.Model('real')
# adding the physical equations of the linear elasticity to the model
model.add_fem_variable('u', displacement)
# for the Steel
model.add_initialized_data('lambda_S', clambda_S)
model.add_initialized_data('mu_S', cmu_S)
model.add_filtered_fem_variable('lambdastar_S', clambdastar_S)
model.add_isotropic_linearized_elasticity_brick(mim,'u','lambda_S','mu_S',STEEL_CANT)
# for the Aluminium
model.add_initialized_data('lambda_Al', clambda_Al)
model.add_initialized_data('mu_Al', cmu_Al)
model.add_initialized_data('lambdastar_Al', clambdastar_Al)
model.add_isotropic_linearized_elasticity_brick(mim,'u','lambda_Al','mu_Al',ALUMINIUM_CANT)

#######################
# Boundary Conditions #
#######################
# adding a homogeneous Dirichlet condition
model.add_Dirichlet_condition_with_multipliers(mim, 'u', elements_degree-1, LEFT_BOUND)
model.add_Dirichlet_condition_with_multipliers(mim, 'u', elements_degree-1, RIGHT_BOUND)
# adding Neuman add_Dirichlet_condition_with_multipliers
vf = 1      												# Vertical force
F = np.zeros(elements_degree)
F[elements_degree-1] = -vf
model.add_initialized_data('Fdata', F)
model.add_source_term_brick(mim, 'u', 'Fdata', TOP_BOUND)

###########
# Solving #
###########
print 'running solve...'
time_start = time.clock()
model.solve('max_res', max_res, 'max_iter', max_iter, 'noisy')
time_elapsed = (time.clock() - time_start)
print "solve done in {} second" .format(time_elapsed)

# visualisation of results
U = model.variable('u')
#VM_S = model.compute_isotropic_linearized_Von_Mises_or_Tresca('u', 'lambdastar', 'mu', von_Mises)
displacement.export_to_vtk('solution_step2.vtk', displacement, U, 'Displacement')