import getfem as gf
import numpy as np
from numpy import sin, cos, pi
NY = 50
ls_deg  = 1
deg     = 1

# Mesh definition, simplices or square or ....
mesh = gf.Mesh("regular simplices", np.linspace(-1.0,1.0, 2*NY+1),
                                         np.linspace(-0.5,0.5,   NY+1))
# boundary selection
flst  = mesh.outer_faces()
fnor  = mesh.normal_of_faces(flst)
tleft = abs(fnor[1,:]+1) < 1e-14
ttop  = abs(fnor[0,:]-1) < 1e-14
fleft = np.compress(tleft, flst, axis=1)
ftop  = np.compress(ttop, flst, axis=1)

# mark it as boundary
DIRICHLET_BOUNDARY_NUM1 = 1
DIRICHLET_BOUNDARY_NUM2 = 2
NEUMANN_BOUNDARY_NUM = 3
mesh.set_region(DIRICHLET_BOUNDARY_NUM1, fleft)
mesh.set_region(DIRICHLET_BOUNDARY_NUM2, ftop)

# define level set on the defined mesh
ls = gf.LevelSet(mesh, ls_deg)
mf_ls = ls.mf()
ULS = mf_ls.eval('(-0.6+sin(pi*4*x)*cos(pi*4*y))', globals(),locals())
ls.set_values(ULS)

# the mesh with LS to reprensent the cut
mesh_ls = gf.MeshLevelSet(mesh)
mesh_ls.add(ls)
print('Adapting the mesh')
mesh_ls.adapt()

mim_in = gf.MeshIm('levelset', mesh_ls,'inside', gf.Integ('IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6),3)'))
mim_in.set_integ(4)
mim_out = gf.MeshIm('levelset', mesh_ls,'outside', gf.Integ('IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6),3)'))
mim_out.set_integ(4)

mf_basic = gf.MeshFem(mesh)
mf_basic.set_fem(gf.Fem("FEM_PK(%d,%d)"%(2,deg)))
mf = gf.MeshFem('levelset', mesh_ls, mf_basic)
mf.set_qdim(2)

print('Integration methods adapted')

# Model
md = gf.Model("real")
md.add_fem_variable("u", mf)
md.add_initialized_data('Force', [0.0, 1.0])
md.add_initialized_data('lambda1', 1)
md.add_initialized_data('mu1', 1)
md.add_initialized_data('lambda2', 2)
md.add_initialized_data('mu2', 2)

md.add_isotropic_linearized_elasticity_brick(mim_in, "u", "lambda1", "mu1")
md.add_isotropic_linearized_elasticity_brick(mim_out, "u", "lambda2", "mu2")
md.add_Dirichlet_condition_with_penalization(mim_in, 'u', 1.E+9, DIRICHLET_BOUNDARY_NUM1)
md.add_source_term_brick(mim_in, 'u', 'Force', DIRICHLET_BOUNDARY_NUM2)

for i in xrange(10): 
    if(i > 0):
      ULS = ULS + 0.0025*i
      ls.set_values(ULS)
      mesh_ls.adapt()
      mim_in.adapt()
      mim_out.adapt()
      print('Integration methods adapted')
    
    md.solve()
    U   = md.variable('u')
    print("iter %g"%(i+1))