[Getfem-commits] r4803 - /trunk/getfem/interface/tests/matlab/demo_plasticity.m

[Yves . Renard]

Author: renard
Date: Sat Nov  8 09:06:33 2014
New Revision: 4803

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

Modified:
    trunk/getfem/interface/tests/matlab/demo_plasticity.m

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=4803&r1=4802&r2=4803&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_plasticity.m       (original)
+++ trunk/getfem/interface/tests/matlab/demo_plasticity.m       Sat Nov  8 
09:06:33 2014
@@ -27,8 +27,10 @@
 
 %%
 
-with_hardening = 1;
+with_hardening = 0;
+bi_material = true;
 test_tangent_matrix = 0;
+do_plot = true;
 
 % Initialize used data
 L = 100;
@@ -39,8 +41,9 @@
 f = [0 -330]';
 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];
 if (with_hardening == 1)
-  f = [20000 0]';
-  t = [0 0.5 0.9 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];
+  f = [0 -330]';
+  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];
+  % t = [0 0.5 0.9 1.2 1.5 1.8 1.5 1.2 0.9 0.5 0.3 0];
 end
 
 % Create the mesh
@@ -55,14 +58,13 @@
 
 % Define used MeshFem
 if (with_hardening == 1)
-  mf_u=gfMeshFem(m,2); set(mf_u, 'fem',gfFem('FEM_PK(2,1)'));
+  mf_u=gfMeshFem(m,2); set(mf_u, 'fem',gfFem('FEM_PK(2,2)'));
 else
   mf_u=gfMeshFem(m,2); set(mf_u, 'fem',gfFem('FEM_PK(2,1)'));
 end
 mf_data=gfMeshFem(m); set(mf_data, 'fem', gfFem('FEM_PK_DISCONTINUOUS(2,0)'));
 % mf_sigma=gfMeshFem(m,4); set(mf_sigma, 
'fem',gfFem('FEM_PK_DISCONTINUOUS(2,1)'));
 mf_sigma=gfMeshFem(m,4); set(mf_sigma, 
'fem',gfFem('FEM_PK_DISCONTINUOUS(2,0)'));
-mf_err=gfMeshFem(m); set(mf_err, 'fem',gfFem('FEM_PK(2,0)'));
 mf_vm = gfMeshFem(m); set(mf_vm, 'fem', gfFem('FEM_PK_DISCONTINUOUS(2,1)'));
 
 % Find the border of the domain
@@ -71,31 +73,28 @@
 pidright=find(abs(P(1,:) - L)<1e-6); % Retrieve index of points which x near 
to L
 fleft =get(m,'faces from pid',pidleft); 
 fright=get(m,'faces from pid',pidright);
+set(m,'boundary',1,fleft); % for Dirichlet condition
+set(m,'boundary',2,fright); % for Neumann condition
 
 % Decomposed the mesh into 2 regions with different values of Lamé coeff
-CV = get(m, 'cvid');
-CVbottom = find(CV <= 250); % Retrieve index of convex located at the bottom
-CVtop = find(CV > 250); % Retrieve index of convex located at the top
+if (bi_material) separation = H/2; else separation = 0; end
+pidtop    = find(P(2,:)>=separation-1E-6); % Retrieve index of points of the 
top part
+pidbottom = find(P(2,:)<=separation+1E-6); % Retrieve index of points of the 
bottom part
+cvidtop   = get(m, 'cvid from pid', pidtop);
+cvidbottom= get(m, 'cvid from pid', pidbottom);
+CVtop     = sort(get(mf_data, 'basic dof from cvid', cvidtop));
+CVbottom  = sort(get(mf_data, 'basic dof from cvid', cvidbottom));
 
 % Definition of Lame coeff
 lambda(CVbottom,1) = 121150; % Steel
 lambda(CVtop,1) = 84605; % Iron
-%lambda = 121150;
 mu(CVbottom,1) = 80769; %Steel
 mu(CVtop,1) = 77839; % Iron
-%mu = 80769;
+% Definition of plastic threshold
 von_mises_threshold(CVbottom) = 7000;
 von_mises_threshold(CVtop) = 8000;
-
-if (with_hardening == 1)
-    lambda(CVtop,1) = 121150;
-    mu(CVtop,1) = 80769;
-    von_mises_threshold(CVtop) = 7000;
-end;
-
-% Assign boundary numbers
-set(m,'boundary',1,fleft); % for Dirichlet condition
-set(m,'boundary',2,fright); % for Neumann condition
+% Definition of hardening parameter
+H = mu(1)/5;
 
 % Create the model
 md = gfModel('real');
@@ -112,7 +111,6 @@
 
 % Declare that von_mises_threshold is a data of the system on mf_data
 set(md, 'add initialized fem data', 'von_mises_threshold', mf_data, 
von_mises_threshold);
-
 
 
 if (with_hardening)
@@ -121,20 +119,17 @@
   mim_data = gf_mesh_im_data(mim, -1, [N, N]);
   gf_model_set(md, 'add im data', 'sigma', mim_data);
    
-  H = mu(1)/2; 
   set(md, 'add initialized data', 'H', [H]);
 
   Is = 'Reshape(Id(meshdim*meshdim),meshdim,meshdim,meshdim,meshdim)';
   IxI = 'Id(meshdim)@Id(meshdim)';
   coeff_long = '((lambda)*(H))/((2*(mu)+(H))*(meshdim*(lambda)+2*(mu)+(H)))';
   B_inv = sprintf('((2*(mu)/(2*(mu)+(H)))*(%s) + (%s)*(%s))', Is, coeff_long, 
IxI);
-  B = sprintf('((1+(H)/(2*(mu)))*(%s) + 
(-(lambda)*(H)/(2*(mu)*(meshdim*lambda+2*(mu))))*(%s))', Is, IxI);
+  B = sprintf('((1+(H)/(2*(mu)))*(%s) + 
(-(lambda)*(H)/(2*(mu)*(meshdim*(lambda)+2*(mu))))*(%s))', Is, IxI);
   ApH = sprintf('((2*(mu)+(H))*(%s) + (lambda)*(%s))', Is, IxI);
   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_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]');
@@ -156,22 +151,18 @@
 VM=zeros(1,get(mf_vm, 'nbdof'));
 
 
-nbstep = size(t,2);
-
-dd = get(mf_err, 'basic dof from cvid');
-
-for step=1:nbstep,
-    if step > 1
-        set(md, 'variable', 'VolumicData', get(mf_data, 
'eval',{f(1,1)*t(step);f(2,1)*t(step)}));
-    end
+
+for step=1:size(t,2),
+    disp(sprintf('step %d / %d, coeff = %g', step, size(t,2), t(step)));
+    set(md, 'variable', 'VolumicData', get(mf_data, 
'eval',{f(1,1)*t(step);f(2,1)*t(step)}));
 
     if (test_tangent_matrix)
       gf_model_get(md, 'test tangent matrix', 1E-8, 10, 0.000001);
     end;
    
     % Solve the system
-    % get(md, 'solve','lsolver', 'superlu', 'noisy', 'lsearch', 'simplest',  
'alpha min', 0.8, 'max_iter', 100, 'max_res', 1e-6);
-    get(md, 'solve', 'noisy', 'max_iter', 80);
+    get(md, 'solve','lsolver', 'superlu', 'noisy', 'lsearch', 'simplest',  
'alpha min', 0.8, 'max_iter', 100, 'max_res', 1e-6);
+    % get(md, 'solve', 'noisy', 'max_iter', 80);
 
     % Retrieve the solution U
     U = get(md, 'variable', 'u', 0);
@@ -199,32 +190,26 @@
       VM = get(md, 'compute elastoplasticity Von Mises or Tresca', 'sigma', 
mf_vm, 'Von Mises');
     end
   
-    figure(2)
-    subplot(2,1,1);
-    gf_plot(mf_vm,VM, 'deformation',U,'deformation_mf',mf_u,'refine', 4, 
'deformation_scale',1, 'disp_options', 0); % 'deformed_mesh', 'on')
-    colorbar;
-    axis([-20 120 -20 40]);
-    % caxis([0 10000]);
-    n = t(step);
-    title(['Von Mises criterion for t = ', num2str(step)]);
-  
-    %ERR=gf_compute(mf_u,U,'error estimate', mim);
-    %E=ERR; E(dd)=ERR;
-    %subplot(2,1,2);
-    %gf_plot(mf_err, E, 'mesh','on', 'refine', 1); 
-    %colorbar;
-    %title('Error estimate');
-
-    %figure(3)
-    subplot(2,1,2);
-    gf_plot(mf_vm,plast, 'deformation',U,'deformation_mf',mf_u,'refine', 4, 
'deformation_scale',1, 'disp_options', 0);  % 'deformed_mesh', 'on')
-    colorbar;
-    axis([-20 120 -20 40]);
-    % caxis([0 10000]);
-    n = t(step);
-    title(['Plastification for t = ', num2str(step)]);
+    if (do_plot)
+      figure(2)
+      subplot(2,1,1);
+      gf_plot(mf_vm,VM, 'deformation',U,'deformation_mf',mf_u,'refine', 4, 
'deformation_scale',1, 'disp_options', 0); % 'deformed_mesh', 'on')
+      colorbar;
+      axis([-20 120 -20 40]);
+      % caxis([0 10000]);
+      n = t(step);
+      title(['Von Mises criterion for t = ', num2str(step)]);
+
+      subplot(2,1,2);
+      gf_plot(mf_vm,plast, 'deformation',U,'deformation_mf',mf_u,'refine', 4, 
'deformation_scale',1, 'disp_options', 0);  % 'deformed_mesh', 'on')
+      colorbar;
+      axis([-20 120 -20 40]);
+      % caxis([0 10000]);
+      n = t(step);
+      title(['Plastification for t = ', num2str(step)]);
     
-    pause(2);
+      pause(0.1);
+    end
 
 end;