/[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 2793 - (show annotations)
Tue Dec 1 06:10:10 2009 UTC (9 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 86160 byte(s)
some new util functions
1
2 ########################################################
3 #
4 # Copyright (c) 2003-2009 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-2009 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-6
918 VERBOSE=False or True
919 A=1.
920 P_max=100
921 NE=2*getMPISizeWorld()*2
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*1.e-3)
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)
996 self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)
997 t+=dt
998 N_t+=1
999
1000 def check(self,N_t,mod,t_ref,v_ref, s_ref,p_ref):
1001 p=mod.getPressure()
1002 p-=integrate(p)/vol(self.dom)
1003 error_p=Lsup(mod.getPressure()-p_ref)/Lsup(p_ref)
1004 error_s=Lsup(mod.getDeviatoricStress()-s_ref)/Lsup(s_ref)
1005 error_v=Lsup(mod.getVelocity()-v_ref)/Lsup(v_ref)
1006 error_t=abs(mod.getTime()-t_ref)/abs(t_ref)
1007 if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v
1008 self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
1009 self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )
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 def tearDown(self):
1013 del self.dom
1014
1015 def test_D2_Fixed_MuNone_LateStart(self):
1016 self.dom = Rectangle(self.NE,self.NE,order=2)
1017 self.mu=None
1018 self.latestart=True
1019 self.runIt()
1020 def test_D2_Fixed_Mu_LateStart(self):
1021 self.dom = Rectangle(self.NE,self.NE,order=2)
1022 self.mu=555.
1023 self.latestart=True
1024 self.runIt()
1025 def test_D2_Fixed_MuNone(self):
1026 self.dom = Rectangle(self.NE,self.NE,order=2)
1027 self.mu=None
1028 self.latestart=False
1029 self.runIt()
1030 def test_D2_Fixed_Mu(self):
1031 self.dom = Rectangle(self.NE,self.NE,order=2)
1032 self.mu=555.
1033 self.latestart=False
1034 self.runIt()
1035 def test_D2_Free_MuNone_LateStart(self):
1036 self.dom = Rectangle(self.NE,self.NE,order=2)
1037 self.mu=None
1038 self.latestart=True
1039 self.runIt(free=0)
1040 def test_D2_Free_Mu_LateStart(self):
1041 self.dom = Rectangle(self.NE,self.NE,order=2)
1042 self.mu=555.
1043 self.latestart=True
1044 self.runIt(free=0)
1045 def test_D2_Free_MuNone(self):
1046 self.dom = Rectangle(self.NE,self.NE,order=2)
1047 self.mu=None
1048 self.latestart=False
1049 self.runIt(free=0)
1050 def test_D2_Free_Mu(self):
1051 self.dom = Rectangle(self.NE,self.NE,order=2)
1052 self.mu=555.
1053 self.latestart=False
1054 self.runIt(free=0)
1055
1056 def test_D3_Fixed_MuNone_LateStart(self):
1057 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1058 self.mu=None
1059 self.latestart=True
1060 self.runIt()
1061 def test_D3_Fixed_Mu_LateStart(self):
1062 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1063 self.mu=555.
1064 self.latestart=True
1065 self.runIt()
1066 def test_D3_Fixed_MuNone(self):
1067 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1068 self.mu=None
1069 self.latestart=False
1070 self.runIt()
1071 def test_D3_Fixed_Mu(self):
1072 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1073 self.mu=555.
1074 self.latestart=False
1075 self.runIt()
1076 def test_D3_Free_MuNone_LateStart(self):
1077 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1078 self.mu=None
1079 self.latestart=True
1080 self.runIt(free=0)
1081 def test_D3_Free_Mu_LateStart(self):
1082 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1083 self.mu=555.
1084 self.latestart=True
1085 self.runIt(free=0)
1086 def test_D3_Free_MuNone(self):
1087 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1088 self.mu=None
1089 self.latestart=False
1090 self.runIt(free=0)
1091 def test_D3_Free_Mu(self):
1092 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1093 self.mu=555.
1094 self.latestart=False
1095 self.runIt(free=0)
1096
1097
1098 class Test_FaultSystem(unittest.TestCase):
1099 EPS=1.e-8
1100 NE=10
1101 def test_Fault_MaxValue(self):
1102 dom=Rectangle(2*self.NE,2*self.NE)
1103 x=dom.getX()
1104 f=FaultSystem(dim=2)
1105 f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
1106 f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
1107
1108 u=x[0]*(1.-x[0])*(1-x[1])
1109 t, loc=f.getMaxValue(u)
1110 p=f.getParametrization(x,t)[0]
1111 m, l=loc(u), loc(p)
1112 self.failUnless( m == 0.25, "wrong max value")
1113 self.failUnless( t == 1, "wrong max tag")
1114 self.failUnless( l == 0., "wrong max location")
1115
1116 u=x[1]*(1.-x[1])*(1-x[0])*x[0]
1117 t, loc=f.getMaxValue(u)
1118 p=f.getParametrization(x,t)[0]
1119 m, l=loc(u), loc(p)
1120 self.failUnless( m == 0.0625, "wrong max value")
1121 self.failUnless( t == 2, "wrong max tag")
1122 self.failUnless( l == 0.5, "wrong max location")
1123
1124 u=x[0]*(1.-x[0])*x[1]
1125 t, loc=f.getMaxValue(u)
1126 p=f.getParametrization(x,t)[0]
1127 m, l=loc(u), loc(p)
1128 self.failUnless( m == 0.25, "wrong max value")
1129 self.failUnless( t == 2, "wrong max tag")
1130 self.failUnless( l == 1.0, "wrong max location")
1131
1132 u=x[1]*(1.-x[1])*x[0]
1133 t, loc=f.getMaxValue(u)
1134 p=f.getParametrization(x,t)[0]
1135 m, l=loc(u), loc(p)
1136 self.failUnless( m == 0.25, "wrong max value")
1137 self.failUnless( t == 2, "wrong max tag")
1138 self.failUnless( l == 0., "wrong max location")
1139
1140 u=x[1]*(1.-x[1])*(1.-x[0])
1141 t, loc=f.getMaxValue(u)
1142 p=f.getParametrization(x,t)[0]
1143 m, l=loc(u), loc(p)
1144 self.failUnless( m == 0.25, "wrong max value")
1145 self.failUnless( t == 1, "wrong max tag")
1146 self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong max location")
1147 def test_Fault_MinValue(self):
1148 dom=Rectangle(2*self.NE,2*self.NE)
1149 x=dom.getX()
1150 f=FaultSystem(dim=2)
1151 f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
1152 f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
1153
1154 u=-x[0]*(1.-x[0])*(1-x[1])
1155 t, loc=f.getMinValue(u)
1156 p=f.getParametrization(x,t)[0]
1157 m, l=loc(u), loc(p)
1158 self.failUnless( m == -0.25, "wrong min value")
1159 self.failUnless( t == 1, "wrong min tag")
1160 self.failUnless( l == 0., "wrong min location")
1161 u=-x[1]*(1.-x[1])*(1-x[0])*x[0]
1162 t, loc=f.getMinValue(u)
1163 p=f.getParametrization(x,t)[0]
1164 m, l=loc(u), loc(p)
1165 self.failUnless( m == -0.0625, "wrong min value")
1166 self.failUnless( t == 2, "wrong min tag")
1167 self.failUnless( l == 0.5, "wrong min location")
1168 u=-x[0]*(1.-x[0])*x[1]
1169 t, loc=f.getMinValue(u)
1170 p=f.getParametrization(x,t)[0]
1171 m, l=loc(u), loc(p)
1172 self.failUnless( m == -0.25, "wrong min value")
1173 self.failUnless( t == 2, "wrong min tag")
1174 self.failUnless( l == 1.0, "wrong min location")
1175 u=-x[1]*(1.-x[1])*x[0]
1176 t, loc=f.getMinValue(u)
1177 p=f.getParametrization(x,t)[0]
1178 m, l=loc(u), loc(p)
1179 self.failUnless( m == -0.25, "wrong min value")
1180 self.failUnless( t == 2, "wrong min tag")
1181 self.failUnless( l == 0., "wrong min location")
1182 u=-x[1]*(1.-x[1])*(1.-x[0])
1183 t, loc=f.getMinValue(u)
1184 p=f.getParametrization(x,t)[0]
1185 m, l=loc(u), loc(p)
1186 self.failUnless( m == -0.25, "wrong min value")
1187 self.failUnless( t == 1, "wrong min tag")
1188 self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong min location")
1189
1190
1191 def test_Fault2D(self):
1192 f=FaultSystem(dim=2)
1193 top1=[ [1.,0.], [1.,1.], [0.,1.] ]
1194 self.failUnlessRaises(ValueError,f.addFault,V0=[1.,0],strikes=[pi/2, pi/2],ls=[1.,1.],tag=1,dips=top1)
1195 f.addFault(V0=[1.,0],strikes=[pi/2, pi],ls=[1.,1.],tag=1)
1196 self.failUnless(f.getDim() == 2, "wrong dimension")
1197 self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
1198 self.failUnless( 2. == f.getTotalLength(1), "length wrong")
1199 self.failUnless( 0. == f.getMediumDepth(1), "depth wrong")
1200 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 range")
1201 self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 range")
1202 self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsets")
1203 segs=f.getTopPolyline(1)
1204 self.failUnless( len(segs) == 3, "wrong number of segments")
1205 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1206 self.failUnless( segs[0].size == 2, "seg 0 has wrong size.")
1207 self.failUnless( numpy.linalg.norm(segs[0]-[1.,0.]) < self.EPS, "wrong vertex. 0 ")
1208 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1209 self.failUnless( segs[1].size == 2, "seg 1 has wrong size.")
1210 self.failUnless( numpy.linalg.norm(segs[1]-[1.,1.]) < self.EPS, "wrong vertex. 1 ")
1211 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1212 self.failUnless( segs[2].size == 2, "seg 2 has wrong size.")
1213 self.failUnless( numpy.linalg.norm(segs[2]-[0.,1.]) < self.EPS, "wrong vertex. 2 ")
1214 c=f.getCenterOnSurface()
1215 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1216 self.failUnless( c.size == 2, "center size is wrong")
1217 self.failUnless( numpy.linalg.norm(c-[2./3.,2./3.]) < self.EPS, "center has wrong coordinates.")
1218 o=f.getOrientationOnSurface()/pi*180.
1219 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1220
1221 top2=[ [10.,0.], [0.,10.] ]
1222 f.addFault(V0=[10.,0],strikes=[3.*pi/4],ls=[14.142135623730951], tag=2, w0_offsets=[0,20], w1_max=20)
1223 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1224 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length")
1225 self.failUnless( 0. == f.getMediumDepth(2), "depth wrong")
1226 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range")
1227 self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range")
1228 self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets")
1229 segs=f.getTopPolyline(2)
1230 self.failUnless( len(segs) == 2, "wrong number of segments")
1231 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1232 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.]) < self.EPS, "wrong vertex. 0 ")
1233 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1234 self.failUnless( numpy.linalg.norm(segs[1]-[0.,10.]) < self.EPS, "wrong vertex. 1 ")
1235 c=f.getCenterOnSurface()
1236 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1237 self.failUnless( c.size == 2, "center size is wrong")
1238 self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.]) < self.EPS, "center has wrong coordinates.")
1239 o=f.getOrientationOnSurface()/pi*180.
1240 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1241
1242 s,d=f.getSideAndDistance([0.,0.], tag=1)
1243 self.failUnless( s<0, "wrong side.")
1244 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1245 s,d=f.getSideAndDistance([0.,2.], tag=1)
1246 self.failUnless( s>0, "wrong side.")
1247 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1248 s,d=f.getSideAndDistance([1.,2.], tag=1)
1249 self.failUnless( s>0, "wrong side.")
1250 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1251 s,d=f.getSideAndDistance([2.,1.], tag=1)
1252 self.failUnless( s>0, "wrong side.")
1253 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1254 s,d=f.getSideAndDistance([2.,0.], tag=1)
1255 self.failUnless( s>0, "wrong side.")
1256 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1257 s,d=f.getSideAndDistance([0.,-1.], tag=1)
1258 self.failUnless( s<0, "wrong side.")
1259 self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
1260 s,d=f.getSideAndDistance([-1.,0], tag=1)
1261 self.failUnless( s<0, "wrong side.")
1262 self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
1263
1264
1265 f.transform(rot=-pi/2., shift=[-1.,-1.])
1266 self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
1267 self.failUnless( 2. == f.getTotalLength(1), "length after transformation wrong")
1268 self.failUnless( 0. == f.getMediumDepth(1), "depth after transformation wrong")
1269 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 after transformation range")
1270 self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
1271 self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsetsafter transformation ")
1272 segs=f.getTopPolyline(1)
1273 self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
1274 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1275 self.failUnless( segs[0].size == 2, "seg 0 has wrong size after transformation.")
1276 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1277 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1278 self.failUnless( segs[1].size == 2, "seg 1 has wrong size after transformation.")
1279 self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1280 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1281 self.failUnless( segs[2].size == 2, "seg 2 has wrong size after transformation.")
1282 self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.]) < self.EPS, "wrong vertex after transformation. 2 ")
1283 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1284 self.failUnless( 0. == f.getMediumDepth(2), "depth wrong after transformation")
1285 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range after transformation")
1286 self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1287 self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1288 segs=f.getTopPolyline(2)
1289 self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
1290 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1291 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9]) < self.EPS, "wrong vertex. 0 after transformation")
1292 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1293 self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.]) < self.EPS, "wrong vertex. 1 after transformation")
1294
1295 c=f.getCenterOnSurface()
1296 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1297 self.failUnless( c.size == 2, "center size is wrong")
1298 self.failUnless( numpy.linalg.norm(c-[7./5.,-7./5.]) < self.EPS, "center has wrong coordinates.")
1299 o=f.getOrientationOnSurface()/pi*180.
1300 self.failUnless( abs(o-45.) < self.EPS, "wrong orientation.")
1301
1302 p=f.getParametrization([-1.,0.],1)
1303 self.failUnless(p[1]==1., "wrong value.")
1304 self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1305 p=f.getParametrization([-0.5,0.],1)
1306 self.failUnless(p[1]==1., "wrong value.")
1307 self.failUnless(abs(p[0]-0.5)<self.EPS* 0.5, "wrong value.")
1308 p=f.getParametrization([0.,0.],1)
1309 self.failUnless(p[1]==1., "wrong value.")
1310 self.failUnless(abs(p[0]-1.)<self.EPS, "wrong value.")
1311 p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-8)
1312 self.failUnless(p[1]==0., "wrong value.")
1313 p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-6)
1314 self.failUnless(p[1]==1., "wrong value.")
1315 self.failUnless(abs(p[0]-1.0000001)<self.EPS, "wrong value.")
1316 p=f.getParametrization([0.,0.5],1)
1317 self.failUnless(p[1]==1., "wrong value.")
1318 self.failUnless(abs(p[0]-1.5)<self.EPS, "wrong value.")
1319 p=f.getParametrization([0,1.],1)
1320 self.failUnless(p[1]==1., "wrong value.")
1321 self.failUnless(abs(p[0]-2.)<self.EPS, "wrong value.")
1322 p=f.getParametrization([1.,1.],1)
1323 self.failUnless(p[1]==0., "wrong value.")
1324 p=f.getParametrization([0,1.11],1)
1325 self.failUnless(p[1]==0., "wrong value.")
1326 p=f.getParametrization([-1,-9.],2)
1327 self.failUnless(p[1]==1., "wrong value.")
1328 self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1329 p=f.getParametrization([9,1],2)
1330 self.failUnless(p[1]==1., "wrong value.")
1331 self.failUnless(abs(p[0]-20.)<self.EPS, "wrong value.")
1332
1333 def test_Fault3D(self):
1334 f=FaultSystem(dim=3)
1335 self.failUnless(f.getDim() == 3, "wrong dimension")
1336
1337 top1=[ [0.,0.,0.], [1., 0., 0.] ]
1338 f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=20.,tag=1)
1339 self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
1340 self.failUnless( 1. == f.getTotalLength(1), "length wrong")
1341 self.failUnless( 20. == f.getMediumDepth(1), "depth wrong")
1342 self.failUnless( (0., 1.) == f.getW0Range(1)," wrong W0 range")
1343 self.failUnless( (-20., 0.) == f.getW1Range(1)," wrong W1 range")
1344 self.failUnless( [0., 1.] == f.getW0Offsets(1)," wrong W0 offsets")
1345 segs=f.getTopPolyline(1)
1346 self.failUnless( len(segs) == 2, "wrong number of segments")
1347 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1348 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1349 self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1350 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1351 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1352 self.failUnless( numpy.linalg.norm(segs[1]-[1.,0.,0]) < self.EPS, "wrong vertex. 1 ")
1353 c=f.getCenterOnSurface()
1354 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1355 self.failUnless( c.size == 3, "center size is wrong")
1356 self.failUnless( numpy.linalg.norm(c-[0.5,0.,0.]) < self.EPS, "center has wrong coordinates.")
1357 o=f.getOrientationOnSurface()/pi*180.
1358 self.failUnless( abs(o) < self.EPS, "wrong orientation.")
1359 d=f.getDips(1)
1360 self.failUnless( len(d) == 1, "wrong number of dips")
1361 self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1362 sn=f.getSegmentNormals(1)
1363 self.failUnless( len(sn) == 1, "wrong number of normals")
1364 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1365 self.failUnless( numpy.linalg.norm(sn[0]-[0, -1., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1366 dv=f.getDepthVectors(1)
1367 self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1368 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1369 self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1370 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1371 self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1372 b=f.getBottomPolyline(1)
1373 self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1374 self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1375 self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1376 self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1377 self.failUnless( numpy.linalg.norm(b[1]-[1., 0., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1378 ds=f.getDepths(1)
1379 self.failUnless( len(ds) == 2, "wrong number of depth")
1380 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1381 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1382
1383 top2=[ [0.,0.,0.], [0., 10., 0.] ]
1384 f.addFault(V0=[0.,0,0],strikes=[pi/2],ls=[10.,], dips=pi/2, depths=20.,tag=2)
1385 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1386 self.failUnless( 10. == f.getTotalLength(2), "length wrong")
1387 self.failUnless( 20. == f.getMediumDepth(2), "depth wrong")
1388 self.failUnless( (0., 10.) == f.getW0Range(2)," wrong W0 range")
1389 self.failUnless( (-20., 0.) == f.getW1Range(2)," wrong W1 range")
1390 self.failUnless( [0., 10.] == f.getW0Offsets(2)," wrong W0 offsets")
1391 segs=f.getTopPolyline(2)
1392 self.failUnless( len(segs) == 2, "wrong number of segments")
1393 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1394 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1395 self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1396 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1397 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1398 self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1399 d=f.getDips(2)
1400 self.failUnless( len(d) == 1, "wrong number of dips")
1401 self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1402 sn=f.getSegmentNormals(2)
1403 self.failUnless( len(sn) == 1, "wrong number of normals")
1404 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1405 self.failUnless( numpy.linalg.norm(sn[0]-[1, 0., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1406 dv=f.getDepthVectors(2)
1407 self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1408 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1409 self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1410 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1411 self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1412 b=f.getBottomPolyline(2)
1413 self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1414 self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1415 self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1416 self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1417 self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1418 ds=f.getDepths(2)
1419 self.failUnless( len(ds) == 2, "wrong number of depth")
1420 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1421 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1422
1423 top2=[ [10.,0.,0.], [0., 10., 0.] ]
1424 f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/2, depths=30.,tag=2)
1425 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1426 self.failUnless( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1427 self.failUnless( 30. == f.getMediumDepth(2), "depth wrong")
1428 self.failUnless( (-30., 0.) == f.getW1Range(2)," wrong W1 range")
1429 segs=f.getTopPolyline(2)
1430 self.failUnless( len(segs) == 2, "wrong number of segments")
1431 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1432 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1433 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1434 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1435 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1436 self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1437 d=f.getDips(2)
1438 self.failUnless( len(d) == 1, "wrong number of dips")
1439 self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1440 sn=f.getSegmentNormals(2)
1441 self.failUnless( len(sn) == 1, "wrong number of normals")
1442 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1443 self.failUnless( numpy.linalg.norm(sn[0]-[0.70710678118654746, 0.70710678118654746, 0.]) < self.EPS, "wrong bottom vertex 1 ")
1444 dv=f.getDepthVectors(2)
1445 self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1446 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1447 self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -30.]) < self.EPS, "wrong depth vector 0 ")
1448 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1449 self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -30.]) < self.EPS, "wrong depth vector 1 ")
1450 b=f.getBottomPolyline(2)
1451 self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1452 self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1453 self.failUnless( numpy.linalg.norm(b[0]-[10., 0., -30.]) < self.EPS, "wrong bottom vertex 0 ")
1454 self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1455 self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -30.]) < self.EPS, "wrong bottom vertex 1 ")
1456 ds=f.getDepths(2)
1457 self.failUnless( len(ds) == 2, "wrong number of depth")
1458 self.failUnless( abs(ds[0]-30.) < self.EPS, "wrong depth at vertex 0 ")
1459 self.failUnless( abs(ds[1]-30.) < self.EPS, "wrong depth at vertex 1 ")
1460
1461 top2=[ [10.,0.,0.], [0., 10., 0.] ]
1462 f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/4, depths=50.,tag=2)
1463 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1464 self.failUnless( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1465 self.failUnless( 50. == f.getMediumDepth(2), "depth wrong")
1466 self.failUnless( (-50., 0.) == f.getW1Range(2)," wrong W1 range")
1467 segs=f.getTopPolyline(2)
1468 self.failUnless( len(segs) == 2, "wrong number of segments")
1469 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1470 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1471 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1472 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1473 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1474 self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1475 d=f.getDips(2)
1476 self.failUnless( len(d) == 1, "wrong number of dips")
1477 self.failUnless( abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1478 sn=f.getSegmentNormals(2)
1479 self.failUnless( len(sn) == 1, "wrong number of normals")
1480 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1481 self.failUnless( numpy.linalg.norm(sn[0]-[0.5,0.5,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1482 dv=f.getDepthVectors(2)
1483 self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1484 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1485 self.failUnless( numpy.linalg.norm(dv[0]-[25., 25., -35.355339059327378]) < self.EPS, "wrong depth vector 0 ")
1486 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1487 self.failUnless( numpy.linalg.norm(dv[1]-[25.,25., -35.355339059327378]) < self.EPS, "wrong depth vector 1 ")
1488 b=f.getBottomPolyline(2)
1489 self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1490 self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1491 self.failUnless( numpy.linalg.norm(b[0]-[35., 25., -35.355339059327378]) < self.EPS, "wrong bottom vertex 0 ")
1492 self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1493 self.failUnless( numpy.linalg.norm(b[1]-[25, 35., -35.355339059327378]) < self.EPS, "wrong bottom vertex 1 ")
1494 ds=f.getDepths(2)
1495 self.failUnless( len(ds) == 2, "wrong number of depth")
1496 self.failUnless( abs(ds[0]-50.) < self.EPS, "wrong depth at vertex 0 ")
1497 self.failUnless( abs(ds[1]-50.) < self.EPS, "wrong depth at vertex 1 ")
1498
1499 top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
1500 f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=1, dips=pi/4, depths=20.)
1501 self.failUnless( 20. == f.getTotalLength(1), "length wrong")
1502 self.failUnless( 20. == f.getMediumDepth(1), "depth wrong")
1503 segs=f.getTopPolyline(1)
1504 self.failUnless( len(segs) == 3, "wrong number of segments")
1505 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1506 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1507 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1508 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1509 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1510 self.failUnless( numpy.linalg.norm(segs[1]-[10.,10.,0.]) < self.EPS, "wrong vertex. 1 ")
1511 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1512 self.failUnless( segs[2].size == 3, "seg 2 has wrong size.")
1513 self.failUnless( numpy.linalg.norm(segs[2]-[0.,10.,0.]) < self.EPS, "wrong vertex. 2 ")
1514 d=f.getDips(1)
1515 self.failUnless( len(d) == 2, "wrong number of dips")
1516 self.failUnless( abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1517 self.failUnless( abs(d[1]-0.78539816339744828) < self.EPS, "wrong dip 0")
1518 ds=f.getDepths(1)
1519 self.failUnless( len(ds) == 3, "wrong number of depth")
1520 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1521 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1522 sn=f.getSegmentNormals(1)
1523 self.failUnless( len(sn) == 2, "wrong number of normals")
1524 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1525 self.failUnless( numpy.linalg.norm(sn[0]-[0.70710678118654746,0.,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1526 self.failUnless( isinstance(sn[1], numpy.ndarray), "wrong class of bottom vertex 0")
1527 self.failUnless( numpy.linalg.norm(sn[1]-[0.,0.70710678118654746,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1528 dv=f.getDepthVectors(1)
1529 self.failUnless( len(dv) == 3, "wrong number of depth vectors.")
1530 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1531 self.failUnless( numpy.linalg.norm(dv[0]-[14.142135623730951, 0., -14.142135623730951]) < self.EPS, "wrong depth vector 0 ")
1532 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1533 self.failUnless( numpy.linalg.norm(dv[1]-[11.547005383792515,11.547005383792515, -11.547005383792515]) < self.EPS, "wrong depth vector 2 ")
1534 self.failUnless( isinstance(dv[2], numpy.ndarray), "wrong class of depth vector 1")
1535 self.failUnless( numpy.linalg.norm(dv[2]-[0.,14.142135623730951, -14.142135623730951]) < self.EPS, "wrong depth vector 2 ")
1536 segs=f.getBottomPolyline(1)
1537 self.failUnless( len(segs) == 3, "wrong number of segments")
1538 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1539 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1540 self.failUnless( numpy.linalg.norm(segs[0]-[24.142135623730951,0.,-14.142135623730951]) < self.EPS, "wrong vertex. 0 ")
1541 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1542 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1543 self.failUnless( numpy.linalg.norm(segs[1]-[21.547005383792515,21.547005383792515, -11.547005383792515]) < self.EPS, "wrong vertex. 1 ")
1544 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1545 self.failUnless( segs[2].size == 3, "seg 2 has wrong size.")
1546 self.failUnless( numpy.linalg.norm(segs[2]-[0., 24.142135623730951, -14.142135623730951]) < self.EPS, "wrong vertex. 2 ")
1547 self.failUnless( abs(0.-f.getW0Range(1)[0]) <=self.EPS," wrong W0 range (0)")
1548 self.failUnless( abs(31.857329272664341-f.getW0Range(1)[1]) <=self.EPS," wrong W0 range (1)")
1549 self.failUnless( abs(-20. - f.getW1Range(1)[0]) <=self.EPS," wrong W1 range (0)")
1550 self.failUnless( abs(0. - f.getW1Range(1)[1]) <=self.EPS," wrong W1 range (1)")
1551 self.failUnless( abs(0.0-f.getW0Offsets(1)[0])<=self.EPS," wrong W0 offsets (0)")
1552 self.failUnless( abs(15.92866463633217-f.getW0Offsets(1)[1])<=self.EPS," wrong W0 offsets (1)")
1553 self.failUnless( abs(31.857329272664341-f.getW0Offsets(1)[2])<=self.EPS," wrong W0 offsets(2)")
1554 #
1555 # ============ fresh start ====================
1556 #
1557 f.addFault(V0=[1.,0,0.],strikes=[pi/2, pi],ls=[1.,1.], dips=pi/2,depths=20,tag=1)
1558 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)
1559 c=f.getCenterOnSurface()
1560 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1561 self.failUnless( c.size == 3, "center size is wrong")
1562 self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.,0.]) < self.EPS, "center has wrong coordinates.")
1563 o=f.getOrientationOnSurface()/pi*180.
1564 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1565
1566 f.transform(rot=-pi/2., shift=[-1.,-1.,0.])
1567 self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
1568 self.failUnless( 2. == f.getTotalLength(1), "length after transformation wrong")
1569 self.failUnless( 20. == f.getMediumDepth(1), "depth after transformation wrong")
1570 rw0=f.getW0Range(1)
1571 self.failUnless( len(rw0) ==2, "wo range has wrong length")
1572 self.failUnless( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
1573 self.failUnless( abs(rw0[1]-2.) < self.EPS,"W0 1 wrong.")
1574 self.failUnless( (-20., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
1575 dips=f.getDips(1)
1576 self.failUnless(len(dips) == 2, "wrong number of dips.")
1577 self.failUnless( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
1578 self.failUnless( abs(dips[1]-1.5707963267948966) <= self.EPS, "wrong dip")
1579 ds=f.getDepths(1)
1580 self.failUnless( len(ds) == 3, "wrong number of depth")
1581 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1582 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1583 self.failUnless( abs(ds[2]-20.) < self.EPS, "wrong depth at vertex 1 ")
1584 segs=f.getTopPolyline(1)
1585 self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
1586 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1587 self.failUnless( segs[0].size == 3, "seg 0 has wrong size after transformation.")
1588 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1589 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1590 self.failUnless( segs[1].size == 3, "seg 1 has wrong size after transformation.")
1591 self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1592 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1593 self.failUnless( segs[2].size == 3, "seg 2 has wrong size after transformation.")
1594 self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.,0.]) < self.EPS, "wrong vertex after transformation. 2 ")
1595 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1596 self.failUnless( 20. == f.getMediumDepth(2), "depth wrong after transformation")
1597 rw0=f.getW0Range(2)
1598 self.failUnless( len(rw0) ==2, "wo range has wrong length")
1599 self.failUnless( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
1600 self.failUnless( abs(rw0[1]-20.) < self.EPS,"W0 1 wrong.")
1601 self.failUnless( (-20., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1602 self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1603 dips=f.getDips(2)
1604 self.failUnless(len(dips) == 1, "wrong number of dips.")
1605 self.failUnless( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
1606 ds=f.getDepths(2)
1607 self.failUnless( len(ds) == 2, "wrong number of depth")
1608 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1609 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1610 segs=f.getTopPolyline(2)
1611 self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
1612 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1613 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1614 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1615 self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.,0.]) < self.EPS, "wrong vertex. 1 after transformation")
1616 #
1617 # ============ fresh start ====================
1618 #
1619 f=FaultSystem(dim=3)
1620
1621 top1=[ [0.,0.,0.], [1., 0., 0.] ]
1622 f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=1.,tag=1)
1623 top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
1624 f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=2, dips=pi/4, depths=20.)
1625
1626 p,m=f.getParametrization([0.3,0.,-0.5],1)
1627 self.failUnless(length(p-[0.3,-0.5]) <= self.EPS, "wrong value.")
1628 self.failUnless(m==1., "wrong value.")
1629
1630 p,m=f.getParametrization([0.5,0.,-0.5],1)
1631 self.failUnless(length(p-[0.5,-0.5]) <= self.EPS, "wrong value.")
1632 self.failUnless(m==1., "wrong value.")
1633
1634 p,m=f.getParametrization([0.25,0.,-0.5],1)
1635 self.failUnless(length(p-[0.25,-0.5]) <= self.EPS, "wrong value.")
1636 self.failUnless(m==1., "wrong value.")
1637
1638 p,m=f.getParametrization([0.5,0.,-0.25],1)
1639 self.failUnless(length(p-[0.5,-0.25]) <= self.EPS, "wrong value.")
1640 self.failUnless(m==1., "wrong value.")
1641
1642 p,m=f.getParametrization([0.001,0.,-0.001],1)
1643 self.failUnless(length(p-[0.001, -0.001]) <= self.EPS, "wrong value.")
1644 self.failUnless(m==1., "wrong value.")
1645
1646 p,m=f.getParametrization([0.001,0.,0.001],1)
1647 self.failUnless(m==0., "wrong value.")
1648
1649 p,m=f.getParametrization([0.999,0.,0.001],1)
1650 self.failUnless(m==0., "wrong value.")
1651
1652 p,m=f.getParametrization([1.001,0.,-0.001],1)
1653 self.failUnless(m==0., "wrong value.")
1654 p,m=f.getParametrization([1.001,0.,-0.1],1)
1655 self.failUnless(m==0., "wrong value.")
1656 p,m=f.getParametrization([1.001,0.,-0.000000001],1)
1657 self.failUnless(m==0., "wrong value.")
1658
1659 p,m=f.getParametrization([0.999,0.,-0.001],1)
1660 self.failUnless(length(p-[0.999, -0.001]) <= self.EPS, "wrong value.")
1661 self.failUnless(m==1., "wrong value.")
1662
1663 p,m=f.getParametrization([ 16.29252873 , 6.46410161 ,-6.29252873], 2, tol=1.e-7)
1664 self.failUnless(m==1., "wrong value.")
1665 self.failUnless(length(p-[4.7785993908996511, -10]) <= self.EPS*10., "wrong value.")
1666 p,m=f.getParametrization([15.77350269, 12.77350269, -5.77350269], 2, tol=1.e-7)
1667 self.failUnless(m==1., "wrong value.")
1668 self.failUnless(length(p-[11.150065245432518, -10]) <= self.EPS*10., "wrong value.")
1669
1670 p,m=f.getParametrization([ 3., 17.0710678, -7.0710678], 2, tol=1.e-7)
1671 self.failUnless(m==1., "wrong value.")
1672 self.failUnless(length(p-[27.078729881764687, -10]) <= self.EPS*10., "wrong value.")
1673 p,m=f.getParametrization([9.30940108, 16.55204176, -6.55204176], 2, tol=1.e-7)
1674 self.failUnless(m==1., "wrong value.")
1675 self.failUnless(length(p-[20.707264027231822, -10]) <= self.EPS*10., "wrong value.")
1676
1677 p,m=f.getParametrization([ 21.54700538, 21.54700538, -11.54700538], 2, tol=1.e-7)
1678 self.failUnless(m==1., "wrong value.")
1679 self.failUnless(length(p-[15.92866463633217, -20]) <= self.EPS*10., "wrong value.")
1680
1681 p,m=f.getParametrization([ 0.,0.,0.], 2, tol=1.e-7)
1682 self.failUnless(m==0., "wrong value.")
1683
1684 p,m=f.getParametrization([ 11.,11.,0.], 2, tol=1.e-7)
1685 self.failUnless(m==0., "wrong value.")
1686
1687
1688 s,d=f.getSideAndDistance([0.,-1.,0.], tag=1)
1689 self.failUnless( s>0, "wrong side.")
1690 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1691 s,d=f.getSideAndDistance([1.,-1.,0.], tag=1)
1692 self.failUnless( s>0, "wrong side.")
1693 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1694 s,d=f.getSideAndDistance([0.,1.,0.], tag=1)
1695 self.failUnless( s<0, "wrong side.")
1696 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1697 s,d=f.getSideAndDistance([1.,1.,0.], tag=1)
1698 self.failUnless( s<0, "wrong side.")
1699 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1700
1701
1702 s,d=f.getSideAndDistance([0.,0.,0.], tag=2)
1703 self.failUnless( s<0, "wrong side.")
1704 self.failUnless( abs(d-10.)<self.EPS, "wrong distance.")
1705 s,d=f.getSideAndDistance([5.,5.,0.], tag=2)
1706 self.failUnless( s<0, "wrong side.")
1707 self.failUnless( abs(d-5.)<self.EPS, "wrong distance.")
1708
1709 s,d=f.getSideAndDistance([10.,10.,-1.], tag=2)
1710 self.failUnless( s<0, "wrong side.")
1711 self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
1712 s,d=f.getSideAndDistance([10.,10.,-2.], tag=2)
1713 self.failUnless( s<0, "wrong side.")
1714 self.failUnless( abs(d-2.*0.70710678118654757)<self.EPS, "wrong distance.")
1715 s,d=f.getSideAndDistance([10.,10.,-3.], tag=2)
1716 self.failUnless( s<0, "wrong side.")
1717 self.failUnless( abs(d-3.*0.70710678118654757)<self.EPS, "wrong distance.")
1718
1719 s,d=f.getSideAndDistance([5.,12.,0], tag=2)
1720 self.failUnless( s>0, "wrong side.")
1721 self.failUnless( abs(d-2*0.70710678118654757)<self.EPS, "wrong distance.")
1722 s,d=f.getSideAndDistance([5.,12.,-1], tag=2)
1723 self.failUnless( s>0, "wrong side.")
1724 self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
1725 s,d=f.getSideAndDistance([5.,12.,-2], tag=2)
1726 # s not checked as it is undefined.
1727 self.failUnless( abs(d)<self.EPS, "wrong distance.")
1728 s,d=f.getSideAndDistance([5.,12.,-3], tag=2)
1729 self.failUnless( s<0, "wrong side.")
1730 self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
1731 s,d=f.getSideAndDistance([5.,12.,-4], tag=2)
1732 self.failUnless( s<0, "wrong side.")
1733 self.failUnless( abs(d-2.*0.70710678118654757)<self.EPS, "wrong distance.")
1734
1735 if __name__ == '__main__':
1736 suite = unittest.TestSuite()
1737 # suite.addTest(unittest.makeSuite(Test_FaultSystem))
1738 # suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
1739 # suite.addTest(unittest.makeSuite(Test_Darcy3D))
1740 # suite.addTest(unittest.makeSuite(Test_Darcy2D))
1741 # suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
1742 # suite.addTest(unittest.makeSuite(Test_Mountains3D))
1743 # suite.addTest(unittest.makeSuite(Test_Mountains2D))
1744 # suite.addTest(unittest.makeSuite(Test_Rheologies))
1745 suite.addTest(Test_IncompressibleIsotropicFlowCartesian("test_D2_Fixed_Mu"))
1746 # suite.addTest(unittest.makeSuite(Test_IncompressibleIsotropicFlowCartesian))
1747 s=unittest.TextTestRunner(verbosity=2).run(suite)
1748 if not s.wasSuccessful(): sys.exit(1)
1749

  ViewVC Help
Powered by ViewVC 1.1.26