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 |
|