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

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

  ViewVC Help
Powered by ViewVC 1.1.26