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

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

  ViewVC Help
Powered by ViewVC 1.1.26