# Contents of /trunk-mpi-branch/escript/py_src/linearPDEs.py

Revision 148 - (show annotations)
Tue Aug 23 01:24:31 2005 UTC (14 years, 1 month ago) by jgs
Original Path: trunk/esys2/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 68376 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-08-23


 1 # $Id$ 2 3 ## @file linearPDEs.py 4 5 """ 6 Functions and classes for linear PDEs 7 """ 8 9 import escript 10 import util 11 import numarray 12 13 14 class IllegalCoefficient(ValueError): 15 """ 16 raised if an illegal coefficient of the general ar particular PDE is requested. 17 """ 18 19 class IllegalCoefficientValue(ValueError): 20 """ 21 raised if an incorrect value for a coefficient is used. 22 """ 23 24 class UndefinedPDEError(ValueError): 25 """ 26 raised if a PDE is not fully defined yet. 27 """ 28 29 def _CompTuple2(t1,t2): 30 """ 31 Compare two tuples 32 33 @param t1 The first tuple 34 @param t2 The second tuple 35 36 """ 37 38 dif=t1[0]+t1[1]-(t2[0]+t2[1]) 39 if dif<0: return 1 40 elif dif>0: return -1 41 else: return 0 42 43 class PDECoefficient: 44 """ 45 A class for PDE coefficients 46 """ 47 # identifier for location of Data objects defining COEFFICIENTS 48 INTERIOR=0 49 BOUNDARY=1 50 CONTACT=2 51 CONTINUOUS=3 52 # identifier in the pattern of COEFFICIENTS: 53 # the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the 54 # number of unknowns. 55 EQUATION=3 56 SOLUTION=4 57 DIM=5 58 # indicator for what is altered if the coefficient is altered: 59 OPERATOR=5 60 RIGHTHANDSIDE=6 61 BOTH=7 62 def __init__(self,where,pattern,altering): 63 """ 64 Initialise a PDE Coefficient type 65 """ 66 self.what=where 67 self.pattern=pattern 68 self.altering=altering 69 self.resetValue() 70 71 def resetValue(self): 72 """ 73 resets coefficient value to default 74 """ 75 self.value=escript.Data() 76 77 def getFunctionSpace(self,domain): 78 """ 79 defines the FunctionSpace of the coefficient on the domain 80 81 @param domain: 82 """ 83 if self.what==self.INTERIOR: return escript.Function(domain) 84 elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain) 85 elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain) 86 elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain) 87 88 def getValue(self): 89 """ 90 returns the value of the coefficient: 91 """ 92 return self.value 93 94 def setValue(self,domain,numEquations=1,numSolutions=1,newValue=None): 95 """ 96 set the value of the coefficient to new value 97 """ 98 if newValue==None: 99 newValue=escript.Data() 100 elif isinstance(newValue,escript.Data): 101 if not newValue.isEmpty(): 102 newValue=escript.Data(newValue,self.getFunctionSpace(domain)) 103 else: 104 newValue=escript.Data(newValue,self.getFunctionSpace(domain)) 105 if not newValue.isEmpty(): 106 if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape(): 107 raise IllegalCoefficientValue,"Expected shape for coefficient %s is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape()) 108 self.value=newValue 109 110 def isAlteringOperator(self): 111 """ 112 return true if the operator of the PDE is changed when the coefficient is changed 113 """ 114 if self.altering==self.OPERATOR or self.altering==self.BOTH: 115 return not None 116 else: 117 return None 118 119 def isAlteringRightHandSide(self): 120 """ 121 return true if the right hand side of the PDE is changed when the coefficient is changed 122 """ 123 if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH: 124 return not None 125 else: 126 return None 127 128 def estimateNumEquationsAndNumSolutions(self,domain,shape=()): 129 """ 130 tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim 131 132 @param shape: 133 @param dim: 134 """ 135 dim=domain.getDim() 136 if len(shape)>0: 137 num=max(shape)+1 138 else: 139 num=1 140 search=[] 141 for u in range(num): 142 for e in range(num): 143 search.append((e,u)) 144 search.sort(_CompTuple2) 145 for item in search: 146 s=self.getShape(domain,item[0],item[1]) 147 if len(s)==0 and len(shape)==0: 148 return (1,1) 149 else: 150 if s==shape: return item 151 return None 152 153 def getShape(self,domain,numEquations=1,numSolutions=1): 154 """ 155 builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim 156 157 @param e: 158 @param u: 159 @param dim: 160 """ 161 dim=domain.getDim() 162 s=() 163 for i in self.pattern: 164 if i==self.EQUATION: 165 if numEquations>1: s=s+(numEquations,) 166 elif i==self.SOLUTION: 167 if numSolutions>1: s=s+(numSolutions,) 168 else: 169 s=s+(dim,) 170 return s 171 172 class LinearPDE: 173 """ 174 Class to define a linear PDE of the form 175 176 \f[ 177 -(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i 178 \f] 179 180 with boundary conditons: 181 182 \f[ 183 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i 184 \f] 185 186 and contact conditions 187 188 \f[ 189 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i 190 \f] 191 192 and constraints: 193 194 \f[ 195 u_i=r_i \quad \mathrm{where} \quad q_i>0 196 \f] 197 198 """ 199 TOL=1.e-13 200 # solver methods 201 UNKNOWN=-1 202 DEFAULT_METHOD=0 203 DIRECT=1 204 CHOLEVSKY=2 205 PCG=3 206 CR=4 207 CGS=5 208 BICGSTAB=6 209 SSOR=7 210 ILU0=8 211 ILUT=9 212 JACOBI=10 213 GMRES=11 214 PRES20=12 215 LUMPING=13 216 # matrix reordering: 217 NO_REORDERING=0 218 MINIMUM_FILL_IN=1 219 NESTED_DISSECTION=2 220 # important keys in the dictonary used to hand over options to the solver library: 221 METHOD_KEY="method" 222 SYMMETRY_KEY="symmetric" 223 TOLERANCE_KEY="tolerance" 224 225 226 def __init__(self,domain,numEquations=None,numSolutions=None,debug=False): 227 """ 228 initializes a new linear PDE 229 230 @param domain: domain of the PDE 231 @type domain: L{Domain} 232 @param numEquations: number of equations. If numEquations==None the number of equations 233 is exracted from the PDE coefficients. 234 @param numSolutions: number of solution components. If numSolutions==None the number of solution components 235 is exracted from the PDE coefficients. 236 @param debug: if True debug informations are printed. 237 238 239 """ 240 # 241 # the coefficients of the general PDE: 242 # 243 self.__COEFFICIENTS_OF_GENEARL_PDE={ 244 "A" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR), 245 "B" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR), 246 "C" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR), 247 "D" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR), 248 "X" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE), 249 "Y" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE), 250 "d" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR), 251 "y" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE), 252 "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR), 253 "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE), 254 "r" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE), 255 "q" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)} 256 257 # COEFFICIENTS can be overwritten by subclasses: 258 self.COEFFICIENTS=self.__COEFFICIENTS_OF_GENEARL_PDE 259 # initialize attributes 260 self.__debug=debug 261 self.__domain=domain 262 self.__numEquations=numEquations 263 self.__numSolutions=numSolutions 264 self.__resetSystem() 265 266 # set some default values: 267 self.__homogeneous_constraint=True 268 self.__row_function_space=escript.Solution(self.__domain) 269 self.__column_function_space=escript.Solution(self.__domain) 270 self.__tolerance=1.e-8 271 self.__solver_method=self.DEFAULT_METHOD 272 self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT_METHOD,False) 273 self.__sym=False 274 275 self.resetCoefficients() 276 self.trace("PDE Coeffients are %s"%str(self.COEFFICIENTS.keys())) 277 # ============================================================================= 278 # general stuff: 279 # ============================================================================= 280 def __str__(self): 281 return ""%id(self) 282 # ============================================================================= 283 # debug : 284 # ============================================================================= 285 def setDebugOn(self): 286 """ 287 switches on debugging 288 """ 289 self.__debug=not None 290 291 def setDebugOff(self): 292 """ 293 switches off debugging 294 """ 295 self.__debug=None 296 297 def trace(self,text): 298 """ 299 print the text message if debugging is swiched on. 300 301 @param name: name of the coefficient enquired. 302 @type name: C{string} 303 """ 304 if self.__debug: print "%s: %s"%(str(self),text) 305 306 # ============================================================================= 307 # some service functions: 308 # ============================================================================= 309 def getDomain(self): 310 """ 311 returns the domain of the PDE 312 313 @return : the domain of the PDE 314 @rtype : C{Domain} 315 316 """ 317 return self.__domain 318 319 def getDim(self): 320 """ 321 returns the spatial dimension of the PDE 322 323 @return : the spatial dimension of the PDE domain 324 @rtype : C{int} 325 """ 326 return self.getDomain().getDim() 327 328 def getNumEquations(self): 329 """ 330 returns the number of equations 331 332 @return : the number of equations 333 @rtype : C{int} 334 @raise UndefinedPDEError: if the number of equations is not be specified yet. 335 """ 336 if self.__numEquations==None: 337 raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations." 338 else: 339 return self.__numEquations 340 341 def getNumSolutions(self): 342 """ 343 returns the number of unknowns 344 345 @return : the number of unknowns 346 @rtype : C{int} 347 @raise UndefinedPDEError: if the number of unknowns is not be specified yet. 348 """ 349 if self.__numSolutions==None: 350 raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions." 351 else: 352 return self.__numSolutions 353 354 def getFunctionSpaceForEquation(self): 355 """ 356 returns the L{escript.FunctionSpace} used to discretize the equation 357 358 @return : representation space of equation 359 @rtype : L{escript.FunctionSpace} 360 361 """ 362 return self.__row_function_space 363 364 def getFunctionSpaceForSolution(self): 365 """ 366 returns the L{escript.FunctionSpace} used to represent the solution 367 368 @return : representation space of solution 369 @rtype : L{escript.FunctionSpace} 370 371 """ 372 return self.__column_function_space 373 374 375 def getOperator(self): 376 """ 377 provides access to the operator of the PDE 378 379 @return : the operator of the PDE 380 @rtype : L{Operator} 381 """ 382 m=self.getSystem()[0] 383 if self.isUsingLumping(): 384 return self.copyConstraint(1./m) 385 else: 386 return m 387 388 def getRightHandSide(self): 389 """ 390 provides access to the right hand side of the PDE 391 392 @return : the right hand side of the PDE 393 @rtype : L{escript.Data} 394 """ 395 r=self.getSystem()[1] 396 if self.isUsingLumping(): 397 return self.copyConstraint(r) 398 else: 399 return r 400 401 def applyOperator(self,u=None): 402 """ 403 applies the operator of the PDE to a given u or the solution of PDE if u is not present. 404 405 @param u: argument of the operator. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None} 406 the current solution is used. 407 @type u: L{escript.Data} or None 408 @return : image of u 409 @rtype : L{escript.Data} 410 """ 411 if u==None: 412 return self.getOperator()*self.getSolution() 413 else: 414 self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution()) 415 416 def getResidual(self,u=None): 417 """ 418 return the residual of u or the current solution if u is not present. 419 420 @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None} 421 the current solution is used. 422 @type u: L{escript.Data} or None 423 @return : residual of u 424 @rtype : L{escript.Data} 425 """ 426 return self.applyOperator(u)-self.getRightHandSide() 427 428 def checkSymmetry(self,verbose=True): 429 """ 430 test the PDE for symmetry. 431 432 433 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed. 434 @type verbose: C{bool} 435 @return: True if the PDE is symmetric. 436 @rtype : C{escript.Data} 437 438 @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered. 439 """ 440 verbose=verbose or self.debug() 441 out=True 442 if self.getNumSolutions()!=self.getNumEquations(): 443 if verbose: print "non-symmetric PDE because of different number of equations and solutions" 444 out=False 445 else: 446 A=self.getCoefficientOfGeneralPDE("A") 447 if not A.isEmpty(): 448 tol=util.Lsup(A)*self.TOL 449 if self.getNumSolutions()>1: 450 for i in range(self.getNumEquations()): 451 for j in range(self.getDim()): 452 for k in range(self.getNumSolutions()): 453 for l in range(self.getDim()): 454 if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol: 455 if verbose: print "non-symmetric PDE because A[%d,%d,%d,%d]!=A[%d,%d,%d,%d]"%(i,j,k,l,k,l,i,j) 456 out=False 457 else: 458 for j in range(self.getDim()): 459 for l in range(self.getDim()): 460 if util.Lsup(A[j,l]-A[l,j])>tol: 461 if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j) 462 out=False 463 B=self.getCoefficientOfGeneralPDE("B") 464 C=self.getCoefficientOfGeneralPDE("C") 465 if B.isEmpty() and not C.isEmpty(): 466 if verbose: print "non-symmetric PDE because B is not present but C is" 467 out=False 468 elif not B.isEmpty() and C.isEmpty(): 469 if verbose: print "non-symmetric PDE because C is not present but B is" 470 out=False 471 elif not B.isEmpty() and not C.isEmpty(): 472 tol=(util.Lsup(B)+util.Lsup(C))*self.TOL/2. 473 if self.getNumSolutions()>1: 474 for i in range(self.getNumEquations()): 475 for j in range(self.getDim()): 476 for k in range(self.getNumSolutions()): 477 if util.Lsup(B[i,j,k]-C[k,i,j])>tol: 478 if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j) 479 out=False 480 else: 481 for j in range(self.getDim()): 482 if util.Lsup(B[j]-C[j])>tol: 483 if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j) 484 out=False 485 if self.getNumSolutions()>1: 486 D=self.getCoefficientOfGeneralPDE("D") 487 if not D.isEmpty(): 488 tol=util.Lsup(D)*self.TOL 489 for i in range(self.getNumEquations()): 490 for k in range(self.getNumSolutions()): 491 if util.Lsup(D[i,k]-D[k,i])>tol: 492 if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i) 493 out=False 494 495 return out 496 497 def getSolution(self,**options): 498 """ 499 returns the solution of the PDE. If the solution is not valid the PDE is solved. 500 501 @return: the solution 502 @rtype: L{escript.Data} 503 @param options: solver options 504 @keyword verbose: 505 @keyword reordering: reordering scheme to be used during elimination 506 @keyword preconditioner: preconditioner method to be used 507 @keyword iter_max: maximum number of iteration steps allowed. 508 @keyword drop_tolerance: 509 @keyword drop_storage: 510 @keyword truncation: 511 @keyword restart: 512 """ 513 if not self.__solution_isValid: 514 mat,f=self.getSystem() 515 if self.isUsingLumping(): 516 self.__solution=self.copyConstraint(f*mat) 517 else: 518 options[self.TOLERANCE_KEY]=self.getTolerance() 519 options[self.METHOD_KEY]=self.getSolverMethod() 520 options[self.SYMMETRY_KEY]=self.isSymmetric() 521 self.trace("PDE is resolved.") 522 self.trace("solver options: %s"%str(options)) 523 self.__solution=mat.solve(f,options) 524 self.__solution_isValid=True 525 return self.__solution 526 527 def getFlux(self,u=None): 528 """ 529 returns the flux J_ij for a given u 530 531 \f[ 532 J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij} 533 \f] 534 535 @param u: argument in the flux. If u is not present or equals L{None} the current solution is used. 536 @type u: L{escript.Data} or None 537 @return : flux 538 @rtype : L{escript.Data} 539 540 """ 541 if u==None: u=self.getSolution() 542 return util.tensormult(self.getCoefficientOfGeneralPDE("A"),util.grad(u))+util.matrixmult(self.getCoefficientOfGeneralPDE("B"),u)-util.self.getCoefficientOfGeneralPDE("X") 543 544 # ============================================================================= 545 # solver settings: 546 # ============================================================================= 547 def setSolverMethod(self,solver=None): 548 """ 549 sets a new solver 550 551 @param solver: sets a new solver method. 552 @type solver: C{int} 553 554 """ 555 if solver==None: solve=self.DEFAULT_METHOD 556 if not solver==self.getSolverMethod(): 557 self.__solver_method=solver 558 self.__checkMatrixType() 559 self.trace("New solver is %s"%self.getSolverMethodName()) 560 561 def getSolverMethodName(self): 562 """ 563 returns the name of the solver currently used 564 565 @return : the name of the solver currently used. 566 @rtype: C{string} 567 """ 568 569 m=self.getSolverMethod() 570 if m==self.DEFAULT_METHOD: return "DEFAULT_METHOD" 571 elif m==self.DIRECT: return "DIRECT" 572 elif m==self.CHOLEVSKY: return "CHOLEVSKY" 573 elif m==self.PCG: return "PCG" 574 elif m==self.CR: return "CR" 575 elif m==self.CGS: return "CGS" 576 elif m==self.BICGSTAB: return "BICGSTAB" 577 elif m==self.SSOR: return "SSOR" 578 elif m==self.GMRES: return "GMRES" 579 elif m==self.PRES20: return "PRES20" 580 elif m==self.LUMPING: return "LUMPING" 581 return None 582 583 584 def getSolverMethod(self): 585 """ 586 returns the solver method 587 588 @return : the solver method currently be used. 589 @rtype : C{int} 590 """ 591 return self.__solver_method 592 593 def isUsingLumping(self): 594 """ 595 checks if matrix lumping is used a solver method 596 597 @return : True is lumping is currently used a solver method. 598 @rtype: C{bool} 599 """ 600 return self.getSolverMethod()==self.LUMPING 601 602 def setTolerance(self,tol=1.e-8): 603 """ 604 resets the tolerance for the solver method to tol where for an appropriate norm |.| 605 606 |self.getResidual()|0: 616 raise ValueException,"Tolerance as to be positive" 617 if tol1: 848 return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True) 849 else: 850 return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True) 851 852 def __getNewSolution(self): 853 """ 854 returns an instance of a new solution 855 """ 856 self.trace("New solution is allocated.") 857 if self.getNumSolutions()>1: 858 return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True) 859 else: 860 return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True) 861 862 def __makeFreshSolution(self): 863 """ 864 makes sure that the solution is instantiated and returns it initialized by zeros 865 """ 866 if self.__solution.isEmpty(): 867 self.__solution=self.__getNewSolution() 868 else: 869 self.__solution*=0 870 self.trace("Solution is reset to zero.") 871 return self.__solution 872 873 def __makeFreshRightHandSide(self): 874 """ 875 makes sure that the right hand side is instantiated and returns it initialized by zeros 876 """ 877 if self.__righthandside.isEmpty(): 878 self.__righthandside=self.__getNewRightHandSide() 879 else: 880 self.__righthandside*=0 881 self.trace("Right hand side is reset to zero.") 882 return self.__righthandside 883 884 def __makeFreshOperator(self): 885 """ 886 makes sure that the operator is instantiated and returns it initialized by zeros 887 """ 888 if self.__operator.isEmpty(): 889 self.__operator=self.__getNewOperator() 890 else: 891 self.__operator.setValue(0.) 892 self.trace("Operator reset to zero") 893 return self.__operator 894 895 def __applyConstraint(self): 896 """ 897 applies the constraints defined by q and r to the system 898 """ 899 if not self.isUsingLumping(): 900 q=self.getCoefficientOfGeneralPDE("q") 901 r=self.getCoefficientOfGeneralPDE("r") 902 if not q.isEmpty() and not self.__operator.isEmpty(): 903 # q is the row and column mask to indicate where constraints are set: 904 row_q=escript.Data(q,self.getFunctionSpaceForEquation()) 905 col_q=escript.Data(q,self.getFunctionSpaceForSolution()) 906 u=self.__getNewSolution() 907 if r.isEmpty(): 908 r_s=self.__getNewSolution() 909 else: 910 r_s=escript.Data(r,self.getFunctionSpaceForSolution()) 911 u.copyWithMask(r_s,col_q) 912 if not self.__righthandside.isEmpty(): 913 self.__righthandside-=self.__operator*u 914 self.__righthandside=self.copyConstraint(self.__righthandside) 915 self.__operator.nullifyRowsAndCols(row_q,col_q,1.) 916 # ============================================================================= 917 # function giving access to coefficients of the general PDE: 918 # ============================================================================= 919 def getCoefficientOfGeneralPDE(self,name): 920 """ 921 return the value of the coefficient name of the general PDE. 922 923 @note This method is called by the assembling routine it can be overwritten 924 to map coefficients of a particular PDE to the general PDE. 925 926 @param name: name of the coefficient requested. 927 @type name: C{string} 928 @return : the value of the coefficient name 929 @rtype : L{escript.Data} 930 @raise IllegalCoefficient: if name is not one of coefficients 931 "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q". 932 """ 933 if self.hasCoefficientOfGeneralPDE(name): 934 return self.getCoefficient(name) 935 else: 936 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name 937 938 def hasCoefficientOfGeneralPDE(self,name): 939 """ 940 checks if name is a the name of a coefficient of the general PDE. 941 942 @param name: name of the coefficient enquired. 943 @type name: C{string} 944 @return : True if name is the name of a coefficient of the general PDE. Otherwise False. 945 @rtype : C{bool} 946 947 """ 948 return self.__COEFFICIENTS_OF_GENEARL_PDE.has_key(name) 949 950 def createCoefficientOfGeneralPDE(self,name): 951 """ 952 returns a new instance of a coefficient for coefficient name of the general PDE 953 954 @param name: name of the coefficient requested. 955 @type name: C{string} 956 @return : a coefficient name initialized to 0. 957 @rtype : L{escript.Data} 958 @raise IllegalCoefficient: if name is not one of coefficients 959 "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q". 960 """ 961 if self.hasCoefficientOfGeneralPDE(name): 962 return escript.Data(0,self.getShapeOfCoefficientOfGeneralPDE(name),self.getFunctionSpaceForCoefficientOfGeneralPDE(name)) 963 else: 964 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name 965 966 def getFunctionSpaceForCoefficientOfGeneralPDE(self,name): 967 """ 968 return the L{escript.FunctionSpace} to be used for coefficient name of the general PDE 969 970 @param name: name of the coefficient enquired. 971 @type name: C{string} 972 @return : the function space to be used for coefficient name 973 @rtype : L{escript.FunctionSpace} 974 @raise IllegalCoefficient: if name is not one of coefficients 975 "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q". 976 """ 977 if self.hasCoefficientOfGeneralPDE(name): 978 return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getFunctionSpace(self.getDomain()) 979 else: 980 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name 981 982 def getShapeOfCoefficientOfGeneralPDE(self,name): 983 """ 984 return the shape of the coefficient name of the general PDE 985 986 @param name: name of the coefficient enquired. 987 @type name: C{string} 988 @return : the shape of the coefficient name 989 @rtype : C{tuple} of C{int} 990 @raise IllegalCoefficient: if name is not one of coefficients 991 "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q". 992 993 """ 994 if self.hasCoefficientOfGeneralPDE(name): 995 return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions()) 996 else: 997 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name 998 999 # ============================================================================= 1000 # functions giving access to coefficients of a particular PDE implementation: 1001 # ============================================================================= 1002 def getCoefficient(self,name): 1003 """ 1004 returns the value of the coefficient name 1005 1006 @param name: name of the coefficient requested. 1007 @type name: C{string} 1008 @return : the value of the coefficient name 1009 @rtype : L{escript.Data} 1010 @raise IllegalCoefficient: if name is not a coefficient of the PDE. 1011 """ 1012 if self.hasCoefficient(name): 1013 return self.COEFFICIENTS[name].getValue() 1014 else: 1015 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name 1016 1017 def hasCoefficient(self,name): 1018 """ 1019 return True if name is the name of a coefficient 1020 1021 @param name: name of the coefficient enquired. 1022 @type name: C{string} 1023 @return : True if name is the name of a coefficient of the general PDE. Otherwise False. 1024 @rtype : C{bool} 1025 """ 1026 return self.COEFFICIENTS.has_key(name) 1027 1028 def createCoefficient(self, name): 1029 """ 1030 create a L{escript.Data} object corresponding to coefficient name 1031 1032 @return : a coefficient name initialized to 0. 1033 @rtype : L{escript.Data} 1034 @raise IllegalCoefficient: if name is not a coefficient of the PDE. 1035 """ 1036 if self.hasCoefficient(name): 1037 return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name)) 1038 else: 1039 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name 1040 1041 def getFunctionSpaceForCoefficient(self,name): 1042 """ 1043 return the L{escript.FunctionSpace} to be used for coefficient name 1044 1045 @param name: name of the coefficient enquired. 1046 @type name: C{string} 1047 @return : the function space to be used for coefficient name 1048 @rtype : L{escript.FunctionSpace} 1049 @raise IllegalCoefficient: if name is not a coefficient of the PDE. 1050 """ 1051 if self.hasCoefficient(name): 1052 return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain()) 1053 else: 1054 raise ValueError,"unknown coefficient %s requested"%name 1055 1056 def getShapeOfCoefficient(self,name): 1057 """ 1058 return the shape of the coefficient name 1059 1060 @param name: name of the coefficient enquired. 1061 @type name: C{string} 1062 @return : the shape of the coefficient name 1063 @rtype : C{tuple} of C{int} 1064 @raise IllegalCoefficient: if name is not a coefficient of the PDE. 1065 """ 1066 if self.hasCoefficient(name): 1067 return self.COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions()) 1068 else: 1069 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name 1070 1071 def resetCoefficients(self): 1072 """ 1073 resets all coefficients to there default values. 1074 """ 1075 for i in self.COEFFICIENTS.iterkeys(): 1076 self.COEFFICIENTS[i].resetValue() 1077 1078 def alteredCoefficient(self,name): 1079 """ 1080 announce that coefficient name has been changed 1081 1082 @param name: name of the coefficient enquired. 1083 @type name: C{string} 1084 @raise IllegalCoefficient: if name is not a coefficient of the PDE. 1085 """ 1086 if self.hasCoefficient(name): 1087 self.trace("Coefficient %s has been altered."%name) 1088 if self.COEFFICIENTS[name].isAlteringOperator(): self.__invalidateOperator() 1089 if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__invalidateRightHandSide() 1090 else: 1091 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name 1092 1093 def copyConstraint(self,u): 1094 """ 1095 copies the constraint into u and returns u. 1096 1097 @param u: a function of rank 0 is a single PDE is solved and of shape (numSolution,) for a system of PDEs 1098 @type u: L{escript.Data} 1099 @return : the input u modified by the constraints. 1100 @rtype : L{escript.Data} 1101 @warning: u is altered if it has the appropriate L{escript.FunctionSpace} 1102 1103 """ 1104 q=self.getCoefficientOfGeneralPDE("q") 1105 r=self.getCoefficientOfGeneralPDE("r") 1106 if not q.isEmpty(): 1107 if u.isEmpty(): u=escript.Data(0.,q.getShape(),q.getFunctionSpace()) 1108 if r.isEmpty(): 1109 r=escript.Data(0,u.getShape(),u.getFunctionSpace()) 1110 else: 1111 r=escript.Data(r,u.getFunctionSpace()) 1112 u.copyWithMask(r,escript.Data(q,u.getFunctionSpace())) 1113 return u 1114 1115 def setValue(self,**coefficients): 1116 """ 1117 sets new values to coefficients 1118 1119 @note This method is called by the assembling routine it can be overwritten 1120 to map coefficients of a particular PDE to the general PDE. 1121 1122 @param name: name of the coefficient requested. 1123 @type name: C{string} 1124 @keyword A: value for coefficient A. 1125 @type A: any type that can be interpreted as L{escript.Data} object on L{escript.Function}. 1126 @keyword B: value for coefficient B 1127 @type B: any type that can be interpreted as L{escript.Data} object on L{escript.Function}. 1128 @keyword C: value for coefficient C 1129 @type C: any type that can be interpreted as L{escript.Data} object on L{escript.Function}. 1130 @keyword D: value for coefficient D 1131 @type D: any type that can be interpreted as L{escript.Data} object on L{escript.Function}. 1132 @keyword X: value for coefficient X 1133 @type X: any type that can be interpreted as L{escript.Data} object on L{escript.Function}. 1134 @keyword Y: value for coefficient Y 1135 @type Y: any type that can be interpreted as L{escript.Data} object on L{escript.Function}. 1136 @keyword d: value for coefficient d 1137 @type d: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnBoundary}. 1138 @keyword y: value for coefficient y 1139 @type y: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnBoundary}. 1140 @keyword d_contact: value for coefficient d_contact 1141 @type d_contact: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnContactOne}. 1142 or L{escript.FunctionOnContactZero}. 1143 @keyword y_contact: value for coefficient y_contact 1144 @type y_contact: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnContactOne}. 1145 or L{escript.FunctionOnContactZero}. 1146 @keyword r: values prescribed to the solution at the locations of constraints 1147 @type r: any type that can be interpreted as L{escript.Data} object on L{escript.Solution} or L{escript.ReducedSolution} 1148 depending of reduced order is used for the solution. 1149 @keyword q: mask for location of constraints 1150 @type q: any type that can be interpreted as L{escript.Data} object on L{escript.Solution} or L{escript.ReducedSolution} 1151 depending of reduced order is used for the representation of the equation. 1152 @raise IllegalCoefficient: if an unknown coefficient keyword is used. 1153 1154 """ 1155 # check if the coefficients are legal: 1156 for i in coefficients.iterkeys(): 1157 if not self.hasCoefficient(i): 1158 raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i 1159 # if the number of unknowns or equations is still unknown we try to estimate them: 1160 if self.__numEquations==None or self.__numSolutions==None: 1161 for i,d in coefficients.iteritems(): 1162 if hasattr(d,"shape"): 1163 s=d.shape 1164 elif hasattr(d,"getShape"): 1165 s=d.getShape() 1166 else: 1167 s=numarray.array(d).shape 1168 if s!=None: 1169 # get number of equations and number of unknowns: 1170 res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s) 1171 if res==None: 1172 raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i) 1173 else: 1174 if self.__numEquations==None: self.__numEquations=res[0] 1175 if self.__numSolutions==None: self.__numSolutions=res[1] 1176 if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations" 1177 if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions" 1178 # now we check the shape of the coefficient if numEquations and numSolutions are set: 1179 for i,d in coefficients.iteritems(): 1180 try: 1181 self.COEFFICIENTS[i].setValue(self.getDomain(),self.getNumEquations(),self.getNumSolutions(),d) 1182 except IllegalCoefficientValue,m: 1183 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m)) 1184 self.alteredCoefficient(i) 1185 1186 # check if the systrem is inhomogeneous: 1187 if len(coefficients)>0 and not self.isUsingLumping(): 1188 q=self.getCoefficientOfGeneralPDE("q") 1189 r=self.getCoefficientOfGeneralPDE("r") 1190 homogeneous_constraint=True 1191 if not q.isEmpty() and not r.isEmpty(): 1192 if util.Lsup(q*r)>=1.e-13*util.Lsup(r): 1193 self.trace("Inhomogeneous constraint detected.") 1194 self.__invalidateSystem() 1195 1196 1197 def getSystem(self): 1198 """ 1199 return the operator and right hand side of the PDE 1200 """ 1201 if not self.__operator_isValid or not self.__righthandside_isValid: 1202 if self.isUsingLumping(): 1203 if not self.__operator_isValid: 1204 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution(): raise TypeError,"Lumped matrix requires same order for equations and unknowns" 1205 if not self.getCoefficientOfGeneralPDE("A").isEmpty(): raise Warning,"Lumped matrix does not allow coefficient A" 1206 if not self.getCoefficientOfGeneralPDE("B").isEmpty(): raise Warning,"Lumped matrix does not allow coefficient B" 1207 if not self.getCoefficientOfGeneralPDE("C").isEmpty(): raise Warning,"Lumped matrix does not allow coefficient C" 1208 mat=self.__getNewOperator() 1209 self.getDomain().addPDEToSystem(mat,escript.Data(), \ 1210 self.getCoefficientOfGeneralPDE("A"), \ 1211 self.getCoefficientOfGeneralPDE("B"), \ 1212 self.getCoefficientOfGeneralPDE("C"), \ 1213 self.getCoefficientOfGeneralPDE("D"), \ 1214 escript.Data(), \ 1215 escript.Data(), \ 1216 self.getCoefficientOfGeneralPDE("d"), \ 1217 escript.Data(),\ 1218 self.getCoefficientOfGeneralPDE("d_contact"), \ 1219 escript.Data()) 1220 self.__operator=1./(mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)) 1221 del mat 1222 self.trace("New lumped operator has been built.") 1223 self.__operator_isValid=True 1224 if not self.__righthandside_isValid: 1225 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \ 1226 self.getCoefficientOfGeneralPDE("X"), \ 1227 self.getCoefficientOfGeneralPDE("Y"),\ 1228 self.getCoefficientOfGeneralPDE("y"),\ 1229 self.getCoefficientOfGeneralPDE("y_contact")) 1230 self.trace("New right hand side as been built.") 1231 self.__righthandside_isValid=True 1232 else: 1233 if not self.__operator_isValid and not self.__righthandside_isValid: 1234 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \ 1235 self.getCoefficientOfGeneralPDE("A"), \ 1236 self.getCoefficientOfGeneralPDE("B"), \ 1237 self.getCoefficientOfGeneralPDE("C"), \ 1238 self.getCoefficientOfGeneralPDE("D"), \ 1239 self.getCoefficientOfGeneralPDE("X"), \ 1240 self.getCoefficientOfGeneralPDE("Y"), \ 1241 self.getCoefficientOfGeneralPDE("d"), \ 1242 self.getCoefficientOfGeneralPDE("y"), \ 1243 self.getCoefficientOfGeneralPDE("d_contact"), \ 1244 self.getCoefficientOfGeneralPDE("y_contact")) 1245 self.__applyConstraint() 1246 self.__righthandside=self.copyConstraint(self.__righthandside) 1247 self.trace("New system has been built.") 1248 self.__operator_isValid=True 1249 self.__righthandside_isValid=True 1250 elif not self.__righthandside_isValid: 1251 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \ 1252 self.getCoefficientOfGeneralPDE("X"), \ 1253 self.getCoefficientOfGeneralPDE("Y"),\ 1254 self.getCoefficientOfGeneralPDE("y"),\ 1255 self.getCoefficientOfGeneralPDE("y_contact")) 1256 self.__righthandside=self.copyConstraint(self.__righthandside) 1257 self.trace("New right hand side has been built.") 1258 self.__righthandside_isValid=True 1259 elif not self.__operator_isValid: 1260 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \ 1261 self.getCoefficientOfGeneralPDE("A"), \ 1262 self.getCoefficientOfGeneralPDE("B"), \ 1263 self.getCoefficientOfGeneralPDE("C"), \ 1264 self.getCoefficientOfGeneralPDE("D"), \ 1265 escript.Data(), \ 1266 escript.Data(), \ 1267 self.getCoefficientOfGeneralPDE("d"), \ 1268 escript.Data(),\ 1269 self.getCoefficientOfGeneralPDE("d_contact"), \ 1270 escript.Data()) 1271 self.__applyConstraint() 1272 self.trace("New operator has been built.") 1273 self.__operator_isValid=True 1274 return (self.__operator,self.__righthandside) 1275 1276 1277 1278 1279 class AdvectivePDE(LinearPDE): 1280 """ 1281 Class to handle a linear PDE dominated by advective terms: 1282 1283 class to define a linear PDE of the form 1284 1285 \f[ 1286 -(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i 1287 \f] 1288 1289 with boundary conditons: 1290 1291 \f[ 1292 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i 1293 \f] 1294 1295 and contact conditions 1296 1297 \f[ 1298 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d^{contact}_{ik}[u_k] = - n_j*X_{ij} + y^{contact}_{i} 1299 \f] 1300 1301 and constraints: 1302 1303 \f[ 1304 u_i=r_i \quad \mathrm{where} \quad q_i>0 1305 \f] 1306 """ 1307 def __init__(self,domain,numEquations=0,numSolutions=0,xi=None,debug=False): 1308 LinearPDE.__init__(self,domain,numEquations,numSolutions,debug) 1309 if xi==None: 1310 self.__xi=AdvectivePDE.ELMAN_RAMAGE 1311 else: 1312 self.__xi=xi 1313 self.__Xi=escript.Data() 1314 1315 def __calculateXi(self,peclet_factor,Z,h): 1316 Z_max=util.Lsup(Z) 1317 if Z_max>0.: 1318 return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.TOL) 1319 else: 1320 return 0. 1321 1322 def setValue(self,**args): 1323 if "A" in args.keys() or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data() 1324 LinearPDE.setValue(**args) 1325 1326 def ELMAN_RAMAGE(P): 1327 """ """ 1328 return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15)) 1329 def SIMPLIFIED_BROOK_HUGHES(P): 1330 """ """ 1331 c=(P-3.).whereNegative() 1332 return P/6.*c+1./2.*(1.-c) 1333 1334 def HALF(P): 1335 """ """ 1336 return escript.Scalar(0.5,P.getFunctionSpace()) 1337 1338 def getXi(self): 1339 if self.__Xi.isEmpty(): 1340 B=self.getCoefficient("B") 1341 C=self.getCoefficient("C") 1342 A=self.getCoefficient("A") 1343 h=self.getDomain().getSize() 1344 self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A")) 1345 if not C.isEmpty() or not B.isEmpty(): 1346 if not C.isEmpty() and not B.isEmpty(): 1347 Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A")) 1348 if self.getNumEquations()>1: 1349 if self.getNumSolutions()>1: 1350 for i in range(self.getNumEquations()): 1351 for k in range(self.getNumSolutions()): 1352 for l in range(self.getDim()): Z2+=(C[i,k,l]-B[i,l,k])**2 1353 else: 1354 for i in range(self.getNumEquations()): 1355 for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2 1356 else: 1357 if self.getNumSolutions()>1: 1358 for k in range(self.getNumSolutions()): 1359 for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2 1360 else: 1361 for l in range(self.getDim()): Z2+=(C[l]-B[l])**2 1362 length_of_Z=util.sqrt(Z2) 1363 elif C.isEmpty(): 1364 length_of_Z=util.length(B) 1365 else: 1366 length_of_Z=util.length(C) 1367 1368 Z_max=util.Lsup(length_of_Z) 1369 if Z_max>0.: 1370 length_of_A=util.length(A) 1371 A_max=util.Lsup(length_of_A) 1372 if A_max>0: 1373 inv_A=1./(length_of_A+A_max*self.TOL) 1374 else: 1375 inv_A=1./self.TOL 1376 peclet_number=length_of_Z*h/2*inv_A 1377 xi=self.__xi(peclet_number) 1378 self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL) 1379 print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z) 1380 return self.__Xi 1381 1382 1383 def getCoefficientOfGeneralPDE(self,name): 1384 """ 1385 return the value of the coefficient name of the general PDE 1386 1387 @param name: 1388 """ 1389 if not self.getNumEquations() == self.getNumSolutions(): 1390 raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal." 1391 1392 if name == "A" : 1393 A=self.getCoefficient("A") 1394 B=self.getCoefficient("B") 1395 C=self.getCoefficient("C") 1396 if B.isEmpty() and C.isEmpty(): 1397 Aout=A 1398 else: 1399 if A.isEmpty(): 1400 Aout=self.createNewCoefficient("A") 1401 else: 1402 Aout=A[:] 1403 Xi=self.getXi() 1404 if self.getNumEquations()>1: 1405 for i in range(self.getNumEquations()): 1406 for j in range(self.getDim()): 1407 for k in range(self.getNumSolutions()): 1408 for l in range(self.getDim()): 1409 if not C.isEmpty() and not B.isEmpty(): 1410 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*(C[p,i,j]-B[p,j,i])*(C[p,k,l]-B[p,l,k]) 1411 elif C.isEmpty(): 1412 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k] 1413 else: 1414 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l] 1415 else: 1416 for j in range(self.getDim()): 1417 for l in range(self.getDim()): 1418 if not C.isEmpty() and not B.isEmpty(): 1419 Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l]) 1420 elif C.isEmpty(): 1421 Aout[j,l]+=Xi*B[j]*B[l] 1422 else: 1423 Aout[j,l]+=Xi*C[j]*C[l] 1424 return Aout 1425 elif name == "B" : 1426 B=self.getCoefficient("B") 1427 C=self.getCoefficient("C") 1428 D=self.getCoefficient("D") 1429 if C.isEmpty() or D.isEmpty(): 1430 Bout=B 1431 else: 1432 Xi=self.getXi() 1433 if B.isEmpty(): 1434 Bout=self.createNewCoefficient("B") 1435 else: 1436 Bout=B[:] 1437 if self.getNumEquations()>1: 1438 for k in range(self.getNumSolutions()): 1439 for p in range(self.getNumEquations()): 1440 tmp=Xi*D[p,k] 1441 for i in range(self.getNumEquations()): 1442 for j in range(self.getDim()): 1443 Bout[i,j,k]+=tmp*C[p,i,j] 1444 else: 1445 tmp=Xi*D 1446 for j in range(self.getDim()): Bout[j]+=tmp*C[j] 1447 return Bout 1448 elif name == "C" : 1449 B=self.getCoefficient("B") 1450 C=self.getCoefficient("C") 1451 D=self.getCoefficient("D") 1452 if B.isEmpty() or D.isEmpty(): 1453 Cout=C 1454 else: 1455 Xi=self.getXi() 1456 if C.isEmpty(): 1457 Cout=self.createNewCoefficient("C") 1458 else: 1459 Cout=C[:] 1460 if self.getNumEquations()>1: 1461 for k in range(self.getNumSolutions()): 1462 for p in range(self.getNumEquations()): 1463 tmp=Xi*D[p,k] 1464 for i in range(self.getNumEquations()): 1465 for l in range(self.getDim()): 1466 Cout[i,k,l]+=tmp*B[p,l,i] 1467 else: 1468 tmp=Xi*D 1469 for j in range(self.getDim()): Cout[j]+=tmp*B[j] 1470 return Cout 1471 elif name == "D" : 1472 return self.getCoefficient("D") 1473 elif name == "X" : 1474 X=self.getCoefficient("X") 1475 Y=self.getCoefficient("Y") 1476 B=self.getCoefficient("B") 1477 C=self.getCoefficient("C") 1478 if Y.isEmpty() or (B.isEmpty() and C.isEmpty()): 1479 Xout=X 1480 else: 1481 if X.isEmpty(): 1482 Xout=self.createNewCoefficient("X") 1483 else: 1484 Xout=X[:] 1485 Xi=self.getXi() 1486 if self.getNumEquations()>1: 1487 for p in range(self.getNumEquations()): 1488 tmp=Xi*Y[p] 1489 for i in range(self.getNumEquations()): 1490 for j in range(self.getDim()): 1491 if not C.isEmpty() and not B.isEmpty(): 1492 Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i]) 1493 elif C.isEmpty(): 1494 Xout[i,j]-=tmp*B[p,j,i] 1495 else: 1496 Xout[i,j]+=tmp*C[p,i,j] 1497 else: 1498 tmp=Xi*Y 1499 for j in range(self.getDim()): 1500 if not C.isEmpty() and not B.isEmpty(): 1501 Xout[j]+=tmp*(C[j]-B[j]) 1502 elif C.isEmpty(): 1503 Xout[j]-=tmp*B[j] 1504 else: 1505 Xout[j]+=tmp*C[j] 1506 return Xout 1507 elif name == "Y" : 1508 return self.getCoefficient("Y") 1509 elif name == "d" : 1510 return self.getCoefficient("d") 1511 elif name == "y" : 1512 return self.getCoefficient("y") 1513 elif name == "d_contact" : 1514 return self.getCoefficient("d_contact") 1515 elif name == "y_contact" : 1516 return self.getCoefficient("y_contact") 1517 elif name == "r" : 1518 return self.getCoefficient("r") 1519 elif name == "q" : 1520 return self.getCoefficient("q") 1521 else: 1522 raise SystemError,"unknown PDE coefficient %s",name 1523 1524 1525 class Poisson(LinearPDE): 1526 """ 1527 Class to define a Poisson equation problem: 1528 1529 class to define a linear PDE of the form 1530 \f[ 1531 -u_{,jj} = f 1532 \f] 1533 1534 with boundary conditons: 1535 1536 \f[ 1537 n_j*u_{,j} = 0 1538 \f] 1539 1540 and constraints: 1541 1542 \f[ 1543 u=0 \quad \mathrm{where} \quad q>0 1544 \f] 1545 """ 1546 1547 def __init__(self,domain,f=escript.Data(),q=escript.Data(),debug=False): 1548 LinearPDE.__init__(self,domain,1,1,debug) 1549 self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE), 1550 "q": PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)} 1551 self.setSymmetryOn() 1552 self.setValue(f,q) 1553 1554 def setValue(self,f=escript.Data(),q=escript.Data()): 1555 """set value of PDE parameters f and q""" 1556 self._LinearPDE__setValue(f=f,q=q) 1557 1558 def getCoefficientOfGeneralPDE(self,name): 1559 """ 1560 return the value of the coefficient name of the general PDE 1561 1562 @param name: 1563 """ 1564 if name == "A" : 1565 return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain())) 1566 elif name == "B" : 1567 return escript.Data() 1568 elif name == "C" : 1569 return escript.Data() 1570 elif name == "D" : 1571 return escript.Data() 1572 elif name == "X" : 1573 return escript.Data() 1574 elif name == "Y" : 1575 return self.getCoefficient("f") 1576 elif name == "d" : 1577 return escript.Data() 1578 elif name == "y" : 1579 return escript.Data() 1580 elif name == "d_contact" : 1581 return escript.Data() 1582 elif name == "y_contact" : 1583 return escript.Data() 1584 elif name == "r" : 1585 return escript.Data() 1586 elif name == "q" : 1587 return self.getCoefficient("q") 1588 else: 1589 raise SystemError,"unknown PDE coefficient %s",name 1590 1591 class LameEquation(LinearPDE): 1592 """ 1593 Class to define a Lame equation problem: 1594 1595 class to define a linear PDE of the form 1596 \f[ 1597 -(\mu (u_{i,j}+u_{j,i}))_{,j} - \lambda u_{j,ji}} = F_i -\sigma_{ij,j} 1598 \f] 1599 1600 with boundary conditons: 1601 1602 \f[ 1603 n_j(\mu (u_{i,j}+u_{j,i})-sigma_{ij}) + n_i\lambda u_{j,j} = f_i 1604 \f] 1605 1606 and constraints: 1607 1608 \f[ 1609 u_i=r_i \quad \mathrm{where} \quad q_i>0 1610 \f] 1611 """ 1612 1613 def __init__(self,domain,f=escript.Data(),q=escript.Data(),debug=False): 1614 LinearPDE.__init__(self,domain,domain.getDim(),domain.getDim(),debug) 1615 self.COEFFICIENTS={ "lame_lambda" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR), 1616 "lame_mu" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR), 1617 "F" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE), 1618 "sigma" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE), 1619 "f" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE), 1620 "r" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH), 1621 "q" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)} 1622 self.setSymmetryOn() 1623 1624 def setValue(self,lame_lambda=escript.Data(),lame_mu=escript.Data(),F=escript.Data(),sigma=escript.Data(),f=escript.Data(),r=escript.Data(),q=escript.Data()): 1625 """set value of PDE parameters""" 1626 self._LinearPDE__setValue(lame_lambda=lame_lambda, \ 1627 lame_mu=lame_mu, \ 1628 F=F, \ 1629 sigma=sigma, \ 1630 f=f, \ 1631 r=r, \ 1632 q=q) 1633 def getCoefficientOfGeneralPDE(self,name): 1634 """ 1635 return the value of the coefficient name of the general PDE 1636 1637 @param name: 1638 """ 1639 if name == "A" : 1640 out =self.createCoefficientOfGeneralPDE("A") 1641 for i in range(self.getDim()): 1642 for j in range(self.getDim()): 1643 out[i,i,j,j] += self.getCoefficient("lame_lambda") 1644 out[i,j,j,i] += self.getCoefficient("lame_mu") 1645 out[i,j,i,j] += self.getCoefficient("lame_mu") 1646 return out 1647 elif name == "B" : 1648 return escript.Data() 1649 elif name == "C" : 1650 return escript.Data() 1651 elif name == "D" : 1652 return escript.Data() 1653 elif name == "X" : 1654 return self.getCoefficient("sigma") 1655 elif name == "Y" : 1656 return self.getCoefficient("F") 1657 elif name == "d" : 1658 return escript.Data() 1659 elif name == "y" : 1660 return self.getCoefficient("f") 1661 elif name == "d_contact" : 1662 return escript.Data() 1663 elif name == "y_contact" : 1664 return escript.Data() 1665 elif name == "r" : 1666 return self.getCoefficient("r") 1667 elif name == "q" : 1668 return self.getCoefficient("q") 1669 else: 1670 raise SystemError,"unknown PDE coefficient %s",name 1671 1672 # $Log$ 1673 # Revision 1.11 2005/08/23 01:24:28 jgs 1674 # Merge of development branch dev-02 back to main trunk on 2005-08-23 1675 # 1676 # Revision 1.10 2005/08/12 01:45:36 jgs 1677 # erge of development branch dev-02 back to main trunk on 2005-08-12 1678 # 1679 # Revision 1.9.2.4 2005/08/22 07:11:09 gross 1680 # some problems with LinearPDEs fixed. 1681 # 1682 # Revision 1.9.2.3 2005/08/18 04:48:48 gross 1683 # the methods SetLumping*() are removed. Lumping is set trough setSolverMethod(LinearPDE.LUMPING) 1684 # 1685 # Revision 1.9.2.2 2005/08/18 04:39:32 gross 1686 # the constants have been removed from util.py as they not needed anymore. PDE related constants are accessed through LinearPDE attributes now 1687 # 1688 # Revision 1.9.2.1 2005/07/29 07:10:27 gross 1689 # new functions in util and a new pde type in linearPDEs 1690 # 1691 # Revision 1.1.2.25 2005/07/28 04:21:09 gross 1692 # Lame equation: (linear elastic, isotropic) added 1693 # 1694 # Revision 1.1.2.24 2005/07/22 06:37:11 gross 1695 # some extensions to modellib and linearPDEs 1696 # 1697 # Revision 1.1.2.23 2005/05/13 00:55:20 cochrane 1698 # Fixed up some docstrings. Moved module-level functions to top of file so 1699 # that epydoc and doxygen can pick them up properly. 1700 # 1701 # Revision 1.1.2.22 2005/05/12 11:41:30 gross 1702 # some basic Models have been added 1703 # 1704 # Revision 1.1.2.21 2005/05/12 07:16:12 cochrane 1705 # Moved ELMAN_RAMAGE, SIMPLIFIED_BROOK_HUGHES, and HALF functions to bottom of 1706 # file so that the AdvectivePDE class is picked up by doxygen. Some 1707 # reformatting of docstrings. Addition of code to make equations come out 1708 # as proper LaTeX. 1709 # 1710 # Revision 1.1.2.20 2005/04/15 07:09:08 gross 1711 # some problems with functionspace and linearPDEs fixed. 1712 # 1713 # Revision 1.1.2.19 2005/03/04 05:27:07 gross 1714 # bug in SystemPattern fixed. 1715 # 1716 # Revision 1.1.2.18 2005/02/08 06:16:45 gross 1717 # Bugs in AdvectivePDE fixed, AdvectiveTest is stable but more testing is needed 1718 # 1719 # Revision 1.1.2.17 2005/02/08 05:56:19 gross 1720 # Reference Number handling added 1721 # 1722 # Revision 1.1.2.16 2005/02/07 04:41:28 gross 1723 # some function exposed to python to make mesh merging running 1724 # 1725 # Revision 1.1.2.15 2005/02/03 00:14:44 gross 1726 # timeseries add and ESySParameter.py renames esysXML.py for consistence 1727 # 1728 # Revision 1.1.2.14 2005/02/01 06:44:10 gross 1729 # new implementation of AdvectivePDE which now also updates right hand side. systems of PDEs are still not working 1730 # 1731 # Revision 1.1.2.13 2005/01/25 00:47:07 gross 1732 # updates in the documentation 1733 # 1734 # Revision 1.1.2.12 2005/01/12 01:28:04 matt 1735 # Added createCoefficient method for linearPDEs. 1736 # 1737 # Revision 1.1.2.11 2005/01/11 01:55:34 gross 1738 # a problem in linearPDE class fixed 1739 # 1740 # Revision 1.1.2.10 2005/01/07 01:13:29 gross 1741 # some bugs in linearPDE fixed 1742 # 1743 # Revision 1.1.2.9 2005/01/06 06:24:58 gross 1744 # some bugs in slicing fixed 1745 # 1746 # Revision 1.1.2.8 2005/01/05 04:21:40 gross 1747 # FunctionSpace checking/matchig in slicing added 1748 # 1749 # Revision 1.1.2.7 2004/12/29 10:03:41 gross 1750 # bug in setValue fixed 1751 # 1752 # Revision 1.1.2.6 2004/12/29 05:29:59 gross 1753 # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data() 1754 # 1755 # Revision 1.1.2.5 2004/12/29 00:18:41 gross 1756 # AdvectivePDE added 1757 # 1758 # Revision 1.1.2.4 2004/12/24 06:05:41 gross 1759 # some changes in linearPDEs to add AdevectivePDE 1760 # 1761 # Revision 1.1.2.3 2004/12/16 00:12:34 gross 1762 # __init__ of LinearPDE does not accept any coefficient anymore 1763 # 1764 # Revision 1.1.2.2 2004/12/14 03:55:01 jgs 1765 # *** empty log message *** 1766 # 1767 # Revision 1.1.2.1 2004/12/12 22:53:47 gross 1768 # linearPDE has been renamed LinearPDE 1769 # 1770 # Revision 1.1.1.1.2.7 2004/12/07 10:13:08 gross 1771 # GMRES added 1772 # 1773 # Revision 1.1.1.1.2.6 2004/12/07 03:19:50 gross 1774 # options for GMRES and PRES20 added 1775 # 1776 # Revision 1.1.1.1.2.5 2004/12/01 06:25:15 gross 1777 # some small changes 1778 # 1779 # Revision 1.1.1.1.2.4 2004/11/24 01:50:21 gross 1780 # Finley solves 4M unknowns now 1781 # 1782 # Revision 1.1.1.1.2.3 2004/11/15 06:05:26 gross 1783 # poisson solver added 1784 # 1785 # Revision 1.1.1.1.2.2 2004/11/12 06:58:15 gross 1786 # a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry 1787 # 1788 # Revision 1.1.1.1.2.1 2004/10/28 22:59:22 gross 1789 # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed 1790 # 1791 # Revision 1.1.1.1 2004/10/26 06:53:56 jgs 1792 # initial import of project esys2 1793 # 1794 # Revision 1.3.2.3 2004/10/26 06:43:48 jgs 1795 # committing Lutz's and Paul's changes to brach jgs 1796 # 1797 # Revision 1.3.4.1 2004/10/20 05:32:51 cochrane 1798 # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form. 1799 # 1800 # Revision 1.3 2004/09/23 00:53:23 jgs 1801 # minor fixes 1802 # 1803 # Revision 1.1 2004/08/28 12:58:06 gross 1804 # SimpleSolve is not running yet: problem with == of functionsspace 1805 # 1806 #

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision