# Diff of /trunk/finley/test/python/run_models.py

revision 2208 by gross, Mon Jan 12 06:37:07 2009 UTC revision 3527 by gross, Tue May 24 06:59:20 2011 UTC
# Line 1  Line 1
1    # -*- coding: utf-8 -*-
2
3  ########################################################  ########################################################
4  #  #
5  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
6  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
7  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
8  #  #
# Line 11  Line 12
12  #  #
13  ########################################################  ########################################################
14
15  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2010 by University of Queensland
16  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
17  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
18  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
22
23  import unittest  import unittest
24  import tempfile  import tempfile
25
26
27
28    VERBOSE=False  and True
29
30  from esys.escript import *  from esys.escript import *
31  from esys.finley import Rectangle  from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian, FaultSystem, DarcyFlow
32  from esys.escript.models import DarcyFlow  from esys.escript.models import Mountains
33    from esys.finley import Rectangle, Brick
34
35    from math import pi
36    import numpy
37  import sys  import sys
38  import os  import os
39    #====================================================================================================================
40  try:  try:
41       FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']       FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
42  except KeyError:  except KeyError:
43       FINLEY_WORKDIR='.'       FINLEY_WORKDIR='.'
44

VERBOSE=False #

from esys.escript import *
from esys.escript.models import StokesProblemCartesian
from esys.finley import Rectangle, Brick

45  #====================================================================================================================  #====================================================================================================================
46  class Test_StokesProblemCartesian2D(unittest.TestCase):  class Test_StokesProblemCartesian2D(unittest.TestCase):
47     def setUp(self):     def setUp(self):
48         NE=6         NE=6
49         self.TOL=1e-3         self.TOL=1e-3
50         self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)         self.domain=Rectangle(NE,NE,order=-1)
51     def tearDown(self):     def tearDown(self):
52         del self.domain         del self.domain
53     def test_PCG_P_0(self):     def test_PCG_P_0(self):
# Line 62  class Test_StokesProblemCartesian2D(unit Line 65  class Test_StokesProblemCartesian2D(unit
65
67         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
68         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
69         sp.setTolerance(self.TOL)         sp.setTolerance(self.TOL)
70         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)         u,p=sp.solve(u0*mask,p0,verbose=VERBOSE,max_iter=100,usePCG=True)
71
72         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
73         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
# Line 88  class Test_StokesProblemCartesian2D(unit Line 91  class Test_StokesProblemCartesian2D(unit
91
93         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
94         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
95         sp.setTolerance(self.TOL*0.2)         sp.setTolerance(self.TOL)
96         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)         u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
97         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
98         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
99         error_p=Lsup(P1*x[0]*x[1]+p)         error_p=Lsup(P1*x[0]*x[1]+p)
saveVTK("d.vtu",p=p, e=P1*x[0]*x[1]+p, p_ref=P1*x[0]*x[1])
100         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")         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.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
102         self.failUnless(error_p<10*self.TOL, "pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
# Line 114  class Test_StokesProblemCartesian2D(unit Line 116  class Test_StokesProblemCartesian2D(unit
116
118         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
119         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
120         sp.setTolerance(self.TOL)         sp.setTolerance(self.TOL)
121         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)         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])         error_v0=Lsup(u[0]-u0[0])
125         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
# Line 140  class Test_StokesProblemCartesian2D(unit Line 143  class Test_StokesProblemCartesian2D(unit
143
145         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
146         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
147         sp.setTolerance(self.TOL)         sp.setTolerance(self.TOL)
148         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=50,useUzawa=False,iter_restart=18)         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])         error_v0=Lsup(u[0]-u0[0])
151         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
# Line 166  class Test_StokesProblemCartesian2D(unit Line 169  class Test_StokesProblemCartesian2D(unit
169
171         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
172         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
173         sp.setTolerance(self.TOL*0.1)         sp.setTolerance(self.TOL)
174         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=20,useUzawa=False)         u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=20,usePCG=False)
175
176         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
177         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
178         error_p=Lsup(P1*x[0]*x[1]+p)         error_p=Lsup(P1*x[0]*x[1]+p)
179
180         self.failUnless(error_p<10*self.TOL, "pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
181         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")         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.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
# Line 192  class Test_StokesProblemCartesian2D(unit Line 196  class Test_StokesProblemCartesian2D(unit
196
198         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
199         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
200         sp.setTolerance(self.TOL)         sp.setTolerance(self.TOL)
201         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=False)         u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)
202
203         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
204         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
# Line 207  class Test_StokesProblemCartesian3D(unit Line 211  class Test_StokesProblemCartesian3D(unit
211     def setUp(self):     def setUp(self):
212         NE=6         NE=6
213         self.TOL=1e-4         self.TOL=1e-4
214         self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)         self.domain=Brick(NE,NE,NE,order=-1)
215     def tearDown(self):     def tearDown(self):
216         del self.domain         del self.domain
217     def test_PCG_P_0(self):     def test_PCG_P_0(self):
# Line 229  class Test_StokesProblemCartesian3D(unit Line 233  class Test_StokesProblemCartesian3D(unit
233
235         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
236         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
237         sp.setTolerance(self.TOL)         sp.setTolerance(self.TOL)
238         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE ,max_iter=100,useUzawa=True)         u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
239
240         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
241         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
# Line 260  class Test_StokesProblemCartesian3D(unit Line 264  class Test_StokesProblemCartesian3D(unit
264
266         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
267         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
268         sp.setTolerance(self.TOL)         sp.setTolerance(self.TOL)
269         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE ,max_iter=100,useUzawa=True)         u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)

270         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
271         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
272         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
# Line 272  class Test_StokesProblemCartesian3D(unit Line 275  class Test_StokesProblemCartesian3D(unit
275         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")         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.")         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.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
278
279     def test_PCG_P_large(self):     def test_PCG_P_large(self):
280         ETA=1.         ETA=1.
281         P1=1000.         P1=1000.
# Line 290  class Test_StokesProblemCartesian3D(unit Line 294  class Test_StokesProblemCartesian3D(unit
294
296         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
297         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
298         sp.setTolerance(self.TOL)         sp.setTolerance(self.TOL)
299         u,p=sp.solve(u0,-p0,show_details=VERBOSE, verbose=VERBOSE ,max_iter=100,useUzawa=True)         u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
300
301         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
302         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
# Line 322  class Test_StokesProblemCartesian3D(unit Line 326  class Test_StokesProblemCartesian3D(unit
326
328         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
329         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
330         sp.setTolerance(self.TOL)         sp.setTolerance(self.TOL)
331         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=False,iter_restart=20)         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])         error_v0=Lsup(u[0]-u0[0])
334         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
335         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
336         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
# print error_p, error_v0,error_v1,error_v2
337         self.failUnless(error_p<10*self.TOL, "pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
338         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")         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.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
# Line 353  class Test_StokesProblemCartesian3D(unit Line 356  class Test_StokesProblemCartesian3D(unit
356
358         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
359         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
360         sp.setTolerance(self.TOL)         sp.setTolerance(self.TOL/10)
361         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE ,max_iter=100,useUzawa=False)         u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)
362
363         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
364         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
365         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
366         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
self.failUnless(error_p<10*self.TOL, "pressure error too large.")
367         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")         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.")         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.")         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):     def test_GMRES_P_large(self):
372         ETA=1.         ETA=1.
373         P1=1000.         P1=1000.
# Line 383  class Test_StokesProblemCartesian3D(unit Line 386  class Test_StokesProblemCartesian3D(unit
386
388         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
389         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
390         sp.setTolerance(self.TOL)         sp.setTolerance(self.TOL)
391         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE ,max_iter=100,useUzawa=False)         u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=False)
392
393         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
394         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
# Line 396  class Test_StokesProblemCartesian3D(unit Line 399  class Test_StokesProblemCartesian3D(unit
399         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")         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.")         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  class Test_Mountains3D(unittest.TestCase):
404      #     def setUp(self):
405      #         NE=16
406      #  p = 1/k * ( 1/2* (fb-f0)/xb* x **2 + f0 * x - ub*x ) +  p0         self.TOL=1e-4
407      #         self.domain=Brick(NE,NE,NE,order=1,useFullElementOrder=True)
408      #  with f = (fb-f0)/xb* x + f0     def tearDown(self):
409      #         del self.domain
410      #    u = f - k * p,x = ub     def test_periodic(self):
411      #         EPS=0.01
412      #  we prescribe pressure at x=x0=0 to p0
413      #         x=self.domain.getX()
414      #  if we prescribe the pressure on the bottom x=xb we set         v = Vector(0.0, Solution(self.domain))
415      #         a0=1
416      #  pb= 1/k * ( 1/2* (fb-f0)* xb + f0 * xb - ub*xb ) +  p0 = 1/k * ((fb+f0)/2 - ub ) * xb  + p0         a1=1
417      #         n0=2
418      #  which leads to ub = (fb+f0)/2-k*(pb-p0)/xb         n1=2
419      #         n2=0.5
420      def rescaleDomain(self):         a2=-(a0*n0+a1*n1)/n2
421          x=self.dom.getX().copy()         v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])* cos(pi*n2*x[2])
422          for i in xrange(self.dom.getDim()):         v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])* cos(pi*n2*x[2])
423               x_inf=inf(x[i])         v[2]=a2*cos(pi*n0*x[0])* cos(pi*n1*x[1])* sin(pi*n2*x[2])
424               x_sup=sup(x[i])
425               if i == self.dom.getDim()-1:         mts=Mountains(self.domain,eps=EPS)
426                  x[i]=-self.WIDTH*(x[i]-x_sup)/(x_inf-x_sup)         mts.setVelocity(v)
427               else:         Z=mts.update()
428                  x[i]=self.WIDTH*(x[i]-x_inf)/(x_sup-x_inf)
429          self.dom.setX(x)         error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[2]-1.)))
430      def getScalarMask(self,include_bottom=True):         self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
431          x=self.dom.getX().copy()
432          x_inf=inf(x[self.dom.getDim()-1])  class Test_Mountains2D(unittest.TestCase):
433          x_sup=sup(x[self.dom.getDim()-1])     def setUp(self):
434          out=whereZero(x[self.dom.getDim()-1]-x_sup)         NE=16
435          if include_bottom: out+=whereZero(x[self.dom.getDim()-1]-x_inf)         self.TOL=1e-4
436          return wherePositive(out)         self.domain=Rectangle(NE,NE,order=1,useFullElementOrder=True)
437      def getVectorMask(self,include_bottom=True):     def tearDown(self):
438          x=self.dom.getX().copy()         del self.domain
439          out=Vector(0.,Solution(self.dom))     def test_periodic(self):
440          for i in xrange(self.dom.getDim()):         EPS=0.01
441               x_inf=inf(x[i])
442               x_sup=sup(x[i])         x=self.domain.getX()
443               if i != self.dom.getDim()-1: out[i]+=whereZero(x[i]-x_sup)         v = Vector(0.0, Solution(self.domain))
444               if i != self.dom.getDim()-1 or include_bottom: out[i]+=whereZero(x[i]-x_inf)         a0=1
445          return wherePositive(out)         n0=1
446           n1=0.5
447      def setSolutionFixedBottom(self, p0, pb, f0, fb, k):         a1=-(a0*n0)/n1
448           d=self.dom.getDim()         v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])
449           x=self.dom.getX()[d-1]         v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])
450           xb=inf(x)
451           u=Vector(0.,Solution(self.dom))+kronecker(d)[d-1]*((f0+fb)/2.-k*(pb-p0)/xb)         H_t=Scalar(0.0, Solution(self.domain))
452           p=1./k*((fb-f0)/(xb*2.)* x**2 - (fb-f0)/2.*x)+(pb-p0)/xb*x +  p0         mts=Mountains(self.domain,eps=EPS)
453           f= ((fb-f0)/xb* x + f0)*kronecker(Function(self.dom))[d-1]         mts.setVelocity(v)
454           return u,p,f         Z=mts.update()
455
456           error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[1]-1.)))
457           self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
458
459
460
461    class Test_Rheologies(unittest.TestCase):
462         """
463         this is the program used to generate the powerlaw tests:
464
465         TAU_Y=100.
466         N=10
467         M=5
468
469         def getE(tau):
470             if tau<=TAU_Y:
471               return 1./(0.5+20*sqrt(tau))
472             else:
473               raise ValueError,"out of range."
474         tau=[ i* (TAU_Y/N) for i in range(N+1)] + [TAU_Y for i in range(M)]
475         e=[ tau[i]/getE(tau[i]) for i in range(N+1)]
476         e+= [ (i+1+j)* (max(e)/N) for j in range(M) ]
477
478         print tau
479         print e
480         """
481         TOL=1.e-8
482         def test_PowerLaw_Init(self):
483             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
484
485             self.failUnlessEqual(pl.getNumMaterials(),3,"num materials is wrong")
490
491             self.failUnlessEqual(pl.getFriction(),None,"initial friction wrong.")
492             self.failUnlessEqual(pl.getTauY(),None,"initial tau y wrong.")
493             pl.setDruckerPragerLaw(tau_Y=10,friction=3)
494             self.failUnlessEqual(pl.getFriction(),3,"friction wrong.")
495             self.failUnlessEqual(pl.getTauY(),10,"tau y wrong.")
496
497             self.failUnlessEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
498             pl.setElasticShearModulus(1000)
499             self.failUnlessEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
500
501             e=pl.getEtaN()
502             self.failUnlessEqual(len(e),3,"initial length of etaN is wrong.")
503             self.failUnlessEqual(e,[None, None, None],"initial etaN are wrong.")
504             self.failUnlessEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
505             self.failUnlessEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
506             self.failUnlessEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
507             self.failUnlessRaises(ValueError, pl.getEtaN, 3)
508
509             self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
510             self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
511             self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
512
513             pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
514             self.failUnlessEqual(pl.getPower(),[1,2,3],"powers are wrong.")
515             self.failUnlessEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
516             self.failUnlessEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
517
518             pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
519             self.failUnlessEqual(pl.getPower(),[4,2,3],"powers are wrong.")
520             self.failUnlessEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
521             self.failUnlessEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
522
523             self.failUnlessRaises(ValueError,pl.getPower,-1)
524             self.failUnlessEqual(pl.getPower(0),4,"power 0 is wrong.")
525             self.failUnlessEqual(pl.getPower(1),2,"power 1 is wrong.")
526             self.failUnlessEqual(pl.getPower(2),3,"power 2 is wrong.")
527             self.failUnlessRaises(ValueError,pl.getPower,3)
528
529             self.failUnlessRaises(ValueError,pl.getTauT,-1)
530             self.failUnlessEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
531             self.failUnlessEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
532             self.failUnlessEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
533             self.failUnlessRaises(ValueError,pl.getTauT,3)
534
535             self.failUnlessRaises(ValueError,pl.getEtaN,-1)
536             self.failUnlessEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
537             self.failUnlessEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
538             self.failUnlessEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
539             self.failUnlessRaises(ValueError,pl.getEtaN,3)
540
541         def checkResult(self,id,gamma_dot_, eta, tau_ref):
542             self.failUnless(eta>=0,"eta needs to be positive (test %s)"%id)
543             error=abs(gamma_dot_*eta-tau_ref)
544             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))
545
546      def testConstF_FixedBottom_smallK(self):       def test_PowerLaw_Linear(self):
547          k=1.e-10           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]
548          mp=self.getScalarMask(include_bottom=True)           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]
550          u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)           pl.setDruckerPragerLaw(tau_Y=100.)
551          p=p_ref*mp           pl.setPowerLaw(eta_N=2.)
552          u=u_ref*mv           pl.setEtaTolerance(self.TOL)
553          df=DarcyFlow(self.dom)           for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
554          df.setValue(g=f,
555                        location_of_fixed_pressure=mp,       def test_PowerLaw_QuadLarge(self):
556                        location_of_fixed_flux=mv,           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]
557                        permeability=Scalar(k,Function(self.dom)))           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]
558          df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)           pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
559          v,p=df.solve(u_ref,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)           pl.setDruckerPragerLaw(tau_Y=100.)
560          self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")           pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2])
561          self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")           pl.setEtaTolerance(self.TOL)
562             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
563      def testConstF_FixedBottom_mediumK(self):
564          k=1.       def test_PowerLaw_QuadSmall(self):
565          mp=self.getScalarMask(include_bottom=True)           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]
566          mv=self.getVectorMask(include_bottom=False)           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]
567          u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)           pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
568          p=p_ref*mp           pl.setDruckerPragerLaw(tau_Y=100.)
569          u=u_ref*mv           pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2])
570          df=DarcyFlow(self.dom)           pl.setEtaTolerance(self.TOL)
571          df.setValue(g=f,           for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
572                        location_of_fixed_pressure=mp,
573                        location_of_fixed_flux=mv,       def test_PowerLaw_CubeLarge(self):
574                        permeability=Scalar(k,Function(self.dom)))           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]
575          df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)           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]
576          v,p=df.solve(u,p,max_iter=100, verbose=VERBOSE ,sub_rtol=self.TOL/200)           pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
577          self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")           pl.setDruckerPragerLaw(tau_Y=100.)
578          self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")           pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3])
579             pl.setEtaTolerance(self.TOL)
580      def testConstF_FixedBottom_largeK(self):           for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
581          k=1.e10
582          mp=self.getScalarMask(include_bottom=True)       def test_PowerLaw_CubeSmall(self):
583          mv=self.getVectorMask(include_bottom=False)           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]
584          u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)           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]
585          p=p_ref*mp           pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
586          u=u_ref*mv           pl.setDruckerPragerLaw(tau_Y=100.)
587          df=DarcyFlow(self.dom)           pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3])
588          df.setValue(g=f,           pl.setEtaTolerance(self.TOL)
589                        location_of_fixed_pressure=mp,           for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
590                        location_of_fixed_flux=mv,
591                        permeability=Scalar(k,Function(self.dom)))       def test_PowerLaw_QuadLarge_CubeLarge(self):
592          df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)           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]
593          v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)           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]
594          self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
595          self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")           pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
596             pl.setDruckerPragerLaw(tau_Y=100.)
597      def testVarioF_FixedBottom_smallK(self):           pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
598          k=1.e-10           pl.setEtaTolerance(self.TOL)
599          mp=self.getScalarMask(include_bottom=True)           for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
601          u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)       def test_PowerLaw_QuadLarge_CubeSmall(self):
602          p=p_ref*mp           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]
603          u=u_ref*mv           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]
604          df=DarcyFlow(self.dom)
605          df.setValue(g=f,           pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
606                        location_of_fixed_pressure=mp,           pl.setDruckerPragerLaw(tau_Y=100.)
607                        location_of_fixed_flux=mv,           pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
608                        permeability=Scalar(k,Function(self.dom)))           pl.setEtaTolerance(self.TOL)
609          df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)           for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
610          v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
611          self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")       def test_PowerLaw_QuadSmall_CubeLarge(self):
612          self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")           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]
613             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]
614      def testVarioF_FixedBottom_mediumK(self):           pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
615          k=1.           pl.setDruckerPragerLaw(tau_Y=100.)
616          mp=self.getScalarMask(include_bottom=True)           pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
618          u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)           for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
619          p=p_ref*mp
620          u=u_ref*mv       def test_PowerLaw_QuadSmall_CubeSmall(self):
621          df=DarcyFlow(self.dom)           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]
622          df.setValue(g=f,           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]
623                        location_of_fixed_pressure=mp,           pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
624                        location_of_fixed_flux=mv,           pl.setDruckerPragerLaw(tau_Y=100.)
625                        permeability=Scalar(k,Function(self.dom)))           pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
626          df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)           pl.setEtaTolerance(self.TOL)
627          v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)           for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
628          self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
629          self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")       def test_PowerLaw_withShear(self):
630             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]
631      def testVarioF_FixedBottom_largeK(self):           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]
632          k=1.e10           pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
635          u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)           pl.setElasticShearModulus(3.)
636          p=p_ref*mp           dt=1./3.
637          u=u_ref*mv           pl.setEtaTolerance(self.TOL)
638          df=DarcyFlow(self.dom)           self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
639          df.setValue(g=f,           for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
640                        location_of_fixed_pressure=mp,
641                        location_of_fixed_flux=mv,  class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase):
642                        permeability=Scalar(k,Function(self.dom)))     TOL=1.e-5
643          df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)     VERBOSE=False # or True
644          v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)     A=1.
645          self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")     P_max=100
646          self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")     NE=2*getMPISizeWorld()
647       tau_Y=10.
648      def testConstF_FreeBottom_smallK(self):     N_dt=10
649          k=1.e-10
650          mp=self.getScalarMask(include_bottom=False)     # material parameter:
652          u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)     tau_2=5.
653          p=p_ref*mp     eta_0=100.
654          u=u_ref*mv     eta_1=50.
655          df=DarcyFlow(self.dom)     eta_2=400.
656          df.setValue(g=f,     N_1=2.
657                      location_of_fixed_pressure=mp,     N_2=3.
658                        location_of_fixed_flux=mv,     def getReference(self, t):
659                        permeability=Scalar(k,Function(self.dom)))
660          df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)        B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
661          v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)        x=self.dom.getX()
662          self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
663          self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")        s_00=min(self.A*t,B)
664          tau=sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)*abs(s_00)
665      def testConstF_FreeBottom_mediumK(self):        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.)
666          k=1.
668          mv=self.getVectorMask(include_bottom=True)        if s_00 <= B and self.mu !=None: alpha+=1./(2*self.mu)*self.A
669          u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)        u_ref=x*alpha
670          p=p_ref*mp        u_ref[self.dom.getDim()-1]=(1.-x[self.dom.getDim()-1])*alpha*(self.dom.getDim()-1)
671          u=u_ref*mv        sigma_ref=kronecker(self.dom)*s_00
672          df=DarcyFlow(self.dom)        sigma_ref[self.dom.getDim()-1,self.dom.getDim()-1]=-s_00*(self.dom.getDim()-1)
673          df.setValue(g=f,
674                        location_of_fixed_pressure=mp,        p_ref=self.P_max
675                        location_of_fixed_flux=mv,        for d in range(self.dom.getDim()): p_ref=p_ref*x[d]
676                        permeability=Scalar(k,Function(self.dom)))        p_ref-=integrate(p_ref)/vol(self.dom)
677          df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)        return u_ref, sigma_ref, p_ref
678          v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
679          self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")     def runIt(self, free=None):
680          self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")        x=self.dom.getX()
681          B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
682      def testConstF_FreeBottom_largeK(self):        dt=B/int(self.N_dt/2)
683          k=1.e10        if self.VERBOSE: print "dt =",dt
684          mp=self.getScalarMask(include_bottom=False)        if self.latestart:
686          u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)        else:
687          p=p_ref*mp            t=0
688          u=u_ref*mv        v,s,p=self.getReference(t)
689          df=DarcyFlow(self.dom)
690          df.setValue(g=f,        mod=IncompressibleIsotropicFlowCartesian(self.dom, stress=s, v=v, p=p, t=t, numMaterials=3, verbose=self.VERBOSE)
691                        location_of_fixed_pressure=mp,        mod.setDruckerPragerLaw(tau_Y=self.tau_Y,friction=None)
692                        location_of_fixed_flux=mv,        mod.setElasticShearModulus(self.mu)
693                        permeability=Scalar(k,Function(self.dom)))        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])
694          df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)        mod.setTolerance(self.TOL)
695          v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)        mod.setEtaTolerance(self.TOL*0.1)
696          self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
697          self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")        BF=Vector(self.P_max,Function(self.dom))
698          for d in range(self.dom.getDim()):
699      def testVarioF_FreeBottom_smallK(self):            for d2 in range(self.dom.getDim()):
700          k=1.e-10                if d!=d2: BF[d]*=x[d2]
702          mv=self.getVectorMask(include_bottom=True)        if free==None:
703          u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)           for d in range(self.dom.getDim()):
705          u=u_ref*mv        else:
706          df=DarcyFlow(self.dom)           for d in range(self.dom.getDim()):
707          df.setValue(g=f,              if d == self.dom.getDim()-1:
709                        location_of_fixed_flux=mv,              else:
712          v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
713          self.failUnless(Lsup(v-u_ref)<self.TOL*25.*Lsup(u_ref), "flux error too big.")  # 25 because of disc. error.        n=self.dom.getNormal()
714          self.failUnless(Lsup(p-p_ref)<self.TOL*25.*Lsup(p_ref), "pressure error too big.")        N_t=0
715          errors=[]
716      def testVarioF_FreeBottom_mediumK(self):        while N_t < self.N_dt:
717          k=1.           t_ref=t+dt
718          mp=self.getScalarMask(include_bottom=False)           v_ref, s_ref,p_ref=self.getReference(t_ref)
719          mv=self.getVectorMask(include_bottom=True)           mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref)
720          u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)           # mod.update(dt, eta_iter_max=10, iter_max=50, verbose=self.VERBOSE, usePCG=True, max_correction_steps=30) new version
721          p=p_ref*mp           mod.update(dt, iter_max=50, verbose=self.VERBOSE, usePCG=True)
722          u=u_ref*mv           self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)
723          df=DarcyFlow(self.dom)           t+=dt
724          df.setValue(g=f,           N_t+=1
725                        location_of_fixed_pressure=mp,
726                        location_of_fixed_flux=mv,     def check(self,N_t,mod,t_ref,v_ref, s_ref,p_ref):
727                        permeability=Scalar(k,Function(self.dom)))           p=mod.getPressure()
728          df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)           p-=integrate(p)/vol(self.dom)
729          v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)           error_p=Lsup(mod.getPressure()-p_ref)/Lsup(p_ref)
730          self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")           error_s=Lsup(mod.getDeviatoricStress()-s_ref)/Lsup(s_ref)
731          self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")           error_v=Lsup(mod.getVelocity()-v_ref)/Lsup(v_ref)
732             error_t=abs(mod.getTime()-t_ref)/abs(t_ref)
733      def testVarioF_FreeBottom_largeK(self):           if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v
734          k=1.e10           self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
735          mp=self.getScalarMask(include_bottom=False)           self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )
736          mv=self.getVectorMask(include_bottom=True)           self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )
737          u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)           self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )
738          p=p_ref*mp     def tearDown(self):
u=u_ref*mv
df=DarcyFlow(self.dom)
df.setValue(g=f,
location_of_fixed_pressure=mp,
location_of_fixed_flux=mv,
permeability=Scalar(k,Function(self.dom)))
df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)
v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")

class Test_Darcy2D(Test_Darcy):
TOL=1e-4
WIDTH=1.
def setUp(self):
NE=40  # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
self.dom = Rectangle(NE,NE)
self.rescaleDomain()
def tearDown(self):
del self.dom
class Test_Darcy3D(Test_Darcy):
TOL=1e-4
WIDTH=1.
def setUp(self):
NE=25  # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
self.dom = Brick(NE,NE,NE)
self.rescaleDomain()
def tearDown(self):
739          del self.dom          del self.dom
740
741       def test_D2_Fixed_MuNone_LateStart(self):
742           self.dom = Rectangle(self.NE,self.NE,order=2)
743           self.mu=None
744           self.latestart=True
745           self.runIt()
746       def test_D2_Fixed_Mu_LateStart(self):
747           self.dom = Rectangle(self.NE,self.NE,order=2)
748           self.mu=555.
749           self.latestart=True
750           self.runIt()
751       def test_D2_Fixed_MuNone(self):
752           self.dom = Rectangle(self.NE,self.NE,order=2)
753           self.mu=None
754           self.latestart=False
755           self.runIt()
756       def test_D2_Fixed_Mu(self):
757           self.dom = Rectangle(self.NE,self.NE,order=2)
758           self.mu=555.
759           self.latestart=False
760           self.runIt()
761       def test_D2_Free_MuNone_LateStart(self):
762           self.dom = Rectangle(self.NE,self.NE,order=2)
763           self.mu=None
764           self.latestart=True
765           self.runIt(free=0)
766       def test_D2_Free_Mu_LateStart(self):
767           self.dom = Rectangle(self.NE,self.NE,order=2)
768           self.mu=555.
769           self.latestart=True
770           self.runIt(free=0)
771       def test_D2_Free_MuNone(self):
772           self.dom = Rectangle(self.NE,self.NE,order=2)
773           self.mu=None
774           self.latestart=False
775           self.runIt(free=0)
776       def test_D2_Free_Mu(self):
777           self.dom = Rectangle(self.NE,self.NE,order=2)
778           self.mu=555.
779           self.latestart=False
780           self.runIt(free=0)
781
782       def test_D3_Fixed_MuNone_LateStart(self):
783           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
784           self.mu=None
785           self.latestart=True
786           self.runIt()
787       def test_D3_Fixed_Mu_LateStart(self):
788           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
789           self.mu=555.
790           self.latestart=True
791           self.runIt()
792       def test_D3_Fixed_MuNone(self):
793           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
794           self.mu=None
795           self.latestart=False
796           self.runIt()
797       def test_D3_Fixed_Mu(self):
798           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
799           self.mu=555.
800           self.latestart=False
801           self.runIt()
802       def test_D3_Free_MuNone_LateStart(self):
803           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
804           self.mu=None
805           self.latestart=True
806           self.runIt(free=0)
807       def test_D3_Free_Mu_LateStart(self):
808           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
809           self.mu=555.
810           self.latestart=True
811           self.runIt(free=0)
812       def test_D3_Free_MuNone(self):
813           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
814           self.mu=None
815           self.latestart=False
816           self.runIt(free=0)
817       def test_D3_Free_Mu(self):
818           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
819           self.mu=555.
820           self.latestart=False
821           self.runIt(free=0)
822
823
824    class Test_FaultSystem(unittest.TestCase):
825       EPS=1.e-8
826       NE=10
827       def test_Fault_MaxValue(self):
828          dom=Rectangle(2*self.NE,2*self.NE)
829          x=dom.getX()
830          f=FaultSystem(dim=2)
831          f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
832          f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
833
834          u=x[0]*(1.-x[0])*(1-x[1])
835          t, loc=f.getMaxValue(u)
836          p=f.getParametrization(x,t)[0]
837          m, l=loc(u), loc(p)
838          self.failUnless(  m == 0.25, "wrong max value")
839          self.failUnless(  t == 1, "wrong max tag")
840          self.failUnless(  l == 0., "wrong max location")
841
842          u=x[1]*(1.-x[1])*(1-x[0])*x[0]
843          t, loc=f.getMaxValue(u)
844          p=f.getParametrization(x,t)[0]
845          m, l=loc(u), loc(p)
846          self.failUnless(  m == 0.0625, "wrong max value")
847          self.failUnless(  t == 2, "wrong max tag")
848          self.failUnless(  l == 0.5, "wrong max location")
849
850          u=x[0]*(1.-x[0])*x[1]
851          t, loc=f.getMaxValue(u)
852          p=f.getParametrization(x,t)[0]
853          m, l=loc(u), loc(p)
854          self.failUnless(  m == 0.25, "wrong max value")
855          self.failUnless(  t == 2, "wrong max tag")
856          self.failUnless(  l == 1.0, "wrong max location")
857
858          u=x[1]*(1.-x[1])*x[0]
859          t, loc=f.getMaxValue(u)
860          p=f.getParametrization(x,t)[0]
861          m, l=loc(u), loc(p)
862          self.failUnless(  m == 0.25, "wrong max value")
863          self.failUnless(  t == 2, "wrong max tag")
864          self.failUnless(  l == 0., "wrong max location")
865
866          u=x[1]*(1.-x[1])*(1.-x[0])
867          t, loc=f.getMaxValue(u)
868          p=f.getParametrization(x,t)[0]
869          m, l=loc(u), loc(p)
870          self.failUnless(  m == 0.25, "wrong max value")
871          self.failUnless(  t == 1, "wrong max tag")
872          self.failUnless(  abs(l-0.70710678118654) <= self.EPS,  "wrong max location")
873       def test_Fault_MinValue(self):
874          dom=Rectangle(2*self.NE,2*self.NE)
875          x=dom.getX()
876          f=FaultSystem(dim=2)
877          f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
878          f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
879
880          u=-x[0]*(1.-x[0])*(1-x[1])
881          t, loc=f.getMinValue(u)
882          p=f.getParametrization(x,t)[0]
883          m, l=loc(u), loc(p)
884          self.failUnless(  m == -0.25, "wrong min value")
885          self.failUnless(  t == 1, "wrong min tag")
886          self.failUnless(  l == 0., "wrong min location")
887          u=-x[1]*(1.-x[1])*(1-x[0])*x[0]
888          t, loc=f.getMinValue(u)
889          p=f.getParametrization(x,t)[0]
890          m, l=loc(u), loc(p)
891          self.failUnless(  m == -0.0625, "wrong min value")
892          self.failUnless(  t == 2, "wrong min tag")
893          self.failUnless(  l == 0.5, "wrong min location")
894          u=-x[0]*(1.-x[0])*x[1]
895          t, loc=f.getMinValue(u)
896          p=f.getParametrization(x,t)[0]
897          m, l=loc(u), loc(p)
898          self.failUnless(  m == -0.25, "wrong min value")
899          self.failUnless(  t == 2, "wrong min tag")
900          self.failUnless(  l == 1.0, "wrong min location")
901          u=-x[1]*(1.-x[1])*x[0]
902          t, loc=f.getMinValue(u)
903          p=f.getParametrization(x,t)[0]
904          m, l=loc(u), loc(p)
905          self.failUnless(  m == -0.25, "wrong min value")
906          self.failUnless(  t == 2, "wrong min tag")
907          self.failUnless(  l == 0., "wrong min location")
908          u=-x[1]*(1.-x[1])*(1.-x[0])
909          t, loc=f.getMinValue(u)
910          p=f.getParametrization(x,t)[0]
911          m, l=loc(u), loc(p)
912          self.failUnless(  m == -0.25, "wrong min value")
913          self.failUnless(  t == 1, "wrong min tag")
914          self.failUnless(  abs(l-0.70710678118654) <= self.EPS,  "wrong min location")
915
916
917       def test_Fault2D(self):
918          f=FaultSystem(dim=2)
919          top1=[ [1.,0.], [1.,1.], [0.,1.] ]
922          self.failUnless(f.getDim() == 2, "wrong dimension")
923          self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
924          self.failUnless(  2. == f.getTotalLength(1), "length wrong")
925          self.failUnless(  0. == f.getMediumDepth(1), "depth wrong")
926          self.failUnless( (0., 2.) ==  f.getW0Range(1)," wrong W0 range")
927          self.failUnless( (0., 0.) ==  f.getW1Range(1)," wrong W1 range")
928          self.failUnless( [0., 1., 2.] ==  f.getW0Offsets(1)," wrong W0 offsets")
929          segs=f.getTopPolyline(1)
930          self.failUnless( len(segs) == 3, "wrong number of segments")
931          self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
932          self.failUnless( segs[0].size == 2, "seg 0 has wrong size.")
933          self.failUnless( numpy.linalg.norm(segs[0]-[1.,0.]) < self.EPS, "wrong vertex. 0 ")
934          self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
935          self.failUnless( segs[1].size == 2, "seg 1 has wrong size.")
936          self.failUnless( numpy.linalg.norm(segs[1]-[1.,1.]) < self.EPS, "wrong vertex. 1 ")
937          self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
938          self.failUnless( segs[2].size == 2, "seg 2 has wrong size.")
939          self.failUnless( numpy.linalg.norm(segs[2]-[0.,1.]) < self.EPS, "wrong vertex. 2 ")
940          c=f.getCenterOnSurface()
941          self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
942          self.failUnless( c.size == 2, "center size is wrong")
943          self.failUnless( numpy.linalg.norm(c-[2./3.,2./3.]) < self.EPS, "center has wrong coordinates.")
944          o=f.getOrientationOnSurface()/pi*180.
945          self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
946
947          top2=[ [10.,0.], [0.,10.] ]
948          f.addFault(V0=[10.,0],strikes=[3.*pi/4],ls=[14.142135623730951], tag=2, w0_offsets=[0,20], w1_max=20)
949          self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
950          self.failUnless(  abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length")
951          self.failUnless(  0. == f.getMediumDepth(2), "depth wrong")
952          self.failUnless( (0., 20.) ==  f.getW0Range(2)," wrong W0 range")
953          self.failUnless( (0., 0.) ==  f.getW1Range(2)," wrong W1 range")
954          self.failUnless( [0., 20.] ==  f.getW0Offsets(2)," wrong W0 offsets")
955          segs=f.getTopPolyline(2)
956          self.failUnless( len(segs) == 2, "wrong number of segments")
957          self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
958          self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.]) < self.EPS, "wrong vertex. 0 ")
959          self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
960          self.failUnless( numpy.linalg.norm(segs[1]-[0.,10.]) < self.EPS, "wrong vertex. 1 ")
961          c=f.getCenterOnSurface()
962          self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
963          self.failUnless( c.size == 2, "center size is wrong")
964          self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.]) < self.EPS, "center has wrong coordinates.")
965          o=f.getOrientationOnSurface()/pi*180.
966          self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
967
968          s,d=f.getSideAndDistance([0.,0.], tag=1)
969          self.failUnless( s<0, "wrong side.")
970          self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
971          s,d=f.getSideAndDistance([0.,2.], tag=1)
972          self.failUnless( s>0, "wrong side.")
973          self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
974          s,d=f.getSideAndDistance([1.,2.], tag=1)
975          self.failUnless( s>0, "wrong side.")
976          self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
977          s,d=f.getSideAndDistance([2.,1.], tag=1)
978          self.failUnless( s>0, "wrong side.")
979          self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
980          s,d=f.getSideAndDistance([2.,0.], tag=1)
981          self.failUnless( s>0, "wrong side.")
982          self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
983          s,d=f.getSideAndDistance([0.,-1.], tag=1)
984          self.failUnless( s<0, "wrong side.")
985          self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
986          s,d=f.getSideAndDistance([-1.,0], tag=1)
987          self.failUnless( s<0, "wrong side.")
988          self.failUnless( abs(d-1.41421356237)<self.EPS, "wrong distance.")
989
990
991          f.transform(rot=-pi/2., shift=[-1.,-1.])
992          self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
993          self.failUnless(  2. == f.getTotalLength(1), "length after transformation wrong")
994          self.failUnless(  0. == f.getMediumDepth(1), "depth after transformation wrong")
995          self.failUnless( (0., 2.) ==  f.getW0Range(1)," wrong W0 after transformation range")
996          self.failUnless( (0., 0.) ==  f.getW1Range(1)," wrong W1 rangeafter transformation ")
997          self.failUnless( [0., 1., 2.] ==  f.getW0Offsets(1)," wrong W0 offsetsafter transformation ")
998          segs=f.getTopPolyline(1)
999          self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
1000          self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1001          self.failUnless( segs[0].size == 2, "seg 0 has wrong size after transformation.")
1002          self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.]) < self.EPS, "wrong vertex. 0  after transformation")
1003          self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex  after transformation1")
1004          self.failUnless( segs[1].size == 2, "seg 1 has wrong size after transformation.")
1005          self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.]) < self.EPS, "wrong vertex.  after transformation1 ")
1006          self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex  after transformation2")
1007          self.failUnless( segs[2].size == 2, "seg 2 has wrong size after transformation.")
1008          self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.]) < self.EPS, "wrong vertex after transformation. 2 ")
1009          self.failUnless(  abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1010          self.failUnless(  0. == f.getMediumDepth(2), "depth wrong after transformation")
1011          self.failUnless( (0., 20.) ==  f.getW0Range(2)," wrong W0 range after transformation")
1012          self.failUnless( (0., 0.) ==  f.getW1Range(2)," wrong W1 range after transformation")
1013          self.failUnless( [0., 20.] ==  f.getW0Offsets(2)," wrong W0 offsets after transformation")
1014          segs=f.getTopPolyline(2)
1015          self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
1016          self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1017          self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9]) < self.EPS, "wrong vertex. 0  after transformation")
1018          self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1019          self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.]) < self.EPS, "wrong vertex. 1  after transformation")
1020
1021          c=f.getCenterOnSurface()
1022          self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1023          self.failUnless( c.size == 2, "center size is wrong")
1024          self.failUnless( numpy.linalg.norm(c-[7./5.,-7./5.]) < self.EPS, "center has wrong coordinates.")
1025          o=f.getOrientationOnSurface()/pi*180.
1026          self.failUnless( abs(o-45.) < self.EPS, "wrong orientation.")
1027
1028          p=f.getParametrization([-1.,0.],1)
1029          self.failUnless(p[1]==1., "wrong value.")
1030          self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1031          p=f.getParametrization([-0.5,0.],1)
1032          self.failUnless(p[1]==1., "wrong value.")
1033          self.failUnless(abs(p[0]-0.5)<self.EPS* 0.5, "wrong value.")
1034          p=f.getParametrization([0.,0.],1)
1035          self.failUnless(p[1]==1., "wrong value.")
1036          self.failUnless(abs(p[0]-1.)<self.EPS, "wrong value.")
1037          p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-8)
1038          self.failUnless(p[1]==0., "wrong value.")
1039          p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-6)
1040          self.failUnless(p[1]==1., "wrong value.")
1041          self.failUnless(abs(p[0]-1.0000001)<self.EPS, "wrong value.")
1042          p=f.getParametrization([0.,0.5],1)
1043          self.failUnless(p[1]==1., "wrong value.")
1044          self.failUnless(abs(p[0]-1.5)<self.EPS, "wrong value.")
1045          p=f.getParametrization([0,1.],1)
1046          self.failUnless(p[1]==1., "wrong value.")
1047          self.failUnless(abs(p[0]-2.)<self.EPS, "wrong value.")
1048          p=f.getParametrization([1.,1.],1)
1049          self.failUnless(p[1]==0., "wrong value.")
1050          p=f.getParametrization([0,1.11],1)
1051          self.failUnless(p[1]==0., "wrong value.")
1052          p=f.getParametrization([-1,-9.],2)
1053          self.failUnless(p[1]==1., "wrong value.")
1054          self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1055          p=f.getParametrization([9,1],2)
1056          self.failUnless(p[1]==1., "wrong value.")
1057          self.failUnless(abs(p[0]-20.)<self.EPS, "wrong value.")
1058
1059       def test_Fault3D(self):
1060          f=FaultSystem(dim=3)
1061          self.failUnless(f.getDim() == 3, "wrong dimension")
1062
1063          top1=[ [0.,0.,0.], [1., 0., 0.] ]
1064          f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=20.,tag=1)
1065          self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
1066          self.failUnless(  1. == f.getTotalLength(1), "length wrong")
1067          self.failUnless(  20. == f.getMediumDepth(1), "depth wrong")
1068          self.failUnless( (0., 1.) ==  f.getW0Range(1)," wrong W0 range")
1069          self.failUnless( (-20., 0.) ==  f.getW1Range(1)," wrong W1 range")
1070          self.failUnless( [0., 1.] ==  f.getW0Offsets(1)," wrong W0 offsets")
1071          segs=f.getTopPolyline(1)
1072          self.failUnless( len(segs) == 2, "wrong number of segments")
1073          self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1074          self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1075          self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1076          self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1077          self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1078          self.failUnless( numpy.linalg.norm(segs[1]-[1.,0.,0]) < self.EPS, "wrong vertex. 1 ")
1079          c=f.getCenterOnSurface()
1080          self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1081          self.failUnless( c.size == 3, "center size is wrong")
1082          self.failUnless( numpy.linalg.norm(c-[0.5,0.,0.]) < self.EPS, "center has wrong coordinates.")
1083          o=f.getOrientationOnSurface()/pi*180.
1084          self.failUnless( abs(o) < self.EPS, "wrong orientation.")
1085          d=f.getDips(1)
1086          self.failUnless( len(d) == 1, "wrong number of dips")
1087          self.failUnless(  abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1088          sn=f.getSegmentNormals(1)
1089          self.failUnless( len(sn) == 1, "wrong number of normals")
1090          self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1091          self.failUnless( numpy.linalg.norm(sn[0]-[0, -1., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1092          dv=f.getDepthVectors(1)
1093          self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1094          self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1095          self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1096          self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1097          self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1098          b=f.getBottomPolyline(1)
1099          self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1100          self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1101          self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1102          self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1103          self.failUnless( numpy.linalg.norm(b[1]-[1., 0., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1104          ds=f.getDepths(1)
1105          self.failUnless( len(ds) == 2, "wrong number of depth")
1106          self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1107          self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1108
1109          top2=[ [0.,0.,0.], [0., 10., 0.] ]
1110          f.addFault(V0=[0.,0,0],strikes=[pi/2],ls=[10.,], dips=pi/2, depths=20.,tag=2)
1111          self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1112          self.failUnless(  10. == f.getTotalLength(2), "length wrong")
1113          self.failUnless(  20. == f.getMediumDepth(2), "depth wrong")
1114          self.failUnless( (0., 10.) ==  f.getW0Range(2)," wrong W0 range")
1115          self.failUnless( (-20., 0.) ==  f.getW1Range(2)," wrong W1 range")
1116          self.failUnless( [0., 10.] ==  f.getW0Offsets(2)," wrong W0 offsets")
1117          segs=f.getTopPolyline(2)
1118          self.failUnless( len(segs) == 2, "wrong number of segments")
1119          self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1120          self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1121          self.failUnless( numpy.linalg.norm(segs[0]-[0.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1122          self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1123          self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1124          self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1125          d=f.getDips(2)
1126          self.failUnless( len(d) == 1, "wrong number of dips")
1127          self.failUnless(  abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1128          sn=f.getSegmentNormals(2)
1129          self.failUnless( len(sn) == 1, "wrong number of normals")
1130          self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1131          self.failUnless( numpy.linalg.norm(sn[0]-[1, 0., 0.]) < self.EPS, "wrong bottom vertex 1 ")
1132          dv=f.getDepthVectors(2)
1133          self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1134          self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1135          self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -20.]) < self.EPS, "wrong depth vector 0 ")
1136          self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1137          self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -20.]) < self.EPS, "wrong depth vector 1 ")
1138          b=f.getBottomPolyline(2)
1139          self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1140          self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1141          self.failUnless( numpy.linalg.norm(b[0]-[0., 0., -20.]) < self.EPS, "wrong bottom vertex 0 ")
1142          self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1143          self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -20.]) < self.EPS, "wrong bottom vertex 1 ")
1144          ds=f.getDepths(2)
1145          self.failUnless( len(ds) == 2, "wrong number of depth")
1146          self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1147          self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1148
1149          top2=[ [10.,0.,0.], [0., 10., 0.] ]
1150          f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/2, depths=30.,tag=2)
1151          self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1152          self.failUnless(  abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1153          self.failUnless(  30. == f.getMediumDepth(2), "depth wrong")
1154          self.failUnless( (-30., 0.) ==  f.getW1Range(2)," wrong W1 range")
1155          segs=f.getTopPolyline(2)
1156          self.failUnless( len(segs) == 2, "wrong number of segments")
1157          self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1158          self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1159          self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1160          self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1161          self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1162          self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1163          d=f.getDips(2)
1164          self.failUnless( len(d) == 1, "wrong number of dips")
1165          self.failUnless(  abs(d[0]-1.5707963267948966) < self.EPS, "wrong dip 0")
1166          sn=f.getSegmentNormals(2)
1167          self.failUnless( len(sn) == 1, "wrong number of normals")
1168          self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1169          self.failUnless( numpy.linalg.norm(sn[0]-[0.70710678118654746, 0.70710678118654746, 0.]) < self.EPS, "wrong bottom vertex 1 ")
1170          dv=f.getDepthVectors(2)
1171          self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1172          self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1173          self.failUnless( numpy.linalg.norm(dv[0]-[0., 0., -30.]) < self.EPS, "wrong depth vector 0 ")
1174          self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1175          self.failUnless( numpy.linalg.norm(dv[1]-[0., 0., -30.]) < self.EPS, "wrong depth vector 1 ")
1176          b=f.getBottomPolyline(2)
1177          self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1178          self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1179          self.failUnless( numpy.linalg.norm(b[0]-[10., 0., -30.]) < self.EPS, "wrong bottom vertex 0 ")
1180          self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1181          self.failUnless( numpy.linalg.norm(b[1]-[0., 10., -30.]) < self.EPS, "wrong bottom vertex 1 ")
1182          ds=f.getDepths(2)
1183          self.failUnless( len(ds) == 2, "wrong number of depth")
1184          self.failUnless( abs(ds[0]-30.) < self.EPS, "wrong depth at vertex 0 ")
1185          self.failUnless( abs(ds[1]-30.) < self.EPS, "wrong depth at vertex 1 ")
1186
1187          top2=[ [10.,0.,0.], [0., 10., 0.] ]
1188          f.addFault(V0=[10.,0,0],strikes=3*pi/4,ls=14.142135623730951, dips=pi/4, depths=50.,tag=2)
1189          self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1190          self.failUnless(  abs(14.142135623730951 - f.getTotalLength(2)) <self.EPS, "length wrong")
1191          self.failUnless(  50. == f.getMediumDepth(2), "depth wrong")
1192          self.failUnless( (-50., 0.) ==  f.getW1Range(2)," wrong W1 range")
1193          segs=f.getTopPolyline(2)
1194          self.failUnless( len(segs) == 2, "wrong number of segments")
1195          self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1196          self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1197          self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1198          self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1199          self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1200          self.failUnless( numpy.linalg.norm(segs[1]-[0., 10., 0]) < self.EPS, "wrong vertex. 1 ")
1201          d=f.getDips(2)
1202          self.failUnless( len(d) == 1, "wrong number of dips")
1203          self.failUnless(  abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1204          sn=f.getSegmentNormals(2)
1205          self.failUnless( len(sn) == 1, "wrong number of normals")
1206          self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1207          self.failUnless( numpy.linalg.norm(sn[0]-[0.5,0.5,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1208          dv=f.getDepthVectors(2)
1209          self.failUnless( len(dv) == 2, "wrong number of depth vectors.")
1210          self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1211          self.failUnless( numpy.linalg.norm(dv[0]-[25., 25., -35.355339059327378]) < self.EPS, "wrong depth vector 0 ")
1212          self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1213          self.failUnless( numpy.linalg.norm(dv[1]-[25.,25., -35.355339059327378]) < self.EPS, "wrong depth vector 1 ")
1214          b=f.getBottomPolyline(2)
1215          self.failUnless( len(b) == 2, "wrong number of bottom vertices")
1216          self.failUnless( isinstance(b[0], numpy.ndarray), "wrong class of bottom vertex 0")
1217          self.failUnless( numpy.linalg.norm(b[0]-[35., 25., -35.355339059327378]) < self.EPS, "wrong bottom vertex 0 ")
1218          self.failUnless( isinstance(b[1], numpy.ndarray), "wrong class of bottom vertex 1")
1219          self.failUnless( numpy.linalg.norm(b[1]-[25, 35., -35.355339059327378]) < self.EPS, "wrong bottom vertex 1 ")
1220          ds=f.getDepths(2)
1221          self.failUnless( len(ds) == 2, "wrong number of depth")
1222          self.failUnless( abs(ds[0]-50.) < self.EPS, "wrong depth at vertex 0 ")
1223          self.failUnless( abs(ds[1]-50.) < self.EPS, "wrong depth at vertex 1 ")
1224
1225          top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
1226          f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=1, dips=pi/4, depths=20.)
1227          self.failUnless(  20. == f.getTotalLength(1), "length wrong")
1228          self.failUnless(  20. == f.getMediumDepth(1), "depth wrong")
1229          segs=f.getTopPolyline(1)
1230          self.failUnless( len(segs) == 3, "wrong number of segments")
1231          self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1232          self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1233          self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.,0.]) < self.EPS, "wrong vertex. 0 ")
1234          self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1235          self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1236          self.failUnless( numpy.linalg.norm(segs[1]-[10.,10.,0.]) < self.EPS, "wrong vertex. 1 ")
1237          self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1238          self.failUnless( segs[2].size == 3, "seg 2 has wrong size.")
1239          self.failUnless( numpy.linalg.norm(segs[2]-[0.,10.,0.]) < self.EPS, "wrong vertex. 2 ")
1240          d=f.getDips(1)
1241          self.failUnless( len(d) == 2, "wrong number of dips")
1242          self.failUnless(  abs(d[0]-0.78539816339744828) < self.EPS, "wrong dip 0")
1243          self.failUnless(  abs(d[1]-0.78539816339744828) < self.EPS, "wrong dip 0")
1244          ds=f.getDepths(1)
1245          self.failUnless( len(ds) == 3, "wrong number of depth")
1246          self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1247          self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1248          sn=f.getSegmentNormals(1)
1249          self.failUnless( len(sn) == 2, "wrong number of normals")
1250          self.failUnless( isinstance(sn[0], numpy.ndarray), "wrong class of bottom vertex 0")
1251          self.failUnless( numpy.linalg.norm(sn[0]-[0.70710678118654746,0.,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1252          self.failUnless( isinstance(sn[1], numpy.ndarray), "wrong class of bottom vertex 0")
1253          self.failUnless( numpy.linalg.norm(sn[1]-[0.,0.70710678118654746,0.70710678118654746]) < self.EPS, "wrong bottom vertex 1 ")
1254          dv=f.getDepthVectors(1)
1255          self.failUnless( len(dv) == 3, "wrong number of depth vectors.")
1256          self.failUnless( isinstance(dv[0], numpy.ndarray), "wrong class of depth vector 0")
1257          self.failUnless( numpy.linalg.norm(dv[0]-[14.142135623730951, 0., -14.142135623730951]) < self.EPS, "wrong depth vector 0 ")
1258          self.failUnless( isinstance(dv[1], numpy.ndarray), "wrong class of depth vector 1")
1259          self.failUnless( numpy.linalg.norm(dv[1]-[11.547005383792515,11.547005383792515, -11.547005383792515]) < self.EPS, "wrong depth vector 2 ")
1260          self.failUnless( isinstance(dv[2], numpy.ndarray), "wrong class of depth vector 1")
1261          self.failUnless( numpy.linalg.norm(dv[2]-[0.,14.142135623730951, -14.142135623730951]) < self.EPS, "wrong depth vector 2 ")
1262          segs=f.getBottomPolyline(1)
1263          self.failUnless( len(segs) == 3, "wrong number of segments")
1264          self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1265          self.failUnless( segs[0].size == 3, "seg 0 has wrong size.")
1266          self.failUnless( numpy.linalg.norm(segs[0]-[24.142135623730951,0.,-14.142135623730951]) < self.EPS, "wrong vertex. 0 ")
1267          self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1268          self.failUnless( segs[1].size == 3, "seg 1 has wrong size.")
1269          self.failUnless( numpy.linalg.norm(segs[1]-[21.547005383792515,21.547005383792515, -11.547005383792515]) < self.EPS, "wrong vertex. 1 ")
1270          self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1271          self.failUnless( segs[2].size == 3, "seg 2 has wrong size.")
1272          self.failUnless( numpy.linalg.norm(segs[2]-[0., 24.142135623730951, -14.142135623730951]) < self.EPS, "wrong vertex. 2 ")
1273          self.failUnless( abs(0.-f.getW0Range(1)[0]) <=self.EPS," wrong W0 range (0)")
1274          self.failUnless( abs(31.857329272664341-f.getW0Range(1)[1]) <=self.EPS," wrong W0 range (1)")
1275          self.failUnless( abs(-20. - f.getW1Range(1)[0]) <=self.EPS," wrong W1 range (0)")
1276          self.failUnless( abs(0. - f.getW1Range(1)[1]) <=self.EPS," wrong W1 range (1)")
1277          self.failUnless( abs(0.0-f.getW0Offsets(1)[0])<=self.EPS," wrong W0 offsets (0)")
1278          self.failUnless( abs(15.92866463633217-f.getW0Offsets(1)[1])<=self.EPS," wrong W0 offsets (1)")
1279          self.failUnless( abs(31.857329272664341-f.getW0Offsets(1)[2])<=self.EPS," wrong W0 offsets(2)")
1280          #
1281          #    ============ fresh start ====================
1282          #
1283          f.addFault(V0=[1.,0,0.],strikes=[pi/2, pi],ls=[1.,1.], dips=pi/2,depths=20,tag=1)
1284          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)
1285          c=f.getCenterOnSurface()
1286          self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1287          self.failUnless( c.size == 3, "center size is wrong")
1288          self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.,0.]) < self.EPS, "center has wrong coordinates.")
1289          o=f.getOrientationOnSurface()/pi*180.
1290          self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1291
1292          f.transform(rot=-pi/2., shift=[-1.,-1.,0.])
1293          self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
1294          self.failUnless(  2. == f.getTotalLength(1), "length after transformation wrong")
1295          self.failUnless(  20. == f.getMediumDepth(1), "depth after transformation wrong")
1296          rw0=f.getW0Range(1)
1297          self.failUnless( len(rw0) ==2, "wo range has wrong length")
1298          self.failUnless( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
1299          self.failUnless( abs(rw0[1]-2.) < self.EPS,"W0 1 wrong.")
1300          self.failUnless( (-20., 0.) ==  f.getW1Range(1)," wrong W1 rangeafter transformation ")
1301          dips=f.getDips(1)
1302          self.failUnless(len(dips) == 2, "wrong number of dips.")
1303          self.failUnless( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
1304          self.failUnless( abs(dips[1]-1.5707963267948966) <= self.EPS, "wrong dip")
1305          ds=f.getDepths(1)
1306          self.failUnless( len(ds) == 3, "wrong number of depth")
1307          self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1308          self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1309          self.failUnless( abs(ds[2]-20.) < self.EPS, "wrong depth at vertex 1 ")
1310          segs=f.getTopPolyline(1)
1311          self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
1312          self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1313          self.failUnless( segs[0].size == 3, "seg 0 has wrong size after transformation.")
1314          self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.,0.]) < self.EPS, "wrong vertex. 0  after transformation")
1315          self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex  after transformation1")
1316          self.failUnless( segs[1].size == 3, "seg 1 has wrong size after transformation.")
1317          self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.,0.]) < self.EPS, "wrong vertex.  after transformation1 ")
1318          self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex  after transformation2")
1319          self.failUnless( segs[2].size == 3, "seg 2 has wrong size after transformation.")
1320          self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.,0.]) < self.EPS, "wrong vertex after transformation. 2 ")
1321          self.failUnless(  abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1322          self.failUnless(  20. == f.getMediumDepth(2), "depth wrong after transformation")
1323          rw0=f.getW0Range(2)
1324          self.failUnless( len(rw0) ==2, "wo range has wrong length")
1325          self.failUnless( abs(rw0[0]) < self.EPS,"W0 0 wrong.")
1326          self.failUnless( abs(rw0[1]-20.) < self.EPS,"W0 1 wrong.")
1327          self.failUnless( (-20., 0.) ==  f.getW1Range(2)," wrong W1 range after transformation")
1328          self.failUnless( [0., 20.] ==  f.getW0Offsets(2)," wrong W0 offsets after transformation")
1329          dips=f.getDips(2)
1330          self.failUnless(len(dips) == 1, "wrong number of dips.")
1331          self.failUnless( abs(dips[0]-1.5707963267948966) <= self.EPS, "wrong dip")
1332          ds=f.getDepths(2)
1333          self.failUnless( len(ds) == 2, "wrong number of depth")
1334          self.failUnless( abs(ds[0]-20.) < self.EPS, "wrong depth at vertex 0 ")
1335          self.failUnless( abs(ds[1]-20.) < self.EPS, "wrong depth at vertex 1 ")
1336          segs=f.getTopPolyline(2)
1337          self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
1338          self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1339          self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9,0.]) < self.EPS, "wrong vertex. 0  after transformation")
1340          self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1341          self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.,0.]) < self.EPS, "wrong vertex. 1  after transformation")
1342          #
1343          #    ============ fresh start ====================
1344          #
1345          f=FaultSystem(dim=3)
1346
1347          top1=[ [0.,0.,0.], [1., 0., 0.] ]
1348          f.addFault(V0=[0.,0,0],strikes=[0.],ls=[1.,], dips=pi/2, depths=1.,tag=1)
1349          top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
1350          f.addFault(V0=[10.,0.,0.],strikes=[pi/2, pi],ls=[10.,10.],tag=2, dips=pi/4, depths=20.)
1351
1352          p,m=f.getParametrization([0.3,0.,-0.5],1)
1353          self.failUnless(length(p-[0.3,-0.5]) <= self.EPS, "wrong value.")
1354          self.failUnless(m==1., "wrong value.")
1355
1356          p,m=f.getParametrization([0.5,0.,-0.5],1)
1357          self.failUnless(length(p-[0.5,-0.5]) <= self.EPS, "wrong value.")
1358          self.failUnless(m==1., "wrong value.")
1359
1360          p,m=f.getParametrization([0.25,0.,-0.5],1)
1361          self.failUnless(length(p-[0.25,-0.5]) <= self.EPS, "wrong value.")
1362          self.failUnless(m==1., "wrong value.")
1363
1364          p,m=f.getParametrization([0.5,0.,-0.25],1)
1365          self.failUnless(length(p-[0.5,-0.25]) <= self.EPS, "wrong value.")
1366          self.failUnless(m==1., "wrong value.")
1367
1368          p,m=f.getParametrization([0.001,0.,-0.001],1)
1369          self.failUnless(length(p-[0.001, -0.001]) <= self.EPS, "wrong value.")
1370          self.failUnless(m==1., "wrong value.")
1371
1372          p,m=f.getParametrization([0.001,0.,0.001],1)
1373          self.failUnless(m==0., "wrong value.")
1374
1375          p,m=f.getParametrization([0.999,0.,0.001],1)
1376          self.failUnless(m==0., "wrong value.")
1377
1378          p,m=f.getParametrization([1.001,0.,-0.001],1)
1379          self.failUnless(m==0., "wrong value.")
1380          p,m=f.getParametrization([1.001,0.,-0.1],1)
1381          self.failUnless(m==0., "wrong value.")
1382          p,m=f.getParametrization([1.001,0.,-0.000000001],1)
1383          self.failUnless(m==0., "wrong value.")
1384
1385          p,m=f.getParametrization([0.999,0.,-0.001],1)
1386          self.failUnless(length(p-[0.999, -0.001]) <= self.EPS, "wrong value.")
1387          self.failUnless(m==1., "wrong value.")
1388
1389          p,m=f.getParametrization([ 16.29252873 , 6.46410161 ,-6.29252873], 2, tol=1.e-7)
1390          self.failUnless(m==1., "wrong value.")
1391          self.failUnless(length(p-[4.7785993908996511, -10]) <= self.EPS*10., "wrong value.")
1392          p,m=f.getParametrization([15.77350269, 12.77350269, -5.77350269], 2, tol=1.e-7)
1393          self.failUnless(m==1., "wrong value.")
1394          self.failUnless(length(p-[11.150065245432518, -10]) <= self.EPS*10., "wrong value.")
1395
1396          p,m=f.getParametrization([  3., 17.0710678, -7.0710678], 2, tol=1.e-7)
1397          self.failUnless(m==1., "wrong value.")
1398          self.failUnless(length(p-[27.078729881764687, -10]) <= self.EPS*10., "wrong value.")
1399          p,m=f.getParametrization([9.30940108, 16.55204176, -6.55204176], 2, tol=1.e-7)
1400          self.failUnless(m==1., "wrong value.")
1401          self.failUnless(length(p-[20.707264027231822, -10]) <= self.EPS*10., "wrong value.")
1402
1403          p,m=f.getParametrization([ 21.54700538,  21.54700538, -11.54700538], 2, tol=1.e-7)
1404          self.failUnless(m==1., "wrong value.")
1405          self.failUnless(length(p-[15.92866463633217, -20]) <= self.EPS*10., "wrong value.")
1406
1407          p,m=f.getParametrization([ 0.,0.,0.], 2, tol=1.e-7)
1408          self.failUnless(m==0., "wrong value.")
1409
1410          p,m=f.getParametrization([ 11.,11.,0.], 2, tol=1.e-7)
1411          self.failUnless(m==0., "wrong value.")
1412
1413
1414          s,d=f.getSideAndDistance([0.,-1.,0.], tag=1)
1415          self.failUnless( s>0, "wrong side.")
1416          self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1417          s,d=f.getSideAndDistance([1.,-1.,0.], tag=1)
1418          self.failUnless( s>0, "wrong side.")
1419          self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1420          s,d=f.getSideAndDistance([0.,1.,0.], tag=1)
1421          self.failUnless( s<0, "wrong side.")
1422          self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1423          s,d=f.getSideAndDistance([1.,1.,0.], tag=1)
1424          self.failUnless( s<0, "wrong side.")
1425          self.failUnless( abs(d-1.)<self.EPS, "wrong distance.")
1426
1427
1428          s,d=f.getSideAndDistance([0.,0.,0.], tag=2)
1429          self.failUnless( s<0, "wrong side.")
1430          self.failUnless( abs(d-10.)<self.EPS, "wrong distance.")
1431          s,d=f.getSideAndDistance([5.,5.,0.], tag=2)
1432          self.failUnless( s<0, "wrong side.")
1433          self.failUnless( abs(d-5.)<self.EPS, "wrong distance.")
1434
1435          s,d=f.getSideAndDistance([10.,10.,-1.], tag=2)
1436          self.failUnless( s<0, "wrong side.")
1437          self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
1438          s,d=f.getSideAndDistance([10.,10.,-2.], tag=2)
1439          self.failUnless( s<0, "wrong side.")
1440          self.failUnless( abs(d-2.*0.70710678118654757)<self.EPS, "wrong distance.")
1441          s,d=f.getSideAndDistance([10.,10.,-3.], tag=2)
1442          self.failUnless( s<0, "wrong side.")
1443          self.failUnless( abs(d-3.*0.70710678118654757)<self.EPS, "wrong distance.")
1444
1445          s,d=f.getSideAndDistance([5.,12.,0], tag=2)
1446          self.failUnless( s>0, "wrong side.")
1447          self.failUnless( abs(d-2*0.70710678118654757)<self.EPS, "wrong distance.")
1448          s,d=f.getSideAndDistance([5.,12.,-1], tag=2)
1449          self.failUnless( s>0, "wrong side.")
1450          self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
1451          s,d=f.getSideAndDistance([5.,12.,-2], tag=2)
1452          # s not checked as it is undefined.
1453          self.failUnless( abs(d)<self.EPS, "wrong distance.")
1454          s,d=f.getSideAndDistance([5.,12.,-3], tag=2)
1455          self.failUnless( s<0, "wrong side.")
1456          self.failUnless( abs(d-0.70710678118654757)<self.EPS, "wrong distance.")
1457          s,d=f.getSideAndDistance([5.,12.,-4], tag=2)
1458          self.failUnless( s<0, "wrong side.")
1459          self.failUnless( abs(d-2.*0.70710678118654757)<self.EPS, "wrong distance.")
1460
1461  if __name__ == '__main__':  if __name__ == '__main__':
1462     suite = unittest.TestSuite()     suite = unittest.TestSuite()