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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 149 - (hide annotations)
Thu Sep 1 03:31:39 2005 UTC (14 years, 6 months ago) by jgs
File MIME type: text/x-python
File size: 5720 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-01

1 jgs 82 # $Id$
2     """tests a variety of functions in connection with contact elements"""
3 jgs 149 from esys.escript import *
4     from esys.escript.linearPDEs import LinearPDE
5     from esys.escript.pdetools import Projector
6     from esys import finley
7 jgs 104
8 jgs 82 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 jgs 147 ms1=finley.Rectangle(numElm,numElm,order,l1=0.5,useElementsOnFace=onElem)
23     ms2=finley.Rectangle(numElm,numElm,order,l1=0.5,useElementsOnFace=onElem)
24 jgs 82 ms2.setX(ms2.getX()+[0,0.5])
25     else:
26 jgs 147 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 jgs 82 ms2.setX(ms2.getX()+[0,0,0.5])
29 jgs 147 return finley.JoinFaces([ms1,ms2])
30 jgs 82
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 jgs 104 mypde=LinearPDE(msh)
38 jgs 147 mypde.setValue(D=1.,Y=(x-0.5).wherePositive())
39     return 2*mypde.getSolution()-1
40 jgs 82
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 jgs 147 c0=FunctionOnContactZero(msh)
54     c1=FunctionOnContactOne(msh)
55 jgs 82 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 jgs 147 error_norm=abs(integrate(g)-0.25)
85 jgs 82 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 jgs 147 error_norm=abs(integrate(g)-1./3.)
93 jgs 82 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 jgs 147 # for red in [False,True]:
116     for red in [False]:
117 jgs 82 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 jgs 147 for i in range(ne):
124     u_ex[i]=(i+1)*n.getX()[0]*char
125 jgs 82 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 jgs 147 d=kronecker(ne)*1.
132 jgs 82 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 jgs 147 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 jgs 82 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