/[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 97 - (show annotations)
Tue Dec 14 05:39:33 2004 UTC (15 years, 1 month ago) by jgs
File MIME type: text/x-python
File size: 5784 byte(s)
*** empty log message ***

1 # $Id$
2 """tests a variety of functions in connection with contact elements"""
3 from escript import *
4 from linearPDEs import LinearPDE
5 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 mypde=LinearPDE(D=1,Y=(x-0.5).whatPositive())
50 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

  ViewVC Help
Powered by ViewVC 1.1.26