/[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 2881 - (show annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 8 months ago) by jfenwick
File MIME type: text/x-python
File size: 86304 byte(s)
Don't panic.
Updating copyright stamps

1
2 ########################################################
3 #
4 # Copyright (c) 2003-2010 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 import unittest
23 import tempfile
24
25
26
27 VERBOSE=False # or True
28
29 from esys.escript import *
30 from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian, FaultSystem, DarcyFlow
31 from esys.escript.models import Mountains
32 from esys.finley import Rectangle, Brick
33
34 from math import pi
35 import numpy
36 import sys
37 import os
38 #====================================================================================================================
39 try:
40 FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
41 except KeyError:
42 FINLEY_WORKDIR='.'
43
44 #====================================================================================================================
45 class Test_StokesProblemCartesian2D(unittest.TestCase):
46 def setUp(self):
47 NE=6
48 self.TOL=1e-3
49 self.domain=Rectangle(NE,NE,order=-1)
50 def tearDown(self):
51 del self.domain
52 def test_PCG_P_0(self):
53 ETA=1.
54 P1=0.
55
56 x=self.domain.getX()
57 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
58 mask=whereZero(x[0]) * [1.,1.] \
59 +whereZero(x[0]-1) * [1.,1.] \
60 +whereZero(x[1]) * [1.,0.] \
61 +whereZero(x[1]-1) * [1.,1.]
62
63 sp=StokesProblemCartesian(self.domain)
64
65 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
66 u0=(1-x[0])*x[0]*[0.,1.]
67 p0=Scalar(-P1,ReducedSolution(self.domain))
68 sp.setTolerance(self.TOL)
69 u,p=sp.solve(u0*mask,p0,verbose=VERBOSE,max_iter=100,usePCG=True)
70
71 error_v0=Lsup(u[0]-u0[0])
72 error_v1=Lsup(u[1]-u0[1])/0.25
73 error_p=Lsup(p+P1*x[0]*x[1])
74 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
75 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
76 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
77
78 def test_PCG_P_small(self):
79 ETA=1.
80 P1=1.
81
82 x=self.domain.getX()
83 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
84 mask=whereZero(x[0]) * [1.,1.] \
85 +whereZero(x[0]-1) * [1.,1.] \
86 +whereZero(x[1]) * [1.,0.] \
87 +whereZero(x[1]-1) * [1.,1.]
88
89 sp=StokesProblemCartesian(self.domain)
90
91 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
92 u0=(1-x[0])*x[0]*[0.,1.]
93 p0=Scalar(-P1,ReducedSolution(self.domain))
94 sp.setTolerance(self.TOL)
95 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
96 error_v0=Lsup(u[0]-u0[0])
97 error_v1=Lsup(u[1]-u0[1])/0.25
98 error_p=Lsup(P1*x[0]*x[1]+p)
99 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
100 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
101 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
102
103 def test_PCG_P_large(self):
104 ETA=1.
105 P1=1000.
106
107 x=self.domain.getX()
108 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
109 mask=whereZero(x[0]) * [1.,1.] \
110 +whereZero(x[0]-1) * [1.,1.] \
111 +whereZero(x[1]) * [1.,0.] \
112 +whereZero(x[1]-1) * [1.,1.]
113
114 sp=StokesProblemCartesian(self.domain)
115
116 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
117 u0=(1-x[0])*x[0]*[0.,1.]
118 p0=Scalar(-P1,ReducedSolution(self.domain))
119 sp.setTolerance(self.TOL)
120 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
121 # u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
122
123 error_v0=Lsup(u[0]-u0[0])
124 error_v1=Lsup(u[1]-u0[1])/0.25
125 error_p=Lsup(P1*x[0]*x[1]+p)/P1
126 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
127 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
128 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
129
130 def test_GMRES_P_0(self):
131 ETA=1.
132 P1=0.
133
134 x=self.domain.getX()
135 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
136 mask=whereZero(x[0]) * [1.,1.] \
137 +whereZero(x[0]-1) * [1.,1.] \
138 +whereZero(x[1]) * [1.,0.] \
139 +whereZero(x[1]-1) * [1.,1.]
140
141 sp=StokesProblemCartesian(self.domain)
142
143 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
144 u0=(1-x[0])*x[0]*[0.,1.]
145 p0=Scalar(-P1,ReducedSolution(self.domain))
146 sp.setTolerance(self.TOL)
147 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=50,usePCG=False,iter_restart=18)
148
149 error_v0=Lsup(u[0]-u0[0])
150 error_v1=Lsup(u[1]-u0[1])/0.25
151 error_p=Lsup(P1*x[0]*x[1]+p)
152 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
153 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
154 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
155
156 def test_GMRES_P_small(self):
157 ETA=1.
158 P1=1.
159
160 x=self.domain.getX()
161 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
162 mask=whereZero(x[0]) * [1.,1.] \
163 +whereZero(x[0]-1) * [1.,1.] \
164 +whereZero(x[1]) * [1.,0.] \
165 +whereZero(x[1]-1) * [1.,1.]
166
167 sp=StokesProblemCartesian(self.domain)
168
169 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
170 u0=(1-x[0])*x[0]*[0.,1.]
171 p0=Scalar(-P1,ReducedSolution(self.domain))
172 sp.setTolerance(self.TOL)
173 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=20,usePCG=False)
174
175 error_v0=Lsup(u[0]-u0[0])
176 error_v1=Lsup(u[1]-u0[1])/0.25
177 error_p=Lsup(P1*x[0]*x[1]+p)
178
179 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
180 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
181 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
182
183 def test_GMRES_P_large(self):
184 ETA=1.
185 P1=1000.
186
187 x=self.domain.getX()
188 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
189 mask=whereZero(x[0]) * [1.,1.] \
190 +whereZero(x[0]-1) * [1.,1.] \
191 +whereZero(x[1]) * [1.,0.] \
192 +whereZero(x[1]-1) * [1.,1.]
193
194 sp=StokesProblemCartesian(self.domain)
195
196 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
197 u0=(1-x[0])*x[0]*[0.,1.]
198 p0=Scalar(-P1,ReducedSolution(self.domain))
199 sp.setTolerance(self.TOL)
200 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)
201
202 error_v0=Lsup(u[0]-u0[0])
203 error_v1=Lsup(u[1]-u0[1])/0.25
204 error_p=Lsup(P1*x[0]*x[1]+p)/P1
205 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
206 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
207 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
208 #====================================================================================================================
209 class Test_StokesProblemCartesian3D(unittest.TestCase):
210 def setUp(self):
211 NE=6
212 self.TOL=1e-4
213 self.domain=Brick(NE,NE,NE,order=-1)
214 def tearDown(self):
215 del self.domain
216 def test_PCG_P_0(self):
217 ETA=1.
218 P1=0.
219
220 x=self.domain.getX()
221 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.]
222 x=self.domain.getX()
223 mask=whereZero(x[0]) * [1.,1.,1.] \
224 +whereZero(x[0]-1) * [1.,1.,1.] \
225 +whereZero(x[1]) * [1.,0.,1.] \
226 +whereZero(x[1]-1) * [1.,1.,1.] \
227 +whereZero(x[2]) * [1.,1.,0.] \
228 +whereZero(x[2]-1) * [1.,1.,1.]
229
230
231 sp=StokesProblemCartesian(self.domain)
232
233 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
234 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
235 p0=Scalar(-P1,ReducedSolution(self.domain))
236 sp.setTolerance(self.TOL)
237 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
238
239 error_v0=Lsup(u[0]-u0[0])
240 error_v1=Lsup(u[1]-u0[1])
241 error_v2=Lsup(u[2]-u0[2])/0.25**2
242 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
243 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
244 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
245 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
246 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
247
248 def test_PCG_P_small(self):
249 ETA=1.
250 P1=1.
251
252 x=self.domain.getX()
253 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.]
254 mask=whereZero(x[0]) * [1.,1.,1.] \
255 +whereZero(x[0]-1) * [1.,1.,1.] \
256 +whereZero(x[1]) * [1.,0.,1.] \
257 +whereZero(x[1]-1) * [1.,1.,1.] \
258 +whereZero(x[2]) * [1.,1.,0.] \
259 +whereZero(x[2]-1) * [1.,1.,1.]
260
261
262 sp=StokesProblemCartesian(self.domain)
263
264 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
265 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
266 p0=Scalar(-P1,ReducedSolution(self.domain))
267 sp.setTolerance(self.TOL)
268 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
269 error_v0=Lsup(u[0]-u0[0])
270 error_v1=Lsup(u[1]-u0[1])
271 error_v2=Lsup(u[2]-u0[2])/0.25**2
272 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
273 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
274 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
275 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
276 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
277
278 def test_PCG_P_large(self):
279 ETA=1.
280 P1=1000.
281
282 x=self.domain.getX()
283 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.]
284 mask=whereZero(x[0]) * [1.,1.,1.] \
285 +whereZero(x[0]-1) * [1.,1.,1.] \
286 +whereZero(x[1]) * [1.,0.,1.] \
287 +whereZero(x[1]-1) * [1.,1.,1.] \
288 +whereZero(x[2]) * [1.,1.,0.] \
289 +whereZero(x[2]-1) * [1.,1.,1.]
290
291
292 sp=StokesProblemCartesian(self.domain)
293
294 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
295 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
296 p0=Scalar(-P1,ReducedSolution(self.domain))
297 sp.setTolerance(self.TOL)
298 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
299
300 error_v0=Lsup(u[0]-u0[0])
301 error_v1=Lsup(u[1]-u0[1])
302 error_v2=Lsup(u[2]-u0[2])/0.25**2
303 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
304 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
305 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
306 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
307 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
308
309 def test_GMRES_P_0(self):
310 ETA=1.
311 P1=0.
312
313 x=self.domain.getX()
314 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.]
315 x=self.domain.getX()
316 mask=whereZero(x[0]) * [1.,1.,1.] \
317 +whereZero(x[0]-1) * [1.,1.,1.] \
318 +whereZero(x[1]) * [1.,1.,1.] \
319 +whereZero(x[1]-1) * [1.,1.,1.] \
320 +whereZero(x[2]) * [1.,1.,0.] \
321 +whereZero(x[2]-1) * [1.,1.,1.]
322
323
324 sp=StokesProblemCartesian(self.domain)
325
326 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
327 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
328 p0=Scalar(-P1,ReducedSolution(self.domain))
329 sp.setTolerance(self.TOL)
330 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False,iter_restart=20)
331
332 error_v0=Lsup(u[0]-u0[0])
333 error_v1=Lsup(u[1]-u0[1])
334 error_v2=Lsup(u[2]-u0[2])/0.25**2
335 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
336 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
337 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
338 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
339 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
340 def test_GMRES_P_small(self):
341 ETA=1.
342 P1=1.
343
344 x=self.domain.getX()
345 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.]
346 mask=whereZero(x[0]) * [1.,1.,1.] \
347 +whereZero(x[0]-1) * [1.,1.,1.] \
348 +whereZero(x[1]) * [1.,1.,1.] \
349 +whereZero(x[1]-1) * [1.,1.,1.] \
350 +whereZero(x[2]) * [1.,1.,0.] \
351 +whereZero(x[2]-1) * [1.,1.,1.]
352
353
354 sp=StokesProblemCartesian(self.domain)
355
356 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
357 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
358 p0=Scalar(-P1,ReducedSolution(self.domain))
359 sp.setTolerance(self.TOL/10)
360 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)
361
362 error_v0=Lsup(u[0]-u0[0])
363 error_v1=Lsup(u[1]-u0[1])
364 error_v2=Lsup(u[2]-u0[2])/0.25**2
365 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
366 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
367 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
368 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
369 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
370 def test_GMRES_P_large(self):
371 ETA=1.
372 P1=1000.
373
374 x=self.domain.getX()
375 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.]
376 mask=whereZero(x[0]) * [1.,1.,1.] \
377 +whereZero(x[0]-1) * [1.,1.,1.] \
378 +whereZero(x[1]) * [1.,0.,1.] \
379 +whereZero(x[1]-1) * [1.,1.,1.] \
380 +whereZero(x[2]) * [1.,1.,0.] \
381 +whereZero(x[2]-1) * [1.,1.,1.]
382
383
384 sp=StokesProblemCartesian(self.domain)
385
386 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
387 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
388 p0=Scalar(-P1,ReducedSolution(self.domain))
389 sp.setTolerance(self.TOL)
390 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=False)
391
392 error_v0=Lsup(u[0]-u0[0])
393 error_v1=Lsup(u[1]-u0[1])
394 error_v2=Lsup(u[2]-u0[2])/0.25**2
395 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
396 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
397 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
398 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
399 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
400 #====================================================================================================================
401 class Test_Darcy(unittest.TestCase):
402 # this is a simple test for the darcy flux problem
403 #
404 #
405 # p = 1/k * ( 1/2* (fb-f0)/xb* x **2 + f0 * x - ub*x ) + p0
406 #
407 # with f = (fb-f0)/xb* x + f0
408 #
409 # u = f - k * p,x = ub
410 #
411 # we prescribe pressure at x=x0=0 to p0
412 #
413 # if we prescribe the pressure on the bottom x=xb we set
414 #
415 # pb= 1/k * ( 1/2* (fb-f0)* xb + f0 * xb - ub*xb ) + p0 = 1/k * ((fb+f0)/2 - ub ) * xb + p0
416 #
417 # which leads to ub = (fb+f0)/2-k*(pb-p0)/xb
418 #
419 def rescaleDomain(self):
420 x=self.dom.getX().copy()
421 for i in xrange(self.dom.getDim()):
422 x_inf=inf(x[i])
423 x_sup=sup(x[i])
424 if i == self.dom.getDim()-1:
425 x[i]=-self.WIDTH*(x[i]-x_sup)/(x_inf-x_sup)
426 else:
427 x[i]=self.WIDTH*(x[i]-x_inf)/(x_sup-x_inf)
428 self.dom.setX(x)
429 def getScalarMask(self,include_bottom=True):
430 x=self.dom.getX().copy()
431 x_inf=inf(x[self.dom.getDim()-1])
432 x_sup=sup(x[self.dom.getDim()-1])
433 out=whereZero(x[self.dom.getDim()-1]-x_sup)
434 if include_bottom: out+=whereZero(x[self.dom.getDim()-1]-x_inf)
435 return wherePositive(out)
436 def getVectorMask(self,include_bottom=True):
437 x=self.dom.getX().copy()
438 out=Vector(0.,Solution(self.dom))
439 for i in xrange(self.dom.getDim()):
440 x_inf=inf(x[i])
441 x_sup=sup(x[i])
442 if i != self.dom.getDim()-1: out[i]+=whereZero(x[i]-x_sup)
443 if i != self.dom.getDim()-1 or include_bottom: out[i]+=whereZero(x[i]-x_inf)
444 return wherePositive(out)
445
446 def setSolutionFixedBottom(self, p0, pb, f0, fb, k):
447 d=self.dom.getDim()
448 x=self.dom.getX()[d-1]
449 xb=inf(x)
450 u=Vector(0.,Solution(self.dom))+kronecker(d)[d-1]*((f0+fb)/2.-k*(pb-p0)/xb)
451 p=1./k*((fb-f0)/(xb*2.)* x**2 - (fb-f0)/2.*x)+(pb-p0)/xb*x + p0
452 f= ((fb-f0)/xb* x + f0)*kronecker(Function(self.dom))[d-1]
453 return u,p,f
454
455 def testConstF_FixedBottom_smallK(self):
456 k=1.e-10
457 mp=self.getScalarMask(include_bottom=True)
458 mv=self.getVectorMask(include_bottom=False)
459 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
460 p=p_ref*mp
461 u=u_ref*mv
462 df=DarcyFlow(self.dom)
463 df.setValue(g=f,
464 location_of_fixed_pressure=mp,
465 location_of_fixed_flux=mv,
466 permeability=Scalar(k,Function(self.dom)))
467 df.setTolerance(rtol=self.TOL)
468 v,p=df.solve(u_ref,p, max_iter=100, verbose=VERBOSE)
469 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
470 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
471 def testConstF_FixedBottom_mediumK(self):
472 k=1.
473 mp=self.getScalarMask(include_bottom=True)
474 mv=self.getVectorMask(include_bottom=False)
475 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
476 p=p_ref*mp
477 u=u_ref*mv
478 df=DarcyFlow(self.dom)
479 df.setValue(g=f,
480 location_of_fixed_pressure=mp,
481 location_of_fixed_flux=mv,
482 permeability=Scalar(k,Function(self.dom)))
483 df.setTolerance(rtol=self.TOL)
484 v,p=df.solve(u,p,max_iter=100, verbose=VERBOSE )
485 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
486 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
487
488 def testConstF_FixedBottom_largeK(self):
489 k=1.e10
490 mp=self.getScalarMask(include_bottom=True)
491 mv=self.getVectorMask(include_bottom=False)
492 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
493 p=p_ref*mp
494 u=u_ref*mv
495 df=DarcyFlow(self.dom)
496 df.setValue(g=f,
497 location_of_fixed_pressure=mp,
498 location_of_fixed_flux=mv,
499 permeability=Scalar(k,Function(self.dom)))
500 df.setTolerance(rtol=self.TOL)
501 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
502 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
503 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
504
505 def testVarioF_FixedBottom_smallK(self):
506 k=1.e-10
507 mp=self.getScalarMask(include_bottom=True)
508 mv=self.getVectorMask(include_bottom=False)
509 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
510 p=p_ref*mp
511 u=u_ref*mv
512 df=DarcyFlow(self.dom)
513 df.setValue(g=f,
514 location_of_fixed_pressure=mp,
515 location_of_fixed_flux=mv,
516 permeability=Scalar(k,Function(self.dom)))
517 df.setTolerance(rtol=self.TOL)
518 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
519 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
520 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
521
522 def testVarioF_FixedBottom_mediumK(self):
523 k=1.
524 mp=self.getScalarMask(include_bottom=True)
525 mv=self.getVectorMask(include_bottom=False)
526 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
527 p=p_ref*mp
528 u=u_ref*mv
529 df=DarcyFlow(self.dom)
530 df.setValue(g=f,
531 location_of_fixed_pressure=mp,
532 location_of_fixed_flux=mv,
533 permeability=Scalar(k,Function(self.dom)))
534 df.setTolerance(rtol=self.TOL)
535 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
536 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
537 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
538
539 def testVarioF_FixedBottom_largeK(self):
540 k=1.e10
541 mp=self.getScalarMask(include_bottom=True)
542 mv=self.getVectorMask(include_bottom=False)
543 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
544 p=p_ref*mp
545 u=u_ref*mv
546 df=DarcyFlow(self.dom)
547 df.setValue(g=f,
548 location_of_fixed_pressure=mp,
549 location_of_fixed_flux=mv,
550 permeability=Scalar(k,Function(self.dom)))
551 df.setTolerance(rtol=self.TOL)
552 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
553 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
554 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
555
556 def testConstF_FreeBottom_smallK(self):
557 k=1.e-10
558 mp=self.getScalarMask(include_bottom=False)
559 mv=self.getVectorMask(include_bottom=True)
560 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
561 p=p_ref*mp
562 u=u_ref*mv
563 df=DarcyFlow(self.dom)
564 df.setValue(g=f,
565 location_of_fixed_pressure=mp,
566 location_of_fixed_flux=mv,
567 permeability=Scalar(k,Function(self.dom)))
568 df.setTolerance(rtol=self.TOL)
569 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
570 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
571 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
572
573 def testConstF_FreeBottom_mediumK(self):
574 k=1.
575 mp=self.getScalarMask(include_bottom=False)
576 mv=self.getVectorMask(include_bottom=True)
577 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
578 p=p_ref*mp
579 u=u_ref*mv
580 df=DarcyFlow(self.dom)
581 df.setValue(g=f,
582 location_of_fixed_pressure=mp,
583 location_of_fixed_flux=mv,
584 permeability=Scalar(k,Function(self.dom)))
585 df.setTolerance(rtol=self.TOL)
586 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
587 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
588 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
589
590 def testConstF_FreeBottom_largeK(self):
591 k=1.e10
592 mp=self.getScalarMask(include_bottom=False)
593 mv=self.getVectorMask(include_bottom=True)
594 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
595 p=p_ref*mp
596 u=u_ref*mv
597 df=DarcyFlow(self.dom)
598 df.setValue(g=f,
599 location_of_fixed_pressure=mp,
600 location_of_fixed_flux=mv,
601 permeability=Scalar(k,Function(self.dom)))
602 df.setTolerance(rtol=self.TOL)
603 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
604 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
605 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
606
607 def testVarioF_FreeBottom_smallK(self):
608 k=1.e-10
609 mp=self.getScalarMask(include_bottom=False)
610 mv=self.getVectorMask(include_bottom=True)
611 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
612 p=p_ref*mp
613 u=u_ref*mv
614 df=DarcyFlow(self.dom)
615 df.setValue(g=f,
616 location_of_fixed_pressure=mp,
617 location_of_fixed_flux=mv,
618 permeability=Scalar(k,Function(self.dom)))
619 df.setTolerance(rtol=self.TOL)
620 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
621 self.failUnless(Lsup(v-u_ref)<self.TOL*25.*Lsup(u_ref), "flux error too big.") # 25 because of disc. error.
622 self.failUnless(Lsup(p-p_ref)<self.TOL*25.*Lsup(p_ref), "pressure error too big.")
623
624 def testVarioF_FreeBottom_mediumK(self):
625 k=1.
626 mp=self.getScalarMask(include_bottom=False)
627 mv=self.getVectorMask(include_bottom=True)
628 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
629 p=p_ref*mp
630 u=u_ref*mv
631 df=DarcyFlow(self.dom)
632 df.setValue(g=f,
633 location_of_fixed_pressure=mp,
634 location_of_fixed_flux=mv,
635 permeability=Scalar(k,Function(self.dom)))
636 df.setTolerance(rtol=self.TOL)
637 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
638 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
639 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
640
641 def testVarioF_FreeBottom_largeK(self):
642 k=1.e10
643 mp=self.getScalarMask(include_bottom=False)
644 mv=self.getVectorMask(include_bottom=True)
645 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
646 p=p_ref*mp
647 u=u_ref*mv
648 df=DarcyFlow(self.dom)
649 df.setValue(g=f,
650 location_of_fixed_pressure=mp,
651 location_of_fixed_flux=mv,
652 permeability=Scalar(k,Function(self.dom)))
653 df.setTolerance(rtol=self.TOL)
654 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
655 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
656 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
657
658 class Test_Darcy2D(Test_Darcy):
659 TOL=1e-4
660 WIDTH=1.
661 def setUp(self):
662 NE=40 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
663 self.dom = Rectangle(NE,NE)
664 self.rescaleDomain()
665 def tearDown(self):
666 del self.dom
667 class Test_Darcy3D(Test_Darcy):
668 TOL=1e-4
669 WIDTH=1.
670 def setUp(self):
671 NE=25 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
672 self.dom = Brick(NE,NE,NE)
673 self.rescaleDomain()
674 def tearDown(self):
675 del self.dom
676
677
678 class Test_Mountains3D(unittest.TestCase):
679 def setUp(self):
680 NE=16
681 self.TOL=1e-4
682 self.domain=Brick(NE,NE,NE,order=1,useFullElementOrder=True)
683 def tearDown(self):
684 del self.domain
685 def test_periodic(self):
686 EPS=0.01
687
688 x=self.domain.getX()
689 v = Vector(0.0, Solution(self.domain))
690 a0=1
691 a1=1
692 n0=2
693 n1=2
694 n2=0.5
695 a2=-(a0*n0+a1*n1)/n2
696 v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])* cos(pi*n2*x[2])
697 v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])* cos(pi*n2*x[2])
698 v[2]=a2*cos(pi*n0*x[0])* cos(pi*n1*x[1])* sin(pi*n2*x[2])
699
700 mts=Mountains(self.domain,eps=EPS)
701 mts.setVelocity(v)
702 Z=mts.update()
703
704 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[2]-1.)))
705 self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
706
707 class Test_Mountains2D(unittest.TestCase):
708 def setUp(self):
709 NE=16
710 self.TOL=1e-4
711 self.domain=Rectangle(NE,NE,order=1,useFullElementOrder=True)
712 def tearDown(self):
713 del self.domain
714 def test_periodic(self):
715 EPS=0.01
716
717 x=self.domain.getX()
718 v = Vector(0.0, Solution(self.domain))
719 a0=1
720 n0=1
721 n1=0.5
722 a1=-(a0*n0)/n1
723 v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])
724 v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])
725
726 H_t=Scalar(0.0, Solution(self.domain))
727 mts=Mountains(self.domain,eps=EPS)
728 mts.setVelocity(v)
729 Z=mts.update()
730
731 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[1]-1.)))
732 self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
733
734
735
736 class Test_Rheologies(unittest.TestCase):
737 """
738 this is the program used to generate the powerlaw tests:
739
740 TAU_Y=100.
741 N=10
742 M=5
743
744 def getE(tau):
745 if tau<=TAU_Y:
746 return 1./(0.5+20*sqrt(tau))
747 else:
748 raise ValueError,"out of range."
749 tau=[ i* (TAU_Y/N) for i in range(N+1)] + [TAU_Y for i in range(M)]
750 e=[ tau[i]/getE(tau[i]) for i in range(N+1)]
751 e+= [ (i+1+j)* (max(e)/N) for j in range(M) ]
752
753 print tau
754 print e
755 """
756 TOL=1.e-8
757 def test_PowerLaw_Init(self):
758 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
759
760 self.failUnlessEqual(pl.getNumMaterials(),3,"num materials is wrong")
761 self.failUnlessEqual(pl.validMaterialId(0),True,"material id 0 not found")
762 self.failUnlessEqual(pl.validMaterialId(1),True,"material id 1 not found")
763 self.failUnlessEqual(pl.validMaterialId(2),True,"material id 2 not found")
764 self.failUnlessEqual(pl.validMaterialId(3),False,"material id 3 not found")
765
766 self.failUnlessEqual(pl.getFriction(),None,"initial friction wrong.")
767 self.failUnlessEqual(pl.getTauY(),None,"initial tau y wrong.")
768 pl.setDruckerPragerLaw(tau_Y=10,friction=3)
769 self.failUnlessEqual(pl.getFriction(),3,"friction wrong.")
770 self.failUnlessEqual(pl.getTauY(),10,"tau y wrong.")
771
772 self.failUnlessEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
773 pl.setElasticShearModulus(1000)
774 self.failUnlessEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
775
776 e=pl.getEtaN()
777 self.failUnlessEqual(len(e),3,"initial length of etaN is wrong.")
778 self.failUnlessEqual(e,[None, None, None],"initial etaN are wrong.")
779 self.failUnlessEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
780 self.failUnlessEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
781 self.failUnlessEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
782 self.failUnlessRaises(ValueError, pl.getEtaN, 3)
783
784 self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
785 self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
786 self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
787
788 pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
789 self.failUnlessEqual(pl.getPower(),[1,2,3],"powers are wrong.")
790 self.failUnlessEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
791 self.failUnlessEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
792
793 pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
794 self.failUnlessEqual(pl.getPower(),[4,2,3],"powers are wrong.")
795 self.failUnlessEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
796 self.failUnlessEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
797
798 self.failUnlessRaises(ValueError,pl.getPower,-1)
799 self.failUnlessEqual(pl.getPower(0),4,"power 0 is wrong.")
800 self.failUnlessEqual(pl.getPower(1),2,"power 1 is wrong.")
801 self.failUnlessEqual(pl.getPower(2),3,"power 2 is wrong.")
802 self.failUnlessRaises(ValueError,pl.getPower,3)
803
804 self.failUnlessRaises(ValueError,pl.getTauT,-1)
805 self.failUnlessEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
806 self.failUnlessEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
807 self.failUnlessEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
808 self.failUnlessRaises(ValueError,pl.getTauT,3)
809
810 self.failUnlessRaises(ValueError,pl.getEtaN,-1)
811 self.failUnlessEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
812 self.failUnlessEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
813 self.failUnlessEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
814 self.failUnlessRaises(ValueError,pl.getEtaN,3)
815
816 def checkResult(self,id,gamma_dot_, eta, tau_ref):
817 self.failUnless(eta>=0,"eta needs to be positive (test %s)"%id)
818 error=abs(gamma_dot_*eta-tau_ref)
819 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))
820
821 def test_PowerLaw_Linear(self):
822 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]
823 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]
824 pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
825 pl.setDruckerPragerLaw(tau_Y=100.)
826 pl.setPowerLaw(eta_N=2.)
827 pl.setEtaTolerance(self.TOL)
828 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
829
830 def test_PowerLaw_QuadLarge(self):
831 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]
832 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]
833 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
834 pl.setDruckerPragerLaw(tau_Y=100.)
835 pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2])
836 pl.setEtaTolerance(self.TOL)
837 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
838
839 def test_PowerLaw_QuadSmall(self):
840 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]
841 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]
842 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
843 pl.setDruckerPragerLaw(tau_Y=100.)
844 pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2])
845 pl.setEtaTolerance(self.TOL)
846 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
847
848 def test_PowerLaw_CubeLarge(self):
849 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]
850 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]
851 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
852 pl.setDruckerPragerLaw(tau_Y=100.)
853 pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3])
854 pl.setEtaTolerance(self.TOL)
855 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
856
857 def test_PowerLaw_CubeSmall(self):
858 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]
859 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]
860 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
861 pl.setDruckerPragerLaw(tau_Y=100.)
862 pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3])
863 pl.setEtaTolerance(self.TOL)
864 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
865
866 def test_PowerLaw_QuadLarge_CubeLarge(self):
867 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]
868 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]
869
870 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
871 pl.setDruckerPragerLaw(tau_Y=100.)
872 pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
873 pl.setEtaTolerance(self.TOL)
874 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
875
876 def test_PowerLaw_QuadLarge_CubeSmall(self):
877 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]
878 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]
879
880 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
881 pl.setDruckerPragerLaw(tau_Y=100.)
882 pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
883 pl.setEtaTolerance(self.TOL)
884 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
885
886 def test_PowerLaw_QuadSmall_CubeLarge(self):
887 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]
888 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]
889 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
890 pl.setDruckerPragerLaw(tau_Y=100.)
891 pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
892 pl.setEtaTolerance(self.TOL)
893 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
894
895 def test_PowerLaw_QuadSmall_CubeSmall(self):
896 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]
897 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]
898 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
899 pl.setDruckerPragerLaw(tau_Y=100.)
900 pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
901 pl.setEtaTolerance(self.TOL)
902 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
903
904 def test_PowerLaw_withShear(self):
905 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]
906 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]
907 pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
908 pl.setDruckerPragerLaw(tau_Y=100.)
909 pl.setPowerLaw(eta_N=2.)
910 pl.setElasticShearModulus(3.)
911 dt=1./3.
912 pl.setEtaTolerance(self.TOL)
913 self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
914 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
915
916 class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase):
917 TOL=1.e-5
918 VERBOSE=False # or True
919 A=1.
920 P_max=100
921 NE=2*getMPISizeWorld()
922 tau_Y=10.
923 N_dt=10
924
925 # material parameter:
926 tau_1=5.
927 tau_2=5.
928 eta_0=100.
929 eta_1=50.
930 eta_2=400.
931 N_1=2.
932 N_2=3.
933 def getReference(self, t):
934
935 B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
936 x=self.dom.getX()
937
938 s_00=min(self.A*t,B)
939 tau=sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)*abs(s_00)
940 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.)
941
942 alpha=0.5*inv_eta*s_00
943 if s_00 <= B and self.mu !=None: alpha+=1./(2*self.mu)*self.A
944 u_ref=x*alpha
945 u_ref[self.dom.getDim()-1]=(1.-x[self.dom.getDim()-1])*alpha*(self.dom.getDim()-1)
946 sigma_ref=kronecker(self.dom)*s_00
947 sigma_ref[self.dom.getDim()-1,self.dom.getDim()-1]=-s_00*(self.dom.getDim()-1)
948
949 p_ref=self.P_max
950 for d in range(self.dom.getDim()): p_ref=p_ref*x[d]
951 p_ref-=integrate(p_ref)/vol(self.dom)
952 return u_ref, sigma_ref, p_ref
953
954 def runIt(self, free=None):
955 x=self.dom.getX()
956 B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
957 dt=B/int(self.N_dt/2)
958 if self.VERBOSE: print "dt =",dt
959 if self.latestart:
960 t=dt
961 else:
962 t=0
963 v,s,p=self.getReference(t)
964
965 mod=IncompressibleIsotropicFlowCartesian(self.dom, stress=s, v=v, p=p, t=t, numMaterials=3, verbose=self.VERBOSE)
966 mod.setDruckerPragerLaw(tau_Y=self.tau_Y,friction=None)
967 mod.setElasticShearModulus(self.mu)
968 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])
969 mod.setTolerance(self.TOL)
970 mod.setEtaTolerance(self.TOL*0.1)
971
972 BF=Vector(self.P_max,Function(self.dom))
973 for d in range(self.dom.getDim()):
974 for d2 in range(self.dom.getDim()):
975 if d!=d2: BF[d]*=x[d2]
976 v_mask=Vector(0,Solution(self.dom))
977 if free==None:
978 for d in range(self.dom.getDim()):
979 v_mask[d]=whereZero(x[d])+whereZero(x[d]-1.)
980 else:
981 for d in range(self.dom.getDim()):
982 if d == self.dom.getDim()-1:
983 v_mask[d]=whereZero(x[d]-1.)
984 else:
985 v_mask[d]=whereZero(x[d])
986 mod.setExternals(F=BF,fixed_v_mask=v_mask)
987
988 n=self.dom.getNormal()
989 N_t=0
990 errors=[]
991 while N_t < self.N_dt:
992 t_ref=t+dt
993 v_ref, s_ref,p_ref=self.getReference(t_ref)
994 mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref)
995 # mod.update(dt, eta_iter_max=10, iter_max=50, verbose=self.VERBOSE, usePCG=True, max_correction_steps=30) new version
996 mod.update(dt, iter_max=50, verbose=self.VERBOSE, usePCG=True)
997 self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)
998 t+=dt
999 N_t+=1
1000
1001 def check(self,N_t,mod,t_ref,v_ref, s_ref,p_ref):
1002 p=mod.getPressure()
1003 p-=integrate(p)/vol(self.dom)
1004 error_p=Lsup(mod.getPressure()-p_ref)/Lsup(p_ref)
1005 error_s=Lsup(mod.getDeviatoricStress()-s_ref)/Lsup(s_ref)
1006 error_v=Lsup(mod.getVelocity()-v_ref)/Lsup(v_ref)
1007 error_t=abs(mod.getTime()-t_ref)/abs(t_ref)
1008 if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v
1009 self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
1010 self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )
1011 self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )
1012 self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )
1013 def tearDown(self):
1014 del self.dom
1015
1016 def test_D2_Fixed_MuNone_LateStart(self):
1017 self.dom = Rectangle(self.NE,self.NE,order=2)
1018 self.mu=None
1019 self.latestart=True
1020 self.runIt()
1021 def test_D2_Fixed_Mu_LateStart(self):
1022 self.dom = Rectangle(self.NE,self.NE,order=2)
1023 self.mu=555.
1024 self.latestart=True
1025 self.runIt()
1026 def test_D2_Fixed_MuNone(self):
1027 self.dom = Rectangle(self.NE,self.NE,order=2)
1028 self.mu=None
1029 self.latestart=False
1030 self.runIt()
1031 def test_D2_Fixed_Mu(self):
1032 self.dom = Rectangle(self.NE,self.NE,order=2)
1033 self.mu=555.
1034 self.latestart=False
1035 self.runIt()
1036 def test_D2_Free_MuNone_LateStart(self):
1037 self.dom = Rectangle(self.NE,self.NE,order=2)
1038 self.mu=None
1039 self.latestart=True
1040 self.runIt(free=0)
1041 def test_D2_Free_Mu_LateStart(self):
1042 self.dom = Rectangle(self.NE,self.NE,order=2)
1043 self.mu=555.
1044 self.latestart=True
1045 self.runIt(free=0)
1046 def test_D2_Free_MuNone(self):
1047 self.dom = Rectangle(self.NE,self.NE,order=2)
1048 self.mu=None
1049 self.latestart=False
1050 self.runIt(free=0)
1051 def test_D2_Free_Mu(self):
1052 self.dom = Rectangle(self.NE,self.NE,order=2)
1053 self.mu=555.
1054 self.latestart=False
1055 self.runIt(free=0)
1056
1057 def test_D3_Fixed_MuNone_LateStart(self):
1058 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1059 self.mu=None
1060 self.latestart=True
1061 self.runIt()
1062 def test_D3_Fixed_Mu_LateStart(self):
1063 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1064 self.mu=555.
1065 self.latestart=True
1066 self.runIt()
1067 def test_D3_Fixed_MuNone(self):
1068 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1069 self.mu=None
1070 self.latestart=False
1071 self.runIt()
1072 def test_D3_Fixed_Mu(self):
1073 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1074 self.mu=555.
1075 self.latestart=False
1076 self.runIt()
1077 def test_D3_Free_MuNone_LateStart(self):
1078 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1079 self.mu=None
1080 self.latestart=True
1081 self.runIt(free=0)
1082 def test_D3_Free_Mu_LateStart(self):
1083 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1084 self.mu=555.
1085 self.latestart=True
1086 self.runIt(free=0)
1087 def test_D3_Free_MuNone(self):
1088 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1089 self.mu=None
1090 self.latestart=False
1091 self.runIt(free=0)
1092 def test_D3_Free_Mu(self):
1093 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1094 self.mu=555.
1095 self.latestart=False
1096 self.runIt(free=0)
1097
1098
1099 class Test_FaultSystem(unittest.TestCase):
1100 EPS=1.e-8
1101 NE=10
1102 def test_Fault_MaxValue(self):
1103 dom=Rectangle(2*self.NE,2*self.NE)
1104 x=dom.getX()
1105 f=FaultSystem(dim=2)
1106 f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
1107 f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
1108
1109 u=x[0]*(1.-x[0])*(1-x[1])
1110 t, loc=f.getMaxValue(u)
1111 p=f.getParametrization(x,t)[0]
1112 m, l=loc(u), loc(p)
1113 self.failUnless( m == 0.25, "wrong max value")
1114 self.failUnless( t == 1, "wrong max tag")
1115 self.failUnless( l == 0., "wrong max location")
1116
1117 u=x[1]*(1.-x[1])*(1-x[0])*x[0]
1118 t, loc=f.getMaxValue(u)
1119 p=f.getParametrization(x,t)[0]
1120 m, l=loc(u), loc(p)
1121 self.failUnless( m == 0.0625, "wrong max value")
1122 self.failUnless( t == 2, "wrong max tag")
1123 self.failUnless( l == 0.5, "wrong max location")
1124
1125 u=x[0]*(1.-x[0])*x[1]
1126 t, loc=f.getMaxValue(u)
1127 p=f.getParametrization(x,t)[0]
1128 m, l=loc(u), loc(p)
1129 self.failUnless( m == 0.25, "wrong max value")
1130 self.failUnless( t == 2, "wrong max tag")
1131 self.failUnless( l == 1.0, "wrong max location")
1132
1133 u=x[1]*(1.-x[1])*x[0]
1134 t, loc=f.getMaxValue(u)
1135 p=f.getParametrization(x,t)[0]
1136 m, l=loc(u), loc(p)
1137 self.failUnless( m == 0.25, "wrong max value")
1138 self.failUnless( t == 2, "wrong max tag")
1139 self.failUnless( l == 0., "wrong max location")
1140
1141 u=x[1]*(1.-x[1])*(1.-x[0])
1142 t, loc=f.getMaxValue(u)
1143 p=f.getParametrization(x,t)[0]
1144 m, l=loc(u), loc(p)
1145 self.failUnless( m == 0.25, "wrong max value")
1146 self.failUnless( t == 1, "wrong max tag")
1147 self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong max location")
1148 def test_Fault_MinValue(self):
1149 dom=Rectangle(2*self.NE,2*self.NE)
1150 x=dom.getX()
1151 f=FaultSystem(dim=2)
1152 f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
1153 f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
1154
1155 u=-x[0]*(1.-x[0])*(1-x[1])
1156 t, loc=f.getMinValue(u)
1157 p=f.getParametrization(x,t)[0]
1158 m, l=loc(u), loc(p)
1159 self.failUnless( m == -0.25, "wrong min value")
1160 self.failUnless( t == 1, "wrong min tag")
1161 self.failUnless( l == 0., "wrong min location")
1162 u=-x[1]*(1.-x[1])*(1-x[0])*x[0]
1163 t, loc=f.getMinValue(u)
1164 p=f.getParametrization(x,t)[0]
1165 m, l=loc(u), loc(p)
1166 self.failUnless( m == -0.0625, "wrong min value")
1167 self.failUnless( t == 2, "wrong min tag")
1168 self.failUnless( l == 0.5, "wrong min location")
1169 u=-x[0]*(1.-x[0])*x[1]
1170 t, loc=f.getMinValue(u)
1171 p=f.getParametrization(x,t)[0]
1172 m, l=loc(u), loc(p)
1173 self.failUnless( m == -0.25, "wrong min value")
1174 self.failUnless( t == 2, "wrong min tag")
1175 self.failUnless( l == 1.0, "wrong min location")
1176 u=-x[1]*(1.-x[1])*x[0]
1177 t, loc=f.getMinValue(u)
1178 p=f.getParametrization(x,t)[0]
1179 m, l=loc(u), loc(p)
1180 self.failUnless( m == -0.25, "wrong min value")
1181 self.failUnless( t == 2, "wrong min tag")
1182 self.failUnless( l == 0., "wrong min location")
1183 u=-x[1]*(1.-x[1])*(1.-x[0])
1184 t, loc=f.getMinValue(u)
1185 p=f.getParametrization(x,t)[0]
1186 m, l=loc(u), loc(p)
1187 self.failUnless( m == -0.25, "wrong min value")
1188 self.failUnless( t == 1, "wrong min tag")
1189 self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong min location")
1190
1191
1192 def test_Fault2D(self):
1193 f=FaultSystem(dim=2)
1194 top1=[ [1.,0.], [1.,1.], [0.,1.] ]
1195 self.failUnlessRaises(ValueError,f.addFault,V0=[1.,0],strikes=[pi/2, pi/2],ls=[1.,1.],tag=1,dips=top1)
1196 f.addFault(V0=[1.,0],strikes=[pi/2, pi],ls=[1.,1.],tag=1)
1197 self.failUnless(f.getDim() == 2, "wrong dimension")
1198 self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
1199 self.failUnless( 2. == f.getTotalLength(1), "length wrong")
1200 self.failUnless( 0. == f.getMediumDepth(1), "depth wrong")
1201 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 range")
1202 self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 range")
1203 self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsets")
1204 segs=f.getTopPolyline(1)
1205 self.failUnless( len(segs) == 3, "wrong number of segments")
1206 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1207 self.failUnless( segs[0].size == 2, "seg 0 has wrong size.")
1208 self.failUnless( numpy.linalg.norm(segs[0]-[1.,0.]) < self.EPS, "wrong vertex. 0 ")
1209 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1210 self.failUnless( segs[1].size == 2, "seg 1 has wrong size.")
1211 self.failUnless( numpy.linalg.norm(segs[1]-[1.,1.]) < self.EPS, "wrong vertex. 1 ")
1212 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1213 self.failUnless( segs[2].size == 2, "seg 2 has wrong size.")
1214 self.failUnless( numpy.linalg.norm(segs[2]-[0.,1.]) < self.EPS, "wrong vertex. 2 ")
1215 c=f.getCenterOnSurface()
1216 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1217 self.failUnless( c.size == 2, "center size is wrong")
1218 self.failUnless( numpy.linalg.norm(c-[2./3.,2./3.]) < self.EPS, "center has wrong coordinates.")
1219 o=f.getOrientationOnSurface()/pi*180.
1220 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1221
1222 top2=[ [10.,0.], [0.,10.] ]
1223 f.addFault(V0=[10.,0],strikes=[3.*pi/4],ls=[14.142135623730951], tag=2, w0_offsets=[0,20], w1_max=20)
1224 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1225 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length")
1226 self.failUnless( 0. == f.getMediumDepth(2), "depth wrong")
1227 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range")
1228 self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range")
1229 self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets")
1230 segs=f.getTopPolyline(2)
1231 self.failUnless( len(segs) == 2, "wrong number of segments")
1232 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1233 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.]) < self.EPS, "wrong vertex. 0 ")
1234 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1235 self.failUnless( numpy.linalg.norm(segs[1]-[0.,10.]) < self.EPS, "wrong vertex. 1 ")
1236 c=f.getCenterOnSurface()
1237 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1238 self.failUnless( c.size == 2, "center size is wrong")
1239 self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.]) < self.EPS, "center has wrong coordinates.")
1240 o=f.getOrientationOnSurface()/pi*180.
1241 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1242
1243 s,d=f.getSideAndDistance([0.,0.], tag=1)
1244 self.failUnless( s<0, "wrong side.")
1245 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1246 s,d=f.getSideAndDistance([0.,2.], tag=1)
1247 self.failUnless( s>0, "wrong side.")
1248 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1249 s,d=f.getSideAndDistance([1.,2.], tag=1)
1250 self.failUnless( s>0, "wrong side.")
1251 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1252 s,d=f.getSideAndDistance([2.,1.], tag=1)
1253 self.failUnless( s>0, "wrong side.")
1254 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1255 s,d=f.getSideAndDistance([2.,0.], tag=1)
1256 self.failUnless( s>0, "wrong side.")
1257 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1258 s,d=f.getSideAndDistance([0.,-1.], tag=1)
1259 self.failUnless( s<0, "wrong side.")
1260 self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
1261 s,d=f.getSideAndDistance([-1.,0], tag=1)
1262 self.failUnless( s<0, "wrong side.")
1263 self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
1264
1265
1266 f.transform(rot=-pi/2., shift=[-1.,-1.])
1267 self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
1268 self.failUnless( 2. == f.getTotalLength(1), "length after transformation wrong")
1269 self.failUnless( 0. == f.getMediumDepth(1), "depth after transformation wrong")
1270 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 after transformation range")
1271 self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
1272 self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsetsafter transformation ")
1273 segs=f.getTopPolyline(1)
1274 self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
1275 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1276 self.failUnless( segs[0].size == 2, "seg 0 has wrong size after transformation.")
1277 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1278 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1279 self.failUnless( segs[1].size == 2, "seg 1 has wrong size after transformation.")
1280 self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1281 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1282 self.failUnless( segs[2].size == 2, "seg 2 has wrong size after transformation.")
1283 self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.]) < self.EPS, "wrong vertex after transformation. 2 ")
1284 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1285 self.failUnless( 0. == f.getMediumDepth(2), "depth wrong after transformation")
1286 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range after transformation")
1287 self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1288 self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1289 segs=f.getTopPolyline(2)
1290 self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
1291 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1292 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9]) < self.EPS, "wrong vertex. 0 after transformation")
1293 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1294 self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.]) < self.EPS, "wrong vertex. 1 after transformation")
1295
1296 c=f.getCenterOnSurface()
1297 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1298 self.failUnless( c.size == 2, "center size is wrong")
1299 self.failUnless( numpy.linalg.norm(c-[7./5.,-7./5.]) < self.EPS, "center has wrong coordinates.")
1300 o=f.getOrientationOnSurface()/pi*180.
1301 self.failUnless( abs(o-45.) < self.EPS, "wrong orientation.")
1302
1303 p=f.getParametrization([-1.,0.],1)
1304 self.failUnless(p[1]==1., "wrong value.")
1305 self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1306 p=f.getParametrization([-0.5,0.],1)
1307 self.failUnless(p[1]==1., "wrong value.")
1308 self.failUnless(abs(p[0]-0.5)<self.EPS* 0.5, "wrong value.")
1309 p=f.getParametrization([0.,0.],1)
1310 self.failUnless(p[1]==1., "wrong value.")
1311 self.failUnless(abs(p[0]-1.)<self.EPS, "wrong value.")
1312 p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-8)
1313 self.failUnless(p[1]==0., "wrong value.")
1314 p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-6)
1315 self.failUnless(p[1]==1., "wrong value.")
1316 self.failUnless(abs(p[0]-1.0000001)<self.EPS, "wrong value.")
1317 p=f.getParametrization([0.,0.5],1)
1318 self.failUnless(p[1]==1., "wrong value.")
1319 self.failUnless(abs(p[0]-1.5)<self.EPS, "wrong value.")
1320 p=f.getParametrization([0,1.],1)
1321 self.failUnless(p[1]==1., "wrong value.")
1322 self.failUnless(abs(p[0]-2.)<self.EPS, "wrong value.")
1323 p=f.getParametrization([1.,1.],1)
1324 self.failUnless(p[1]==0., "wrong value.")
1325 p=f.getParametrization([0,1.11],1)
1326 self.failUnless(p[1]==0., "wrong value.")
1327 p=f.getParametrization([-1,-9.],2)
1328 self.failUnless(p[1]==1., "wrong value.")
1329 self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1330 p=f.getParametrization([9,1],2)
1331 self.failUnless(p[1]==1., "wrong value.")
1332 self.failUnless(abs(p[0]-20.)<self.EPS, "wrong value.")
1333
1334 def test_Fault3D(self):
1335 f=FaultSystem(dim=3)
1336 self.failUnless(f.getDim() == 3, "wrong dimension")
1337
1338 top1=[ [0.,0.,0.], [1., 0., 0.] ]
1339 f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=20.,tag=1)
1340 self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
1341 self.failUnless( 1. == f.getTotalLength(1), "length wrong")
1342 self.failUnless( 20. == f.getMediumDepth(1), "depth wrong")
1343 self.failUnless( (0., 1.) == f.getW0Range(1)," wrong W0 range")
1344 self.failUnless( (-20., 0.) == f.getW1Range(1)," wrong W1 range")
1345 self.failUnless( [0., 1.] == f.getW0Offsets(1)," wrong W0 offsets")
1346 segs=f.getTopPolyline(1)
1347 self.failUnless( len(segs) == 2, "wrong number of segments")
1348 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1349 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1350 self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1351 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1352 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1353 self.failUnless( numpy.linalg.norm(segs[1]-[1.,0.,0]) < self.EPS, "wrong vertex. 1 ")
1354 c=f.getCenterOnSurface()
1355 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1356 self.failUnless( c.size == 3, "center size is wrong")
1357 self.failUnless( numpy.linalg.norm(c-[0.5,0.,0.]) < self.EPS, "center has wrong coordinates.")
1358 o=f.getOrientationOnSurface()/pi*180.
1359 self.failUnless( abs(o) < self.EPS, "wrong orientation.")
1360 d=f.getDips(1)
1361 self.failUnless( len(d) == 1, "wrong number of dips")
1362 self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1363 sn=f.getSegmentNormals(1)
1364 self.failUnless( len(sn) == 1, "wrong number of normals")
1365 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1366 self.failUnless( numpy.linalg.norm(sn[0]-[0, -1., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1367 dv=f.getDepthVectors(1)
1368 self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1369 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1370 self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1371 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1372 self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1373 b=f.getBottomPolyline(1)
1374 self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1375 self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1376 self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1377 self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1378 self.failUnless( numpy.linalg.norm(b[1]-[1., 0., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1379 ds=f.getDepths(1)
1380 self.failUnless( len(ds) == 2, "wrong number of depth")
1381 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1382 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1383
1384 top2=[ [0.,0.,0.], [0., 10., 0.] ]
1385 f.addFault(V0=[0.,0,0],strikes=[pi/2],ls=[10.,], dips=pi/2, depths=20.,tag=2)
1386 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1387 self.failUnless( 10. == f.getTotalLength(2), "length wrong")
1388 self.failUnless( 20. == f.getMediumDepth(2), "depth wrong")
1389 self.failUnless( (0., 10.) == f.getW0Range(2)," wrong W0 range")
1390 self.failUnless( (-20., 0.) == f.getW1Range(2)," wrong W1 range")
1391 self.failUnless( [0., 10.] == f.getW0Offsets(2)," wrong W0 offsets")
1392 segs=f.getTopPolyline(2)
1393 self.failUnless( len(segs) == 2, "wrong number of segments")
1394 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1395 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1396 self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1397 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1398 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1399 self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1400 d=f.getDips(2)
1401 self.failUnless( len(d) == 1, "wrong number of dips")
1402 self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1403 sn=f.getSegmentNormals(2)
1404 self.failUnless( len(sn) == 1, "wrong number of normals")
1405 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1406 self.failUnless( numpy.linalg.norm(sn[0]-[1, 0., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1407 dv=f.getDepthVectors(2)
1408 self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1409 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1410 self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1411 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1412 self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1413 b=f.getBottomPolyline(2)
1414 self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1415 self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1416 self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1417 self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1418 self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1419 ds=f.getDepths(2)
1420 self.failUnless( len(ds) == 2, "wrong number of depth")
1421 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1422 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1423
1424 top2=[ [10.,0.,0.], [0., 10., 0.] ]
1425 f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/2, depths=30.,tag=2)
1426 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1427 self.failUnless( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1428 self.failUnless( 30. == f.getMediumDepth(2), "depth wrong")
1429 self.failUnless( (-30., 0.) == f.getW1Range(2)," wrong W1 range")
1430 segs=f.getTopPolyline(2)
1431 self.failUnless( len(segs) == 2, "wrong number of segments")
1432 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1433 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1434 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1435 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1436 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1437 self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1438 d=f.getDips(2)
1439 self.failUnless( len(d) == 1, "wrong number of dips")
1440 self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1441 sn=f.getSegmentNormals(2)
1442 self.failUnless( len(sn) == 1, "wrong number of normals")
1443 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1444 self.failUnless( numpy.linalg.norm(sn[0]-[0.70710678118654746, 0.70710678118654746, 0.]) < self.EPS, "wrong bottom vertex 1 ")
1445 dv=f.getDepthVectors(2)
1446 self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1447 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1448 self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -30.]) < self.EPS, "wrong depth vector 0 ")
1449 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1450 self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -30.]) < self.EPS, "wrong depth vector 1 ")
1451 b=f.getBottomPolyline(2)
1452 self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1453 self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1454 self.failUnless( numpy.linalg.norm(b[0]-[10., 0., -30.]) < self.EPS, "wrong bottom vertex 0 ")
1455 self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1456 self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -30.]) < self.EPS, "wrong bottom vertex 1 ")
1457 ds=f.getDepths(2)
1458 self.failUnless( len(ds) == 2, "wrong number of depth")
1459 self.failUnless( abs(ds[0]-30.) < self.EPS, "wrong depth at vertex 0 ")
1460 self.failUnless( abs(ds[1]-30.) < self.EPS, "wrong depth at vertex 1 ")
1461
1462 top2=[ [10.,0.,0.], [0., 10., 0.] ]
1463 f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/4, depths=50.,tag=2)
1464 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1465 self.failUnless( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1466 self.failUnless( 50. == f.getMediumDepth(2), "depth wrong")
1467 self.failUnless( (-50., 0.) == f.getW1Range(2)," wrong W1 range")
1468 segs=f.getTopPolyline(2)
1469 self.failUnless( len(segs) == 2, "wrong number of segments")
1470 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1471 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1472 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1473 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1474 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1475 self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1476 d=f.getDips(2)
1477 self.failUnless( len(d) == 1, "wrong number of dips")
1478 self.failUnless( abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1479 sn=f.getSegmentNormals(2)
1480 self.failUnless( len(sn) == 1, "wrong number of normals")
1481 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1482 self.failUnless( numpy.linalg.norm(sn[0]-[0.5,0.5,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1483 dv=f.getDepthVectors(2)
1484 self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1485 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1486 self.failUnless( numpy.linalg.norm(dv[0]-[25., 25., -35.355339059327378]) < self.EPS, "wrong depth vector 0 ")
1487 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1488 self.failUnless( numpy.linalg.norm(dv[1]-[25.,25., -35.355339059327378]) < self.EPS, "wrong depth vector 1 ")
1489 b=f.getBottomPolyline(2)
1490 self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1491 self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1492 self.failUnless( numpy.linalg.norm(b[0]-[35., 25., -35.355339059327378]) < self.EPS, "wrong bottom vertex 0 ")
1493 self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1494 self.failUnless( numpy.linalg.norm(b[1]-[25, 35., -35.355339059327378]) < self.EPS, "wrong bottom vertex 1 ")
1495 ds=f.getDepths(2)
1496 self.failUnless( len(ds) == 2, "wrong number of depth")
1497 self.failUnless( abs(ds[0]-50.) < self.EPS, "wrong depth at vertex 0 ")
1498 self.failUnless( abs(ds[1]-50.) < self.EPS, "wrong depth at vertex 1 ")
1499
1500 top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
1501 f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=1, dips=pi/4, depths=20.)
1502 self.failUnless( 20. == f.getTotalLength(1), "length wrong")
1503 self.failUnless( 20. == f.getMediumDepth(1), "depth wrong")
1504 segs=f.getTopPolyline(1)
1505 self.failUnless( len(segs) == 3, "wrong number of segments")
1506 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1507 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1508 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1509 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1510 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1511 self.failUnless( numpy.linalg.norm(segs[1]-[10.,10.,0.]) < self.EPS, "wrong vertex. 1 ")
1512 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1513 self.failUnless( segs[2].size == 3, "seg 2 has wrong size.")
1514 self.failUnless( numpy.linalg.norm(segs[2]-[0.,10.,0.]) < self.EPS, "wrong vertex. 2 ")
1515 d=f.getDips(1)
1516 self.failUnless( len(d) == 2, "wrong number of dips")
1517 self.failUnless( abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1518 self.failUnless( abs(d[1]-0.78539816339744828) < self.EPS, "wrong dip 0")
1519 ds=f.getDepths(1)
1520 self.failUnless( len(ds) == 3, "wrong number of depth")
1521 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1522 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1523 sn=f.getSegmentNormals(1)
1524 self.failUnless( len(sn) == 2, "wrong number of normals")
1525 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1526 self.failUnless( numpy.linalg.norm(sn[0]-[0.70710678118654746,0.,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1527 self.failUnless( isinstance(sn[1], numpy.ndarray), "wrong class of bottom vertex 0")
1528 self.failUnless( numpy.linalg.norm(sn[1]-[0.,0.70710678118654746,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1529 dv=f.getDepthVectors(1)
1530 self.failUnless( len(dv) == 3, "wrong number of depth vectors.")
1531 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1532 self.failUnless( numpy.linalg.norm(dv[0]-[14.142135623730951, 0., -14.142135623730951]) < self.EPS, "wrong depth vector 0 ")
1533 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1534 self.failUnless( numpy.linalg.norm(dv[1]-[11.547005383792515,11.547005383792515, -11.547005383792515]) < self.EPS, "wrong depth vector 2 ")
1535 self.failUnless( isinstance(dv[2], numpy.ndarray), "wrong class of depth vector 1")
1536 self.failUnless( numpy.linalg.norm(dv[2]-[0.,14.142135623730951, -14.142135623730951]) < self.EPS, "wrong depth vector 2 ")
1537 segs=f.getBottomPolyline(1)
1538 self.failUnless( len(segs) == 3, "wrong number of segments")
1539 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1540 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1541 self.failUnless( numpy.linalg.norm(segs[0]-[24.142135623730951,0.,-14.142135623730951]) < self.EPS, "wrong vertex. 0 ")
1542 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1543 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1544 self.failUnless( numpy.linalg.norm(segs[1]-[21.547005383792515,21.547005383792515, -11.547005383792515]) < self.EPS, "wrong vertex. 1 ")
1545 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1546 self.failUnless( segs[2].size == 3, "seg 2 has wrong size.")
1547 self.failUnless( numpy.linalg.norm(segs[2]-[0., 24.142135623730951, -14.142135623730951]) < self.EPS, "wrong vertex. 2 ")
1548 self.failUnless( abs(0.-f.getW0Range(1)[0]) <=self.EPS," wrong W0 range (0)")
1549 self.failUnless( abs(31.857329272664341-f.getW0Range(1)[1]) <=self.EPS," wrong W0 range (1)")
1550 self.failUnless( abs(-20. - f.getW1Range(1)[0]) <=self.EPS," wrong W1 range (0)")
1551 self.failUnless( abs(0. - f.getW1Range(1)[1]) <=self.EPS," wrong W1 range (1)")
1552 self.failUnless( abs(0.0-f.getW0Offsets(1)[0])<=self.EPS," wrong W0 offsets (0)")
1553 self.failUnless( abs(15.92866463633217-f.getW0Offsets(1)[1])<=self.EPS," wrong W0 offsets (1)")
1554 self.failUnless( abs(31.857329272664341-f.getW0Offsets(1)[2])<=self.EPS," wrong W0 offsets(2)")
1555 #
1556 # ============ fresh start ====================
1557 #
1558 f.addFault(V0=[1.,0,0.],strikes=[pi/2, pi],ls=[1.,1.], dips=pi/2,depths=20,tag=1)
1559 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)
1560 c=f.getCenterOnSurface()
1561 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1562 self.failUnless( c.size == 3, "center size is wrong")
1563 self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.,0.]) < self.EPS, "center has wrong coordinates.")
1564 o=f.getOrientationOnSurface()/pi*180.
1565 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1566
1567 f.transform(rot=-pi/2., shift=[-1.,-1.,0.])
1568 self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
1569 self.failUnless( 2. == f.getTotalLength(1), "length after transformation wrong")
1570 self.failUnless( 20. == f.getMediumDepth(1), "depth after transformation wrong")
1571 rw0=f.getW0Range(1)
1572 self.failUnless( len(rw0) ==2, "wo range has wrong length")
1573 self.failUnless( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
1574 self.failUnless( abs(rw0[1]-2.) < self.EPS,"W0 1 wrong.")
1575 self.failUnless( (-20., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
1576 dips=f.getDips(1)
1577 self.failUnless(len(dips) == 2, "wrong number of dips.")
1578 self.failUnless( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
1579 self.failUnless( abs(dips[1]-1.5707963267948966) <= self.EPS, "wrong dip")
1580 ds=f.getDepths(1)
1581 self.failUnless( len(ds) == 3, "wrong number of depth")
1582 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1583 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1584 self.failUnless( abs(ds[2]-20.) < self.EPS, "wrong depth at vertex 1 ")
1585 segs=f.getTopPolyline(1)
1586 self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
1587 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1588 self.failUnless( segs[0].size == 3, "seg 0 has wrong size after transformation.")
1589 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1590 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1591 self.failUnless( segs[1].size == 3, "seg 1 has wrong size after transformation.")
1592 self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1593 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1594 self.failUnless( segs[2].size == 3, "seg 2 has wrong size after transformation.")
1595 self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.,0.]) < self.EPS, "wrong vertex after transformation. 2 ")
1596 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1597 self.failUnless( 20. == f.getMediumDepth(2), "depth wrong after transformation")
1598 rw0=f.getW0Range(2)
1599 self.failUnless( len(rw0) ==2, "wo range has wrong length")
1600 self.failUnless( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
1601 self.failUnless( abs(rw0[1]-20.) < self.EPS,"W0 1 wrong.")
1602 self.failUnless( (-20., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1603 self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1604 dips=f.getDips(2)
1605 self.failUnless(len(dips) == 1, "wrong number of dips.")
1606 self.failUnless( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
1607 ds=f.getDepths(2)
1608 self.failUnless( len(ds) == 2, "wrong number of depth")
1609 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1610 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1611 segs=f.getTopPolyline(2)
1612 self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
1613 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1614 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1615 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1616 self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.,0.]) < self.EPS, "wrong vertex. 1 after transformation")
1617 #
1618 # ============ fresh start ====================
1619 #
1620 f=FaultSystem(dim=3)
1621
1622 top1=[ [0.,0.,0.], [1., 0., 0.] ]
1623 f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=1.,tag=1)
1624 top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
1625 f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=2, dips=pi/4, depths=20.)
1626
1627 p,m=f.getParametrization([0.3,0.,-0.5],1)
1628 self.failUnless(length(p-[0.3,-0.5]) <= self.EPS, "wrong value.")
1629 self.failUnless(m==1., "wrong value.")
1630
1631 p,m=f.getParametrization([0.5,0.,-0.5],1)
1632 self.failUnless(length(p-[0.5,-0.5]) <= self.EPS, "wrong value.")
1633 self.failUnless(m==1., "wrong value.")
1634
1635 p,m=f.getParametrization([0.25,0.,-0.5],1)
1636 self.failUnless(length(p-[0.25,-0.5]) <= self.EPS, "wrong value.")
1637 self.failUnless(m==1., "wrong value.")
1638
1639 p,m=f.getParametrization([0.5,0.,-0.25],1)
1640 self.failUnless(length(p-[0.5,-0.25]) <= self.EPS, "wrong value.")
1641 self.failUnless(m==1., "wrong value.")
1642
1643 p,m=f.getParametrization([0.001,0.,-0.001],1)
1644 self.failUnless(length(p-[0.001, -0.001]) <= self.EPS, "wrong value.")
1645 self.failUnless(m==1., "wrong value.")
1646
1647 p,m=f.getParametrization([0.001,0.,0.001],1)
1648 self.failUnless(m==0., "wrong value.")
1649
1650 p,m=f.getParametrization([0.999,0.,0.001],1)
1651 self.failUnless(m==0., "wrong value.")
1652
1653 p,m=f.getParametrization([1.001,0.,-0.001],1)
1654 self.failUnless(m==0., "wrong value.")
1655 p,m=f.getParametrization([1.001,0.,-0.1],1)
1656 self.failUnless(m==0., "wrong value.")
1657 p,m=f.getParametrization([1.001,0.,-0.000000001],1)
1658 self.failUnless(m==0., "wrong value.")
1659
1660 p,m=f.getParametrization([0.999,0.,-0.001],1)
1661 self.failUnless(length(p-[0.999, -0.001]) <= self.EPS, "wrong value.")
1662 self.failUnless(m==1., "wrong value.")
1663
1664 p,m=f.getParametrization([ 16.29252873 , 6.46410161 ,-6.29252873], 2, tol=1.e-7)
1665 self.failUnless(m==1., "wrong value.")
1666 self.failUnless(length(p-[4.7785993908996511, -10]) <= self.EPS*10., "wrong value.")
1667 p,m=f.getParametrization([15.77350269, 12.77350269, -5.77350269], 2, tol=1.e-7)
1668 self.failUnless(m==1., "wrong value.")
1669 self.failUnless(length(p-[11.150065245432518, -10]) <= self.EPS*10., "wrong value.")
1670
1671 p,m=f.getParametrization([ 3., 17.0710678, -7.0710678], 2, tol=1.e-7)
1672 self.failUnless(m==1., "wrong value.")
1673 self.failUnless(length(p-[27.078729881764687, -10]) <= self.EPS*10., "wrong value.")
1674 p,m=f.getParametrization([9.30940108, 16.55204176, -6.55204176], 2, tol=1.e-7)
1675 self.failUnless(m==1., "wrong value.")
1676 self.failUnless(length(p-[20.707264027231822, -10]) <= self.EPS*10., "wrong value.")
1677
1678 p,m=f.getParametrization([ 21.54700538, 21.54700538, -11.54700538], 2, tol=1.e-7)
1679 self.failUnless(m==1., "wrong value.")
1680 self.failUnless(length(p-[15.92866463633217, -20]) <= self.EPS*10., "wrong value.")
1681
1682 p,m=f.getParametrization([ 0.,0.,0.], 2, tol=1.e-7)
1683 self.failUnless(m==0., "wrong value.")
1684
1685 p,m=f.getParametrization([ 11.,11.,0.], 2, tol=1.e-7)
1686 self.failUnless(m==0., "wrong value.")
1687
1688
1689 s,d=f.getSideAndDistance([0.,-1.,0.], tag=1)
1690 self.failUnless( s>0, "wrong side.")
1691 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1692 s,d=f.getSideAndDistance([1.,-1.,0.], tag=1)
1693 self.failUnless( s>0, "wrong side.")
1694 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1695 s,d=f.getSideAndDistance([0.,1.,0.], tag=1)
1696 self.failUnless( s<0, "wrong side.")
1697 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1698 s,d=f.getSideAndDistance([1.,1.,0.], tag=1)
1699 self.failUnless( s<0, "wrong side.")
1700 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1701
1702
1703 s,d=f.getSideAndDistance([0.,0.,0.], tag=2)
1704 self.failUnless( s<0, "wrong side.")
1705 self.failUnless( abs(d-10.)<self.EPS, "wrong distance.")
1706 s,d=f.getSideAndDistance([5.,5.,0.], tag=2)
1707 self.failUnless( s<0, "wrong side.")
1708 self.failUnless( abs(d-5.)<self.EPS, "wrong distance.")
1709
1710 s,d=f.getSideAndDistance([10.,10.,-1.], tag=2)
1711 self.failUnless( s<0, "wrong side.")
1712 self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
1713 s,d=f.getSideAndDistance([10.,10.,-2.], tag=2)
1714 self.failUnless( s<0, "wrong side.")
1715 self.failUnless( abs(d-2.*0.70710678118654757)<self.EPS, "wrong distance.")
1716 s,d=f.getSideAndDistance([10.,10.,-3.], tag=2)
1717 self.failUnless( s<0, "wrong side.")
1718 self.failUnless( abs(d-3.*0.70710678118654757)<self.EPS, "wrong distance.")
1719
1720 s,d=f.getSideAndDistance([5.,12.,0], tag=2)
1721 self.failUnless( s>0, "wrong side.")
1722 self.failUnless( abs(d-2*0.70710678118654757)<self.EPS, "wrong distance.")
1723 s,d=f.getSideAndDistance([5.,12.,-1], tag=2)
1724 self.failUnless( s>0, "wrong side.")
1725 self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
1726 s,d=f.getSideAndDistance([5.,12.,-2], tag=2)
1727 # s not checked as it is undefined.
1728 self.failUnless( abs(d)<self.EPS, "wrong distance.")
1729 s,d=f.getSideAndDistance([5.,12.,-3], tag=2)
1730 self.failUnless( s<0, "wrong side.")
1731 self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
1732 s,d=f.getSideAndDistance([5.,12.,-4], tag=2)
1733 self.failUnless( s<0, "wrong side.")
1734 self.failUnless( abs(d-2.*0.70710678118654757)<self.EPS, "wrong distance.")
1735
1736 if __name__ == '__main__':
1737 suite = unittest.TestSuite()
1738 suite.addTest(unittest.makeSuite(Test_FaultSystem))
1739 suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
1740 suite.addTest(unittest.makeSuite(Test_Darcy3D))
1741 suite.addTest(unittest.makeSuite(Test_Darcy2D))
1742 suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
1743 # suite.addTest(Test_StokesProblemCartesian3D("test_PCG_P_large"))
1744 suite.addTest(unittest.makeSuite(Test_Mountains3D))
1745 suite.addTest(unittest.makeSuite(Test_Mountains2D))
1746 suite.addTest(unittest.makeSuite(Test_Rheologies))
1747 # suite.addTest(Test_IncompressibleIsotropicFlowCartesian("test_D2_Fixed_MuNone"))
1748 suite.addTest(unittest.makeSuite(Test_IncompressibleIsotropicFlowCartesian))
1749 s=unittest.TextTestRunner(verbosity=2).run(suite)
1750 if not s.wasSuccessful(): sys.exit(1)
1751

  ViewVC Help
Powered by ViewVC 1.1.26