/[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 3527 - (hide annotations)
Tue May 24 06:59:20 2011 UTC (8 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 73991 byte(s)
modifications to the coal seam gas 
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 3502 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 gross 1471
403 artak 2234 class Test_Mountains3D(unittest.TestCase):
404     def setUp(self):
405     NE=16
406     self.TOL=1e-4
407     self.domain=Brick(NE,NE,NE,order=1,useFullElementOrder=True)
408     def tearDown(self):
409     del self.domain
410     def test_periodic(self):
411     EPS=0.01
412 gross 2100
413 artak 2234 x=self.domain.getX()
414     v = Vector(0.0, Solution(self.domain))
415     a0=1
416     a1=1
417     n0=2
418     n1=2
419     n2=0.5
420     a2=-(a0*n0+a1*n1)/n2
421     v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])* cos(pi*n2*x[2])
422     v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])* cos(pi*n2*x[2])
423     v[2]=a2*cos(pi*n0*x[0])* cos(pi*n1*x[1])* sin(pi*n2*x[2])
424    
425 gross 2563 mts=Mountains(self.domain,eps=EPS)
426     mts.setVelocity(v)
427     Z=mts.update()
428 artak 2234
429 gross 2564 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[2]-1.)))
430 artak 2234 self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
431    
432     class Test_Mountains2D(unittest.TestCase):
433     def setUp(self):
434     NE=16
435     self.TOL=1e-4
436     self.domain=Rectangle(NE,NE,order=1,useFullElementOrder=True)
437     def tearDown(self):
438     del self.domain
439     def test_periodic(self):
440     EPS=0.01
441    
442     x=self.domain.getX()
443     v = Vector(0.0, Solution(self.domain))
444     a0=1
445     n0=1
446     n1=0.5
447     a1=-(a0*n0)/n1
448     v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])
449     v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])
450    
451     H_t=Scalar(0.0, Solution(self.domain))
452 gross 2563 mts=Mountains(self.domain,eps=EPS)
453     mts.setVelocity(v)
454     Z=mts.update()
455 artak 2234
456 gross 2564 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[1]-1.)))
457 artak 2234 self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
458    
459 gross 2301
460    
461     class Test_Rheologies(unittest.TestCase):
462     """
463     this is the program used to generate the powerlaw tests:
464    
465     TAU_Y=100.
466     N=10
467     M=5
468    
469     def getE(tau):
470     if tau<=TAU_Y:
471     return 1./(0.5+20*sqrt(tau))
472     else:
473     raise ValueError,"out of range."
474     tau=[ i* (TAU_Y/N) for i in range(N+1)] + [TAU_Y for i in range(M)]
475     e=[ tau[i]/getE(tau[i]) for i in range(N+1)]
476     e+= [ (i+1+j)* (max(e)/N) for j in range(M) ]
477    
478     print tau
479     print e
480     """
481     TOL=1.e-8
482     def test_PowerLaw_Init(self):
483     pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
484    
485     self.failUnlessEqual(pl.getNumMaterials(),3,"num materials is wrong")
486     self.failUnlessEqual(pl.validMaterialId(0),True,"material id 0 not found")
487     self.failUnlessEqual(pl.validMaterialId(1),True,"material id 1 not found")
488     self.failUnlessEqual(pl.validMaterialId(2),True,"material id 2 not found")
489     self.failUnlessEqual(pl.validMaterialId(3),False,"material id 3 not found")
490    
491     self.failUnlessEqual(pl.getFriction(),None,"initial friction wrong.")
492     self.failUnlessEqual(pl.getTauY(),None,"initial tau y wrong.")
493     pl.setDruckerPragerLaw(tau_Y=10,friction=3)
494     self.failUnlessEqual(pl.getFriction(),3,"friction wrong.")
495     self.failUnlessEqual(pl.getTauY(),10,"tau y wrong.")
496    
497     self.failUnlessEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
498     pl.setElasticShearModulus(1000)
499     self.failUnlessEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
500    
501     e=pl.getEtaN()
502     self.failUnlessEqual(len(e),3,"initial length of etaN is wrong.")
503     self.failUnlessEqual(e,[None, None, None],"initial etaN are wrong.")
504     self.failUnlessEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
505     self.failUnlessEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
506     self.failUnlessEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
507     self.failUnlessRaises(ValueError, pl.getEtaN, 3)
508    
509     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
510     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
511     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
512    
513     pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
514     self.failUnlessEqual(pl.getPower(),[1,2,3],"powers are wrong.")
515     self.failUnlessEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
516     self.failUnlessEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
517    
518     pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
519     self.failUnlessEqual(pl.getPower(),[4,2,3],"powers are wrong.")
520     self.failUnlessEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
521     self.failUnlessEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
522    
523     self.failUnlessRaises(ValueError,pl.getPower,-1)
524     self.failUnlessEqual(pl.getPower(0),4,"power 0 is wrong.")
525     self.failUnlessEqual(pl.getPower(1),2,"power 1 is wrong.")
526     self.failUnlessEqual(pl.getPower(2),3,"power 2 is wrong.")
527     self.failUnlessRaises(ValueError,pl.getPower,3)
528    
529     self.failUnlessRaises(ValueError,pl.getTauT,-1)
530     self.failUnlessEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
531     self.failUnlessEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
532     self.failUnlessEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
533     self.failUnlessRaises(ValueError,pl.getTauT,3)
534    
535     self.failUnlessRaises(ValueError,pl.getEtaN,-1)
536     self.failUnlessEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
537     self.failUnlessEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
538     self.failUnlessEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
539     self.failUnlessRaises(ValueError,pl.getEtaN,3)
540    
541     def checkResult(self,id,gamma_dot_, eta, tau_ref):
542     self.failUnless(eta>=0,"eta needs to be positive (test %s)"%id)
543     error=abs(gamma_dot_*eta-tau_ref)
544     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))
545    
546     def test_PowerLaw_Linear(self):
547     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]
548     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]
549     pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
550     pl.setDruckerPragerLaw(tau_Y=100.)
551     pl.setPowerLaw(eta_N=2.)
552     pl.setEtaTolerance(self.TOL)
553     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
554    
555     def test_PowerLaw_QuadLarge(self):
556     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]
557 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]
558 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
559     pl.setDruckerPragerLaw(tau_Y=100.)
560     pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2])
561     pl.setEtaTolerance(self.TOL)
562     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
563    
564     def test_PowerLaw_QuadSmall(self):
565     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]
566 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]
567 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
568     pl.setDruckerPragerLaw(tau_Y=100.)
569     pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2])
570     pl.setEtaTolerance(self.TOL)
571     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
572    
573     def test_PowerLaw_CubeLarge(self):
574     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]
575 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]
576 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
577     pl.setDruckerPragerLaw(tau_Y=100.)
578     pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3])
579     pl.setEtaTolerance(self.TOL)
580     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
581    
582     def test_PowerLaw_CubeSmall(self):
583     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]
584 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]
585 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
586     pl.setDruckerPragerLaw(tau_Y=100.)
587     pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3])
588     pl.setEtaTolerance(self.TOL)
589     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
590    
591     def test_PowerLaw_QuadLarge_CubeLarge(self):
592     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]
593 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]
594    
595 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
596     pl.setDruckerPragerLaw(tau_Y=100.)
597     pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
598     pl.setEtaTolerance(self.TOL)
599     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
600    
601     def test_PowerLaw_QuadLarge_CubeSmall(self):
602     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]
603 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]
604    
605 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
606     pl.setDruckerPragerLaw(tau_Y=100.)
607     pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
608     pl.setEtaTolerance(self.TOL)
609     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
610    
611     def test_PowerLaw_QuadSmall_CubeLarge(self):
612     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]
613 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]
614 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
615     pl.setDruckerPragerLaw(tau_Y=100.)
616     pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
617     pl.setEtaTolerance(self.TOL)
618     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
619    
620     def test_PowerLaw_QuadSmall_CubeSmall(self):
621     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]
622 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]
623 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
624     pl.setDruckerPragerLaw(tau_Y=100.)
625     pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
626     pl.setEtaTolerance(self.TOL)
627     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
628    
629     def test_PowerLaw_withShear(self):
630     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]
631     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]
632     pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
633     pl.setDruckerPragerLaw(tau_Y=100.)
634     pl.setPowerLaw(eta_N=2.)
635     pl.setElasticShearModulus(3.)
636     dt=1./3.
637     pl.setEtaTolerance(self.TOL)
638     self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
639     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
640    
641 gross 2432 class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase):
642 gross 2850 TOL=1.e-5
643     VERBOSE=False # or True
644 gross 2432 A=1.
645     P_max=100
646 gross 2794 NE=2*getMPISizeWorld()
647 gross 2432 tau_Y=10.
648     N_dt=10
649    
650     # material parameter:
651     tau_1=5.
652     tau_2=5.
653     eta_0=100.
654     eta_1=50.
655     eta_2=400.
656     N_1=2.
657     N_2=3.
658     def getReference(self, t):
659    
660     B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
661     x=self.dom.getX()
662    
663     s_00=min(self.A*t,B)
664     tau=sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)*abs(s_00)
665 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.)
666 gross 2432
667     alpha=0.5*inv_eta*s_00
668     if s_00 <= B and self.mu !=None: alpha+=1./(2*self.mu)*self.A
669     u_ref=x*alpha
670     u_ref[self.dom.getDim()-1]=(1.-x[self.dom.getDim()-1])*alpha*(self.dom.getDim()-1)
671     sigma_ref=kronecker(self.dom)*s_00
672     sigma_ref[self.dom.getDim()-1,self.dom.getDim()-1]=-s_00*(self.dom.getDim()-1)
673    
674     p_ref=self.P_max
675     for d in range(self.dom.getDim()): p_ref=p_ref*x[d]
676     p_ref-=integrate(p_ref)/vol(self.dom)
677     return u_ref, sigma_ref, p_ref
678    
679     def runIt(self, free=None):
680     x=self.dom.getX()
681     B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
682     dt=B/int(self.N_dt/2)
683     if self.VERBOSE: print "dt =",dt
684     if self.latestart:
685     t=dt
686     else:
687     t=0
688     v,s,p=self.getReference(t)
689    
690     mod=IncompressibleIsotropicFlowCartesian(self.dom, stress=s, v=v, p=p, t=t, numMaterials=3, verbose=self.VERBOSE)
691     mod.setDruckerPragerLaw(tau_Y=self.tau_Y,friction=None)
692     mod.setElasticShearModulus(self.mu)
693     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])
694     mod.setTolerance(self.TOL)
695 gross 2850 mod.setEtaTolerance(self.TOL*0.1)
696 gross 2432
697     BF=Vector(self.P_max,Function(self.dom))
698     for d in range(self.dom.getDim()):
699     for d2 in range(self.dom.getDim()):
700     if d!=d2: BF[d]*=x[d2]
701     v_mask=Vector(0,Solution(self.dom))
702     if free==None:
703     for d in range(self.dom.getDim()):
704     v_mask[d]=whereZero(x[d])+whereZero(x[d]-1.)
705     else:
706     for d in range(self.dom.getDim()):
707     if d == self.dom.getDim()-1:
708     v_mask[d]=whereZero(x[d]-1.)
709     else:
710     v_mask[d]=whereZero(x[d])
711     mod.setExternals(F=BF,fixed_v_mask=v_mask)
712    
713     n=self.dom.getNormal()
714     N_t=0
715     errors=[]
716     while N_t < self.N_dt:
717     t_ref=t+dt
718     v_ref, s_ref,p_ref=self.getReference(t_ref)
719     mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref)
720 gross 2794 # mod.update(dt, eta_iter_max=10, iter_max=50, verbose=self.VERBOSE, usePCG=True, max_correction_steps=30) new version
721     mod.update(dt, iter_max=50, verbose=self.VERBOSE, usePCG=True)
722 gross 2432 self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)
723     t+=dt
724     N_t+=1
725    
726     def check(self,N_t,mod,t_ref,v_ref, s_ref,p_ref):
727     p=mod.getPressure()
728     p-=integrate(p)/vol(self.dom)
729     error_p=Lsup(mod.getPressure()-p_ref)/Lsup(p_ref)
730     error_s=Lsup(mod.getDeviatoricStress()-s_ref)/Lsup(s_ref)
731     error_v=Lsup(mod.getVelocity()-v_ref)/Lsup(v_ref)
732     error_t=abs(mod.getTime()-t_ref)/abs(t_ref)
733     if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v
734 gross 2850 self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
735     self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )
736     self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )
737     self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )
738 gross 2432 def tearDown(self):
739     del self.dom
740    
741     def test_D2_Fixed_MuNone_LateStart(self):
742     self.dom = Rectangle(self.NE,self.NE,order=2)
743     self.mu=None
744     self.latestart=True
745     self.runIt()
746     def test_D2_Fixed_Mu_LateStart(self):
747     self.dom = Rectangle(self.NE,self.NE,order=2)
748     self.mu=555.
749     self.latestart=True
750     self.runIt()
751     def test_D2_Fixed_MuNone(self):
752     self.dom = Rectangle(self.NE,self.NE,order=2)
753     self.mu=None
754     self.latestart=False
755     self.runIt()
756     def test_D2_Fixed_Mu(self):
757     self.dom = Rectangle(self.NE,self.NE,order=2)
758     self.mu=555.
759     self.latestart=False
760     self.runIt()
761     def test_D2_Free_MuNone_LateStart(self):
762     self.dom = Rectangle(self.NE,self.NE,order=2)
763     self.mu=None
764     self.latestart=True
765     self.runIt(free=0)
766     def test_D2_Free_Mu_LateStart(self):
767     self.dom = Rectangle(self.NE,self.NE,order=2)
768     self.mu=555.
769     self.latestart=True
770     self.runIt(free=0)
771     def test_D2_Free_MuNone(self):
772     self.dom = Rectangle(self.NE,self.NE,order=2)
773     self.mu=None
774     self.latestart=False
775     self.runIt(free=0)
776     def test_D2_Free_Mu(self):
777     self.dom = Rectangle(self.NE,self.NE,order=2)
778     self.mu=555.
779     self.latestart=False
780     self.runIt(free=0)
781    
782     def test_D3_Fixed_MuNone_LateStart(self):
783     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
784     self.mu=None
785     self.latestart=True
786     self.runIt()
787     def test_D3_Fixed_Mu_LateStart(self):
788     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
789     self.mu=555.
790     self.latestart=True
791     self.runIt()
792     def test_D3_Fixed_MuNone(self):
793     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
794     self.mu=None
795     self.latestart=False
796     self.runIt()
797     def test_D3_Fixed_Mu(self):
798     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
799     self.mu=555.
800     self.latestart=False
801     self.runIt()
802     def test_D3_Free_MuNone_LateStart(self):
803     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
804     self.mu=None
805     self.latestart=True
806     self.runIt(free=0)
807     def test_D3_Free_Mu_LateStart(self):
808     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
809     self.mu=555.
810     self.latestart=True
811     self.runIt(free=0)
812     def test_D3_Free_MuNone(self):
813     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
814     self.mu=None
815     self.latestart=False
816     self.runIt(free=0)
817     def test_D3_Free_Mu(self):
818     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
819     self.mu=555.
820     self.latestart=False
821     self.runIt(free=0)
822    
823 gross 2647
824     class Test_FaultSystem(unittest.TestCase):
825     EPS=1.e-8
826     NE=10
827 gross 2663 def test_Fault_MaxValue(self):
828 gross 2647 dom=Rectangle(2*self.NE,2*self.NE)
829     x=dom.getX()
830     f=FaultSystem(dim=2)
831 gross 2663 f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
832     f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
833 gross 2647
834 gross 2676 u=x[0]*(1.-x[0])*(1-x[1])
835     t, loc=f.getMaxValue(u)
836     p=f.getParametrization(x,t)[0]
837     m, l=loc(u), loc(p)
838 gross 2647 self.failUnless( m == 0.25, "wrong max value")
839     self.failUnless( t == 1, "wrong max tag")
840     self.failUnless( l == 0., "wrong max location")
841 gross 2676
842     u=x[1]*(1.-x[1])*(1-x[0])*x[0]
843     t, loc=f.getMaxValue(u)
844     p=f.getParametrization(x,t)[0]
845     m, l=loc(u), loc(p)
846 gross 2647 self.failUnless( m == 0.0625, "wrong max value")
847     self.failUnless( t == 2, "wrong max tag")
848     self.failUnless( l == 0.5, "wrong max location")
849 gross 2676
850     u=x[0]*(1.-x[0])*x[1]
851     t, loc=f.getMaxValue(u)
852     p=f.getParametrization(x,t)[0]
853     m, l=loc(u), loc(p)
854 gross 2647 self.failUnless( m == 0.25, "wrong max value")
855     self.failUnless( t == 2, "wrong max tag")
856     self.failUnless( l == 1.0, "wrong max location")
857 gross 2676
858     u=x[1]*(1.-x[1])*x[0]
859     t, loc=f.getMaxValue(u)
860     p=f.getParametrization(x,t)[0]
861     m, l=loc(u), loc(p)
862 gross 2647 self.failUnless( m == 0.25, "wrong max value")
863     self.failUnless( t == 2, "wrong max tag")
864     self.failUnless( l == 0., "wrong max location")
865 gross 2676
866     u=x[1]*(1.-x[1])*(1.-x[0])
867     t, loc=f.getMaxValue(u)
868     p=f.getParametrization(x,t)[0]
869     m, l=loc(u), loc(p)
870 gross 2647 self.failUnless( m == 0.25, "wrong max value")
871     self.failUnless( t == 1, "wrong max tag")
872     self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong max location")
873 gross 2675 def test_Fault_MinValue(self):
874     dom=Rectangle(2*self.NE,2*self.NE)
875     x=dom.getX()
876     f=FaultSystem(dim=2)
877     f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
878     f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
879 gross 2647
880 gross 2676 u=-x[0]*(1.-x[0])*(1-x[1])
881     t, loc=f.getMinValue(u)
882     p=f.getParametrization(x,t)[0]
883     m, l=loc(u), loc(p)
884 gross 2675 self.failUnless( m == -0.25, "wrong min value")
885     self.failUnless( t == 1, "wrong min tag")
886     self.failUnless( l == 0., "wrong min location")
887 gross 2676 u=-x[1]*(1.-x[1])*(1-x[0])*x[0]
888     t, loc=f.getMinValue(u)
889     p=f.getParametrization(x,t)[0]
890     m, l=loc(u), loc(p)
891 gross 2675 self.failUnless( m == -0.0625, "wrong min value")
892     self.failUnless( t == 2, "wrong min tag")
893     self.failUnless( l == 0.5, "wrong min location")
894 gross 2676 u=-x[0]*(1.-x[0])*x[1]
895     t, loc=f.getMinValue(u)
896     p=f.getParametrization(x,t)[0]
897     m, l=loc(u), loc(p)
898 gross 2675 self.failUnless( m == -0.25, "wrong min value")
899     self.failUnless( t == 2, "wrong min tag")
900     self.failUnless( l == 1.0, "wrong min location")
901 gross 2676 u=-x[1]*(1.-x[1])*x[0]
902     t, loc=f.getMinValue(u)
903     p=f.getParametrization(x,t)[0]
904     m, l=loc(u), loc(p)
905 gross 2675 self.failUnless( m == -0.25, "wrong min value")
906     self.failUnless( t == 2, "wrong min tag")
907     self.failUnless( l == 0., "wrong min location")
908 gross 2676 u=-x[1]*(1.-x[1])*(1.-x[0])
909     t, loc=f.getMinValue(u)
910     p=f.getParametrization(x,t)[0]
911     m, l=loc(u), loc(p)
912 gross 2675 self.failUnless( m == -0.25, "wrong min value")
913     self.failUnless( t == 1, "wrong min tag")
914     self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong min location")
915    
916 gross 2647
917 gross 2663 def test_Fault2D(self):
918 gross 2647 f=FaultSystem(dim=2)
919     top1=[ [1.,0.], [1.,1.], [0.,1.] ]
920 gross 2663 self.failUnlessRaises(ValueError,f.addFault,V0=[1.,0],strikes=[pi/2, pi/2],ls=[1.,1.],tag=1,dips=top1)
921     f.addFault(V0=[1.,0],strikes=[pi/2, pi],ls=[1.,1.],tag=1)
922 gross 2647 self.failUnless(f.getDim() == 2, "wrong dimension")
923     self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
924 gross 2663 self.failUnless( 2. == f.getTotalLength(1), "length wrong")
925     self.failUnless( 0. == f.getMediumDepth(1), "depth wrong")
926 gross 2647 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 range")
927     self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 range")
928     self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsets")
929 gross 2663 segs=f.getTopPolyline(1)
930 gross 2647 self.failUnless( len(segs) == 3, "wrong number of segments")
931     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
932     self.failUnless( segs[0].size == 2, "seg 0 has wrong size.")
933     self.failUnless( numpy.linalg.norm(segs[0]-[1.,0.]) < self.EPS, "wrong vertex. 0 ")
934     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
935     self.failUnless( segs[1].size == 2, "seg 1 has wrong size.")
936     self.failUnless( numpy.linalg.norm(segs[1]-[1.,1.]) < self.EPS, "wrong vertex. 1 ")
937     self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
938     self.failUnless( segs[2].size == 2, "seg 2 has wrong size.")
939     self.failUnless( numpy.linalg.norm(segs[2]-[0.,1.]) < self.EPS, "wrong vertex. 2 ")
940     c=f.getCenterOnSurface()
941     self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
942     self.failUnless( c.size == 2, "center size is wrong")
943     self.failUnless( numpy.linalg.norm(c-[2./3.,2./3.]) < self.EPS, "center has wrong coordinates.")
944     o=f.getOrientationOnSurface()/pi*180.
945 gross 2663 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
946 gross 2647
947     top2=[ [10.,0.], [0.,10.] ]
948 gross 2663 f.addFault(V0=[10.,0],strikes=[3.*pi/4],ls=[14.142135623730951], tag=2, w0_offsets=[0,20], w1_max=20)
949 gross 2647 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
950 gross 2663 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length")
951     self.failUnless( 0. == f.getMediumDepth(2), "depth wrong")
952 gross 2647 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range")
953     self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range")
954     self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets")
955 gross 2663 segs=f.getTopPolyline(2)
956 gross 2647 self.failUnless( len(segs) == 2, "wrong number of segments")
957     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
958     self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.]) < self.EPS, "wrong vertex. 0 ")
959     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
960     self.failUnless( numpy.linalg.norm(segs[1]-[0.,10.]) < self.EPS, "wrong vertex. 1 ")
961     c=f.getCenterOnSurface()
962     self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
963     self.failUnless( c.size == 2, "center size is wrong")
964     self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.]) < self.EPS, "center has wrong coordinates.")
965     o=f.getOrientationOnSurface()/pi*180.
966 gross 2663 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
967 gross 2647
968 gross 2654 s,d=f.getSideAndDistance([0.,0.], tag=1)
969     self.failUnless( s<0, "wrong side.")
970     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
971     s,d=f.getSideAndDistance([0.,2.], tag=1)
972     self.failUnless( s>0, "wrong side.")
973     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
974     s,d=f.getSideAndDistance([1.,2.], tag=1)
975     self.failUnless( s>0, "wrong side.")
976     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
977     s,d=f.getSideAndDistance([2.,1.], tag=1)
978     self.failUnless( s>0, "wrong side.")
979     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
980     s,d=f.getSideAndDistance([2.,0.], tag=1)
981     self.failUnless( s>0, "wrong side.")
982     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
983     s,d=f.getSideAndDistance([0.,-1.], tag=1)
984     self.failUnless( s<0, "wrong side.")
985     self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
986     s,d=f.getSideAndDistance([-1.,0], tag=1)
987     self.failUnless( s<0, "wrong side.")
988     self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
989    
990    
991 gross 2647 f.transform(rot=-pi/2., shift=[-1.,-1.])
992     self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
993 gross 2663 self.failUnless( 2. == f.getTotalLength(1), "length after transformation wrong")
994     self.failUnless( 0. == f.getMediumDepth(1), "depth after transformation wrong")
995 gross 2647 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 after transformation range")
996     self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
997     self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsetsafter transformation ")
998 gross 2663 segs=f.getTopPolyline(1)
999 gross 2647 self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
1000     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1001     self.failUnless( segs[0].size == 2, "seg 0 has wrong size after transformation.")
1002     self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1003     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1004     self.failUnless( segs[1].size == 2, "seg 1 has wrong size after transformation.")
1005     self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1006     self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1007     self.failUnless( segs[2].size == 2, "seg 2 has wrong size after transformation.")
1008     self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.]) < self.EPS, "wrong vertex after transformation. 2 ")
1009 gross 2663 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1010     self.failUnless( 0. == f.getMediumDepth(2), "depth wrong after transformation")
1011 gross 2647 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range after transformation")
1012     self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1013     self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1014 gross 2663 segs=f.getTopPolyline(2)
1015 gross 2647 self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
1016     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1017     self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9]) < self.EPS, "wrong vertex. 0 after transformation")
1018     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1019     self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.]) < self.EPS, "wrong vertex. 1 after transformation")
1020    
1021     c=f.getCenterOnSurface()
1022     self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1023     self.failUnless( c.size == 2, "center size is wrong")
1024     self.failUnless( numpy.linalg.norm(c-[7./5.,-7./5.]) < self.EPS, "center has wrong coordinates.")
1025     o=f.getOrientationOnSurface()/pi*180.
1026 gross 2663 self.failUnless( abs(o-45.) < self.EPS, "wrong orientation.")
1027 gross 2647
1028     p=f.getParametrization([-1.,0.],1)
1029     self.failUnless(p[1]==1., "wrong value.")
1030     self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1031     p=f.getParametrization([-0.5,0.],1)
1032     self.failUnless(p[1]==1., "wrong value.")
1033     self.failUnless(abs(p[0]-0.5)<self.EPS* 0.5, "wrong value.")
1034     p=f.getParametrization([0.,0.],1)
1035     self.failUnless(p[1]==1., "wrong value.")
1036     self.failUnless(abs(p[0]-1.)<self.EPS, "wrong value.")
1037     p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-8)
1038     self.failUnless(p[1]==0., "wrong value.")
1039     p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-6)
1040     self.failUnless(p[1]==1., "wrong value.")
1041     self.failUnless(abs(p[0]-1.0000001)<self.EPS, "wrong value.")
1042     p=f.getParametrization([0.,0.5],1)
1043     self.failUnless(p[1]==1., "wrong value.")
1044     self.failUnless(abs(p[0]-1.5)<self.EPS, "wrong value.")
1045     p=f.getParametrization([0,1.],1)
1046     self.failUnless(p[1]==1., "wrong value.")
1047     self.failUnless(abs(p[0]-2.)<self.EPS, "wrong value.")
1048     p=f.getParametrization([1.,1.],1)
1049     self.failUnless(p[1]==0., "wrong value.")
1050     p=f.getParametrization([0,1.11],1)
1051     self.failUnless(p[1]==0., "wrong value.")
1052     p=f.getParametrization([-1,-9.],2)
1053     self.failUnless(p[1]==1., "wrong value.")
1054     self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1055     p=f.getParametrization([9,1],2)
1056     self.failUnless(p[1]==1., "wrong value.")
1057     self.failUnless(abs(p[0]-20.)<self.EPS, "wrong value.")
1058    
1059 gross 2663 def test_Fault3D(self):
1060     f=FaultSystem(dim=3)
1061     self.failUnless(f.getDim() == 3, "wrong dimension")
1062 gross 2647
1063 gross 2663 top1=[ [0.,0.,0.], [1., 0., 0.] ]
1064     f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=20.,tag=1)
1065     self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
1066     self.failUnless( 1. == f.getTotalLength(1), "length wrong")
1067     self.failUnless( 20. == f.getMediumDepth(1), "depth wrong")
1068     self.failUnless( (0., 1.) == f.getW0Range(1)," wrong W0 range")
1069     self.failUnless( (-20., 0.) == f.getW1Range(1)," wrong W1 range")
1070     self.failUnless( [0., 1.] == f.getW0Offsets(1)," wrong W0 offsets")
1071     segs=f.getTopPolyline(1)
1072     self.failUnless( len(segs) == 2, "wrong number of segments")
1073     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1074     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1075     self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1076     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1077     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1078     self.failUnless( numpy.linalg.norm(segs[1]-[1.,0.,0]) < self.EPS, "wrong vertex. 1 ")
1079     c=f.getCenterOnSurface()
1080     self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1081     self.failUnless( c.size == 3, "center size is wrong")
1082     self.failUnless( numpy.linalg.norm(c-[0.5,0.,0.]) < self.EPS, "center has wrong coordinates.")
1083     o=f.getOrientationOnSurface()/pi*180.
1084     self.failUnless( abs(o) < self.EPS, "wrong orientation.")
1085     d=f.getDips(1)
1086     self.failUnless( len(d) == 1, "wrong number of dips")
1087     self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1088     sn=f.getSegmentNormals(1)
1089     self.failUnless( len(sn) == 1, "wrong number of normals")
1090     self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1091     self.failUnless( numpy.linalg.norm(sn[0]-[0, -1., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1092     dv=f.getDepthVectors(1)
1093     self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1094     self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1095     self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1096     self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1097     self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1098     b=f.getBottomPolyline(1)
1099     self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1100     self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1101     self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1102     self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1103     self.failUnless( numpy.linalg.norm(b[1]-[1., 0., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1104     ds=f.getDepths(1)
1105     self.failUnless( len(ds) == 2, "wrong number of depth")
1106     self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1107     self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1108 gross 2647
1109 gross 2663 top2=[ [0.,0.,0.], [0., 10., 0.] ]
1110     f.addFault(V0=[0.,0,0],strikes=[pi/2],ls=[10.,], dips=pi/2, depths=20.,tag=2)
1111     self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1112     self.failUnless( 10. == f.getTotalLength(2), "length wrong")
1113     self.failUnless( 20. == f.getMediumDepth(2), "depth wrong")
1114     self.failUnless( (0., 10.) == f.getW0Range(2)," wrong W0 range")
1115     self.failUnless( (-20., 0.) == f.getW1Range(2)," wrong W1 range")
1116     self.failUnless( [0., 10.] == f.getW0Offsets(2)," wrong W0 offsets")
1117     segs=f.getTopPolyline(2)
1118     self.failUnless( len(segs) == 2, "wrong number of segments")
1119     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1120     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1121     self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1122     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1123     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1124     self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1125     d=f.getDips(2)
1126     self.failUnless( len(d) == 1, "wrong number of dips")
1127     self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1128     sn=f.getSegmentNormals(2)
1129     self.failUnless( len(sn) == 1, "wrong number of normals")
1130     self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1131     self.failUnless( numpy.linalg.norm(sn[0]-[1, 0., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1132     dv=f.getDepthVectors(2)
1133     self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1134     self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1135     self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1136     self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1137     self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1138     b=f.getBottomPolyline(2)
1139     self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1140     self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1141     self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1142     self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1143     self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1144     ds=f.getDepths(2)
1145     self.failUnless( len(ds) == 2, "wrong number of depth")
1146     self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1147     self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1148 gross 2647
1149 gross 2663 top2=[ [10.,0.,0.], [0., 10., 0.] ]
1150     f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/2, depths=30.,tag=2)
1151     self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1152     self.failUnless( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1153     self.failUnless( 30. == f.getMediumDepth(2), "depth wrong")
1154     self.failUnless( (-30., 0.) == f.getW1Range(2)," wrong W1 range")
1155     segs=f.getTopPolyline(2)
1156     self.failUnless( len(segs) == 2, "wrong number of segments")
1157     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1158     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1159     self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1160     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1161     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1162     self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1163     d=f.getDips(2)
1164     self.failUnless( len(d) == 1, "wrong number of dips")
1165     self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1166     sn=f.getSegmentNormals(2)
1167     self.failUnless( len(sn) == 1, "wrong number of normals")
1168     self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1169     self.failUnless( numpy.linalg.norm(sn[0]-[0.70710678118654746, 0.70710678118654746, 0.]) < self.EPS, "wrong bottom vertex 1 ")
1170     dv=f.getDepthVectors(2)
1171     self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1172     self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1173     self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -30.]) < self.EPS, "wrong depth vector 0 ")
1174     self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1175     self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -30.]) < self.EPS, "wrong depth vector 1 ")
1176     b=f.getBottomPolyline(2)
1177     self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1178     self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1179     self.failUnless( numpy.linalg.norm(b[0]-[10., 0., -30.]) < self.EPS, "wrong bottom vertex 0 ")
1180     self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1181     self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -30.]) < self.EPS, "wrong bottom vertex 1 ")
1182     ds=f.getDepths(2)
1183     self.failUnless( len(ds) == 2, "wrong number of depth")
1184     self.failUnless( abs(ds[0]-30.) < self.EPS, "wrong depth at vertex 0 ")
1185     self.failUnless( abs(ds[1]-30.) < self.EPS, "wrong depth at vertex 1 ")
1186    
1187     top2=[ [10.,0.,0.], [0., 10., 0.] ]
1188     f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/4, depths=50.,tag=2)
1189     self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1190     self.failUnless( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1191     self.failUnless( 50. == f.getMediumDepth(2), "depth wrong")
1192     self.failUnless( (-50., 0.) == f.getW1Range(2)," wrong W1 range")
1193     segs=f.getTopPolyline(2)
1194     self.failUnless( len(segs) == 2, "wrong number of segments")
1195     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1196     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1197     self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1198     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1199     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1200     self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1201     d=f.getDips(2)
1202     self.failUnless( len(d) == 1, "wrong number of dips")
1203     self.failUnless( abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1204     sn=f.getSegmentNormals(2)
1205     self.failUnless( len(sn) == 1, "wrong number of normals")
1206     self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1207     self.failUnless( numpy.linalg.norm(sn[0]-[0.5,0.5,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1208     dv=f.getDepthVectors(2)
1209     self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1210     self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1211     self.failUnless( numpy.linalg.norm(dv[0]-[25., 25., -35.355339059327378]) < self.EPS, "wrong depth vector 0 ")
1212     self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1213     self.failUnless( numpy.linalg.norm(dv[1]-[25.,25., -35.355339059327378]) < self.EPS, "wrong depth vector 1 ")
1214     b=f.getBottomPolyline(2)
1215     self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1216     self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1217     self.failUnless( numpy.linalg.norm(b[0]-[35., 25., -35.355339059327378]) < self.EPS, "wrong bottom vertex 0 ")
1218     self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1219     self.failUnless( numpy.linalg.norm(b[1]-[25, 35., -35.355339059327378]) < self.EPS, "wrong bottom vertex 1 ")
1220     ds=f.getDepths(2)
1221     self.failUnless( len(ds) == 2, "wrong number of depth")
1222     self.failUnless( abs(ds[0]-50.) < self.EPS, "wrong depth at vertex 0 ")
1223     self.failUnless( abs(ds[1]-50.) < self.EPS, "wrong depth at vertex 1 ")
1224    
1225     top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
1226     f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=1, dips=pi/4, depths=20.)
1227     self.failUnless( 20. == f.getTotalLength(1), "length wrong")
1228     self.failUnless( 20. == f.getMediumDepth(1), "depth wrong")
1229     segs=f.getTopPolyline(1)
1230     self.failUnless( len(segs) == 3, "wrong number of segments")
1231     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1232     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1233     self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1234     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1235     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1236     self.failUnless( numpy.linalg.norm(segs[1]-[10.,10.,0.]) < self.EPS, "wrong vertex. 1 ")
1237     self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1238     self.failUnless( segs[2].size == 3, "seg 2 has wrong size.")
1239     self.failUnless( numpy.linalg.norm(segs[2]-[0.,10.,0.]) < self.EPS, "wrong vertex. 2 ")
1240     d=f.getDips(1)
1241     self.failUnless( len(d) == 2, "wrong number of dips")
1242     self.failUnless( abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1243     self.failUnless( abs(d[1]-0.78539816339744828) < self.EPS, "wrong dip 0")
1244     ds=f.getDepths(1)
1245     self.failUnless( len(ds) == 3, "wrong number of depth")
1246     self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1247     self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1248     sn=f.getSegmentNormals(1)
1249     self.failUnless( len(sn) == 2, "wrong number of normals")
1250     self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1251     self.failUnless( numpy.linalg.norm(sn[0]-[0.70710678118654746,0.,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1252     self.failUnless( isinstance(sn[1], numpy.ndarray), "wrong class of bottom vertex 0")
1253     self.failUnless( numpy.linalg.norm(sn[1]-[0.,0.70710678118654746,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1254     dv=f.getDepthVectors(1)
1255     self.failUnless( len(dv) == 3, "wrong number of depth vectors.")
1256     self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1257     self.failUnless( numpy.linalg.norm(dv[0]-[14.142135623730951, 0., -14.142135623730951]) < self.EPS, "wrong depth vector 0 ")
1258     self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1259     self.failUnless( numpy.linalg.norm(dv[1]-[11.547005383792515,11.547005383792515, -11.547005383792515]) < self.EPS, "wrong depth vector 2 ")
1260     self.failUnless( isinstance(dv[2], numpy.ndarray), "wrong class of depth vector 1")
1261     self.failUnless( numpy.linalg.norm(dv[2]-[0.,14.142135623730951, -14.142135623730951]) < self.EPS, "wrong depth vector 2 ")
1262     segs=f.getBottomPolyline(1)
1263     self.failUnless( len(segs) == 3, "wrong number of segments")
1264     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1265     self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1266     self.failUnless( numpy.linalg.norm(segs[0]-[24.142135623730951,0.,-14.142135623730951]) < self.EPS, "wrong vertex. 0 ")
1267     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1268     self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1269     self.failUnless( numpy.linalg.norm(segs[1]-[21.547005383792515,21.547005383792515, -11.547005383792515]) < self.EPS, "wrong vertex. 1 ")
1270     self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1271     self.failUnless( segs[2].size == 3, "seg 2 has wrong size.")
1272     self.failUnless( numpy.linalg.norm(segs[2]-[0., 24.142135623730951, -14.142135623730951]) < self.EPS, "wrong vertex. 2 ")
1273     self.failUnless( abs(0.-f.getW0Range(1)[0]) <=self.EPS," wrong W0 range (0)")
1274     self.failUnless( abs(31.857329272664341-f.getW0Range(1)[1]) <=self.EPS," wrong W0 range (1)")
1275     self.failUnless( abs(-20. - f.getW1Range(1)[0]) <=self.EPS," wrong W1 range (0)")
1276     self.failUnless( abs(0. - f.getW1Range(1)[1]) <=self.EPS," wrong W1 range (1)")
1277     self.failUnless( abs(0.0-f.getW0Offsets(1)[0])<=self.EPS," wrong W0 offsets (0)")
1278     self.failUnless( abs(15.92866463633217-f.getW0Offsets(1)[1])<=self.EPS," wrong W0 offsets (1)")
1279     self.failUnless( abs(31.857329272664341-f.getW0Offsets(1)[2])<=self.EPS," wrong W0 offsets(2)")
1280     #
1281     # ============ fresh start ====================
1282     #
1283     f.addFault(V0=[1.,0,0.],strikes=[pi/2, pi],ls=[1.,1.], dips=pi/2,depths=20,tag=1)
1284     f.addFault(V0=[10.,0,0],strikes=[3.*pi/4],ls=[14.142135623730951], tag=2, w0_offsets=[0,20], w1_max=20, dips=pi/2,depths=20)
1285     c=f.getCenterOnSurface()
1286     self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1287     self.failUnless( c.size == 3, "center size is wrong")
1288     self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.,0.]) < self.EPS, "center has wrong coordinates.")
1289     o=f.getOrientationOnSurface()/pi*180.
1290     self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1291    
1292     f.transform(rot=-pi/2., shift=[-1.,-1.,0.])
1293     self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
1294     self.failUnless( 2. == f.getTotalLength(1), "length after transformation wrong")
1295     self.failUnless( 20. == f.getMediumDepth(1), "depth after transformation wrong")
1296     rw0=f.getW0Range(1)
1297     self.failUnless( len(rw0) ==2, "wo range has wrong length")
1298     self.failUnless( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
1299     self.failUnless( abs(rw0[1]-2.) < self.EPS,"W0 1 wrong.")
1300     self.failUnless( (-20., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
1301     dips=f.getDips(1)
1302     self.failUnless(len(dips) == 2, "wrong number of dips.")
1303     self.failUnless( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
1304     self.failUnless( abs(dips[1]-1.5707963267948966) <= self.EPS, "wrong dip")
1305     ds=f.getDepths(1)
1306     self.failUnless( len(ds) == 3, "wrong number of depth")
1307     self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1308     self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1309     self.failUnless( abs(ds[2]-20.) < self.EPS, "wrong depth at vertex 1 ")
1310     segs=f.getTopPolyline(1)
1311     self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
1312     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1313     self.failUnless( segs[0].size == 3, "seg 0 has wrong size after transformation.")
1314     self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1315     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1316     self.failUnless( segs[1].size == 3, "seg 1 has wrong size after transformation.")
1317     self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1318     self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1319     self.failUnless( segs[2].size == 3, "seg 2 has wrong size after transformation.")
1320     self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.,0.]) < self.EPS, "wrong vertex after transformation. 2 ")
1321     self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1322     self.failUnless( 20. == f.getMediumDepth(2), "depth wrong after transformation")
1323     rw0=f.getW0Range(2)
1324     self.failUnless( len(rw0) ==2, "wo range has wrong length")
1325     self.failUnless( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
1326     self.failUnless( abs(rw0[1]-20.) < self.EPS,"W0 1 wrong.")
1327     self.failUnless( (-20., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1328     self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1329     dips=f.getDips(2)
1330     self.failUnless(len(dips) == 1, "wrong number of dips.")
1331     self.failUnless( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
1332     ds=f.getDepths(2)
1333     self.failUnless( len(ds) == 2, "wrong number of depth")
1334     self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1335     self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1336     segs=f.getTopPolyline(2)
1337     self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
1338     self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1339     self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1340     self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1341     self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.,0.]) < self.EPS, "wrong vertex. 1 after transformation")
1342     #
1343     # ============ fresh start ====================
1344     #
1345     f=FaultSystem(dim=3)
1346    
1347     top1=[ [0.,0.,0.], [1., 0., 0.] ]
1348     f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=1.,tag=1)
1349     top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
1350     f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=2, dips=pi/4, depths=20.)
1351    
1352     p,m=f.getParametrization([0.3,0.,-0.5],1)
1353     self.failUnless(length(p-[0.3,-0.5]) <= self.EPS, "wrong value.")
1354     self.failUnless(m==1., "wrong value.")
1355    
1356     p,m=f.getParametrization([0.5,0.,-0.5],1)
1357     self.failUnless(length(p-[0.5,-0.5]) <= self.EPS, "wrong value.")
1358     self.failUnless(m==1., "wrong value.")
1359    
1360     p,m=f.getParametrization([0.25,0.,-0.5],1)
1361     self.failUnless(length(p-[0.25,-0.5]) <= self.EPS, "wrong value.")
1362     self.failUnless(m==1., "wrong value.")
1363    
1364     p,m=f.getParametrization([0.5,0.,-0.25],1)
1365     self.failUnless(length(p-[0.5,-0.25]) <= self.EPS, "wrong value.")
1366     self.failUnless(m==1., "wrong value.")
1367    
1368     p,m=f.getParametrization([0.001,0.,-0.001],1)
1369     self.failUnless(length(p-[0.001, -0.001]) <= self.EPS, "wrong value.")
1370     self.failUnless(m==1., "wrong value.")
1371    
1372     p,m=f.getParametrization([0.001,0.,0.001],1)
1373     self.failUnless(m==0., "wrong value.")
1374    
1375     p,m=f.getParametrization([0.999,0.,0.001],1)
1376     self.failUnless(m==0., "wrong value.")
1377    
1378     p,m=f.getParametrization([1.001,0.,-0.001],1)
1379     self.failUnless(m==0., "wrong value.")
1380     p,m=f.getParametrization([1.001,0.,-0.1],1)
1381     self.failUnless(m==0., "wrong value.")
1382     p,m=f.getParametrization([1.001,0.,-0.000000001],1)
1383     self.failUnless(m==0., "wrong value.")
1384    
1385     p,m=f.getParametrization([0.999,0.,-0.001],1)
1386     self.failUnless(length(p-[0.999, -0.001]) <= self.EPS, "wrong value.")
1387     self.failUnless(m==1., "wrong value.")
1388    
1389     p,m=f.getParametrization([ 16.29252873 , 6.46410161 ,-6.29252873], 2, tol=1.e-7)
1390     self.failUnless(m==1., "wrong value.")
1391     self.failUnless(length(p-[4.7785993908996511, -10]) <= self.EPS*10., "wrong value.")
1392     p,m=f.getParametrization([15.77350269, 12.77350269, -5.77350269], 2, tol=1.e-7)
1393     self.failUnless(m==1., "wrong value.")
1394     self.failUnless(length(p-[11.150065245432518, -10]) <= self.EPS*10., "wrong value.")
1395    
1396     p,m=f.getParametrization([ 3., 17.0710678, -7.0710678], 2, tol=1.e-7)
1397     self.failUnless(m==1., "wrong value.")
1398     self.failUnless(length(p-[27.078729881764687, -10]) <= self.EPS*10., "wrong value.")
1399     p,m=f.getParametrization([9.30940108, 16.55204176, -6.55204176], 2, tol=1.e-7)
1400     self.failUnless(m==1., "wrong value.")
1401     self.failUnless(length(p-[20.707264027231822, -10]) <= self.EPS*10., "wrong value.")
1402    
1403     p,m=f.getParametrization([ 21.54700538, 21.54700538, -11.54700538], 2, tol=1.e-7)
1404     self.failUnless(m==1., "wrong value.")
1405     self.failUnless(length(p-[15.92866463633217, -20]) <= self.EPS*10., "wrong value.")
1406    
1407     p,m=f.getParametrization([ 0.,0.,0.], 2, tol=1.e-7)
1408     self.failUnless(m==0., "wrong value.")
1409    
1410     p,m=f.getParametrization([ 11.,11.,0.], 2, tol=1.e-7)
1411     self.failUnless(m==0., "wrong value.")
1412    
1413    
1414     s,d=f.getSideAndDistance([0.,-1.,0.], tag=1)
1415     self.failUnless( s>0, "wrong side.")
1416     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1417     s,d=f.getSideAndDistance([1.,-1.,0.], tag=1)
1418     self.failUnless( s>0, "wrong side.")
1419     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1420     s,d=f.getSideAndDistance([0.,1.,0.], tag=1)
1421     self.failUnless( s<0, "wrong side.")
1422     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1423     s,d=f.getSideAndDistance([1.,1.,0.], tag=1)
1424     self.failUnless( s<0, "wrong side.")
1425     self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1426    
1427    
1428     s,d=f.getSideAndDistance([0.,0.,0.], tag=2)
1429     self.failUnless( s<0, "wrong side.")
1430     self.failUnless( abs(d-10.)<self.EPS, "wrong distance.")
1431     s,d=f.getSideAndDistance([5.,5.,0.], tag=2)
1432     self.failUnless( s<0, "wrong side.")
1433     self.failUnless( abs(d-5.)<self.EPS, "wrong distance.")
1434    
1435     s,d=f.getSideAndDistance([10.,10.,-1.], tag=2)
1436     self.failUnless( s<0, "wrong side.")
1437     self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
1438     s,d=f.getSideAndDistance([10.,10.,-2.], tag=2)
1439     self.failUnless( s<0, "wrong side.")
1440     self.failUnless( abs(d-2.*0.70710678118654757)<self.EPS, "wrong distance.")
1441     s,d=f.getSideAndDistance([10.,10.,-3.], tag=2)
1442     self.failUnless( s<0, "wrong side.")
1443     self.failUnless( abs(d-3.*0.70710678118654757)<self.EPS, "wrong distance.")
1444    
1445     s,d=f.getSideAndDistance([5.,12.,0], tag=2)
1446     self.failUnless( s>0, "wrong side.")
1447     self.failUnless( abs(d-2*0.70710678118654757)<self.EPS, "wrong distance.")
1448     s,d=f.getSideAndDistance([5.,12.,-1], tag=2)
1449     self.failUnless( s>0, "wrong side.")
1450     self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
1451     s,d=f.getSideAndDistance([5.,12.,-2], tag=2)
1452     # s not checked as it is undefined.
1453     self.failUnless( abs(d)<self.EPS, "wrong distance.")
1454     s,d=f.getSideAndDistance([5.,12.,-3], tag=2)
1455     self.failUnless( s<0, "wrong side.")
1456     self.failUnless( abs(d-0.7071