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

Legend:
Removed from v.2234  
changed lines
  Added in v.3527

  ViewVC Help
Powered by ViewVC 1.1.26