/[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 108 by jgs, Thu Jan 27 06:21:59 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 172  class LinearPDE: Line 191  class LinearPDE:
191     SSOR=util.SSOR     SSOR=util.SSOR
192     GMRES=util.GMRES     GMRES=util.GMRES
193     PRES20=util.PRES20     PRES20=util.PRES20
194       ILU0=util.ILU0
195       JACOBI=util.JACOBI
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 220  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 = getFunctionSpaceOfCoefficient(name))                           what = getFunctionSpaceForCoefficient(name))
249    
250     def __del__(self):     def __del__(self):
251       pass       pass
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 301  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 312  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 325  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 353  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 362  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 371  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 387  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 404  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 426  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 435  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 447  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 472  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 482  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 492  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 504  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 514  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 524  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 534  class LinearPDE: Line 555  class LinearPDE:
555          self.setReducedOrderForEquationOff()          self.setReducedOrderForEquationOff()
556                                                                                                                                                                                                                                                                                                                        
557     # ==== initialization =====================================================================     # ==== initialization =====================================================================
558     def __makeNewOperator(self):     def __getNewOperator(self):
559         """         """
        @brief  
560         """         """
561         return self.getDomain().newOperator( \         return self.getDomain().newOperator( \
562                             self.getNumEquations(), \                             self.getNumEquations(), \
# Line 545  class LinearPDE: Line 565  class LinearPDE:
565                             self.getFunctionSpaceForSolution(), \                             self.getFunctionSpaceForSolution(), \
566                             self.__matrix_type)                             self.__matrix_type)
567    
568     def __makeNewRightHandSide(self):     def __makeFreshRightHandSide(self):
569         """         """
        @brief  
570         """         """
571         return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)         if self.debug() : print "PDE Debug: New right hand side allocated"
572           if self.getNumEquations()>1:
573               self.__righthandside=escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
574           else:
575               self.__righthandside=escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
576           return self.__righthandside
577    
578     def __makeNewSolution(self):     def __getNewSolution(self):
579         """         """
        @brief  
580         """         """
581         return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)         if self.debug() : print "PDE Debug: New right hand side allocated"
582           if self.getNumSolutions()>1:
583               return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
584           else:
585               return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
586    
587     def __getFreshOperator(self):     def __makeFreshOperator(self):
588         """         """
        @brief  
589         """         """
590         if self.__operator.isEmpty():         if self.__operator.isEmpty():
591             self.__operator=self.__makeNewOperator()             self.__operator=self.__getNewOperator()
592             if self.debug() : print "PDE Debug: New operator allocated"             if self.debug() : print "PDE Debug: New operator allocated"
593         else:         else:
594             self.__operator.setValue(0.)             self.__operator.setValue(0.)
# Line 570  class LinearPDE: Line 596  class LinearPDE:
596             if self.debug() : print "PDE Debug: Operator reset to zero"             if self.debug() : print "PDE Debug: Operator reset to zero"
597         return self.__operator         return self.__operator
598    
    def __getFreshRightHandSide(self):  
        """  
        @brief  
        """  
        if self.__righthandside.isEmpty():  
            self.__righthandside=self.__makeNewRightHandSide()  
            if self.debug() : print "PDE Debug: New right hand side allocated"  
        else:  
            print "fix self.__righthandside*=0"  
            self.__righthandside*=0.  
            if self.debug() : print "PDE Debug: Right hand side reset to zero"  
        return  self.__righthandside  
   
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 607  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 617  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 656  class LinearPDE: Line 669  class LinearPDE:
669                 for i in range(self.getNumEquations()):                 for i in range(self.getNumEquations()):
670                     for j in range(self.getDim()):                     for j in range(self.getDim()):
671                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
672                           if util.Lsup(B[i,j,k]-C[i,j,k])>tol:                           if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
673                                if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,i,j,k)                                if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
674                                out=False                                out=False
675              else:              else:
676                 for j in range(self.getDim()):                 for j in range(self.getDim()):
# Line 678  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
695    
696              J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}         \f[
697           J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}
698         @param u argument of the operator         \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 759  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 776  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 785  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 794  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 803  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 811  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 822  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 836  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 844  class LinearPDE: Line 857  class LinearPDE:
857            # q is the row and column mask to indicate where constraints are set:            # q is the row and column mask to indicate where constraints are set:
858            row_q=escript.Data(q,self.getFunctionSpaceForEquation())            row_q=escript.Data(q,self.getFunctionSpaceForEquation())
859            col_q=escript.Data(q,self.getFunctionSpaceForSolution())            col_q=escript.Data(q,self.getFunctionSpaceForSolution())
860            u=self.__makeNewSolution()            u=self.__getNewSolution()
861            if r.isEmpty():            if r.isEmpty():
862               r_s=self.__makeNewSolution()               r_s=self.__getNewSolution()
863            else:            else:
864               r_s=escript.Data(r,self.getFunctionSpaceForSolution())               r_s=escript.Data(r,self.getFunctionSpaceForSolution())
865            u.copyWithMask(r_s,col_q)            u.copyWithMask(r_s,col_q)
# Line 858  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 872  class LinearPDE: Line 885  class LinearPDE:
885                   if not self.getCoefficientOfPDE("C").isEmpty():                   if not self.getCoefficientOfPDE("C").isEmpty():
886                            raise Warning,"Lumped matrix does not allow coefficient C"                            raise Warning,"Lumped matrix does not allow coefficient C"
887                   if self.debug() : print "PDE Debug: New lumped operator is built."                   if self.debug() : print "PDE Debug: New lumped operator is built."
888                   mat=self.__makeNewOperator()                   mat=self.__getNewOperator()
889                   self.getDomain().addPDEToSystem(mat,escript.Data(), \                   self.getDomain().addPDEToSystem(mat,escript.Data(), \
890                             self.getCoefficientOfPDE("A"), \                             self.getCoefficientOfPDE("A"), \
891                             self.getCoefficientOfPDE("B"), \                             self.getCoefficientOfPDE("B"), \
# Line 889  class LinearPDE: Line 902  class LinearPDE:
902                   self.__operator_isValid=True                   self.__operator_isValid=True
903                if not self.__righthandside_isValid:                if not self.__righthandside_isValid:
904                   if self.debug() : print "PDE Debug: New right hand side is built."                   if self.debug() : print "PDE Debug: New right hand side is built."
905                   self.getDomain().addPDEToRHS(self.__getFreshRightHandSide(), \                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
906                                 self.getCoefficientOfPDE("X"), \                                 self.getCoefficientOfPDE("X"), \
907                                 self.getCoefficientOfPDE("Y"),\                                 self.getCoefficientOfPDE("Y"),\
908                                 self.getCoefficientOfPDE("y"),\                                 self.getCoefficientOfPDE("y"),\
# Line 899  class LinearPDE: Line 912  class LinearPDE:
912            else:            else:
913               if not self.__operator_isValid and not self.__righthandside_isValid:               if not self.__operator_isValid and not self.__righthandside_isValid:
914                   if self.debug() : print "PDE Debug: New system is built."                   if self.debug() : print "PDE Debug: New system is built."
915                   self.getDomain().addPDEToSystem(self.__getFreshOperator(),self.__getFreshRightHandSide(), \                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
916                                 self.getCoefficientOfPDE("A"), \                                 self.getCoefficientOfPDE("A"), \
917                                 self.getCoefficientOfPDE("B"), \                                 self.getCoefficientOfPDE("B"), \
918                                 self.getCoefficientOfPDE("C"), \                                 self.getCoefficientOfPDE("C"), \
# Line 916  class LinearPDE: Line 929  class LinearPDE:
929                   self.__righthandside_isValid=True                   self.__righthandside_isValid=True
930               elif not self.__righthandside_isValid:               elif not self.__righthandside_isValid:
931                   if self.debug() : print "PDE Debug: New right hand side is built."                   if self.debug() : print "PDE Debug: New right hand side is built."
932                   self.getDomain().addPDEToRHS(self.__getFreshRightHandSide(), \                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
933                                 self.getCoefficientOfPDE("X"), \                                 self.getCoefficientOfPDE("X"), \
934                                 self.getCoefficientOfPDE("Y"),\                                 self.getCoefficientOfPDE("Y"),\
935                                 self.getCoefficientOfPDE("y"),\                                 self.getCoefficientOfPDE("y"),\
# Line 925  class LinearPDE: Line 938  class LinearPDE:
938                   self.__righthandside_isValid=True                   self.__righthandside_isValid=True
939               elif not self.__operator_isValid:               elif not self.__operator_isValid:
940                   if self.debug() : print "PDE Debug: New operator is built."                   if self.debug() : print "PDE Debug: New operator is built."
941                   self.getDomain().addPDEToSystem(self.__getFreshOperator(),escript.Data(), \                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
942                              self.getCoefficientOfPDE("A"), \                              self.getCoefficientOfPDE("A"), \
943                              self.getCoefficientOfPDE("B"), \                              self.getCoefficientOfPDE("B"), \
944                              self.getCoefficientOfPDE("C"), \                              self.getCoefficientOfPDE("C"), \
# Line 941  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 970  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 993  class LinearPDE:
993             self.__solution_isValid=True             self.__solution_isValid=True
994         return self.__solution         return self.__solution
995    
996    
997    
998    def ELMAN_RAMAGE(P):
999         """   """
1000         return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
1001    def SIMPLIFIED_BROOK_HUGHES(P):
1002         """   """
1003         c=(P-3.).whereNegative()
1004         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]
   
         n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i  
1019    
1020      and contact conditions     with boundary conditons:
1021    
1022          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[
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      and constraints:     and contact conditions
1027    
1028           u_i=r_i where q_i>0     \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      The PDE is solved by stabilizing the advective terms using SUPG approach:     and constraints:
1033    
1034         A_{ijkl}<-A_{ijkl}+0.5*h*(xi(b_{ik})*B_{ijk}*B_{ilk}/length(B_{i:k})^2)+0.5*h*xi_{c_{ik}}*(C_{ikj}*C_{ikl}/length(C_{ik:})^2)     \f[
1035       u_i=r_i \quad \mathrm{where} \quad q_i>0
1036      where       \f]
   
            b_{ik}=length(B_{i:k})*h/2/length(A_{i:k:})  
            c_{ik}=length(C_{i:k})*h/2/length(A_{i:k:})  
   
                       alpha/3        alpha<3  
            xi(alpha)=          for                  approximating cotanh(alpha)-1/alpha  
                        1             alpha>=3  
1037     """     """
1038     def __getXi(self,alpha):     def __init__(self,domain,numEquations=0,numSolutions=0,xi=ELMAN_RAMAGE):
1039           c=alpha-3.        LinearPDE.__init__(self,domain,numEquations,numSolutions)
1040           return c*c.whereNegative()/3.+1.        self.__xi=xi
1041          self.__Xi=escript.Data()
1042     def __getUpdateVector(self,V,hover2,alphaByU):  
1043       v=util.length(V)     def __calculateXi(self,peclet_factor,Z,h):
1044       v_max=util.Lsup(v)         Z_max=util.Lsup(Z)
1045       if v_max>0:         if Z_max>0.:
1046           V/=v+v_max*self.TOL            return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.TOL)
1047           alpha=alphaByU*v         else:
1048           A_bar=v*hover2*self.__getXi(alpha)            return 0.
          print "-------------"  
          print "@ max alpha ",util.Lsup(alpha)  
          print "-------------"  
      else:  
          A_bar=1.  
      return V,A_bar  
1049    
1050     def __getAlphaByU(self,A,hover2):     def setValue(self,**args):
1051        a=util.length(A)         if "A" in args.keys()   or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data()
1052        a_max=util.Lsup(a)         self._LinearPDE__setValue(**args)
1053        if a_max>0:            
1054           return hover2/(a+a_max*self.TOL)     def getXi(self):
1055        else:        if self.__Xi.isEmpty():
1056           return 1./self.TOL           B=self.getCoefficient("B")
1057             C=self.getCoefficient("C")
1058             A=self.getCoefficient("A")
1059             h=self.getDomain().getSize()
1060             self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
1061             if not C.isEmpty() or not B.isEmpty():
1062                if not C.isEmpty() and not B.isEmpty():
1063                    Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
1064                    if self.getNumEquations()>1:
1065                       if self.getNumSolutions()>1:
1066                          for i in range(self.getNumEquations()):
1067                             for k in range(self.getNumSolutions()):
1068                                for l in range(self.getDim()): Z2+=(C[i,k,l]-B[i,l,k])**2
1069                       else:
1070                          for i in range(self.getNumEquations()):
1071                             for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2
1072                    else:
1073                       if self.getNumSolutions()>1:
1074                          for k in range(self.getNumSolutions()):
1075                             for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2
1076                       else:
1077                          for l in range(self.getDim()): Z2+=(C[l]-B[l])**2
1078                    length_of_Z=util.sqrt(Z2)
1079                elif C.isEmpty():
1080                  length_of_Z=util.length(B)
1081                else:
1082                  length_of_Z=util.length(C)
1083    
1084                Z_max=util.Lsup(length_of_Z)
1085                if Z_max>0.:
1086                   length_of_A=util.length(A)
1087                   A_max=util.Lsup(length_of_A)
1088                   if A_max>0:
1089                        inv_A=1./(length_of_A+A_max*self.TOL)
1090                   else:
1091                        inv_A=1./self.TOL
1092                   peclet_number=length_of_Z*h/2*inv_A
1093                   xi=self.__xi(peclet_number)
1094                   self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL)
1095                   print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z)
1096          return self.__Xi
1097          
1098    
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():
1106              raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."
1107    
1108       if name == "A" :       if name == "A" :
1109           A=self.getCoefficient("A")           A=self.getCoefficient("A")
1110           B=self.getCoefficient("B")           B=self.getCoefficient("B")
1111           C=self.getCoefficient("C")           C=self.getCoefficient("C")
1112           if not B.isEmpty() or not C.isEmpty():           if B.isEmpty() and C.isEmpty():
1113               if A.isEmpty():              Aout=A
1114                   A=self.createNewCoefficient("A")           else:
1115               else:              if A.isEmpty():
1116                   A=A[:]                 Aout=self.createNewCoefficient("A")
1117               hover2=self.getDomain().getSize()/2.              else:
1118               if self.getNumEquations()>1:                 Aout=A[:]
1119                  if self.getNumSolutions()>1:              Xi=self.getXi()
1120                     for i in range(self.getNumEquations()):              if self.getNumEquations()>1:
1121                    for i in range(self.getNumEquations()):
1122                       for j in range(self.getDim()):
1123                        for k in range(self.getNumSolutions()):                        for k in range(self.getNumSolutions()):
1124                           alphaByU=self.__getAlphaByU(A[i,:,k,:],hover2)                           for l in range(self.getDim()):
1125                           if not B.isEmpty():                              if not C.isEmpty() and not B.isEmpty():
1126                               b_sub,f=self.__getUpdateVector(B[i,:,k],hover2,alphaByU)                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*(C[p,i,j]-B[p,j,i])*(C[p,k,l]-B[p,l,k])
1127                               for j in range(self.getDim()):                              elif C.isEmpty():
1128                                  for l in range(self.getDim()):                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]
1129                                     A[i,j,k,l]+=f*b_sub[j]*b_sub[l]                              else:
1130                           if not C.isEmpty():                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
1131                               c_sub,f=self.__getUpdateVector(C[i,k,:],hover2,alphaByU)              else:
1132                               for j in range(self.getDim()):                  for j in range(self.getDim()):
1133                                  for l in range(self.getDim()):                     for l in range(self.getDim()):
1134                                     A[i,j,k,l]+=f*c_sub[j]*c_sub[l]                        if not C.isEmpty() and not B.isEmpty():
1135                  else:                              Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])
1136                     for i in range(self.getNumEquations()):                        elif C.isEmpty():
1137                        alphaByU=self.__getAlphaByU(A[i,:,:],hover2)                            Aout[j,l]+=Xi*B[j]*B[l]
1138                        if not B.isEmpty():                        else:
1139                            b_sub,f=self.__getUpdateVector(B[i,:],hover2,alphaByU)                            Aout[j,l]+=Xi*C[j]*C[l]
1140                            for j in range(self.getDim()):           return Aout
                              for l in range(self.getDim()):  
                                  A[i,j,l]+=f*b_sub[j]*b_sub[l]  
                       if not C.isEmpty():  
                            c_sub,f=self.__getUpdateVector(C[i,:],hover2,alphaByU)  
                            for j in range(self.getDim()):  
                               for l in range(self.getDim()):  
                                  A[i,j,l]+=f*c_sub[j]*c_sub[l]  
              else:  
                 if self.getNumSolutions()>1:  
                    for k in range(self.getNumSolutions()):  
                       alphaByU=self.__getAlphaByU(A[:,k,:],hover2)  
                       if not B.isEmpty():  
                          b_sub,f=self.__getUpdateVector(B[:,k],hover2,alphaByU)  
                          for j in range(self.getDim()):  
                             for l in range(self.getDim()):  
                                    A[j,k,l]+=f*b_sub[j]*b_sub[l]  
                       if not C.isEmpty():  
                          c_sub,f=self.__getUpdateVector(C[k,:],hover2,alphaByU)  
                          for j in range(self.getDim()):  
                             for l in range(self.getDim()):  
                                A[j,k,l]+=f*c_sub[j]*c_sub[l]  
                 else:    
                    alphaByU=self.__getAlphaByU(A[:,:],hover2)  
                    if not B.isEmpty():  
                        b_sub,f=self.__getUpdateVector(B[:],hover2,alphaByU)  
                        for j in range(self.getDim()):  
                           for l in range(self.getDim()):  
                              A[j,l]+=f*b_sub[j]*b_sub[l]  
                    if not C.isEmpty():  
                       c_sub,f=self.__getUpdateVector(C[:],hover2,alphaByU)  
                       for j in range(self.getDim()):  
                           for l in range(self.getDim()):  
                              A[j,l]+=f*c_sub[j]*c_sub[l]  
          return A  
1141       elif name == "B" :       elif name == "B" :
1142           return self.getCoefficient("B")           B=self.getCoefficient("B")
1143             C=self.getCoefficient("C")
1144             D=self.getCoefficient("D")
1145             if C.isEmpty() or D.isEmpty():
1146                Bout=B
1147             else:
1148                Xi=self.getXi()
1149                if B.isEmpty():
1150                    Bout=self.createNewCoefficient("B")
1151                else:
1152                    Bout=B[:]
1153                if self.getNumEquations()>1:
1154                   for k in range(self.getNumSolutions()):
1155                      for p in range(self.getNumEquations()):
1156                         tmp=Xi*D[p,k]
1157                         for i in range(self.getNumEquations()):
1158                            for j in range(self.getDim()):
1159                               Bout[i,j,k]+=tmp*C[p,i,j]
1160                else:
1161                   tmp=Xi*D
1162                   for j in range(self.getDim()): Bout[j]+=tmp*C[j]
1163             return Bout
1164       elif name == "C" :       elif name == "C" :
1165           return self.getCoefficient("C")           B=self.getCoefficient("B")
1166             C=self.getCoefficient("C")
1167             D=self.getCoefficient("D")
1168             if B.isEmpty() or D.isEmpty():
1169                Cout=C
1170             else:
1171                Xi=self.getXi()
1172                if C.isEmpty():
1173                    Cout=self.createNewCoefficient("C")
1174                else:
1175                    Cout=C[:]
1176                if self.getNumEquations()>1:
1177                   for k in range(self.getNumSolutions()):
1178                       for p in range(self.getNumEquations()):
1179                          tmp=Xi*D[p,k]
1180                          for i in range(self.getNumEquations()):
1181                            for l in range(self.getDim()):
1182                                     Cout[i,k,l]+=tmp*B[p,l,i]
1183                else:
1184                   tmp=Xi*D
1185                   for j in range(self.getDim()): Cout[j]+=tmp*B[j]
1186             return Cout
1187       elif name == "D" :       elif name == "D" :
1188           return self.getCoefficient("D")           return self.getCoefficient("D")
1189       elif name == "X" :       elif name == "X" :
1190           return self.getCoefficient("X")           X=self.getCoefficient("X")
1191             Y=self.getCoefficient("Y")
1192             B=self.getCoefficient("B")
1193             C=self.getCoefficient("C")
1194             if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):
1195                Xout=X
1196             else:
1197                if X.isEmpty():
1198                    Xout=self.createNewCoefficient("X")
1199                else:
1200                    Xout=X[:]
1201                Xi=self.getXi()
1202                if self.getNumEquations()>1:
1203                     for p in range(self.getNumEquations()):
1204                        tmp=Xi*Y[p]
1205                        for i in range(self.getNumEquations()):
1206                           for j in range(self.getDim()):
1207                              if not C.isEmpty() and not B.isEmpty():
1208                                 Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
1209                              elif C.isEmpty():
1210                                 Xout[i,j]-=tmp*B[p,j,i]
1211                              else:
1212                                 Xout[i,j]+=tmp*C[p,i,j]
1213                else:
1214                     tmp=Xi*Y
1215                     for j in range(self.getDim()):
1216                        if not C.isEmpty() and not B.isEmpty():
1217                           Xout[j]+=tmp*(C[j]-B[j])
1218                        elif C.isEmpty():
1219                           Xout[j]-=tmp*B[j]
1220                        else:
1221                           Xout[j]+=tmp*C[j]
1222             return Xout
1223       elif name == "Y" :       elif name == "Y" :
1224           return self.getCoefficient("Y")           return self.getCoefficient("Y")
1225       elif name == "d" :       elif name == "d" :
# Line 1138  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 1163  class Poisson(LinearPDE): Line 1269  class Poisson(LinearPDE):
1269         self.setValue(f,q)         self.setValue(f,q)
1270    
1271     def setValue(self,f=escript.Data(),q=escript.Data()):     def setValue(self,f=escript.Data(),q=escript.Data()):
1272         self._setValue(f=f,q=q)         self._LinearPDE__setValue(f=f,q=q)
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 1196  class Poisson(LinearPDE): Line 1303  class Poisson(LinearPDE):
1303           return self.getCoefficient("q")           return self.getCoefficient("q")
1304       else:       else:
1305           raise SystemError,"unknown PDE coefficient %s",name           raise SystemError,"unknown PDE coefficient %s",name
1306    
1307    # $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
1312    # 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
1328    # some problems with functionspace and linearPDEs fixed.
1329    #
1330    # Revision 1.1.2.19  2005/03/04 05:27:07  gross
1331    # bug in SystemPattern fixed.
1332    #
1333    # Revision 1.1.2.18  2005/02/08 06:16:45  gross
1334    # Bugs in AdvectivePDE fixed, AdvectiveTest is stable but more testing is needed
1335    #
1336    # Revision 1.1.2.17  2005/02/08 05:56:19  gross
1337    # Reference Number handling added
1338    #
1339    # Revision 1.1.2.16  2005/02/07 04:41:28  gross
1340    # some function exposed to python to make mesh merging running
1341    #
1342    # Revision 1.1.2.15  2005/02/03 00:14:44  gross
1343    # timeseries add and ESySParameter.py renames esysXML.py for consistence
1344    #
1345    # Revision 1.1.2.14  2005/02/01 06:44:10  gross
1346    # new implementation of AdvectivePDE which now also updates right hand side. systems of PDEs are still not working
1347    #
1348    # Revision 1.1.2.13  2005/01/25 00:47:07  gross
1349    # updates in the documentation
1350    #
1351    # Revision 1.1.2.12  2005/01/12 01:28:04  matt
1352    # Added createCoefficient method for linearPDEs.
1353    #
1354    # Revision 1.1.2.11  2005/01/11 01:55:34  gross
1355    # a problem in linearPDE class fixed
1356    #
1357    # Revision 1.1.2.10  2005/01/07 01:13:29  gross
1358    # some bugs in linearPDE fixed
1359    #
1360    # Revision 1.1.2.9  2005/01/06 06:24:58  gross
1361    # some bugs in slicing fixed
1362    #
1363    # Revision 1.1.2.8  2005/01/05 04:21:40  gross
1364    # FunctionSpace checking/matchig in slicing added
1365    #
1366    # Revision 1.1.2.7  2004/12/29 10:03:41  gross
1367    # bug in setValue fixed
1368    #
1369    # Revision 1.1.2.6  2004/12/29 05:29:59  gross
1370    # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()
1371    #
1372    # Revision 1.1.2.5  2004/12/29 00:18:41  gross
1373    # AdvectivePDE added
1374    #
1375    # Revision 1.1.2.4  2004/12/24 06:05:41  gross
1376    # some changes in linearPDEs to add AdevectivePDE
1377    #
1378    # Revision 1.1.2.3  2004/12/16 00:12:34  gross
1379    # __init__ of LinearPDE does not accept any coefficient anymore
1380    #
1381    # Revision 1.1.2.2  2004/12/14 03:55:01  jgs
1382    # *** empty log message ***
1383    #
1384    # Revision 1.1.2.1  2004/12/12 22:53:47  gross
1385    # linearPDE has been renamed LinearPDE
1386    #
1387    # Revision 1.1.1.1.2.7  2004/12/07 10:13:08  gross
1388    # GMRES added
1389    #
1390    # Revision 1.1.1.1.2.6  2004/12/07 03:19:50  gross
1391    # options for GMRES and PRES20 added
1392    #
1393    # Revision 1.1.1.1.2.5  2004/12/01 06:25:15  gross
1394    # some small changes
1395    #
1396    # Revision 1.1.1.1.2.4  2004/11/24 01:50:21  gross
1397    # Finley solves 4M unknowns now
1398    #
1399    # Revision 1.1.1.1.2.3  2004/11/15 06:05:26  gross
1400    # poisson solver added
1401    #
1402    # Revision 1.1.1.1.2.2  2004/11/12 06:58:15  gross
1403    # 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
1404    #
1405    # Revision 1.1.1.1.2.1  2004/10/28 22:59:22  gross
1406    # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed
1407    #
1408    # Revision 1.1.1.1  2004/10/26 06:53:56  jgs
1409    # initial import of project esys2
1410    #
1411    # Revision 1.3.2.3  2004/10/26 06:43:48  jgs
1412    # committing Lutz's and Paul's changes to brach jgs
1413    #
1414    # Revision 1.3.4.1  2004/10/20 05:32:51  cochrane
1415    # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.
1416    #
1417    # Revision 1.3  2004/09/23 00:53:23  jgs
1418    # minor fixes
1419    #
1420    # Revision 1.1  2004/08/28 12:58:06  gross
1421    # SimpleSolve is not running yet: problem with == of functionsspace
1422    #
1423    #

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

  ViewVC Help
Powered by ViewVC 1.1.26