# -*- coding: utf-8 -*-
import getfem as gf
import numpy as np
#clear all;
#% parameters
d = 2;                 #% dimension (cannot be changed for the moment)
E = 1;                 #% Young's modulus
nu = 0.3;              #% Poisson ratio
cmu = E/(2*(1+nu));    #% Lame coefficient

clambda = 2*cmu*nu/(1-nu);   #% Lame coefficient, plane stress

NSize = 0.05; #NX = 30;

m = gf.Mesh('cartesian',np.arange(0,1+NSize,NSize),np.arange(0,1+NSize,NSize)); 

#% create a mesh_fem of for a field of dimension d (i.e. a vector field)
mf = gf.MeshFem(m,d);#mf = gf_mesh_fem(m,d);
#% assign the Q2 fem to all convexes of the mesh_fem,
gf.MeshFem.set(mf, 'fem', gf.Fem('FEM_QK(2,2)'));

#% Integration which will be used
mim = gf.MeshIm(m, gf.Integ('IM_GAUSS_PARALLELEPIPED(2,4)')); 

md=gf.Model('real');#md=gf_model('real');
md.add_fem_variable('u', mf)
md.add_initialized_data('cmu', [cmu]);
md.add_initialized_data('clambda', [clambda]);
md.add_linear_term(mim,'(clambda*Trace(Grad_u)*Id(qdim(u)) + cmu*(Grad_u+Grad_u’)):Grad_Test_u')
#md.add_linear_term('add linear term', mim,'(clambda*Div_Test_u*Id(qdim(u)) + cmu*(Grad_Test_u''+Grad_Test_u)):Grad_Test2_u');
#%md.add_isotropic_linearized_elasticity_brick(mim, 'u', 'clambda', 'cmu');
print 'test ok'