Added build folder

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: 97675b57fb48f4bcf85cca2fef8e548c
+tags: 645f666f9bcd5a90fca523b33c5a78b7

doc/_build/html/_downloads/2D_elasticity.py

+#
+# ..    # gedit: set fileencoding=utf8 :
+#
+# .. _LinearElasticity2D:
+#
+#
+# =========================
+#  2D linear elasticity
+# =========================
+#
+#
+# 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 displacement-based finite element formulation. The corresponding
+# file can be obtained from :download:`2D_elasticity.py`.
+#
+# .. seealso::
+#
+#  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
+# :math:`L\times  H`. Geometrical parameters and mesh density are first defined
+# and the rectangular domain is  generated using the ``RectangleMesh`` function.
+# We also choose a criss-crossed structured mesh::
+
+from __future__ import print_function
+from fenics import *
+
+L = 25.
+H = 1.
+Nx = 250
+Ny = 10
+mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny, "crossed")
+
+
+# Constitutive relation
+# ---------------------
+#
+# We now define the material parameters which are here given in terms of a Young's
+# modulus :math:`E` and a Poisson coefficient :math:`\nu`. In the following, we will
+# need to define the constitutive relation between the stress tensor :math:`\boldsymbol{\sigma}`
+# and the strain tensor :math:`\boldsymbol{\varepsilon}`. Let us recall
+# that the general expression of the linear elastic isotropic constitutive relation
+# for a 3D medium is given by:
+#
+# .. math::
+#  \boldsymbol{\sigma} = \lambda \text{tr}(\boldsymbol{\varepsilon})\mathbf{1} + 2\mu\boldsymbol{\varepsilon}
+#  :label: constitutive_3D
+#
+# for a natural (no prestress) initial state where the Lamé coefficients are given by:
+#
+# .. math::
+#  \lambda = \dfrac{E\nu}{(1+\nu)(1-2\nu)}, \quad \mu = \dfrac{E}{2(1+\nu)}
+#  :label: Lame_coeff
+#
+# In this demo, we consider a 2D model either in plane strain or in plane stress conditions.
+# Irrespective of this choice, we will work only with a 2D displacement vector :math:`\boldsymbol{u}=(u_x,u_y)`
+# and will subsequently define the strain operator ``eps`` as follows::
+
+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:
+#
+# .. math::
+#  \boldsymbol{\varepsilon} = \begin{bmatrix} \varepsilon_{xx} & \varepsilon_{xy} & 0\\
+#  \varepsilon_{xy} & \varepsilon_{yy} & 0 \\ 0 & 0 & 0\end{bmatrix}
+#
+# so that the 2x2 plane part of the stress tensor is defined in the same way as for the 3D case
+# (the out-of-plane stress component being given by :math:`\sigma_{zz}=\lambda(\varepsilon_{xx}+\varepsilon_{yy})`.
+#
+# In the plane stress case, an out-of-plane strain component :math:`\varepsilon_{zz}`
+# must be considered so that :math:`\sigma_{zz}=0`. Using this condition in the
+# 3D constitutive relation, one has :math:`\varepsilon_{zz}=-\dfrac{\lambda}{\lambda+2\mu}(\varepsilon_{xx}+\varepsilon_{yy})`.
+# Injecting into :eq:`constitutive_3D`, we have for the 2D plane stress relation:
+#
+# .. math::
+#  \boldsymbol{\sigma} = \lambda^* \text{tr}(\boldsymbol{\varepsilon})\mathbf{1} + 2\mu\boldsymbol{\varepsilon}
+#
+# where :math:`\boldsymbol{\sigma}, \boldsymbol{\varepsilon}, \mathbf{1}` are 2D tensors and with
+# :math:`\lambda^* = \dfrac{2\lambda\mu}{\lambda+2\mu}`. Hence, the 2D constitutive relation
+# is identical to the plane strain case by changing only the value of the Lamé coefficient :math:`\lambda`.
+# We can then have::
+
+E = Constant(1e5)
+nu = Constant(0.3)
+model = "plane_stress"
+
+mu = E/2/(1+nu)
+lmbda = E*nu/(1+nu)/(1-2*nu)
+if model == "plane_stress":
+    lmbda = 2*mu*lmbda/(lmbda+2*mu)
+
+def sigma(v):
+    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
+
+# .. note::
+#  Note that we used the variable name ``lmbda`` to avoid any confusion with the
+#  lambda functions of Python
+#
+#  We also used an intrinsic formulation of the constitutive relation. Example of
+#  constitutive relation implemented with a matrix/vector engineering notation
+#  will be provided in the :ref:`OrthotropicElasticity` example.
+#
+#
+# Variational formulation
+# -----------------------
+#
+# For this example, we consider a continuous polynomial interpolation of degree 2
+# and a uniformly distributed loading :math:`\boldsymbol{f}=(0,-f)` corresponding
+# to the beam self-weight. The continuum mechanics variational formulation (obtained
+# from the virtual work principle) is given by:
+#
+# .. math::
+#  \text{Find } \boldsymbol{u}\in V \text{ s.t. } \int_{\Omega}
+#  \boldsymbol{\sigma}(\boldsymbol{u}):\boldsymbol{\varepsilon}(\boldsymbol{v}) d\Omega
+#  = \int_{\Omega} \boldsymbol{f}\cdot\boldsymbol{v}  d\Omega \quad \forall\boldsymbol{v} \in V
+#
+# which translates into the following FEniCS code::
+
+rho_g = 1e-3
+f = Constant((0,-rho_g))
+
+V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
+du = TrialFunction(V)
+u_ = TestFunction(V)
+a = inner(sigma(du), eps(u_))*dx
+l = inner(f, u_)*dx
+
+
+# Resolution
+# ----------
+#
+# Fixed displacements are imposed on the left part of the beam, the ``solve``
+# function is then called and solution is plotted by deforming the mesh::
+
+def left(x, on_boundary):
+    return near(x[0],0.)
+
+bc = DirichletBC(V, Constant((0.,0.)), left)
+
+u = Function(V, name="Displacement")
+solve(a == l, u, bc)
+
+plot(1e3*u, mode="displacement")
+
+# The (amplified) solution should look like this:
+#
+# .. image:: cantilever_deformed.png
+#    :scale: 15%
+#
+#
+# Validation and post-processing
+# ------------------------------
+#
+# The maximal deflection is compared against the analytical solution from
+# Euler-Bernoulli beam theory which is here :math:`w_{beam} = \dfrac{qL^4}{8EI}`::
+
+print("Maximal deflection:", -u(L,H/2.)[1])
+print("Beam theory deflection:", float(3*rho_g*L**4/2/E/H**3))
+
+# One finds :math:`w_{FE} = 5.8638\text{e-3}` against :math:`w_{beam} = 5.8594\text{e-3}`
+# 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.
+# VTK-based 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

doc/_build/html/_downloads/axisymmetric_elasticity.py

+
+# coding: utf-8
+
+# # Axisymmetric formulation for elastic structures of revolution
+# 
+# In this numerical tour, we will deal with axisymmetric problems of elastic solids. We will consider a solid of revolution around a fixed axis $(Oz)$, the loading, boundary conditions and material properties being also invariant with respect to a rotation along the symmetry axis. The solid cross-section in a plane $\theta=\text{cst}$ will be represented by a two-dimensional domain $\omega$ for which the first spatial variable (`x[0]` in FEniCS) will represent the radial coordinate $r$ whereas the second spatial variable will denote the axial variable $z$.
+# 
+# ## Problem position
+# 
+# We will investigate here the case of a hollow hemisphere of inner (resp. outer) radius $R_i$ (resp. $R_e$). Due to the revolution symmetry, the 2D cross-section corresponds to a quarter of a hollow cylinder.
+
+# In[50]:
+
+
+from __future__ import print_function
+from dolfin import *
+from mshr import *
+import matplotlib.pyplot as plt
+get_ipython().magic(u'matplotlib notebook')
+
+Re = 11.
+Ri = 9.
+rect = Rectangle(Point(0., 0.), Point(Re, Re))
+domain = Circle(Point(0., 0.), Re, 100) - Circle(Point(0., 0.), Ri, 100)
+domain = domain - Rectangle(Point(0., -Re), Point(-Re, Re))                 - Rectangle(Point(0., 0.), Point(Re, -Re))
+
+mesh = generate_mesh(domain, 40)
+plot(mesh)
+
+class Bottom(SubDomain):
+    def inside(self, x, on_boundary):
+        return near(x[1], 0) and on_boundary
+class Left(SubDomain):
+    def inside(self, x, on_boundary):
+        return near(x[0], 0) and on_boundary
+class Outer(SubDomain):
+    def inside(self, x, on_boundary):
+        return near(sqrt(x[0]**2+x[1]**2), Re, 1e-1) and on_boundary  
+    
+facets = MeshFunction("size_t", mesh, 1)
+facets.set_all(0)
+Bottom().mark(facets, 1)
+Left().mark(facets, 2)
+Outer().mark(facets, 3)
+ds = Measure("ds", subdomain_data=facets)
+
+
+# ## Definition of axisymmetric strains
+# 
+# For axisymmetric conditions, the unkown displacement field is of the form:
+# 
+# \begin{equation}
+# \boldsymbol{u} = u_r(r,z)\boldsymbol{e}_r + u_z(r,z)\boldsymbol{e}_z
+# \end{equation}
+# 
+# As a result, we will work with a standard VectorFunctionSpace of dimension 2. The associated strain components are however given by:
+# 
+# \begin{equation}
+# \boldsymbol{\varepsilon} = \begin{bmatrix} \partial_r u_r & 0 & (\partial_z u_r + \partial_r u_z)/2 \\ 
+# 0 & u_r/r & 0 \\
+# (\partial_z u_r + \partial_r u_z)/2 & 0 & \partial_z u_z\end{bmatrix}_{(\boldsymbol{e}_r,\boldsymbol{e}_\theta,\boldsymbol{e}_z)}
+# \end{equation}
+# 
+# The previous relation involves explicitly the radial variable $r$, which can be obtained from the SpatialCoordinate `x[0]`, the strain-displacement relation is then defined explicitly in the `eps` function.
+# 
+# > **Note**: we could also express the strain components in the form of a vector of size 4 in alternative of the 3D tensor representation implemented below.
+
+# In[43]:
+
+
+x = SpatialCoordinate(mesh)
+
+def eps(v):
+    return sym(as_tensor([[v[0].dx(0), 0, v[0].dx(1)],
+                          [0, v[0]/x[0], 0],
+                          [v[1].dx(0), 0, v[1].dx(1)]]))
+
+E = Constant(1e5)
+nu = Constant(0.3)
+mu = E/2/(1+nu)
+lmbda = E*nu/(1+nu)/(1-2*nu)
+def sigma(v):
+    return lmbda*tr(eps(v))*Identity(3) + 2.0*mu*eps(v)
+
+
+# ## Resolution
+# 
+# The rest of the formulation is similar to the 2D elastic case with a small difference in the integration measure. Indeed, the virtual work principle reads as:
+# \begin{equation}
+#  \text{Find } \boldsymbol{u}\in V \text{ s.t. } \int_{\Omega}
+#  \boldsymbol{\sigma}(\boldsymbol{u}):\boldsymbol{\varepsilon}(\boldsymbol{v}) d\Omega
+#  = \int_{\partial \Omega_T} \boldsymbol{T}\cdot\boldsymbol{v}  dS \quad \forall\boldsymbol{v} \in V
+#  \end{equation}
+# where $\boldsymbol{T}$ is the imposed traction on some part $\partial \Omega_T$ of the domain boundary.
+# 
+# In axisymmetric conditions, the full 3D domain $\Omega$ can be decomposed as $\Omega = \omega \times [0;2\pi]$ where the interval represents the $\theta$ variable. The integration measures therefore reduce to $d\Omega = d\omega\cdot(rd\theta)$ and $dS = ds\cdot(rd\theta)$ where $dS$ is the surface integration measure on the 3D domain $\Omega$ and $ds$ its counterpart on the cross-section boundary $\partial \omega$. Exploiting the invariance of all fields with respect to $\theta$, the previous virtual work principle is reformulated on the cross-section only as follows:
+# 
+# \begin{equation}
+#  \text{Find } \boldsymbol{u}\in V \text{ s.t. } \int_{\omega}
+#  \boldsymbol{\sigma}(\boldsymbol{u}):\boldsymbol{\varepsilon}(\boldsymbol{v}) rd\omega
+#  = \int_{\partial \omega_T} \boldsymbol{T}\cdot\boldsymbol{v}  rds \quad \forall\boldsymbol{v} \in V
+#  \end{equation}
+# where the $2\pi$ constants arising from the integration on $\theta$ have been cancelled on both sides. As a result, the bilinear and linear form are similar to the plane 2D case with the exception of the additional $r$ term in the integration measures.
+# 
+# The final formulation is therefore pretty straightforward. Since a uniform pressure loading is applied on the outer boundary, we will also need the exterior normal vector to define the work of external forces form.
+
+# In[44]:
+
+
+n = FacetNormal(mesh)
+p = Constant(10.)
+
+V = VectorFunctionSpace(mesh, 'CG', degree=2)
+du = TrialFunction(V)
+u_ = TestFunction(V)
+a = inner(sigma(du), eps(u_))*x[0]*dx
+l = inner(-p*n, u_)*x[0]*ds(3)
+
+u = Function(V, name="Displacement")
+
+
+# First, smooth contact conditions are assumed on both $r=0$ (`Left`) and $z=0$ (`Bottom`) boundaries. For this specific case, the solution corresponds to the classical hollow sphere under external pressure with a purely radial displacement:
+# \begin{equation}
+# u_r(r) = -\dfrac{R_e^3}{R_e^3-R_i^3}\left((1 − 2\nu)r + (1 + \nu)\dfrac{R_i^3}{2r^2}\right)\dfrac{p}{E},
+# \quad u_z=0
+# \end{equation}
+
+# In[45]:
+
+
+bcs = [DirichletBC(V.sub(1), Constant(0), facets, 1),
+       DirichletBC(V.sub(0), Constant(0), facets, 2)]
+solve(a == l, u, bcs)
+print("Inwards radial displacement at (r=Re, theta=0): {:1.7f} (FE) {:1.7f} (Exact)".format(-u(Re, 0.)[0], float(Re**3/(Re**3-Ri**3)*((1-2*nu)*Re+(1+nu)*Ri**3/2/Re**2)*p/E)))
+print("Inwards radial displacement at (r=Ri, theta=0): {:1.7f} (FE) {:1.7f} (Exact)".format(-u(Ri, 0.)[0], float(Re**3/(Re**3-Ri**3)*((1-2*nu)*Ri+(1+nu)*Ri/2)*p/E)))
+
+
+# In[46]:
+
+
+plt.figure()
+plot(100*u, mode="displacement")
+plt.show()
+
+
+# The second loading case corresponds to a fully clamped condition on $z=0$, the vertical boundary remaining in smooth contact.
+
+# In[48]:
+
+
+bcs = [DirichletBC(V, Constant((0., 0.)), facets, 1),
+       DirichletBC(V.sub(0), Constant(0), facets, 2)]
+solve(a == l, u, bcs)
+plt.figure()
+plot(mesh, linewidth=0.2)
+plot(200*u, mode="displacement")
+plt.show()
+

doc/_build/html/_downloads/cantilever_modal.py

+#
+# .. _ModalAnalysis:
+#
+# ==========================================
+# Modal analysis of an elastic structure
+# ==========================================
+#
+# -------------
+# Introduction
+# -------------