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

Contents of /trunk/finley/test/python/run_models.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3527 - (show annotations)
Tue May 24 06:59:20 2011 UTC (8 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 73991 byte(s)
modifications to the coal seam gas 
1 # -*- coding: utf-8 -*-
2
3 ########################################################
4 #
5 # Copyright (c) 2003-2010 by University of Queensland
6 # Earth Systems Science Computational Center (ESSCC)
7 # http://www.uq.edu.au/esscc
8 #
9 # 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 #
13 ########################################################
14
15 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
16 Earth Systems Science Computational Center (ESSCC)
17 http://www.uq.edu.au/esscc
18 Primary Business: Queensland, Australia"""
19 __license__="""Licensed under the Open Software License version 3.0
20 http://www.opensource.org/licenses/osl-3.0.php"""
21 __url__="https://launchpad.net/escript-finley"
22
23 import unittest
24 import tempfile
25
26
27
28 VERBOSE=False and True
29
30 from esys.escript import *
31 from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian, FaultSystem, DarcyFlow
32 from esys.escript.models import Mountains
33 from esys.finley import Rectangle, Brick
34
35 from math import pi
36 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
45 #====================================================================================================================
46 class Test_StokesProblemCartesian2D(unittest.TestCase):
47 def setUp(self):
48 NE=6
49 self.TOL=1e-3
50 self.domain=Rectangle(NE,NE,order=-1)
51 def tearDown(self):
52 del self.domain
53 def test_PCG_P_0(self):
54 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 p0=Scalar(-P1,ReducedSolution(self.domain))
69 sp.setTolerance(self.TOL)
70 u,p=sp.solve(u0*mask,p0,verbose=VERBOSE,max_iter=100,usePCG=True)
71
72 error_v0=Lsup(u[0]-u0[0])
73 error_v1=Lsup(u[1]-u0[1])/0.25
74 error_p=Lsup(p+P1*x[0]*x[1])
75 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
76 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
77 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
78
79 def test_PCG_P_small(self):
80 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 p0=Scalar(-P1,ReducedSolution(self.domain))
95 sp.setTolerance(self.TOL)
96 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
97 error_v0=Lsup(u[0]-u0[0])
98 error_v1=Lsup(u[1]-u0[1])/0.25
99 error_p=Lsup(P1*x[0]*x[1]+p)
100 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
101 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
102 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
103
104 def test_PCG_P_large(self):
105 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 p0=Scalar(-P1,ReducedSolution(self.domain))
120 sp.setTolerance(self.TOL)
121 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
122 # u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
123
124 error_v0=Lsup(u[0]-u0[0])
125 error_v1=Lsup(u[1]-u0[1])/0.25
126 error_p=Lsup(P1*x[0]*x[1]+p)/P1
127 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
128 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
129 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
130
131 def test_GMRES_P_0(self):
132 ETA=1.
133 P1=0.
134
135 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 p0=Scalar(-P1,ReducedSolution(self.domain))
147 sp.setTolerance(self.TOL)
148 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=50,usePCG=False,iter_restart=18)
149
150 error_v0=Lsup(u[0]-u0[0])
151 error_v1=Lsup(u[1]-u0[1])/0.25
152 error_p=Lsup(P1*x[0]*x[1]+p)
153 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
154 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
155 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
156
157 def test_GMRES_P_small(self):
158 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 p0=Scalar(-P1,ReducedSolution(self.domain))
173 sp.setTolerance(self.TOL)
174 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=20,usePCG=False)
175
176 error_v0=Lsup(u[0]-u0[0])
177 error_v1=Lsup(u[1]-u0[1])/0.25
178 error_p=Lsup(P1*x[0]*x[1]+p)
179
180 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
181 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
182 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
183
184 def test_GMRES_P_large(self):
185 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 p0=Scalar(-P1,ReducedSolution(self.domain))
200 sp.setTolerance(self.TOL)
201 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)
202
203 error_v0=Lsup(u[0]-u0[0])
204 error_v1=Lsup(u[1]-u0[1])/0.25
205 error_p=Lsup(P1*x[0]*x[1]+p)/P1
206 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
207 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
208 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
209 #====================================================================================================================
210 class Test_StokesProblemCartesian3D(unittest.TestCase):
211 def setUp(self):
212 NE=6
213 self.TOL=1e-4
214 self.domain=Brick(NE,NE,NE,order=-1)
215 def tearDown(self):
216 del self.domain
217 def test_PCG_P_0(self):
218 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 +whereZero(x[1]) * [1.,0.,1.] \
227 +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 p0=Scalar(-P1,ReducedSolution(self.domain))
237 sp.setTolerance(self.TOL)
238 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
239
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 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
244 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
245 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
246 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
247 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
248
249 def test_PCG_P_small(self):
250 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 +whereZero(x[1]) * [1.,0.,1.] \
258 +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 p0=Scalar(-P1,ReducedSolution(self.domain))
268 sp.setTolerance(self.TOL)
269 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
270 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 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
274 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
275 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
276 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
277 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
278
279 def test_PCG_P_large(self):
280 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 +whereZero(x[1]) * [1.,0.,1.] \
288 +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 p0=Scalar(-P1,ReducedSolution(self.domain))
298 sp.setTolerance(self.TOL)
299 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
300
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 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
305 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
306 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
307 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
308 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
309
310 def test_GMRES_P_0(self):
311 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 p0=Scalar(-P1,ReducedSolution(self.domain))
330 sp.setTolerance(self.TOL)
331 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False,iter_restart=20)
332
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 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
337 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
338 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
339 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
340 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
341 def test_GMRES_P_small(self):
342 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 p0=Scalar(-P1,ReducedSolution(self.domain))
360 sp.setTolerance(self.TOL/10)
361 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)
362
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 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
367 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
368 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
369 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
370 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
371 def test_GMRES_P_large(self):
372 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 +whereZero(x[1]) * [1.,0.,1.] \
380 +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 p0=Scalar(-P1,ReducedSolution(self.domain))
390 sp.setTolerance(self.TOL)
391 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=False)
392
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 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
397 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
398 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
399 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
400 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
401 #====================================================================================================================
402
403 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
413 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 mts=Mountains(self.domain,eps=EPS)
426 mts.setVelocity(v)
427 Z=mts.update()
428
429 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[2]-1.)))
430 self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
431
432 class Test_Mountains2D(unittest.TestCase):
433 def setUp(self):
434 NE=16
435 self.TOL=1e-4
436 self.domain=Rectangle(NE,NE,order=1,useFullElementOrder=True)
437 def tearDown(self):
438 del self.domain
439 def test_periodic(self):
440 EPS=0.01
441
442 x=self.domain.getX()
443 v = Vector(0.0, Solution(self.domain))
444 a0=1
445 n0=1
446 n1=0.5
447 a1=-(a0*n0)/n1
448 v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])
449 v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])
450
451 H_t=Scalar(0.0, Solution(self.domain))
452 mts=Mountains(self.domain,eps=EPS)
453 mts.setVelocity(v)
454 Z=mts.update()
455
456 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[1]-1.)))
457 self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
458
459
460
461 class Test_Rheologies(unittest.TestCase):
462 """
463 this is the program used to generate the powerlaw tests:
464
465 TAU_Y=100.
466 N=10
467 M=5
468
469 def getE(tau):
470 if tau<=TAU_Y:
471 return 1./(0.5+20*sqrt(tau))
472 else:
473 raise ValueError,"out of range."
474 tau=[ i* (TAU_Y/N) for i in range(N+1)] + [TAU_Y for i in range(M)]
475 e=[ tau[i]/getE(tau[i]) for i in range(N+1)]
476 e+= [ (i+1+j)* (max(e)/N) for j in range(M) ]
477
478 print tau
479 print e
480 """
481 TOL=1.e-8
482 def test_PowerLaw_Init(self):
483 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
484
485 self.failUnlessEqual(pl.getNumMaterials(),3,"num materials is wrong")
486 self.failUnlessEqual(pl.validMaterialId(0),True,"material id 0 not found")
487 self.failUnlessEqual(pl.validMaterialId(1),True,"material id 1 not found")
488 self.failUnlessEqual(pl.validMaterialId(2),True,"material id 2 not found")
489 self.failUnlessEqual(pl.validMaterialId(3),False,"material id 3 not found")
490
491 self.failUnlessEqual(pl.getFriction(),None,"initial friction wrong.")
492 self.failUnlessEqual(pl.getTauY(),None,"initial tau y wrong.")
493 pl.setDruckerPragerLaw(tau_Y=10,friction=3)
494 self.failUnlessEqual(pl.getFriction(),3,"friction wrong.")
495 self.failUnlessEqual(pl.getTauY(),10,"tau y wrong.")
496
497 self.failUnlessEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
498 pl.setElasticShearModulus(1000)
499 self.failUnlessEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
500
501 e=pl.getEtaN()
502 self.failUnlessEqual(len(e),3,"initial length of etaN is wrong.")
503 self.failUnlessEqual(e,[None, None, None],"initial etaN are wrong.")
504 self.failUnlessEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
505 self.failUnlessEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
506 self.failUnlessEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
507 self.failUnlessRaises(ValueError, pl.getEtaN, 3)
508
509 self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
510 self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
511 self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
512
513 pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
514 self.failUnlessEqual(pl.getPower(),[1,2,3],"powers are wrong.")
515 self.failUnlessEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
516 self.failUnlessEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
517
518 pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
519 self.failUnlessEqual(pl.getPower(),[4,2,3],"powers are wrong.")
520 self.failUnlessEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
521 self.failUnlessEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
522
523 self.failUnlessRaises(ValueError,pl.getPower,-1)
524 self.failUnlessEqual(pl.getPower(0),4,"power 0 is wrong.")
525 self.failUnlessEqual(pl.getPower(1),2,"power 1 is wrong.")
526 self.failUnlessEqual(pl.getPower(2),3,"power 2 is wrong.")
527 self.failUnlessRaises(ValueError,pl.getPower,3)
528
529 self.failUnlessRaises(ValueError,pl.getTauT,-1)
530 self.failUnlessEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
531 self.failUnlessEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
532 self.failUnlessEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
533 self.failUnlessRaises(ValueError,pl.getTauT,3)
534
535 self.failUnlessRaises(ValueError,pl.getEtaN,-1)
536 self.failUnlessEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
537 self.failUnlessEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
538 self.failUnlessEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
539 self.failUnlessRaises(ValueError,pl.getEtaN,3)
540
541 def checkResult(self,id,gamma_dot_, eta, tau_ref):
542 self.failUnless(eta>=0,"eta needs to be positive (test %s)"%id)
543 error=abs(gamma_dot_*eta-tau_ref)
544 self.failUnless(error<=self.TOL*tau_ref,"eta is wrong: error = gamma_dot_*eta-tau_ref = %s * %s - %s = %s (test %s)"%(gamma_dot_,eta,tau_ref,error,id))
545
546 def test_PowerLaw_Linear(self):
547 taus= [0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
548 gamma_dot_s=[0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0, 75.0]
549 pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
550 pl.setDruckerPragerLaw(tau_Y=100.)
551 pl.setPowerLaw(eta_N=2.)
552 pl.setEtaTolerance(self.TOL)
553 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
554
555 def test_PowerLaw_QuadLarge(self):
556 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
557 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 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
559 pl.setDruckerPragerLaw(tau_Y=100.)
560 pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2])
561 pl.setEtaTolerance(self.TOL)
562 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
563
564 def test_PowerLaw_QuadSmall(self):
565 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
566 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 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
568 pl.setDruckerPragerLaw(tau_Y=100.)
569 pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2])
570 pl.setEtaTolerance(self.TOL)
571 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
572
573 def test_PowerLaw_CubeLarge(self):
574 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
575 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 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
577 pl.setDruckerPragerLaw(tau_Y=100.)
578 pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3])
579 pl.setEtaTolerance(self.TOL)
580 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
581
582 def test_PowerLaw_CubeSmall(self):
583 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
584 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 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
586 pl.setDruckerPragerLaw(tau_Y=100.)
587 pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3])
588 pl.setEtaTolerance(self.TOL)
589 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
590
591 def test_PowerLaw_QuadLarge_CubeLarge(self):
592 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
593 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 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
596 pl.setDruckerPragerLaw(tau_Y=100.)
597 pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
598 pl.setEtaTolerance(self.TOL)
599 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
600
601 def test_PowerLaw_QuadLarge_CubeSmall(self):
602 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
603 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 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
606 pl.setDruckerPragerLaw(tau_Y=100.)
607 pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
608 pl.setEtaTolerance(self.TOL)
609 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
610
611 def test_PowerLaw_QuadSmall_CubeLarge(self):
612 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
613 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 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
615 pl.setDruckerPragerLaw(tau_Y=100.)
616 pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
617 pl.setEtaTolerance(self.TOL)
618 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
619
620 def test_PowerLaw_QuadSmall_CubeSmall(self):
621 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
622 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 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
624 pl.setDruckerPragerLaw(tau_Y=100.)
625 pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
626 pl.setEtaTolerance(self.TOL)
627 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
628
629 def test_PowerLaw_withShear(self):
630 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
631 gamma_dot_s=[0.0, 15.0, 30.0, 45.0, 60.0, 75.0, 90.0, 105.0, 120.0, 135.0, 150.0, 165.0, 180.0, 195.0, 210.0, 225.0]
632 pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
633 pl.setDruckerPragerLaw(tau_Y=100.)
634 pl.setPowerLaw(eta_N=2.)
635 pl.setElasticShearModulus(3.)
636 dt=1./3.
637 pl.setEtaTolerance(self.TOL)
638 self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
639 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
640
641 class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase):
642 TOL=1.e-5
643 VERBOSE=False # or True
644 A=1.
645 P_max=100
646 NE=2*getMPISizeWorld()
647 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 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
667 alpha=0.5*inv_eta*s_00
668 if s_00 <= B and self.mu !=None: alpha+=1./(2*self.mu)*self.A
669 u_ref=x*alpha
670 u_ref[self.dom.getDim()-1]=(1.-x[self.dom.getDim()-1])*alpha*(self.dom.getDim()-1)
671 sigma_ref=kronecker(self.dom)*s_00
672 sigma_ref[self.dom.getDim()-1,self.dom.getDim()-1]=-s_00*(self.dom.getDim()-1)
673
674 p_ref=self.P_max
675 for d in range(self.dom.getDim()): p_ref=p_ref*x[d]
676 p_ref-=integrate(p_ref)/vol(self.dom)
677 return u_ref, sigma_ref, p_ref
678
679 def runIt(self, free=None):
680 x=self.dom.getX()
681 B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
682 dt=B/int(self.N_dt/2)
683 if self.VERBOSE: print "dt =",dt
684 if self.latestart:
685 t=dt
686 else:
687 t=0
688 v,s,p=self.getReference(t)
689
690 mod=IncompressibleIsotropicFlowCartesian(self.dom, stress=s, v=v, p=p, t=t, numMaterials=3, verbose=self.VERBOSE)
691 mod.setDruckerPragerLaw(tau_Y=self.tau_Y,friction=None)
692 mod.setElasticShearModulus(self.mu)
693 mod.setPowerLaws([self.eta_0, self.eta_1, self.eta_2], [ 1., self.tau_1, self.tau_2], [1.,self.N_1,self.N_2])
694 mod.setTolerance(self.TOL)
695 mod.setEtaTolerance(self.TOL*0.1)
696
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 # 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 self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)
723 t+=dt
724 N_t+=1
725
726 def check(self,N_t,mod,t_ref,v_ref, s_ref,p_ref):
727 p=mod.getPressure()
728 p-=integrate(p)/vol(self.dom)
729 error_p=Lsup(mod.getPressure()-p_ref)/Lsup(p_ref)
730 error_s=Lsup(mod.getDeviatoricStress()-s_ref)/Lsup(s_ref)
731 error_v=Lsup(mod.getVelocity()-v_ref)/Lsup(v_ref)
732 error_t=abs(mod.getTime()-t_ref)/abs(t_ref)
733 if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v
734 self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
735 self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )
736 self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )
737 self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )
738 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
824 class Test_FaultSystem(unittest.TestCase):
825 EPS=1.e-8
826 NE=10
827 def test_Fault_MaxValue(self):
828 dom=Rectangle(2*self.NE,2*self.NE)
829 x=dom.getX()
830 f=FaultSystem(dim=2)
831 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
834 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 self.failUnless( m == 0.25, "wrong max value")
839 self.failUnless( t == 1, "wrong max tag")
840 self.failUnless( l == 0., "wrong max location")
841
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 self.failUnless( m == 0.0625, "wrong max value")
847 self.failUnless( t == 2, "wrong max tag")
848 self.failUnless( l == 0.5, "wrong max location")
849
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 self.failUnless( m == 0.25, "wrong max value")
855 self.failUnless( t == 2, "wrong max tag")
856 self.failUnless( l == 1.0, "wrong max location")
857
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 self.failUnless( m == 0.25, "wrong max value")
863 self.failUnless( t == 2, "wrong max tag")
864 self.failUnless( l == 0., "wrong max location")
865
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 self.failUnless( m == 0.25, "wrong max value")
871 self.failUnless( t == 1, "wrong max tag")
872 self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong max location")
873 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
880 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 self.failUnless( m == -0.25, "wrong min value")
885 self.failUnless( t == 1, "wrong min tag")
886 self.failUnless( l == 0., "wrong min location")
887 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 self.failUnless( m == -0.0625, "wrong min value")
892 self.failUnless( t == 2, "wrong min tag")
893 self.failUnless( l == 0.5, "wrong min location")
894 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 self.failUnless( m == -0.25, "wrong min value")
899 self.failUnless( t == 2, "wrong min tag")
900 self.failUnless( l == 1.0, "wrong min location")
901 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 self.failUnless( m == -0.25, "wrong min value")
906 self.failUnless( t == 2, "wrong min tag")
907 self.failUnless( l == 0., "wrong min location")
908 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 self.failUnless( m == -0.25, "wrong min value")
913 self.failUnless( t == 1, "wrong min tag")
914 self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong min location")
915
916
917 def test_Fault2D(self):
918 f=FaultSystem(dim=2)
919 top1=[ [1.,0.], [1.,1.], [0.,1.] ]
920 self.failUnlessRaises(ValueError,f.addFault,V0=[1.,0],strikes=[pi/2, pi/2],ls=[1.,1.],tag=1,dips=top1)
921 f.addFault(V0=[1.,0],strikes=[pi/2, pi],ls=[1.,1.],tag=1)
922 self.failUnless(f.getDim() == 2, "wrong dimension")
923 self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
924 self.failUnless( 2. == f.getTotalLength(1), "length wrong")
925 self.failUnless( 0. == f.getMediumDepth(1), "depth wrong")
926 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 range")
927 self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 range")
928 self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsets")
929 segs=f.getTopPolyline(1)
930 self.failUnless( len(segs) == 3, "wrong number of segments")
931 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
932 self.failUnless( segs[0].size == 2, "seg 0 has wrong size.")
933 self.failUnless( numpy.linalg.norm(segs[0]-[1.,0.]) < self.EPS, "wrong vertex. 0 ")
934 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
935 self.failUnless( segs[1].size == 2, "seg 1 has wrong size.")
936 self.failUnless( numpy.linalg.norm(segs[1]-[1.,1.]) < self.EPS, "wrong vertex. 1 ")
937 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
938 self.failUnless( segs[2].size == 2, "seg 2 has wrong size.")
939 self.failUnless( numpy.linalg.norm(segs[2]-[0.,1.]) < self.EPS, "wrong vertex. 2 ")
940 c=f.getCenterOnSurface()
941 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
942 self.failUnless( c.size == 2, "center size is wrong")
943 self.failUnless( numpy.linalg.norm(c-[2./3.,2./3.]) < self.EPS, "center has wrong coordinates.")
944 o=f.getOrientationOnSurface()/pi*180.
945 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
946
947 top2=[ [10.,0.], [0.,10.] ]
948 f.addFault(V0=[10.,0],strikes=[3.*pi/4],ls=[14.142135623730951], tag=2, w0_offsets=[0,20], w1_max=20)
949 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
950 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length")
951 self.failUnless( 0. == f.getMediumDepth(2), "depth wrong")
952 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range")
953 self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range")
954 self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets")
955 segs=f.getTopPolyline(2)
956 self.failUnless( len(segs) == 2, "wrong number of segments")
957 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
958 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.]) < self.EPS, "wrong vertex. 0 ")
959 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
960 self.failUnless( numpy.linalg.norm(segs[1]-[0.,10.]) < self.EPS, "wrong vertex. 1 ")
961 c=f.getCenterOnSurface()
962 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
963 self.failUnless( c.size == 2, "center size is wrong")
964 self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.]) < self.EPS, "center has wrong coordinates.")
965 o=f.getOrientationOnSurface()/pi*180.
966 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
967
968 s,d=f.getSideAndDistance([0.,0.], tag=1)
969 self.failUnless( s<0, "wrong side.")
970 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
971 s,d=f.getSideAndDistance([0.,2.], tag=1)
972 self.failUnless( s>0, "wrong side.")
973 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
974 s,d=f.getSideAndDistance([1.,2.], tag=1)
975 self.failUnless( s>0, "wrong side.")
976 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
977 s,d=f.getSideAndDistance([2.,1.], tag=1)
978 self.failUnless( s>0, "wrong side.")
979 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
980 s,d=f.getSideAndDistance([2.,0.], tag=1)
981 self.failUnless( s>0, "wrong side.")
982 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
983 s,d=f.getSideAndDistance([0.,-1.], tag=1)
984 self.failUnless( s<0, "wrong side.")
985 self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
986 s,d=f.getSideAndDistance([-1.,0], tag=1)
987 self.failUnless( s<0, "wrong side.")
988 self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
989
990
991 f.transform(rot=-pi/2., shift=[-1.,-1.])
992 self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
993 self.failUnless( 2. == f.getTotalLength(1), "length after transformation wrong")
994 self.failUnless( 0. == f.getMediumDepth(1), "depth after transformation wrong")
995 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 after transformation range")
996 self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
997 self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsetsafter transformation ")
998 segs=f.getTopPolyline(1)
999 self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
1000 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1001 self.failUnless( segs[0].size == 2, "seg 0 has wrong size after transformation.")
1002 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1003 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1004 self.failUnless( segs[1].size == 2, "seg 1 has wrong size after transformation.")
1005 self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1006 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1007 self.failUnless( segs[2].size == 2, "seg 2 has wrong size after transformation.")
1008 self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.]) < self.EPS, "wrong vertex after transformation. 2 ")
1009 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1010 self.failUnless( 0. == f.getMediumDepth(2), "depth wrong after transformation")
1011 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range after transformation")
1012 self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1013 self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1014 segs=f.getTopPolyline(2)
1015 self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
1016 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1017 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9]) < self.EPS, "wrong vertex. 0 after transformation")
1018 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1019 self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.]) < self.EPS, "wrong vertex. 1 after transformation")
1020
1021 c=f.getCenterOnSurface()
1022 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1023 self.failUnless( c.size == 2, "center size is wrong")
1024 self.failUnless( numpy.linalg.norm(c-[7./5.,-7./5.]) < self.EPS, "center has wrong coordinates.")
1025 o=f.getOrientationOnSurface()/pi*180.
1026 self.failUnless( abs(o-45.) < self.EPS, "wrong orientation.")
1027
1028 p=f.getParametrization([-1.,0.],1)
1029 self.failUnless(p[1]==1., "wrong value.")
1030 self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1031 p=f.getParametrization([-0.5,0.],1)
1032 self.failUnless(p[1]==1., "wrong value.")
1033 self.failUnless(abs(p[0]-0.5)<self.EPS* 0.5, "wrong value.")
1034 p=f.getParametrization([0.,0.],1)
1035 self.failUnless(p[1]==1., "wrong value.")
1036 self.failUnless(abs(p[0]-1.)<self.EPS, "wrong value.")
1037 p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-8)
1038 self.failUnless(p[1]==0., "wrong value.")
1039 p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-6)
1040 self.failUnless(p[1]==1., "wrong value.")
1041 self.failUnless(abs(p[0]-1.0000001)<self.EPS, "wrong value.")
1042 p=f.getParametrization([0.,0.5],1)
1043 self.failUnless(p[1]==1., "wrong value.")
1044 self.failUnless(abs(p[0]-1.5)<self.EPS, "wrong value.")
1045 p=f.getParametrization([0,1.],1)
1046 self.failUnless(p[1]==1., "wrong value.")
1047 self.failUnless(abs(p[0]-2.)<self.EPS, "wrong value.")
1048 p=f.getParametrization([1.,1.],1)
1049 self.failUnless(p[1]==0., "wrong value.")
1050 p=f.getParametrization([0,1.11],1)
1051 self.failUnless(p[1]==0., "wrong value.")
1052 p=f.getParametrization([-1,-9.],2)
1053 self.failUnless(p[1]==1., "wrong value.")
1054 self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1055 p=f.getParametrization([9,1],2)
1056 self.failUnless(p[1]==1., "wrong value.")
1057 self.failUnless(abs(p[0]-20.)<self.EPS, "wrong value.")
1058
1059 def test_Fault3D(self):
1060 f=FaultSystem(dim=3)
1061 self.failUnless(f.getDim() == 3, "wrong dimension")
1062
1063 top1=[ [0.,0.,0.], [1., 0., 0.] ]
1064 f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=20.,tag=1)
1065 self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
1066 self.failUnless( 1. == f.getTotalLength(1), "length wrong")
1067 self.failUnless( 20. == f.getMediumDepth(1), "depth wrong")
1068 self.failUnless( (0., 1.) == f.getW0Range(1)," wrong W0 range")
1069 self.failUnless( (-20., 0.) == f.getW1Range(1)," wrong W1 range")
1070 self.failUnless( [0., 1.] == f.getW0Offsets(1)," wrong W0 offsets")
1071 segs=f.getTopPolyline(1)
1072 self.failUnless( len(segs) == 2, "wrong number of segments")
1073 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1074 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1075 self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1076 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1077 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1078 self.failUnless( numpy.linalg.norm(segs[1]-[1.,0.,0]) < self.EPS, "wrong vertex. 1 ")
1079 c=f.getCenterOnSurface()
1080 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1081 self.failUnless( c.size == 3, "center size is wrong")
1082 self.failUnless( numpy.linalg.norm(c-[0.5,0.,0.]) < self.EPS, "center has wrong coordinates.")
1083 o=f.getOrientationOnSurface()/pi*180.
1084 self.failUnless( abs(o) < self.EPS, "wrong orientation.")
1085 d=f.getDips(1)
1086 self.failUnless( len(d) == 1, "wrong number of dips")
1087 self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1088 sn=f.getSegmentNormals(1)
1089 self.failUnless( len(sn) == 1, "wrong number of normals")
1090 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1091 self.failUnless( numpy.linalg.norm(sn[0]-[0, -1., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1092 dv=f.getDepthVectors(1)
1093 self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1094 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1095 self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1096 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1097 self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1098 b=f.getBottomPolyline(1)
1099 self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1100 self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1101 self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1102 self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1103 self.failUnless( numpy.linalg.norm(b[1]-[1., 0., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1104 ds=f.getDepths(1)
1105 self.failUnless( len(ds) == 2, "wrong number of depth")
1106 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1107 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1108
1109 top2=[ [0.,0.,0.], [0., 10., 0.] ]
1110 f.addFault(V0=[0.,0,0],strikes=[pi/2],ls=[10.,], dips=pi/2, depths=20.,tag=2)
1111 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1112 self.failUnless( 10. == f.getTotalLength(2), "length wrong")
1113 self.failUnless( 20. == f.getMediumDepth(2), "depth wrong")
1114 self.failUnless( (0., 10.) == f.getW0Range(2)," wrong W0 range")
1115 self.failUnless( (-20., 0.) == f.getW1Range(2)," wrong W1 range")
1116 self.failUnless( [0., 10.] == f.getW0Offsets(2)," wrong W0 offsets")
1117 segs=f.getTopPolyline(2)
1118 self.failUnless( len(segs) == 2, "wrong number of segments")
1119 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1120 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1121 self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1122 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1123 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1124 self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1125 d=f.getDips(2)
1126 self.failUnless( len(d) == 1, "wrong number of dips")
1127 self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1128 sn=f.getSegmentNormals(2)
1129 self.failUnless( len(sn) == 1, "wrong number of normals")
1130 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1131 self.failUnless( numpy.linalg.norm(sn[0]-[1, 0., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1132 dv=f.getDepthVectors(2)
1133 self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1134 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1135 self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1136 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1137 self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1138 b=f.getBottomPolyline(2)
1139 self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1140 self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1141 self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1142 self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1143 self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1144 ds=f.getDepths(2)
1145 self.failUnless( len(ds) == 2, "wrong number of depth")
1146 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1147 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1148
1149 top2=[ [10.,0.,0.], [0., 10., 0.] ]
1150 f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/2, depths=30.,tag=2)
1151 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1152 self.failUnless( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1153 self.failUnless( 30. == f.getMediumDepth(2), "depth wrong")
1154 self.failUnless( (-30., 0.) == f.getW1Range(2)," wrong W1 range")
1155 segs=f.getTopPolyline(2)
1156 self.failUnless( len(segs) == 2, "wrong number of segments")
1157 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1158 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1159 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1160 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1161 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1162 self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1163 d=f.getDips(2)
1164 self.failUnless( len(d) == 1, "wrong number of dips")
1165 self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1166 sn=f.getSegmentNormals(2)
1167 self.failUnless( len(sn) == 1, "wrong number of normals")
1168 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1169 self.failUnless( numpy.linalg.norm(sn[0]-[0.70710678118654746, 0.70710678118654746, 0.]) < self.EPS, "wrong bottom vertex 1 ")
1170 dv=f.getDepthVectors(2)
1171 self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1172 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1173 self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -30.]) < self.EPS, "wrong depth vector 0 ")
1174 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1175 self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -30.]) < self.EPS, "wrong depth vector 1 ")
1176 b=f.getBottomPolyline(2)
1177 self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1178 self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1179 self.failUnless( numpy.linalg.norm(b[0]-[10., 0., -30.]) < self.EPS, "wrong bottom vertex 0 ")
1180 self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1181 self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -30.]) < self.EPS, "wrong bottom vertex 1 ")
1182 ds=f.getDepths(2)
1183 self.failUnless( len(ds) == 2, "wrong number of depth")
1184 self.failUnless( abs(ds[0]-30.) < self.EPS, "wrong depth at vertex 0 ")
1185 self.failUnless( abs(ds[1]-30.) < self.EPS, "wrong depth at vertex 1 ")
1186
1187 top2=[ [10.,0.,0.], [0., 10., 0.] ]
1188 f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/4, depths=50.,tag=2)
1189 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1190 self.failUnless( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1191 self.failUnless( 50. == f.getMediumDepth(2), "depth wrong")
1192 self.failUnless( (-50., 0.) == f.getW1Range(2)," wrong W1 range")
1193 segs=f.getTopPolyline(2)
1194 self.failUnless( len(segs) == 2, "wrong number of segments")
1195 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1196 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1197 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1198 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1199 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1200 self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1201 d=f.getDips(2)
1202 self.failUnless( len(d) == 1, "wrong number of dips")
1203 self.failUnless( abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1204 sn=f.getSegmentNormals(2)
1205 self.failUnless( len(sn) == 1, "wrong number of normals")
1206 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1207 self.failUnless( numpy.linalg.norm(sn[0]-[0.5,0.5,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1208 dv=f.getDepthVectors(2)
1209 self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1210 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1211 self.failUnless( numpy.linalg.norm(dv[0]-[25., 25., -35.355339059327378]) < self.EPS, "wrong depth vector 0 ")
1212 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1213 self.failUnless( numpy.linalg.norm(dv[1]-[25.,25., -35.355339059327378]) < self.EPS, "wrong depth vector 1 ")
1214 b=f.getBottomPolyline(2)
1215 self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1216 self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1217 self.failUnless( numpy.linalg.norm(b[0]-[35., 25., -35.355339059327378]) < self.EPS, "wrong bottom vertex 0 ")
1218 self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1219 self.failUnless( numpy.linalg.norm(b[1]-[25, 35., -35.355339059327378]) < self.EPS, "wrong bottom vertex 1 ")
1220 ds=f.getDepths(2)
1221 self.failUnless( len(ds) == 2, "wrong number of depth")
1222 self.failUnless( abs(ds[0]-50.) < self.EPS, "wrong depth at vertex 0 ")
1223 self.failUnless( abs(ds[1]-50.) < self.EPS, "wrong depth at vertex 1 ")
1224
1225 top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
1226 f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=1, dips=pi/4, depths=20.)
1227 self.failUnless( 20. == f.getTotalLength(1), "length wrong")
1228 self.failUnless( 20. == f.getMediumDepth(1), "depth wrong")
1229 segs=f.getTopPolyline(1)
1230 self.failUnless( len(segs) == 3, "wrong number of segments")
1231 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1232 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1233 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1234 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1235 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1236 self.failUnless( numpy.linalg.norm(segs[1]-[10.,10.,0.]) < self.EPS, "wrong vertex. 1 ")
1237 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1238 self.failUnless( segs[2].size == 3, "seg 2 has wrong size.")
1239 self.failUnless( numpy.linalg.norm(segs[2]-[0.,10.,0.]) < self.EPS, "wrong vertex. 2 ")
1240 d=f.getDips(1)
1241 self.failUnless( len(d) == 2, "wrong number of dips")
1242 self.failUnless( abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1243 self.failUnless( abs(d[1]-0.78539816339744828) < self.EPS, "wrong dip 0")
1244 ds=f.getDepths(1)
1245 self.failUnless( len(ds) == 3, "wrong number of depth")
1246 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1247 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1248 sn=f.getSegmentNormals(1)
1249 self.failUnless( len(sn) == 2, "wrong number of normals")
1250 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1251 self.failUnless( numpy.linalg.norm(sn[0]-[0.70710678118654746,0.,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1252 self.failUnless( isinstance(sn[1], numpy.ndarray), "wrong class of bottom vertex 0")
1253 self.failUnless( numpy.linalg.norm(sn[1]-[0.,0.70710678118654746,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1254 dv=f.getDepthVectors(1)
1255 self.failUnless( len(dv) == 3, "wrong number of depth vectors.")
1256 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1257 self.failUnless( numpy.linalg.norm(dv[0]-[14.142135623730951, 0., -14.142135623730951]) < self.EPS, "wrong depth vector 0 ")
1258 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1259 self.failUnless( numpy.linalg.norm(dv[1]-[11.547005383792515,11.547005383792515, -11.547005383792515]) < self.EPS, "wrong depth vector 2 ")
1260 self.failUnless( isinstance(dv[2], numpy.ndarray), "wrong class of depth vector 1")
1261 self.failUnless( numpy.linalg.norm(dv[2]-[0.,14.142135623730951, -14.142135623730951]) < self.EPS, "wrong depth vector 2 ")
1262 segs=f.getBottomPolyline(1)
1263 self.failUnless( len(segs) == 3, "wrong number of segments")
1264 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1265 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1266 self.failUnless( numpy.linalg.norm(segs[0]-[24.142135623730951,0.,-14.142135623730951]) < self.EPS, "wrong vertex. 0 ")
1267 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1268 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1269 self.failUnless( numpy.linalg.norm(segs[1]-[21.547005383792515,21.547005383792515, -11.547005383792515]) < self.EPS, "wrong vertex. 1 ")
1270 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1271 self.failUnless( segs[2].size == 3, "seg 2 has wrong size.")
1272 self.failUnless( numpy.linalg.norm(segs[2]-[0., 24.142135623730951, -14.142135623730951]) < self.EPS, "wrong vertex. 2 ")
1273 self.failUnless( abs(0.-f.getW0Range(1)[0]) <=self.EPS," wrong W0 range (0)")
1274 self.failUnless( abs(31.857329272664341-f.getW0Range(1)[1]) <=self.EPS," wrong W0 range (1)")
1275 self.failUnless( abs(-20. - f.getW1Range(1)[0]) <=self.EPS," wrong W1 range (0)")
1276 self.failUnless( abs(0. - f.getW1Range(1)[1]) <=self.EPS," wrong W1 range (1)")
1277 self.failUnless( abs(0.0-f.getW0Offsets(1)[0])<=self.EPS," wrong W0 offsets (0)")
1278 self.failUnless( abs(15.92866463633217-f.getW0Offsets(1)[1])<=self.EPS," wrong W0 offsets (1)")
1279 self.failUnless( abs(31.857329272664341-f.getW0Offsets(1)[2])<=self.EPS," wrong W0 offsets(2)")
1280 #
1281 # ============ fresh start ====================
1282 #
1283 f.addFault(V0=[1.,0,0.],strikes=[pi/2, pi],ls=[1.,1.], dips=pi/2,depths=20,tag=1)
1284 f.addFault(V0=[10.,0,0],strikes=[3.*pi/4],ls=[14.142135623730951], tag=2, w0_offsets=[0,20], w1_max=20, dips=pi/2,depths=20)
1285 c=f.getCenterOnSurface()
1286 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1287 self.failUnless( c.size == 3, "center size is wrong")
1288 self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.,0.]) < self.EPS, "center has wrong coordinates.")
1289 o=f.getOrientationOnSurface()/pi*180.
1290 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1291
1292 f.transform(rot=-pi/2., shift=[-1.,-1.,0.])
1293 self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
1294 self.failUnless( 2. == f.getTotalLength(1), "length after transformation wrong")
1295 self.failUnless( 20. == f.getMediumDepth(1), "depth after transformation wrong")
1296 rw0=f.getW0Range(1)
1297 self.failUnless( len(rw0) ==2, "wo range has wrong length")
1298 self.failUnless( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
1299 self.failUnless( abs(rw0[1]-2.) < self.EPS,"W0 1 wrong.")
1300 self.failUnless( (-20., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
1301 dips=f.getDips(1)
1302 self.failUnless(len(dips) == 2, "wrong number of dips.")
1303 self.failUnless( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
1304 self.failUnless( abs(dips[1]-1.5707963267948966) <= self.EPS, "wrong dip")
1305 ds=f.getDepths(1)
1306 self.failUnless( len(ds) == 3, "wrong number of depth")
1307 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1308 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1309 self.failUnless( abs(ds[2]-20.) < self.EPS, "wrong depth at vertex 1 ")
1310 segs=f.getTopPolyline(1)
1311 self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
1312 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1313 self.failUnless( segs[0].size == 3, "seg 0 has wrong size after transformation.")
1314 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1315 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1316 self.failUnless( segs[1].size == 3, "seg 1 has wrong size after transformation.")
1317 self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1318 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1319 self.failUnless( segs[2].size == 3, "seg 2 has wrong size after transformation.")
1320 self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.,0.]) < self.EPS, "wrong vertex after transformation. 2 ")
1321 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1322 self.failUnless( 20. == f.getMediumDepth(2), "depth wrong after transformation")
1323 rw0=f.getW0Range(2)
1324 self.failUnless( len(rw0) ==2, "wo range has wrong length")
1325 self.failUnless( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
1326 self.failUnless( abs(rw0[1]-20.) < self.EPS,"W0 1 wrong.")
1327 self.failUnless( (-20., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1328 self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1329 dips=f.getDips(2)
1330 self.failUnless(len(dips) == 1, "wrong number of dips.")
1331 self.failUnless( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
1332 ds=f.getDepths(2)
1333 self.failUnless( len(ds) == 2, "wrong number of depth")
1334 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1335 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1336 segs=f.getTopPolyline(2)
1337 self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
1338 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1339 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1340 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1341 self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.,0.]) < self.EPS, "wrong vertex. 1 after transformation")
1342 #
1343 # ============ fresh start ====================
1344 #
1345 f=FaultSystem(dim=3)
1346
1347 top1=[ [0.,0.,0.], [1., 0., 0.] ]
1348 f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=1.,tag=1)
1349 top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
1350 f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=2, dips=pi/4, depths=20.)
1351
1352 p,m=f.getParametrization([0.3,0.,-0.5],1)
1353 self.failUnless(length(p-[0.3,-0.5]) <= self.EPS, "wrong value.")
1354 self.failUnless(m==1., "wrong value.")
1355
1356 p,m=f.getParametrization([0.5,0.,-0.5],1)
1357 self.failUnless(length(p-[0.5,-0.5]) <= self.EPS, "wrong value.")
1358 self.failUnless(m==1., "wrong value.")
1359
1360 p,m=f.getParametrization([0.25,0.,-0.5],1)
1361 self.failUnless(length(p-[0.25,-0.5]) <= self.EPS, "wrong value.")
1362 self.failUnless(m==1., "wrong value.")
1363
1364 p,m=f.getParametrization([0.5,0.,-0.25],1)
1365 self.failUnless(length(p-[0.5,-0.25]) <= self.EPS, "wrong value.")
1366 self.failUnless(m==1., "wrong value.")
1367
1368 p,m=f.getParametrization([0.001,0.,-0.001],1)
1369 self.failUnless(length(p-[0.001, -0.001]) <= self.EPS, "wrong value.")
1370 self.failUnless(m==1., "wrong value.")
1371
1372 p,m=f.getParametrization([0.001,0.,0.001],1)
1373 self.failUnless(m==0., "wrong value.")
1374
1375 p,m=f.getParametrization([0.999,0.,0.001],1)
1376 self.failUnless(m==0., "wrong value.")
1377
1378 p,m=f.getParametrization([1.001,0.,-0.001],1)
1379 self.failUnless(m==0., "wrong value.")
1380 p,m=f.getParametrization([1.001,0.,-0.1],1)
1381 self.failUnless(m==0., "wrong value.")
1382 p,m=f.getParametrization([1.001,0.,-0.000000001],1)
1383 self.failUnless(m==0., "wrong value.")
1384
1385 p,m=f.getParametrization([0.999,0.,-0.001],1)
1386 self.failUnless(length(p-[0.999, -0.001]) <= self.EPS, "wrong value.")
1387 self.failUnless(m==1., "wrong value.")
1388
1389 p,m=f.getParametrization([ 16.29252873 , 6.46410161 ,-6.29252873], 2, tol=1.e-7)
1390 self.failUnless(m==1., "wrong value.")
1391 self.failUnless(length(p-[4.7785993908996511, -10]) <= self.EPS*10., "wrong value.")
1392 p,m=f.getParametrization([15.77350269, 12.77350269, -5.77350269], 2, tol=1.e-7)
1393 self.failUnless(m==1., "wrong value.")
1394 self.failUnless(length(p-[11.150065245432518, -10]) <= self.EPS*10., "wrong value.")
1395
1396 p,m=f.getParametrization([ 3., 17.0710678, -7.0710678], 2, tol=1.e-7)
1397 self.failUnless(m==1., "wrong value.")
1398 self.failUnless(length(p-[27.078729881764687, -10]) <= self.EPS*10., "wrong value.")
1399 p,m=f.getParametrization([9.30940108, 16.55204176, -6.55204176], 2, tol=1.e-7)
1400 self.failUnless(m==1., "wrong value.")
1401 self.failUnless(length(p-[20.707264027231822, -10]) <= self.EPS*10., "wrong value.")
1402
1403 p,m=f.getParametrization([ 21.54700538, 21.54700538, -11.54700538], 2, tol=1.e-7)
1404 self.failUnless(m==1., "wrong value.")
1405 self.failUnless(length(p-[15.92866463633217, -20]) <= self.EPS*10., "wrong value.")
1406
1407 p,m=f.getParametrization([ 0.,0.,0.], 2, tol=1.e-7)
1408 self.failUnless(m==0., "wrong value.")
1409
1410 p,m=f.getParametrization([ 11.,11.,0.], 2, tol=1.e-7)
1411 self.failUnless(m==0., "wrong value.")
1412
1413
1414 s,d=f.getSideAndDistance([0.,-1.,0.], tag=1)
1415 self.failUnless( s>0, "wrong side.")
1416 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1417 s,d=f.getSideAndDistance([1.,-1.,0.], tag=1)
1418 self.failUnless( s>0, "wrong side.")
1419 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1420 s,d=f.getSideAndDistance([0.,1.,0.], tag=1)
1421 self.failUnless( s<0, "wrong side.")
1422 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1423 s,d=f.getSideAndDistance([1.,1.,0.], tag=1)
1424 self.failUnless( s<0, "wrong side.")
1425 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1426
1427
1428 s,d=f.getSideAndDistance([0.,0.,0.], tag=2)
1429 self.failUnless( s<0, "wrong side.")
1430 self.failUnless( abs(d-10.)<self.EPS, "wrong distance.")
1431 s,d=f.getSideAndDistance([5.,5.,0.], tag=2)
1432 self.failUnless( s<0, "wrong side.")
1433 self.failUnless( abs(d-5.)<self.EPS, "wrong distance.")
1434
1435 s,d=f.getSideAndDistance([10.,10.,-1.], tag=2)
1436 self.failUnless( s<0, "wrong side.")
1437 self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
1438 s,d=f.getSideAndDistance([10.,10.,-2.], tag=2)
1439 self.failUnless( s<0, "wrong side.")
1440 self.failUnless( abs(d-2.*0.70710678118654757)<self.EPS, "wrong distance.")
1441 s,d=f.getSideAndDistance([10.,10.,-3.], tag=2)
1442 self.failUnless( s<0, "wrong side.")
1443 self.failUnless( abs(d-3.*0.70710678118654757)<self.EPS, "wrong distance.")
1444
1445 s,d=f.getSideAndDistance([5.,12.,0], tag=2)
1446 self.failUnless( s>0, "wrong side.")
1447 self.failUnless( abs(d-2*0.70710678118654757)<self.EPS, "wrong distance.")
1448 s,d=f.getSideAndDistance([5.,12.,-1], tag=2)
1449 self.failUnless( s>0, "wrong side.")
1450 self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
1451 s,d=f.getSideAndDistance([5.,12.,-2], tag=2)
1452 # s not checked as it is undefined.
1453 self.failUnless( abs(d)<self.EPS, "wrong distance.")
1454 s,d=f.getSideAndDistance([5.,12.,-3], tag=2)
1455 self.failUnless( s<0, "wrong side.")
1456 self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
1457 s,d=f.getSideAndDistance([5.,12.,-4], tag=2)
1458 self.failUnless( s<0, "wrong side.")
1459 self.failUnless( abs(d-2.*0.70710678118654757)<self.EPS, "wrong distance.")
1460
1461 if __name__ == '__main__':
1462 suite = unittest.TestSuite()
1463 # suite.addTest(unittest.makeSuite(Test_FaultSystem))
1464 # suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
1465 suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
1466 suite.addTest(unittest.makeSuite(Test_Mountains3D))
1467 suite.addTest(unittest.makeSuite(Test_Mountains2D))
1468 suite.addTest(unittest.makeSuite(Test_Rheologies))
1469 suite.addTest(unittest.makeSuite(Test_IncompressibleIsotropicFlowCartesian))
1470 s=unittest.TextTestRunner(verbosity=2).run(suite)
1471 if not s.wasSuccessful(): sys.exit(1)

  ViewVC Help
Powered by ViewVC 1.1.26