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

[Yves . Renard]

Author: renard
Date: Thu May 26 17:55:11 2016
New Revision: 5328

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

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

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=5328&r1=5327&r2=5328&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   
Thu May 26 17:55:11 2016
@@ -33,14 +33,14 @@
 import sys
 
 NT = 128; NX = 256; option = 3; Hi = 0; Hk = 0; load_type = 1; theta = 0.5;
-LX=100.; order = 2;
+LX=100.; order = 2; do_export = 1;
 resultspath = './exported_solutions'
 
 
 def call_test_plasticity():
     return os.system("python test_small_strain_plasticity.py option"+
-        ("=%d NX=%d NT=%d Hi=%g Hk=%g theta=%g resultspath=%s load_type=%d 
order=%d"
-         % (option, NX, NT, Hi, Hk, theta, resultspath, load_type, order)))
+        ("=%d NX=%d NT=%d Hi=%g Hk=%g theta=%g resultspath=%s load_type=%d 
order=%d do_export=%d"
+         % (option, NX, NT, Hi, Hk, theta, resultspath, load_type, order, 
do_export)))
   
 
 
@@ -52,7 +52,7 @@
 # Computation of the reference solution if necessary
 refname_U  = resultspath+'/ref_perfect_plasticity_U.dat'
 refname_mf = resultspath+'/ref_perfect_plasticity_mf.mf'
-NT = 256; NX = 256; option = 4; Hi = 0; Hk = 0; load_type = 1;
+NT = 256; NX = 256; option = 3; Hi = 0; Hk = 0; load_type = 1; theta = 0.5; 
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)
@@ -61,7 +61,8 @@
   filename = resultspath+('/U_%d.dat' % (NT))
   os.system("cp %s %s" % (filename, refname_U))
 
-hrange=[LX/16, LX/22.5, LX/32, LX/45, LX/64, LX/90., LX/128]
+# hrange=[LX/16, LX/22.5, LX/32, LX/45, LX/64, LX/90., LX/128.]
+hrange=[LX/16, LX/22.5, LX/32, LX/45, LX/64]
 
 for order in [1, 2]:
 
@@ -159,7 +160,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; theta = 
0.5; LX=100.; order = 2;
+NT = 256; NX = 256; option = 3; Hi = 12000; Hk = 12000; load_type = 2; theta = 
0.5; LX=100.; order = 2; do_plot = 0;
 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.py
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/contrib/test_plasticity/test_small_strain_plasticity.py?rev=5328&r1=5327&r2=5328&view=diff
==============================================================================
--- trunk/getfem/contrib/test_plasticity/test_small_strain_plasticity.py        
(original)
+++ trunk/getfem/contrib/test_plasticity/test_small_strain_plasticity.py        
Thu May 26 17:55:11 2016
@@ -46,10 +46,11 @@
                # 2 : horizontal
 
 constraint_at_np1 = True
+trapezoidal = False  # trapezoidal or generalized mid-point  (for option 2 
only for the moment ...)
                
 bi_material = False
 test_tangent_matrix = False
-do_plot = True
+do_export = True
 
 resultspath = './exported_solutions'
 
@@ -99,13 +100,16 @@
   if (a[0:12] == 'resultspath='):
       resultspath=a[12:]; print 'resultspath set to %s from argv' % resultspath
       continue
+  if (a[0:10] == 'do_export='):
+      do_export=a[10:]; print 'do_export set to %s from argv' % do_export
+      continue
   print "Unknow argument '%s' from the command line, exiting" % a; exit(1);
 
 NY = int(np.ceil(NX * LY / (2 * LX))*2)
 DT = T/NT
 
 
-if (do_plot):
+if (do_export):
     if (not os.path.exists(resultspath)):
         os.makedirs(resultspath)
     print('You can vizualize the optimization steps by launching')
@@ -224,20 +228,26 @@
                      +')*Id(meshdim) + 2*mu*(Sym(Grad_u)-'+Epnp1+'))')
         sigma_theta = sigma_np1;
     else:
-        Etheta = '(Sym(theta*Grad_u+(1-theta)*Grad_Previous_u))'
-        Eptheta = 
'((Epn+2*mu*theta*xi*Deviator('+Etheta+'))/(1+2*mu*theta*xi))'
-        Epnp1 = '(('+Eptheta+'-(1-theta)*Epn)/theta)'
+        En = '(Sym(Grad_Previous_u))'
+        Enp1 =  '(Sym(Grad_u))'
+        Etheta = '(theta*'+Enp1+'+(1-theta)*'+En+')'
+        if (trapezoidal):
+            md.add_fem_variable('Previous_xi', mf_xi)
+            Epnp1='((Epn+(1-theta)*2*mu*Previous_xi*(Deviator('+En+')-Epn) + 
2*mu*theta*xi*Deviator('+Enp1+'))/(1+2*mu*theta*xi))'
+            Eptheta='(theta*'+Epnp1+'+(1-theta)*'+Epn+')'
+        else:
+          
Eptheta='((Epn+2*mu*theta*xi*Deviator('+Etheta+'))/(1+2*mu*theta*xi))'
+          Epnp1 = '(('+Eptheta+'-(1-theta)*Epn)/theta)'
         sigma_np1 = ('(lambda*Trace(Sym(Grad_u)-'+Epnp1
                      +')*Id(meshdim) + 2*mu*(Sym(Grad_u)-'+Epnp1+'))')
         sigma_theta = ('(lambda*Trace('+Etheta+'-'+Eptheta
                        +')*Id(meshdim) + 2*mu*('+Etheta+'-'+Eptheta+'))')
     
-    if (constraint_at_np1):
+    if (constraint_at_np1 or trapezoidal):
       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+
     #        '))*Test_xi'
@@ -433,13 +443,13 @@
    
     # Solve the system
     md.solve('noisy', 'lsearch', 'simplest',  'alpha min', 0.4, 'max_iter', 50,
-             'max_res', 1e-6)
+             'max_res', 1e-6, 'lsolver', 'mumps')
     # md.solve('noisy', 'max_iter', 80)
 
     U = md.variable('u')
 
     # Compute Von Mises and plastic part for graphical post-treatment
-    if (do_plot):
+    if (do_export):
         if (option == 1):
             sigma_np1 = 'sigma';
         else:
@@ -480,6 +490,8 @@
         NewEpn = md.interpolation(Epnp1, mim_data)
         md.set_variable('Epn', NewEpn)
         md.set_variable('Previous_u', U)
+        if (trapezoidal):
+            md.set_variable('Previous_xi', md.variable('xi'))
         
     if (option == 3 or option == 4):
         NewEpn = md.interpolation(Epnp1, mim_data)
@@ -504,7 +516,7 @@
         md.set_variable('Previous_u', U)
       
        
-    if (do_plot):
+    if (do_export):
         # Export Von Mises and plastic part
         filename = resultspath+('/von_mises_%d.vtk' % (step))
         mf_vm.export_to_vtk(filename, mf_vm, VM.reshape((len(VM))),

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=5328&r1=5327&r2=5328&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    
Thu May 26 17:55:11 2016
@@ -131,7 +131,7 @@
 
 which results from the local flow rule integration  (the pair 
:math:`(\varepsilon^p_{n+\theta}, \alpha_{n+\theta}) = ({\mathscr 
E}^p(u_{n+\theta},  \varepsilon^p_{n}, \alpha_n), {\mathscr A}(u_{n+\theta}, 
\varepsilon^p_{n}, \alpha_n))` is the solution to equations :eq:`flowrule1`, 
:eq:`flowrule2` and  :eq:`flowrule3`). Both these maps and their tangent moduli 
(usually called consistent tangent moduli) are then used in the global solve of 
the problem with a Newton method and for :math:`u_{n+1}` the unique remaining 
variable. The advantage of the return mapping strategy is that the unique 
variable of the global solve is the displacement :math:`u_{n+1}`. A nonlinear 
solve on each Gauss point is often necessary which is usualy performed with a 
local Newton method.
 
-In |gf| we propose both the return mapping trategy and also an alternative 
strategy developped below which is mainly inspired from  [PO-NI2016]_,  
[SE-PO-WO2015]_ and [HA-WO2009]_ and allow more simple tangent moduli. It 
consists in keeping (a multiple of) :math:`\Delta \gamma` as an additional 
unknown with respect to :math:`u_{n+1}`. As we will see, this will allow a more 
generic treatment of the yield functions, the price for the simplicity being 
this additional unknown scalar field.
+In |gf| we propose both the return mapping strategy and also an alternative 
strategy developped below which is mainly inspired from  [PO-NI2016]_,  
[SE-PO-WO2015]_ and [HA-WO2009]_ and allow more simple tangent moduli. It 
consists in keeping (a multiple of) :math:`\Delta \gamma` as an additional 
unknown with respect to :math:`u_{n+1}`. As we will see, this will allow a more 
generic treatment of the yield functions, the price for the simplicity being 
this additional unknown scalar field.
 
 First, we consider an additional (and optional) given function 
:math:`\alpha(\sigma_{n+\theta}, A_{n+\theta}) > 0` whose interest will appear 
later on (it will allow simple local inverses) and the new unknown scalar field
 
@@ -288,7 +288,7 @@
 
 .. math:: \dot{\varepsilon}^p \in  \partial I_K(\sigma),
 
-where :math:`K = \left\{ \sigma : \|\mbox{Dev}(\sigma)\| \le 
\sqrt{\frac{2}{3}}\sigma_y\right\}` is the set of admissible stres tensors, 
:math:`I_K` is the indicator function of this set and :math:`\partial I_K` its 
sub-differential (normal cone to :math:`K`). This can be equivalently written
+where :math:`K = \left\{ \sigma : \|\mbox{Dev}(\sigma)\| \le 
\sqrt{\frac{2}{3}}\sigma_y\right\}` is the set of admissible stress tensors, 
:math:`I_K` is the indicator function of this set and :math:`\partial I_K` its 
sub-differential (normal cone to :math:`K`). This can be equivalently written
 
 .. math:: \sigma\in K, ~~~(\tau - \sigma):\dot{\varepsilon}^p \le 0 ~~~ 
\forall \tau \in K.
 

Modified: trunk/getfem/src/dal_static_stored_objects.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/dal_static_stored_objects.cc?rev=5328&r1=5327&r2=5328&view=diff
==============================================================================
--- trunk/getfem/src/dal_static_stored_objects.cc       (original)
+++ trunk/getfem/src/dal_static_stored_objects.cc       Thu May 26 17:55:11 2016
@@ -341,8 +341,8 @@
   {
       getfem::omp_guard lock;
       GMM_NOPERATION(lock);
-      stored_object_tab& stored_objects
-        = dal::singleton<stored_object_tab>::instance();
+      // stored_object_tab& stored_objects
+      //   = dal::singleton<stored_object_tab>::instance();
 
       if (dal_static_stored_tab_valid__) {