From 9c1a45c885255f15790935a73588682544cb86ac Mon Sep 17 00:00:00 2001
From: Jeremy Bleyer
Date: Thu, 19 Apr 2018 18:49:23 +0200
Subject: [PATCH] Finished 2D elasticity

doc/2D_elasticity.rst  34 
examples/elasticity/2D_elasticity.py.rst  26 ++++++++++++++++
toolbox_list  9 ++++++
3 files changed, 32 insertions(+), 37 deletions()
delete mode 100644 doc/2D_elasticity.rst
diff git a/doc/2D_elasticity.rst b/doc/2D_elasticity.rst
deleted file mode 100644
index e8dc3fe..0000000
 a/doc/2D_elasticity.rst
+++ /dev/null
@@ 1,34 +0,0 @@
=========================
 2D linear elasticity
=========================

In this first numerical tour, we will show how to compute a static solution for a 2D isotropic linear elastic medium, either in plane stress or in plane strain, using FEniCS in a tradtional displacementbased finite element formulation.

We will illustrate this on the case of a cantilever beam modeled as a 2D medium of dimensions L x H. We first define the geometrical parameters and mesh density which will be used. The mesh is generated using the RectangleMesh function and we choose a crossed configuration for the mesh structure.

We now define the material parameters which are here given in terms of a Young's modulus and a Poisson coefficient. In the following, we will need to define the constitutive relation relating \sigma as a function of \varepsilon. Let us recall that the general expression for a 3D medium of the linear elastic isotropic constitutive relation is given by :

here we will consider a 2D model either in plane strain or in plane stress. Irrespective of this choice, we will work only with a 2D displacement vector u and will subsequently define the strain operator \epsilon as follows

.. codeblock:: python

 def eps(v):
 return sym(grad(v))

which computes the 2x2 plane components of the symmetrized gradient tensor of any 2D vectorial field.
In the plane strain case, the full 3D strain tensor is defined as follows, so that the 2x2 plane part of the stress tensor is defined in the same way as for the 3D case.

In the plane stress case, an outofplane \epsilon_{zz}

Hence, the 2D constitutive relation can be defined as follows, by changing only the value of the LamÃ© coefficient \lambda.


.. codeblock:: python

 def sigma(v):
 dim = 2
 if plane_stress:
 lmbda = lmbda_plane_stress
 else:
 lmbda = lmbda_plane_strain
 return lmbda*tr(eps(v))*Identity(dim) + 2.0*mu*eps(v)
diff git a/examples/elasticity/2D_elasticity.py.rst b/examples/elasticity/2D_elasticity.py.rst
index 5344063..d08865e 100644
 a/examples/elasticity/2D_elasticity.py.rst
+++ b/examples/elasticity/2D_elasticity.py.rst
@@ 14,7 +14,8 @@ Introduction
In this first numerical tour, we will show how to compute a small strain solution for
a 2D isotropic linear elastic medium, either in plane stress or in plane strain,
in a tradtional displacementbased finite element formulation. Extension to 3D
+in a tradtional displacementbased finite element formulation. The corresponding
+file can be obtained from :download:`2D_elasticity.py`.Extension to 3D
is straightforward and an example can be found in the :ref:`ModalAnalysis` example.
We consider here the case of a cantilever beam modeled as a 2D medium of dimensions
@@ 160,4 +161,25 @@ EulerBernoulli beam theory which is here :math:`w_{beam} = \dfrac{qL^4}{8EI}`::
print("Beam theory deflection:", float(3*rho_g*L**4/2/E/H**3))
One finds :math:`w_{FE} = 5.8638\text{e3}` against :math:`w_{beam} = 5.8594\text{e3}`
that is a 0.07% difference.
\ No newline at end of file
+that is a 0.07% difference.
+
+
+The stress tensor must be projected on an appropriate function space in order to
+evaluate pointwise values or export it for Paraview vizualisation. Here we choose
+to describe it as a (2D) tensor and project it onto a piecewise constant function
+space::
+
+ Vsig = TensorFunctionSpace(mesh, "DG", degree=0)
+ sig = Function(Vsig, name="Stress")
+ sig.assign(project(sigma(u), Vsig))
+ print("Stress at (0,H):", sig(0, H))
+
+Fields can be exported in a suitable format for vizualisation using Paraview.
+VTKbased extensions (.pvd,.vtu) are not suited for multiple fields and parallel
+writing/reading. Prefered output format is now .xdmf::
+
+ file_results = XDMFFile("elasticity_results.xdmf")
+ file_results.parameters["flush_output"] = True
+ file_results.parameters["functions_share_mesh"] = True
+ file_results.write(u, 0.)
+ file_results.write(sig, 0.)
\ No newline at end of file
diff git a/toolbox_list b/toolbox_list
index cbec9dd..56d946c 100644
 a/toolbox_list
+++ b/toolbox_list
@@ 1,3 +1,10 @@
+#############
+Introduction
+#############
+
+Workflow : geometry with Gmsh, dolfinconvert, export Paraview
+
+
###############
Linear problems
###############
@@ 15,7 +22,7 @@ mixed formulation using HuWashizu
> volumetric locking
Elastodynamics
* time integration Newmark scheme, theta
+* time integration Newmark scheme, theta (lumped mass matrix, efficient solving)
Poroelasticity

GitLab