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

  ViewVC Help
Powered by ViewVC 1.1.26