Merge branch 'tangent-blocks'

demos/finite_strain_elastoplasticity.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Finite-strain elastoplasticity within the logarithmic strain framework\n",
+    "\n",
+    "This demo is dedicated to the resolution of a finite-strain elastoplastic problem using the logarithmic strain framework proposed in <cite data-cite=\"miehe_anisotropic_2002\">(Miehe et al., 2002)</cite>. \n",
+    "\n",
+    "## Logarithmic strains \n",
+    "\n",
+    "This framework expresses constitutive relations between the Hencky strain measure $\\boldsymbol{H} = \\dfrac{1}{2}\\log (\\boldsymbol{F}^T\\cdot\\boldsymbol{F})$ and its dual stress measure $\\boldsymbol{T}$. This approach makes it possible to extend classical small strain constitutive relations to a finite-strain setting. In particular, the total (Hencky) strain can be split **additively** into many contributions (elastic, plastic, thermal, swelling, etc.). Its trace is also linked with the volume change $J=\\exp(\\operatorname{tr}(\\boldsymbol{H}))$. As a result, the deformation gradient $\\boldsymbol{F}$ is used for expressing the Hencky strain $\\boldsymbol{H}$, a small-strain constitutive law is then written for the $(\\boldsymbol{H},\\boldsymbol{T})$-pair and the dual stress $\\boldsymbol{T}$ is then post-processed to an appropriate stress measure such as the Cauchy stress $\\boldsymbol{\\sigma}$ or Piola-Kirchhoff stresses.\n",
+    "\n",
+    "## MFront implementation\n",
+    "\n",
+    "<font color=red> TODO:\n",
+    "* Write the gallery example comments and refer to it?\n",
+    "\n",
+    "* or use standard brick with comments here ?\n",
+    "</font> \n",
+    "\n",
+    "\n",
+    "## FEniCS implementation\n",
+    "\n",
+    "We define a box mesh representing half of a beam oriented along the $x$-direction. The beam will be fully clamped on its left side and symmetry conditions will be imposed on its right extremity. The loading consists of a uniform self-weight."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "from dolfin import *\n",
+    "import mfront_wrapper as mf\n",
+    "import numpy as np\n",
+    "import ufl\n",
+    "\n",
+    "length, width, height = 1., 0.04, 0.1\n",
+    "nx, ny, nz = 30, 5, 10\n",
+    "mesh = BoxMesh(Point(0, -width/2, -height/2.), Point(length, width/2, height/2.), nx, ny, nz)\n",
+    "\n",
+    "V = VectorFunctionSpace(mesh, \"CG\", 2)\n",
+    "u = Function(V, name=\"Displacement\")\n",
+    "\n",
+    "def left(x, on_boundary):\n",
+    "    return near(x[0], 0) and on_boundary\n",
+    "def right(x, on_boundary):\n",
+    "    return near(x[0], length) and on_boundary\n",
+    "\n",
+    "bc = [DirichletBC(V, Constant((0.,)*3), left),\n",
+    "      DirichletBC(V.sub(0), Constant(0.), right)]\n",
+    "\n",
+    "selfweight = Expression((\"0\", \"0\", \"-t*qmax\"), t=0., qmax = 50e6, degree=0)\n",
+    "\n",
+    "file_results = XDMFFile(\"results/finite_strain_plasticity.xdmf\")\n",
+    "file_results.parameters[\"flush_output\"] = True\n",
+    "file_results.parameters[\"functions_share_mesh\"] = True"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The `MFrontNonlinearMaterial` instance is loaded from the MFront `LogarithmicStrainPlasticity` behaviour. This behaviour is a finite-strain behaviour (`material.finite_strain=True`) which relies on a kinematic description using the total deformation gradient $\\boldsymbol{F}$. By default, a MFront behaviour always returns the Cauchy stress as the stress measure after integration. However, the stress variable dual to the deformation gradient is the first Piloa-Kirchhoff (PK1) stress. An internal option of the MGIS interface is therefore used in the finite-strain context to return the PK1 stress as the \"flux\" associated to the \"gradient\" $\\boldsymbol{F}$. Both quantities are non-symmetric tensors, aranged as a 9-dimensional vector in 3D following [MFront conventions on tensors](http://tfel.sourceforge.net/tensors.html)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "StandardFiniteStrainBehaviour\n",
+      "F_CAUCHY\n",
+      "['DeformationGradient'] [9]\n",
+      "['FirstPiolaKirchhoffStress'] [9]\n"
+     ]
+    }
+   ],
+   "source": [
+    "material = mf.MFrontNonlinearMaterial(\"../materials/src/libBehaviour.so\",\n",
+    "                                      \"LogarithmicStrainPlasticity\")\n",
+    "print(material.behaviour.getBehaviourType())\n",
+    "print(material.behaviour.getKinematic())\n",
+    "print(material.get_gradient_names(), material.get_gradient_sizes())\n",
+    "print(material.get_flux_names(), material.get_flux_sizes())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The `MFrontNonlinearProblem` instance must therefore register the deformation gradient as `Identity(3)+grad(u)`. This again done automatically since `\"DeformationGradient\"` is a predefined gradient. The following message will be shown upon calling `solve`:\n",
+    "```\n",
+    "Automatic registration of 'DeformationGradient' as I + (grad(Displacement)).\n",
+    "```\n",
+    "The loading is then defined and, as for the [small-strain elastoplasticity example](small_strain_elastoplasticity.ipynb), state variables include the `ElasticStrain` and `EquivalentPlasticStrain` since the same behaviour is used as in the small-strain case with the only difference that the total strain is now given by the Hencky strain measure. In particular, the `ElasticStrain` is still a symmetric tensor (vector of dimension 6). Note that it has not been explicitly defined as a state variable in the MFront behaviour file since this is done automatically when using the `IsotropicPlasticMisesFlow` parser.\n",
+    "\n",
+    "Finally, we setup a few parameters of the Newton non-linear solver."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "(6,)\n"
+     ]
+    }
+   ],
+   "source": [
+    "problem = mf.MFrontNonlinearProblem(u, material, bcs=bc)\n",
+    "problem.set_loading(dot(selfweight, u)*dx)\n",
+    "\n",
+    "p = problem.get_state_variable(\"EquivalentPlasticStrain\")\n",
+    "epsel = problem.get_state_variable(\"ElasticStrain\")\n",
+    "print(ufl.shape(epsel))\n",
+    "\n",
+    "prm = problem.solver.parameters\n",
+    "prm[\"absolute_tolerance\"] = 1e-6\n",
+    "prm[\"relative_tolerance\"] = 1e-6\n",
+    "prm[\"linear_solver\"] = \"mumps\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "During the load incrementation, we monitor the evolution of the vertical downwards displacement at the middle of the right extremity.\n",
+    "\n",
+    "This simulation is a bit heavy to run so we suggest running it in parallel:\n",
+    "```bash\n",
+    "mpirun -np 16 python3\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Increment  1\n",
+      "Automatic registration of 'DeformationGradient' as I + (grad(Displacement)).\n",
+      "\n",
+      "Automatic registration of 'Temperature' as a Constant value = 293.15.\n",
+      "\n",
+      "Increment  2\n",
+      "Increment  3\n",
+      "Increment  4\n"
+     ]
+    }
+   ],
+   "source": [
+    "P0 = FunctionSpace(mesh, \"DG\", 0)\n",
+    "p_avg = Function(P0, name=\"Plastic strain\")\n",
+    "\n",
+    "Nincr = 30\n",
+    "load_steps = np.linspace(0., 1., Nincr+1)\n",
+    "results = np.zeros((Nincr+1, 3))\n",
+    "for (i, t) in enumerate(load_steps[1:]):\n",
+    "    selfweight.t = t\n",
+    "    print(\"Increment \", i+1)\n",
+    "    problem.solve(u.vector())\n",
+    "\n",
+    "    results[i+1, 0] = -u(length, 0, 0)[2]\n",
+    "    results[i+1, 1] = t\n",
+    "\n",
+    "    file_results.write(u, t)\n",
+    "    p_avg.assign(project(p, P0))\n",
+    "    file_results.write(p_avg, t)\n",
+    "\n",
+    "plt.figure()\n",
+    "plt.plot(results[:, 0], results[:, 1], \"-o\")\n",
+    "plt.xlabel(r\"Displacement\")\n",
+    "plt.ylabel(\"Load\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The load-displacement curve exhibits a classical elastoplastic behaviour rapidly followed by a stiffening behaviour due to membrane catenary effects. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.9"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

demos/logarithmic_strain_plasticity.py

+import matplotlib.pyplot as plt
 from dolfin import *
 import mfront_wrapper as mf
 import numpy as np
@@ -17,63 +18,47 @@ def right(x, on_boundary):
 
 bc = [DirichletBC(V, Constant((0.,)*3), left),
       DirichletBC(V.sub(0), Constant(0.), right)]
-facets = MeshFunction("size_t", mesh, 2)
-AutoSubDomain(right).mark(facets, 1)
-ds = Measure("ds", subdomain_data=facets)
-
-# Building function with vx=1 on left boundary to compute reaction
-v = Function(V)
-bcv = DirichletBC(V.sub(0), Constant(1.), left)
-bcv.apply(v.vector())
 
-Vpost = FunctionSpace(mesh, "CG", 1)
-file_results = XDMFFile("results/beam_GL_plasticity.xdmf")
+file_results = XDMFFile("results/finie_strain_plasticity.xdmf")
 file_results.parameters["flush_output"] = True
 file_results.parameters["functions_share_mesh"] = True
 
-
 selfweight = Expression(("0", "0", "-t*qmax"), t=0., qmax = 50e6, degree=0)
+
 material = mf.MFrontNonlinearMaterial("../materials/src/libBehaviour.so",
                                       "LogarithmicStrainPlasticity")
-problem = mf.MFrontNonlinearProblem(u, material)
+problem = mf.MFrontNonlinearProblem(u, material, bcs=bc)
 problem.set_loading(dot(selfweight, u)*dx)
-problem.bc = bc
 
-p = problem.get_state_variable(name="EquivalentPlasticStrain")
-assert (ufl.shape(p) == ())
-epsel = problem.get_state_variable(name="ElasticStrain")
+p = problem.get_state_variable("EquivalentPlasticStrain")
+epsel = problem.get_state_variable("ElasticStrain")
 print(ufl.shape(epsel))
 
 prm = problem.solver.parameters
 prm["absolute_tolerance"] = 1e-6
 prm["relative_tolerance"] = 1e-6
+prm["linear_solver"] = "mumps"
 
 P0 = FunctionSpace(mesh, "DG", 0)
 p_avg = Function(P0, name="Plastic strain")
-Nitermax, tol = 20, 1e-4  # parameters of the Newton-Raphson procedure
+
 Nincr = 30
 load_steps = np.linspace(0., 1., Nincr+1)
-file_results.write(u, 0)
-file_results.write(p_avg, 0)
 results = np.zeros((Nincr+1, 3))
 for (i, t) in enumerate(load_steps[1:]):
     selfweight.t = t
 
-    problem.solve(u.vector())
-    results[i+1, 0] = assemble(-u[2]*ds(1))/(width*height)
+    # problem.solve(u.vector())
+
+    results[i+1, 0] = u(length, 0, 0)
     results[i+1, 1] = t
-    results[i+1, 2] = assemble(action(-problem.L, v))
 
     file_results.write(u, t)
-#    file_results.write(p_avg, t)
+    p_avg.assign(project(p, V0))
+    file_results.write(p_avg, t)
 
-import matplotlib.pyplot as plt
 plt.figure()
 plt.plot(results[:, 0], results[:, 1], "-o")
 plt.xlabel(r"Displacement")
 plt.ylabel("Load")
-plt.figure()
-plt.plot(results[:, 1], results[:, 2], "-o")
-plt.xlabel(r"Load")
-plt.ylabel("Horizontal Reaction")
-plt.show()
+

demos/multiphase_model.py

+#!/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 = 0.5
+height = 1.
+mesh = RectangleMesh(Point(0., 0.), Point(width, height), 100, 100)
+# mesh = BoxMesh(Point(0., 0., 0.), Point(width, height, height), 100, 2, 2)
+
+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": 10.,
+            "MatrixPoissonRatio": 0.45,
+            "FiberYoungModulus": 10000,
+            "FiberPoissonRatio": 0.3,
+            "FiberVolumeFraction": 0.01,
+            "Size": 1/16.}
+
+material = mf.MFrontNonlinearMaterial("../materials/src/libBehaviour.so",
+                                      "MultiphaseModel",
+                                      hypothesis="plane_strain",
+                                      material_properties=mat_prop)
+
+problem = mf.MFrontNonlinearProblem(u, material, quadrature_degree=0, bcs=bc)
+problem.register_gradient("MatrixStrain", sym(grad(u1)))
+problem.register_gradient("FiberStrain", sym(grad(u2)))
+problem.register_gradient("RelativeDisplacement", diag(u2-u1))
+# E = Function(problem.W0)
+# E.interpolate(Expression("20e3", degree=1))
+# mat_prop["MatrixYoungModulus"] = E
+# problem.material.update_material_properties(mat_prop)
+
+problem.solve(u.vector())
+
+x = np.linspace(0, 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^\textrm{matrix}$ EF")
+plt.plot(x, np.array([u(xi, height/2)[2] for xi in x]), "or", label=r"$u_x^\textrm{fiber}$ EF")
+
+E1 = mat_prop["MatrixYoungModulus"]
+nu1 = mat_prop["MatrixPoissonRatio"]
+E2 = mat_prop["FiberYoungModulus"]
+nu2 = mat_prop["FiberPoissonRatio"]
+eta = mat_prop["FiberVolumeFraction"]
+s = mat_prop["Size"]
+mu1 = E1/2/(1+nu1)
+mu2 = E2/2/(1+nu2)
+lmbda1 = E1*nu1/(1+nu1)/(1-2*nu1)
+lmbda2 = E2*nu2/(1+nu2)/(1-2*nu2)
+kappa = 12*mu1*mu2/((1-eta)*mu2+eta*mu1)/s**2
+M1 = lmbda1+2*mu1
+M2 = lmbda2+2*mu2
+Er = eta*E2/(1-nu2**2)
+ahom = M1+Er
+l = sqrt(M1*Er/ahom/kappa)
+um = lmbda1/ahom*(x+l*(Er/M1)*np.sinh(x/l)/np.cosh(width/l))
+ur = lmbda1/ahom*(x-l*np.sinh(x/l)/np.cosh(width/l))
+plt.plot(x, um, '--b', label=r"$u_x^\textrm{matrix}$ ref.")
+plt.plot(x, ur, '--r', label=r"$u_x^\textrm{fiber}$ ref.")
+plt.legend()
+plt.show()
\ No newline at end of file

demos/non_linear_heat_equation.py

+#!/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
+import meshio
+
+fname = "../meshes/rod3D.msh"
+msh = meshio.read("../meshes/rod3D.msh")
+for cell in msh.cells:
+    if cell.type == "triangle":
+        triangle_cells = cell.data
+    elif  cell.type == "tetra":
+        tetra_cells = cell.data
+
+for key in msh.cell_data_dict["gmsh:physical"].keys():
+    if key == "triangle":
+        triangle_data = msh.cell_data_dict["gmsh:physical"][key]
+    elif key == "tetra":
+        tetra_data = msh.cell_data_dict["gmsh:physical"][key]
+tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells},
+                         cell_data={"name_to_read":[tetra_data]})
+# triangle_mesh =meshio.Mesh(points=msh.points,
+#                            cells=[("triangle", triangle_cells)],
+#                            cell_data={"name_to_read":[triangle_data]})
+triangle_mesh =meshio.Mesh(points=msh.points,
+                           cells=[("triangle", triangle_cells)],
+                           cell_data={"name_to_read":[triangle_data]})
+meshio.write("mesh.xdmf", tetra_mesh)
+meshio.write("mf.xdmf", triangle_mesh)
+
+# remove blank spaces in field_date keys
+tags = dict([(key.strip(), value) for (key, value) in msh.field_data.items()])
+
+
+mesh = Mesh()
+mvc = MeshValueCollection("size_t", mesh, 3)
+with XDMFFile("mesh.xdmf") as infile:
+    infile.read(mesh)
+    infile.read(mvc, "name_to_read")
+cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
+mvc = MeshValueCollection("size_t", mesh, 2)
+with XDMFFile("mf.xdmf") as infile:
+    infile.read(mvc, "name_to_read")
+mf = MeshFunction("size_t", mesh, 2)
+
+def get_tag(tag):
+    return tags[tag][0]
+
+dx = Measure("dx", domain=mesh, subdomain_data=cf)
+ds = Measure("ds", domain=mesh, subdomain_data=mf)
+
+V = FunctionSpace(mesh, "CG", 2)
+T = Function(V, name="Temperature")
+
+Text = 1e3
+v = Function(V)
+bc = DirichletBC(V, Constant(Tl), mf, get_tag("00SOO"))
+bc[0].apply(v.vector())
+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)))
\ No newline at end of file

demos/nonlinear_heat_transfer_3D.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Stationnary non-linear heat transfer: 3D problem and performance comparisons\n",
+    "\n",
+    "## Description of the non-linear constitutive heat transfer law\n",
+    "\n",
+    "This example is a direct continuation of the [previous 2D example on non-linear heat transfer](nonlinear_heat_transfer.ipynb). The present computations will use the same behaviour `StationaryHeatTransfer.mfront`.\n",
+    "\n",
+    "![3D fuel rod example](supplementary_files/fuel_rod_solution.png)\n",
+    "\n",
+    "## FEniCS implementation\n",
+    "\n",
+    "We now consider a portion of nuclear fuel rod (Uranium Dioxide $\\text{UO}_2$) subject to an external imposed temperature $T_{ext}=1000\\text{ K}$ and uniform volumetric heat source $r=300 \\text{ MW/m}^3$. From the steady state heat balance equation $\\operatorname{div}\\boldsymbol{j} = r$, the variational formulation is now:\n",
+    "\n",
+    "\\begin{equation}\n",
+    "F(\\widehat{T}) = \\int_\\Omega \\boldsymbol{j}(T,\\nabla T)\\cdot\\nabla \\widehat{T}\\,\\text{dx} + \\int_\\Omega r \\widehat{T} \\,\\text{dx}=0 \\quad \\forall \\widehat{T}\n",
+    "\\end{equation}\n",
+    "\n",
+    "which fits the general default format of a `MFrontNonlinearProblem`:\n",
+    "\\begin{equation}\n",
+    "F(\\widehat{u}) = \\sum_i \\int_\\Omega \\boldsymbol{\\sigma}_i(\\boldsymbol{u})\\cdot \\boldsymbol{g}_i(\\widehat{u})\\,\\text{dx} -L(\\widehat{u}) =0 \\quad \\forall \\widehat{u}\n",
+    "\\end{equation}\n",
+    "\n",
+    "where $(\\boldsymbol{\\sigma}_i,\\boldsymbol{g}_i)$ are pairs of dual flux/gradient and here the external loading form $L$ is given by $-\\int_\\Omega r \\widehat{T} \\,\\text{dx}$. Compared to the previous example, we just add this source term using the `set_loading` method. Here we use a quadratic interpolation for the temperature field and external temperature is imposed on the surface numbered 12. Finally, we also rely on automatic registration of the gradient and external state variables as explained in the first demo."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from dolfin import *\n",
+    "import mfront_wrapper as mf\n",
+    "from time import time\n",
+    "\n",
+    "mesh = Mesh()\n",
+    "mvc = MeshValueCollection(\"size_t\", mesh, 3)\n",
+    "with XDMFFile(\"mesh.xdmf\") as infile:\n",
+    "    infile.read(mesh)\n",
+    "mvc = MeshValueCollection(\"size_t\", mesh, 2)\n",
+    "with XDMFFile(\"mf.xdmf\") as infile:\n",
+    "    infile.read(mvc, \"name_to_read\")\n",
+    "facets = cpp.mesh.MeshFunctionSizet(mesh, mvc)\n",
+    "\n",
+    "V = FunctionSpace(mesh, \"CG\", 2)\n",
+    "T = Function(V, name=\"Temperature\")\n",
+    "T_ = TestFunction(V)\n",
+    "dT = TrialFunction(V)\n",
+    "T0 = Constant(300.)\n",
+    "T.interpolate(T0)\n",
+    "\n",
+    "Text = Constant(1e3)\n",
+    "bc = DirichletBC(V, Text, facets, 12)\n",
+    "\n",
+    "r = Constant(3e8)\n",
+    "\n",
+    "quad_deg = 2\n",
+    "material = mf.MFrontNonlinearMaterial(\"../materials/src/libBehaviour.so\",\n",
+    "                                      \"StationaryHeatTransfer\")\n",
+    "problem = mf.MFrontNonlinearProblem(T, material, quadrature_degree=quad_deg, bcs=bc)\n",
+    "problem.set_loading(-r*T*dx)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The `solve` method computing time is monitored:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+