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

Revision 3528 - (show annotations)
Tue May 24 07:02:08 2011 UTC (9 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 73988 byte(s)

 1 # -*- coding: utf-8 -*- 2 3 ######################################################## 4 # 5 # Copyright (c) 2003-2010 by University of Queensland 6 # Earth Systems Science Computational Center (ESSCC) 7 8 # 9 # Primary Business: Queensland, Australia 10 # Licensed under the Open Software License version 3.0 11 12 # 13 ######################################################## 14 15 __copyright__="""Copyright (c) 2003-2010 by University of Queensland 16 Earth Systems Science Computational Center (ESSCC) 17 http://www.uq.edu.au/esscc 18 Primary Business: Queensland, Australia""" 19 __license__="""Licensed under the Open Software License version 3.0 20 21 __url__= 22 23 import unittest 24 import tempfile 25 26 27 28 VERBOSE=False and True 29 30 from esys.escript import * 31 from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian, FaultSystem, DarcyFlow 32 from esys.escript.models import Mountains 33 from esys.finley import Rectangle, Brick 34 35 from math import pi 36 import numpy 37 import sys 38 import os 39 #==================================================================================================================== 40 try: 41 FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR'] 42 except KeyError: 43 FINLEY_WORKDIR='.' 44 45 #==================================================================================================================== 46 class Test_StokesProblemCartesian2D(unittest.TestCase): 47 def setUp(self): 48 NE=6 49 self.TOL=1e-3 50 self.domain=Rectangle(NE,NE,order=-1) 51 def tearDown(self): 52 del self.domain 53 def test_PCG_P_0(self): 54 ETA=1. 55 P1=0. 56 57 x=self.domain.getX() 58 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.] 59 mask=whereZero(x[0]) * [1.,1.] \ 60 +whereZero(x[0]-1) * [1.,1.] \ 61 +whereZero(x[1]) * [1.,0.] \ 62 +whereZero(x[1]-1) * [1.,1.] 63 64 sp=StokesProblemCartesian(self.domain) 65 66 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 67 u0=(1-x[0])*x[0]*[0.,1.] 68 p0=Scalar(-P1,ReducedSolution(self.domain)) 69 sp.setTolerance(self.TOL) 70 u,p=sp.solve(u0*mask,p0,verbose=VERBOSE,max_iter=100,usePCG=True) 71 72 error_v0=Lsup(u[0]-u0[0]) 73 error_v1=Lsup(u[1]-u0[1])/0.25 74 error_p=Lsup(p+P1*x[0]*x[1]) 75 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 76 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 77 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 78 79 def test_PCG_P_small(self): 80 ETA=1. 81 P1=1. 82 83 x=self.domain.getX() 84 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.] 85 mask=whereZero(x[0]) * [1.,1.] \ 86 +whereZero(x[0]-1) * [1.,1.] \ 87 +whereZero(x[1]) * [1.,0.] \ 88 +whereZero(x[1]-1) * [1.,1.] 89 90 sp=StokesProblemCartesian(self.domain) 91 92 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 93 u0=(1-x[0])*x[0]*[0.,1.] 94 p0=Scalar(-P1,ReducedSolution(self.domain)) 95 sp.setTolerance(self.TOL) 96 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True) 97 error_v0=Lsup(u[0]-u0[0]) 98 error_v1=Lsup(u[1]-u0[1])/0.25 99 error_p=Lsup(P1*x[0]*x[1]+p) 100 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 101 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 102 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 103 104 def test_PCG_P_large(self): 105 ETA=1. 106 P1=1000. 107 108 x=self.domain.getX() 109 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.] 110 mask=whereZero(x[0]) * [1.,1.] \ 111 +whereZero(x[0]-1) * [1.,1.] \ 112 +whereZero(x[1]) * [1.,0.] \ 113 +whereZero(x[1]-1) * [1.,1.] 114 115 sp=StokesProblemCartesian(self.domain) 116 117 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 118 u0=(1-x[0])*x[0]*[0.,1.] 119 p0=Scalar(-P1,ReducedSolution(self.domain)) 120 sp.setTolerance(self.TOL) 121 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True) 122 # u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True) 123 124 error_v0=Lsup(u[0]-u0[0]) 125 error_v1=Lsup(u[1]-u0[1])/0.25 126 error_p=Lsup(P1*x[0]*x[1]+p)/P1 127 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 128 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 129 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 130 131 def test_GMRES_P_0(self): 132 ETA=1. 133 P1=0. 134 135 x=self.domain.getX() 136 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.] 137 mask=whereZero(x[0]) * [1.,1.] \ 138 +whereZero(x[0]-1) * [1.,1.] \ 139 +whereZero(x[1]) * [1.,0.] \ 140 +whereZero(x[1]-1) * [1.,1.] 141 142 sp=StokesProblemCartesian(self.domain) 143 144 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 145 u0=(1-x[0])*x[0]*[0.,1.] 146 p0=Scalar(-P1,ReducedSolution(self.domain)) 147 sp.setTolerance(self.TOL) 148 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=50,usePCG=False,iter_restart=18) 149 150 error_v0=Lsup(u[0]-u0[0]) 151 error_v1=Lsup(u[1]-u0[1])/0.25 152 error_p=Lsup(P1*x[0]*x[1]+p) 153 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 154 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 155 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 156 157 def test_GMRES_P_small(self): 158 ETA=1. 159 P1=1. 160 161 x=self.domain.getX() 162 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.] 163 mask=whereZero(x[0]) * [1.,1.] \ 164 +whereZero(x[0]-1) * [1.,1.] \ 165 +whereZero(x[1]) * [1.,0.] \ 166 +whereZero(x[1]-1) * [1.,1.] 167 168 sp=StokesProblemCartesian(self.domain) 169 170 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 171 u0=(1-x[0])*x[0]*[0.,1.] 172 p0=Scalar(-P1,ReducedSolution(self.domain)) 173 sp.setTolerance(self.TOL) 174 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=20,usePCG=False) 175 176 error_v0=Lsup(u[0]-u0[0]) 177 error_v1=Lsup(u[1]-u0[1])/0.25 178 error_p=Lsup(P1*x[0]*x[1]+p) 179 180 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 181 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 182 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 183 184 def test_GMRES_P_large(self): 185 ETA=1. 186 P1=1000. 187 188 x=self.domain.getX() 189 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.] 190 mask=whereZero(x[0]) * [1.,1.] \ 191 +whereZero(x[0]-1) * [1.,1.] \ 192 +whereZero(x[1]) * [1.,0.] \ 193 +whereZero(x[1]-1) * [1.,1.] 194 195 sp=StokesProblemCartesian(self.domain) 196 197 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 198 u0=(1-x[0])*x[0]*[0.,1.] 199 p0=Scalar(-P1,ReducedSolution(self.domain)) 200 sp.setTolerance(self.TOL) 201 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False) 202 203 error_v0=Lsup(u[0]-u0[0]) 204 error_v1=Lsup(u[1]-u0[1])/0.25 205 error_p=Lsup(P1*x[0]*x[1]+p)/P1 206 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 207 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 208 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 209 #==================================================================================================================== 210 class Test_StokesProblemCartesian3D(unittest.TestCase): 211 def setUp(self): 212 NE=6 213 self.TOL=1e-4 214 self.domain=Brick(NE,NE,NE,order=-1) 215 def tearDown(self): 216 del self.domain 217 def test_PCG_P_0(self): 218 ETA=1. 219 P1=0. 220 221 x=self.domain.getX() 222 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.] 223 x=self.domain.getX() 224 mask=whereZero(x[0]) * [1.,1.,1.] \ 225 +whereZero(x[0]-1) * [1.,1.,1.] \ 226 +whereZero(x[1]) * [1.,0.,1.] \ 227 +whereZero(x[1]-1) * [1.,1.,1.] \ 228 +whereZero(x[2]) * [1.,1.,0.] \ 229 +whereZero(x[2]-1) * [1.,1.,1.] 230 231 232 sp=StokesProblemCartesian(self.domain) 233 234 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 235 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.] 236 p0=Scalar(-P1,ReducedSolution(self.domain)) 237 sp.setTolerance(self.TOL) 238 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True) 239 240 error_v0=Lsup(u[0]-u0[0]) 241 error_v1=Lsup(u[1]-u0[1]) 242 error_v2=Lsup(u[2]-u0[2])/0.25**2 243 error_p=Lsup(P1*x[0]*x[1]*x[2]+p) 244 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 245 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 246 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 247 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.") 248 249 def test_PCG_P_small(self): 250 ETA=1. 251 P1=1. 252 253 x=self.domain.getX() 254 F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.] 255 mask=whereZero(x[0]) * [1.,1.,1.] \ 256 +whereZero(x[0]-1) * [1.,1.,1.] \ 257 +whereZero(x[1]) * [1.,0.,1.] \ 258 +whereZero(x[1]-1) * [1.,1.,1.] \ 259 +whereZero(x[2]) * [1.,1.,0.] \ 260 +whereZero(x[2]-1) * [1.,1.,1.] 261 262 263 sp=StokesProblemCartesian(self.domain) 264 265 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 266 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.] 267 p0=Scalar(-P1,ReducedSolution(self.domain)) 268 sp.setTolerance(self.TOL) 269 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True) 270 error_v0=Lsup(u[0]-u0[0]) 271 error_v1=Lsup(u[1]-u0[1]) 272 error_v2=Lsup(u[2]-u0[2])/0.25**2 273 error_p=Lsup(P1*x[0]*x[1]*x[2]+p) 274 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 275 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 276 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 277 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.") 278 279 def test_PCG_P_large(self): 280 ETA=1. 281 P1=1000. 282 283 x=self.domain.getX() 284 F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.] 285 mask=whereZero(x[0]) * [1.,1.,1.] \ 286 +whereZero(x[0]-1) * [1.,1.,1.] \ 287 +whereZero(x[1]) * [1.,0.,1.] \ 288 +whereZero(x[1]-1) * [1.,1.,1.] \ 289 +whereZero(x[2]) * [1.,1.,0.] \ 290 +whereZero(x[2]-1) * [1.,1.,1.] 291 292 293 sp=StokesProblemCartesian(self.domain) 294 295 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 296 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.] 297 p0=Scalar(-P1,ReducedSolution(self.domain)) 298 sp.setTolerance(self.TOL) 299 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True) 300 301 error_v0=Lsup(u[0]-u0[0]) 302 error_v1=Lsup(u[1]-u0[1]) 303 error_v2=Lsup(u[2]-u0[2])/0.25**2 304 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1 305 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 306 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 307 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 308 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.") 309 310 def test_GMRES_P_0(self): 311 ETA=1. 312 P1=0. 313 314 x=self.domain.getX() 315 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.] 316 x=self.domain.getX() 317 mask=whereZero(x[0]) * [1.,1.,1.] \ 318 +whereZero(x[0]-1) * [1.,1.,1.] \ 319 +whereZero(x[1]) * [1.,1.,1.] \ 320 +whereZero(x[1]-1) * [1.,1.,1.] \ 321 +whereZero(x[2]) * [1.,1.,0.] \ 322 +whereZero(x[2]-1) * [1.,1.,1.] 323 324 325 sp=StokesProblemCartesian(self.domain) 326 327 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 328 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.] 329 p0=Scalar(-P1,ReducedSolution(self.domain)) 330 sp.setTolerance(self.TOL) 331 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False,iter_restart=20) 332 333 error_v0=Lsup(u[0]-u0[0]) 334 error_v1=Lsup(u[1]-u0[1]) 335 error_v2=Lsup(u[2]-u0[2])/0.25**2 336 error_p=Lsup(P1*x[0]*x[1]*x[2]+p) 337 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 338 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 339 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 340 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.") 341 def test_GMRES_P_small(self): 342 ETA=1. 343 P1=1. 344 345 x=self.domain.getX() 346 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.] 347 mask=whereZero(x[0]) * [1.,1.,1.] \ 348 +whereZero(x[0]-1) * [1.,1.,1.] \ 349 +whereZero(x[1]) * [1.,1.,1.] \ 350 +whereZero(x[1]-1) * [1.,1.,1.] \ 351 +whereZero(x[2]) * [1.,1.,0.] \ 352 +whereZero(x[2]-1) * [1.,1.,1.] 353 354 355 sp=StokesProblemCartesian(self.domain) 356 357 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 358 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.] 359 p0=Scalar(-P1,ReducedSolution(self.domain)) 360 sp.setTolerance(self.TOL/10) 361 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False) 362 363 error_v0=Lsup(u[0]-u0[0]) 364 error_v1=Lsup(u[1]-u0[1]) 365 error_v2=Lsup(u[2]-u0[2])/0.25**2 366 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1 367 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 368 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 369 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.") 370 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 371 def test_GMRES_P_large(self): 372 ETA=1. 373 P1=1000. 374 375 x=self.domain.getX() 376 F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.] 377 mask=whereZero(x[0]) * [1.,1.,1.] \ 378 +whereZero(x[0]-1) * [1.,1.,1.] \ 379 +whereZero(x[1]) * [1.,0.,1.] \ 380 +whereZero(x[1]-1) * [1.,1.,1.] \ 381 +whereZero(x[2]) * [1.,1.,0.] \ 382 +whereZero(x[2]-1) * [1.,1.,1.] 383 384 385 sp=StokesProblemCartesian(self.domain) 386 387 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) 388 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.] 389 p0=Scalar(-P1,ReducedSolution(self.domain)) 390 sp.setTolerance(self.TOL) 391 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=False) 392 393 error_v0=Lsup(u[0]-u0[0]) 394 error_v1=Lsup(u[1]-u0[1]) 395 error_v2=Lsup(u[2]-u0[2])/0.25**2 396 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1 397 self.failUnless(error_p<10*self.TOL, "pressure error too large.") 398 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.") 399 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.") 400 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.") 401 #==================================================================================================================== 402 403 class Test_Mountains3D(unittest.TestCase): 404 def setUp(self): 405 NE=16 406 self.TOL=1e-4 407 self.domain=Brick(NE,NE,NE,order=1,useFullElementOrder=True) 408 def tearDown(self): 409 del self.domain 410 def test_periodic(self): 411 EPS=0.01 412 413 x=self.domain.getX() 414 v = Vector(0.0, Solution(self.domain)) 415 a0=1 416 a1=1 417 n0=2 418 n1=2 419 n2=0.5 420 a2=-(a0*n0+a1*n1)/n2 421 v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])* cos(pi*n2*x[2]) 422 v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])* cos(pi*n2*x[2]) 423 v[2]=a2*cos(pi*n0*x[0])* cos(pi*n1*x[1])* sin(pi*n2*x[2]) 424 425 mts=Mountains(self.domain,eps=EPS) 426 mts.setVelocity(v) 427 Z=mts.update() 428 429 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[2]-1.))) 430 self.failUnless(error_int=0,"eta needs to be positive (test %s)"%id) 543 error=abs(gamma_dot_*eta-tau_ref) 544 self.failUnless(error<=self.TOL*tau_ref,"eta is wrong: error = gamma_dot_*eta-tau_ref = %s * %s - %s = %s (test %s)"%(gamma_dot_,eta,tau_ref,error,id)) 545 546 def test_PowerLaw_Linear(self): 547 taus= [0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] 548 gamma_dot_s=[0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0, 75.0] 549 pl=PowerLaw(numMaterials=1,verbose=VERBOSE) 550 pl.setDruckerPragerLaw(tau_Y=100.) 551 pl.setPowerLaw(eta_N=2.) 552 pl.setEtaTolerance(self.TOL) 553 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 554 555 def test_PowerLaw_QuadLarge(self): 556 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] 557 gamma_dot_s=[0.0, 405.0, 1610.0, 3615.0, 6420.0, 10025.0, 14430.0, 19635.0, 25640.0, 32445.0, 40050.0, 44055.0, 48060.0, 52065.0, 56070.0, 60075.0] 558 pl=PowerLaw(numMaterials=2,verbose=VERBOSE) 559 pl.setDruckerPragerLaw(tau_Y=100.) 560 pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2]) 561 pl.setEtaTolerance(self.TOL) 562 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 563 564 def test_PowerLaw_QuadSmall(self): 565 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] 566 gamma_dot_s=[0.0, 5.4, 11.6, 18.6, 26.4, 35.0, 44.4, 54.6, 65.6, 77.4, 90.0, 99.0, 108.0, 117.0, 126.0, 135.0] 567 pl=PowerLaw(numMaterials=2,verbose=VERBOSE) 568 pl.setDruckerPragerLaw(tau_Y=100.) 569 pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2]) 570 pl.setEtaTolerance(self.TOL) 571 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 572 573 def test_PowerLaw_CubeLarge(self): 574 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] 575 gamma_dot_s=[0.0, 8.90625, 41.25, 120.46875, 270.0, 513.28125, 873.75, 1374.84375, 2040.0, 2892.65625, 3956.25, 4351.875, 4747.5, 5143.125, 5538.75, 5934.375] 576 pl=PowerLaw(numMaterials=2,verbose=VERBOSE) 577 pl.setDruckerPragerLaw(tau_Y=100.) 578 pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3]) 579 pl.setEtaTolerance(self.TOL) 580 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 581 582 def test_PowerLaw_CubeSmall(self): 583 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] 584 gamma_dot_s=[0.0, 5.0390625, 10.3125, 16.0546875, 22.5, 29.8828125, 38.4375, 48.3984375, 60.0, 73.4765625, 89.0625, 97.96875, 106.875, 115.78125, 124.6875, 133.59375] 585 pl=PowerLaw(numMaterials=2,verbose=VERBOSE) 586 pl.setDruckerPragerLaw(tau_Y=100.) 587 pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3]) 588 pl.setEtaTolerance(self.TOL) 589 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 590 591 def test_PowerLaw_QuadLarge_CubeLarge(self): 592 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] 593 gamma_dot_s=[0.0, 408.90625, 1641.25, 3720.46875, 6670.0, 10513.28125, 15273.75, 20974.84375, 27640.000000000004, 35292.65625, 43956.25, 48351.875, 52747.5, 57143.125, 61538.75, 65934.375] 594 595 pl=PowerLaw(numMaterials=3,verbose=VERBOSE) 596 pl.setDruckerPragerLaw(tau_Y=100.) 597 pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3]) 598 pl.setEtaTolerance(self.TOL) 599 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 600 601 def test_PowerLaw_QuadLarge_CubeSmall(self): 602 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] 603 gamma_dot_s=[0.0, 405.0390625, 1610.3125, 3616.0546875, 6422.5, 10029.8828125, 14438.4375, 19648.3984375, 25660.0, 32473.4765625, 40089.0625, 44097.96875, 48106.875, 52115.78125, 56124.6875, 60133.59375] 604 605 pl=PowerLaw(numMaterials=3,verbose=VERBOSE) 606 pl.setDruckerPragerLaw(tau_Y=100.) 607 pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3]) 608 pl.setEtaTolerance(self.TOL) 609 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 610 611 def test_PowerLaw_QuadSmall_CubeLarge(self): 612 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] 613 gamma_dot_s=[0.0, 9.30625, 42.85, 124.06875, 276.4, 523.28125, 888.15, 1394.44375, 2065.6, 2925.05625, 3996.25, 4395.875, 4795.5, 5195.125, 5594.75, 5994.375] 614 pl=PowerLaw(numMaterials=3,verbose=VERBOSE) 615 pl.setDruckerPragerLaw(tau_Y=100.) 616 pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3]) 617 pl.setEtaTolerance(self.TOL) 618 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 619 620 def test_PowerLaw_QuadSmall_CubeSmall(self): 621 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] 622 gamma_dot_s=[0.0, 5.4390625, 11.9125, 19.6546875, 28.9, 39.8828125, 52.8375, 67.9984375, 85.6, 105.8765625, 129.0625, 141.96875, 154.875, 167.78125, 180.6875, 193.59375] 623 pl=PowerLaw(numMaterials=3,verbose=VERBOSE) 624 pl.setDruckerPragerLaw(tau_Y=100.) 625 pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3]) 626 pl.setEtaTolerance(self.TOL) 627 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i]) 628 629 def test_PowerLaw_withShear(self): 630 taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0] 631 gamma_dot_s=[0.0, 15.0, 30.0, 45.0, 60.0, 75.0, 90.0, 105.0, 120.0, 135.0, 150.0, 165.0, 180.0, 195.0, 210.0, 225.0] 632 pl=PowerLaw(numMaterials=1,verbose=VERBOSE) 633 pl.setDruckerPragerLaw(tau_Y=100.) 634 pl.setPowerLaw(eta_N=2.) 635 pl.setElasticShearModulus(3.) 636 dt=1./3. 637 pl.setEtaTolerance(self.TOL) 638 self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0]) 639 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i]) 640 641 class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase): 642 TOL=1.e-5 643 VERBOSE=False # or True 644 A=1. 645 P_max=100 646 NE=2*getMPISizeWorld() 647 tau_Y=10. 648 N_dt=10 649 650 # material parameter: 651 tau_1=5. 652 tau_2=5. 653 eta_0=100. 654 eta_1=50. 655 eta_2=400. 656 N_1=2. 657 N_2=3. 658 def getReference(self, t): 659 660 B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5) 661 x=self.dom.getX() 662 663 s_00=min(self.A*t,B) 664 tau=sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)*abs(s_00) 665 inv_eta= 1./self.eta_0 + 1./self.eta_1*(tau/self.tau_1)**(self.N_1-1.) + 1./self.eta_2*(tau/self.tau_2)**(self.N_2-1.) 666 667 alpha=0.5*inv_eta*s_00 668 if s_00 <= B and self.mu !=None: alpha+=1./(2*self.mu)*self.A 669 u_ref=x*alpha 670 u_ref[self.dom.getDim()-1]=(1.-x[self.dom.getDim()-1])*alpha*(self.dom.getDim()-1) 671 sigma_ref=kronecker(self.dom)*s_00 672 sigma_ref[self.dom.getDim()-1,self.dom.getDim()-1]=-s_00*(self.dom.getDim()-1) 673 674 p_ref=self.P_max 675 for d in range(self.dom.getDim()): p_ref=p_ref*x[d] 676 p_ref-=integrate(p_ref)/vol(self.dom) 677 return u_ref, sigma_ref, p_ref 678 679 def runIt(self, free=None): 680 x=self.dom.getX() 681 B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5) 682 dt=B/int(self.N_dt/2) 683 if self.VERBOSE: print "dt =",dt 684 if self.latestart: 685 t=dt 686 else: 687 t=0 688 v,s,p=self.getReference(t) 689 690 mod=IncompressibleIsotropicFlowCartesian(self.dom, stress=s, v=v, p=p, t=t, numMaterials=3, verbose=self.VERBOSE) 691 mod.setDruckerPragerLaw(tau_Y=self.tau_Y,friction=None) 692 mod.setElasticShearModulus(self.mu) 693 mod.setPowerLaws([self.eta_0, self.eta_1, self.eta_2], [ 1., self.tau_1, self.tau_2], [1.,self.N_1,self.N_2]) 694 mod.setTolerance(self.TOL) 695 mod.setEtaTolerance(self.TOL*0.1) 696 697 BF=Vector(self.P_max,Function(self.dom)) 698 for d in range(self.dom.getDim()): 699 for d2 in range(self.dom.getDim()): 700 if d!=d2: BF[d]*=x[d2] 701 v_mask=Vector(0,Solution(self.dom)) 702 if free==None: 703 for d in range(self.dom.getDim()): 704 v_mask[d]=whereZero(x[d])+whereZero(x[d]-1.) 705 else: 706 for d in range(self.dom.getDim()): 707 if d == self.dom.getDim()-1: 708 v_mask[d]=whereZero(x[d]-1.) 709 else: 710 v_mask[d]=whereZero(x[d]) 711 mod.setExternals(F=BF,fixed_v_mask=v_mask) 712 713 n=self.dom.getNormal() 714 N_t=0 715 errors=[] 716 while N_t < self.N_dt: 717 t_ref=t+dt 718 v_ref, s_ref,p_ref=self.getReference(t_ref) 719 mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref) 720 # mod.update(dt, eta_iter_max=10, iter_max=50, verbose=self.VERBOSE, usePCG=True, max_correction_steps=30) new version 721 mod.update(dt, iter_max=50, verbose=self.VERBOSE, usePCG=True) 722 self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref) 723 t+=dt 724 N_t+=1 725 726 def check(self,N_t,mod,t_ref,v_ref, s_ref,p_ref): 727 p=mod.getPressure() 728 p-=integrate(p)/vol(self.dom) 729 error_p=Lsup(mod.getPressure()-p_ref)/Lsup(p_ref) 730 error_s=Lsup(mod.getDeviatoricStress()-s_ref)/Lsup(s_ref) 731 error_v=Lsup(mod.getVelocity()-v_ref)/Lsup(v_ref) 732 error_t=abs(mod.getTime()-t_ref)/abs(t_ref) 733 if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v 734 self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) ) 735 self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) ) 736 self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) ) 737 self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) ) 738 def tearDown(self): 739 del self.dom 740 741 def test_D2_Fixed_MuNone_LateStart(self): 742 self.dom = Rectangle(self.NE,self.NE,order=2) 743 self.mu=None 744 self.latestart=True 745 self.runIt() 746 def test_D2_Fixed_Mu_LateStart(self): 747 self.dom = Rectangle(self.NE,self.NE,order=2) 748 self.mu=555. 749 self.latestart=True 750 self.runIt() 751 def test_D2_Fixed_MuNone(self): 752 self.dom = Rectangle(self.NE,self.NE,order=2) 753 self.mu=None 754 self.latestart=False 755 self.runIt() 756 def test_D2_Fixed_Mu(self): 757 self.dom = Rectangle(self.NE,self.NE,order=2) 758 self.mu=555. 759 self.latestart=False 760 self.runIt() 761 def test_D2_Free_MuNone_LateStart(self): 762 self.dom = Rectangle(self.NE,self.NE,order=2) 763 self.mu=None 764 self.latestart=True 765 self.runIt(free=0) 766 def test_D2_Free_Mu_LateStart(self): 767 self.dom = Rectangle(self.NE,self.NE,order=2) 768 self.mu=555. 769 self.latestart=True 770 self.runIt(free=0) 771 def test_D2_Free_MuNone(self): 772 self.dom = Rectangle(self.NE,self.NE,order=2) 773 self.mu=None 774 self.latestart=False 775 self.runIt(free=0) 776 def test_D2_Free_Mu(self): 777 self.dom = Rectangle(self.NE,self.NE,order=2) 778 self.mu=555. 779 self.latestart=False 780 self.runIt(free=0) 781 782 def test_D3_Fixed_MuNone_LateStart(self): 783 self.dom = Brick(self.NE,self.NE,self.NE,order=2) 784 self.mu=None 785 self.latestart=True 786 self.runIt() 787 def test_D3_Fixed_Mu_LateStart(self): 788 self.dom = Brick(self.NE,self.NE,self.NE,order=2) 789 self.mu=555. 790 self.latestart=True 791 self.runIt() 792 def test_D3_Fixed_MuNone(self): 793 self.dom = Brick(self.NE,self.NE,self.NE,order=2) 794 self.mu=None 795 self.latestart=False 796 self.runIt() 797 def test_D3_Fixed_Mu(self): 798 self.dom = Brick(self.NE,self.NE,self.NE,order=2) 799 self.mu=555. 800 self.latestart=False 801 self.runIt() 802 def test_D3_Free_MuNone_LateStart(self): 803 self.dom = Brick(self.NE,self.NE,self.NE,order=2) 804 self.mu=None 805 self.latestart=True 806 self.runIt(free=0) 807 def test_D3_Free_Mu_LateStart(self): 808 self.dom = Brick(self.NE,self.NE,self.NE,order=2) 809 self.mu=555. 810 self.latestart=True 811 self.runIt(free=0) 812 def test_D3_Free_MuNone(self): 813 self.dom = Brick(self.NE,self.NE,self.NE,order=2) 814 self.mu=None 815 self.latestart=False 816 self.runIt(free=0) 817 def test_D3_Free_Mu(self): 818 self.dom = Brick(self.NE,self.NE,self.NE,order=2) 819 self.mu=555. 820 self.latestart=False 821 self.runIt(free=0) 822 823 824 class Test_FaultSystem(unittest.TestCase): 825 EPS=1.e-8 826 NE=10 827 def test_Fault_MaxValue(self): 828 dom=Rectangle(2*self.NE,2*self.NE) 829 x=dom.getX() 830 f=FaultSystem(dim=2) 831 f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1) 832 f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2) 833 834 u=x[0]*(1.-x[0])*(1-x[1]) 835 t, loc=f.getMaxValue(u) 836 p=f.getParametrization(x,t)[0] 837 m, l=loc(u), loc(p) 838 self.failUnless( m == 0.25, "wrong max value") 839 self.failUnless( t == 1, "wrong max tag") 840 self.failUnless( l == 0., "wrong max location") 841 842 u=x[1]*(1.-x[1])*(1-x[0])*x[0] 843 t, loc=f.getMaxValue(u) 844 p=f.getParametrization(x,t)[0] 845 m, l=loc(u), loc(p) 846 self.failUnless( m == 0.0625, "wrong max value") 847 self.failUnless( t == 2, "wrong max tag") 848 self.failUnless( l == 0.5, "wrong max location") 849 850 u=x[0]*(1.-x[0])*x[1] 851 t, loc=f.getMaxValue(u) 852 p=f.getParametrization(x,t)[0] 853 m, l=loc(u), loc(p) 854 self.failUnless( m == 0.25, "wrong max value") 855 self.failUnless( t == 2, "wrong max tag") 856 self.failUnless( l == 1.0, "wrong max location") 857 858 u=x[1]*(1.-x[1])*x[0] 859 t, loc=f.getMaxValue(u) 860 p=f.getParametrization(x,t)[0] 861 m, l=loc(u), loc(p) 862 self.failUnless( m == 0.25, "wrong max value") 863 self.failUnless( t == 2, "wrong max tag") 864 self.failUnless( l == 0., "wrong max location") 865 866 u=x[1]*(1.-x[1])*(1.-x[0]) 867 t, loc=f.getMaxValue(u) 868 p=f.getParametrization(x,t)[0] 869 m, l=loc(u), loc(p) 870 self.failUnless( m == 0.25, "wrong max value") 871 self.failUnless( t == 1, "wrong max tag") 872 self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong max location") 873 def test_Fault_MinValue(self): 874 dom=Rectangle(2*self.NE,2*self.NE) 875 x=dom.getX() 876 f=FaultSystem(dim=2) 877 f.addFault(V0=[0.5,0.], strikes=[3.*pi/4], ls=[0.70710678118654757], tag=1) 878 f.addFault(V0=[1.,0.5], strikes=[pi, pi/2], ls=[0.5,0.5], tag=2) 879 880 u=-x[0]*(1.-x[0])*(1-x[1]) 881 t, loc=f.getMinValue(u) 882 p=f.getParametrization(x,t)[0] 883 m, l=loc(u), loc(p) 884 self.failUnless( m == -0.25, "wrong min value") 885 self.failUnless( t == 1, "wrong min tag") 886 self.failUnless( l == 0., "wrong min location") 887 u=-x[1]*(1.-x[1])*(1-x[0])*x[0] 888 t, loc=f.getMinValue(u) 889 p=f.getParametrization(x,t)[0] 890 m, l=loc(u), loc(p) 891 self.failUnless( m == -0.0625, "wrong min value") 892 self.failUnless( t == 2, "wrong min tag") 893 self.failUnless( l == 0.5, "wrong min location") 894 u=-x[0]*(1.-x[0])*x[1] 895 t, loc=f.getMinValue(u) 896 p=f.getParametrization(x,t)[0] 897 m, l=loc(u), loc(p) 898 self.failUnless( m == -0.25, "wrong min value") 899 self.failUnless( t == 2, "wrong min tag") 900 self.failUnless( l == 1.0, "wrong min location") 901 u=-x[1]*(1.-x[1])*x[0] 902 t, loc=f.getMinValue(u) 903 p=f.getParametrization(x,t)[0] 904 m, l=loc(u), loc(p) 905 self.failUnless( m == -0.25, "wrong min value") 906 self.failUnless( t == 2, "wrong min tag") 907 self.failUnless( l == 0., "wrong min location") 908 u=-x[1]*(1.-x[1])*(1.-x[0]) 909 t, loc=f.getMinValue(u) 910 p=f.getParametrization(x,t)[0] 911 m, l=loc(u), loc(p) 912 self.failUnless( m == -0.25, "wrong min value") 913 self.failUnless( t == 1, "wrong min tag") 914 self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong min location") 915 916 917 def test_Fault2D(self): 918 f=FaultSystem(dim=2) 919 top1=[ [1.,0.], [1.,1.], [0.,1.] ] 920 self.failUnlessRaises(ValueError,f.addFault,V0=[1.,0],strikes=[pi/2, pi/2],ls=[1.,1.],tag=1,dips=top1) 921 f.addFault(V0=[1.,0],strikes=[pi/2, pi],ls=[1.,1.],tag=1) 922 self.failUnless(f.getDim() == 2, "wrong dimension") 923 self.failUnless( [ 1 ] == f.getTags(), "tags wrong") 924 self.failUnless( 2. == f.getTotalLength(1), "length wrong") 925 self.failUnless( 0. == f.getMediumDepth(1), "depth wrong") 926 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 range") 927 self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 range") 928 self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsets") 929 segs=f.getTopPolyline(1) 930 self.failUnless( len(segs) == 3, "wrong number of segments") 931 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0") 932 self.failUnless( segs[0].size == 2, "seg 0 has wrong size.") 933 self.failUnless( numpy.linalg.norm(segs[0]-[1.,0.]) < self.EPS, "wrong vertex. 0 ") 934 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1") 935 self.failUnless( segs[1].size == 2, "seg 1 has wrong size.") 936 self.failUnless( numpy.linalg.norm(segs[1]-[1.,1.]) < self.EPS, "wrong vertex. 1 ") 937 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2") 938 self.failUnless( segs[2].size == 2, "seg 2 has wrong size.") 939 self.failUnless( numpy.linalg.norm(segs[2]-[0.,1.]) < self.EPS, "wrong vertex. 2 ") 940 c=f.getCenterOnSurface() 941 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class") 942 self.failUnless( c.size == 2, "center size is wrong") 943 self.failUnless( numpy.linalg.norm(c-[2./3.,2./3.]) < self.EPS, "center has wrong coordinates.") 944 o=f.getOrientationOnSurface()/pi*180. 945 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.") 946 947 top2=[ [10.,0.], [0.,10.] ] 948 f.addFault(V0=[10.,0],strikes=[3.*pi/4],ls=[14.142135623730951], tag=2, w0_offsets=[0,20], w1_max=20) 949 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong") 950 self.failUnless( abs(f.getTotalLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length") 951 self.failUnless( 0. == f.getMediumDepth(2), "depth wrong") 952 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range") 953 self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range") 954 self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets") 955 segs=f.getTopPolyline(2) 956 self.failUnless( len(segs) == 2, "wrong number of segments") 957 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0") 958 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.]) < self.EPS, "wrong vertex. 0 ") 959 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1") 960 self.failUnless( numpy.linalg.norm(segs[1]-[0.,10.]) < self.EPS, "wrong vertex. 1 ") 961 c=f.getCenterOnSurface() 962 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class") 963 self.failUnless( c.size == 2, "center size is wrong") 964 self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.]) < self.EPS, "center has wrong coordinates.") 965 o=f.getOrientationOnSurface()/pi*180. 966 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.") 967 968 s,d=f.getSideAndDistance([0.,0.], tag=1) 969 self.failUnless( s<0, "wrong side.") 970 self.failUnless( abs(d-1.)0, "wrong side.") 973 self.failUnless( abs(d-1.)0, "wrong side.") 976 self.failUnless( abs(d-1.)0, "wrong side.") 979 self.failUnless( abs(d-1.)0, "wrong side.") 982 self.failUnless( abs(d-1.)0, "wrong side.") 1416 self.failUnless( abs(d-1.)0, "wrong side.") 1419 self.failUnless( abs(d-1.)0, "wrong side.") 1447 self.failUnless( abs(d-2*0.70710678118654757)0, "wrong side.") 1450 self.failUnless( abs(d-0.70710678118654757)

 ViewVC Help Powered by ViewVC 1.1.26