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