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