[Getfem-commits] (no subject)

[Yves Renard]

branch: devel-yves-ctrl-newton
commit f62d8dd528c0acd1f896de8a262e5337c7cef795
Author: Yves Renard <address@hidden>
Date:   Fri Oct 13 10:05:53 2017 +0200

    minor changes
---
 interface/tests/matlab/demo_nonlinear_elasticity.m | 23 +++++++++++-----------
 src/getfem/getfem_model_solvers.h                  | 17 +++++++++-------
 tests/nonlinear_elastostatic.param                 |  4 ++--
 3 files changed, 24 insertions(+), 20 deletions(-)

diff --git a/interface/tests/matlab/demo_nonlinear_elasticity.m 
b/interface/tests/matlab/demo_nonlinear_elasticity.m
index 709b631..d0d1376 100644
--- a/interface/tests/matlab/demo_nonlinear_elasticity.m
+++ b/interface/tests/matlab/demo_nonlinear_elasticity.m
@@ -26,12 +26,12 @@ colormap(r);
 dirichlet_version = 2; % 1 = simplification, 2 = penalisation
 drawing = true;
 test_tangent_matrix = false;
-incompressible = true;
+incompressible = false;
 
-% lawname = 'Ciarlet Geymonat';
-% params = [1;1;0.25];
-lawname = 'SaintVenant Kirchhoff';
-params = [1;1];
+lawname = 'Ciarlet Geymonat';
+params = [1;1;0.25];
+% lawname = 'SaintVenant Kirchhoff';
+% params = [1;1];
 if (incompressible)
     lawname = 'Incompressible Mooney Rivlin';
     params = [1;1];
@@ -41,9 +41,9 @@ N1=2; N2=4; h=20;
 m=gf_mesh('cartesian',(0:N1)/N1 - .5, (0:N2)/N2*h, ((0:N1)/N1 - .5)*3);
 mfu=gf_mesh_fem(m,3);     % mesh-fem supporting a 3D-vector field
 % the mesh_im stores the integration methods for each tetrahedron
-mim=gf_mesh_im(m,gf_Integ('IM_GAUSS_PARALLELEPIPED(3,4)'));
+mim=gf_mesh_im(m,gf_integ('IM_GAUSS_PARALLELEPIPED(3,4)'));
 % we choose a P2 fem for the main unknown
-gf_mesh_fem_set(mfu, 'fem',gf_Fem('FEM_QK(3,2)'));
+gf_mesh_fem_set(mfu, 'fem',gf_fem('FEM_QK(3,2)'));
 mfdu=gf_mesh_fem(m,1);
 % the material is homogeneous, hence we use a P0 fem for the data
 if (dirichlet_version == 1)
@@ -103,9 +103,10 @@ gf_model_set(md, 'add finite strain elasticity brick', 
mim, lawname, 'u',  'para
 %      
'((Id(meshdim)+Grad_u)*(Saint_Venant_Kirchhoff_sigma(Grad_u,params))):Grad_Test_u');
 
     
-if (incompressible || true)
-  mfp = gf_Mesh_Fem(m,1); 
-  gf_mesh_fem_set(mfp, 'classical discontinuous fem', 1);
+if (incompressible)
+  mfp = gf_mesh_fem(m,1); 
+  % gf_mesh_fem_set(mfp, 'classical discontinuous fem', 1);
+  gf_mesh_fem_set(mfp, 'classical fem', 1);
   gf_model_set(md, 'add fem variable', 'p', mfp);
   % gf_model_set(md, 'add nonlinear incompressibility brick',  mim, 'u', 'p');
   gf_model_set(md, 'add finite strain incompressibility brick',  mim, 'u', 
'p');
@@ -197,7 +198,7 @@ for step=1:nbstep,
     end
     
     gf_model_set(md, 'variable', 'DirichletData', R);
-    gf_model_get(md, 'solve', 'very noisy', 'max_iter', 100, 'max_res', 1e-5, 
'lsearch', 'simplest');
+    gf_model_get(md, 'solve', 'noisy', 'max_iter', 100, 'max_res', 1e-5); %  , 
'lsearch', 'simplest');
       
     if (test_tangent_matrix)
       gf_model_get(md, 'test tangent matrix', 1E-8, 10, 0.0001);
diff --git a/src/getfem/getfem_model_solvers.h 
b/src/getfem/getfem_model_solvers.h
index c39f6d7..278e7ba 100644
--- a/src/getfem/getfem_model_solvers.h
+++ b/src/getfem/getfem_model_solvers.h
@@ -480,14 +480,14 @@ namespace getfem {
        gmm::add(dr, pb.state_vector());
        pb.compute_residual();
        R res = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha));
-       R dec = R(1)/R(2);
+       R dec = R(1)/R(2), coeff2 = coeff * R(1.5);
        
        while (dec > R(1E-5) && res >= res0 * coeff) {
          gmm::add(gmm::scaled(dr, -dec), pb.state_vector());
          pb.compute_residual();
          R res2 = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha));
-         if (res2 < res*R(0.95)) {
-           dec /= R(2); res = res2;
+         if (res2 < res*R(0.95) || res2 >= res0 * coeff2) {
+           dec /= R(2); res = res2; coeff2 *= R(1.5);
          } else {
            gmm::add(gmm::scaled(dr, dec), pb.state_vector());
            break;
@@ -496,12 +496,15 @@ namespace getfem {
        dec *= R(2);
 
        nit++;
-       coeff = std::max(R(1), coeff*R(0.93));
-       if (res > minres && nit > 4) {
+       coeff = std::max(R(1.05), coeff*R(0.93));
+       bool near_end = (iter.get_iteration() > iter.get_maxiter()/2);
+       bool cut = (alpha < R(1)) && near_end;
+       if ((res > minres && nit > 4) || cut) {
          stagn++;
 
-         if (stagn > 10 && alpha > alpha0 + 5E-2) {
+         if ((stagn > 10 && alpha > alpha0 + R(5E-2)) || cut) {
            alpha = (alpha + alpha0) / R(2);
+           if (near_end) alpha = R(1);
            gmm::copy(Xi, pb.state_vector());
            pb.compute_residual();
            res = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha));
@@ -529,7 +532,7 @@ namespace getfem {
          alpha0 = a;
          gmm::copy(pb.state_vector(), Xi);
          nit = 0; stagn = 0; coeff = R(2);
-       } else if (alpha == 1) return;
+       } else return;
       }
     }
   }
diff --git a/tests/nonlinear_elastostatic.param 
b/tests/nonlinear_elastostatic.param
index 9914d71..cf3630d 100644
--- a/tests/nonlinear_elastostatic.param
+++ b/tests/nonlinear_elastostatic.param
@@ -26,7 +26,7 @@ LZ = 2.0;             % size in Z.
 P1 = 1.;               % First elastic coefficient.
 P2 = 1.;               % Second elastic coefficient.
 P3 = 0.5;              % Third elastic coefficient.
-LAW = 3;     % 0 : SaintVenant-Kirchhoff
+LAW = 2;     % 0 : SaintVenant-Kirchhoff
              % 1 : SaintVenant-Kirchhoff+incompressibility
              % 2 : Ciarlet-Geymonat
              % 3 : Mooney-Rivlin (+incompressibility)
@@ -81,7 +81,7 @@ INTEGRATION = 'IM_TETRAHEDRON(6)'
 %INTEGRATION = 'IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(7), 3)';
 
 RESIDUAL = 1E-6;       % residu for iterative solvers.
-MAXITER = 100000;
+MAXITER = 500;
 DIRICHLET_VERSION=0;
 NBSTEP = 40;
 %%%%%   saving parameters                                             %%%%%