/[escript]/trunk/finley/test/python/run_models.py
ViewVC logotype

Contents of /trunk/finley/test/python/run_models.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6657 - (show annotations)
Thu Mar 29 04:43:17 2018 UTC (17 months, 3 weeks ago) by uqagarro
File MIME type: text/x-python
File size: 79161 byte(s)
Added test for IncompressibleIsotropicFlowCartesianOnFinley
with higher NE


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

  ViewVC Help
Powered by ViewVC 1.1.26