/[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 82 - (hide annotations)
Tue Oct 26 06:53:54 2004 UTC (15 years, 11 months ago) by jgs
File MIME type: text/x-python
File size: 5751 byte(s)
Initial revision

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