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

  ViewVC Help
Powered by ViewVC 1.1.26