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

Revision 97 - (hide annotations)
Tue Dec 14 05:39:33 2004 UTC (15 years, 9 months ago) by jgs
File MIME type: text/x-python
File size: 5784 byte(s)
```*** empty log message ***

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

## Properties

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