set contact element

Ning

Gmsh does not support contact elements so one need to apply some tricks to get this going. There is also not much support to create meshes with contact elements in esys-escript. They have been introduced for a very specific project which had it own means to create meshes with contact elements.
There is a function function esys.finley.JoinFaces which merges meshes and creates contact elements across contacting faces but the meshes on the touching faces need to match. This function needs not very reliable and may require more work to be used with gmsh files but I would be happy to help you to get the function going if you want to try it.

Lutz

Hello Lutz,

Thanks for your reply and kindness to help. One of my friends sent me a script to generate simple meshes with contact elements, which can be imported into escript. I have worked around with it and succeeded to reproduce the results for a sliding block test (Oden & Pires, 1984), but I was struggling with those 'negative-positive signs' in the PDE coefficients. I attached the scripts here. Hope it could be useful to others.

Ning

#### sliding_block.py ####
import numpy as np
from esys.escript import *
from esys.escript.linearPDEs import LinearPDE
from esys.finley import ReadMesh
from esys.weipa import saveVTK

# material parameter
E = 1000.; nu = .3 # Young's modulus and Poisson's ratio
mu = .5; kt = E * 10; kn = kt * 10 # frictional coefficient and contact stiffness
lam1=E*nu/((1+nu)*(1-2*nu)); lam2=.5*E/(1+nu);
# dimension
lx = 4.; ly = 2.2; base = .2

mydomain=ReadMesh('mesh.fly')
mypde=LinearPDE(mydomain)

k=kronecker(mydomain)
kxk=outer(k,k)
C=lam2*(swap_axes(kxk,0,3)+swap_axes(kxk,1,3)) + lam1*kxk

x=mydomain.getX()
fx = Function(mydomain).getX()
bx = FunctionOnBoundary(mydomain).getX()
cx = FunctionOnContactOne(mydomain).getX()

Dbc = whereZero(x[1])*[1,1] # fixed bottom
top = whereZero(bx[1]-ly) * wherePositive(bx[0]-.2) * whereNegative(bx[0]-3.8) # where surcharge applies
Nbc = whereZero(bx[0]-lx)*[60.,0]+top*[0.,-200.]

fric_surf = whereZero(cx[1]-base) * wherePositive(cx[0]-.2) * whereNegative(cx[0]-3.8)

normal = FunctionOnContactOne(mydomain).getNormal()
tangential = normal[0]*[0,1] - normal[1]*[1,0]

d_contact = Tensor([[kt,0],[0,kn]],FunctionOnContactOne(mydomain))

S = Tensor4(C,Function(mydomain))
S = S + whereNegative(fx[1]-base) * S * 99 # make the base 100 times stiffer

mypde.setValue(A=S, y=Nbc, q=Dbc, d_contact=d_contact)
u = mypde.getSolution()

err = 1.0
while (err > 0.0001):
   g = grad(u)
   D = symmetric(g)
   stress = tensor_mult(S,D)
   jmp = jump(u)
   y_contact = matrix_mult(d_contact,jmp)
   fn = -inner(y_contact, normal)
   fn_mu = fn * mu * fric_surf
   ft = inner(y_contact, tangential)
   ft = maximum(minimum(abs(ft),fn_mu),0) * sign(ft)
   y_contact = -ft*tangential+fn*normal # struggling with the signs here
   mypde.setValue(X=-stress,y_contact=y_contact)
   du = mypde.getSolution()
   u += du
   err = L2(du)/L2(u) # any better criterion
   print err

mydomain.setX(x+u)
saveVTK('mesh.vtu',mydomain)

import matplotlib.pyplot as plt

fric = np.array(ft.toListOfTuples())
fnorm = np.array(fn.toListOfTuples())
x_coord = np.array(FunctionOnContactOne(mydomain).getX()[0].toListOfTuples())
xmsk = (x_coord > 0.2) * (x_coord < 3.8)

plt.plot(x_coord[xmsk],fric[xmsk],'rP',ms=8)
plt.plot(x_coord[xmsk],fnorm[xmsk],'kP',ms=8)

plt.show()

##### gen_mesh.py ######
__author__="Shu-Gang Ai, <email address hidden>"
__institution__="Tongji University"

class Meshing:

 def run(self, params):

  lengthX=params['lengthX']; numberX=params['numberX']; meshSizeX=params['meshSizeX'];
  lengthY1=params['lengthY1']; numberY1=params['numberY1']; meshSizeY1=params['meshSizeY1'];
  lengthY2=params['lengthY2']; numberY2=params['numberY2']; meshSizeY2=params['meshSizeY2'];

  rm_points={}
  rm_rectangles={}
  rm_line2={}
  rm_line2contacts={}

  for j in range(numberY1+1):
   for i in range(numberX+1):
    rm_points[int(j*(numberX+1)+(i+1))]=( \
     float(meshSizeX*i), \
     float(meshSizeY1*j))

  for j in range(numberY1+1,numberY1+numberY2+2):
   for i in range(numberX+1):
    rm_points[int(j*(numberX+1)+(i+1))]=( \
     float(meshSizeX*i), \
     float(lengthY1+meshSizeY2*(j-numberY1-1)))

  for j in range(0,numberY1):
   for i in range(0,numberX):
    rm_rectangles[int(j*numberX+i+1)]=( \
     int((numberX+1)*j+i+1), \
     int((numberX+1)*j+i+2), \
     int((numberX+1)*(j+1)+i+2), \
     int((numberX+1)*(j+1)+i+1))

  for j in range(numberY1+1,numberY1+numberY2+1):
   for i in range(0,numberX):
    rm_rectangles[int((j-1)*numberX+i+1)]=( \
     int((numberX+1)*j+i+1), \
     int((numberX+1)*j+i+2), \
     int((numberX+1)*(j+1)+i+2), \
     int((numberX+1)*(j+1)+i+1))

  for i in range(0,numberX):
   rm_line2contacts[int(numberX*(numberY1+numberY2)+i+1)]=( \
    int((numberX+1)*(numberY1+1)+i+2), \
    int((numberX+1)*(numberY1+1)+i+1), \
    int((numberX+1)*numberY1+i+2), \
    int((numberX+1)*numberY1+i+1))

  # line2-bottom
  for i in range(numberX):
   rm_line2[int(numberX*(numberY1+numberY2)+numberX+i+1)]=( \
    int(i+1), \
    int(i+2))

  # line2-right-1
  for i in range(numberY1):
   rm_line2[int(numberX*(numberY1+numberY2)+numberX+numberX+i+1)]=( \
    int((numberX+1)*(i+1)), \
    int((numberX+1)*(i+2)))

  # line2-right-2
  for i in range(numberY1+1,numberY1+numberY2+1):
   rm_line2[int(numberX*(numberY1+numberY2)+numberX+numberX+i)]=( \
    int((numberX+1)*(i+1)), \
    int((numberX+1)*(i+2)))

  # line2-top
  for i in range(numberX):
   rm_line2[int(numberX*(numberY1+numberY2)+numberX*2+(numberY1+numberY2)+i+1)]=( \
    int(((numberY1+numberY2)+2)*(numberX+1)-i), \
    int(((numberY1+numberY2)+2)*(numberX+1)-i-1))

  # line2-left-2
  for i in range(numberY2):
   rm_line2[int(numberX*(numberY1+numberY2)+numberX*3+(numberY1+numberY2)+i+1)]=( \
    int((numberX+1)*((numberY1+numberY2)+1-i)+1), \
    int((numberX+1)*((numberY1+numberY2)-i)+1))

  # line2-left-1
  for i in range(numberY2+1,(numberY1+numberY2)+1):
   rm_line2[int(numberX*(numberY1+numberY2)+numberX*3+(numberY1+numberY2)+i)]=( \
    int((numberX+1)*((numberY1+numberY2)+1-i)+1), \
    int((numberX+1)*((numberY1+numberY2)-i)+1))

  fh=open(params['paramID'] + '.fly', 'w')
  fh.write('mesh\n')

  fh.write('2D-nodes ' + str(len(rm_points)) + '\n')
  for (no, pt) in rm_points.items():
   fh.write('%d %d %d %f %f\n' % (no, no, 0, pt[0], pt[1]))

  fh.write('Rec4 ' + str(len(rm_rectangles)) + '\n')
  for (no, rec) in rm_rectangles.items():
   if no <= (numberX*numberY1):
    fh.write('%d %d %d %d %d %d\n' % (no, 1, rec[0], rec[1], rec[2], rec[3]))
   else:
    fh.write('%d %d %d %d %d %d\n' % (no, 2, rec[0], rec[1], rec[2], rec[3]))

  fh.write('Line2 ' + str(len(rm_line2)) + '\n')
  for (no, ln) in rm_line2.items():
   fh.write('%d %d %d %d\n' % (no, 3, ln[0], ln[1]))

  fh.write('Line2_Contact ' + str(len(rm_line2contacts)) + '\n')
  for (no, ce) in rm_line2contacts.items():
   fh.write('%d %d %d %d %d %d\n' % (no, 4, ce[0], ce[1], ce[2], ce[3]))

  fh.write('Point1 0\n')
  fh.write('\n\n\n')
  fh.close()

if __name__ == '__main__':

 Parameters={}

 lengthX=4.0; numberX=20; meshSizeX=lengthX/numberX;
 lengthY1=0.2; numberY1=1; meshSizeY1=lengthY1/numberY1;
 lengthY2=2.0; numberY2=10; meshSizeY2=lengthY2/numberY2;

 Parameters['lengthX']=lengthX; Parameters['numberX']=numberX; Parameters['meshSizeX']=meshSizeX
 Parameters['lengthY1']=lengthY1; Parameters['numberY1']=numberY1; Parameters['meshSizeY1']=meshSizeY1
 Parameters['lengthY2']=lengthY2; Parameters['numberY2']=numberY2; Parameters['meshSizeY2']=meshSizeY2

 Parameters['paramID']='mesh'
 params=Parameters.copy()
 meshing=Meshing()
 meshing.run(params)
 del meshing

Thanks for making this available.