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

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

  ViewVC Help
Powered by ViewVC 1.1.26