/[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 617 - (show annotations)
Wed Mar 22 02:58:17 2006 UTC (14 years, 4 months ago) by elspeth
File MIME type: text/x-python
File size: 6010 byte(s)
More copyright.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26