/[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 3071 - (show annotations)
Wed Jul 21 05:37:30 2010 UTC (9 years, 2 months ago) by gross
File MIME type: text/x-python
File size: 86497 byte(s)
new darcy flux
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 and 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 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
523 self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
524
525 def testVarioF_FixedBottom_mediumK(self):
526 k=1.
527 mp=self.getScalarMask(include_bottom=True)
528 mv=self.getVectorMask(include_bottom=False)
529 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
530 p=p_ref*mp
531 u=u_ref*mv
532 df=DarcyFlow(self.dom)
533 df.setValue(g=f,
534 location_of_fixed_pressure=mp,
535 location_of_fixed_flux=mv,
536 permeability=Scalar(k,Function(self.dom)))
537 df.setTolerance(rtol=self.TOL)
538 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
539 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
540 self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
541
542 def testVarioF_FixedBottom_largeK(self):
543 k=1.e10
544 mp=self.getScalarMask(include_bottom=True)
545 mv=self.getVectorMask(include_bottom=False)
546 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
547 p=p_ref*mp
548 u=u_ref*mv
549 df=DarcyFlow(self.dom)
550 df.setValue(g=f,
551 location_of_fixed_pressure=mp,
552 location_of_fixed_flux=mv,
553 permeability=Scalar(k,Function(self.dom)))
554 df.setTolerance(rtol=self.TOL)
555 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
556 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
557 self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
558
559 def testConstF_FreeBottom_smallK(self):
560 k=1.e-10
561 mp=self.getScalarMask(include_bottom=False)
562 mv=self.getVectorMask(include_bottom=True)
563 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
564 p=p_ref*mp
565 u=u_ref*mv
566 df=DarcyFlow(self.dom)
567 df.setValue(g=f,
568 location_of_fixed_pressure=mp,
569 location_of_fixed_flux=mv,
570 permeability=Scalar(k,Function(self.dom)))
571 df.setTolerance(rtol=self.TOL)
572 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
573
574
575 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
576 self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
577
578 def testConstF_FreeBottom_mediumK(self):
579 k=1.
580 mp=self.getScalarMask(include_bottom=False)
581 mv=self.getVectorMask(include_bottom=True)
582 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
583 p=p_ref*mp
584 u=u_ref*mv
585 df=DarcyFlow(self.dom)
586 df.setValue(g=f,
587 location_of_fixed_pressure=mp,
588 location_of_fixed_flux=mv,
589 permeability=Scalar(k,Function(self.dom)))
590 df.setTolerance(rtol=self.TOL)
591 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
592 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
593 self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
594
595 def testConstF_FreeBottom_largeK(self):
596 k=1.e10
597 mp=self.getScalarMask(include_bottom=False)
598 mv=self.getVectorMask(include_bottom=True)
599 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
600 p=p_ref*mp
601 u=u_ref*mv
602 df=DarcyFlow(self.dom)
603 df.setValue(g=f,
604 location_of_fixed_pressure=mp,
605 location_of_fixed_flux=mv,
606 permeability=Scalar(k,Function(self.dom)))
607 df.setTolerance(rtol=self.TOL)
608 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
609 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
610 self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
611
612 def testVarioF_FreeBottom_smallK(self):
613 k=1.e-10
614 mp=self.getScalarMask(include_bottom=False)
615 mv=self.getVectorMask(include_bottom=True)
616 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
617 p=p_ref*mp
618 u=u_ref*mv
619 df=DarcyFlow(self.dom)
620 df.setValue(g=f,
621 location_of_fixed_pressure=mp,
622 location_of_fixed_flux=mv,
623 permeability=Scalar(k,Function(self.dom)))
624 df.setTolerance(rtol=self.TOL)
625 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
626 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
627 self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
628
629 def testVarioF_FreeBottom_mediumK(self):
630 k=1.
631 mp=self.getScalarMask(include_bottom=False)
632 mv=self.getVectorMask(include_bottom=True)
633 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
634 p=p_ref*mp
635 u=u_ref*mv
636 df=DarcyFlow(self.dom)
637 df.setValue(g=f,
638 location_of_fixed_pressure=mp,
639 location_of_fixed_flux=mv,
640 permeability=Scalar(k,Function(self.dom)))
641 df.setTolerance(rtol=self.TOL)
642 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
643 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
644 self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
645
646 def testVarioF_FreeBottom_largeK(self):
647 k=1.e10
648 mp=self.getScalarMask(include_bottom=False)
649 mv=self.getVectorMask(include_bottom=True)
650 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
651 p=p_ref*mp
652 u=u_ref*mv
653 df=DarcyFlow(self.dom)
654 df.setValue(g=f,
655 location_of_fixed_pressure=mp,
656 location_of_fixed_flux=mv,
657 permeability=Scalar(k,Function(self.dom)))
658 df.setTolerance(rtol=self.TOL)
659 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
660 self.failUnless(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
661 self.failUnless(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
662
663 class Test_Darcy2D(Test_Darcy):
664 TOL=1e-6
665 TEST_TOL=2.e-3
666 WIDTH=1.
667 def setUp(self):
668 NE=40 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
669 self.dom = Rectangle(NE/2,NE)
670 self.rescaleDomain()
671 def tearDown(self):
672 del self.dom
673 class Test_Darcy3D(Test_Darcy):
674 TOL=1e-6
675 WIDTH=1.
676 TEST_TOL=4.e-3
677 def setUp(self):
678 NE=25 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
679 self.dom = Brick(NE,NE,NE)
680 self.rescaleDomain()
681 def tearDown(self):
682 del self.dom
683
684
685 class Test_Mountains3D(unittest.TestCase):
686 def setUp(self):
687 NE=16
688 self.TOL=1e-4
689 self.domain=Brick(NE,NE,NE,order=1,useFullElementOrder=True)
690 def tearDown(self):
691 del self.domain
692 def test_periodic(self):
693 EPS=0.01
694
695 x=self.domain.getX()
696 v = Vector(0.0, Solution(self.domain))
697 a0=1
698 a1=1
699 n0=2
700 n1=2
701 n2=0.5
702 a2=-(a0*n0+a1*n1)/n2
703 v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])* cos(pi*n2*x[2])
704 v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])* cos(pi*n2*x[2])
705 v[2]=a2*cos(pi*n0*x[0])* cos(pi*n1*x[1])* sin(pi*n2*x[2])
706
707 mts=Mountains(self.domain,eps=EPS)
708 mts.setVelocity(v)
709 Z=mts.update()
710
711 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[2]-1.)))
712 self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
713
714 class Test_Mountains2D(unittest.TestCase):
715 def setUp(self):
716 NE=16
717 self.TOL=1e-4
718 self.domain=Rectangle(NE,NE,order=1,useFullElementOrder=True)
719 def tearDown(self):
720 del self.domain
721 def test_periodic(self):
722 EPS=0.01
723
724 x=self.domain.getX()
725 v = Vector(0.0, Solution(self.domain))
726 a0=1
727 n0=1
728 n1=0.5
729 a1=-(a0*n0)/n1
730 v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])
731 v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])
732
733 H_t=Scalar(0.0, Solution(self.domain))
734 mts=Mountains(self.domain,eps=EPS)
735 mts.setVelocity(v)
736 Z=mts.update()
737
738 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[1]-1.)))
739 self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
740
741
742
743 class Test_Rheologies(unittest.TestCase):
744 """
745 this is the program used to generate the powerlaw tests:
746
747 TAU_Y=100.
748 N=10
749 M=5
750
751 def getE(tau):
752 if tau<=TAU_Y:
753 return 1./(0.5+20*sqrt(tau))
754 else:
755 raise ValueError,"out of range."
756 tau=[ i* (TAU_Y/N) for i in range(N+1)] + [TAU_Y for i in range(M)]
757 e=[ tau[i]/getE(tau[i]) for i in range(N+1)]
758 e+= [ (i+1+j)* (max(e)/N) for j in range(M) ]
759
760 print tau
761 print e
762 """
763 TOL=1.e-8
764 def test_PowerLaw_Init(self):
765 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
766
767 self.failUnlessEqual(pl.getNumMaterials(),3,"num materials is wrong")
768 self.failUnlessEqual(pl.validMaterialId(0),True,"material id 0 not found")
769 self.failUnlessEqual(pl.validMaterialId(1),True,"material id 1 not found")
770 self.failUnlessEqual(pl.validMaterialId(2),True,"material id 2 not found")
771 self.failUnlessEqual(pl.validMaterialId(3),False,"material id 3 not found")
772
773 self.failUnlessEqual(pl.getFriction(),None,"initial friction wrong.")
774 self.failUnlessEqual(pl.getTauY(),None,"initial tau y wrong.")
775 pl.setDruckerPragerLaw(tau_Y=10,friction=3)
776 self.failUnlessEqual(pl.getFriction(),3,"friction wrong.")
777 self.failUnlessEqual(pl.getTauY(),10,"tau y wrong.")
778
779 self.failUnlessEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
780 pl.setElasticShearModulus(1000)
781 self.failUnlessEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
782
783 e=pl.getEtaN()
784 self.failUnlessEqual(len(e),3,"initial length of etaN is wrong.")
785 self.failUnlessEqual(e,[None, None, None],"initial etaN are wrong.")
786 self.failUnlessEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
787 self.failUnlessEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
788 self.failUnlessEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
789 self.failUnlessRaises(ValueError, pl.getEtaN, 3)
790
791 self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
792 self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
793 self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
794
795 pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
796 self.failUnlessEqual(pl.getPower(),[1,2,3],"powers are wrong.")
797 self.failUnlessEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
798 self.failUnlessEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
799
800 pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
801 self.failUnlessEqual(pl.getPower(),[4,2,3],"powers are wrong.")
802 self.failUnlessEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
803 self.failUnlessEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
804
805 self.failUnlessRaises(ValueError,pl.getPower,-1)
806 self.failUnlessEqual(pl.getPower(0),4,"power 0 is wrong.")
807 self.failUnlessEqual(pl.getPower(1),2,"power 1 is wrong.")
808 self.failUnlessEqual(pl.getPower(2),3,"power 2 is wrong.")
809 self.failUnlessRaises(ValueError,pl.getPower,3)
810
811 self.failUnlessRaises(ValueError,pl.getTauT,-1)
812 self.failUnlessEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
813 self.failUnlessEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
814 self.failUnlessEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
815 self.failUnlessRaises(ValueError,pl.getTauT,3)
816
817 self.failUnlessRaises(ValueError,pl.getEtaN,-1)
818 self.failUnlessEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
819 self.failUnlessEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
820 self.failUnlessEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
821 self.failUnlessRaises(ValueError,pl.getEtaN,3)
822
823 def checkResult(self,id,gamma_dot_, eta, tau_ref):
824 self.failUnless(eta>=0,"eta needs to be positive (test %s)"%id)
825 error=abs(gamma_dot_*eta-tau_ref)
826 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))
827
828 def test_PowerLaw_Linear(self):
829 taus= [0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
830 gamma_dot_s=[0.0, 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]
831 pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
832 pl.setDruckerPragerLaw(tau_Y=100.)
833 pl.setPowerLaw(eta_N=2.)
834 pl.setEtaTolerance(self.TOL)
835 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
836
837 def test_PowerLaw_QuadLarge(self):
838 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
839 gamma_dot_s=[0.0, 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]
840 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
841 pl.setDruckerPragerLaw(tau_Y=100.)
842 pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2])
843 pl.setEtaTolerance(self.TOL)
844 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
845
846 def test_PowerLaw_QuadSmall(self):
847 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
848 gamma_dot_s=[0.0, 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]
849 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
850 pl.setDruckerPragerLaw(tau_Y=100.)
851 pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2])
852 pl.setEtaTolerance(self.TOL)
853 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
854
855 def test_PowerLaw_CubeLarge(self):
856 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
857 gamma_dot_s=[0.0, 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]
858 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
859 pl.setDruckerPragerLaw(tau_Y=100.)
860 pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3])
861 pl.setEtaTolerance(self.TOL)
862 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
863
864 def test_PowerLaw_CubeSmall(self):
865 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
866 gamma_dot_s=[0.0, 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]
867 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
868 pl.setDruckerPragerLaw(tau_Y=100.)
869 pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3])
870 pl.setEtaTolerance(self.TOL)
871 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
872
873 def test_PowerLaw_QuadLarge_CubeLarge(self):
874 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]
875 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]
876
877 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
878 pl.setDruckerPragerLaw(tau_Y=100.)
879 pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
880 pl.setEtaTolerance(self.TOL)
881 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
882
883 def test_PowerLaw_QuadLarge_CubeSmall(self):
884 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]
885 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]
886
887 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
888 pl.setDruckerPragerLaw(tau_Y=100.)
889 pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
890 pl.setEtaTolerance(self.TOL)
891 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
892
893 def test_PowerLaw_QuadSmall_CubeLarge(self):
894 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
895 gamma_dot_s=[0.0, 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]
896 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
897 pl.setDruckerPragerLaw(tau_Y=100.)
898 pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
899 pl.setEtaTolerance(self.TOL)
900 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
901
902 def test_PowerLaw_QuadSmall_CubeSmall(self):
903 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
904 gamma_dot_s=[0.0, 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]
905 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
906 pl.setDruckerPragerLaw(tau_Y=100.)
907 pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
908 pl.setEtaTolerance(self.TOL)
909 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
910
911 def test_PowerLaw_withShear(self):
912 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]
913 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]
914 pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
915 pl.setDruckerPragerLaw(tau_Y=100.)
916 pl.setPowerLaw(eta_N=2.)
917 pl.setElasticShearModulus(3.)
918 dt=1./3.
919 pl.setEtaTolerance(self.TOL)
920 self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
921 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
922
923 class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase):
924 TOL=1.e-5
925 VERBOSE=False # or True
926 A=1.
927 P_max=100
928 NE=2*getMPISizeWorld()
929 tau_Y=10.
930 N_dt=10
931
932 # material parameter:
933 tau_1=5.
934 tau_2=5.
935 eta_0=100.
936 eta_1=50.
937 eta_2=400.
938 N_1=2.
939 N_2=3.
940 def getReference(self, t):
941
942 B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
943 x=self.dom.getX()
944
945 s_00=min(self.A*t,B)
946 tau=sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)*abs(s_00)
947 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.)
948
949 alpha=0.5*inv_eta*s_00
950 if s_00 <= B and self.mu !=None: alpha+=1./(2*self.mu)*self.A
951 u_ref=x*alpha
952 u_ref[self.dom.getDim()-1]=(1.-x[self.dom.getDim()-1])*alpha*(self.dom.getDim()-1)
953 sigma_ref=kronecker(self.dom)*s_00
954 sigma_ref[self.dom.getDim()-1,self.dom.getDim()-1]=-s_00*(self.dom.getDim()-1)
955
956 p_ref=self.P_max
957 for d in range(self.dom.getDim()): p_ref=p_ref*x[d]
958 p_ref-=integrate(p_ref)/vol(self.dom)
959 return u_ref, sigma_ref, p_ref
960
961 def runIt(self, free=None):
962 x=self.dom.getX()
963 B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
964 dt=B/int(self.N_dt/2)
965 if self.VERBOSE: print "dt =",dt
966 if self.latestart:
967 t=dt
968 else:
969 t=0
970 v,s,p=self.getReference(t)
971
972 mod=IncompressibleIsotropicFlowCartesian(self.dom, stress=s, v=v, p=p, t=t, numMaterials=3, verbose=self.VERBOSE)
973 mod.setDruckerPragerLaw(tau_Y=self.tau_Y,friction=None)
974 mod.setElasticShearModulus(self.mu)
975 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])
976 mod.setTolerance(self.TOL)
977 mod.setEtaTolerance(self.TOL*0.1)
978
979 BF=Vector(self.P_max,Function(self.dom))
980 for d in range(self.dom.getDim()):
981 for d2 in range(self.dom.getDim()):
982 if d!=d2: BF[d]*=x[d2]
983 v_mask=Vector(0,Solution(self.dom))
984 if free==None:
985 for d in range(self.dom.getDim()):
986 v_mask[d]=whereZero(x[d])+whereZero(x[d]-1.)
987 else:
988 for d in range(self.dom.getDim()):
989 if d == self.dom.getDim()-1:
990 v_mask[d]=whereZero(x[d]-1.)
991 else:
992 v_mask[d]=whereZero(x[d])
993 mod.setExternals(F=BF,fixed_v_mask=v_mask)
994
995 n=self.dom.getNormal()
996 N_t=0
997 errors=[]
998 while N_t < self.N_dt:
999 t_ref=t+dt
1000 v_ref, s_ref,p_ref=self.getReference(t_ref)
1001 mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref)
1002 # mod.update(dt, eta_iter_max=10, iter_max=50, verbose=self.VERBOSE, usePCG=True, max_correction_steps=30) new version
1003 mod.update(dt, iter_max=50, verbose=self.VERBOSE, usePCG=True)
1004 self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)
1005 t+=dt
1006 N_t+=1
1007
1008 def check(self,N_t,mod,t_ref,v_ref, s_ref,p_ref):
1009 p=mod.getPressure()
1010 p-=integrate(p)/vol(self.dom)
1011 error_p=Lsup(mod.getPressure()-p_ref)/Lsup(p_ref)
1012 error_s=Lsup(mod.getDeviatoricStress()-s_ref)/Lsup(s_ref)
1013 error_v=Lsup(mod.getVelocity()-v_ref)/Lsup(v_ref)
1014 error_t=abs(mod.getTime()-t_ref)/abs(t_ref)
1015 if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v
1016 self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
1017 self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )
1018 self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )
1019 self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )
1020 def tearDown(self):
1021 del self.dom
1022
1023 def test_D2_Fixed_MuNone_LateStart(self):
1024 self.dom = Rectangle(self.NE,self.NE,order=2)
1025 self.mu=None
1026 self.latestart=True
1027 self.runIt()
1028 def test_D2_Fixed_Mu_LateStart(self):
1029 self.dom = Rectangle(self.NE,self.NE,order=2)
1030 self.mu=555.
1031 self.latestart=True
1032 self.runIt()
1033 def test_D2_Fixed_MuNone(self):
1034 self.dom = Rectangle(self.NE,self.NE,order=2)
1035 self.mu=None
1036 self.latestart=False
1037 self.runIt()
1038 def test_D2_Fixed_Mu(self):
1039 self.dom = Rectangle(self.NE,self.NE,order=2)
1040 self.mu=555.
1041 self.latestart=False
1042 self.runIt()
1043 def test_D2_Free_MuNone_LateStart(self):
1044 self.dom = Rectangle(self.NE,self.NE,order=2)
1045 self.mu=None
1046 self.latestart=True
1047 self.runIt(free=0)
1048 def test_D2_Free_Mu_LateStart(self):
1049 self.dom = Rectangle(self.NE,self.NE,order=2)
1050 self.mu=555.
1051 self.latestart=True
1052 self.runIt(free=0)
1053 def test_D2_Free_MuNone(self):
1054 self.dom = Rectangle(self.NE,self.NE,order=2)
1055 self.mu=None
1056 self.latestart=False
1057 self.runIt(free=0)
1058 def test_D2_Free_Mu(self):
1059 self.dom = Rectangle(self.NE,self.NE,order=2)
1060 self.mu=555.
1061 self.latestart=False
1062 self.runIt(free=0)
1063
1064 def test_D3_Fixed_MuNone_LateStart(self):
1065 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1066 self.mu=None
1067 self.latestart=True
1068 self.runIt()
1069 def test_D3_Fixed_Mu_LateStart(self):
1070 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1071 self.mu=555.
1072 self.latestart=True
1073 self.runIt()
1074 def test_D3_Fixed_MuNone(self):
1075 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1076 self.mu=None
1077 self.latestart=False
1078 self.runIt()
1079 def test_D3_Fixed_Mu(self):
1080 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1081 self.mu=555.
1082 self.latestart=False
1083 self.runIt()
1084 def test_D3_Free_MuNone_LateStart(self):
1085 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1086 self.mu=None
1087 self.latestart=True
1088 self.runIt(free=0)
1089 def test_D3_Free_Mu_LateStart(self):
1090 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1091 self.mu=555.
1092 self.latestart=True
1093 self.runIt(free=0)
1094 def test_D3_Free_MuNone(self):
1095 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1096 self.mu=None
1097 self.latestart=False
1098 self.runIt(free=0)
1099 def test_D3_Free_Mu(self):
1100 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1101 self.mu=555.
1102 self.latestart=False
1103 self.runIt(free=0)
1104
1105
1106 class Test_FaultSystem(unittest.TestCase):
1107 EPS=1.e-8
1108 NE=10
1109 def test_Fault_MaxValue(self):
1110 dom=Rectangle(2*self.NE,2*self.NE)
1111 x=dom.getX()
1112 f=FaultSystem(dim=2)
1113 f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
1114 f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
1115
1116 u=x[0]*(1.-x[0])*(1-x[1])
1117 t, loc=f.getMaxValue(u)
1118 p=f.getParametrization(x,t)[0]
1119 m, l=loc(u), loc(p)
1120 self.failUnless( m == 0.25, "wrong max value")
1121 self.failUnless( t == 1, "wrong max tag")
1122 self.failUnless( l == 0., "wrong max location")
1123
1124 u=x[1]*(1.-x[1])*(1-x[0])*x[0]
1125 t, loc=f.getMaxValue(u)
1126 p=f.getParametrization(x,t)[0]
1127 m, l=loc(u), loc(p)
1128 self.failUnless( m == 0.0625, "wrong max value")
1129 self.failUnless( t == 2, "wrong max tag")
1130 self.failUnless( l == 0.5, "wrong max location")
1131
1132 u=x[0]*(1.-x[0])*x[1]
1133 t, loc=f.getMaxValue(u)
1134 p=f.getParametrization(x,t)[0]
1135 m, l=loc(u), loc(p)
1136 self.failUnless( m == 0.25, "wrong max value")
1137 self.failUnless( t == 2, "wrong max tag")
1138 self.failUnless( l == 1.0, "wrong max location")
1139
1140 u=x[1]*(1.-x[1])*x[0]
1141 t, loc=f.getMaxValue(u)
1142 p=f.getParametrization(x,t)[0]
1143 m, l=loc(u), loc(p)
1144 self.failUnless( m == 0.25, "wrong max value")
1145 self.failUnless( t == 2, "wrong max tag")
1146 self.failUnless( l == 0., "wrong max location")
1147
1148 u=x[1]*(1.-x[1])*(1.-x[0])
1149 t, loc=f.getMaxValue(u)
1150 p=f.getParametrization(x,t)[0]
1151 m, l=loc(u), loc(p)
1152 self.failUnless( m == 0.25, "wrong max value")
1153 self.failUnless( t == 1, "wrong max tag")
1154 self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong max location")
1155 def test_Fault_MinValue(self):
1156 dom=Rectangle(2*self.NE,2*self.NE)
1157 x=dom.getX()
1158 f=FaultSystem(dim=2)
1159 f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
1160 f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
1161
1162 u=-x[0]*(1.-x[0])*(1-x[1])
1163 t, loc=f.getMinValue(u)
1164 p=f.getParametrization(x,t)[0]
1165 m, l=loc(u), loc(p)
1166 self.failUnless( m == -0.25, "wrong min value")
1167 self.failUnless( t == 1, "wrong min tag")
1168 self.failUnless( l == 0., "wrong min location")
1169 u=-x[1]*(1.-x[1])*(1-x[0])*x[0]
1170 t, loc=f.getMinValue(u)
1171 p=f.getParametrization(x,t)[0]
1172 m, l=loc(u), loc(p)
1173 self.failUnless( m == -0.0625, "wrong min value")
1174 self.failUnless( t == 2, "wrong min tag")
1175 self.failUnless( l == 0.5, "wrong min location")
1176 u=-x[0]*(1.-x[0])*x[1]
1177 t, loc=f.getMinValue(u)
1178 p=f.getParametrization(x,t)[0]
1179 m, l=loc(u), loc(p)
1180 self.failUnless( m == -0.25, "wrong min value")
1181 self.failUnless( t == 2, "wrong min tag")
1182 self.failUnless( l == 1.0, "wrong min location")
1183 u=-x[1]*(1.-x[1])*x[0]
1184 t, loc=f.getMinValue(u)
1185 p=f.getParametrization(x,t)[0]
1186 m, l=loc(u), loc(p)
1187 self.failUnless( m == -0.25, "wrong min value")
1188 self.failUnless( t == 2, "wrong min tag")
1189 self.failUnless( l == 0., "wrong min location")
1190 u=-x[1]*(1.-x[1])*(1.-x[0])
1191 t, loc=f.getMinValue(u)
1192 p=f.getParametrization(x,t)[0]
1193 m, l=loc(u), loc(p)
1194 self.failUnless( m == -0.25, "wrong min value")
1195 self.failUnless( t == 1, "wrong min tag")
1196 self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong min location")
1197
1198
1199 def test_Fault2D(self):
1200 f=FaultSystem(dim=2)
1201 top1=[ [1.,0.], [1.,1.], [0.,1.] ]
1202 self.failUnlessRaises(ValueError,f.addFault,V0=[1.,0],strikes=[pi/2, pi/2],ls=[1.,1.],tag=1,dips=top1)
1203 f.addFault(V0=[1.,0],strikes=[pi/2, pi],ls=[1.,1.],tag=1)
1204 self.failUnless(f.getDim() == 2, "wrong dimension")
1205 self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
1206 self.failUnless( 2. == f.getTotalLength(1), "length wrong")
1207 self.failUnless( 0. == f.getMediumDepth(1), "depth wrong")
1208 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 range")
1209 self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 range")
1210 self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsets")
1211 segs=f.getTopPolyline(1)
1212 self.failUnless( len(segs) == 3, "wrong number of segments")
1213 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1214 self.failUnless( segs[0].size == 2, "seg 0 has wrong size.")
1215 self.failUnless( numpy.linalg.norm(segs[0]-[1.,0.]) < self.EPS, "wrong vertex. 0 ")
1216 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1217 self.failUnless( segs[1].size == 2, "seg 1 has wrong size.")
1218 self.failUnless( numpy.linalg.norm(segs[1]-[1.,1.]) < self.EPS, "wrong vertex. 1 ")
1219 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1220 self.failUnless( segs[2].size == 2, "seg 2 has wrong size.")
1221 self.failUnless( numpy.linalg.norm(segs[2]-[0.,1.]) < self.EPS, "wrong vertex. 2 ")
1222 c=f.getCenterOnSurface()
1223 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1224 self.failUnless( c.size == 2, "center size is wrong")
1225 self.failUnless( numpy.linalg.norm(c-[2./3.,2./3.]) < self.EPS, "center has wrong coordinates.")
1226 o=f.getOrientationOnSurface()/pi*180.
1227 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1228
1229 top2=[ [10.,0.], [0.,10.] ]
1230 f.addFault(V0=[10.,0],strikes=[3.*pi/4],ls=[14.142135623730951], tag=2, w0_offsets=[0,20], w1_max=20)
1231 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1232 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length")
1233 self.failUnless( 0. == f.getMediumDepth(2), "depth wrong")
1234 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range")
1235 self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range")
1236 self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets")
1237 segs=f.getTopPolyline(2)
1238 self.failUnless( len(segs) == 2, "wrong number of segments")
1239 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1240 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.]) < self.EPS, "wrong vertex. 0 ")
1241 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1242 self.failUnless( numpy.linalg.norm(segs[1]-[0.,10.]) < self.EPS, "wrong vertex. 1 ")
1243 c=f.getCenterOnSurface()
1244 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1245 self.failUnless( c.size == 2, "center size is wrong")
1246 self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.]) < self.EPS, "center has wrong coordinates.")
1247 o=f.getOrientationOnSurface()/pi*180.
1248 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1249
1250 s,d=f.getSideAndDistance([0.,0.], tag=1)
1251 self.failUnless( s<0, "wrong side.")
1252 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1253 s,d=f.getSideAndDistance([0.,2.], tag=1)
1254 self.failUnless( s>0, "wrong side.")
1255 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1256 s,d=f.getSideAndDistance([1.,2.], tag=1)
1257 self.failUnless( s>0, "wrong side.")
1258 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1259 s,d=f.getSideAndDistance([2.,1.], tag=1)
1260 self.failUnless( s>0, "wrong side.")
1261 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1262 s,d=f.getSideAndDistance([2.,0.], tag=1)
1263 self.failUnless( s>0, "wrong side.")
1264 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1265 s,d=f.getSideAndDistance([0.,-1.], tag=1)
1266 self.failUnless( s<0, "wrong side.")
1267 self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
1268 s,d=f.getSideAndDistance([-1.,0], tag=1)
1269 self.failUnless( s<0, "wrong side.")
1270 self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
1271
1272
1273 f.transform(rot=-pi/2., shift=[-1.,-1.])
1274 self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
1275 self.failUnless( 2. == f.getTotalLength(1), "length after transformation wrong")
1276 self.failUnless( 0. == f.getMediumDepth(1), "depth after transformation wrong")
1277 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 after transformation range")
1278 self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
1279 self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsetsafter transformation ")
1280 segs=f.getTopPolyline(1)
1281 self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
1282 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1283 self.failUnless( segs[0].size == 2, "seg 0 has wrong size after transformation.")
1284 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1285 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1286 self.failUnless( segs[1].size == 2, "seg 1 has wrong size after transformation.")
1287 self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1288 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1289 self.failUnless( segs[2].size == 2, "seg 2 has wrong size after transformation.")
1290 self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.]) < self.EPS, "wrong vertex after transformation. 2 ")
1291 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1292 self.failUnless( 0. == f.getMediumDepth(2), "depth wrong after transformation")
1293 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range after transformation")
1294 self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1295 self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1296 segs=f.getTopPolyline(2)
1297 self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
1298 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1299 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9]) < self.EPS, "wrong vertex. 0 after transformation")
1300 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1301 self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.]) < self.EPS, "wrong vertex. 1 after transformation")
1302
1303 c=f.getCenterOnSurface()
1304 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1305 self.failUnless( c.size == 2, "center size is wrong")
1306 self.failUnless( numpy.linalg.norm(c-[7./5.,-7./5.]) < self.EPS, "center has wrong coordinates.")
1307 o=f.getOrientationOnSurface()/pi*180.
1308 self.failUnless( abs(o-45.) < self.EPS, "wrong orientation.")
1309
1310 p=f.getParametrization([-1.,0.],1)
1311 self.failUnless(p[1]==1., "wrong value.")
1312 self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1313 p=f.getParametrization([-0.5,0.],1)
1314 self.failUnless(p[1]==1., "wrong value.")
1315 self.failUnless(abs(p[0]-0.5)<self.EPS* 0.5, "wrong value.")
1316 p=f.getParametrization([0.,0.],1)
1317 self.failUnless(p[1]==1., "wrong value.")
1318 self.failUnless(abs(p[0]-1.)<self.EPS, "wrong value.")
1319 p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-8)
1320 self.failUnless(p[1]==0., "wrong value.")
1321 p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-6)
1322 self.failUnless(p[1]==1., "wrong value.")
1323 self.failUnless(abs(p[0]-1.0000001)<self.EPS, "wrong value.")
1324 p=f.getParametrization([0.,0.5],1)
1325 self.failUnless(p[1]==1., "wrong value.")
1326 self.failUnless(abs(p[0]-1.5)<self.EPS, "wrong value.")
1327 p=f.getParametrization([0,1.],1)
1328 self.failUnless(p[1]==1., "wrong value.")
1329 self.failUnless(abs(p[0]-2.)<self.EPS, "wrong value.")
1330 p=f.getParametrization([1.,1.],1)
1331 self.failUnless(p[1]==0., "wrong value.")
1332 p=f.getParametrization([0,1.11],1)
1333 self.failUnless(p[1]==0., "wrong value.")
1334 p=f.getParametrization([-1,-9.],2)
1335 self.failUnless(p[1]==1., "wrong value.")
1336 self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1337 p=f.getParametrization([9,1],2)
1338 self.failUnless(p[1]==1., "wrong value.")
1339 self.failUnless(abs(p[0]-20.)<self.EPS, "wrong value.")
1340
1341 def test_Fault3D(self):
1342 f=FaultSystem(dim=3)
1343 self.failUnless(f.getDim() == 3, "wrong dimension")
1344
1345 top1=[ [0.,0.,0.], [1., 0., 0.] ]
1346 f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=20.,tag=1)
1347 self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
1348 self.failUnless( 1. == f.getTotalLength(1), "length wrong")
1349 self.failUnless( 20. == f.getMediumDepth(1), "depth wrong")
1350 self.failUnless( (0., 1.) == f.getW0Range(1)," wrong W0 range")
1351 self.failUnless( (-20., 0.) == f.getW1Range(1)," wrong W1 range")
1352 self.failUnless( [0., 1.] == f.getW0Offsets(1)," wrong W0 offsets")
1353 segs=f.getTopPolyline(1)
1354 self.failUnless( len(segs) == 2, "wrong number of segments")
1355 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1356 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1357 self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1358 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1359 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1360 self.failUnless( numpy.linalg.norm(segs[1]-[1.,0.,0]) < self.EPS, "wrong vertex. 1 ")
1361 c=f.getCenterOnSurface()
1362 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1363 self.failUnless( c.size == 3, "center size is wrong")
1364 self.failUnless( numpy.linalg.norm(c-[0.5,0.,0.]) < self.EPS, "center has wrong coordinates.")
1365 o=f.getOrientationOnSurface()/pi*180.
1366 self.failUnless( abs(o) < self.EPS, "wrong orientation.")
1367 d=f.getDips(1)
1368 self.failUnless( len(d) == 1, "wrong number of dips")
1369 self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1370 sn=f.getSegmentNormals(1)
1371 self.failUnless( len(sn) == 1, "wrong number of normals")
1372 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1373 self.failUnless( numpy.linalg.norm(sn[0]-[0, -1., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1374 dv=f.getDepthVectors(1)
1375 self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1376 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1377 self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1378 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1379 self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1380 b=f.getBottomPolyline(1)
1381 self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1382 self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1383 self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1384 self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1385 self.failUnless( numpy.linalg.norm(b[1]-[1., 0., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1386 ds=f.getDepths(1)
1387 self.failUnless( len(ds) == 2, "wrong number of depth")
1388 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1389 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1390
1391 top2=[ [0.,0.,0.], [0., 10., 0.] ]
1392 f.addFault(V0=[0.,0,0],strikes=[pi/2],ls=[10.,], dips=pi/2, depths=20.,tag=2)
1393 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1394 self.failUnless( 10. == f.getTotalLength(2), "length wrong")
1395 self.failUnless( 20. == f.getMediumDepth(2), "depth wrong")
1396 self.failUnless( (0., 10.) == f.getW0Range(2)," wrong W0 range")
1397 self.failUnless( (-20., 0.) == f.getW1Range(2)," wrong W1 range")
1398 self.failUnless( [0., 10.] == f.getW0Offsets(2)," wrong W0 offsets")
1399 segs=f.getTopPolyline(2)
1400 self.failUnless( len(segs) == 2, "wrong number of segments")
1401 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1402 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1403 self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1404 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1405 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1406 self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1407 d=f.getDips(2)
1408 self.failUnless( len(d) == 1, "wrong number of dips")
1409 self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1410 sn=f.getSegmentNormals(2)
1411 self.failUnless( len(sn) == 1, "wrong number of normals")
1412 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1413 self.failUnless( numpy.linalg.norm(sn[0]-[1, 0., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1414 dv=f.getDepthVectors(2)
1415 self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1416 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1417 self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1418 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1419 self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1420 b=f.getBottomPolyline(2)
1421 self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1422 self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1423 self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1424 self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1425 self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1426 ds=f.getDepths(2)
1427 self.failUnless( len(ds) == 2, "wrong number of depth")
1428 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1429 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1430
1431 top2=[ [10.,0.,0.], [0., 10., 0.] ]
1432 f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/2, depths=30.,tag=2)
1433 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1434 self.failUnless( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1435 self.failUnless( 30. == f.getMediumDepth(2), "depth wrong")
1436 self.failUnless( (-30., 0.) == f.getW1Range(2)," wrong W1 range")
1437 segs=f.getTopPolyline(2)
1438 self.failUnless( len(segs) == 2, "wrong number of segments")
1439 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1440 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1441 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1442 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1443 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1444 self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1445 d=f.getDips(2)
1446 self.failUnless( len(d) == 1, "wrong number of dips")
1447 self.failUnless( abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1448 sn=f.getSegmentNormals(2)
1449 self.failUnless( len(sn) == 1, "wrong number of normals")
1450 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1451 self.failUnless( numpy.linalg.norm(sn[0]-[0.70710678118654746, 0.70710678118654746, 0.]) < self.EPS, "wrong bottom vertex 1 ")
1452 dv=f.getDepthVectors(2)
1453 self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1454 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1455 self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -30.]) < self.EPS, "wrong depth vector 0 ")
1456 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1457 self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -30.]) < self.EPS, "wrong depth vector 1 ")
1458 b=f.getBottomPolyline(2)
1459 self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1460 self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1461 self.failUnless( numpy.linalg.norm(b[0]-[10., 0., -30.]) < self.EPS, "wrong bottom vertex 0 ")
1462 self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1463 self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -30.]) < self.EPS, "wrong bottom vertex 1 ")
1464 ds=f.getDepths(2)
1465 self.failUnless( len(ds) == 2, "wrong number of depth")
1466 self.failUnless( abs(ds[0]-30.) < self.EPS, "wrong depth at vertex 0 ")
1467 self.failUnless( abs(ds[1]-30.) < self.EPS, "wrong depth at vertex 1 ")
1468
1469 top2=[ [10.,0.,0.], [0., 10., 0.] ]
1470 f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/4, depths=50.,tag=2)
1471 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1472 self.failUnless( abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1473 self.failUnless( 50. == f.getMediumDepth(2), "depth wrong")
1474 self.failUnless( (-50., 0.) == f.getW1Range(2)," wrong W1 range")
1475 segs=f.getTopPolyline(2)
1476 self.failUnless( len(segs) == 2, "wrong number of segments")
1477 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1478 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1479 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1480 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1481 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1482 self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1483 d=f.getDips(2)
1484 self.failUnless( len(d) == 1, "wrong number of dips")
1485 self.failUnless( abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1486 sn=f.getSegmentNormals(2)
1487 self.failUnless( len(sn) == 1, "wrong number of normals")
1488 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1489 self.failUnless( numpy.linalg.norm(sn[0]-[0.5,0.5,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1490 dv=f.getDepthVectors(2)
1491 self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1492 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1493 self.failUnless( numpy.linalg.norm(dv[0]-[25., 25., -35.355339059327378]) < self.EPS, "wrong depth vector 0 ")
1494 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1495 self.failUnless( numpy.linalg.norm(dv[1]-[25.,25., -35.355339059327378]) < self.EPS, "wrong depth vector 1 ")
1496 b=f.getBottomPolyline(2)
1497 self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1498 self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1499 self.failUnless( numpy.linalg.norm(b[0]-[35., 25., -35.355339059327378]) < self.EPS, "wrong bottom vertex 0 ")
1500 self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1501 self.failUnless( numpy.linalg.norm(b[1]-[25, 35., -35.355339059327378]) < self.EPS, "wrong bottom vertex 1 ")
1502 ds=f.getDepths(2)
1503 self.failUnless( len(ds) == 2, "wrong number of depth")
1504 self.failUnless( abs(ds[0]-50.) < self.EPS, "wrong depth at vertex 0 ")
1505 self.failUnless( abs(ds[1]-50.) < self.EPS, "wrong depth at vertex 1 ")
1506
1507 top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
1508 f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=1, dips=pi/4, depths=20.)
1509 self.failUnless( 20. == f.getTotalLength(1), "length wrong")
1510 self.failUnless( 20. == f.getMediumDepth(1), "depth wrong")
1511 segs=f.getTopPolyline(1)
1512 self.failUnless( len(segs) == 3, "wrong number of segments")
1513 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1514 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1515 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1516 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1517 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1518 self.failUnless( numpy.linalg.norm(segs[1]-[10.,10.,0.]) < self.EPS, "wrong vertex. 1 ")
1519 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1520 self.failUnless( segs[2].size == 3, "seg 2 has wrong size.")
1521 self.failUnless( numpy.linalg.norm(segs[2]-[0.,10.,0.]) < self.EPS, "wrong vertex. 2 ")
1522 d=f.getDips(1)
1523 self.failUnless( len(d) == 2, "wrong number of dips")
1524 self.failUnless( abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1525 self.failUnless( abs(d[1]-0.78539816339744828) < self.EPS, "wrong dip 0")
1526 ds=f.getDepths(1)
1527 self.failUnless( len(ds) == 3, "wrong number of depth")
1528 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1529 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1530 sn=f.getSegmentNormals(1)
1531 self.failUnless( len(sn) == 2, "wrong number of normals")
1532 self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1533 self.failUnless( numpy.linalg.norm(sn[0]-[0.70710678118654746,0.,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1534 self.failUnless( isinstance(sn[1], numpy.ndarray), "wrong class of bottom vertex 0")
1535 self.failUnless( numpy.linalg.norm(sn[1]-[0.,0.70710678118654746,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1536 dv=f.getDepthVectors(1)
1537 self.failUnless( len(dv) == 3, "wrong number of depth vectors.")
1538 self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1539 self.failUnless( numpy.linalg.norm(dv[0]-[14.142135623730951, 0., -14.142135623730951]) < self.EPS, "wrong depth vector 0 ")
1540 self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1541 self.failUnless( numpy.linalg.norm(dv[1]-[11.547005383792515,11.547005383792515, -11.547005383792515]) < self.EPS, "wrong depth vector 2 ")
1542 self.failUnless( isinstance(dv[2], numpy.ndarray), "wrong class of depth vector 1")
1543 self.failUnless( numpy.linalg.norm(dv[2]-[0.,14.142135623730951, -14.142135623730951]) < self.EPS, "wrong depth vector 2 ")
1544 segs=f.getBottomPolyline(1)
1545 self.failUnless( len(segs) == 3, "wrong number of segments")
1546 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1547 self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1548 self.failUnless( numpy.linalg.norm(segs[0]-[24.142135623730951,0.,-14.142135623730951]) < self.EPS, "wrong vertex. 0 ")
1549 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1550 self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1551 self.failUnless( numpy.linalg.norm(segs[1]-[21.547005383792515,21.547005383792515, -11.547005383792515]) < self.EPS, "wrong vertex. 1 ")
1552 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1553 self.failUnless( segs[2].size == 3, "seg 2 has wrong size.")
1554 self.failUnless( numpy.linalg.norm(segs[2]-[0., 24.142135623730951, -14.142135623730951]) < self.EPS, "wrong vertex. 2 ")
1555 self.failUnless( abs(0.-f.getW0Range(1)[0]) <=self.EPS," wrong W0 range (0)")
1556 self.failUnless( abs(31.857329272664341-f.getW0Range(1)[1]) <=self.EPS," wrong W0 range (1)")
1557 self.failUnless( abs(-20. - f.getW1Range(1)[0]) <=self.EPS," wrong W1 range (0)")
1558 self.failUnless( abs(0. - f.getW1Range(1)[1]) <=self.EPS," wrong W1 range (1)")
1559 self.failUnless( abs(0.0-f.getW0Offsets(1)[0])<=self.EPS," wrong W0 offsets (0)")
1560 self.failUnless( abs(15.92866463633217-f.getW0Offsets(1)[1])<=self.EPS," wrong W0 offsets (1)")
1561 self.failUnless( abs(31.857329272664341-f.getW0Offsets(1)[2])<=self.EPS," wrong W0 offsets(2)")
1562 #
1563 # ============ fresh start ====================
1564 #
1565 f.addFault(V0=[1.,0,0.],strikes=[pi/2, pi],ls=[1.,1.], dips=pi/2,depths=20,tag=1)
1566 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)
1567 c=f.getCenterOnSurface()
1568 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1569 self.failUnless( c.size == 3, "center size is wrong")
1570 self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.,0.]) < self.EPS, "center has wrong coordinates.")
1571 o=f.getOrientationOnSurface()/pi*180.
1572 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1573
1574 f.transform(rot=-pi/2., shift=[-1.,-1.,0.])
1575 self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
1576 self.failUnless( 2. == f.getTotalLength(1), "length after transformation wrong")
1577 self.failUnless( 20. == f.getMediumDepth(1), "depth after transformation wrong")
1578 rw0=f.getW0Range(1)
1579 self.failUnless( len(rw0) ==2, "wo range has wrong length")
1580 self.failUnless( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
1581 self.failUnless( abs(rw0[1]-2.) < self.EPS,"W0 1 wrong.")
1582 self.failUnless( (-20., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
1583 dips=f.getDips(1)
1584 self.failUnless(len(dips) == 2, "wrong number of dips.")
1585 self.failUnless( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
1586 self.failUnless( abs(dips[1]-1.5707963267948966) <= self.EPS, "wrong dip")
1587 ds=f.getDepths(1)
1588 self.failUnless( len(ds) == 3, "wrong number of depth")
1589 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1590 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1591 self.failUnless( abs(ds[2]-20.) < self.EPS, "wrong depth at vertex 1 ")
1592 segs=f.getTopPolyline(1)
1593 self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
1594 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1595 self.failUnless( segs[0].size == 3, "seg 0 has wrong size after transformation.")
1596 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1597 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1598 self.failUnless( segs[1].size == 3, "seg 1 has wrong size after transformation.")
1599 self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1600 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1601 self.failUnless( segs[2].size == 3, "seg 2 has wrong size after transformation.")
1602 self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.,0.]) < self.EPS, "wrong vertex after transformation. 2 ")
1603 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1604 self.failUnless( 20. == f.getMediumDepth(2), "depth wrong after transformation")
1605 rw0=f.getW0Range(2)
1606 self.failUnless( len(rw0) ==2, "wo range has wrong length")
1607 self.failUnless( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
1608 self.failUnless( abs(rw0[1]-20.) < self.EPS,"W0 1 wrong.")
1609 self.failUnless( (-20., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1610 self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1611 dips=f.getDips(2)
1612 self.failUnless(len(dips) == 1, "wrong number of dips.")
1613 self.failUnless( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
1614 ds=f.getDepths(2)
1615 self.failUnless( len(ds) == 2, "wrong number of depth")
1616 self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1617 self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1618 segs=f.getTopPolyline(2)
1619 self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
1620 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1621 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1622 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1623 self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.,0.]) < self.EPS, "wrong vertex. 1 after transformation")
1624 #
1625 # ============ fresh start ====================
1626 #
1627 f=FaultSystem(dim=3)
1628
1629 top1=[ [0.,0.,0.], [1., 0., 0.] ]
1630 f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=1.,tag=1)
1631 top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
1632 f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=2, dips=pi/4, depths=20.)
1633
1634 p,m=f.getParametrization([0.3,0.,-0.5],1)
1635 self.failUnless(length(p-[0.3,-0.5]) <= self.EPS, "wrong value.")
1636 self.failUnless(m==1., "wrong value.")
1637
1638 p,m=f.getParametrization([0.5,0.,-0.5],1)
1639 self.failUnless(length(p-[0.5,-0.5]) <= self.EPS, "wrong value.")
1640 self.failUnless(m==1., "wrong value.")
1641
1642 p,m=f.getParametrization([0.25,0.,-0.5],1)
1643 self.failUnless(length(p-[0.25,-0.5]) <= self.EPS, "wrong value.")
1644 self.failUnless(m==1., "wrong value.")
1645
1646 p,m=f.getParametrization([0.5,0.,-0.25],1)
1647 self.failUnless(length(p-[0.5,-0.25]) <= self.EPS, "wrong value.")
1648 self.failUnless(m==1., "wrong value.")
1649
1650 p,m=f.getParametrization([0.001,0.,-0.001],1)
1651 self.failUnless(length(p-[0.001, -0.001]) <= self.EPS, "wrong value.")
1652 self.failUnless(m==1., "wrong value.")
1653
1654 p,m=f.getParametrization([0.001,0.,0.001],1)
1655 self.failUnless(m==0., "wrong value.")
1656
1657 p,m=f.getParametrization([0.999,0.,0.001],1)
1658 self.failUnless(m==0., "wrong value.")
1659
1660 p,m=f.getParametrization([1.001,0.,-0.001],1)
1661 self.failUnless(m==0., "wrong value.")
1662 p,m=f.getParametrization([1.001,0.,-0.1],1)
1663 self.failUnless(m==0., "wrong value.")
1664 p,m=f.getParametrization([1.001,0.,-0.000000001],1)
1665 self.failUnless(m==0., "wrong value.")
1666
1667 p,m=f.getParametrization([0.999,0.,-0.001],1)
1668 self.failUnless(length(p-[0.999, -0.001]) <= self.EPS, "wrong value.")
1669 self.failUnless(m==1., "wrong value.")
1670
1671 p,m=f.getParametrization([ 16.29252873 , 6.46410161 ,-6.29252873], 2, tol=1.e-7)
1672 self.failUnless(m==1., "wrong value.")
1673 self.failUnless(length(p-[4.7785993908996511, -10]) <= self.EPS*10., "wrong value.")
1674 p,m=f.getParametrization([15.77350269, 12.77350269, -5.77350269], 2, tol=1.e-7)
1675 self.failUnless(m==1., "wrong value.")
1676 self.failUnless(length(p-[11.150065245432518, -10]) <= self.EPS*10., "wrong value.")
1677
1678 p,m=f.getParametrization([ 3., 17.0710678, -7.0710678], 2, tol=1.e-7)
1679 self.failUnless(m==1., "wrong value.")
1680 self.failUnless(length(p-[27.078729881764687, -10]) <= self.EPS*10., "wrong value.")
1681 p,m=f.getParametrization([9.30940108, 16.55204176, -6.55204176], 2, tol=1.e-7)
1682 self.failUnless(m==1., "wrong value.")
1683 self.failUnless(length(p-[20.707264027231822, -10]) <= self.EPS*10., "wrong value.")
1684
1685 p,m=f.getParametrization([ 21.54700538, 21.54700538, -11.54700538], 2, tol=1.e-7)
1686 self.failUnless(m==1., "wrong value.")
1687 self.failUnless(length(p-[15.92866463633217, -20]) <= self.EPS*10., "wrong value.")
1688
1689 p,m=f.getParametrization([ 0.,0.,0.], 2, tol=1.e-7)
1690 self.failUnless(m==0., "wrong value.")
1691
1692 p,m=f.getParametrization([ 11.,11.,0.], 2, tol=1.e-7)
1693 self.failUnless(m==0., "wrong value.")
1694
1695
1696 s,d=f.getSideAndDistance([0.,-1.,0.], tag=1)
1697 self.failUnless( s>0, "wrong side.")
1698 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1699 s,d=f.getSideAndDistance([1.,-1.,0.], tag=1)
1700 self.failUnless( s>0, "wrong side.")
1701 self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
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
1709
1710 s,d=f.getSideAndDistance([0.,0.,0.], tag=2)
1711 self.failUnless( s<0, "wrong side.")
1712 self.failUnless( abs(d-10.)<self.EPS, "wrong distance.")
1713 s,d=f.getSideAndDistance([5.,5.,0.], tag=2)
1714 self.failUnless( s<0, "wrong side.")
1715 self.failUnless( abs(d-5.)<self.EPS, "wrong distance.")
1716
1717 s,d=f.getSideAndDistance([10.,10.,-1.], tag=2)
1718 self.failUnless( s<0, "wrong side.")
1719 self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
1720 s,d=f.getSideAndDistance([10.,10.,-2.], tag=2)
1721 self.failUnless( s<0, "wrong side.")
1722 self.failUnless( abs(d-2.*0.70710678118654757)<self.EPS, "wrong distance.")
1723 s,d=f.getSideAndDistance([10.,10.,-3.], tag=2)
1724 self.failUnless( s<0, "wrong side.")
1725 self.failUnless( abs(d-3.*0.70710678118654757)<self.EPS, "wrong distance.")
1726
1727 s,d=f.getSideAndDistance([5.,12.,0], tag=2)
1728 self.failUnless( s>0, "wrong side.")
1729 self.failUnless( abs(d-2*0.70710678118654757)<self.EPS, "wrong distance.")
1730 s,d=f.getSideAndDistance([5.,12.,-1], tag=2)
1731 self.failUnless( s>0, "wrong side.")
1732 self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
1733 s,d=f.getSideAndDistance([5.,12.,-2], tag=2)
1734 # s not checked as it is undefined.
1735 self.failUnless( abs(d)<self.EPS, "wrong distance.")
1736 s,d=f.getSideAndDistance([5.,12.,-3], 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.,-4], tag=2)
1740 self.failUnless( s<0, "wrong side.")
1741 self.failUnless( abs(d-2.*0.70710678118654757)<self.EPS, "wrong distance.")
1742
1743 if __name__ == '__main__':
1744 suite = unittest.TestSuite()
1745 suite.addTest(unittest.makeSuite(Test_FaultSystem))
1746 suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
1747 suite.addTest(unittest.makeSuite(Test_Darcy3D))
1748 suite.addTest(unittest.makeSuite(Test_Darcy2D))
1749 # suite.addTest(Test_Darcy2D("testVarioF_FixedBottom_smallK"))
1750 suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
1751 suite.addTest(Test_StokesProblemCartesian3D("test_PCG_P_large"))
1752 suite.addTest(unittest.makeSuite(Test_Mountains3D))
1753 suite.addTest(unittest.makeSuite(Test_Mountains2D))
1754 suite.addTest(unittest.makeSuite(Test_Rheologies))
1755 ## suite.addTest(Test_IncompressibleIsotropicFlowCartesian("test_D2_Fixed_MuNone"))
1756 suite.addTest(unittest.makeSuite(Test_IncompressibleIsotropicFlowCartesian))
1757 s=unittest.TextTestRunner(verbosity=2).run(suite)
1758 if not s.wasSuccessful(): sys.exit(1)
1759

  ViewVC Help
Powered by ViewVC 1.1.26