# -*- coding: utf-8 -*-
"""
Spyder Editor

This is a temporary script file.
"""
import fenics as fe
import matplotlib.pyplot as plt
import numpy as np
import dolfin
import os

# --------------------
# Functions and classes
# --------------------
#def bottom(x, on_boundary):
#   return (on_boundary and fe.near(x[1], 0.0))

# Strain function
def epsilon(u):
    return fe.sym(fe.grad(u))

# Stress function
def sigma(u):
    return lambda_*fe.div(u)*fe.Identity(2) + 2*mu*epsilon(u)


# --------------------
# Parameters
# --------------------

# Density
rho = fe.Constant(200.0)

# Young's modulus and Poisson's ratio
E = 135e9
nu = 0.35

# Lame's constants
lambda_ = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

l_x, l_y = 5.0, 5.0  # Domain dimensions
n_x, n_y = 20, 20  # Number of elements
lc = l_x/n_x
r = 2 # rayon du trou 

# Load
f = 600
F = fe.Constant((0.0, f))

# Model type
model = "plane_strain"
if model == "plane_stress":
    lambda_ = 2*mu*lambda_/(lambda_+2*mu)
#%% --------------------
# Geometry
# --------------------
#on adapte la géometrie aux paramètres
f1 = open("quartvierge.geo",'r')
toto = f1.read()
f1.close
toto = toto.replace("lx =","lx = %f;" % l_x)
toto = toto.replace("ly =","ly = %f;" % l_y)
toto = toto.replace("lc =","lc = %f;" % lc)
toto = toto.replace("r =","r = %f;" % r)
f2 = open("quart.geo",'w')
f2.writelines(toto)
f2.close()
# Commandes pour générer le maillage en appelant GMSH
# 
command = 'gmsh -2 quart.geo -format msh2'
os.system(command)
#print ("Le maillage est fait!\n")
command = 'dolfin-convert quart.msh geometry.xml'
os.system(command)
#print ("Finished meshing!\n")
# demonstration du maillage
mesh = fe.Mesh("geometry.xml")
plt.figure(1)
fe.plot(mesh)
plt.show()
#%%

# Definition des limites 
left = dolfin.CompiledSubDomain("on_boundary && near(x[0],0)")
right = dolfin.CompiledSubDomain("on_boundary && near(x[0], Lx)", Lx = l_x)
top = dolfin.CompiledSubDomain("on_boundary && near(x[1], Ly)", Ly = l_y)
bottom = dolfin.CompiledSubDomain("on_boundary && near(x[1], 0)")

boundary_indices = {"left": 0, "right": 1, "top": 2, "bottom": 3}
boundary_markers = dolfin.MeshFunction("size_t", mesh, dim=1, value=0)
left.mark(boundary_markers, boundary_indices["left"])
right.mark(boundary_markers, boundary_indices["right"])
right.mark(boundary_markers, boundary_indices["right"])
top.mark(boundary_markers, boundary_indices["top"])
bottom.mark(boundary_markers, boundary_indices["bottom"])

ds = dolfin.ds(domain=mesh, subdomain_data=boundary_markers)
dx = dolfin.dx(domain=mesh)

# --------------------
# Function spaces
# --------------------
V = fe.VectorFunctionSpace(mesh, "CG", 1)
u_tr = fe.TrialFunction(V)
u_test = fe.TestFunction(V)

# --------------------
# Boundary conditions
# --------------------
bc_l = fe.DirichletBC(V.sub(1), 0.0, left)

bc_b = fe.DirichletBC(V.sub(0), 0.0, bottom)

#bc = fe.DirichletBC(V, fe.Constant((0.0, 0.0)), bottom)

# --------------------
# Weak form
# --------------------
a = fe.inner(sigma(u_tr), epsilon(u_test))*fe.dx
l = fe.inner(F, u_test)*ds(boundary_indices["top"])

# --------------------
# Solver
# --------------------
u = fe.Function(V)
A_ass, L_ass = fe.assemble_system(a, l)

fe.solve(A_ass, u.vector(), L_ass)

print("jusqu'ici tout va bien")

#%% --------------------
# Post-process
# --------------------

plt.figure(2)
fe.plot(u.sub(0), title="Ux")
plt.show()
plt.figure(3)
fe.plot(u.sub(1), title="Uy")
plt.show()
plt.figure(4)
dolfin.plot(u, mode="displacement", title="displacement")
plt.show()
plt.figure(5)
stress = sigma(u)
fe.plot(stress[1, 1], title="Stress")
plt.show()