/[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 2793 - (hide annotations)
Tue Dec 1 06:10:10 2009 UTC (9 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 86160 byte(s)
some new util functions
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 2793 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     TOL=1.e-6
918 gross 2793 VERBOSE=False or True
919 gross 2432 A=1.
920     P_max=100
921 gross 2793 NE=2*getMPISizeWorld()*2
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 2793 mod.setEtaTolerance(self.TOL*1.e-3)
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 2793 mod.update(dt, eta_iter_max=10, iter_max=50, verbose=self.VERBOSE, usePCG=True, max_correction_steps=30)
996 gross 2432 self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)
997     t+=dt
998     N_t+=1
999    
1000     def check(self,N_t,mod,t_ref,v_ref, s_ref,p_ref):
1001     p=mod.getPressure()
1002     p-=integrate(p)/vol(self.dom)
1003     error_p=Lsup(mod.getPressure()-p_ref)/Lsup(p_ref)
1004     error_s=Lsup(mod.getDeviatoricStress()-s_ref)/Lsup(s_ref)
1005     error_v=Lsup(mod.getVelocity()-v_ref)/Lsup(v_ref)
1006     error_t=abs(mod.getTime()-t_ref)/abs(t_ref)
1007     if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v
1008     self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
1009     self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )
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     def tearDown(self):
1013     del self.dom
1014    
1015     def test_D2_Fixed_MuNone_LateStart(self):
1016     self.dom = Rectangle(self.NE,self.NE,order=2)
1017     self.mu=None
1018     self.latestart=True
1019     self.runIt()
1020     def test_D2_Fixed_Mu_LateStart(self):
1021     self.dom = Rectangle(self.NE,self.NE,order=2)
1022     self.mu=555.
1023     self.latestart=True
1024     self.runIt()
1025     def test_D2_Fixed_MuNone(self):
1026     self.dom = Rectangle(self.NE,self.NE,order=2)
1027     self.mu=None
1028     self.latestart=False
1029     self.runIt()
1030     def test_D2_Fixed_Mu(self):
1031     self.dom = Rectangle(self.NE,self.NE,order=2)
1032     self.mu=555.
1033     self.latestart=False
1034     self.runIt()
1035     def test_D2_Free_MuNone_LateStart(self):
1036     self.dom = Rectangle(self.NE,self.NE,order=2)
1037     self.mu=None
1038     self.latestart=True
1039     self.runIt(free=0)
1040     def test_D2_Free_Mu_LateStart(self):
1041     self.dom = Rectangle(self.NE,self.NE,order=2)
1042     self.mu=555.
1043     self.latestart=True
1044     self.runIt(free=0)
1045     def test_D2_Free_MuNone(self):
1046     self.dom = Rectangle(self.NE,self.NE,order=2)
1047     self.mu=None
1048     self.latestart=False
1049     self.runIt(free=0)
1050     def test_D2_Free_Mu(self):
1051     self.dom = Rectangle(self.NE,self.NE,order=2)
1052     self.mu=555.
1053     self.latestart=False
1054     self.runIt(free=0)
1055    
1056     def test_D3_Fixed_MuNone_LateStart(self):
1057     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1058     self.mu=None
1059     self.latestart=True
1060     self.runIt()
1061     def test_D3_Fixed_Mu_LateStart(self):
1062     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1063     self.mu=555.
1064     self.latestart=True
1065     self.runIt()
1066     def test_D3_Fixed_MuNone(self):
1067     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1068     self.mu=None
1069     self.latestart=False
1070     self.runIt()
1071     def test_D3_Fixed_Mu(self):
1072     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1073     self.mu=555.
1074     self.latestart=False
1075     self.runIt()
1076     def test_D3_Free_MuNone_LateStart(self):
1077     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1078     self.mu=None
1079     self.latestart=True
1080     self.runIt(free=0)
1081     def test_D3_Free_Mu_LateStart(self):
1082     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1083     self.mu=555.
1084     self.latestart=True
1085     self.runIt(free=0)
1086     def test_D3_Free_MuNone(self):
1087     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1088     self.mu=None
1089     self.latestart=False
1090     self.runIt(free=0)
1091     def test_D3_Free_Mu(self):
1092     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1093     self.mu=555.
1094     self.latestart=False
1095     self.runIt(free=0)
1096    
1097 gross 2647
1098     class Test_FaultSystem(unittest.TestCase):
1099     EPS=1.e-8
1100     NE=10
1101 gross 2663 def test_Fault_MaxValue(self):
1102 gross 2647 dom=Rectangle(2*self.NE,2*self.NE)
1103     x=dom.getX()
1104     f=FaultSystem(dim=2)
1105 gross 2663 f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
1106     f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
1107 gross 2647
1108 gross 2676 u=x[0]*(1.-x[0])*(1-x[1])
1109     t, loc=f.getMaxValue(u)
1110     p=f.getParametrization(x,t)[0]
1111     m, l=loc(u), loc(p)
1112 gross 2647 self.failUnless( m == 0.25, "wrong max value")
1113     self.failUnless( t == 1, "wrong max tag")
1114     self.failUnless( l == 0., "wrong max location")
1115 gross 2676
1116     u=x[1]*(1.-x[1])*(1-x[0])*x[0]
1117     t, loc=f.getMaxValue(u)
1118     p=f.getParametrization(x,t)[0]
1119     m, l=loc(u), loc(p)
1120 gross 2647 self.failUnless( m == 0.0625, "wrong max value")
1121     self.failUnless( t == 2, "wrong max tag")
1122     self.failUnless( l == 0.5, "wrong max location")
1123 gross 2676
1124     u=x[0]*(1.-x[0])*x[1]
1125     t, loc=f.getMaxValue(u)
1126     p=f.getParametrization(x,t)[0]
1127     m, l=loc(u), loc(p)
1128 gross 2647 self.failUnless( m == 0.25, "wrong max value")
1129     self.failUnless( t == 2, "wrong max tag")
1130     self.failUnless( l == 1.0, "wrong max location")
1131 gross 2676
1132     u=x[1]*(1.-x[1])*x[0]
1133     t, loc=f.getMaxValue(u)
1134     p=f.getParametrization(x,t)[0]
1135     m, l=loc(u), loc(p)
1136 gross 2647 self.failUnless( m == 0.25, "wrong max value")
1137     self.failUnless( t == 2, "wrong max tag")
1138     self.failUnless( l == 0., "wrong max location")
1139 gross 2676
1140     u=x[1]*(1.-x[1])*(1.-x[0])
1141     t, loc=f.getMaxValue(u)
1142     p=f.getParametrization(x,t)[0]
1143     m, l=loc(u), loc(p)
1144 gross 2647 self.failUnless( m == 0.25, "wrong max value")
1145     self.failUnless( t == 1, "wrong max tag")
1146     self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong max location")
1147 gross 2675 def test_Fault_MinValue(self):
1148     dom=Rectangle(2*self.NE,2*self.NE)
1149     x=dom.getX()
1150     f=FaultSystem(dim=2)
1151     f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
1152     f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
1153 gross 2647
1154 gross 2676 u=-x[0]*(1.-x[0])*(1-x[1])
1155     t, loc=f.getMinValue(u)
1156     p=f.getParametrization(x,t)[0]
1157     m, l=loc(u), loc(p)
1158 gross 2675 self.failUnless( m == -0.25, "wrong min value")
1159     self.failUnless( t == 1, "wrong min tag")
1160     self.failUnless( l == 0., "wrong min location")
1161 gross 2676 u=-x[1]*(1.-x[1])*(1-x[0])*x[0]
1162     t, loc=f.getMinValue(u)
1163     p=f.getParametrization(x,t)[0]
1164     m, l=loc(u), loc(p)
1165 gross 2675 self.failUnless( m == -0.0625, "wrong min value")
1166     self.failUnless( t == 2, "wrong min tag")
1167     self.failUnless( l == 0.5, "wrong min location")
1168 gross 2676 u=-x[0]*(1.-x[0])*x[1]
1169     t, loc=f.getMinValue(u)
1170     p=f.getParametrization(x,t)[0]
1171     m, l=loc(u), loc(p)
1172 gross 2675 self.failUnless( m == -0.25, "wrong min value")
1173     self.failUnless( t == 2, "wrong min tag")
1174     self.failUnless( l == 1.0, "wrong min location")
1175 gross 2676 u=-x[1]*(1.-x[1])*x[0]
1176     t, loc=f.getMinValue(u)
1177     p=f.getParametrization(x,t)[0]
1178     m, l=loc(u), loc(p)
1179 gross 2675 self.failUnless( m == -0.25, "wrong min value")
1180     self.failUnless( t == 2, "wrong min tag")
1181     self.failUnless( l == 0., "wrong min location")
1182 gross 2676 u=-x[1]*(1.-x[1])*(1.-x[0])
1183     t, loc=f.getMinValue(u)
1184     p=f.getParametrization(x,t)[0]
1185     m, l=loc(u), loc(p)
1186 gross 2675 self.failUnless( m == -0.25, "wrong min value")
1187     self.failUnless( t == 1, "wrong min tag")
1188     self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong min location")
1189    
1190 gross 2647
1191 gross 2663 def test_Fault2D(self):
1192 gross 2647 f=FaultSystem(dim=2)
1193     top1=[ [1.,0.], [1.,1.], [0.,1.] ]
1194 gross 2663 self.failUnlessRaises(ValueError,f.addFault,V0=[1.,0],strikes=[pi/2, pi/2],ls=[1.,1.],tag=1,dips=top1)
1195     f.addFault(V0=[1.,0],strikes=[pi/2, pi],ls=[1.,1.],tag=1)
1196 gross 2647 self.failUnless(f.getDim() == 2, "wrong dimension")
1197     self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
1198 gross 2663 self.failUnless( 2. == f.getTotalLength(1), "length wrong")
1199     self.failUnless( 0. == f.getMediumDepth(1), "depth wrong")
1200 gross 2647 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 range")
1201     self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 range")
1202     self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsets")
1203 gross 2663 segs=f.getTopPolyline(1)
1204 gross 2647 self.failUnless( len(segs) == 3, "wrong number of segments")
1205     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1206     self.failUnless( segs[0].size == 2, "seg 0 has wrong size.")
1207     self.failUnless( numpy.linalg.norm(segs[0]-[1.,0.]) < self.EPS, "wrong vertex. 0 ")
1208     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1209     self.failUnless( segs[1].size == 2, "seg 1 has wrong size.")
1210     self.failUnless( numpy.linalg.norm(segs[1]-[1.,1.]) < self.EPS, "wrong vertex. 1 ")
1211     self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1212     self.failUnless( segs[2].size == 2, "seg 2 has wrong size.")
1213     self.failUnless( numpy.linalg.norm(segs[2]-[0.,1.]) < self.EPS, "wrong vertex. 2 ")
1214     c=f.getCenterOnSurface()
1215     self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1216     self.failUnless( c.size == 2, "center size is wrong")
1217     self.failUnless( numpy.linalg.norm(c-[2./3.,2./3.]) < self.EPS, "center has wrong coordinates.")
1218     o=f.getOrientationOnSurface()/pi*180.
1219 gross 2663 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1220 gross 2647
1221     top2=[ [10.,0.], [0.,10.] ]
1222 gross 2663 f.addFault(V0=[10.,0],strikes=[3.*pi/4],ls=[14.142135623730951], tag=2, w0_offsets=[0,20], w1_max=20)
1223 gross 2647 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1224 gross 2663 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length")
1225     self.failUnless( 0. == f.getMediumDepth(2), "depth wrong")
1226 gross 2647 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range")
1227     self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range")
1228     self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets")
1229 gross 2663 segs=f.getTopPolyline(2)
1230 gross 2647 self.failUnless( len(segs) == 2, "wrong number of segments")
1231     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1232     self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.]) < self.EPS, "wrong vertex. 0 ")
1233     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1234     self.failUnless( numpy.linalg.norm(segs[1]-[0.,10.]) < self.EPS, "wrong vertex. 1 ")
1235     c=f.getCenterOnSurface()
1236     self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1237     self.failUnless( c.size == 2, "center size is wrong")
1238     self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.]) < self.EPS, "center has wrong coordinates.")
1239     o=f.getOrientationOnSurface()/pi*180.
1240 gross 2663 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1241 gross 2647
1242 gross 2654 s,d=f.getSideAndDistance([0.,0.], tag=1)
1243     self.failUnless( s<0, "wrong side.")
1244     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1245     s,d=f.getSideAndDistance([0.,2.], tag=1)
1246     self.failUnless( s>0, "wrong side.")
1247     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1248     s,d=f.getSideAndDistance([1.,2.], tag=1)
1249     self.failUnless( s>0, "wrong side.")
1250     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1251     s,d=f.getSideAndDistance([2.,1.], tag=1)
1252     self.failUnless( s>0, "wrong side.")
1253     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1254     s,d=f.getSideAndDistance([2.,0.], tag=1)
1255     self.failUnless( s>0, "wrong side.")
1256     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1257     s,d=f.getSideAndDistance([0.,-1.], tag=1)
1258     self.failUnless( s<0, "wrong side.")
1259     self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
1260     s,d=f.getSideAndDistance([-1.,0], tag=1)
1261     self.failUnless( s<0, "wrong side.")
1262     self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
1263    
1264    
1265 gross 2647 f.transform(rot=-pi/2., shift=[-1.,-1.])
1266     self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
1267 gross 2663 self.failUnless( 2. == f.getTotalLength(1), "length after transformation wrong")
1268     self.failUnless( 0. == f.getMediumDepth(1), "depth after transformation wrong")
1269 gross 2647 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 after transformation range")
1270     self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
1271     self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsetsafter transformation ")
1272 gross 2663 segs=f.getTopPolyline(1)
1273 gross 2647 self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
1274     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1275     self.failUnless( segs[0].size == 2, "seg 0 has wrong size after transformation.")
1276     self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1277     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1278     self.failUnless( segs[1].size == 2, "seg 1 has wrong size after transformation.")
1279     self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1280     self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1281     self.failUnless( segs[2].size == 2, "seg 2 has wrong size after transformation.")
1282     self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.]) < self.EPS, "wrong vertex after transformation. 2 ")
1283 gross 2663 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1284     self.failUnless( 0. == f.getMediumDepth(2), "depth wrong after transformation")
1285 gross 2647 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range after transformation")
1286     self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1287     self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1288 gross 2663 segs=f.getTopPolyline(2)
1289 gross 2647 self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
1290     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1291     self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9]) < self.EPS, "wrong vertex. 0 after transformation")
1292     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1293     self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.]) < self.EPS, "wrong vertex. 1 after transformation")
1294    
1295     c=f.getCenterOnSurface()
1296     self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1297     self.failUnless( c.size == 2, "center size is wrong")
1298     self.failUnless( numpy.linalg.norm(c-[7./5.,-7./5.]) < self.EPS, "center has wrong coordinates.")
1299     o=f.getOrientationOnSurface()/pi*180.
1300 gross 2663 self.failUnless( abs(o-45.) < self.EPS, "wrong orientation.")
1301 gross 2647
1302     p=f.getParametrization([-1.,0.],1)
1303     self.failUnless(p[1]==1., "wrong value.")
1304     self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1305     p=f.getParametrization([-0.5,0.],1)
1306     self.failUnless(p[1]==1., "wrong value.")
1307     self.failUnless(abs(p[0]-0.5)<self.EPS* 0.5, "wrong value.")
1308     p=f.getParametrization([0.,0.],1)
1309     self.failUnless(p[1]==1., "wrong value.")
1310     self.failUnless(abs(p[0]-1.)<self.EPS, "wrong value.")
1311     p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-8)
1312     self.failUnless(p[1]==0., "wrong value.")
1313     p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-6)
1314     self.failUnless(p[1]==1., "wrong value.")
1315     self.failUnless(abs(p[0]-1.0000001)<self.EPS, "wrong value.")
1316     p=f.getParametrization([0.,0.5],1)
1317     self.failUnless(p[1]==1., "wrong value.")
1318     self.failUnless(abs(p[0]-1.5)<self.EPS, "wrong value.")
1319     p=f.getParametrization([0,1.],1)
1320     self.failUnless(p[1]==1., "wrong value.")
1321     self.failUnless(abs(p[0]-2.)<self.EPS, "wrong value.")
1322     p=f.getParametrization([1.,1.],1)
1323     self.failUnless(p[1]==0., "wrong value.")
1324     p=f.getParametrization([0,1.11],1)
1325     self.failUnless(p[1]==0., "wrong value.")
1326     p=f.getParametrization([-1,-9.],2)
1327     self.failUnless(p[1]==1., "wrong value.")
1328     self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1329     p=f.getParametrization([9,1],2)
1330     self.failUnless(p[1]==1., "wrong value.")
1331     self.failUnless(abs(p[0]-20.)<self.EPS, "wrong value.")
1332    
1333 gross 2663 def test_Fault3D(self):
1334     f=FaultSystem(dim=3)
1335     self.failUnless(f.getDim() == 3, "wrong dimension")
1336 gross 2647
1337 gross 2663 top1=[ [0.,0.,0.], [1., 0., 0.] ]
1338     f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=20.,tag=1)
1339     self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
1340     self.failUnless( 1. == f.getTotalLength(1), "length wrong")
1341     self.failUnless( 20. == f.getMediumDepth(1), "depth wrong")
1342     self.failUnless( (0., 1.) == f.getW0Range(1)," wrong W0 range")
1343     self.failUnless( (-20., 0.) == f.getW1Range(1)," wrong W1 range")
1344     self.failUnless( [0., 1.] == f.getW0Offsets(1)," wrong W0 offsets")
1345     segs=f.getTopPolyline(1)
1346     self.failUnless( len(segs) == 2, "wrong number of segments")
1347     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1348     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1349     self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1350     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1351     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1352     self.failUnless( numpy.linalg.norm(segs[1]-[1.,0.,0]) < self.EPS, "wrong vertex. 1 ")
1353     c=f.getCenterOnSurface()
1354     self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1355     self.failUnless( c.size == 3, "center size is wrong")
1356     self.failUnless( numpy.linalg.norm(c-[0.5,0.,0.]) < self.EPS, "center has wrong coordinates.")
1357     o=f.getOrientationOnSurface()/pi*180.
1358     self.failUnless( abs(o) < self.EPS, "wrong orientation.")
1359     d=f.getDips(1)
1360     self.failUnless( len(d) == 1, "wrong number of dips")
1361     self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1362     sn=f.getSegmentNormals(1)
1363     self.failUnless( len(sn) == 1, "wrong number of normals")
1364     self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1365     self.failUnless( numpy.linalg.norm(sn[0]-[0, -1., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1366     dv=f.getDepthVectors(1)
1367     self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1368     self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1369     self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1370     self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1371     self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1372     b=f.getBottomPolyline(1)
1373     self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1374     self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1375     self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1376     self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1377     self.failUnless( numpy.linalg.norm(b[1]-[1., 0., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1378     ds=f.getDepths(1)
1379     self.failUnless( len(ds) == 2, "wrong number of depth")
1380     self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1381     self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1382 gross 2647
1383 gross 2663 top2=[ [0.,0.,0.], [0., 10., 0.] ]
1384     f.addFault(V0=[0.,0,0],strikes=[pi/2],ls=[10.,], dips=pi/2, depths=20.,tag=2)
1385     self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1386     self.failUnless( 10. == f.getTotalLength(2), "length wrong")
1387     self.failUnless( 20. == f.getMediumDepth(2), "depth wrong")
1388     self.failUnless( (0., 10.) == f.getW0Range(2)," wrong W0 range")
1389     self.failUnless( (-20., 0.) == f.getW1Range(2)," wrong W1 range")
1390     self.failUnless( [0., 10.] == f.getW0Offsets(2)," wrong W0 offsets")
1391     segs=f.getTopPolyline(2)
1392     self.failUnless( len(segs) == 2, "wrong number of segments")
1393     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1394     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1395     self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1396     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1397     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1398     self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1399     d=f.getDips(2)
1400     self.failUnless( len(d) == 1, "wrong number of dips")
1401     self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1402     sn=f.getSegmentNormals(2)
1403     self.failUnless( len(sn) == 1, "wrong number of normals")
1404     self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1405     self.failUnless( numpy.linalg.norm(sn[0]-[1, 0., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1406     dv=f.getDepthVectors(2)
1407     self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1408     self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1409     self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1410     self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1411     self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1412     b=f.getBottomPolyline(2)
1413     self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1414     self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1415     self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1416     self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1417     self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1418     ds=f.getDepths(2)
1419     self.failUnless( len(ds) == 2, "wrong number of depth")
1420     self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1421     self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1422 gross 2647
1423 gross 2663 top2=[ [10.,0.,0.], [0., 10., 0.] ]
1424     f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/2, depths=30.,tag=2)
1425     self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1426     self.failUnless( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1427     self.failUnless( 30. == f.getMediumDepth(2), "depth wrong")
1428     self.failUnless( (-30., 0.) == f.getW1Range(2)," wrong W1 range")
1429     segs=f.getTopPolyline(2)
1430     self.failUnless( len(segs) == 2, "wrong number of segments")
1431     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1432     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1433     self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1434     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1435     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1436     self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1437     d=f.getDips(2)
1438     self.failUnless( len(d) == 1, "wrong number of dips")
1439     self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1440     sn=f.getSegmentNormals(2)
1441     self.failUnless( len(sn) == 1, "wrong number of normals")
1442     self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1443     self.failUnless( numpy.linalg.norm(sn[0]-[0.70710678118654746, 0.70710678118654746, 0.]) < self.EPS, "wrong bottom vertex 1 ")
1444     dv=f.getDepthVectors(2)
1445     self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1446     self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1447     self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -30.]) < self.EPS, "wrong depth vector 0 ")
1448     self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1449     self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -30.]) < self.EPS, "wrong depth vector 1 ")
1450     b=f.getBottomPolyline(2)
1451     self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1452     self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1453     self.failUnless( numpy.linalg.norm(b[0]-[10., 0., -30.]) < self.EPS, "wrong bottom vertex 0 ")
1454     self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1455     self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -30.]) < self.EPS, "wrong bottom vertex 1 ")
1456     ds=f.getDepths(2)
1457     self.failUnless( len(ds) == 2, "wrong number of depth")
1458     self.failUnless( abs(ds[0]-30.) < self.EPS, "wrong depth at vertex 0 ")
1459     self.failUnless( abs(ds[1]-30.) < self.EPS, "wrong depth at vertex 1 ")
1460    
1461     top2=[ [10.,0.,0.], [0., 10., 0.] ]
1462     f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/4, depths=50.,tag=2)
1463     self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1464     self.failUnless( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1465     self.failUnless( 50. == f.getMediumDepth(2), "depth wrong")
1466     self.failUnless( (-50., 0.) == f.getW1Range(2)," wrong W1 range")
1467     segs=f.getTopPolyline(2)
1468     self.failUnless( len(segs) == 2, "wrong number of segments")
1469     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1470     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1471     self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1472     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1473     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1474     self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1475     d=f.getDips(2)
1476     self.failUnless( len(d) == 1, "wrong number of dips")
1477     self.failUnless( abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1478     sn=f.getSegmentNormals(2)
1479     self.failUnless( len(sn) == 1, "wrong number of normals")
1480     self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1481     self.failUnless( numpy.linalg.norm(sn[0]-[0.5,0.5,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1482     dv=f.getDepthVectors(2)
1483     self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1484     self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1485     self.failUnless( numpy.linalg.norm(dv[0]-[25., 25., -35.355339059327378]) < self.EPS, "wrong depth vector 0 ")
1486     self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1487     self.failUnless( numpy.linalg.norm(dv[1]-[25.,25., -35.355339059327378]) < self.EPS, "wrong depth vector 1 ")
1488     b=f.getBottomPolyline(2)
1489     self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1490     self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1491     self.failUnless( numpy.linalg.norm(b[0]-[35., 25., -35.355339059327378]) < self.EPS, "wrong bottom vertex 0 ")
1492     self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1493     self.failUnless( numpy.linalg.norm(b[1]-[25, 35., -35.355339059327378]) < self.EPS, "wrong bottom vertex 1 ")
1494     ds=f.getDepths(2)
1495     self.failUnless( len(ds) == 2, "wrong number of depth")
1496     self.failUnless( abs(ds[0]-50.) < self.EPS, "wrong depth at vertex 0 ")
1497     self.failUnless( abs(ds[1]-50.) < self.EPS, "wrong depth at vertex 1 ")
1498    
1499     top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
1500     f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=1, dips=pi/4, depths=20.)
1501     self.failUnless( 20. == f.getTotalLength(1), "length wrong")
1502     self.failUnless( 20. == f.getMediumDepth(1), "depth wrong")
1503     segs=f.getTopPolyline(1)
1504     self.failUnless( len(segs) == 3, "wrong number of segments")
1505     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1506     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1507     self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1508     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1509     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1510     self.failUnless( numpy.linalg.norm(segs[1]-[10.,10.,0.]) < self.EPS, "wrong vertex. 1 ")
1511     self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1512     self.failUnless( segs[2].size