[Getfem-commits] r4800 - in /trunk/getfem: interface/tests/matlab/ src/

[Yves . Renard]

Author: renard
Date: Wed Nov  5 17:01:57 2014
New Revision: 4800

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

Modified:
    trunk/getfem/interface/tests/matlab/demo_laplacian.m
    trunk/getfem/interface/tests/matlab/demo_plasticity.m
    trunk/getfem/src/getfem_plasticity.cc

Modified: trunk/getfem/interface/tests/matlab/demo_laplacian.m
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/demo_laplacian.m?rev=4800&r1=4799&r2=4800&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_laplacian.m        (original)
+++ trunk/getfem/interface/tests/matlab/demo_laplacian.m        Wed Nov  5 
17:01:57 2014
@@ -21,6 +21,8 @@
 gamma0 = 0.001;  % Nitsche's method parameter gamma0 (gamma = gamma0*h)
 r = 1e8;         % Penalization parameter
 draw = true;
+quadrangles = true;
+K = 2;           % Degree of the finite element method
 
 asize =  size(who('automatic_var654'));
 if (asize(1)) draw = false; end;
@@ -28,17 +30,24 @@
 % trace on;
 gf_workspace('clear all');
 NX = 20;
-m = gf_mesh('cartesian',[0:1/NX:1],[0:1/NX:1]);
-%m=gf_mesh('import','structured','GT="GT_QK(2,1)";SIZES=[1,1];NOISED=1;NSUBDIV=[1,1];')
+if (quadrangles)
+  m = gf_mesh('cartesian',[0:1/NX:1],[0:1/NX:1]);
+else
+  
m=gf_mesh('import','structured',sprintf('GT="GT_PK(2,1)";SIZES=[1,1];NOISED=0;NSUBDIV=[%d,%d];',
 NX, NX))
+end
 
 % create a mesh_fem of for a field of dimension 1 (i.e. a scalar field)
 mf = gf_mesh_fem(m,1);
-% assign the Q2 fem to all convexes of the mesh_fem,
-gf_mesh_fem_set(mf,'fem',gf_fem('FEM_QK(2,2)'));
+% assign the QK or PK fem to all convexes of the mesh_fem, and define an
+% integration method
+if (quadrangles)
+  gf_mesh_fem_set(mf,'fem',gf_fem(sprintf('FEM_QK(2,%d)', K)));
+  mim = gf_mesh_im(m, gf_integ('IM_GAUSS_PARALLELEPIPED(2,6)'));
+else
+  gf_mesh_fem_set(mf,'fem',gf_fem(sprintf('FEM_PK(2,%d)', K)));
+  mim = gf_mesh_im(m, gf_integ('IM_TRIANGLE(6)'));
+end
 
-% Integration which will be used
-mim = gf_mesh_im(m, gf_integ('IM_GAUSS_PARALLELEPIPED(2,4)'));
-%mim = gf_mesh_im(m, 
gf_integ('IM_STRUCTURED_COMPOSITE(IM_GAUSS_PARALLELEPIPED(2,5),4)'));
 % detect the border of the mesh
 border = gf_mesh_get(m,'outer faces');
 % mark it as boundary #1
@@ -49,9 +58,11 @@
 end
 
 % interpolate the exact solution
-Uexact = gf_mesh_fem_get(mf, 'eval', { 'y.*(y-1).*x.*(x-1)+x.^5' });
-% its second derivative
-F      = gf_mesh_fem_get(mf, 'eval', { '-(2*(x.^2+y.^2)-2*x-2*y+20*x.^3)' });
+% Uexact = gf_mesh_fem_get(mf, 'eval', { '10*y.*(y-1).*x.*(x-1)+10*x.^5' });
+Uexact = gf_mesh_fem_get(mf, 'eval', { 'x.*sin(2*pi*x).*sin(2*pi*y)' });
+% Opposite of its laplacian
+% F      = gf_mesh_fem_get(mf, 'eval', { 
'-(20*(x.^2+y.^2)-20*x-20*y+200*x.^3)' });
+F      = gf_mesh_fem_get(mf, 'eval', { '4*pi*(2*pi*x.*sin(2*pi*x) - 
cos(2*pi*x)).*sin(2*pi*y)' });
 
 md=gf_model('real');
 gf_model_set(md, 'add fem variable', 'u', mf);
@@ -74,18 +85,18 @@
 U = gf_model_get(md, 'variable', 'u');
 
 if (draw)
-  subplot(2,1,1); gf_plot(mf,U,'mesh','on','contour',.01:.01:.1); 
+  subplot(2,1,1); gf_plot(mf,U,'mesh','off'); 
   colorbar; title('computed solution');
 
-  subplot(2,1,2); gf_plot(mf,U-Uexact,'mesh','on'); 
+  subplot(2,1,2); gf_plot(mf,Uexact-U,'mesh','on'); 
   colorbar;title('difference with exact solution');
 end
 
-err = gf_compute(mf, U-Uexact, 'H1 norm', mim);
+err = gf_compute(mf, Uexact-U, 'H1 norm', mim);
 
 disp(sprintf('H1 norm of error: %g', err));
 
-if (err > 2E-5)
+if (err > 2E-4)
    error('Laplacian test: error to big');
 end
 

Modified: trunk/getfem/interface/tests/matlab/demo_plasticity.m
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/demo_plasticity.m?rev=4800&r1=4799&r2=4800&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_plasticity.m       (original)
+++ trunk/getfem/interface/tests/matlab/demo_plasticity.m       Wed Nov  5 
17:01:57 2014
@@ -28,14 +28,14 @@
 %%
 
 new_test = 0;
-test_tangent_matrix = 1;
+test_tangent_matrix = 0;
 
 % Initialize used data
 L = 100;
 H = 20;
 
 % f = [0 -330]';
-f = [10000 0]';
+f = [1 0]';
 t = [0 0.9032 1 1.1 1.3 1.5 1.7 1.74 1.7 1.5 1.3 1.1 1 0.9032 0.7 0.5 0.3 0.1 
0];
 
 % Create the mesh
@@ -70,14 +70,20 @@
 CVtop = find(CV > 250); % Retrieve index of convex located at the top
 
 % Definition of Lame coeff
-lambda(CVbottom,1) = 121150; % Stell
-lambda(CVtop,1) = 84605; % Iron
+lambda(CVbottom,1) = 12.1150; % Steel
+lambda(CVtop,1) = 8.4605; % Iron
 %lambda = 121150;
-mu(CVbottom,1) = 80769; %Stell
-mu(CVtop,1) = 77839; % Iron
+mu(CVbottom,1) = 8.0769; %Steel
+mu(CVtop,1) = 7.7839; % Iron
 %mu = 80769;
-von_mises_threshold(CVbottom) = 7000;
-von_mises_threshold(CVtop) = 8000;
+von_mises_threshold(CVbottom) = 0.7000;
+von_mises_threshold(CVtop) = 0.8000;
+
+if (new_test == 1)
+    lambda(CVtop,1) = 12.1150;
+    mu(CVtop,1) = 8.0769;
+    von_mises_threshold(CVtop) = 0.7000;
+end;
 
 % Assign boundary numbers
 set(m,'boundary',1,fleft); % for Dirichlet condition
@@ -110,16 +116,18 @@
   H = mu(1)/2; 
   set(md, 'add initialized data', 'H', [H]);
 
-  coeff_long = '(2*lambda*H)/((2*mu+H)*(2*meshdim*lambda+4*mu+2*H))';
-  B_inv = 
sprintf('((2*mu/(2*mu+H))*Reshape(Id(meshdim*meshdim),meshdim,meshdim,meshdim,meshdim)
 + (%s)*(Id(meshdim)@Id(meshdim)))', coeff_long);
-  B = 
'((1+H/(2*mu))*Reshape(Id(meshdim*meshdim),meshdim,meshdim,meshdim,meshdim) + 
(-lambda*H/(2*mu*(meshdim*lambda+2*mu)))*(Id(meshdim)@Id(meshdim)))';
-  ApH = 
'((2*mu+H)*Reshape(Id(meshdim*meshdim),meshdim,meshdim,meshdim,meshdim) + 
(lambda)*(Id(meshdim)@Id(meshdim)))';
+  coeff_long = '((lambda)*(H))/((2*(mu)+(H))*(meshdim*(lambda)+2*(mu)+(H)))';
+  B_inv = 
sprintf('((2*(mu)/(2*(mu)+(H)))*Reshape(Id(meshdim*meshdim),meshdim,meshdim,meshdim,meshdim)
 + (%s)*(Id(meshdim)@Id(meshdim)))', coeff_long);
+  B = 
'((1+(H)/(2*(mu)))*Reshape(Id(meshdim*meshdim),meshdim,meshdim,meshdim,meshdim) 
+ (-(lambda)*(H)/(2*(mu)*(meshdim*lambda+2*(mu))))*(Id(meshdim)@Id(meshdim)))';
+  ApH = 
'((2*(mu)+(H))*Reshape(Id(meshdim*meshdim),meshdim,meshdim,meshdim,meshdim) + 
(lambda)*(Id(meshdim)@Id(meshdim)))';
   Enp1 = '((Grad_u+Grad_u'')/2)';
   En = '((Grad_Previous_u+Grad_Previous_u'')/2)';
-
-  expr_sigma = strcat('(', B_inv, '*(Von_Mises_projection((H*', Enp1, ')+(', 
ApH, '*(',Enp1,'+',En,')) + (', B, '*sigma), von_mises_threshold) + H*', Enp1, 
'))');
+  expr_sigma = strcat('(', B_inv, '*(Von_Mises_projection(((H)*', Enp1, ')+(', 
ApH, '*(',Enp1,'+',En,')) + (', B, '*sigma), von_mises_threshold) + H*', Enp1, 
'))');
+  
+  % expr_sigma2 = 'Von_Mises_projection(Grad_u+Grad_u'', von_mises_threshold)';
   
   gf_model_set(md, 'add nonlinear generic assembly brick', mim, 
strcat(expr_sigma, ':Grad_Test_u'));
+  % gf_model_set(md, 'add finite strain elasticity brick', mim, 'u', 
'SaintVenant Kirchhoff', '[lambda; mu]');
 else
     
   % Declare that sigma is a data of the system on mf_sigma
@@ -148,9 +156,9 @@
     end
 
     if (test_tangent_matrix)
-      gf_model_get(md, 'test tangent matrix', 1E-8, 10, 0.0001);
+      gf_model_get(md, 'test tangent matrix', 1E-8, 10, 0.00000001);
     end;
-    
+   
     
     % Solve the system
     get(md, 'solve','lsolver', 'superlu', 'lsearch', 'simplest',  'alpha min', 
0.8, 'very noisy', 'max_iter', 100, 'max_res', 1e-6);

Modified: trunk/getfem/src/getfem_plasticity.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_plasticity.cc?rev=4800&r1=4799&r2=4800&view=diff
==============================================================================
--- trunk/getfem/src/getfem_plasticity.cc       (original)
+++ trunk/getfem/src/getfem_plasticity.cc       Wed Nov  5 17:01:57 2014
@@ -885,16 +885,14 @@
       base_matrix tau(N, N), tau_D(N, N);
       gmm::copy(args[0]->as_vector(), tau.as_vector());
       scalar_type s = (*(args[1]))[0];
-      
-      
+
       scalar_type tau_m = gmm::mat_trace(tau) / scalar_type(N);
       gmm::copy(tau, tau_D);
       for (size_type i = 0; i < N; ++i) tau_D(i,i) -= tau_m;
 
       scalar_type norm_tau_D = gmm::mat_euclidean_norm(tau_D);
       
-      if (norm_tau_D != scalar_type(0))
-        gmm::scale(tau_D, std::min(norm_tau_D, s) / norm_tau_D);
+      if (norm_tau_D > s) gmm::scale(tau_D, s / norm_tau_D);
 
       for (size_type i = 0; i < N; ++i) tau_D(i,i) += tau_m;
 
@@ -918,23 +916,21 @@
 
       switch(nder) {
       case 1: 
-        if (norm_tau_D < s) {
+        if (norm_tau_D <= s) {
           gmm::clear(result.as_vector());
           for (size_type i = 0; i < N; ++i)
             for (size_type j = 0; j < N; ++j)
               result(i,j,i,j) = scalar_type(1);
         } else {
-          if (norm_tau_D != scalar_type(0)) {
-            for (size_type i = 0; i < N; ++i)
-              for (size_type j = 0; j < N; ++j)
-                for (size_type m = 0; m < N; ++m)
-                  for (size_type n = 0; n < N; ++n)
-                    result(i,j,m,n)
-                      = s * (-tau_D(i,j) * tau_D(m,n)
-                             + (i == m && j == n) ? scalar_type(1) : 
scalar_type(0)
-                             - (i == j && m == n) ? 
scalar_type(1)/scalar_type(N)
-                                                    : scalar_type(0)) / 
norm_tau_D;
-          } else gmm::clear(result.as_vector());
+          for (size_type i = 0; i < N; ++i)
+            for (size_type j = 0; j < N; ++j)
+              for (size_type m = 0; m < N; ++m)
+                for (size_type n = 0; n < N; ++n)
+                  result(i,j,m,n)
+                    = s * (-tau_D(i,j) * tau_D(m,n)
+                           + ((i == m && j == n) ? scalar_type(1) : 
scalar_type(0))
+                           - ((i == j && m == n) ? 
scalar_type(1)/scalar_type(N)
+                              : scalar_type(0))) / norm_tau_D;
           for (size_type i = 0; i < N; ++i)
             for (size_type j = 0; j < N; ++j)
               result(i,i,j,j) += scalar_type(1)/scalar_type(N);