/[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 3501 - (show annotations)
Wed Apr 13 04:07:53 2011 UTC (8 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 86914 byte(s)


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