/[escript]/trunk/finley/test/python/ContactTest.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 155 - (show annotations)
Wed Nov 9 02:02:19 2005 UTC (14 years, 3 months ago) by jgs
File MIME type: text/x-python
File size: 5720 byte(s)
move all directories from trunk/esys2 into trunk and remove esys2

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26