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

Revision 82 - (hide annotations)
Tue Oct 26 06:53:54 2004 UTC (15 years, 11 months ago) by jgs
File MIME type: text/x-python
File size: 5751 byte(s)
```Initial revision

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

Properties

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