/[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 3911 - (hide annotations)
Thu Jun 14 01:01:03 2012 UTC (7 years, 3 months ago) by jfenwick
File MIME type: text/x-python
File size: 73848 byte(s)
Copyright changes
1 gross 3071 # -*- coding: utf-8 -*-
2 ksteube 1809
3     ########################################################
4 gross 1414 #
5 jfenwick 3911 # Copyright (c) 2003-2012 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 3911 __copyright__="""Copyright (c) 2003-2012 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 jfenwick 3901 self.domain=Rectangle(NE,NE,order=-1,useElementsOnFace=0)
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 jfenwick 3551 self.assertTrue(error_p<10*self.TOL, "pressure error too large.")
76     self.assertTrue(error_v0<10*self.TOL, "0-velocity error too large.")
77     self.assertTrue(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 jfenwick 3551 self.assertTrue(error_v0<10*self.TOL, "0-velocity error too large.")
101     self.assertTrue(error_v1<10*self.TOL, "1-velocity error too large.")
102     self.assertTrue(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 jfenwick 3551 self.assertTrue(error_p<10*self.TOL, "pressure error too large.")
128     self.assertTrue(error_v0<10*self.TOL, "0-velocity error too large.")
129     self.assertTrue(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 jfenwick 3551 self.assertTrue(error_p<10*self.TOL, "pressure error too large.")
154     self.assertTrue(error_v0<10*self.TOL, "0-velocity error too large.")
155     self.assertTrue(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 jfenwick 3551 self.assertTrue(error_p<10*self.TOL, "pressure error too large.")
181     self.assertTrue(error_v0<10*self.TOL, "0-velocity error too large.")
182     self.assertTrue(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 jfenwick 3551 self.assertTrue(error_p<10*self.TOL, "pressure error too large.")
207     self.assertTrue(error_v0<10*self.TOL, "0-velocity error too large.")
208     self.assertTrue(error_v1<10*self.TOL, "1-velocity error too large.")
209 gross 2100 #====================================================================================================================
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 jfenwick 3901 self.domain=Brick(NE,NE,NE,order=-1,useElementsOnFace=0)
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 jfenwick 3551 self.assertTrue(error_p<10*self.TOL, "pressure error too large.")
245     self.assertTrue(error_v0<10*self.TOL, "0-velocity error too large.")
246     self.assertTrue(error_v1<10*self.TOL, "1-velocity error too large.")
247     self.assertTrue(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 jfenwick 3551 self.assertTrue(error_p<10*self.TOL, "pressure error too large.")
275     self.assertTrue(error_v0<10*self.TOL, "0-velocity error too large.")
276     self.assertTrue(error_v1<10*self.TOL, "1-velocity error too large.")
277     self.assertTrue(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 jfenwick 3551 self.assertTrue(error_p<10*self.TOL, "pressure error too large.")
306     self.assertTrue(error_v0<10*self.TOL, "0-velocity error too large.")
307     self.assertTrue(error_v1<10*self.TOL, "1-velocity error too large.")
308     self.assertTrue(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 jfenwick 3551 self.assertTrue(error_p<10*self.TOL, "pressure error too large.")
338     self.assertTrue(error_v0<10*self.TOL, "0-velocity error too large.")
339     self.assertTrue(error_v1<10*self.TOL, "1-velocity error too large.")
340     self.assertTrue(error_v2<10*self.TOL, "2-velocity error too large.")
341 gross 2100 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 jfenwick 3551 self.assertTrue(error_v0<10*self.TOL, "0-velocity error too large.")
368     self.assertTrue(error_v1<10*self.TOL, "1-velocity error too large.")
369     self.assertTrue(error_v2<10*self.TOL, "2-velocity error too large.")
370     self.assertTrue(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 jfenwick 3551 self.assertTrue(error_p<10*self.TOL, "pressure error too large.")
398     self.assertTrue(error_v0<10*self.TOL, "0-velocity error too large.")
399     self.assertTrue(error_v1<10*self.TOL, "1-velocity error too large.")
400     self.assertTrue(error_v2<10*self.TOL, "2-velocity error too large.")
401 gross 2100 #====================================================================================================================
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 jfenwick 3551 self.assertTrue(error_int<self.TOL, "Boundary intergral is too large.")
431 artak 2234
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 jfenwick 3551 self.assertTrue(error_int<self.TOL, "Boundary intergral is too large.")
458 artak 2234
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 jfenwick 3551 self.assertEqual(pl.getNumMaterials(),3,"num materials is wrong")
486     self.assertEqual(pl.validMaterialId(0),True,"material id 0 not found")
487     self.assertEqual(pl.validMaterialId(1),True,"material id 1 not found")
488     self.assertEqual(pl.validMaterialId(2),True,"material id 2 not found")
489     self.assertEqual(pl.validMaterialId(3),False,"material id 3 not found")
490 gross 2301
491 jfenwick 3551 self.assertEqual(pl.getFriction(),None,"initial friction wrong.")
492     self.assertEqual(pl.getTauY(),None,"initial tau y wrong.")
493 gross 2301 pl.setDruckerPragerLaw(tau_Y=10,friction=3)
494 jfenwick 3551 self.assertEqual(pl.getFriction(),3,"friction wrong.")
495     self.assertEqual(pl.getTauY(),10,"tau y wrong.")
496 gross 2301
497 jfenwick 3551 self.assertEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
498 gross 2301 pl.setElasticShearModulus(1000)
499 jfenwick 3551 self.assertEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
500 gross 2301
501     e=pl.getEtaN()
502 jfenwick 3551 self.assertEqual(len(e),3,"initial length of etaN is wrong.")
503     self.assertEqual(e,[None, None, None],"initial etaN are wrong.")
504     self.assertEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
505     self.assertEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
506     self.assertEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
507     self.assertRaises(ValueError, pl.getEtaN, 3)
508 gross 2301
509 jfenwick 3551 self.assertRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
510     self.assertRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
511     self.assertRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
512 gross 2301
513     pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
514 jfenwick 3551 self.assertEqual(pl.getPower(),[1,2,3],"powers are wrong.")
515     self.assertEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
516     self.assertEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
517 gross 2301
518     pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
519 jfenwick 3551 self.assertEqual(pl.getPower(),[4,2,3],"powers are wrong.")
520     self.assertEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
521     self.assertEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
522 gross 2301
523 jfenwick 3551 self.assertRaises(ValueError,pl.getPower,-1)
524     self.assertEqual(pl.getPower(0),4,"power 0 is wrong.")
525     self.assertEqual(pl.getPower(1),2,"power 1 is wrong.")
526     self.assertEqual(pl.getPower(2),3,"power 2 is wrong.")
527     self.assertRaises(ValueError,pl.getPower,3)
528 gross 2301
529 jfenwick 3551 self.assertRaises(ValueError,pl.getTauT,-1)
530     self.assertEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
531     self.assertEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
532     self.assertEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
533     self.assertRaises(ValueError,pl.getTauT,3)
534 gross 2301
535 jfenwick 3551 self.assertRaises(ValueError,pl.getEtaN,-1)
536     self.assertEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
537     self.assertEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
538     self.assertEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
539     self.assertRaises(ValueError,pl.getEtaN,3)
540 gross 2301
541     def checkResult(self,id,gamma_dot_, eta, tau_ref):
542 jfenwick 3551 self.assertTrue(eta>=0,"eta needs to be positive (test %s)"%id)
543 gross 2301 error=abs(gamma_dot_*eta-tau_ref)
544 jfenwick 3551 self.assertTrue(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 gross 2301
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 jfenwick 3772 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
554 gross 2301
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 jfenwick 3772 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
563 gross 2301
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 jfenwick 3772 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
572 gross 2301
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 jfenwick 3772 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
581 gross 2301
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 jfenwick 3772 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
590 gross 2301
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 jfenwick 3772 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
600 gross 2301
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 jfenwick 3772 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
610 gross 2301
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 jfenwick 3772 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
619 gross 2301
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 jfenwick 3772 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
628 gross 2301
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 jfenwick 3551 self.assertRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
639 jfenwick 3772 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
640 gross 2301
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 jfenwick 3772 if self.VERBOSE: print("dt =",dt)
684 gross 2432 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 jfenwick 3772 if self.VERBOSE: print("time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v)
734 jfenwick 3551 self.assertTrue( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
735     self.assertTrue( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )
736     self.assertTrue( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )
737     self.assertTrue( 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 jfenwick 3551 self.assertTrue( m == 0.25, "wrong max value")
839     self.assertTrue( t == 1, "wrong max tag")
840     self.assertTrue( 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 jfenwick 3551 self.assertTrue( m == 0.0625, "wrong max value")
847     self.assertTrue( t == 2, "wrong max tag")
848     self.assertTrue( 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 jfenwick 3551 self.assertTrue( m == 0.25, "wrong max value")
855     self.assertTrue( t == 2, "wrong max tag")
856     self.assertTrue( 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 jfenwick 3551 self.assertTrue( m == 0.25, "wrong max value")
863     self.assertTrue( t == 2, "wrong max tag")
864     self.assertTrue( 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 jfenwick 3551 self.assertTrue( m == 0.25, "wrong max value")
871     self.assertTrue( t == 1, "wrong max tag")
872     self.assertTrue( 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 jfenwick 3551 self.assertTrue( m == -0.25, "wrong min value")
885     self.assertTrue( t == 1, "wrong min tag")
886     self.assertTrue( 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 jfenwick 3551 self.assertTrue( m == -0.0625, "wrong min value")
892     self.assertTrue( t == 2, "wrong min tag")
893     self.assertTrue( 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 jfenwick 3551 self.assertTrue( m == -0.25, "wrong min value")
899     self.assertTrue( t == 2, "wrong min tag")
900     self.assertTrue( 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 jfenwick 3551 self.assertTrue( m == -0.25, "wrong min value")
906     self.assertTrue( t == 2, "wrong min tag")
907     self.assertTrue( 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 jfenwick 3551 self.assertTrue( m == -0.25, "wrong min value")
913     self.assertTrue( t == 1, "wrong min tag")
914     self.assertTrue( abs(l-0.70710678118654) <= self.EPS, "wrong min location")
915 gross 2675
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 jfenwick 3551 self.assertRaises(ValueError,f.addFault,V0=[1.,0],strikes=[pi/2, pi/2],ls=[1.,1.],tag=1,dips=top1)
921 gross 2663 f.addFault(V0=[1.,0],strikes=[pi/2, pi],ls=[1.,1.],tag=1)
922 jfenwick 3551 self.assertTrue(f.getDim() == 2, "wrong dimension")
923     self.assertTrue( [ 1 ] == f.getTags(), "tags wrong")
924     self.assertTrue( 2. == f.getTotalLength(1), "length wrong")
925     self.assertTrue( 0. == f.getMediumDepth(1), "depth wrong")
926     self.assertTrue( (0., 2.) == f.getW0Range(1)," wrong W0 range")
927     self.assertTrue( (0., 0.) == f.getW1Range(1)," wrong W1 range")
928     self.assertTrue( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsets")
929 gross 2663 segs=f.getTopPolyline(1)
930 jfenwick 3551 self.assertTrue( len(segs) == 3, "wrong number of segments")
931     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
932     self.assertTrue( segs[0].size == 2, "seg 0 has wrong size.")
933     self.assertTrue( numpy.linalg.norm(segs[0]-[1.,0.]) < self.EPS, "wrong vertex. 0 ")
934     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
935     self.assertTrue( segs[1].size == 2, "seg 1 has wrong size.")
936     self.assertTrue( numpy.linalg.norm(segs[1]-[1.,1.]) < self.EPS, "wrong vertex. 1 ")
937     self.assertTrue( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
938     self.assertTrue( segs[2].size == 2, "seg 2 has wrong size.")
939     self.assertTrue( numpy.linalg.norm(segs[2]-[0.,1.]) < self.EPS, "wrong vertex. 2 ")
940 gross 2647 c=f.getCenterOnSurface()
941 jfenwick 3551 self.assertTrue( isinstance(c, numpy.ndarray), "center has wrong class")
942     self.assertTrue( c.size == 2, "center size is wrong")
943     self.assertTrue( numpy.linalg.norm(c-[2./3.,2./3.]) < self.EPS, "center has wrong coordinates.")
944 gross 2647 o=f.getOrientationOnSurface()/pi*180.
945 jfenwick 3551 self.assertTrue( 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 jfenwick 3551 self.assertTrue( [ 1, 2 ] == f.getTags(), "tags wrong")
950     self.assertTrue( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length")
951     self.assertTrue( 0. == f.getMediumDepth(2), "depth wrong")
952     self.assertTrue( (0., 20.) == f.getW0Range(2)," wrong W0 range")
953     self.assertTrue( (0., 0.) == f.getW1Range(2)," wrong W1 range")
954     self.assertTrue( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets")
955 gross 2663 segs=f.getTopPolyline(2)
956 jfenwick 3551 self.assertTrue( len(segs) == 2, "wrong number of segments")
957     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
958     self.assertTrue( numpy.linalg.norm(segs[0]-[10.,0.]) < self.EPS, "wrong vertex. 0 ")
959     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
960     self.assertTrue( numpy.linalg.norm(segs[1]-[0.,10.]) < self.EPS, "wrong vertex. 1 ")
961 gross 2647 c=f.getCenterOnSurface()
962 jfenwick 3551 self.assertTrue( isinstance(c, numpy.ndarray), "center has wrong class")
963     self.assertTrue( c.size == 2, "center size is wrong")
964     self.assertTrue( numpy.linalg.norm(c-[12./5.,12./5.]) < self.EPS, "center has wrong coordinates.")
965 gross 2647 o=f.getOrientationOnSurface()/pi*180.
966 jfenwick 3551 self.assertTrue( abs(o+45.) < self.EPS, "wrong orientation.")
967 gross 2647
968 gross 2654 s,d=f.getSideAndDistance([0.,0.], tag=1)
969 jfenwick 3551 self.assertTrue( s<0, "wrong side.")
970     self.assertTrue( abs(d-1.)<self.EPS, "wrong distance.")
971 gross 2654 s,d=f.getSideAndDistance([0.,2.], tag=1)
972 jfenwick 3551 self.assertTrue( s>0, "wrong side.")
973     self.assertTrue( abs(d-1.)<self.EPS, "wrong distance.")
974 gross 2654 s,d=f.getSideAndDistance([1.,2.], tag=1)
975 jfenwick 3551 self.assertTrue( s>0, "wrong side.")
976     self.assertTrue( abs(d-1.)<self.EPS, "wrong distance.")
977 gross 2654 s,d=f.getSideAndDistance([2.,1.], tag=1)
978 jfenwick 3551 self.assertTrue( s>0, "wrong side.")
979     self.assertTrue( abs(d-1.)<self.EPS, "wrong distance.")
980 gross 2654 s,d=f.getSideAndDistance([2.,0.], tag=1)
981 jfenwick 3551 self.assertTrue( s>0, "wrong side.")
982     self.assertTrue( abs(d-1.)<self.EPS, "wrong distance.")
983 gross 2654 s,d=f.getSideAndDistance([0.,-1.], tag=1)
984 jfenwick 3551 self.assertTrue( s<0, "wrong side.")
985     self.assertTrue( abs(d-1.41421356237)<self.EPS, "wrong distance.")
986 gross 2654 s,d=f.getSideAndDistance([-1.,0], tag=1)
987 jfenwick 3551 self.assertTrue( s<0, "wrong side.")
988     self.assertTrue( abs(d-1.41421356237)<self.EPS, "wrong distance.")
989 gross 2654
990    
991 gross 2647 f.transform(rot=-pi/2., shift=[-1.,-1.])
992 jfenwick 3551 self.assertTrue( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
993     self.assertTrue( 2. == f.getTotalLength(1), "length after transformation wrong")
994     self.assertTrue( 0. == f.getMediumDepth(1), "depth after transformation wrong")
995     self.assertTrue( (0., 2.) == f.getW0Range(1)," wrong W0 after transformation range")
996     self.assertTrue( (0., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
997     self.assertTrue( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsetsafter transformation ")
998 gross 2663 segs=f.getTopPolyline(1)
999 jfenwick 3551 self.assertTrue( len(segs) == 3, "wrong number of segmentsafter transformation ")
1000     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1001     self.assertTrue( segs[0].size == 2, "seg 0 has wrong size after transformation.")
1002     self.assertTrue( numpy.linalg.norm(segs[0]-[-1.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1003     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1004     self.assertTrue( segs[1].size == 2, "seg 1 has wrong size after transformation.")
1005     self.assertTrue( numpy.linalg.norm(segs[1]-[0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1006     self.assertTrue( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1007     self.assertTrue( segs[2].size == 2, "seg 2 has wrong size after transformation.")
1008     self.assertTrue( numpy.linalg.norm(segs[2]-[0., 1.]) < self.EPS, "wrong vertex after transformation. 2 ")
1009     self.assertTrue( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1010     self.assertTrue( 0. == f.getMediumDepth(2), "depth wrong after transformation")
1011     self.assertTrue( (0., 20.) == f.getW0Range(2)," wrong W0 range after transformation")
1012     self.assertTrue( (0., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1013     self.assertTrue( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1014 gross 2663 segs=f.getTopPolyline(2)
1015 jfenwick 3551 self.assertTrue( len(segs) == 2, "wrong number of segments after transformation")
1016     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1017     self.assertTrue( numpy.linalg.norm(segs[0]-[-1.,-9]) < self.EPS, "wrong vertex. 0 after transformation")
1018     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1019     self.assertTrue( numpy.linalg.norm(segs[1]-[9.,1.]) < self.EPS, "wrong vertex. 1 after transformation")
1020 gross 2647
1021     c=f.getCenterOnSurface()
1022 jfenwick 3551 self.assertTrue( isinstance(c, numpy.ndarray), "center has wrong class")
1023     self.assertTrue( c.size == 2, "center size is wrong")
1024     self.assertTrue( numpy.linalg.norm(c-[7./5.,-7./5.]) < self.EPS, "center has wrong coordinates.")
1025 gross 2647 o=f.getOrientationOnSurface()/pi*180.
1026 jfenwick 3551 self.assertTrue( abs(o-45.) < self.EPS, "wrong orientation.")
1027 gross 2647
1028     p=f.getParametrization([-1.,0.],1)
1029 jfenwick 3551 self.assertTrue(p[1]==1., "wrong value.")
1030     self.assertTrue(abs(p[0])<self.EPS, "wrong value.")
1031 gross 2647 p=f.getParametrization([-0.5,0.],1)
1032 jfenwick 3551 self.assertTrue(p[1]==1., "wrong value.")
1033     self.assertTrue(abs(p[0]-0.5)<self.EPS* 0.5, "wrong value.")
1034 gross 2647 p=f.getParametrization([0.,0.],1)
1035 jfenwick 3551 self.assertTrue(p[1]==1., "wrong value.")
1036     self.assertTrue(abs(p[0]-1.)<self.EPS, "wrong value.")
1037 gross 2647 p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-8)
1038 jfenwick 3551 self.assertTrue(p[1]==0., "wrong value.")
1039 gross 2647 p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-6)
1040 jfenwick 3551 self.assertTrue(p[1]==1., "wrong value.")
1041     self.assertTrue(abs(p[0]-1.0000001)<self.EPS, "wrong value.")
1042 gross 2647 p=f.getParametrization([0.,0.5],1)
1043 jfenwick 3551 self.assertTrue(p[1]==1., "wrong value.")
1044     self.assertTrue(abs(p[0]-1.5)<self.EPS, "wrong value.")
1045 gross 2647 p=f.getParametrization([0,1.],1)
1046 jfenwick 3551 self.assertTrue(p[1]==1., "wrong value.")
1047     self.assertTrue(abs(p[0]-2.)<self.EPS, "wrong value.")
1048 gross 2647 p=f.getParametrization([1.,1.],1)
1049 jfenwick 3551 self.assertTrue(p[1]==0., "wrong value.")
1050 gross 2647 p=f.getParametrization([0,1.11],1)
1051 jfenwick 3551 self.assertTrue(p[1]==0., "wrong value.")
1052 gross 2647 p=f.getParametrization([-1,-9.],2)
1053 jfenwick 3551 self.assertTrue(p[1]==1., "wrong value.")
1054     self.assertTrue(abs(p[0])<self.EPS, "wrong value.")
1055 gross 2647 p=f.getParametrization([9,1],2)
1056 jfenwick 3551 self.assertTrue(p[1]==1., "wrong value.")
1057     self.assertTrue(abs(p[0]-20.)<self.EPS, "wrong value.")
1058 gross 2647
1059 gross 2663 def test_Fault3D(self):
1060     f=FaultSystem(dim=3)
1061 jfenwick 3551 self.assertTrue(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 jfenwick 3551 self.assertTrue( [ 1 ] == f.getTags(), "tags wrong")
1066     self.assertTrue( 1. == f.getTotalLength(1), "length wrong")
1067     self.assertTrue( 20. == f.getMediumDepth(1), "depth wrong")
1068     self.assertTrue( (0., 1.) == f.getW0Range(1)," wrong W0 range")
1069     self.assertTrue( (-20., 0.) == f.getW1Range(1)," wrong W1 range")
1070     self.assertTrue( [0., 1.] == f.getW0Offsets(1)," wrong W0 offsets")
1071 gross 2663 segs=f.getTopPolyline(1)
1072 jfenwick 3551 self.assertTrue( len(segs) == 2, "wrong number of segments")
1073     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1074     self.assertTrue( segs[0].size == 3, "seg 0 has wrong size.")
1075     self.assertTrue( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1076     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1077     self.assertTrue( segs[1].size == 3, "seg 1 has wrong size.")
1078     self.assertTrue( numpy.linalg.norm(segs[1]-[1.,0.,0]) < self.EPS, "wrong vertex. 1 ")
1079 gross 2663 c=f.getCenterOnSurface()
1080 jfenwick 3551 self.assertTrue( isinstance(c, numpy.ndarray), "center has wrong class")
1081     self.assertTrue( c.size == 3, "center size is wrong")
1082     self.assertTrue( numpy.linalg.norm(c-[0.5,0.,0.]) < self.EPS, "center has wrong coordinates.")
1083 gross 2663 o=f.getOrientationOnSurface()/pi*180.
1084 jfenwick 3551 self.assertTrue( abs(o) < self.EPS, "wrong orientation.")
1085 gross 2663 d=f.getDips(1)
1086 jfenwick 3551 self.assertTrue( len(d) == 1, "wrong number of dips")
1087     self.assertTrue( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1088 gross 2663 sn=f.getSegmentNormals(1)
1089 jfenwick 3551 self.assertTrue( len(sn) == 1, "wrong number of normals")
1090     self.assertTrue( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1091     self.assertTrue( numpy.linalg.norm(sn[0]-[0, -1., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1092 gross 2663 dv=f.getDepthVectors(1)
1093 jfenwick 3551 self.assertTrue( len(dv) == 2, "wrong number of depth vectors.")
1094     self.assertTrue( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1095     self.assertTrue( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1096     self.assertTrue( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1097     self.assertTrue( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1098 gross 2663 b=f.getBottomPolyline(1)
1099 jfenwick 3551 self.assertTrue( len(b) == 2, "wrong number of bottom vertices")
1100     self.assertTrue( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1101     self.assertTrue( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1102     self.assertTrue( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1103     self.assertTrue( numpy.linalg.norm(b[1]-[1., 0., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1104 gross 2663 ds=f.getDepths(1)
1105 jfenwick 3551 self.assertTrue( len(ds) == 2, "wrong number of depth")
1106     self.assertTrue( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1107     self.assertTrue( 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 jfenwick 3551 self.assertTrue( [ 1, 2 ] == f.getTags(), "tags wrong")
1112     self.assertTrue( 10. == f.getTotalLength(2), "length wrong")
1113     self.assertTrue( 20. == f.getMediumDepth(2), "depth wrong")
1114     self.assertTrue( (0., 10.) == f.getW0Range(2)," wrong W0 range")
1115     self.assertTrue( (-20., 0.) == f.getW1Range(2)," wrong W1 range")
1116     self.assertTrue( [0., 10.] == f.getW0Offsets(2)," wrong W0 offsets")
1117 gross 2663 segs=f.getTopPolyline(2)
1118 jfenwick 3551 self.assertTrue( len(segs) == 2, "wrong number of segments")
1119     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1120     self.assertTrue( segs[0].size == 3, "seg 0 has wrong size.")
1121     self.assertTrue( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1122     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1123     self.assertTrue( segs[1].size == 3, "seg 1 has wrong size.")
1124     self.assertTrue( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1125 gross 2663 d=f.getDips(2)
1126 jfenwick 3551 self.assertTrue( len(d) == 1, "wrong number of dips")
1127     self.assertTrue( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1128 gross 2663 sn=f.getSegmentNormals(2)
1129 jfenwick 3551 self.assertTrue( len(sn) == 1, "wrong number of normals")
1130     self.assertTrue( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1131     self.assertTrue( numpy.linalg.norm(sn[0]-[1, 0., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1132 gross 2663 dv=f.getDepthVectors(2)
1133 jfenwick 3551 self.assertTrue( len(dv) == 2, "wrong number of depth vectors.")
1134     self.assertTrue( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1135     self.assertTrue( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1136     self.assertTrue( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1137     self.assertTrue( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1138 gross 2663 b=f.getBottomPolyline(2)
1139 jfenwick 3551 self.assertTrue( len(b) == 2, "wrong number of bottom vertices")
1140     self.assertTrue( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1141     self.assertTrue( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1142     self.assertTrue( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1143     self.assertTrue( numpy.linalg.norm(b[1]-[0., 10., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1144 gross 2663 ds=f.getDepths(2)
1145 jfenwick 3551 self.assertTrue( len(ds) == 2, "wrong number of depth")
1146     self.assertTrue( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1147     self.assertTrue( 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 jfenwick 3551 self.assertTrue( [ 1, 2 ] == f.getTags(), "tags wrong")
1152     self.assertTrue( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1153     self.assertTrue( 30. == f.getMediumDepth(2), "depth wrong")
1154     self.assertTrue( (-30., 0.) == f.getW1Range(2)," wrong W1 range")
1155 gross 2663 segs=f.getTopPolyline(2)
1156 jfenwick 3551 self.assertTrue( len(segs) == 2, "wrong number of segments")
1157     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1158     self.assertTrue( segs[0].size == 3, "seg 0 has wrong size.")
1159     self.assertTrue( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1160     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1161     self.assertTrue( segs[1].size == 3, "seg 1 has wrong size.")
1162     self.assertTrue( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1163 gross 2663 d=f.getDips(2)
1164 jfenwick 3551 self.assertTrue( len(d) == 1, "wrong number of dips")
1165     self.assertTrue( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1166 gross 2663 sn=f.getSegmentNormals(2)
1167 jfenwick 3551 self.assertTrue( len(sn) == 1, "wrong number of normals")
1168     self.assertTrue( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1169     self.assertTrue( numpy.linalg.norm(sn[0]-[0.70710678118654746, 0.70710678118654746, 0.]) < self.EPS, "wrong bottom vertex 1 ")
1170 gross 2663 dv=f.getDepthVectors(2)
1171 jfenwick 3551 self.assertTrue( len(dv) == 2, "wrong number of depth vectors.")
1172     self.assertTrue( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1173     self.assertTrue( numpy.linalg.norm(dv[0]-[0., 0., -30.]) < self.EPS, "wrong depth vector 0 ")
1174     self.assertTrue( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1175     self.assertTrue( numpy.linalg.norm(dv[1]-[0., 0., -30.]) < self.EPS, "wrong depth vector 1 ")
1176 gross 2663 b=f.getBottomPolyline(2)
1177 jfenwick 3551 self.assertTrue( len(b) == 2, "wrong number of bottom vertices")
1178     self.assertTrue( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1179     self.assertTrue( numpy.linalg.norm(b[0]-[10., 0., -30.]) < self.EPS, "wrong bottom vertex 0 ")
1180     self.assertTrue( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1181     self.assertTrue( numpy.linalg.norm(b[1]-[0., 10., -30.]) < self.EPS, "wrong bottom vertex 1 ")
1182 gross 2663 ds=f.getDepths(2)
1183 jfenwick 3551 self.assertTrue( len(ds) == 2, "wrong number of depth")
1184     self.assertTrue( abs(ds[0]-30.) < self.EPS, "wrong depth at vertex 0 ")
1185     self.assertTrue( abs(ds[1]-30.) < self.EPS, "wrong depth at vertex 1 ")
1186 gross 2663
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 jfenwick 3551 self.assertTrue( [ 1, 2 ] == f.getTags(), "tags wrong")
1190     self.assertTrue( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1191     self.assertTrue( 50. == f.getMediumDepth(2), "depth wrong")
1192     self.assertTrue( (-50., 0.) == f.getW1Range(2)," wrong W1 range")
1193 gross 2663 segs=f.getTopPolyline(2)
1194 jfenwick 3551 self.assertTrue( len(segs) == 2, "wrong number of segments")
1195     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1196     self.assertTrue( segs[0].size == 3, "seg 0 has wrong size.")
1197     self.assertTrue( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1198     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1199     self.assertTrue( segs[1].size == 3, "seg 1 has wrong size.")
1200     self.assertTrue( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1201 gross 2663 d=f.getDips(2)
1202 jfenwick 3551 self.assertTrue( len(d) == 1, "wrong number of dips")
1203     self.assertTrue( abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1204 gross 2663 sn=f.getSegmentNormals(2)
1205 jfenwick 3551 self.assertTrue( len(sn) == 1, "wrong number of normals")
1206     self.assertTrue( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1207     self.assertTrue( numpy.linalg.norm(sn[0]-[0.5,0.5,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1208 gross 2663 dv=f.getDepthVectors(2)
1209 jfenwick 3551 self.assertTrue( len(dv) == 2, "wrong number of depth vectors.")
1210     self.assertTrue( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1211     self.assertTrue( numpy.linalg.norm(dv[0]-[25., 25., -35.355339059327378]) < self.EPS, "wrong depth vector 0 ")
1212     self.assertTrue( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1213     self.assertTrue( numpy.linalg.norm(dv[1]-[25.,25., -35.355339059327378]) < self.EPS, "wrong depth vector 1 ")
1214 gross 2663 b=f.getBottomPolyline(2)
1215 jfenwick 3551 self.assertTrue( len(b) == 2, "wrong number of bottom vertices")
1216     self.assertTrue( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1217     self.assertTrue( numpy.linalg.norm(b[0]-[35., 25., -35.355339059327378]) < self.EPS, "wrong bottom vertex 0 ")
1218     self.assertTrue( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1219     self.assertTrue( numpy.linalg.norm(b[1]-[25, 35., -35.355339059327378]) < self.EPS, "wrong bottom vertex 1 ")
1220 gross 2663 ds=f.getDepths(2)
1221 jfenwick 3551 self.assertTrue( len(ds) == 2, "wrong number of depth")
1222     self.assertTrue( abs(ds[0]-50.) < self.EPS, "wrong depth at vertex 0 ")
1223     self.assertTrue( abs(ds[1]-50.) < self.EPS, "wrong depth at vertex 1 ")
1224 gross 2663
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 jfenwick 3551 self.assertTrue( 20. == f.getTotalLength(1), "length wrong")
1228     self.assertTrue( 20. == f.getMediumDepth(1), "depth wrong")
1229 gross 2663 segs=f.getTopPolyline(1)
1230 jfenwick 3551 self.assertTrue( len(segs) == 3, "wrong number of segments")
1231     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1232     self.assertTrue( segs[0].size == 3, "seg 0 has wrong size.")
1233     self.assertTrue( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1234     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1235     self.assertTrue( segs[1].size == 3, "seg 1 has wrong size.")
1236     self.assertTrue( numpy.linalg.norm(segs[1]-[10.,10.,0.]) < self.EPS, "wrong vertex. 1 ")
1237     self.assertTrue( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1238     self.assertTrue( segs[2].size == 3, "seg 2 has wrong size.")
1239     self.assertTrue( numpy.linalg.norm(segs[2]-[0.,10.,0.]) < self.EPS, "wrong vertex. 2 ")
1240 gross 2663 d=f.getDips(1)
1241 jfenwick 3551 self.assertTrue( len(d) == 2, "wrong number of dips")
1242     self.assertTrue( abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1243     self.assertTrue( abs(d[1]-0.78539816339744828) < self.EPS, "wrong dip 0")
1244 gross 2663 ds=f.getDepths(1)
1245 jfenwick 3551 self.assertTrue( len(ds) == 3, "wrong number of depth")
1246     self.assertTrue( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1247     self.assertTrue( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1248 gross 2663 sn=f.getSegmentNormals(1)
1249 jfenwick 3551 self.assertTrue( len(sn) == 2, "wrong number of normals")
1250     self.assertTrue( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1251     self.assertTrue( numpy.linalg.norm(sn[0]-[0.70710678118654746,0.,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1252     self.assertTrue( isinstance(sn[1], numpy.ndarray), "wrong class of bottom vertex 0")
1253     self.assertTrue( numpy.linalg.norm(sn[1]-[0.,0.70710678118654746,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1254 gross 2663 dv=f.getDepthVectors(1)
1255 jfenwick 3551 self.assertTrue( len(dv) == 3, "wrong number of depth vectors.")
1256     self.assertTrue( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1257     self.assertTrue( numpy.linalg.norm(dv[0]-[14.142135623730951, 0., -14.142135623730951]) < self.EPS, "wrong depth vector 0 ")
1258     self.assertTrue( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1259     self.assertTrue( numpy.linalg.norm(dv[1]-[11.547005383792515,11.547005383792515, -11.547005383792515]) < self.EPS, "wrong depth vector 2 ")
1260     self.assertTrue( isinstance(dv[2], numpy.ndarray), "wrong class of depth vector 1")
1261     self.assertTrue( numpy.linalg.norm(dv[2]-[0.,14.142135623730951, -14.142135623730951]) < self.EPS, "wrong depth vector 2 ")
1262 gross 2663 segs=f.getBottomPolyline(1)
1263 jfenwick 3551 self.assertTrue( len(segs) == 3, "wrong number of segments")
1264     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1265     self.assertTrue( segs[0].size == 3, "seg 0 has wrong size.")
1266     self.assertTrue( numpy.linalg.norm(segs[0]-[24.142135623730951,0.,-14.142135623730951]) < self.EPS, "wrong vertex. 0 ")
1267     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1268     self.assertTrue( segs[1].size == 3, "seg 1 has wrong size.")
1269     self.assertTrue( numpy.linalg.norm(segs[1]-[21.547005383792515,21.547005383792515, -11.547005383792515]) < self.EPS, "wrong vertex. 1 ")
1270     self.assertTrue( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1271     self.assertTrue( segs[2].size == 3, "seg 2 has wrong size.")
1272     self.assertTrue( numpy.linalg.norm(segs[2]-[0., 24.142135623730951, -14.142135623730951]) < self.EPS, "wrong vertex. 2 ")
1273     self.assertTrue( abs(0.-f.getW0Range(1)[0]) <=self.EPS," wrong W0 range (0)")
1274     self.assertTrue( abs(31.857329272664341-f.getW0Range(1)[1]) <=self.EPS," wrong W0 range (1)")
1275     self.assertTrue( abs(-20. - f.getW1Range(1)[0]) <=self.EPS," wrong W1 range (0)")
1276     self.assertTrue( abs(0. - f.getW1Range(1)[1]) <=self.EPS," wrong W1 range (1)")
1277     self.assertTrue( abs(0.0-f.getW0Offsets(1)[0])<=self.EPS," wrong W0 offsets (0)")
1278     self.assertTrue( abs(15.92866463633217-f.getW0Offsets(1)[1])<=self.EPS," wrong W0 offsets (1)")
1279     self.assertTrue( abs(31.857329272664341-f.getW0Offsets(1)[2])<=self.EPS," wrong W0 offsets(2)")
1280 gross 2663 #
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 jfenwick 3551 self.assertTrue( isinstance(c, numpy.ndarray), "center has wrong class")
1287     self.assertTrue( c.size == 3, "center size is wrong")
1288     self.assertTrue( numpy.linalg.norm(c-[12./5.,12./5.,0.]) < self.EPS, "center has wrong coordinates.")
1289 gross 2663 o=f.getOrientationOnSurface()/pi*180.
1290 jfenwick 3551 self.assertTrue( abs(o+45.) < self.EPS, "wrong orientation.")
1291 gross 2663
1292     f.transform(rot=-pi/2., shift=[-1.,-1.,0.])
1293 jfenwick 3551 self.assertTrue( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
1294     self.assertTrue( 2. == f.getTotalLength(1), "length after transformation wrong")
1295     self.assertTrue( 20. == f.getMediumDepth(1), "depth after transformation wrong")
1296 gross 2663 rw0=f.getW0Range(1)
1297 jfenwick 3551 self.assertTrue( len(rw0) ==2, "wo range has wrong length")
1298     self.assertTrue( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
1299     self.assertTrue( abs(rw0[1]-2.) < self.EPS,"W0 1 wrong.")
1300     self.assertTrue( (-20., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
1301 gross 2663 dips=f.getDips(1)
1302 jfenwick 3551 self.assertTrue(len(dips) == 2, "wrong number of dips.")
1303     self.assertTrue( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
1304     self.assertTrue( abs(dips[1]-1.5707963267948966) <= self.EPS, "wrong dip")
1305 gross 2663 ds=f.getDepths(1)
1306 jfenwick 3551 self.assertTrue( len(ds) == 3, "wrong number of depth")
1307     self.assertTrue( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1308     self.assertTrue( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1309     self.assertTrue( abs(ds[2]-20.) < self.EPS, "wrong depth at vertex 1 ")
1310 gross 2663 segs=f.getTopPolyline(1)
1311 jfenwick 3551 self.assertTrue( len(segs) == 3, "wrong number of segmentsafter transformation ")
1312     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1313     self.assertTrue( segs[0].size == 3, "seg 0 has wrong size after transformation.")
1314     self.assertTrue( numpy.linalg.norm(segs[0]-[-1.,0.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1315     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1316     self.assertTrue( segs[1].size == 3, "seg 1 has wrong size after transformation.")
1317     self.assertTrue( numpy.linalg.norm(segs[1]-[0.,0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1318     self.assertTrue( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1319     self.assertTrue( segs[2].size == 3, "seg 2 has wrong size after transformation.")
1320     self.assertTrue( numpy.linalg.norm(segs[2]-[0., 1.,0.]) < self.EPS, "wrong vertex after transformation. 2 ")
1321     self.assertTrue( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1322     self.assertTrue( 20. == f.getMediumDepth(2), "depth wrong after transformation")
1323 gross 2663 rw0=f.getW0Range(2)
1324 jfenwick 3551 self.assertTrue( len(rw0) ==2, "wo range has wrong length")
1325     self.assertTrue( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
1326     self.assertTrue( abs(rw0[1]-20.) < self.EPS,"W0 1 wrong.")
1327     self.assertTrue( (-20., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1328     self.assertTrue( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1329 gross 2663 dips=f.getDips(2)
1330 jfenwick 3551 self.assertTrue(len(dips) == 1, "wrong number of dips.")
1331     self.assertTrue( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
1332 gross 2663 ds=f.getDepths(2)
1333 jfenwick 3551 self.assertTrue( len(ds) == 2, "wrong number of depth")
1334     self.assertTrue( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1335     self.assertTrue( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1336 gross 2663 segs=f.getTopPolyline(2)
1337 jfenwick 3551 self.assertTrue( len(segs) == 2, "wrong number of segments after transformation")
1338     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1339     self.assertTrue( numpy.linalg.norm(segs[0]-[-1.,-9,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1340     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1341     self.assertTrue( numpy.linalg.norm(segs[1]-[9.,1.,0.]) < self.EPS, "wrong vertex. 1 after transformation")
1342 gross 2663 #
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 jfenwick 3551 self.assertTrue(length(p-[0.3,-0.5]) <= self.EPS, "wrong value.")
1354     self.assertTrue(m==1., "wrong value.")
1355 gross 2663
1356     p,m=f.getParametrization([0.5,0.,-0.5],1)
1357 jfenwick 3551 self.assertTrue(length(p-[0.5,-0.5]) <= self.EPS, "wrong value.")
1358     self.assertTrue(m==1., "wrong value.")
1359 gross 2663
1360     p,m=f.getParametrization([0.25,0.,-0.5],1)
1361 jfenwick 3551 self.assertTrue(length(p-[0.25,-0.5]) <= self.EPS, "wrong value.")
1362     self.assertTrue(m==1., "wrong value.")
1363 gross 2663
1364     p,m=f.getParametrization([0.5,0.,-0.25],1)
1365 jfenwick 3551 self.assertTrue(length(p-[0.5,-0.25]) <= self.EPS, "wrong value.")
1366     self.assertTrue(m==1., "wrong value.")
1367 gross 2663
1368     p,m=f.getParametrization([0.001,0.,-0.001],1)
1369 jfenwick 3551 self.assertTrue(length(p-[0.001, -0.001]) <= self.EPS, "wrong value.")
1370     self.assertTrue(m==1., "wrong value.")
1371 gross 2663
1372     p,m=f.getParametrization([0.001,0.,0.001],1)
1373 jfenwick 3551 self.assertTrue(m==0., "wrong value.")
1374 gross 2663
1375     p,m=f.getParametrization([0.999,0.,0.001],1)
1376 jfenwick 3551 self.assertTrue(m==0., "wrong value.")
1377 gross 2663
1378     p,m=f.getParametrization([1.001,0.,-0.001],1)
1379 jfenwick 3551 self.assertTrue(m==0., "wrong value.")
1380 gross 2663 p,m=f.getParametrization([1.001,0.,-0.1],1)
1381 jfenwick 3551 self.assertTrue(m==0., "wrong value.")
1382 gross 2663 p,m=f.getParametrization([1.001,0.,-0.000000001],1)
1383 jfenwick 3551 self.assertTrue(m==0., "wrong value.")
1384 gross 2663
1385     p,m=f.getParametrization([0.999,0.,-0.001],1)
1386 jfenwick 3551 self.assertTrue(length(p-[0.999, -0.001]) <= self.EPS, "wrong value.")
1387     self.assertTrue(m==1., "wrong value.")
1388 gross 2663
1389     p,m=f.getParametrization([ 16.29252873 , 6.46410161 ,-6.29252873], 2, tol=1.e-7)
1390 jfenwick 3551 self.assertTrue(m==1., "wrong value.")
1391     self.assertTrue(length(p-[4.7785993908996511, -10]) <= self.EPS*10., "wrong value.")
1392 gross 2663 p,m=f.getParametrization([15.77350269, 12.77350269, -5.77350269], 2, tol=1.e-7)
1393 jfenwick 3551 self.assertTrue(m==1., "wrong value.")
1394     self.assertTrue(length(p-[11.150065245432518, -10]) <= self.EPS*10., "wrong value.")
1395 gross 2663
1396     p,m=f.getParametrization([ 3., 17.0710678, -7.0710678], 2, tol=1.e-7)
1397 jfenwick 3551 self.assertTrue(m==1., "wrong value.")
1398     self.assertTrue(length(p-[27.078729881764687, -10]) <= self.EPS*10., "wrong value.")
1399 gross 2663 p,m=f.getParametrization([9.30940108, 16.55204176, -6.55204176], 2, tol=1.e-7)
1400 jfenwick 3551 self.assertTrue(m==1., "wrong value.")
1401     self.assertTrue(length(p-[20.707264027231822, -10]) <= self.EPS*10., "wrong value.")
1402 gross 2663
1403     p,m=f.getParametrization([ 21.54700538, 21.54700538, -11.54700538], 2, tol=1.e-7)
1404 jfenwick 3551 self.assertTrue(m==1., "wrong value.")
1405     self.assertTrue(length(p-[15.92866463633217, -20]) <= self.EPS*10., "wrong value.")
1406 gross 2663
1407     p,m=f.getParametrization([ 0.,0.,0.], 2, tol=1.e-7)
1408 jfenwick 3551 self.assertTrue(m==0., "wrong value.")
1409 gross 2663
1410     p,m=f.getParametrization([ 11.,11.,0.], 2, tol=1.e-7)
1411 jfenwick 3551 self.assertTrue(m==0., "wrong value.")
1412 gross 2663
1413    
1414     s,d=f.getSideAndDistance([0.,-1.,0.], tag=1)
1415 jfenwick 3551 self.assertTrue( s>0, "wrong side.")
1416     self.assertTrue( abs(d-1.)<self.EPS, "wrong distance.")