branch: master
commit 19f83737c2b2d3049fac0253059ac5516d1d631b
Author: Yves Renard <yves.ren...@insa-lyon.fr>
Date:   Wed Feb 14 14:07:22 2018 +0100

    fix inconsistencies for small strain plasticity brick
---
 .gitignore                                        |  1 +
 interface/src/gf_model_get.cc                     |  4 +-
 interface/src/gf_model_set.cc                     |  2 +-
 interface/tests/python/demo_dynamic_contact_1D.py | 98 +++++++++++++++--------
 interface/tests/python/demo_plasticity.py         |  2 +-
 src/getfem_plasticity.cc                          | 10 +--
 6 files changed, 74 insertions(+), 43 deletions(-)

diff --git a/.gitignore b/.gitignore
index 59920d6..351184f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -114,6 +114,7 @@ Makefile
 /interface/src/scilab/src/c/libsp_get.so
 /interface/src/scilab/src/c/loader.sce
 /interface/src/scilab/unloader.sce
+/interface/tests/python/*.pyc
 /doc/sphinx/tools/
 /doc/sphinx/source/scilab/cmdref*
 /doc/sphinx/source/python/cmdref*
diff --git a/interface/src/gf_model_get.cc b/interface/src/gf_model_get.cc
index ba0d365..e73d1a5 100644
--- a/interface/src/gf_model_get.cc
+++ b/interface/src/gf_model_get.cc
@@ -913,7 +913,7 @@ void gf_model_get(getfemint::mexargs_in& m_in,
       `compute_small_strain_elastoplasticity_Von_Mises`.
       @*/
     sub_command
-      ("small strain elastoplasticity next iter", 10, 15, 0, 0,
+      ("small strain elastoplasticity next iter", 3, 15, 0, 0,
        getfem::mesh_im *mim = to_meshim_object(in.pop());
        std::string lawname = in.pop().to_string();
        filter_lawname(lawname);
@@ -986,7 +986,7 @@ void gf_model_get(getfemint::mexargs_in& m_in,
       before any call of this function.
       @*/
     sub_command         
-      ("small strain elastoplasticity Von Mises", 11, 16, 0, 0,
+      ("small strain elastoplasticity Von Mises", 4, 16, 0, 0,
        getfem::mesh_im *mim = to_meshim_object(in.pop());
        const getfem::mesh_fem *mf_vm = to_meshfem_object(in.pop());
        std::string lawname = in.pop().to_string();
diff --git a/interface/src/gf_model_set.cc b/interface/src/gf_model_set.cc
index 7e60db8..8254899 100644
--- a/interface/src/gf_model_set.cc
+++ b/interface/src/gf_model_set.cc
@@ -1968,7 +1968,7 @@ void gf_model_set(getfemint::mexargs_in& m_in,
       strain and plastic multiplier).
       @*/
     sub_command
-      ("add small strain elastoplasticity brick", 10, 15, 0, 1,
+      ("add small strain elastoplasticity brick", 3, 15, 0, 1,
        getfem::mesh_im *mim = to_meshim_object(in.pop());
        std::string lawname = in.pop().to_string();
        filter_lawname(lawname);
diff --git a/interface/tests/python/demo_dynamic_contact_1D.py 
b/interface/tests/python/demo_dynamic_contact_1D.py
index 12515b4..567f4c3 100644
--- a/interface/tests/python/demo_dynamic_contact_1D.py
+++ b/interface/tests/python/demo_dynamic_contact_1D.py
@@ -34,24 +34,25 @@ import time, os, sys
 
 # Numerical parameters
 NX = 20               # Number of elements
-T = 9                 # Simulation time
-dt = 0.001            # Time step
-u_degree = 2          # Degree of the finite element method for u
+T = 12                # Simulation time
+dt = 0.002            # Time step
+u_degree = 1          # Degree of the finite element method for u
 
-gamma0_N = 1          # Nitsche parameter gamma
-theta_N =  0          # Nitsche parameter theta
+gamma0_N = 5.         # Nitsche parameter gamma
+theta_N =  0.         # Nitsche parameter theta
 
-gamma0_P = 10          # Penalization parameter
+gamma0_P = 1.         # Penalization parameter
 
 beta = 0.             # Newmark parameter beta
 gamma = 0.5           # Newmark parameter gamma
 
 e_PS = 0.             # Restitution coefficient for Paoli-Schatzman scheme
 
-version = 3           # 0 = pure Signorini contact
+version = 4           # 0 = pure Signorini contact
                       # 1 = pure Signorini contact with Paoli-Schatzman scheme
                       # 2 = penalized contact
                       # 3 = Nitsche's method
+                      # 4 = Taylor-Flanagan method
 
 lump_mass_matrix = 0  # 0 = standard mass matrix
                       # 1 = basic lumped mass matrix
@@ -69,10 +70,19 @@ output_directory = './expe_num'
 root_filename = 'dyn1d'
 do_export_in_files = False;
 
-
 # Read optional parameters on the command line 
 for i in range(1,len(sys.argv)): exec(sys.argv[i])
 
+print "Begin experiment for",
+if    (version == 0): print "Pure Signorini contact",
+elif  (version == 1): print "Paoli-Schatzman scheme",
+elif  (version == 2): print "Penalized contact",
+elif  (version == 3): print "Nitsche's method",
+elif  (version == 4): print "Taylor-Flanagan method",
+print " in P%d, with NX = %d, dt = %g" % (u_degree,NX, dt)
+
+if (version == 4 and beta != 0): print 'Incompatibility'; exit(1)
+
 # Deduced parameters
 h = 1./NX
 TT = np.arange(0, T+dt, dt)
@@ -83,7 +93,6 @@ if (version == 3): dt_max_approx = min(dt_max_approx, 
2*h/(gamma0_N))
 print 'Approximative dt_max for CFL :', dt_max_approx
 if (beta == 0 and dt > dt_max_approx): print 'Time step too large'; exit(1)
 
-
 # Exact solution. The solution is periodic of period 3
 # Return the displacement (d=0), x derivative (d=1) or time derivative (d=2)
 def uExact(x, t, d = 0):
@@ -195,7 +204,8 @@ K_N_C = gf.Spmat('copy', K_N);  # Newmark matrix with 
effective contact
 if (version == 0):     # Pure Signorini contact
     for i in range(0, u_degree+1): K_N_C[0,i] = 0.
     K_N_C[0,0] = 1;
-elif (version == 1):   # Pure Signorini contact with Paoli-Schatzman scheme
+elif (version == 1 or version == 4):
+    # Pure Signorini contact with Paoli-Schatzman scheme
     for i in range(0, u_degree+1): K_N_C[0,i] = 0.
     K_N_C[0,0] = 1;
 elif (version == 2):   # Penalized contact
@@ -260,8 +270,22 @@ for nit in range(0, NT):
                         B0 = np.copy(B); B0[0] = -e_PS*Um1[0];
                         U1=linsolve(K_N_C, B0)
                         F=M.mult(U1)-B; s1 = -F[0]/(dt*dt)
-                        V1 += (U1-U1_0)/dt
+                        V1 += (U1-U1_0)/dt; 
                 Fc[0] = 0;
+            elif (version == 4): # Pure Signorini with Taylor-Flanagan scheme
+                if (mass_matrix_type >= 1):
+                    U1[0] = max(0., -a) / K[0,0];
+                    s1 = min(0., -a);
+                else:
+                    s1 = 0;
+                    if (U1[0] < 0):
+                        F *= 0.; F[0] = 1;
+                        F2=linsolve(M, F);
+                        s1 = V1[0]/(F2[0]*dt);
+                        F *= 0.; F[0] = -s1;
+                        F2=linsolve(M, F);
+                        U1 += dt*dt*F2;
+                        V1 += dt*F2;
             elif (version == 2): # Penalized contact
                 if (mass_matrix_type >= 1):
                     U1[0] = -a / K[0,0];
@@ -329,6 +353,7 @@ for nit in range(0, NT):
     # Compute the Energy
     MMV0 = M.mult(V0); KU0 = K.mult(U0)
     E = (np.vdot(MMV0, V0) + np.vdot(KU0, U0))*0.5
+    E_org = E
     if (version == 2): # Energy stored in the penalized contact
         E += (gamma0_P/h)*pow(min(0., U0[0]),2)/2
     if (version == 3): # Nitsche Energy
@@ -347,31 +372,36 @@ for nit in range(0, NT):
     store_UL2[nit] = gf.compute_L2_norm(mfu, U0-Uex, mim)
     store_UH1[nit] = gf.compute_H1_norm(mfu, U0-Uex, mim)
 
-    print ("Time %6f"% t), "/", T, " Energy :", store_E[nit]
+    
         
     # Draw the approximated and exact solutions
-    if (do_inter_plot and t >= tplot):
-        tplot = t + dtplot; UUex = np.copy(Xdraw)
-        plt.figure(1); plt.rc('text', usetex=True)
-        plt.subplot(311) # Displacement
-        for i in range(0,Ndraw): UUex[i] = uExact(Xdraw[i], t)
-        UU = md.interpolation('u', Xdraw, m)
-        plt.cla(); plt.axis([0.,1.,-0.6,0.6])
-        plt.plot(Xdraw, UUex, 'red'); plt.plot(Xdraw, UU, 'blue')
-        plt.ylabel('u'); plt.xlabel('x')
-        plt.subplot(312) # Velocity
-        for i in range(0,Ndraw): UUex[i] = uExact(Xdraw[i], t, 2)
-        UU = md.interpolation('v', Xdraw, m)
-        plt.cla(); plt.axis([0.,1.,-0.6,0.6])
-        plt.plot(Xdraw, UUex, 'red'); plt.plot(Xdraw, UU, 'blue')
-        plt.ylabel('v'); plt.xlabel('x')
-        plt.subplot(313) # Stress
-        for i in range(0,Ndraw): UUex[i] = uExact(Xdraw[i], t, 1)
-        UU = md.interpolation('Grad_u', Xdraw, m)
-        plt.cla(); plt.axis([0.,1.,-0.6,0.6])
-        plt.plot(Xdraw, UUex, 'red'); plt.plot(Xdraw, UU, 'blue')
-        plt.ylabel('$\partial_x u$'); plt.xlabel('x')
-        plt.pause(0.01); plt.show(0)
+    if (t >= tplot-(1e-10)):
+        tplot += dtplot;
+        print ("Time %3f"% t), "/", T,
+        print (" Energy %7f" % E), (" Mech energy %7f" % E_org)
+        
+        if (do_inter_plot):
+            UUex = np.copy(Xdraw)
+            plt.figure(1); plt.rc('text', usetex=True)
+            plt.subplot(311) # Displacement
+            for i in range(0,Ndraw): UUex[i] = uExact(Xdraw[i], t)
+            UU = md.interpolation('u', Xdraw, m)
+            plt.cla(); plt.axis([0.,1.,-0.6,0.6])
+            plt.plot(Xdraw, UUex, 'red'); plt.plot(Xdraw, UU, 'blue')
+            plt.ylabel('u'); plt.xlabel('x')
+            plt.subplot(312) # Velocity
+            for i in range(0,Ndraw): UUex[i] = uExact(Xdraw[i], t, 2)
+            UU = md.interpolation('v', Xdraw, m)
+            plt.cla(); plt.axis([0.,1.,-0.6,0.6])
+            plt.plot(Xdraw, UUex, 'red'); plt.plot(Xdraw, UU, 'blue')
+            plt.ylabel('v'); plt.xlabel('x')
+            plt.subplot(313) # Stress
+            for i in range(0,Ndraw): UUex[i] = uExact(Xdraw[i], t, 1)
+            UU = md.interpolation('Grad_u', Xdraw, m)
+            plt.cla(); plt.axis([0.,1.,-0.6,0.6])
+            plt.plot(Xdraw, UUex, 'red'); plt.plot(Xdraw, UU, 'blue')
+            plt.ylabel('$\partial_x u$'); plt.xlabel('x')
+            plt.pause(0.01); plt.show(0)
     
 
     
diff --git a/interface/tests/python/demo_plasticity.py 
b/interface/tests/python/demo_plasticity.py
index c2baa7a..2307fcc 100644
--- a/interface/tests/python/demo_plasticity.py
+++ b/interface/tests/python/demo_plasticity.py
@@ -32,7 +32,7 @@ import getfem as gf
 import numpy as np
 
 
-with_graphics=True
+with_graphics=False
 try:
     import getfem_tvtk
 except:
diff --git a/src/getfem_plasticity.cc b/src/getfem_plasticity.cc
index 820ed62..8519549 100644
--- a/src/getfem_plasticity.cc
+++ b/src/getfem_plasticity.cc
@@ -770,8 +770,8 @@ namespace getfem {
     compcond = ga_substitute
       ("((mu)*xi-pos_part((mu)*xi+100*(fbound)/(mu)))", dict);
     von_mises = ga_substitute
-      ("sqrt(3/2)*sqrt(Norm_sqr((2*(mu))*(Dev_En)-(2*(mu)+(Hk))*(Epn))"
-       "+sqr(2*(mu)*Trace(En)/3-(2*(mu)+(Hk))*Trace(Epn)))", dict);
+      ("sqrt(3/2)*sqrt(Norm_sqr((2*(mu))*(Dev_En)-(2*(mu))*(Epn))"
+       "+sqr(2*(mu)*Trace(En)/3-(2*(mu))*Trace(Epn)))", dict);
   }
 
 
@@ -842,8 +842,8 @@ namespace getfem {
     xi_np1 = ga_substitute
       ("pos_part(sqrt(3/2)*Norm(B)/(sigma_y)-1/(2*(mu)))/((theta)*(dt))", 
dict);
     von_mises = ga_substitute
-      ("sqrt(3/2)*sqrt(Norm_sqr((2*(mu))*(Dev_En)-(2*(mu)+(Hk))*(Epn))"
-       "+sqr(2*(mu)*Trace(En)/3-(2*(mu)+(Hk))*Trace(Epn)))", dict);
+      ("sqrt(3/2)*sqrt(Norm_sqr((2*(mu))*(Dev_En)-(2*(mu))*(Epn))"
+       "+sqr(2*(mu)*Trace(En)/3-(2*(mu))*Trace(Epn)))", dict);
   }
 
   // Assembly strings for isotropic elastoplasticity with Von-Mises
@@ -1385,7 +1385,7 @@ namespace getfem {
       std::string dum1, dum2, dum3, dum4, dum7;
       build_isotropic_perfect_elastoplasticity_expressions_generic
         (md, lawname, unknowns_type, varnames, params,
-         dum1, dum2, dum3, dum4, sigma_after, von_mises, dum7);
+         dum1, dum2, dum3, dum4, sigma_after, von_mises, dum7);      
     }
 
     size_type n_ep = 2; // Index of the plastic strain variable

Reply via email to