/[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 100 - (show annotations)
Wed Dec 15 03:48:48 2004 UTC (15 years, 9 months ago) by jgs
File MIME type: text/x-python
File size: 5751 byte(s)
*** empty log message ***

1 # $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

  ViewVC Help
Powered by ViewVC 1.1.26