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

Legend:
Removed from v.1562  
changed lines
  Added in v.2960

  ViewVC Help
Powered by ViewVC 1.1.26