/[escript]/trunk/finley/test/python/run_models.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2438 by gross, Tue May 26 02:26:52 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"""
# Line 22  __url__="https://launchpad.net/escript-f Line 23  __url__="https://launchpad.net/escript-f
23  import unittest  import unittest
24  import tempfile  import tempfile
25                
 from esys.escript import *  
 from esys.finley import Rectangle  
 from esys.escript.models import DarcyFlow  
 import sys  
 import os  
 try:  
      FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']  
 except KeyError:  
      FINLEY_WORKDIR='.'  
26    
27    
28  VERBOSE=False # or True  VERBOSE=False  and True
 DETAIL_VERBOSE=False  
29    
30  from esys.escript import *  from esys.escript import *
31  from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian  from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian, FaultSystem, DarcyFlow
32    from esys.escript.models import Mountains
33  from esys.finley import Rectangle, Brick  from esys.finley import Rectangle, Brick
34    
 from esys.escript.models import Mountains  
35  from math import pi  from math import pi
36    import numpy
37    import sys
38    import os
39    #====================================================================================================================
40    try:
41         FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
42    except KeyError:
43         FINLEY_WORKDIR='.'
44    
45  #====================================================================================================================  #====================================================================================================================
46  class Test_StokesProblemCartesian2D(unittest.TestCase):  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 66  class Test_StokesProblemCartesian2D(unit Line 65  class Test_StokesProblemCartesian2D(unit
65                
66         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
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=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=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 92  class Test_StokesProblemCartesian2D(unit Line 91  class Test_StokesProblemCartesian2D(unit
91                
92         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
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=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=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 118  class Test_StokesProblemCartesian2D(unit Line 116  class Test_StokesProblemCartesian2D(unit
116                
117         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
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=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=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 144  class Test_StokesProblemCartesian2D(unit Line 143  class Test_StokesProblemCartesian2D(unit
143                
144         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
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=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=50,usePCG=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 170  class Test_StokesProblemCartesian2D(unit Line 169  class Test_StokesProblemCartesian2D(unit
169                
170         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
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=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=20,usePCG=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 196  class Test_StokesProblemCartesian2D(unit Line 196  class Test_StokesProblemCartesian2D(unit
196                
197         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
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=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=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 211  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 233  class Test_StokesProblemCartesian3D(unit Line 233  class Test_StokesProblemCartesian3D(unit
233                
234         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         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.]         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=DETAIL_VERBOSE, verbose=VERBOSE ,max_iter=100,usePCG=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 264  class Test_StokesProblemCartesian3D(unit Line 264  class Test_StokesProblemCartesian3D(unit
264                
265         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         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.]         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*0.1)         sp.setTolerance(self.TOL)
269         u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE ,max_iter=100,usePCG=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 275  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 293  class Test_StokesProblemCartesian3D(unit Line 294  class Test_StokesProblemCartesian3D(unit
294                
295         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         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.]         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=DETAIL_VERBOSE, verbose=VERBOSE ,max_iter=100,usePCG=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 325  class Test_StokesProblemCartesian3D(unit Line 326  class Test_StokesProblemCartesian3D(unit
326                
327         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         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.]         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=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=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])
# Line 355  class Test_StokesProblemCartesian3D(unit Line 356  class Test_StokesProblemCartesian3D(unit
356                
357         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         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.]         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*0.1)         sp.setTolerance(self.TOL/10)
361         u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE and False, verbose=VERBOSE,max_iter=100,usePCG=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])
# Line 385  class Test_StokesProblemCartesian3D(unit Line 386  class Test_StokesProblemCartesian3D(unit
386                
387         sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)         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.]         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=DETAIL_VERBOSE, verbose=VERBOSE ,max_iter=100,usePCG=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 398  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  #====================================================================================================================  #====================================================================================================================
 class Test_Darcy(unittest.TestCase):  
     # this is a simple test for the darcy flux problem  
     #  
     #  
     #  p = 1/k * ( 1/2* (fb-f0)/xb* x **2 + f0 * x - ub*x ) +  p0  
     #  
     #  with f = (fb-f0)/xb* x + f0  
     #  
     #    u = f - k * p,x = ub  
     #  
     #  we prescribe pressure at x=x0=0 to p0  
     #  
     #  if we prescribe the pressure on the bottom x=xb we set  
     #  
     #  pb= 1/k * ( 1/2* (fb-f0)* xb + f0 * xb - ub*xb ) +  p0 = 1/k * ((fb+f0)/2 - ub ) * xb  + p0  
     #  
     #  which leads to ub = (fb+f0)/2-k*(pb-p0)/xb  
     #  
     def rescaleDomain(self):  
         x=self.dom.getX().copy()  
         for i in xrange(self.dom.getDim()):  
              x_inf=inf(x[i])  
              x_sup=sup(x[i])  
              if i == self.dom.getDim()-1:  
                 x[i]=-self.WIDTH*(x[i]-x_sup)/(x_inf-x_sup)  
              else:  
                 x[i]=self.WIDTH*(x[i]-x_inf)/(x_sup-x_inf)  
         self.dom.setX(x)  
     def getScalarMask(self,include_bottom=True):  
         x=self.dom.getX().copy()  
         x_inf=inf(x[self.dom.getDim()-1])  
         x_sup=sup(x[self.dom.getDim()-1])  
         out=whereZero(x[self.dom.getDim()-1]-x_sup)  
         if include_bottom: out+=whereZero(x[self.dom.getDim()-1]-x_inf)  
         return wherePositive(out)  
     def getVectorMask(self,include_bottom=True):  
         x=self.dom.getX().copy()  
         out=Vector(0.,Solution(self.dom))  
         for i in xrange(self.dom.getDim()):  
              x_inf=inf(x[i])  
              x_sup=sup(x[i])  
              if i != self.dom.getDim()-1: out[i]+=whereZero(x[i]-x_sup)  
              if i != self.dom.getDim()-1 or include_bottom: out[i]+=whereZero(x[i]-x_inf)  
         return wherePositive(out)  
   
     def setSolutionFixedBottom(self, p0, pb, f0, fb, k):  
          d=self.dom.getDim()  
          x=self.dom.getX()[d-1]  
          xb=inf(x)  
          u=Vector(0.,Solution(self.dom))+kronecker(d)[d-1]*((f0+fb)/2.-k*(pb-p0)/xb)  
          p=1./k*((fb-f0)/(xb*2.)* x**2 - (fb-f0)/2.*x)+(pb-p0)/xb*x +  p0  
          f= ((fb-f0)/xb* x + f0)*kronecker(Function(self.dom))[d-1]  
          return u,p,f  
           
     def testConstF_FixedBottom_smallK(self):  
         k=1.e-10  
         mp=self.getScalarMask(include_bottom=True)  
         mv=self.getVectorMask(include_bottom=False)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)  
         p=p_ref*mp  
         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)  
         df.setSubProblemTolerance()  
         v,p=df.solve(u_ref,p, max_iter=100, verbose=VERBOSE)  
         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.")  
     def testConstF_FixedBottom_mediumK(self):  
         k=1.  
         mp=self.getScalarMask(include_bottom=True)  
         mv=self.getVectorMask(include_bottom=False)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)  
         p=p_ref*mp  
         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)  
         df.setSubProblemTolerance()  
         v,p=df.solve(u,p,max_iter=100, verbose=VERBOSE )  
         self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")  
         self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")  
   
     def testConstF_FixedBottom_largeK(self):  
         k=1.e10  
         mp=self.getScalarMask(include_bottom=True)  
         mv=self.getVectorMask(include_bottom=False)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)  
         p=p_ref*mp  
         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)  
         df.setSubProblemTolerance()  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         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.")  
   
     def testVarioF_FixedBottom_smallK(self):  
         k=1.e-10  
         mp=self.getScalarMask(include_bottom=True)  
         mv=self.getVectorMask(include_bottom=False)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)  
         p=p_ref*mp  
         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)  
         df.setSubProblemTolerance()  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         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.")  
   
     def testVarioF_FixedBottom_mediumK(self):  
         k=1.  
         mp=self.getScalarMask(include_bottom=True)  
         mv=self.getVectorMask(include_bottom=False)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)  
         p=p_ref*mp  
         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)  
         df.setSubProblemTolerance()  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         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.")  
   
     def testVarioF_FixedBottom_largeK(self):  
         k=1.e10  
         mp=self.getScalarMask(include_bottom=True)  
         mv=self.getVectorMask(include_bottom=False)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)  
         p=p_ref*mp  
         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)  
         df.setSubProblemTolerance()  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         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.")  
   
     def testConstF_FreeBottom_smallK(self):  
         k=1.e-10  
         mp=self.getScalarMask(include_bottom=False)  
         mv=self.getVectorMask(include_bottom=True)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)  
         p=p_ref*mp  
         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)  
         df.setSubProblemTolerance()  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         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.")  
   
     def testConstF_FreeBottom_mediumK(self):  
         k=1.  
         mp=self.getScalarMask(include_bottom=False)  
         mv=self.getVectorMask(include_bottom=True)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)  
         p=p_ref*mp  
         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)  
         df.setSubProblemTolerance()  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         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.")  
   
     def testConstF_FreeBottom_largeK(self):  
         k=1.e10  
         mp=self.getScalarMask(include_bottom=False)  
         mv=self.getVectorMask(include_bottom=True)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)  
         p=p_ref*mp  
         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)  
         df.setSubProblemTolerance()  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         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.")  
   
     def testVarioF_FreeBottom_smallK(self):  
         k=1.e-10  
         mp=self.getScalarMask(include_bottom=False)  
         mv=self.getVectorMask(include_bottom=True)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)  
         p=p_ref*mp  
         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)  
         df.setSubProblemTolerance()  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         self.failUnless(Lsup(v-u_ref)<self.TOL*25.*Lsup(u_ref), "flux error too big.")  # 25 because of disc. error.  
         self.failUnless(Lsup(p-p_ref)<self.TOL*25.*Lsup(p_ref), "pressure error too big.")  
   
     def testVarioF_FreeBottom_mediumK(self):  
         k=1.  
         mp=self.getScalarMask(include_bottom=False)  
         mv=self.getVectorMask(include_bottom=True)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)  
         p=p_ref*mp  
         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)  
         df.setSubProblemTolerance()  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         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.")  
   
     def testVarioF_FreeBottom_largeK(self):  
         k=1.e10  
         mp=self.getScalarMask(include_bottom=False)  
         mv=self.getVectorMask(include_bottom=True)  
         u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)  
         p=p_ref*mp  
         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)  
         df.setSubProblemTolerance()  
         v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)  
         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):  
         del self.dom  
   
402    
403  class Test_Mountains3D(unittest.TestCase):  class Test_Mountains3D(unittest.TestCase):
404     def setUp(self):     def setUp(self):
# Line 709  class Test_Mountains3D(unittest.TestCase Line 422  class Test_Mountains3D(unittest.TestCase
422         v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])* cos(pi*n2*x[2])         v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])* cos(pi*n2*x[2])
423         v[2]=a2*cos(pi*n0*x[0])* cos(pi*n1*x[1])* sin(pi*n2*x[2])         v[2]=a2*cos(pi*n0*x[0])* cos(pi*n1*x[1])* sin(pi*n2*x[2])
424    
425         H_t=Scalar(0.0, Solution(self.domain))         mts=Mountains(self.domain,eps=EPS)
426         mts=Mountains(self.domain,v,eps=EPS,z=1)         mts.setVelocity(v)
427         u,Z=mts.update(u=v,H_t=H_t)         Z=mts.update()
428                
429         error_int=integrate(Z)         error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[2]-1.)))
430         self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")         self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
431    
432  class Test_Mountains2D(unittest.TestCase):  class Test_Mountains2D(unittest.TestCase):
# Line 736  class Test_Mountains2D(unittest.TestCase Line 449  class Test_Mountains2D(unittest.TestCase
449         v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])         v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])
450                
451         H_t=Scalar(0.0, Solution(self.domain))         H_t=Scalar(0.0, Solution(self.domain))
452         mts=Mountains(self.domain,v,eps=EPS,z=1)         mts=Mountains(self.domain,eps=EPS)
453         u,Z=mts.update(u=v,H_t=H_t)         mts.setVelocity(v)
454           Z=mts.update()
455                
456         error_int=integrate(Z)         error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[1]-1.)))
457         self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")         self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
458                
459    
# Line 925  class Test_Rheologies(unittest.TestCase) Line 639  class Test_Rheologies(unittest.TestCase)
639           for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])           for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
640    
641  class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase):  class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase):
642     TOL=1.e-6     TOL=1.e-5
643     VERBOSE=False     VERBOSE=False # or True
644     A=1.     A=1.
645     P_max=100     P_max=100
646     NE=2*getMPISizeWorld()     NE=2*getMPISizeWorld()
# Line 978  class Test_IncompressibleIsotropicFlowCa Line 692  class Test_IncompressibleIsotropicFlowCa
692        mod.setElasticShearModulus(self.mu)        mod.setElasticShearModulus(self.mu)
693        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])        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        mod.setTolerance(self.TOL)        mod.setTolerance(self.TOL)
695        mod.setFlowSubTolerance(self.TOL**2)        mod.setEtaTolerance(self.TOL*0.1)
       mod.setFlowTolerance(self.TOL)  
696    
697        BF=Vector(self.P_max,Function(self.dom))        BF=Vector(self.P_max,Function(self.dom))
698        for d in range(self.dom.getDim()):        for d in range(self.dom.getDim()):
# Line 1004  class Test_IncompressibleIsotropicFlowCa Line 717  class Test_IncompressibleIsotropicFlowCa
717           t_ref=t+dt           t_ref=t+dt
718           v_ref, s_ref,p_ref=self.getReference(t_ref)           v_ref, s_ref,p_ref=self.getReference(t_ref)
719           mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref)           mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref)
720           mod.update(dt, iter_max=100, inner_iter_max=20, verbose=self.VERBOSE, usePCG=True)           # mod.update(dt, eta_iter_max=10, iter_max=50, verbose=self.VERBOSE, usePCG=True, max_correction_steps=30) new version
721             mod.update(dt, iter_max=50, verbose=self.VERBOSE, usePCG=True)
722           self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)           self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)
723           t+=dt           t+=dt
724           N_t+=1           N_t+=1
# Line 1018  class Test_IncompressibleIsotropicFlowCa Line 732  class Test_IncompressibleIsotropicFlowCa
732           error_t=abs(mod.getTime()-t_ref)/abs(t_ref)           error_t=abs(mod.getTime()-t_ref)/abs(t_ref)
733           if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v           if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v
734           self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )           self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
          self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )  
735           self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )           self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )
736           self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )           self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )
737             self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )
738     def tearDown(self):     def tearDown(self):
739          del self.dom          del self.dom
740    
# Line 1106  class Test_IncompressibleIsotropicFlowCa Line 820  class Test_IncompressibleIsotropicFlowCa
820         self.latestart=False         self.latestart=False
821         self.runIt(free=0)         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.] ]
920          self.failUnlessRaises(ValueError,f.addFault,V0=[1.,0],strikes=[pi/2, pi/2],ls=[1.,1.],tag=1,dips=top1)
921          f.addFault(V0=[1.,0],strikes=[pi/2, pi],ls=[1.,1.],tag=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()
1463     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))     # suite.addTest(unittest.makeSuite(Test_FaultSystem))
1464     suite.addTest(unittest.makeSuite(Test_Darcy3D))     # suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
    suite.addTest(unittest.makeSuite(Test_Darcy2D))  
1465     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
1466     suite.addTest(unittest.makeSuite(Test_Mountains3D))     suite.addTest(unittest.makeSuite(Test_Mountains3D))
1467     suite.addTest(unittest.makeSuite(Test_Mountains2D))     suite.addTest(unittest.makeSuite(Test_Mountains2D))
# Line 1118  if __name__ == '__main__': Line 1469  if __name__ == '__main__':
1469     suite.addTest(unittest.makeSuite(Test_IncompressibleIsotropicFlowCartesian))     suite.addTest(unittest.makeSuite(Test_IncompressibleIsotropicFlowCartesian))
1470     s=unittest.TextTestRunner(verbosity=2).run(suite)     s=unittest.TextTestRunner(verbosity=2).run(suite)
1471     if not s.wasSuccessful(): sys.exit(1)     if not s.wasSuccessful(): sys.exit(1)
   

Legend:
Removed from v.2438  
changed lines
  Added in v.3527

  ViewVC Help
Powered by ViewVC 1.1.26