#!/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

length = 30e-3
width = 5.4e-3
mesh = RectangleMesh(Point(0., 0.), Point(length, width), 100, 10)

V = FunctionSpace(mesh, "CG", 1)
T = Function(V, name="Temperature")

def left(x, on_boundary):
    return near(x[1], 0) and on_boundary
def right(x, on_boundary):
    return near(x[0], length) and on_boundary

Tl = 300
Tr = 800

bc = [DirichletBC(V, Constant(Tl), left),
      DirichletBC(V, Constant(Tr), right)]

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

dt = Constant(1e-3)
theta = Constant(0.5)

material = mf.MFrontNonlinearMaterial("../materials/src/libBehaviour.so",
                                      "HeatTransfer",
                                      hypothesis="plane_strain")
problem = mf.MFrontNonlinearProblem(T, material, quadrature_degree=0)
problem.bc = bc
j = mf.ThermalFlux(T)
j.initialize_functions(mesh, 0)
problem.flux = j
H = mf.Flux(mf.Var(T), name="Enthalpy")
H.initialize_functions(mesh, 0)
problem.state_variables.append(H)

H0 = H.previous()
j0 = j.previous()
problem.L = (problem.u_*(H.function - H0)-dt*(theta*dot(grad(problem.u_), j.function)+
                 -(1-theta)*dot(grad(problem.u_), j0)))*problem.dx
problem.a = (problem.u_*H.tangents[0]*problem.du -
             dt*theta*dot(grad(problem.u_), j.tangents[0]*grad(problem.du) +
                          j.tangents[1]*problem.du))*problem.dx

print(problem.material.data_manager.K.shape)

T.interpolate(Constant(Tl))

problem.solve(T.vector())

x = np.linspace(0, length)
import matplotlib.pyplot as plt
plt.plot(x, np.array([T(xi, width/2) for xi in x]))
#plt.figure()
#plot(project(flux.function[0], FunctionSpace(mesh, "DG",0)))