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

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

  ViewVC Help
Powered by ViewVC 1.1.26