Contents of /trunk/esys2/escript/py_src/linearPDEs.py

Revision 147 - (show annotations)
Fri Aug 12 01:45:47 2005 UTC (15 years, 2 months ago) by jgs
File MIME type: text/x-python
File size: 55333 byte(s)
erge of development branch dev-02 back to main trunk on 2005-08-12


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

Properties

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