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

[Yves . Renard]

Author: renard
Date: Thu Oct  2 16:46:03 2014
New Revision: 4768

URL: http://svn.gna.org/viewcvs/getfem?rev=4768&view=rev
Log:
Adding the projection on the plastic addmissible stresses for a Von Mises 
criterium in the generic assembly language

Added:
    trunk/getfem/interface/tests/matlab/demo_plasticity.m
      - copied, changed from r4766, 
trunk/getfem/interface/tests/matlab/test_plasticity_new_brick.m
Removed:
    trunk/getfem/interface/tests/matlab/test_plasticity_new_brick.m
Modified:
    trunk/getfem/src/getfem_plasticity.cc

Copied: trunk/getfem/interface/tests/matlab/demo_plasticity.m (from r4766, 
trunk/getfem/interface/tests/matlab/test_plasticity_new_brick.m)
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/demo_plasticity.m?p2=trunk/getfem/interface/tests/matlab/demo_plasticity.m&p1=trunk/getfem/interface/tests/matlab/test_plasticity_new_brick.m&r1=4766&r2=4768&rev=4768&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/test_plasticity_new_brick.m     
(original)
+++ trunk/getfem/interface/tests/matlab/demo_plasticity.m       Thu Oct  2 
16:46:03 2014
@@ -102,7 +102,28 @@
 set(md, 'add fem data', 'sigma', mf_sigma);
 
 % Add plasticity brick on u
-set(md, 'add elastoplasticity brick', mim, 'VM', 'u', 'lambda', 'mu', 
'von_mises_threshold', 'sigma');
+% set(md, 'add elastoplasticity brick', mim, 'VM', 'u', 'lambda', 'mu', 
'von_mises_threshold', 'sigma');
+
+new_test = 0;
+if (new_test) 
+
+  % Remain to define: Previous_u, sigma
+    
+  H = mu(CVtop)/2; 
+  set(md, 'add initialized data', 'H', [H]);
+
+  coeff_long = 'lambda*H / (..)';
+  B_inv = strcat('(2*mu/(2*mu+H))*Reshape(Id(N*N),N,N,N,N) + 
(',coeff_long,')*(Id(N)@Id(N))');
+  B = '(1+H/(2*mu))*Reshape(Id(N*N),N,N,N,N) + 
(-lambda*H/(2*mu*(3*lambda+2*mu)))*(Id(N)@Id(N))';
+  ApH = '(2*mu+H)*Reshape(Id(N*N),N,N,N,N) + (lambda)*(Id(N)@Id(N))';
+  Enp1 = '((Grad_u+Grad_u'')/2)';
+  En = '((Grad_Previous_u+Grad_Previous_u'')/2)';
+
+  assline = strcat('(', B_inv, '*(Von_Mises_projection(H*', Epn1, '+', ApH, 
'*(',Enp1,'+',En,') + ', B, '*sigma) + H*', Enp1, ')):Grad_Test_u');
+  gf_model_set(md, 'add nonlinear generic assembly brick', mim, assline);
+else
+  set(md, 'add elastoplasticity brick', mim, 'VM', 'u', 'lambda', 'mu', 
'von_mises_threshold', 'sigma');
+end
 
 % Add homogeneous Dirichlet condition to u on the left hand side of the domain
 set(md, 'add Dirichlet condition with multipliers', mim, 'u', mf_u, 1);

Removed: trunk/getfem/interface/tests/matlab/test_plasticity_new_brick.m
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/test_plasticity_new_brick.m?rev=4767&view=auto
==============================================================================
--- trunk/getfem/interface/tests/matlab/test_plasticity_new_brick.m     
(original)
+++ trunk/getfem/interface/tests/matlab/test_plasticity_new_brick.m     
(removed)
@@ -1,173 +0,0 @@
-% Copyright (C) 2010-2012 Amandine Cottaz, Yves Renard.
-%
-% This file is a part of GETFEM++
-%
-% Getfem++  is  free software;  you  can  redistribute  it  and/or modify it
-% under  the  terms  of the  GNU  Lesser General Public License as published
-% by  the  Free Software Foundation;  either version 3 of the License,  or
-% (at your option) any later version along with the GCC Runtime Library
-% Exception either version 3.1 or (at your option) any later version.
-% This program  is  distributed  in  the  hope  that it will be useful,  but
-% WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
-% or  FITNESS  FOR  A PARTICULAR PURPOSE.  See the GNU Lesser General Public
-% License and GCC Runtime Library Exception for more details.
-% You  should  have received a copy of the GNU Lesser General Public License
-% along  with  this program;  if not, write to the Free Software Foundation,
-% Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
-
-clear all
-gf_workspace('clear all');
-clc
-
-%%
-
-% We try to compute a plasticity problem with a Von Mises crierion
-% For convenience we consider an homogenous Dirichlet condition on the left
-% of the domain and an easy computed Neumann Condition on the right
-
-%%
-
-
-
-% Initialize used data
-L = 100;
-H = 20;
-
-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];
-
-% Create the mesh
-m = gfMesh('triangles grid', [0:2:L], [0:1:H]);
-
-% Plotting
-% gf_plot_mesh(m, 'vertices', 'on', 'convexes', 'on');
-
-% Define used MeshIm
-mim=gfMeshIm(m);  set(mim, 'integ', gfInteg('IM_TRIANGLE(6)')); % Gauss 
methods on triangles
-
-% Define used MeshFem
-% mf_u=gfMeshFem(m,2); set(mf_u, 'fem',gfFem('FEM_PK(2,2)'));
-mf_u=gfMeshFem(m,2); set(mf_u, 'fem',gfFem('FEM_PK(2,1)'));
-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)'));
-mf_pl = gfMeshFem(m); set(mf_pl, 'fem', gfFem('FEM_PK_DISCONTINUOUS(2,1)'));
-
-% Find the border of the domain
-P=get(m, 'pts');
-pidleft=find(abs(P(1,:))<1e-6); % Retrieve index of points which x near to 0
-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);
-
-% 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
-
-% Definition of Lame coeff
-lambda(CVbottom) = 121150; % Stell
-lambda(CVtop) = 84605; % Iron
-%lambda = 121150;
-mu(CVbottom) = 80769; %Stell
-mu(CVtop) = 77839; % Iron
-%mu = 80769;
-von_mises_threshold(CVbottom) = 7000;
-von_mises_threshold(CVtop) = 8000;
-
-% Assign boundary numbers
-set(m,'boundary',1,fleft); % for Dirichlet condition
-set(m,'boundary',2,fright); % for Neumann condition
-
-% Create the model
-md = gfModel('real');
-
-% Declare that u is the unknown of the system on mf_u
-% 2 is the number of version of the data stored, for the time integration 
scheme 
-set(md, 'add fem variable', 'u', mf_u, 2);
-
-% Declare that lambda is a data of the system on mf_data
-set(md, 'add initialized fem data', 'lambda', mf_data, lambda);
-
-% Declare that mu is a data of the system on mf_data
-set(md, 'add initialized fem data', 'mu', mf_data, mu);
-
-% 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);
-
-% Declare that sigma is a data of the system on mf_sigma
-% 2 is the number of version of the data stored, for the time integration 
scheme
-set(md, 'add fem data', 'sigma', mf_sigma);
-
-% Add plasticity brick on u
-set(md, 'add elastoplasticity brick', mim, 'VM', 'u', 'lambda', 'mu', 
'von_mises_threshold', 'sigma');
-
-% Add homogeneous Dirichlet condition to u on the left hand side of the domain
-set(md, 'add Dirichlet condition with multipliers', mim, 'u', mf_u, 1);
-
-% Add a source term to the system
-set(md,'add initialized fem data', 'VolumicData', mf_data, get(mf_data, 
'eval',{f(1,1);f(2,1)*t(1)}));
-set(md, 'add source term brick', mim, 'u', 'VolumicData', 2);
-
-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);f(2,1)*t(step)}));
-    end
-
-    % Solve the system
-    get(md, 'solve','lsolver', 'superlu', 'lsearch', 'simplest',  'alpha min', 
0.8, 'very noisy', 'max_iter', 100, 'max_res', 1e-6);
-
-    % Retrieve the solution U
-    U = get(md, 'variable', 'u', 0);
-    
-    % Compute new plasticity constraints used to compute 
-    % the Von Mises or Tresca stress
-    get(md, 'elastoplasticity next iter', mim, 'u', 'VM', 'lambda', 'mu', 
'von_mises_threshold', 'sigma');
-    plast = get(md, 'compute plastic part', mim, mf_pl, 'u', 'VM', 'lambda', 
'mu', 'von_mises_threshold', 'sigma');
-      
-    % Compute Von Mises or Tresca stress
-    VM = get(md, 'compute elastoplasticity Von Mises or Tresca', 'sigma', 
mf_vm, 'Von Mises');
-    max(abs(VM));
-  
-    figure(2)
-    subplot(2,1,1);
-    gf_plot(mf_vm,VM, 'deformation',U,'deformation_mf',mf_u,'refine', 4, 
'deformation_scale',1); % 'deformed_mesh', 'on')
-    colorbar;
-    caxis([0 10000]);
-    n = t(step);
-    title(['Von Mises criterion for t = ', num2str(n)]);
-  
-    %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)
-    gf_plot(mf_pl,plast, 'deformation',U,'deformation_mf',mf_u,'refine', 4, 
'deformation_scale',1);  % 'deformed_mesh', 'on')
-    colorbar;
-    caxis([0 10000]);
-    n = t(step);
-    title(['Plastification for t = ', num2str(n)]);
-    
-    % pause();
-
-end;
-
-
-
-
-
-
-
-

Modified: trunk/getfem/src/getfem_plasticity.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_plasticity.cc?rev=4768&r1=4767&r2=4768&view=diff
==============================================================================
--- trunk/getfem/src/getfem_plasticity.cc       (original)
+++ trunk/getfem/src/getfem_plasticity.cc       Thu Oct  2 16:46:03 2014
@@ -689,13 +689,13 @@
   //
   //=========================================================================
 
-  // static void ga_init_scalar(bgeot::multi_index &mi) { mi.resize(0); }
-  static void ga_init_vector(bgeot::multi_index &mi, size_type N)
-  { mi.resize(1); mi[0] = N; }
+  static void ga_init_scalar(bgeot::multi_index &mi) { mi.resize(0); }
+  //  static void ga_init_vector(bgeot::multi_index &mi, size_type N)
+  // { mi.resize(1); mi[0] = N; }
   static void ga_init_matrix(bgeot::multi_index &mi, size_type M, size_type N)
   { mi.resize(2); mi[0] = M; mi[1] = N; }
   static void ga_init_square_matrix(bgeot::multi_index &mi, size_type N)
-  { mi.resize(2); mi[0] = mi[1] = N; }
+    { mi.resize(2); mi[0] = mi[1] = N; }
 
 
   bool expm(const base_matrix &a_, base_matrix &aexp, scalar_type tol=1e-15) {
@@ -863,6 +863,98 @@
     }
   };
 
+  //=================================================================
+  // Von Mises projection
+  //=================================================================
+
+
+  struct Von_Mises_projection_operator : public ga_nonlinear_operator {
+    bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
+      if (args.size() != 2 || args[0]->sizes().size() <= 2
+          || args[1]->size() != 1) return false;
+      size_type N = (args[0]->sizes().size() == 2) ?  args[0]->sizes()[0] : 1;
+      if (args[0]->sizes().size() == 2 && args[0]->sizes()[1] != N) return 
false;
+      if (args[0]->sizes().size() != 2 && args[0]->size() != 1)  return false;
+      if (N > 1) ga_init_square_matrix(sizes, N); else ga_init_scalar(sizes);
+      return true;
+    }
+
+    // Value:
+    void value(const arg_list &args, base_tensor &result) const {
+      size_type N = (args[0]->sizes().size() == 2) ? args[0]->sizes()[0] : 1;
+      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);
+
+      for (size_type i = 0; i < N; ++i) tau_D(i,i) += tau_m;
+
+      gmm::copy(tau_D.as_vector(), result.as_vector());
+    }
+
+    // Derivative:
+    void derivative(const arg_list &args, size_type nder,
+                    base_tensor &result) const {
+      size_type N = (args[0]->sizes().size() == 2) ? args[0]->sizes()[0] : 1;
+      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, scalar_type(1)/norm_tau_D);
+
+      switch(nder) {
+      case 1: 
+        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)
+              result(i,i,j,j) += scalar_type(1)/scalar_type(N);
+        }
+        break;
+      case 2:
+        if (norm_tau_D < s)
+          gmm::clear(result.as_vector());
+        else
+          gmm::copy(tau_D.as_vector(), result.as_vector());
+        break;
+      }
+    }
+
+    // Second derivative : not implemented
+    void second_derivative(const arg_list &, size_type, size_type,
+                           base_tensor &) const {
+      GMM_ASSERT1(false, "Sorry, second derivative not implemented");
+    }
+  };
 
   static bool init_predef_operators(void) {
 
@@ -871,6 +963,8 @@
     
     PREDEF_OPERATORS.add_method("expm",
                                 new matrix_exponential_operator());
+    PREDEF_OPERATORS.add_method("Von_Mises_projection",
+                                new Von_Mises_projection_operator());
     return true;
    }