# Annotation of /trunk/escript/py_src/linearPDEs.py

Revision 142 - (hide annotations)
Mon Jul 25 05:28:20 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: 54182 byte(s)
Merge of development branch back to main trunk on 2005-07-25


 1 jgs 102 # $Id$ 2 3 ## @file linearPDEs.py 4 5 """ 6 jgs 122 Functions and classes for linear PDEs 7 jgs 102 """ 8 9 import escript 10 import util 11 import numarray 12 13 14 def _CompTuple2(t1,t2): 15 """ 16 jgs 122 Compare two tuples 17 jgs 102 18 jgs 122 \param t1 The first tuple 19 \param t2 The second tuple 20 jgs 102 """ 21 jgs 122 22 jgs 102 dif=t1[0]+t1[1]-(t2[0]+t2[1]) 23 if dif<0: return 1 24 elif dif>0: return -1 25 else: return 0 26 27 jgs 122 def ELMAN_RAMAGE(P): 28 return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15)) 29 30 def SIMPLIFIED_BROOK_HUGHES(P): 31 c=(P-3.).whereNegative() 32 return P/6.*c+1./2.*(1.-c) 33 34 def HALF(P): 35 return escript.Scalar(0.5,P.getFunctionSpace()) 36 37 jgs 108 class PDECoefficient: 38 jgs 102 """ 39 jgs 122 A class for PDE coefficients 40 jgs 102 """ 41 jgs 108 # identifier for location of Data objects defining COEFFICIENTS 42 jgs 102 INTERIOR=0 43 BOUNDARY=1 44 CONTACT=2 45 CONTINUOUS=3 46 jgs 108 # identifier in the pattern of COEFFICIENTS: 47 jgs 102 # the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the 48 # number of unknowns. 49 EQUATION=3 50 SOLUTION=4 51 DIM=5 52 # indicator for what is altered if the coefficient is altered: 53 OPERATOR=5 54 RIGHTHANDSIDE=6 55 BOTH=7 56 def __init__(self,where,pattern,altering): 57 """ 58 jgs 122 Initialise a PDE Coefficient type 59 jgs 102 """ 60 self.what=where 61 self.pattern=pattern 62 self.altering=altering 63 jgs 108 self.resetValue() 64 jgs 102 65 jgs 108 def resetValue(self): 66 """ 67 jgs 122 resets coefficient value to default 68 jgs 108 """ 69 self.value=escript.Data() 70 71 jgs 102 def getFunctionSpace(self,domain): 72 """ 73 jgs 122 defines the FunctionSpace of the coefficient on the domain 74 jgs 102 75 jgs 122 @param domain: 76 jgs 102 """ 77 if self.what==self.INTERIOR: return escript.Function(domain) 78 elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain) 79 elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain) 80 elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain) 81 82 jgs 108 def getValue(self): 83 """ 84 jgs 122 returns the value of the coefficient: 85 jgs 108 """ 86 return self.value 87 88 def setValue(self,newValue): 89 """ 90 jgs 122 set the value of the coefficient to new value 91 jgs 108 """ 92 self.value=newValue 93 94 jgs 102 def isAlteringOperator(self): 95 """ 96 jgs 122 return true if the operator of the PDE is changed when the coefficient is changed 97 jgs 102 """ 98 if self.altering==self.OPERATOR or self.altering==self.BOTH: 99 return not None 100 else: 101 return None 102 103 def isAlteringRightHandSide(self): 104 """ 105 jgs 122 return true if the right hand side of the PDE is changed when the coefficient is changed 106 jgs 102 """ 107 if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH: 108 return not None 109 else: 110 return None 111 112 def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3): 113 """ 114 jgs 122 tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim 115 jgs 102 116 jgs 122 @param shape: 117 @param dim: 118 jgs 102 """ 119 if len(shape)>0: 120 num=max(shape)+1 121 else: 122 num=1 123 search=[] 124 for u in range(num): 125 for e in range(num): 126 search.append((e,u)) 127 search.sort(_CompTuple2) 128 for item in search: 129 s=self.buildShape(item[0],item[1],dim) 130 if len(s)==0 and len(shape)==0: 131 return (1,1) 132 else: 133 if s==shape: return item 134 return None 135 136 def buildShape(self,e=1,u=1,dim=3): 137 """ 138 jgs 122 builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim 139 jgs 102 140 jgs 122 @param e: 141 @param u: 142 @param dim: 143 jgs 102 """ 144 s=() 145 for i in self.pattern: 146 if i==self.EQUATION: 147 if e>1: s=s+(e,) 148 elif i==self.SOLUTION: 149 if u>1: s=s+(u,) 150 else: 151 s=s+(dim,) 152 return s 153 154 class LinearPDE: 155 """ 156 jgs 122 Class to handle a linear PDE 157 jgs 102 158 class to define a linear PDE of the form 159 160 jgs 122 \f[ 161 jgs 102 -(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 162 jgs 122 \f] 163 jgs 102 164 jgs 122 with boundary conditons: 165 jgs 102 166 jgs 122 \f[ 167 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i 168 \f] 169 jgs 102 170 jgs 122 and contact conditions 171 jgs 102 172 jgs 122 \f[ 173 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i 174 \f] 175 jgs 102 176 jgs 122 and constraints: 177 jgs 102 178 jgs 122 \f[ 179 u_i=r_i \quad \mathrm{where} \quad q_i>0 180 \f] 181 jgs 102 182 """ 183 jgs 108 TOL=1.e-13 184 jgs 102 DEFAULT_METHOD=util.DEFAULT_METHOD 185 DIRECT=util.DIRECT 186 CHOLEVSKY=util.CHOLEVSKY 187 PCG=util.PCG 188 CR=util.CR 189 CGS=util.CGS 190 BICGSTAB=util.BICGSTAB 191 SSOR=util.SSOR 192 GMRES=util.GMRES 193 PRES20=util.PRES20 194 jgs 114 ILU0=util.ILU0 195 JACOBI=util.JACOBI 196 jgs 102 197 jgs 104 def __init__(self,domain,numEquations=0,numSolutions=0): 198 jgs 102 """ 199 jgs 122 initializes a new linear PDE. 200 jgs 102 201 jgs 122 @param args: 202 jgs 102 """ 203 jgs 108 # COEFFICIENTS can be overwritten by subclasses: 204 self.COEFFICIENTS={ 205 "A" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR), 206 "B" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR), 207 "C" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR), 208 "D" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR), 209 "X" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE), 210 "Y" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE), 211 "d" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR), 212 "y" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE), 213 "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR), 214 "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE), 215 "r" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE), 216 "q" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)} 217 jgs 102 218 # initialize attributes 219 self.__debug=None 220 jgs 104 self.__domain=domain 221 self.__numEquations=numEquations 222 self.__numSolutions=numSolutions 223 jgs 102 self.cleanCoefficients() 224 225 self.__operator=escript.Operator() 226 self.__operator_isValid=False 227 self.__righthandside=escript.Data() 228 self.__righthandside_isValid=False 229 self.__solution=escript.Data() 230 self.__solution_isValid=False 231 232 # set some default values: 233 self.__homogeneous_constraint=True 234 self.__row_function_space=escript.Solution(self.__domain) 235 self.__column_function_space=escript.Solution(self.__domain) 236 self.__tolerance=1.e-8 237 self.__solver_method=util.DEFAULT_METHOD 238 self.__matrix_type=self.__domain.getSystemMatrixTypeId(util.DEFAULT_METHOD,False) 239 self.__sym=False 240 self.__lumping=False 241 242 jgs 108 def createCoefficient(self, name): 243 """ 244 jgs 122 create a data object corresponding to coefficient name 245 @param name: 246 jgs 108 """ 247 return escript.Data(shape = getShapeOfCoefficient(name), \ 248 jgs 110 what = getFunctionSpaceForCoefficient(name)) 249 jgs 108 250 def __del__(self): 251 pass 252 253 jgs 102 def getCoefficient(self,name): 254 """ 255 jgs 122 return the value of the parameter name 256 jgs 102 257 jgs 122 @param name: 258 jgs 102 """ 259 jgs 108 return self.COEFFICIENTS[name].getValue() 260 jgs 102 261 jgs 108 def getCoefficientOfPDE(self,name): 262 """ 263 jgs 122 return the value of the coefficient name of the general PDE. 264 This method is called by the assembling routine it can be 265 overwritten to map coefficients of a particualr PDE to the general PDE. 266 267 @param name: 268 jgs 108 """ 269 return self.getCoefficient(name) 270 271 def hasCoefficient(self,name): 272 jgs 102 """ 273 jgs 122 return true if name is the name of a coefficient 274 jgs 102 275 jgs 122 @param name: 276 jgs 102 """ 277 jgs 108 return self.COEFFICIENTS.has_key(name) 278 jgs 102 279 jgs 108 def getFunctionSpaceForEquation(self): 280 """ 281 jgs 122 return true if the test functions should use reduced order 282 jgs 108 """ 283 return self.__row_function_space 284 285 def getFunctionSpaceForSolution(self): 286 """ 287 jgs 122 return true if the interpolation of the solution should use reduced order 288 jgs 108 """ 289 return self.__column_function_space 290 291 def setValue(self,**coefficients): 292 jgs 102 """ 293 jgs 122 sets new values to coefficients 294 jgs 102 295 jgs 122 @param coefficients: 296 jgs 102 """ 297 jgs 121 self.__setValue(**coefficients) 298 jgs 102 299 300 def cleanCoefficients(self): 301 """ 302 jgs 122 resets all coefficients to default values. 303 jgs 102 """ 304 jgs 108 for i in self.COEFFICIENTS.iterkeys(): 305 self.COEFFICIENTS[i].resetValue() 306 jgs 102 307 jgs 108 def createNewCoefficient(self,name): 308 """ 309 jgs 122 returns a new coefficient appropriate for coefficient name: 310 jgs 108 """ 311 return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name)) 312 313 314 jgs 102 def getShapeOfCoefficient(self,name): 315 """ 316 jgs 122 return the shape of the coefficient name 317 jgs 102 318 jgs 122 @param name: 319 jgs 102 """ 320 if self.hasCoefficient(name): 321 jgs 108 return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim()) 322 jgs 102 else: 323 raise ValueError,"Solution coefficient %s requested"%name 324 325 jgs 108 def getFunctionSpaceForCoefficient(self,name): 326 jgs 102 """ 327 jgs 122 return the atoms of the coefficient name 328 jgs 102 329 jgs 122 @param name: 330 jgs 102 """ 331 if self.hasCoefficient(name): 332 jgs 108 return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain()) 333 jgs 102 else: 334 raise ValueError,"Solution coefficient %s requested"%name 335 336 def alteredCoefficient(self,name): 337 """ 338 jgs 122 announce that coefficient name has been changed 339 jgs 102 340 jgs 122 @param name: 341 jgs 102 """ 342 if self.hasCoefficient(name): 343 jgs 108 if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator() 344 if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__rebuildRightHandSide() 345 jgs 102 else: 346 jgs 108 raise ValueError,"unknown coefficient %s requested"%name 347 jgs 102 348 # ===== debug ============================================================== 349 def setDebugOn(self): 350 """ 351 """ 352 self.__debug=not None 353 354 def setDebugOff(self): 355 """ 356 """ 357 self.__debug=None 358 359 def debug(self): 360 """ 361 jgs 122 returns true if the PDE is in the debug mode 362 jgs 102 """ 363 return self.__debug 364 365 #===== Lumping =========================== 366 def setLumpingOn(self): 367 """ 368 jgs 122 indicates to use matrix lumping 369 jgs 102 """ 370 if not self.isUsingLumping(): 371 if self.debug() : print "PDE Debug: lumping is set on" 372 self.__rebuildOperator() 373 self.__lumping=True 374 375 def setLumpingOff(self): 376 """ 377 jgs 122 switches off matrix lumping 378 jgs 102 """ 379 if self.isUsingLumping(): 380 if self.debug() : print "PDE Debug: lumping is set off" 381 self.__rebuildOperator() 382 self.__lumping=False 383 384 def setLumping(self,flag=False): 385 """ 386 jgs 122 set the matrix lumping flag to flag 387 jgs 102 """ 388 if flag: 389 self.setLumpingOn() 390 else: 391 self.setLumpingOff() 392 393 def isUsingLumping(self): 394 """ 395 jgs 122 396 jgs 102 """ 397 return self.__lumping 398 399 #============ method business ========================================================= 400 def setSolverMethod(self,solver=util.DEFAULT_METHOD): 401 """ 402 jgs 122 sets a new solver 403 jgs 102 """ 404 if not solver==self.getSolverMethod(): 405 self.__solver_method=solver 406 if self.debug() : print "PDE Debug: New solver is %s"%solver 407 self.__checkMatrixType() 408 409 def getSolverMethod(self): 410 """ 411 jgs 122 returns the solver method 412 jgs 102 """ 413 return self.__solver_method 414 415 #============ tolerance business ========================================================= 416 def setTolerance(self,tol=1.e-8): 417 """ 418 jgs 122 resets the tolerance to tol. 419 jgs 102 """ 420 if not tol>0: 421 raise ValueException,"Tolerance as to be positive" 422 if tol1: 573 self.__righthandside=escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True) 574 else: 575 self.__righthandside=escript.Data(0.,(),self.getFunctionSpaceForEquation(),True) 576 return self.__righthandside 577 jgs 102 578 jgs 121 def __getNewSolution(self): 579 jgs 102 """ 580 """ 581 jgs 121 if self.debug() : print "PDE Debug: New right hand side allocated" 582 if self.getNumSolutions()>1: 583 return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True) 584 else: 585 return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True) 586 jgs 102 587 jgs 121 def __makeFreshOperator(self): 588 jgs 102 """ 589 """ 590 if self.__operator.isEmpty(): 591 jgs 121 self.__operator=self.__getNewOperator() 592 jgs 102 if self.debug() : print "PDE Debug: New operator allocated" 593 else: 594 self.__operator.setValue(0.) 595 jgs 108 self.__operator.resetSolver() 596 jgs 102 if self.debug() : print "PDE Debug: Operator reset to zero" 597 return self.__operator 598 599 jgs 108 #============ some serivice functions ===================================================== 600 def getDomain(self): 601 """ 602 jgs 122 returns the domain of the PDE 603 jgs 108 """ 604 return self.__domain 605 606 def getDim(self): 607 """ 608 jgs 122 returns the spatial dimension of the PDE 609 jgs 108 """ 610 return self.getDomain().getDim() 611 612 def getNumEquations(self): 613 """ 614 jgs 122 returns the number of equations 615 jgs 108 """ 616 if self.__numEquations>0: 617 return self.__numEquations 618 else: 619 raise ValueError,"Number of equations is undefined. Please specify argument numEquations." 620 621 def getNumSolutions(self): 622 """ 623 jgs 122 returns the number of unknowns 624 jgs 108 """ 625 if self.__numSolutions>0: 626 return self.__numSolutions 627 else: 628 raise ValueError,"Number of solution is undefined. Please specify argument numSolutions." 629 630 631 def checkSymmetry(self,verbose=True): 632 """ 633 jgs 122 returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered. 634 jgs 108 """ 635 verbose=verbose or self.debug() 636 out=True 637 if self.getNumSolutions()!=self.getNumEquations(): 638 if verbose: print "non-symmetric PDE because of different number of equations and solutions" 639 out=False 640 else: 641 A=self.getCoefficientOfPDE("A") 642 if not A.isEmpty(): 643 tol=util.Lsup(A)*self.TOL 644 if self.getNumSolutions()>1: 645 for i in range(self.getNumEquations()): 646 for j in range(self.getDim()): 647 for k in range(self.getNumSolutions()): 648 for l in range(self.getDim()): 649 if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol: 650 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) 651 out=False 652 else: 653 for j in range(self.getDim()): 654 for l in range(self.getDim()): 655 if util.Lsup(A[j,l]-A[l,j])>tol: 656 if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j) 657 out=False 658 B=self.getCoefficientOfPDE("B") 659 C=self.getCoefficientOfPDE("C") 660 if B.isEmpty() and not C.isEmpty(): 661 if verbose: print "non-symmetric PDE because B is not present but C is" 662 out=False 663 elif not B.isEmpty() and C.isEmpty(): 664 if verbose: print "non-symmetric PDE because C is not present but B is" 665 out=False 666 elif not B.isEmpty() and not C.isEmpty(): 667 tol=(util.Lsup(B)+util.Lsup(C))*self.TOL/2. 668 if self.getNumSolutions()>1: 669 for i in range(self.getNumEquations()): 670 for j in range(self.getDim()): 671 for k in range(self.getNumSolutions()): 672 jgs 110 if util.Lsup(B[i,j,k]-C[k,i,j])>tol: 673 if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j) 674 jgs 108 out=False 675 else: 676 for j in range(self.getDim()): 677 if util.Lsup(B[j]-C[j])>tol: 678 if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j) 679 out=False 680 if self.getNumSolutions()>1: 681 D=self.getCoefficientOfPDE("D") 682 if not D.isEmpty(): 683 tol=util.Lsup(D)*self.TOL 684 for i in range(self.getNumEquations()): 685 for k in range(self.getNumSolutions()): 686 if util.Lsup(D[i,k]-D[k,i])>tol: 687 if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i) 688 out=False 689 690 return out 691 692 def getFlux(self,u): 693 """ 694 jgs 122 returns the flux J_ij for a given u 695 jgs 108 696 jgs 122 \f[ 697 J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij} 698 \f] 699 jgs 108 700 jgs 122 @param u: argument of the operator 701 jgs 108 """ 702 raise SystemError,"getFlux is not implemented yet" 703 return None 704 705 def applyOperator(self,u): 706 """ 707 jgs 122 applies the operator of the PDE to a given solution u in weak from 708 jgs 108 709 jgs 122 @param u: argument of the operator 710 jgs 108 """ 711 return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution()) 712 713 def getResidual(self,u): 714 """ 715 jgs 122 return the residual of u in the weak from 716 jgs 108 717 jgs 122 @param u: 718 jgs 108 """ 719 return self.applyOperator(u)-self.getRightHandSide() 720 721 jgs 121 def __setValue(self,**coefficients): 722 jgs 108 """ 723 jgs 122 sets new values to coefficient 724 jgs 108 725 jgs 122 @param coefficients: 726 jgs 108 """ 727 # check if the coefficients are legal: 728 for i in coefficients.iterkeys(): 729 if not self.hasCoefficient(i): 730 raise ValueError,"Attempt to set unknown coefficient %s"%i 731 # if the number of unknowns or equations is still unknown we try to estimate them: 732 if self.__numEquations<1 or self.__numSolutions<1: 733 for i,d in coefficients.iteritems(): 734 if hasattr(d,"shape"): 735 s=d.shape 736 elif hasattr(d,"getShape"): 737 s=d.getShape() 738 else: 739 s=numarray.array(d).shape 740 if s!=None: 741 # get number of equations and number of unknowns: 742 res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(s,self.getDim()) 743 if res==None: 744 raise ValueError,"Illegal shape %s of coefficient %s"%(s,i) 745 else: 746 if self.__numEquations<1: self.__numEquations=res[0] 747 if self.__numSolutions<1: self.__numSolutions=res[1] 748 if self.__numEquations<1: raise ValueError,"unidententified number of equations" 749 if self.__numSolutions<1: raise ValueError,"unidententified number of solutions" 750 # now we check the shape of the coefficient if numEquations and numSolutions are set: 751 for i,d in coefficients.iteritems(): 752 if d==None: 753 d2=escript.Data() 754 elif isinstance(d,escript.Data): 755 if d.isEmpty(): 756 d2=d 757 else: 758 d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i)) 759 else: 760 d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i)) 761 if not d2.isEmpty(): 762 if not self.getShapeOfCoefficient(i)==d2.getShape(): 763 raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),d2.getShape()) 764 # overwrite new values: 765 if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i 766 self.COEFFICIENTS[i].setValue(d2) 767 self.alteredCoefficient(i) 768 769 # reset the HomogeneousConstraintFlag: 770 self.__setHomogeneousConstraintFlag() 771 if len(coefficients)>0 and not self.isUsingLumping() and not self.__homogeneous_constraint: self.__rebuildSystem() 772 773 def __setHomogeneousConstraintFlag(self): 774 """ 775 jgs 122 checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly. 776 jgs 108 """ 777 self.__homogeneous_constraint=True 778 q=self.getCoefficientOfPDE("q") 779 r=self.getCoefficientOfPDE("r") 780 if not q.isEmpty() and not r.isEmpty(): 781 if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False 782 if self.debug(): 783 if self.__homogeneous_constraint: 784 print "PDE Debug: Constraints are homogeneous." 785 else: 786 print "PDE Debug: Constraints are inhomogeneous." 787 788 789 jgs 102 # ==== rebuild switches ===================================================================== 790 def __rebuildSolution(self,deep=False): 791 """ 792 jgs 122 indicates the PDE has to be reolved if the solution is requested 793 jgs 102 """ 794 if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved." 795 self.__solution_isValid=False 796 jgs 108 if deep: self.__solution=escript.Data() 797 jgs 102 798 799 def __rebuildOperator(self,deep=False): 800 """ 801 jgs 122 indicates the operator has to be rebuilt next time it is used 802 jgs 102 """ 803 if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt." 804 self.__rebuildSolution(deep) 805 self.__operator_isValid=False 806 if deep: self.__operator=escript.Operator() 807 808 def __rebuildRightHandSide(self,deep=False): 809 """ 810 jgs 122 indicates the right hand side has to be rebuild next time it is used 811 jgs 102 """ 812 if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt." 813 self.__rebuildSolution(deep) 814 self.__righthandside_isValid=False 815 if deep: self.__righthandside=escript.Data() 816 817 def __rebuildSystem(self,deep=False): 818 """ 819 jgs 122 annonced that all coefficient name has been changed 820 jgs 102 """ 821 self.__rebuildSolution(deep) 822 self.__rebuildOperator(deep) 823 self.__rebuildRightHandSide(deep) 824 825 def __checkMatrixType(self): 826 """ 827 jgs 122 reassess the matrix type and, if needed, initiates an operator rebuild 828 jgs 102 """ 829 new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric()) 830 if not new_matrix_type==self.__matrix_type: 831 if self.debug() : print "PDE Debug: Matrix type is now %d."%new_matrix_type 832 self.__matrix_type=new_matrix_type 833 self.__rebuildOperator(deep=True) 834 835 #============ assembling ======================================================= 836 jgs 108 def __copyConstraint(self): 837 jgs 102 """ 838 jgs 122 copies the constrint condition into u 839 jgs 102 """ 840 jgs 108 if not self.__righthandside.isEmpty(): 841 q=self.getCoefficientOfPDE("q") 842 r=self.getCoefficientOfPDE("r") 843 if not q.isEmpty(): 844 if r.isEmpty(): 845 r2=escript.Data(0,self.__righthandside.getShape(),self.__righthandside.getFunctionSpace()) 846 else: 847 r2=escript.Data(r,self.__righthandside.getFunctionSpace()) 848 self.__righthandside.copyWithMask(r2,escript.Data(q,self.__righthandside.getFunctionSpace())) 849 jgs 102 850 jgs 108 def __applyConstraint(self): 851 jgs 102 """ 852 jgs 122 applies the constraints defined by q and r to the system 853 jgs 102 """ 854 jgs 108 q=self.getCoefficientOfPDE("q") 855 r=self.getCoefficientOfPDE("r") 856 jgs 102 if not q.isEmpty() and not self.__operator.isEmpty(): 857 # q is the row and column mask to indicate where constraints are set: 858 row_q=escript.Data(q,self.getFunctionSpaceForEquation()) 859 col_q=escript.Data(q,self.getFunctionSpaceForSolution()) 860 jgs 121 u=self.__getNewSolution() 861 jgs 102 if r.isEmpty(): 862 jgs 121 r_s=self.__getNewSolution() 863 jgs 102 else: 864 r_s=escript.Data(r,self.getFunctionSpaceForSolution()) 865 u.copyWithMask(r_s,col_q) 866 jgs 108 if self.isUsingLumping(): 867 self.__operator.copyWithMask(escript.Data(1,q.getShape(),self.getFunctionSpaceForEquation()),row_q) 868 else: 869 if not self.__righthandside.isEmpty(): self.__righthandside-=self.__operator*u 870 self.__operator.nullifyRowsAndCols(row_q,col_q,1.) 871 jgs 102 872 jgs 108 def getSystem(self): 873 """ 874 jgs 122 return the operator and right hand side of the PDE 875 jgs 108 """ 876 if not self.__operator_isValid or not self.__righthandside_isValid: 877 if self.isUsingLumping(): 878 if not self.__operator_isValid: 879 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution(): 880 raise TypeError,"Lumped matrix requires same order for equations and unknowns" 881 if not self.getCoefficientOfPDE("A").isEmpty(): 882 raise Warning,"Lumped matrix does not allow coefficient A" 883 if not self.getCoefficientOfPDE("B").isEmpty(): 884 raise Warning,"Lumped matrix does not allow coefficient B" 885 if not self.getCoefficientOfPDE("C").isEmpty(): 886 raise Warning,"Lumped matrix does not allow coefficient C" 887 if self.debug() : print "PDE Debug: New lumped operator is built." 888 jgs 121 mat=self.__getNewOperator() 889 jgs 108 self.getDomain().addPDEToSystem(mat,escript.Data(), \ 890 self.getCoefficientOfPDE("A"), \ 891 self.getCoefficientOfPDE("B"), \ 892 self.getCoefficientOfPDE("C"), \ 893 self.getCoefficientOfPDE("D"), \ 894 escript.Data(), \ 895 escript.Data(), \ 896 self.getCoefficientOfPDE("d"), \ 897 escript.Data(),\ 898 self.getCoefficientOfPDE("d_contact"), \ 899 escript.Data()) 900 self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True) 901 self.__applyConstraint() 902 self.__operator_isValid=True 903 if not self.__righthandside_isValid: 904 if self.debug() : print "PDE Debug: New right hand side is built." 905 jgs 121 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \ 906 jgs 108 self.getCoefficientOfPDE("X"), \ 907 self.getCoefficientOfPDE("Y"),\ 908 self.getCoefficientOfPDE("y"),\ 909 self.getCoefficientOfPDE("y_contact")) 910 self.__copyConstraint() 911 self.__righthandside_isValid=True 912 else: 913 if not self.__operator_isValid and not self.__righthandside_isValid: 914 if self.debug() : print "PDE Debug: New system is built." 915 jgs 121 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \ 916 jgs 108 self.getCoefficientOfPDE("A"), \ 917 self.getCoefficientOfPDE("B"), \ 918 self.getCoefficientOfPDE("C"), \ 919 self.getCoefficientOfPDE("D"), \ 920 self.getCoefficientOfPDE("X"), \ 921 self.getCoefficientOfPDE("Y"), \ 922 self.getCoefficientOfPDE("d"), \ 923 self.getCoefficientOfPDE("y"), \ 924 self.getCoefficientOfPDE("d_contact"), \ 925 self.getCoefficientOfPDE("y_contact")) 926 self.__applyConstraint() 927 self.__copyConstraint() 928 self.__operator_isValid=True 929 self.__righthandside_isValid=True 930 elif not self.__righthandside_isValid: 931 if self.debug() : print "PDE Debug: New right hand side is built." 932 jgs 121 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \ 933 jgs 108 self.getCoefficientOfPDE("X"), \ 934 self.getCoefficientOfPDE("Y"),\ 935 self.getCoefficientOfPDE("y"),\ 936 self.getCoefficientOfPDE("y_contact")) 937 self.__copyConstraint() 938 self.__righthandside_isValid=True 939 elif not self.__operator_isValid: 940 if self.debug() : print "PDE Debug: New operator is built." 941 jgs 121 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \ 942 jgs 108 self.getCoefficientOfPDE("A"), \ 943 self.getCoefficientOfPDE("B"), \ 944 self.getCoefficientOfPDE("C"), \ 945 self.getCoefficientOfPDE("D"), \ 946 escript.Data(), \ 947 escript.Data(), \ 948 self.getCoefficientOfPDE("d"), \ 949 escript.Data(),\ 950 self.getCoefficientOfPDE("d_contact"), \ 951 escript.Data()) 952 self.__applyConstraint() 953 self.__operator_isValid=True 954 return (self.__operator,self.__righthandside) 955 jgs 102 def getOperator(self): 956 """ 957 jgs 122 returns the operator of the PDE 958 jgs 102 """ 959 jgs 108 return self.getSystem()[0] 960 jgs 102 961 jgs 108 def getRightHandSide(self): 962 jgs 102 """ 963 jgs 122 returns the right hand side of the PDE 964 jgs 102 """ 965 jgs 108 return self.getSystem()[1] 966 jgs 102 967 def solve(self,**options): 968 """ 969 jgs 122 solve the PDE 970 jgs 102 971 jgs 122 @param options: 972 jgs 102 """ 973 mat,f=self.getSystem() 974 if self.isUsingLumping(): 975 out=f/mat 976 else: 977 options[util.TOLERANCE_KEY]=self.getTolerance() 978 options[util.METHOD_KEY]=self.getSolverMethod() 979 options[util.SYMMETRY_KEY]=self.isSymmetric() 980 if self.debug() : print "PDE Debug: solver options: ",options 981 out=mat.solve(f,options) 982 return out 983 984 def getSolution(self,**options): 985 """ 986 jgs 122 returns the solution of the PDE 987 jgs 102 988 jgs 122 @param options: 989 jgs 102 """ 990 if not self.__solution_isValid: 991 if self.debug() : print "PDE Debug: PDE is resolved." 992 self.__solution=self.solve(**options) 993 self.__solution_isValid=True 994 return self.__solution 995 996 jgs 110 997 998 jgs 122 def ELMAN_RAMAGE(P): 999 """ """ 1000 return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15)) 1001 jgs 110 def SIMPLIFIED_BROOK_HUGHES(P): 1002 jgs 122 """ """ 1003 c=(P-3.).whereNegative() 1004 return P/6.*c+1./2.*(1.-c) 1005 jgs 110 1006 jgs 122 def HALF(P): 1007 """ """ 1008 return escript.Scalar(0.5,P.getFunctionSpace()) 1009 jgs 110 1010 jgs 108 class AdvectivePDE(LinearPDE): 1011 """ 1012 jgs 122 Class to handle a linear PDE dominated by advective terms: 1013 jgs 108 1014 class to define a linear PDE of the form 1015 jgs 104 1016 jgs 122 \f[ 1017 -(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 1018 \f] 1019 jgs 102 1020 jgs 122 with boundary conditons: 1021 jgs 102 1022 jgs 122 \f[ 1023 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i 1024 \f] 1025 jgs 102 1026 jgs 122 and contact conditions 1027 jgs 102 1028 jgs 122 \f[ 1029 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d^{contact}_{ik}[u_k] = - n_j*X_{ij} + y^{contact}_{i} 1030 \f] 1031 jgs 102 1032 jgs 122 and constraints: 1033 jgs 102 1034 jgs 122 \f[ 1035 u_i=r_i \quad \mathrm{where} \quad q_i>0 1036 \f] 1037 jgs 110 """ 1038 def __init__(self,domain,numEquations=0,numSolutions=0,xi=ELMAN_RAMAGE): 1039 LinearPDE.__init__(self,domain,numEquations,numSolutions) 1040 self.__xi=xi 1041 self.__Xi=escript.Data() 1042 jgs 102 1043 jgs 110 def __calculateXi(self,peclet_factor,Z,h): 1044 Z_max=util.Lsup(Z) 1045 if Z_max>0.: 1046 return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.TOL) 1047 else: 1048 return 0. 1049 jgs 102 1050 jgs 110 def setValue(self,**args): 1051 if "A" in args.keys() or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data() 1052 jgs 121 self._LinearPDE__setValue(**args) 1053 jgs 110 1054 def getXi(self): 1055 if self.__Xi.isEmpty(): 1056 B=self.getCoefficient("B") 1057 C=self.getCoefficient("C") 1058 A=self.getCoefficient("A") 1059 h=self.getDomain().getSize() 1060 self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A")) 1061 if not C.isEmpty() or not B.isEmpty(): 1062 if not C.isEmpty() and not B.isEmpty(): 1063 Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A")) 1064 if self.getNumEquations()>1: 1065 if self.getNumSolutions()>1: 1066 for i in range(self.getNumEquations()): 1067 for k in range(self.getNumSolutions()): 1068 for l in range(self.getDim()): Z2+=(C[i,k,l]-B[i,l,k])**2 1069 else: 1070 for i in range(self.getNumEquations()): 1071 for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2 1072 else: 1073 if self.getNumSolutions()>1: 1074 for k in range(self.getNumSolutions()): 1075 for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2 1076 else: 1077 for l in range(self.getDim()): Z2+=(C[l]-B[l])**2 1078 length_of_Z=util.sqrt(Z2) 1079 elif C.isEmpty(): 1080 length_of_Z=util.length(B) 1081 else: 1082 length_of_Z=util.length(C) 1083 jgs 102 1084 jgs 110 Z_max=util.Lsup(length_of_Z) 1085 if Z_max>0.: 1086 length_of_A=util.length(A) 1087 A_max=util.Lsup(length_of_A) 1088 if A_max>0: 1089 inv_A=1./(length_of_A+A_max*self.TOL) 1090 else: 1091 inv_A=1./self.TOL 1092 peclet_number=length_of_Z*h/2*inv_A 1093 xi=self.__xi(peclet_number) 1094 self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL) 1095 print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z) 1096 return self.__Xi 1097 1098 jgs 102 1099 jgs 108 def getCoefficientOfPDE(self,name): 1100 """ 1101 jgs 122 return the value of the coefficient name of the general PDE 1102 1103 @param name: 1104 jgs 108 """ 1105 jgs 110 if not self.getNumEquations() == self.getNumSolutions(): 1106 raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal." 1107 1108 jgs 108 if name == "A" : 1109 A=self.getCoefficient("A") 1110 B=self.getCoefficient("B") 1111 C=self.getCoefficient("C") 1112 jgs 110 if B.isEmpty() and C.isEmpty(): 1113 Aout=A 1114 else: 1115 if A.isEmpty(): 1116 Aout=self.createNewCoefficient("A") 1117 else: 1118 Aout=A[:] 1119 Xi=self.getXi() 1120 if self.getNumEquations()>1: 1121 for i in range(self.getNumEquations()): 1122 for j in range(self.getDim()): 1123 jgs 108 for k in range(self.getNumSolutions()): 1124 jgs 110 for l in range(self.getDim()): 1125 if not C.isEmpty() and not B.isEmpty(): 1126 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]) 1127 elif C.isEmpty(): 1128 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k] 1129 else: 1130 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l] 1131 else: 1132 for j in range(self.getDim()): 1133 for l in range(self.getDim()): 1134 if not C.isEmpty() and not B.isEmpty(): 1135 Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l]) 1136 elif C.isEmpty(): 1137 Aout[j,l]+=Xi*B[j]*B[l] 1138 else: 1139 Aout[j,l]+=Xi*C[j]*C[l] 1140 return Aout 1141 jgs 108 elif name == "B" : 1142 jgs 110 B=self.getCoefficient("B") 1143 C=self.getCoefficient("C") 1144 D=self.getCoefficient("D") 1145 if C.isEmpty() or D.isEmpty(): 1146 Bout=B 1147 else: 1148 Xi=self.getXi() 1149 if B.isEmpty(): 1150 Bout=self.createNewCoefficient("B") 1151 else: 1152 Bout=B[:] 1153 if self.getNumEquations()>1: 1154 for k in range(self.getNumSolutions()): 1155 for p in range(self.getNumEquations()): 1156 tmp=Xi*D[p,k] 1157 for i in range(self.getNumEquations()): 1158 for j in range(self.getDim()): 1159 Bout[i,j,k]+=tmp*C[p,i,j] 1160 else: 1161 tmp=Xi*D 1162 for j in range(self.getDim()): Bout[j]+=tmp*C[j] 1163 return Bout 1164 jgs 108 elif name == "C" : 1165 jgs 110 B=self.getCoefficient("B") 1166 C=self.getCoefficient("C") 1167 D=self.getCoefficient("D") 1168 if B.isEmpty() or D.isEmpty(): 1169 Cout=C 1170 else: 1171 Xi=self.getXi() 1172 if C.isEmpty(): 1173 Cout=self.createNewCoefficient("C") 1174 else: 1175 Cout=C[:] 1176 if self.getNumEquations()>1: 1177 for k in range(self.getNumSolutions()): 1178 for p in range(self.getNumEquations()): 1179 tmp=Xi*D[p,k] 1180 for i in range(self.getNumEquations()): 1181 for l in range(self.getDim()): 1182 Cout[i,k,l]+=tmp*B[p,l,i] 1183 else: 1184 tmp=Xi*D 1185 for j in range(self.getDim()): Cout[j]+=tmp*B[j] 1186 return Cout 1187 jgs 108 elif name == "D" : 1188 return self.getCoefficient("D") 1189 elif name == "X" : 1190 jgs 110 X=self.getCoefficient("X") 1191 Y=self.getCoefficient("Y") 1192 B=self.getCoefficient("B") 1193 C=self.getCoefficient("C") 1194 if Y.isEmpty() or (B.isEmpty() and C.isEmpty()): 1195 Xout=X 1196 else: 1197 if X.isEmpty(): 1198 Xout=self.createNewCoefficient("X") 1199 else: 1200 Xout=X[:] 1201 Xi=self.getXi() 1202 if self.getNumEquations()>1: 1203 for p in range(self.getNumEquations()): 1204 tmp=Xi*Y[p] 1205 for i in range(self.getNumEquations()): 1206 for j in range(self.getDim()): 1207 if not C.isEmpty() and not B.isEmpty(): 1208 Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i]) 1209 elif C.isEmpty(): 1210 Xout[i,j]-=tmp*B[p,j,i] 1211 else: 1212 Xout[i,j]+=tmp*C[p,i,j] 1213 else: 1214 tmp=Xi*Y 1215 for j in range(self.getDim()): 1216 if not C.isEmpty() and not B.isEmpty(): 1217 Xout[j]+=tmp*(C[j]-B[j]) 1218 elif C.isEmpty(): 1219 Xout[j]-=tmp*B[j] 1220 else: 1221 Xout[j]+=tmp*C[j] 1222 return Xout 1223 jgs 108 elif name == "Y" : 1224 return self.getCoefficient("Y") 1225 elif name == "d" : 1226 return self.getCoefficient("d") 1227 elif name == "y" : 1228 return self.getCoefficient("y") 1229 elif name == "d_contact" : 1230 return self.getCoefficient("d_contact") 1231 elif name == "y_contact" : 1232 return self.getCoefficient("y_contact") 1233 elif name == "r" : 1234 return self.getCoefficient("r") 1235 elif name == "q" : 1236 return self.getCoefficient("q") 1237 else: 1238 raise SystemError,"unknown PDE coefficient %s",name 1239 1240 1241 jgs 102 class Poisson(LinearPDE): 1242 """ 1243 jgs 142 Class to define a Poisson equation problem: 1244 jgs 122 1245 jgs 102 class to define a linear PDE of the form 1246 jgs 122 \f[ 1247 -u_{,jj} = f 1248 \f] 1249 1250 with boundary conditons: 1251 1252 \f[ 1253 n_j*u_{,j} = 0 1254 \f] 1255 1256 and constraints: 1257 1258 \f[ 1259 u=0 \quad \mathrm{where} \quad q>0 1260 \f] 1261 jgs 102 """ 1262 1263 jgs 108 def __init__(self,domain,f=escript.Data(),q=escript.Data()): 1264 LinearPDE.__init__(self,domain,1,1) 1265 self.COEFFICIENTS={ 1266 "f" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE), 1267 "q" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)} 1268 jgs 102 self.setSymmetryOn() 1269 self.setValue(f,q) 1270 1271 def setValue(self,f=escript.Data(),q=escript.Data()): 1272 jgs 142 """set value of PDE parameters f and q""" 1273 jgs 121 self._LinearPDE__setValue(f=f,q=q) 1274 jgs 108 1275 def getCoefficientOfPDE(self,name): 1276 """ 1277 jgs 122 return the value of the coefficient name of the general PDE 1278 1279 @param name: 1280 jgs 108 """ 1281 if name == "A" : 1282 return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain())) 1283 elif name == "B" : 1284 return escript.Data() 1285 elif name == "C" : 1286 return escript.Data() 1287 elif name == "D" : 1288 return escript.Data() 1289 elif name == "X" : 1290 return escript.Data() 1291 elif name == "Y" : 1292 return self.getCoefficient("f") 1293 elif name == "d" : 1294 return escript.Data() 1295 elif name == "y" : 1296 return escript.Data() 1297 elif name == "d_contact" : 1298 return escript.Data() 1299 elif name == "y_contact" : 1300 return escript.Data() 1301 elif name == "r" : 1302 return escript.Data() 1303 elif name == "q" : 1304 return self.getCoefficient("q") 1305 else: 1306 raise SystemError,"unknown PDE coefficient %s",name 1307 jgs 121 1308 jgs 142 class LameEquation(LinearPDE): 1309 """ 1310 Class to define a Lame equation problem: 1311 1312 class to define a linear PDE of the form 1313 \f[ 1314 -(\lambda (u_{i,j}+u_{j,i}))_{,j} - \mu u_{j,ji}} = F_i -\sigma_{ij,j} 1315 \f] 1316 1317 with boundary conditons: 1318 1319 \f[ 1320 n_j(\lambda(u_{i,j}+u_{j,i})-sigma_{ij}) + n_i\mu u_{j,j} = f_i 1321 \f] 1322 1323 and constraints: 1324 1325 \f[ 1326 u_i=r_i \quad \mathrm{where} \quad q_i>0 1327 \f] 1328 """ 1329 1330 def __init__(self,domain,f=escript.Data(),q=escript.Data()): 1331 LinearPDE.__init__(self,domain,domain.getDim(),domain.getDim()) 1332 self.COEFFICIENTS={ 1333 "lame_lambda" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR), 1334 "lame_mu" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR), 1335 "F" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE), 1336 "sigma" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.EQUATION),PDECoefficient.RIGHTHANDSIDE), 1337 "f" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE), 1338 "r" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH), 1339 "q" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)} 1340 self.setSymmetryOn() 1341 1342 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()): 1343 """set value of PDE parameters""" 1344 self._LinearPDE__setValue(lame_lambda=lame_lambda, \ 1345 lame_mu=lame_mu, \ 1346 F=F, \ 1347 sigma=sigma, \ 1348 f=f, \ 1349 r=r, \ 1350 q=q) 1351 def getCoefficientOfPDE(self,name): 1352 """ 1353 return the value of the coefficient name of the general PDE 1354 1355 @param name: 1356 """ 1357 if name == "A" : 1358 A =self.createNewCoefficient("A") 1359 for i in range(self.getDim()): 1360 for j in range(self.getDim()): 1361 out[i,i,j,j] += self.getCoefficient("lame_mu") 1362 out[i,j,j,i] += self.getCoefficient("lame_lambda") 1363 out[i,j,i,j] += self.getCoefficient("lame_lambda") 1364 return out 1365 elif name == "B" : 1366 return escript.Data() 1367 elif name == "C" : 1368 return escript.Data() 1369 elif name == "D" : 1370 return escript.Data() 1371 elif name == "X" : 1372 return self.getCoefficient("sigma") 1373 elif name == "Y" : 1374 return self.getCoefficient("F") 1375 elif name == "d" : 1376 return escript.Data() 1377 elif name == "y" : 1378 return self.getCoefficient("f") 1379 elif name == "d_contact" : 1380 return escript.Data() 1381 elif name == "y_contact" : 1382 return escript.Data() 1383 elif name == "r" : 1384 return self.getCoefficient("r") 1385 elif name == "q" : 1386 return self.getCoefficient("q") 1387 else: 1388 raise SystemError,"unknown PDE coefficient %s",name 1389 1390 jgs 121 # $Log$ 1391 jgs 142 # Revision 1.9 2005/07/25 05:28:13 jgs 1392 # Merge of development branch back to main trunk on 2005-07-25 1393 # 1394 jgs 122 # Revision 1.8 2005/06/09 05:37:59 jgs 1395 # Merge of development branch back to main trunk on 2005-06-09 1396 # 1397 jgs 121 # Revision 1.7 2005/05/06 04:26:10 jgs 1398 # Merge of development branch back to main trunk on 2005-05-06 1399 # 1400 jgs 142 # Revision 1.1.2.24 2005/07/22 06:37:11 gross 1401 # some extensions to modellib and linearPDEs 1402 # 1403 jgs 122 # Revision 1.1.2.23 2005/05/13 00:55:20 cochrane 1404 # Fixed up some docstrings. Moved module-level functions to top of file so 1405 # that epydoc and doxygen can pick them up properly. 1406 # 1407 # Revision 1.1.2.22 2005/05/12 11:41:30 gross 1408 # some basic Models have been added 1409 # 1410 # Revision 1.1.2.21 2005/05/12 07:16:12 cochrane 1411 # Moved ELMAN_RAMAGE, SIMPLIFIED_BROOK_HUGHES, and HALF functions to bottom of 1412 # file so that the AdvectivePDE class is picked up by doxygen. Some 1413 # reformatting of docstrings. Addition of code to make equations come out 1414 # as proper LaTeX. 1415 # 1416 jgs 121 # Revision 1.1.2.20 2005/04/15 07:09:08 gross 1417 # some problems with functionspace and linearPDEs fixed. 1418 # 1419 # Revision 1.1.2.19 2005/03/04 05:27:07 gross 1420 # bug in SystemPattern fixed. 1421 # 1422 # Revision 1.1.2.18 2005/02/08 06:16:45 gross 1423 # Bugs in AdvectivePDE fixed, AdvectiveTest is stable but more testing is needed 1424 # 1425 # Revision 1.1.2.17 2005/02/08 05:56:19 gross 1426 # Reference Number handling added 1427 # 1428 # Revision 1.1.2.16 2005/02/07 04:41:28 gross 1429 # some function exposed to python to make mesh merging running 1430 # 1431 # Revision 1.1.2.15 2005/02/03 00:14:44 gross 1432 # timeseries add and ESySParameter.py renames esysXML.py for consistence 1433 # 1434 # Revision 1.1.2.14 2005/02/01 06:44:10 gross 1435 # new implementation of AdvectivePDE which now also updates right hand side. systems of PDEs are still not working 1436 # 1437 # Revision 1.1.2.13 2005/01/25 00:47:07 gross 1438 # updates in the documentation 1439 # 1440 # Revision 1.1.2.12 2005/01/12 01:28:04 matt 1441 # Added createCoefficient method for linearPDEs. 1442 # 1443 # Revision 1.1.2.11 2005/01/11 01:55:34 gross 1444 # a problem in linearPDE class fixed 1445 # 1446 # Revision 1.1.2.10 2005/01/07 01:13:29 gross 1447 # some bugs in linearPDE fixed 1448 # 1449 # Revision 1.1.2.9 2005/01/06 06:24:58 gross 1450 # some bugs in slicing fixed 1451 # 1452 # Revision 1.1.2.8 2005/01/05 04:21:40 gross 1453 # FunctionSpace checking/matchig in slicing added 1454 # 1455 # Revision 1.1.2.7 2004/12/29 10:03:41 gross 1456 # bug in setValue fixed 1457 # 1458 # Revision 1.1.2.6 2004/12/29 05:29:59 gross 1459 # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data() 1460 # 1461 # Revision 1.1.2.5 2004/12/29 00:18:41 gross 1462 # AdvectivePDE added 1463 # 1464 # Revision 1.1.2.4 2004/12/24 06:05:41 gross 1465 # some changes in linearPDEs to add AdevectivePDE 1466 # 1467 # Revision 1.1.2.3 2004/12/16 00:12:34 gross 1468 # __init__ of LinearPDE does not accept any coefficient anymore 1469 # 1470 # Revision 1.1.2.2 2004/12/14 03:55:01 jgs 1471 # *** empty log message *** 1472 # 1473 # Revision 1.1.2.1 2004/12/12 22:53:47 gross 1474 # linearPDE has been renamed LinearPDE 1475 # 1476 # Revision 1.1.1.1.2.7 2004/12/07 10:13:08 gross 1477 # GMRES added 1478 # 1479 # Revision 1.1.1.1.2.6 2004/12/07 03:19:50 gross 1480 # options for GMRES and PRES20 added 1481 # 1482 # Revision 1.1.1.1.2.5 2004/12/01 06:25:15 gross 1483 # some small changes 1484 # 1485 # Revision 1.1.1.1.2.4 2004/11/24 01:50:21 gross 1486 # Finley solves 4M unknowns now 1487 # 1488 # Revision 1.1.1.1.2.3 2004/11/15 06:05:26 gross 1489 # poisson solver added 1490 # 1491 # Revision 1.1.1.1.2.2 2004/11/12 06:58:15 gross 1492 # 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 1493 # 1494 # Revision 1.1.1.1.2.1 2004/10/28 22:59:22 gross 1495 # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed 1496 # 1497 # Revision 1.1.1.1 2004/10/26 06:53:56 jgs 1498 # initial import of project esys2 1499 # 1500 # Revision 1.3.2.3 2004/10/26 06:43:48 jgs 1501 # committing Lutz's and Paul's changes to brach jgs 1502 # 1503 # Revision 1.3.4.1 2004/10/20 05:32:51 cochrane 1504 # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form. 1505 # 1506 # Revision 1.3 2004/09/23 00:53:23 jgs 1507 # minor fixes 1508 # 1509 # Revision 1.1 2004/08/28 12:58:06 gross 1510 # SimpleSolve is not running yet: problem with == of functionsspace 1511 # 1512 #

## Properties

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