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

Revision 2675 - (show annotations)
Mon Sep 21 02:32:17 2009 UTC (12 years ago) by gross
File MIME type: text/x-python
File size: 85292 byte(s)
```getMinValue added to the fault system.
```
 1 2 ######################################################## 3 # 4 # Copyright (c) 2003-2009 by University of Queensland 5 # Earth Systems Science Computational Center (ESSCC) 6 7 # 8 # Primary Business: Queensland, Australia 9 # Licensed under the Open Software License version 3.0 10 11 # 12 ######################################################## 13 14 __copyright__="""Copyright (c) 2003-2009 by University of Queensland 15 Earth Systems Science Computational Center (ESSCC) 16 http://www.uq.edu.au/esscc 17 Primary Business: Queensland, Australia""" 18 __license__="""Licensed under the Open Software License version 3.0 19 20 __url__= 21 22 import unittest 23 import tempfile 24 25 26 27 VERBOSE=False # or True 28 29 from esys.escript import * 30 from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian, FaultSystem, DarcyFlow 31 from esys.escript.models import Mountains 32 from esys.finley import Rectangle, Brick 33 34 from math import pi 35 import numpy 36 import sys 37 import os 38 #==================================================================================================================== 39 try: 40 FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR'] 41 except KeyError: 42 FINLEY_WORKDIR='.' 43 44 #==================================================================================================================== 45 class Test_StokesProblemCartesian2D(unittest.TestCase): 46 def setUp(self): 47 NE=6 48 self.TOL=1e-3 49 self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True) 50 def tearDown(self): 51 del self.domain 52 def test_PCG_P_0(self): 53 ETA=1. 54 P1=0. 55 56 x=self.domain.getX() 57 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.] 58 mask=whereZero(x[0]) * [1.,1.] \ 59 +whereZero(x[0]-1) * [1.,1.] \ 60 +whereZero(x[1]) * [1.,0.] \ 61 +whereZero(x[1]-1) * [1.,1.] 62 63 sp=StokesProblemCartesian(self.domain) 64 65 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 66 u0=(1-x[0])*x[0]*[0.,1.] 67 p0=Scalar(-P1,ReducedSolution(self.domain)) 68 sp.setTolerance(self.TOL) 69 u,p=sp.solve(u0,p0,verbose=VERBOSE,max_iter=100,usePCG=True) 70 71 error_v0=Lsup(u[0]-u0[0]) 72 error_v1=Lsup(u[1]-u0[1])/0.25 73 error_p=Lsup(p+P1*x[0]*x[1]) 74 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 75 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 76 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 77 78 def test_PCG_P_small(self): 79 ETA=1. 80 P1=1. 81 82 x=self.domain.getX() 83 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.] 84 mask=whereZero(x[0]) * [1.,1.] \ 85 +whereZero(x[0]-1) * [1.,1.] \ 86 +whereZero(x[1]) * [1.,0.] \ 87 +whereZero(x[1]-1) * [1.,1.] 88 89 sp=StokesProblemCartesian(self.domain) 90 91 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 92 u0=(1-x[0])*x[0]*[0.,1.] 93 p0=Scalar(-P1,ReducedSolution(self.domain)) 94 sp.setTolerance(self.TOL*0.2) 95 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True) 96 error_v0=Lsup(u[0]-u0[0]) 97 error_v1=Lsup(u[1]-u0[1])/0.25 98 error_p=Lsup(P1*x[0]*x[1]+p) 99 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 100 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 101 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 102 103 def test_PCG_P_large(self): 104 ETA=1. 105 P1=1000. 106 107 x=self.domain.getX() 108 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.] 109 mask=whereZero(x[0]) * [1.,1.] \ 110 +whereZero(x[0]-1) * [1.,1.] \ 111 +whereZero(x[1]) * [1.,0.] \ 112 +whereZero(x[1]-1) * [1.,1.] 113 114 sp=StokesProblemCartesian(self.domain) 115 116 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 117 u0=(1-x[0])*x[0]*[0.,1.] 118 p0=Scalar(-P1,ReducedSolution(self.domain)) 119 sp.setTolerance(self.TOL) 120 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True) 121 122 error_v0=Lsup(u[0]-u0[0]) 123 error_v1=Lsup(u[1]-u0[1])/0.25 124 error_p=Lsup(P1*x[0]*x[1]+p)/P1 125 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 126 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 127 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 128 129 def test_GMRES_P_0(self): 130 ETA=1. 131 P1=0. 132 133 x=self.domain.getX() 134 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.] 135 mask=whereZero(x[0]) * [1.,1.] \ 136 +whereZero(x[0]-1) * [1.,1.] \ 137 +whereZero(x[1]) * [1.,0.] \ 138 +whereZero(x[1]-1) * [1.,1.] 139 140 sp=StokesProblemCartesian(self.domain) 141 142 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 143 u0=(1-x[0])*x[0]*[0.,1.] 144 p0=Scalar(-P1,ReducedSolution(self.domain)) 145 sp.setTolerance(self.TOL) 146 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=50,usePCG=False,iter_restart=18) 147 148 error_v0=Lsup(u[0]-u0[0]) 149 error_v1=Lsup(u[1]-u0[1])/0.25 150 error_p=Lsup(P1*x[0]*x[1]+p) 151 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 152 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 153 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 154 155 def test_GMRES_P_small(self): 156 ETA=1. 157 P1=1. 158 159 x=self.domain.getX() 160 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.] 161 mask=whereZero(x[0]) * [1.,1.] \ 162 +whereZero(x[0]-1) * [1.,1.] \ 163 +whereZero(x[1]) * [1.,0.] \ 164 +whereZero(x[1]-1) * [1.,1.] 165 166 sp=StokesProblemCartesian(self.domain) 167 168 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 169 u0=(1-x[0])*x[0]*[0.,1.] 170 p0=Scalar(-P1,ReducedSolution(self.domain)) 171 sp.setTolerance(self.TOL*0.1) 172 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=20,usePCG=False) 173 174 error_v0=Lsup(u[0]-u0[0]) 175 error_v1=Lsup(u[1]-u0[1])/0.25 176 error_p=Lsup(P1*x[0]*x[1]+p) 177 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 178 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 179 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 180 181 def test_GMRES_P_large(self): 182 ETA=1. 183 P1=1000. 184 185 x=self.domain.getX() 186 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.] 187 mask=whereZero(x[0]) * [1.,1.] \ 188 +whereZero(x[0]-1) * [1.,1.] \ 189 +whereZero(x[1]) * [1.,0.] \ 190 +whereZero(x[1]-1) * [1.,1.] 191 192 sp=StokesProblemCartesian(self.domain) 193 194 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 195 u0=(1-x[0])*x[0]*[0.,1.] 196 p0=Scalar(-P1,ReducedSolution(self.domain)) 197 sp.setTolerance(self.TOL) 198 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False) 199 200 error_v0=Lsup(u[0]-u0[0]) 201 error_v1=Lsup(u[1]-u0[1])/0.25 202 error_p=Lsup(P1*x[0]*x[1]+p)/P1 203 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 204 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 205 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 206 #==================================================================================================================== 207 class Test_StokesProblemCartesian3D(unittest.TestCase): 208 def setUp(self): 209 NE=6 210 self.TOL=1e-4 211 self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True) 212 def tearDown(self): 213 del self.domain 214 def test_PCG_P_0(self): 215 ETA=1. 216 P1=0. 217 218 x=self.domain.getX() 219 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.] 220 x=self.domain.getX() 221 mask=whereZero(x[0]) * [1.,1.,1.] \ 222 +whereZero(x[0]-1) * [1.,1.,1.] \ 223 +whereZero(x[1]) * [1.,0.,1.] \ 224 +whereZero(x[1]-1) * [1.,1.,1.] \ 225 +whereZero(x[2]) * [1.,1.,0.] \ 226 +whereZero(x[2]-1) * [1.,1.,1.] 227 228 229 sp=StokesProblemCartesian(self.domain) 230 231 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 232 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.] 233 p0=Scalar(-P1,ReducedSolution(self.domain)) 234 sp.setTolerance(self.TOL) 235 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True) 236 237 error_v0=Lsup(u[0]-u0[0]) 238 error_v1=Lsup(u[1]-u0[1]) 239 error_v2=Lsup(u[2]-u0[2])/0.25**2 240 error_p=Lsup(P1*x[0]*x[1]*x[2]+p) 241 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 242 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 243 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 244 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.") 245 246 def test_PCG_P_small(self): 247 ETA=1. 248 P1=1. 249 250 x=self.domain.getX() 251 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.] 252 mask=whereZero(x[0]) * [1.,1.,1.] \ 253 +whereZero(x[0]-1) * [1.,1.,1.] \ 254 +whereZero(x[1]) * [1.,0.,1.] \ 255 +whereZero(x[1]-1) * [1.,1.,1.] \ 256 +whereZero(x[2]) * [1.,1.,0.] \ 257 +whereZero(x[2]-1) * [1.,1.,1.] 258 259 260 sp=StokesProblemCartesian(self.domain) 261 262 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 263 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.] 264 p0=Scalar(-P1,ReducedSolution(self.domain)) 265 sp.setTolerance(self.TOL*0.1) 266 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True) 267 error_v0=Lsup(u[0]-u0[0]) 268 error_v1=Lsup(u[1]-u0[1]) 269 error_v2=Lsup(u[2]-u0[2])/0.25**2 270 error_p=Lsup(P1*x[0]*x[1]*x[2]+p) 271 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 272 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 273 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 274 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.") 275 276 def test_PCG_P_large(self): 277 ETA=1. 278 P1=1000. 279 280 x=self.domain.getX() 281 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.] 282 mask=whereZero(x[0]) * [1.,1.,1.] \ 283 +whereZero(x[0]-1) * [1.,1.,1.] \ 284 +whereZero(x[1]) * [1.,0.,1.] \ 285 +whereZero(x[1]-1) * [1.,1.,1.] \ 286 +whereZero(x[2]) * [1.,1.,0.] \ 287 +whereZero(x[2]-1) * [1.,1.,1.] 288 289 290 sp=StokesProblemCartesian(self.domain) 291 292 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 293 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.] 294 p0=Scalar(-P1,ReducedSolution(self.domain)) 295 sp.setTolerance(self.TOL) 296 u,p=sp.solve(u0,-p0, verbose=VERBOSE ,max_iter=100,usePCG=True) 297 298 error_v0=Lsup(u[0]-u0[0]) 299 error_v1=Lsup(u[1]-u0[1]) 300 error_v2=Lsup(u[2]-u0[2])/0.25**2 301 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1 302 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 303 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 304 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 305 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.") 306 307 def test_GMRES_P_0(self): 308 ETA=1. 309 P1=0. 310 311 x=self.domain.getX() 312 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.] 313 x=self.domain.getX() 314 mask=whereZero(x[0]) * [1.,1.,1.] \ 315 +whereZero(x[0]-1) * [1.,1.,1.] \ 316 +whereZero(x[1]) * [1.,1.,1.] \ 317 +whereZero(x[1]-1) * [1.,1.,1.] \ 318 +whereZero(x[2]) * [1.,1.,0.] \ 319 +whereZero(x[2]-1) * [1.,1.,1.] 320 321 322 sp=StokesProblemCartesian(self.domain) 323 324 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 325 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.] 326 p0=Scalar(-P1,ReducedSolution(self.domain)) 327 sp.setTolerance(self.TOL) 328 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False,iter_restart=20) 329 330 error_v0=Lsup(u[0]-u0[0]) 331 error_v1=Lsup(u[1]-u0[1]) 332 error_v2=Lsup(u[2]-u0[2])/0.25**2 333 error_p=Lsup(P1*x[0]*x[1]*x[2]+p) 334 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 335 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 336 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 337 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.") 338 def test_GMRES_P_small(self): 339 ETA=1. 340 P1=1. 341 342 x=self.domain.getX() 343 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.] 344 mask=whereZero(x[0]) * [1.,1.,1.] \ 345 +whereZero(x[0]-1) * [1.,1.,1.] \ 346 +whereZero(x[1]) * [1.,1.,1.] \ 347 +whereZero(x[1]-1) * [1.,1.,1.] \ 348 +whereZero(x[2]) * [1.,1.,0.] \ 349 +whereZero(x[2]-1) * [1.,1.,1.] 350 351 352 sp=StokesProblemCartesian(self.domain) 353 354 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 355 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.] 356 p0=Scalar(-P1,ReducedSolution(self.domain)) 357 sp.setTolerance(self.TOL*0.1) 358 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False) 359 360 error_v0=Lsup(u[0]-u0[0]) 361 error_v1=Lsup(u[1]-u0[1]) 362 error_v2=Lsup(u[2]-u0[2])/0.25**2 363 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1 364 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 365 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 366 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.") 367 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 368 def test_GMRES_P_large(self): 369 ETA=1. 370 P1=1000. 371 372 x=self.domain.getX() 373 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.] 374 mask=whereZero(x[0]) * [1.,1.,1.] \ 375 +whereZero(x[0]-1) * [1.,1.,1.] \ 376 +whereZero(x[1]) * [1.,0.,1.] \ 377 +whereZero(x[1]-1) * [1.,1.,1.] \ 378 +whereZero(x[2]) * [1.,1.,0.] \ 379 +whereZero(x[2]-1) * [1.,1.,1.] 380 381 382 sp=StokesProblemCartesian(self.domain) 383 384 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 385 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.] 386 p0=Scalar(-P1,ReducedSolution(self.domain)) 387 sp.setTolerance(self.TOL) 388 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=False) 389 390 error_v0=Lsup(u[0]-u0[0]) 391 error_v1=Lsup(u[1]-u0[1]) 392 error_v2=Lsup(u[2]-u0[2])/0.25**2 393 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1 394 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 395 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 396 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 397 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.") 398 #==================================================================================================================== 399 class Test_Darcy(unittest.TestCase): 400 # this is a simple test for the darcy flux problem 401 # 402 # 403 # p = 1/k * ( 1/2* (fb-f0)/xb* x **2 + f0 * x - ub*x ) + p0 404 # 405 # with f = (fb-f0)/xb* x + f0 406 # 407 # u = f - k * p,x = ub 408 # 409 # we prescribe pressure at x=x0=0 to p0 410 # 411 # if we prescribe the pressure on the bottom x=xb we set 412 # 413 # pb= 1/k * ( 1/2* (fb-f0)* xb + f0 * xb - ub*xb ) + p0 = 1/k * ((fb+f0)/2 - ub ) * xb + p0 414 # 415 # which leads to ub = (fb+f0)/2-k*(pb-p0)/xb 416 # 417 def rescaleDomain(self): 418 x=self.dom.getX().copy() 419 for i in xrange(self.dom.getDim()): 420 x_inf=inf(x[i]) 421 x_sup=sup(x[i]) 422 if i == self.dom.getDim()-1: 423 x[i]=-self.WIDTH*(x[i]-x_sup)/(x_inf-x_sup) 424 else: 425 x[i]=self.WIDTH*(x[i]-x_inf)/(x_sup-x_inf) 426 self.dom.setX(x) 427 def getScalarMask(self,include_bottom=True): 428 x=self.dom.getX().copy() 429 x_inf=inf(x[self.dom.getDim()-1]) 430 x_sup=sup(x[self.dom.getDim()-1]) 431 out=whereZero(x[self.dom.getDim()-1]-x_sup) 432 if include_bottom: out+=whereZero(x[self.dom.getDim()-1]-x_inf) 433 return wherePositive(out) 434 def getVectorMask(self,include_bottom=True): 435 x=self.dom.getX().copy() 436 out=Vector(0.,Solution(self.dom)) 437 for i in xrange(self.dom.getDim()): 438 x_inf=inf(x[i]) 439 x_sup=sup(x[i]) 440 if i != self.dom.getDim()-1: out[i]+=whereZero(x[i]-x_sup) 441 if i != self.dom.getDim()-1 or include_bottom: out[i]+=whereZero(x[i]-x_inf) 442 return wherePositive(out) 443 444 def setSolutionFixedBottom(self, p0, pb, f0, fb, k): 445 d=self.dom.getDim() 446 x=self.dom.getX()[d-1] 447 xb=inf(x) 448 u=Vector(0.,Solution(self.dom))+kronecker(d)[d-1]*((f0+fb)/2.-k*(pb-p0)/xb) 449 p=1./k*((fb-f0)/(xb*2.)* x**2 - (fb-f0)/2.*x)+(pb-p0)/xb*x + p0 450 f= ((fb-f0)/xb* x + f0)*kronecker(Function(self.dom))[d-1] 451 return u,p,f 452 453 def testConstF_FixedBottom_smallK(self): 454 k=1.e-10 455 mp=self.getScalarMask(include_bottom=True) 456 mv=self.getVectorMask(include_bottom=False) 457 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k) 458 p=p_ref*mp 459 u=u_ref*mv 460 df=DarcyFlow(self.dom) 461 df.setValue(g=f, 462 location_of_fixed_pressure=mp, 463 location_of_fixed_flux=mv, 464 permeability=Scalar(k,Function(self.dom))) 465 df.setTolerance(rtol=self.TOL) 466 v,p=df.solve(u_ref,p, max_iter=100, verbose=VERBOSE) 467 self.failUnless(Lsup(v-u_ref)=0,"eta needs to be positive (test %s)"%id) 816 error=abs(gamma_dot_*eta-tau_ref) 817 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)) 818 819 def test_PowerLaw_Linear(self): 820 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] 821 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] 822 pl=PowerLaw(numMaterials=1,verbose=VERBOSE) 823 pl.setDruckerPragerLaw(tau_Y=100.) 824 pl.setPowerLaw(eta_N=2.) 825 pl.setEtaTolerance(self.TOL) 826 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 827 828 def test_PowerLaw_QuadLarge(self): 829 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] 830 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] 831 pl=PowerLaw(numMaterials=2,verbose=VERBOSE) 832 pl.setDruckerPragerLaw(tau_Y=100.) 833 pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2]) 834 pl.setEtaTolerance(self.TOL) 835 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 836 837 def test_PowerLaw_QuadSmall(self): 838 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] 839 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] 840 pl=PowerLaw(numMaterials=2,verbose=VERBOSE) 841 pl.setDruckerPragerLaw(tau_Y=100.) 842 pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2]) 843 pl.setEtaTolerance(self.TOL) 844 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 845 846 def test_PowerLaw_CubeLarge(self): 847 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] 848 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] 849 pl=PowerLaw(numMaterials=2,verbose=VERBOSE) 850 pl.setDruckerPragerLaw(tau_Y=100.) 851 pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3]) 852 pl.setEtaTolerance(self.TOL) 853 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 854 855 def test_PowerLaw_CubeSmall(self): 856 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] 857 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] 858 pl=PowerLaw(numMaterials=2,verbose=VERBOSE) 859 pl.setDruckerPragerLaw(tau_Y=100.) 860 pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3]) 861 pl.setEtaTolerance(self.TOL) 862 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 863 864 def test_PowerLaw_QuadLarge_CubeLarge(self): 865 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] 866 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] 867 868 pl=PowerLaw(numMaterials=3,verbose=VERBOSE) 869 pl.setDruckerPragerLaw(tau_Y=100.) 870 pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3]) 871 pl.setEtaTolerance(self.TOL) 872 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 873 874 def test_PowerLaw_QuadLarge_CubeSmall(self): 875 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] 876 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] 877 878 pl=PowerLaw(numMaterials=3,verbose=VERBOSE) 879 pl.setDruckerPragerLaw(tau_Y=100.) 880 pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3]) 881 pl.setEtaTolerance(self.TOL) 882 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 883 884 def test_PowerLaw_QuadSmall_CubeLarge(self): 885 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] 886 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] 887 pl=PowerLaw(numMaterials=3,verbose=VERBOSE) 888 pl.setDruckerPragerLaw(tau_Y=100.) 889 pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3]) 890 pl.setEtaTolerance(self.TOL) 891 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 892 893 def test_PowerLaw_QuadSmall_CubeSmall(self): 894 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] 895 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] 896 pl=PowerLaw(numMaterials=3,verbose=VERBOSE) 897 pl.setDruckerPragerLaw(tau_Y=100.) 898 pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3]) 899 pl.setEtaTolerance(self.TOL) 900 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 901 902 def test_PowerLaw_withShear(self): 903 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] 904 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] 905 pl=PowerLaw(numMaterials=1,verbose=VERBOSE) 906 pl.setDruckerPragerLaw(tau_Y=100.) 907 pl.setPowerLaw(eta_N=2.) 908 pl.setElasticShearModulus(3.) 909 dt=1./3. 910 pl.setEtaTolerance(self.TOL) 911 self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0]) 912 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i]) 913 914 class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase): 915 TOL=1.e-6 916 VERBOSE=False 917 A=1. 918 P_max=100 919 NE=2*getMPISizeWorld() 920 tau_Y=10. 921 N_dt=10 922 923 # material parameter: 924 tau_1=5. 925 tau_2=5. 926 eta_0=100. 927 eta_1=50. 928 eta_2=400. 929 N_1=2. 930 N_2=3. 931 def getReference(self, t): 932 933 B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5) 934 x=self.dom.getX() 935 936 s_00=min(self.A*t,B) 937 tau=sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)*abs(s_00) 938 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.) 939 940 alpha=0.5*inv_eta*s_00 941 if s_00 <= B and self.mu !=None: alpha+=1./(2*self.mu)*self.A 942 u_ref=x*alpha 943 u_ref[self.dom.getDim()-1]=(1.-x[self.dom.getDim()-1])*alpha*(self.dom.getDim()-1) 944 sigma_ref=kronecker(self.dom)*s_00 945 sigma_ref[self.dom.getDim()-1,self.dom.getDim()-1]=-s_00*(self.dom.getDim()-1) 946 947 p_ref=self.P_max 948 for d in range(self.dom.getDim()): p_ref=p_ref*x[d] 949 p_ref-=integrate(p_ref)/vol(self.dom) 950 return u_ref, sigma_ref, p_ref 951 952 def runIt(self, free=None): 953 x=self.dom.getX() 954 B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5) 955 dt=B/int(self.N_dt/2) 956 if self.VERBOSE: print "dt =",dt 957 if self.latestart: 958 t=dt 959 else: 960 t=0 961 v,s,p=self.getReference(t) 962 963 mod=IncompressibleIsotropicFlowCartesian(self.dom, stress=s, v=v, p=p, t=t, numMaterials=3, verbose=self.VERBOSE) 964 mod.setDruckerPragerLaw(tau_Y=self.tau_Y,friction=None) 965 mod.setElasticShearModulus(self.mu) 966 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]) 967 mod.setTolerance(self.TOL) 968 mod.setFlowTolerance(self.TOL) 969 970 BF=Vector(self.P_max,Function(self.dom)) 971 for d in range(self.dom.getDim()): 972 for d2 in range(self.dom.getDim()): 973 if d!=d2: BF[d]*=x[d2] 974 v_mask=Vector(0,Solution(self.dom)) 975 if free==None: 976 for d in range(self.dom.getDim()): 977 v_mask[d]=whereZero(x[d])+whereZero(x[d]-1.) 978 else: 979 for d in range(self.dom.getDim()): 980 if d == self.dom.getDim()-1: 981 v_mask[d]=whereZero(x[d]-1.) 982 else: 983 v_mask[d]=whereZero(x[d]) 984 mod.setExternals(F=BF,fixed_v_mask=v_mask) 985 986 n=self.dom.getNormal() 987 N_t=0 988 errors=[] 989 while N_t < self.N_dt: 990 t_ref=t+dt 991 v_ref, s_ref,p_ref=self.getReference(t_ref) 992 mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref) 993 mod.update(dt, iter_max=100, inner_iter_max=20, verbose=self.VERBOSE, usePCG=True) 994 self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref) 995 t+=dt 996 N_t+=1 997 998 def check(self,N_t,mod,t_ref,v_ref, s_ref,p_ref): 999 p=mod.getPressure() 1000 p-=integrate(p)/vol(self.dom) 1001 error_p=Lsup(mod.getPressure()-p_ref)/Lsup(p_ref) 1002 error_s=Lsup(mod.getDeviatoricStress()-s_ref)/Lsup(s_ref) 1003 error_v=Lsup(mod.getVelocity()-v_ref)/Lsup(v_ref) 1004 error_t=abs(mod.getTime()-t_ref)/abs(t_ref) 1005 if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v 1006 self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) ) 1007 self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) ) 1008 self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) ) 1009 self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) ) 1010 def tearDown(self): 1011 del self.dom 1012 1013 def test_D2_Fixed_MuNone_LateStart(self): 1014 self.dom = Rectangle(self.NE,self.NE,order=2) 1015 self.mu=None 1016 self.latestart=True 1017 self.runIt() 1018 def test_D2_Fixed_Mu_LateStart(self): 1019 self.dom = Rectangle(self.NE,self.NE,order=2) 1020 self.mu=555. 1021 self.latestart=True 1022 self.runIt() 1023 def test_D2_Fixed_MuNone(self): 1024 self.dom = Rectangle(self.NE,self.NE,order=2) 1025 self.mu=None 1026 self.latestart=False 1027 self.runIt() 1028 def test_D2_Fixed_Mu(self): 1029 self.dom = Rectangle(self.NE,self.NE,order=2) 1030 self.mu=555. 1031 self.latestart=False 1032 self.runIt() 1033 def test_D2_Free_MuNone_LateStart(self): 1034 self.dom = Rectangle(self.NE,self.NE,order=2) 1035 self.mu=None 1036 self.latestart=True 1037 self.runIt(free=0) 1038 def test_D2_Free_Mu_LateStart(self): 1039 self.dom = Rectangle(self.NE,self.NE,order=2) 1040 self.mu=555. 1041 self.latestart=True 1042 self.runIt(free=0) 1043 def test_D2_Free_MuNone(self): 1044 self.dom = Rectangle(self.NE,self.NE,order=2) 1045 self.mu=None 1046 self.latestart=False 1047 self.runIt(free=0) 1048 def test_D2_Free_Mu(self): 1049 self.dom = Rectangle(self.NE,self.NE,order=2) 1050 self.mu=555. 1051 self.latestart=False 1052 self.runIt(free=0) 1053 1054 def test_D3_Fixed_MuNone_LateStart(self): 1055 self.dom = Brick(self.NE,self.NE,self.NE,order=2) 1056 self.mu=None 1057 self.latestart=True 1058 self.runIt() 1059 def test_D3_Fixed_Mu_LateStart(self): 1060 self.dom = Brick(self.NE,self.NE,self.NE,order=2) 1061 self.mu=555. 1062 self.latestart=True 1063 self.runIt() 1064 def test_D3_Fixed_MuNone(self): 1065 self.dom = Brick(self.NE,self.NE,self.NE,order=2) 1066 self.mu=None 1067 self.latestart=False 1068 self.runIt() 1069 def test_D3_Fixed_Mu(self): 1070 self.dom = Brick(self.NE,self.NE,self.NE,order=2) 1071 self.mu=555. 1072 self.latestart=False 1073 self.runIt() 1074 def test_D3_Free_MuNone_LateStart(self): 1075 self.dom = Brick(self.NE,self.NE,self.NE,order=2) 1076 self.mu=None 1077 self.latestart=True 1078 self.runIt(free=0) 1079 def test_D3_Free_Mu_LateStart(self): 1080 self.dom = Brick(self.NE,self.NE,self.NE,order=2) 1081 self.mu=555. 1082 self.latestart=True 1083 self.runIt(free=0) 1084 def test_D3_Free_MuNone(self): 1085 self.dom = Brick(self.NE,self.NE,self.NE,order=2) 1086 self.mu=None 1087 self.latestart=False 1088 self.runIt(free=0) 1089 def test_D3_Free_Mu(self): 1090 self.dom = Brick(self.NE,self.NE,self.NE,order=2) 1091 self.mu=555. 1092 self.latestart=False 1093 self.runIt(free=0) 1094 1095 1096 class Test_FaultSystem(unittest.TestCase): 1097 EPS=1.e-8 1098 NE=10 1099 def test_Fault_MaxValue(self): 1100 dom=Rectangle(2*self.NE,2*self.NE) 1101 x=dom.getX() 1102 f=FaultSystem(dim=2) 1103 f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1) 1104 f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2) 1105 1106 m, t, l=f.getMaxValue(x[0]*(1.-x[0])*(1-x[1])) 1107 self.failUnless( m == 0.25, "wrong max value") 1108 self.failUnless( t == 1, "wrong max tag") 1109 self.failUnless( l == 0., "wrong max location") 1110 m, t, l=f.getMaxValue(x[1]*(1.-x[1])*(1-x[0])*x[0]) 1111 self.failUnless( m == 0.0625, "wrong max value") 1112 self.failUnless( t == 2, "wrong max tag") 1113 self.failUnless( l == 0.5, "wrong max location") 1114 m, t, l=f.getMaxValue(x[0]*(1.-x[0])*x[1]) 1115 self.failUnless( m == 0.25, "wrong max value") 1116 self.failUnless( t == 2, "wrong max tag") 1117 self.failUnless( l == 1.0, "wrong max location") 1118 m, t, l= f.getMaxValue(x[1]*(1.-x[1])*x[0]) 1119 self.failUnless( m == 0.25, "wrong max value") 1120 self.failUnless( t == 2, "wrong max tag") 1121 self.failUnless( l == 0., "wrong max location") 1122 m, t, l= f.getMaxValue(x[1]*(1.-x[1])*(1.-x[0])) 1123 self.failUnless( m == 0.25, "wrong max value") 1124 self.failUnless( t == 1, "wrong max tag") 1125 self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong max location") 1126 def test_Fault_MinValue(self): 1127 dom=Rectangle(2*self.NE,2*self.NE) 1128 x=dom.getX() 1129 f=FaultSystem(dim=2) 1130 f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1) 1131 f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2) 1132 1133 m, t, l=f.getMinValue(-x[0]*(1.-x[0])*(1-x[1])) 1134 self.failUnless( m == -0.25, "wrong min value") 1135 self.failUnless( t == 1, "wrong min tag") 1136 self.failUnless( l == 0., "wrong min location") 1137 m, t, l=f.getMinValue(-x[1]*(1.-x[1])*(1-x[0])*x[0]) 1138 self.failUnless( m == -0.0625, "wrong min value") 1139 self.failUnless( t == 2, "wrong min tag") 1140 self.failUnless( l == 0.5, "wrong min location") 1141 m, t, l=f.getMinValue(-x[0]*(1.-x[0])*x[1]) 1142 self.failUnless( m == -0.25, "wrong min value") 1143 self.failUnless( t == 2, "wrong min tag") 1144 self.failUnless( l == 1.0, "wrong min location") 1145 m, t, l= f.getMinValue(-x[1]*(1.-x[1])*x[0]) 1146 self.failUnless( m == -0.25, "wrong min value") 1147 self.failUnless( t == 2, "wrong min tag") 1148 self.failUnless( l == 0., "wrong min location") 1149 m, t, l= f.getMinValue(-x[1]*(1.-x[1])*(1.-x[0])) 1150 self.failUnless( m == -0.25, "wrong min value") 1151 self.failUnless( t == 1, "wrong min tag") 1152 self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong min location") 1153 1154 1155 def test_Fault2D(self): 1156 f=FaultSystem(dim=2) 1157 top1=[ [1.,0.], [1.,1.], [0.,1.] ] 1158 self.failUnlessRaises(ValueError,f.addFault,V0=[1.,0],strikes=[pi/2, pi/2],ls=[1.,1.],tag=1,dips=top1) 1159 f.addFault(V0=[1.,0],strikes=[pi/2, pi],ls=[1.,1.],tag=1) 1160 self.failUnless(f.getDim() == 2, "wrong dimension") 1161 self.failUnless( [ 1 ] == f.getTags(), "tags wrong") 1162 self.failUnless( 2. == f.getTotalLength(1), "length wrong") 1163 self.failUnless( 0. == f.getMediumDepth(1), "depth wrong") 1164 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 range") 1165 self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 range") 1166 self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsets") 1167 segs=f.getTopPolyline(1) 1168 self.failUnless( len(segs) == 3, "wrong number of segments") 1169 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0") 1170 self.failUnless( segs[0].size == 2, "seg 0 has wrong size.") 1171 self.failUnless( numpy.linalg.norm(segs[0]-[1.,0.]) < self.EPS, "wrong vertex. 0 ") 1172 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1") 1173 self.failUnless( segs[1].size == 2, "seg 1 has wrong size.") 1174 self.failUnless( numpy.linalg.norm(segs[1]-[1.,1.]) < self.EPS, "wrong vertex. 1 ") 1175 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2") 1176 self.failUnless( segs[2].size == 2, "seg 2 has wrong size.") 1177 self.failUnless( numpy.linalg.norm(segs[2]-[0.,1.]) < self.EPS, "wrong vertex. 2 ") 1178 c=f.getCenterOnSurface() 1179 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class") 1180 self.failUnless( c.size == 2, "center size is wrong") 1181 self.failUnless( numpy.linalg.norm(c-[2./3.,2./3.]) < self.EPS, "center has wrong coordinates.") 1182 o=f.getOrientationOnSurface()/pi*180. 1183 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.") 1184 1185 top2=[ [10.,0.], [0.,10.] ] 1186 f.addFault(V0=[10.,0],strikes=[3.*pi/4],ls=[14.142135623730951], tag=2, w0_offsets=[0,20], w1_max=20) 1187 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong") 1188 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length") 1189 self.failUnless( 0. == f.getMediumDepth(2), "depth wrong") 1190 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range") 1191 self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range") 1192 self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets") 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( numpy.linalg.norm(segs[0]-[10.,0.]) < self.EPS, "wrong vertex. 0 ") 1197 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1") 1198 self.failUnless( numpy.linalg.norm(segs[1]-[0.,10.]) < self.EPS, "wrong vertex. 1 ") 1199 c=f.getCenterOnSurface() 1200 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class") 1201 self.failUnless( c.size == 2, "center size is wrong") 1202 self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.]) < self.EPS, "center has wrong coordinates.") 1203 o=f.getOrientationOnSurface()/pi*180. 1204 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.") 1205 1206 s,d=f.getSideAndDistance([0.,0.], tag=1) 1207 self.failUnless( s<0, "wrong side.") 1208 self.failUnless( abs(d-1.)0, "wrong side.") 1211 self.failUnless( abs(d-1.)0, "wrong side.") 1214 self.failUnless( abs(d-1.)0, "wrong side.") 1217 self.failUnless( abs(d-1.)0, "wrong side.") 1220 self.failUnless( abs(d-1.)0, "wrong side.") 1654 self.failUnless( abs(d-1.)0, "wrong side.") 1657 self.failUnless( abs(d-1.)0, "wrong side.") 1685 self.failUnless( abs(d-2*0.70710678118654757)0, "wrong side.") 1688 self.failUnless( abs(d-0.70710678118654757)

 ViewVC Help Powered by ViewVC 1.1.26