# Diff of /trunk/finley/test/python/run_models.py

revision 1562 by gross, Wed May 21 13:04:40 2008 UTC revision 3502 by gross, Thu Apr 28 05:06:24 2011 UTC
# Line 1  Line 1
1  #######################################################  # -*- coding: utf-8 -*-
2  #
3  #           Copyright 2003-2007 by ACceSS MNRF  ########################################################
#       Copyright 2007 by University of Queensland
4  #  #
5  #                http://esscc.uq.edu.au  # Copyright (c) 2003-2010 by University of Queensland
6  #        Primary Business: Queensland, Australia  # Earth Systems Science Computational Center (ESSCC)
7  #  Licensed under the Open Software License version 3.0  # http://www.uq.edu.au/esscc
8  #  #
9  #######################################################  # Primary Business: Queensland, Australia
12  #  #
13    ########################################################
14
16                      http://www.access.edu.au  Earth Systems Science Computational Center (ESSCC)
17                  Primary Business: Queensland, Australia"""  http://www.uq.edu.au/esscc
22
23  import unittest  import unittest
24  import tempfile  import tempfile
25
26
27
28    VERBOSE=False  and True
29
30  from esys.escript import *  from esys.escript import *
31  from esys.finley import Rectangle  from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian, FaultSystem, DarcyFlow
32    from esys.escript.models import Mountains
33    from esys.finley import Rectangle, Brick
34
35    from math import pi
36    import numpy
37  import sys  import sys
38  import os  import os
39    #====================================================================================================================
40  try:  try:
41       FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']       FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
42  except KeyError:  except KeyError:
43       FINLEY_WORKDIR='.'       FINLEY_WORKDIR='.'
44
45    #====================================================================================================================
46  NE=6  class Test_StokesProblemCartesian2D(unittest.TestCase):
TOL=1.e-5
VERBOSE=False or True

from esys.escript import *
from esys.escript.models import StokesProblemCartesian
from esys.finley import Rectangle, Brick
class Test_Simple2DModels(unittest.TestCase):
47     def setUp(self):     def setUp(self):
48         self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)         NE=6
49           self.TOL=1e-3
50           self.domain=Rectangle(NE,NE,order=-1)
51     def tearDown(self):     def tearDown(self):
52         del self.domain         del self.domain
53     def test_StokesProblemCartesian_PCG_P_0(self):     def test_PCG_P_0(self):
54         ETA=1.         ETA=1.
55         P1=0.         P1=0.
56
# Line 56  class Test_Simple2DModels(unittest.TestC Line 65  class Test_Simple2DModels(unittest.TestC
65
67         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
68         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
69         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL)
71
72         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
73         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
74         zz=P1*x[0]*x[1]-p         error_p=Lsup(p+P1*x[0]*x[1])
75         error_p=Lsup(zz-integrate(zz))         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
76         # print error_p, error_v0,error_v1         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
77         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
self.failUnless(error_v0<TOL,"0-velocity error too large.")
self.failUnless(error_v1<TOL,"1-velocity error too large.")
78
79     def test_StokesProblemCartesian_PCG_P_small(self):     def test_PCG_P_small(self):
80         ETA=1.         ETA=1.
81         P1=1.         P1=1.
82
# Line 83  class Test_Simple2DModels(unittest.TestC Line 91  class Test_Simple2DModels(unittest.TestC
91
93         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
94         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
95         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL)
96                 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
97         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
98         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
99         zz=P1*x[0]*x[1]-p         error_p=Lsup(P1*x[0]*x[1]+p)
100         error_p=Lsup(zz-integrate(zz))         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
101         # print error_p, error_v0,error_v1         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
102         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
self.failUnless(error_v0<TOL,"0-velocity error too large.")
self.failUnless(error_v1<TOL,"1-velocity error too large.")
103
104     def test_StokesProblemCartesian_PCG_P_large(self):     def test_PCG_P_large(self):
105         ETA=1.         ETA=1.
106         P1=1000.         P1=1000.
107
# Line 110  class Test_Simple2DModels(unittest.TestC Line 116  class Test_Simple2DModels(unittest.TestC
116
118         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
119         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
120         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL)
121           u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
122           # u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
123
124         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
125         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
126         zz=P1*x[0]*x[1]-p         error_p=Lsup(P1*x[0]*x[1]+p)/P1
127         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
128         # print error_p, error_v0,error_v1         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
129         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
self.failUnless(error_v0<TOL,"0-velocity error too large.")
self.failUnless(error_v1<TOL,"1-velocity error too large.")
130
131     def test_StokesProblemCartesian_GMRES_P_0(self):     def test_GMRES_P_0(self):
132         ETA=1.         ETA=1.
133         P1=0.         P1=0.
134
# Line 137  class Test_Simple2DModels(unittest.TestC Line 143  class Test_Simple2DModels(unittest.TestC
143
145         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
146         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
147         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=20)         sp.setTolerance(self.TOL)
148           u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=50,usePCG=False,iter_restart=18)
149
150         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
151         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
152         zz=P1*x[0]*x[1]-p         error_p=Lsup(P1*x[0]*x[1]+p)
153         error_p=Lsup(zz-integrate(zz))         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
154           self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
155        # print error_p, error_v0,error_v1         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
self.failUnless(error_p<TOL,"pressure error too large.")
self.failUnless(error_v0<TOL,"0-velocity error too large.")
self.failUnless(error_v1<TOL,"1-velocity error too large.")
156
157     def test_StokesProblemCartesian_GMRES_P_small(self):     def test_GMRES_P_small(self):
158         ETA=1.         ETA=1.
159         P1=1.         P1=1.
160
# Line 165  class Test_Simple2DModels(unittest.TestC Line 169  class Test_Simple2DModels(unittest.TestC
169
171         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
172         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
173         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")         sp.setTolerance(self.TOL)
174           u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=20,usePCG=False)
175
176         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
177         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
178         zz=P1*x[0]*x[1]-p         error_p=Lsup(P1*x[0]*x[1]+p)
error_p=Lsup(zz-integrate(zz))
# print error_p, error_v0,error_v1
self.failUnless(error_p<TOL,"pressure error too large.")
self.failUnless(error_v0<TOL,"0-velocity error too large.")
self.failUnless(error_v1<TOL,"1-velocity error too large.")
179
180     def test_StokesProblemCartesian_GMRES_P_large(self):         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
181         ETA=1.         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
182         P1=1000.         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
183
184         x=self.domain.getX()     def test_GMRES_P_large(self):
F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
+whereZero(x[0]-1)  * [1.,1.] \
+whereZero(x[1])    * [1.,0.] \
+whereZero(x[1]-1)  * [1.,1.]

sp=StokesProblemCartesian(self.domain)

u0=(1-x[0])*x[0]*[0.,1.]
p0=Scalar(P1,ReducedSolution(self.domain))
u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")

error_v0=Lsup(u[0]-u0[0])
error_v1=Lsup(u[1]-u0[1])/0.25
zz=P1*x[0]*x[1]-p
error_p=Lsup(zz-integrate(zz))/P1
# print error_p, error_v0,error_v1
self.failUnless(error_p<TOL,"pressure error too large.")
self.failUnless(error_v0<TOL,"0-velocity error too large.")
self.failUnless(error_v1<TOL,"1-velocity error too large.")

def test_StokesProblemCartesian_MINRES_P_0(self):
ETA=1.
P1=0.

x=self.domain.getX()
F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
+whereZero(x[0]-1)  * [1.,1.] \
+whereZero(x[1])    * [1.,0.] \
+whereZero(x[1]-1)  * [1.,1.]

sp=StokesProblemCartesian(self.domain)

u0=(1-x[0])*x[0]*[0.,1.]
p0=Scalar(P1,ReducedSolution(self.domain))
u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")

error_v0=Lsup(u[0]-u0[0])
error_v1=Lsup(u[1]-u0[1])/0.25
zz=P1*x[0]*x[1]-p
error_p=Lsup(zz-integrate(zz))
# print error_p, error_v0,error_v1
self.failUnless(error_p<TOL,"pressure error too large.")
self.failUnless(error_v0<TOL,"0-velocity error too large.")
self.failUnless(error_v1<TOL,"1-velocity error too large.")

def test_StokesProblemCartesian_MINRES_P_small(self):
ETA=1.
P1=1.

x=self.domain.getX()
F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
+whereZero(x[0]-1)  * [1.,1.] \
+whereZero(x[1])    * [1.,0.] \
+whereZero(x[1]-1)  * [1.,1.]

sp=StokesProblemCartesian(self.domain)

u0=(1-x[0])*x[0]*[0.,1.]
p0=Scalar(P1,ReducedSolution(self.domain))
u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")

error_v0=Lsup(u[0]-u0[0])
error_v1=Lsup(u[1]-u0[1])/0.25
zz=P1*x[0]*x[1]-p
error_p=Lsup(zz-integrate(zz))/P1
# print error_p, error_v0,error_v1
self.failUnless(error_p<TOL,"pressure error too large.")
self.failUnless(error_v0<TOL,"0-velocity error too large.")
self.failUnless(error_v1<TOL,"1-velocity error too large.")

#   def test_StokesProblemCartesian_MINRES_P_large(self):
#       ETA=1.
#       P1=1000.
#
#       x=self.domain.getX()
#       F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
#              +whereZero(x[0]-1)  * [1.,1.] \
#              +whereZero(x[1])    * [1.,0.] \
#              +whereZero(x[1]-1)  * [1.,1.]

#       sp=StokesProblemCartesian(self.domain)

#       u0=(1-x[0])*x[0]*[0.,1.]
#       p0=Scalar(P1,ReducedSolution(self.domain))
#       u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")

#       error_v0=Lsup(u[0]-u0[0])
#       error_v1=Lsup(u[1]-u0[1])/0.25
#       zz=P1*x[0]*x[1]-p
#       error_p=Lsup(zz-integrate(zz))/P1
# print error_p, error_v0,error_v1
#       self.failUnless(error_p<TOL,"pressure error too large.")
#       self.failUnless(error_v0<TOL,"0-velocity error too large.")
#       self.failUnless(error_v1<TOL,"1-velocity error too large.")

#   def test_StokesProblemCartesian_TFQMR_P_0(self):
#       ETA=1.
#       P1=0.

#       x=self.domain.getX()
#       F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
#              +whereZero(x[0]-1)  * [1.,1.] \
#              +whereZero(x[1])    * [1.,0.] \
#              +whereZero(x[1]-1)  * [1.,1.]

#       sp=StokesProblemCartesian(self.domain)

#       u0=(1-x[0])*x[0]*[0.,1.]
#       p0=Scalar(P1,ReducedSolution(self.domain))
#       u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")

#       error_v0=Lsup(u[0]-u0[0])
#       error_v1=Lsup(u[1]-u0[1])/0.25
#       zz=P1*x[0]*x[1]-p
#       error_p=Lsup(zz-integrate(zz))
# print error_p, error_v0,error_v1
#       self.failUnless(error_p<TOL,"pressure error too large.")
#       self.failUnless(error_v0<TOL,"0-velocity error too large.")
#       self.failUnless(error_v1<TOL,"1-velocity error too large.")

def test_StokesProblemCartesian_TFQMR_P_small(self):
ETA=1.
P1=1.

x=self.domain.getX()
F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
+whereZero(x[0]-1)  * [1.,1.] \
+whereZero(x[1])    * [1.,0.] \
+whereZero(x[1]-1)  * [1.,1.]

sp=StokesProblemCartesian(self.domain)

u0=(1-x[0])*x[0]*[0.,1.]
p0=Scalar(P1,ReducedSolution(self.domain))
u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")

error_v0=Lsup(u[0]-u0[0])
error_v1=Lsup(u[1]-u0[1])/0.25
zz=P1*x[0]*x[1]-p
error_p=Lsup(zz-integrate(zz))/P1
# print error_p, error_v0,error_v1
self.failUnless(error_p<TOL,"pressure error too large.")
self.failUnless(error_v0<TOL,"0-velocity error too large.")
self.failUnless(error_v1<TOL,"1-velocity error too large.")

def test_StokesProblemCartesian_TFQMR_P_large(self):
185         ETA=1.         ETA=1.
186         P1=1000.         P1=1000.
187
# Line 355  class Test_Simple2DModels(unittest.TestC Line 196  class Test_Simple2DModels(unittest.TestC
196
198         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
199         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
200         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")         sp.setTolerance(self.TOL)
201           u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)
202
203         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
204         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
205         zz=P1*x[0]*x[1]-p         error_p=Lsup(P1*x[0]*x[1]+p)/P1
206         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
207         # print error_p, error_v0,error_v1         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
208         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
209         self.failUnless(error_v0<TOL,"0-velocity error too large.")  #====================================================================================================================
210         self.failUnless(error_v1<TOL,"1-velocity error too large.")  class Test_StokesProblemCartesian3D(unittest.TestCase):

class Test_Simple3DModels(unittest.TestCase):
211     def setUp(self):     def setUp(self):
212         self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)         NE=6
213           self.TOL=1e-4
214           self.domain=Brick(NE,NE,NE,order=-1)
215     def tearDown(self):     def tearDown(self):
216         del self.domain         del self.domain
217     def test_StokesProblemCartesian_PCG_P_0(self):     def test_PCG_P_0(self):
218         ETA=1.         ETA=1.
219         P1=0.         P1=0.
220
# Line 383  class Test_Simple3DModels(unittest.TestC Line 223  class Test_Simple3DModels(unittest.TestC
223         x=self.domain.getX()         x=self.domain.getX()
225                +whereZero(x[0]-1)  * [1.,1.,1.] \                +whereZero(x[0]-1)  * [1.,1.,1.] \
226                +whereZero(x[1])    * [1.,1.,1.] \                +whereZero(x[1])    * [1.,0.,1.] \
227                +whereZero(x[1]-1)  * [1.,1.,1.] \                +whereZero(x[1]-1)  * [1.,1.,1.] \
228                +whereZero(x[2])    * [1.,1.,0.] \                +whereZero(x[2])    * [1.,1.,0.] \
229                +whereZero(x[2]-1)  * [1.,1.,1.]                +whereZero(x[2]-1)  * [1.,1.,1.]
# Line 393  class Test_Simple3DModels(unittest.TestC Line 233  class Test_Simple3DModels(unittest.TestC
233
235         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
236         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
237         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL)
238           u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
239
240         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
241         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
242         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
243         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
244         error_p=Lsup(zz-integrate(zz))         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
245         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
246         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
247         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
self.failUnless(error_v1<TOL,"1-velocity error too large.")
self.failUnless(error_v2<TOL,"2-velocity error too large.")
248
249     def test_StokesProblemCartesian_PCG_P_small(self):     def test_PCG_P_small(self):
250         ETA=1.         ETA=1.
251         P1=1.         P1=1.
252
# Line 415  class Test_Simple3DModels(unittest.TestC Line 254  class Test_Simple3DModels(unittest.TestC
254         F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]         F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
256                +whereZero(x[0]-1)  * [1.,1.,1.] \                +whereZero(x[0]-1)  * [1.,1.,1.] \
257                +whereZero(x[1])    * [1.,1.,1.] \                +whereZero(x[1])    * [1.,0.,1.] \
258                +whereZero(x[1]-1)  * [1.,1.,1.] \                +whereZero(x[1]-1)  * [1.,1.,1.] \
259                +whereZero(x[2])    * [1.,1.,0.] \                +whereZero(x[2])    * [1.,1.,0.] \
260                +whereZero(x[2]-1)  * [1.,1.,1.]                +whereZero(x[2]-1)  * [1.,1.,1.]
# Line 425  class Test_Simple3DModels(unittest.TestC Line 264  class Test_Simple3DModels(unittest.TestC
264
266         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
267         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
268         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL)
269                 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
270         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
271         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
272         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
273         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
274         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
275         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
276         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
277         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
278         self.failUnless(error_v1<TOL,"1-velocity error too large.")
279         self.failUnless(error_v2<TOL,"2-velocity error too large.")     def test_PCG_P_large(self):
def test_StokesProblemCartesian_PCG_P_large(self):
280         ETA=1.         ETA=1.
281         P1=1000.         P1=1000.
282
# Line 446  class Test_Simple3DModels(unittest.TestC Line 284  class Test_Simple3DModels(unittest.TestC
284         F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]         F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
286                +whereZero(x[0]-1)  * [1.,1.,1.] \                +whereZero(x[0]-1)  * [1.,1.,1.] \
287                +whereZero(x[1])    * [1.,1.,1.] \                +whereZero(x[1])    * [1.,0.,1.] \
288                +whereZero(x[1]-1)  * [1.,1.,1.] \                +whereZero(x[1]-1)  * [1.,1.,1.] \
289                +whereZero(x[2])    * [1.,1.,0.] \                +whereZero(x[2])    * [1.,1.,0.] \
290                +whereZero(x[2]-1)  * [1.,1.,1.]                +whereZero(x[2]-1)  * [1.,1.,1.]
# Line 456  class Test_Simple3DModels(unittest.TestC Line 294  class Test_Simple3DModels(unittest.TestC
294
296         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
297         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
298         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")         sp.setTolerance(self.TOL)
299           u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
300
301         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
302         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
303         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
304         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
305         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
306         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
307         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
308         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
self.failUnless(error_v1<TOL,"1-velocity error too large.")
self.failUnless(error_v2<TOL,"2-velocity error too large.")
309
310     def test_StokesProblemCartesian_GMRES_P_0(self):     def test_GMRES_P_0(self):
311         ETA=1.         ETA=1.
312         P1=0.         P1=0.
313
# Line 489  class Test_Simple3DModels(unittest.TestC Line 326  class Test_Simple3DModels(unittest.TestC
326
328         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
329         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
330         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=20)         sp.setTolerance(self.TOL)
331           u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False,iter_restart=20)
332
333         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
334         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
335         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
336         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
337         error_p=Lsup(zz-integrate(zz))         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
338         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
339         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
340         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
341         self.failUnless(error_v1<TOL,"1-velocity error too large.")     def test_GMRES_P_small(self):
self.failUnless(error_v2<TOL,"2-velocity error too large.")
def test_StokesProblemCartesian_GMRES_P_small(self):
342         ETA=1.         ETA=1.
343         P1=1.         P1=1.
344
# Line 520  class Test_Simple3DModels(unittest.TestC Line 356  class Test_Simple3DModels(unittest.TestC
356
358         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
359         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
360         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")         sp.setTolerance(self.TOL/10)
361           u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)
362
363         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
364         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
365         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
366         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
367         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
368         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
369         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
370         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
371         self.failUnless(error_v1<TOL,"1-velocity error too large.")     def test_GMRES_P_large(self):
self.failUnless(error_v2<TOL,"2-velocity error too large.")
def test_StokesProblemCartesian_GMRES_P_large(self):
372         ETA=1.         ETA=1.
373         P1=1000.         P1=1000.
374
# Line 541  class Test_Simple3DModels(unittest.TestC Line 376  class Test_Simple3DModels(unittest.TestC
376         F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]         F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
378                +whereZero(x[0]-1)  * [1.,1.,1.] \                +whereZero(x[0]-1)  * [1.,1.,1.] \
379                +whereZero(x[1])    * [1.,1.,1.] \                +whereZero(x[1])    * [1.,0.,1.] \
+whereZero(x[1]-1)  * [1.,1.,1.] \
+whereZero(x[2])    * [1.,1.,0.] \
+whereZero(x[2]-1)  * [1.,1.,1.]

sp=StokesProblemCartesian(self.domain)

u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
p0=Scalar(P1,ReducedSolution(self.domain))
u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")

error_v0=Lsup(u[0]-u0[0])
error_v1=Lsup(u[1]-u0[1])
error_v2=Lsup(u[2]-u0[2])/0.25**2
zz=P1*x[0]*x[1]*x[2]-p
error_p=Lsup(zz-integrate(zz))/P1
# print error_p, error_v0,error_v1,error_v2
self.failUnless(error_p<TOL,"pressure error too large.")
self.failUnless(error_v0<TOL,"0-velocity error too large.")
self.failUnless(error_v1<TOL,"1-velocity error too large.")
self.failUnless(error_v2<TOL,"2-velocity error too large.")

def test_StokesProblemCartesian_MINRES_P_0(self):
ETA=1.
P1=0.

x=self.domain.getX()
F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
x=self.domain.getX()
+whereZero(x[0]-1)  * [1.,1.,1.] \
+whereZero(x[1])    * [1.,1.,1.] \
+whereZero(x[1]-1)  * [1.,1.,1.] \
+whereZero(x[2])    * [1.,1.,0.] \
+whereZero(x[2]-1)  * [1.,1.,1.]

sp=StokesProblemCartesian(self.domain)

u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
p0=Scalar(P1,ReducedSolution(self.domain))
u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")

error_v0=Lsup(u[0]-u0[0])
error_v1=Lsup(u[1]-u0[1])
error_v2=Lsup(u[2]-u0[2])/0.25**2
zz=P1*x[0]*x[1]*x[2]-p
error_p=Lsup(zz-integrate(zz))
# print error_p, error_v0,error_v1,error_v2
self.failUnless(error_p<TOL,"pressure error too large.")
self.failUnless(error_v0<TOL,"0-velocity error too large.")
self.failUnless(error_v1<TOL,"1-velocity error too large.")
self.failUnless(error_v2<TOL,"2-velocity error too large.")

def test_StokesProblemCartesian_MINRES_P_small(self):
ETA=1.
P1=1.

x=self.domain.getX()
F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
+whereZero(x[0]-1)  * [1.,1.,1.] \
+whereZero(x[1])    * [1.,1.,1.] \
380                +whereZero(x[1]-1)  * [1.,1.,1.] \                +whereZero(x[1]-1)  * [1.,1.,1.] \
381                +whereZero(x[2])    * [1.,1.,0.] \                +whereZero(x[2])    * [1.,1.,0.] \
382                +whereZero(x[2]-1)  * [1.,1.,1.]                +whereZero(x[2]-1)  * [1.,1.,1.]
# Line 616  class Test_Simple3DModels(unittest.TestC Line 386  class Test_Simple3DModels(unittest.TestC
386
388         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
389         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(-P1,ReducedSolution(self.domain))
390         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")         sp.setTolerance(self.TOL)
391           u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=False)
392
393         error_v0=Lsup(u[0]-u0[0])         error_v0=Lsup(u[0]-u0[0])
394         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
395         error_v2=Lsup(u[2]-u0[2])/0.25**2         error_v2=Lsup(u[2]-u0[2])/0.25**2
396         zz=P1*x[0]*x[1]*x[2]-p         error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
397         error_p=Lsup(zz-integrate(zz))/P1         self.failUnless(error_p<10*self.TOL, "pressure error too large.")
398         # print error_p, error_v0,error_v1,error_v2         self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
399         self.failUnless(error_p<TOL,"pressure error too large.")         self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
400         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
401         self.failUnless(error_v1<TOL,"1-velocity error too large.")  #====================================================================================================================
self.failUnless(error_v2<TOL,"2-velocity error too large.")

#   def test_StokesProblemCartesian_MINRES_P_large(self):
#       ETA=1.
#       P1=1000.

#       x=self.domain.getX()
#       F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
#              +whereZero(x[0]-1)  * [1.,1.,1.] \
#              +whereZero(x[1])    * [1.,1.,1.] \
#              +whereZero(x[1]-1)  * [1.,1.,1.] \
#              +whereZero(x[2])    * [1.,1.,0.] \
#              +whereZero(x[2]-1)  * [1.,1.,1.]

#       sp=StokesProblemCartesian(self.domain)

#       u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
#       p0=Scalar(P1,ReducedSolution(self.domain))
#       u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")

#       error_v0=Lsup(u[0]-u0[0])
#       error_v1=Lsup(u[1]-u0[1])
#       error_v2=Lsup(u[2]-u0[2])/0.25**2
#       zz=P1*x[0]*x[1]*x[2]-p
#       error_p=Lsup(zz-integrate(zz))/P1
# print error_p, error_v0,error_v1,error_v2
#       self.failUnless(error_p<TOL,"pressure error too large.")
#       self.failUnless(error_v0<TOL,"0-velocity error too large.")
#       self.failUnless(error_v1<TOL,"1-velocity error too large.")
#       self.failUnless(error_v2<TOL,"2-velocity error too large.")

#   def test_StokesProblemCartesian_TFQMR_P_0(self):
#       ETA=1.
#       P1=0.

#       x=self.domain.getX()
#       F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
#       x=self.domain.getX()
#              +whereZero(x[0]-1)  * [1.,1.,1.] \
#              +whereZero(x[1])    * [1.,1.,1.] \
#              +whereZero(x[1]-1)  * [1.,1.,1.] \
#              +whereZero(x[2])    * [1.,1.,0.] \
#              +whereZero(x[2]-1)  * [1.,1.,1.]

#       sp=StokesProblemCartesian(self.domain)

#       u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
#       p0=Scalar(P1,ReducedSolution(self.domain))
#       u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")

#       error_v0=Lsup(u[0]-u0[0])
#       error_v1=Lsup(u[1]-u0[1])
#       error_v2=Lsup(u[2]-u0[2])/0.25**2
#       zz=P1*x[0]*x[1]*x[2]-p
#       error_p=Lsup(zz-integrate(zz))
# print error_p, error_v0,error_v1,error_v2
#       self.failUnless(error_p<TOL,"pressure error too large.")
#       self.failUnless(error_v0<TOL,"0-velocity error too large.")
#       self.failUnless(error_v1<TOL,"1-velocity error too large.")
#       self.failUnless(error_v2<TOL,"2-velocity error too large.")
402
403     def test_StokesProblemCartesian_TFQMR_P_small(self):  class Test_Mountains3D(unittest.TestCase):
404         ETA=1.     def setUp(self):
405         P1=1.         NE=16
406           self.TOL=1e-4
407           self.domain=Brick(NE,NE,NE,order=1,useFullElementOrder=True)
408       def tearDown(self):
409           del self.domain
410       def test_periodic(self):
411           EPS=0.01
412
413         x=self.domain.getX()         x=self.domain.getX()
414         F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]         v = Vector(0.0, Solution(self.domain))
415         mask=whereZero(x[0])    * [1.,1.,1.] \         a0=1
416                +whereZero(x[0]-1)  * [1.,1.,1.] \         a1=1
417                +whereZero(x[1])    * [1.,1.,1.] \         n0=2
418                +whereZero(x[1]-1)  * [1.,1.,1.] \         n1=2
419                +whereZero(x[2])    * [1.,1.,0.] \         n2=0.5
420                +whereZero(x[2]-1)  * [1.,1.,1.]         a2=-(a0*n0+a1*n1)/n2
421           v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])* cos(pi*n2*x[2])
422           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])
424
425           mts=Mountains(self.domain,eps=EPS)
426           mts.setVelocity(v)
427           Z=mts.update()
428
429                 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[2]-1.)))
430         sp=StokesProblemCartesian(self.domain)         self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")

u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
p0=Scalar(P1,ReducedSolution(self.domain))
u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")

error_v0=Lsup(u[0]-u0[0])
error_v1=Lsup(u[1]-u0[1])
error_v2=Lsup(u[2]-u0[2])/0.25**2
zz=P1*x[0]*x[1]*x[2]-p
error_p=Lsup(zz-integrate(zz))/P1
# print error_p, error_v0,error_v1,error_v2
self.failUnless(error_p<TOL,"pressure error too large.")
self.failUnless(error_v0<TOL,"0-velocity error too large.")
self.failUnless(error_v1<TOL,"1-velocity error too large.")
self.failUnless(error_v2<TOL,"2-velocity error too large.")
431
432     def test_StokesProblemCartesian_TFQMR_P_large(self):  class Test_Mountains2D(unittest.TestCase):
433         ETA=1.     def setUp(self):
434         P1=1000.         NE=16
435           self.TOL=1e-4
436           self.domain=Rectangle(NE,NE,order=1,useFullElementOrder=True)
437       def tearDown(self):
438           del self.domain
439       def test_periodic(self):
440           EPS=0.01
441
442         x=self.domain.getX()         x=self.domain.getX()
443         F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]         v = Vector(0.0, Solution(self.domain))
444         mask=whereZero(x[0])    * [1.,1.,1.] \         a0=1
445                +whereZero(x[0]-1)  * [1.,1.,1.] \         n0=1
446                +whereZero(x[1])    * [1.,1.,1.] \         n1=0.5
447                +whereZero(x[1]-1)  * [1.,1.,1.] \         a1=-(a0*n0)/n1
448                +whereZero(x[2])    * [1.,1.,0.] \         v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])
449                +whereZero(x[2]-1)  * [1.,1.,1.]         v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])
450
451                 H_t=Scalar(0.0, Solution(self.domain))
452         sp=StokesProblemCartesian(self.domain)         mts=Mountains(self.domain,eps=EPS)
453                 mts.setVelocity(v)
455         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
456         p0=Scalar(P1,ReducedSolution(self.domain))         error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[1]-1.)))
457         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")         self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
458
459         error_v0=Lsup(u[0]-u0[0])
460         error_v1=Lsup(u[1]-u0[1])
461         error_v2=Lsup(u[2]-u0[2])/0.25**2  class Test_Rheologies(unittest.TestCase):
462         zz=P1*x[0]*x[1]*x[2]-p       """
463         error_p=Lsup(zz-integrate(zz))/P1       this is the program used to generate the powerlaw tests:
464         # print error_p, error_v0,error_v1,error_v2
465         self.failUnless(error_p<TOL,"pressure error too large.")       TAU_Y=100.
466         self.failUnless(error_v0<TOL,"0-velocity error too large.")       N=10
467         self.failUnless(error_v1<TOL,"1-velocity error too large.")       M=5
468         self.failUnless(error_v2<TOL,"2-velocity error too large.")
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")
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
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
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
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
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
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
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]
702          if free==None:
703             for d in range(self.dom.getDim()):
705          else:
706             for d in range(self.dom.getDim()):
707                if d == self.dom.getDim()-1:
709                else:
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)
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)
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.] ]
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.] ]
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.] ]
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.] ]
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.] ]
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.] ]
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] ]
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          #
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.] ]
1349          top1=[ [10.,0.,0], [10.,10.,0], [0.,10.,0] ]
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()