/[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 3501 - (hide annotations)
Wed Apr 13 04:07:53 2011 UTC (8 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 86914 byte(s)


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