/[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 2647 by gross, Fri Sep 4 05:25:25 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-2009 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-2009 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 24  import tempfile Line 25  import tempfile
25                
26    
27    
28  VERBOSE=False # or True  VERBOSE=False  and True
29    
30  from esys.escript import *  from esys.escript import *
31  from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian, FaultSystem, DarcyFlow  from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian, FaultSystem, DarcyFlow
# Line 46  class Test_StokesProblemCartesian2D(unit Line 47  class Test_StokesProblemCartesian2D(unit
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 67  class Test_StokesProblemCartesian2D(unit
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,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 91  class Test_StokesProblemCartesian2D(unit Line 92  class Test_StokesProblemCartesian2D(unit
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, 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
# Line 118  class Test_StokesProblemCartesian2D(unit Line 119  class Test_StokesProblemCartesian2D(unit
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, 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 168  class Test_StokesProblemCartesian2D(unit Line 170  class Test_StokesProblemCartesian2D(unit
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, 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 208  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 262  class Test_StokesProblemCartesian3D(unit Line 265  class Test_StokesProblemCartesian3D(unit
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, 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])
# Line 293  class Test_StokesProblemCartesian3D(unit Line 296  class Test_StokesProblemCartesian3D(unit
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, 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 354  class Test_StokesProblemCartesian3D(unit Line 357  class Test_StokesProblemCartesian3D(unit
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, 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])
# 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  #====================================================================================================================  #====================================================================================================================
 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)  
         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)  
         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)  
         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)  
         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)  
         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)  
         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)  
         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)  
         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)  
         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)  
         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)  
         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)  
         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 912  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 965  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.setFlowTolerance(self.TOL)        mod.setEtaTolerance(self.TOL*0.1)
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 990  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 1004  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 1096  class Test_IncompressibleIsotropicFlowCa Line 824  class Test_IncompressibleIsotropicFlowCa
824  class Test_FaultSystem(unittest.TestCase):  class Test_FaultSystem(unittest.TestCase):
825     EPS=1.e-8     EPS=1.e-8
826     NE=10     NE=10
827     def test_Fault2D_MaxValue(self):     def test_Fault_MaxValue(self):
828        dom=Rectangle(2*self.NE,2*self.NE)        dom=Rectangle(2*self.NE,2*self.NE)
829        x=dom.getX()        x=dom.getX()
830        f=FaultSystem(dim=2)        f=FaultSystem(dim=2)
831        f.addFault([[0.5,0.],[0.,0.5]], tag=1)        f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1)
832        f.addFault([[1.,0.5],[0.5,0.5],[0.5,1.]], tag=2)        f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2)
833    
834        m, t, l=f.getMaxValue(x[0]*(1.-x[0])*(1-x[1]))        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")        self.failUnless(  m == 0.25, "wrong max value")
839        self.failUnless(  t == 1, "wrong max tag")        self.failUnless(  t == 1, "wrong max tag")
840        self.failUnless(  l == 0., "wrong max location")        self.failUnless(  l == 0., "wrong max location")
841        m, t, l=f.getMaxValue(x[1]*(1.-x[1])*(1-x[0])*x[0])  
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")        self.failUnless(  m == 0.0625, "wrong max value")
847        self.failUnless(  t == 2, "wrong max tag")        self.failUnless(  t == 2, "wrong max tag")
848        self.failUnless(  l == 0.5, "wrong max location")        self.failUnless(  l == 0.5, "wrong max location")
849        m, t, l=f.getMaxValue(x[0]*(1.-x[0])*x[1])  
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")        self.failUnless(  m == 0.25, "wrong max value")
855        self.failUnless(  t == 2, "wrong max tag")        self.failUnless(  t == 2, "wrong max tag")
856        self.failUnless(  l == 1.0, "wrong max location")        self.failUnless(  l == 1.0, "wrong max location")
857        m, t, l= f.getMaxValue(x[1]*(1.-x[1])*x[0])  
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")        self.failUnless(  m == 0.25, "wrong max value")
863        self.failUnless(  t == 2, "wrong max tag")        self.failUnless(  t == 2, "wrong max tag")
864        self.failUnless(  l == 0., "wrong max location")        self.failUnless(  l == 0., "wrong max location")
865        m, t, l= f.getMaxValue(x[1]*(1.-x[1])*(1.-x[0]))  
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")        self.failUnless(  m == 0.25, "wrong max value")
871        self.failUnless(  t == 1, "wrong max tag")        self.failUnless(  t == 1, "wrong max tag")
872        self.failUnless(  abs(l-0.70710678118654) <= self.EPS,  "wrong max location")        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 xtest_Fault2D_TwoFaults(self):     def test_Fault2D(self):
918        f=FaultSystem(dim=2)        f=FaultSystem(dim=2)
919        top1=[ [1.,0.], [1.,1.], [0.,1.] ]        top1=[ [1.,0.], [1.,1.], [0.,1.] ]
920        self.failUnlessRaises(ValueError,f.addFault,top=top1,tag=1,bottom=top1)        self.failUnlessRaises(ValueError,f.addFault,V0=[1.,0],strikes=[pi/2, pi/2],ls=[1.,1.],tag=1,dips=top1)
921        f.addFault(top=top1,tag=1)        f.addFault(V0=[1.,0],strikes=[pi/2, pi],ls=[1.,1.],tag=1)
922        self.failUnless(f.getDim() == 2, "wrong dimension")        self.failUnless(f.getDim() == 2, "wrong dimension")
923        self.failUnless( [ 1 ] == f.getTags(), "tags wrong")        self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
924        self.failUnless(  2. == f.getLength(1), "length wrong")        self.failUnless(  2. == f.getTotalLength(1), "length wrong")
925        self.failUnless(  0. == f.getDepth(1), "depth wrong")        self.failUnless(  0. == f.getMediumDepth(1), "depth wrong")
926        self.failUnless( (0., 2.) ==  f.getW0Range(1)," wrong W0 range")        self.failUnless( (0., 2.) ==  f.getW0Range(1)," wrong W0 range")
927        self.failUnless( (0., 0.) ==  f.getW1Range(1)," wrong W1 range")        self.failUnless( (0., 0.) ==  f.getW1Range(1)," wrong W1 range")
928        self.failUnless( [0., 1., 2.] ==  f.getW0Offsets(1)," wrong W0 offsets")        self.failUnless( [0., 1., 2.] ==  f.getW0Offsets(1)," wrong W0 offsets")
929        segs=f.getFaultSegments(1)[0]        segs=f.getTopPolyline(1)
930        self.failUnless( len(segs) == 3, "wrong number of segments")        self.failUnless( len(segs) == 3, "wrong number of segments")
931        self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")        self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
932        self.failUnless( segs[0].size == 2, "seg 0 has wrong size.")        self.failUnless( segs[0].size == 2, "seg 0 has wrong size.")
# Line 1153  class Test_FaultSystem(unittest.TestCase Line 942  class Test_FaultSystem(unittest.TestCase
942        self.failUnless( c.size == 2, "center size is wrong")        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.")        self.failUnless( numpy.linalg.norm(c-[2./3.,2./3.]) < self.EPS, "center has wrong coordinates.")
944        o=f.getOrientationOnSurface()/pi*180.        o=f.getOrientationOnSurface()/pi*180.
945        self.failUnless( abs(o-45.) < self.EPS, "wrong orientation.")        self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
946    
947        top2=[ [10.,0.], [0.,10.] ]        top2=[ [10.,0.], [0.,10.] ]
948        f.addFault(top2, tag=2, w0_offsets=[0,20], w1_max=20)        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")        self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
950        self.failUnless(  abs(f.getLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length")        self.failUnless(  abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length")
951        self.failUnless(  0. == f.getDepth(2), "depth wrong")        self.failUnless(  0. == f.getMediumDepth(2), "depth wrong")
952        self.failUnless( (0., 20.) ==  f.getW0Range(2)," wrong W0 range")        self.failUnless( (0., 20.) ==  f.getW0Range(2)," wrong W0 range")
953        self.failUnless( (0., 0.) ==  f.getW1Range(2)," wrong W1 range")        self.failUnless( (0., 0.) ==  f.getW1Range(2)," wrong W1 range")
954        self.failUnless( [0., 20.] ==  f.getW0Offsets(2)," wrong W0 offsets")        self.failUnless( [0., 20.] ==  f.getW0Offsets(2)," wrong W0 offsets")
955        segs=f.getFaultSegments(2)[0]        segs=f.getTopPolyline(2)
956        self.failUnless( len(segs) == 2, "wrong number of segments")        self.failUnless( len(segs) == 2, "wrong number of segments")
957        self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")        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 ")        self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.]) < self.EPS, "wrong vertex. 0 ")
# Line 1174  class Test_FaultSystem(unittest.TestCase Line 963  class Test_FaultSystem(unittest.TestCase
963        self.failUnless( c.size == 2, "center size is wrong")        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.")        self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.]) < self.EPS, "center has wrong coordinates.")
965        o=f.getOrientationOnSurface()/pi*180.        o=f.getOrientationOnSurface()/pi*180.
966        self.failUnless( abs(o-45.) < self.EPS, "wrong orientation.")        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.])        f.transform(rot=-pi/2., shift=[-1.,-1.])
992        self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")        self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
993        self.failUnless(  2. == f.getLength(1), "length after transformation wrong")        self.failUnless(  2. == f.getTotalLength(1), "length after transformation wrong")
994        self.failUnless(  0. == f.getDepth(1), "depth after transformation wrong")        self.failUnless(  0. == f.getMediumDepth(1), "depth after transformation wrong")
995        self.failUnless( (0., 2.) ==  f.getW0Range(1)," wrong W0 after transformation range")        self.failUnless( (0., 2.) ==  f.getW0Range(1)," wrong W0 after transformation range")
996        self.failUnless( (0., 0.) ==  f.getW1Range(1)," wrong W1 rangeafter transformation ")        self.failUnless( (0., 0.) ==  f.getW1Range(1)," wrong W1 rangeafter transformation ")
997        self.failUnless( [0., 1., 2.] ==  f.getW0Offsets(1)," wrong W0 offsetsafter transformation ")        self.failUnless( [0., 1., 2.] ==  f.getW0Offsets(1)," wrong W0 offsetsafter transformation ")
998        segs=f.getFaultSegments(1)[0]        segs=f.getTopPolyline(1)
999        self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")        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")        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.")        self.failUnless( segs[0].size == 2, "seg 0 has wrong size after transformation.")
# Line 1194  class Test_FaultSystem(unittest.TestCase Line 1006  class Test_FaultSystem(unittest.TestCase
1006        self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex  after transformation2")        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.")        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 ")        self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.]) < self.EPS, "wrong vertex after transformation. 2 ")
1009        self.failUnless(  abs(f.getLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")        self.failUnless(  abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1010        self.failUnless(  0. == f.getDepth(2), "depth wrong after transformation")        self.failUnless(  0. == f.getMediumDepth(2), "depth wrong after transformation")
1011        self.failUnless( (0., 20.) ==  f.getW0Range(2)," wrong W0 range after transformation")        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")        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")        self.failUnless( [0., 20.] ==  f.getW0Offsets(2)," wrong W0 offsets after transformation")
1014        segs=f.getFaultSegments(2)[0]        segs=f.getTopPolyline(2)
1015        self.failUnless( len(segs) == 2, "wrong number of segments after transformation")        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")        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")        self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9]) < self.EPS, "wrong vertex. 0  after transformation")
# Line 1211  class Test_FaultSystem(unittest.TestCase Line 1023  class Test_FaultSystem(unittest.TestCase
1023        self.failUnless( c.size == 2, "center size is wrong")        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.")        self.failUnless( numpy.linalg.norm(c-[7./5.,-7./5.]) < self.EPS, "center has wrong coordinates.")
1025        o=f.getOrientationOnSurface()/pi*180.        o=f.getOrientationOnSurface()/pi*180.
1026        self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")        self.failUnless( abs(o-45.) < self.EPS, "wrong orientation.")
1027    
1028        p=f.getParametrization([-1.,0.],1)        p=f.getParametrization([-1.,0.],1)
1029        self.failUnless(p[1]==1., "wrong value.")        self.failUnless(p[1]==1., "wrong value.")
# Line 1244  class Test_FaultSystem(unittest.TestCase Line 1056  class Test_FaultSystem(unittest.TestCase
1056        self.failUnless(p[1]==1., "wrong value.")        self.failUnless(p[1]==1., "wrong value.")
1057        self.failUnless(abs(p[0]-20.)<self.EPS, "wrong value.")        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_FaultSystem))     # suite.addTest(unittest.makeSuite(Test_FaultSystem))
1464     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))     # suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
    suite.addTest(unittest.makeSuite(Test_Darcy3D))  
    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 1260  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.2647  
changed lines
  Added in v.3527

  ViewVC Help
Powered by ViewVC 1.1.26