vonMises_plasticity.py.html 57.2 KB
Newer Older
Jeremy BLEYER's avatar
Jeremy BLEYER committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

<!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>Elasto-plastic analysis of a 2D von Mises material &#8212; Numerical tours of continuum mechanics using FEniCS master 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:     'master',
        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" />
Jeremy BLEYER's avatar
Jeremy BLEYER committed
27
    <link rel="next" title="Beams and plates" href="../../beams_and_plates.html" />
Jeremy BLEYER's avatar
Jeremy BLEYER committed
28
    <link rel="prev" title="Linear viscoelasticity" href="../viscoelasticity/linear_viscoelasticity.html" /> 
Jeremy BLEYER's avatar
Jeremy BLEYER committed
29 30 31 32 33 34 35
  </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</a></div>
        <div class="rel" role="navigation" aria-label="related navigation">
Jeremy BLEYER's avatar
Jeremy BLEYER committed
36
          <a href="../viscoelasticity/linear_viscoelasticity.html" title="Linear viscoelasticity"
Jeremy BLEYER's avatar
Jeremy BLEYER committed
37
             accesskey="P">previous</a> |
Jeremy BLEYER's avatar
Jeremy BLEYER committed
38
          <a href="../../beams_and_plates.html" title="Beams and plates"
Jeremy BLEYER's avatar
Jeremy BLEYER committed
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430
             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="elasto-plastic-analysis-of-a-2d-von-mises-material">
<span id="vonmisesplasticity"></span><h1>Elasto-plastic analysis of a 2D von Mises material<a class="headerlink" href="#elasto-plastic-analysis-of-a-2d-von-mises-material" title="Permalink to this headline"></a></h1>
<div class="section" id="introduction">
<h2>Introduction<a class="headerlink" href="#introduction" title="Permalink to this headline"></a></h2>
<p>This example is concerned with the incremental analysis of an elasto-plastic
von Mises material. The structure response is computed using an iterative
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,
the return mapping procedure is completely analytical (with linear isotropic
hardening). The corresponding sources can be obtained from <a class="reference download internal" href="../../_downloads/vonMises_plasticity.zip" download=""><code class="xref download docutils literal"><span class="pre">vonMises_plasticity.zip</span></code></a>.
Another implementation of von Mises plasticity can also be found at
<a class="reference external" href="https://bitbucket.org/fenics-apps/fenics-solid-mechanics">https://bitbucket.org/fenics-apps/fenics-solid-mechanics</a>.</p>
<p>We point out that the 2D nature of the problem will impose keeping
track of the out-of-plane <span class="math">\(\varepsilon_{zz}^p\)</span> plastic strain and dealing with
representations of stress/strain states including the <span class="math">\(zz\)</span> component. Note
also that are not concerned with the issue of volumetric locking induced by
incompressible plastic deformations since quadratic triangles in 2D is enough
to mitigate the locking phenomenon.</p>
<p>The plastic strain evolution during the cylinder expansion will look like this:</p>
<a class="reference internal image-reference" href="../../_images/plastic_strain.gif"><img alt="../../_images/plastic_strain.gif" src="../../_images/plastic_strain.gif" style="width: 725.6px; height: 628.0px;" /></a>
</div>
<div class="section" id="problem-position">
<h2>Problem position<a class="headerlink" href="#problem-position" title="Permalink to this headline"></a></h2>
<p>In FEniCS 2017.2, the FEniCS Form Compiler <code class="docutils literal"><span class="pre">ffc</span></code> now uses <code class="docutils literal"><span class="pre">uflacs</span></code> as a default
representation instead of the old <code class="docutils literal"><span class="pre">quadrature</span></code> representation. However, using
<code class="docutils literal"><span class="pre">Quadrature</span></code> elements generates some bugs in this representation and we therefore
need to revert to the old representation. Deprecation warning messages are also disabled.
See <a class="reference external" href="https://www.allanswered.com/post/lknbq/assemble-quadrature-representation-vs-uflacs/">this post</a>
and the corresponding <a class="reference external" href="https://bitbucket.org/fenics-project/ffc/issues/146/uflacs-generates-undefined-variable-for">issue</a>
for more information.:</p>
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="kn">from</span> <span class="nn">__future__</span> <span class="k">import</span> <span class="n">print_function</span>
<span class="kn">from</span> <span class="nn">fenics</span> <span class="k">import</span> <span class="o">*</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="n">parameters</span><span class="p">[</span><span class="s2">&quot;form_compiler&quot;</span><span class="p">][</span><span class="s2">&quot;representation&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="s1">&#39;quadrature&#39;</span>
<span class="kn">import</span> <span class="nn">warnings</span>
<span class="kn">from</span> <span class="nn">ffc.quadrature.deprecation</span> <span class="k">import</span> <span class="n">QuadratureRepresentationDeprecationWarning</span>
<span class="n">warnings</span><span class="o">.</span><span class="n">simplefilter</span><span class="p">(</span><span class="s2">&quot;once&quot;</span><span class="p">,</span> <span class="n">QuadratureRepresentationDeprecationWarning</span><span class="p">)</span>
</pre></div>
</div>
<p>The material is represented by an isotropic elasto-plastic von Mises yield condition
of uniaxial strength <span class="math">\(\sigma_0\)</span> and with isotropic hardening of modulus <span class="math">\(H\)</span>.
The yield condition is thus given by:</p>
<div class="math">
\[f(\boldsymbol{\sigma}) = \sqrt{\frac{3}{2}\boldsymbol{s}:\boldsymbol{s}} - \sigma_0 -Hp \leq 0\]</div>
<p>where <span class="math">\(p\)</span> is the cumulated equivalent plastic strain. The hardening modulus
can also be related to a tangent elastic modulus <span class="math">\(E_t = \dfrac{EH}{E+H}\)</span>.</p>
<p>The considered problem is that of a plane strain hollow cylinder of internal (resp. external)
radius <span class="math">\(R_i\)</span> (resp. <span class="math">\(R_e\)</span>) under internal uniform pressure <span class="math">\(q\)</span>.
Only a quarter of cylinder is generated using <code class="docutils literal"><span class="pre">Gmsh</span></code> and converted to <code class="docutils literal"><span class="pre">.xml</span></code> format.</p>
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="c1"># elastic parameters</span>
<span class="n">E</span> <span class="o">=</span> <span class="n">Constant</span><span class="p">(</span><span class="mf">70e3</span><span class="p">)</span>
<span class="n">nu</span> <span class="o">=</span> <span class="n">Constant</span><span class="p">(</span><span class="mf">0.3</span><span class="p">)</span>
<span class="n">lmbda</span> <span class="o">=</span> <span class="n">E</span><span class="o">*</span><span class="n">nu</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">+</span><span class="n">nu</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="mi">2</span><span class="o">*</span><span class="n">nu</span><span class="p">)</span>
<span class="n">mu</span> <span class="o">=</span> <span class="n">E</span><span class="o">/</span><span class="mf">2.</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">+</span><span class="n">nu</span><span class="p">)</span>
<span class="n">sig0</span> <span class="o">=</span> <span class="n">Constant</span><span class="p">(</span><span class="mf">250.</span><span class="p">)</span>  <span class="c1"># yield strength</span>
<span class="n">Et</span> <span class="o">=</span> <span class="n">E</span><span class="o">/</span><span class="mf">100.</span>  <span class="c1"># tangent modulus</span>
<span class="n">H</span> <span class="o">=</span> <span class="n">E</span><span class="o">*</span><span class="n">Et</span><span class="o">/</span><span class="p">(</span><span class="n">E</span><span class="o">-</span><span class="n">Et</span><span class="p">)</span>  <span class="c1"># hardening modulus</span>

<span class="n">Re</span><span class="p">,</span> <span class="n">Ri</span> <span class="o">=</span> <span class="mf">1.3</span><span class="p">,</span> <span class="mf">1.</span>   <span class="c1"># external/internal radius</span>
<span class="n">mesh</span> <span class="o">=</span> <span class="n">Mesh</span><span class="p">(</span><span class="s2">&quot;thick_cylinder.xml&quot;</span><span class="p">)</span>
<span class="n">facets</span> <span class="o">=</span> <span class="n">MeshFunction</span><span class="p">(</span><span class="s2">&quot;size_t&quot;</span><span class="p">,</span> <span class="n">mesh</span><span class="p">,</span> <span class="s2">&quot;thick_cylinder_facet_region.xml&quot;</span><span class="p">)</span>
<span class="n">ds</span> <span class="o">=</span> <span class="n">Measure</span><span class="p">(</span><span class="s1">&#39;ds&#39;</span><span class="p">)[</span><span class="n">facets</span><span class="p">]</span>
</pre></div>
</div>
<p>Function spaces will involve a standard CG space for the displacement whereas internal
state variables such as plastic strains will be represented using a <strong>Quadrature</strong> element.
This choice will make it possible to express the complex non-linear material constitutive
equation at the Gauss points only, without involving any interpolation of non-linear
expressions throughout the element. It will ensure an optimal convergence rate
for the Newton-Raphson method. See Chapter 26 of the <a class="reference external" href="https://fenicsproject.org/book/">FEniCS book</a>.
We will need Quadrature elements for 4-dimensional vectors and scalars, the number
of Gauss points will be determined by the required degree of the Quadrature element
(e.g. <code class="docutils literal"><span class="pre">degree=1</span></code> yields only 1 Gauss point, <code class="docutils literal"><span class="pre">degree=2</span></code> yields 3 Gauss points and
<code class="docutils literal"><span class="pre">degree=3</span></code> yields 6 Gauss points (note that this is suboptimal)):</p>
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">deg_u</span> <span class="o">=</span> <span class="mi">2</span>
<span class="n">deg_stress</span> <span class="o">=</span> <span class="mi">2</span>
<span class="n">V</span> <span class="o">=</span> <span class="n">VectorFunctionSpace</span><span class="p">(</span><span class="n">mesh</span><span class="p">,</span> <span class="s2">&quot;CG&quot;</span><span class="p">,</span> <span class="n">deg_u</span><span class="p">)</span>
<span class="n">We</span> <span class="o">=</span> <span class="n">VectorElement</span><span class="p">(</span><span class="s2">&quot;Quadrature&quot;</span><span class="p">,</span> <span class="n">mesh</span><span class="o">.</span><span class="n">ufl_cell</span><span class="p">(),</span> <span class="n">degree</span><span class="o">=</span><span class="n">deg_stress</span><span class="p">,</span> <span class="n">dim</span><span class="o">=</span><span class="mi">4</span><span class="p">,</span> <span class="n">quad_scheme</span><span class="o">=</span><span class="s1">&#39;default&#39;</span><span class="p">)</span>
<span class="n">W</span> <span class="o">=</span> <span class="n">FunctionSpace</span><span class="p">(</span><span class="n">mesh</span><span class="p">,</span> <span class="n">We</span><span class="p">)</span>
<span class="n">W0e</span> <span class="o">=</span> <span class="n">FiniteElement</span><span class="p">(</span><span class="s2">&quot;Quadrature&quot;</span><span class="p">,</span> <span class="n">mesh</span><span class="o">.</span><span class="n">ufl_cell</span><span class="p">(),</span> <span class="n">degree</span><span class="o">=</span><span class="n">deg_stress</span><span class="p">,</span> <span class="n">quad_scheme</span><span class="o">=</span><span class="s1">&#39;default&#39;</span><span class="p">)</span>
<span class="n">W0</span> <span class="o">=</span> <span class="n">FunctionSpace</span><span class="p">(</span><span class="n">mesh</span><span class="p">,</span> <span class="n">W0e</span><span class="p">)</span>
</pre></div>
</div>
<div class="admonition note">
<p class="first admonition-title">Note</p>
<p class="last">In older versions, it was possible to define Quadrature function spaces directly
using <code class="docutils literal"><span class="pre">FunctionSpace(mesh,</span> <span class="pre">&quot;Quadrature&quot;,</span> <span class="pre">1)</span></code>. This is no longer the case since
FEniCS 2016.1 (see <a class="reference external" href="https://bitbucket.org/fenics-project/dolfin/issues/757/functionspace-mesh-quadrature-1-broken-in">this issue</a>). Instead, Quadrature elements must first be defined
by specifying the associated degree and quadrature scheme before defining the
associated FunctionSpace.</p>
</div>
<p>Various functions are defined to keep track of the current internal state and
currently computed increments:</p>
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">sig</span> <span class="o">=</span> <span class="n">Function</span><span class="p">(</span><span class="n">W</span><span class="p">)</span>
<span class="n">sig_old</span> <span class="o">=</span> <span class="n">Function</span><span class="p">(</span><span class="n">W</span><span class="p">)</span>
<span class="n">n_elas</span> <span class="o">=</span> <span class="n">Function</span><span class="p">(</span><span class="n">W</span><span class="p">)</span>
<span class="n">beta</span> <span class="o">=</span> <span class="n">Function</span><span class="p">(</span><span class="n">W0</span><span class="p">)</span>
<span class="n">p</span> <span class="o">=</span> <span class="n">Function</span><span class="p">(</span><span class="n">W0</span><span class="p">,</span> <span class="n">name</span><span class="o">=</span><span class="s2">&quot;Cumulative plastic strain&quot;</span><span class="p">)</span>
<span class="n">u</span> <span class="o">=</span> <span class="n">Function</span><span class="p">(</span><span class="n">V</span><span class="p">,</span> <span class="n">name</span><span class="o">=</span><span class="s2">&quot;Total displacement&quot;</span><span class="p">)</span>
<span class="n">du</span> <span class="o">=</span> <span class="n">Function</span><span class="p">(</span><span class="n">V</span><span class="p">,</span> <span class="n">name</span><span class="o">=</span><span class="s2">&quot;Iteration correction&quot;</span><span class="p">)</span>
<span class="n">Du</span> <span class="o">=</span> <span class="n">Function</span><span class="p">(</span><span class="n">V</span><span class="p">,</span> <span class="n">name</span><span class="o">=</span><span class="s2">&quot;Current increment&quot;</span><span class="p">)</span>
<span class="n">v</span> <span class="o">=</span> <span class="n">TrialFunction</span><span class="p">(</span><span class="n">V</span><span class="p">)</span>
<span class="n">u_</span> <span class="o">=</span> <span class="n">TestFunction</span><span class="p">(</span><span class="n">V</span><span class="p">)</span>
</pre></div>
</div>
<p>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
internal boundary (numbered 4). It will be progressively increased from 0 to
<span class="math">\(q_{lim}=\dfrac{2}{\sqrt{3}}\sigma_0\log\left(\dfrac{R_e}{R_i}\right)\)</span>
which is the analytical collapse load for a perfectly-plastic material (no hardening):</p>
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">bc</span> <span class="o">=</span> <span class="p">[</span><span class="n">DirichletBC</span><span class="p">(</span><span class="n">V</span><span class="o">.</span><span class="n">sub</span><span class="p">(</span><span class="mi">1</span><span class="p">),</span> <span class="mi">0</span><span class="p">,</span> <span class="n">facets</span><span class="p">,</span> <span class="mi">1</span><span class="p">),</span> <span class="n">DirichletBC</span><span class="p">(</span><span class="n">V</span><span class="o">.</span><span class="n">sub</span><span class="p">(</span><span class="mi">0</span><span class="p">),</span> <span class="mi">0</span><span class="p">,</span> <span class="n">facets</span><span class="p">,</span> <span class="mi">3</span><span class="p">)]</span>


<span class="n">n</span> <span class="o">=</span> <span class="n">FacetNormal</span><span class="p">(</span><span class="n">mesh</span><span class="p">)</span>
<span class="n">q_lim</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="mi">2</span><span class="o">/</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">3</span><span class="p">)</span><span class="o">*</span><span class="n">ln</span><span class="p">(</span><span class="n">Re</span><span class="o">/</span><span class="n">Ri</span><span class="p">)</span><span class="o">*</span><span class="n">sig0</span><span class="p">)</span>
<span class="n">loading</span> <span class="o">=</span> <span class="n">Expression</span><span class="p">(</span><span class="s2">&quot;-q*t&quot;</span><span class="p">,</span> <span class="n">q</span><span class="o">=</span><span class="n">q_lim</span><span class="p">,</span> <span class="n">t</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">degree</span><span class="o">=</span><span class="mi">2</span><span class="p">)</span>

<span class="k">def</span> <span class="nf">F_ext</span><span class="p">(</span><span class="n">v</span><span class="p">):</span>
    <span class="k">return</span> <span class="n">loading</span><span class="o">*</span><span class="n">dot</span><span class="p">(</span><span class="n">n</span><span class="p">,</span> <span class="n">v</span><span class="p">)</span><span class="o">*</span><span class="n">ds</span><span class="p">(</span><span class="mi">4</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="section" id="constitutive-relation-update">
<h2>Constitutive relation update<a class="headerlink" href="#constitutive-relation-update" title="Permalink to this headline"></a></h2>
<p>Before writing the variational form, we now define some useful functions which
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
with isotropic hardening and follow notations from <a class="reference internal" href="#bon2014" id="id1">[BON2014]</a>. First, the strain
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
involve out-of-plane plastic strains. The elastic consitutive relation is also defined
and a function <code class="docutils literal"><span class="pre">as_3D_tensor</span></code> will enable to represent a 4 dimensional vector
containing respectively <span class="math">\(xx, yy, zz\)</span> and <span class="math">\(xy\)</span> components as a 3D tensor:</p>
<div class="highlight-default"><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="n">e</span> <span class="o">=</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>
    <span class="k">return</span> <span class="n">as_tensor</span><span class="p">([[</span><span class="n">e</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span> <span class="n">e</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="mi">0</span><span class="p">],</span>
                      <span class="p">[</span><span class="n">e</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="n">e</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="mi">0</span><span class="p">],</span>
                      <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">]])</span>
<span class="k">def</span> <span class="nf">sigma</span><span class="p">(</span><span class="n">eps_el</span><span class="p">):</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_el</span><span class="p">)</span><span class="o">*</span><span class="n">Identity</span><span class="p">(</span><span class="mi">3</span><span class="p">)</span> <span class="o">+</span> <span class="mi">2</span><span class="o">*</span><span class="n">mu</span><span class="o">*</span><span class="n">eps_el</span>
<span class="k">def</span> <span class="nf">as_3D_tensor</span><span class="p">(</span><span class="n">X</span><span class="p">):</span>
    <span class="k">return</span> <span class="n">as_tensor</span><span class="p">([[</span><span class="n">X</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span> <span class="n">X</span><span class="p">[</span><span class="mi">3</span><span class="p">],</span> <span class="mi">0</span><span class="p">],</span>
                      <span class="p">[</span><span class="n">X</span><span class="p">[</span><span class="mi">3</span><span class="p">],</span> <span class="n">X</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span> <span class="mi">0</span><span class="p">],</span>
                      <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="n">X</span><span class="p">[</span><span class="mi">2</span><span class="p">]]])</span>
</pre></div>
</div>
<p>The return mapping procedure consists in finding a new stress <span class="math">\(\boldsymbol{\sigma}_{n+1}\)</span>
and internal variable <span class="math">\(p_{n+1}\)</span> state verifying the current plasticity condition
from a previous stress <span class="math">\(\boldsymbol{\sigma}_{n}\)</span> and internal variable <span class="math">\(p_n\)</span> state  and
an increment of total deformation <span class="math">\(\Delta \boldsymbol{\varepsilon}\)</span>. An elastic
trial stress <span class="math">\(\boldsymbol{\sigma}_{\text{elas}} = \boldsymbol{\sigma}_{n} + \mathbf{C}\Delta \boldsymbol{\varepsilon}\)</span>
is first computed. The plasticity criterion is then evaluated with the previous plastic strain
<span class="math">\(f_{\text{elas}} = \sigma^{eq}_{\text{elas}} - \sigma_0 - H p_n\)</span> where
<span class="math">\(\sigma^{eq}_{\text{elas}} = \sqrt{\frac{3}{2}\boldsymbol{s}:\boldsymbol{s}}\)</span>
with the deviatoric elastic stress <span class="math">\(\boldsymbol{s} = \operatorname{dev}\boldsymbol{\sigma}_{\text{elas}}\)</span>.
If <span class="math">\(f_{\text{elas}} &lt; 0\)</span>, no plasticity occurs during this time increment and
<span class="math">\(\Delta p,\Delta  \boldsymbol{\varepsilon}^p =0\)</span>.</p>
<p>Otherwise, plasticity
occurs and the increment of plastic strain is given by <span class="math">\(\Delta p = \dfrac{f_{\text{elas}}}{3\mu+H}\)</span>.
Hence, both elastic and plastic evolution can be accounted for by defining the
plastic strain increment as follows:</p>
<div class="math">
\[\Delta p = \dfrac{\langle f_{\text{elas}}\rangle_+}{3\mu+H}\]</div>
<p>where <span class="math">\(\langle \star \rangle_+\)</span> represents the positive part of <span class="math">\(\star\)</span>
and is obtained by function <code class="docutils literal"><span class="pre">ppos</span></code>. Plastic evolution also requires the computation
of the normal vector to the final yield surface given by
<span class="math">\(\boldsymbol{n}_{\text{elas}} = \boldsymbol{s}/\sigma_{\text{elas}}^{eq}\)</span>. In the following,
this vector must be zero in case of elastic evolution. Hence, we multiply it by
<span class="math">\(\dfrac{\langle f_{\text{elas}}\rangle_+}{ f_{\text{elas}}}\)</span> to tackle
both cases in a single expression. The final stress state is corrected by the
plastic strain as follows <span class="math">\(\boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}_{\text{elas}} -
\beta \boldsymbol{s}\)</span> with <span class="math">\(\beta = \dfrac{3\mu}{\sigma_{\text{elas}}^{eq}}\Delta p\)</span>.
It can be observed that the last term vanishes in case of elastic evolution so
that the final stress is indeed the elastic predictor.</p>
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">ppos</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="p">(</span><span class="n">x</span><span class="o">+</span><span class="nb">abs</span><span class="p">(</span><span class="n">x</span><span class="p">))</span><span class="o">/</span><span class="mf">2.</span>
<span class="k">def</span> <span class="nf">proj_sig</span><span class="p">(</span><span class="n">deps</span><span class="p">,</span> <span class="n">old_sig</span><span class="p">,</span> <span class="n">old_p</span><span class="p">):</span>
    <span class="n">sig_n</span> <span class="o">=</span> <span class="n">as_3D_tensor</span><span class="p">(</span><span class="n">old_sig</span><span class="p">)</span>
    <span class="n">sig_elas</span> <span class="o">=</span> <span class="n">sig_n</span> <span class="o">+</span> <span class="n">sigma</span><span class="p">(</span><span class="n">deps</span><span class="p">)</span>
    <span class="n">s</span> <span class="o">=</span> <span class="n">dev</span><span class="p">(</span><span class="n">sig_elas</span><span class="p">)</span>
    <span class="n">sig_eq</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="mi">3</span><span class="o">/</span><span class="mf">2.</span><span class="o">*</span><span class="n">inner</span><span class="p">(</span><span class="n">s</span><span class="p">,</span> <span class="n">s</span><span class="p">))</span>
    <span class="n">f_elas</span> <span class="o">=</span> <span class="n">sig_eq</span> <span class="o">-</span> <span class="n">sig0</span> <span class="o">-</span> <span class="n">H</span><span class="o">*</span><span class="n">old_p</span>
    <span class="n">dp</span> <span class="o">=</span> <span class="n">ppos</span><span class="p">(</span><span class="n">f_elas</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="mi">3</span><span class="o">*</span><span class="n">mu</span><span class="o">+</span><span class="n">H</span><span class="p">)</span>
    <span class="n">n_elas</span> <span class="o">=</span> <span class="n">s</span><span class="o">/</span><span class="n">sig_eq</span><span class="o">*</span><span class="n">ppos</span><span class="p">(</span><span class="n">f_elas</span><span class="p">)</span><span class="o">/</span><span class="n">f_elas</span>
    <span class="n">beta</span> <span class="o">=</span> <span class="mi">3</span><span class="o">*</span><span class="n">mu</span><span class="o">*</span><span class="n">dp</span><span class="o">/</span><span class="n">sig_eq</span>
    <span class="n">new_sig</span> <span class="o">=</span> <span class="n">sig_elas</span><span class="o">-</span><span class="n">beta</span><span class="o">*</span><span class="n">s</span>
    <span class="k">return</span> <span class="n">as_vector</span><span class="p">([</span><span class="n">new_sig</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span> <span class="n">new_sig</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="n">new_sig</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">],</span> <span class="n">new_sig</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">]]),</span> \
           <span class="n">as_vector</span><span class="p">([</span><span class="n">n_elas</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span> <span class="n">n_elas</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="n">n_elas</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">],</span> <span class="n">n_elas</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">]]),</span> \
           <span class="n">beta</span><span class="p">,</span> <span class="n">dp</span>
</pre></div>
</div>
<div class="admonition note">
<p class="first admonition-title">Note</p>
<p class="last">We could have used conditionals to write more explicitly the difference
between elastic and plastic evolution.</p>
</div>
<p>In order to use a Newton-Raphson procedure to resolve global equilibrium, we also
need to derive the algorithmic consistent tangent matrix given by:</p>
<div class="math">
\[\mathbf{C}_{\text{tang}}^{\text{alg}} = \mathbf{C} - 3\mu\left(\dfrac{3\mu}{3\mu+H}-\beta\right)
\boldsymbol{n}_{\text{elas}} \otimes \boldsymbol{n}_{\text{elas}} - 2\mu\beta\mathbf{DEV}\]</div>
<p>where <span class="math">\(\mathbf{DEV}\)</span> is the 4th-order tensor associated with the deviatoric
operator (note that <span class="math">\(\mathbf{C}_{\text{tang}}^{\text{alg}}=\mathbf{C}\)</span> for
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
the <span class="math">\(\beta\)</span> parameter related to the plastic strains. We instead define a function
computing the tangent stress <span class="math">\(\boldsymbol{\sigma}_\text{tang} = \mathbf{C}_{\text{tang}}^{\text{alg}}
\boldsymbol{\varepsilon}\)</span> as follows:</p>
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="k">def</span> <span class="nf">sigma_tang</span><span class="p">(</span><span class="n">e</span><span class="p">):</span>
    <span class="n">N_elas</span> <span class="o">=</span> <span class="n">as_3D_tensor</span><span class="p">(</span><span class="n">n_elas</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">sigma</span><span class="p">(</span><span class="n">e</span><span class="p">)</span> <span class="o">-</span> <span class="mi">3</span><span class="o">*</span><span class="n">mu</span><span class="o">*</span><span class="p">(</span><span class="mi">3</span><span class="o">*</span><span class="n">mu</span><span class="o">/</span><span class="p">(</span><span class="mi">3</span><span class="o">*</span><span class="n">mu</span><span class="o">+</span><span class="n">H</span><span class="p">)</span><span class="o">-</span><span class="n">beta</span><span class="p">)</span><span class="o">*</span><span class="n">inner</span><span class="p">(</span><span class="n">N_elas</span><span class="p">,</span> <span class="n">e</span><span class="p">)</span><span class="o">*</span><span class="n">N_elas</span><span class="o">-</span><span class="mi">2</span><span class="o">*</span><span class="n">mu</span><span class="o">*</span><span class="n">beta</span><span class="o">*</span><span class="n">dev</span><span class="p">(</span><span class="n">e</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="section" id="global-problem-and-newton-raphson-procedure">
<h2>Global problem and Newton-Raphson procedure<a class="headerlink" href="#global-problem-and-newton-raphson-procedure" title="Permalink to this headline"></a></h2>
<p>We are now in position to derive the global problem with its associated
Newton-Raphson procedure. Each iteration will require establishing equilibrium
by driving to zero the residual between the internal forces associated with the current
stress state <code class="docutils literal"><span class="pre">sig</span></code> and the external force vector. Because we use Quadrature
elements a custom integration measure must be defined to match the quadrature
degree and scheme used by the Quadrature elements:</p>
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">metadata</span> <span class="o">=</span> <span class="p">{</span><span class="s2">&quot;quadrature_degree&quot;</span><span class="p">:</span> <span class="n">deg_stress</span><span class="p">,</span> <span class="s2">&quot;quadrature_scheme&quot;</span><span class="p">:</span> <span class="s2">&quot;default&quot;</span><span class="p">}</span>
<span class="n">dxm</span> <span class="o">=</span> <span class="n">dx</span><span class="p">(</span><span class="n">metadata</span><span class="o">=</span><span class="n">metadata</span><span class="p">)</span>

<span class="n">a_Newton</span> <span class="o">=</span> <span class="n">inner</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="n">sigma_tang</span><span class="p">(</span><span class="n">eps</span><span class="p">(</span><span class="n">u_</span><span class="p">)))</span><span class="o">*</span><span class="n">dxm</span>
<span class="n">res</span> <span class="o">=</span> <span class="o">-</span><span class="n">inner</span><span class="p">(</span><span class="n">eps</span><span class="p">(</span><span class="n">u_</span><span class="p">),</span> <span class="n">as_3D_tensor</span><span class="p">(</span><span class="n">sig</span><span class="p">))</span><span class="o">*</span><span class="n">dxm</span> <span class="o">+</span> <span class="n">F_ext</span><span class="p">(</span><span class="n">u_</span><span class="p">)</span>
</pre></div>
</div>
<p>The consitutive update defined earlier will perform nonlinear operations on
the stress and strain tensors. These nonlinear expressions must then be projected
back onto the associated Quadrature spaces. Since these fields are defined locally
in each cell (in fact only at their associated Gauss point), this projection can
be performed locally. For this reason, we define a <code class="docutils literal"><span class="pre">local_project</span></code> function
that use the <code class="docutils literal"><span class="pre">LocalSolver</span></code> to gain in efficiency (see also <a class="reference internal" href="../../tips_and_tricks.html#tipstricksprojection"><span class="std std-ref">Efficient projection on DG or Quadrature spaces</span></a>)
for more details:</p>
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="k">def</span> <span class="nf">local_project</span><span class="p">(</span><span class="n">v</span><span class="p">,</span> <span class="n">V</span><span class="p">,</span> <span class="n">u</span><span class="o">=</span><span class="kc">None</span><span class="p">):</span>
    <span class="n">dv</span> <span class="o">=</span> <span class="n">TrialFunction</span><span class="p">(</span><span class="n">V</span><span class="p">)</span>
    <span class="n">v_</span> <span class="o">=</span> <span class="n">TestFunction</span><span class="p">(</span><span class="n">V</span><span class="p">)</span>
    <span class="n">a_proj</span> <span class="o">=</span> <span class="n">inner</span><span class="p">(</span><span class="n">dv</span><span class="p">,</span> <span class="n">v_</span><span class="p">)</span><span class="o">*</span><span class="n">dxm</span>
    <span class="n">b_proj</span> <span class="o">=</span> <span class="n">inner</span><span class="p">(</span><span class="n">v</span><span class="p">,</span> <span class="n">v_</span><span class="p">)</span><span class="o">*</span><span class="n">dxm</span>
    <span class="n">solver</span> <span class="o">=</span> <span class="n">LocalSolver</span><span class="p">(</span><span class="n">a_proj</span><span class="p">,</span> <span class="n">b_proj</span><span class="p">)</span>
    <span class="n">solver</span><span class="o">.</span><span class="n">factorize</span><span class="p">()</span>
    <span class="k">if</span> <span class="n">u</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
        <span class="n">u</span> <span class="o">=</span> <span class="n">Function</span><span class="p">(</span><span class="n">V</span><span class="p">)</span>
        <span class="n">solver</span><span class="o">.</span><span class="n">solve_local_rhs</span><span class="p">(</span><span class="n">u</span><span class="p">)</span>
        <span class="k">return</span> <span class="n">u</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">solver</span><span class="o">.</span><span class="n">solve_local_rhs</span><span class="p">(</span><span class="n">u</span><span class="p">)</span>
        <span class="k">return</span>
</pre></div>
</div>
<div class="admonition note">
<p class="first admonition-title">Note</p>
<p class="last">We could have used the standard <code class="docutils literal"><span class="pre">project</span></code> if we are not interested in optimizing
the code. However, the use of Quadrature elements would have required telling
<code class="docutils literal"><span class="pre">project</span></code> to use an appropriate integration measure to solve the global <span class="math">\(L^2\)</span>
projection that occurs under the hood. This would have needed either redefining
explicitly the projection associated forms (as we just did) or specifiying the
appropriate quadrature degree to the form compiler as follows
<code class="code docutils literal"><span class="pre">project(sig_,</span> <span class="pre">W,</span> <span class="pre">form_compiler_parameters={&quot;quadrature_degree&quot;:deg_stress})</span></code></p>
</div>
<p>Before defining the Newton-Raphson loop, we set up the output file and appropriate
FunctionSpace (here piecewise constant) and Function for output of the equivalent
plastic strain since XDMF output does not handle Quadrature elements:</p>
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">file_results</span> <span class="o">=</span> <span class="n">XDMFFile</span><span class="p">(</span><span class="s2">&quot;plasticity_results.xdmf&quot;</span><span class="p">)</span>
<span class="n">file_results</span><span class="o">.</span><span class="n">parameters</span><span class="p">[</span><span class="s2">&quot;flush_output&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="kc">True</span>
<span class="n">file_results</span><span class="o">.</span><span class="n">parameters</span><span class="p">[</span><span class="s2">&quot;functions_share_mesh&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="kc">True</span>
<span class="n">P0</span> <span class="o">=</span> <span class="n">FunctionSpace</span><span class="p">(</span><span class="n">mesh</span><span class="p">,</span> <span class="s2">&quot;DG&quot;</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span>
<span class="n">p_avg</span> <span class="o">=</span> <span class="n">Function</span><span class="p">(</span><span class="n">P0</span><span class="p">,</span> <span class="n">name</span><span class="o">=</span><span class="s2">&quot;Plastic strain&quot;</span><span class="p">)</span>
</pre></div>
</div>
<p>We now define the global Newton-Raphson loop. We will discretize the applied
loading using <code class="docutils literal"><span class="pre">Nincr</span></code> increments from 0 up to 1.1 (we exclude zero from
the list of load steps). A nonlinear discretization is adopted to refine the
load steps during the plastic evolution phase. At each time increment, the
system is assembled and the residual norm is computed. The incremental displacement
<code class="docutils literal"><span class="pre">Du</span></code> is initialized to zero and the inner iteration loop performing the constitutive
update is initiated. Inside this loop, corrections <code class="docutils literal"><span class="pre">du</span></code> to the displacement
increment are computed by solving the Newton system and the return mapping
update is performed using the current total strain increment <code class="docutils literal"><span class="pre">deps</span></code>. The resulting
quantities are then projected onto their appropriate FunctionSpaces. The Newton
system and residuals are reassembled and this procedure continues until the residual
norm falls below a given tolerance. After convergence of the iteration loop,
the total displacement, stress and plastic strain states are updated</p>
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">Nitermax</span><span class="p">,</span> <span class="n">tol</span> <span class="o">=</span> <span class="mi">200</span><span class="p">,</span> <span class="mf">1e-8</span>  <span class="c1"># parameters of the Newton-Raphson procedure</span>
<span class="n">Nincr</span> <span class="o">=</span> <span class="mi">20</span>
<span class="n">load_steps</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="mf">1.1</span><span class="p">,</span> <span class="n">Nincr</span><span class="o">+</span><span class="mi">1</span><span class="p">)[</span><span class="mi">1</span><span class="p">:]</span><span class="o">**</span><span class="mf">0.5</span>
<span class="n">results</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">Nincr</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">))</span>
<span class="k">for</span> <span class="p">(</span><span class="n">i</span><span class="p">,</span> <span class="n">t</span><span class="p">)</span> <span class="ow">in</span> <span class="nb">enumerate</span><span class="p">(</span><span class="n">load_steps</span><span class="p">):</span>
    <span class="n">loading</span><span class="o">.</span><span class="n">t</span> <span class="o">=</span> <span class="n">t</span>
    <span class="n">A</span><span class="p">,</span> <span class="n">Res</span> <span class="o">=</span> <span class="n">assemble_system</span><span class="p">(</span><span class="n">a_Newton</span><span class="p">,</span> <span class="n">res</span><span class="p">,</span> <span class="n">bc</span><span class="p">)</span>
    <span class="n">nRes0</span> <span class="o">=</span> <span class="n">Res</span><span class="o">.</span><span class="n">norm</span><span class="p">(</span><span class="s2">&quot;l2&quot;</span><span class="p">)</span>
    <span class="n">nRes</span> <span class="o">=</span> <span class="n">nRes0</span>
    <span class="n">Du</span><span class="o">.</span><span class="n">interpolate</span><span class="p">(</span><span class="n">Constant</span><span class="p">((</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">)))</span>
    <span class="nb">print</span><span class="p">(</span><span class="s2">&quot;Increment:&quot;</span><span class="p">,</span> <span class="nb">str</span><span class="p">(</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">))</span>
    <span class="n">niter</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="k">while</span> <span class="n">nRes</span><span class="o">/</span><span class="n">nRes0</span> <span class="o">&gt;</span> <span class="n">tol</span> <span class="ow">and</span> <span class="n">niter</span> <span class="o">&lt;</span> <span class="n">Nitermax</span><span class="p">:</span>
        <span class="n">solve</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">du</span><span class="o">.</span><span class="n">vector</span><span class="p">(),</span> <span class="n">Res</span><span class="p">,</span> <span class="s2">&quot;mumps&quot;</span><span class="p">)</span>
        <span class="n">Du</span><span class="o">.</span><span class="n">assign</span><span class="p">(</span><span class="n">Du</span><span class="o">+</span><span class="n">du</span><span class="p">)</span>
        <span class="n">deps</span> <span class="o">=</span> <span class="n">eps</span><span class="p">(</span><span class="n">Du</span><span class="p">)</span>
        <span class="n">sig_</span><span class="p">,</span> <span class="n">n_elas_</span><span class="p">,</span> <span class="n">beta_</span><span class="p">,</span> <span class="n">dp_</span> <span class="o">=</span> <span class="n">proj_sig</span><span class="p">(</span><span class="n">deps</span><span class="p">,</span> <span class="n">sig_old</span><span class="p">,</span> <span class="n">p</span><span class="p">)</span>
        <span class="n">local_project</span><span class="p">(</span><span class="n">sig_</span><span class="p">,</span> <span class="n">W</span><span class="p">,</span> <span class="n">sig</span><span class="p">)</span>
        <span class="n">local_project</span><span class="p">(</span><span class="n">n_elas_</span><span class="p">,</span> <span class="n">W</span><span class="p">,</span> <span class="n">n_elas</span><span class="p">)</span>
        <span class="n">local_project</span><span class="p">(</span><span class="n">beta_</span><span class="p">,</span> <span class="n">W0</span><span class="p">,</span> <span class="n">beta</span><span class="p">)</span>
        <span class="n">A</span><span class="p">,</span> <span class="n">Res</span> <span class="o">=</span> <span class="n">assemble_system</span><span class="p">(</span><span class="n">a_Newton</span><span class="p">,</span> <span class="n">res</span><span class="p">,</span> <span class="n">bc</span><span class="p">)</span>
        <span class="n">nRes</span> <span class="o">=</span> <span class="n">Res</span><span class="o">.</span><span class="n">norm</span><span class="p">(</span><span class="s2">&quot;l2&quot;</span><span class="p">)</span>
        <span class="nb">print</span><span class="p">(</span><span class="s2">&quot;    Residual:&quot;</span><span class="p">,</span> <span class="n">nRes</span><span class="p">)</span>
        <span class="n">niter</span> <span class="o">+=</span> <span class="mi">1</span>
    <span class="n">u</span><span class="o">.</span><span class="n">assign</span><span class="p">(</span><span class="n">u</span><span class="o">+</span><span class="n">Du</span><span class="p">)</span>
    <span class="n">sig_old</span><span class="o">.</span><span class="n">assign</span><span class="p">(</span><span class="n">sig</span><span class="p">)</span>
    <span class="n">p</span><span class="o">.</span><span class="n">assign</span><span class="p">(</span><span class="n">p</span><span class="o">+</span><span class="n">local_project</span><span class="p">(</span><span class="n">dp_</span><span class="p">,</span> <span class="n">W0</span><span class="p">))</span>
</pre></div>
</div>
</div>
<div class="section" id="post-processing">
<h2>Post-processing<a class="headerlink" href="#post-processing" title="Permalink to this headline"></a></h2>
<p>Inside the incremental loop, the displacement and plastic strains are exported
at each time increment, the plastic strain must first be projected onto the
previously defined DG FunctionSpace. We also monitor the value of the cylinder
displacement on the inner boundary. The load-displacement curve is then plotted:</p>
<div class="highlight-default"><div class="highlight"><pre><span></span>    <span class="n">file_results</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">)</span>
    <span class="n">p_avg</span><span class="o">.</span><span class="n">assign</span><span class="p">(</span><span class="n">project</span><span class="p">(</span><span class="n">p</span><span class="p">,</span> <span class="n">P0</span><span class="p">))</span>
    <span class="n">file_results</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="n">p_avg</span><span class="p">,</span> <span class="n">t</span><span class="p">)</span>
    <span class="n">results</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span> <span class="p">:]</span> <span class="o">=</span> <span class="p">(</span><span class="n">u</span><span class="p">(</span><span class="n">Ri</span><span class="p">,</span> <span class="mi">0</span><span class="p">)[</span><span class="mi">0</span><span class="p">],</span> <span class="n">t</span><span class="p">)</span>

<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">results</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">],</span> <span class="n">results</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">],</span> <span class="s2">&quot;-o&quot;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="s2">&quot;Displacement of inner boundary&quot;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="sa">r</span><span class="s2">&quot;Applied pressure $q/q_</span><span class="si">{lim}</span><span class="s2">$&quot;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
<p>The load-displacement curve looks as follows:</p>
<a class="reference internal image-reference" href="../../_images/cylinder_expansion_load_displ.png"><img alt="../../_images/cylinder_expansion_load_displ.png" src="../../_images/cylinder_expansion_load_displ.png" style="width: 678.4px; height: 454.0px;" /></a>
<p>It can also be checked that the analytical limit load is also well reproduced
when considering a zero hardening modulus.</p>
</div>
<div class="section" id="references">
<h2>References<a class="headerlink" href="#references" title="Permalink to this headline"></a></h2>
<table class="docutils citation" frame="void" id="bon2014" rules="none">
<colgroup><col class="label" /><col /></colgroup>
<tbody valign="top">
<tr><td class="label"><a class="fn-backref" href="#id1">[BON2014]</a></td><td>Marc Bonnet, Attilio Frangi, Christian Rey.
<em>The finite element method in solid mechanics.</em> McGraw Hill Education, pp.365, 2014</td></tr>
</tbody>
</table>
</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><ul>
<li class="toctree-l2"><a class="reference internal" href="../../intro.html#what-is-it-about">What is it about ?</a></li>
<li class="toctree-l2"><a class="reference internal" href="../../intro.html#how-do-i-get-started">How do I get started ?</a></li>
<li class="toctree-l2"><a class="reference internal" href="../../intro.html#about-the-author">About the author</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="../../linear_problems.html">Linear problems in solid mechanics</a><ul>
<li class="toctree-l2"><a class="reference internal" href="../elasticity/2D_elasticity.py.html">2D linear elasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="../elasticity/orthotropic_elasticity.py.html">Orthotropic linear elasticity</a></li>
<li class="toctree-l2"><a class="reference internal" href="../elasticity/axisymmetric_elasticity.html">Axisymmetric formulation for elastic structures of revolution</a></li>
<li class="toctree-l2"><a class="reference internal" href="../thermoelasticity/thermoelasticity.html">Linear thermoelasticity (weak coupling)</a></li>
<li class="toctree-l2"><a class="reference internal" href="../thermoelasticity/thermoelasticity_transient.html">Thermo-elastic evolution problem (full coupling)</a></li>
<li class="toctree-l2"><a class="reference internal" href="../modal_analysis_dynamics/cantilever_modal.py.html">Modal analysis of an elastic structure</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="../../homogenization.html">Homogenization of heterogeneous materials</a><ul>
<li class="toctree-l2"><a class="reference internal" href="../periodic_homog_elas/periodic_homog_elas.html">Periodic homogenization of linear elastic materials</a></li>
</ul>
</li>
<li class="toctree-l1 current"><a class="reference internal" href="../../nonlinear_problems.html">Nonlinear problems in solid mechanics</a><ul class="current">
Jeremy BLEYER's avatar
Jeremy BLEYER committed
431
<li class="toctree-l2"><a class="reference internal" href="../viscoelasticity/linear_viscoelasticity.html">Linear viscoelasticity</a></li>
Jeremy BLEYER's avatar
Jeremy BLEYER committed
432 433 434
<li class="toctree-l2 current"><a class="current reference internal" href="#">Elasto-plastic analysis of a 2D von Mises material</a></li>
</ul>
</li>
Jeremy BLEYER's avatar
Jeremy BLEYER committed
435 436
<li class="toctree-l1"><a class="reference internal" href="../../beams_and_plates.html">Beams and plates</a><ul>
<li class="toctree-l2"><a class="reference internal" href="../timoshenko/beam_buckling.html">Eulerian buckling of a beam</a></li>
Jeremy BLEYER's avatar
Jeremy BLEYER committed
437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464
<li class="toctree-l2"><a class="reference internal" href="../reissner_mindlin/reissner_mindlin_quads.py.html">Reissner-Mindlin plate with Quadrilaterals</a></li>
<li class="toctree-l2"><a class="reference internal" href="../reissner_mindlin/reissner_mindlin_dg.py.html">Reissner-Mindlin plate with a Discontinuous-Galerkin approach</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="../../tips_and_tricks.html">Tips and Tricks</a><ul>
<li class="toctree-l2"><a class="reference internal" href="../../tips_and_tricks.html#efficient-projection-on-dg-or-quadrature-spaces">Efficient projection on DG or Quadrature spaces</a></li>
</ul>
</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">
Jeremy BLEYER's avatar
Jeremy BLEYER committed
465
            <a href="../viscoelasticity/linear_viscoelasticity.html" title="Linear viscoelasticity"
Jeremy BLEYER's avatar
Jeremy BLEYER committed
466
              >previous</a> |
Jeremy BLEYER's avatar
Jeremy BLEYER committed
467
            <a href="../../beams_and_plates.html" title="Beams and plates"
Jeremy BLEYER's avatar
Jeremy BLEYER committed
468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491
              >next</a> |
            <a href="../../genindex.html" title="General Index"
              >index</a>
          </div>
          <div role="note" aria-label="source link">
              <br/>
              <a href="../../_sources/demo/2D_plasticity/vonMises_plasticity.py.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>