[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] r5022 - in /trunk/getfem: interface/tests/python/demo_w
From: |
Yves . Renard |
Subject: |
[Getfem-commits] r5022 - in /trunk/getfem: interface/tests/python/demo_wheel_contact.py src/getfem_models.cc |
Date: |
Tue, 02 Jun 2015 15:40:59 -0000 |
Author: renard
Date: Tue Jun 2 17:40:58 2015
New Revision: 5022
URL: http://svn.gna.org/viewcvs/getfem?rev=5022&view=rev
Log:
Passing linear elasticity brick to the high level generic assembly
Modified:
trunk/getfem/interface/tests/python/demo_wheel_contact.py
trunk/getfem/src/getfem_models.cc
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=5022&r1=5021&r2=5022&view=diff
==============================================================================
--- trunk/getfem/interface/tests/python/demo_wheel_contact.py (original)
+++ trunk/getfem/interface/tests/python/demo_wheel_contact.py Tue Jun 2
17:40:58 2015
@@ -20,13 +20,17 @@
# Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
#
############################################################################
+#
+# Contact of a deformable 'wheel' onto a plane deformable obstacle.
+#
+############################################################################
import getfem as gf
import numpy as np
-
-# Contact of a deformable 'wheel' onto a plane deformable obstacle.
+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
#
# Physical parameters
@@ -42,7 +46,7 @@
#
# 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
@@ -54,7 +58,6 @@
mo3 = gf.MesherObject('set minus', mo1, mo2)
print ('Meshes generation')
-gf.util('trace level', 1) # No trace for mesh generation nor for assembly
mesh1 = gf.Mesh('generate', mo3, h, 2)
mesh2 =
gf.Mesh('import','structured','GT="GT_PK(2,1)";SIZES=[30,10];NOISED=0;NSUBDIV=[%d,%d];'
% (int(30/h)+1, int(10/h)+1));
mesh2.translate([-15.,-10.])
@@ -98,8 +101,9 @@
mfvm1.set_classical_discontinuous_fem(elements_degree)
mfvm2 = gf.MeshFem(mesh2, 1)
mfvm2.set_classical_discontinuous_fem(elements_degree)
-mim1 = gf.MeshIm(mesh1, pow(elements_degree,2))
-mim2 = gf.MeshIm(mesh2, pow(elements_degree,2))
+mim1 = gf.MeshIm(mesh1, 8) # Order 8 seems to be a minimum. Why ?
+mim1c = gf.MeshIm(mesh1, gf.Integ('IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(4),2)'))
+mim2 = gf.MeshIm(mesh2, 4)
#
@@ -120,32 +124,30 @@
# 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(mim1, 'lambda1*(Test_u1.[0;1])'
+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(mim1, '-(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);
+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
-
-# md.add_initialized_data('DData', [0., -1.0])
-# md.add_Dirichlet_condition_with_multipliers(mim1, 'u1', elements_degree-1,
HOLE_BOUND, 'DData');
-
-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)
+if (Dirichlet_condition):
+ 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)
#
# Model solve
#
print 'Solve problem with ', md.nbdof(), ' dofs'
-print 'alpha_D = ', md.variable('alpha_D')[0]
-md.test_tangent_matrix(1e-4, 10, 0.001);
md.solve('max_res', 1E-9, 'max_iter', 20, 'noisy') # , 'lsearch', 'simplest',
'alpha min', 0.8)
-print 'alpha_D = ', md.variable('alpha_D')[0]
-md.test_tangent_matrix(1e-4, 10, 0.001);
+if not(Dirichlet_condition):
+ print 'alpha_D = ', md.variable('alpha_D')[0]
#
# Solution export
Modified: trunk/getfem/src/getfem_models.cc
URL:
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_models.cc?rev=5022&r1=5021&r2=5022&view=diff
==============================================================================
--- trunk/getfem/src/getfem_models.cc (original)
+++ trunk/getfem/src/getfem_models.cc Tue Jun 2 17:40:58 2015
@@ -6113,7 +6113,7 @@
+ test_varname;
-#if 1 // Old brick
+#if 0 // Old brick
pbrick pbr = new iso_lin_elasticity_brick(dataname3.size() ? expr1:expr2);
model::termlist tl;
tl.push_back(model::term_description(varname, varname, true));
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] r5022 - in /trunk/getfem: interface/tests/python/demo_wheel_contact.py src/getfem_models.cc,
Yves . Renard <=