#!/usr/bin/env python
# coding: utf-8

# # Stationnary non-linear heat transfer
# 
# ## Description of the non-linear constitutive heat transfer law
# 
# ## FEniCS implementation
# 
# We consider a rectanglar domain with imposed temperatures `Tl` (resp. `Tr`) on the left (resp. right boundaries). We want to solve for the temperature field `T` inside the domain using a $P^1$-interpolation. We initialize the temperature at value `Tl` throughout the domain. We finally load the material library with a `plane_strain` hypothesis.

# In[1]:


import matplotlib.pyplot as plt
from dolfin import *
import mfront_wrapper as mf

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[0], 0) and on_boundary
def right(x, on_boundary):
    return near(x[0], length) and on_boundary

Tl = 300
Tr = 800
T.interpolate(Constant(Tl))

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

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


# The MFront behaviour implicitly declares the temperature as an external state variable called `"Temperature"`. We must therefore associate this external state variable to a known mechanical field. This can be achieved explicitly using the `register_external_state_variable` method. In the present case, this can be donc automatically since the name of the unknown temperature field matches the [TFEL Glossary](http://tfel.sourceforge.net/glossary.html) name `"Temperature"`. In this case, the following message is printed:
# ```
# Automatic registration of 'Temperature' as an external state variable.
# ```
# For problems in which the temperature only acts as a parameter (no jacobian blocks with respect to the temperature), the temperature can be automatically registered as a constant value ($293.15 \text{ K}$ by default) or to any other (`dolfin.Constant` or `float`) value using the `register_external_state_variable` method.
# 
# In the FEniCS interface, we instantiate the main mechanical unknown, here the temperature field `T` which has to be named `"Temperature"` in order to match MFront's predefined name. Using another name than this will later result in an error saying:
# ```
# ValueError: 'Temperature' could not be associated with a registered gradient or a known state variable.
# ```
# 
# The MFront behaviour declares the field `"TemperatureGradient"` as a Gradient variable, with its associated Flux called `"HeatFlux"`. We can check that the `material` object retrieves MFront's gradient and flux names, as well as the different tangent operator blocks which have been defined, namely `dj_ddgT` and `dj_ddT` in the present case:

# In[2]:


print(material.get_gradient_names())
print(material.get_flux_names())
print(["d{}_d{}".format(*t) for t in material.get_tangent_block_names()])


# When defining the non-linear problem, we will specify the boundary conditions and the requested quadrature degree which will control the number of quadrature points used in each cell to compute the non-linear constitutive law. Here, we specify a quadrature of degree 2 (i.e. 3 Gauss points for a triangular element). Finally, we need to associate to MFront gradient object the corresponding UFL expression as a function of the unknown field `T`. To do so, we use the `register_gradient` method linking MFront `"TemperatureGradient"` object to the UFL expression `grad(T)`. Doing so, the corresponding non-linear variational problem will be automatically be built:
# 
# \begin{equation}
# F(\widehat{T}) = \int_\Omega \boldsymbol{j}\cdot \nabla \widehat{T} \text{dx} = 0 \quad \forall \widehat{T}
# \end{equation}

# In[3]:


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


# From the two tangent operator blocks `dj_ddgT` and `dj_ddT`, it will automatically be deduced that the heat flux $\boldsymbol{j}$ is a function of both the temperature gradient $\boldsymbol{g}=\nabla T$ and the temperature itself i.e. $\boldsymbol{j}=\boldsymbol{j}(\boldsymbol{g}, T)$. The following tangent bilinear form will therefore be used when solving the above non-linear problem:
# 
# \begin{equation}
# J(\widehat{T},T^*) = \int_{\Omega} \nabla \widehat{T}\cdot\left(\dfrac{\partial \boldsymbol{j}}{\partial \boldsymbol{g}}\cdot \nabla T^*+\dfrac{\partial \boldsymbol{j}}{\partial T}\cdot T^*\right) \text{dx}
# \end{equation}
# 
# Similarly to the case of external state variables, common gradient expressions for some [TFEL Glossary](http://tfel.sourceforge.net/glossary.html) names have been already predefined which avoid calling explicitly the `register_gradient` method. Predefined expressions can be obtained from:

# In[4]:


mf.list_predefined_gradients()


# We can see that the name `"Temperature Gradient"` is in fact a predefined gradient. Omitting calling the `register_gradient` method will in this case print the following message upon calling `solve`:
# ```
# Automatic registration of 'TemperatureGradient' as grad(Temperature).
# ```
# meaning that a predefined gradient name has been found and registered as the UFL expression $\nabla T$.
# 
# We finally solve the non-linear problem using a default Newton non-linear solver. The `solve` method returns the number of Newton iterations (4 in the present case) and converged status .

# In[5]:


problem.solve(T.vector())


# We finally check that the thermal conductivity coefficient $k$, computed from the ratio between the horizontal heat flux and temperature gradient matches the temperature-dependent expressions implemented in the MFront behaviour.

# In[7]:


j = problem.fluxes["HeatFlux"].function
g = problem.gradients["TemperatureGradient"].function
k_gauss = j.vector().get_local()[::2]/g.vector().get_local()[::2]
T_gauss = problem.state_variables["external"]["Temperature"].function.vector().get_local()
A = 0.0375;
B = 2.165e-4;
k_ref = 1/(A + B*T_gauss)
plt.plot(T_gauss, k_gauss, 'o', label="FE")
plt.plot(T_gauss, k_ref, '.', label="ref")
plt.xlabel(r"Temperature $T\: (K)$")
plt.ylabel(r"Thermal conductivity $k\: (W.m^{-1}.K^{-1})$")
plt.legend()
plt.show()


# In[ ]: