Commit 98142cc2 authored by Jeremy BLEYER's avatar Jeremy BLEYER

Added DG Reissner Mindlin example

parent 5349a067
# 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
<!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
=========================
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)
=========================
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)
.. 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:
.. toctree::
:maxdepth: 2
intro
2D_elasticity
demo/reissner_mindlin/reissner_mindlin.py.rst
Indices and tables
==================
* :ref:`genindex`
* :ref:`modindex`
* :ref:`search`
=========================
Introduction
=========================
These numerical tours will introduce you to a wide variety of topics in Computational Continuum Mechanics using the finite element software FEniCS.
You can find instructions on how to install FEniCS on the FEniCS project website http://fenicsproject.org. In the following numerical tours, we will use the Python interface for the different FEniCS scripts.
FEniCS is also distributed along with an important number of documented or undocumented examples, some of them will be revisited in these tours but do not hesitate to look for other interesting examples.
In particular, we advise you to go through the documentation and tutorials to get started with FEniCS.
--------------------
A few words about FEniCS
--------------------
Contrary to most finite element softwares in the computational mechanics community, FEniCS is a multi-purpose software that allows for the discretization and numerical resolution of a wide range of problems governed by partial differential equations.
The formulation of a problem in FEniCS relies on the definition of discrete functional spaces depending on a mesh. For instance, we can define :
.. code-block:: python
# A continuous Galerkin ("CG") scalar function space interpolated with polynomials of degree 2
Vh = FunctionSpace(mesh,"CG",2)
# A discontinous Galerkin ("DG") vectorial function space (the dimension is fixed by the mesh) interpolated with polynomials of degree 1
Wh = VectorFunctionSpace(mesh,"DG",1)
The next step is to define test and trial functions which will then be used to define forms (linear, bilinear or non-linear) representing the weak variational formulation of the problem which we aim at solving.
=========================
Introduction
=========================
These numerical tours will introduce you to a wide variety of topics in Computational Continuum Mechanics using the finite element software FEniCS.
You can find instructions on how to install FEniCS on the FEniCS project website http://fenicsproject.org. In the following numerical tours, we will use the Python interface for the different FEniCS scripts.
FEniCS is also distributed along with an important number of documented or undocumented examples, some of them will be revisited in these tours but do not hesitate to look for other interesting examples.
In particular, we advise you to go through the documentation and tutorials to get started with FEniCS.
--------------------
A few words about FEniCS
--------------------
Contrary to most finite element softwares in the computational mechanics community, FEniCS is a multi-purpose software that allows for the discretization and numerical resolution of a wide range of problems governed by partial differential equations.
The formulation of a problem in FEniCS relies on the definition of discrete functional spaces depending on a mesh. For instance, we can define :
.. code-block:: python
# A continuous Galerkin ("CG") scalar function space interpolated with polynomials of degree 2
Vh = FunctionSpace(mesh,"CG",2)
# A discontinous Galerkin ("DG") vectorial function space (the dimension is fixed by the mesh) interpolated with polynomials of degree 1
Wh = VectorFunctionSpace(mesh,"DG",1)
The next step is to define test and trial functions which will then be used to define forms (linear, bilinear or non-linear) representing the weak variational formulation of the problem which we aim at solving.
Efficient assignement of vector-valued functions
------------------------------------------------
Suppose that you have a vector-valued (dimension 2) function ``v`` defined on space ``VV`` and two scalar valued functions ``u0`` and ``u1`` defined on space ``V`` that you want to assign as the components of ``v``.
Projection of ``as_vector([u0,u1])`` onto ``VV`` or interpolation of Expression is too slow for big meshes. An efficient way to perform this is to use the ``FunctionAssigner`` between space ``VV`` and spaces ``[V,V]`` as follows :
assigner = FunctionAssigner(VV, [V, V])
assigner.assign(vv, [u0, u1])
Efficient projection on DG or Quadrature spaces
------------------------------------------------
For projecting a Function on a DG or Quadrature space, that is a space with no coupling between elements, the projection can be performed element-wise. For this purpose, using the LocalSolver is much more faster than performing a global projection :
def local_project(v,V):
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv,v_)*dx(metadata=metadata)
b_proj = inner(v,v_)*dx(metadata=metadata)
solver = LocalSolver(a_proj,b_proj)
solver.factorize()
u = Function(V)
solver.solve_local_rhs(u)
return u
Local factorizations can be cached if projection is performed many times.
/*
* agogo.css_t
* ~~~~~~~~~~~
*
* Sphinx stylesheet -- agogo theme.
*
* :copyright: Copyright 2007-2017 by the Sphinx team, see AUTHORS.
* :license: BSD, see LICENSE for details.
*
*/
* {
margin: 0px;
padding: 0px;
}
body {
font-family: "Verdana", Arial, sans-serif;
line-height: 1.4em;
color: black;
background-color: #eeeeec;
}
/* Page layout */
div.header, div.content, div.footer {
width: 70em;
margin-left: auto;
margin-right: auto;
}
div.header-wrapper {
background: #555573 url(bgtop.png) top left repeat-x;
border-bottom: 3px solid #2e3436;
}
/* Default body styles */
a {
color: #ce5c00;
}
div.bodywrapper a, div.footer a {
text-decoration: underline;
}
.clearer {
clear: both;
}
.left {
float: left;
}
.right {
float: right;
}
.line-block {
display: block;
margin-top: 1em;
margin-bottom: 1em;
}
.line-block .line-block {
margin-top: 0;
margin-bottom: 0;
margin-left: 1.5em;
}
h1, h2, h3, h4 {
font-family: "Georgia", "Times New Roman", serif;
font-weight: normal;
color: #3465a4;
margin-bottom: .8em;
}
h1 {
color: #204a87;
}
h2 {
padding-bottom: .5em;
border-bottom: 1px solid #3465a4;
}
a.headerlink {
visibility: hidden;
color: #dddddd;
padding-left: .3em;
}
h1:hover > a.headerlink,
h2:hover > a.headerlink,
h3:hover > a.headerlink,
h4:hover > a.headerlink,
h5:hover > a.headerlink,
h6:hover > a.headerlink,
dt:hover > a.headerlink,
caption:hover > a.headerlink,
p.caption:hover > a.headerlink,
div.code-block-caption:hover > a.headerlink {
visibility: visible;
}
img {
border: 0;
}
div.admonition {
margin-top: 10px;
margin-bottom: 10px;
padding: 2px 7px 1px 7px;
border-left: 0.2em solid black;
}
p.admonition-title {
margin: 0px 10px 5px 0px;
font-weight: bold;
}
dt:target, .highlighted {
background-color: #fbe54e;
}
/* Header */
div.header {
padding-top: 10px;
padding-bottom: 10px;
}
div.header .headertitle {
font-family: "Georgia", "Times New Roman", serif;
font-weight: normal;
font-size: 180%;
letter-spacing: .08em;
margin-bottom: .8em;
}
div.header .headertitle a {
color: white;
}
div.header div.rel {
margin-top: 1em;
}
div.header div.rel a {
color: #fcaf3e;
letter-spacing: .1em;
text-transform: uppercase;
}
p.logo {
float: right;
}
img.logo {
border: 0;
}
/* Content */
div.content-wrapper {
background-color: white;
padding-top: 20px;
padding-bottom: 20px;
}
div.document {
width: 50em;
float: left;
}
div.body {
padding-right: 2em;
text-align: justify;
}
div.document h1 {
line-height: 120%;
}
div.document ul {
margin: 1.5em;
list-style-type: square;
}
div.document dd {
margin-left: 1.2em;
margin-top: .4em;
margin-bottom: 1em;
}
div.document .section {
margin-top: 1.7em;
}
div.document .section:first-child {
margin-top: 0px;
}
div.document div.highlight {
padding: 3px;
background-color: #eeeeec;
border-top: 2px solid #dddddd;
border-bottom: 2px solid #dddddd;
margin-top: .8em;
margin-bottom: .8em;
}
div.document div.literal-block-wrapper {
margin-top: .8em;
margin-bottom: .8em;
}
div.document div.literal-block-wrapper div.highlight {
margin: 0;
}
div.document div.code-block-caption span.caption-number {
padding: 0.1em 0.3em;
font-style: italic;
}
div.document div.code-block-caption span.caption-text {
}
div.document h2 {
margin-top: .7em;
}
div.document p {
margin-bottom: .5em;
}