/[escript]/branches/symbolic_from_3470/escript/py_src/linearPDEs.py
ViewVC logotype

Diff of /branches/symbolic_from_3470/escript/py_src/linearPDEs.py

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

revision 2168 by gross, Wed Nov 26 08:13:00 2008 UTC revision 2169 by caltinay, Wed Dec 17 03:08:58 2008 UTC
# Line 21  __url__="http://www.uq.edu.au/esscc/escr Line 21  __url__="http://www.uq.edu.au/esscc/escr
21    
22  """  """
23  The module provides an interface to define and solve linear partial  The module provides an interface to define and solve linear partial
24  differential equations (PDEs) and Transport problems within L{escript}. L{linearPDEs} does not provide any  differential equations (PDEs) and Transport problems within L{escript}.
25  solver capabilities in itself but hands the PDE over to  L{linearPDEs} does not provide any solver capabilities in itself but hands the
26  the PDE solver library defined through the L{Domain<escript.Domain>} of the PDE.  PDE over to the PDE solver library defined through the L{Domain<escript.Domain>}
27  The general interface is provided through the L{LinearPDE} class.  of the PDE. The general interface is provided through the L{LinearPDE} class.
28  L{TransportProblem} provides an interface to initial value problems dominated by its advective terms.  L{TransportProblem} provides an interface to initial value problems dominated
29    by its advective terms.
30    
31  @var __author__: name of author  @var __author__: name of author
32  @var __copyright__: copyrights  @var __copyright__: copyrights
# Line 45  __author__="Lutz Gross, l.gross@uq.edu.a Line 46  __author__="Lutz Gross, l.gross@uq.edu.a
46    
47  class IllegalCoefficient(ValueError):  class IllegalCoefficient(ValueError):
48     """     """
49     raised if an illegal coefficient of the general ar particular PDE is requested.     Exception that is raised if an illegal coefficient of the general or
50       particular PDE is requested.
51     """     """
52     pass     pass
53    
54  class IllegalCoefficientValue(ValueError):  class IllegalCoefficientValue(ValueError):
55     """     """
56     raised if an incorrect value for a coefficient is used.     Exception that is raised if an incorrect value for a coefficient is used.
57     """     """
58     pass     pass
59    
60  class IllegalCoefficientFunctionSpace(ValueError):  class IllegalCoefficientFunctionSpace(ValueError):
61     """     """
62     raised if an incorrect function space for a coefficient is used.     Exception that is raised if an incorrect function space for a coefficient
63       is used.
64     """     """
65    
66  class UndefinedPDEError(ValueError):  class UndefinedPDEError(ValueError):
67     """     """
68     raised if a PDE is not fully defined yet.     Exception that is raised if a PDE is not fully defined yet.
69     """     """
70     pass     pass
71    
72  class PDECoef(object):  class PDECoef(object):
73      """      """
74      A class for describing a PDE coefficient      A class for describing a PDE coefficient.
75    
76      @cvar INTERIOR: indicator that coefficient is defined on the interior of the PDE domain      @cvar INTERIOR: indicator that coefficient is defined on the interior of
77      @cvar BOUNDARY: indicator that coefficient is defined on the boundary of the PDE domain                      the PDE domain
78      @cvar CONTACT: indicator that coefficient is defined on the contact region within the PDE domain      @cvar BOUNDARY: indicator that coefficient is defined on the boundary of
79      @cvar INTERIOR_REDUCED: indicator that coefficient is defined on the interior of the PDE domain using a reduced integration order                      the PDE domain
80      @cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the boundary of the PDE domain using a reduced integration order      @cvar CONTACT: indicator that coefficient is defined on the contact region
81      @cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact region within the PDE domain using a reduced integration order                     within the PDE domain
82      @cvar SOLUTION: indicator that coefficient is defined trough a solution of the PDE      @cvar INTERIOR_REDUCED: indicator that coefficient is defined on the
83      @cvar REDUCED: indicator that coefficient is defined trough a reduced solution of the PDE                              interior of the PDE domain using a reduced
84      @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is defined by the number PDE equations                              integration order
85      @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is defined by the number PDE solutions      @cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the
86      @cvar BY_DIM: indicator that the dimension of the coefficient shape is defined by the spatial dimension                              boundary of the PDE domain using a reduced
87      @cvar OPERATOR: indicator that the the coefficient alters the operator of the PDE                              integration order
88      @cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right hand side of the PDE      @cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact
89      @cvar BOTH: indicator that the the coefficient alters the operator as well as the right hand side of the PDE                             region within the PDE domain using a reduced
90                               integration order
91        @cvar SOLUTION: indicator that coefficient is defined through a solution of
92                        the PDE
93        @cvar REDUCED: indicator that coefficient is defined through a reduced
94                       solution of the PDE
95        @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is
96                           defined by the number of PDE equations
97        @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is
98                           defined by the number of PDE solutions
99        @cvar BY_DIM: indicator that the dimension of the coefficient shape is
100                      defined by the spatial dimension
101        @cvar OPERATOR: indicator that the the coefficient alters the operator of
102                        the PDE
103        @cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right
104                             hand side of the PDE
105        @cvar BOTH: indicator that the the coefficient alters the operator as well
106                    as the right hand side of the PDE
107    
108      """      """
109      INTERIOR=0      INTERIOR=0
# Line 103  class PDECoef(object): Line 123  class PDECoef(object):
123    
124      def __init__(self, where, pattern, altering):      def __init__(self, where, pattern, altering):
125         """         """
126         Initialise a PDE Coefficient type         Initialises a PDE coefficient type.
127    
128         @param where: describes where the coefficient lives         @param where: describes where the coefficient lives
129         @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION}, L{REDUCED},         @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION},
130                             L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED}, L{CONTACT_REDUCED}.                      L{REDUCED}, L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED},
131         @param pattern: describes the shape of the coefficient and how the shape is build for a given                      L{CONTACT_REDUCED}
132                spatial dimension and numbers of equation and solution in then PDE. For instance,         @param pattern: describes the shape of the coefficient and how the shape
133                (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) descrbes a rank 3 coefficient which                         is built for a given spatial dimension and numbers of
134                is instanciated as shape (3,2,2) in case of a three equations and two solution components                         equations and solutions in then PDE. For instance,
135                on a 2-dimensional domain. In the case of single equation and a single solution component                         (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) describes a
136                the shape compoments marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped. In this case                         rank 3 coefficient which is instantiated as shape (3,2,2)
137                the example would be read as (2,).                         in case of three equations and two solution components
138                           on a 2-dimensional domain. In the case of single equation
139                           and a single solution component the shape components
140                           marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped.
141                           In this case the example would be read as (2,).
142         @type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM}         @type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM}
143         @param altering: indicates what part of the PDE is altered if the coefficiennt is altered         @param altering: indicates what part of the PDE is altered if the
144                            coefficient is altered
145         @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}         @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}
146         """         """
147         super(PDECoef, self).__init__()         super(PDECoef, self).__init__()
# Line 127  class PDECoef(object): Line 152  class PDECoef(object):
152    
153      def resetValue(self):      def resetValue(self):
154         """         """
155         resets coefficient value to default         Resets the coefficient value to the default.
156         """         """
157         self.value=escript.Data()         self.value=escript.Data()
158    
159      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
160         """         """
161         defines the L{FunctionSpace<escript.FunctionSpace>} of the coefficient         Returns the L{FunctionSpace<escript.FunctionSpace>} of the coefficient.
162    
163         @param domain: domain on which the PDE uses the coefficient         @param domain: domain on which the PDE uses the coefficient
164         @type domain: L{Domain<escript.Domain>}         @type domain: L{Domain<escript.Domain>}
165         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation         @param reducedEquationOrder: True to indicate that reduced order is used
166                                        to represent the equation
167         @type reducedEquationOrder: C{bool}         @type reducedEquationOrder: C{bool}
168         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution         @param reducedSolutionOrder: True to indicate that reduced order is used
169                                        to represent the solution
170         @type reducedSolutionOrder: C{bool}         @type reducedSolutionOrder: C{bool}
171         @return:  L{FunctionSpace<escript.FunctionSpace>} of the coefficient         @return: L{FunctionSpace<escript.FunctionSpace>} of the coefficient
172         @rtype:  L{FunctionSpace<escript.FunctionSpace>}         @rtype: L{FunctionSpace<escript.FunctionSpace>}
173         """         """
174         if self.what==self.INTERIOR:         if self.what==self.INTERIOR:
175              return escript.Function(domain)              return escript.Function(domain)
# Line 166  class PDECoef(object): Line 193  class PDECoef(object):
193    
194      def getValue(self):      def getValue(self):
195         """         """
196         returns the value of the coefficient         Returns the value of the coefficient.
197    
198         @return:  value of the coefficient         @return: value of the coefficient
199         @rtype:  L{Data<escript.Data>}         @rtype: L{Data<escript.Data>}
200         """         """
201         return self.value         return self.value
202    
203      def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):      def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
204         """         """
205         set the value of the coefficient to a new value         Sets the value of the coefficient to a new value.
206    
207         @param domain: domain on which the PDE uses the coefficient         @param domain: domain on which the PDE uses the coefficient
208         @type domain: L{Domain<escript.Domain>}         @type domain: L{Domain<escript.Domain>}
# Line 183  class PDECoef(object): Line 210  class PDECoef(object):
210         @type numEquations: C{int}         @type numEquations: C{int}
211         @param numSolutions: number of components of the PDE solution         @param numSolutions: number of components of the PDE solution
212         @type numSolutions: C{int}         @type numSolutions: C{int}
213         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation         @param reducedEquationOrder: True to indicate that reduced order is used
214                                        to represent the equation
215         @type reducedEquationOrder: C{bool}         @type reducedEquationOrder: C{bool}
216         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution         @param reducedSolutionOrder: True to indicate that reduced order is used
217                                        to represent the solution
218         @type reducedSolutionOrder: C{bool}         @type reducedSolutionOrder: C{bool}
219         @param newValue: number of components of the PDE solution         @param newValue: number of components of the PDE solution
220         @type newValue: any object that can be converted into a L{Data<escript.Data>} object with the appropriate shape and L{FunctionSpace<escript.FunctionSpace>}         @type newValue: any object that can be converted into a
221         @raise IllegalCoefficientValue: if the shape of the assigned value does not match the shape of the coefficient                         L{Data<escript.Data>} object with the appropriate shape
222         @raise IllegalCoefficientFunctionSpace: if unable to interploate value to appropriate function space                         and L{FunctionSpace<escript.FunctionSpace>}
223           @raise IllegalCoefficientValue: if the shape of the assigned value does
224                                           not match the shape of the coefficient
225           @raise IllegalCoefficientFunctionSpace: if unable to interpolate value
226                                                   to appropriate function space
227         """         """
228         if newValue==None:         if newValue==None:
229             newValue=escript.Data()             newValue=escript.Data()
# Line 210  class PDECoef(object): Line 243  class PDECoef(object):
243    
244      def isAlteringOperator(self):      def isAlteringOperator(self):
245          """          """
246          checks if the coefficient alters the operator of the PDE          Checks if the coefficient alters the operator of the PDE.
247    
248          @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
249          @rtype:  C{bool}                   coefficient is changed
250      """          @rtype: C{bool}
251            """
252          if self.altering==self.OPERATOR or self.altering==self.BOTH:          if self.altering==self.OPERATOR or self.altering==self.BOTH:
253              return not None              return not None
254          else:          else:
# Line 222  class PDECoef(object): Line 256  class PDECoef(object):
256    
257      def isAlteringRightHandSide(self):      def isAlteringRightHandSide(self):
258          """          """
259          checks if the coefficeint alters the right hand side of the PDE          Checks if the coefficient alters the right hand side of the PDE.
260    
261      @rtype:  C{bool}          @rtype: C{bool}
262          @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
263      """                   coefficient is changed, C{None} otherwise.
264            """
265          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
266              return not None              return not None
267          else:          else:
# Line 234  class PDECoef(object): Line 269  class PDECoef(object):
269    
270      def estimateNumEquationsAndNumSolutions(self,domain,shape=()):      def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
271         """         """
272         tries to estimate the number of equations and number of solutions if the coefficient has the given shape         Tries to estimate the number of equations and number of solutions if
273           the coefficient has the given shape.
274    
275         @param domain: domain on which the PDE uses the coefficient         @param domain: domain on which the PDE uses the coefficient
276         @type domain: L{Domain<escript.Domain>}         @type domain: L{Domain<escript.Domain>}
277         @param shape: suggested shape of the coefficient         @param shape: suggested shape of the coefficient
278         @type shape: C{tuple} of C{int} values         @type shape: C{tuple} of C{int} values
279         @return: the number of equations and number of solutions of the PDE is the coefficient has shape s.         @return: the number of equations and number of solutions of the PDE if
280                   If no appropriate numbers could be identified, C{None} is returned                  the coefficient has given shape. If no appropriate numbers
281                    could be identified, C{None} is returned
282         @rtype: C{tuple} of two C{int} values or C{None}         @rtype: C{tuple} of two C{int} values or C{None}
283         """         """
284         dim=domain.getDim()         dim=domain.getDim()
# Line 277  class PDECoef(object): Line 314  class PDECoef(object):
314               else:               else:
315                   if s==shape: return (None,u)                   if s==shape: return (None,u)
316         return None         return None
317    
318      def definesNumSolutions(self):      def definesNumSolutions(self):
319         """         """
320         checks if the coefficient allows to estimate the number of solution components         Checks if the coefficient allows to estimate the number of solution
321           components.
322    
323         @return: True if the coefficient allows an estimate of the number of solution components         @return: True if the coefficient allows an estimate of the number of
324                    solution components, False otherwise
325         @rtype: C{bool}         @rtype: C{bool}
326         """         """
327         for i in self.pattern:         for i in self.pattern:
# Line 290  class PDECoef(object): Line 330  class PDECoef(object):
330    
331      def definesNumEquation(self):      def definesNumEquation(self):
332         """         """
333         checks if the coefficient allows to estimate the number of equations         Checks if the coefficient allows to estimate the number of equations.
334    
335         @return: True if the coefficient allows an estimate of the number of equations         @return: True if the coefficient allows an estimate of the number of
336                    equations, False otherwise
337         @rtype: C{bool}         @rtype: C{bool}
338         """         """
339         for i in self.pattern:         for i in self.pattern:
# Line 301  class PDECoef(object): Line 342  class PDECoef(object):
342    
343      def __CompTuple2(self,t1,t2):      def __CompTuple2(self,t1,t2):
344        """        """
345        Compare two tuples of possible number of equations and number of solutions        Compares two tuples of possible number of equations and number of
346          solutions.
       @param t1: The first tuple  
       @param t2: The second tuple  
347    
348          @param t1: the first tuple
349          @param t2: the second tuple
350          @return: 0, 1, or -1
351        """        """
352    
353        dif=t1[0]+t1[1]-(t2[0]+t2[1])        dif=t1[0]+t1[1]-(t2[0]+t2[1])
# Line 315  class PDECoef(object): Line 357  class PDECoef(object):
357    
358      def getShape(self,domain,numEquations=1,numSolutions=1):      def getShape(self,domain,numEquations=1,numSolutions=1):
359         """         """
360         builds the required shape of the coefficient         Builds the required shape of the coefficient.
361    
362         @param domain: domain on which the PDE uses the coefficient         @param domain: domain on which the PDE uses the coefficient
363         @type domain: L{Domain<escript.Domain>}         @type domain: L{Domain<escript.Domain>}
# Line 339  class PDECoef(object): Line 381  class PDECoef(object):
381    
382  class LinearProblem(object):  class LinearProblem(object):
383     """     """
384     This is the base class is to define a general linear PDE-type problem for an u     This is the base class to define a general linear PDE-type problem for
385     for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.     for an unknown function M{u} on a given domain defined through a
386     The problem can be given as a single equations or as a system of equations.     L{Domain<escript.Domain>} object. The problem can be given as a single
387       equation or as a system of equations.
388     The class assumes that some sort of assembling process is required to form a problem of the form  
389       The class assumes that some sort of assembling process is required to form
390     M{L u=f}     a problem of the form
391      
392     where M{L} is an operator and M{f} is the right hand side. This operator problem will be solved to get the     M{L u=f}
    unknown M{u}.  
393    
394     @cvar DEFAULT: The default method used to solve the system of linear equations     where M{L} is an operator and M{f} is the right hand side. This operator
395       problem will be solved to get the unknown M{u}.
396    
397       @cvar DEFAULT: The default method used to solve the system of linear
398                      equations
399     @cvar DIRECT: The direct solver based on LDU factorization     @cvar DIRECT: The direct solver based on LDU factorization
400     @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)     @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be
401     @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)                      applied for symmetric PDEs)
402       @cvar PCG: The preconditioned conjugate gradient method (can only be applied
403                  for symmetric PDEs)
404     @cvar CR: The conjugate residual method     @cvar CR: The conjugate residual method
405     @cvar CGS: The conjugate gardient square method     @cvar CGS: The conjugate gradient square method
406     @cvar BICGSTAB: The stabilized BiConjugate Gradient method.     @cvar BICGSTAB: The stabilized BiConjugate Gradient method
407     @cvar TFQMR: Transport Free Quasi Minimal Residual method.     @cvar TFQMR: Transport Free Quasi Minimal Residual method
408     @cvar MINRES: Minimum residual method.     @cvar MINRES: Minimum residual method
409     @cvar SSOR: The symmetric overrealaxtion method     @cvar SSOR: The symmetric overrelaxation method
410     @cvar ILU0: The incomplete LU factorization preconditioner  with no fill in     @cvar ILU0: The incomplete LU factorization preconditioner with no fill-in
411     @cvar ILUT: The incomplete LU factorization preconditioner with will in     @cvar ILUT: The incomplete LU factorization preconditioner with fill-in
412     @cvar JACOBI: The Jacobi preconditioner     @cvar JACOBI: The Jacobi preconditioner
413     @cvar GMRES: The Gram-Schmidt minimum residual method     @cvar GMRES: The Gram-Schmidt minimum residual method
414     @cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals     @cvar PRES20: Special GMRES with restart after 20 steps and truncation after
415     @cvar LUMPING: Matrix lumping.                   5 residuals
416       @cvar LUMPING: Matrix lumping
417     @cvar NO_REORDERING: No matrix reordering allowed     @cvar NO_REORDERING: No matrix reordering allowed
418     @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization     @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
419     @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization     @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during
420                                factorization
421     @cvar PASO: PASO solver package     @cvar PASO: PASO solver package
422     @cvar SCSL: SGI SCSL solver library     @cvar SCSL: SGI SCSL solver library
423     @cvar MKL: Intel's MKL solver library     @cvar MKL: Intel's MKL solver library
424     @cvar UMFPACK: the UMFPACK library     @cvar UMFPACK: The UMFPACK library
425     @cvar TRILINOS: the TRILINOS parallel solver class library from Sandia Natl Labs     @cvar TRILINOS: The TRILINOS parallel solver class library from Sandia Natl
426                       Labs
427     @cvar ITERATIVE: The default iterative solver     @cvar ITERATIVE: The default iterative solver
428     @cvar AMG: algebraic multi grid     @cvar AMG: Algebraic Multi Grid
429     @cvar RILU: recursive ILU     @cvar RILU: Recursive ILU
430     @cvar GS: Gauss-Seidel solver     @cvar GS: Gauss-Seidel solver
431    
432     """     """
# Line 420  class LinearProblem(object): Line 470  class LinearProblem(object):
470    
471     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
472       """       """
473       initializes a linear problem       Initializes a linear problem.
474    
475       @param domain: domain of the PDE       @param domain: domain of the PDE
476       @type domain: L{Domain<escript.Domain>}       @type domain: L{Domain<escript.Domain>}
477       @param numEquations: number of equations. If numEquations==None the number of equations       @param numEquations: number of equations. If C{None} the number of
478                            is extracted from the coefficients.                            equations is extracted from the coefficients.
479       @param numSolutions: number of solution components. If  numSolutions==None the number of solution components       @param numSolutions: number of solution components. If C{None} the number
480                            is extracted from the coefficients.                            of solution components is extracted from the
481       @param debug: if True debug informations are printed.                            coefficients.
482         @param debug: if True debug information is printed
483    
484       """       """
485       super(LinearProblem, self).__init__()       super(LinearProblem, self).__init__()
# Line 452  class LinearProblem(object): Line 503  class LinearProblem(object):
503       self.resetAllCoefficients()       self.resetAllCoefficients()
504       self.__system_type=None       self.__system_type=None
505       self.updateSystemType()       self.updateSystemType()
506     # =============================================================================     # ==========================================================================
507     #    general stuff:     #    general stuff:
508     # =============================================================================     # ==========================================================================
509     def __str__(self):     def __str__(self):
510       """       """
511       returns string representation of the PDE       Returns a string representation of the PDE.
512    
513       @return: a simple representation of the PDE       @return: a simple representation of the PDE
514       @rtype: C{str}       @rtype: C{str}
515       """       """
516       return "<LinearProblem %d>"%id(self)       return "<LinearProblem %d>"%id(self)
517     # =============================================================================     # ==========================================================================
518     #    debug :     #    debug :
519     # =============================================================================     # ==========================================================================
520     def setDebugOn(self):     def setDebugOn(self):
521       """       """
522       switches on debugging       Switches debug output on.
523       """       """
524       self.__debug=not None       self.__debug=not None
525    
526     def setDebugOff(self):     def setDebugOff(self):
527       """       """
528       switches off debugging       Switches debug output off.
529       """       """
530       self.__debug=None       self.__debug=None
531    
532     def setDebug(self, flag):     def setDebug(self, flag):
533       """       """
534       switches debugging to on if C{flag} is true outherwise debugging is set to off       Switches debug output on if C{flag} is True otherwise it is switched off.
535    
536       @param flag: desired debug status       @param flag: desired debug status
537       @type flag: C{bool}       @type flag: C{bool}
538       """       """
# Line 492  class LinearProblem(object): Line 543  class LinearProblem(object):
543    
544     def trace(self,text):     def trace(self,text):
545       """       """
546       print the text message if debugging is swiched on.       Prints the text message if debug mode is switched on.
547       @param text: message  
548         @param text: message to be printed
549       @type text: C{string}       @type text: C{string}
550       """       """
551       if self.__debug: print "%s: %s"%(str(self),text)       if self.__debug: print "%s: %s"%(str(self),text)
552    
553     # =============================================================================     # ==========================================================================
554     # some service functions:     # some service functions:
555     # =============================================================================     # ==========================================================================
556     def introduceCoefficients(self,**coeff):     def introduceCoefficients(self,**coeff):
557         """         """
558         introduces a new coefficient into the problem.         Introduces new coefficients into the problem.
559    
560         use         Use::
561    
562         self.introduceCoefficients(A=PDECoef(...),B=PDECoef(...))         p.introduceCoefficients(A=PDECoef(...), B=PDECoef(...))
563    
564         to introduce the coefficients M{A} ans M{B}.         to introduce the coefficients M{A} ans M{B}.
565         """         """
# Line 520  class LinearProblem(object): Line 572  class LinearProblem(object):
572    
573     def getDomain(self):     def getDomain(self):
574       """       """
575       returns the domain of the PDE       Returns the domain of the PDE.
576    
577       @return: the domain of the PDE       @return: the domain of the PDE
578       @rtype: L{Domain<escript.Domain>}       @rtype: L{Domain<escript.Domain>}
# Line 529  class LinearProblem(object): Line 581  class LinearProblem(object):
581    
582     def getDim(self):     def getDim(self):
583       """       """
584       returns the spatial dimension of the PDE       Returns the spatial dimension of the PDE.
585    
586       @return: the spatial dimension of the PDE domain       @return: the spatial dimension of the PDE domain
587       @rtype: C{int}       @rtype: C{int}
# Line 538  class LinearProblem(object): Line 590  class LinearProblem(object):
590    
591     def getNumEquations(self):     def getNumEquations(self):
592       """       """
593       returns the number of equations       Returns the number of equations.
594    
595       @return: the number of equations       @return: the number of equations
596       @rtype: C{int}       @rtype: C{int}
597       @raise UndefinedPDEError: if the number of equations is not be specified yet.       @raise UndefinedPDEError: if the number of equations is not specified yet
598       """       """
599       if self.__numEquations==None:       if self.__numEquations==None:
600           if self.__numSolutions==None:           if self.__numSolutions==None:
# Line 553  class LinearProblem(object): Line 605  class LinearProblem(object):
605    
606     def getNumSolutions(self):     def getNumSolutions(self):
607       """       """
608       returns the number of unknowns       Returns the number of unknowns.
609    
610       @return: the number of unknowns       @return: the number of unknowns
611       @rtype: C{int}       @rtype: C{int}
612       @raise UndefinedPDEError: if the number of unknowns is not be specified yet.       @raise UndefinedPDEError: if the number of unknowns is not specified yet
613       """       """
614       if self.__numSolutions==None:       if self.__numSolutions==None:
615          if self.__numEquations==None:          if self.__numEquations==None:
# Line 568  class LinearProblem(object): Line 620  class LinearProblem(object):
620    
621     def reduceEquationOrder(self):     def reduceEquationOrder(self):
622       """       """
623       return status for order reduction for equation       Returns the status of order reduction for the equation.
624    
625       @return: return True is reduced interpolation order is used for the represenation of the equation       @return: True if reduced interpolation order is used for the
626                  representation of the equation, False otherwise
627       @rtype: L{bool}       @rtype: L{bool}
628       """       """
629       return self.__reduce_equation_order       return self.__reduce_equation_order
630    
631     def reduceSolutionOrder(self):     def reduceSolutionOrder(self):
632       """       """
633       return status for order reduction for the solution       Returns the status of order reduction for the solution.
634    
635       @return: return True is reduced interpolation order is used for the represenation of the solution       @return: True if reduced interpolation order is used for the
636                  representation of the solution, False otherwise
637       @rtype: L{bool}       @rtype: L{bool}
638       """       """
639       return self.__reduce_solution_order       return self.__reduce_solution_order
640    
641     def getFunctionSpaceForEquation(self):     def getFunctionSpaceForEquation(self):
642       """       """
643       returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation       Returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize
644         the equation.
645    
646       @return: representation space of equation       @return: representation space of equation
647       @rtype: L{FunctionSpace<escript.FunctionSpace>}       @rtype: L{FunctionSpace<escript.FunctionSpace>}
# Line 598  class LinearProblem(object): Line 653  class LinearProblem(object):
653    
654     def getFunctionSpaceForSolution(self):     def getFunctionSpaceForSolution(self):
655       """       """
656       returns the L{FunctionSpace<escript.FunctionSpace>} used to represent the solution       Returns the L{FunctionSpace<escript.FunctionSpace>} used to represent
657         the solution.
658    
659       @return: representation space of solution       @return: representation space of solution
660       @rtype: L{FunctionSpace<escript.FunctionSpace>}       @rtype: L{FunctionSpace<escript.FunctionSpace>}
# Line 608  class LinearProblem(object): Line 664  class LinearProblem(object):
664       else:       else:
665           return escript.Solution(self.getDomain())           return escript.Solution(self.getDomain())
666    
667     # =============================================================================     # ==========================================================================
668     #   solver settings:     #   solver settings:
669     # =============================================================================     # ==========================================================================
670     def setSolverMethod(self,solver=None,preconditioner=None):     def setSolverMethod(self,solver=None,preconditioner=None):
671         """         """
672         sets a new solver         Sets a new solver method and/or preconditioner.
673    
674         @param solver: sets a new solver method.         @param solver: new solver method to use
675         @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY}, L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES}, L{TFQMR}, L{MINRES}, L{PRES20}, L{LUMPING}, L{AMG}         @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY},
676         @param preconditioner: sets a new solver method.                       L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES},
677         @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU},  L{GS}                       L{TFQMR}, L{MINRES}, L{PRES20}, L{LUMPING}, L{AMG}
678           @param preconditioner: new preconditioner to use
679           @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},
680                                 L{SSOR}, L{RILU}, L{GS}
681    
682         @note: the solver method choosen may not be available by the selected package.         @note: the solver method chosen may not be available by the selected
683                  package.
684         """         """
685         if solver==None: solver=self.__solver_method         if solver==None: solver=self.__solver_method
686         if preconditioner==None: preconditioner=self.__preconditioner         if preconditioner==None: preconditioner=self.__preconditioner
# Line 632  class LinearProblem(object): Line 692  class LinearProblem(object):
692             self.updateSystemType()             self.updateSystemType()
693             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
694    
695     def getSolverMethodName(self):     def getSolverMethodName(self):
696         """         """
697         returns the name of the solver currently used         Returns the name of the solver currently used.
698    
699         @return: the name of the solver currently used.         @return: the name of the solver currently used
700         @rtype: C{string}         @rtype: C{string}
701         """         """
702    
# Line 675  class LinearProblem(object): Line 735  class LinearProblem(object):
735         else : method="unknown"         else : method="unknown"
736         return "%s solver of %s package"%(method,package)         return "%s solver of %s package"%(method,package)
737    
738     def getSolverMethod(self):     def getSolverMethod(self):
739         """         """
740         returns the solver method and preconditioner used         Returns the solver method and preconditioner used.
741    
742         @return: the solver method currently be used.         @return: the solver method currently used.
743         @rtype: C{tuple} of C{int}         @rtype: C{tuple} of C{int}
744         """         """
745         return self.__solver_method,self.__preconditioner         return self.__solver_method,self.__preconditioner
746    
747     def setSolverPackage(self,package=None):     def setSolverPackage(self,package=None):
748         """         """
749         sets a new solver package         Sets a new solver package.
750    
751         @param package: sets a new solver method.         @param package: the new solver package
752         @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}, L{TRILINOS}         @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK},
753                          L{TRILINOS}
754         """         """
755         if package==None: package=self.DEFAULT         if package==None: package=self.DEFAULT
756         if not package==self.getSolverPackage():         if not package==self.getSolverPackage():
# Line 697  class LinearProblem(object): Line 758  class LinearProblem(object):
758             self.updateSystemType()             self.updateSystemType()
759             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
760    
761     def getSolverPackage(self):     def getSolverPackage(self):
762         """         """
763         returns the package of the solver         Returns the package of the solver.
764    
765         @return: the solver package currently being used.         @return: the solver package currently being used
766         @rtype: C{int}         @rtype: C{int}
767         """         """
768         return self.__solver_package         return self.__solver_package
769    
770     def isUsingLumping(self):     def isUsingLumping(self):
771        """        """
772        checks if matrix lumping is used as a solver method        Checks if matrix lumping is the current solver method.
773    
774        @return: True is lumping is currently used a solver method.        @return: True if the current solver method is lumping
775        @rtype: C{bool}        @rtype: C{bool}
776        """        """
777        return self.getSolverMethod()[0]==self.LUMPING        return self.getSolverMethod()[0]==self.LUMPING
778    
779     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
780         """         """
781         resets the tolerance for the solver method to tol.         Resets the tolerance for the solver method to C{tol}.
782    
783         @param tol: new tolerance for the solver. If the tol is lower then the current tolerence         @param tol: new tolerance for the solver. If C{tol} is lower than the
784                     the system will be resolved.                     current tolerence the system will be resolved.
785         @type tol: positive C{float}         @type tol: positive C{float}
786         @raise ValueError: if tolerance is not positive.         @raise ValueError: if tolerance is not positive
787         """         """
788         if not tol>0:         if not tol>0:
789             raise ValueError,"Tolerance as to be positive"             raise ValueError,"Tolerance has to be positive"
790         if tol<self.getTolerance(): self.invalidateSolution()         if tol<self.getTolerance(): self.invalidateSolution()
791         self.trace("New tolerance %e"%tol)         self.trace("New tolerance %e"%tol)
792         self.__tolerance=tol         self.__tolerance=tol
793         return         return
794    
795     def getTolerance(self):     def getTolerance(self):
796         """         """
797         returns the tolerance set for the solution         Returns the tolerance currently set for the solution.
798    
799         @return: tolerance currently used.         @return: tolerance currently used
800         @rtype: C{float}         @rtype: C{float}
801         """         """
802         return self.__tolerance         return self.__tolerance
803    
804     # =============================================================================     # ==========================================================================
805     #    symmetry  flag:     #    symmetry  flag:
806     # =============================================================================     # ==========================================================================
807     def isSymmetric(self):     def isSymmetric(self):
808        """        """
809        checks if symmetry is indicated.        Checks if symmetry is indicated.
810    
811        @return: True is a symmetric PDE is indicated, otherwise False is returned        @return: True if a symmetric PDE is indicated, False otherwise
812        @rtype: C{bool}        @rtype: C{bool}
813        """        """
814        return self.__sym        return self.__sym
815    
816     def setSymmetryOn(self):     def setSymmetryOn(self):
817        """        """
818        sets the symmetry flag.        Sets the symmetry flag.
819        """        """
820        if not self.isSymmetric():        if not self.isSymmetric():
821           self.trace("PDE is set to be symmetric")           self.trace("PDE is set to be symmetric")
822           self.__sym=True           self.__sym=True
823           self.updateSystemType()           self.updateSystemType()
824    
825     def setSymmetryOff(self):     def setSymmetryOff(self):
826        """        """
827        removes the symmetry flag.        Clears the symmetry flag.
828        """        """
829        if self.isSymmetric():        if self.isSymmetric():
830           self.trace("PDE is set to be unsymmetric")           self.trace("PDE is set to be nonsymmetric")
831           self.__sym=False           self.__sym=False
832           self.updateSystemType()           self.updateSystemType()
833    
834     def setSymmetryTo(self,flag=False):     def setSymmetryTo(self,flag=False):
835        """        """
836        sets the symmetry flag to flag        Sets the symmetry flag to C{flag}.
837    
838        @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.        @param flag: If True, the symmetry flag is set otherwise reset.
839        @type flag: C{bool}        @type flag: C{bool}
840        """        """
841        if flag:        if flag:
# Line 782  class LinearProblem(object): Line 843  class LinearProblem(object):
843        else:        else:
844           self.setSymmetryOff()           self.setSymmetryOff()
845    
846     # =============================================================================     # ==========================================================================
847     # function space handling for the equation as well as the solution     # function space handling for the equation as well as the solution
848     # =============================================================================     # ==========================================================================
849     def setReducedOrderOn(self):     def setReducedOrderOn(self):
850       """       """
851       switches on reduced order for solution and equation representation       Switches reduced order on for solution and equation representation.
852    
853       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @raise RuntimeError: if order reduction is altered after a coefficient has
854                              been set
855       """       """
856       self.setReducedOrderForSolutionOn()       self.setReducedOrderForSolutionOn()
857       self.setReducedOrderForEquationOn()       self.setReducedOrderForEquationOn()
858    
859     def setReducedOrderOff(self):     def setReducedOrderOff(self):
860       """       """
861       switches off reduced order for solution and equation representation       Switches reduced order off for solution and equation representation
862    
863       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @raise RuntimeError: if order reduction is altered after a coefficient has
864                              been set
865       """       """
866       self.setReducedOrderForSolutionOff()       self.setReducedOrderForSolutionOff()
867       self.setReducedOrderForEquationOff()       self.setReducedOrderForEquationOff()
868    
869     def setReducedOrderTo(self,flag=False):     def setReducedOrderTo(self,flag=False):
870       """       """
871       sets order reduction for both solution and equation representation according to flag.       Sets order reduction state for both solution and equation representation
872       @param flag: if flag is True, the order reduction is switched on for both  solution and equation representation, otherwise or       according to flag.
873                    if flag is not present order reduction is switched off  
874         @param flag: if True, the order reduction is switched on for both solution
875                      and equation representation, otherwise or if flag is not
876                      present order reduction is switched off
877       @type flag: C{bool}       @type flag: C{bool}
878       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @raise RuntimeError: if order reduction is altered after a coefficient has
879                              been set
880       """       """
881       self.setReducedOrderForSolutionTo(flag)       self.setReducedOrderForSolutionTo(flag)
882       self.setReducedOrderForEquationTo(flag)       self.setReducedOrderForEquationTo(flag)
883    
884    
885     def setReducedOrderForSolutionOn(self):     def setReducedOrderForSolutionOn(self):
886       """       """
887       switches on reduced order for solution representation       Switches reduced order on for solution representation.
888    
889       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @raise RuntimeError: if order reduction is altered after a coefficient has
890                              been set
891       """       """
892       if not self.__reduce_solution_order:       if not self.__reduce_solution_order:
893           if self.__altered_coefficients:           if self.__altered_coefficients:
894                raise RuntimeError,"order cannot be altered after coefficients have been defined."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
895           self.trace("Reduced order is used to solution representation.")           self.trace("Reduced order is used for solution representation.")
896           self.__reduce_solution_order=True           self.__reduce_solution_order=True
897           self.initializeSystem()           self.initializeSystem()
898    
899     def setReducedOrderForSolutionOff(self):     def setReducedOrderForSolutionOff(self):
900       """       """
901       switches off reduced order for solution representation       Switches reduced order off for solution representation
902    
903       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @raise RuntimeError: if order reduction is altered after a coefficient has
904                              been set.
905       """       """
906       if self.__reduce_solution_order:       if self.__reduce_solution_order:
907           if self.__altered_coefficients:           if self.__altered_coefficients:
# Line 841  class LinearProblem(object): Line 910  class LinearProblem(object):
910           self.__reduce_solution_order=False           self.__reduce_solution_order=False
911           self.initializeSystem()           self.initializeSystem()
912    
913     def setReducedOrderForSolutionTo(self,flag=False):     def setReducedOrderForSolutionTo(self,flag=False):
914       """       """
915       sets order for test functions according to flag       Sets order reduction state for solution representation according to flag.
916    
917       @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or       @param flag: if flag is True, the order reduction is switched on for
918                    if flag is not present order reduction is switched off                    solution representation, otherwise or if flag is not present
919                      order reduction is switched off
920       @type flag: C{bool}       @type flag: C{bool}
921       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @raise RuntimeError: if order reduction is altered after a coefficient has
922                              been set
923       """       """
924       if flag:       if flag:
925          self.setReducedOrderForSolutionOn()          self.setReducedOrderForSolutionOn()
926       else:       else:
927          self.setReducedOrderForSolutionOff()          self.setReducedOrderForSolutionOff()
928    
929     def setReducedOrderForEquationOn(self):     def setReducedOrderForEquationOn(self):
930       """       """
931       switches on reduced order for equation representation       Switches reduced order on for equation representation.
932    
933       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @raise RuntimeError: if order reduction is altered after a coefficient has
934                              been set
935       """       """
936       if not self.__reduce_equation_order:       if not self.__reduce_equation_order:
937           if self.__altered_coefficients:           if self.__altered_coefficients:
# Line 868  class LinearProblem(object): Line 940  class LinearProblem(object):
940           self.__reduce_equation_order=True           self.__reduce_equation_order=True
941           self.initializeSystem()           self.initializeSystem()
942    
943     def setReducedOrderForEquationOff(self):     def setReducedOrderForEquationOff(self):
944       """       """
945       switches off reduced order for equation representation       Switches reduced order off for equation representation.
946    
947       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @raise RuntimeError: if order reduction is altered after a coefficient has
948                              been set
949       """       """
950       if self.__reduce_equation_order:       if self.__reduce_equation_order:
951           if self.__altered_coefficients:           if self.__altered_coefficients:
# Line 881  class LinearProblem(object): Line 954  class LinearProblem(object):
954           self.__reduce_equation_order=False           self.__reduce_equation_order=False
955           self.initializeSystem()           self.initializeSystem()
956    
957     def setReducedOrderForEquationTo(self,flag=False):     def setReducedOrderForEquationTo(self,flag=False):
958       """       """
959       sets order for test functions according to flag       Sets order reduction state for equation representation according to flag.
960    
961       @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or       @param flag: if flag is True, the order reduction is switched on for
962                    if flag is not present order reduction is switched off                    equation representation, otherwise or if flag is not present
963                      order reduction is switched off
964       @type flag: C{bool}       @type flag: C{bool}
965       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @raise RuntimeError: if order reduction is altered after a coefficient has
966                              been set
967       """       """
968       if flag:       if flag:
969          self.setReducedOrderForEquationOn()          self.setReducedOrderForEquationOn()
970       else:       else:
971          self.setReducedOrderForEquationOff()          self.setReducedOrderForEquationOff()
972    
973     def updateSystemType(self):     def updateSystemType(self):
974       """       """
975       reasses the system type and, if a new matrix is needed, resets the system.       Reassesses the system type and, if a new matrix is needed, resets the
976         system.
977       """       """
978       new_system_type=self.getRequiredSystemType()       new_system_type=self.getRequiredSystemType()
979       if not new_system_type==self.__system_type:       if not new_system_type==self.__system_type:
980           self.trace("Matrix type is now %d."%new_system_type)           self.trace("Matrix type is now %d."%new_system_type)
981           self.__system_type=new_system_type           self.__system_type=new_system_type
982           self.initializeSystem()           self.initializeSystem()
983    
984     def getSystemType(self):     def getSystemType(self):
985        """        """
986        return the current system type        Returns the current system type.
987        """        """
988        return self.__system_type        return self.__system_type
989    
990     def checkSymmetricTensor(self,name,verbose=True):     def checkSymmetricTensor(self,name,verbose=True):
991        """        """
992        tests a coefficient for symmetry        Tests a coefficient for symmetry.
993    
994        @param name: name of the coefficient        @param name: name of the coefficient
995        @type name: C{str}        @type name: C{str}
996        @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.        @param verbose: if set to True or not present a report on coefficients
997                          which break the symmetry is printed.
998        @type verbose: C{bool}        @type verbose: C{bool}
999        @return:  True if coefficient C{name} is symmetric.        @return: True if coefficient C{name} is symmetric
1000        @rtype: C{bool}        @rtype: C{bool}
1001        """        """
1002        A=self.getCoefficient(name)        A=self.getCoefficient(name)
# Line 957  class LinearProblem(object): Line 1035  class LinearProblem(object):
1035    
1036     def checkReciprocalSymmetry(self,name0,name1,verbose=True):     def checkReciprocalSymmetry(self,name0,name1,verbose=True):
1037        """        """
1038        test two coefficients for resciprocal symmetry        Tests two coefficients for reciprocal symmetry.
1039    
1040        @param name0: name of the first coefficient        @param name0: name of the first coefficient
1041        @type name0: C{str}        @type name0: C{str}
1042        @param name1: name of the second coefficient        @param name1: name of the second coefficient
1043        @type name1: C{str}        @type name1: C{str}
1044        @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.        @param verbose: if set to True or not present a report on coefficients
1045                          which break the symmetry is printed
1046        @type verbose: C{bool}        @type verbose: C{bool}
1047        @return:  True if coefficients C{name0} and C{name1}  are resciprocally symmetric.        @return: True if coefficients C{name0} and C{name1} are reciprocally
1048                   symmetric.
1049        @rtype: C{bool}        @rtype: C{bool}
1050        """        """
1051        B=self.getCoefficient(name0)        B=self.getCoefficient(name0)
# Line 1012  class LinearProblem(object): Line 1092  class LinearProblem(object):
1092                   raise ValueError,"Cannot check rank %s of %s and %s."%(len(sB),name0,name1)                   raise ValueError,"Cannot check rank %s of %s and %s."%(len(sB),name0,name1)
1093        return out        return out
1094    
1095     def getCoefficient(self,name):     def getCoefficient(self,name):
1096       """       """
1097       returns the value of the coefficient name       Returns the value of the coefficient C{name}.
1098    
1099       @param name: name of the coefficient requested.       @param name: name of the coefficient requested
1100       @type name: C{string}       @type name: C{string}
1101       @return: the value of the coefficient name       @return: the value of the coefficient
1102       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
1103       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1104       """       """
1105       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1106           return self.__COEFFICIENTS[name].getValue()           return self.__COEFFICIENTS[name].getValue()
1107       else:       else:
1108          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1109    
1110     def hasCoefficient(self,name):     def hasCoefficient(self,name):
1111       """       """
1112       return True if name is the name of a coefficient       Returns True if C{name} is the name of a coefficient.
1113    
1114       @param name: name of the coefficient enquired.       @param name: name of the coefficient enquired
1115       @type name: C{string}       @type name: C{string}
1116       @return: True if name is the name of a coefficient of the general PDE. Otherwise False.       @return: True if C{name} is the name of a coefficient of the general PDE,
1117                  False otherwise
1118       @rtype: C{bool}       @rtype: C{bool}
1119       """       """
1120       return self.__COEFFICIENTS.has_key(name)       return self.__COEFFICIENTS.has_key(name)
1121    
1122     def createCoefficient(self, name):     def createCoefficient(self, name):
1123       """       """
1124       create a L{Data<escript.Data>} object corresponding to coefficient name       Creates a L{Data<escript.Data>} object corresponding to coefficient
1125         C{name}.
1126    
1127       @return: a coefficient name initialized to 0.       @return: the coefficient C{name} initialized to 0
1128       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
1129       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1130       """       """
1131       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1132          return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))          return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1133       else:       else:
1134          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1135    
1136     def getFunctionSpaceForCoefficient(self,name):     def getFunctionSpaceForCoefficient(self,name):
1137       """       """
1138       return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name       Returns the L{FunctionSpace<escript.FunctionSpace>} to be used for
1139         coefficient C{name}.
1140    
1141       @param name: name of the coefficient enquired.       @param name: name of the coefficient enquired
1142       @type name: C{string}       @type name: C{string}
1143       @return: the function space to be used for coefficient name       @return: the function space to be used for coefficient C{name}
1144       @rtype: L{FunctionSpace<escript.FunctionSpace>}       @rtype: L{FunctionSpace<escript.FunctionSpace>}
1145       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1146       """       """
1147       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1148          return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())          return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1149       else:       else:
1150          raise ValueError,"unknown coefficient %s requested"%name          raise ValueError,"unknown coefficient %s requested"%name
1151    
1152     def getShapeOfCoefficient(self,name):     def getShapeOfCoefficient(self,name):
1153       """       """
1154       return the shape of the coefficient name       Returns the shape of the coefficient C{name}.
1155    
1156       @param name: name of the coefficient enquired.       @param name: name of the coefficient enquired
1157       @type name: C{string}       @type name: C{string}
1158       @return: the shape of the coefficient name       @return: the shape of the coefficient C{name}
1159       @rtype: C{tuple} of C{int}       @rtype: C{tuple} of C{int}
1160       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1161       """       """
1162       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1163          return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())          return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1164       else:       else:
1165          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1166    
1167     def resetAllCoefficients(self):     def resetAllCoefficients(self):
1168       """       """
1169       resets all coefficients to there default values.       Resets all coefficients to their default values.
1170       """       """
1171       for i in self.__COEFFICIENTS.iterkeys():       for i in self.__COEFFICIENTS.iterkeys():
1172           self.__COEFFICIENTS[i].resetValue()           self.__COEFFICIENTS[i].resetValue()
1173    
1174     def alteredCoefficient(self,name):     def alteredCoefficient(self,name):
1175       """       """
1176       announce that coefficient name has been changed       Announces that coefficient C{name} has been changed.
1177    
1178       @param name: name of the coefficient enquired.       @param name: name of the coefficient affected
1179       @type name: C{string}       @type name: C{string}
1180       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1181       @note: if name is q or r, the method will not trigger a rebuilt of the system as constraints are applied to the solved system.       @note: if C{name} is q or r, the method will not trigger a rebuild of the
1182                system as constraints are applied to the solved system.
1183       """       """
1184       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1185          self.trace("Coefficient %s has been altered."%name)          self.trace("Coefficient %s has been altered."%name)
# Line 1105  class LinearProblem(object): Line 1189  class LinearProblem(object):
1189       else:       else:
1190          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1191    
1192     def validSolution(self):     def validSolution(self):
1193         """         """
1194         markes the solution as valid         Marks the solution as valid.
1195         """         """
1196         self.__is_solution_valid=True         self.__is_solution_valid=True
1197    
1198     def invalidateSolution(self):     def invalidateSolution(self):
1199         """         """
1200         indicates the PDE has to be resolved if the solution is requested         Indicates the PDE has to be resolved if the solution is requested.
1201         """         """
1202         self.trace("System will be resolved.")         self.trace("System will be resolved.")
1203         self.__is_solution_valid=False         self.__is_solution_valid=False
1204    
1205     def isSolutionValid(self):     def isSolutionValid(self):
1206         """         """
1207         returns True is the solution is still valid.         Returns True if the solution is still valid.
1208         """         """
1209         return self.__is_solution_valid         return self.__is_solution_valid
1210    
1211     def validOperator(self):     def validOperator(self):
1212         """         """
1213         markes the  operator as valid         Marks the operator as valid.
1214         """         """
1215         self.__is_operator_valid=True         self.__is_operator_valid=True
1216    
1217     def invalidateOperator(self):     def invalidateOperator(self):
1218         """         """
1219         indicates the operator has to be rebuilt next time it is used         Indicates the operator has to be rebuilt next time it is used.
1220         """         """
1221         self.trace("Operator will be rebuilt.")         self.trace("Operator will be rebuilt.")
1222         self.invalidateSolution()         self.invalidateSolution()
# Line 1140  class LinearProblem(object): Line 1224  class LinearProblem(object):
1224    
1225     def isOperatorValid(self):     def isOperatorValid(self):
1226         """         """
1227         returns True is the operator is still valid.         Returns True if the operator is still valid.
1228         """         """
1229         return self.__is_operator_valid         return self.__is_operator_valid
1230    
1231     def validRightHandSide(self):     def validRightHandSide(self):
1232         """         """
1233         markes the right hand side as valid         Marks the right hand side as valid.
1234         """         """
1235         self.__is_RHS_valid=True         self.__is_RHS_valid=True
1236    
1237     def invalidateRightHandSide(self):     def invalidateRightHandSide(self):
1238         """         """
1239         indicates the right hand side has to be rebuild next time it is used         Indicates the right hand side has to be rebuilt next time it is used.
1240         """         """
1241         if self.isRightHandSideValid(): self.trace("Right hand side has to be rebuilt.")         if self.isRightHandSideValid(): self.trace("Right hand side has to be rebuilt.")
1242         self.invalidateSolution()         self.invalidateSolution()
# Line 1160  class LinearProblem(object): Line 1244  class LinearProblem(object):
1244    
1245     def isRightHandSideValid(self):     def isRightHandSideValid(self):
1246         """         """
1247         returns True is the operator is still valid.         Returns True if the operator is still valid.
1248         """         """
1249         return self.__is_RHS_valid         return self.__is_RHS_valid
1250    
1251     def invalidateSystem(self):     def invalidateSystem(self):
1252         """         """
1253         annonced that everthing has to be rebuild:         Announces that everything has to be rebuilt.
1254         """         """
1255         if self.isRightHandSideValid(): self.trace("System has to be rebuilt.")         if self.isRightHandSideValid(): self.trace("System has to be rebuilt.")
1256         self.invalidateSolution()         self.invalidateSolution()
# Line 1175  class LinearProblem(object): Line 1259  class LinearProblem(object):
1259    
1260     def isSystemValid(self):     def isSystemValid(self):
1261         """         """
1262         returns true if the system (including solution) is still vaild:         Returns True if the system (including solution) is still vaild.
1263         """         """
1264         return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()         return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()
1265    
1266     def initializeSystem(self):     def initializeSystem(self):
1267         """         """
1268         resets the system         Resets the system clearing the operator, right hand side and solution.
1269         """         """
1270         self.trace("New System has been created.")         self.trace("New System has been created.")
1271         self.__operator=escript.Operator()         self.__operator=escript.Operator()
# Line 1189  class LinearProblem(object): Line 1273  class LinearProblem(object):
1273         self.__solution=escript.Data()         self.__solution=escript.Data()
1274         self.invalidateSystem()         self.invalidateSystem()
1275    
1276     def getOperator(self):     def getOperator(self):
1277       """       """
1278       returns the operator of the linear problem       Returns the operator of the linear problem.
1279    
1280       @return: the operator of the problem       @return: the operator of the problem
1281       """       """
1282       return self.getSystem()[0]       return self.getSystem()[0]
1283    
1284     def getRightHandSide(self):     def getRightHandSide(self):
1285       """       """
1286       returns the right hand side of the linear problem       Returns the right hand side of the linear problem.
1287    
1288       @return: the right hand side of the problem       @return: the right hand side of the problem
1289       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
1290       """       """
1291       return self.getSystem()[1]       return self.getSystem()[1]
1292    
1293     def createRightHandSide(self):     def createRightHandSide(self):
1294         """         """
1295         returns an instance of a new right hand side         Returns an instance of a new right hand side.
1296         """         """
1297         self.trace("New right hand side is allocated.")         self.trace("New right hand side is allocated.")
1298         if self.getNumEquations()>1:         if self.getNumEquations()>1:
# Line 1216  class LinearProblem(object): Line 1300  class LinearProblem(object):
1300         else:         else:
1301             return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)             return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1302    
1303     def createSolution(self):     def createSolution(self):
1304         """         """
1305         returns an instance of a new solution         Returns an instance of a new solution.
1306         """         """
1307         self.trace("New solution is allocated.")         self.trace("New solution is allocated.")
1308         if self.getNumSolutions()>1:         if self.getNumSolutions()>1:
# Line 1226  class LinearProblem(object): Line 1310  class LinearProblem(object):
1310         else:         else:
1311             return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)             return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1312    
1313     def resetSolution(self):     def resetSolution(self):
1314         """         """
1315         sets the solution to zero         Sets the solution to zero.
1316         """         """
1317         if self.__solution.isEmpty():         if self.__solution.isEmpty():
1318             self.__solution=self.createSolution()             self.__solution=self.createSolution()
1319         else:         else:
1320             self.__solution.setToZero()             self.__solution.setToZero()
1321             self.trace("Solution is reset to zero.")             self.trace("Solution is reset to zero.")
1322    
1323     def setSolution(self,u):     def setSolution(self,u):
1324         """         """
1325         sets the solution assuming that makes the sytem valid.         Sets the solution assuming that makes the system valid.
1326         """         """
1327         self.__solution=u         self.__solution=u
1328         self.validSolution()         self.validSolution()
1329    
1330     def getCurrentSolution(self):     def getCurrentSolution(self):
1331         """         """
1332         returns the solution in its current state         Returns the solution in its current state.
1333         """         """
1334         if self.__solution.isEmpty(): self.__solution=self.createSolution()         if self.__solution.isEmpty(): self.__solution=self.createSolution()
1335         return self.__solution         return self.__solution
1336    
1337     def resetRightHandSide(self):     def resetRightHandSide(self):
1338         """         """
1339         sets the right hand side to zero         Sets the right hand side to zero.
1340         """         """
1341         if self.__righthandside.isEmpty():         if self.__righthandside.isEmpty():
1342             self.__righthandside=self.createRightHandSide()             self.__righthandside=self.createRightHandSide()
# Line 1261  class LinearProblem(object): Line 1346  class LinearProblem(object):
1346    
1347     def getCurrentRightHandSide(self):     def getCurrentRightHandSide(self):
1348         """         """
1349         returns the right hand side  in its current state         Returns the right hand side in its current state.
1350         """         """
1351         if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()         if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
1352         return self.__righthandside         return self.__righthandside
           
1353    
1354     def resetOperator(self):     def resetOperator(self):
1355         """         """
1356         makes sure that the operator is instantiated and returns it initialized by zeros         Makes sure that the operator is instantiated and returns it initialized
1357           with zeros.
1358         """         """
1359         if self.__operator.isEmpty():         if self.__operator.isEmpty():
1360             if self.isUsingLumping():             if self.isUsingLumping():
# Line 1285  class LinearProblem(object): Line 1370  class LinearProblem(object):
1370    
1371     def getCurrentOperator(self):     def getCurrentOperator(self):
1372         """         """
1373         returns the operator in its current state         Returns the operator in its current state.
1374         """         """
1375         if self.__operator.isEmpty(): self.resetOperator()         if self.__operator.isEmpty(): self.resetOperator()
1376         return self.__operator         return self.__operator
1377    
1378     def setValue(self,**coefficients):     def setValue(self,**coefficients):
1379        """        """
1380        sets new values to coefficients        Sets new values to coefficients.
1381    
1382        @raise IllegalCoefficient: if an unknown coefficient keyword is used.        @raise IllegalCoefficient: if an unknown coefficient keyword is used
1383        """        """
1384        # check if the coefficients are  legal:        # check if the coefficients are  legal:
1385        for i in coefficients.iterkeys():        for i in coefficients.iterkeys():
# Line 1317  class LinearProblem(object): Line 1402  class LinearProblem(object):
1402                  else:                  else:
1403                      if self.__numEquations==None: self.__numEquations=res[0]                      if self.__numEquations==None: self.__numEquations=res[0]
1404                      if self.__numSolutions==None: self.__numSolutions=res[1]                      if self.__numSolutions==None: self.__numSolutions=res[1]
1405        if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"        if self.__numEquations==None: raise UndefinedPDEError,"unidentified number of equations"
1406        if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions"        if self.__numSolutions==None: raise UndefinedPDEError,"unidentified number of solutions"
1407        # now we check the shape of the coefficient if numEquations and numSolutions are set:        # now we check the shape of the coefficient if numEquations and numSolutions are set:
1408        for i,d in coefficients.iteritems():        for i,d in coefficients.iteritems():
1409          try:          try:
1410             self.__COEFFICIENTS[i].setValue(self.getDomain(),             self.__COEFFICIENTS[i].setValue(self.getDomain(),
1411                                           self.getNumEquations(),self.getNumSolutions(),                       self.getNumEquations(),self.getNumSolutions(),
1412                                           self.reduceEquationOrder(),self.reduceSolutionOrder(),d)                       self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1413             self.alteredCoefficient(i)             self.alteredCoefficient(i)
1414          except IllegalCoefficientFunctionSpace,m:          except IllegalCoefficientFunctionSpace,m:
1415              # if the function space is wrong then we try the reduced version:              # if the function space is wrong then we try the reduced version:
# Line 1345  class LinearProblem(object): Line 1430  class LinearProblem(object):
1430             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1431        self.__altered_coefficients=True        self.__altered_coefficients=True
1432    
1433     # ====================================================================================     # ==========================================================================
1434     # methods that are typically overwritten when implementing a particular linear problem     # methods that are typically overwritten when implementing a particular
1435     # ====================================================================================     # linear problem
1436       # ==========================================================================
1437     def getRequiredSystemType(self):     def getRequiredSystemType(self):
1438        """        """
1439        returns the system type which needs to be used by the current set up.        Returns the system type which needs to be used by the current set up.
1440    
1441        @note: Typically this method is overwritten when implementing a particular linear problem.        @note: Typically this method is overwritten when implementing a
1442                 particular linear problem.
1443        """        """
1444        return 0        return 0
1445    
1446     def createOperator(self):     def createOperator(self):
1447         """         """
1448         returns an instance of a new operator         Returns an instance of a new operator.
1449    
1450         @note: This method is overwritten when implementing a particular linear problem.         @note: This method is overwritten when implementing a particular
1451                  linear problem.
1452         """         """
1453         return escript.Operator()         return escript.Operator()
1454    
1455     def checkSymmetry(self,verbose=True):     def checkSymmetry(self,verbose=True):
1456        """        """
1457        test the PDE for symmetry.        Tests the PDE for symmetry.
1458    
1459        @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.        @param verbose: if set to True or not present a report on coefficients
1460                          which break the symmetry is printed
1461        @type verbose: C{bool}        @type verbose: C{bool}
1462        @return:  True if the problem is symmetric.        @return: True if the problem is symmetric
1463        @rtype: C{bool}        @rtype: C{bool}
1464        @note: Typically this method is overwritten when implementing a particular linear problem.        @note: Typically this method is overwritten when implementing a
1465                 particular linear problem.
1466        """        """
1467        out=True        out=True
1468        return out        return out
1469    
1470     def getSolution(self,**options):     def getSolution(self,**options):
1471         """         """
1472         returns the solution of the problem         Returns the solution of the problem.
1473    
1474         @return: the solution         @return: the solution
1475         @rtype: L{Data<escript.Data>}         @rtype: L{Data<escript.Data>}
1476    
1477         @note: This method is overwritten when implementing a particular linear problem.         @note: This method is overwritten when implementing a particular
1478                  linear problem.
1479         """         """
1480         return self.getCurrentSolution()         return self.getCurrentSolution()
1481    
1482     def getSystem(self):     def getSystem(self):
1483         """         """
1484         return the operator and right hand side of the PDE         Returns the operator and right hand side of the PDE.
1485    
1486         @return: the discrete version of the PDE         @return: the discrete version of the PDE
1487         @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.         @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1488    
1489         @note: This method is overwritten when implementing a particular linear problem.         @note: This method is overwritten when implementing a particular
1490                  linear problem.
1491         """         """
1492         return (self.getCurrentOperator(), self.getCurrentRightHandSide())         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1493  #=====================  #=====================
# Line 1402  class LinearProblem(object): Line 1495  class LinearProblem(object):
1495  class LinearPDE(LinearProblem):  class LinearPDE(LinearProblem):
1496     """     """
1497     This class is used to define a general linear, steady, second order PDE     This class is used to define a general linear, steady, second order PDE
1498     for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.     for an unknown function M{u} on a given domain defined through a
1499       L{Domain<escript.Domain>} object.
1500    
1501     For a single PDE with a solution with a single component the linear PDE is defined in the following form:     For a single PDE having a solution with a single component the linear PDE
1502       is defined in the following form:
1503    
1504     M{-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)=-grad(X+X_reduced)[j,j]+(Y+Y_reduced)}     M{-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)=-grad(X+X_reduced)[j,j]+(Y+Y_reduced)}
1505    
1506       where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's
1507     where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,     summation convention, ie. summation over indexes appearing twice in a term
1508     ie. summation over indexes appearing twice in a term of a sum is performed, is used.     of a sum performed, is used.
1509     The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified through L{Data<escript.Data>} objects in the     The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified
1510     L{Function<escript.Function>} and the coefficients M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through L{Data<escript.Data>} objects in the     through L{Data<escript.Data>} objects in L{Function<escript.Function>} and
1511     L{ReducedFunction<escript.ReducedFunction>}. It is also allowd to use objects that can be converted into     the coefficients M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced},
1512     such L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B}, M{C}, M{X},     M{X_reduced} and M{Y_reduced} have to be specified through
1513     M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and M{D}, M{D_reduced} M{Y} and M{Y_reduced} are scalar.     L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}.
1514       It is also allowed to use objects that can be converted into such
1515       L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B},
1516       M{C}, M{X}, M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and
1517       M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are scalar.
1518    
1519     The following natural boundary conditions are considered:     The following natural boundary conditions are considered:
1520    
1521     M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u)+(d+d_reduced)*u=n[j]*(X[j]+X_reduced[j])+y}     M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u)+(d+d_reduced)*u=n[j]*(X[j]+X_reduced[j])+y}
1522    
1523     where M{n} is the outer normal field. Notice that the coefficients M{A}, M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the PDE. The coefficients M{d} and M{y} and are each a scalar in the L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients M{d_reduced} and M{y_reduced} and are each a scalar in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.     where M{n} is the outer normal field. Notice that the coefficients M{A},
1524       M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the
1525       PDE. The coefficients M{d} and M{y} are each a scalar in
1526       L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients
1527       M{d_reduced} and M{y_reduced} are each a scalar in
1528       L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1529    
1530     Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form     Constraints for the solution prescribe the value of the solution at certain
1531       locations in the domain. They have the form
1532    
1533     M{u=r}  where M{q>0}     M{u=r} where M{q>0}
1534    
1535     M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.     M{r} and M{q} are each scalar where M{q} is the characteristic function
1536     The constraints override any other condition set by the PDE or the boundary condition.     defining where the constraint is applied. The constraints override any
1537       other condition set by the PDE or the boundary condition.
1538    
1539     The PDE is symmetrical if     The PDE is symmetrical if
1540    
1541     M{A[i,j]=A[j,i]}  and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]}  and M{B_reduced[j]=C_reduced[j]}     M{A[i,j]=A[j,i]}  and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]}
1542       and M{B_reduced[j]=C_reduced[j]}
1543    
1544     For a system of PDEs and a solution with several components the PDE has the form     For a system of PDEs and a solution with several components the PDE has the
1545       form
1546    
1547     M{-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k] =-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] }     M{-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k] =-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] }
1548    
1549     M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and M{C_reduced} are each of rank three, M{D}, M{D_reduced}, M{X_reduced} and M{X} are each a rank two and M{Y} and M{Y_reduced} are of rank one.     M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and
1550       M{C_reduced} are each of rank three, M{D}, M{D_reduced}, M{X_reduced} and
1551       M{X} are each of rank two and M{Y} and M{Y_reduced} are of rank one.
1552     The natural boundary conditions take the form:     The natural boundary conditions take the form:
1553    
1554     M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])+(d[i,k]+d_reduced[i,k])*u[k]=n[j]*(X[i,j]+X_reduced[i,j])+y[i]+y_reduced[i]}     M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])+(d[i,k]+d_reduced[i,k])*u[k]=n[j]*(X[i,j]+X_reduced[i,j])+y[i]+y_reduced[i]}
1555    
1556       The coefficient M{d} is of rank two and M{y} is of rank one both in
1557       L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients
1558       M{d_reduced} is of rank two and M{y_reduced} is of rank one both in
1559       L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1560    
1561     The coefficient M{d} is a rank two and M{y} is a rank one both in the L{FunctionOnBoundary<escript.FunctionOnBoundary>}. Constraints take the form and the coefficients M{d_reduced} is a rank two and M{y_reduced} is a rank one both in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.     Constraints take the form
   
    Constraints take the form  
1562    
1563     M{u[i]=r[i]}  where  M{q[i]>0}     M{u[i]=r[i]}  where  M{q[i]>0}
1564    
1565     M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.     M{r} and M{q} are each rank one. Notice that at some locations not
1566       necessarily all components must have a constraint.
1567    
1568     The system of PDEs is symmetrical if     The system of PDEs is symmetrical if
1569    
1570          - M{A[i,j,k,l]=A[k,l,i,j]}        - M{A[i,j,k,l]=A[k,l,i,j]}
1571          - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}        - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
1572          - M{B[i,j,k]=C[k,i,j]}        - M{B[i,j,k]=C[k,i,j]}
1573          - M{B_reduced[i,j,k]=C_reduced[k,i,j]}        - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
1574          - M{D[i,k]=D[i,k]}        - M{D[i,k]=D[i,k]}
1575          - M{D_reduced[i,k]=D_reduced[i,k]}        - M{D_reduced[i,k]=D_reduced[i,k]}
1576          - M{d[i,k]=d[k,i]}        - M{d[i,k]=d[k,i]}
1577          - M{d_reduced[i,k]=d_reduced[k,i]}        - M{d_reduced[i,k]=d_reduced[k,i]}
1578    
1579     L{LinearPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the     L{LinearPDE} also supports solution discontinuities over a contact region
1580     discontinuity we are using the generalised flux M{J} which is in the case of a systems of PDEs and several components of the solution     in the domain. To specify the conditions across the discontinuity we are
1581     defined as     using the generalised flux M{J} which, in the case of a system of PDEs
1582       and several components of the solution, is defined as
1583    
1584     M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]-X[i,j]-X_reduced[i,j]}     M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]-X[i,j]-X_reduced[i,j]}
1585    
1586     For the case of single solution component and single PDE M{J} is defined     For the case of single solution component and single PDE M{J} is defined as
1587    
1588     M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}     M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}
1589    
1590     In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1     In the context of discontinuities M{n} denotes the normal on the
1591     calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of PDEs     discontinuity pointing from side 0 towards side 1 calculated from
1592     the contact condition takes the form     L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1593       For a system of PDEs the contact condition takes the form
1594    
1595     M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]}     M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]}
1596    
1597     where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the discontinuity, respectively. M{jump(u)}, which is the difference     where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the
1598     of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calculated by     discontinuity, respectively. M{jump(u)}, which is the difference of the
1599     L{jump<util.jump>}.     solution at side 1 and at side 0, denotes the jump of M{u} across
1600     The coefficient M{d_contact} is a rank two and M{y_contact} is a rank one both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.     discontinuity along the normal calculated by L{jump<util.jump>}.
1601     The coefficient M{d_contact_reduced} is a rank two and M{y_contact_reduced} is a rank one both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.     The coefficient M{d_contact} is of rank two and M{y_contact} is of rank one
1602     In case of a single PDE and a single component solution the contact condition takes the form     both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or
1603       L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1604       The coefficient M{d_contact_reduced} is of rank two and M{y_contact_reduced}
1605       is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}
1606       or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
1607       In case of a single PDE and a single component solution the contact
1608       condition takes the form
1609    
1610     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
1611    
1612     In this case the coefficient M{d_contact} and M{y_contact} are each scalar both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}     In this case the coefficient M{d_contact} and M{y_contact} are each scalar
1613       both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or
1614     typical usage:     L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient
1615       M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in
1616        p=LinearPDE(dom)     L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or
1617        p.setValue(A=kronecker(dom),D=1,Y=0.5)     L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
1618        u=p.getSolution()  
1619       Typical usage::
1620    
1621           p = LinearPDE(dom)
1622           p.setValue(A=kronecker(dom), D=1, Y=0.5)
1623           u = p.getSolution()
1624    
1625     """     """
1626    
   
1627     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1628       """       """
1629       initializes a new linear PDE       Initializes a new linear PDE.
1630    
1631       @param domain: domain of the PDE       @param domain: domain of the PDE
1632       @type domain: L{Domain<escript.Domain>}       @type domain: L{Domain<escript.Domain>}
1633       @param numEquations: number of equations. If numEquations==None the number of equations       @param numEquations: number of equations. If C{None} the number of
1634                            is extracted from the PDE coefficients.                            equations is extracted from the PDE coefficients.
1635       @param numSolutions: number of solution components. If  numSolutions==None the number of solution components       @param numSolutions: number of solution components. If C{None} the number
1636                            is extracted from the PDE coefficients.                            of solution components is extracted from the PDE
1637       @param debug: if True debug informations are printed.                            coefficients.
1638         @param debug: if True debug information is printed
1639    
1640       """       """
1641       super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)       super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
# Line 1543  class LinearPDE(LinearProblem): Line 1668  class LinearPDE(LinearProblem):
1668    
1669     def __str__(self):     def __str__(self):
1670       """       """
1671       returns string representation of the PDE       Returns the string representation of the PDE.
1672    
1673       @return: a simple representation of the PDE       @return: a simple representation of the PDE
1674       @rtype: C{str}       @rtype: C{str}
# Line 1552  class LinearPDE(LinearProblem): Line 1677  class LinearPDE(LinearProblem):
1677    
1678     def getRequiredSystemType(self):     def getRequiredSystemType(self):
1679        """        """
1680        returns the system type which needs to be used by the current set up.        Returns the system type which needs to be used by the current set up.
1681        """        """
1682        return self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())        return self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())
1683    
1684     def checkSymmetry(self,verbose=True):     def checkSymmetry(self,verbose=True):
1685        """        """
1686        test the PDE for symmetry.        Tests the PDE for symmetry.
1687    
1688        @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.        @param verbose: if set to True or not present a report on coefficients
1689                          which break the symmetry is printed.
1690        @type verbose: C{bool}        @type verbose: C{bool}
1691        @return:  True if the PDE is symmetric.        @return: True if the PDE is symmetric
1692        @rtype: L{bool}        @rtype: L{bool}
1693        @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.        @note: This is a very expensive operation. It should be used for
1694                 degugging only! The symmetry flag is not altered.
1695        """        """
1696        out=True        out=True
1697        out=out and self.checkSymmetricTensor("A", verbose)        out=out and self.checkSymmetricTensor("A", verbose)
# Line 1579  class LinearPDE(LinearProblem): Line 1706  class LinearPDE(LinearProblem):
1706        out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)        out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
1707        return out        return out
1708    
1709     def createOperator(self):     def createOperator(self):
1710         """         """
1711         returns an instance of a new operator         Returns an instance of a new operator.
1712         """         """
1713         self.trace("New operator is allocated.")         self.trace("New operator is allocated.")
1714         return self.getDomain().newOperator( \         return self.getDomain().newOperator( \
# Line 1591  class LinearPDE(LinearProblem): Line 1718  class LinearPDE(LinearProblem):
1718                             self.getFunctionSpaceForSolution(), \                             self.getFunctionSpaceForSolution(), \
1719                             self.getSystemType())                             self.getSystemType())
1720    
1721     def getSolution(self,**options):     def getSolution(self,**options):
1722         """         """
1723         returns the solution of the PDE         Returns the solution of the PDE.
1724    
1725         @return: the solution         @return: the solution
1726         @rtype: L{Data<escript.Data>}         @rtype: L{Data<escript.Data>}
1727         @param options: solver options         @param options: solver options
1728         @keyword verbose: True to get some information during PDE solution         @keyword verbose: True to get some information during PDE solution
1729         @type verbose: C{bool}         @type verbose: C{bool}
1730         @keyword reordering: reordering scheme to be used during elimination. Allowed values are         @keyword reordering: reordering scheme to be used during elimination.
1731                              L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}                              Allowed values are L{NO_REORDERING},
1732         @keyword iter_max: maximum number of iteration steps allowed.                              L{MINIMUM_FILL_IN} and L{NESTED_DISSECTION}
1733         @keyword drop_tolerance: threshold for drupping in L{ILUT}         @keyword iter_max: maximum number of iteration steps allowed
1734           @keyword drop_tolerance: threshold for dropping in L{ILUT}
1735         @keyword drop_storage: maximum of allowed memory in L{ILUT}         @keyword drop_storage: maximum of allowed memory in L{ILUT}
1736         @keyword truncation: maximum number of residuals in L{GMRES}         @keyword truncation: maximum number of residuals in L{GMRES}
1737         @keyword restart: restart cycle length in L{GMRES}         @keyword restart: restart cycle length in L{GMRES}
# Line 1623  class LinearPDE(LinearProblem): Line 1751  class LinearPDE(LinearProblem):
1751               self.setSolution(mat.solve(f,options))               self.setSolution(mat.solve(f,options))
1752         return self.getCurrentSolution()         return self.getCurrentSolution()
1753    
1754     def getSystem(self):     def getSystem(self):
1755         """         """
1756         return the operator and right hand side of the PDE         Returns the operator and right hand side of the PDE.
1757    
1758         @return: the discrete version of the PDE         @return: the discrete version of the PDE
1759         @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.         @rtype: C{tuple} of L{Operator,<escript.Operator>} and
1760                   L{Data<escript.Data>}
1761         """         """
1762         if not self.isOperatorValid() or not self.isRightHandSideValid():         if not self.isOperatorValid() or not self.isRightHandSideValid():
1763            if self.isUsingLumping():            if self.isUsingLumping():
1764                if not self.isOperatorValid():                if not self.isOperatorValid():
1765                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1766                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1767                   if not self.getCoefficient("A").isEmpty():                   if not self.getCoefficient("A").isEmpty():
1768                        raise ValueError,"coefficient A in lumped matrix may not be present."                        raise ValueError,"coefficient A in lumped matrix may not be present."
1769                   if not self.getCoefficient("B").isEmpty():                   if not self.getCoefficient("B").isEmpty():
1770                        raise ValueError,"coefficient B in lumped matrix may not be present."                        raise ValueError,"coefficient B in lumped matrix may not be present."
# Line 1643  class LinearPDE(LinearProblem): Line 1772  class LinearPDE(LinearProblem):
1772                        raise ValueError,"coefficient C in lumped matrix may not be present."                        raise ValueError,"coefficient C in lumped matrix may not be present."
1773                   if not self.getCoefficient("d_contact").isEmpty():                   if not self.getCoefficient("d_contact").isEmpty():
1774                        raise ValueError,"coefficient d_contact in lumped matrix may not be present."                        raise ValueError,"coefficient d_contact in lumped matrix may not be present."
1775                   if not self.getCoefficient("A_reduced").isEmpty():                   if not self.getCoefficient("A_reduced").isEmpty():
1776                        raise ValueError,"coefficient A_reduced in lumped matrix may not be present."                        raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
1777                   if not self.getCoefficient("B_reduced").isEmpty():                   if not self.getCoefficient("B_reduced").isEmpty():
1778                        raise ValueError,"coefficient B_reduced in lumped matrix may not be present."                        raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
# Line 1669  class LinearPDE(LinearProblem): Line 1798  class LinearPDE(LinearProblem):
1798                          d_times_e=d                          d_times_e=d
1799                   else:                   else:
1800                      d_times_e=escript.Data()                      d_times_e=escript.Data()
1801          
1802                   if not D_reduced.isEmpty():                   if not D_reduced.isEmpty():
1803                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
1804                          D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))                          D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))
# Line 1715  class LinearPDE(LinearProblem): Line 1844  class LinearPDE(LinearProblem):
1844                                 self.getCoefficient("Y_reduced"),\                                 self.getCoefficient("Y_reduced"),\
1845                                 self.getCoefficient("y_reduced"),\                                 self.getCoefficient("y_reduced"),\
1846                                 self.getCoefficient("y_contact_reduced"))                                 self.getCoefficient("y_contact_reduced"))
1847                   self.trace("New right hand side as been built.")                   self.trace("New right hand side has been built.")
1848                   self.validRightHandSide()                   self.validRightHandSide()
1849                self.insertConstraint(rhs_only=False)                self.insertConstraint(rhs_only=False)
1850                self.validOperator()                self.validOperator()
# Line 1764  class LinearPDE(LinearProblem): Line 1893  class LinearPDE(LinearProblem):
1893                                 self.getCoefficient("Y_reduced"),\                                 self.getCoefficient("Y_reduced"),\
1894                                 self.getCoefficient("y_reduced"),\                                 self.getCoefficient("y_reduced"),\
1895                                 self.getCoefficient("y_contact_reduced"))                                 self.getCoefficient("y_contact_reduced"))
1896                   self.insertConstraint(rhs_only=True)                   self.insertConstraint(rhs_only=True)
1897                   self.trace("New right hand side has been built.")                   self.trace("New right hand side has been built.")
1898                   self.validRightHandSide()                   self.validRightHandSide()
1899               elif not self.isOperatorValid():               elif not self.isOperatorValid():
# Line 1797  class LinearPDE(LinearProblem): Line 1926  class LinearPDE(LinearProblem):
1926                   self.validOperator()                   self.validOperator()
1927         return (self.getCurrentOperator(), self.getCurrentRightHandSide())         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1928    
1929     def insertConstraint(self, rhs_only=False):     def insertConstraint(self, rhs_only=False):
1930        """        """
1931        applies the constraints defined by q and r to the PDE        Applies the constraints defined by q and r to the PDE.
1932          
1933        @param rhs_only: if True only the right hand side is altered by the constraint.        @param rhs_only: if True only the right hand side is altered by the
1934                           constraint
1935        @type rhs_only: C{bool}        @type rhs_only: C{bool}
1936        """        """
1937        q=self.getCoefficient("q")        q=self.getCoefficient("q")
# Line 1826  class LinearPDE(LinearProblem): Line 1956  class LinearPDE(LinearProblem):
1956                   operator.nullifyRowsAndCols(row_q,col_q,1.)                   operator.nullifyRowsAndCols(row_q,col_q,1.)
1957           righthandside.copyWithMask(r_s,q)           righthandside.copyWithMask(r_s,q)
1958    
1959     def setValue(self,**coefficients):     def setValue(self,**coefficients):
1960        """        """
1961        sets new values to coefficients        Sets new values to coefficients.
1962    
1963        @param coefficients: new values assigned to coefficients        @param coefficients: new values assigned to coefficients
1964        @keyword A: value for coefficient A.        @keyword A: value for coefficient C{A}
1965        @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @type A: any type that can be cast to a L{Data<escript.Data>} object on
1966        @keyword A_reduced: value for coefficient A_reduced.                 L{Function<escript.Function>}
1967        @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.        @keyword A_reduced: value for coefficient C{A_reduced}
1968        @keyword B: value for coefficient B        @type A_reduced: any type that can be cast to a L{Data<escript.Data>}
1969        @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.                         object on L{ReducedFunction<escript.ReducedFunction>}
1970        @keyword B_reduced: value for coefficient B_reduced        @keyword B: value for coefficient C{B}
1971        @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.        @type B: any type that can be cast to a L{Data<escript.Data>} object on
1972        @keyword C: value for coefficient C                 L{Function<escript.Function>}
1973        @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @keyword B_reduced: value for coefficient C{B_reduced}
1974        @keyword C_reduced: value for coefficient C_reduced        @type B_reduced: any type that can be cast to a L{Data<escript.Data>}
1975        @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.                         object on L{ReducedFunction<escript.ReducedFunction>}
1976        @keyword D: value for coefficient D        @keyword C: value for coefficient C{C}
1977        @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @type C: any type that can be cast to a L{Data<escript.Data>} object on
1978        @keyword D_reduced: value for coefficient D_reduced                 L{Function<escript.Function>}
1979        @type D_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.        @keyword C_reduced: value for coefficient C{C_reduced}
1980        @keyword X: value for coefficient X        @type C_reduced: any type that can be cast to a L{Data<escript.Data>}
1981        @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.                         object on L{ReducedFunction<escript.ReducedFunction>}
1982        @keyword X_reduced: value for coefficient X_reduced        @keyword D: value for coefficient C{D}
1983        @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.        @type D: any type that can be cast to a L{Data<escript.Data>} object on
1984        @keyword Y: value for coefficient Y                 L{Function<escript.Function>}
1985        @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @keyword D_reduced: value for coefficient C{D_reduced}
1986        @keyword Y_reduced: value for coefficient Y_reduced        @type D_reduced: any type that can be cast to a L{Data<escript.Data>}
1987        @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.                         object on L{ReducedFunction<escript.ReducedFunction>}
1988        @keyword d: value for coefficient d        @keyword X: value for coefficient C{X}
1989        @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.        @type X: any type that can be cast to a L{Data<escript.Data>} object on
1990        @keyword d_reduced: value for coefficient d_reduced                 L{Function<escript.Function>}
1991        @type d_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.        @keyword X_reduced: value for coefficient C{X_reduced}
1992        @keyword y: value for coefficient y        @type X_reduced: any type that can be cast to a L{Data<escript.Data>}
1993        @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.                         object on L{ReducedFunction<escript.ReducedFunction>}
1994        @keyword d_contact: value for coefficient d_contact        @keyword Y: value for coefficient C{Y}
1995        @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.        @type Y: any type that can be cast to a L{Data<escript.Data>} object on
1996        @keyword d_contact_reduced: value for coefficient d_contact_reduced                 L{Function<escript.Function>}
1997        @type d_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>} or  L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}.        @keyword Y_reduced: value for coefficient C{Y_reduced}
1998        @keyword y_contact: value for coefficient y_contact        @type Y_reduced: any type that can be cast to a L{Data<escript.Data>}
1999        @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.                         object on L{ReducedFunction<escript.Function>}
2000        @keyword y_contact_reduced: value for coefficient y_contact_reduced        @keyword d: value for coefficient C{d}
2001        @type y_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}.        @type d: any type that can be cast to a L{Data<escript.Data>} object on
2002        @keyword r: values prescribed to the solution at the locations of constraints                 L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2003        @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}        @keyword d_reduced: value for coefficient C{d_reduced}
2004                 depending of reduced order is used for the solution.        @type d_reduced: any type that can be cast to a L{Data<escript.Data>}
2005                           object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}
2006          @keyword y: value for coefficient C{y}
2007          @type y: any type that can be cast to a L{Data<escript.Data>} object on
2008                   L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2009          @keyword d_contact: value for coefficient C{d_contact}
2010          @type d_contact: any type that can be cast to a L{Data<escript.Data>}
2011                           object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}
2012                           or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
2013          @keyword d_contact_reduced: value for coefficient C{d_contact_reduced}
2014          @type d_contact_reduced: any type that can be cast to a L{Data<escript.Data>}
2015                                   object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}
2016                                   or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}
2017          @keyword y_contact: value for coefficient C{y_contact}
2018          @type y_contact: any type that can be cast to a L{Data<escript.Data>}
2019                           object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}
2020                           or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
2021          @keyword y_contact_reduced: value for coefficient C{y_contact_reduced}
2022          @type y_contact_reduced: any type that can be cast to a L{Data<escript.Data>}
2023                                   object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>}
2024                                   or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}
2025          @keyword r: values prescribed to the solution at the locations of
2026                      constraints
2027          @type r: any type that can be cast to a L{Data<escript.Data>} object on
2028                   L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2029                   depending on whether reduced order is used for the solution
2030        @keyword q: mask for location of constraints        @keyword q: mask for location of constraints
2031        @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}        @type q: any type that can be cast to a L{Data<escript.Data>} object on
2032                 depending of reduced order is used for the representation of the equation.                 L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2033        @raise IllegalCoefficient: if an unknown coefficient keyword is used.                 depending on whether reduced order is used for the
2034                   representation of the equation
2035          @raise IllegalCoefficient: if an unknown coefficient keyword is used
2036        """        """
2037        super(LinearPDE,self).setValue(**coefficients)        super(LinearPDE,self).setValue(**coefficients)
2038        # check if the systrem is inhomogeneous:        # check if the systrem is inhomogeneous:
# Line 1888  class LinearPDE(LinearProblem): Line 2045  class LinearPDE(LinearProblem):
2045                 self.invalidateSystem()                 self.invalidateSystem()
2046    
2047    
2048     def getResidual(self,u=None):     def getResidual(self,u=None):
2049       """       """
2050       return the residual of u or the current solution if u is not present.       Returns the residual of u or the current solution if u is not present.
2051    
2052       @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}       @param u: argument in the residual calculation. It must be representable
2053                 the current solution is used.                 in L{self.getFunctionSpaceForSolution()}. If u is not present
2054                   or equals C{None} the current solution is used.
2055       @type u: L{Data<escript.Data>} or None       @type u: L{Data<escript.Data>} or None
2056       @return: residual of u       @return: residual of u
2057       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
# Line 1903  class LinearPDE(LinearProblem): Line 2061  class LinearPDE(LinearProblem):
2061       else:       else:
2062          return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())-self.getRightHandSide()          return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())-self.getRightHandSide()
2063    
2064     def getFlux(self,u=None):     def getFlux(self,u=None):
2065       """       """
2066       returns the flux M{J} for a given M{u}       Returns the flux M{J} for a given M{u}.
2067    
2068       M{J[i,j]=(A[i,j,k,l]+A_reduced[A[i,j,k,l]]*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])u[k]-X[i,j]-X_reduced[i,j]}       M{J[i,j]=(A[i,j,k,l]+A_reduced[A[i,j,k,l]]*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])u[k]-X[i,j]-X_reduced[i,j]}
2069    
# Line 1913  class LinearPDE(LinearProblem): Line 2071  class LinearPDE(LinearProblem):
2071    
2072       M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}       M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}
2073    
2074       @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.       @param u: argument in the flux. If u is not present or equals L{None} the
2075                   current solution is used.
2076       @type u: L{Data<escript.Data>} or None       @type u: L{Data<escript.Data>} or None
2077       @return: flux       @return: flux
2078       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
# Line 1929  class LinearPDE(LinearProblem): Line 2088  class LinearPDE(LinearProblem):
2088    
2089  class Poisson(LinearPDE):  class Poisson(LinearPDE):
2090     """     """
2091     Class to define a Poisson equation problem, which is genear L{LinearPDE} of the form     Class to define a Poisson equation problem. This is generally a
2092       L{LinearPDE} of the form
2093    
2094     M{-grad(grad(u)[j])[j] = f}     M{-grad(grad(u)[j])[j] = f}
2095    
2096     with natural boundary conditons     with natural boundary conditions
2097    
2098     M{n[j]*grad(u)[j] = 0 }     M{n[j]*grad(u)[j] = 0 }
2099    
# Line 1945  class Poisson(LinearPDE): Line 2105  class Poisson(LinearPDE):
2105    
2106     def __init__(self,domain,debug=False):     def __init__(self,domain,debug=False):
2107       """       """
2108       initializes a new Poisson equation       Initializes a new Poisson equation.
2109    
2110       @param domain: domain of the PDE       @param domain: domain of the PDE
2111       @type domain: L{Domain<escript.Domain>}       @type domain: L{Domain<escript.Domain>}
2112       @param debug: if True debug informations are printed.       @param debug: if True debug information is printed
2113    
2114       """       """
2115       super(Poisson, self).__init__(domain,1,1,debug)       super(Poisson, self).__init__(domain,1,1,debug)
# Line 1960  class Poisson(LinearPDE): Line 2120  class Poisson(LinearPDE):
2120    
2121     def setValue(self,**coefficients):     def setValue(self,**coefficients):
2122       """       """
2123       sets new values to coefficients       Sets new values to coefficients.
2124    
2125       @param coefficients: new values assigned to coefficients       @param coefficients: new values assigned to coefficients
2126       @keyword f: value for right hand side M{f}       @keyword f: value for right hand side M{f}
2127       @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.       @type f: any type that can be cast to a L{Scalar<escript.Scalar>} object
2128                  on L{Function<escript.Function>}
2129       @keyword q: mask for location of constraints       @keyword q: mask for location of constraints
2130       @type q: any type that can be casted to rank zeo L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}       @type q: any type that can be cast to a rank zero L{Data<escript.Data>}
2131                 depending of reduced order is used for the representation of the equation.                object on L{Solution<escript.Solution>} or
2132       @raise IllegalCoefficient: if an unknown coefficient keyword is used.                L{ReducedSolution<escript.ReducedSolution>} depending on whether
2133                  reduced order is used for the representation of the equation
2134         @raise IllegalCoefficient: if an unknown coefficient keyword is used
2135       """       """
2136       super(Poisson, self).setValue(**coefficients)       super(Poisson, self).setValue(**coefficients)
2137    
2138    
2139     def getCoefficient(self,name):     def getCoefficient(self,name):
2140       """       """
2141       return the value of the coefficient name of the general PDE       Returns the value of the coefficient C{name} of the general PDE.
2142       @param name: name of the coefficient requested.  
2143         @param name: name of the coefficient requested
2144       @type name: C{string}       @type name: C{string}
2145       @return: the value of the coefficient  name       @return: the value of the coefficient C{name}
2146       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
2147       @raise IllegalCoefficient: invalid coefficient name       @raise IllegalCoefficient: invalid coefficient name
2148       @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.       @note: This method is called by the assembling routine to map the Poisson
2149                equation onto the general PDE.
2150       """       """
2151       if name == "A" :       if name == "A" :
2152           return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))           return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
# Line 1990  class Poisson(LinearPDE): Line 2155  class Poisson(LinearPDE):
2155       elif name == "Y_reduced" :       elif name == "Y_reduced" :
2156           return self.getCoefficient("f_reduced")           return self.getCoefficient("f_reduced")
2157       else:       else:
2158          return super(Poisson, self).getCoefficient(name)           return super(Poisson, self).getCoefficient(name)
2159    
2160  class Helmholtz(LinearPDE):  class Helmholtz(LinearPDE):
2161     """     """
2162     Class to define a Helmhotz equation problem, which is genear L{LinearPDE} of the form     Class to define a Helmholtz equation problem. This is generally a
2163       L{LinearPDE} of the form
2164    
2165     M{S{omega}*u - grad(k*grad(u)[j])[j] = f}     M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
2166    
2167     with natural boundary conditons     with natural boundary conditions
2168    
2169     M{k*n[j]*grad(u)[j] = g- S{alpha}u }     M{k*n[j]*grad(u)[j] = g- S{alpha}u }
2170    
# Line 2010  class Helmholtz(LinearPDE): Line 2176  class Helmholtz(LinearPDE):
2176    
2177     def __init__(self,domain,debug=False):     def __init__(self,domain,debug=False):
2178       """       """
2179       initializes a new Poisson equation       Initializes a new Helmholtz equation.
2180    
2181       @param domain: domain of the PDE       @param domain: domain of the PDE
2182       @type domain: L{Domain<escript.Domain>}       @type domain: L{Domain<escript.Domain>}
2183       @param debug: if True debug informations are printed.       @param debug: if True debug information is printed
2184    
2185       """       """
2186       super(Helmholtz, self).__init__(domain,1,1,debug)       super(Helmholtz, self).__init__(domain,1,1,debug)
# Line 2028  class Helmholtz(LinearPDE): Line 2194  class Helmholtz(LinearPDE):
2194                          g_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))                          g_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2195       self.setSymmetryOn()       self.setSymmetryOn()
2196    
2197     def setValue(self,**coefficients):     def setValue(self,**coefficients):
2198       """       """
2199       sets new values to coefficients       Sets new values to coefficients.
2200    
2201       @param coefficients: new values assigned to coefficients       @param coefficients: new values assigned to coefficients
2202       @keyword omega: value for coefficient M{S{omega}}       @keyword omega: value for coefficient M{S{omega}}
2203       @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.       @type omega: any type that can be cast to a L{Scalar<escript.Scalar>}
2204       @keyword k: value for coefficeint M{k}                    object on L{Function<escript.Function>}
2205       @type k: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.       @keyword k: value for coefficient M{k}
2206         @type k: any type that can be cast to a L{Scalar<escript.Scalar>} object
2207                  on L{Function<escript.Function>}
2208       @keyword f: value for right hand side M{f}       @keyword f: value for right hand side M{f}
2209       @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.       @type f: any type that can be cast to a L{Scalar<escript.Scalar>} object
2210                  on L{Function<escript.Function>}
2211       @keyword alpha: value for right hand side M{S{alpha}}       @keyword alpha: value for right hand side M{S{alpha}}
2212       @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.       @type alpha: any type that can be cast to a L{Scalar<escript.Scalar>}
2213                      object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2214       @keyword g: value for right hand side M{g}       @keyword g: value for right hand side M{g}
2215       @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.       @type g: any type that can be cast to a L{Scalar<escript.Scalar>} object
2216       @keyword r: prescribed values M{r} for the solution in constraints.                on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2217       @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}       @keyword r: prescribed values M{r} for the solution in constraints
2218                 depending of reduced order is used for the representation of the equation.       @type r: any type that can be cast to a L{Scalar<escript.Scalar>} object
2219       @keyword q: mask for location of constraints                on L{Solution<escript.Solution>} or
2220       @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}                L{ReducedSolution<escript.ReducedSolution>} depending on whether
2221                 depending of reduced order is used for the representation of the equation.                reduced order is used for the representation of the equation
2222       @raise IllegalCoefficient: if an unknown coefficient keyword is used.       @keyword q: mask for the location of constraints
2223         @type q: any type that can be cast to a L{Scalar<escript.Scalar>} object
2224                  on L{Solution<escript.Solution>} or
2225                  L{ReducedSolution<escript.ReducedSolution>} depending on whether
2226                  reduced order is used for the representation of the equation
2227         @raise IllegalCoefficient: if an unknown coefficient keyword is used
2228       """       """
2229       super(Helmholtz, self).setValue(**coefficients)       super(Helmholtz, self).setValue(**coefficients)
2230    
2231     def getCoefficient(self,name):     def getCoefficient(self,name):
2232       """       """
2233       return the value of the coefficient name of the general PDE       Returns the value of the coefficient C{name} of the general PDE.
2234    
2235       @param name: name of the coefficient requested.       @param name: name of the coefficient requested
2236       @type name: C{string}       @type name: C{string}
2237       @return: the value of the coefficient  name       @return: the value of the coefficient C{name}
2238       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
2239       @raise IllegalCoefficient: invalid name       @raise IllegalCoefficient: invalid name
2240       """       """
# Line 2082  class Helmholtz(LinearPDE): Line 2257  class Helmholtz(LinearPDE):
2257    
2258  class LameEquation(LinearPDE):  class LameEquation(LinearPDE):
2259     """     """
2260     Class to define a Lame equation problem:     Class to define a Lame equation problem. This problem is defined as:
2261    
2262     M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[k])[k])[j] = F_i -grad(S{sigma}[ij])[j] }     M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[k])[k])[j] = F_i -grad(S{sigma}[ij])[j] }
2263    
2264     with natural boundary conditons:     with natural boundary conditions:
2265    
2266     M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) + S{lambda}*grad(u[k])[k]) = f_i +n[j]*S{sigma}[ij] }     M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) + S{lambda}*grad(u[k])[k]) = f_i +n[j]*S{sigma}[ij] }
2267    
# Line 2097  class LameEquation(LinearPDE): Line 2272  class LameEquation(LinearPDE):
2272     """     """
2273    
2274     def __init__(self,domain,debug=False):     def __init__(self,domain,debug=False):
2275          """
2276          Initializes a new Lame equation.
2277    
2278          @param domain: domain of the PDE
2279          @type domain: L{Domain<escript.Domain>}
2280          @param debug: if True debug information is printed
2281    
2282          """
2283        super(LameEquation, self).__init__(domain,\        super(LameEquation, self).__init__(domain,\
2284                                           domain.getDim(),domain.getDim(),debug)                                           domain.getDim(),domain.getDim(),debug)
2285        self.introduceCoefficients(lame_lambda=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),        self.introduceCoefficients(lame_lambda=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
# Line 2106  class LameEquation(LinearPDE): Line 2289  class LameEquation(LinearPDE):
2289                                   f=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))                                   f=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2290        self.setSymmetryOn()        self.setSymmetryOn()
2291    
2292     def setValues(self,**coefficients):     def setValues(self,**coefficients):
2293       """       """
2294       sets new values to coefficients       Sets new values to coefficients.
2295    
2296       @param coefficients: new values assigned to coefficients       @param coefficients: new values assigned to coefficients
2297       @keyword lame_mu: value for coefficient M{S{mu}}       @keyword lame_mu: value for coefficient M{S{mu}}
2298       @type lame_mu: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.       @type lame_mu: any type that can be cast to a L{Scalar<escript.Scalar>}
2299                        object on L{Function<escript.Function>}
2300       @keyword lame_lambda: value for coefficient M{S{lambda}}       @keyword lame_lambda: value for coefficient M{S{lambda}}
2301       @type lame_lambda: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.       @type lame_lambda: any type that can be cast to a L{Scalar<escript.Scalar>}
2302                            object on L{Function<escript.Function>}
2303       @keyword F: value for internal force M{F}       @keyword F: value for internal force M{F}
2304       @type F: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}       @type F: any type that can be cast to a L{Vector<escript.Vector>} object
2305                  on L{Function<escript.Function>}
2306       @keyword sigma: value for initial stress M{S{sigma}}       @keyword sigma: value for initial stress M{S{sigma}}
2307       @type sigma: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}       @type sigma: any type that can be cast to a L{Tensor<escript.Tensor>}
2308       @keyword f: value for extrenal force M{f}                    object on L{Function<escript.Function>}
2309       @type f: any type that can be casted to L{Vector<escript.Vector>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}       @keyword f: value for external force M{f}
2310       @keyword r: prescribed values M{r} for the solution in constraints.       @type f: any type that can be cast to a L{Vector<escript.Vector>} object
2311       @type r: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}                on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2312                 depending of reduced order is used for the representation of the equation.       @keyword r: prescribed values M{r} for the solution in constraints
2313       @keyword q: mask for location of constraints       @type r: any type that can be cast to a L{Vector<escript.Vector>} object
2314       @type q: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}                on L{Solution<escript.Solution>} or
2315                 depending of reduced order is used for the representation of the equation.                L{ReducedSolution<escript.ReducedSolution>} depending on whether
2316       @raise IllegalCoefficient: if an unknown coefficient keyword is used.                reduced order is used for the representation of the equation
2317         @keyword q: mask for the location of constraints
2318         @type q: any type that can be cast to a L{Vector<escript.Vector>} object
2319                  on L{Solution<escript.Solution>} or
2320                  L{ReducedSolution<escript.ReducedSolution>} depending on whether
2321                  reduced order is used for the representation of the equation
2322         @raise IllegalCoefficient: if an unknown coefficient keyword is used
2323       """       """
2324       super(LameEquation, self).setValues(**coefficients)       super(LameEquation, self).setValues(**coefficients)
2325    
2326     def getCoefficient(self,name):     def getCoefficient(self,name):
2327       """       """
2328       return the value of the coefficient name of the general PDE       Returns the value of the coefficient C{name} of the general PDE.
2329    
2330       @param name: name of the coefficient requested.       @param name: name of the coefficient requested
2331       @type name: C{string}       @type name: C{string}
2332       @return: the value of the coefficient  name       @return: the value of the coefficient C{name}
2333       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
2334       @raise IllegalCoefficient: invalid coefficient name       @raise IllegalCoefficient: invalid coefficient name
2335       """       """
# Line 2160  class LameEquation(LinearPDE): Line 2352  class LameEquation(LinearPDE):
2352    
2353  def LinearSinglePDE(domain,debug=False):  def LinearSinglePDE(domain,debug=False):
2354     """     """
2355     defines a single linear PDEs     Defines a single linear PDE.
2356    
2357     @param domain: domain of the PDE     @param domain: domain of the PDE
2358     @type domain: L{Domain<escript.Domain>}     @type domain: L{Domain<escript.Domain>}
2359     @param debug: if True debug informations are printed.     @param debug: if True debug information is printed
2360     @rtype: L{LinearPDE}     @rtype: L{LinearPDE}
2361     """     """
2362     return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)     return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
2363    
2364  def LinearPDESystem(domain,debug=False):  def LinearPDESystem(domain,debug=False):
2365     """     """
2366     defines a system of linear PDEs     Defines a system of linear PDEs.
2367    
2368     @param domain: domain of the PDE     @param domain: domain of the PDEs
2369     @type domain: L{Domain<escript.Domain>}     @type domain: L{Domain<escript.Domain>}
2370     @param debug: if True debug informations are printed.     @param debug: if True debug information is printed
2371     @rtype: L{LinearPDE}     @rtype: L{LinearPDE}
2372     """     """
2373     return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)     return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
2374    
2375    
2376  class TransportPDE(LinearProblem):  class TransportPDE(LinearProblem):
2377     """     """
2378     This class is used to define a transport problem given by a general linear, time dependend, second order PDE     This class is used to define a transport problem given by a general linear,
2379     for an unknown, non-negative, time-dependend function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.     time dependent, second order PDE for an unknown, non-negative,
2380       time-dependent function M{u} on a given domain defined through a
2381       L{Domain<escript.Domain>} object.
2382    
2383     For a single equation with a solution with a single component the transport problem is defined in the following form:     For a single equation with a solution with a single component the transport
2384       problem is defined in the following form:
2385    
2386     M{(M+M_reduced)*u_t=-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)-grad(X+X_reduced)[j,j]+(Y+Y_reduced)}     M{(M+M_reduced)*u_t=-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)-grad(X+X_reduced)[j,j]+(Y+Y_reduced)}
2387    
2388       where M{u_t} denotes the time derivative of M{u} and M{grad(F)} denotes the
2389     where M{u_t} denotes the time derivative of M{u} and M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,     spatial derivative of M{F}. Einstein's summation convention,  ie. summation
2390     ie. summation over indexes appearing twice in a term of a sum is performed, is used.     over indexes appearing twice in a term of a sum performed, is used.
2391     The coefficients M{M}, M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified through L{Data<escript.Data>} objects in the     The coefficients M{M}, M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be
2392     L{Function<escript.Function>} and the coefficients M{M_reduced}, M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through L{Data<escript.Data>} objects in the     specified through L{Data<escript.Data>} objects in L{Function<escript.Function>}
2393     L{ReducedFunction<escript.ReducedFunction>}. It is also allowed to use objects that can be converted into     and the coefficients M{M_reduced}, M{A_reduced}, M{B_reduced}, M{C_reduced},
2394     such L{Data<escript.Data>} objects. M{M} and M{M_reduced} are scalar, M{A} and M{A_reduced} are rank two,     M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through
2395     M{B}, M{C}, M{X}, M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are scalar.     L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}.
2396       It is also allowed to use objects that can be converted into such
2397       L{Data<escript.Data>} objects. M{M} and M{M_reduced} are scalar, M{A} and
2398       M{A_reduced} are rank two, M{B}, M{C}, M{X}, M{B_reduced}, M{C_reduced} and
2399       M{X_reduced} are rank one and M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are
2400       scalar.
2401    
2402     The following natural boundary conditions are considered:     The following natural boundary conditions are considered:
2403    
2404     M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u+X[j]+X_reduced[j])+(d+d_reduced)*u+y+y_reduced=(m+m_reduced)*u_t}     M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u+X[j]+X_reduced[j])+(d+d_reduced)*u+y+y_reduced=(m+m_reduced)*u_t}
2405    
2406     where M{n} is the outer normal field. Notice that the coefficients M{A}, M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the transport problem. The coefficients M{m}, M{d} and M{y} and are each a scalar in the L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients M{m_reduced}, M{d_reduced} and M{y_reduced} and are each a scalar in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.     where M{n} is the outer normal field. Notice that the coefficients M{A},
2407       M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the
2408     Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form     transport problem. The coefficients M{m}, M{d} and M{y} are each a scalar in
2409       L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients
2410     M{u_t=r}  where M{q>0}     M{m_reduced}, M{d_reduced} and M{y_reduced} are each a scalar in
2411       L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2412     M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.  
2413     The constraints override any other condition set by the transport problem or the boundary condition.     Constraints for the solution prescribing the value of the solution at
2414       certain locations in the domain have the form
2415    
2416       M{u_t=r} where M{q>0}
2417    
2418       M{r} and M{q} are each scalar where M{q} is the characteristic function
2419       defining where the constraint is applied. The constraints override any other
2420       condition set by the transport problem or the boundary condition.
2421    
2422     The transport problem is symmetrical if     The transport problem is symmetrical if
2423    
2424     M{A[i,j]=A[j,i]}  and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]}  and M{B_reduced[j]=C_reduced[j]}     M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]} and
2425       M{B_reduced[j]=C_reduced[j]}
2426    
2427     For a system and a solution with several components the transport problem has the form     For a system and a solution with several components the transport problem
2428       has the form
2429    
2430     M{(M[i,k]+M_reduced[i,k])*u[k]_t=-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k]-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] }     M{(M[i,k]+M_reduced[i,k])*u[k]_t=-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k]-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] }
2431    
2432     M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and M{C_reduced} are each of rank three, M{M}, M{M_reduced}, M{D}, M{D_reduced}, M{X_reduced} and M{X} are each a rank two and M{Y} and M{Y_reduced} are of rank one. The natural boundary conditions take the form:     M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and
2433       M{C_reduced} are each of rank three, M{M}, M{M_reduced}, M{D}, M{D_reduced},
2434       M{X_reduced} and M{X} are each of rank two and M{Y} and M{Y_reduced} are of
2435       rank one. The natural boundary conditions take the form:
2436    
2437     M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j])+(d[i,k]+d_reduced[i,k])*u[k]+y[i]+y_reduced[i]= (m[i,k]+m_reduced[i,k])*u[k]_t}     M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j])+(d[i,k]+d_reduced[i,k])*u[k]+y[i]+y_reduced[i]= (m[i,k]+m_reduced[i,k])*u[k]_t}
2438    
2439     The coefficient M{d} and M{m} are of rank two and M{y} are of rank one with all in the L{FunctionOnBoundary<escript.FunctionOnBoundary>}. Constraints take the form and the coefficients M{d_reduced} and M{m_reduced} are of rank two and M{y_reduced} is of rank one with all in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.     The coefficient M{d} and M{m} are of rank two and M{y} is of rank one with
2440       all in L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients
2441       M{d_reduced} and M{m_reduced} are of rank two and M{y_reduced} is of rank
2442       one all in L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2443    
2444     Constraints take the form     Constraints take the form
2445    
2446     M{u[i]_t=r[i]}  where  M{q[i]>0}     M{u[i]_t=r[i]} where M{q[i]>0}
2447    
2448     M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.     M{r} and M{q} are each rank one. Notice that at some locations not
2449       necessarily all components must have a constraint.
2450    
2451     The transport problem is symmetrical if     The transport problem is symmetrical if
2452    
2453          - M{M[i,k]=M[i,k]}        - M{M[i,k]=M[i,k]}
2454          - M{M_reduced[i,k]=M_reduced[i,k]}        - M{M_reduced[i,k]=M_reduced[i,k]}
2455          - M{A[i,j,k,l]=A[k,l,i,j]}        - M{A[i,j,k,l]=A[k,l,i,j]}
2456          - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}        - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
2457          - M{B[i,j,k]=C[k,i,j]}        - M{B[i,j,k]=C[k,i,j]}
2458          - M{B_reduced[i,j,k]=C_reduced[k,i,j]}        - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
2459          - M{D[i,k]=D[i,k]}        - M{D[i,k]=D[i,k]}
2460          - M{D_reduced[i,k]=D_reduced[i,k]}        - M{D_reduced[i,k]=D_reduced[i,k]}
2461          - M{m[i,k]=m[k,i]}        - M{m[i,k]=m[k,i]}
2462          - M{m_reduced[i,k]=m_reduced[k,i]}        - M{m_reduced[i,k]=m_reduced[k,i]}
2463          - M{d[i,k]=d[k,i]}        - M{d[i,k]=d[k,i]}
2464          - M{d_reduced[i,k]=d_reduced[k,i]}        - M{d_reduced[i,k]=d_reduced[k,i]}
2465    
2466     L{TransportPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the     L{TransportPDE} also supports solution discontinuities over a contact region
2467     discontinuity we are using the generalised flux M{J} which is in the case of a systems of PDEs and several components of the solution     in the domain. To specify the conditions across the discontinuity we are
2468     defined as     using the generalised flux M{J} which, in the case of a system of PDEs and
2469       several components of the solution, is defined as
2470    
2471     M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j]}     M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j]}
2472    
2473     For the case of single solution component and single PDE M{J} is defined     For the case of single solution component and single PDE M{J} is defined as
2474    
2475     M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]}     M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]}
2476    
2477     In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1     In the context of discontinuities M{n} denotes the normal on the
2478     calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of transport problems the contact condition takes the form     discontinuity pointing from side 0 towards side 1 calculated from
2479       L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
2480       For a system of transport problems the contact condition takes the form
2481    
2482     M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]}     M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]}
2483    
2484     where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the discontinuity, respectively. M{jump(u)}, which is the difference     where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the
2485     of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calculated by     discontinuity, respectively. M{jump(u)}, which is the difference of the
2486     L{jump<util.jump>}.     solution at side 1 and at side 0, denotes the jump of M{u} across
2487     The coefficient M{d_contact} is a rank two and M{y_contact} is a rank one both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.     discontinuity along the normal calculated by L{jump<util.jump>}.
2488     The coefficient M{d_contact_reduced} is a rank two and M{y_contact_reduced} is a rank one both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.     The coefficient M{d_contact} is of rank two and M{y_contact} is of rank one
2489     In case of a single PDE and a single component solution the contact condition takes the form     both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
2490       The coefficient M{d_contact_reduced} is of rank two and M{y_contact_reduced}
2491       is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
2492       In case of a single PDE and a single component solution the contact
2493       condition takes the form
2494    
2495     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
2496    
2497     In this case the coefficient M{d_contact} and M{y_contact} are each scalar both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}     In this case the coefficient M{d_contact} and M{y_contact} are each scalar
2498       both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or
2499     typical usage:     L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient
2500       M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in
2501        p=TransportPDE(dom)     L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or
2502        p.setValue(M=1.,C=[-1.,0.])     L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
2503        p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)  
2504        t=0     Typical usage::
2505        dt=0.1  
2506        while (t<1.):         p = TransportPDE(dom)
2507            u=p.solve(dt)         p.setValue(M=1., C=[-1.,0.])
2508           p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)
2509           t = 0
2510           dt = 0.1
2511           while (t < 1.):
2512               u = p.solve(dt)
2513    
2514     """     """
2515     def __init__(self,domain,numEquations=None,numSolutions=None, theta=0.5,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None, theta=0.5,debug=False):
2516       """       """
2517       initializes a linear problem       Initializes a transport problem.
2518    
2519       @param domain: domain of the PDE       @param domain: domain of the PDE
2520       @type domain: L{Domain<escript.Domain>}       @type domain: L{Domain<escript.Domain>}
2521       @param numEquations: number of equations. If numEquations==None the number of equations       @param numEquations: number of equations. If C{None} the number of
2522                            is extracted from the coefficients.                            equations is extracted from the coefficients.
2523       @param numSolutions: number of solution components. If  numSolutions==None the number of solution components       @param numSolutions: number of solution components. If C{None} the number
2524                            is extracted from the coefficients.                            of solution components is extracted from the
2525       @param debug: if True debug informations are printed.                            coefficients.
2526       @param theta: time stepping control: =1 backward Euler, =0 forward Euler, =0.5 Crank-Nicholson scheme, U{see here<http://en.wikipedia.org/wiki/Crank-Nicolson_method>}.       @param debug: if True debug information is printed
2527         @param theta: time stepping control:
2528            - 1.0: backward Euler,
2529            - 0.0: forward Euler,
2530            - 0.5: Crank-Nicholson scheme, U{see here<http://en.wikipedia.org/wiki/Crank-Nicolson_method>}
2531    
2532       """       """
2533       if theta<0 or theta>1:       if theta<0 or theta>1:
# Line 2335  class TransportPDE(LinearProblem): Line 2568  class TransportPDE(LinearProblem):
2568    
2569     def __str__(self):     def __str__(self):
2570       """       """
2571       returns string representation of the transport problem       Returns the string representation of the transport problem.
2572    
2573       @return: a simple representation of the transport problem       @return: a simple representation of the transport problem
2574       @rtype: C{str}       @rtype: C{str}
# Line 2344  class TransportPDE(LinearProblem): Line 2577  class TransportPDE(LinearProblem):
2577    
2578     def getTheta(self):     def getTheta(self):
2579        """        """
2580        returns the time stepping control parameter        Returns the time stepping control parameter.
2581        @rtype: float        @rtype: float
2582        """        """
2583        return self.__theta        return self.__theta
# Line 2352  class TransportPDE(LinearProblem): Line 2585  class TransportPDE(LinearProblem):
2585    
2586     def checkSymmetry(self,verbose=True):     def checkSymmetry(self,verbose=True):
2587        """        """
2588        test the transport problem for symmetry.        Tests the transport problem for symmetry.
2589    
2590        @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.        @param verbose: if set to True or not present a report on coefficients
2591                          which break the symmetry is printed.
2592        @type verbose: C{bool}        @type verbose: C{bool}
2593        @return:  True if the PDE is symmetric.        @return:  True if the PDE is symmetric
2594        @rtype: C{bool}        @rtype: C{bool}
2595        @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.        @note: This is a very expensive operation. It should be used for
2596                 degugging only! The symmetry flag is not altered.
2597        """        """
2598        out=True        out=True
2599        out=out and self.checkSymmetricTensor("M", verbose)        out=out and self.checkSymmetricTensor("M", verbose)
# Line 2377  class TransportPDE(LinearProblem): Line 2612  class TransportPDE(LinearProblem):
2612        out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)        out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
2613        return out        return out
2614    
2615     def setValue(self,**coefficients):     def setValue(self,**coefficients):
2616        """        """
2617        sets new values to coefficients        Sets new values to coefficients.
2618    
2619        @param coefficients: new values assigned to coefficients        @param coefficients: new values assigned to coefficients
2620        @keyword M: value for coefficient M.        @keyword M: value for coefficient C{M}
2621        @type M: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @type M: any type that can be cast to a L{Data<escript.Data>} object on
2622        @keyword M_reduced: value for coefficient M_reduced.                 L{Function<escript.Function>}
2623        @type M_reduced: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.ReducedFunction>}.        @keyword M_reduced: value for coefficient C{M_reduced}
2624        @keyword A: value for coefficient A.        @type M_reduced: any type that can be cast to a L{Data<escript.Data>}
2625        @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.                         object on L{Function<escript.ReducedFunction>}
2626        @keyword A_reduced: value for coefficient A_reduced.        @keyword A: value for coefficient C{A}
2627        @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.        @type A: any type that can be cast to a L{Data<escript.Data>} object on
2628        @keyword B: value for coefficient B                 L{Function<escript.Function>}
2629        @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @keyword A_reduced: value for coefficient C{A_reduced}
2630        @keyword B_reduced: value for coefficient B_reduced        @type A_reduced: any type that can be cast to a L{Data<escript.Data>}
2631        @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.                         object on L{ReducedFunction<escript.ReducedFunction>}
2632        @keyword C: value for coefficient C        @keyword B: value for coefficient C{B}
2633        @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @type B: any type that can be cast to a L{Data<escript.Data>} object on
2634        @keyword C_reduced: value for coefficient C_reduced                 L{Function<escript.Function>}
2635        @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.        @keyword B_reduced: value for coefficient C{B_reduced}
2636        @keyword D: value for coefficient D        @type B_reduced: any type that can be cast to a L{Data<escript.Data>}
2637        @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.                         object on L{ReducedFunction<escript.ReducedFunction>}
2638        @keyword D_reduced: value for coefficient D_reduced        @keyword C: value for coefficient C{C}
2639        @type D_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.        @type C: any type that can be cast to a L{Data<escript.Data>} object on
2640        @keyword X: value for coefficient X                 L{Function<escript.Function>}
2641        @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @keyword C_reduced: value for coefficient C{C_reduced}
2642        @keyword X_reduced: value for coefficient X_reduced        @type C_reduced: any type that can be cast to a L{Data<escript.Data>}
2643        @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.                         object on L{ReducedFunction<escript.ReducedFunction>}
2644        @keyword Y: value for coefficient Y        @keyword D: value for coefficient C{D}
2645        @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @type D: any type that can be cast to a L{Data<escript.Data>} object on
2646        @keyword Y_reduced: value for coefficient Y_reduced                 L{Function<escript.Function>}
2647        @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.        @keyword D_reduced: value for coefficient C{D_reduced}
2648        @keyword m: value for coefficient m        @type D_reduced: any type that can be cast to a L{Data<escript.Data>}
2649        @type m: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.                         object on L{ReducedFunction<escript.ReducedFunction>}
2650        @keyword m_reduced: value for coefficient m_reduced        @keyword X: value for coefficient C{X}
2651        @type m_reduced: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.        @type X: any type that can be cast to a L{Data<escript.Data>} object on
2652        @keyword d: value for coefficient d                 L{Function<escript.Function>}
2653        @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.        @keyword X_reduced: value for coefficient C{X_reduced}
2654        @keyword d_reduced: value for coefficient d_reduced        @type X_reduced: any type that can be cast to a L{Data<escript.Data>}
2655        @type d_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.                         object on L{ReducedFunction<escript.ReducedFunction>}
2656        @keyword y: value for coefficient y        @keyword Y: value for coefficient C{Y}
2657        @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.        @type Y: any type that can be cast to a L{Data<escript.Data>} object on
2658        @keyword d_contact: value for coefficient d_contact                 L{Function<escript.Function>}
2659        @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.        @keyword Y_reduced: value for coefficient C{Y_reduced}
2660        @keyword d_contact_reduced: value for coefficient d_contact_reduced        @type Y_reduced: any type that can be cast to a L{Data<escript.Data>}
2661        @type d_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>} or  L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}.                         object on L{ReducedFunction<escript.Function>}
2662        @keyword y_contact: value for coefficient y_contact        @keyword m: value for coefficient C{m}
2663        @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.        @type m: any type that can be cast to a L{Data<escript.Data>} object on
2664        @keyword y_contact_reduced: value for coefficient y_contact_reduced                 L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2665        @type y_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}.        @keyword m_reduced: value for coefficient C{m_reduced}
2666          @type m_reduced: any type that can be cast to a L{Data<escript.Data>}
2667                           object on L{FunctionOnBoundary<escript.ReducedFunctionOnBoundary>}
2668          @keyword d: value for coefficient C{d}
2669          @type d: any type that can be cast to a L{Data<escript.Data>} object on
2670                   L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2671          @keyword d_reduced: value for coefficient C{d_reduced}
2672          @type d_reduced: any type that can be cast to a L{Data<escript.Data>}
2673                           object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}
2674          @keyword y: value for coefficient C{y}
2675          @type y: any type that can be cast to a L{Data<escript.Data>} object on
2676                   L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2677          @keyword d_contact: value for coefficient C{d_contact}
2678          @type d_contact: any type that can be cast to a L{Data<escript.Data>}
2679                           object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
2680          @keyword d_contact_reduced: value for coefficient C{d_contact_reduced}
2681          @type d_contact_reduced: any type that can be cast to a L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}
2682          @keyword y_contact: value for coefficient C{y_contact}
2683          @type y_contact: any type that can be cast to a L{Data<escript.Data>}
2684                           object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
2685          @keyword y_contact_reduced: value for coefficient C{y_contact_reduced}
2686          @type y_contact_reduced: any type that can be cast to a L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}
2687        @keyword r: values prescribed to the solution at the locations of constraints        @keyword r: values prescribed to the solution at the locations of constraints
2688        @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}        @type r: any type that can be cast to a L{Data<escript.Data>} object on
2689                 depending of reduced order is used for the solution.                 L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2690        @keyword q: mask for location of constraints                 depending on whether reduced order is used for the solution
2691        @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}        @keyword q: mask for the location of constraints
2692                 depending of reduced order is used for the representation of the equation.        @type q: any type that can be cast to a L{Data<escript.Data>} object on
2693        @raise IllegalCoefficient: if an unknown coefficient keyword is used.                 L{Solution<escript.Solution>} or
2694                   L{ReducedSolution<escript.ReducedSolution>} depending on whether
2695                   reduced order is used for the representation of the equation
2696          @raise IllegalCoefficient: if an unknown coefficient keyword is used
2697        """        """
2698        super(TransportPDE,self).setValue(**coefficients)        super(TransportPDE,self).setValue(**coefficients)
2699    
2700       def createOperator(self):
    def createOperator(self):  
2701         """         """
2702         returns an instance of a new transport operator         Returns an instance of a new transport operator.
2703         """         """
2704         self.trace("New Transport problem is allocated.")         self.trace("New Transport problem is allocated.")
2705         return self.getDomain().newTransportProblem( \         return self.getDomain().newTransportProblem( \
# Line 2452  class TransportPDE(LinearProblem): Line 2710  class TransportPDE(LinearProblem):
2710    
2711     def setInitialSolution(self,u):     def setInitialSolution(self,u):
2712         """         """
2713         sets the initial solution         Sets the initial solution.
2714        
2715         @param u: new initial solution         @param u: new initial solution
2716         @type u: any object that can ne interpolated to a L{Data<escript.Data>} object on on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}.         @type u: any object that can be interpolated to a L{Data<escript.Data>}
2717         @note: C{u} must be non negative.                  object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2718           @note: C{u} must be non-negative
2719         """         """
2720         u2=util.interpolate(u,self.getFunctionSpaceForSolution())         u2=util.interpolate(u,self.getFunctionSpaceForSolution())
2721         if self.getNumSolutions() == 1:         if self.getNumSolutions() == 1:
# Line 2466  class TransportPDE(LinearProblem): Line 2725  class TransportPDE(LinearProblem):
2725            if u2.getShape()!=(self.getNumSolutions(),):            if u2.getShape()!=(self.getNumSolutions(),):
2726                raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)                raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)
2727         self.getOperator().setInitialValue(u2)         self.getOperator().setInitialValue(u2)
2728    
2729     def getRequiredSystemType(self):     def getRequiredSystemType(self):
2730        """        """
2731        returns the system type which needs to be used by the current set up.        Returns the system type which needs to be used by the current set up.
2732    
2733        @return: a code to indicate the type of transport problem scheme used.        @return: a code to indicate the type of transport problem scheme used
2734        @rtype: C{float}        @rtype: C{float}
2735        """        """
2736        return self.getDomain().getTransportTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())        return self.getDomain().getTransportTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())
2737    
2738     def getUnlimitedTimeStepSize(self):     def getUnlimitedTimeStepSize(self):
2739        """        """
2740        returns the value returned by the C{getSafeTimeStepSize} method to indicate no limit on the save time step size.        Returns the value returned by the C{getSafeTimeStepSize} method to
2741          indicate no limit on the safe time step size.
2742    
2743         @return: the value used to indicate that no limit is set to the time step size         @return: the value used to indicate that no limit is set to the time
2744                    step size
2745         @rtype: C{float}         @rtype: C{float}
2746         @note: Typically the biggest positive float is returned.         @note: Typically the biggest positive float is returned
2747        """        """
2748        return self.getOperator().getUnlimitedTimeStepSize()        return self.getOperator().getUnlimitedTimeStepSize()
2749    
2750     def getSafeTimeStepSize(self):     def getSafeTimeStepSize(self):
2751         """         """
2752         returns a safe time step size to do the next time step.         Returns a safe time step size to do the next time step.
2753    
2754         @return: safe time step size         @return: safe time step size
2755         @rtype: C{float}         @rtype: C{float}
2756         @note: If not C{getSafeTimeStepSize()} < C{getUnlimitedTimeStepSize()} any time step size can be used.         @note: If not C{getSafeTimeStepSize()} < C{getUnlimitedTimeStepSize()}
2757                  any time step size can be used.
2758         """         """
2759         return self.getOperator().getSafeTimeStepSize()         return self.getOperator().getSafeTimeStepSize()
2760     #====================================================================     #====================================================================
2761     def getSolution(self,dt,**options):     def getSolution(self,dt,**options):
2762         """         """
2763         returns the solution of the problem         Returns the solution of the problem.
2764    
2765         @return: the solution         @return: the solution
2766         @rtype: L{Data<escript.Data>}         @rtype: L{Data<escript.Data>}
2767    
        @note: This method is overwritten when implementing a particular linear problem.  
2768         """         """
2769         if dt<=0:         if dt<=0:
2770             raise ValueError,"step size needs to be positive."             raise ValueError,"step size needs to be positive."
# Line 2510  class TransportPDE(LinearProblem): Line 2773  class TransportPDE(LinearProblem):
2773         self.validSolution()         self.validSolution()
2774         return self.getCurrentSolution()         return self.getCurrentSolution()
2775    
2776     def getSystem(self):     def getSystem(self):
2777         """         """
2778         return the operator and right hand side of the PDE         Returns the operator and right hand side of the PDE.
2779    
2780         @return: the discrete version of the PDE         @return: the discrete version of the PDE
2781         @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.         @rtype: C{tuple} of L{Operator,<escript.Operator>} and
2782                   L{Data<escript.Data>}
2783    
        @note: This method is overwritten when implementing a particular linear problem.  
2784         """         """
2785         if not self.isOperatorValid() or not self.isRightHandSideValid():         if not self.isOperatorValid() or not self.isRightHandSideValid():
2786            self.resetRightHandSide()            self.resetRightHandSide()
# Line 2557  class TransportPDE(LinearProblem): Line 2820  class TransportPDE(LinearProblem):
2820            self.validOperator()            self.validOperator()
2821            self.validRightHandSide()            self.validRightHandSide()
2822         return (self.getCurrentOperator(), self.getCurrentRightHandSide())         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
2823    

Legend:
Removed from v.2168  
changed lines
  Added in v.2169

  ViewVC Help
Powered by ViewVC 1.1.26