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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26