# Contents of /trunk/esys2/finley/test/python/ContactTest.py

Revision 104 - (show annotations)
Fri Dec 17 07:43:12 2004 UTC (14 years, 4 months ago) by jgs
File MIME type: text/x-python
File size: 5793 byte(s)
```*** empty log message ***

```
 1 # \$Id\$ 2 3 """tests a variety of functions in connection with contact elements""" 4 5 from esys.escript import * 6 from esys.linearPDEs import LinearPDE 7 import esys.finley 8 9 import math 10 numElements=4 11 numEquations=2 12 13 seed=1 14 rand_factor=sqrt(2.) 15 def randomNum(): 16 global seed 17 seed+=1 18 s=rand_factor*seed 19 return s-int(s) 20 21 22 def mkMesh(dim,order,numElm,onElem): 23 if dim==2: 24 ms1=Rectangle(numElm,numElm,order,l1=0.5,useElementsOnFace=onElem) 25 ms2=Rectangle(numElm,numElm,order,l1=0.5,useElementsOnFace=onElem) 26 ms2.setX(ms2.getX()+[0,0.5]) 27 else: 28 ms1=Brick(numElm,numElm,numElm,order,l2=0.5,useElementsOnFace=onElem) 29 ms2=Brick(numElm,numElm,numElm,order,l2=0.5,useElementsOnFace=onElem) 30 ms2.setX(ms2.getX()+[0,0,0.5]) 31 return JoinFaces([ms1,ms2]) 32 33 34 def mkCharateristicFunction(msh): 35 """returns a scalar function on nodes which is -1 below and 1 above the fault""" 36 options = { 37 "verbose" : True, 38 "reordering" : NO_REORDERING, 39 "tolerance" : 1.E-13, 40 "final_residual" : 0., 41 "iterative" : True, 42 "iterative_method" : BICGSTAB, 43 "preconditioner" : JACOBI, 44 "iter_max" : 4000, 45 "iter" : 0, 46 "drop_tolerance" : 1.10, 47 "drop_storage" : 2. 48 } 49 e=Function(msh) 50 d=msh.getDim() 51 x=e.getX()[d-1] 52 mypde=LinearPDE(msh) 53 mypde.setValue(D=1,Y=(x-0.5).whatPositive()) 54 return 2*mypde.getSolution(**options)-1 55 56 max_error_text="" 57 max_error=0. 58 for dim in [2,3]: 59 for order in [1,2]: 60 for onElem in [False,True]: 61 if onElem: 62 onElem_text=",elements on face" 63 else: 64 onElem_text="" 65 case="Contact: %dD, order %d%s"%(dim,order,onElem_text) 66 print case 67 msh=mkMesh(dim,order,numElements,onElem) 68 c0=FunctionOnContact0(msh) 69 c1=FunctionOnContact1(msh) 70 n=ContinuousFunction(msh) 71 # 72 # check the normal on the fault 0: 73 # 74 refNormal=Vector(0,what=c0) 75 if dim==3: 76 refNormal.setTaggedValue(100,[0,0,1]) 77 else: 78 refNormal.setTaggedValue(10,[0,1]) 79 error_norm=Lsup(c0.getNormal()-refNormal) 80 text=" %s : error normals 0= %e"%(case,error_norm) 81 print "@@@ %s"%text 82 if error_norm>max_error: max_error_text,max_error=text,error_norm 83 # 84 # check the normal on the fault 1: 85 # 86 refNormal=Vector(0,what=c1) 87 if dim==3: 88 refNormal.setTaggedValue(100,[0,0,-1]) 89 else: 90 refNormal.setTaggedValue(10,[0,-1]) 91 error_norm=Lsup(c1.getNormal()-refNormal) 92 text=" %s : error normals 1= %e"%(case,error_norm) 93 print "@@@ %s"%text 94 if error_norm>max_error: max_error_text,max_error=text,error_norm 95 # 96 # integration on 0: 97 # 98 g=c0.getX()[dim-1]**2 99 error_norm=abs(g.integrate()-0.25) 100 text=" %s : error integrate 0= %e"%(case,error_norm) 101 print "@@@ %s"%text 102 if error_norm>max_error: max_error_text,max_error=text,error_norm 103 # 104 # integration on 1: 105 # 106 g=c1.getX()[0]**2 107 error_norm=abs(g.integrate()-1./3.) 108 text=" %s : error integrate 1= %e"%(case,error_norm) 109 print "@@@ %s"%text 110 if error_norm>max_error: max_error_text,max_error=text,error_norm 111 # 112 # gradient: 113 # 114 if onElem: 115 x=n.getX() 116 char=(x[dim-1]-0.5).whereNegative() 117 foo=2*x[dim-1]*char+(10*x[dim-1]-4)*(1-char) 118 119 error_norm=Lsup(foo.grad(c0)[dim-1]-2) 120 text=" %s : error gradient on 0= %e"%(case,error_norm) 121 print "@@@ %s"%text 122 if error_norm>max_error: max_error_text,max_error=text,error_norm 123 124 error_norm=Lsup(foo.grad(c1)[dim-1]-10) 125 text=" %s : error gradient on 1= %e"%(case,error_norm) 126 print "@@@ %s"%text 127 if error_norm>max_error: max_error_text,max_error=text,error_norm 128 129 char=mkCharateristicFunction(msh) 130 for red in [False,True]: 131 if red: 132 case_red=",reduced" 133 else: 134 case_red="" 135 for ne in range(1,numEquations+1): 136 u_ex=Data(0,shape=(ne,),what=n) 137 for i in range(ne): 138 u_ex[i]=char*(i+1)*n.getX()[0] 139 u0=u_ex.interpolate(c0) 140 u1=u_ex.interpolate(c1) 141 if ne==1: 142 d=randomNum() 143 y=d*(u1-u0) 144 else: 145 d=Id(ne) 146 for i in range(ne): 147 for j in range(ne): 148 d[i,j]=randomNum() 149 y=Data(0,shape=(ne,),what=c0) 150 for i in range(ne): 151 for j in range(ne): 152 y[i]+=d[i,j]*(u1[j]-u0[j]) 153 mypde=linearPDE(d_contact=d,y_contact=y) 154 mypde.setReducedOrderForEquationsTo(red) 155 error_norm=Lsup(mypde.getResidual(u_ex))/Lsup(mypde.getRightHandSide()) 156 text=" %s%s, assembling %d equations: error = %e"%(case,case_red,ne,error_norm) 157 print "@@@ %s"%text 158 if error_norm>max_error: max_error_text,max_error=text,error_norm 159 160 161 print "maximal error for ",max_error_text 162

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision