/[escript]/trunk/dudley/test/python/run_models.py
ViewVC logotype

Annotation of /trunk/dudley/test/python/run_models.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3774 - (hide annotations)
Wed Jan 18 06:29:34 2012 UTC (7 years, 9 months ago) by jfenwick
File MIME type: text/x-python
File size: 77744 byte(s)
dudley, pasowrap, pycad

1 gross 3071 # -*- coding: utf-8 -*-
2 ksteube 1809
3     ########################################################
4 gross 1414 #
5 jfenwick 2881 # Copyright (c) 2003-2010 by University of Queensland
6 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
7     # http://www.uq.edu.au/esscc
8 gross 1414 #
9 ksteube 1809 # Primary Business: Queensland, Australia
10     # Licensed under the Open Software License version 3.0
11     # http://www.opensource.org/licenses/osl-3.0.php
12 gross 1414 #
13 ksteube 1809 ########################################################
14 gross 1414
15 jfenwick 2881 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
16 ksteube 1809 Earth Systems Science Computational Center (ESSCC)
17     http://www.uq.edu.au/esscc
18     Primary Business: Queensland, Australia"""
19 gross 1414 __license__="""Licensed under the Open Software License version 3.0
20 ksteube 1809 http://www.opensource.org/licenses/osl-3.0.php"""
21 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
22 ksteube 1809
23 gross 1414 import unittest
24     import tempfile
25    
26    
27    
28 gross 3071 VERBOSE=False and True
29 gross 1414
30     from esys.escript import *
31 gross 2647 from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian, FaultSystem, DarcyFlow
32     from esys.escript.models import Mountains
33 jfenwick 3085 from esys.dudley 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 jfenwick 3085 DUDLEY_WORKDIR=os.environ['DUDLEY_WORKDIR']
42 gross 2647 except KeyError:
43 jfenwick 3085 DUDLEY_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 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 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 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     class Test_Darcy(unittest.TestCase):
403     # this is a simple test for the darcy flux problem
404     #
405     #
406     # p = 1/k * ( 1/2* (fb-f0)/xb* x **2 + f0 * x - ub*x ) + p0
407     #
408     # with f = (fb-f0)/xb* x + f0
409     #
410     # u = f - k * p,x = ub
411     #
412     # we prescribe pressure at x=x0=0 to p0
413     #
414     # if we prescribe the pressure on the bottom x=xb we set
415     #
416     # pb= 1/k * ( 1/2* (fb-f0)* xb + f0 * xb - ub*xb ) + p0 = 1/k * ((fb+f0)/2 - ub ) * xb + p0
417     #
418     # which leads to ub = (fb+f0)/2-k*(pb-p0)/xb
419     #
420     def rescaleDomain(self):
421     x=self.dom.getX().copy()
422 jfenwick 3774 for i in range(self.dom.getDim()):
423 gross 2100 x_inf=inf(x[i])
424     x_sup=sup(x[i])
425     if i == self.dom.getDim()-1:
426     x[i]=-self.WIDTH*(x[i]-x_sup)/(x_inf-x_sup)
427     else:
428     x[i]=self.WIDTH*(x[i]-x_inf)/(x_sup-x_inf)
429     self.dom.setX(x)
430     def getScalarMask(self,include_bottom=True):
431     x=self.dom.getX().copy()
432     x_inf=inf(x[self.dom.getDim()-1])
433     x_sup=sup(x[self.dom.getDim()-1])
434     out=whereZero(x[self.dom.getDim()-1]-x_sup)
435     if include_bottom: out+=whereZero(x[self.dom.getDim()-1]-x_inf)
436     return wherePositive(out)
437     def getVectorMask(self,include_bottom=True):
438     x=self.dom.getX().copy()
439     out=Vector(0.,Solution(self.dom))
440 jfenwick 3774 for i in range(self.dom.getDim()):
441 gross 2100 x_inf=inf(x[i])
442     x_sup=sup(x[i])
443     if i != self.dom.getDim()-1: out[i]+=whereZero(x[i]-x_sup)
444     if i != self.dom.getDim()-1 or include_bottom: out[i]+=whereZero(x[i]-x_inf)
445     return wherePositive(out)
446 gross 1471
447 gross 2100 def setSolutionFixedBottom(self, p0, pb, f0, fb, k):
448     d=self.dom.getDim()
449     x=self.dom.getX()[d-1]
450     xb=inf(x)
451     u=Vector(0.,Solution(self.dom))+kronecker(d)[d-1]*((f0+fb)/2.-k*(pb-p0)/xb)
452     p=1./k*((fb-f0)/(xb*2.)* x**2 - (fb-f0)/2.*x)+(pb-p0)/xb*x + p0
453     f= ((fb-f0)/xb* x + f0)*kronecker(Function(self.dom))[d-1]
454     return u,p,f
455    
456     def testConstF_FixedBottom_smallK(self):
457     k=1.e-10
458     mp=self.getScalarMask(include_bottom=True)
459     mv=self.getVectorMask(include_bottom=False)
460     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
461     p=p_ref*mp
462     u=u_ref*mv
463     df=DarcyFlow(self.dom)
464     df.setValue(g=f,
465     location_of_fixed_pressure=mp,
466     location_of_fixed_flux=mv,
467     permeability=Scalar(k,Function(self.dom)))
468 gross 2264 df.setTolerance(rtol=self.TOL)
469     v,p=df.solve(u_ref,p, max_iter=100, verbose=VERBOSE)
470 jfenwick 3551 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
471     self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
472 gross 2100 def testConstF_FixedBottom_mediumK(self):
473     k=1.
474     mp=self.getScalarMask(include_bottom=True)
475     mv=self.getVectorMask(include_bottom=False)
476     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
477     p=p_ref*mp
478     u=u_ref*mv
479     df=DarcyFlow(self.dom)
480     df.setValue(g=f,
481     location_of_fixed_pressure=mp,
482     location_of_fixed_flux=mv,
483     permeability=Scalar(k,Function(self.dom)))
484 gross 2264 df.setTolerance(rtol=self.TOL)
485     v,p=df.solve(u,p,max_iter=100, verbose=VERBOSE )
486 jfenwick 3551 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
487     self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
488 artak 1518
489 gross 2100 def testConstF_FixedBottom_largeK(self):
490     k=1.e10
491     mp=self.getScalarMask(include_bottom=True)
492     mv=self.getVectorMask(include_bottom=False)
493     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
494     p=p_ref*mp
495     u=u_ref*mv
496     df=DarcyFlow(self.dom)
497     df.setValue(g=f,
498     location_of_fixed_pressure=mp,
499     location_of_fixed_flux=mv,
500     permeability=Scalar(k,Function(self.dom)))
501 gross 2264 df.setTolerance(rtol=self.TOL)
502     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
503 jfenwick 3551 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
504     self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
505 artak 1518
506 gross 2100 def testVarioF_FixedBottom_smallK(self):
507     k=1.e-10
508     mp=self.getScalarMask(include_bottom=True)
509     mv=self.getVectorMask(include_bottom=False)
510     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
511     p=p_ref*mp
512     u=u_ref*mv
513     df=DarcyFlow(self.dom)
514 gross 3071 #df.getSolverOptionsPressure().setVerbosityOn()
515 gross 2100 df.setValue(g=f,
516     location_of_fixed_pressure=mp,
517     location_of_fixed_flux=mv,
518     permeability=Scalar(k,Function(self.dom)))
519 gross 2264 df.setTolerance(rtol=self.TOL)
520     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
521 gross 3071
522 jfenwick 3551 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
523     self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
524 artak 1518
525 gross 2100 def testVarioF_FixedBottom_mediumK(self):
526     k=1.
527     mp=self.getScalarMask(include_bottom=True)
528     mv=self.getVectorMask(include_bottom=False)
529     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
530     p=p_ref*mp
531     u=u_ref*mv
532     df=DarcyFlow(self.dom)
533     df.setValue(g=f,
534     location_of_fixed_pressure=mp,
535     location_of_fixed_flux=mv,
536     permeability=Scalar(k,Function(self.dom)))
537 gross 2264 df.setTolerance(rtol=self.TOL)
538     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
539 jfenwick 3551 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
540     self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
541 artak 1518
542 gross 2100 def testVarioF_FixedBottom_largeK(self):
543     k=1.e10
544     mp=self.getScalarMask(include_bottom=True)
545     mv=self.getVectorMask(include_bottom=False)
546     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
547     p=p_ref*mp
548     u=u_ref*mv
549     df=DarcyFlow(self.dom)
550     df.setValue(g=f,
551     location_of_fixed_pressure=mp,
552     location_of_fixed_flux=mv,
553     permeability=Scalar(k,Function(self.dom)))
554 gross 2264 df.setTolerance(rtol=self.TOL)
555     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
556 jfenwick 3551 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
557     self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
558 artak 1518
559 gross 2100 def testConstF_FreeBottom_smallK(self):
560     k=1.e-10
561     mp=self.getScalarMask(include_bottom=False)
562     mv=self.getVectorMask(include_bottom=True)
563     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
564     p=p_ref*mp
565     u=u_ref*mv
566     df=DarcyFlow(self.dom)
567     df.setValue(g=f,
568 gross 2208 location_of_fixed_pressure=mp,
569 gross 2100 location_of_fixed_flux=mv,
570     permeability=Scalar(k,Function(self.dom)))
571 gross 2264 df.setTolerance(rtol=self.TOL)
572     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
573 gross 3071
574    
575 jfenwick 3551 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
576     self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
577 artak 1518
578 gross 2100 def testConstF_FreeBottom_mediumK(self):
579     k=1.
580     mp=self.getScalarMask(include_bottom=False)
581     mv=self.getVectorMask(include_bottom=True)
582     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
583     p=p_ref*mp
584     u=u_ref*mv
585     df=DarcyFlow(self.dom)
586     df.setValue(g=f,
587     location_of_fixed_pressure=mp,
588     location_of_fixed_flux=mv,
589     permeability=Scalar(k,Function(self.dom)))
590 gross 2264 df.setTolerance(rtol=self.TOL)
591     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
592 jfenwick 3551 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
593     self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
594 artak 1518
595 gross 2100 def testConstF_FreeBottom_largeK(self):
596     k=1.e10
597     mp=self.getScalarMask(include_bottom=False)
598     mv=self.getVectorMask(include_bottom=True)
599     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
600     p=p_ref*mp
601     u=u_ref*mv
602     df=DarcyFlow(self.dom)
603     df.setValue(g=f,
604     location_of_fixed_pressure=mp,
605     location_of_fixed_flux=mv,
606     permeability=Scalar(k,Function(self.dom)))
607 gross 2264 df.setTolerance(rtol=self.TOL)
608     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
609 jfenwick 3551 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
610     self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
611 artak 1518
612 gross 2100 def testVarioF_FreeBottom_smallK(self):
613     k=1.e-10
614     mp=self.getScalarMask(include_bottom=False)
615     mv=self.getVectorMask(include_bottom=True)
616     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
617     p=p_ref*mp
618     u=u_ref*mv
619     df=DarcyFlow(self.dom)
620     df.setValue(g=f,
621     location_of_fixed_pressure=mp,
622     location_of_fixed_flux=mv,
623     permeability=Scalar(k,Function(self.dom)))
624 gross 2264 df.setTolerance(rtol=self.TOL)
625     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
626 jfenwick 3551 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
627     self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
628 artak 1518
629 gross 2100 def testVarioF_FreeBottom_mediumK(self):
630     k=1.
631     mp=self.getScalarMask(include_bottom=False)
632     mv=self.getVectorMask(include_bottom=True)
633     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
634     p=p_ref*mp
635     u=u_ref*mv
636     df=DarcyFlow(self.dom)
637     df.setValue(g=f,
638     location_of_fixed_pressure=mp,
639     location_of_fixed_flux=mv,
640     permeability=Scalar(k,Function(self.dom)))
641 gross 2264 df.setTolerance(rtol=self.TOL)
642     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
643 jfenwick 3551 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
644     self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
645 artak 1518
646 gross 2100 def testVarioF_FreeBottom_largeK(self):
647     k=1.e10
648     mp=self.getScalarMask(include_bottom=False)
649     mv=self.getVectorMask(include_bottom=True)
650     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
651     p=p_ref*mp
652     u=u_ref*mv
653     df=DarcyFlow(self.dom)
654     df.setValue(g=f,
655     location_of_fixed_pressure=mp,
656     location_of_fixed_flux=mv,
657     permeability=Scalar(k,Function(self.dom)))
658 gross 2264 df.setTolerance(rtol=self.TOL)
659     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
660 jfenwick 3551 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
661     self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
662 artak 1518
663 gross 2100 class Test_Darcy2D(Test_Darcy):
664 gross 2960 TOL=1e-6
665     TEST_TOL=2.e-3
666 gross 2100 WIDTH=1.
667     def setUp(self):
668 gross 2208 NE=40 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
669 gross 3071 self.dom = Rectangle(NE/2,NE)
670 gross 2100 self.rescaleDomain()
671     def tearDown(self):
672     del self.dom
673     class Test_Darcy3D(Test_Darcy):
674 gross 2960 TOL=1e-6
675 gross 2100 WIDTH=1.
676 gross 2960 TEST_TOL=4.e-3
677 gross 2100 def setUp(self):
678     NE=25 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
679     self.dom = Brick(NE,NE,NE)
680     self.rescaleDomain()
681     def tearDown(self):
682     del self.dom
683 artak 1518
684 gross 2301 class Test_Rheologies(unittest.TestCase):
685     """
686     this is the program used to generate the powerlaw tests:
687    
688     TAU_Y=100.
689     N=10
690     M=5
691    
692     def getE(tau):
693     if tau<=TAU_Y:
694     return 1./(0.5+20*sqrt(tau))
695     else:
696     raise ValueError,"out of range."
697     tau=[ i* (TAU_Y/N) for i in range(N+1)] + [TAU_Y for i in range(M)]
698     e=[ tau[i]/getE(tau[i]) for i in range(N+1)]
699     e+= [ (i+1+j)* (max(e)/N) for j in range(M) ]
700    
701     print tau
702     print e
703     """
704     TOL=1.e-8
705     def test_PowerLaw_Init(self):
706     pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
707    
708 jfenwick 3551 self.assertEqual(pl.getNumMaterials(),3,"num materials is wrong")
709     self.assertEqual(pl.validMaterialId(0),True,"material id 0 not found")
710     self.assertEqual(pl.validMaterialId(1),True,"material id 1 not found")
711     self.assertEqual(pl.validMaterialId(2),True,"material id 2 not found")
712     self.assertEqual(pl.validMaterialId(3),False,"material id 3 not found")
713 gross 2301
714 jfenwick 3551 self.assertEqual(pl.getFriction(),None,"initial friction wrong.")
715     self.assertEqual(pl.getTauY(),None,"initial tau y wrong.")
716 gross 2301 pl.setDruckerPragerLaw(tau_Y=10,friction=3)
717 jfenwick 3551 self.assertEqual(pl.getFriction(),3,"friction wrong.")
718     self.assertEqual(pl.getTauY(),10,"tau y wrong.")
719 gross 2301
720 jfenwick 3551 self.assertEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
721 gross 2301 pl.setElasticShearModulus(1000)
722 jfenwick 3551 self.assertEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
723 gross 2301
724     e=pl.getEtaN()
725 jfenwick 3551 self.assertEqual(len(e),3,"initial length of etaN is wrong.")
726     self.assertEqual(e,[None, None, None],"initial etaN are wrong.")
727     self.assertEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
728     self.assertEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
729     self.assertEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
730     self.assertRaises(ValueError, pl.getEtaN, 3)
731 gross 2301
732 jfenwick 3551 self.assertRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
733     self.assertRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
734     self.assertRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
735 gross 2301
736     pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
737 jfenwick 3551 self.assertEqual(pl.getPower(),[1,2,3],"powers are wrong.")
738     self.assertEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
739     self.assertEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
740 gross 2301
741     pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
742 jfenwick 3551 self.assertEqual(pl.getPower(),[4,2,3],"powers are wrong.")
743     self.assertEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
744     self.assertEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
745 gross 2301
746 jfenwick 3551 self.assertRaises(ValueError,pl.getPower,-1)
747     self.assertEqual(pl.getPower(0),4,"power 0 is wrong.")
748     self.assertEqual(pl.getPower(1),2,"power 1 is wrong.")
749     self.assertEqual(pl.getPower(2),3,"power 2 is wrong.")
750     self.assertRaises(ValueError,pl.getPower,3)
751 gross 2301
752 jfenwick 3551 self.assertRaises(ValueError,pl.getTauT,-1)
753     self.assertEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
754     self.assertEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
755     self.assertEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
756     self.assertRaises(ValueError,pl.getTauT,3)
757 gross 2301
758 jfenwick 3551 self.assertRaises(ValueError,pl.getEtaN,-1)
759     self.assertEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
760     self.assertEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
761     self.assertEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
762     self.assertRaises(ValueError,pl.getEtaN,3)
763 gross 2301
764     def checkResult(self,id,gamma_dot_, eta, tau_ref):
765 jfenwick 3551 self.assertTrue(eta>=0,"eta needs to be positive (test %s)"%id)
766 gross 2301 error=abs(gamma_dot_*eta-tau_ref)
767 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))
768 gross 2301
769     def test_PowerLaw_Linear(self):
770     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]
771     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]
772     pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
773     pl.setDruckerPragerLaw(tau_Y=100.)
774     pl.setPowerLaw(eta_N=2.)
775     pl.setEtaTolerance(self.TOL)
776 jfenwick 3774 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
777 gross 2301
778     def test_PowerLaw_QuadLarge(self):
779     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]
780 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]
781 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
782     pl.setDruckerPragerLaw(tau_Y=100.)
783     pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2])
784     pl.setEtaTolerance(self.TOL)
785 jfenwick 3774 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
786 gross 2301
787     def test_PowerLaw_QuadSmall(self):
788     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]
789 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]
790 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
791     pl.setDruckerPragerLaw(tau_Y=100.)
792     pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2])
793     pl.setEtaTolerance(self.TOL)
794 jfenwick 3774 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
795 gross 2301
796     def test_PowerLaw_CubeLarge(self):
797     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]
798 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]
799 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
800     pl.setDruckerPragerLaw(tau_Y=100.)
801     pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3])
802     pl.setEtaTolerance(self.TOL)
803 jfenwick 3774 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
804 gross 2301
805     def test_PowerLaw_CubeSmall(self):
806     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]
807 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]
808 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
809     pl.setDruckerPragerLaw(tau_Y=100.)
810     pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3])
811     pl.setEtaTolerance(self.TOL)
812 jfenwick 3774 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
813 gross 2301
814     def test_PowerLaw_QuadLarge_CubeLarge(self):
815     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]
816 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]
817    
818 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
819     pl.setDruckerPragerLaw(tau_Y=100.)
820     pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
821     pl.setEtaTolerance(self.TOL)
822 jfenwick 3774 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
823 gross 2301
824     def test_PowerLaw_QuadLarge_CubeSmall(self):
825     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]
826 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]
827    
828 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
829     pl.setDruckerPragerLaw(tau_Y=100.)
830     pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
831     pl.setEtaTolerance(self.TOL)
832 jfenwick 3774 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
833 gross 2301
834     def test_PowerLaw_QuadSmall_CubeLarge(self):
835     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
836 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]
837 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
838     pl.setDruckerPragerLaw(tau_Y=100.)
839     pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
840     pl.setEtaTolerance(self.TOL)
841 jfenwick 3774 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
842 gross 2301
843     def test_PowerLaw_QuadSmall_CubeSmall(self):
844     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
845 gross 2438 gamma_dot_s=[0.0, 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]
846 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
847     pl.setDruckerPragerLaw(tau_Y=100.)
848     pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
849     pl.setEtaTolerance(self.TOL)
850 jfenwick 3774 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
851 gross 2301
852     def test_PowerLaw_withShear(self):
853     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
854     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]
855     pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
856     pl.setDruckerPragerLaw(tau_Y=100.)
857     pl.setPowerLaw(eta_N=2.)
858     pl.setElasticShearModulus(3.)
859     dt=1./3.
860     pl.setEtaTolerance(self.TOL)
861 jfenwick 3551 self.assertRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
862 jfenwick 3774 for i in range(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
863 gross 2301
864 gross 2432
865 gross 2647 class Test_FaultSystem(unittest.TestCase):
866     EPS=1.e-8
867     NE=10
868 gross 2663 def test_Fault_MaxValue(self):
869 gross 2647 dom=Rectangle(2*self.NE,2*self.NE)
870     x=dom.getX()
871     f=FaultSystem(dim=2)
872 gross 2663 f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
873     f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
874 gross 2647
875 gross 2676 u=x[0]*(1.-x[0])*(1-x[1])
876     t, loc=f.getMaxValue(u)
877     p=f.getParametrization(x,t)[0]
878     m, l=loc(u), loc(p)
879 jfenwick 3551 self.assertTrue( m == 0.25, "wrong max value")
880     self.assertTrue( t == 1, "wrong max tag")
881     self.assertTrue( l == 0., "wrong max location")
882 gross 2676
883     u=x[1]*(1.-x[1])*(1-x[0])*x[0]
884     t, loc=f.getMaxValue(u)
885     p=f.getParametrization(x,t)[0]
886     m, l=loc(u), loc(p)
887 jfenwick 3551 self.assertTrue( m == 0.0625, "wrong max value")
888     self.assertTrue( t == 2, "wrong max tag")
889     self.assertTrue( l == 0.5, "wrong max location")
890 gross 2676
891     u=x[0]*(1.-x[0])*x[1]
892     t, loc=f.getMaxValue(u)
893     p=f.getParametrization(x,t)[0]
894     m, l=loc(u), loc(p)
895 jfenwick 3551 self.assertTrue( m == 0.25, "wrong max value")
896     self.assertTrue( t == 2, "wrong max tag")
897     self.assertTrue( l == 1.0, "wrong max location")
898 gross 2676
899     u=x[1]*(1.-x[1])*x[0]
900     t, loc=f.getMaxValue(u)
901     p=f.getParametrization(x,t)[0]
902     m, l=loc(u), loc(p)
903 jfenwick 3551 self.assertTrue( m == 0.25, "wrong max value")
904     self.assertTrue( t == 2, "wrong max tag")
905     self.assertTrue( l == 0., "wrong max location")
906 gross 2676
907     u=x[1]*(1.-x[1])*(1.-x[0])
908     t, loc=f.getMaxValue(u)
909     p=f.getParametrization(x,t)[0]
910     m, l=loc(u), loc(p)
911 jfenwick 3551 self.assertTrue( m == 0.25, "wrong max value")
912     self.assertTrue( t == 1, "wrong max tag")
913     self.assertTrue( abs(l-0.70710678118654) <= self.EPS, "wrong max location")
914 gross 2675 def test_Fault_MinValue(self):
915     dom=Rectangle(2*self.NE,2*self.NE)
916     x=dom.getX()
917     f=FaultSystem(dim=2)
918     f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
919     f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
920 gross 2647
921 gross 2676 u=-x[0]*(1.-x[0])*(1-x[1])
922     t, loc=f.getMinValue(u)
923     p=f.getParametrization(x,t)[0]
924     m, l=loc(u), loc(p)
925 jfenwick 3551 self.assertTrue( m == -0.25, "wrong min value")
926     self.assertTrue( t == 1, "wrong min tag")
927     self.assertTrue( l == 0., "wrong min location")
928 gross 2676 u=-x[1]*(1.-x[1])*(1-x[0])*x[0]
929     t, loc=f.getMinValue(u)
930     p=f.getParametrization(x,t)[0]
931     m, l=loc(u), loc(p)
932 jfenwick 3551 self.assertTrue( m == -0.0625, "wrong min value")
933     self.assertTrue( t == 2, "wrong min tag")
934     self.assertTrue( l == 0.5, "wrong min location")
935 gross 2676 u=-x[0]*(1.-x[0])*x[1]
936     t, loc=f.getMinValue(u)
937     p=f.getParametrization(x,t)[0]
938     m, l=loc(u), loc(p)
939 jfenwick 3551 self.assertTrue( m == -0.25, "wrong min value")
940     self.assertTrue( t == 2, "wrong min tag")
941     self.assertTrue( l == 1.0, "wrong min location")
942 gross 2676 u=-x[1]*(1.-x[1])*x[0]
943     t, loc=f.getMinValue(u)
944     p=f.getParametrization(x,t)[0]
945     m, l=loc(u), loc(p)
946 jfenwick 3551 self.assertTrue( m == -0.25, "wrong min value")
947     self.assertTrue( t == 2, "wrong min tag")
948     self.assertTrue( l == 0., "wrong min location")
949 gross 2676 u=-x[1]*(1.-x[1])*(1.-x[0])
950     t, loc=f.getMinValue(u)
951     p=f.getParametrization(x,t)[0]
952     m, l=loc(u), loc(p)
953 jfenwick 3551 self.assertTrue( m == -0.25, "wrong min value")
954     self.assertTrue( t == 1, "wrong min tag")
955     self.assertTrue( abs(l-0.70710678118654) <= self.EPS, "wrong min location")
956 gross 2675
957 gross 2647
958 gross 2663 def test_Fault2D(self):
959 gross 2647 f=FaultSystem(dim=2)
960     top1=[ [1.,0.], [1.,1.], [0.,1.] ]
961 jfenwick 3551 self.assertRaises(ValueError,f.addFault,V0=[1.,0],strikes=[pi/2, pi/2],ls=[1.,1.],tag=1,dips=top1)
962 gross 2663 f.addFault(V0=[1.,0],strikes=[pi/2, pi],ls=[1.,1.],tag=1)
963 jfenwick 3551 self.assertTrue(f.getDim() == 2, "wrong dimension")
964     self.assertTrue( [ 1 ] == f.getTags(), "tags wrong")
965     self.assertTrue( 2. == f.getTotalLength(1), "length wrong")
966     self.assertTrue( 0. == f.getMediumDepth(1), "depth wrong")
967     self.assertTrue( (0., 2.) == f.getW0Range(1)," wrong W0 range")
968     self.assertTrue( (0., 0.) == f.getW1Range(1)," wrong W1 range")
969     self.assertTrue( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsets")
970 gross 2663 segs=f.getTopPolyline(1)
971 jfenwick 3551 self.assertTrue( len(segs) == 3, "wrong number of segments")
972     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
973     self.assertTrue( segs[0].size == 2, "seg 0 has wrong size.")
974     self.assertTrue( numpy.linalg.norm(segs[0]-[1.,0.]) < self.EPS, "wrong vertex. 0 ")
975     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
976     self.assertTrue( segs[1].size == 2, "seg 1 has wrong size.")
977     self.assertTrue( numpy.linalg.norm(segs[1]-[1.,1.]) < self.EPS, "wrong vertex. 1 ")
978     self.assertTrue( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
979     self.assertTrue( segs[2].size == 2, "seg 2 has wrong size.")
980     self.assertTrue( numpy.linalg.norm(segs[2]-[0.,1.]) < self.EPS, "wrong vertex. 2 ")
981 gross 2647 c=f.getCenterOnSurface()
982 jfenwick 3551 self.assertTrue( isinstance(c, numpy.ndarray), "center has wrong class")
983     self.assertTrue( c.size == 2, "center size is wrong")
984     self.assertTrue( numpy.linalg.norm(c-[2./3.,2./3.]) < self.EPS, "center has wrong coordinates.")
985 gross 2647 o=f.getOrientationOnSurface()/pi*180.
986 jfenwick 3551 self.assertTrue( abs(o+45.) < self.EPS, "wrong orientation.")
987 gross 2647
988     top2=[ [10.,0.], [0.,10.] ]
989 gross 2663 f.addFault(V0=[10.,0],strikes=[3.*pi/4],ls=[14.142135623730951], tag=2, w0_offsets=[0,20], w1_max=20)
990 jfenwick 3551 self.assertTrue( [ 1, 2 ] == f.getTags(), "tags wrong")
991     self.assertTrue( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length")
992     self.assertTrue( 0. == f.getMediumDepth(2), "depth wrong")
993     self.assertTrue( (0., 20.) == f.getW0Range(2)," wrong W0 range")
994     self.assertTrue( (0., 0.) == f.getW1Range(2)," wrong W1 range")
995     self.assertTrue( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets")
996 gross 2663 segs=f.getTopPolyline(2)
997 jfenwick 3551 self.assertTrue( len(segs) == 2, "wrong number of segments")
998     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
999     self.assertTrue( numpy.linalg.norm(segs[0]-[10.,0.]) < self.EPS, "wrong vertex. 0 ")
1000     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1001     self.assertTrue( numpy.linalg.norm(segs[1]-[0.,10.]) < self.EPS, "wrong vertex. 1 ")
1002 gross 2647 c=f.getCenterOnSurface()
1003 jfenwick 3551 self.assertTrue( isinstance(c, numpy.ndarray), "center has wrong class")
1004     self.assertTrue( c.size == 2, "center size is wrong")
1005     self.assertTrue( numpy.linalg.norm(c-[12./5.,12./5.]) < self.EPS, "center has wrong coordinates.")
1006 gross 2647 o=f.getOrientationOnSurface()/pi*180.
1007 jfenwick 3551 self.assertTrue( abs(o+45.) < self.EPS, "wrong orientation.")
1008 gross 2647
1009 gross 2654 s,d=f.getSideAndDistance([0.,0.], tag=1)
1010 jfenwick 3551 self.assertTrue( s<0, "wrong side.")
1011     self.assertTrue( abs(d-1.)<self.EPS, "wrong distance.")
1012 gross 2654 s,d=f.getSideAndDistance([0.,2.], tag=1)
1013 jfenwick 3551 self.assertTrue( s>0, "wrong side.")
1014     self.assertTrue( abs(d-1.)<self.EPS, "wrong distance.")
1015 gross 2654 s,d=f.getSideAndDistance([1.,2.], tag=1)
1016 jfenwick 3551 self.assertTrue( s>0, "wrong side.")
1017     self.assertTrue( abs(d-1.)<self.EPS, "wrong distance.")
1018 gross 2654 s,d=f.getSideAndDistance([2.,1.], tag=1)
1019 jfenwick 3551 self.assertTrue( s>0, "wrong side.")
1020     self.assertTrue( abs(d-1.)<self.EPS, "wrong distance.")
1021 gross 2654 s,d=f.getSideAndDistance([2.,0.], tag=1)
1022 jfenwick 3551 self.assertTrue( s>0, "wrong side.")
1023     self.assertTrue( abs(d-1.)<self.EPS, "wrong distance.")
1024 gross 2654 s,d=f.getSideAndDistance([0.,-1.], tag=1)
1025 jfenwick 3551 self.assertTrue( s<0, "wrong side.")
1026     self.assertTrue( abs(d-1.41421356237)<self.EPS, "wrong distance.")
1027 gross 2654 s,d=f.getSideAndDistance([-1.,0], tag=1)
1028 jfenwick 3551 self.assertTrue( s<0, "wrong side.")
1029     self.assertTrue( abs(d-1.41421356237)<self.EPS, "wrong distance.")
1030 gross 2654
1031    
1032 gross 2647 f.transform(rot=-pi/2., shift=[-1.,-1.])
1033 jfenwick 3551 self.assertTrue( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
1034     self.assertTrue( 2. == f.getTotalLength(1), "length after transformation wrong")
1035     self.assertTrue( 0. == f.getMediumDepth(1), "depth after transformation wrong")
1036     self.assertTrue( (0., 2.) == f.getW0Range(1)," wrong W0 after transformation range")
1037     self.assertTrue( (0., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
1038     self.assertTrue( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsetsafter transformation ")
1039 gross 2663 segs=f.getTopPolyline(1)
1040 jfenwick 3551 self.assertTrue( len(segs) == 3, "wrong number of segmentsafter transformation ")
1041     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1042     self.assertTrue( segs[0].size == 2, "seg 0 has wrong size after transformation.")
1043     self.assertTrue( numpy.linalg.norm(segs[0]-[-1.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1044     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1045     self.assertTrue( segs[1].size == 2, "seg 1 has wrong size after transformation.")
1046     self.assertTrue( numpy.linalg.norm(segs[1]-[0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1047     self.assertTrue( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1048     self.assertTrue( segs[2].size == 2, "seg 2 has wrong size after transformation.")
1049     self.assertTrue( numpy.linalg.norm(segs[2]-[0., 1.]) < self.EPS, "wrong vertex after transformation. 2 ")
1050     self.assertTrue( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1051     self.assertTrue( 0. == f.getMediumDepth(2), "depth wrong after transformation")
1052     self.assertTrue( (0., 20.) == f.getW0Range(2)," wrong W0 range after transformation")
1053     self.assertTrue( (0., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1054     self.assertTrue( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1055 gross 2663 segs=f.getTopPolyline(2)
1056 jfenwick 3551 self.assertTrue( len(segs) == 2, "wrong number of segments after transformation")
1057     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1058     self.assertTrue( numpy.linalg.norm(segs[0]-[-1.,-9]) < self.EPS, "wrong vertex. 0 after transformation")
1059     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1060     self.assertTrue( numpy.linalg.norm(segs[1]-[9.,1.]) < self.EPS, "wrong vertex. 1 after transformation")
1061 gross 2647
1062     c=f.getCenterOnSurface()
1063 jfenwick 3551 self.assertTrue( isinstance(c, numpy.ndarray), "center has wrong class")
1064     self.assertTrue( c.size == 2, "center size is wrong")
1065     self.assertTrue( numpy.linalg.norm(c-[7./5.,-7./5.]) < self.EPS, "center has wrong coordinates.")
1066 gross 2647 o=f.getOrientationOnSurface()/pi*180.
1067 jfenwick 3551 self.assertTrue( abs(o-45.) < self.EPS, "wrong orientation.")
1068 gross 2647
1069     p=f.getParametrization([-1.,0.],1)
1070 jfenwick 3551 self.assertTrue(p[1]==1., "wrong value.")
1071     self.assertTrue(abs(p[0])<self.EPS, "wrong value.")
1072 gross 2647 p=f.getParametrization([-0.5,0.],1)
1073 jfenwick 3551 self.assertTrue(p[1]==1., "wrong value.")
1074     self.assertTrue(abs(p[0]-0.5)<self.EPS* 0.5, "wrong value.")
1075 gross 2647 p=f.getParametrization([0.,0.],1)
1076 jfenwick 3551 self.assertTrue(p[1]==1., "wrong value.")
1077     self.assertTrue(abs(p[0]-1.)<self.EPS, "wrong value.")
1078 gross 2647 p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-8)
1079 jfenwick 3551 self.assertTrue(p[1]==0., "wrong value.")
1080 gross 2647 p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-6)
1081 jfenwick 3551 self.assertTrue(p[1]==1., "wrong value.")
1082     self.assertTrue(abs(p[0]-1.0000001)<self.EPS, "wrong value.")
1083 gross 2647 p=f.getParametrization([0.,0.5],1)
1084 jfenwick 3551 self.assertTrue(p[1]==1., "wrong value.")
1085     self.assertTrue(abs(p[0]-1.5)<self.EPS, "wrong value.")
1086 gross 2647 p=f.getParametrization([0,1.],1)
1087 jfenwick 3551 self.assertTrue(p[1]==1., "wrong value.")
1088     self.assertTrue(abs(p[0]-2.)<self.EPS, "wrong value.")
1089 gross 2647 p=f.getParametrization([1.,1.],1)
1090 jfenwick 3551 self.assertTrue(p[1]==0., "wrong value.")
1091 gross 2647 p=f.getParametrization([0,1.11],1)
1092 jfenwick 3551 self.assertTrue(p[1]==0., "wrong value.")
1093 gross 2647 p=f.getParametrization([-1,-9.],2)
1094 jfenwick 3551 self.assertTrue(p[1]==1., "wrong value.")
1095     self.assertTrue(abs(p[0])<self.EPS, "wrong value.")
1096 gross 2647 p=f.getParametrization([9,1],2)
1097 jfenwick 3551 self.assertTrue(p[1]==1., "wrong value.")
1098     self.assertTrue(abs(p[0]-20.)<self.EPS, "wrong value.")
1099 gross 2647
1100 gross 2663 def test_Fault3D(self):
1101     f=FaultSystem(dim=3)
1102 jfenwick 3551 self.assertTrue(f.getDim() == 3, "wrong dimension")
1103 gross 2647
1104 gross 2663 top1=[ [0.,0.,0.], [1., 0., 0.] ]
1105     f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=20.,tag=1)
1106 jfenwick 3551 self.assertTrue( [ 1 ] == f.getTags(), "tags wrong")
1107     self.assertTrue( 1. == f.getTotalLength(1), "length wrong")
1108     self.assertTrue( 20. == f.getMediumDepth(1), "depth wrong")
1109     self.assertTrue( (0., 1.) == f.getW0Range(1)," wrong W0 range")
1110     self.assertTrue( (-20., 0.) == f.getW1Range(1)," wrong W1 range")
1111     self.assertTrue( [0., 1.] == f.getW0Offsets(1)," wrong W0 offsets")
1112 gross 2663 segs=f.getTopPolyline(1)
1113 jfenwick 3551 self.assertTrue( len(segs) == 2, "wrong number of segments")
1114     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1115     self.assertTrue( segs[0].size == 3, "seg 0 has wrong size.")
1116     self.assertTrue( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1117     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1118     self.assertTrue( segs[1].size == 3, "seg 1 has wrong size.")
1119     self.assertTrue( numpy.linalg.norm(segs[1]-[1.,0.,0]) < self.EPS, "wrong vertex. 1 ")
1120 gross 2663 c=f.getCenterOnSurface()
1121 jfenwick 3551 self.assertTrue( isinstance(c, numpy.ndarray), "center has wrong class")
1122     self.assertTrue( c.size == 3, "center size is wrong")
1123     self.assertTrue( numpy.linalg.norm(c-[0.5,0.,0.]) < self.EPS, "center has wrong coordinates.")
1124 gross 2663 o=f.getOrientationOnSurface()/pi*180.
1125 jfenwick 3551 self.assertTrue( abs(o) < self.EPS, "wrong orientation.")
1126 gross 2663 d=f.getDips(1)
1127 jfenwick 3551 self.assertTrue( len(d) == 1, "wrong number of dips")
1128     self.assertTrue( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1129 gross 2663 sn=f.getSegmentNormals(1)
1130 jfenwick 3551 self.assertTrue( len(sn) == 1, "wrong number of normals")
1131     self.assertTrue( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1132     self.assertTrue( numpy.linalg.norm(sn[0]-[0, -1., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1133 gross 2663 dv=f.getDepthVectors(1)
1134 jfenwick 3551 self.assertTrue( len(dv) == 2, "wrong number of depth vectors.")
1135     self.assertTrue( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1136     self.assertTrue( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1137     self.assertTrue( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1138     self.assertTrue( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1139 gross 2663 b=f.getBottomPolyline(1)
1140 jfenwick 3551 self.assertTrue( len(b) == 2, "wrong number of bottom vertices")
1141     self.assertTrue( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1142     self.assertTrue( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1143     self.assertTrue( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1144     self.assertTrue( numpy.linalg.norm(b[1]-[1., 0., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1145 gross 2663 ds=f.getDepths(1)
1146 jfenwick 3551 self.assertTrue( len(ds) == 2, "wrong number of depth")
1147     self.assertTrue( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1148     self.assertTrue( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1149 gross 2647
1150 gross 2663 top2=[ [0.,0.,0.], [0., 10., 0.] ]
1151     f.addFault(V0=[0.,0,0],strikes=[pi/2],ls=[10.,], dips=pi/2, depths=20.,tag=2)
1152 jfenwick 3551 self.assertTrue( [ 1, 2 ] == f.getTags(), "tags wrong")
1153     self.assertTrue( 10. == f.getTotalLength(2), "length wrong")
1154     self.assertTrue( 20. == f.getMediumDepth(2), "depth wrong")
1155     self.assertTrue( (0., 10.) == f.getW0Range(2)," wrong W0 range")
1156     self.assertTrue( (-20., 0.) == f.getW1Range(2)," wrong W1 range")
1157     self.assertTrue( [0., 10.] == f.getW0Offsets(2)," wrong W0 offsets")
1158 gross 2663 segs=f.getTopPolyline(2)
1159 jfenwick 3551 self.assertTrue( len(segs) == 2, "wrong number of segments")
1160     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1161     self.assertTrue( segs[0].size == 3, "seg 0 has wrong size.")
1162     self.assertTrue( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1163     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1164     self.assertTrue( segs[1].size == 3, "seg 1 has wrong size.")
1165     self.assertTrue( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1166 gross 2663 d=f.getDips(2)
1167 jfenwick 3551 self.assertTrue( len(d) == 1, "wrong number of dips")
1168     self.assertTrue( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1169 gross 2663 sn=f.getSegmentNormals(2)
1170 jfenwick 3551 self.assertTrue( len(sn) == 1, "wrong number of normals")
1171     self.assertTrue( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1172     self.assertTrue( numpy.linalg.norm(sn[0]-[1, 0., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1173 gross 2663 dv=f.getDepthVectors(2)
1174 jfenwick 3551 self.assertTrue( len(dv) == 2, "wrong number of depth vectors.")
1175     self.assertTrue( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1176     self.assertTrue( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1177     self.assertTrue( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1178     self.assertTrue( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1179 gross 2663 b=f.getBottomPolyline(2)
1180 jfenwick 3551 self.assertTrue( len(b) == 2, "wrong number of bottom vertices")
1181     self.assertTrue( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1182     self.assertTrue( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1183     self.assertTrue( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1184     self.assertTrue( numpy.linalg.norm(b[1]-[0., 10., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1185 gross 2663 ds=f.getDepths(2)
1186 jfenwick 3551 self.assertTrue( len(ds) == 2, "wrong number of depth")
1187     self.assertTrue( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1188     self.assertTrue( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1189 gross 2647
1190 gross 2663 top2=[ [10.,0.,0.], [0., 10., 0.] ]
1191     f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/2, depths=30.,tag=2)
1192 jfenwick 3551 self.assertTrue( [ 1, 2 ] == f.getTags(), "tags wrong")
1193     self.assertTrue( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1194     self.assertTrue( 30. == f.getMediumDepth(2), "depth wrong")
1195     self.assertTrue( (-30., 0.) == f.getW1Range(2)," wrong W1 range")
1196 gross 2663 segs=f.getTopPolyline(2)
1197 jfenwick 3551 self.assertTrue( len(segs) == 2, "wrong number of segments")
1198     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1199     self.assertTrue( segs[0].size == 3, "seg 0 has wrong size.")
1200     self.assertTrue( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1201     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1202     self.assertTrue( segs[1].size == 3, "seg 1 has wrong size.")
1203     self.assertTrue( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1204 gross 2663 d=f.getDips(2)
1205 jfenwick 3551 self.assertTrue( len(d) == 1, "wrong number of dips")
1206     self.assertTrue( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1207 gross 2663 sn=f.getSegmentNormals(2)
1208 jfenwick 3551 self.assertTrue( len(sn) == 1, "wrong number of normals")
1209     self.assertTrue( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1210     self.assertTrue( numpy.linalg.norm(sn[0]-[0.70710678118654746, 0.70710678118654746, 0.]) < self.EPS, "wrong bottom vertex 1 ")
1211 gross 2663 dv=f.getDepthVectors(2)
1212 jfenwick 3551 self.assertTrue( len(dv) == 2, "wrong number of depth vectors.")
1213     self.assertTrue( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1214     self.assertTrue( numpy.linalg.norm(dv[0]-[0., 0., -30.]) < self.EPS, "wrong depth vector 0 ")
1215     self.assertTrue( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1216     self.assertTrue( numpy.linalg.norm(dv[1]-[0., 0., -30.]) < self.EPS, "wrong depth vector 1 ")
1217 gross 2663 b=f.getBottomPolyline(2)
1218 jfenwick 3551 self.assertTrue( len(b) == 2, "wrong number of bottom vertices")
1219     self.assertTrue( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1220     self.assertTrue( numpy.linalg.norm(b[0]-[10., 0., -30.]) < self.EPS, "wrong bottom vertex 0 ")
1221     self.assertTrue( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1222     self.assertTrue( numpy.linalg.norm(b[1]-[0., 10., -30.]) < self.EPS, "wrong bottom vertex 1 ")
1223 gross 2663 ds=f.getDepths(2)
1224 jfenwick 3551 self.assertTrue( len(ds) == 2, "wrong number of depth")
1225     self.assertTrue( abs(ds[0]-30.) < self.EPS, "wrong depth at vertex 0 ")
1226     self.assertTrue( abs(ds[1]-30.) < self.EPS, "wrong depth at vertex 1 ")
1227 gross 2663
1228     top2=[ [10.,0.,0.], [0., 10., 0.] ]
1229     f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/4, depths=50.,tag=2)
1230 jfenwick 3551 self.assertTrue( [ 1, 2 ] == f.getTags(), "tags wrong")
1231     self.assertTrue( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1232     self.assertTrue( 50. == f.getMediumDepth(2), "depth wrong")
1233     self.assertTrue( (-50., 0.) == f.getW1Range(2)," wrong W1 range")
1234 gross 2663 segs=f.getTopPolyline(2)
1235 jfenwick 3551 self.assertTrue( len(segs) == 2, "wrong number of segments")
1236     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1237     self.assertTrue( segs[0].size == 3, "seg 0 has wrong size.")
1238     self.assertTrue( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1239     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1240     self.assertTrue( segs[1].size == 3, "seg 1 has wrong size.")
1241     self.assertTrue( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1242 gross 2663 d=f.getDips(2)
1243 jfenwick 3551 self.assertTrue( len(d) == 1, "wrong number of dips")
1244     self.assertTrue( abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1245 gross 2663 sn=f.getSegmentNormals(2)
1246 jfenwick 3551 self.assertTrue( len(sn) == 1, "wrong number of normals")
1247     self.assertTrue( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1248     self.assertTrue( numpy.linalg.norm(sn[0]-[0.5,0.5,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1249 gross 2663 dv=f.getDepthVectors(2)
1250 jfenwick 3551 self.assertTrue( len(dv) == 2, "wrong number of depth vectors.")
1251     self.assertTrue( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1252     self.assertTrue( numpy.linalg.norm(dv[0]-[25., 25., -35.355339059327378]) < self.EPS, "wrong depth vector 0 ")
1253     self.assertTrue( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1254     self.assertTrue( numpy.linalg.norm(dv[1]-[25.,25., -35.355339059327378]) < self.EPS, "wrong depth vector 1 ")
1255 gross 2663 b=f.getBottomPolyline(2)
1256 jfenwick 3551 self.assertTrue( len(b) == 2, "wrong number of bottom vertices")
1257     self.assertTrue( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1258     self.assertTrue( numpy.linalg.norm(b[0]-[35., 25., -35.355339059327378]) < self.EPS, "wrong bottom vertex 0 ")
1259     self.assertTrue( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1260     self.assertTrue( numpy.linalg.norm(b[1]-[25, 35., -35.355339059327378]) < self.EPS, "wrong bottom vertex 1 ")
1261 gross 2663 ds=f.getDepths(2)
1262 jfenwick 3551 self.assertTrue( len(ds) == 2, "wrong number of depth")
1263     self.assertTrue( abs(ds[0]-50.) < self.EPS, "wrong depth at vertex 0 ")
1264     self.assertTrue( abs(ds[1]-50.) < self.EPS, "wrong depth at vertex 1 ")
1265 gross 2663
1266     top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
1267     f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=1, dips=pi/4, depths=20.)
1268 jfenwick 3551 self.assertTrue( 20. == f.getTotalLength(1), "length wrong")
1269     self.assertTrue( 20. == f.getMediumDepth(1), "depth wrong")
1270 gross 2663 segs=f.getTopPolyline(1)
1271 jfenwick 3551 self.assertTrue( len(segs) == 3, "wrong number of segments")
1272     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1273     self.assertTrue( segs[0].size == 3, "seg 0 has wrong size.")
1274     self.assertTrue( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1275     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1276     self.assertTrue( segs[1].size == 3, "seg 1 has wrong size.")
1277     self.assertTrue( numpy.linalg.norm(segs[1]-[10.,10.,0.]) < self.EPS, "wrong vertex. 1 ")
1278     self.assertTrue( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1279     self.assertTrue( segs[2].size == 3, "seg 2 has wrong size.")
1280     self.assertTrue( numpy.linalg.norm(segs[2]-[0.,10.,0.]) < self.EPS, "wrong vertex. 2 ")
1281 gross 2663 d=f.getDips(1)
1282 jfenwick 3551 self.assertTrue( len(d) == 2, "wrong number of dips")
1283     self.assertTrue( abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1284     self.assertTrue( abs(d[1]-0.78539816339744828) < self.EPS, "wrong dip 0")
1285 gross 2663 ds=f.getDepths(1)
1286 jfenwick 3551 self.assertTrue( len(ds) == 3, "wrong number of depth")
1287     self.assertTrue( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1288     self.assertTrue( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1289 gross 2663 sn=f.getSegmentNormals(1)
1290 jfenwick 3551 self.assertTrue( len(sn) == 2, "wrong number of normals")
1291     self.assertTrue( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1292     self.assertTrue( numpy.linalg.norm(sn[0]-[0.70710678118654746,0.,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1293     self.assertTrue( isinstance(sn[1], numpy.ndarray), "wrong class of bottom vertex 0")
1294     self.assertTrue( numpy.linalg.norm(sn[1]-[0.,0.70710678118654746,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1295 gross 2663 dv=f.getDepthVectors(1)
1296 jfenwick 3551 self.assertTrue( len(dv) == 3, "wrong number of depth vectors.")
1297     self.assertTrue( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1298     self.assertTrue( numpy.linalg.norm(dv[0]-[14.142135623730951, 0., -14.142135623730951]) < self.EPS, "wrong depth vector 0 ")
1299     self.assertTrue( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1300     self.assertTrue( numpy.linalg.norm(dv[1]-[11.547005383792515,11.547005383792515, -11.547005383792515]) < self.EPS, "wrong depth vector 2 ")
1301     self.assertTrue( isinstance(dv[2], numpy.ndarray), "wrong class of depth vector 1")
1302     self.assertTrue( numpy.linalg.norm(dv[2]-[0.,14.142135623730951, -14.142135623730951]) < self.EPS, "wrong depth vector 2 ")
1303 gross 2663 segs=f.getBottomPolyline(1)
1304 jfenwick 3551 self.assertTrue( len(segs) == 3, "wrong number of segments")
1305     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1306     self.assertTrue( segs[0].size == 3, "seg 0 has wrong size.")
1307     self.assertTrue( numpy.linalg.norm(segs[0]-[24.142135623730951,0.,-14.142135623730951]) < self.EPS, "wrong vertex. 0 ")
1308     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1309     self.assertTrue( segs[1].size == 3, "seg 1 has wrong size.")
1310     self.assertTrue( numpy.linalg.norm(segs[1]-[21.547005383792515,21.547005383792515, -11.547005383792515]) < self.EPS, "wrong vertex. 1 ")
1311     self.assertTrue( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1312     self.assertTrue( segs[2].size == 3, "seg 2 has wrong size.")
1313     self.assertTrue( numpy.linalg.norm(segs[2]-[0., 24.142135623730951, -14.142135623730951]) < self.EPS, "wrong vertex. 2 ")
1314     self.assertTrue( abs(0.-f.getW0Range(1)[0]) <=self.EPS," wrong W0 range (0)")
1315     self.assertTrue( abs(31.857329272664341-f.getW0Range(1)[1]) <=self.EPS," wrong W0 range (1)")
1316     self.assertTrue( abs(-20. - f.getW1Range(1)[0]) <=self.EPS," wrong W1 range (0)")
1317     self.assertTrue( abs(0. - f.getW1Range(1)[1]) <=self.EPS," wrong W1 range (1)")
1318     self.assertTrue( abs(0.0-f.getW0Offsets(1)[0])<=self.EPS," wrong W0 offsets (0)")
1319     self.assertTrue( abs(15.92866463633217-f.getW0Offsets(1)[1])<=self.EPS," wrong W0 offsets (1)")
1320     self.assertTrue( abs(31.857329272664341-f.getW0Offsets(1)[2])<=self.EPS," wrong W0 offsets(2)")
1321 gross 2663 #
1322     # ============ fresh start ====================
1323     #
1324     f.addFault(V0=[1.,0,0.],strikes=[pi/2, pi],ls=[1.,1.], dips=pi/2,depths=20,tag=1)
1325     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)
1326     c=f.getCenterOnSurface()
1327 jfenwick 3551 self.assertTrue( isinstance(c, numpy.ndarray), "center has wrong class")
1328     self.assertTrue( c.size == 3, "center size is wrong")
1329     self.assertTrue( numpy.linalg.norm(c-[12./5.,12./5.,0.]) < self.EPS, "center has wrong coordinates.")
1330 gross 2663 o=f.getOrientationOnSurface()/pi*180.
1331 jfenwick 3551 self.assertTrue( abs(o+45.) < self.EPS, "wrong orientation.")
1332 gross 2663
1333     f.transform(rot=-pi/2., shift=[-1.,-1.,0.])
1334 jfenwick 3551 self.assertTrue( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
1335     self.assertTrue( 2. == f.getTotalLength(1), "length after transformation wrong")
1336     self.assertTrue( 20. == f.getMediumDepth(1), "depth after transformation wrong")
1337 gross 2663 rw0=f.getW0Range(1)
1338 jfenwick 3551 self.assertTrue( len(rw0) ==2, "wo range has wrong length")
1339     self.assertTrue( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
1340     self.assertTrue( abs(rw0[1]-2.) < self.EPS,"W0 1 wrong.")
1341     self.assertTrue( (-20., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
1342 gross 2663 dips=f.getDips(1)
1343 jfenwick 3551 self.assertTrue(len(dips) == 2, "wrong number of dips.")
1344     self.assertTrue( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
1345     self.assertTrue( abs(dips[1]-1.5707963267948966) <= self.EPS, "wrong dip")
1346 gross 2663 ds=f.getDepths(1)
1347 jfenwick 3551 self.assertTrue( len(ds) == 3, "wrong number of depth")
1348     self.assertTrue( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1349     self.assertTrue( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1350     self.assertTrue( abs(ds[2]-20.) < self.EPS, "wrong depth at vertex 1 ")
1351 gross 2663 segs=f.getTopPolyline(1)
1352 jfenwick 3551 self.assertTrue( len(segs) == 3, "wrong number of segmentsafter transformation ")
1353     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1354     self.assertTrue( segs[0].size == 3, "seg 0 has wrong size after transformation.")
1355     self.assertTrue( numpy.linalg.norm(segs[0]-[-1.,0.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1356     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1357     self.assertTrue( segs[1].size == 3, "seg 1 has wrong size after transformation.")
1358     self.assertTrue( numpy.linalg.norm(segs[1]-[0.,0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1359     self.assertTrue( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1360     self.assertTrue( segs[2].size == 3, "seg 2 has wrong size after transformation.")
1361     self.assertTrue( numpy.linalg.norm(segs[2]-[0., 1.,0.]) < self.EPS, "wrong vertex after transformation. 2 ")
1362     self.assertTrue( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1363     self.assertTrue( 20. == f.getMediumDepth(2), "depth wrong after transformation")
1364 gross 2663 rw0=f.getW0Range(2)
1365 jfenwick 3551 self.assertTrue( len(rw0) ==2, "wo range has wrong length")
1366     self.assertTrue( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
1367     self.assertTrue( abs(rw0[1]-20.) < self.EPS,"W0 1 wrong.")
1368     self.assertTrue( (-20., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1369     self.assertTrue( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1370 gross 2663 dips=f.getDips(2)
1371 jfenwick 3551 self.assertTrue(len(dips) == 1, "wrong number of dips.")
1372     self.assertTrue( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
1373 gross 2663 ds=f.getDepths(2)
1374 jfenwick 3551 self.assertTrue( len(ds) == 2, "wrong number of depth")
1375     self.assertTrue( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1376     self.assertTrue( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1377 gross 2663 segs=f.getTopPolyline(2)
1378 jfenwick 3551 self.assertTrue( len(segs) == 2, "wrong number of segments after transformation")
1379     self.assertTrue( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1380     self.assertTrue( numpy.linalg.norm(segs[0]-[-1.,-9,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1381     self.assertTrue( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1382     self.assertTrue( numpy.linalg.norm(segs[1]-[9.,1.,0.]) < self.EPS, "wrong vertex. 1 after transformation")
1383 gross 2663 #
1384     # ============ fresh start ====================
1385     #
1386     f=FaultSystem(dim=3)
1387    
1388     top1=[ [0.,0.,0.], [1., 0., 0.] ]
1389     f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=1.,tag=1)
1390     top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
1391     f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=2, dips=pi/4, depths=20.)
1392    
1393     p,m=f.getParametrization([0.3,0.,-0.5],1)
1394 jfenwick 3551 self.assertTrue(length(p-[0.3,-0.5]) <= self.EPS, "wrong value.")
1395     self.assertTrue(m==1., "wrong value.")
1396 gross 2663
1397     p,m=f.getParametrization([0.5,0.,-0.5],1)
1398 jfenwick 3551 self.assertTrue(length(p-[0.5,-0.5]) <= self.EPS, "wrong value.")
1399     self.assertTrue(m==1., "wrong value.")
1400 gross 2663
1401     p,m=f.getParametrization([0.25,0.,-0.5],1)
1402 jfenwick 3551 self.assertTrue(length(p-[0.25,-0.5]) <= self.EPS, "wrong value.")
1403     self.assertTrue(m==1., "wrong value.")
1404 gross 2663
1405     p,m=f.getParametrization([0.5,0.,-0.25],1)
1406 jfenwick 3551 self.assertTrue(length(p-[0.5,-0.25]) <= self.EPS, "wrong value.")
1407     self.assertTrue(m==1., "wrong value.")
1408 gross 2663
1409     p,m=f.getParametrization([0.001,0.