#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Feb  1 08:49:11 2019

@author: bleyerj
"""
from dolfin import *
import mfront_wrapper as mf
import numpy as np
from ufl import diag

width = 1.
height = 1.
mesh = RectangleMesh(Point(0., 0.), Point(width, height), 10, 10)

Ve = VectorElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, MixedElement([Ve, Ve]))
u = Function(V, name="Displacements")
(u1, u2) = split(u)

def bottom(x, on_boundary):
    return near(x[1], 0) and on_boundary
def top(x, on_boundary):
    return near(x[1], height) and on_boundary
def left(x, on_boundary):
    return near(x[0], 0) and on_boundary


bc = [DirichletBC(V.sub(0).sub(0), Constant(0), left),
      DirichletBC(V.sub(1).sub(0), Constant(0), left),
      DirichletBC(V.sub(0).sub(1), Constant(0), bottom),
      DirichletBC(V.sub(1).sub(1), Constant(0), bottom),
      DirichletBC(V.sub(0).sub(1), Constant(-1), top),
      DirichletBC(V.sub(1).sub(1), Constant(-1), top)]

facets = MeshFunction("size_t", mesh, 1)
ds = Measure("ds", subdomain_data=facets)


mat_prop = {"MatrixYoungModulus": 70e3,
            "MatrixPoissonRatio": 0.2,
            "FiberYoungModulus": 200e3,
            "FiberPoissonRatio": 0.3,
            "FiberVolumeFraction": 0.1,
            "Size": 0.1}

material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
                                      "MultiphaseModel",
                                      hypothesis="plane_strain",
                                      material_properties=mat_prop)

problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=0)
problem.bc = bc
problem.register_gradient("MatrixStrain", sym(grad(u1)), symmetric=True)
problem.register_gradient("InclusionStrain", sym(grad(u2)), symmetric=True)
problem.register_gradient("RelativeDisplacement", diag(u2-u1), symmetric=True)

problem.solve(u.vector())

x = np.linspace(width/2, width, 21)
import matplotlib.pyplot as plt
plt.figure()
plt.plot(x, np.array([u(xi, height/2)[0] for xi in x]), "ob", label=r"$u_x$ (matrix, size={})".format(mat_prop["Size"]))
plt.plot(x, np.array([u(xi, height/2)[2] for xi in x]), "or", label=r"$u_x$ (fiber, size={})".format(mat_prop["Size"]))

# mat_prop["Size"] = 0.02

# problem2 = mf.MFrontNonlinearProblem(u, material, quadrature_degree=0)
# problem2.bc = bc
# problem2.register_gradient("MatrixStrain", sym(grad(u1)), symmetric=True)
# problem2.register_gradient("InclusionStrain", sym(grad(u2)), symmetric=True)
# problem2.register_gradient("RelativeDisplacement", diag(u2-u1), symmetric=True)
# u = Function(V)
# problem2.solve(u.vector())

# plt.plot(x, np.array([u(xi, height/2)[0] for xi in x]), "sb", markerfacecolor="None", label=r"$u_x$ (matrix, size={})".format(mat_prop["Size"]))
# plt.plot(x, np.array([u(xi, height/2)[2] for xi in x]), "sr", markerfacecolor="None", label=r"$u_x$ (fiber, size={})".format(mat_prop["Size"]))

plt.legend()
plt.show()