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

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

  ViewVC Help
Powered by ViewVC 1.1.26