Skip to content

GitLab

M
mfront-wrapper
Commit 03f2794a authored by Jeremy BLEYER's avatar Jeremy BLEYER
1 Browse files
Options

Adds a few notebook demos

parent a6efe243
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 10056 additions and 252 deletions
README.md 100644 → 100755
View file @ 03f2794a
File mode changed from 100644 to 100755
demos/finite_strain_elastoplasticity.ipynb 0 → 100644
View file @ 03f2794a
{
"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 100644 → 100755
View file @ 03f2794a
......@@ -18,29 +18,19 @@ 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, bcs=bc)
problem.set_loading(dot(selfweight, u)*dx)
p = problem.get_state_variable("EquivalentPlasticStrain")
assert (ufl.shape(p) == ())
epsel = problem.get_state_variable("ElasticStrain")
print(ufl.shape(epsel))
......@@ -51,29 +41,24 @@ 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)
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/non_linear_heat_equation.py
View file @ 03f2794a
......@@ -8,25 +8,60 @@ Created on Fri Feb 1 08:49:11 2019
from dolfin import *
import mfront_wrapper as mf
import numpy as np
import meshio
length = 30e-3
width = 5.4e-3
mesh = RectangleMesh(Point(0., 0.), Point(length, width), 100, 10)
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()])
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
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)
Tl = 300
Tr = 800
def get_tag(tag):
return tags[tag][0]
bc = [DirichletBC(V, Constant(Tl), left),
DirichletBC(V, Constant(Tr), right)]
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)
......
demos/nonlinear_heat_transfer.ipynb 100644 → 100755
View file @ 03f2794a
......@@ -8,6 +8,189 @@
"\n",
"## Description of the non-linear constitutive heat transfer law\n",
"\n",
"The thermal material is described by the following non linear Fourier\n",
"Law:\n",
"\n",
"$$\n",
"\\mathbf{j}=-k\\left(T\\right)\\,\\mathbf{\\nabla} T\n",
"$$\n",
"\n",
"where $\\mathbf{j}$ is the heat flux and $\\mathbf{\\nabla} T$ is the\n",
"temperature gradient.\n",
"\n",
"Expression of the thermal conductivity\n",
"--------------------------------------\n",
"\n",
"The thermal conductivity is assumed to be given by:\n",
"\n",
"$$\n",
"k\\left(T\\right)={\\displaystyle \\frac{\\displaystyle 1}{\\displaystyle A+B\\,T}}\n",
"$$\n",
"\n",
"This expression accounts for the phononic contribution to the thermal\n",
"conductivity.\n",
"\n",
"Derivatives\n",
"-----------\n",
"\n",
"As discussed below, the consistent linearisation of the heat transfer\n",
"equilibrium requires to compute:\n",
"\n",
"- the derivative\n",
" ${\\displaystyle \\frac{\\displaystyle \\partial \\mathbf{j}}{\\displaystyle \\partial \\mathbf{\\nabla} T}}$\n",
" of the heat flux with respect to the temperature gradient.\n",
" ${\\displaystyle \\frac{\\displaystyle \\partial \\mathbf{j}}{\\displaystyle \\partial \\mathbf{\\nabla} T}}$\n",
" is given by: $$\n",
" {\\displaystyle \\frac{\\displaystyle \\partial \\mathbf{j}}{\\displaystyle \\partial \\mathbf{\\nabla} T}}=-k\\left(T\\right)\\,\\matrix{I}\n",
" $$\n",
"- the derivative\n",
" ${\\displaystyle \\frac{\\displaystyle \\partial \\mathbf{j}}{\\displaystyle \\partial T}}$\n",
" of the heat flux with respect to the temperature.\n",
" ${\\displaystyle \\frac{\\displaystyle \\partial \\mathbf{j}}{\\displaystyle \\partial T}}$\n",
" is given by: $$\n",
" {\\displaystyle \\frac{\\displaystyle \\partial \\mathbf{j}}{\\displaystyle \\partial T}}=-{\\displaystyle \\frac{\\displaystyle \\partial k\\left(T\\right)}{\\displaystyle \\partial T}}\\,\\mathbf{\\nabla} T=B\\,k^{2}\\,\\mathbf{\\nabla} T\n",
" $$\n",
"\n",
"`MFront`’ implementation\n",
"========================\n",
"\n",
"Choice of the the domain specific language\n",
"------------------------------------------\n",
"\n",
"Every `MFront` file is handled by a domain specific language (DSL), which\n",
"aims at providing the most suitable abstraction for a particular choice\n",
"of behaviour and integration algorithm. See `mfront mfront --list-dsl`\n",
"for a list of the available DSLs.\n",
"\n",
"The name of DSL’s handling generic behaviours ends with\n",
"`GenericBehaviour`. The first part of a DSL’s name is related to the\n",
"integration algorithm used.\n",
"\n",
"In the case of this non linear transfer behaviour, the heat flux is\n",
"explicitly computed from the temperature and the temperature gradient.\n",
"The `DefaultGenericBehaviour` is the most suitable choice:\n",
"\n",
"``` cxx\n",
"@DSL DefaultGenericBehaviour;\n",
"```\n",
"\n",
"Some metadata\n",
"-------------\n",
"\n",
"The following lines define the name of the behaviour, the name of the\n",
"author and the date of its writing:\n",
"\n",
"``` cxx\n",
"@Behaviour StationaryHeatTransfer;\n",
"@Author Thomas Helfer;\n",
"@Date 15/02/2019;\n",
"```\n",
"\n",
"Gradients and fluxes\n",
"--------------------\n",
"\n",
"Generic behaviours relate pairs of gradients and fluxes. Gradients and\n",
"fluxes are declared independently but the first declared gradient is\n",
"assumed to be conjugated with the first declared fluxes and so on…\n",
"\n",
"The temperature gradient is declared as follows (note that Unicode characters are supported):\n",
"\n",
"``` cxx\n",
"@Gradient TemperatureGradient ∇T;\n",
"∇T.setGlossaryName(\"TemperatureGradient\");\n",
"```\n",
"\n",
"Note that we associated to `∇T` the glossary name `TemperatureGradient`.\n",
"This is helpful for the calling code.\n",
"\n",
"After this declaration, the following variables will be defined:\n",
"\n",
"- The temperature gradient `∇T` at the beginning of the time step.\n",
"- The increment of the temperature gradient `Δ∇T` over the time step.\n",
"\n",
"The heat flux is then declared as follows:\n",
"\n",
"``` cxx\n",
"@Flux HeatFlux j;\n",
"j.setGlossaryName(\"HeatFlux\");\n",
"```\n",
"\n",
"In the following code blocks, `j` will be the heat flux at the end of\n",
"the time step.\n",
"\n",
"Tangent operator blocks\n",
"-----------------------\n",
"\n",
"By default, the derivatives of the gradients with respect to the fluxes\n",
"are declared. Thus the variable `∂j∕∂Δ∇T` is automatically declared.\n",
"\n",
"However, as discussed in the next section, the consistent linearisation\n",
"of the thermal equilibrium requires to return the derivate of the heat\n",
"flux with respect to the increment of the temperature (or equivalently\n",
"with respect to the temperature at the end of the time step).\n",
"\n",
"``` cxx\n",
"@AdditionalTangentOperatorBlock ∂j∕∂ΔT;\n",
"```\n",
"\n",
"Parameters\n",
"----------\n",
"\n",
"The `A` and `B` coefficients that appears in the definition of the\n",
"thermal conductivity are declared as parameters:\n",
"\n",
"``` cxx\n",
"@Parameter real A = 0.0375;\n",
"@Parameter real B = 2.165e-4;\n",
"```\n",
"\n",
"Parameters are stored globally and can be modified from the calling\n",
"solver or from `python` in the case of the coupling with `FEniCS`\n",
"discussed below.\n",
"\n",
"Local variable\n",
"--------------\n",
"\n",
"A local variable is accessible in each code blocks.\n",
"\n",
"Here, we declare the thermal conductivity `k` as a local variable in\n",
"order to be able to compute its value during the behaviour integration\n",
"and to reuse this value when computing the tangent operator.\n",
"\n",
"``` cxx\n",
"@LocalVariable thermalconductivity k;\n",
"```\n",
"\n",
"Integration of the behaviour\n",
"----------------------------\n",
"\n",
"The behaviour integration is straightforward: one starts to compute the\n",
"temperature at the end of the time step, then we compute the thermal\n",
"conductivity (at the end of the time step) and the heat flux using the\n",
"temperature gradient (at the end of the time step).\n",
"\n",
"``` cxx\n",
"@Integrator{\n",
" // temperature at the end of the time step\n",
" const auto T_ = T + ΔT;\n",
" // thermal conductivity\n",
" k = 1 / (A + B ⋅ T_);\n",
" // heat flux\n",
" j = -k ⋅ (∇T + Δ∇T);\n",
"} // end of @Integrator\n",
"```\n",
"\n",
"Tangent operator\n",
"----------------\n",
"\n",
"The computation of the tangent operator blocks is equally simple:\n",
"\n",
"``` cxx\n",
"@TangentOperator {\n",
" ∂j∕∂Δ∇T = -k ⋅ tmatrix<N, N, real>::Id();\n",
" ∂j∕∂ΔT = B ⋅ k ⋅ k ⋅ (∇T + Δ∇T);\n",
"} // end of @TangentOperator \n",
"```\n",
"## FEniCS implementation\n",
"\n",
"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."
......@@ -211,10 +394,10 @@
"source": [
"j = problem.fluxes[\"HeatFlux\"].function\n",
"g = problem.gradients[\"TemperatureGradient\"].function\n",
"k_gauss = j.vector().get_local()[::2]/g.vector().get_local()[::2]\n",
"k_gauss = -j.vector().get_local()[::2]/g.vector().get_local()[::2]\n",
"T_gauss = problem.state_variables[\"external\"][\"Temperature\"].function.vector().get_local()\n",
"A = 0.0375;\n",
"B = 2.165e-4;\n",
"A = material.get_parameter(\"A\");\n",
"B = material.get_parameter(\"B\");\n",
"k_ref = 1/(A + B*T_gauss)\n",
"plt.plot(T_gauss, k_gauss, 'o', label=\"FE\")\n",
"plt.plot(T_gauss, k_ref, '.', label=\"ref\")\n",
......@@ -248,7 +431,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
"version": "3.6.9"
}
},
"nbformat": 4,
......
demos/nonlinear_heat_transfer.py deleted 100644 → 0
View file @ a6efe243
#!/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")