Commit 5e56a751 authored by Jeremy BLEYER's avatar Jeremy BLEYER

Added periodic homog example + changes in process files for including zip and notebooks

parent 7513fac1
...@@ -24,8 +24,9 @@ sys.path.insert(0, os.path.abspath('..')) ...@@ -24,8 +24,9 @@ sys.path.insert(0, os.path.abspath('..'))
#sys.path.insert(0, os.path.abspath('../examples</')) #sys.path.insert(0, os.path.abspath('../examples</'))
sys.path.append(os.getcwd()) sys.path.append(os.getcwd())
import rstprocess import process_sources
rstprocess.process() process_sources.process()
# -- General configuration ------------------------------------------------ # -- General configuration ------------------------------------------------
# If your documentation needs a minimal Sphinx version, state it here. # If your documentation needs a minimal Sphinx version, state it here.
...@@ -35,13 +36,17 @@ rstprocess.process() ...@@ -35,13 +36,17 @@ rstprocess.process()
# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
# ones. # ones.
extensions = ['sphinx.ext.autodoc', 'sphinx.ext.doctest', 'sphinx.ext.intersphinx', extensions = ['sphinx.ext.autodoc', 'sphinx.ext.doctest', 'sphinx.ext.intersphinx',
'sphinx.ext.mathjax', 'sphinx.ext.napoleon'] 'sphinx.ext.mathjax', 'sphinx.ext.napoleon', 'nbsphinx']
# Add any paths that contain templates here, relative to this directory. # Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates'] templates_path = ['_templates']
source_parsers = {
'.md': 'recommonmark.parser.CommonMarkParser',
}
# The suffix of source filenames. # The suffix of source filenames.
source_suffix = '.rst'#,'.txt'] source_suffix = ['.rst','.md']#,'.txt']
# The encoding of source files. # The encoding of source files.
#source_encoding = 'utf-8-sig' #source_encoding = 'utf-8-sig'
...@@ -58,9 +63,9 @@ copyright = u'2016, Jeremy Bleyer' ...@@ -58,9 +63,9 @@ copyright = u'2016, Jeremy Bleyer'
# built documents. # built documents.
# #
# The short X.Y version. # The short X.Y version.
version = '1.0' version = 'master'
# The full version, including alpha/beta/rc tags. # The full version, including alpha/beta/rc tags.
release = '1.0' release = 'master'
# The language for content autogenerated by Sphinx. Refer to documentation # The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages. # for a list of supported languages.
...@@ -74,7 +79,7 @@ release = '1.0' ...@@ -74,7 +79,7 @@ release = '1.0'
# List of patterns, relative to source directory, that match files and # List of patterns, relative to source directory, that match files and
# directories to ignore when looking for source files. # directories to ignore when looking for source files.
exclude_patterns = ['_build'] exclude_patterns = ['_build', '**.ipynb_checkpoints']
# The reST default role (used for this markup: `text`) to use for all # The reST default role (used for this markup: `text`) to use for all
# documents. # documents.
...@@ -186,6 +191,7 @@ html_static_path = ['_static'] ...@@ -186,6 +191,7 @@ html_static_path = ['_static']
# Output file base name for HTML help builder. # Output file base name for HTML help builder.
htmlhelp_basename = 'NumericaltoursofcontinuummechanicsusingFEniCSdoc' htmlhelp_basename = 'NumericaltoursofcontinuummechanicsusingFEniCSdoc'
html_short_title = "Numerical tours of continuum mechanics using FEniCS"
# -- Options for LaTeX output --------------------------------------------- # -- Options for LaTeX output ---------------------------------------------
...@@ -204,7 +210,7 @@ latex_elements = { ...@@ -204,7 +210,7 @@ latex_elements = {
# (source start file, target name, title, # (source start file, target name, title,
# author, documentclass [howto, manual, or own class]). # author, documentclass [howto, manual, or own class]).
latex_documents = [ latex_documents = [
('index', 'NumericaltoursofcontinuummechanicsusingFEniCS.tex', u'Numerical tours of continuum mechanics using FEniCS Documentation', ('index', 'NumericaltoursofcontinuummechanicsusingFEniCS.tex', u'Numerical tours of continuum mechanics using FEniCS',
u'Jeremy Bleyer', 'manual'), u'Jeremy Bleyer', 'manual'),
] ]
...@@ -216,6 +222,7 @@ latex_documents = [ ...@@ -216,6 +222,7 @@ latex_documents = [
# not chapters. # not chapters.
#latex_use_parts = False #latex_use_parts = False
# If true, show page references after internal links. # If true, show page references after internal links.
#latex_show_pagerefs = False #latex_show_pagerefs = False
...@@ -234,7 +241,7 @@ latex_documents = [ ...@@ -234,7 +241,7 @@ latex_documents = [
# One entry per manual page. List of tuples # One entry per manual page. List of tuples
# (source start file, name, description, authors, manual section). # (source start file, name, description, authors, manual section).
man_pages = [ man_pages = [
('index', 'numericaltoursofcontinuummechanicsusingfenics', u'Numerical tours of continuum mechanics using FEniCS Documentation', ('index', 'numericaltoursofcontinuummechanicsusingfenics', u'Numerical tours of continuum mechanics using FEniCS',
[u'Jeremy Bleyer'], 1) [u'Jeremy Bleyer'], 1)
] ]
...@@ -248,7 +255,7 @@ man_pages = [ ...@@ -248,7 +255,7 @@ man_pages = [
# (source start file, target name, title, author, # (source start file, target name, title, author,
# dir menu entry, description, category) # dir menu entry, description, category)
texinfo_documents = [ texinfo_documents = [
('index', 'NumericaltoursofcontinuummechanicsusingFEniCS', u'Numerical tours of continuum mechanics using FEniCS Documentation', ('index', 'NumericaltoursofcontinuummechanicsusingFEniCS', u'Numerical tours of continuum mechanics using FEniCS',
u'Jeremy Bleyer', 'NumericaltoursofcontinuummechanicsusingFEniCS', 'One line description of project.', u'Jeremy Bleyer', 'NumericaltoursofcontinuummechanicsusingFEniCS', 'One line description of project.',
'Miscellaneous'), 'Miscellaneous'),
] ]
...@@ -264,3 +271,19 @@ texinfo_documents = [ ...@@ -264,3 +271,19 @@ texinfo_documents = [
# If true, do not generate a @detailmenu in the "Top" node's menu. # If true, do not generate a @detailmenu in the "Top" node's menu.
#texinfo_no_detailmenu = False #texinfo_no_detailmenu = False
#nbsphinx_prolog = """
#Go there: https://example.org/notebooks/{{ env.doc2path(env.docname, base=None) }}
#
#----
#"""
#{% set docname = env.doc2path(env.docname, base='doc') %}
nbsphinx_prolog = r"""
{% set docname = env.docname.split("/")[-1] %}
.. role:: raw-html(raw)
:format: html
.. nbinfo::
The corresponding files can be obtained from:
- Jupyter Notebook: :download:`{{docname+".ipynb"}}`
- Python script: :download:`{{docname+".py"}}`"""
\ No newline at end of file
.. 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.
Homogenization of heterogeneous materials
=========================================
Contents:
.. toctree::
:maxdepth: 1
demo/periodic_homog_elas/periodic_homog_elas.ipynb
...@@ -13,6 +13,7 @@ Contents: ...@@ -13,6 +13,7 @@ Contents:
intro intro
linear_problems linear_problems
homogenization
nonlinear_problems nonlinear_problems
demo/reissner_mindlin/reissner_mindlin.rst demo/reissner_mindlin/reissner_mindlin.rst
tips_and_tricks tips_and_tricks
......
...@@ -10,7 +10,7 @@ Contents: ...@@ -10,7 +10,7 @@ Contents:
.. toctree:: .. toctree::
:maxdepth: 1 :maxdepth: 1
demo/nonlinear_materials/vonMises_plasticity.py.rst demo/2D_plasticity/vonMises_plasticity.py.rst
...@@ -19,9 +19,11 @@ ...@@ -19,9 +19,11 @@
import sys import sys
import os import os
import shutil import shutil
import zipfile
# extensions of files which must be copied with the demo when building docs # extensions of files which must be copied with the demo when building docs
file_extensions = [".png",".gif", ".geo", ".xml", ".msh", ".pdf"] file_extensions = [".png",".gif", ".geo", ".xml", ".msh", ".pdf"]
zipfile_extensions = [".geo", ".msh", ".xml", ".py", ".ipynb"]
def process(): def process():
"""Copy demo rst files (C++ and Python) from the DOLFIN source tree """Copy demo rst files (C++ and Python) from the DOLFIN source tree
...@@ -36,7 +38,7 @@ def process(): ...@@ -36,7 +38,7 @@ def process():
# Directories to scan # Directories to scan
subdirs = ["../examples"] subdirs = ["../examples"]
retval = os.getcwd()
# Iterate over subdirectories containing demos # Iterate over subdirectories containing demos
for subdir in subdirs: for subdir in subdirs:
...@@ -55,28 +57,51 @@ def process(): ...@@ -55,28 +57,51 @@ def process():
if not os.path.isdir(version_path): if not os.path.isdir(version_path):
continue continue
# Build list of rst files in demo source directory
rst_files = [f for f in os.listdir(version_path) if os.path.splitext(f)[1] == ".rst" ]
# Create directory in documentation tree for demo # Create directory in documentation tree for demo
demo_dir = os.path.join('./demo/', demo, version) demo_dir = os.path.join('./demo/', demo)
if not os.path.exists(demo_dir): if not os.path.exists(demo_dir):
os.makedirs(demo_dir) os.makedirs(demo_dir)
ipynb_files = [f for f in os.listdir(version_path) if f.endswith("ipynb") ]
for f in ipynb_files:
ff = os.path.join(path, f)
shutil.copy(ff, demo_dir)
ret = os.system("jupyter-nbconvert --to script %s --output %s" % (ff, os.path.join("../../doc",demo_dir, os.path.splitext(f)[0])))
if not ret == 0:
raise RuntimeError("Unable to convert ipynb file to a .py ({})".format(f))
# Build list of rst files in demo source directory
rst_files = [f for f in os.listdir(version_path) if f.endswith("rst") ]
# Copy rst files into documentation demo directory and process with Pylit # Copy rst files into documentation demo directory and process with Pylit
for f in rst_files: for f in rst_files:
shutil.copy(os.path.join(version_path, f), demo_dir) shutil.copy(os.path.join(version_path, f), demo_dir)
# Run pylit on cpp.rst and py.rst files (file with 'double extensions') # Run pylit on cpp.rst and py.rst files (file with 'double extensions')
if os.path.splitext(os.path.splitext(f)[0])[1] in (".py"): if "py" in f.split("."):
rst_file = os.path.join(demo_dir, f) rst_file = os.path.join(demo_dir, f)
command = pylit_parser + " " + rst_file command = pylit_parser + " " + rst_file
ret = os.system("python "+command) ret = os.system("python "+command)
if not ret == 0: if not ret == 0:
raise RuntimeError("Unable to convert rst file to a .py ({})".format(f)) raise RuntimeError("Unable to convert rst file to a .py ({})".format(f))
tocopy_files = [f for f in os.listdir(version_path) if os.path.splitext(f)[1] in file_extensions ] tocopy_files = [f for f in os.listdir(version_path) if os.path.splitext(f)[1] in file_extensions ]
for f in tocopy_files: for f in tocopy_files:
source = os.path.join(version_path, f) source = os.path.join(version_path, f)
print("Copying {} to {}".format(source, demo_dir)) print("Copying {} to {}".format(source, demo_dir))
shutil.copy(source, demo_dir) shutil.copy(source, demo_dir)
for f in ipynb_files + rst_files:
cwd = os.getcwd()
os.chdir(demo_dir)
files = os.listdir(".")
zip_files = []
for ff in files:
if any([ff.endswith(ext) for ext in zipfile_extensions]):
zip_files.append(ff)
with zipfile.ZipFile(f.split(".")[0]+".zip", 'w') as myzip:
for a in zip_files:
myzip.write(a, compress_type=zipfile.ZIP_DEFLATED)
os.chdir(cwd)
...@@ -2,28 +2,28 @@ ...@@ -2,28 +2,28 @@
.. _vonMisesPlasticity: .. _vonMisesPlasticity:
================================================== ==================================================
Elasto-plastic analysis of a 2D von Mises material Elasto-plastic analysis of a 2D von Mises material
================================================== ==================================================
------------- -------------
Introduction Introduction
------------- -------------
This example is concerned with the incremental analysis of an elasto-plastic This example is concerned with the incremental analysis of an elasto-plastic
von Mises material. The structure response is computed using an iterative von Mises material. The structure response is computed using an iterative
predictor-corrector return mapping algorithm embedded in a Newton-Raphson global predictor-corrector return mapping algorithm embedded in a Newton-Raphson global
loop for restoring equilibrium. Due to the simple expression of the von Mises criterion, loop for restoring equilibrium. Due to the simple expression of the von Mises criterion,
the return mapping procedure is completely analytical (with linear isotropic the return mapping procedure is completely analytical (with linear isotropic
hardening). The corresponding file can be obtained from :download:`vonMises_plasticity.py`. hardening). The corresponding sources can be obtained from :download:`vonMises_plasticity.zip`.
Another implementation of von Mises plasticity can also be found at Another implementation of von Mises plasticity can also be found at
https://bitbucket.org/fenics-apps/fenics-solid-mechanics. https://bitbucket.org/fenics-apps/fenics-solid-mechanics.
We point out that the 2D nature of the problem will impose keeping We point out that the 2D nature of the problem will impose keeping
track of the out-of-plane :math:`\varepsilon_{zz}^p` plastic strain and dealing with track of the out-of-plane :math:`\varepsilon_{zz}^p` plastic strain and dealing with
representations of stress/strain states including the :math:`zz` component. Note representations of stress/strain states including the :math:`zz` component. Note
also that are not concerned with the issue of volumetric locking induced by also that are not concerned with the issue of volumetric locking induced by
incompressible plastic deformations since quadratic triangles in 2D is enough incompressible plastic deformations since quadratic triangles in 2D is enough
to mitigate the locking phenomenon. to mitigate the locking phenomenon.
The plastic strain evolution during the cylinder expansion will look like this: The plastic strain evolution during the cylinder expansion will look like this:
...@@ -35,12 +35,12 @@ The plastic strain evolution during the cylinder expansion will look like this: ...@@ -35,12 +35,12 @@ The plastic strain evolution during the cylinder expansion will look like this:
Problem position Problem position
----------------- -----------------
In FEniCS 2017.2, the FEniCS Form Compiler ``ffc`` now uses ``uflacs`` as a default In FEniCS 2017.2, the FEniCS Form Compiler ``ffc`` now uses ``uflacs`` as a default
representation instead of the old ``quadrature`` representation. However, using representation instead of the old ``quadrature`` representation. However, using
``Quadrature`` elements generates some bugs in this representation and we therefore ``Quadrature`` elements generates some bugs in this representation and we therefore
need to revert to the old representation. Deprecation warning messages are also disabled. need to revert to the old representation. Deprecation warning messages are also disabled.
See `this post <https://www.allanswered.com/post/lknbq/assemble-quadrature-representation-vs-uflacs/>`_ See `this post <https://www.allanswered.com/post/lknbq/assemble-quadrature-representation-vs-uflacs/>`_
and the corresponding `issue <https://bitbucket.org/fenics-project/ffc/issues/146/uflacs-generates-undefined-variable-for>`_ and the corresponding `issue <https://bitbucket.org/fenics-project/ffc/issues/146/uflacs-generates-undefined-variable-for>`_
for more information.:: for more information.::
from __future__ import print_function from __future__ import print_function
...@@ -64,9 +64,9 @@ can also be related to a tangent elastic modulus :math:`E_t = \dfrac{EH}{E+H}`. ...@@ -64,9 +64,9 @@ can also be related to a tangent elastic modulus :math:`E_t = \dfrac{EH}{E+H}`.
The considered problem is that of a plane strain hollow cylinder of internal (resp. external) The considered problem is that of a plane strain hollow cylinder of internal (resp. external)
radius :math:`R_i` (resp. :math:`R_e`) under internal uniform pressure :math:`q`. radius :math:`R_i` (resp. :math:`R_e`) under internal uniform pressure :math:`q`.
Only a quarter of cylinder is generated using ``Gmsh`` and converted to ``.xml`` format. :: Only a quarter of cylinder is generated using ``Gmsh`` and converted to ``.xml`` format. ::
# elastic parameters # elastic parameters
E = Constant(70e3) E = Constant(70e3)
nu = Constant(0.3) nu = Constant(0.3)
...@@ -80,7 +80,7 @@ Only a quarter of cylinder is generated using ``Gmsh`` and converted to ``.xml`` ...@@ -80,7 +80,7 @@ Only a quarter of cylinder is generated using ``Gmsh`` and converted to ``.xml``
mesh = Mesh("thick_cylinder.xml") mesh = Mesh("thick_cylinder.xml")
facets = MeshFunction("size_t", mesh, "thick_cylinder_facet_region.xml") facets = MeshFunction("size_t", mesh, "thick_cylinder_facet_region.xml")
ds = Measure('ds')[facets] ds = Measure('ds')[facets]
Function spaces will involve a standard CG space for the displacement whereas internal Function spaces will involve a standard CG space for the displacement whereas internal
state variables such as plastic strains will be represented using a **Quadrature** element. state variables such as plastic strains will be represented using a **Quadrature** element.
This choice will make it possible to express the complex non-linear material constitutive This choice will make it possible to express the complex non-linear material constitutive
...@@ -99,16 +99,16 @@ of Gauss points will be determined by the required degree of the Quadrature elem ...@@ -99,16 +99,16 @@ of Gauss points will be determined by the required degree of the Quadrature elem
W = FunctionSpace(mesh, We) W = FunctionSpace(mesh, We)
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default') W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')
W0 = FunctionSpace(mesh, W0e) W0 = FunctionSpace(mesh, W0e)
.. note:: .. note::
In older versions, it was possible to define Quadrature function spaces directly In older versions, it was possible to define Quadrature function spaces directly
using ``FunctionSpace(mesh, "Quadrature", 1)``. This is no longer the case since using ``FunctionSpace(mesh, "Quadrature", 1)``. This is no longer the case since
FEniCS 2016.1 (see `this issue <https://bitbucket.org/fenics-project/dolfin/issues/757/functionspace-mesh-quadrature-1-broken-in>`_). Instead, Quadrature elements must first be defined FEniCS 2016.1 (see `this issue <https://bitbucket.org/fenics-project/dolfin/issues/757/functionspace-mesh-quadrature-1-broken-in>`_). Instead, Quadrature elements must first be defined
by specifying the associated degree and quadrature scheme before defining the by specifying the associated degree and quadrature scheme before defining the
associated FunctionSpace. associated FunctionSpace.
Various functions are defined to keep track of the current internal state and Various functions are defined to keep track of the current internal state and
currently computed increments:: currently computed increments::
...@@ -124,21 +124,21 @@ currently computed increments:: ...@@ -124,21 +124,21 @@ currently computed increments::
u_ = TestFunction(V) u_ = TestFunction(V)
Boundary conditions correspond to symmetry conditions on the bottom and left Boundary conditions correspond to symmetry conditions on the bottom and left
parts (resp. numbered 1 and 3). Loading consists of a uniform pressure on the parts (resp. numbered 1 and 3). Loading consists of a uniform pressure on the
internal boundary (numbered 4). It will be progressively increased from 0 to internal boundary (numbered 4). It will be progressively increased from 0 to
:math:`q_{lim}=\dfrac{2}{\sqrt{3}}\sigma_0\log\left(\dfrac{R_e}{R_i}\right)` :math:`q_{lim}=\dfrac{2}{\sqrt{3}}\sigma_0\log\left(\dfrac{R_e}{R_i}\right)`
which is the analytical collapse load for a perfectly-plastic material (no hardening):: which is the analytical collapse load for a perfectly-plastic material (no hardening)::
bc = [DirichletBC(V.sub(1), 0, facets, 1), DirichletBC(V.sub(0), 0, facets, 3)] bc = [DirichletBC(V.sub(1), 0, facets, 1), DirichletBC(V.sub(0), 0, facets, 3)]
n = FacetNormal(mesh) n = FacetNormal(mesh)
q_lim = float(2/sqrt(3)*ln(Re/Ri)*sig0) q_lim = float(2/sqrt(3)*ln(Re/Ri)*sig0)
loading = Expression("-q*t", q=q_lim, t=0, degree=2) loading = Expression("-q*t", q=q_lim, t=0, degree=2)
def F_ext(v): def F_ext(v):
return loading*dot(n, v)*ds(4) return loading*dot(n, v)*ds(4)
----------------------------- -----------------------------
Constitutive relation update Constitutive relation update
----------------------------- -----------------------------
...@@ -147,7 +147,7 @@ Before writing the variational form, we now define some useful functions which ...@@ -147,7 +147,7 @@ Before writing the variational form, we now define some useful functions which
will enable performing the constitutive relation update using a return mapping will enable performing the constitutive relation update using a return mapping
procedure. This step is quite classical in FEM plasticity for a von Mises criterion procedure. This step is quite classical in FEM plasticity for a von Mises criterion
with isotropic hardening and follow notations from [BON2014]_. First, the strain with isotropic hardening and follow notations from [BON2014]_. First, the strain
tensor will be represented in a 3D fashion by appending zeros on the out-of-plane tensor will be represented in a 3D fashion by appending zeros on the out-of-plane
components since, even if the problem is 2D, the plastic constitutive relation will components since, even if the problem is 2D, the plastic constitutive relation will
involve out-of-plane plastic strains. The elastic consitutive relation is also defined involve out-of-plane plastic strains. The elastic consitutive relation is also defined
and a function ``as_3D_tensor`` will enable to represent a 4 dimensional vector and a function ``as_3D_tensor`` will enable to represent a 4 dimensional vector
...@@ -166,35 +166,35 @@ containing respectively :math:`xx, yy, zz` and :math:`xy` components as a 3D ten ...@@ -166,35 +166,35 @@ containing respectively :math:`xx, yy, zz` and :math:`xy` components as a 3D ten
[0, 0, X[2]]]) [0, 0, X[2]]])
The return mapping procedure consists in finding a new stress :math:`\boldsymbol{\sigma}_{n+1}` The return mapping procedure consists in finding a new stress :math:`\boldsymbol{\sigma}_{n+1}`
and internal variable :math:`p_{n+1}` state verifying the current plasticity condition and internal variable :math:`p_{n+1}` state verifying the current plasticity condition
from a previous stress :math:`\boldsymbol{\sigma}_{n}` and internal variable :math:`p_n` state and from a previous stress :math:`\boldsymbol{\sigma}_{n}` and internal variable :math:`p_n` state and
an increment of total deformation :math:`\Delta \boldsymbol{\varepsilon}`. An elastic an increment of total deformation :math:`\Delta \boldsymbol{\varepsilon}`. An elastic
trial stress :math:`\boldsymbol{\sigma}_{\text{elas}} = \boldsymbol{\sigma}_{n} + \mathbf{C}\Delta \boldsymbol{\varepsilon}` trial stress :math:`\boldsymbol{\sigma}_{\text{elas}} = \boldsymbol{\sigma}_{n} + \mathbf{C}\Delta \boldsymbol{\varepsilon}`
is first computed. The plasticity criterion is then evaluated with the previous plastic strain is first computed. The plasticity criterion is then evaluated with the previous plastic strain
:math:`f_{\text{elas}} = \sigma^{eq}_{\text{elas}} - \sigma_0 - H p_n` where :math:`f_{\text{elas}} = \sigma^{eq}_{\text{elas}} - \sigma_0 - H p_n` where
:math:`\sigma^{eq}_{\text{elas}} = \sqrt{\frac{3}{2}\boldsymbol{s}:\boldsymbol{s}}` :math:`\sigma^{eq}_{\text{elas}} = \sqrt{\frac{3}{2}\boldsymbol{s}:\boldsymbol{s}}`
with the deviatoric elastic stress :math:`\boldsymbol{s} = \operatorname{dev}\boldsymbol{\sigma}_{\text{elas}}`. with the deviatoric elastic stress :math:`\boldsymbol{s} = \operatorname{dev}\boldsymbol{\sigma}_{\text{elas}}`.
If :math:`f_{\text{elas}} < 0`, no plasticity occurs during this time increment and If :math:`f_{\text{elas}} < 0`, no plasticity occurs during this time increment and
:math:`\Delta p,\Delta \boldsymbol{\varepsilon}^p =0`. :math:`\Delta p,\Delta \boldsymbol{\varepsilon}^p =0`.
Otherwise, plasticity Otherwise, plasticity
occurs and the increment of plastic strain is given by :math:`\Delta p = \dfrac{f_{\text{elas}}}{3\mu+H}`. occurs and the increment of plastic strain is given by :math:`\Delta p = \dfrac{f_{\text{elas}}}{3\mu+H}`.
Hence, both elastic and plastic evolution can be accounted for by defining the Hence, both elastic and plastic evolution can be accounted for by defining the
plastic strain increment as follows: plastic strain increment as follows:
.. math:: .. math::
\Delta p = \dfrac{\langle f_{\text{elas}}\rangle_+}{3\mu+H} \Delta p = \dfrac{\langle f_{\text{elas}}\rangle_+}{3\mu+H}
where :math:`\langle \star \rangle_+` represents the positive part of :math:`\star` where :math:`\langle \star \rangle_+` represents the positive part of :math:`\star`
and is obtained by function ``ppos``. Plastic evolution also requires the computation and is obtained by function ``ppos``. Plastic evolution also requires the computation
of the normal vector to the final yield surface given by of the normal vector to the final yield surface given by
:math:`\boldsymbol{n}_{\text{elas}} = \boldsymbol{s}/\sigma_{\text{elas}}^{eq}`. In the following, :math:`\boldsymbol{n}_{\text{elas}} = \boldsymbol{s}/\sigma_{\text{elas}}^{eq}`. In the following,
this vector must be zero in case of elastic evolution. Hence, we multiply it by this vector must be zero in case of elastic evolution. Hence, we multiply it by
:math:`\dfrac{\langle f_{\text{elas}}\rangle_+}{ f_{\text{elas}}}` to tackle :math:`\dfrac{\langle f_{\text{elas}}\rangle_+}{ f_{\text{elas}}}` to tackle
both cases in a single expression. The final stress state is corrected by the both cases in a single expression. The final stress state is corrected by the
plastic strain as follows :math:`\boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}_{\text{elas}} - plastic strain as follows :math:`\boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}_{\text{elas}} -
\beta \boldsymbol{s}` with :math:`\beta = \dfrac{3\mu}{\sigma_{\text{elas}}^{eq}}\Delta p`. \beta \boldsymbol{s}` with :math:`\beta = \dfrac{3\mu}{\sigma_{\text{elas}}^{eq}}\Delta p`.
It can be observed that the last term vanishes in case of elastic evolution so It can be observed that the last term vanishes in case of elastic evolution so
that the final stress is indeed the elastic predictor. :: that the final stress is indeed the elastic predictor. ::
ppos = lambda x: (x+abs(x))/2. ppos = lambda x: (x+abs(x))/2.
...@@ -202,7 +202,7 @@ that the final stress is indeed the elastic predictor. :: ...@@ -202,7 +202,7 @@ that the final stress is indeed the elastic predictor. ::
sig_n = as_3D_tensor(old_sig) sig_n = as_3D_tensor(old_sig)
sig_elas = sig_n + sigma(deps) sig_elas = sig_n + sigma(deps)
s = dev(sig_elas) s = dev(sig_elas)
sig_eq = sqrt(3/2.*inner(s, s)) sig_eq = sqrt(3/2.*inner(s, s))
f_elas = sig_eq - sig0 - H*old_p f_elas = sig_eq - sig0 - H*old_p
dp = ppos(f_elas)/(3*mu+H) dp = ppos(f_elas)/(3*mu+H)
n_elas = s/sig_eq*ppos(f_elas)/f_elas n_elas = s/sig_eq*ppos(f_elas)/f_elas
...@@ -212,7 +212,7 @@ that the final stress is indeed the elastic predictor. :: ...@@ -212,7 +212,7 @@ that the final stress is indeed the elastic predictor. ::
as_vector([n_elas[0, 0], n_elas[1, 1], n_elas[2, 2], n_elas[0, 1]]), \ as_vector([n_elas[0, 0], n_elas[1, 1], n_elas[2, 2], n_elas[0, 1]]), \
beta, dp beta, dp
.. note:: .. note::
We could have used conditionals to write more explicitly the difference We could have used conditionals to write more explicitly the difference
between elastic and plastic evolution. between elastic and plastic evolution.
...@@ -224,10 +224,10 @@ need to derive the algorithmic consistent tangent matrix given by: ...@@ -224,10 +224,10 @@ need to derive the algorithmic consistent tangent matrix given by:
\boldsymbol{n}_{\text{elas}} \otimes \boldsymbol{n}_{\text{elas}} - 2\mu\beta\mathbf{DEV} \boldsymbol{n}_{\text{elas}} \otimes \boldsymbol{n}_{\text{elas}} - 2\mu\beta\mathbf{DEV}
where :math:`\mathbf{DEV}` is the 4th-order tensor associated with the deviatoric where :math:`\mathbf{DEV}` is the 4th-order tensor associated with the deviatoric
operator (note that :math:`\mathbf{C}_{\text{tang}}^{\text{alg}}=\mathbf{C}` for operator (note that :math:`\mathbf{C}_{\text{tang}}^{\text{alg}}=\mathbf{C}` for
elastic evolution). Contrary to what is done in the FEniCS book, we do not store it as the components elastic evolution). Contrary to what is done in the FEniCS book, we do not store it as the components
of a 4th-order tensor but it will suffice keeping track of the normal vector and of a 4th-order tensor but it will suffice keeping track of the normal vector and
the :math:`\beta` parameter related to the plastic strains. We instead define a function the :math:`\beta` parameter related to the plastic strains. We instead define a function
computing the tangent stress :math:`\boldsymbol{\sigma}_\text{tang} = \mathbf{C}_{\text{tang}}^{\text{alg}} computing the tangent stress :math:`\boldsymbol{\sigma}_\text{tang} = \mathbf{C}_{\text{tang}}^{\text{alg}}
\boldsymbol{\varepsilon}` as follows:: \boldsymbol{\varepsilon}` as follows::
...@@ -242,26 +242,26 @@ Global problem and Newton-Raphson procedure ...@@ -242,26 +242,26 @@ Global problem and Newton-Raphson procedure
We are now in position to derive the global problem with its associated We are now in position to derive the global problem with its associated
Newton-Raphson procedure. Each iteration will require establishing equilibrium