First commit with ReissnerMindlin

2D_elasticity.txt~

+=========================
+ 2D elasticity in Fenics
+=========================
+
+:math:`\int_{\Omega}\underline{\underline{\sigma}}:\underline{\underline{\varepsilon}}[v]d\Omega`

doc/2D_elasticity.rst

+=========================
+ 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 displacement-based 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 
+
+.. code-block:: 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 out-of-plane \epsilon_{zz}  
+
+Hence, the 2D constitutive relation can be defined as follows, by changing only the value of the Lamé coefficient \lambda. 
+
+
+.. code-block:: 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)

doc/2D_elasticity.rst~

+=========================
+ 2D elasticity in Fenics
+=========================
+

doc/Makefile

+# Makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line.
+SPHINXOPTS    =
+SPHINXBUILD   = sphinx-build
+PAPER         =
+BUILDDIR      = _build
+
+# User-friendly check for sphinx-build
+ifeq ($(shell which $(SPHINXBUILD) >/dev/null 2>&1; echo $$?), 1)
+$(error The '$(SPHINXBUILD)' command was not found. Make sure you have Sphinx installed, then set the SPHINXBUILD environment variable to point to the full path of the '$(SPHINXBUILD)' executable. Alternatively you can add the directory with the executable to your PATH. If you don't have Sphinx installed, grab it from http://sphinx-doc.org/)
+endif
+
+# Internal variables.
+PAPEROPT_a4     = -D latex_paper_size=a4
+PAPEROPT_letter = -D latex_paper_size=letter
+ALLSPHINXOPTS   = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) .
+# the i18n builder cannot share the environment and doctrees with the others
+I18NSPHINXOPTS  = $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) .
+
+.PHONY: help clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest gettext
+
+help:
+	@echo "Please use \`make <target>' where <target> is one of"
+	@echo "  html       to make standalone HTML files"
+	@echo "  dirhtml    to make HTML files named index.html in directories"
+	@echo "  singlehtml to make a single large HTML file"
+	@echo "  pickle     to make pickle files"
+	@echo "  json       to make JSON files"
+	@echo "  htmlhelp   to make HTML files and a HTML help project"
+	@echo "  qthelp     to make HTML files and a qthelp project"
+	@echo "  devhelp    to make HTML files and a Devhelp project"
+	@echo "  epub       to make an epub"
+	@echo "  latex      to make LaTeX files, you can set PAPER=a4 or PAPER=letter"
+	@echo "  latexpdf   to make LaTeX files and run them through pdflatex"
+	@echo "  latexpdfja to make LaTeX files and run them through platex/dvipdfmx"
+	@echo "  text       to make text files"
+	@echo "  man        to make manual pages"
+	@echo "  texinfo    to make Texinfo files"
+	@echo "  info       to make Texinfo files and run them through makeinfo"
+	@echo "  gettext    to make PO message catalogs"
+	@echo "  changes    to make an overview of all changed/added/deprecated items"
+	@echo "  xml        to make Docutils-native XML files"
+	@echo "  pseudoxml  to make pseudoxml-XML files for display purposes"
+	@echo "  linkcheck  to check all external links for integrity"
+	@echo "  doctest    to run all doctests embedded in the documentation (if enabled)"
+
+clean:
+	rm -rf $(BUILDDIR)/*
+
+html:
+	$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
+	@echo
+	@echo "Build finished. The HTML pages are in $(BUILDDIR)/html."
+
+dirhtml:
+	$(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml
+	@echo
+	@echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml."
+
+singlehtml:
+	$(SPHINXBUILD) -b singlehtml $(ALLSPHINXOPTS) $(BUILDDIR)/singlehtml
+	@echo
+	@echo "Build finished. The HTML page is in $(BUILDDIR)/singlehtml."
+
+pickle:
+	$(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle
+	@echo
+	@echo "Build finished; now you can process the pickle files."
+
+json:
+	$(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json
+	@echo
+	@echo "Build finished; now you can process the JSON files."
+
+htmlhelp:
+	$(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp
+	@echo
+	@echo "Build finished; now you can run HTML Help Workshop with the" \
+	      ".hhp project file in $(BUILDDIR)/htmlhelp."
+
+qthelp:
+	$(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp
+	@echo
+	@echo "Build finished; now you can run "qcollectiongenerator" with the" \
+	      ".qhcp project file in $(BUILDDIR)/qthelp, like this:"
+	@echo "# qcollectiongenerator $(BUILDDIR)/qthelp/NumericaltoursofcontinuummechanicsusingFEniCS.qhcp"
+	@echo "To view the help file:"
+	@echo "# assistant -collectionFile $(BUILDDIR)/qthelp/NumericaltoursofcontinuummechanicsusingFEniCS.qhc"
+
+devhelp:
+	$(SPHINXBUILD) -b devhelp $(ALLSPHINXOPTS) $(BUILDDIR)/devhelp
+	@echo
+	@echo "Build finished."
+	@echo "To view the help file:"
+	@echo "# mkdir -p $$HOME/.local/share/devhelp/NumericaltoursofcontinuummechanicsusingFEniCS"
+	@echo "# ln -s $(BUILDDIR)/devhelp $$HOME/.local/share/devhelp/NumericaltoursofcontinuummechanicsusingFEniCS"
+	@echo "# devhelp"
+
+epub:
+	$(SPHINXBUILD) -b epub $(ALLSPHINXOPTS) $(BUILDDIR)/epub
+	@echo
+	@echo "Build finished. The epub file is in $(BUILDDIR)/epub."
+
+latex:
+	$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
+	@echo
+	@echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex."
+	@echo "Run \`make' in that directory to run these through (pdf)latex" \
+	      "(use \`make latexpdf' here to do that automatically)."
+
+latexpdf:
+	$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
+	@echo "Running LaTeX files through pdflatex..."
+	$(MAKE) -C $(BUILDDIR)/latex all-pdf
+	@echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex."
+
+latexpdfja:
+	$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
+	@echo "Running LaTeX files through platex and dvipdfmx..."
+	$(MAKE) -C $(BUILDDIR)/latex all-pdf-ja
+	@echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex."
+
+text:
+	$(SPHINXBUILD) -b text $(ALLSPHINXOPTS) $(BUILDDIR)/text
+	@echo
+	@echo "Build finished. The text files are in $(BUILDDIR)/text."
+
+man:
+	$(SPHINXBUILD) -b man $(ALLSPHINXOPTS) $(BUILDDIR)/man
+	@echo
+	@echo "Build finished. The manual pages are in $(BUILDDIR)/man."
+
+texinfo:
+	$(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo
+	@echo
+	@echo "Build finished. The Texinfo files are in $(BUILDDIR)/texinfo."
+	@echo "Run \`make' in that directory to run these through makeinfo" \
+	      "(use \`make info' here to do that automatically)."
+
+info:
+	$(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo
+	@echo "Running Texinfo files through makeinfo..."
+	make -C $(BUILDDIR)/texinfo info
+	@echo "makeinfo finished; the Info files are in $(BUILDDIR)/texinfo."
+
+gettext:
+	$(SPHINXBUILD) -b gettext $(I18NSPHINXOPTS) $(BUILDDIR)/locale
+	@echo
+	@echo "Build finished. The message catalogs are in $(BUILDDIR)/locale."
+
+changes:
+	$(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) $(BUILDDIR)/changes
+	@echo
+	@echo "The overview file is in $(BUILDDIR)/changes."
+
+linkcheck:
+	$(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) $(BUILDDIR)/linkcheck
+	@echo
+	@echo "Link check complete; look for any errors in the above output " \
+	      "or in $(BUILDDIR)/linkcheck/output.txt."
+
+doctest:
+	$(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest
+	@echo "Testing of doctests in the sources finished, look at the " \
+	      "results in $(BUILDDIR)/doctest/output.txt."
+
+xml:
+	$(SPHINXBUILD) -b xml $(ALLSPHINXOPTS) $(BUILDDIR)/xml
+	@echo
+	@echo "Build finished. The XML files are in $(BUILDDIR)/xml."
+
+pseudoxml:
+	$(SPHINXBUILD) -b pseudoxml $(ALLSPHINXOPTS) $(BUILDDIR)/pseudoxml
+	@echo
+	@echo "Build finished. The pseudo-XML files are in $(BUILDDIR)/pseudoxml."

doc/_build/html/.buildinfo

+# Sphinx build info version 1
+# This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done.
+config: 8034634c671452947be2ed52157c7e13
+tags: 645f666f9bcd5a90fca523b33c5a78b7

doc/_build/html/2D_elasticity.html

+
+<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
+  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
+
+<html xmlns="http://www.w3.org/1999/xhtml">
+  <head>
+    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
+    <title>2D linear elasticity &#8212; Numerical tours of continuum mechanics using FEniCS 1.0 documentation</title>
+    <link rel="stylesheet" href="_static/agogo.css" type="text/css" />
+    <link rel="stylesheet" href="_static/pygments.css" type="text/css" />
+    <script type="text/javascript">
+      var DOCUMENTATION_OPTIONS = {
+        URL_ROOT:    './',
+        VERSION:     '1.0',
+        COLLAPSE_INDEX: false,
+        FILE_SUFFIX: '.html',
+        HAS_SOURCE:  true,
+        SOURCELINK_SUFFIX: '.txt'
+      };
+    </script>
+    <script type="text/javascript" src="_static/jquery.js"></script>
+    <script type="text/javascript" src="_static/underscore.js"></script>
+    <script type="text/javascript" src="_static/doctools.js"></script>
+    <script type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
+    <link rel="index" title="Index" href="genindex.html" />
+    <link rel="search" title="Search" href="search.html" />
+    <link rel="next" title="Reissner-Mindlin plate with Quadrilaterals" href="demo/reissner_mindlin/reissner_mindlin.py.html" />
+    <link rel="prev" title="Introduction" href="intro.html" /> 
+  </head>
+  <body>
+    <div class="header-wrapper" role="banner">
+      <div class="header">
+        <div class="headertitle"><a
+          href="index.html">Numerical tours of continuum mechanics using FEniCS 1.0 documentation</a></div>
+        <div class="rel" role="navigation" aria-label="related navigation">
+          <a href="intro.html" title="Introduction"
+             accesskey="P">previous</a> |
+          <a href="demo/reissner_mindlin/reissner_mindlin.py.html" title="Reissner-Mindlin plate with Quadrilaterals"
+             accesskey="N">next</a> |
+          <a href="genindex.html" title="General Index"
+             accesskey="I">index</a>
+        </div>
+       </div>
+    </div>
+
+    <div class="content-wrapper">
+      <div class="content">
+        <div class="document">
+            
+      <div class="documentwrapper">
+        <div class="bodywrapper">
+          <div class="body" role="main">
+            
+  <div class="section" id="d-linear-elasticity">
+<h1>2D linear elasticity<a class="headerlink" href="#d-linear-elasticity" title="Permalink to this headline">¶</a></h1>
+<p>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 displacement-based finite element formulation.</p>
+<p>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.</p>
+<p>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 :</p>
+<p>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</p>
+<div class="highlight-python"><div class="highlight"><pre><span></span><span class="k">def</span> <span class="nf">eps</span><span class="p">(</span><span class="n">v</span><span class="p">):</span>
+     <span class="k">return</span> <span class="n">sym</span><span class="p">(</span><span class="n">grad</span><span class="p">(</span><span class="n">v</span><span class="p">))</span>
+</pre></div>
+</div>
+<p>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.</p>
+<p>In the plane stress case, an out-of-plane epsilon_{zz}</p>
+<p>Hence, the 2D constitutive relation can be defined as follows, by changing only the value of the Lamé coefficient lambda.</p>
+<div class="highlight-python"><div class="highlight"><pre><span></span><span class="k">def</span> <span class="nf">sigma</span><span class="p">(</span><span class="n">v</span><span class="p">):</span>
+    <span class="n">dim</span> <span class="o">=</span> <span class="mi">2</span>
+    <span class="k">if</span> <span class="n">plane_stress</span><span class="p">:</span>
+        <span class="n">lmbda</span> <span class="o">=</span> <span class="n">lmbda_plane_stress</span>
+    <span class="k">else</span><span class="p">:</span>
+        <span class="n">lmbda</span> <span class="o">=</span> <span class="n">lmbda_plane_strain</span>
+    <span class="k">return</span> <span class="n">lmbda</span><span class="o">*</span><span class="n">tr</span><span class="p">(</span><span class="n">eps</span><span class="p">(</span><span class="n">v</span><span class="p">))</span><span class="o">*</span><span class="n">Identity</span><span class="p">(</span><span class="n">dim</span><span class="p">)</span> <span class="o">+</span> <span class="mf">2.0</span><span class="o">*</span><span class="n">mu</span><span class="o">*</span><span class="n">eps</span><span class="p">(</span><span class="n">v</span><span class="p">)</span>
+</pre></div>
+</div>
+</div>
+
+
+          </div>
+        </div>
+      </div>
+        </div>
+        <div class="sidebar">
+          <h3>Table Of Contents</h3>
+          <ul class="current">
+<li class="toctree-l1"><a class="reference internal" href="intro.html">Introduction</a></li>
+<li class="toctree-l1 current"><a class="current reference internal" href="#">2D linear elasticity</a></li>
+<li class="toctree-l1"><a class="reference internal" href="demo/reissner_mindlin/reissner_mindlin.py.html">Reissner-Mindlin plate with Quadrilaterals</a></li>
+</ul>
+
+          <div role="search">
+            <h3 style="margin-top: 1.5em;">Search</h3>
+            <form class="search" action="search.html" method="get">
+                <input type="text" name="q" />
+                <input type="submit" value="Go" />
+                <input type="hidden" name="check_keywords" value="yes" />
+                <input type="hidden" name="area" value="default" />
+            </form>
+          </div>
+        </div>
+        <div class="clearer"></div>
+      </div>
+    </div>
+
+    <div class="footer-wrapper">
+      <div class="footer">
+        <div class="left">
+          <div role="navigation" aria-label="related navigaton">
+            <a href="intro.html" title="Introduction"
+              >previous</a> |
+            <a href="demo/reissner_mindlin/reissner_mindlin.py.html" title="Reissner-Mindlin plate with Quadrilaterals"
+              >next</a> |
+            <a href="genindex.html" title="General Index"
+              >index</a>
+          </div>
+          <div role="note" aria-label="source link">
+              <br/>
+              <a href="_sources/2D_elasticity.rst.txt"
+                rel="nofollow">Show Source</a>
+          </div>
+        </div>
+
+        <div class="right">
+          
+    <div class="footer" role="contentinfo">
+        &#169; Copyright 2016, Jeremy Bleyer.
+      Created using <a href="http://sphinx-doc.org/">Sphinx</a> 1.6.5.
+    </div>
+        </div>
+        <div class="clearer"></div>
+      </div>
+    </div>
+
+  </body>
+</html>
\ No newline at end of file

doc/_build/html/_sources/2D_elasticity.rst.txt

+=========================
+ 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 displacement-based 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 
+
+.. code-block:: 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 out-of-plane \epsilon_{zz}  
+
+Hence, the 2D constitutive relation can be defined as follows, by changing only the value of the Lamé coefficient \lambda. 
+
+
+.. code-block:: 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)

doc/_build/html/_sources/2D_elasticity.txt

+=========================
+ 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 displacement-based finite element formulation.
+
+We will illustrate this on the case of a cantilever beam modeled as 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 
+
+.. code-block:: 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 out-of-plane \epsilon_{zz}  
+
+Hence, the 2D constitutive relation can be defined as follows, by changing only the value of the Lamé coefficient \lambda. 
+
+
+.. code-block:: python
+
+    def sigma(v):
+	dim = v.geometric_dimension()
+	if plane_stress:
+	    lmbda = lmbda_plane_stress
+	else:
+	    lmbda = lmbda_plane_strain
+	return lmbda*tr(eps(v))*Identity(dim) + 2.0*mu*eps(v)

doc/_build/html/_sources/demo/reissner_mindlin/reissner_mindlin.py.rst.txt

+
+.. _ReissnerMindlin:
+
+==========================================
+Reissner-Mindlin plate with Quadrilaterals
+==========================================
+
+
+This program solves the Reissner-Mindlin plate equations on the unit
+square with uniform transverse loading and fully clamped boundary conditions.
+
+It uses quadrilateral cells and selective reduced integration (SRI) to
+remove shear-locking issues in the thin plate limit. Both linear and 
+quadratic interpolation are considered for the transverse deflection 
+:math:`w` and rotation :math:`\underline{\theta}`. 
+
+The solution for :math:`w` in this demo will look as follows:
+
+.. image:: clamped_40x40.png
+   :scale: 40 %
+
+
+Material parameters for isotropic linear elastic behavior are first defined::
+
+ from fenics import *
+
+ E = Constant(1e3)
+ nu = Constant(0.3)
+  
+Plate bending stiffness :math:`D=\dfrac{Eh^3}{12(1-\nu^2)}` and shear stiffness :math:`F = \kappa Gh`
+with a shear correction factor :math:`\kappa = 5/6` for a homogeneous plate
+of thickness :math:`h`::
+
+ thick = Constant(1e-2)
+ D = E*thick**3/(1-nu**2)/12.
+ F = E/2/(1+nu)*thick*5./6.
+
+The uniform loading :math:`f` is scaled by the plate thickness so that the deflection converges to a
+constant value in the thin plate Love-Kirchhoff limit::
+
+ f = Constant(-thick**3)
+
+The unit square mesh is divided in :math:`N\times N` quadrilaterals::
+
+ N = 10
+ mesh = UnitSquareMesh.create(N, N, CellType.Type_quadrilateral)
+
+Continuous interpolation using of degree :math:`d=\texttt{deg}` is chosen for both deflection and rotation::
+
+ deg = 1
+ We = FiniteElement("Lagrange", mesh.ufl_cell(), deg)
+ Te = VectorElement("Lagrange", mesh.ufl_cell(), deg)
+ V = FunctionSpace(mesh,MixedElement([We,Te]))
+
+Clamped boundary conditions on the lateral boundary are defined as::
+
+ def border(x, on_boundary):
+     return on_boundary
+      
+ bc =  [DirichletBC(V,Constant((0.,0.,0.)), border)]
+  
+
+Some useful functions for implementing generalized constitutive relations are now
+defined::
+
+ def strain2voigt(eps):
+     return as_vector([eps[0,0],eps[1,1],2*eps[0,1]])
+ def voigt2stress(S):
+     return as_tensor([[S[0],S[2]],[S[2],S[1]]])
+ def curv(u):
+     (w,theta) = split(u)
+     return sym(grad(theta))
+ def shear_strain(u):
+     (w,theta) = split(u)
+     return theta-grad(w)
+ def bending_moment(u):
+     DD = as_tensor([[D,nu*D,0],[nu*D,D,0],[0,0,D*(1-nu)/2.]])
+     return voigt2stress(dot(DD,strain2voigt(curv(u))))
+ def shear_force(u):
+     return F*shear_strain(u)
+
+
+The contribution of shear forces to the total energy is under-integrated using
+a custom quadrature rule of degree :math:`2d-2` i.e. for linear (:math:`d=1`) 
+quadrilaterals, the shear energy is integrated as if it were constant (1 Gauss point instead of 2x2)
+and for quadratic (:math:`d=2`) quadrilaterals, as if it were quadratic (2x2 Gauss points instead of 3x3)::
+
+ u = Function(V)
+ u_ = TestFunction(V)
+ du = TrialFunction(V)
+
+ dx_shear = dx(metadata={"quadrature_degree":2*deg-2})
+ 
+ L = f*u_[0]*dx
+ a = inner(bending_moment(u_),curv(du))*dx + dot(shear_force(u_),shear_strain(du))*dx_shear
+  
+
+We then solve for the solution and export the relevant fields to XDMF files ::
+
+ solve(a == L, u, bc)
+ 
+ (w,theta) = split(u)
+  
+ Vw = FunctionSpace(mesh,We)
+ Vt = FunctionSpace(mesh,Te)
+ ww = Function(Vw, name="Deflection")
+ tt = Function(Vt, name="Rotation")
+ ww.assign(project(w, Vw))
+ tt.assign(project(theta, Vt))
+  
+ file_results = XDMFFile("RM_results.xdmf")
+ file_results.parameters["flush_output"] = True
+ file_results.parameters["functions_share_mesh"] = True
+ file_results.write(ww, 0.)
+ file_results.write(tt, 0.)
+ 
+The solution is compared to the Kirchhoff analytical solution::
+
+ print "Kirchhoff deflection:", -1.265319087e-3*float(f/D)
+ print "Reissner-Mindlin FE deflection:", min(ww.vector().get_local()) # point evaluation for quads
+                                                                       # is not implemented in fenics 2017.2
\ No newline at end of file

doc/_build/html/_sources/index.rst.txt

+.. Numerical tours of continuum mechanics using FEniCS documentation master file, created by
+   sphinx-quickstart on Wed Jun  8 21:25:10 2016.
+   You can adapt this file completely to your liking, but it should at least
+   contain the root `toctree` directive.
+
+Welcome to Numerical tours of Computational Continuum Mechanics using FEniCS
+============================================================================
+
+Contents:
+