[Getfem-commits] r5024 - in /trunk/getfem: doc/sphinx/source/tutorial/wheel.rst interface/tests/python/demo_wheel_contact.py

[Yves . Renard]

Author: renard
Date: Wed Jun  3 13:43:15 2015
New Revision: 5024

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

Modified:
    trunk/getfem/doc/sphinx/source/tutorial/wheel.rst
    trunk/getfem/interface/tests/python/demo_wheel_contact.py

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=5024&r1=5023&r2=5024&view=diff
==============================================================================
--- trunk/getfem/doc/sphinx/source/tutorial/wheel.rst   (original)
+++ trunk/getfem/doc/sphinx/source/tutorial/wheel.rst   Wed Jun  3 13:43:15 2015
@@ -199,7 +199,7 @@
 
 .. math::
 
-  \int_{\Gamma_D} \lambda_D.\delta_{u^1} + (u^1 - (0, 
-\alpha_D)).\delta_{\lambda_D} + (\lambda_D.(0,1) - F)\delta_{\alpha_D} d\Gamma 
= 0,
+  \cdots + \int_{\Gamma_D} -\lambda_D.\delta_{u^1} + ((0, \alpha_D) - 
u^1).\delta_{\lambda_D} + (\lambda_D.(0,1) + F)\delta_{\alpha_D} d\Gamma = 0,
 
 where :math:`\Gamma_D` is the rim boundary, :math:`F` is the applied density 
of force.
 
@@ -209,8 +209,15 @@
 
   md.add_filtered_fem_variable('lambda_D', mflambda, HOLE_BOUND)
   md.add_initialized_data('F', [applied_force/(8*2*np.pi)])
-  md.add_linear_generic_assembly_brick(mim1, 'lambda_D.Test_u1 + (u1 - 
alpha_D*[0;1]).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', HOLE_BOUND)
+
+For more robustness, a small penalization on :math:`alpha_D` can be added
+
+.. code-block:: python
+
+  md.add_linear_generic_assembly_brick(mim1, '1E-6*alpha_D*Test_alpha_D');
+
 
 Note that the fixed size variable `alpha_D` is linked to each points of the 
rim boundary. This means that the line of the tangent matrix corresponding to 
`alpha_D` may have a lot of nonzero components. This is why such a use of fixed 
size variable have to be done with care.
 

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=5024&r1=5023&r2=5024&view=diff
==============================================================================
--- trunk/getfem/interface/tests/python/demo_wheel_contact.py   (original)
+++ trunk/getfem/interface/tests/python/demo_wheel_contact.py   Wed Jun  3 
13:43:15 2015
@@ -39,14 +39,14 @@
 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)
+clambdastar = 2*clambda*cmu/(clambda+2*cmu) # Lame coefficient for Plane stress
 applied_force = 1E7              # Force at the hole boundary (N)
 
 
 #
 # Numerical parameters
 #
-h = 1.                   # Approximate mesh size
+h = 1                    # Approximate mesh size
 elements_degree = 2      # Degree of the finite element methods
 gamma0 = 1./E;           # Augmentation parameter for the augmented Lagrangian 
 
@@ -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, 6)
+mim1 = gf.MeshIm(mesh1, 4)
 mim1c = gf.MeshIm(mesh1, gf.Integ('IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(4),2)'))
 mim2 = gf.MeshIm(mesh2, 4)
 
@@ -124,7 +124,7 @@
 
 # Contact condition (augmented Lagrangian)
 md.add_initialized_data('gamma0', [gamma0])
-md.add_interpolate_transformation_from_expression('Proj1', mesh1, mesh2, 
'[X(1);0.]')
+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)
@@ -151,6 +151,7 @@
 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]
+# print 'Contact multiplier ', md.variable('lambda1')
 
 #
 # Solution export