1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
|
#!/usr/bin/env python
# -*- coding: UTF8 -*-
# Python GetFEM++ interface
#
# Copyright (C) 2011 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 2.1 of the License, 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 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.
#
############################################################################
""" Static equilibrium of an elastic solid in contact with a rigid foundation
This program is used to check that python-getfem is working. This is also
a good example of use of GetFEM++.
"""
import getfem as gf
import numpy as np
# Import the mesh : disc
# m = gf.Mesh('load', '../../../tests/meshes/disc_P2_h4.mesh')
#m = gf.Mesh('load', '../../../tests/meshes/disc_P2_h2.mesh')
# m = gf.Mesh('load', '../../../tests/meshes/disc_P2_h1.mesh')
# m = gf.Mesh('load', '../../../tests/meshes/disc_P2_h0_5.mesh')
# m = gf.Mesh('load', '../../../tests/meshes/disc_P2_h0_3.mesh')
# Import the mesh : sphere
# m = gf.Mesh('load', '../../../tests/meshes/sphere_with_quadratic_tetra_8_elts.mesh')
# m = gf.Mesh('load', '../../../tests/meshes/sphere_with_quadratic_tetra_80_elts.mesh')
m = gf.Mesh('load', '../../../tests/meshes/sphere_with_quadratic_tetra_400_elts.mesh')
# m = gf.Mesh('load', '../../../tests/meshes/sphere_with_quadratic_tetra_2000_elts.mesh')
# m = gf.Mesh('load', '../../../tests/meshes/sphere_with_quadratic_tetra_16000_elts.mesh')
d = m.dim() # Mesh dimension
# Parameters of the model
clambda = 1. # Lame coefficient
cmu = 1. # Lame coefficient
friction_coeff = 0.4 # coefficient of friction
vertical_force = 0.05 # Volumic load in the vertical direction
r = 10. # Augmentation parameter
condition_type = 0 # 0 = Explicitely kill horizontal rigid displacements
# 1 = Kill rigid displacements using a global penalization
# 2 = Add a Dirichlet condition on the top of the structure
penalty_parameter = 1E-6 # Penalization coefficient for the global penalization
if d == 2:
cpoints = [0, 0] # constrained points for 2d
cunitv = [1, 0] # corresponding constrained directions for 2d
else:
cpoints = [0, 0, 0, 0, 0, 0, 5, 0, 5] # constrained points for 3d
cunitv = [1, 0, 0, 0, 1, 0, 0, 1, 0] # corresponding constrained directions for 3d
niter = 100 # Maximum number of iterations for Newton's algorithm.
version = 13 # 1 : frictionless contact and the basic contact brick
# 2 : contact with 'static' Coulomb friction and basic contact brick
# 3 : frictionless contact and the contact with a
# rigid obstacle brick
# 4 : contact with 'static' Coulomb friction and the contact with a
# rigid obstacle brick
# 5 : frictionless contact and the integral brick
# Newton and Alart-Curnier augmented lagrangian,
# unsymmetric version
# 6 : frictionless contact and the integral brick
# Newton and Alart-Curnier augmented lagrangian, symmetric
# version.
# 7 : frictionless contact and the integral brick
# Newton and Alart-Curnier augmented lagrangian,
# unsymmetric version with an additional augmentation.
# 8 : frictionless contact and the integral brick
# New unsymmetric method.
# 9 : frictionless contact and the integral brick : Uzawa
# on the Lagrangian augmented by the penalization term.
# 10 : contact with 'static' Coulomb friction and the integral brick
# Newton and Alart-Curnier augmented lagrangian,
# unsymmetric version.
# 11 : contact with 'static' Coulomb friction and the integral brick
# Newton and Alart-Curnier augmented lagrangian,
# nearly symmetric version.
# 12 : contact with 'static' Coulomb friction and the integral brick
# Newton and Alart-Curnier augmented lagrangian,
# unsymmetric version with an additional augmentation.
# 13 : contact with 'static' Coulomb friction and the integral brick
# New unsymmetric method.
# 14 : contact with 'static' Coulomb friction and the integral brick : Uzawa
# on the Lagrangian augmented by the penalization term.
# 15 : penalized contact with 'static' Coulomb friction (r is the penalization
# coefficient).
# Signed distance representing the obstacle
if d == 2:
obstacle = 'y'
else:
obstacle = 'z'
# Selection of the contact and Dirichlet boundaries
GAMMAC = 1
GAMMAD = 2
border = m.outer_faces()
normals = m.normal_of_faces(border)
contact_boundary = border[:,np.nonzero(normals[d-1] < -0.01)[0]]
m.set_region(GAMMAC, contact_boundary)
contact_boundary = border[:,np.nonzero(normals[d-1] > 0.01)[0]]
m.set_region(GAMMAD, contact_boundary)
# Finite element methods
u_degree = 2
lambda_degree = 2
mfu = gf.MeshFem(m, d)
mfu.set_classical_fem(u_degree)
mfd = gf.MeshFem(m, 1)
mfd.set_classical_fem(u_degree)
mflambda = gf.MeshFem(m, 1) # used only by version 5 to 13
mflambda.set_classical_fem(lambda_degree)
mfvm = gf.MeshFem(m, 1)
mfvm.set_classical_discontinuous_fem(u_degree-1)
# Integration method
mim = gf.MeshIm(m, 4)
if d == 2:
mim_friction = gf.MeshIm(m,
gf.Integ('IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(4),4)'))
else:
mim_friction = gf.MeshIm(m,
gf.Integ('IM_STRUCTURED_COMPOSITE(IM_TETRAHEDRON(5),4)'))
# Volumic density of force
nbdofd = mfd.nbdof()
nbdofu = mfu.nbdof()
F = np.zeros(nbdofd*d)
F[d-1:nbdofd*d:d] = -vertical_force;
# Elasticity model
md = gf.Model('real')
md.add_fem_variable('u', mfu)
md.add_initialized_data('cmu', [cmu])
md.add_initialized_data('clambda', [clambda])
md.add_isotropic_linearized_elasticity_brick(mim, 'u', 'clambda', 'cmu')
md.add_initialized_fem_data('volumicload', mfd, F)
md.add_source_term_brick(mim, 'u', 'volumicload')
if condition_type == 2:
Ddata = np.zeros(d)
Ddata[d-1] = -5
md.add_initialized_data('Ddata', Ddata)
md.add_Dirichlet_condition_with_multipliers(mim, 'u', u_degree, GAMMAD, 'Ddata')
elif condition_type == 0:
md.add_initialized_data('cpoints', cpoints)
md.add_initialized_data('cunitv', cunitv)
md.add_pointwise_constraints_with_multipliers('u', 'cpoints', 'cunitv')
elif condition_type == 1:
# Small penalty term to avoid rigid motion (should be replaced by an
# explicit treatment of the rigid motion with a constraint matrix)
md.add_initialized_data('penalty_param', [penalty_parameter])
md.add_mass_brick(mim, 'u', 'penalty_param')
# The contact condition
cdof = mfu.dof_on_region(GAMMAC)
nbc = cdof.shape[0] / d
solved = False
if version == 1 or version == 2: # defining the matrices BN and BT by hand
contact_dof = cdof[d-1:nbc*d:d]
contact_nodes = mfu.basic_dof_nodes(contact_dof)
BN = gf.Spmat('empty', nbc, nbdofu)
ngap = np.zeros(nbc)
for i in range(nbc):
BN[i, contact_dof[i]] = -1.
ngap[i] = contact_nodes[d-1, i]
if version == 2:
BT = gf.Spmat('empty', nbc*(d-1), nbdofu)
for i in range(nbc):
for j in range(d-1):
BT[j+i*(d-1), contact_dof[i]-d+j+1] = 1.0
md.add_variable('lambda_n', nbc)
md.add_initialized_data('r', [r])
if version == 2:
md.add_variable('lambda_t', nbc * (d-1))
md.add_initialized_data('friction_coeff', [friction_coeff])
md.add_initialized_data('ngap', ngap)
md.add_initialized_data('alpha', np.ones(nbc))
if version == 1:
md.add_basic_contact_brick('u', 'lambda_n', 'r', BN, 'ngap', 'alpha', 1)
else:
md.add_basic_contact_brick('u', 'lambda_n', 'lambda_t', 'r', BN, BT, 'friction_coeff', 'ngap', 'alpha', 1);
elif version == 3 or version == 4: # BN and BT defined by the contact brick
md.add_variable('lambda_n', nbc)
md.add_initialized_data('r', [r])
if version == 3:
md.add_nodal_contact_with_rigid_obstacle_brick(mim, 'u', 'lambda_n', 'r', GAMMAC, obstacle, 1);
else:
md.add_variable('lambda_t', nbc*(d-1))
md.add_initialized_data('friction_coeff', [friction_coeff])
md.add_nodal_contact_with_rigid_obstacle_brick(mim, 'u', 'lambda_n', 'lambda_t', 'r',
'friction_coeff', GAMMAC, obstacle, 1)
elif version >= 5 and version <= 8: # The integral version, Newton
ldof = mflambda.dof_on_region(GAMMAC)
mflambda_partial = gf.MeshFem('partial', mflambda, ldof)
md.add_fem_variable('lambda_n', mflambda_partial)
md.add_initialized_data('r', [r])
OBS = mfd.eval(obstacle)
md.add_initialized_fem_data('obstacle', mfd, OBS)
md.add_integral_contact_with_rigid_obstacle_brick(mim_friction, 'u', 'lambda_n',
'obstacle', 'r', GAMMAC, version-4);
elif version == 9: # The integral version, Uzawa on the augmented Lagrangian
ldof = mflambda.dof_on_region(GAMMAC)
mflambda_partial = gf.MeshFem('partial', mflambda, ldof)
nbc = mflambda_partial.nbdof()
M = gf.asm_mass_matrix(mim, mflambda_partial, mflambda_partial, GAMMAC)
lambda_n = np.zeros(nbc)
md.add_initialized_fem_data('lambda_n', mflambda_partial, lambda_n)
md.add_initialized_data('r', [r])
OBS = mfd.eval(obstacle) # np.array([mfd.eval(obstacle)])
md.add_initialized_fem_data('obstacle', mfd, OBS)
md.add_penalized_contact_with_rigid_obstacle_brick \
(mim_friction, 'u', 'obstacle', 'r', GAMMAC, 2, 'lambda_n')
for ii in range(100):
print 'iteration %d' % (ii+1)
md.solve('max_res', 1E-9, 'max_iter', niter)
U = md.get('variable', 'u')
lambda_n_old = lambda_n
sol = gf.linsolve_superlu(M, gf.asm_integral_contact_Uzawa_projection(GAMMAC, mim_friction, mfu, U, mflambda_partial, lambda_n, mfd, OBS, r))
lambda_n = sol[0].transpose()
md.set_variable('lambda_n', lambda_n)
difff = max(abs(lambda_n-lambda_n_old))[0]/max(abs(lambda_n))[0]
print 'diff : %g' % difff
if difff < penalty_parameter:
break
solved = True
elif version >= 10 and version <= 13: # The integral version with friction, Newton
mflambda.set_qdim(d);
ldof = mflambda.dof_on_region(GAMMAC)
mflambda_partial = gf.MeshFem('partial', mflambda, ldof)
md.add_fem_variable('lambda', mflambda_partial)
md.add_initialized_data('r', [r])
md.add_initialized_data('friction_coeff', [friction_coeff])
OBS = mfd.eval(obstacle)
md.add_initialized_fem_data('obstacle', mfd, OBS)
md.add_integral_contact_with_rigid_obstacle_brick \
(mim_friction, 'u', 'lambda', 'obstacle', 'r', 'friction_coeff', GAMMAC, version-9)
elif version == 14: # The integral version, Uzawa on the augmented Lagrangian with friction
mflambda.set_qdim(d)
ldof = mflambda.dof_on_region(GAMMAC)
mflambda_partial = gf.MeshFem('partial', mflambda, ldof)
nbc = mflambda_partial.nbdof()
md.add_initialized_data('friction_coeff', [friction_coeff])
M = gf.asm_mass_matrix(mim, mflambda_partial, mflambda_partial, GAMMAC)
lambda_nt = np.zeros(nbc)
md.add_initialized_fem_data('lambda', mflambda_partial, lambda_nt)
md.add_initialized_data('r', [r])
OBS = mfd.eval(obstacle)
md.add_initialized_fem_data('obstacle', mfd, OBS)
md.add_penalized_contact_with_rigid_obstacle_brick \
(mim_friction, 'u', 'obstacle', 'r', 'friction_coeff', GAMMAC, 2, 'lambda')
for ii in range(100):
print 'iteration %d' % (ii+1)
md.solve('max_res', 1E-9, 'max_iter', niter)
U = md.get('variable', 'u')
lambda_nt_old = lambda_nt
sol = gf.linsolve_superlu(M,
gf.asm_integral_contact_Uzawa_projection(
GAMMAC, mim_friction, mfu, U, mflambda_partial, lambda_nt, mfd, OBS, r, friction_coeff))
lambda_nt = sol[0].transpose()
md.set_variable('lambda', lambda_nt)
difff = max(abs(lambda_nt-lambda_nt_old))[0]/max(abs(lambda_nt))[0]
print 'diff : %g' % difff
if difff < penalty_parameter:
break
solved = True
elif version == 15:
md.add_initialized_data('r', [r])
md.add_initialized_data('friction_coeff', [friction_coeff])
OBS = mfd.eval(obstacle)
md.add_initialized_fem_data('obstacle', mfd, OBS);
md.add_penalized_contact_with_rigid_obstacle_brick \
(mim_friction, 'u', 'obstacle', 'r', 'friction_coeff', GAMMAC)
else:
print 'Inexistent version'
# Solve the problem
if not solved:
md.solve('max_res', 1E-9, 'very noisy', 'max_iter', niter, 'lsearch', 'default') #, 'with pseudo potential')
U = md.get('variable', 'u')
# LAMBDA = md.get('variable', 'lambda_n')
VM = md.compute_isotropic_linearized_Von_Mises_or_Tresca('u', 'clambda', 'cmu', mfvm)
mfd.export_to_vtk('static_contact.vtk', 'ascii', mfvm, VM, 'Von Mises Stress', mfu, U, 'Displacement')
|
