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



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

problem = mf.MFrontNonlinearProblem(T, material, quadrature_degree=0)
problem.bc = bc
# problem.register_gradient("TemperatureGradient", grad(T))

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)))