[Getfem-commits] r5323 - in /trunk/getfem: contrib/test_plasticity/ doc/sphinx/source/userdoc/

[Yves . Renard]

Author: renard
Date: Sat May  7 10:23:38 2016
New Revision: 5323

URL: http://svn.gna.org/viewcvs/getfem?rev=5323&view=rev
Log:
minor modifications

Modified:
    trunk/getfem/contrib/test_plasticity/conv_test_small_strain_plasticity.py
    trunk/getfem/contrib/test_plasticity/test_small_strain_plasticity.m
    trunk/getfem/contrib/test_plasticity/test_small_strain_plasticity.py
    trunk/getfem/doc/sphinx/source/userdoc/model_plasticity_small_strain.rst

Modified: 
trunk/getfem/contrib/test_plasticity/conv_test_small_strain_plasticity.py
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/contrib/test_plasticity/conv_test_small_strain_plasticity.py?rev=5323&r1=5322&r2=5323&view=diff
==============================================================================
--- trunk/getfem/contrib/test_plasticity/conv_test_small_strain_plasticity.py   
(original)
+++ trunk/getfem/contrib/test_plasticity/conv_test_small_strain_plasticity.py   
Sat May  7 10:23:38 2016
@@ -159,7 +159,7 @@
 # Computation of the reference solution if necessary
 refname_U  = resultspath+'/ref_hardening_plasticity_U.dat'
 refname_mf = resultspath+'/ref_hardening_plasticity_mf.mf'
-NT = 256; NX = 256; option = 4; Hi = 12000; Hk = 12000; load_type = 2;
+NT = 256; NX = 256; option = 4; Hi = 12000; Hk = 12000; load_type = 2; theta = 
0.5; LX=100.; order = 2;
 if (not(os.path.exists(refname_U)) or not(os.path.isfile(refname_U))):
   if (call_test_plasticity() != 0):
       print ('Error in the computation of the reference solution'); exit(1)

Modified: trunk/getfem/contrib/test_plasticity/test_small_strain_plasticity.m
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/contrib/test_plasticity/test_small_strain_plasticity.m?rev=5323&r1=5322&r2=5323&view=diff
==============================================================================
--- trunk/getfem/contrib/test_plasticity/test_small_strain_plasticity.m 
(original)
+++ trunk/getfem/contrib/test_plasticity/test_small_strain_plasticity.m Sat May 
 7 10:23:38 2016
@@ -64,7 +64,7 @@
 
 % Numerica parameters
 T = 10;
-NT = 200;
+NT = 40;
 LX = 100;
 LY = 20;
 NX = 40;
@@ -102,7 +102,8 @@
   mf_sigma=gfMeshFem(m,2,2); set(mf_sigma, 
'fem',gfFem('FEM_PK_DISCONTINUOUS(2,0)'));
 end
 % mf_xi = gfMeshFem(m); set(mf_xi, 'fem', gfFem('FEM_PK(2,2)'));
-mf_xi = gfMeshFem(m); set(mf_xi, 'fem', gfFem('FEM_PK_DISCONTINUOUS(2,2)'));
+mf_xi = gfMeshFem(m); set(mf_xi, 'fem', gfFem('FEM_PK_DISCONTINUOUS(2,1)'));
+mf_delta = gfMeshFem(m); set(mf_delta, 'fem', 
gfFem('FEM_PK_DISCONTINUOUS(2,1)'));
 mf_data=gfMeshFem(m); set(mf_data, 'fem', gfFem('FEM_PK_DISCONTINUOUS(2,0)'));
 mf_vm = gfMeshFem(m); set(mf_vm, 'fem', gfFem('FEM_PK_DISCONTINUOUS(2,1)'));
 
@@ -261,34 +262,33 @@
     
   case 6
     set(md, 'add fem variable', 'xi', mf_xi);
-    set(md, 'add fem variable', 'delta', mf_xi);
     set(md, 'add initialized data', 'theta', [theta]);
     set(md, 'add initialized data', 'r1', [1e-8]);
     set(md, 'add initialized data', 'r2', [1]);
     set(md, 'add im data', 'Epn', mim_data);
     set(md, 'add initialized data', 'c1', [0]); % [7.5*3]);
     set(md, 'add initialized data', 'c2', [Hk]);
-    set(md, 'add initialized data', 'c3', [0.1]); % [0.03]);
-    
-    
-    
-    if (1)
+    set(md, 'add initialized data', 'c3', [0.15]); % [0.03]);
+    
+    
+    
+    if (1) % Version with two multipliers
+        set(md, 'add fem variable', 'delta', mf_delta);
         Etheta = '(Sym(theta*Grad_u+(1-theta)*Grad_Previous_u))';
         Btheta = strcat('(Epn+theta*xi*2*mu*Deviator(',Etheta,'))');
-        Eptheta = strcat('((',Btheta,')/(1+(2*mu+c2+delta)*theta*xi))'); % 
version sans c1
-        % Eptheta = 
strcat('(',Btheta,'*pos_part(1-theta*xi*c1/(Norm(',Btheta,')+1E-6))/(1+(2*mu+c2+delta)*theta*xi))');
+        Eptheta = strcat('((',Btheta,')/(1+(2*mu+c2+delta)*theta*xi))'); % 
version without c1
+        % Eptheta = 
strcat('(',Btheta,'*pos_part(1-theta*xi*c1/(Norm(',Btheta,')+1E-6))/(1+(2*mu+c2)*theta*xi
 + theta*delta))');
         Epnp1 = strcat('((', Eptheta, ' - (1-theta)*Epn)/theta)');
         sigma_np1 = strcat('(lambda*Trace(Sym(Grad_u)-',Epnp1, ')*Id(meshdim) 
+ 2*mu*(Sym(Grad_u)-', Epnp1,'))');
         sigma_theta = strcat('(lambda*Trace(',Etheta,'-',Eptheta, 
')*Id(meshdim) + 2*mu*(',Etheta,'-', Eptheta,'))');
     
         
-        fbound = 
strcat('(Norm(2*mu*Deviator(',Etheta,')-(2*mu+c2+delta)*',Eptheta,') - 
sqrt(2/3)*von_mises_threshold)');
+        fbound = 
strcat('(Norm(2*mu*Deviator(',Etheta,')-(2*mu+c2+delta)*',Eptheta,') - 
sqrt(2/3)*von_mises_threshold)'); % version without c1
         
         
         % fbound = 
strcat('(Norm(2*mu*Deviator(',Etheta,')-(2*mu+c2+delta)*',Eptheta,'-c1*Normalized_reg(',Eptheta,',1E-6))
 - sqrt(2/3)*von_mises_threshold)');
         fbound_delta = strcat('(Norm(',Eptheta,')-c3)');
-        % fbound = 
strcat('(Norm(2*mu*Deviator(',Etheta,')-(2*mu+c2)*',Eptheta,') - 
sqrt(2/3)*von_mises_threshold)');
-        expr = strcat(sigma_np1, ':Grad_Test_u + (1/r1)*(xi - 
pos_part(xi+r1*',fbound,'))*Test_xi - (1/r2)*(delta - 
pos_part(delta+r2*',fbound_delta,'))*Test_delta');
+        expr = strcat(sigma_np1, ':Grad_Test_u + (10/r1)*(xi - 
pos_part(xi+r1*',fbound,'))*Test_xi - (100/r2)*(delta - 
pos_part(delta+r2*',fbound_delta,'))*Test_delta');
         gf_model_set(md, 'add nonlinear generic assembly brick', mim, expr);
     
     else
@@ -296,7 +296,7 @@
         Etheta = '(Sym(theta*Grad_u+(1-theta)*Grad_Previous_u))';
         Btheta = strcat('(Epn+theta*xi*2*mu*Deviator(',Etheta,'))');
         Eptheta = strcat('((',Btheta,')*min(c3/(max(Norm(',Btheta,'), c3/2)), 
pos_part(1-theta*xi*c1/(Norm(',Btheta,')+0.001))/(1+(2*mu+c2)*theta*xi)))');
-        % Eptheta = strcat('(',Btheta,'*min(c3/(max(Norm(',Btheta,'), c3/2)), 
1/(1+(2*mu+c2)*theta*xi)))');
+        Eptheta = strcat('(',Btheta,'*min(c3/(max(Norm(',Btheta,'),c3/2)), 
1/(1+(2*mu+c2)*theta*xi)))'); % version without c1
         % Eptheta = strcat('(',Btheta,'*min(c3/(Norm(',Btheta,')+1e-10), 
1/(1+(2*mu+c2)*theta*xi)))');
         Epnp1 = strcat('((', Eptheta, ' - (1-theta)*Epn)/theta)');
         sigma_np1 = strcat('(lambda*Trace(Sym(Grad_u)-',Epnp1, ')*Id(meshdim) 
+ 2*mu*(Sym(Grad_u)-', Epnp1,'))');
@@ -328,13 +328,12 @@
     end;
     
     if (option == 6)
-       set(md, 'variable', 'delta', zeros(1, get(mf_xi, 'nbdof')));
+       set(md, 'variable', 'delta', zeros(1, get(mf_delta, 'nbdof')));
        set(md, 'variable', 'xi', zeros(1, get(mf_xi, 'nbdof')));
     end
    
     % Solve the system
-    get(md, 'solve', 'noisy', 'max_iter', 50, 'lsearch', 'simplest',  'alpha 
min', 0.1, 'max_res', 1e-6);
-    % get(md, 'solve', 'noisy', 'max_iter', 80);
+    get(md, 'solve', 'noisy', 'max_iter', 50, 'lsearch', 'simplest',  'alpha 
min', 0.5, 'max_res', 1e-6);
     
     if (option == 6)
        delta = get(md, 'variable', 'delta');
@@ -428,7 +427,8 @@
       title(['Von Mises criterion for t = ', num2str(step)]);
 
       subplot(3,1,2);
-      gf_plot(mf_vm,plast, 'deformation',U,'deformation_mf',mf_u,'refine', 4, 
'deformation_scale',1, 'disp_options', 0);  % 'deformed_mesh', 'on')
+      % gf_plot(mf_vm,plast, 'deformation',U,'deformation_mf',mf_u,'refine', 
4, 'deformation_scale',1, 'disp_options', 0);  % 'deformed_mesh', 'on')
+      gf_plot(mf_vm,plast,'refine', 4, 'disp_options', 0);  % 'deformed_mesh', 
'on')
       colorbar;
       axis([-20 120 -20 40]);
       % caxis([0 10000]);

Modified: trunk/getfem/contrib/test_plasticity/test_small_strain_plasticity.py
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/contrib/test_plasticity/test_small_strain_plasticity.py?rev=5323&r1=5322&r2=5323&view=diff
==============================================================================
--- trunk/getfem/contrib/test_plasticity/test_small_strain_plasticity.py        
(original)
+++ trunk/getfem/contrib/test_plasticity/test_small_strain_plasticity.py        
Sat May  7 10:23:38 2016
@@ -44,6 +44,8 @@
             
 load_type = 1  # 1 : vertical
                # 2 : horizontal
+
+constraint_at_np1 = True
                
 bi_material = False
 test_tangent_matrix = False
@@ -71,7 +73,7 @@
 LX = 40.
 LY = 20.
 NX = 40
-theta = 1.; # Parameter for the generalized mid point scheme.
+theta = 0.5; # Parameter for the generalized mid point scheme.
 order = 2;
 
 # Arguments from the command line if any
@@ -230,7 +232,11 @@
         sigma_theta = ('(lambda*Trace('+Etheta+'-'+Eptheta
                        +')*Id(meshdim) + 2*mu*('+Etheta+'-'+Eptheta+'))')
     
-    fbound = '(Norm(Deviator('+sigma_theta+'))-sqrt(2/3)*von_mises_threshold)'
+    if (constraint_at_np1):
+      fbound = '(Norm(Deviator('+sigma_np1+'))-sqrt(2/3)*von_mises_threshold)'
+    else:
+      fbound = 
'(Norm(Deviator('+sigma_theta+'))-sqrt(2/3)*von_mises_threshold)'
+    
     # fbound = '(Norm('+Eptheta+'-Epn)-theta*xi*von_mises_threshold)'
     expr = 
sigma_np1+':Grad_Test_u+(1/r)*(xi-pos_part(xi+r*'+fbound+'))*Test_xi'
     # expr = sigma_np1+':Grad_Test_u+('+fbound+'+pos_part(-xi/r-'+fbound+
@@ -273,8 +279,12 @@
     
     # fbound = ('(Norm(Deviator('+sigma_theta+')-Hk*'+Eptheta
     #           +') - von_mises_threshold - Hi*'+alpha_theta+')')
-    fbound = ('(Norm(2*mu*Deviator('+Etheta+')-(2*mu+Hk)*'+Eptheta
-              +') - sqrt(2/3)*(von_mises_threshold + Hi*'+alpha_theta+'))')
+    if (constraint_at_np1):
+      fbound = ('(Norm(2*mu*Deviator(Sym(Grad_u))-(2*mu+Hk)*'+Epnp1
+                +') - sqrt(2/3)*(von_mises_threshold + Hi*'+alpha_np1+'))')
+    else:
+      fbound = ('(Norm(2*mu*Deviator('+Etheta+')-(2*mu+Hk)*'+Eptheta
+                +') - sqrt(2/3)*(von_mises_threshold + Hi*'+alpha_theta+'))')
     expr = (sigma_np1+':Grad_Test_u + (1/r)*(xi - pos_part(xi+r*'+fbound
             +'))*Test_xi')
     # expr = (sigma_np1+':Grad_Test_u + ('+fbound+' + pos_part(-xi/r-'+fbound

Modified: 
trunk/getfem/doc/sphinx/source/userdoc/model_plasticity_small_strain.rst
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/doc/sphinx/source/userdoc/model_plasticity_small_strain.rst?rev=5323&r1=5322&r2=5323&view=diff
==============================================================================
--- trunk/getfem/doc/sphinx/source/userdoc/model_plasticity_small_strain.rst    
(original)
+++ trunk/getfem/doc/sphinx/source/userdoc/model_plasticity_small_strain.rst    
Sat May  7 10:23:38 2016
@@ -13,9 +13,7 @@
 Small strain plasticity
 -----------------------
 
-Work in progress. Not available for the moment ...
-
-A framework for the approximation of plasticity models in |gf|.
+A framework for the approximation of plasticity models in |gf|. See in 
:file:`src/getfem_plasticity.cc` and :file:`interface/src/gf_model_set.cc` for 
the brick implementation and to extend the implementation to new plasticity 
models.
 
 
 Theoretical background
@@ -175,9 +173,6 @@
 .. math:: \ds \int_{\Omega} (f(\sigma_{n+\theta} + (-f(\sigma_{n+\theta}, 
A_{n+\theta}) - \Delta \xi/r)_+ , A_{n+\theta}) ) \lambda dx = 0,   \forall 
\lambda
 
 
-pb : need of :math:`A_{n+\theta}` 
-
-
 
 Plane strain approximation
 ==========================
@@ -246,6 +241,7 @@
 Moreover
 
 .. math:: \|\mbox{Dev}(\sigma)\| = \left(\|\bar{\sigma}\|^2 - 
\Frac{1}{3}(\mbox{tr}(\bar{\sigma}))^2\right)^{1/2}.
+   :label: plane_stress_dev
 
 Note that in the case where isochoric plastic strain is assumed, one still has
 
@@ -326,30 +322,31 @@
 
 The plane strain approximation has the same expression replacing the 3D strain 
tensors by the in-plane ones :math:`\bar{\varepsilon}^p` and  
:math:`\bar{\varepsilon}(u_{n+\theta})`.
 
-.. math:: \bar{{\mathscr E}}^p(\bar{u}_{n+\theta}, \theta \Delta \xi, 
\bar{\varepsilon}^p_{n}) = \Frac{1}{1+2\mu\theta\Delta 
\xi}(\bar{\varepsilon}^p_{n} + 2\mu\theta\Delta \xi 
\mbox{Dev}^*(\bar{\varepsilon}(\bar{u}_{n+\theta}))),
-
-where :math:`\mbox{Dev}^*(\bar{\varepsilon}) = \bar{\varepsilon} - 
\Frac{\mbox{tr}(\bar{\varepsilon})}{3} \bar{I}` is still the 3D deviator.
+.. math:: \bar{\tilde{\mathscr E}}^p(\bar{u}_{n+\theta}, \theta \Delta \xi, 
\bar{\varepsilon}^p_{n}) = \Frac{1}{1+2\mu\theta\Delta 
\xi}(\bar{\varepsilon}^p_{n} + 2\mu\theta\Delta \xi 
\overline{\mbox{Dev}}(\bar{\varepsilon}(\bar{u}_{n+\theta}))),
+
+where :math:`\overline{\mbox{Dev}}(\bar{\varepsilon}) = \bar{\varepsilon} - 
\Frac{\mbox{tr}(\bar{\varepsilon})}{3} \bar{I}` is still the 3D deviator.
 
 Moreover, for the yield condition, 
 
-.. math:: \mbox{Dev}(\sigma) = 2\mu\mbox{Dev}(\varepsilon(u) - \varepsilon^p) 
= 2\mu\left(\varepsilon(u) - \varepsilon^p - 
\Frac{\mbox{tr}(\bar{\varepsilon}(u)) - \mbox{tr}(\bar{\varepsilon}^p)}{3} 
I\right)
-
-.. math:: \begin{array}{rcl} \|\mbox{Dev}(\sigma)\| &=& 
2\mu\sqrt{\left\|\bar{\varepsilon}(u) - \bar{\varepsilon}^p - 
\Frac{\mbox{tr}(\bar{\varepsilon}(u)) - \mbox{tr}(\bar{\varepsilon}^p)}{3} 
\bar{I}\right\|^2 + \Frac{(\mbox{tr}(\bar{\varepsilon}(u)) - 
\mbox{tr}(\bar{\varepsilon}^p))^2}{9}} \\ &=& \sqrt{\left\|\bar{\sigma} - 
\Frac{3\lambda+2\mu}{6(\lambda+\mu)}\mbox{tr}(\bar{\sigma})\bar{I} \right\|^2 + 
\Frac{\mu^2}{9(\lambda+\mu)^2}\mbox{tr}(\bar{\sigma})^2 } \end{array}
+.. math:: \|\mbox{Dev}(\sigma)\|^2 = 
4\mu^2\left(\|\overline{\mbox{Dev}}\bar{\varepsilon}(u) - 
\bar{\varepsilon}^p\|^2 + \left(\Frac{\mbox{tr}(\bar{\varepsilon}(u))}{3} 
-\mbox{tr}(\bar{\varepsilon}^p) \right)^2\right)
+
+And for the closest point projection approach,
+
+.. math:: \bar{\mathscr E}^p(\bar{u}_{n+\theta}, \bar{\varepsilon}^p_{n}) = 
\bar{\varepsilon}^p_{n} + \left( 1 - 
\sqrt{\frac{2}{3}}\Frac{\sigma_y}{2\mu\|B\|}\right)_+ \bar{B}
+
+with :math:`\bar{B} = 
\overline{\mbox{Dev}}(\bar{\varepsilon}(u_{n+\theta}))-\bar{\varepsilon}^p_{n}` 
and :math:`\|B\|^2 = \|\overline{\mbox{Dev}}(\bar{\varepsilon}(u_{n+\theta})) - 
\bar{\varepsilon}^p_n\|^2 + 
\left(\Frac{\mbox{tr}(\bar{\varepsilon}(u_{n+\theta}))}{3} 
-\mbox{tr}(\bar{\varepsilon}^p_n) \right)^2`.
 
 **Plane stress approximation**
 
-For plane stress approximation, we use :eq:`plane_stress_iso` which gives
-
-.. math::  \bar{\varepsilon}^p_{n+\theta} - \bar{\varepsilon}^p_{n} = \theta 
\Delta \xi \mbox{Dev}^*(\bar{\sigma}_{n+\theta}) =  \theta \Delta \xi 
\mbox{Dev}^*(\lambda^*\mbox{tr}(\bar{\varepsilon}^e_{n+\theta})\bar{I} + 2\mu 
\bar{\varepsilon}^e_{n+\theta}) = \theta \Delta 
\xi\left(\Frac{\lambda^*-2\mu}{3}\mbox{tr}(\bar{\varepsilon}^e_{n+\theta})\bar{I}
 + 2\mu\bar{\varepsilon}^e_{n+\theta}\right)
-
-thus with :math:`\beta = \Frac{\lambda^*-2\mu}{3}` one has
-
-.. math::  (1+2\mu\theta \Delta \xi)\bar{\varepsilon}^p_{n+\theta} + 
\beta\theta \Delta \xi \mbox{tr}(\bar{\varepsilon}^p_{n+\theta})\bar{I} = 
\bar{\varepsilon}^p_{n} + \theta \Delta 
\xi\left(\beta\mbox{tr}(\bar{\varepsilon}(u_{n+\theta}))\bar{I} + 
2\mu\bar{\varepsilon}(u_{n+\theta})\right)
-
-By inverting this relation we find for :math:`A = \bar{\varepsilon}^p_{n} + 
\theta \Delta \xi\left(\beta\mbox{tr}(\bar{\varepsilon}(u_{n+\theta}))\bar{I} + 
2\mu\bar{\varepsilon}(u_{n+\theta})\right)`
-
-.. math::  \bar{{\mathscr E}}^p(\bar{u}_{n+\theta}, \theta \Delta \xi, 
\bar{\varepsilon}^p_{n}) = \Frac{1}{1+2\mu\theta\Delta \xi} A - \left( 
\Frac{\beta\theta\Delta \xi}{(1+2\mu\theta\Delta 
\xi)(1+(2\mu+2\beta)\theta\Delta \xi)} \right) \mbox{tr}(A)\bar{I}
-
+For plane stress approximation, using :eq:`plane_stress_iso` we deduce from 
the expression of the 3D case
+
+.. math::  \bar{\varepsilon}^p_{n+\theta} = \Frac{1}{1+2\mu\theta\Delta 
\xi}\left(\bar{\varepsilon}^p_{n} +2\mu\theta\Delta 
\xi\left(\bar{\varepsilon}(u_{n+\theta}) - 
\Frac{2\mu}{3(\lambda+2\mu)}(\mbox{tr}(\bar{\varepsilon}(u_{n+\theta})) - 
\mbox{tr}(\bar{\varepsilon}_{n+\theta}^p))\bar{I}\right) \right),
+
+since :math:`\mbox{Dev}(\varepsilon(u)) = \varepsilon(u) - 
\Frac{2\mu}{3(\lambda+2\mu)}(\mbox{tr}(\bar{\varepsilon}(u)) - 
\mbox{tr}(\bar{\varepsilon}^p))`. Of course, this relation still has to be 
inverted. Denoting :math:`\alpha = 1+2\mu\theta\Delta \xi`, :math:`\beta = 
\Frac{4\mu^2\theta\Delta \xi}{3\lambda+6\mu}` and :math:`C = 
\bar{\varepsilon}^p_{n} +2\mu\theta\Delta 
\xi\left(\bar{\varepsilon}(u_{n+\theta}) - 
\Frac{2\mu}{3(\lambda+2\mu)}(\mbox{tr}(\bar{\varepsilon}(u_{n+\theta}))))\bar{I}\right)`
 one obtains
+
+.. math:: \bar{\varepsilon}^p_{n+\theta} = \Frac{\beta 
\mbox{tr}(C)}{\alpha(\alpha-2\beta)}\bar{I} + \Frac{1}{\alpha}C.
+
+Moreover, for the yield condition, expression :eq:`label: plane_stress_dev` 
can be used.
 
 Isotropic elastoplasticity with linear isotropic and kinematic hardening and 
Von-Mises criterion
 
================================================================================================
@@ -411,10 +408,6 @@
 
 which complete the expression.
 
-**Plane strain approximation**
-
-
-The plane strain approximation has the same expression replacing the 3D strain 
tensors by the in-plane ones :math:`\bar{\varepsilon}^p` and  
:math:`\bar{\varepsilon}(u_{n+\theta})`.
 
 Souza-Auricchio elastoplasticity law (for shape memory alloys)
 ==============================================================
@@ -461,11 +454,6 @@
 .. or using :eq:`souza_auri_comp`
 .. .. math:: \|\tilde{\mathscr E}^p(u_{n+\theta}, \theta \Delta \xi, 
\varepsilon^p_{n}) - \varepsilon^p_{n}\| \le \theta\Delta 
\xi\sqrt{\frac{2}{3}}\sigma_{y}.
 .. stupid ? Yes a priori ! 
-
-**Plane strain approximation**
-
-
-The plane strain approximation has the same expression replacing the 3D strain 
tensors by the in-plane ones :math:`\bar{\varepsilon}^p` and  
:math:`\bar{\varepsilon}(u_{n+\theta})`.