#!/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')