[Getfem-commits] r5023 - in /trunk/getfem: doc/sphinx/source/tutorial/ interface/tests/matlab/ interface/tests/python/ src/

[Yves . Renard]

Author: renard
Date: Tue Jun  2 20:06:17 2015
New Revision: 5023

URL: http://svn.gna.org/viewcvs/getfem?rev=5023&view=rev
Log:
stabilizating demo_wheel_contact with a small penalization. adding matlab demo

Added:
    trunk/getfem/interface/tests/matlab/demo_wheel_contact.m
Modified:
    trunk/getfem/doc/sphinx/source/tutorial/wheel.rst
    trunk/getfem/interface/tests/matlab/Makefile.am
    trunk/getfem/interface/tests/python/demo_wheel_contact.py
    trunk/getfem/src/getfem_generic_assembly.cc

Modified: trunk/getfem/doc/sphinx/source/tutorial/wheel.rst
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/doc/sphinx/source/tutorial/wheel.rst?rev=5023&r1=5022&r2=5023&view=diff
==============================================================================
--- trunk/getfem/doc/sphinx/source/tutorial/wheel.rst   (original)
+++ trunk/getfem/doc/sphinx/source/tutorial/wheel.rst   Tue Jun  2 20:06:17 2015
@@ -151,7 +151,7 @@
 
   X \longmapsto (X(1), 0)
 
-where :math`X` is the vector of coordinates of the point. We add this 
transformation to the model with the command
+where :math:`X` is the vector of coordinates of the point. We add this 
transformation to the model with the command
 
 .. code-block:: python
 
@@ -164,7 +164,7 @@
 .. math::
 
   & \cdots + \ds \int_{\Gamma_c} \lambda_N(X) 
(\delta_{u^1}(X)-\delta_{u^2}(\Pi(X)))\cdot n d\Gamma \\
-  & -  \ds \int_{\Gamma_c} \left(\lambda_N(X) + \left(\lambda_N(X) - 
\Frac{1}{h_T\gamma_0}((X + u^1(X))\cdot n - (\Pi(X) - u^2(\Pi(X)))\cdot 
n\right)_-\right)\delta_{\lambda_N}(X) d\Gamma = 0 ~~~~ \forall 
\delta_{\lambda_N}, \forall \delta_{u^1}, \forall \delta_{u^2},
+  & -  \ds \int_{\Gamma_c} \left(\lambda_N(X) + \left(\lambda_N(X) + 
\Frac{1}{h_T\gamma_0}((X + u^1(X))\cdot n - (\Pi(X) - u^2(\Pi(X)))\cdot 
n\right)_-\right)\delta_{\lambda_N}(X) d\Gamma = 0 ~~~~ \forall 
\delta_{\lambda_N}, \forall \delta_{u^1}, \forall \delta_{u^2},
 
 where :math:`\Gamma_c` is the slave contact boundary, :math:`\lambda_N` is the 
contact multiplier (contact pressure), :math:`h_T` is the radius of the 
element, :math:`\Pi` is the transformation, `n` is the outward normal vector to 
the master contact boundary (here :math:`n = (0,1)`), :math:`\gamma_0` is an 
augmentation parameter, :math:`(\cdot)_-:I\hspace{-0.2em}R\rightarrow 
I\hspace{-0.2em}R_+` is the negative part and :math:`\delta_{\lambda_N}, 
\delta_{u^1}, \delta_{u^2}` are the test  functions corresponding to 
:math:`\lambda_N, u^1, u^2`, respectively.
 

Modified: trunk/getfem/interface/tests/matlab/Makefile.am
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/Makefile.am?rev=5023&r1=5022&r2=5023&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/Makefile.am     (original)
+++ trunk/getfem/interface/tests/matlab/Makefile.am     Tue Jun  2 20:06:17 2015
@@ -53,6 +53,7 @@
        demo_contact_fictitious_domain_nitsche.m \
        demo_static_contact.m \
        demo_large_sliding_contact.m \
+       demo_wheel_contact.m \
        demo_wave2D.m \
        demo_wave2D_alt.m \
        demo_wave_equation.m \

Added: trunk/getfem/interface/tests/matlab/demo_wheel_contact.m
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/demo_wheel_contact.m?rev=5023&view=auto
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_wheel_contact.m    (added)
+++ trunk/getfem/interface/tests/matlab/demo_wheel_contact.m    Tue Jun  2 
20:06:17 2015
@@ -0,0 +1,165 @@
+% Copyright (C) 2015-2015 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.
+
+%
+% Contact of a deformable 'wheel' onto a plane deformable obstacle.
+%
+
+view_mesh = true;
+Dirichlet_version = false; % A Dirichlet condition instead of the applied load
+
+%
+% Physical parameters
+%
+E = 21E6;                         % Young Modulus (N/cm^2)
+nu = 0.3;                         % Poisson ratio
+clambda = E*nu/((1+nu)*(1-2*nu)); % First Lame coefficient (N/cm^2)
+cmu = E/(2*(1+nu));               % Second Lame coefficient (N/cm^2)
+clambdastar = 2*clambda*cmu/(clambda+2*cmu); % Lame coefficient for Plane 
stress (N/cm^2)
+applied_force = 1E7;              % Force at the hole boundary (N)
+
+
+%
+% Numerical parameters
+%
+h = 4;                   % Approximate mesh size
+elements_degree = 2;     % Degree of the finite element methods
+gamma0 = 1./E;           % Augmentation parameter for the augmented Lagrangian 
+
+%
+% Mesh generation. Meshes can also been imported from several formats.
+%
+mo1 = gfMesherObject('ball', [0., 15.], 15.);
+mo2 = gfMesherObject('ball', [0., 15.], 8.);
+mo3 = gfMesherObject('set minus', mo1, mo2);
+
+disp('Meshes generation');
+gf_util('trace level', 1)   % No trace for mesh generation nor for assembly
+mesh1 = gfMesh('generate', mo3, h, 2);
+mesh2 = gfMesh('import','structured', 
sprintf('GT="GT_PK(2,1)";SIZES=[30,10];NOISED=0;NSUBDIV=[%d,%d];', 
floor(30/h)+1, floor(10/h)+1));
+mesh2.translate([-15.;-10.]);
+
+
+if (view_mesh)
+  figure(2);
+  gf_plot_mesh(mesh1);
+  hold on;
+  gf_plot_mesh(mesh2);
+  hold off;
+  pause(1);
+end
+
+%
+% Boundary selection
+%
+fb1 = mesh1.outer_faces_in_box([-8.1, 6.9], [8.1, 23.1]);  % Boundary of the 
hole
+fb2 = mesh1.outer_faces_with_direction([0., -1.], pi/4.5);  % Contact boundary 
of the wheel
+fb3 = mesh2.outer_faces_with_direction([0., -1.], 0.01);   % Bottom boundary 
of the foundation
+
+HOLE_BOUND=1; CONTACT_BOUND=2; BOTTOM_BOUND=3;
+
+mesh1.set_region(HOLE_BOUND, fb1);
+mesh1.set_region(CONTACT_BOUND, fb2);
+mesh1.region_subtract(CONTACT_BOUND, HOLE_BOUND);
+mesh2.set_region(BOTTOM_BOUND, fb3);
+
+
+%
+% Definition of finite elements methods and integration method
+%
+
+mfu1 = gfMeshFem(mesh1, 2);
+mfu1.set_classical_fem(elements_degree);
+mflambda = gfMeshFem(mesh1, 2);
+mflambda.set_classical_fem(elements_degree-1);
+mflambda_C = gfMeshFem(mesh1, 1);
+mflambda_C.set_classical_fem(elements_degree-1);
+mfu2 = gfMeshFem(mesh2, 2);
+mfu2.set_classical_fem(elements_degree);
+mfvm1 = gfMeshFem(mesh1, 1);
+mfvm1.set_classical_discontinuous_fem(elements_degree);
+mfvm2 = gfMeshFem(mesh2, 1);
+mfvm2.set_classical_discontinuous_fem(elements_degree);
+mim1  = gfMeshIm(mesh1, 6);
+mim1c = gfMeshIm(mesh1, gfInteg('IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(4),4)'));
+mim2  = gfMeshIm(mesh2, 4);
+
+
+%
+% Model definition
+%
+
+md=gfModel('real');
+md.add_fem_variable('u1', mfu1);       % Displacement of the structure 1
+md.add_fem_variable('u2', mfu2);       % Displacement of the structure 2
+
+md.add_initialized_data('cmu', [cmu]);
+md.add_initialized_data('clambdastar', [clambdastar]);
+md.add_isotropic_linearized_elasticity_brick(mim1, 'u1', 'clambdastar', 'cmu');
+md.add_isotropic_linearized_elasticity_brick(mim2, 'u2', 'clambdastar', 'cmu');
+md.add_Dirichlet_condition_with_multipliers(mim2, 'u2', elements_degree-1, 
BOTTOM_BOUND);
+
+
+
+% Contact condition (augmented Lagrangian)
+md.add_initialized_data('gamma0', [gamma0]);
+md.add_interpolate_transformation_from_expression('Proj1', mesh1, mesh2, 
'[X(1);0]');
+md.add_filtered_fem_variable('lambda1', mflambda_C, CONTACT_BOUND);
+md.add_nonlinear_generic_assembly_brick(mim1c, 
'lambda1*(Test_u1.[0;1])-lambda1*(Interpolate(Test_u2,Proj1).[0;1])', 
CONTACT_BOUND);
+md.add_nonlinear_generic_assembly_brick(mim1c, 
'-(gamma0*element_size)*(lambda1 + 
neg_part(lambda1+(1/(gamma0*element_size))*((u1-Interpolate(u2,Proj1)+X-Interpolate(X,Proj1)).[0;1])))*Test_lambda1',
 CONTACT_BOUND);
+
+% Prescribed force in the hole (or Dirichlet condition)
+if (Dirichlet_version)
+  md.add_initialized_data('DData', [0., -1.0]);
+  md.add_Dirichlet_condition_with_multipliers(mim1, 'u1', elements_degree-1, 
HOLE_BOUND, 'DData');
+else
+  md.add_filtered_fem_variable('lambda_D', mflambda, HOLE_BOUND);
+  md.add_initialized_data('F', [applied_force/(8.*2.*pi)]);
+  md.add_variable('alpha_D', 1);
+  md.add_linear_generic_assembly_brick(mim1, '-lambda_D.Test_u1 + 
(alpha_D*[0;1] - u1).Test_lambda_D + (lambda_D.[0;1] + F)*Test_alpha_D + 
1E-6*alpha_D*Test_alpha_D', HOLE_BOUND);
+  % The small penalization 1E-6*alpha_D*Test_alpha_D seems necessary to
+  % have a convergence in all cases. Why ?
+end
+
+%
+% Model solve
+%
+
+disp(sprintf('Solve problem with %d dofs', md.nbdof()));
+md.solve('max_res', 1E-9, 'max_iter', 40, 'noisy'); % , 'lsearch', 'simplest', 
 'alpha min', 0.8)
+if (~Dirichlet_version)
+  disp(sprintf('alpha_D = %g', md.variable('alpha_D')));
+end
+% md.test_tangent_matrix(1e-4, 10, 0.001);
+md.variable('lambda1')
+
+%
+% Solution draw
+%  
+U1 = md.variable('u1');
+U2 = md.variable('u2');
+VM1 = md.compute_isotropic_linearized_Von_Mises_or_Tresca('u1', 'clambdastar', 
'cmu', mfvm1);
+VM2 = md.compute_isotropic_linearized_Von_Mises_or_Tresca('u2', 'clambdastar', 
'cmu', mfvm2);
+
+figure(1);
+gf_plot(mfvm1, VM1, 'deformed_mesh', 'on', 'deformation', U1, 
'deformation_mf', mfu1, 'deformation_scale', 1, 'refine', 8);
+hold on;
+gf_plot(mfvm2, VM2, 'deformed_mesh', 'on', 'deformation', U2, 
'deformation_mf', mfu2, 'deformation_scale', 1, 'refine', 8);
+xlabel('x'); ylabel('y');
+title('Deformed configuration');
+hold off;
+

Modified: trunk/getfem/interface/tests/python/demo_wheel_contact.py
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/python/demo_wheel_contact.py?rev=5023&r1=5022&r2=5023&view=diff
==============================================================================
--- trunk/getfem/interface/tests/python/demo_wheel_contact.py   (original)
+++ trunk/getfem/interface/tests/python/demo_wheel_contact.py   Tue Jun  2 
20:06:17 2015
@@ -30,7 +30,7 @@
 gf.util('trace level', 1)    # No trace for mesh generation nor for assembly
 
 export_mesh = True;
-Dirichlet_condition = False; # Use a dirichlet condition instead of a global 
load
+Dirichlet_version = False; # Use a dirichlet condition instead of a global load
 
 #
 # Physical parameters
@@ -101,7 +101,7 @@
 mfvm1.set_classical_discontinuous_fem(elements_degree)
 mfvm2 = gf.MeshFem(mesh2, 1)
 mfvm2.set_classical_discontinuous_fem(elements_degree)
-mim1 = gf.MeshIm(mesh1, 8) # Order 8 seems to be a minimum. Why ?
+mim1 = gf.MeshIm(mesh1, 6)
 mim1c = gf.MeshIm(mesh1, gf.Integ('IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(4),2)'))
 mim2 = gf.MeshIm(mesh2, 4)
 
@@ -131,22 +131,25 @@
 md.add_nonlinear_generic_assembly_brick(mim1c, 
'-(gamma0*element_size)*(lambda1 + 
neg_part(lambda1+(1/(gamma0*element_size))*((u1-Interpolate(u2,Proj1)+X-Interpolate(X,Proj1)).[0;1])))*Test_lambda1',
 CONTACT_BOUND);
 
 # Prescribed force in the hole
-if (Dirichlet_condition):
+if (Dirichlet_version):
   md.add_initialized_data('DData', [0., -1.0])
   md.add_Dirichlet_condition_with_multipliers(mim1, 'u1', elements_degree-1, 
HOLE_BOUND, 'DData');
 else:
   md.add_filtered_fem_variable('lambda_D', mflambda, HOLE_BOUND)
   md.add_initialized_data('F', [applied_force/(8*2*np.pi)])
   md.add_variable('alpha_D', 1)
-  md.add_linear_generic_assembly_brick(mim1, 'lambda_D.Test_u1 + 
(alpha_D*[0;1] - u1).Test_lambda_D + (lambda_D.[0;1] - F)*Test_alpha_D', 
HOLE_BOUND)
+  md.add_linear_generic_assembly_brick(mim1, '-lambda_D.Test_u1 + 
(alpha_D*[0;1] - u1).Test_lambda_D + (lambda_D.[0;1] + F)*Test_alpha_D + 
1E-6*alpha_D*Test_alpha_D', HOLE_BOUND)
+  # The small penalization 1E-6*alpha_D*Test_alpha_D seems necessary to have
+  # a convergence in all cases. Why ?
+
 
 #
 # Model solve
 #
 
 print 'Solve problem with ', md.nbdof(), ' dofs'
-md.solve('max_res', 1E-9, 'max_iter', 20, 'noisy') # , 'lsearch', 'simplest',  
'alpha min', 0.8)
-if not(Dirichlet_condition):
+md.solve('max_res', 1E-9, 'max_iter', 40, 'noisy') # , 'lsearch', 'simplest',  
'alpha min', 0.8)
+if not(Dirichlet_version):
   print 'alpha_D = ', md.variable('alpha_D')[0]
 
 #

Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=5023&r1=5022&r2=5023&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Tue Jun  2 20:06:17 2015
@@ -10307,6 +10307,12 @@
         bool converged = true;
         bool is_in = gic.invert(P, P_ref, converged);
 
+        // cout << "cv = " << cv << " P = " << P << " P_ref = " << P_ref
+        // cout << " is_in = " << int(is_in) << endl;
+        // for (size_type iii = 0;
+        //      iii < target_mesh.points_of_convex(cv).size(); ++iii)
+        //  cout << target_mesh.points_of_convex(cv)[iii] << endl;
+        
         if (is_in && converged) {
           face_num = short_type(-1); // Should detect potential faces ?
           ret_type = 1;