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

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

  ViewVC Help
Powered by ViewVC 1.1.26