← Back to branch summary

~vcs-imports/escript-finley/trunk

« back to all changes in this revision

Viewing changes to dudley/test/python/run_models.py

  • Committer: jfenwick
  • Date: 2010-10-11 01:48:14 UTC
  • Revision ID: svn-v4:77569008-7704-0410-b7a0-a92fef0b09fd:trunk:3259
Merging dudley and scons updates from branches

Show diffs side-by-side

added added

removed removed

Lines of Context:
dudley/test/python/run_models.py
 
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.dudley import Rectangle, Brick
 
34
 
 
35
from math import pi
 
36
import numpy
 
37
import sys
 
38
import os
 
39
#====================================================================================================================
 
40
try:
 
41
     DUDLEY_WORKDIR=os.environ['DUDLEY_WORKDIR']
 
42
except KeyError:
 
43
     DUDLEY_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
class Test_Rheologies(unittest.TestCase):
 
685
     """
 
686
     this is the program used to generate the powerlaw tests:
 
687
 
 
688
     TAU_Y=100.
 
689
     N=10
 
690
     M=5
 
691
 
 
692
     def getE(tau):
 
693
         if tau<=TAU_Y:
 
694
           return 1./(0.5+20*sqrt(tau))
 
695
         else:
 
696
           raise ValueError,"out of range."
 
697
     tau=[ i* (TAU_Y/N) for i in range(N+1)] + [TAU_Y for i in range(M)]
 
698
     e=[ tau[i]/getE(tau[i]) for i in range(N+1)]
 
699
     e+= [ (i+1+j)* (max(e)/N) for j in range(M) ]
 
700
 
 
701
     print tau
 
702
     print e
 
703
     """
 
704
     TOL=1.e-8
 
705
     def test_PowerLaw_Init(self):
 
706
         pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
 
707
 
 
708
         self.failUnlessEqual(pl.getNumMaterials(),3,"num materials is wrong")
 
709
         self.failUnlessEqual(pl.validMaterialId(0),True,"material id 0 not found")
 
710
         self.failUnlessEqual(pl.validMaterialId(1),True,"material id 1 not found")
 
711
         self.failUnlessEqual(pl.validMaterialId(2),True,"material id 2 not found")
 
712
         self.failUnlessEqual(pl.validMaterialId(3),False,"material id 3 not found")
 
713
 
 
714
         self.failUnlessEqual(pl.getFriction(),None,"initial friction wrong.")
 
715
         self.failUnlessEqual(pl.getTauY(),None,"initial tau y wrong.")
 
716
         pl.setDruckerPragerLaw(tau_Y=10,friction=3)
 
717
         self.failUnlessEqual(pl.getFriction(),3,"friction wrong.")
 
718
         self.failUnlessEqual(pl.getTauY(),10,"tau y wrong.")
 
719
 
 
720
         self.failUnlessEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
 
721
         pl.setElasticShearModulus(1000)
 
722
         self.failUnlessEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
 
723
 
 
724
         e=pl.getEtaN()
 
725
         self.failUnlessEqual(len(e),3,"initial length of etaN is wrong.")
 
726
         self.failUnlessEqual(e,[None, None, None],"initial etaN are wrong.")
 
727
         self.failUnlessEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
 
728
         self.failUnlessEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
 
729
         self.failUnlessEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
 
730
         self.failUnlessRaises(ValueError, pl.getEtaN, 3)
 
731
 
 
732
         self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
 
733
         self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
 
734
         self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
 
735
 
 
736
         pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
 
737
         self.failUnlessEqual(pl.getPower(),[1,2,3],"powers are wrong.")
 
738
         self.failUnlessEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
 
739
         self.failUnlessEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
 
740
 
 
741
         pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
 
742
         self.failUnlessEqual(pl.getPower(),[4,2,3],"powers are wrong.")
 
743
         self.failUnlessEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
 
744
         self.failUnlessEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
 
745
 
 
746
         self.failUnlessRaises(ValueError,pl.getPower,-1)
 
747
         self.failUnlessEqual(pl.getPower(0),4,"power 0 is wrong.")
 
748
         self.failUnlessEqual(pl.getPower(1),2,"power 1 is wrong.")
 
749
         self.failUnlessEqual(pl.getPower(2),3,"power 2 is wrong.")
 
750
         self.failUnlessRaises(ValueError,pl.getPower,3)
 
751
 
 
752
         self.failUnlessRaises(ValueError,pl.getTauT,-1)
 
753
         self.failUnlessEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
 
754
         self.failUnlessEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
 
755
         self.failUnlessEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
 
756
         self.failUnlessRaises(ValueError,pl.getTauT,3)
 
757
 
 
758
         self.failUnlessRaises(ValueError,pl.getEtaN,-1)
 
759
         self.failUnlessEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
 
760
         self.failUnlessEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
 
761
         self.failUnlessEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
 
762
         self.failUnlessRaises(ValueError,pl.getEtaN,3)
 
763
 
 
764
     def checkResult(self,id,gamma_dot_, eta, tau_ref):
 
765
         self.failUnless(eta>=0,"eta needs to be positive (test %s)"%id)
 
766
         error=abs(gamma_dot_*eta-tau_ref)
 
767
         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))
 
768
        
 
769
     def test_PowerLaw_Linear(self):
 
770
         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]
 
771
         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]
 
772
         pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
 
773
         pl.setDruckerPragerLaw(tau_Y=100.)
 
774
         pl.setPowerLaw(eta_N=2.)
 
775
         pl.setEtaTolerance(self.TOL)
 
776
         for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
 
777
        
 
778
     def test_PowerLaw_QuadLarge(self):
 
779
         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]
 
780
         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]
 
781
         pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
 
782
         pl.setDruckerPragerLaw(tau_Y=100.)
 
783
         pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2])
 
784
         pl.setEtaTolerance(self.TOL)
 
785
         for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
 
786
 
 
787
     def test_PowerLaw_QuadSmall(self):
 
788
         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]
 
789
         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]
 
790
         pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
 
791
         pl.setDruckerPragerLaw(tau_Y=100.)
 
792
         pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2])
 
793
         pl.setEtaTolerance(self.TOL)
 
794
         for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
 
795
 
 
796
     def test_PowerLaw_CubeLarge(self):
 
797
         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]
 
798
         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]
 
799
         pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
 
800
         pl.setDruckerPragerLaw(tau_Y=100.)
 
801
         pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3])
 
802
         pl.setEtaTolerance(self.TOL)
 
803
         for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
 
804
 
 
805
     def test_PowerLaw_CubeSmall(self):
 
806
         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]
 
807
         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]
 
808
         pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
 
809
         pl.setDruckerPragerLaw(tau_Y=100.)
 
810
         pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3])
 
811
         pl.setEtaTolerance(self.TOL)
 
812
         for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
 
813
 
 
814
     def test_PowerLaw_QuadLarge_CubeLarge(self):
 
815
         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]
 
816
         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]
 
817
 
 
818
         pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
 
819
         pl.setDruckerPragerLaw(tau_Y=100.)
 
820
         pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
 
821
         pl.setEtaTolerance(self.TOL)
 
822
         for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
 
823
 
 
824
     def test_PowerLaw_QuadLarge_CubeSmall(self):
 
825
         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]
 
826
         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]
 
827
 
 
828
         pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
 
829
         pl.setDruckerPragerLaw(tau_Y=100.)
 
830
         pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
 
831
         pl.setEtaTolerance(self.TOL)
 
832
         for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
 
833
 
 
834
     def test_PowerLaw_QuadSmall_CubeLarge(self):
 
835
         taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
 
836
         gamma_dot_s=[0.0, 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]
 
837
         pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
 
838
         pl.setDruckerPragerLaw(tau_Y=100.)
 
839
         pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
 
840
         pl.setEtaTolerance(self.TOL)
 
841
         for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
 
842
 
 
843
     def test_PowerLaw_QuadSmall_CubeSmall(self):
 
844
         taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
 
845
         gamma_dot_s=[0.0, 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]
 
846
         pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
 
847
         pl.setDruckerPragerLaw(tau_Y=100.)
 
848
         pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
 
849
         pl.setEtaTolerance(self.TOL)
 
850
         for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
 
851
 
 
852
     def test_PowerLaw_withShear(self):
 
853
         taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
 
854
         gamma_dot_s=[0.0, 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]
 
855
         pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
 
856
         pl.setDruckerPragerLaw(tau_Y=100.)
 
857
         pl.setPowerLaw(eta_N=2.)
 
858
         pl.setElasticShearModulus(3.)
 
859
         dt=1./3.
 
860
         pl.setEtaTolerance(self.TOL)
 
861
         self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
 
862
         for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
 
863
 
 
864
 
 
865
class Test_FaultSystem(unittest.TestCase):
 
866
   EPS=1.e-8
 
867
   NE=10
 
868
   def test_Fault_MaxValue(self):
 
869
      dom=Rectangle(2*self.NE,2*self.NE)
 
870
      x=dom.getX()
 
871
      f=FaultSystem(dim=2)
 
872
      f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
 
873
      f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
 
874
 
 
875
      u=x[0]*(1.-x[0])*(1-x[1])
 
876
      t, loc=f.getMaxValue(u)
 
877
      p=f.getParametrization(x,t)[0]
 
878
      m, l=loc(u), loc(p)
 
879
      self.failUnless(  m == 0.25, "wrong max value")
 
880
      self.failUnless(  t == 1, "wrong max tag")
 
881
      self.failUnless(  l == 0., "wrong max location")
 
882
 
 
883
      u=x[1]*(1.-x[1])*(1-x[0])*x[0]
 
884
      t, loc=f.getMaxValue(u)
 
885
      p=f.getParametrization(x,t)[0]
 
886
      m, l=loc(u), loc(p)
 
887
      self.failUnless(  m == 0.0625, "wrong max value")
 
888
      self.failUnless(  t == 2, "wrong max tag")
 
889
      self.failUnless(  l == 0.5, "wrong max location")
 
890
 
 
891
      u=x[0]*(1.-x[0])*x[1]
 
892
      t, loc=f.getMaxValue(u)
 
893
      p=f.getParametrization(x,t)[0]
 
894
      m, l=loc(u), loc(p)
 
895
      self.failUnless(  m == 0.25, "wrong max value")
 
896
      self.failUnless(  t == 2, "wrong max tag")
 
897
      self.failUnless(  l == 1.0, "wrong max location")
 
898
 
 
899
      u=x[1]*(1.-x[1])*x[0]
 
900
      t, loc=f.getMaxValue(u)
 
901
      p=f.getParametrization(x,t)[0]
 
902
      m, l=loc(u), loc(p)
 
903
      self.failUnless(  m == 0.25, "wrong max value")
 
904
      self.failUnless(  t == 2, "wrong max tag")
 
905
      self.failUnless(  l == 0., "wrong max location")
 
906
 
 
907
      u=x[1]*(1.-x[1])*(1.-x[0])
 
908
      t, loc=f.getMaxValue(u)
 
909
      p=f.getParametrization(x,t)[0]
 
910
      m, l=loc(u), loc(p)
 
911
      self.failUnless(  m == 0.25, "wrong max value")
 
912
      self.failUnless(  t == 1, "wrong max tag")
 
913
      self.failUnless(  abs(l-0.70710678118654) <= self.EPS,  "wrong max location")
 
914
   def test_Fault_MinValue(self):
 
915
      dom=Rectangle(2*self.NE,2*self.NE)
 
916
      x=dom.getX()
 
917
      f=FaultSystem(dim=2)
 
918
      f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
 
919
      f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
 
920
 
 
921
      u=-x[0]*(1.-x[0])*(1-x[1])
 
922
      t, loc=f.getMinValue(u)
 
923
      p=f.getParametrization(x,t)[0]
 
924
      m, l=loc(u), loc(p)
 
925
      self.failUnless(  m == -0.25, "wrong min value")
 
926
      self.failUnless(  t == 1, "wrong min tag")
 
927
      self.failUnless(  l == 0., "wrong min location")
 
928
      u=-x[1]*(1.-x[1])*(1-x[0])*x[0]
 
929
      t, loc=f.getMinValue(u)
 
930
      p=f.getParametrization(x,t)[0]
 
931
      m, l=loc(u), loc(p)
 
932
      self.failUnless(  m == -0.0625, "wrong min value")
 
933
      self.failUnless(  t == 2, "wrong min tag")
 
934
      self.failUnless(  l == 0.5, "wrong min location")
 
935
      u=-x[0]*(1.-x[0])*x[1]
 
936
      t, loc=f.getMinValue(u)
 
937
      p=f.getParametrization(x,t)[0]
 
938
      m, l=loc(u), loc(p)
 
939
      self.failUnless(  m == -0.25, "wrong min value")
 
940
      self.failUnless(  t == 2, "wrong min tag")
 
941
      self.failUnless(  l == 1.0, "wrong min location")
 
942
      u=-x[1]*(1.-x[1])*x[0]
 
943
      t, loc=f.getMinValue(u)
 
944
      p=f.getParametrization(x,t)[0]
 
945
      m, l=loc(u), loc(p)
 
946
      self.failUnless(  m == -0.25, "wrong min value")
 
947
      self.failUnless(  t == 2, "wrong min tag")
 
948
      self.failUnless(  l == 0., "wrong min location")
 
949
      u=-x[1]*(1.-x[1])*(1.-x[0])
 
950
      t, loc=f.getMinValue(u)
 
951
      p=f.getParametrization(x,t)[0]
 
952
      m, l=loc(u), loc(p)
 
953
      self.failUnless(  m == -0.25, "wrong min value")
 
954
      self.failUnless(  t == 1, "wrong min tag")
 
955
      self.failUnless(  abs(l-0.70710678118654) <= self.EPS,  "wrong min location")
 
956
 
 
957
      
 
958
   def test_Fault2D(self):
 
959
      f=FaultSystem(dim=2)
 
960
      top1=[ [1.,0.], [1.,1.], [0.,1.] ]
 
961
      self.failUnlessRaises(ValueError,f.addFault,V0=[1.,0],strikes=[pi/2, pi/2],ls=[1.,1.],tag=1,dips=top1)
 
962
      f.addFault(V0=[1.,0],strikes=[pi/2, pi],ls=[1.,1.],tag=1)
 
963
      self.failUnless(f.getDim() == 2, "wrong dimension")
 
964
      self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
 
965
      self.failUnless(  2. == f.getTotalLength(1), "length wrong")
 
966
      self.failUnless(  0. == f.getMediumDepth(1), "depth wrong")
 
967
      self.failUnless( (0., 2.) ==  f.getW0Range(1)," wrong W0 range")
 
968
      self.failUnless( (0., 0.) ==  f.getW1Range(1)," wrong W1 range")
 
969
      self.failUnless( [0., 1., 2.] ==  f.getW0Offsets(1)," wrong W0 offsets")
 
970
      segs=f.getTopPolyline(1)
 
971
      self.failUnless( len(segs) == 3, "wrong number of segments")
 
972
      self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
 
973
      self.failUnless( segs[0].size == 2, "seg 0 has wrong size.")
 
974
      self.failUnless( numpy.linalg.norm(segs[0]-[1.,0.]) < self.EPS, "wrong vertex. 0 ")
 
975
      self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
 
976
      self.failUnless( segs[1].size == 2, "seg 1 has wrong size.")
 
977
      self.failUnless( numpy.linalg.norm(segs[1]-[1.,1.]) < self.EPS, "wrong vertex. 1 ")
 
978
      self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
 
979
      self.failUnless( segs[2].size == 2, "seg 2 has wrong size.")
 
980
      self.failUnless( numpy.linalg.norm(segs[2]-[0.,1.]) < self.EPS, "wrong vertex. 2 ")
 
981
      c=f.getCenterOnSurface()
 
982
      self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
 
983
      self.failUnless( c.size == 2, "center size is wrong")
 
984
      self.failUnless( numpy.linalg.norm(c-[2./3.,2./3.]) < self.EPS, "center has wrong coordinates.")
 
985
      o=f.getOrientationOnSurface()/pi*180.
 
986
      self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
 
987
 
 
988
      top2=[ [10.,0.], [0.,10.] ]
 
989
      f.addFault(V0=[10.,0],strikes=[3.*pi/4],ls=[14.142135623730951], tag=2, w0_offsets=[0,20], w1_max=20)
 
990
      self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
 
991
      self.failUnless(  abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length")
 
992
      self.failUnless(  0. == f.getMediumDepth(2), "depth wrong")
 
993
      self.failUnless( (0., 20.) ==  f.getW0Range(2)," wrong W0 range")
 
994
      self.failUnless( (0., 0.) ==  f.getW1Range(2)," wrong W1 range")
 
995
      self.failUnless( [0., 20.] ==  f.getW0Offsets(2)," wrong W0 offsets")
 
996
      segs=f.getTopPolyline(2)
 
997
      self.failUnless( len(segs) == 2, "wrong number of segments")
 
998
      self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
 
999
      self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.]) < self.EPS, "wrong vertex. 0 ")
 
1000
      self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
 
1001
      self.failUnless( numpy.linalg.norm(segs[1]-[0.,10.]) < self.EPS, "wrong vertex. 1 ")
 
1002
      c=f.getCenterOnSurface()
 
1003
      self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
 
1004
      self.failUnless( c.size == 2, "center size is wrong")
 
1005
      self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.]) < self.EPS, "center has wrong coordinates.")
 
1006
      o=f.getOrientationOnSurface()/pi*180.
 
1007
      self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
 
1008
 
 
1009
      s,d=f.getSideAndDistance([0.,0.], tag=1)
 
1010
      self.failUnless( s<0, "wrong side.")
 
1011
      self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
 
1012
      s,d=f.getSideAndDistance([0.,2.], tag=1)
 
1013
      self.failUnless( s>0, "wrong side.")
 
1014
      self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
 
1015
      s,d=f.getSideAndDistance([1.,2.], tag=1)
 
1016
      self.failUnless( s>0, "wrong side.")
 
1017
      self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
 
1018
      s,d=f.getSideAndDistance([2.,1.], tag=1)
 
1019
      self.failUnless( s>0, "wrong side.")
 
1020
      self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
 
1021
      s,d=f.getSideAndDistance([2.,0.], tag=1)
 
1022
      self.failUnless( s>0, "wrong side.")
 
1023
      self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
 
1024
      s,d=f.getSideAndDistance([0.,-1.], tag=1)
 
1025
      self.failUnless( s<0, "wrong side.")
 
1026
      self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
 
1027
      s,d=f.getSideAndDistance([-1.,0], tag=1)
 
1028
      self.failUnless( s<0, "wrong side.")
 
1029
      self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
 
1030
 
 
1031
 
 
1032
      f.transform(rot=-pi/2., shift=[-1.,-1.])
 
1033
      self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
 
1034
      self.failUnless(  2. == f.getTotalLength(1), "length after transformation wrong")
 
1035
      self.failUnless(  0. == f.getMediumDepth(1), "depth after transformation wrong")
 
1036
      self.failUnless( (0., 2.) ==  f.getW0Range(1)," wrong W0 after transformation range")
 
1037
      self.failUnless( (0., 0.) ==  f.getW1Range(1)," wrong W1 rangeafter transformation ")
 
1038
      self.failUnless( [0., 1., 2.] ==  f.getW0Offsets(1)," wrong W0 offsetsafter transformation ")
 
1039
      segs=f.getTopPolyline(1)
 
1040
      self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
 
1041
      self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
 
1042
      self.failUnless( segs[0].size == 2, "seg 0 has wrong size after transformation.")
 
1043
      self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.]) < self.EPS, "wrong vertex. 0  after transformation")
 
1044
      self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex  after transformation1")
 
1045
      self.failUnless( segs[1].size == 2, "seg 1 has wrong size after transformation.")
 
1046
      self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.]) < self.EPS, "wrong vertex.  after transformation1 ")
 
1047
      self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex  after transformation2")
 
1048
      self.failUnless( segs[2].size == 2, "seg 2 has wrong size after transformation.")
 
1049
      self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.]) < self.EPS, "wrong vertex after transformation. 2 ")
 
1050
      self.failUnless(  abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
 
1051
      self.failUnless(  0. == f.getMediumDepth(2), "depth wrong after transformation")
 
1052
      self.failUnless( (0., 20.) ==  f.getW0Range(2)," wrong W0 range after transformation")
 
1053
      self.failUnless( (0., 0.) ==  f.getW1Range(2)," wrong W1 range after transformation")
 
1054
      self.failUnless( [0., 20.] ==  f.getW0Offsets(2)," wrong W0 offsets after transformation")
 
1055
      segs=f.getTopPolyline(2)
 
1056
      self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
 
1057
      self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
 
1058
      self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9]) < self.EPS, "wrong vertex. 0  after transformation")
 
1059
      self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
 
1060
      self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.]) < self.EPS, "wrong vertex. 1  after transformation")
 
1061
 
 
1062
      c=f.getCenterOnSurface()
 
1063
      self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
 
1064
      self.failUnless( c.size == 2, "center size is wrong")
 
1065
      self.failUnless( numpy.linalg.norm(c-[7./5.,-7./5.]) < self.EPS, "center has wrong coordinates.")
 
1066
      o=f.getOrientationOnSurface()/pi*180.
 
1067
      self.failUnless( abs(o-45.) < self.EPS, "wrong orientation.")
 
1068
 
 
1069
      p=f.getParametrization([-1.,0.],1)
 
1070
      self.failUnless(p[1]==1., "wrong value.")
 
1071
      self.failUnless(abs(p[0])<self.EPS, "wrong value.")
 
1072
      p=f.getParametrization([-0.5,0.],1)
 
1073
      self.failUnless(p[1]==1., "wrong value.")
 
1074
      self.failUnless(abs(p[0]-0.5)<self.EPS* 0.5, "wrong value.")
 
1075
      p=f.getParametrization([0.,0.],1)
 
1076
      self.failUnless(p[1]==1., "wrong value.")
 
1077
      self.failUnless(abs(p[0]-1.)<self.EPS, "wrong value.")
 
1078
      p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-8)
 
1079
      self.failUnless(p[1]==0., "wrong value.")
 
1080
      p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-6)
 
1081
      self.failUnless(p[1]==1., "wrong value.")
 
1082
      self.failUnless(abs(p[0]-1.0000001)<self.EPS, "wrong value.")
 
1083
      p=f.getParametrization([0.,0.5],1)
 
1084
      self.failUnless(p[1]==1., "wrong value.")
 
1085
      self.failUnless(abs(p[0]-1.5)<self.EPS, "wrong value.")
 
1086
      p=f.getParametrization([0,1.],1)
 
1087
      self.failUnless(p[1]==1., "wrong value.")
 
1088
      self.failUnless(abs(p[0]-2.)<self.EPS, "wrong value.")
 
1089
      p=f.getParametrization([1.,1.],1)
 
1090
      self.failUnless(p[1]==0., "wrong value.")
 
1091
      p=f.getParametrization([0,1.11],1)
 
1092
      self.failUnless(p[1]==0., "wrong value.")
 
1093
      p=f.getParametrization([-1,-9.],2)
 
1094
      self.failUnless(p[1]==1., "wrong value.")
 
1095
      self.failUnless(abs(p[0])<self.EPS, "wrong value.")
 
1096
      p=f.getParametrization([9,1],2)
 
1097
      self.failUnless(p[1]==1., "wrong value.")
 
1098
      self.failUnless(abs(p[0]-20.)<self.EPS, "wrong value.")
 
1099
 
 
1100
   def test_Fault3D(self):
 
1101
      f=FaultSystem(dim=3)
 
1102
      self.failUnless(f.getDim() == 3, "wrong dimension")
 
1103
 
 
1104
      top1=[ [0.,0.,0.], [1., 0., 0.] ]
 
1105
      f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=20.,tag=1)
 
1106
      self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
 
1107
      self.failUnless(  1. == f.getTotalLength(1), "length wrong")
 
1108
      self.failUnless(  20. == f.getMediumDepth(1), "depth wrong")
 
1109
      self.failUnless( (0., 1.) ==  f.getW0Range(1)," wrong W0 range")
 
1110
      self.failUnless( (-20., 0.) ==  f.getW1Range(1)," wrong W1 range")
 
1111
      self.failUnless( [0., 1.] ==  f.getW0Offsets(1)," wrong W0 offsets")
 
1112
      segs=f.getTopPolyline(1)
 
1113
      self.failUnless( len(segs) == 2, "wrong number of segments")
 
1114
      self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
 
1115
      self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
 
1116
      self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
 
1117
      self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
 
1118
      self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
 
1119
      self.failUnless( numpy.linalg.norm(segs[1]-[1.,0.,0]) < self.EPS, "wrong vertex. 1 ")
 
1120
      c=f.getCenterOnSurface()
 
1121
      self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
 
1122
      self.failUnless( c.size == 3, "center size is wrong")
 
1123
      self.failUnless( numpy.linalg.norm(c-[0.5,0.,0.]) < self.EPS, "center has wrong coordinates.")
 
1124
      o=f.getOrientationOnSurface()/pi*180.
 
1125
      self.failUnless( abs(o) < self.EPS, "wrong orientation.")
 
1126
      d=f.getDips(1)
 
1127
      self.failUnless( len(d) == 1, "wrong number of dips")
 
1128
      self.failUnless(  abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
 
1129
      sn=f.getSegmentNormals(1)
 
1130
      self.failUnless( len(sn) == 1, "wrong number of normals")
 
1131
      self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
 
1132
      self.failUnless( numpy.linalg.norm(sn[0]-[0, -1., 0.]) < self.EPS, "wrong bottom vertex 1 ")
 
1133
      dv=f.getDepthVectors(1)
 
1134
      self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
 
1135
      self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
 
1136
      self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
 
1137
      self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
 
1138
      self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
 
1139
      b=f.getBottomPolyline(1)
 
1140
      self.failUnless( len(b) == 2, "wrong number of bottom vertices")
 
1141
      self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
 
1142
      self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
 
1143
      self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
 
1144
      self.failUnless( numpy.linalg.norm(b[1]-[1., 0., -20.]) < self.EPS, "wrong bottom vertex 1 ")
 
1145
      ds=f.getDepths(1)
 
1146
      self.failUnless( len(ds) == 2, "wrong number of depth")
 
1147
      self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
 
1148
      self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
 
1149
 
 
1150
      top2=[ [0.,0.,0.], [0., 10., 0.] ]
 
1151
      f.addFault(V0=[0.,0,0],strikes=[pi/2],ls=[10.,], dips=pi/2, depths=20.,tag=2)
 
1152
      self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
 
1153
      self.failUnless(  10. == f.getTotalLength(2), "length wrong")
 
1154
      self.failUnless(  20. == f.getMediumDepth(2), "depth wrong")
 
1155
      self.failUnless( (0., 10.) ==  f.getW0Range(2)," wrong W0 range")
 
1156
      self.failUnless( (-20., 0.) ==  f.getW1Range(2)," wrong W1 range")
 
1157
      self.failUnless( [0., 10.] ==  f.getW0Offsets(2)," wrong W0 offsets")
 
1158
      segs=f.getTopPolyline(2)
 
1159
      self.failUnless( len(segs) == 2, "wrong number of segments")
 
1160
      self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
 
1161
      self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
 
1162
      self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
 
1163
      self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
 
1164
      self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
 
1165
      self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
 
1166
      d=f.getDips(2)
 
1167
      self.failUnless( len(d) == 1, "wrong number of dips")
 
1168
      self.failUnless(  abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
 
1169
      sn=f.getSegmentNormals(2)
 
1170
      self.failUnless( len(sn) == 1, "wrong number of normals")
 
1171
      self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
 
1172
      self.failUnless( numpy.linalg.norm(sn[0]-[1, 0., 0.]) < self.EPS, "wrong bottom vertex 1 ")
 
1173
      dv=f.getDepthVectors(2)
 
1174
      self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
 
1175
      self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
 
1176
      self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
 
1177
      self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
 
1178
      self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
 
1179
      b=f.getBottomPolyline(2)
 
1180
      self.failUnless( len(b) == 2, "wrong number of bottom vertices")
 
1181
      self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
 
1182
      self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
 
1183
      self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
 
1184
      self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -20.]) < self.EPS, "wrong bottom vertex 1 ")
 
1185
      ds=f.getDepths(2)
 
1186
      self.failUnless( len(ds) == 2, "wrong number of depth")
 
1187
      self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
 
1188
      self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
 
1189
 
 
1190
      top2=[ [10.,0.,0.], [0., 10., 0.] ]
 
1191
      f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/2, depths=30.,tag=2)
 
1192
      self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
 
1193
      self.failUnless(  abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
 
1194
      self.failUnless(  30. == f.getMediumDepth(2), "depth wrong")
 
1195
      self.failUnless( (-30., 0.) ==  f.getW1Range(2)," wrong W1 range")
 
1196
      segs=f.getTopPolyline(2)
 
1197
      self.failUnless( len(segs) == 2, "wrong number of segments")
 
1198
      self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
 
1199
      self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
 
1200
      self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
 
1201
      self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
 
1202
      self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
 
1203
      self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
 
1204
      d=f.getDips(2)
 
1205
      self.failUnless( len(d) == 1, "wrong number of dips")
 
1206
      self.failUnless(  abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
 
1207
      sn=f.getSegmentNormals(2)
 
1208
      self.failUnless( len(sn) == 1, "wrong number of normals")
 
1209
      self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
 
1210
      self.failUnless( numpy.linalg.norm(sn[0]-[0.70710678118654746, 0.70710678118654746, 0.]) < self.EPS, "wrong bottom vertex 1 ")
 
1211
      dv=f.getDepthVectors(2)
 
1212
      self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
 
1213
      self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
 
1214
      self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -30.]) < self.EPS, "wrong depth vector 0 ")
 
1215
      self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
 
1216
      self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -30.]) < self.EPS, "wrong depth vector 1 ")
 
1217
      b=f.getBottomPolyline(2)
 
1218
      self.failUnless( len(b) == 2, "wrong number of bottom vertices")
 
1219
      self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
 
1220
      self.failUnless( numpy.linalg.norm(b[0]-[10., 0., -30.]) < self.EPS, "wrong bottom vertex 0 ")
 
1221
      self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
 
1222
      self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -30.]) < self.EPS, "wrong bottom vertex 1 ")
 
1223
      ds=f.getDepths(2)
 
1224
      self.failUnless( len(ds) == 2, "wrong number of depth")
 
1225
      self.failUnless( abs(ds[0]-30.) < self.EPS, "wrong depth at vertex 0 ")
 
1226
      self.failUnless( abs(ds[1]-30.) < self.EPS, "wrong depth at vertex 1 ")
 
1227
 
 
1228
      top2=[ [10.,0.,0.], [0., 10., 0.] ]
 
1229
      f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/4, depths=50.,tag=2)
 
1230
      self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
 
1231
      self.failUnless(  abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
 
1232
      self.failUnless(  50. == f.getMediumDepth(2), "depth wrong")
 
1233
      self.failUnless( (-50., 0.) ==  f.getW1Range(2)," wrong W1 range")
 
1234
      segs=f.getTopPolyline(2)
 
1235
      self.failUnless( len(segs) == 2, "wrong number of segments")
 
1236
      self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
 
1237
      self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
 
1238
      self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
 
1239
      self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
 
1240
      self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
 
1241
      self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
 
1242
      d=f.getDips(2)
 
1243
      self.failUnless( len(d) == 1, "wrong number of dips")
 
1244
      self.failUnless(  abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
 
1245
      sn=f.getSegmentNormals(2)
 
1246
      self.failUnless( len(sn) == 1, "wrong number of normals")
 
1247
      self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
 
1248
      self.failUnless( numpy.linalg.norm(sn[0]-[0.5,0.5,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
 
1249
      dv=f.getDepthVectors(2)
 
1250
      self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
 
1251
      self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
 
1252
      self.failUnless( numpy.linalg.norm(dv[0]-[25., 25., -35.355339059327378]) < self.EPS, "wrong depth vector 0 ")
 
1253
      self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
 
1254
      self.failUnless( numpy.linalg.norm(dv[1]-[25.,25., -35.355339059327378]) < self.EPS, "wrong depth vector 1 ")
 
1255
      b=f.getBottomPolyline(2)
 
1256
      self.failUnless( len(b) == 2, "wrong number of bottom vertices")
 
1257
      self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
 
1258
      self.failUnless( numpy.linalg.norm(b[0]-[35., 25., -35.355339059327378]) < self.EPS, "wrong bottom vertex 0 ")
 
1259
      self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
 
1260
      self.failUnless( numpy.linalg.norm(b[1]-[25, 35., -35.355339059327378]) < self.EPS, "wrong bottom vertex 1 ")
 
1261
      ds=f.getDepths(2)
 
1262
      self.failUnless( len(ds) == 2, "wrong number of depth")
 
1263
      self.failUnless( abs(ds[0]-50.) < self.EPS, "wrong depth at vertex 0 ")
 
1264
      self.failUnless( abs(ds[1]-50.) < self.EPS, "wrong depth at vertex 1 ")
 
1265
 
 
1266
      top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
 
1267
      f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=1, dips=pi/4, depths=20.)
 
1268
      self.failUnless(  20. == f.getTotalLength(1), "length wrong")
 
1269
      self.failUnless(  20. == f.getMediumDepth(1), "depth wrong")
 
1270
      segs=f.getTopPolyline(1)
 
1271
      self.failUnless( len(segs) == 3, "wrong number of segments")
 
1272
      self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
 
1273
      self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
 
1274
      self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
 
1275
      self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
 
1276
      self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
 
1277
      self.failUnless( numpy.linalg.norm(segs[1]-[10.,10.,0.]) < self.EPS, "wrong vertex. 1 ")
 
1278
      self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
 
1279
      self.failUnless( segs[2].size == 3, "seg 2 has wrong size.")
 
1280
      self.failUnless( numpy.linalg.norm(segs[2]-[0.,10.,0.]) < self.EPS, "wrong vertex. 2 ")
 
1281
      d=f.getDips(1)
 
1282
      self.failUnless( len(d) == 2, "wrong number of dips")
 
1283
      self.failUnless(  abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
 
1284
      self.failUnless(  abs(d[1]-0.78539816339744828) < self.EPS, "wrong dip 0")
 
1285
      ds=f.getDepths(1)
 
1286
      self.failUnless( len(ds) == 3, "wrong number of depth")
 
1287
      self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
 
1288
      self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
 
1289
      sn=f.getSegmentNormals(1)
 
1290
      self.failUnless( len(sn) == 2, "wrong number of normals")
 
1291
      self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
 
1292
      self.failUnless( numpy.linalg.norm(sn[0]-[0.70710678118654746,0.,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
 
1293
      self.failUnless( isinstance(sn[1], numpy.ndarray), "wrong class of bottom vertex 0")
 
1294
      self.failUnless( numpy.linalg.norm(sn[1]-[0.,0.70710678118654746,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
 
1295
      dv=f.getDepthVectors(1)
 
1296
      self.failUnless( len(dv) == 3, "wrong number of depth vectors.") 
 
1297
      self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
 
1298
      self.failUnless( numpy.linalg.norm(dv[0]-[14.142135623730951, 0., -14.142135623730951]) < self.EPS, "wrong depth vector 0 ")
 
1299
      self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
 
1300
      self.failUnless( numpy.linalg.norm(dv[1]-[11.547005383792515,11.547005383792515, -11.547005383792515]) < self.EPS, "wrong depth vector 2 ")
 
1301
      self.failUnless( isinstance(dv[2], numpy.ndarray), "wrong class of depth vector 1")
 
1302
      self.failUnless( numpy.linalg.norm(dv[2]-[0.,14.142135623730951, -14.142135623730951]) < self.EPS, "wrong depth vector 2 ")
 
1303
      segs=f.getBottomPolyline(1)
 
1304
      self.failUnless( len(segs) == 3, "wrong number of segments")
 
1305
      self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
 
1306
      self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
 
1307
      self.failUnless( numpy.linalg.norm(segs[0]-[24.142135623730951,0.,-14.142135623730951]) < self.EPS, "wrong vertex. 0 ")
 
1308
      self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
 
1309
      self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
 
1310
      self.failUnless( numpy.linalg.norm(segs[1]-[21.547005383792515,21.547005383792515, -11.547005383792515]) < self.EPS, "wrong vertex. 1 ")
 
1311
      self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
 
1312
      self.failUnless( segs[2].size == 3, "seg 2 has wrong size.")
 
1313
      self.failUnless( numpy.linalg.norm(segs[2]-[0., 24.142135623730951, -14.142135623730951]) < self.EPS, "wrong vertex. 2 ")
 
1314
      self.failUnless( abs(0.-f.getW0Range(1)[0]) <=self.EPS," wrong W0 range (0)")
 
1315
      self.failUnless( abs(31.857329272664341-f.getW0Range(1)[1]) <=self.EPS," wrong W0 range (1)")
 
1316
      self.failUnless( abs(-20. - f.getW1Range(1)[0]) <=self.EPS," wrong W1 range (0)")
 
1317
      self.failUnless( abs(0. - f.getW1Range(1)[1]) <=self.EPS," wrong W1 range (1)")
 
1318
      self.failUnless( abs(0.0-f.getW0Offsets(1)[0])<=self.EPS," wrong W0 offsets (0)")
 
1319
      self.failUnless( abs(15.92866463633217-f.getW0Offsets(1)[1])<=self.EPS," wrong W0 offsets (1)")
 
1320
      self.failUnless( abs(31.857329272664341-f.getW0Offsets(1)[2])<=self.EPS," wrong W0 offsets(2)")
 
1321
      #
 
1322
      #    ============ fresh start ====================
 
1323
      #
 
1324
      f.addFault(V0=[1.,0,0.],strikes=[pi/2, pi],ls=[1.,1.], dips=pi/2,depths=20,tag=1)
 
1325
      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)
 
1326
      c=f.getCenterOnSurface()
 
1327
      self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
 
1328
      self.failUnless( c.size == 3, "center size is wrong")
 
1329
      self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.,0.]) < self.EPS, "center has wrong coordinates.")
 
1330
      o=f.getOrientationOnSurface()/pi*180.
 
1331
      self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
 
1332
 
 
1333
      f.transform(rot=-pi/2., shift=[-1.,-1.,0.])
 
1334
      self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
 
1335
      self.failUnless(  2. == f.getTotalLength(1), "length after transformation wrong")
 
1336
      self.failUnless(  20. == f.getMediumDepth(1), "depth after transformation wrong")
 
1337
      rw0=f.getW0Range(1)
 
1338
      self.failUnless( len(rw0) ==2, "wo range has wrong length")
 
1339
      self.failUnless( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
 
1340
      self.failUnless( abs(rw0[1]-2.) < self.EPS,"W0 1 wrong.")
 
1341
      self.failUnless( (-20., 0.) ==  f.getW1Range(1)," wrong W1 rangeafter transformation ")
 
1342
      dips=f.getDips(1)
 
1343
      self.failUnless(len(dips) == 2, "wrong number of dips.")
 
1344
      self.failUnless( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
 
1345
      self.failUnless( abs(dips[1]-1.5707963267948966) <= self.EPS, "wrong dip")
 
1346
      ds=f.getDepths(1)
 
1347
      self.failUnless( len(ds) == 3, "wrong number of depth")
 
1348
      self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
 
1349
      self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
 
1350
      self.failUnless( abs(ds[2]-20.) < self.EPS, "wrong depth at vertex 1 ")
 
1351
      segs=f.getTopPolyline(1)
 
1352
      self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
 
1353
      self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
 
1354
      self.failUnless( segs[0].size == 3, "seg 0 has wrong size after transformation.")
 
1355
      self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.,0.]) < self.EPS, "wrong vertex. 0  after transformation")
 
1356
      self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex  after transformation1")
 
1357
      self.failUnless( segs[1].size == 3, "seg 1 has wrong size after transformation.")
 
1358
      self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.,0.]) < self.EPS, "wrong vertex.  after transformation1 ")
 
1359
      self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex  after transformation2")
 
1360
      self.failUnless( segs[2].size == 3, "seg 2 has wrong size after transformation.")
 
1361
      self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.,0.]) < self.EPS, "wrong vertex after transformation. 2 ")
 
1362
      self.failUnless(  abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
 
1363
      self.failUnless(  20. == f.getMediumDepth(2), "depth wrong after transformation")
 
1364
      rw0=f.getW0Range(2)
 
1365
      self.failUnless( len(rw0) ==2, "wo range has wrong length")
 
1366
      self.failUnless( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
 
1367
      self.failUnless( abs(rw0[1]-20.) < self.EPS,"W0 1 wrong.")
 
1368
      self.failUnless( (-20., 0.) ==  f.getW1Range(2)," wrong W1 range after transformation")
 
1369
      self.failUnless( [0., 20.] ==  f.getW0Offsets(2)," wrong W0 offsets after transformation")
 
1370
      dips=f.getDips(2)
 
1371
      self.failUnless(len(dips) == 1, "wrong number of dips.")
 
1372
      self.failUnless( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
 
1373
      ds=f.getDepths(2)
 
1374
      self.failUnless( len(ds) == 2, "wrong number of depth")
 
1375
      self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
 
1376
      self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
 
1377
      segs=f.getTopPolyline(2)
 
1378
      self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
 
1379
      self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
 
1380
      self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9,0.]) < self.EPS, "wrong vertex. 0  after transformation")
 
1381
      self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
 
1382
      self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.,0.]) < self.EPS, "wrong vertex. 1  after transformation")
 
1383
      #
 
1384
      #    ============ fresh start ====================
 
1385
      #
 
1386
      f=FaultSystem(dim=3)
 
1387
 
 
1388
      top1=[ [0.,0.,0.], [1., 0., 0.] ]
 
1389
      f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=1.,tag=1)
 
1390
      top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
 
1391
      f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=2, dips=pi/4, depths=20.)
 
1392
 
 
1393
      p,m=f.getParametrization([0.3,0.,-0.5],1)
 
1394
      self.failUnless(length(p-[0.3,-0.5]) <= self.EPS, "wrong value.")
 
1395
      self.failUnless(m==1., "wrong value.")
 
1396
 
 
1397
      p,m=f.getParametrization([0.5,0.,-0.5],1)
 
1398
      self.failUnless(length(p-[0.5,-0.5]) <= self.EPS, "wrong value.")
 
1399
      self.failUnless(m==1., "wrong value.")
 
1400
 
 
1401
      p,m=f.getParametrization([0.25,0.,-0.5],1)
 
1402
      self.failUnless(length(p-[0.25,-0.5]) <= self.EPS, "wrong value.")
 
1403
      self.failUnless(m==1., "wrong value.")
 
1404
 
 
1405
      p,m=f.getParametrization([0.5,0.,-0.25],1)
 
1406
      self.failUnless(length(p-[0.5,-0.25]) <= self.EPS, "wrong value.")
 
1407
      self.failUnless(m==1., "wrong value.")
 
1408
 
 
1409
      p,m=f.getParametrization([0.001,0.,-0.001],1)
 
1410
      self.failUnless(length(p-[0.001, -0.001]) <= self.EPS, "wrong value.")
 
1411
      self.failUnless(m==1., "wrong value.")
 
1412
 
 
1413
      p,m=f.getParametrization([0.001,0.,0.001],1)
 
1414
      self.failUnless(m==0., "wrong value.")
 
1415
 
 
1416
      p,m=f.getParametrization([0.999,0.,0.001],1)
 
1417
      self.failUnless(m==0., "wrong value.")
 
1418
 
 
1419
      p,m=f.getParametrization([1.001,0.,-0.001],1)
 
1420
      self.failUnless(m==0., "wrong value.")
 
1421
      p,m=f.getParametrization([1.001,0.,-0.1],1)
 
1422
      self.failUnless(m==0., "wrong value.")
 
1423
      p,m=f.getParametrization([1.001,0.,-0.000000001],1)
 
1424
      self.failUnless(m==0., "wrong value.")
 
1425
 
 
1426
      p,m=f.getParametrization([0.999,0.,-0.001],1)
 
1427
      self.failUnless(length(p-[0.999, -0.001]) <= self.EPS, "wrong value.")
 
1428
      self.failUnless(m==1., "wrong value.")
 
1429
 
 
1430
      p,m=f.getParametrization([ 16.29252873 , 6.46410161 ,-6.29252873], 2, tol=1.e-7)
 
1431
      self.failUnless(m==1., "wrong value.")
 
1432
      self.failUnless(length(p-[4.7785993908996511, -10]) <= self.EPS*10., "wrong value.")
 
1433
      p,m=f.getParametrization([15.77350269, 12.77350269, -5.77350269], 2, tol=1.e-7)
 
1434
      self.failUnless(m==1., "wrong value.")
 
1435
      self.failUnless(length(p-[11.150065245432518, -10]) <= self.EPS*10., "wrong value.")
 
1436
 
 
1437
      p,m=f.getParametrization([  3., 17.0710678, -7.0710678], 2, tol=1.e-7)
 
1438
      self.failUnless(m==1., "wrong value.")
 
1439
      self.failUnless(length(p-[27.078729881764687, -10]) <= self.EPS*10., "wrong value.")
 
1440
      p,m=f.getParametrization([9.30940108, 16.55204176, -6.55204176], 2, tol=1.e-7)
 
1441
      self.failUnless(m==1., "wrong value.")
 
1442
      self.failUnless(length(p-[20.707264027231822, -10]) <= self.EPS*10., "wrong value.")
 
1443
 
 
1444
      p,m=f.getParametrization([ 21.54700538,  21.54700538, -11.54700538], 2, tol=1.e-7)
 
1445
      self.failUnless(m==1., "wrong value.")
 
1446
      self.failUnless(length(p-[15.92866463633217, -20]) <= self.EPS*10., "wrong value.")
 
1447
 
 
1448
      p,m=f.getParametrization([ 0.,0.,0.], 2, tol=1.e-7)
 
1449
      self.failUnless(m==0., "wrong value.")
 
1450
 
 
1451
      p,m=f.getParametrization([ 11.,11.,0.], 2, tol=1.e-7)
 
1452
      self.failUnless(m==0., "wrong value.")
 
1453
 
 
1454
 
 
1455
      s,d=f.getSideAndDistance([0.,-1.,0.], tag=1)
 
1456
      self.failUnless( s>0, "wrong side.")
 
1457
      self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
 
1458
      s,d=f.getSideAndDistance([1.,-1.,0.], tag=1)
 
1459
      self.failUnless( s>0, "wrong side.")
 
1460
      self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
 
1461
      s,d=f.getSideAndDistance([0.,1.,0.], tag=1)
 
1462
      self.failUnless( s<0, "wrong side.")
 
1463
      self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
 
1464
      s,d=f.getSideAndDistance([1.,1.,0.], tag=1)
 
1465
      self.failUnless( s<0, "wrong side.")
 
1466
      self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
 
1467
 
 
1468
    
 
1469
      s,d=f.getSideAndDistance([0.,0.,0.], tag=2)
 
1470
      self.failUnless( s<0, "wrong side.")
 
1471
      self.failUnless( abs(d-10.)<self.EPS, "wrong distance.")
 
1472
      s,d=f.getSideAndDistance([5.,5.,0.], tag=2)
 
1473
      self.failUnless( s<0, "wrong side.")
 
1474
      self.failUnless( abs(d-5.)<self.EPS, "wrong distance.")
 
1475
 
 
1476
      s,d=f.getSideAndDistance([10.,10.,-1.], tag=2)
 
1477
      self.failUnless( s<0, "wrong side.")
 
1478
      self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
 
1479
      s,d=f.getSideAndDistance([10.,10.,-2.], tag=2)
 
1480
      self.failUnless( s<0, "wrong side.")
 
1481
      self.failUnless( abs(d-2.*0.70710678118654757)<self.EPS, "wrong distance.")
 
1482
      s,d=f.getSideAndDistance([10.,10.,-3.], tag=2)
 
1483
      self.failUnless( s<0, "wrong side.")
 
1484
      self.failUnless( abs(d-3.*0.70710678118654757)<self.EPS, "wrong distance.")
 
1485
 
 
1486
      s,d=f.getSideAndDistance([5.,12.,0], tag=2)
 
1487
      self.failUnless( s>0, "wrong side.")
 
1488
      self.failUnless( abs(d-2*0.70710678118654757)<self.EPS, "wrong distance.")
 
1489
      s,d=f.getSideAndDistance([5.,12.,-1], tag=2)
 
1490
      self.failUnless( s>0, "wrong side.")
 
1491
      self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
 
1492
      s,d=f.getSideAndDistance([5.,12.,-2], tag=2)
 
1493
      # s not checked as it is undefined.
 
1494
      self.failUnless( abs(d)<self.EPS, "wrong distance.")
 
1495
      s,d=f.getSideAndDistance([5.,12.,-3], tag=2)
 
1496
      self.failUnless( s<0, "wrong side.")
 
1497
      self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
 
1498
      s,d=f.getSideAndDistance([5.,12.,-4], tag=2)
 
1499
      self.failUnless( s<0, "wrong side.")
 
1500
      self.failUnless( abs(d-2.*0.70710678118654757)<self.EPS, "wrong distance.")
 
1501
 
 
1502
if __name__ == '__main__':
 
1503
   suite = unittest.TestSuite()
 
1504
   suite.addTest(unittest.makeSuite(Test_FaultSystem))
 
1505
#   suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
 
1506
   suite.addTest(unittest.makeSuite(Test_Darcy3D))
 
1507
   suite.addTest(unittest.makeSuite(Test_Darcy2D))
 
1508
#   suite.addTest(Test_StokesProblemCartesian2D("test_PCG_P_small"))
 
1509
#   suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
 
1510
   suite.addTest(unittest.makeSuite(Test_Rheologies))
 
1511
   s=unittest.TextTestRunner(verbosity=2).run(suite)
 
1512
   if not s.wasSuccessful(): sys.exit(1)
 
1513