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

Annotation of /trunk/finley/test/python/run_models.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2850 - (hide annotations)
Fri Jan 15 07:14:25 2010 UTC (9 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 86304 byte(s)
the powerlaw is working again. tests are switched on again.
1 ksteube 1809
2     ########################################################
3 gross 1414 #
4 jfenwick 2548 # Copyright (c) 2003-2009 by University of Queensland
5 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 gross 1414 #
8 ksteube 1809 # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11 gross 1414 #
12 ksteube 1809 ########################################################
13 gross 1414
14 jfenwick 2549 __copyright__="""Copyright (c) 2003-2009 by University of Queensland
15 ksteube 1809 Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18 gross 1414 __license__="""Licensed under the Open Software License version 3.0
19 ksteube 1809 http://www.opensource.org/licenses/osl-3.0.php"""
20 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1809
22 gross 1414 import unittest
23     import tempfile
24    
25    
26    
27 gross 2850 VERBOSE=False # or True
28 gross 1414
29     from esys.escript import *
30 gross 2647 from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian, FaultSystem, DarcyFlow
31     from esys.escript.models import Mountains
32 gross 1414 from esys.finley import Rectangle, Brick
33 gross 2100
34 artak 2234 from math import pi
35 gross 2647 import numpy
36     import sys
37     import os
38     #====================================================================================================================
39     try:
40     FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
41     except KeyError:
42     FINLEY_WORKDIR='.'
43 artak 2234
44 gross 2100 #====================================================================================================================
45     class Test_StokesProblemCartesian2D(unittest.TestCase):
46 gross 1414 def setUp(self):
47 gross 2100 NE=6
48 gross 2208 self.TOL=1e-3
49 gross 2793 self.domain=Rectangle(NE,NE,order=-1)
50 gross 1414 def tearDown(self):
51     del self.domain
52 gross 2100 def test_PCG_P_0(self):
53 gross 1414 ETA=1.
54     P1=0.
55    
56     x=self.domain.getX()
57     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
58     mask=whereZero(x[0]) * [1.,1.] \
59     +whereZero(x[0]-1) * [1.,1.] \
60     +whereZero(x[1]) * [1.,0.] \
61     +whereZero(x[1]-1) * [1.,1.]
62    
63     sp=StokesProblemCartesian(self.domain)
64    
65     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
66     u0=(1-x[0])*x[0]*[0.,1.]
67 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
68 gross 2100 sp.setTolerance(self.TOL)
69 gross 2793 u,p=sp.solve(u0*mask,p0,verbose=VERBOSE,max_iter=100,usePCG=True)
70 gross 1414
71     error_v0=Lsup(u[0]-u0[0])
72     error_v1=Lsup(u[1]-u0[1])/0.25
73 gross 2100 error_p=Lsup(p+P1*x[0]*x[1])
74     self.failUnless(error_p<10*self.TOL, "pressure error too large.")
75     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
76     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
77 gross 1414
78 gross 2100 def test_PCG_P_small(self):
79 gross 1414 ETA=1.
80     P1=1.
81    
82     x=self.domain.getX()
83     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
84     mask=whereZero(x[0]) * [1.,1.] \
85     +whereZero(x[0]-1) * [1.,1.] \
86     +whereZero(x[1]) * [1.,0.] \
87     +whereZero(x[1]-1) * [1.,1.]
88    
89     sp=StokesProblemCartesian(self.domain)
90    
91     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
92     u0=(1-x[0])*x[0]*[0.,1.]
93 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
94 gross 2719 sp.setTolerance(self.TOL)
95 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
96 gross 1414 error_v0=Lsup(u[0]-u0[0])
97     error_v1=Lsup(u[1]-u0[1])/0.25
98 gross 2100 error_p=Lsup(P1*x[0]*x[1]+p)
99     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
100     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
101 gross 2156 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
102 gross 1414
103 gross 2100 def test_PCG_P_large(self):
104 gross 1414 ETA=1.
105     P1=1000.
106    
107     x=self.domain.getX()
108     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
109     mask=whereZero(x[0]) * [1.,1.] \
110     +whereZero(x[0]-1) * [1.,1.] \
111     +whereZero(x[1]) * [1.,0.] \
112     +whereZero(x[1]-1) * [1.,1.]
113    
114     sp=StokesProblemCartesian(self.domain)
115    
116     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
117     u0=(1-x[0])*x[0]*[0.,1.]
118 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
119 gross 2100 sp.setTolerance(self.TOL)
120 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
121 gross 2719 # u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
122 gross 1414
123     error_v0=Lsup(u[0]-u0[0])
124     error_v1=Lsup(u[1]-u0[1])/0.25
125 gross 2100 error_p=Lsup(P1*x[0]*x[1]+p)/P1
126     self.failUnless(error_p<10*self.TOL, "pressure error too large.")
127     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
128     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
129 gross 1414
130 gross 2100 def test_GMRES_P_0(self):
131 gross 1471 ETA=1.
132     P1=0.
133 gross 1414
134 gross 1471 x=self.domain.getX()
135     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
136     mask=whereZero(x[0]) * [1.,1.] \
137     +whereZero(x[0]-1) * [1.,1.] \
138     +whereZero(x[1]) * [1.,0.] \
139     +whereZero(x[1]-1) * [1.,1.]
140    
141     sp=StokesProblemCartesian(self.domain)
142    
143     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
144     u0=(1-x[0])*x[0]*[0.,1.]
145 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
146 gross 2100 sp.setTolerance(self.TOL)
147 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=50,usePCG=False,iter_restart=18)
148 gross 1471
149     error_v0=Lsup(u[0]-u0[0])
150     error_v1=Lsup(u[1]-u0[1])/0.25
151 gross 2156 error_p=Lsup(P1*x[0]*x[1]+p)
152 gross 2100 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
153     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
154     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
155 artak 1518
156 gross 2100 def test_GMRES_P_small(self):
157 gross 1471 ETA=1.
158     P1=1.
159    
160     x=self.domain.getX()
161     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
162     mask=whereZero(x[0]) * [1.,1.] \
163     +whereZero(x[0]-1) * [1.,1.] \
164     +whereZero(x[1]) * [1.,0.] \
165     +whereZero(x[1]-1) * [1.,1.]
166    
167     sp=StokesProblemCartesian(self.domain)
168    
169     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
170     u0=(1-x[0])*x[0]*[0.,1.]
171 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
172 gross 2719 sp.setTolerance(self.TOL)
173 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=20,usePCG=False)
174 gross 1471
175     error_v0=Lsup(u[0]-u0[0])
176     error_v1=Lsup(u[1]-u0[1])/0.25
177 gross 2156 error_p=Lsup(P1*x[0]*x[1]+p)
178 gross 2719
179 gross 2100 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
180     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
181     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
182 gross 1471
183 gross 2100 def test_GMRES_P_large(self):
184 gross 1471 ETA=1.
185     P1=1000.
186    
187     x=self.domain.getX()
188     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
189     mask=whereZero(x[0]) * [1.,1.] \
190     +whereZero(x[0]-1) * [1.,1.] \
191     +whereZero(x[1]) * [1.,0.] \
192     +whereZero(x[1]-1) * [1.,1.]
193    
194     sp=StokesProblemCartesian(self.domain)
195    
196     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
197     u0=(1-x[0])*x[0]*[0.,1.]
198 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
199 gross 2100 sp.setTolerance(self.TOL)
200 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)
201 gross 1471
202     error_v0=Lsup(u[0]-u0[0])
203     error_v1=Lsup(u[1]-u0[1])/0.25
204 gross 2156 error_p=Lsup(P1*x[0]*x[1]+p)/P1
205 gross 2100 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
206     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
207     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
208     #====================================================================================================================
209     class Test_StokesProblemCartesian3D(unittest.TestCase):
210 gross 1414 def setUp(self):
211 gross 2100 NE=6
212 gross 2208 self.TOL=1e-4
213 gross 2793 self.domain=Brick(NE,NE,NE,order=-1)
214 gross 1414 def tearDown(self):
215     del self.domain
216 gross 2100 def test_PCG_P_0(self):
217 gross 1414 ETA=1.
218     P1=0.
219    
220     x=self.domain.getX()
221     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
222     x=self.domain.getX()
223     mask=whereZero(x[0]) * [1.,1.,1.] \
224     +whereZero(x[0]-1) * [1.,1.,1.] \
225 gross 2156 +whereZero(x[1]) * [1.,0.,1.] \
226 gross 1414 +whereZero(x[1]-1) * [1.,1.,1.] \
227     +whereZero(x[2]) * [1.,1.,0.] \
228     +whereZero(x[2]-1) * [1.,1.,1.]
229    
230    
231     sp=StokesProblemCartesian(self.domain)
232    
233     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
234     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
235 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
236 gross 2156 sp.setTolerance(self.TOL)
237 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
238 gross 1414
239     error_v0=Lsup(u[0]-u0[0])
240     error_v1=Lsup(u[1]-u0[1])
241     error_v2=Lsup(u[2]-u0[2])/0.25**2
242 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
243 gross 2100 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
244     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
245     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
246     self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
247 artak 1518
248 gross 2100 def test_PCG_P_small(self):
249 gross 1414 ETA=1.
250     P1=1.
251    
252     x=self.domain.getX()
253     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
254     mask=whereZero(x[0]) * [1.,1.,1.] \
255     +whereZero(x[0]-1) * [1.,1.,1.] \
256 gross 2156 +whereZero(x[1]) * [1.,0.,1.] \
257 gross 1414 +whereZero(x[1]-1) * [1.,1.,1.] \
258     +whereZero(x[2]) * [1.,1.,0.] \
259     +whereZero(x[2]-1) * [1.,1.,1.]
260    
261    
262     sp=StokesProblemCartesian(self.domain)
263    
264     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
265     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
266 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
267 gross 2719 sp.setTolerance(self.TOL)
268 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
269 gross 1414 error_v0=Lsup(u[0]-u0[0])
270     error_v1=Lsup(u[1]-u0[1])
271     error_v2=Lsup(u[2]-u0[2])/0.25**2
272 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
273 gross 2100 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
274     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
275     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
276     self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
277 gross 2445
278 gross 2100 def test_PCG_P_large(self):
279 gross 1414 ETA=1.
280     P1=1000.
281    
282     x=self.domain.getX()
283     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
284     mask=whereZero(x[0]) * [1.,1.,1.] \
285     +whereZero(x[0]-1) * [1.,1.,1.] \
286 gross 2156 +whereZero(x[1]) * [1.,0.,1.] \
287 gross 1414 +whereZero(x[1]-1) * [1.,1.,1.] \
288     +whereZero(x[2]) * [1.,1.,0.] \
289     +whereZero(x[2]-1) * [1.,1.,1.]
290    
291    
292     sp=StokesProblemCartesian(self.domain)
293    
294     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
295     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
296 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
297 gross 2156 sp.setTolerance(self.TOL)
298 gross 2719 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
299 gross 1414
300     error_v0=Lsup(u[0]-u0[0])
301     error_v1=Lsup(u[1]-u0[1])
302     error_v2=Lsup(u[2]-u0[2])/0.25**2
303 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
304 gross 2100 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
305     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
306     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
307     self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
308 gross 1414
309 gross 2100 def test_GMRES_P_0(self):
310 gross 1471 ETA=1.
311     P1=0.
312    
313     x=self.domain.getX()
314     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
315     x=self.domain.getX()
316     mask=whereZero(x[0]) * [1.,1.,1.] \
317     +whereZero(x[0]-1) * [1.,1.,1.] \
318     +whereZero(x[1]) * [1.,1.,1.] \
319     +whereZero(x[1]-1) * [1.,1.,1.] \
320     +whereZero(x[2]) * [1.,1.,0.] \
321     +whereZero(x[2]-1) * [1.,1.,1.]
322    
323    
324     sp=StokesProblemCartesian(self.domain)
325    
326     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
327     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
328 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
329 gross 2156 sp.setTolerance(self.TOL)
330 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False,iter_restart=20)
331 gross 1471
332     error_v0=Lsup(u[0]-u0[0])
333     error_v1=Lsup(u[1]-u0[1])
334     error_v2=Lsup(u[2]-u0[2])/0.25**2
335 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
336 gross 2100 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
337     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
338     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
339     self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
340     def test_GMRES_P_small(self):
341 gross 1471 ETA=1.
342     P1=1.
343    
344     x=self.domain.getX()
345     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
346     mask=whereZero(x[0]) * [1.,1.,1.] \
347     +whereZero(x[0]-1) * [1.,1.,1.] \
348     +whereZero(x[1]) * [1.,1.,1.] \
349     +whereZero(x[1]-1) * [1.,1.,1.] \
350     +whereZero(x[2]) * [1.,1.,0.] \
351     +whereZero(x[2]-1) * [1.,1.,1.]
352    
353    
354     sp=StokesProblemCartesian(self.domain)
355    
356     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
357     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
358 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
359 gross 2719 sp.setTolerance(self.TOL/10)
360 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)
361 gross 1471
362     error_v0=Lsup(u[0]-u0[0])
363     error_v1=Lsup(u[1]-u0[1])
364     error_v2=Lsup(u[2]-u0[2])/0.25**2
365 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
366 gross 2100 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
367     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
368     self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
369 gross 2251 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
370 gross 2100 def test_GMRES_P_large(self):
371 gross 1471 ETA=1.
372     P1=1000.
373    
374     x=self.domain.getX()
375     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
376     mask=whereZero(x[0]) * [1.,1.,1.] \
377     +whereZero(x[0]-1) * [1.,1.,1.] \
378 gross 2156 +whereZero(x[1]) * [1.,0.,1.] \
379 gross 1471 +whereZero(x[1]-1) * [1.,1.,1.] \
380     +whereZero(x[2]) * [1.,1.,0.] \
381     +whereZero(x[2]-1) * [1.,1.,1.]
382    
383    
384     sp=StokesProblemCartesian(self.domain)
385    
386     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
387     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
388 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
389 gross 2156 sp.setTolerance(self.TOL)
390 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=False)
391 gross 1471
392     error_v0=Lsup(u[0]-u0[0])
393     error_v1=Lsup(u[1]-u0[1])
394     error_v2=Lsup(u[2]-u0[2])/0.25**2
395 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
396 gross 2100 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
397     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
398     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
399     self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
400     #====================================================================================================================
401     class Test_Darcy(unittest.TestCase):
402     # this is a simple test for the darcy flux problem
403     #
404     #
405     # p = 1/k * ( 1/2* (fb-f0)/xb* x **2 + f0 * x - ub*x ) + p0
406     #
407     # with f = (fb-f0)/xb* x + f0
408     #
409     # u = f - k * p,x = ub
410     #
411     # we prescribe pressure at x=x0=0 to p0
412     #
413     # if we prescribe the pressure on the bottom x=xb we set
414     #
415     # pb= 1/k * ( 1/2* (fb-f0)* xb + f0 * xb - ub*xb ) + p0 = 1/k * ((fb+f0)/2 - ub ) * xb + p0
416     #
417     # which leads to ub = (fb+f0)/2-k*(pb-p0)/xb
418     #
419     def rescaleDomain(self):
420     x=self.dom.getX().copy()
421     for i in xrange(self.dom.getDim()):
422     x_inf=inf(x[i])
423     x_sup=sup(x[i])
424     if i == self.dom.getDim()-1:
425     x[i]=-self.WIDTH*(x[i]-x_sup)/(x_inf-x_sup)
426     else:
427     x[i]=self.WIDTH*(x[i]-x_inf)/(x_sup-x_inf)
428     self.dom.setX(x)
429     def getScalarMask(self,include_bottom=True):
430     x=self.dom.getX().copy()
431     x_inf=inf(x[self.dom.getDim()-1])
432     x_sup=sup(x[self.dom.getDim()-1])
433     out=whereZero(x[self.dom.getDim()-1]-x_sup)
434     if include_bottom: out+=whereZero(x[self.dom.getDim()-1]-x_inf)
435     return wherePositive(out)
436     def getVectorMask(self,include_bottom=True):
437     x=self.dom.getX().copy()
438     out=Vector(0.,Solution(self.dom))
439     for i in xrange(self.dom.getDim()):
440     x_inf=inf(x[i])
441     x_sup=sup(x[i])
442     if i != self.dom.getDim()-1: out[i]+=whereZero(x[i]-x_sup)
443     if i != self.dom.getDim()-1 or include_bottom: out[i]+=whereZero(x[i]-x_inf)
444     return wherePositive(out)
445 gross 1471
446 gross 2100 def setSolutionFixedBottom(self, p0, pb, f0, fb, k):
447     d=self.dom.getDim()
448     x=self.dom.getX()[d-1]
449     xb=inf(x)
450     u=Vector(0.,Solution(self.dom))+kronecker(d)[d-1]*((f0+fb)/2.-k*(pb-p0)/xb)
451     p=1./k*((fb-f0)/(xb*2.)* x**2 - (fb-f0)/2.*x)+(pb-p0)/xb*x + p0
452     f= ((fb-f0)/xb* x + f0)*kronecker(Function(self.dom))[d-1]
453     return u,p,f
454    
455     def testConstF_FixedBottom_smallK(self):
456     k=1.e-10
457     mp=self.getScalarMask(include_bottom=True)
458     mv=self.getVectorMask(include_bottom=False)
459     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
460     p=p_ref*mp
461     u=u_ref*mv
462     df=DarcyFlow(self.dom)
463     df.setValue(g=f,
464     location_of_fixed_pressure=mp,
465     location_of_fixed_flux=mv,
466     permeability=Scalar(k,Function(self.dom)))
467 gross 2264 df.setTolerance(rtol=self.TOL)
468     v,p=df.solve(u_ref,p, max_iter=100, verbose=VERBOSE)
469 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
470     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
471     def testConstF_FixedBottom_mediumK(self):
472     k=1.
473     mp=self.getScalarMask(include_bottom=True)
474     mv=self.getVectorMask(include_bottom=False)
475     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
476     p=p_ref*mp
477     u=u_ref*mv
478     df=DarcyFlow(self.dom)
479     df.setValue(g=f,
480     location_of_fixed_pressure=mp,
481     location_of_fixed_flux=mv,
482     permeability=Scalar(k,Function(self.dom)))
483 gross 2264 df.setTolerance(rtol=self.TOL)
484     v,p=df.solve(u,p,max_iter=100, verbose=VERBOSE )
485 gross 2208 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
486 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
487 artak 1518
488 gross 2100 def testConstF_FixedBottom_largeK(self):
489     k=1.e10
490     mp=self.getScalarMask(include_bottom=True)
491     mv=self.getVectorMask(include_bottom=False)
492     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
493     p=p_ref*mp
494     u=u_ref*mv
495     df=DarcyFlow(self.dom)
496     df.setValue(g=f,
497     location_of_fixed_pressure=mp,
498     location_of_fixed_flux=mv,
499     permeability=Scalar(k,Function(self.dom)))
500 gross 2264 df.setTolerance(rtol=self.TOL)
501     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
502 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
503     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
504 artak 1518
505 gross 2100 def testVarioF_FixedBottom_smallK(self):
506     k=1.e-10
507     mp=self.getScalarMask(include_bottom=True)
508     mv=self.getVectorMask(include_bottom=False)
509     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
510     p=p_ref*mp
511     u=u_ref*mv
512     df=DarcyFlow(self.dom)
513     df.setValue(g=f,
514     location_of_fixed_pressure=mp,
515     location_of_fixed_flux=mv,
516     permeability=Scalar(k,Function(self.dom)))
517 gross 2264 df.setTolerance(rtol=self.TOL)
518     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
519 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
520     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
521 artak 1518
522 gross 2100 def testVarioF_FixedBottom_mediumK(self):
523     k=1.
524     mp=self.getScalarMask(include_bottom=True)
525     mv=self.getVectorMask(include_bottom=False)
526     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
527     p=p_ref*mp
528     u=u_ref*mv
529     df=DarcyFlow(self.dom)
530     df.setValue(g=f,
531     location_of_fixed_pressure=mp,
532     location_of_fixed_flux=mv,
533     permeability=Scalar(k,Function(self.dom)))
534 gross 2264 df.setTolerance(rtol=self.TOL)
535     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
536 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
537     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
538 artak 1518
539 gross 2100 def testVarioF_FixedBottom_largeK(self):
540     k=1.e10
541     mp=self.getScalarMask(include_bottom=True)
542     mv=self.getVectorMask(include_bottom=False)
543     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
544     p=p_ref*mp
545     u=u_ref*mv
546     df=DarcyFlow(self.dom)
547     df.setValue(g=f,
548     location_of_fixed_pressure=mp,
549     location_of_fixed_flux=mv,
550     permeability=Scalar(k,Function(self.dom)))
551 gross 2264 df.setTolerance(rtol=self.TOL)
552     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
553 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
554     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
555 artak 1518
556 gross 2100 def testConstF_FreeBottom_smallK(self):
557     k=1.e-10
558     mp=self.getScalarMask(include_bottom=False)
559     mv=self.getVectorMask(include_bottom=True)
560     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
561     p=p_ref*mp
562     u=u_ref*mv
563     df=DarcyFlow(self.dom)
564     df.setValue(g=f,
565 gross 2208 location_of_fixed_pressure=mp,
566 gross 2100 location_of_fixed_flux=mv,
567     permeability=Scalar(k,Function(self.dom)))
568 gross 2264 df.setTolerance(rtol=self.TOL)
569     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
570 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
571     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
572 artak 1518
573 gross 2100 def testConstF_FreeBottom_mediumK(self):
574     k=1.
575     mp=self.getScalarMask(include_bottom=False)
576     mv=self.getVectorMask(include_bottom=True)
577     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
578     p=p_ref*mp
579     u=u_ref*mv
580     df=DarcyFlow(self.dom)
581     df.setValue(g=f,
582     location_of_fixed_pressure=mp,
583     location_of_fixed_flux=mv,
584     permeability=Scalar(k,Function(self.dom)))
585 gross 2264 df.setTolerance(rtol=self.TOL)
586     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
587 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
588     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
589 artak 1518
590 gross 2100 def testConstF_FreeBottom_largeK(self):
591     k=1.e10
592     mp=self.getScalarMask(include_bottom=False)
593     mv=self.getVectorMask(include_bottom=True)
594     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
595     p=p_ref*mp
596     u=u_ref*mv
597     df=DarcyFlow(self.dom)
598     df.setValue(g=f,
599     location_of_fixed_pressure=mp,
600     location_of_fixed_flux=mv,
601     permeability=Scalar(k,Function(self.dom)))
602 gross 2264 df.setTolerance(rtol=self.TOL)
603     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
604 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
605     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
606 artak 1518
607 gross 2100 def testVarioF_FreeBottom_smallK(self):
608     k=1.e-10
609     mp=self.getScalarMask(include_bottom=False)
610     mv=self.getVectorMask(include_bottom=True)
611     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
612     p=p_ref*mp
613     u=u_ref*mv
614     df=DarcyFlow(self.dom)
615     df.setValue(g=f,
616     location_of_fixed_pressure=mp,
617     location_of_fixed_flux=mv,
618     permeability=Scalar(k,Function(self.dom)))
619 gross 2264 df.setTolerance(rtol=self.TOL)
620     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
621 gross 2208 self.failUnless(Lsup(v-u_ref)<self.TOL*25.*Lsup(u_ref), "flux error too big.") # 25 because of disc. error.
622     self.failUnless(Lsup(p-p_ref)<self.TOL*25.*Lsup(p_ref), "pressure error too big.")
623 artak 1518
624 gross 2100 def testVarioF_FreeBottom_mediumK(self):
625     k=1.
626     mp=self.getScalarMask(include_bottom=False)
627     mv=self.getVectorMask(include_bottom=True)
628     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
629     p=p_ref*mp
630     u=u_ref*mv
631     df=DarcyFlow(self.dom)
632     df.setValue(g=f,
633     location_of_fixed_pressure=mp,
634     location_of_fixed_flux=mv,
635     permeability=Scalar(k,Function(self.dom)))
636 gross 2264 df.setTolerance(rtol=self.TOL)
637     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
638 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
639     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
640 artak 1518
641 gross 2100 def testVarioF_FreeBottom_largeK(self):
642     k=1.e10
643     mp=self.getScalarMask(include_bottom=False)
644     mv=self.getVectorMask(include_bottom=True)
645     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
646     p=p_ref*mp
647     u=u_ref*mv
648     df=DarcyFlow(self.dom)
649     df.setValue(g=f,
650     location_of_fixed_pressure=mp,
651     location_of_fixed_flux=mv,
652     permeability=Scalar(k,Function(self.dom)))
653 gross 2264 df.setTolerance(rtol=self.TOL)
654     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
655 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
656     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
657 artak 1518
658 gross 2100 class Test_Darcy2D(Test_Darcy):
659 gross 2208 TOL=1e-4
660 gross 2100 WIDTH=1.
661     def setUp(self):
662 gross 2208 NE=40 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
663 gross 2100 self.dom = Rectangle(NE,NE)
664     self.rescaleDomain()
665     def tearDown(self):
666     del self.dom
667     class Test_Darcy3D(Test_Darcy):
668     TOL=1e-4
669     WIDTH=1.
670     def setUp(self):
671     NE=25 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
672     self.dom = Brick(NE,NE,NE)
673     self.rescaleDomain()
674     def tearDown(self):
675     del self.dom
676 artak 1518
677 gross 2100
678 artak 2234 class Test_Mountains3D(unittest.TestCase):
679     def setUp(self):
680     NE=16
681     self.TOL=1e-4
682     self.domain=Brick(NE,NE,NE,order=1,useFullElementOrder=True)
683     def tearDown(self):
684     del self.domain
685     def test_periodic(self):
686     EPS=0.01
687 gross 2100
688 artak 2234 x=self.domain.getX()
689     v = Vector(0.0, Solution(self.domain))
690     a0=1
691     a1=1
692     n0=2
693     n1=2
694     n2=0.5
695     a2=-(a0*n0+a1*n1)/n2
696     v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])* cos(pi*n2*x[2])
697     v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])* cos(pi*n2*x[2])
698     v[2]=a2*cos(pi*n0*x[0])* cos(pi*n1*x[1])* sin(pi*n2*x[2])
699    
700 gross 2563 mts=Mountains(self.domain,eps=EPS)
701     mts.setVelocity(v)
702     Z=mts.update()
703 artak 2234
704 gross 2564 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[2]-1.)))
705 artak 2234 self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
706    
707     class Test_Mountains2D(unittest.TestCase):
708     def setUp(self):
709     NE=16
710     self.TOL=1e-4
711     self.domain=Rectangle(NE,NE,order=1,useFullElementOrder=True)
712     def tearDown(self):
713     del self.domain
714     def test_periodic(self):
715     EPS=0.01
716    
717     x=self.domain.getX()
718     v = Vector(0.0, Solution(self.domain))
719     a0=1
720     n0=1
721     n1=0.5
722     a1=-(a0*n0)/n1
723     v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])
724     v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])
725    
726     H_t=Scalar(0.0, Solution(self.domain))
727 gross 2563 mts=Mountains(self.domain,eps=EPS)
728     mts.setVelocity(v)
729     Z=mts.update()
730 artak 2234
731 gross 2564 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[1]-1.)))
732 artak 2234 self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
733    
734 gross 2301
735    
736     class Test_Rheologies(unittest.TestCase):
737     """
738     this is the program used to generate the powerlaw tests:
739    
740     TAU_Y=100.
741     N=10
742     M=5
743    
744     def getE(tau):
745     if tau<=TAU_Y:
746     return 1./(0.5+20*sqrt(tau))
747     else:
748     raise ValueError,"out of range."
749     tau=[ i* (TAU_Y/N) for i in range(N+1)] + [TAU_Y for i in range(M)]
750     e=[ tau[i]/getE(tau[i]) for i in range(N+1)]
751     e+= [ (i+1+j)* (max(e)/N) for j in range(M) ]
752    
753     print tau
754     print e
755     """
756     TOL=1.e-8
757     def test_PowerLaw_Init(self):
758     pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
759    
760     self.failUnlessEqual(pl.getNumMaterials(),3,"num materials is wrong")
761     self.failUnlessEqual(pl.validMaterialId(0),True,"material id 0 not found")
762     self.failUnlessEqual(pl.validMaterialId(1),True,"material id 1 not found")
763     self.failUnlessEqual(pl.validMaterialId(2),True,"material id 2 not found")
764     self.failUnlessEqual(pl.validMaterialId(3),False,"material id 3 not found")
765    
766     self.failUnlessEqual(pl.getFriction(),None,"initial friction wrong.")
767     self.failUnlessEqual(pl.getTauY(),None,"initial tau y wrong.")
768     pl.setDruckerPragerLaw(tau_Y=10,friction=3)
769     self.failUnlessEqual(pl.getFriction(),3,"friction wrong.")
770     self.failUnlessEqual(pl.getTauY(),10,"tau y wrong.")
771    
772     self.failUnlessEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
773     pl.setElasticShearModulus(1000)
774     self.failUnlessEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
775    
776     e=pl.getEtaN()
777     self.failUnlessEqual(len(e),3,"initial length of etaN is wrong.")
778     self.failUnlessEqual(e,[None, None, None],"initial etaN are wrong.")
779     self.failUnlessEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
780     self.failUnlessEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
781     self.failUnlessEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
782     self.failUnlessRaises(ValueError, pl.getEtaN, 3)
783    
784     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
785     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
786     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
787    
788     pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
789     self.failUnlessEqual(pl.getPower(),[1,2,3],"powers are wrong.")
790     self.failUnlessEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
791     self.failUnlessEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
792    
793     pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
794     self.failUnlessEqual(pl.getPower(),[4,2,3],"powers are wrong.")
795     self.failUnlessEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
796     self.failUnlessEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
797    
798     self.failUnlessRaises(ValueError,pl.getPower,-1)
799     self.failUnlessEqual(pl.getPower(0),4,"power 0 is wrong.")
800     self.failUnlessEqual(pl.getPower(1),2,"power 1 is wrong.")
801     self.failUnlessEqual(pl.getPower(2),3,"power 2 is wrong.")
802     self.failUnlessRaises(ValueError,pl.getPower,3)
803    
804     self.failUnlessRaises(ValueError,pl.getTauT,-1)
805     self.failUnlessEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
806     self.failUnlessEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
807     self.failUnlessEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
808     self.failUnlessRaises(ValueError,pl.getTauT,3)
809    
810     self.failUnlessRaises(ValueError,pl.getEtaN,-1)
811     self.failUnlessEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
812     self.failUnlessEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
813     self.failUnlessEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
814     self.failUnlessRaises(ValueError,pl.getEtaN,3)
815    
816     def checkResult(self,id,gamma_dot_, eta, tau_ref):
817     self.failUnless(eta>=0,"eta needs to be positive (test %s)"%id)
818     error=abs(gamma_dot_*eta-tau_ref)
819     self.failUnless(error<=self.TOL*tau_ref,"eta is wrong: error = gamma_dot_*eta-tau_ref = %s * %s - %s = %s (test %s)"%(gamma_dot_,eta,tau_ref,error,id))
820    
821     def test_PowerLaw_Linear(self):
822     taus= [0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
823     gamma_dot_s=[0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0, 75.0]
824     pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
825     pl.setDruckerPragerLaw(tau_Y=100.)
826     pl.setPowerLaw(eta_N=2.)
827     pl.setEtaTolerance(self.TOL)
828     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
829    
830     def test_PowerLaw_QuadLarge(self):
831     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
832 gross 2438 gamma_dot_s=[0.0, 405.0, 1610.0, 3615.0, 6420.0, 10025.0, 14430.0, 19635.0, 25640.0, 32445.0, 40050.0, 44055.0, 48060.0, 52065.0, 56070.0, 60075.0]
833 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
834     pl.setDruckerPragerLaw(tau_Y=100.)
835     pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2])
836     pl.setEtaTolerance(self.TOL)
837     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
838    
839     def test_PowerLaw_QuadSmall(self):
840     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
841 gross 2438 gamma_dot_s=[0.0, 5.4, 11.6, 18.6, 26.4, 35.0, 44.4, 54.6, 65.6, 77.4, 90.0, 99.0, 108.0, 117.0, 126.0, 135.0]
842 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
843     pl.setDruckerPragerLaw(tau_Y=100.)
844     pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2])
845     pl.setEtaTolerance(self.TOL)
846     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
847    
848     def test_PowerLaw_CubeLarge(self):
849     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
850 gross 2438 gamma_dot_s=[0.0, 8.90625, 41.25, 120.46875, 270.0, 513.28125, 873.75, 1374.84375, 2040.0, 2892.65625, 3956.25, 4351.875, 4747.5, 5143.125, 5538.75, 5934.375]
851 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
852     pl.setDruckerPragerLaw(tau_Y=100.)
853     pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3])
854     pl.setEtaTolerance(self.TOL)
855     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
856    
857     def test_PowerLaw_CubeSmall(self):
858     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
859 gross 2438 gamma_dot_s=[0.0, 5.0390625, 10.3125, 16.0546875, 22.5, 29.8828125, 38.4375, 48.3984375, 60.0, 73.4765625, 89.0625, 97.96875, 106.875, 115.78125, 124.6875, 133.59375]
860 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
861     pl.setDruckerPragerLaw(tau_Y=100.)
862     pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3])
863     pl.setEtaTolerance(self.TOL)
864     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
865    
866     def test_PowerLaw_QuadLarge_CubeLarge(self):
867     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
868 gross 2438 gamma_dot_s=[0.0, 408.90625, 1641.25, 3720.46875, 6670.0, 10513.28125, 15273.75, 20974.84375, 27640.000000000004, 35292.65625, 43956.25, 48351.875, 52747.5, 57143.125, 61538.75, 65934.375]
869    
870 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
871     pl.setDruckerPragerLaw(tau_Y=100.)
872     pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
873     pl.setEtaTolerance(self.TOL)
874     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
875    
876     def test_PowerLaw_QuadLarge_CubeSmall(self):
877     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
878 gross 2438 gamma_dot_s=[0.0, 405.0390625, 1610.3125, 3616.0546875, 6422.5, 10029.8828125, 14438.4375, 19648.3984375, 25660.0, 32473.4765625, 40089.0625, 44097.96875, 48106.875, 52115.78125, 56124.6875, 60133.59375]
879    
880 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
881     pl.setDruckerPragerLaw(tau_Y=100.)
882     pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
883     pl.setEtaTolerance(self.TOL)
884     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
885    
886     def test_PowerLaw_QuadSmall_CubeLarge(self):
887     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
888 gross 2438 gamma_dot_s=[0.0, 9.30625, 42.85, 124.06875, 276.4, 523.28125, 888.15, 1394.44375, 2065.6, 2925.05625, 3996.25, 4395.875, 4795.5, 5195.125, 5594.75, 5994.375]
889 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
890     pl.setDruckerPragerLaw(tau_Y=100.)
891     pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
892     pl.setEtaTolerance(self.TOL)
893     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
894    
895     def test_PowerLaw_QuadSmall_CubeSmall(self):
896     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
897 gross 2438 gamma_dot_s=[0.0, 5.4390625, 11.9125, 19.6546875, 28.9, 39.8828125, 52.8375, 67.9984375, 85.6, 105.8765625, 129.0625, 141.96875, 154.875, 167.78125, 180.6875, 193.59375]
898 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
899     pl.setDruckerPragerLaw(tau_Y=100.)
900     pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
901     pl.setEtaTolerance(self.TOL)
902     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
903    
904     def test_PowerLaw_withShear(self):
905     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
906     gamma_dot_s=[0.0, 15.0, 30.0, 45.0, 60.0, 75.0, 90.0, 105.0, 120.0, 135.0, 150.0, 165.0, 180.0, 195.0, 210.0, 225.0]
907     pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
908     pl.setDruckerPragerLaw(tau_Y=100.)
909     pl.setPowerLaw(eta_N=2.)
910     pl.setElasticShearModulus(3.)
911     dt=1./3.
912     pl.setEtaTolerance(self.TOL)
913     self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
914     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
915    
916 gross 2432 class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase):
917 gross 2850 TOL=1.e-5
918     VERBOSE=False # or True
919 gross 2432 A=1.
920     P_max=100
921 gross 2794 NE=2*getMPISizeWorld()
922 gross 2432 tau_Y=10.
923     N_dt=10
924    
925     # material parameter:
926     tau_1=5.
927     tau_2=5.
928     eta_0=100.
929     eta_1=50.
930     eta_2=400.
931     N_1=2.
932     N_2=3.
933     def getReference(self, t):
934    
935     B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
936     x=self.dom.getX()
937    
938     s_00=min(self.A*t,B)
939     tau=sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)*abs(s_00)
940 gross 2438 inv_eta= 1./self.eta_0 + 1./self.eta_1*(tau/self.tau_1)**(self.N_1-1.) + 1./self.eta_2*(tau/self.tau_2)**(self.N_2-1.)
941 gross 2432
942     alpha=0.5*inv_eta*s_00
943     if s_00 <= B and self.mu !=None: alpha+=1./(2*self.mu)*self.A
944     u_ref=x*alpha
945     u_ref[self.dom.getDim()-1]=(1.-x[self.dom.getDim()-1])*alpha*(self.dom.getDim()-1)
946     sigma_ref=kronecker(self.dom)*s_00
947     sigma_ref[self.dom.getDim()-1,self.dom.getDim()-1]=-s_00*(self.dom.getDim()-1)
948    
949     p_ref=self.P_max
950     for d in range(self.dom.getDim()): p_ref=p_ref*x[d]
951     p_ref-=integrate(p_ref)/vol(self.dom)
952     return u_ref, sigma_ref, p_ref
953    
954     def runIt(self, free=None):
955     x=self.dom.getX()
956     B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
957     dt=B/int(self.N_dt/2)
958     if self.VERBOSE: print "dt =",dt
959     if self.latestart:
960     t=dt
961     else:
962     t=0
963     v,s,p=self.getReference(t)
964    
965     mod=IncompressibleIsotropicFlowCartesian(self.dom, stress=s, v=v, p=p, t=t, numMaterials=3, verbose=self.VERBOSE)
966     mod.setDruckerPragerLaw(tau_Y=self.tau_Y,friction=None)
967     mod.setElasticShearModulus(self.mu)
968     mod.setPowerLaws([self.eta_0, self.eta_1, self.eta_2], [ 1., self.tau_1, self.tau_2], [1.,self.N_1,self.N_2])
969     mod.setTolerance(self.TOL)
970 gross 2850 mod.setEtaTolerance(self.TOL*0.1)
971 gross 2432
972     BF=Vector(self.P_max,Function(self.dom))
973     for d in range(self.dom.getDim()):
974     for d2 in range(self.dom.getDim()):
975     if d!=d2: BF[d]*=x[d2]
976     v_mask=Vector(0,Solution(self.dom))
977     if free==None:
978     for d in range(self.dom.getDim()):
979     v_mask[d]=whereZero(x[d])+whereZero(x[d]-1.)
980     else:
981     for d in range(self.dom.getDim()):
982     if d == self.dom.getDim()-1:
983     v_mask[d]=whereZero(x[d]-1.)
984     else:
985     v_mask[d]=whereZero(x[d])
986     mod.setExternals(F=BF,fixed_v_mask=v_mask)
987    
988     n=self.dom.getNormal()
989     N_t=0
990     errors=[]
991     while N_t < self.N_dt:
992     t_ref=t+dt
993     v_ref, s_ref,p_ref=self.getReference(t_ref)
994     mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref)
995 gross 2794 # mod.update(dt, eta_iter_max=10, iter_max=50, verbose=self.VERBOSE, usePCG=True, max_correction_steps=30) new version
996     mod.update(dt, iter_max=50, verbose=self.VERBOSE, usePCG=True)
997 gross 2432 self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)
998     t+=dt
999     N_t+=1
1000    
1001     def check(self,N_t,mod,t_ref,v_ref, s_ref,p_ref):
1002     p=mod.getPressure()
1003     p-=integrate(p)/vol(self.dom)
1004     error_p=Lsup(mod.getPressure()-p_ref)/Lsup(p_ref)
1005     error_s=Lsup(mod.getDeviatoricStress()-s_ref)/Lsup(s_ref)
1006     error_v=Lsup(mod.getVelocity()-v_ref)/Lsup(v_ref)
1007     error_t=abs(mod.getTime()-t_ref)/abs(t_ref)
1008     if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v
1009 gross 2850 self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
1010     self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )
1011     self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )
1012     self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )
1013 gross 2432 def tearDown(self):
1014     del self.dom
1015    
1016     def test_D2_Fixed_MuNone_LateStart(self):
1017     self.dom = Rectangle(self.NE,self.NE,order=2)
1018     self.mu=None
1019     self.latestart=True
1020     self.runIt()
1021     def test_D2_Fixed_Mu_LateStart(self):
1022     self.dom = Rectangle(self.NE,self.NE,order=2)
1023     self.mu=555.
1024     self.latestart=True
1025     self.runIt()
1026     def test_D2_Fixed_MuNone(self):
1027     self.dom = Rectangle(self.NE,self.NE,order=2)
1028     self.mu=None
1029     self.latestart=False
1030     self.runIt()
1031     def test_D2_Fixed_Mu(self):
1032     self.dom = Rectangle(self.NE,self.NE,order=2)
1033     self.mu=555.
1034     self.latestart=False
1035     self.runIt()
1036     def test_D2_Free_MuNone_LateStart(self):
1037     self.dom = Rectangle(self.NE,self.NE,order=2)
1038     self.mu=None
1039     self.latestart=True
1040     self.runIt(free=0)
1041     def test_D2_Free_Mu_LateStart(self):
1042     self.dom = Rectangle(self.NE,self.NE,order=2)
1043     self.mu=555.
1044     self.latestart=True
1045     self.runIt(free=0)
1046     def test_D2_Free_MuNone(self):
1047     self.dom = Rectangle(self.NE,self.NE,order=2)
1048     self.mu=None
1049     self.latestart=False
1050     self.runIt(free=0)
1051     def test_D2_Free_Mu(self):
1052     self.dom = Rectangle(self.NE,self.NE,order=2)
1053     self.mu=555.
1054     self.latestart=False
1055     self.runIt(free=0)
1056    
1057     def test_D3_Fixed_MuNone_LateStart(self):
1058     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1059     self.mu=None
1060     self.latestart=True
1061     self.runIt()
1062     def test_D3_Fixed_Mu_LateStart(self):
1063     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1064     self.mu=555.
1065     self.latestart=True
1066     self.runIt()
1067     def test_D3_Fixed_MuNone(self):
1068     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1069     self.mu=None
1070     self.latestart=False
1071     self.runIt()
1072     def test_D3_Fixed_Mu(self):
1073     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1074     self.mu=555.
1075     self.latestart=False
1076     self.runIt()
1077     def test_D3_Free_MuNone_LateStart(self):
1078     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1079     self.mu=None
1080     self.latestart=True
1081     self.runIt(free=0)
1082     def test_D3_Free_Mu_LateStart(self):
1083     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1084     self.mu=555.
1085     self.latestart=True
1086     self.runIt(free=0)
1087     def test_D3_Free_MuNone(self):
1088     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1089     self.mu=None
1090     self.latestart=False
1091     self.runIt(free=0)
1092     def test_D3_Free_Mu(self):
1093     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1094     self.mu=555.
1095     self.latestart=False
1096     self.runIt(free=0)
1097    
1098 gross 2647
1099     class Test_FaultSystem(unittest.TestCase):
1100     EPS=1.e-8
1101     NE=10
1102 gross 2663 def test_Fault_MaxValue(self):
1103 gross 2647 dom=Rectangle(2*self.NE,2*self.NE)
1104     x=dom.getX()
1105     f=FaultSystem(dim=2)
1106 gross 2663 f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
1107     f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
1108 gross 2647
1109 gross 2676 u=x[0]*(1.-x[0])*(1-x[1])
1110     t, loc=f.getMaxValue(u)
1111     p=f.getParametrization(x,t)[0]
1112     m, l=loc(u), loc(p)
1113 gross 2647 self.failUnless( m == 0.25, "wrong max value")
1114     self.failUnless( t == 1, "wrong max tag")
1115     self.failUnless( l == 0., "wrong max location")
1116 gross 2676
1117     u=x[1]*(1.-x[1])*(1-x[0])*x[0]
1118     t, loc=f.getMaxValue(u)
1119     p=f.getParametrization(x,t)[0]
1120     m, l=loc(u), loc(p)
1121 gross 2647 self.failUnless( m == 0.0625, "wrong max value")
1122     self.failUnless( t == 2, "wrong max tag")
1123     self.failUnless( l == 0.5, "wrong max location")
1124 gross 2676
1125     u=x[0]*(1.-x[0])*x[1]
1126     t, loc=f.getMaxValue(u)
1127     p=f.getParametrization(x,t)[0]
1128     m, l=loc(u), loc(p)
1129 gross 2647 self.failUnless( m == 0.25, "wrong max value")
1130     self.failUnless( t == 2, "wrong max tag")
1131     self.failUnless( l == 1.0, "wrong max location")
1132 gross 2676
1133     u=x[1]*(1.-x[1])*x[0]
1134     t, loc=f.getMaxValue(u)
1135     p=f.getParametrization(x,t)[0]
1136     m, l=loc(u), loc(p)
1137 gross 2647 self.failUnless( m == 0.25, "wrong max value")
1138     self.failUnless( t == 2, "wrong max tag")
1139     self.failUnless( l == 0., "wrong max location")
1140 gross 2676
1141     u=x[1]*(1.-x[1])*(1.-x[0])
1142     t, loc=f.getMaxValue(u)
1143     p=f.getParametrization(x,t)[0]
1144     m, l=loc(u), loc(p)
1145 gross 2647 self.failUnless( m == 0.25, "wrong max value")
1146     self.failUnless( t == 1, "wrong max tag")
1147     self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong max location")
1148 gross 2675 def test_Fault_MinValue(self):
1149     dom=Rectangle(2*self.NE,2*self.NE)
1150     x=dom.getX()
1151     f=FaultSystem(dim=2)
1152     f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
1153     f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
1154 gross 2647
1155 gross 2676 u=-x[0]*(1.-x[0])*(1-x[1])
1156     t, loc=f.getMinValue(u)
1157     p=f.getParametrization(x,t)[0]
1158     m, l=loc(u), loc(p)
1159 gross 2675 self.failUnless( m == -0.25, "wrong min value")
1160     self.failUnless( t == 1, "wrong min tag")
1161     self.failUnless( l == 0., "wrong min location")
1162 gross 2676 u=-x[1]*(1.-x[1])*(1-x[0])*x[0]
1163     t, loc=f.getMinValue(u)
1164     p=f.getParametrization(x,t)[0]
1165     m, l=loc(u), loc(p)
1166 gross 2675 self.failUnless( m == -0.0625, "wrong min value")
1167     self.failUnless( t == 2, "wrong min tag")
1168     self.failUnless( l == 0.5, "wrong min location")
1169 gross 2676 u=-x[0]*(1.-x[0])*x[1]
1170     t, loc=f.getMinValue(u)
1171     p=f.getParametrization(x,t)[0]
1172     m, l=loc(u), loc(p)
1173 gross 2675 self.failUnless( m == -0.25, "wrong min value")
1174     self.failUnless( t == 2, "wrong min tag")
1175     self.failUnless( l == 1.0, "wrong min location")
1176 gross 2676 u=-x[1]*(1.-x[1])*x[0]
1177     t, loc=f.getMinValue(u)
1178     p=f.getParametrization(x,t)[0]
1179     m, l=loc(u), loc(p)
1180 gross 2675 self.failUnless( m == -0.25, "wrong min value")
1181     self.failUnless( t == 2, "wrong min tag")
1182     self.failUnless( l == 0., "wrong min location")
1183 gross 2676 u=-x[1]*(1.-x[1])*(1.-x[0])
1184     t, loc=f.getMinValue(u)
1185     p=f.getParametrization(x,t)[0]
1186     m, l=loc(u), loc(p)
1187 gross 2675 self.failUnless( m == -0.25, "wrong min value")
1188     self.failUnless( t == 1, "wrong min tag")
1189     self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong min location")
1190    
1191 gross 2647
1192 gross 2663 def test_Fault2D(self):
1193 gross 2647 f=FaultSystem(dim=2)
1194     top1=[ [1.,0.], [1.,1.], [0.,1.] ]
1195 gross 2663 self.failUnlessRaises(ValueError,f.addFault,V0=[1.,0],strikes=[pi/2, pi/2],ls=[1.,1.],tag=1,dips=top1)
1196     f.addFault(V0=[1.,0],strikes=[pi/2, pi],ls=[1.,1.],tag=1)
1197 gross 2647 self.failUnless(f.getDim() == 2, "wrong dimension")
1198     self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
1199 gross 2663 self.failUnless( 2. == f.getTotalLength(1), "length wrong")
1200     self.failUnless( 0. == f.getMediumDepth(1), "depth wrong")
1201 gross 2647 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 range")
1202     self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 range")
1203     self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsets")
1204 gross 2663 segs=f.getTopPolyline(1)
1205 gross 2647 self.failUnless( len(segs) == 3, "wrong number of segments")
1206     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1207     self.failUnless( segs[0].size == 2, "seg 0 has wrong size.")
1208     self.failUnless( numpy.linalg.norm(segs[0]-[1.,0.]) < self.EPS, "wrong vertex. 0 ")
1209     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1210     self.failUnless( segs[1].size == 2, "seg 1 has wrong size.")
1211     self.failUnless( numpy.linalg.norm(segs[1]-[1.,1.]) < self.EPS, "wrong vertex. 1 ")
1212     self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1213     self.failUnless( segs[2].size == 2, "seg 2 has wrong size.")
1214     self.failUnless( numpy.linalg.norm(segs[2]-[0.,1.]) < self.EPS, "wrong vertex. 2 ")
1215     c=f.getCenterOnSurface()
1216     self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1217     self.failUnless( c.size == 2, "center size is wrong")
1218     self.failUnless( numpy.linalg.norm(c-[2./3.,2./3.]) < self.EPS, "center has wrong coordinates.")
1219     o=f.getOrientationOnSurface()/pi*180.
1220 gross 2663 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1221 gross 2647
1222     top2=[ [10.,0.], [0.,10.] ]
1223 gross 2663 f.addFault(V0=[10.,0],strikes=[3.*pi/4],ls=[14.142135623730951], tag=2, w0_offsets=[0,20], w1_max=20)
1224 gross 2647 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1225 gross 2663 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length")
1226     self.failUnless( 0. == f.getMediumDepth(2), "depth wrong")
1227 gross 2647 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range")
1228     self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range")
1229     self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets")
1230 gross 2663 segs=f.getTopPolyline(2)
1231 gross 2647 self.failUnless( len(segs) == 2, "wrong number of segments")
1232     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1233     self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.]) < self.EPS, "wrong vertex. 0 ")
1234     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1235     self.failUnless( numpy.linalg.norm(segs[1]-[0.,10.]) < self.EPS, "wrong vertex. 1 ")
1236     c=f.getCenterOnSurface()
1237     self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1238     self.failUnless( c.size == 2, "center size is wrong")
1239     self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.]) < self.EPS, "center has wrong coordinates.")
1240     o=f.getOrientationOnSurface()/pi*180.
1241 gross 2663 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1242 gross 2647
1243 gross 2654 s,d=f.getSideAndDistance([0.,0.], tag=1)
1244     self.failUnless( s<0, "wrong side.")
1245     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1246     s,d=f.getSideAndDistance([0.,2.], tag=1)
1247     self.failUnless( s>0, "wrong side.")
1248     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1249     s,d=f.getSideAndDistance([1.,2.], tag=1)
1250     self.failUnless( s>0, "wrong side.")
1251     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1252     s,d=f.getSideAndDistance([2.,1.], tag=1)
1253     self.failUnless( s>0, "wrong side.")
1254     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1255     s,d=f.getSideAndDistance([2.,0.], tag=1)
1256     self.failUnless( s>0, "wrong side.")
1257     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1258     s,d=f.getSideAndDistance([0.,-1.], tag=1)
1259     self.failUnless( s<0, "wrong side.")
1260     self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
1261     s,d=f.getSideAndDistance([-1.,0], tag=1)
1262     self.failUnless( s<0, "wrong side.")
1263     self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
1264    
1265    
1266 gross 2647 f.transform(rot=-pi/2., shift=[-1.,-1.])
1267     self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
1268 gross 2663 self.failUnless( 2. == f.getTotalLength(1), "length after transformation wrong")
1269     self.failUnless( 0. == f.getMediumDepth(1), "depth after transformation wrong")
1270 gross 2647 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 after transformation range")
1271     self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
1272     self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsetsafter transformation ")
1273 gross 2663 segs=f.getTopPolyline(1)
1274 gross 2647 self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
1275     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1276     self.failUnless( segs[0].size == 2, "seg 0 has wrong size after transformation.")
1277     self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1278     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1279     self.failUnless( segs[1].size == 2, "seg 1 has wrong size after transformation.")
1280     self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1281     self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1282     self.failUnless( segs[2].size == 2, "seg 2 has wrong size after transformation.")
1283     self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.]) < self.EPS, "wrong vertex after transformation. 2 ")
1284 gross 2663 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1285     self.failUnless( 0. == f.getMediumDepth(2), "depth wrong after transformation")
1286 gross 2647 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range after transformation")
1287     self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1288     self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1289 gross 2663 segs=f.getTopPolyline(2)
1290 gross 2647 self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
1291     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1292     self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9]) < self.EPS, "wrong vertex. 0 after transformation")
1293     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1294     self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.]) < self.EPS, "wrong vertex. 1 after transformation")
1295    
1296     c=f.getCenterOnSurface()
1297     self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1298     self.failUnless( c.size == 2, "center size is wrong")
1299     self.failUnless( numpy.linalg.norm(c-[7./5.,-7./5.]) < self.EPS, "center has wrong coordinates.")
1300     o=f.getOrientationOnSurface()/pi*180.
1301 gross 2663 self.failUnless( abs(o-45.) < self.EPS, "wrong orientation.")
1302 gross 2647
1303     p=f.getParametrization([-1.,0.],1)
1304     self.failUnless(p[1]==1., "wrong value.")
1305     self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1306     p=f.getParametrization([-0.5,0.],1)
1307     self.failUnless(p[1]==1., "wrong value.")
1308     self.failUnless(abs(p[0]-0.5)<self.EPS* 0.5, "wrong value.")
1309     p=f.getParametrization([0.,0.],1)
1310     self.failUnless(p[1]==1., "wrong value.")
1311     self.failUnless(abs(p[0]-1.)<self.EPS, "wrong value.")
1312     p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-8)
1313     self.failUnless(p[1]==0., "wrong value.")
1314     p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-6)
1315     self.failUnless(p[1]==1., "wrong value.")
1316     self.failUnless(abs(p[0]-1.0000001)<self.EPS, "wrong value.")
1317     p=f.getParametrization([0.,0.5],1)
1318     self.failUnless(p[1]==1., "wrong value.")
1319     self.failUnless(abs(p[0]-1.5)<self.EPS, "wrong value.")
1320     p=f.getParametrization([0,1.],1)
1321     self.failUnless(p[1]==1., "wrong value.")
1322     self.failUnless(abs(p[0]-2.)<self.EPS, "wrong value.")
1323     p=f.getParametrization([1.,1.],1)
1324     self.failUnless(p[1]==0., "wrong value.")
1325     p=f.getParametrization([0,1.11],1)
1326     self.failUnless(p[1]==0., "wrong value.")
1327     p=f.getParametrization([-1,-9.],2)
1328     self.failUnless(p[1]==1., "wrong value.")
1329     self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1330     p=f.getParametrization([9,1],2)
1331     self.failUnless(p[1]==1., "wrong value.")
1332     self.failUnless(abs(p[0]-20.)<self.EPS, "wrong value.")
1333    
1334 gross 2663 def test_Fault3D(self):
1335     f=FaultSystem(dim=3)
1336     self.failUnless(f.getDim() == 3, "wrong dimension")
1337 gross 2647
1338 gross 2663 top1=[ [0.,0.,0.], [1., 0., 0.] ]
1339     f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=20.,tag=1)
1340     self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
1341     self.failUnless( 1. == f.getTotalLength(1), "length wrong")
1342     self.failUnless( 20. == f.getMediumDepth(1), "depth wrong")
1343     self.failUnless( (0., 1.) == f.getW0Range(1)," wrong W0 range")
1344     self.failUnless( (-20., 0.) == f.getW1Range(1)," wrong W1 range")
1345     self.failUnless( [0., 1.] == f.getW0Offsets(1)," wrong W0 offsets")
1346     segs=f.getTopPolyline(1)
1347     self.failUnless( len(segs) == 2, "wrong number of segments")
1348     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1349     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1350     self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1351     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1352     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1353     self.failUnless( numpy.linalg.norm(segs[1]-[1.,0.,0]) < self.EPS, "wrong vertex. 1 ")
1354     c=f.getCenterOnSurface()
1355     self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1356     self.failUnless( c.size == 3, "center size is wrong")
1357     self.failUnless( numpy.linalg.norm(c-[0.5,0.,0.]) < self.EPS, "center has wrong coordinates.")
1358     o=f.getOrientationOnSurface()/pi*180.
1359     self.failUnless( abs(o) < self.EPS, "wrong orientation.")
1360     d=f.getDips(1)
1361     self.failUnless( len(d) == 1, "wrong number of dips")
1362     self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1363     sn=f.getSegmentNormals(1)
1364     self.failUnless( len(sn) == 1, "wrong number of normals")
1365     self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1366     self.failUnless( numpy.linalg.norm(sn[0]-[0, -1., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1367     dv=f.getDepthVectors(1)
1368     self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1369     self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1370     self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1371     self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1372     self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1373     b=f.getBottomPolyline(1)
1374     self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1375     self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1376     self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1377     self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1378     self.failUnless( numpy.linalg.norm(b[1]-[1., 0., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1379     ds=f.getDepths(1)
1380     self.failUnless( len(ds) == 2, "wrong number of depth")
1381     self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1382     self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1383 gross 2647
1384 gross 2663 top2=[ [0.,0.,0.], [0., 10., 0.] ]
1385     f.addFault(V0=[0.,0,0],strikes=[pi/2],ls=[10.,], dips=pi/2, depths=20.,tag=2)
1386     self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1387     self.failUnless( 10. == f.getTotalLength(2), "length wrong")
1388     self.failUnless( 20. == f.getMediumDepth(2), "depth wrong")
1389     self.failUnless( (0., 10.) == f.getW0Range(2)," wrong W0 range")
1390     self.failUnless( (-20., 0.) == f.getW1Range(2)," wrong W1 range")
1391     self.failUnless( [0., 10.] == f.getW0Offsets(2)," wrong W0 offsets")
1392     segs=f.getTopPolyline(2)
1393     self.failUnless( len(segs) == 2, "wrong number of segments")
1394     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1395     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1396     self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1397     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1398     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1399     self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1400     d=f.getDips(2)
1401     self.failUnless( len(d) == 1, "wrong number of dips")
1402     self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1403     sn=f.getSegmentNormals(2)
1404     self.failUnless( len(sn) == 1, "wrong number of normals")
1405     self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1406     self.failUnless( numpy.linalg.norm(sn[0]-[1, 0., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1407     dv=f.getDepthVectors(2)
1408     self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1409     self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1410     self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1411     self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1412     self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1413     b=f.getBottomPolyline(2)
1414     self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1415     self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1416     self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1417     self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1418     self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1419     ds=f.getDepths(2)
1420     self.failUnless( len(ds) == 2, "wrong number of depth")
1421     self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1422     self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1423 gross 2647
1424 gross 2663 top2=[ [10.,0.,0.], [0., 10., 0.] ]
1425     f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/2, depths=30.,tag=2)
1426     self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1427     self.failUnless( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1428     self.failUnless( 30. == f.getMediumDepth(2), "depth wrong")
1429     self.failUnless( (-30., 0.) == f.getW1Range(2)," wrong W1 range")
1430     segs=f.getTopPolyline(2)
1431     self.failUnless( len(segs) == 2, "wrong number of segments")
1432     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1433     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1434     self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1435     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1436     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1437     self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1438     d=f.getDips(2)
1439     self.failUnless( len(d) == 1, "wrong number of dips")
1440     self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1441     sn=f.getSegmentNormals(2)
1442     self.failUnless( len(sn) == 1, "wrong number of normals")
1443     self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1444     self.failUnless( numpy.linalg.norm(sn[0]-[0.70710678118654746, 0.70710678118654746, 0.]) < self.EPS, "wrong bottom vertex 1 ")
1445     dv=f.getDepthVectors(2)
1446     self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1447     self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1448     self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -30.]) < self.EPS, "wrong depth vector 0 ")
1449     self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1450     self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -30.]) < self.EPS, "wrong depth vector 1 ")
1451     b=f.getBottomPolyline(2)
1452     self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1453     self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1454     self.failUnless( numpy.linalg.norm(b[0]-[10., 0., -30.]) < self.EPS, "wrong bottom vertex 0 ")
1455     self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1456     self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -30.]) < self.EPS, "wrong bottom vertex 1 ")
1457     ds=f.getDepths(2)
1458     self.failUnless( len(ds) == 2, "wrong number of depth")
1459     self.failUnless( abs(ds[0]-30.) < self.EPS, "wrong depth at vertex 0 ")
1460     self.failUnless( abs(ds[1]-30.) < self.EPS, "wrong depth at vertex 1 ")
1461    
1462     top2=[ [10.,0.,0.], [0., 10., 0.] ]
1463     f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/4, depths=50.,tag=2)
1464     self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1465     self.failUnless( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1466     self.failUnless( 50. == f.getMediumDepth(2), "depth wrong")
1467     self.failUnless( (-50., 0.) == f.getW1Range(2)," wrong W1 range")
1468     segs=f.getTopPolyline(2)
1469     self.failUnless( len(segs) == 2, "wrong number of segments")
1470     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1471     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1472     self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1473     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1474     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1475     self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1476     d=f.getDips(2)
1477     self.failUnless( len(d) == 1, "wrong number of dips")
1478     self.failUnless( abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1479     sn=f.getSegmentNormals(2)
1480     self.failUnless( len(sn) == 1, "wrong number of normals")
1481     self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1482     self.failUnless( numpy.linalg.norm(sn[0]-[0.5,0.5,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1483     dv=f.getDepthVectors(2)
1484     self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1485     self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1486     self.failUnless( numpy.linalg.norm(dv[0]-[25., 25., -35.355339059327378]) < self.EPS, "wrong depth vector 0 ")
1487     self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1488     self.failUnless( numpy.linalg.norm(dv[1]-[25.,25., -35.355339059327378]) < self.EPS, "wrong depth vector 1 ")
1489     b=f.getBottomPolyline(2)
1490     self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1491     self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1492     self.failUnless( numpy.linalg.norm(b[0]-[35., 25., -35.355339059327378]) < self.EPS, "wrong bottom vertex 0 ")
1493     self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1494     self.failUnless( numpy.linalg.norm(b[1]-[25, 35., -35.355339059327378]) < self.EPS, "wrong bottom vertex 1 ")
1495     ds=f.getDepths(2)
1496     self.failUnless( len(ds) == 2, "wrong number of depth")
1497     self.failUnless( abs(ds[0]-50.) < self.EPS, "wrong depth at vertex 0 ")
1498     self.failUnless( abs(ds[1]-50.) < self.EPS, "wrong depth at vertex 1 ")
1499    
1500     top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
1501     f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=1, dips=pi/4, depths=20.)
1502     self.failUnless( 20. == f.getTotalLength(1), "length wrong")
1503     self.failUnless( 20. == f.getMediumDepth(1), "depth wrong")
1504     segs=f.getTopPolyline(1)
1505     self.failUnless( len(segs) == 3, "wrong number of segments")
1506     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1507     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1508     self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1509     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1510     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1511     self.failUnless( numpy.linalg.norm(segs[1]-[10.,10.,0.]) < self.EPS, "wrong vertex. 1 ")
1512     self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1513     self.failUnless