/[escript]/trunk/esys2/escript/py_src/linearPDEs.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 121 by jgs, Fri May 6 04:26:16 2005 UTC revision 122 by jgs, Thu Jun 9 05:38:05 2005 UTC
# Line 3  Line 3 
3  ## @file linearPDEs.py  ## @file linearPDEs.py
4    
5  """  """
6  @brief Functions and classes for linear PDEs  Functions and classes for linear PDEs
7  """  """
8    
9  import escript  import escript
# Line 13  import numarray Line 13  import numarray
13    
14  def _CompTuple2(t1,t2):  def _CompTuple2(t1,t2):
15     """     """
16     @brief     Compare two tuples
17    
18     @param t1     \param t1 The first tuple
19     @param t2     \param t2 The second tuple
20     """     """
21    
22     dif=t1[0]+t1[1]-(t2[0]+t2[1])     dif=t1[0]+t1[1]-(t2[0]+t2[1])
23     if dif<0: return 1     if dif<0: return 1
24     elif dif>0: return -1     elif dif>0: return -1
25     else: return 0     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:  class PDECoefficient:
38      """      """
39      @brief      A class for PDE coefficients
40      """      """
41      # identifier for location of Data objects defining COEFFICIENTS      # identifier for location of Data objects defining COEFFICIENTS
42      INTERIOR=0      INTERIOR=0
# Line 44  class PDECoefficient: Line 55  class PDECoefficient:
55      BOTH=7      BOTH=7
56      def __init__(self,where,pattern,altering):      def __init__(self,where,pattern,altering):
57         """         """
58         @brief Initialise a PDE Coefficient type         Initialise a PDE Coefficient type
59         """         """
60         self.what=where         self.what=where
61         self.pattern=pattern         self.pattern=pattern
# Line 53  class PDECoefficient: Line 64  class PDECoefficient:
64    
65      def resetValue(self):      def resetValue(self):
66         """         """
67         @brief resets coefficient value to default         resets coefficient value to default
68         """         """
69         self.value=escript.Data()         self.value=escript.Data()
70    
71      def getFunctionSpace(self,domain):      def getFunctionSpace(self,domain):
72         """         """
73         @brief defines the FunctionSpace of the coefficient on the domain         defines the FunctionSpace of the coefficient on the domain
74    
75         @param domain         @param domain:
76         """         """
77         if self.what==self.INTERIOR: return escript.Function(domain)         if self.what==self.INTERIOR: return escript.Function(domain)
78         elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)         elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)
# Line 70  class PDECoefficient: Line 81  class PDECoefficient:
81    
82      def getValue(self):      def getValue(self):
83         """         """
84         @brief returns the value of the coefficient:         returns the value of the coefficient:
85         """         """
86         return self.value         return self.value
87            
88      def setValue(self,newValue):      def setValue(self,newValue):
89         """         """
90         @brief set the value of the coefficient to new value         set the value of the coefficient to new value
91         """         """
92         self.value=newValue         self.value=newValue
93            
94      def isAlteringOperator(self):      def isAlteringOperator(self):
95          """          """
96      @brief return true if the operator of the PDE is changed when the coefficient is changed      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:          if self.altering==self.OPERATOR or self.altering==self.BOTH:
99              return not None              return not None
# Line 91  class PDECoefficient: Line 102  class PDECoefficient:
102    
103      def isAlteringRightHandSide(self):      def isAlteringRightHandSide(self):
104          """          """
105      @brief return true if the right hand side of the PDE is changed when the coefficient is changed      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:          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
108              return not None              return not None
# Line 100  class PDECoefficient: Line 111  class PDECoefficient:
111    
112      def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3):      def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3):
113         """         """
114         @brief tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim         tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim
115    
116         @param shape         @param shape:
117         @param dim         @param dim:
118         """         """
119         if len(shape)>0:         if len(shape)>0:
120             num=max(shape)+1             num=max(shape)+1
# Line 124  class PDECoefficient: Line 135  class PDECoefficient:
135    
136      def buildShape(self,e=1,u=1,dim=3):      def buildShape(self,e=1,u=1,dim=3):
137          """          """
138      @brief builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim      builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim
139    
140      @param e      @param e:
141      @param u      @param u:
142      @param dim      @param dim:
143      """      """
144          s=()          s=()
145          for i in self.pattern:          for i in self.pattern:
# Line 142  class PDECoefficient: Line 153  class PDECoefficient:
153    
154  class LinearPDE:  class LinearPDE:
155     """     """
156     @brief Class to handel a linear PDE     Class to handle a linear PDE
157        
158     class to define a linear PDE of the form     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       -(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:     with boundary conditons:
165    
166          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i     \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     and contact conditions
171    
172          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i     \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:     and constraints:
177    
178           u_i=r_i where q_i>0     \f[
179       u_i=r_i \quad \mathrm{where} \quad q_i>0
180       \f]
181    
182     """     """
183     TOL=1.e-13     TOL=1.e-13
# Line 177  class LinearPDE: Line 196  class LinearPDE:
196    
197     def __init__(self,domain,numEquations=0,numSolutions=0):     def __init__(self,domain,numEquations=0,numSolutions=0):
198       """       """
199       @brief initializes a new linear PDE.       initializes a new linear PDE.
200    
201       @param args       @param args:
202       """       """
203       # COEFFICIENTS can be overwritten by subclasses:       # COEFFICIENTS can be overwritten by subclasses:
204       self.COEFFICIENTS={       self.COEFFICIENTS={
# Line 222  class LinearPDE: Line 241  class LinearPDE:
241    
242     def createCoefficient(self, name):     def createCoefficient(self, name):
243       """       """
244       @brief create a data object corresponding to coefficient name       create a data object corresponding to coefficient name
245       @param name       @param name:
246       """       """
247       return escript.Data(shape = getShapeOfCoefficient(name), \       return escript.Data(shape = getShapeOfCoefficient(name), \
248                           what = getFunctionSpaceForCoefficient(name))                           what = getFunctionSpaceForCoefficient(name))
# Line 233  class LinearPDE: Line 252  class LinearPDE:
252    
253     def getCoefficient(self,name):     def getCoefficient(self,name):
254       """       """
255       @brief return the value of the parameter name       return the value of the parameter name
256    
257       @param name       @param name:
258       """       """
259       return self.COEFFICIENTS[name].getValue()       return self.COEFFICIENTS[name].getValue()
260    
261     def getCoefficientOfPDE(self,name):     def getCoefficientOfPDE(self,name):
262       """       """
263       @brief return the value of the coefficient name of the general PDE. This method is called by the assembling routine       return the value of the coefficient name of the general PDE.
264              it can be overwritten to map coefficients of a particualr PDE to the general PDE.       This method is called by the assembling routine it can be
265       @param name       overwritten to map coefficients of a particualr PDE to the general PDE.
266    
267         @param name:
268       """       """
269       return self.getCoefficient(name)       return self.getCoefficient(name)
270    
271     def hasCoefficient(self,name):     def hasCoefficient(self,name):
272        """        """
273        @brief return true if name is the name of a coefficient        return true if name is the name of a coefficient
274    
275        @param name        @param name:
276        """        """
277        return self.COEFFICIENTS.has_key(name)        return self.COEFFICIENTS.has_key(name)
278    
279     def getFunctionSpaceForEquation(self):     def getFunctionSpaceForEquation(self):
280       """       """
281       @brief return true if the test functions should use reduced order       return true if the test functions should use reduced order
282       """       """
283       return self.__row_function_space       return self.__row_function_space
284    
285     def getFunctionSpaceForSolution(self):     def getFunctionSpaceForSolution(self):
286       """       """
287       @brief return true if the interpolation of the solution should use reduced order       return true if the interpolation of the solution should use reduced order
288       """       """
289       return self.__column_function_space       return self.__column_function_space
290    
291     def setValue(self,**coefficients):     def setValue(self,**coefficients):
292        """        """
293        @brief sets new values to coefficients        sets new values to coefficients
294    
295        @param coefficients        @param coefficients:
296        """        """
297        self.__setValue(**coefficients)        self.__setValue(**coefficients)
298                
299    
300     def cleanCoefficients(self):     def cleanCoefficients(self):
301       """       """
302       @brief resets all coefficients to default values.       resets all coefficients to default values.
303       """       """
304       for i in self.COEFFICIENTS.iterkeys():       for i in self.COEFFICIENTS.iterkeys():
305           self.COEFFICIENTS[i].resetValue()           self.COEFFICIENTS[i].resetValue()
306    
307     def createNewCoefficient(self,name):     def createNewCoefficient(self,name):
308       """       """
309       @brief returns a new coefficient appropriate for coefficient name:       returns a new coefficient appropriate for coefficient name:
310       """       """
311       return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))       return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
312                
313    
314     def getShapeOfCoefficient(self,name):     def getShapeOfCoefficient(self,name):
315       """       """
316       @brief return the shape of the coefficient name       return the shape of the coefficient name
317    
318       @param name       @param name:
319       """       """
320       if self.hasCoefficient(name):       if self.hasCoefficient(name):
321          return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())          return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())
# Line 303  class LinearPDE: Line 324  class LinearPDE:
324    
325     def getFunctionSpaceForCoefficient(self,name):     def getFunctionSpaceForCoefficient(self,name):
326       """       """
327       @brief return the atoms of the coefficient name       return the atoms of the coefficient name
328    
329       @param name       @param name:
330       """       """
331       if self.hasCoefficient(name):       if self.hasCoefficient(name):
332          return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())          return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
# Line 314  class LinearPDE: Line 335  class LinearPDE:
335    
336     def alteredCoefficient(self,name):     def alteredCoefficient(self,name):
337       """       """
338       @brief annonced that coefficient name has been changed       announce that coefficient name has been changed
339    
340       @param name       @param name:
341       """       """
342       if self.hasCoefficient(name):       if self.hasCoefficient(name):
343          if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator()          if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator()
# Line 327  class LinearPDE: Line 348  class LinearPDE:
348     # ===== debug ==============================================================     # ===== debug ==============================================================
349     def setDebugOn(self):     def setDebugOn(self):
350         """         """
        @brief  
351         """         """
352         self.__debug=not None         self.__debug=not None
353    
354     def setDebugOff(self):     def setDebugOff(self):
355         """         """
        @brief  
356         """         """
357         self.__debug=None         self.__debug=None
358    
359     def debug(self):     def debug(self):
360         """         """
361         @brief returns true if the PDE is in the debug mode         returns true if the PDE is in the debug mode
362         """         """
363         return self.__debug         return self.__debug
364    
365     #===== Lumping ===========================     #===== Lumping ===========================
366     def setLumpingOn(self):     def setLumpingOn(self):
367        """        """
368        @brief indicates to use matrix lumping        indicates to use matrix lumping
369        """        """
370        if not self.isUsingLumping():        if not self.isUsingLumping():
371           if self.debug() : print "PDE Debug: lumping is set on"           if self.debug() : print "PDE Debug: lumping is set on"
# Line 355  class LinearPDE: Line 374  class LinearPDE:
374    
375     def setLumpingOff(self):     def setLumpingOff(self):
376        """        """
377        @brief switches off matrix lumping        switches off matrix lumping
378        """        """
379        if self.isUsingLumping():        if self.isUsingLumping():
380           if self.debug() : print "PDE Debug: lumping is set off"           if self.debug() : print "PDE Debug: lumping is set off"
# Line 364  class LinearPDE: Line 383  class LinearPDE:
383    
384     def setLumping(self,flag=False):     def setLumping(self,flag=False):
385        """        """
386        @brief set the matrix lumping flag to flag        set the matrix lumping flag to flag
387        """        """
388        if flag:        if flag:
389           self.setLumpingOn()           self.setLumpingOn()
# Line 373  class LinearPDE: Line 392  class LinearPDE:
392    
393     def isUsingLumping(self):     def isUsingLumping(self):
394        """        """
395        @brief        
396        """        """
397        return self.__lumping        return self.__lumping
398    
399     #============ method business =========================================================     #============ method business =========================================================
400     def setSolverMethod(self,solver=util.DEFAULT_METHOD):     def setSolverMethod(self,solver=util.DEFAULT_METHOD):
401         """         """
402         @brief sets a new solver         sets a new solver
403         """         """
404         if not solver==self.getSolverMethod():         if not solver==self.getSolverMethod():
405             self.__solver_method=solver             self.__solver_method=solver
# Line 389  class LinearPDE: Line 408  class LinearPDE:
408    
409     def getSolverMethod(self):     def getSolverMethod(self):
410         """         """
411         @brief returns the solver method         returns the solver method
412         """         """
413         return self.__solver_method         return self.__solver_method
414    
415     #============ tolerance business =========================================================     #============ tolerance business =========================================================
416     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
417         """         """
418         @brief resets the tolerance to tol.         resets the tolerance to tol.
419         """         """
420         if not tol>0:         if not tol>0:
421             raise ValueException,"Tolerance as to be positive"             raise ValueException,"Tolerance as to be positive"
# Line 406  class LinearPDE: Line 425  class LinearPDE:
425         return         return
426     def getTolerance(self):     def getTolerance(self):
427         """         """
428         @brief returns the tolerance set for the solution         returns the tolerance set for the solution
429         """         """
430         return self.__tolerance         return self.__tolerance
431    
432     #===== symmetry  flag ==========================     #===== symmetry  flag ==========================
433     def isSymmetric(self):     def isSymmetric(self):
434        """        """
435        @brief returns true is the operator is considered to be symmetric        returns true is the operator is considered to be symmetric
436        """        """
437        return self.__sym        return self.__sym
438    
439     def setSymmetryOn(self):     def setSymmetryOn(self):
440        """        """
441        @brief sets the symmetry flag to true        sets the symmetry flag to true
442        """        """
443        if not self.isSymmetric():        if not self.isSymmetric():
444           if self.debug() : print "PDE Debug: Operator is set to be symmetric"           if self.debug() : print "PDE Debug: Operator is set to be symmetric"
# Line 428  class LinearPDE: Line 447  class LinearPDE:
447    
448     def setSymmetryOff(self):     def setSymmetryOff(self):
449        """        """
450        @brief sets the symmetry flag to false        sets the symmetry flag to false
451        """        """
452        if self.isSymmetric():        if self.isSymmetric():
453           if self.debug() : print "PDE Debug: Operator is set to be unsymmetric"           if self.debug() : print "PDE Debug: Operator is set to be unsymmetric"
# Line 437  class LinearPDE: Line 456  class LinearPDE:
456    
457     def setSymmetryTo(self,flag=False):     def setSymmetryTo(self,flag=False):
458       """       """
459       @brief sets the symmetry flag to flag       sets the symmetry flag to flag
460    
461       @param flag       @param flag:
462       """       """
463       if flag:       if flag:
464          self.setSymmetryOn()          self.setSymmetryOn()
# Line 449  class LinearPDE: Line 468  class LinearPDE:
468     #===== order reduction ==========================     #===== order reduction ==========================
469     def setReducedOrderOn(self):     def setReducedOrderOn(self):
470       """       """
471       @brief switches to on reduced order       switches to on reduced order
472       """       """
473       self.setReducedOrderForSolutionOn()       self.setReducedOrderForSolutionOn()
474       self.setReducedOrderForEquationOn()       self.setReducedOrderForEquationOn()
475    
476     def setReducedOrderOff(self):     def setReducedOrderOff(self):
477       """       """
478       @brief switches to full order       switches to full order
479       """       """
480       self.setReducedOrderForSolutionOff()       self.setReducedOrderForSolutionOff()
481       self.setReducedOrderForEquationOff()       self.setReducedOrderForEquationOff()
482    
483     def setReducedOrderTo(self,flag=False):     def setReducedOrderTo(self,flag=False):
484       """       """
485       @brief sets order according to flag       sets order according to flag
486    
487       @param flag       @param flag:
488       """       """
489       self.setReducedOrderForSolutionTo(flag)       self.setReducedOrderForSolutionTo(flag)
490       self.setReducedOrderForEquationTo(flag)       self.setReducedOrderForEquationTo(flag)
# Line 474  class LinearPDE: Line 493  class LinearPDE:
493     #===== order reduction solution ==========================     #===== order reduction solution ==========================
494     def setReducedOrderForSolutionOn(self):     def setReducedOrderForSolutionOn(self):
495       """       """
496       @brief switches to reduced order to interpolate solution       switches to reduced order to interpolate solution
497       """       """
498       new_fs=escript.ReducedSolution(self.getDomain())       new_fs=escript.ReducedSolution(self.getDomain())
499       if self.getFunctionSpaceForSolution()!=new_fs:       if self.getFunctionSpaceForSolution()!=new_fs:
# Line 484  class LinearPDE: Line 503  class LinearPDE:
503    
504     def setReducedOrderForSolutionOff(self):     def setReducedOrderForSolutionOff(self):
505       """       """
506       @brief switches to full order to interpolate solution       switches to full order to interpolate solution
507       """       """
508       new_fs=escript.Solution(self.getDomain())       new_fs=escript.Solution(self.getDomain())
509       if self.getFunctionSpaceForSolution()!=new_fs:       if self.getFunctionSpaceForSolution()!=new_fs:
# Line 494  class LinearPDE: Line 513  class LinearPDE:
513    
514     def setReducedOrderForSolutionTo(self,flag=False):     def setReducedOrderForSolutionTo(self,flag=False):
515       """       """
516       @brief sets order for test functions according to flag       sets order for test functions according to flag
517    
518       @param flag       @param flag:
519       """       """
520       if flag:       if flag:
521          self.setReducedOrderForSolutionOn()          self.setReducedOrderForSolutionOn()
# Line 506  class LinearPDE: Line 525  class LinearPDE:
525     #===== order reduction equation ==========================     #===== order reduction equation ==========================
526     def setReducedOrderForEquationOn(self):     def setReducedOrderForEquationOn(self):
527       """       """
528       @brief switches to reduced order for test functions       switches to reduced order for test functions
529       """       """
530       new_fs=escript.ReducedSolution(self.getDomain())       new_fs=escript.ReducedSolution(self.getDomain())
531       if self.getFunctionSpaceForEquation()!=new_fs:       if self.getFunctionSpaceForEquation()!=new_fs:
# Line 516  class LinearPDE: Line 535  class LinearPDE:
535    
536     def setReducedOrderForEquationOff(self):     def setReducedOrderForEquationOff(self):
537       """       """
538       @brief switches to full order for test functions       switches to full order for test functions
539       """       """
540       new_fs=escript.Solution(self.getDomain())       new_fs=escript.Solution(self.getDomain())
541       if self.getFunctionSpaceForEquation()!=new_fs:       if self.getFunctionSpaceForEquation()!=new_fs:
# Line 526  class LinearPDE: Line 545  class LinearPDE:
545    
546     def setReducedOrderForEquationTo(self,flag=False):     def setReducedOrderForEquationTo(self,flag=False):
547       """       """
548       @brief sets order for test functions according to flag       sets order for test functions according to flag
549    
550       @param flag       @param flag:
551       """       """
552       if flag:       if flag:
553          self.setReducedOrderForEquationOn()          self.setReducedOrderForEquationOn()
# Line 538  class LinearPDE: Line 557  class LinearPDE:
557     # ==== initialization =====================================================================     # ==== initialization =====================================================================
558     def __getNewOperator(self):     def __getNewOperator(self):
559         """         """
        @brief  
560         """         """
561         return self.getDomain().newOperator( \         return self.getDomain().newOperator( \
562                             self.getNumEquations(), \                             self.getNumEquations(), \
# Line 549  class LinearPDE: Line 567  class LinearPDE:
567    
568     def __makeFreshRightHandSide(self):     def __makeFreshRightHandSide(self):
569         """         """
        @brief  
570         """         """
571         if self.debug() : print "PDE Debug: New right hand side allocated"         if self.debug() : print "PDE Debug: New right hand side allocated"
572         if self.getNumEquations()>1:         if self.getNumEquations()>1:
# Line 560  class LinearPDE: Line 577  class LinearPDE:
577    
578     def __getNewSolution(self):     def __getNewSolution(self):
579         """         """
        @brief  
580         """         """
581         if self.debug() : print "PDE Debug: New right hand side allocated"         if self.debug() : print "PDE Debug: New right hand side allocated"
582         if self.getNumSolutions()>1:         if self.getNumSolutions()>1:
# Line 570  class LinearPDE: Line 586  class LinearPDE:
586    
587     def __makeFreshOperator(self):     def __makeFreshOperator(self):
588         """         """
        @brief  
589         """         """
590         if self.__operator.isEmpty():         if self.__operator.isEmpty():
591             self.__operator=self.__getNewOperator()             self.__operator=self.__getNewOperator()
# Line 584  class LinearPDE: Line 599  class LinearPDE:
599     #============ some serivice functions  =====================================================     #============ some serivice functions  =====================================================
600     def getDomain(self):     def getDomain(self):
601       """       """
602       @brief returns the domain of the PDE       returns the domain of the PDE
603       """       """
604       return self.__domain       return self.__domain
605    
606     def getDim(self):     def getDim(self):
607       """       """
608       @brief returns the spatial dimension of the PDE       returns the spatial dimension of the PDE
609       """       """
610       return self.getDomain().getDim()       return self.getDomain().getDim()
611    
612     def getNumEquations(self):     def getNumEquations(self):
613       """       """
614       @brief returns the number of equations       returns the number of equations
615       """       """
616       if self.__numEquations>0:       if self.__numEquations>0:
617           return self.__numEquations           return self.__numEquations
# Line 605  class LinearPDE: Line 620  class LinearPDE:
620    
621     def getNumSolutions(self):     def getNumSolutions(self):
622       """       """
623       @brief returns the number of unknowns       returns the number of unknowns
624       """       """
625       if self.__numSolutions>0:       if self.__numSolutions>0:
626          return self.__numSolutions          return self.__numSolutions
# Line 615  class LinearPDE: Line 630  class LinearPDE:
630    
631     def checkSymmetry(self,verbose=True):     def checkSymmetry(self,verbose=True):
632        """        """
633        @brief returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.        returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.
634        """        """
635        verbose=verbose or self.debug()        verbose=verbose or self.debug()
636        out=True        out=True
# Line 676  class LinearPDE: Line 691  class LinearPDE:
691    
692     def getFlux(self,u):     def getFlux(self,u):
693         """         """
694         @brief returns the flux J_ij for a given u         returns the flux J_ij for a given u
   
             J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}  
695    
696         @param u argument of the operator         \f[
697           J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}
698           \f]
699    
700           @param u: argument of the operator
701         """         """
702         raise SystemError,"getFlux is not implemented yet"         raise SystemError,"getFlux is not implemented yet"
703         return None         return None
704    
705     def applyOperator(self,u):     def applyOperator(self,u):
706         """         """
707         @brief applies the operator of the PDE to a given solution u in weak from         applies the operator of the PDE to a given solution u in weak from
   
        @param u argument of the operator  
708    
709           @param u: argument of the operator
710         """         """
711         return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())         return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
712                                                                                                                                                                                                                                                                                                                        
713     def getResidual(self,u):     def getResidual(self,u):
714         """         """
715         @brief return the residual of u in the weak from         return the residual of u in the weak from
716    
717         @param u         @param u:
718         """         """
719         return self.applyOperator(u)-self.getRightHandSide()         return self.applyOperator(u)-self.getRightHandSide()
720    
721     def __setValue(self,**coefficients):     def __setValue(self,**coefficients):
722        """        """
723        @brief sets new values to coefficient        sets new values to coefficient
724    
725        @param coefficients        @param coefficients:
726        """        """
727        # check if the coefficients are  legal:        # check if the coefficients are  legal:
728        for i in coefficients.iterkeys():        for i in coefficients.iterkeys():
# Line 757  class LinearPDE: Line 772  class LinearPDE:
772    
773     def __setHomogeneousConstraintFlag(self):     def __setHomogeneousConstraintFlag(self):
774        """        """
775        @brief checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.        checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.
776        """        """
777        self.__homogeneous_constraint=True        self.__homogeneous_constraint=True
778        q=self.getCoefficientOfPDE("q")        q=self.getCoefficientOfPDE("q")
# Line 774  class LinearPDE: Line 789  class LinearPDE:
789     # ==== rebuild switches =====================================================================     # ==== rebuild switches =====================================================================
790     def __rebuildSolution(self,deep=False):     def __rebuildSolution(self,deep=False):
791         """         """
792         @brief indicates the PDE has to be reolved if the solution is requested         indicates the PDE has to be reolved if the solution is requested
793         """         """
794         if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."         if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."
795         self.__solution_isValid=False         self.__solution_isValid=False
# Line 783  class LinearPDE: Line 798  class LinearPDE:
798    
799     def __rebuildOperator(self,deep=False):     def __rebuildOperator(self,deep=False):
800         """         """
801         @brief indicates the operator has to be rebuilt next time it is used         indicates the operator has to be rebuilt next time it is used
802         """         """
803         if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."         if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."
804         self.__rebuildSolution(deep)         self.__rebuildSolution(deep)
# Line 792  class LinearPDE: Line 807  class LinearPDE:
807    
808     def __rebuildRightHandSide(self,deep=False):     def __rebuildRightHandSide(self,deep=False):
809         """         """
810         @brief indicates the right hand side has to be rebuild next time it is used         indicates the right hand side has to be rebuild next time it is used
811         """         """
812         if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."         if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."
813         self.__rebuildSolution(deep)         self.__rebuildSolution(deep)
# Line 801  class LinearPDE: Line 816  class LinearPDE:
816    
817     def __rebuildSystem(self,deep=False):     def __rebuildSystem(self,deep=False):
818         """         """
819         @brief annonced that all coefficient name has been changed         annonced that all coefficient name has been changed
820         """         """
821         self.__rebuildSolution(deep)         self.__rebuildSolution(deep)
822         self.__rebuildOperator(deep)         self.__rebuildOperator(deep)
# Line 809  class LinearPDE: Line 824  class LinearPDE:
824        
825     def __checkMatrixType(self):     def __checkMatrixType(self):
826       """       """
827       @brief reassess the matrix type and, if needed, initiates an operator rebuild       reassess the matrix type and, if needed, initiates an operator rebuild
828       """       """
829       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())
830       if not new_matrix_type==self.__matrix_type:       if not new_matrix_type==self.__matrix_type:
# Line 820  class LinearPDE: Line 835  class LinearPDE:
835     #============ assembling =======================================================     #============ assembling =======================================================
836     def __copyConstraint(self):     def __copyConstraint(self):
837        """        """
838        @brief copies the constrint condition into u        copies the constrint condition into u
839        """        """
840        if not self.__righthandside.isEmpty():        if not self.__righthandside.isEmpty():
841           q=self.getCoefficientOfPDE("q")           q=self.getCoefficientOfPDE("q")
# Line 834  class LinearPDE: Line 849  class LinearPDE:
849    
850     def __applyConstraint(self):     def __applyConstraint(self):
851         """         """
852         @brief applies the constraints defined by q and r to the system         applies the constraints defined by q and r to the system
853         """         """
854         q=self.getCoefficientOfPDE("q")         q=self.getCoefficientOfPDE("q")
855         r=self.getCoefficientOfPDE("r")         r=self.getCoefficientOfPDE("r")
# Line 856  class LinearPDE: Line 871  class LinearPDE:
871    
872     def getSystem(self):     def getSystem(self):
873         """         """
874         @brief return the operator and right hand side of the PDE         return the operator and right hand side of the PDE
875         """         """
876         if not self.__operator_isValid or not self.__righthandside_isValid:         if not self.__operator_isValid or not self.__righthandside_isValid:
877            if self.isUsingLumping():            if self.isUsingLumping():
# Line 939  class LinearPDE: Line 954  class LinearPDE:
954         return (self.__operator,self.__righthandside)         return (self.__operator,self.__righthandside)
955     def getOperator(self):     def getOperator(self):
956         """         """
957         @brief returns the operator of the PDE         returns the operator of the PDE
958         """         """
959         return self.getSystem()[0]         return self.getSystem()[0]
960    
961     def getRightHandSide(self):     def getRightHandSide(self):
962         """         """
963         @brief returns the right hand side of the PDE         returns the right hand side of the PDE
964         """         """
965         return self.getSystem()[1]         return self.getSystem()[1]
966    
967     def solve(self,**options):     def solve(self,**options):
968        """        """
969        @brief solve the PDE        solve the PDE
970    
971        @param options        @param options:
972        """        """
973        mat,f=self.getSystem()        mat,f=self.getSystem()
974        if self.isUsingLumping():        if self.isUsingLumping():
# Line 968  class LinearPDE: Line 983  class LinearPDE:
983    
984     def getSolution(self,**options):     def getSolution(self,**options):
985         """         """
986         @brief returns the solution of the PDE         returns the solution of the PDE
987    
988         @param options         @param options:
989         """         """
990         if not self.__solution_isValid:         if not self.__solution_isValid:
991             if self.debug() : print "PDE Debug: PDE is resolved."             if self.debug() : print "PDE Debug: PDE is resolved."
# Line 980  class LinearPDE: Line 995  class LinearPDE:
995    
996    
997    
998  def ELMAN_RAMAGE(P): return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))  def ELMAN_RAMAGE(P):
999         """   """
1000         return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
1001  def SIMPLIFIED_BROOK_HUGHES(P):  def SIMPLIFIED_BROOK_HUGHES(P):
1002           c=(P-3.).whereNegative()       """   """
1003           return P/6.*c+1./2.*(1.-c)       c=(P-3.).whereNegative()
1004  def HALF(P): return escript.Scalar(0.5,P.getFunctionSpace())       return P/6.*c+1./2.*(1.-c)
1005    
1006    def HALF(P):
1007        """ """
1008        return escript.Scalar(0.5,P.getFunctionSpace())
1009    
1010  class AdvectivePDE(LinearPDE):  class AdvectivePDE(LinearPDE):
1011     """     """
1012     @brief Class to handel a linear PDE domineated by advective terms:     Class to handle a linear PDE dominated by advective terms:
1013        
1014     class to define a linear PDE of the form     class to define a linear PDE of the form
1015    
1016       -(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     \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       with boundary conditons:     \f]
1019    
1020          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i     with boundary conditons:
1021    
1022      and contact conditions     \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    
1026          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i     and contact conditions
1027    
1028      and constraints:     \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    
1032           u_i=r_i where q_i>0     and constraints:
1033    
1034       \f[
1035       u_i=r_i \quad \mathrm{where} \quad q_i>0
1036       \f]
1037     """     """
1038     def __init__(self,domain,numEquations=0,numSolutions=0,xi=ELMAN_RAMAGE):     def __init__(self,domain,numEquations=0,numSolutions=0,xi=ELMAN_RAMAGE):
1039        LinearPDE.__init__(self,domain,numEquations,numSolutions)        LinearPDE.__init__(self,domain,numEquations,numSolutions)
# Line 1071  class AdvectivePDE(LinearPDE): Line 1098  class AdvectivePDE(LinearPDE):
1098    
1099     def getCoefficientOfPDE(self,name):     def getCoefficientOfPDE(self,name):
1100       """       """
1101       @brief return the value of the coefficient name of the general PDE       return the value of the coefficient name of the general PDE
1102       @param name  
1103         @param name:
1104       """       """
1105       if not self.getNumEquations() == self.getNumSolutions():       if not self.getNumEquations() == self.getNumSolutions():
1106            raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."            raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."
# Line 1212  class AdvectivePDE(LinearPDE): Line 1240  class AdvectivePDE(LinearPDE):
1240    
1241  class Poisson(LinearPDE):  class Poisson(LinearPDE):
1242     """     """
1243     @brief Class to define a Poisson equstion problem:     Class to define a Poisson equstion problem:
1244                                                                                                                                                                
1245     class to define a linear PDE of the form     class to define a linear PDE of the form
1246                                                                                                                                                                   \f[
1247          -u_{,jj} = f     -u_{,jj} = f
1248                                                                                                                                                                   \f]
1249       with boundary conditons:  
1250                                                                                                                                                                   with boundary conditons:
1251          n_j*u_{,j} = 0  
1252                                                                                                                                                                   \f[
1253      and constraints:     n_j*u_{,j} = 0
1254                                                                                                                                                                   \f]
1255           u=0 where q>0  
1256                                                                                                                                                                   and constraints:
1257    
1258       \f[
1259       u=0 \quad \mathrm{where} \quad q>0
1260       \f]
1261     """     """
1262    
1263     def __init__(self,domain,f=escript.Data(),q=escript.Data()):     def __init__(self,domain,f=escript.Data(),q=escript.Data()):
# Line 1241  class Poisson(LinearPDE): Line 1273  class Poisson(LinearPDE):
1273    
1274     def getCoefficientOfPDE(self,name):     def getCoefficientOfPDE(self,name):
1275       """       """
1276       @brief return the value of the coefficient name of the general PDE       return the value of the coefficient name of the general PDE
1277       @param name  
1278         @param name:
1279       """       """
1280       if name == "A" :       if name == "A" :
1281           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))
# Line 1272  class Poisson(LinearPDE): Line 1305  class Poisson(LinearPDE):
1305           raise SystemError,"unknown PDE coefficient %s",name           raise SystemError,"unknown PDE coefficient %s",name
1306    
1307  # $Log$  # $Log$
1308    # Revision 1.8  2005/06/09 05:37:59  jgs
1309    # Merge of development branch back to main trunk on 2005-06-09
1310    #
1311  # Revision 1.7  2005/05/06 04:26:10  jgs  # Revision 1.7  2005/05/06 04:26:10  jgs
1312  # Merge of development branch back to main trunk on 2005-05-06  # Merge of development branch back to main trunk on 2005-05-06
1313  #  #
1314    # Revision 1.1.2.23  2005/05/13 00:55:20  cochrane
1315    # Fixed up some docstrings.  Moved module-level functions to top of file so
1316    # that epydoc and doxygen can pick them up properly.
1317    #
1318    # Revision 1.1.2.22  2005/05/12 11:41:30  gross
1319    # some basic Models have been added
1320    #
1321    # Revision 1.1.2.21  2005/05/12 07:16:12  cochrane
1322    # Moved ELMAN_RAMAGE, SIMPLIFIED_BROOK_HUGHES, and HALF functions to bottom of
1323    # file so that the AdvectivePDE class is picked up by doxygen.  Some
1324    # reformatting of docstrings.  Addition of code to make equations come out
1325    # as proper LaTeX.
1326    #
1327  # Revision 1.1.2.20  2005/04/15 07:09:08  gross  # Revision 1.1.2.20  2005/04/15 07:09:08  gross
1328  # some problems with functionspace and linearPDEs fixed.  # some problems with functionspace and linearPDEs fixed.
1329  #  #

Legend:
Removed from v.121  
changed lines
  Added in v.122

  ViewVC Help
Powered by ViewVC 1.1.26