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

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