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

  ViewVC Help
Powered by ViewVC 1.1.26