/[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 2960 - (hide annotations)
Tue Mar 2 07:54:11 2010 UTC (9 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 86341 byte(s)
new DarcyFlux solver
1 ksteube 1809
2     ########################################################
3 gross 1414 #
4 jfenwick 2881 # Copyright (c) 2003-2010 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 2881 __copyright__="""Copyright (c) 2003-2010 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 2960 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 2960 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
470     self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
471 gross 2100 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 2960 self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
486     self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*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 2960 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
503     self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*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 2960 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
520     self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*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 2960 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
537     self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*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 2960 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
554     self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*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 2960 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
571     self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*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 2960 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
588     self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*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 2960 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
605     self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*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 2960 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
622     self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*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 2960 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
639     self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*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 2960 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
656     self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
657 artak 1518
658 gross 2100 class Test_Darcy2D(Test_Darcy):
659 gross 2960 TOL=1e-6
660     TEST_TOL=2.e-3
661 gross 2100 WIDTH=1.
662     def setUp(self):
663 gross 2208 NE=40 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
664 gross 2100 self.dom = Rectangle(NE,NE)
665     self.rescaleDomain()
666     def tearDown(self):
667     del self.dom
668     class Test_Darcy3D(Test_Darcy):
669 gross 2960 TOL=1e-6
670 gross 2100 WIDTH=1.
671 gross 2960 TEST_TOL=4.e-3
672 gross 2100 def setUp(self):
673     NE=25 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
674     self.dom = Brick(NE,NE,NE)
675     self.rescaleDomain()
676     def tearDown(self):
677     del self.dom
678 artak 1518
679 gross 2100
680 artak 2234 class Test_Mountains3D(unittest.TestCase):
681     def setUp(self):
682     NE=16
683     self.TOL=1e-4
684     self.domain=Brick(NE,NE,NE,order=1,useFullElementOrder=True)
685     def tearDown(self):
686     del self.domain
687     def test_periodic(self):
688     EPS=0.01
689 gross 2100
690 artak 2234 x=self.domain.getX()
691     v = Vector(0.0, Solution(self.domain))
692     a0=1
693     a1=1
694     n0=2
695     n1=2
696     n2=0.5
697     a2=-(a0*n0+a1*n1)/n2
698     v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])* cos(pi*n2*x[2])
699     v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])* cos(pi*n2*x[2])
700     v[2]=a2*cos(pi*n0*x[0])* cos(pi*n1*x[1])* sin(pi*n2*x[2])
701    
702 gross 2563 mts=Mountains(self.domain,eps=EPS)
703     mts.setVelocity(v)
704     Z=mts.update()
705 artak 2234
706 gross 2564 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[2]-1.)))
707 artak 2234 self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
708    
709     class Test_Mountains2D(unittest.TestCase):
710     def setUp(self):
711     NE=16
712     self.TOL=1e-4
713     self.domain=Rectangle(NE,NE,order=1,useFullElementOrder=True)
714     def tearDown(self):
715     del self.domain
716     def test_periodic(self):
717     EPS=0.01
718    
719     x=self.domain.getX()
720     v = Vector(0.0, Solution(self.domain))
721     a0=1
722     n0=1
723     n1=0.5
724     a1=-(a0*n0)/n1
725     v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])
726     v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])
727    
728     H_t=Scalar(0.0, Solution(self.domain))
729 gross 2563 mts=Mountains(self.domain,eps=EPS)
730     mts.setVelocity(v)
731     Z=mts.update()
732 artak 2234
733 gross 2564 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[1]-1.)))
734 artak 2234 self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
735    
736 gross 2301
737    
738     class Test_Rheologies(unittest.TestCase):
739     """
740     this is the program used to generate the powerlaw tests:
741    
742     TAU_Y=100.
743     N=10
744     M=5
745    
746     def getE(tau):
747     if tau<=TAU_Y:
748     return 1./(0.5+20*sqrt(tau))
749     else:
750     raise ValueError,"out of range."
751     tau=[ i* (TAU_Y/N) for i in range(N+1)] + [TAU_Y for i in range(M)]
752     e=[ tau[i]/getE(tau[i]) for i in range(N+1)]
753     e+= [ (i+1+j)* (max(e)/N) for j in range(M) ]
754    
755     print tau
756     print e
757     """
758     TOL=1.e-8
759     def test_PowerLaw_Init(self):
760     pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
761    
762     self.failUnlessEqual(pl.getNumMaterials(),3,"num materials is wrong")
763     self.failUnlessEqual(pl.validMaterialId(0),True,"material id 0 not found")
764     self.failUnlessEqual(pl.validMaterialId(1),True,"material id 1 not found")
765     self.failUnlessEqual(pl.validMaterialId(2),True,"material id 2 not found")
766     self.failUnlessEqual(pl.validMaterialId(3),False,"material id 3 not found")
767    
768     self.failUnlessEqual(pl.getFriction(),None,"initial friction wrong.")
769     self.failUnlessEqual(pl.getTauY(),None,"initial tau y wrong.")
770     pl.setDruckerPragerLaw(tau_Y=10,friction=3)
771     self.failUnlessEqual(pl.getFriction(),3,"friction wrong.")
772     self.failUnlessEqual(pl.getTauY(),10,"tau y wrong.")
773    
774     self.failUnlessEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
775     pl.setElasticShearModulus(1000)
776     self.failUnlessEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
777    
778     e=pl.getEtaN()
779     self.failUnlessEqual(len(e),3,"initial length of etaN is wrong.")
780     self.failUnlessEqual(e,[None, None, None],"initial etaN are wrong.")
781     self.failUnlessEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
782     self.failUnlessEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
783     self.failUnlessEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
784     self.failUnlessRaises(ValueError, pl.getEtaN, 3)
785    
786     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
787     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
788     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
789    
790     pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
791     self.failUnlessEqual(pl.getPower(),[1,2,3],"powers are wrong.")
792     self.failUnlessEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
793     self.failUnlessEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
794    
795     pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
796     self.failUnlessEqual(pl.getPower(),[4,2,3],"powers are wrong.")
797     self.failUnlessEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
798     self.failUnlessEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
799    
800     self.failUnlessRaises(ValueError,pl.getPower,-1)
801     self.failUnlessEqual(pl.getPower(0),4,"power 0 is wrong.")
802     self.failUnlessEqual(pl.getPower(1),2,"power 1 is wrong.")
803     self.failUnlessEqual(pl.getPower(2),3,"power 2 is wrong.")
804     self.failUnlessRaises(ValueError,pl.getPower,3)
805    
806     self.failUnlessRaises(ValueError,pl.getTauT,-1)
807     self.failUnlessEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
808     self.failUnlessEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
809     self.failUnlessEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
810     self.failUnlessRaises(ValueError,pl.getTauT,3)
811    
812     self.failUnlessRaises(ValueError,pl.getEtaN,-1)
813     self.failUnlessEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
814     self.failUnlessEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
815     self.failUnlessEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
816     self.failUnlessRaises(ValueError,pl.getEtaN,3)
817    
818     def checkResult(self,id,gamma_dot_, eta, tau_ref):
819     self.failUnless(eta>=0,"eta needs to be positive (test %s)"%id)
820     error=abs(gamma_dot_*eta-tau_ref)
821     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))
822    
823     def test_PowerLaw_Linear(self):
824     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]
825     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]
826     pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
827     pl.setDruckerPragerLaw(tau_Y=100.)
828     pl.setPowerLaw(eta_N=2.)
829     pl.setEtaTolerance(self.TOL)
830     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
831    
832     def test_PowerLaw_QuadLarge(self):
833     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]
834 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]
835 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
836     pl.setDruckerPragerLaw(tau_Y=100.)
837     pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2])
838     pl.setEtaTolerance(self.TOL)
839     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
840    
841     def test_PowerLaw_QuadSmall(self):
842     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]
843 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]
844 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
845     pl.setDruckerPragerLaw(tau_Y=100.)
846     pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2])
847     pl.setEtaTolerance(self.TOL)
848     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
849    
850     def test_PowerLaw_CubeLarge(self):
851     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]
852 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]
853 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
854     pl.setDruckerPragerLaw(tau_Y=100.)
855     pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3])
856     pl.setEtaTolerance(self.TOL)
857     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
858    
859     def test_PowerLaw_CubeSmall(self):
860     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]
861 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]
862 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
863     pl.setDruckerPragerLaw(tau_Y=100.)
864     pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3])
865     pl.setEtaTolerance(self.TOL)
866     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
867    
868     def test_PowerLaw_QuadLarge_CubeLarge(self):
869     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]
870 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]
871    
872 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
873     pl.setDruckerPragerLaw(tau_Y=100.)
874     pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
875     pl.setEtaTolerance(self.TOL)
876     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
877    
878     def test_PowerLaw_QuadLarge_CubeSmall(self):
879     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]
880 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]
881    
882 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
883     pl.setDruckerPragerLaw(tau_Y=100.)
884     pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
885     pl.setEtaTolerance(self.TOL)
886     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
887    
888     def test_PowerLaw_QuadSmall_CubeLarge(self):
889     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]
890 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]
891 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
892     pl.setDruckerPragerLaw(tau_Y=100.)
893     pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
894     pl.setEtaTolerance(self.TOL)
895     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
896    
897     def test_PowerLaw_QuadSmall_CubeSmall(self):
898     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]
899 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]
900 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
901     pl.setDruckerPragerLaw(tau_Y=100.)
902     pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
903     pl.setEtaTolerance(self.TOL)
904     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
905    
906     def test_PowerLaw_withShear(self):
907     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]
908     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]
909     pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
910     pl.setDruckerPragerLaw(tau_Y=100.)
911     pl.setPowerLaw(eta_N=2.)
912     pl.setElasticShearModulus(3.)
913     dt=1./3.
914     pl.setEtaTolerance(self.TOL)
915     self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
916     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
917    
918 gross 2432 class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase):
919 gross 2850 TOL=1.e-5
920     VERBOSE=False # or True
921 gross 2432 A=1.
922     P_max=100
923 gross 2794 NE=2*getMPISizeWorld()
924 gross 2432 tau_Y=10.
925     N_dt=10
926    
927     # material parameter:
928     tau_1=5.
929     tau_2=5.
930     eta_0=100.
931     eta_1=50.
932     eta_2=400.
933     N_1=2.
934     N_2=3.
935     def getReference(self, t):
936    
937     B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
938     x=self.dom.getX()
939    
940     s_00=min(self.A*t,B)
941     tau=sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)*abs(s_00)
942 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.)
943 gross 2432
944     alpha=0.5*inv_eta*s_00
945     if s_00 <= B and self.mu !=None: alpha+=1./(2*self.mu)*self.A
946     u_ref=x*alpha
947     u_ref[self.dom.getDim()-1]=(1.-x[self.dom.getDim()-1])*alpha*(self.dom.getDim()-1)
948     sigma_ref=kronecker(self.dom)*s_00
949     sigma_ref[self.dom.getDim()-1,self.dom.getDim()-1]=-s_00*(self.dom.getDim()-1)
950    
951     p_ref=self.P_max
952     for d in range(self.dom.getDim()): p_ref=p_ref*x[d]
953     p_ref-=integrate(p_ref)/vol(self.dom)
954     return u_ref, sigma_ref, p_ref
955    
956     def runIt(self, free=None):
957     x=self.dom.getX()
958     B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
959     dt=B/int(self.N_dt/2)
960     if self.VERBOSE: print "dt =",dt
961     if self.latestart:
962     t=dt
963     else:
964     t=0
965     v,s,p=self.getReference(t)
966    
967     mod=IncompressibleIsotropicFlowCartesian(self.dom, stress=s, v=v, p=p, t=t, numMaterials=3, verbose=self.VERBOSE)
968     mod.setDruckerPragerLaw(tau_Y=self.tau_Y,friction=None)
969     mod.setElasticShearModulus(self.mu)
970     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])
971     mod.setTolerance(self.TOL)
972 gross 2850 mod.setEtaTolerance(self.TOL*0.1)
973 gross 2432
974     BF=Vector(self.P_max,Function(self.dom))
975     for d in range(self.dom.getDim()):
976     for d2 in range(self.dom.getDim()):
977     if d!=d2: BF[d]*=x[d2]
978     v_mask=Vector(0,Solution(self.dom))
979     if free==None:
980     for d in range(self.dom.getDim()):
981     v_mask[d]=whereZero(x[d])+whereZero(x[d]-1.)
982     else:
983     for d in range(self.dom.getDim()):
984     if d == self.dom.getDim()-1:
985     v_mask[d]=whereZero(x[d]-1.)
986     else:
987     v_mask[d]=whereZero(x[d])
988     mod.setExternals(F=BF,fixed_v_mask=v_mask)
989    
990     n=self.dom.getNormal()
991     N_t=0
992     errors=[]
993     while N_t < self.N_dt:
994     t_ref=t+dt
995     v_ref, s_ref,p_ref=self.getReference(t_ref)
996     mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref)
997 gross 2794 # mod.update(dt, eta_iter_max=10, iter_max=50, verbose=self.VERBOSE, usePCG=True, max_correction_steps=30) new version
998     mod.update(dt, iter_max=50, verbose=self.VERBOSE, usePCG=True)
999 gross 2432 self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)
1000     t+=dt
1001     N_t+=1
1002    
1003     def check(self,N_t,mod,t_ref,v_ref, s_ref,p_ref):
1004     p=mod.getPressure()
1005     p-=integrate(p)/vol(self.dom)
1006     error_p=Lsup(mod.getPressure()-p_ref)/Lsup(p_ref)
1007     error_s=Lsup(mod.getDeviatoricStress()-s_ref)/Lsup(s_ref)
1008     error_v=Lsup(mod.getVelocity()-v_ref)/Lsup(v_ref)
1009     error_t=abs(mod.getTime()-t_ref)/abs(t_ref)
1010     if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v
1011 gross 2850 self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
1012     self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )
1013     self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )
1014     self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )
1015 gross 2432 def tearDown(self):
1016     del self.dom
1017    
1018     def test_D2_Fixed_MuNone_LateStart(self):
1019     self.dom = Rectangle(self.NE,self.NE,order=2)
1020     self.mu=None
1021     self.latestart=True
1022     self.runIt()
1023     def test_D2_Fixed_Mu_LateStart(self):
1024     self.dom = Rectangle(self.NE,self.NE,order=2)
1025     self.mu=555.
1026     self.latestart=True
1027     self.runIt()
1028     def test_D2_Fixed_MuNone(self):
1029     self.dom = Rectangle(self.NE,self.NE,order=2)
1030     self.mu=None
1031     self.latestart=False
1032     self.runIt()
1033     def test_D2_Fixed_Mu(self):
1034     self.dom = Rectangle(self.NE,self.NE,order=2)
1035     self.mu=555.
1036     self.latestart=False
1037     self.runIt()
1038     def test_D2_Free_MuNone_LateStart(self):
1039     self.dom = Rectangle(self.NE,self.NE,order=2)
1040     self.mu=None
1041     self.latestart=True
1042     self.runIt(free=0)
1043     def test_D2_Free_Mu_LateStart(self):
1044     self.dom = Rectangle(self.NE,self.NE,order=2)
1045     self.mu=555.
1046     self.latestart=True
1047     self.runIt(free=0)
1048     def test_D2_Free_MuNone(self):
1049     self.dom = Rectangle(self.NE,self.NE,order=2)
1050     self.mu=None
1051     self.latestart=False
1052     self.runIt(free=0)
1053     def test_D2_Free_Mu(self):
1054     self.dom = Rectangle(self.NE,self.NE,order=2)
1055     self.mu=555.
1056     self.latestart=False
1057     self.runIt(free=0)
1058    
1059     def test_D3_Fixed_MuNone_LateStart(self):
1060     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1061     self.mu=None
1062     self.latestart=True
1063     self.runIt()
1064     def test_D3_Fixed_Mu_LateStart(self):
1065     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1066     self.mu=555.
1067     self.latestart=True
1068     self.runIt()
1069     def test_D3_Fixed_MuNone(self):
1070     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1071     self.mu=None
1072     self.latestart=False
1073     self.runIt()
1074     def test_D3_Fixed_Mu(self):
1075     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1076     self.mu=555.
1077     self.latestart=False
1078     self.runIt()
1079     def test_D3_Free_MuNone_LateStart(self):
1080     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1081     self.mu=None
1082     self.latestart=True
1083     self.runIt(free=0)
1084     def test_D3_Free_Mu_LateStart(self):
1085     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1086     self.mu=555.
1087     self.latestart=True
1088     self.runIt(free=0)
1089     def test_D3_Free_MuNone(self):
1090     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1091     self.mu=None
1092     self.latestart=False
1093     self.runIt(free=0)
1094     def test_D3_Free_Mu(self):
1095     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1096     self.mu=555.
1097     self.latestart=False
1098     self.runIt(free=0)
1099    
1100 gross 2647
1101     class Test_FaultSystem(unittest.TestCase):
1102     EPS=1.e-8
1103     NE=10
1104 gross 2663 def test_Fault_MaxValue(self):
1105 gross 2647 dom=Rectangle(2*self.NE,2*self.NE)
1106     x=dom.getX()
1107     f=FaultSystem(dim=2)
1108 gross 2663 f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
1109     f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
1110 gross 2647
1111 gross 2676 u=x[0]*(1.-x[0])*(1-x[1])
1112     t, loc=f.getMaxValue(u)
1113     p=f.getParametrization(x,t)[0]
1114     m, l=loc(u), loc(p)
1115 gross 2647 self.failUnless( m == 0.25, "wrong max value")
1116     self.failUnless( t == 1, "wrong max tag")
1117     self.failUnless( l == 0., "wrong max location")
1118 gross 2676
1119     u=x[1]*(1.-x[1])*(1-x[0])*x[0]
1120     t, loc=f.getMaxValue(u)
1121     p=f.getParametrization(x,t)[0]
1122     m, l=loc(u), loc(p)
1123 gross 2647 self.failUnless( m == 0.0625, "wrong max value")
1124     self.failUnless( t == 2, "wrong max tag")
1125     self.failUnless( l == 0.5, "wrong max location")
1126 gross 2676
1127     u=x[0]*(1.-x[0])*x[1]
1128     t, loc=f.getMaxValue(u)
1129     p=f.getParametrization(x,t)[0]
1130     m, l=loc(u), loc(p)
1131 gross 2647 self.failUnless( m == 0.25, "wrong max value")
1132     self.failUnless( t == 2, "wrong max tag")
1133     self.failUnless( l == 1.0, "wrong max location")
1134 gross 2676
1135     u=x[1]*(1.-x[1])*x[0]
1136     t, loc=f.getMaxValue(u)
1137     p=f.getParametrization(x,t)[0]
1138     m, l=loc(u), loc(p)
1139 gross 2647 self.failUnless( m == 0.25, "wrong max value")
1140     self.failUnless( t == 2, "wrong max tag")
1141     self.failUnless( l == 0., "wrong max location")
1142 gross 2676
1143     u=x[1]*(1.-x[1])*(1.-x[0])
1144     t, loc=f.getMaxValue(u)
1145     p=f.getParametrization(x,t)[0]
1146     m, l=loc(u), loc(p)
1147 gross 2647 self.failUnless( m == 0.25, "wrong max value")
1148     self.failUnless( t == 1, "wrong max tag")
1149     self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong max location")
1150 gross 2675 def test_Fault_MinValue(self):
1151     dom=Rectangle(2*self.NE,2*self.NE)
1152     x=dom.getX()
1153     f=FaultSystem(dim=2)
1154     f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
1155     f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
1156 gross 2647
1157 gross 2676 u=-x[0]*(1.-x[0])*(1-x[1])
1158     t, loc=f.getMinValue(u)
1159     p=f.getParametrization(x,t)[0]
1160     m, l=loc(u), loc(p)
1161 gross 2675 self.failUnless( m == -0.25, "wrong min value")
1162     self.failUnless( t == 1, "wrong min tag")
1163     self.failUnless( l == 0., "wrong min location")
1164 gross 2676 u=-x[1]*(1.-x[1])*(1-x[0])*x[0]
1165     t, loc=f.getMinValue(u)
1166     p=f.getParametrization(x,t)[0]
1167     m, l=loc(u), loc(p)
1168 gross 2675 self.failUnless( m == -0.0625, "wrong min value")
1169     self.failUnless( t == 2, "wrong min tag")
1170     self.failUnless( l == 0.5, "wrong min location")
1171 gross 2676 u=-x[0]*(1.-x[0])*x[1]
1172     t, loc=f.getMinValue(u)
1173     p=f.getParametrization(x,t)[0]
1174     m, l=loc(u), loc(p)
1175 gross 2675 self.failUnless( m == -0.25, "wrong min value")
1176     self.failUnless( t == 2, "wrong min tag")
1177     self.failUnless( l == 1.0, "wrong min location")
1178 gross 2676 u=-x[1]*(1.-x[1])*x[0]
1179     t, loc=f.getMinValue(u)
1180     p=f.getParametrization(x,t)[0]
1181     m, l=loc(u), loc(p)
1182 gross 2675 self.failUnless( m == -0.25, "wrong min value")
1183     self.failUnless( t == 2, "wrong min tag")
1184     self.failUnless( l == 0., "wrong min location")
1185 gross 2676 u=-x[1]*(1.-x[1])*(1.-x[0])
1186     t, loc=f.getMinValue(u)
1187     p=f.getParametrization(x,t)[0]
1188     m, l=loc(u), loc(p)
1189 gross 2675 self.failUnless( m == -0.25, "wrong min value")
1190     self.failUnless( t == 1, "wrong min tag")
1191     self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong min location")
1192    
1193 gross 2647
1194 gross 2663 def test_Fault2D(self):
1195 gross 2647 f=FaultSystem(dim=2)
1196     top1=[ [1.,0.], [1.,1.], [0.,1.] ]
1197 gross 2663 self.failUnlessRaises(ValueError,f.addFault,V0=[1.,0],strikes=[pi/2, pi/2],ls=[1.,1.],tag=1,dips=top1)
1198     f.addFault(V0=[1.,0],strikes=[pi/2, pi],ls=[1.,1.],tag=1)
1199 gross 2647 self.failUnless(f.getDim() == 2, "wrong dimension")
1200     self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
1201 gross 2663 self.failUnless( 2. == f.getTotalLength(1), "length wrong")
1202     self.failUnless( 0. == f.getMediumDepth(1), "depth wrong")
1203 gross 2647 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 range")
1204     self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 range")
1205     self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsets")
1206 gross 2663 segs=f.getTopPolyline(1)
1207 gross 2647 self.failUnless( len(segs) == 3, "wrong number of segments")
1208     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1209     self.failUnless( segs[0].size == 2, "seg 0 has wrong size.")
1210     self.failUnless( numpy.linalg.norm(segs[0]-[1.,0.]) < self.EPS, "wrong vertex. 0 ")
1211     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1212     self.failUnless( segs[1].size == 2, "seg 1 has wrong size.")
1213     self.failUnless( numpy.linalg.norm(segs[1]-[1.,1.]) < self.EPS, "wrong vertex. 1 ")
1214     self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1215     self.failUnless( segs[2].size == 2, "seg 2 has wrong size.")
1216     self.failUnless( numpy.linalg.norm(segs[2]-[0.,1.]) < self.EPS, "wrong vertex. 2 ")
1217     c=f.getCenterOnSurface()
1218     self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1219     self.failUnless( c.size == 2, "center size is wrong")
1220     self.failUnless( numpy.linalg.norm(c-[2./3.,2./3.]) < self.EPS, "center has wrong coordinates.")
1221     o=f.getOrientationOnSurface()/pi*180.
1222 gross 2663 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1223 gross 2647
1224     top2=[ [10.,0.], [0.,10.] ]
1225 gross 2663 f.addFault(V0=[10.,0],strikes=[3.*pi/4],ls=[14.142135623730951], tag=2, w0_offsets=[0,20], w1_max=20)
1226 gross 2647 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1227 gross 2663 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length")
1228     self.failUnless( 0. == f.getMediumDepth(2), "depth wrong")
1229 gross 2647 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range")
1230     self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range")
1231     self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets")
1232 gross 2663 segs=f.getTopPolyline(2)
1233 gross 2647 self.failUnless( len(segs) == 2, "wrong number of segments")
1234     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1235     self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.]) < self.EPS, "wrong vertex. 0 ")
1236     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1237     self.failUnless( numpy.linalg.norm(segs[1]-[0.,10.]) < self.EPS, "wrong vertex. 1 ")
1238     c=f.getCenterOnSurface()
1239     self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1240     self.failUnless( c.size == 2, "center size is wrong")
1241     self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.]) < self.EPS, "center has wrong coordinates.")
1242     o=f.getOrientationOnSurface()/pi*180.
1243 gross 2663 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1244 gross 2647
1245 gross 2654 s,d=f.getSideAndDistance([0.,0.], tag=1)
1246     self.failUnless( s<0, "wrong side.")
1247     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1248     s,d=f.getSideAndDistance([0.,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([1.,2.], 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.,1.], tag=1)
1255     self.failUnless( s>0, "wrong side.")
1256     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1257     s,d=f.getSideAndDistance([2.,0.], tag=1)
1258     self.failUnless( s>0, "wrong side.")
1259     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1260     s,d=f.getSideAndDistance([0.,-1.], tag=1)
1261     self.failUnless( s<0, "wrong side.")
1262     self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
1263     s,d=f.getSideAndDistance([-1.,0], tag=1)
1264     self.failUnless( s<0, "wrong side.")
1265     self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
1266    
1267    
1268 gross 2647 f.transform(rot=-pi/2., shift=[-1.,-1.])
1269     self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
1270 gross 2663 self.failUnless( 2. == f.getTotalLength(1), "length after transformation wrong")
1271     self.failUnless( 0. == f.getMediumDepth(1), "depth after transformation wrong")
1272 gross 2647 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 after transformation range")
1273     self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
1274     self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsetsafter transformation ")
1275 gross 2663 segs=f.getTopPolyline(1)
1276 gross 2647 self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
1277     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1278     self.failUnless( segs[0].size == 2, "seg 0 has wrong size after transformation.")
1279     self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1280     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1281     self.failUnless( segs[1].size == 2, "seg 1 has wrong size after transformation.")
1282     self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1283     self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1284     self.failUnless( segs[2].size == 2, "seg 2 has wrong size after transformation.")
1285     self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.]) < self.EPS, "wrong vertex after transformation. 2 ")
1286 gross 2663 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1287     self.failUnless( 0. == f.getMediumDepth(2), "depth wrong after transformation")
1288 gross 2647 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range after transformation")
1289     self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1290     self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1291 gross 2663 segs=f.getTopPolyline(2)
1292 gross 2647 self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
1293     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1294     self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9]) < self.EPS, "wrong vertex. 0 after transformation")
1295     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1296     self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.]) < self.EPS, "wrong vertex. 1 after transformation")
1297    
1298     c=f.getCenterOnSurface()
1299     self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1300     self.failUnless( c.size == 2, "center size is wrong")
1301     self.failUnless( numpy.linalg.norm(c-[7./5.,-7./5.]) < self.EPS, "center has wrong coordinates.")
1302     o=f.getOrientationOnSurface()/pi*180.
1303 gross 2663 self.failUnless( abs(o-45.) < self.EPS, "wrong orientation.")
1304 gross 2647
1305     p=f.getParametrization([-1.,0.],1)
1306     self.failUnless(p[1]==1., "wrong value.")
1307     self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1308     p=f.getParametrization([-0.5,0.],1)
1309     self.failUnless(p[1]==1., "wrong value.")
1310     self.failUnless(abs(p[0]-0.5)<self.EPS* 0.5, "wrong value.")
1311     p=f.getParametrization([0.,0.],1)
1312     self.failUnless(p[1]==1., "wrong value.")
1313     self.failUnless(abs(p[0]-1.)<self.EPS, "wrong value.")
1314     p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-8)
1315     self.failUnless(p[1]==0., "wrong value.")
1316     p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-6)
1317     self.failUnless(p[1]==1., "wrong value.")
1318     self.failUnless(abs(p[0]-1.0000001)<self.EPS, "wrong value.")
1319     p=f.getParametrization([0.,0.5],1)
1320     self.failUnless(p[1]==1., "wrong value.")
1321     self.failUnless(abs(p[0]-1.5)<self.EPS, "wrong value.")
1322     p=f.getParametrization([0,1.],1)
1323     self.failUnless(p[1]==1., "wrong value.")
1324     self.failUnless(abs(p[0]-2.)<self.EPS, "wrong value.")
1325     p=f.getParametrization([1.,1.],1)
1326     self.failUnless(p[1]==0., "wrong value.")
1327     p=f.getParametrization([0,1.11],1)
1328     self.failUnless(p[1]==0., "wrong value.")
1329     p=f.getParametrization([-1,-9.],2)
1330     self.failUnless(p[1]==1., "wrong value.")
1331     self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1332     p=f.getParametrization([9,1],2)
1333     self.failUnless(p[1]==1., "wrong value.")
1334     self.failUnless(abs(p[0]-20.)<self.EPS, "wrong value.")
1335    
1336 gross 2663 def test_Fault3D(self):
1337     f=FaultSystem(dim=3)
1338     self.failUnless(f.getDim() == 3, "wrong dimension")
1339 gross 2647
1340 gross 2663 top1=[ [0.,0.,0.], [1., 0., 0.] ]
1341     f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=20.,tag=1)
1342     self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
1343     self.failUnless( 1. == f.getTotalLength(1), "length wrong")
1344     self.failUnless( 20. == f.getMediumDepth(1), "depth wrong")
1345     self.failUnless( (0., 1.) == f.getW0Range(1)," wrong W0 range")
1346     self.failUnless( (-20., 0.) == f.getW1Range(1)," wrong W1 range")
1347     self.failUnless( [0., 1.] == f.getW0Offsets(1)," wrong W0 offsets")
1348     segs=f.getTopPolyline(1)
1349     self.failUnless( len(segs) == 2, "wrong number of segments")
1350     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1351     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1352     self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1353     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1354     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1355     self.failUnless( numpy.linalg.norm(segs[1]-[1.,0.,0]) < self.EPS, "wrong vertex. 1 ")
1356     c=f.getCenterOnSurface()
1357     self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1358     self.failUnless( c.size == 3, "center size is wrong")
1359     self.failUnless( numpy.linalg.norm(c-[0.5,0.,0.]) < self.EPS, "center has wrong coordinates.")
1360     o=f.getOrientationOnSurface()/pi*180.
1361     self.failUnless( abs(o) < self.EPS, "wrong orientation.")
1362     d=f.getDips(1)
1363     self.failUnless( len(d) == 1, "wrong number of dips")
1364     self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1365     sn=f.getSegmentNormals(1)
1366     self.failUnless( len(sn) == 1, "wrong number of normals")
1367     self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1368     self.failUnless( numpy.linalg.norm(sn[0]-[0, -1., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1369     dv=f.getDepthVectors(1)
1370     self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1371     self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1372     self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1373     self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1374     self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1375     b=f.getBottomPolyline(1)
1376     self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1377     self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1378     self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1379     self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1380     self.failUnless( numpy.linalg.norm(b[1]-[1., 0., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1381     ds=f.getDepths(1)
1382     self.failUnless( len(ds) == 2, "wrong number of depth")
1383     self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1384     self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1385 gross 2647
1386 gross 2663 top2=[ [0.,0.,0.], [0., 10., 0.] ]
1387     f.addFault(V0=[0.,0,0],strikes=[pi/2],ls=[10.,], dips=pi/2, depths=20.,tag=2)
1388     self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1389     self.failUnless( 10. == f.getTotalLength(2), "length wrong")
1390     self.failUnless( 20. == f.getMediumDepth(2), "depth wrong")
1391     self.failUnless( (0., 10.) == f.getW0Range(2)," wrong W0 range")
1392     self.failUnless( (-20., 0.) == f.getW1Range(2)," wrong W1 range")
1393     self.failUnless( [0., 10.] == f.getW0Offsets(2)," wrong W0 offsets")
1394     segs=f.getTopPolyline(2)
1395     self.failUnless( len(segs) == 2, "wrong number of segments")
1396     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1397     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1398     self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1399     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1400     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1401     self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1402     d=f.getDips(2)
1403     self.failUnless( len(d) == 1, "wrong number of dips")
1404     self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1405     sn=f.getSegmentNormals(2)
1406     self.failUnless( len(sn) == 1, "wrong number of normals")
1407     self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1408     self.failUnless( numpy.linalg.norm(sn[0]-[1, 0., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1409     dv=f.getDepthVectors(2)
1410     self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1411     self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1412     self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1413     self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1414     self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1415     b=f.getBottomPolyline(2)
1416     self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1417     self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1418     self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1419     self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1420     self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1421     ds=f.getDepths(2)
1422     self.failUnless( len(ds) == 2, "wrong number of depth")
1423     self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1424     self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1425 gross 2647
1426 gross 2663 top2=[ [10.,0.,0.], [0., 10., 0.] ]
1427     f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/2, depths=30.,tag=2)
1428     self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1429     self.failUnless( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1430     self.failUnless( 30. == f.getMediumDepth(2), "depth wrong")
1431     self.failUnless( (-30., 0.) == f.getW1Range(2)," wrong W1 range")
1432     segs=f.getTopPolyline(2)
1433     self.failUnless( len(segs) == 2, "wrong number of segments")
1434     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1435     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1436     self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1437     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1438     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1439     self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1440     d=f.getDips(2)
1441     self.failUnless( len(d) == 1, "wrong number of dips")
1442     self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1443     sn=f.getSegmentNormals(2)
1444     self.failUnless( len(sn) == 1, "wrong number of normals")
1445     self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1446     self.failUnless( numpy.linalg.norm(sn[0]-[0.70710678118654746, 0.70710678118654746, 0.]) < self.EPS, "wrong bottom vertex 1 ")
1447     dv=f.getDepthVectors(2)
1448     self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1449     self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1450     self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -30.]) < self.EPS, "wrong depth vector 0 ")
1451     self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1452     self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -30.]) < self.EPS, "wrong depth vector 1 ")
1453     b=f.getBottomPolyline(2)
1454     self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1455     self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1456     self.failUnless( numpy.linalg.norm(b[0]-[10., 0., -30.]) < self.EPS, "wrong bottom vertex 0 ")
1457     self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1458     self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -30.]) < self.EPS, "wrong bottom vertex 1 ")
1459     ds=f.getDepths(2)
1460     self.failUnless( len(ds) == 2, "wrong number of depth")
1461     self.failUnless( abs(ds[0]-30.) < self.EPS, "wrong depth at vertex 0 ")
1462     self.failUnless( abs(ds[1]-30.) < self.EPS, "wrong depth at vertex 1 ")
1463    
1464     top2=[ [10.,0.,0.], [0., 10., 0.] ]
1465     f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/4, depths=50.,tag=2)
1466     self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1467     self.failUnless( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1468     self.failUnless( 50. == f.getMediumDepth(2), "depth wrong")
1469     self.failUnless( (-50., 0.) == f.getW1Range(2)," wrong W1 range")
1470     segs=f.getTopPolyline(2)
1471     self.failUnless( len(segs) == 2, "wrong number of segments")
1472     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1473     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1474     self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1475     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1476     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1477     self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1478     d=f.getDips(2)
1479     self.failUnless( len(d) == 1, "wrong number of dips")
1480     self.failUnless( abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1481     sn=f.getSegmentNormals(2)
1482     self.failUnless( len(sn) == 1, "wrong number of normals")
1483     self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1484     self.failUnless( numpy.linalg.norm(sn[0]-[0.5,0.5,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1485     dv=f.getDepthVectors(2)
1486     self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1487     self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1488     self.failUnless( numpy.linalg.norm(dv[0]-[25., 25., -35.355339059327378]) < self.EPS, "wrong depth vector 0 ")
1489     self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1490     self.failUnless( numpy.linalg.norm(dv[1]-[25.,25., -35.355339059327378]) < self.EPS, "wrong depth vector 1 ")
1491     b=f.getBottomPolyline(2)
1492     self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1493     self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1494     self.failUnless( numpy.linalg.norm(b[0]-[35., 25., -35.355339059327378]) < self.EPS, "wrong bottom vertex 0 ")
1495     self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1496     self.failUnless( numpy.linalg.norm(b[1]-[25, 35., -35.355339059327378]) < self.EPS, "wrong bottom vertex 1 ")
1497     ds=f.getDepths(2)
1498     self.failUnless( len(ds) == 2, "wrong number of depth")
1499     self.failUnless( abs(ds[0]-50.) < self.EPS, "wrong depth at vertex 0 ")
1500     self.failUnless( abs(ds[1]-50.) < self.EPS, "wrong depth at vertex 1 ")
1501    
1502     top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
1503     f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=1, dips=pi/4, depths=20.)
1504     self.failUnless( 20. == f.getTotalLength(1), "length wrong")
1505     self.failUnless( 20. == f.getMediumDepth(1), "depth wrong")
1506     segs=f.getTopPolyline(1)
1507     self.failUnless( len(segs) == 3, "wrong number of segments")
1508     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1509     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1510     self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1511     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1512     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1513     self.failUnless( numpy.linalg.norm(segs[1]-[10.,10.,0.]) < self.EPS, "wrong vertex. 1 ")
1514     self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1515     self.failUnless( segs[2].size == 3, "seg 2 has wrong size.")
1516