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