Hello Yves,
Than you very much for the reply. I thought and tried that rg1 and rg2
should be boundary numbers as you explained (and in the documents). Just
in case, I am attaching my code with a comment to highlight the line
issues an error. It is line number 200. Let me know if there is anything
unclear about this revised code.
Best regards,
Hojae
On 12/23/2016 2:43 PM, Yves Renard wrote:
Dear Hojae,
The best is that you send the lines of your code that are not working.
rg1 and rg2 are boundary numbers so that they should be some integers.
Yves.
----- Original Message -----
From: "Hojae Yi" <[email protected]>
To: "getfem-users" <[email protected]>
Sent: Thursday, December 22, 2016 8:25:34 PM
Subject: [Getfem-users] nodal contact problem
Hello,
First of all, thank you so much for a great FEA package!
I am currently using GetFEM 5.1 using Python. I have been toying with
provided 'demo_wheel_contact.py' and trying to apply
add_nodal_contact_between_nonmatching_meshes_brick using two contact
surfaces. However, when I run, it complains that rg1 (and rg2) should be
a string. Python Interface reference explains that rg1 and rg2 as
integers while User Document explains them as vectors.
Can someone kindly point me where I can find further information on how
these rg1 and rg2 should be defined? Let me know if any further details
are needed.
Thank you very much!
Hojae
_______________________________________________
Getfem-users mailing list
[email protected]
https://mail.gna.org/listinfo/getfem-users
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Python GetFEM++ interface
#
# 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.
#
############################################################################
import getfem as gf
import numpy as np
gf.util('trace level', 3) # No trace for mesh generation nor for assembly
export_mesh = False;
Dirichlet_version = False; #True # Use a dirichlet condition instead of a
global load
#
# Physical parameters
#
#E = 2.1E5 # Young's Modulus (N/m^2)
approximately 30 psi
E = 2.1E10 # Young's Modulus (N/m^2)
approximately 30 psi + rubber (0.01 ~ 0.1 GPa)
nu = 0.3 # Poisson ratio
clambda = E*nu/((1+nu)*(1-2*nu)) # First Lame coefficient
(N/m^2)
cmu = E/(2*(1+nu)) # Second Lame coefficient
(N/m^2)
clambdastar = 2*clambda*cmu/(clambda+2*cmu) # Lame coefficient for
Plane stress
applied_force = 9.8E4 # Force at the hole
boundary (N) approximately 0.1 ton
E_f = 17.E9 # Concrete
nu_f = 0.15 # Poisson's ratio for the
foundation
clambda_f = E_f*nu_f/((1+nu_f)*(1-2*nu_f)) # First Lame coefficient
(N/m^2)
cmu_f = E_f/(2*(1+nu_f)) # Second Lame coefficient
(N/m^2)
clambdastar_f = 2*clambda_f*cmu_f/(clambda_f+2*cmu_f) # Lame coefficient for
Plane stress
friction_coeff = 0.4 ## coefficient of friction
#
# Numerical parameters
#
r = 5.0 # Augmentation Parameter
h = 0.05 # Approximate mesh size
elements_degree = 2 # Degree of the finite element methods
gamma0 = 1./E; # Augmentation parameter for the augmented Lagrangian
version = 1 # 1 : frictionless contact and the basic contact brick
# 2 : contact with 'static' Coulomb friction and basic contact
brick
#
# Mesh generation. Meshes can also been imported from several formats.
#
## set up base geometries
mo1 = gf.MesherObject('ball', [0., 1.5], 1.5) ## wheel's outer surface; 'ball',
vec center, scalar radius
mo2 = gf.MesherObject('ball', [0., 1.5], 0.8) ## wheel's inner surface;
mo3 = gf.MesherObject('set minus', mo1, mo2) ## whell with hole at the center
dimension = 2 ## Dimension of the problem
print ('Meshes are being generated...')
## wheel
mesh1 = gf.Mesh('generate', mo3, h, 2)
## bottom foundation
mesh2_length = 3
mesh2_height = 1
mesh2 =
gf.Mesh('import','structured','GT="GT_PK(2,1)";SIZES=[%d,%d];NOISED=0;NSUBDIV=[%d,%d];'
% (mesh2_length,mesh2_height,int(mesh2_length/h)+1, int(mesh2_height/h)+1));
## move mesh2's origion to [-1.5, -1.0] to align with the wheel center with the
foundation
mesh2.translate([-1.5,-1.0])
if (export_mesh):
mesh1.export_to_vtk('mesh1.vtk')
mesh2.export_to_vtk('mesh2.vtk')
print ('\nYou can view the meshes for instance with')
print ('mayavi2 -d mesh1.vtk -f ExtractEdges -m Surface -d mesh2.vtk -f
ExtractEdges -m Surface \n')
#
# Boundary selection
#
fb1 = mesh1.outer_faces_in_box([-0.81, 0.69], [0.81, 2.31]) # Boundary of the
inner hole
fb2 = mesh1.outer_faces_with_direction([0., -1.], np.pi/4) # Contact boundary
of the wheel ## increased contact boundary
fb3 = mesh2.outer_faces_with_direction([0., -1.], 0.01) # Bottom boundary of
the foundation; faces of bottom horizontal line
fb4 = mesh2.outer_faces_with_direction([0., 1.], np.pi/4) # Contact boundary of
the foundation
## define boundary region numbers
HOLE_BOUND=1; CONTACT_BOUND_WHEEL=2; BOTTOM_BOUND=3; CONTACT_BOUND_FND = 4
mesh1.set_region(HOLE_BOUND, fb1)
mesh1.set_region(CONTACT_BOUND_WHEEL, fb2)
mesh1.region_subtract(CONTACT_BOUND_WHEEL, HOLE_BOUND)
mesh2.set_region(BOTTOM_BOUND, fb3)
mesh2.set_region(CONTACT_BOUND_FND, fb4)
#
# Definition of finite elements methods and integration method
#
mfu1 = gf.MeshFem(mesh1, 2) ## on Mesh1 (wheel), qdim=2 for a vector field of
size 2
mfu1.set_classical_fem(elements_degree) ## fem with order of 'elements_degree'=2
mfu2 = gf.MeshFem(mesh2, 2)
mfu2.set_classical_fem(elements_degree)
mflambda = gf.MeshFem(mesh1, 2) ##
mflambda.set_classical_fem(elements_degree-1) ## fem of order 1 =
elements_degree - 1
mflambda_C = gf.MeshFem(mesh1, 1) ## on Mesh 1 (wheel), qdim=2 for a scalar
field
mflambda_C.set_classical_fem(elements_degree-1)
mfvm1 = gf.MeshFem(mesh1, 1) ## on mesh1 (wheel), qdim = 1 for a scalar field
mfvm1.set_classical_discontinuous_fem(elements_degree)
mfvm2 = gf.MeshFem(mesh2, 1) ## on mesh2 (foundation), qdim = 2 for a scalar
field
mfvm2.set_classical_discontinuous_fem(elements_degree)
## mim1, mim2 are two integration methods on the wheel and the foundation
mim1 = gf.MeshIm(mesh1, 4)
mim1c = gf.MeshIm(mesh1, gf.Integ('IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(4),2)'))
mim2 = gf.MeshIm(mesh2, 4)
## various mesh descriptors
cdof_wheel = mfu1.dof_on_region(CONTACT_BOUND_WHEEL)
cdof_fnd = mfu2.dof_on_region(CONTACT_BOUND_FND)
nbc_wheel = cdof_wheel.shape[0] / dimension
nbc_fnd = cdof_fnd.shape[0] / dimension
nbdofu1 = mfu1.nbdof()
nbdofu2 = mfu2.nbdof()
#
# Model definition
#
md=gf.Model('real');
md.add_fem_variable('u1', mfu1) # Displacement of the structure 1; wheel
md.add_fem_variable('u2', mfu2) # Displacement of the structure 2;
foundation
md.add_initialized_data('cmu', [cmu]) ## First Lame Coefficient for wheel
md.add_initialized_data('clambda', [clambda]) ## Second Lame Coefficient for
wheel
md.add_initialized_data('clambdastar', [clambdastar]) ## Second Lame
Coefficient for wheel with plane strain
md.add_initialized_data('cmu_f', [cmu_f]) ## First Lame Coefficient for
foundation
md.add_initialized_data('clambda_f', [clambda_f]) ## Second Lame Coefficient
for foundation
md.add_initialized_data('clambdastar_f', [clambdastar_f]) ## Second Lame
Coefficient for foundation for plane strain
md.add_isotropic_linearized_elasticity_brick(mim1, 'u1', 'clambda', 'cmu') ##
Linear Elasticity using Lame constants on MeshIm (mim1) and a fem variable
('u1')
md.add_isotropic_linearized_elasticity_brick(mim2, 'u2', 'clambda_f', 'cmu_f')
## Linear Elasticity using Lame constants on MeshIm (mim1) and a fem variable
('u1') for plane strain case
bottom_bv = mfu2.eval('[0,0]') ### Bottom boundary prescribed value (fixed)
md.add_initialized_fem_data('bbv', mfu2, bottom_bv)
#md.add_Dirichlet_condition_with_multipliers(mim2, 'u2', elements_degree-1,
BOTTOM_BOUND) ## Dirichlet boundary condition with multipliers on the fem
variable 'u2' on 'BOTTOM_BOUND' boundary. mult_description = elements_degree-1
= 1
## Above Dirichlet condition implies zero displacement at the bottom of the
foundation
md.add_Dirichlet_condition_with_multipliers(mim2, 'u2', elements_degree-1,
BOTTOM_BOUND, 'bbv') ## or use prescribed value with the given Lagrange
multiplier
### Conatac Conditions:
if version == 91:
# Contact condition (augmented Lagrangian)
md.add_initialized_data('gamma0', [gamma0]) ## Augmentation parameter for
the augmented Lagrangian defined as gamma0 = 1./E
md.add_interpolate_transformation_from_expression('Proj1', mesh1, mesh2,
'[X(1);0]') ## Project mesh1's coordinates to mesh2's coordinate
md.add_filtered_fem_variable('lambda1', mflambda_C, CONTACT_BOUND_WHEEL) #
mflambda_C is a contact multiplier MeshFem
## An integral contact condition using augmented Lagrangian formulations
using the high-level generic assembly bricks
md.add_nonlinear_generic_assembly_brick(mim1c, 'lambda1*(Test_u1.[0;1])'
'-
lambda1*(Interpolate(Test_u2,Proj1).[0;1])', CONTACT_BOUND_WHEEL)
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_WHEEL);
elif version == 1 or version == 2: # defining the matrices BN (B in normal
direction) and BT (T in traction direction) by hand ## Boundary normal/traction
## for the wheel
contact_wheel_dof = cdof_wheel[dimension-1:nbc_wheel*dimension:dimension]
## DOF in y directions
contact_wheel_nodes = mfu1.basic_dof_nodes(contact_wheel_dof) ## get
location of basic degrees of freedom.
BN_wheel = gf.Spmat('empty', nbc_wheel, nbdofu1)
ngap_wheel = np.zeros(nbc_wheel)
for i in range(nbc_wheel):
BN_wheel[i, contact_wheel_dof[i]] = -1. ## set contact dof's as a
negative one; contact??
ngap_wheel[i] = contact_wheel_nodes[dimension-1, i] ## gap defined at the
master contact ??
## for the foundation
contact_fnd_dof = cdof_fnd[dimension-1:nbc_fnd*dimension:dimension] ## DOF
in y directions
contact_fnd_nodes = mfu2.basic_dof_nodes(contact_fnd_dof) ## get location
of basic degrees of freedom.
BN_fnd = gf.Spmat('empty', nbc_fnd, nbdofu2)
#ngap_fnd = np.zeros(nbc_fnd)
for i in range(nbc_fnd):
BN_fnd[i, contact_fnd_dof[i]] = -1. ## set contact dof's as a negative
one; contact??
## ngap_fnd[i] = contact_fnd_nodes[dimension-1, i]
## model variables and data to account for the contact
md.add_variable('lambda_n', nbc_wheel) ## Number of contact boundary dof
for normal contact`
md.add_initialized_data('r', [r]) ## Augmentation parameter
md.add_initialized_data('f_coeff', 0.0)
md.add_initialized_data('ngap_wheel', ngap_wheel) ## number of Gap nodes
(dof)
md.add_initialized_data('alpha', np.ones(nbc_wheel)) ## an optional
homogenization parameter for the augmentation parameter
if version == 1: ## THIS SHOULD BE REPLACED WITH
add_*_contact_between_nonmatching_meshes_brick() ## 2016-12-20
###
### BELOW LINE COMPLAINS THAT ARGUMENT 01 MUST BE A STRING
###
md.add_nodal_contact_between_nonmatching_meshes_brick(mim1, mim2, 'u1',
'u2', 'lambda_n', 'lambda_t', 'r', mesh1.pid_in_regions(CONTACT_BOUND_WHEEL),
mesh2.pid_in_regions(CONTACT_BOUND_FND), True, False, 1)
#md.add_multiplier('lambda_n', mfu1, 'u1', mim1, CONTACT_BOUND_WHEEL)
#md.add_integral_contact_between_nonmatching_meshes_brick(mim1, 'u1',
'u2', 'lambda_n', 'r', CONTACT_BOUND_WHEEL, CONTACT_BOUND_FND) ## Frictionless
contact ## NOT WORKING PROPERLY
#md.add_integral_contact_between_nonmatching_meshes_brick(mim1, 'u1',
'u2', 'lambda_n', 'r', 'f_coeff', CONTACT_BOUND_WHEEL, CONTACT_BOUND_FND) ##
Frictional contact
elif version == 2:
BT = gf.Spmat('empty', nbc_wheel*(dimension-1), nbdofu1)
for i in range(nbc_wheel):
for j in range(dimension-1):
BT[j+i*(dimension-1), contact_wheel_dof[i]-dimension+j+1] = 1.0
md.add_variable('lambda_t', nbc_wheel * (dimension-1)) ## number of
contact boundary for traction contact
md.add_initialized_data('friction_coeff', [friction_coeff])
md.add_basic_contact_brick('u1', 'lambda_n', 'lambda_t', 'r', BN, BT,
'friction_coeff', 'ngap', 'alpha', 1);
# Prescribed force in the hole
if (Dirichlet_version): ## prescribed vertical displacement
md.add_initialized_data('DData', [0., -0.1]) ## vertical displacement of 0.1
m deformation in the negative y direction
md.add_Dirichlet_condition_with_multipliers(mim1, 'u1', elements_degree-1,
HOLE_BOUND, 'DData'); ## vertical displacement BC on wheel's inner rim
else: ## prescribed force on the rim considering the rigidity of rim
## mflambda is meshfem to approximate a multiplier to take into account the
rigidity of the rim
md.add_filtered_fem_variable('lambda_D', mflambda, HOLE_BOUND) ## add a
variable to the model linked to a MeshFem, mflambda: The variable is filtered
in the sense that only the dof on the region are considered
md.add_initialized_data('F', [applied_force/(0.8*2*np.pi)]) ## applied_force
= 9.8E4 = 1 ton over the diameter of the rim
## add a variable to the model of constant sizes. 'sizes'=1 is either an
integer (for a scalar or vector variable) or a vector of dimensions for a
tensor variable.
## name = 'alpha_D' is th evariable name
md.add_variable('alpha_D', 1) ## Unknown vertical displacement of wheel rim
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-9*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-5, 'max_iter', 100, 'very_noisy' , 'lsearch',
'simplest', 'alpha min', 0.8)
#HY#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 (md.variable_list())
#print 'Contact multiplier ', md.variable('lambda1')
#
# Solution export
#
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)
if (Dirichlet_version):
mfvm1.export_to_vtk('displacement_with_von_mises1_DisplacementBC.vtk', mfvm1,
VM1, 'Von Mises Stresses', mfu1, U1, 'Displacements')
mfvm2.export_to_vtk('displacement_with_von_mises2_DisplacementBC.vtk', mfvm2,
VM2, 'Von Mises Stresses', mfu2, U2, 'Displacements')
else:
mfvm1.export_to_vtk('displacement_with_von_mises1_ForceBC.vtk', mfvm1, VM1,
'Von Mises Stresses', mfu1, U1, 'Displacements')
mfvm2.export_to_vtk('displacement_with_von_mises2_ForceBC.vtk', mfvm2, VM2,
'Von Mises Stresses', mfu2, U2, 'Displacements')
print ('You can view solutions with for instance:\nmayavi2 -d
displacement_with_von_mises1.vtk -f WarpVector -m Surface -d
displacement_with_von_mises2.vtk -f WarpVector -m Surface')
_______________________________________________
Getfem-users mailing list
[email protected]
https://mail.gna.org/listinfo/getfem-users
