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

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

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

trunk/esys2/escript/py_src/linearPDEs.py revision 108 by jgs, Thu Jan 27 06:21:59 2005 UTC trunk/escript/py_src/linearPDEs.py revision 969 by ksteube, Tue Feb 13 23:02:23 2007 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
   
 ## @file linearPDEs.py  
   
2  """  """
3  @brief Functions and classes for linear PDEs  The module provides an interface to define and solve linear partial
4    differential equations (PDEs) within L{escript}. L{linearPDEs} does not provide any
5    solver capabilities in itself but hands the PDE over to
6    the PDE solver library defined through the L{Domain<escript.Domain>} of the PDE.
7    The general interface is provided through the L{LinearPDE} class. The
8    L{AdvectivePDE} which is derived from the L{LinearPDE} class
9    provides an interface to PDE dominated by its advective terms. The L{Poisson},
10    L{Helmholtz}, L{LameEquation}, L{AdvectivePDE}
11    classs which are also derived form the L{LinearPDE} class should be used
12    to define of solve these sepecial PDEs.
13    
14    @var __author__: name of author
15    @var __copyright__: copyrights
16    @var __license__: licence agreement
17    @var __url__: url entry point on documentation
18    @var __version__: version
19    @var __date__: date of the version
20  """  """
21    
22  import escript  import escript
23  import util  import util
24  import numarray  import numarray
25    
26    __author__="Lutz Gross, l.gross@uq.edu.au"
27    __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
28                        http://www.access.edu.au
29                    Primary Business: Queensland, Australia"""
30    __license__="""Licensed under the Open Software License version 3.0
31                 http://www.opensource.org/licenses/osl-3.0.php"""
32    __url__="http://www.iservo.edu.au/esys"
33    __version__="$Revision$"
34    __date__="$Date$"
35    
36    
37    class IllegalCoefficient(ValueError):
38       """
39       raised if an illegal coefficient of the general ar particular PDE is requested.
40       """
41    
42  def _CompTuple2(t1,t2):  class IllegalCoefficientValue(ValueError):
43       """
44       raised if an incorrect value for a coefficient is used.
45     """     """
    @brief  
46    
47     @param t1  class UndefinedPDEError(ValueError):
48     @param t2     """
49       raised if a PDE is not fully defined yet.
50     """     """
    dif=t1[0]+t1[1]-(t2[0]+t2[1])  
    if dif<0: return 1  
    elif dif>0: return -1  
    else: return 0  
51    
52  class PDECoefficient:  class PDECoefficient(object):
53      """      """
54      @brief      A class for describing a PDE coefficient
55    
56        @cvar INTERIOR: indicator that coefficient is defined on the interior of the PDE domain
57        @cvar BOUNDARY: indicator that coefficient is defined on the boundary of the PDE domain
58        @cvar CONTACT: indicator that coefficient is defined on the contact region within the PDE domain
59        @cvar SOLUTION: indicator that coefficient is defined trough a solution of the PDE
60        @cvar REDUCED: indicator that coefficient is defined trough a reduced solution of the PDE
61        @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is defined by the number PDE equations
62        @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is defined by the number PDE solutions
63        @cvar BY_DIM: indicator that the dimension of the coefficient shape is defined by the spatial dimension
64        @cvar OPERATOR: indicator that the the coefficient alters the operator of the PDE
65        @cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right hand side of the PDE
66        @cvar BOTH: indicator that the the coefficient alters the operator as well as the right hand side of the PDE
67    
68      """      """
     # identifier for location of Data objects defining COEFFICIENTS  
69      INTERIOR=0      INTERIOR=0
70      BOUNDARY=1      BOUNDARY=1
71      CONTACT=2      CONTACT=2
72      CONTINUOUS=3      SOLUTION=3
73      # identifier in the pattern of COEFFICIENTS:      REDUCED=4
74      # the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the      BY_EQUATION=5
75      # number of unknowns.      BY_SOLUTION=6
76      EQUATION=3      BY_DIM=7
77      SOLUTION=4      OPERATOR=10
78      DIM=5      RIGHTHANDSIDE=11
79      # indicator for what is altered if the coefficient is altered:      BOTH=12
80      OPERATOR=5  
     RIGHTHANDSIDE=6  
     BOTH=7  
81      def __init__(self,where,pattern,altering):      def __init__(self,where,pattern,altering):
82         """         """
83         @brief Initialise a PDE Coefficient type         Initialise a PDE Coefficient type
84    
85           @param where: describes where the coefficient lives
86           @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION}, L{REDUCED}
87           @param pattern: describes the shape of the coefficient and how the shape is build for a given
88                  spatial dimension and numbers of equation and solution in then PDE. For instance,
89                  (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) descrbes a rank 3 coefficient which
90                  is instanciated as shape (3,2,2) in case of a three equations and two solution components
91                  on a 2-dimensional domain. In the case of single equation and a single solution component
92                  the shape compoments marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped. In this case
93                  the example would be read as (2,).
94           @type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM}
95           @param altering: indicates what part of the PDE is altered if the coefficiennt is altered
96           @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}
97    
98         """         """
99           super(PDECoefficient, self).__init__()
100         self.what=where         self.what=where
101         self.pattern=pattern         self.pattern=pattern
102         self.altering=altering         self.altering=altering
# Line 53  class PDECoefficient: Line 104  class PDECoefficient:
104    
105      def resetValue(self):      def resetValue(self):
106         """         """
107         @brief resets coefficient value to default         resets coefficient value to default
108         """         """
109         self.value=escript.Data()         self.value=escript.Data()
110    
111      def getFunctionSpace(self,domain):      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
112         """         """
113         @brief defines the FunctionSpace of the coefficient on the domain         defines the L{FunctionSpace<escript.FunctionSpace>} of the coefficient
114    
115         @param domain         @param domain: domain on which the PDE uses the coefficient
116         """         @type domain: L{Domain<escript.Domain>}
117         if self.what==self.INTERIOR: return escript.Function(domain)         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
118         elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)         @type reducedEquationOrder: C{bool}
119         elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
120         elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)         @type reducedSolutionOrder: C{bool}
121           @return:  L{FunctionSpace<escript.FunctionSpace>} of the coefficient
122           @rtype:  L{FunctionSpace<escript.FunctionSpace>}
123           """
124           if self.what==self.INTERIOR:
125                return escript.Function(domain)
126           elif self.what==self.BOUNDARY:
127                return escript.FunctionOnBoundary(domain)
128           elif self.what==self.CONTACT:
129                return escript.FunctionOnContactZero(domain)
130           elif self.what==self.SOLUTION:
131                if reducedEquationOrder and reducedSolutionOrder:
132                    return escript.ReducedSolution(domain)
133                else:
134                    return escript.Solution(domain)
135           elif self.what==self.REDUCED:
136                return escript.ReducedSolution(domain)
137    
138      def getValue(self):      def getValue(self):
139         """         """
140         @brief returns the value of the coefficient:         returns the value of the coefficient
141    
142           @return:  value of the coefficient
143           @rtype:  L{Data<escript.Data>}
144         """         """
145         return self.value         return self.value
146        
147      def setValue(self,newValue):      def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
        """  
        @brief set the value of the coefficient to new value  
148         """         """
149           set the value of the coefficient to a new value
150    
151           @param domain: domain on which the PDE uses the coefficient
152           @type domain: L{Domain<escript.Domain>}
153           @param numEquations: number of equations of the PDE
154           @type numEquations: C{int}
155           @param numSolutions: number of components of the PDE solution
156           @type numSolutions: C{int}
157           @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
158           @type reducedEquationOrder: C{bool}
159           @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
160           @type reducedSolutionOrder: C{bool}
161           @param newValue: number of components of the PDE solution
162           @type newValue: any object that can be converted into a L{Data<escript.Data>} object with the appropriate shape and L{FunctionSpace<escript.FunctionSpace>}
163           @raise IllegalCoefficientValue: if the shape of the assigned value does not match the shape of the coefficient
164           """
165           if newValue==None:
166               newValue=escript.Data()
167           elif isinstance(newValue,escript.Data):
168               if not newValue.isEmpty():
169                  try:
170                     newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
171                  except:
172                     raise IllegalCoefficientValue,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)
173           else:
174               newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
175           if not newValue.isEmpty():
176               if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
177                   raise IllegalCoefficientValue,"Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())
178         self.value=newValue         self.value=newValue
179        
180      def isAlteringOperator(self):      def isAlteringOperator(self):
181          """          """
182      @brief return true if the operator of the PDE is changed when the coefficient is changed          checks if the coefficient alters the operator of the PDE
183    
184            @return:  True if the operator of the PDE is changed when the coefficient is changed
185            @rtype:  C{bool}
186      """      """
187          if self.altering==self.OPERATOR or self.altering==self.BOTH:          if self.altering==self.OPERATOR or self.altering==self.BOTH:
188              return not None              return not None
# Line 91  class PDECoefficient: Line 191  class PDECoefficient:
191    
192      def isAlteringRightHandSide(self):      def isAlteringRightHandSide(self):
193          """          """
194      @brief return true if the right hand side of the PDE is changed when the coefficient is changed          checks if the coefficeint alters the right hand side of the PDE
195    
196        @rtype:  C{bool}
197            @return:  True if the right hand side of the PDE is changed when the coefficient is changed
198      """      """
199          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
200              return not None              return not None
201          else:          else:
202              return None              return None
203    
204      def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3):      def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
205         """         """
206         @brief tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim         tries to estimate the number of equations and number of solutions if the coefficient has the given shape
207    
208         @param shape         @param domain: domain on which the PDE uses the coefficient
209         @param dim         @type domain: L{Domain<escript.Domain>}
210           @param shape: suggested shape of the coefficient
211           @type shape: C{tuple} of C{int} values
212           @return: the number of equations and number of solutions of the PDE is the coefficient has shape s.
213                     If no appropriate numbers could be identified, C{None} is returned
214           @rtype: C{tuple} of two C{int} values or C{None}
215         """         """
216           dim=domain.getDim()
217         if len(shape)>0:         if len(shape)>0:
218             num=max(shape)+1             num=max(shape)+1
219         else:         else:
220             num=1             num=1
221         search=[]         search=[]
222         for u in range(num):         if self.definesNumEquation() and self.definesNumSolutions():
223            for e in range(num):            for u in range(num):
224               search.append((e,u))               for e in range(num):
225         search.sort(_CompTuple2)                  search.append((e,u))
226         for item in search:            search.sort(self.__CompTuple2)
227               s=self.buildShape(item[0],item[1],dim)            for item in search:
228                 s=self.getShape(domain,item[0],item[1])
229               if len(s)==0 and len(shape)==0:               if len(s)==0 and len(shape)==0:
230                   return (1,1)                   return (1,1)
231               else:               else:
232                   if s==shape: return item                   if s==shape: return item
233           elif self.definesNumEquation():
234              for e in range(num,0,-1):
235                 s=self.getShape(domain,e,0)
236                 if len(s)==0 and len(shape)==0:
237                     return (1,None)
238                 else:
239                     if s==shape: return (e,None)
240    
241           elif self.definesNumSolutions():
242              for u in range(num,0,-1):
243                 s=self.getShape(domain,0,u)
244                 if len(s)==0 and len(shape)==0:
245                     return (None,1)
246                 else:
247                     if s==shape: return (None,u)
248         return None         return None
249        def definesNumSolutions(self):
250           """
251           checks if the coefficient allows to estimate the number of solution components
252    
253      def buildShape(self,e=1,u=1,dim=3):         @return: True if the coefficient allows an estimate of the number of solution components
254          """         @rtype: C{bool}
255      @brief builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim         """
256           for i in self.pattern:
257                 if i==self.BY_SOLUTION: return True
258           return False
259    
260      @param e      def definesNumEquation(self):
261      @param u         """
262      @param dim         checks if the coefficient allows to estimate the number of equations
263      """  
264          s=()         @return: True if the coefficient allows an estimate of the number of equations
265          for i in self.pattern:         @rtype: C{bool}
266               if i==self.EQUATION:         """
267                  if e>1: s=s+(e,)         for i in self.pattern:
268               elif i==self.SOLUTION:               if i==self.BY_EQUATION: return True
269                  if u>1: s=s+(u,)         return False
270    
271        def __CompTuple2(self,t1,t2):
272          """
273          Compare two tuples of possible number of equations and number of solutions
274    
275          @param t1: The first tuple
276          @param t2: The second tuple
277    
278          """
279    
280          dif=t1[0]+t1[1]-(t2[0]+t2[1])
281          if dif<0: return 1
282          elif dif>0: return -1
283          else: return 0
284    
285        def getShape(self,domain,numEquations=1,numSolutions=1):
286           """
287           builds the required shape of the coefficient
288    
289           @param domain: domain on which the PDE uses the coefficient
290           @type domain: L{Domain<escript.Domain>}
291           @param numEquations: number of equations of the PDE
292           @type numEquations: C{int}
293           @param numSolutions: number of components of the PDE solution
294           @type numSolutions: C{int}
295           @return: shape of the coefficient
296           @rtype: C{tuple} of C{int} values
297           """
298           dim=domain.getDim()
299           s=()
300           for i in self.pattern:
301                 if i==self.BY_EQUATION:
302                    if numEquations>1: s=s+(numEquations,)
303                 elif i==self.BY_SOLUTION:
304                    if numSolutions>1: s=s+(numSolutions,)
305               else:               else:
306                  s=s+(dim,)                  s=s+(dim,)
307          return s         return s
308    
309  class LinearPDE:  class LinearPDE(object):
310     """     """
311     @brief Class to handel a linear PDE     This class is used to define a general linear, steady, second order PDE
312         for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
    class to define a linear PDE of the form  
313    
314       -(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i     For a single PDE with a solution with a single component the linear PDE is defined in the following form:
315    
316       with boundary conditons:     M{-grad(A[j,l]*grad(u)[l]+B[j]u)[j]+C[l]*grad(u)[l]+D*u =-grad(X)[j,j]+Y}
317    
318          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i     where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,
319       ie. summation over indexes appearing twice in a term of a sum is performed, is used.
320       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
321       L{Function<escript.Function>} on the PDE or objects that can be converted into such L{Data<escript.Data>} objects.
322       M{A} is a rank two, M{B}, M{C} and M{X} are rank one and M{D} and M{Y} are scalar.
323    
324      and contact conditions     The following natural boundary conditions are considered:
325    
326          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i     M{n[j]*(A[i,j]*grad(u)[l]+B[j]*u)+d*u=n[j]*X[j]+y}
327    
328      and constraints:     where M{n} is the outer normal field calculated by L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
329       Notice that the coefficients M{A}, M{B} and M{X} are defined in the PDE. The coefficients M{d} and M{y} are
330       each a scalar in the L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
331    
          u_i=r_i where q_i>0  
332    
333     """     Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
    TOL=1.e-13  
    DEFAULT_METHOD=util.DEFAULT_METHOD  
    DIRECT=util.DIRECT  
    CHOLEVSKY=util.CHOLEVSKY  
    PCG=util.PCG  
    CR=util.CR  
    CGS=util.CGS  
    BICGSTAB=util.BICGSTAB  
    SSOR=util.SSOR  
    GMRES=util.GMRES  
    PRES20=util.PRES20  
334    
335     def __init__(self,domain,numEquations=0,numSolutions=0):     M{u=r}  where M{q>0}
      """  
      @brief initializes a new linear PDE.  
336    
337       @param args     M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
338       """     The constraints override any other condition set by the PDE or the boundary condition.
339       # COEFFICIENTS can be overwritten by subclasses:  
340       self.COEFFICIENTS={     The PDE is symmetrical if
341         "A"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),  
342         "B"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),     M{A[i,j]=A[j,i]}  and M{B[j]=C[j]}
343         "C"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),  
344         "D"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),     For a system of PDEs and a solution with several components the PDE has the form
345         "X"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),  
346         "Y"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),     M{-grad(A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k])[j]+C[i,k,l]*grad(u[k])[l]+D[i,k]*u[k] =-grad(X[i,j])[j]+Y[i] }
347         "d"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),  
348         "y"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),     M{A} is a ramk four, M{B} and M{C} are each a rank three, M{D} and M{X} are each a rank two and M{Y} is a rank one.
349         "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),     The natural boundary conditions take the form:
350         "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
351         "r"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),     M{n[j]*(A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k])+d[i,k]*u[k]=n[j]*X[i,j]+y[i]}
352         "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)}  
353    
354       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
355    
356    
357       M{u[i]=r[i]}  where  M{q[i]>0}
358    
359       M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
360    
361       The system of PDEs is symmetrical if
362    
363            - M{A[i,j,k,l]=A[k,l,i,j]}
364            - M{B[i,j,k]=C[k,i,j]}
365            - M{D[i,k]=D[i,k]}
366            - M{d[i,k]=d[k,i]}
367    
368       L{LinearPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
369       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
370       defined as
371    
372       M{J[i,j]=A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k]-X[i,j]}
373    
374       For the case of single solution component and single PDE M{J} is defined
375    
376       M{J_{j}=A[i,j]*grad(u)[j]+B[i]*u-X[i]}
377    
378       In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
379       calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of PDEs
380       the contact condition takes the form
381    
382       M{n[j]*J0[i,j]=n[j]*J1[i,j]=y_contact[i]- d_contact[i,k]*jump(u)[k]}
383    
384       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
385       of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calcualted by
386       L{jump<util.jump>}.
387       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>}.
388       In case of a single PDE and a single component solution the contact condition takes the form
389    
390       M{n[j]*J0_{j}=n[j]*J1_{j}=y_contact-d_contact*jump(u)}
391    
392       In this case the the coefficient M{d_contact} and M{y_contact} are eaach scalar
393       both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
394    
395       @cvar DEFAULT: The default method used to solve the system of linear equations
396       @cvar DIRECT: The direct solver based on LDU factorization
397       @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
398       @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
399       @cvar CR: The conjugate residual method
400       @cvar CGS: The conjugate gardient square method
401       @cvar BICGSTAB: The stabilized BiConjugate Gradient method.
402       @cvar SSOR: The symmetric overrealaxtion method
403       @cvar ILU0: The incomplete LU factorization preconditioner  with no fill in
404       @cvar ILUT: The incomplete LU factorization preconditioner with will in
405       @cvar JACOBI: The Jacobi preconditioner
406       @cvar GMRES: The Gram-Schmidt minimum residual method
407       @cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals
408       @cvar LUMPING: Matrix lumping.
409       @cvar NO_REORDERING: No matrix reordering allowed
410       @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
411       @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
412       @cvar PASO: PASO solver package
413       @cvar SCSL: SGI SCSL solver library
414       @cvar MKL: Intel's MKL solver library
415       @cvar UMFPACK: the UMFPACK library
416       @cvar TRILINOS: the TRILINOS parallel solver class library from Sandia Natl Labs
417       @cvar ITERATIVE: The default iterative solver
418       @cvar AMG: algebraic multi grid
419       @cvar RILU: recursive ILU
420    
421       """
422       DEFAULT= 0
423       DIRECT= 1
424       CHOLEVSKY= 2
425       PCG= 3
426       CR= 4
427       CGS= 5
428       BICGSTAB= 6
429       SSOR= 7
430       ILU0= 8
431       ILUT= 9
432       JACOBI= 10
433       GMRES= 11
434       PRES20= 12
435       LUMPING= 13
436       NO_REORDERING= 17
437       MINIMUM_FILL_IN= 18
438       NESTED_DISSECTION= 19
439       SCSL= 14
440       MKL= 15
441       UMFPACK= 16
442       ITERATIVE= 20
443       PASO= 21
444       AMG= 22
445       RILU = 23
446       TRILINOS = 24
447    
448       SMALL_TOLERANCE=1.e-13
449       __PACKAGE_KEY="package"
450       __METHOD_KEY="method"
451       __SYMMETRY_KEY="symmetric"
452       __TOLERANCE_KEY="tolerance"
453       __PRECONDITIONER_KEY="preconditioner"
454    
455    
456       def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
457         """
458         initializes a new linear PDE
459    
460         @param domain: domain of the PDE
461         @type domain: L{Domain<escript.Domain>}
462         @param numEquations: number of equations. If numEquations==None the number of equations
463                              is exracted from the PDE coefficients.
464         @param numSolutions: number of solution components. If  numSolutions==None the number of solution components
465                              is exracted from the PDE coefficients.
466         @param debug: if True debug informations are printed.
467    
468         """
469         super(LinearPDE, self).__init__()
470         #
471         #   the coefficients of the general PDE:
472         #
473         self.__COEFFICIENTS_OF_GENEARL_PDE={
474           "A"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
475           "B"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
476           "C"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
477           "D"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
478           "X"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
479           "Y"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
480           "d"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
481           "y"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
482           "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
483           "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
484           "r"         : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.RIGHTHANDSIDE),
485           "q"         : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.BOTH)}
486    
487         # COEFFICIENTS can be overwritten by subclasses:
488         self.COEFFICIENTS=self.__COEFFICIENTS_OF_GENEARL_PDE
489         self.__altered_coefficients=False
490       # initialize attributes       # initialize attributes
491       self.__debug=None       self.__debug=debug
492       self.__domain=domain       self.__domain=domain
493       self.__numEquations=numEquations       self.__numEquations=numEquations
494       self.__numSolutions=numSolutions       self.__numSolutions=numSolutions
495       self.cleanCoefficients()       self.__resetSystem()
   
      self.__operator=escript.Operator()  
      self.__operator_isValid=False  
      self.__righthandside=escript.Data()  
      self.__righthandside_isValid=False  
      self.__solution=escript.Data()  
      self.__solution_isValid=False  
496    
497       # set some default values:       # set some default values:
498       self.__homogeneous_constraint=True       self.__reduce_equation_order=False
499       self.__row_function_space=escript.Solution(self.__domain)       self.__reduce_solution_order=False
      self.__column_function_space=escript.Solution(self.__domain)  
500       self.__tolerance=1.e-8       self.__tolerance=1.e-8
501       self.__solver_method=util.DEFAULT_METHOD       self.__solver_method=self.DEFAULT
502       self.__matrix_type=self.__domain.getSystemMatrixTypeId(util.DEFAULT_METHOD,False)       self.__solver_package=self.DEFAULT
503         self.__preconditioner=self.DEFAULT
504         self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)
505       self.__sym=False       self.__sym=False
      self.__lumping=False  
506    
507     def createCoefficient(self, name):       self.resetCoefficients()
508         self.trace("PDE Coeffients are %s"%str(self.COEFFICIENTS.keys()))
509       # =============================================================================
510       #    general stuff:
511       # =============================================================================
512       def __str__(self):
513         """
514         returns string representation of the PDE
515    
516         @return: a simple representation of the PDE
517         @rtype: C{str}
518         """
519         return "<LinearPDE %d>"%id(self)
520       # =============================================================================
521       #    debug :
522       # =============================================================================
523       def setDebugOn(self):
524       """       """
525       @brief create a data object corresponding to coefficient name       switches on debugging
      @param name  
526       """       """
527       return escript.Data(shape = getShapeOfCoefficient(name), \       self.__debug=not None
                          what = getFunctionSpaceOfCoefficient(name))  
528    
529     def __del__(self):     def setDebugOff(self):
      pass  
   
    def getCoefficient(self,name):  
530       """       """
531       @brief return the value of the parameter name       switches off debugging
532         """
533         self.__debug=None
534    
535       @param name     def trace(self,text):
536       """       """
537       return self.COEFFICIENTS[name].getValue()       print the text message if debugging is swiched on.
538         @param text: message
539         @type text: C{string}
540         """
541         if self.__debug: print "%s: %s"%(str(self),text)
542    
543     def getCoefficientOfPDE(self,name):     # =============================================================================
544       # some service functions:
545       # =============================================================================
546       def getDomain(self):
547       """       """
548       @brief return the value of the coefficient name of the general PDE. This method is called by the assembling routine       returns the domain of the PDE
549              it can be overwritten to map coefficients of a particualr PDE to the general PDE.  
550       @param name       @return: the domain of the PDE
551         @rtype: L{Domain<escript.Domain>}
552       """       """
553       return self.getCoefficient(name)       return self.__domain
554    
555     def hasCoefficient(self,name):     def getDim(self):
556        """       """
557        @brief return true if name is the name of a coefficient       returns the spatial dimension of the PDE
558    
559        @param name       @return: the spatial dimension of the PDE domain
560        """       @rtype: C{int}
561        return self.COEFFICIENTS.has_key(name)       """
562         return self.getDomain().getDim()
563    
564     def getFunctionSpaceForEquation(self):     def getNumEquations(self):
565       """       """
566       @brief return true if the test functions should use reduced order       returns the number of equations
567    
568         @return: the number of equations
569         @rtype: C{int}
570         @raise UndefinedPDEError: if the number of equations is not be specified yet.
571       """       """
572       return self.__row_function_space       if self.__numEquations==None:
573             raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
574         else:
575             return self.__numEquations
576    
577     def getFunctionSpaceForSolution(self):     def getNumSolutions(self):
578       """       """
579       @brief return true if the interpolation of the solution should use reduced order       returns the number of unknowns
580    
581         @return: the number of unknowns
582         @rtype: C{int}
583         @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
584       """       """
585       return self.__column_function_space       if self.__numSolutions==None:
586            raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
587         else:
588            return self.__numSolutions
589    
590     def setValue(self,**coefficients):     def reduceEquationOrder(self):
591        """       """
592        @brief sets new values to coefficients       return status for order reduction for equation
593    
594        @param coefficients       @return: return True is reduced interpolation order is used for the represenation of the equation
595        """       @rtype: L{bool}
596        self._setValue(**coefficients)       """
597               return self.__reduce_equation_order
598    
599     def cleanCoefficients(self):     def reduceSolutionOrder(self):
600       """       """
601       @brief resets all coefficients to default values.       return status for order reduction for the solution
602    
603         @return: return True is reduced interpolation order is used for the represenation of the solution
604         @rtype: L{bool}
605       """       """
606       for i in self.COEFFICIENTS.iterkeys():       return self.__reduce_solution_order
          self.COEFFICIENTS[i].resetValue()  
607    
608     def createNewCoefficient(self,name):     def getFunctionSpaceForEquation(self):
609       """       """
610       @brief returns a new coefficient appropriate for coefficient name:       returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation
611    
612         @return: representation space of equation
613         @rtype: L{FunctionSpace<escript.FunctionSpace>}
614       """       """
615       return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))       if self.reduceEquationOrder():
616                   return escript.ReducedSolution(self.getDomain())
617         else:
618             return escript.Solution(self.getDomain())
619    
620     def getShapeOfCoefficient(self,name):     def getFunctionSpaceForSolution(self):
621       """       """
622       @brief return the shape of the coefficient name       returns the L{FunctionSpace<escript.FunctionSpace>} used to represent the solution
623    
624       @param name       @return: representation space of solution
625         @rtype: L{FunctionSpace<escript.FunctionSpace>}
626       """       """
627       if self.hasCoefficient(name):       if self.reduceSolutionOrder():
628          return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())           return escript.ReducedSolution(self.getDomain())
629       else:       else:
630          raise ValueError,"Solution coefficient %s requested"%name           return escript.Solution(self.getDomain())
631    
632     def getFunctionSpaceForCoefficient(self,name):  
633       def getOperator(self):
634       """       """
635       @brief return the atoms of the coefficient name       provides access to the operator of the PDE
636    
637       @param name       @return: the operator of the PDE
638         @rtype: L{Operator<escript.Operator>}
639       """       """
640       if self.hasCoefficient(name):       m=self.getSystem()[0]
641          return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())       if self.isUsingLumping():
642             return self.copyConstraint(1./m)
643       else:       else:
644          raise ValueError,"Solution coefficient %s requested"%name           return m
645    
646     def alteredCoefficient(self,name):     def getRightHandSide(self):
647       """       """
648       @brief annonced that coefficient name has been changed       provides access to the right hand side of the PDE
649         @return: the right hand side of the PDE
650         @rtype: L{Data<escript.Data>}
651         """
652         r=self.getSystem()[1]
653         if self.isUsingLumping():
654             return self.copyConstraint(r)
655         else:
656             return r
657    
658       @param name     def applyOperator(self,u=None):
659       """       """
660       if self.hasCoefficient(name):       applies the operator of the PDE to a given u or the solution of PDE if u is not present.
661          if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator()  
662          if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()       @param u: argument of the operator. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
663                   the current solution is used.
664         @type u: L{Data<escript.Data>} or None
665         @return: image of u
666         @rtype: L{Data<escript.Data>}
667         """
668         if u==None:
669            return self.getOperator()*self.getSolution()
670       else:       else:
671          raise ValueError,"unknown coefficient %s requested"%name          return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
672    
673     # ===== debug ==============================================================     def getResidual(self,u=None):
674     def setDebugOn(self):       """
675         return the residual of u or the current solution if u is not present.
676    
677         @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
678                   the current solution is used.
679         @type u: L{Data<escript.Data>} or None
680         @return: residual of u
681         @rtype: L{Data<escript.Data>}
682         """
683         return self.applyOperator(u)-self.getRightHandSide()
684    
685       def checkSymmetry(self,verbose=True):
686          """
687          test the PDE for symmetry.
688    
689          @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
690          @type verbose: C{bool}
691          @return:  True if the PDE is symmetric.
692          @rtype: L{Data<escript.Data>}
693          @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
694          """
695          verbose=verbose or self.__debug
696          out=True
697          if self.getNumSolutions()!=self.getNumEquations():
698             if verbose: print "non-symmetric PDE because of different number of equations and solutions"
699             out=False
700          else:
701             A=self.getCoefficientOfGeneralPDE("A")
702             if not A.isEmpty():
703                tol=util.Lsup(A)*self.SMALL_TOLERANCE
704                if self.getNumSolutions()>1:
705                   for i in range(self.getNumEquations()):
706                      for j in range(self.getDim()):
707                         for k in range(self.getNumSolutions()):
708                            for l in range(self.getDim()):
709                                if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
710                                   if verbose: print "non-symmetric PDE because A[%d,%d,%d,%d]!=A[%d,%d,%d,%d]"%(i,j,k,l,k,l,i,j)
711                                   out=False
712                else:
713                   for j in range(self.getDim()):
714                      for l in range(self.getDim()):
715                         if util.Lsup(A[j,l]-A[l,j])>tol:
716                            if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)
717                            out=False
718             B=self.getCoefficientOfGeneralPDE("B")
719             C=self.getCoefficientOfGeneralPDE("C")
720             if B.isEmpty() and not C.isEmpty():
721                if verbose: print "non-symmetric PDE because B is not present but C is"
722                out=False
723             elif not B.isEmpty() and C.isEmpty():
724                if verbose: print "non-symmetric PDE because C is not present but B is"
725                out=False
726             elif not B.isEmpty() and not C.isEmpty():
727                tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
728                if self.getNumSolutions()>1:
729                   for i in range(self.getNumEquations()):
730                       for j in range(self.getDim()):
731                          for k in range(self.getNumSolutions()):
732                             if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
733                                  if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
734                                  out=False
735                else:
736                   for j in range(self.getDim()):
737                      if util.Lsup(B[j]-C[j])>tol:
738                         if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)
739                         out=False
740             if self.getNumSolutions()>1:
741               D=self.getCoefficientOfGeneralPDE("D")
742               if not D.isEmpty():
743                 tol=util.Lsup(D)*self.SMALL_TOLERANCE
744                 for i in range(self.getNumEquations()):
745                    for k in range(self.getNumSolutions()):
746                      if util.Lsup(D[i,k]-D[k,i])>tol:
747                          if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)
748                          out=False
749               d=self.getCoefficientOfGeneralPDE("d")
750               if not d.isEmpty():
751                 tol=util.Lsup(d)*self.SMALL_TOLERANCE
752                 for i in range(self.getNumEquations()):
753                    for k in range(self.getNumSolutions()):
754                      if util.Lsup(d[i,k]-d[k,i])>tol:
755                          if verbose: print "non-symmetric PDE because d[%d,%d]!=d[%d,%d]"%(i,k,k,i)
756                          out=False
757               d_contact=self.getCoefficientOfGeneralPDE("d_contact")
758               if not d_contact.isEmpty():
759                 tol=util.Lsup(d_contact)*self.SMALL_TOLERANCE
760                 for i in range(self.getNumEquations()):
761                    for k in range(self.getNumSolutions()):
762                      if util.Lsup(d_contact[i,k]-d_contact[k,i])>tol:
763                          if verbose: print "non-symmetric PDE because d_contact[%d,%d]!=d_contact[%d,%d]"%(i,k,k,i)
764                          out=False
765          return out
766    
767       def getSolution(self,**options):
768         """         """
769         @brief         returns the solution of the PDE. If the solution is not valid the PDE is solved.
770    
771           @return: the solution
772           @rtype: L{Data<escript.Data>}
773           @param options: solver options
774           @keyword verbose: True to get some information during PDE solution
775           @type verbose: C{bool}
776           @keyword reordering: reordering scheme to be used during elimination. Allowed values are
777                                L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}
778           @keyword iter_max: maximum number of iteration steps allowed.
779           @keyword drop_tolerance: threshold for drupping in L{ILUT}
780           @keyword drop_storage: maximum of allowed memory in L{ILUT}
781           @keyword truncation: maximum number of residuals in L{GMRES}
782           @keyword restart: restart cycle length in L{GMRES}
783         """         """
784         self.__debug=not None         if not self.__solution_isValid:
785              mat,f=self.getSystem()
786              if self.isUsingLumping():
787                 self.__solution=self.copyConstraint(f*mat)
788              else:
789                 options[self.__TOLERANCE_KEY]=self.getTolerance()
790                 options[self.__METHOD_KEY]=self.getSolverMethod()[0]
791                 options[self.__PRECONDITIONER_KEY]=self.getSolverMethod()[1]
792                 options[self.__PACKAGE_KEY]=self.getSolverPackage()
793                 options[self.__SYMMETRY_KEY]=self.isSymmetric()
794                 self.trace("PDE is resolved.")
795                 self.trace("solver options: %s"%str(options))
796                 self.__solution=mat.solve(f,options)
797              self.__solution_isValid=True
798           return self.__solution
799    
800     def setDebugOff(self):     def getFlux(self,u=None):
801         """
802         returns the flux M{J} for a given M{u}
803    
804         M{J[i,j]=A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]u[k]-X[i,j]}
805    
806         or
807    
808         M{J[j]=A[i,j]*grad(u)[l]+B[j]u-X[j]}
809    
810         @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.
811         @type u: L{Data<escript.Data>} or None
812         @return: flux
813         @rtype: L{Data<escript.Data>}
814         """
815         if u==None: u=self.getSolution()
816         return util.tensormult(self.getCoefficientOfGeneralPDE("A"),util.grad(u))+util.matrixmult(self.getCoefficientOfGeneralPDE("B"),u)-util.self.getCoefficientOfGeneralPDE("X")
817       # =============================================================================
818       #   solver settings:
819       # =============================================================================
820       def setSolverMethod(self,solver=None,preconditioner=None):
821         """         """
822         @brief         sets a new solver
823    
824           @param solver: sets a new solver method.
825           @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{PRES20}, L{LUMPING}, L{AMG}
826           @param preconditioner: sets a new solver method.
827           @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}
828         """         """
829         self.__debug=None         if solver==None: solve=self.DEFAULT
830           if preconditioner==None: preconditioner=self.DEFAULT
831           if not (solver,preconditioner)==self.getSolverMethod():
832               self.__solver_method=solver
833               self.__preconditioner=preconditioner
834               self.__checkMatrixType()
835               self.trace("New solver is %s"%self.getSolverMethodName())
836    
837     def debug(self):     def getSolverMethodName(self):
838         """         """
839         @brief returns true if the PDE is in the debug mode         returns the name of the solver currently used
840    
841           @return: the name of the solver currently used.
842           @rtype: C{string}
843         """         """
        return self.__debug  
844    
845     #===== Lumping ===========================         m=self.getSolverMethod()
846     def setLumpingOn(self):         p=self.getSolverPackage()
847        """         method=""
848        @brief indicates to use matrix lumping         if m[0]==self.DEFAULT: method="DEFAULT"
849        """         elif m[0]==self.DIRECT: method= "DIRECT"
850        if not self.isUsingLumping():         elif m[0]==self.ITERATIVE: method= "ITERATIVE"
851           if self.debug() : print "PDE Debug: lumping is set on"         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
852           self.__rebuildOperator()         elif m[0]==self.PCG: method= "PCG"
853           self.__lumping=True         elif m[0]==self.CR: method= "CR"
854           elif m[0]==self.CGS: method= "CGS"
855           elif m[0]==self.BICGSTAB: method= "BICGSTAB"
856           elif m[0]==self.SSOR: method= "SSOR"
857           elif m[0]==self.GMRES: method= "GMRES"
858           elif m[0]==self.PRES20: method= "PRES20"
859           elif m[0]==self.LUMPING: method= "LUMPING"
860           elif m[0]==self.AMG: method= "AMG"
861           if m[1]==self.DEFAULT: method+="+DEFAULT"
862           elif m[1]==self.JACOBI: method+= "+JACOBI"
863           elif m[1]==self.ILU0: method+= "+ILU0"
864           elif m[1]==self.ILUT: method+= "+ILUT"
865           elif m[1]==self.SSOR: method+= "+SSOR"
866           elif m[1]==self.AMG: method+= "+AMG"
867           elif m[1]==self.RILU: method+= "+RILU"
868           if p==self.DEFAULT: package="DEFAULT"
869           elif p==self.PASO: package= "PASO"
870           elif p==self.MKL: package= "MKL"
871           elif p==self.SCSL: package= "SCSL"
872           elif p==self.UMFPACK: package= "UMFPACK"
873           elif p==self.TRILINOS: package= "TRILINOS"
874           else : method="unknown"
875           return "%s solver of %s package"%(method,package)
876    
    def setLumpingOff(self):  
       """  
       @brief switches off matrix lumping  
       """  
       if self.isUsingLumping():  
          if self.debug() : print "PDE Debug: lumping is set off"  
          self.__rebuildOperator()  
          self.__lumping=False  
877    
878     def setLumping(self,flag=False):     def getSolverMethod(self):
879        """         """
880        @brief set the matrix lumping flag to flag         returns the solver method
       """  
       if flag:  
          self.setLumpingOn()  
       else:  
          self.setLumpingOff()  
881    
882     def isUsingLumping(self):         @return: the solver method currently be used.
883        """         @rtype: C{int}
884        @brief         """
885        """         return self.__solver_method,self.__preconditioner
       return self.__lumping  
886    
887     #============ method business =========================================================     def setSolverPackage(self,package=None):
    def setSolverMethod(self,solver=util.DEFAULT_METHOD):  
888         """         """
889         @brief sets a new solver         sets a new solver package
890    
891           @param package: sets a new solver method.
892           @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}, L{TRILINOS}
893         """         """
894         if not solver==self.getSolverMethod():         if package==None: package=self.DEFAULT
895             self.__solver_method=solver         if not package==self.getSolverPackage():
896             if self.debug() : print "PDE Debug: New solver is %s"%solver             self.__solver_package=package
897             self.__checkMatrixType()             self.__checkMatrixType()
898               self.trace("New solver is %s"%self.getSolverMethodName())
899    
900     def getSolverMethod(self):     def getSolverPackage(self):
901         """         """
902         @brief returns the solver method         returns the package of the solver
903    
904           @return: the solver package currently being used.
905           @rtype: C{int}
906         """         """
907         return self.__solver_method         return self.__solver_package
908    
909       def isUsingLumping(self):
910          """
911          checks if matrix lumping is used a solver method
912    
913          @return: True is lumping is currently used a solver method.
914          @rtype: C{bool}
915          """
916          return self.getSolverMethod()[0]==self.LUMPING
917    
    #============ tolerance business =========================================================  
918     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
919         """         """
920         @brief resets the tolerance to tol.         resets the tolerance for the solver method to tol where for an appropriate norm M{|.|}
921    
922           M{|L{getResidual}()|<tol*|L{getRightHandSide}()|}
923    
924           defines the stopping criterion.
925    
926           @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
927                       the system will be resolved.
928           @type tol: positive C{float}
929           @raise ValueError: if tolerance is not positive.
930         """         """
931         if not tol>0:         if not tol>0:
932             raise ValueException,"Tolerance as to be positive"             raise ValueError,"Tolerance as to be positive"
933         if tol<self.getTolerance(): self.__rebuildSolution()         if tol<self.getTolerance(): self.__invalidateSolution()
934         if self.debug() : print "PDE Debug: New tolerance %e",tol         self.trace("New tolerance %e"%tol)
935         self.__tolerance=tol         self.__tolerance=tol
936         return         return
937    
938     def getTolerance(self):     def getTolerance(self):
939         """         """
940         @brief returns the tolerance set for the solution         returns the tolerance set for the solution
941    
942           @return: tolerance currently used.
943           @rtype: C{float}
944         """         """
945         return self.__tolerance         return self.__tolerance
946    
947     #===== symmetry  flag ==========================     # =============================================================================
948       #    symmetry  flag:
949       # =============================================================================
950     def isSymmetric(self):     def isSymmetric(self):
951        """        """
952        @brief returns true is the operator is considered to be symmetric        checks if symmetry is indicated.
953    
954          @return: True is a symmetric PDE is indicated, otherwise False is returned
955          @rtype: C{bool}
956        """        """
957        return self.__sym        return self.__sym
958    
959     def setSymmetryOn(self):     def setSymmetryOn(self):
960        """        """
961        @brief sets the symmetry flag to true        sets the symmetry flag.
962        """        """
963        if not self.isSymmetric():        if not self.isSymmetric():
964           if self.debug() : print "PDE Debug: Operator is set to be symmetric"           self.trace("PDE is set to be symmetric")
965           self.__sym=True           self.__sym=True
966           self.__checkMatrixType()           self.__checkMatrixType()
967    
968     def setSymmetryOff(self):     def setSymmetryOff(self):
969        """        """
970        @brief sets the symmetry flag to false        removes the symmetry flag.
971        """        """
972        if self.isSymmetric():        if self.isSymmetric():
973           if self.debug() : print "PDE Debug: Operator is set to be unsymmetric"           self.trace("PDE is set to be unsymmetric")
974           self.__sym=False           self.__sym=False
975           self.__checkMatrixType()           self.__checkMatrixType()
976    
977     def setSymmetryTo(self,flag=False):     def setSymmetryTo(self,flag=False):
978       """        """
979       @brief sets the symmetry flag to flag        sets the symmetry flag to flag
980    
981       @param flag        @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.
982       """        @type flag: C{bool}
983       if flag:        """
984          self.setSymmetryOn()        if flag:
985       else:           self.setSymmetryOn()
986          self.setSymmetryOff()        else:
987             self.setSymmetryOff()
988    
989     #===== order reduction ==========================     # =============================================================================
990       # function space handling for the equation as well as the solution
991       # =============================================================================
992     def setReducedOrderOn(self):     def setReducedOrderOn(self):
993       """       """
994       @brief switches to on reduced order       switches on reduced order for solution and equation representation
995    
996         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
997       """       """
998       self.setReducedOrderForSolutionOn()       self.setReducedOrderForSolutionOn()
999       self.setReducedOrderForEquationOn()       self.setReducedOrderForEquationOn()
1000    
1001     def setReducedOrderOff(self):     def setReducedOrderOff(self):
1002       """       """
1003       @brief switches to full order       switches off reduced order for solution and equation representation
1004    
1005         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1006       """       """
1007       self.setReducedOrderForSolutionOff()       self.setReducedOrderForSolutionOff()
1008       self.setReducedOrderForEquationOff()       self.setReducedOrderForEquationOff()
1009    
1010     def setReducedOrderTo(self,flag=False):     def setReducedOrderTo(self,flag=False):
1011       """       """
1012       @brief sets order according to flag       sets order reduction for both solution and equation representation according to flag.
1013         @param flag: if flag is True, the order reduction is switched on for both  solution and equation representation, otherwise or
1014       @param flag                    if flag is not present order reduction is switched off
1015         @type flag: C{bool}
1016         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1017       """       """
1018       self.setReducedOrderForSolutionTo(flag)       self.setReducedOrderForSolutionTo(flag)
1019       self.setReducedOrderForEquationTo(flag)       self.setReducedOrderForEquationTo(flag)
                                                                                                                                                             
1020    
1021     #===== order reduction solution ==========================  
1022     def setReducedOrderForSolutionOn(self):     def setReducedOrderForSolutionOn(self):
1023       """       """
1024       @brief switches to reduced order to interpolate solution       switches on reduced order for solution representation
1025    
1026         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1027       """       """
1028       new_fs=escript.ReducedSolution(self.getDomain())       if not self.__reduce_solution_order:
1029       if self.getFunctionSpaceForSolution()!=new_fs:           if self.__altered_coefficients:
1030           if self.debug() : print "PDE Debug: Reduced order is used to interpolate solution."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
1031           self.__column_function_space=new_fs           self.trace("Reduced order is used to solution representation.")
1032           self.__rebuildSystem(deep=True)           self.__reduce_solution_order=True
1033             self.__resetSystem()
1034    
1035     def setReducedOrderForSolutionOff(self):     def setReducedOrderForSolutionOff(self):
1036       """       """
1037       @brief switches to full order to interpolate solution       switches off reduced order for solution representation
1038    
1039         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1040       """       """
1041       new_fs=escript.Solution(self.getDomain())       if self.__reduce_solution_order:
1042       if self.getFunctionSpaceForSolution()!=new_fs:           if self.__altered_coefficients:
1043           if self.debug() : print "PDE Debug: Full order is used to interpolate solution."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
1044           self.__column_function_space=new_fs           self.trace("Full order is used to interpolate solution.")
1045           self.__rebuildSystem(deep=True)           self.__reduce_solution_order=False
1046             self.__resetSystem()
1047    
1048     def setReducedOrderForSolutionTo(self,flag=False):     def setReducedOrderForSolutionTo(self,flag=False):
1049       """       """
1050       @brief sets order for test functions according to flag       sets order for test functions according to flag
1051    
1052       @param flag       @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or
1053                      if flag is not present order reduction is switched off
1054         @type flag: C{bool}
1055         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1056       """       """
1057       if flag:       if flag:
1058          self.setReducedOrderForSolutionOn()          self.setReducedOrderForSolutionOn()
1059       else:       else:
1060          self.setReducedOrderForSolutionOff()          self.setReducedOrderForSolutionOff()
1061                                                                                                                                                              
    #===== order reduction equation ==========================  
1062     def setReducedOrderForEquationOn(self):     def setReducedOrderForEquationOn(self):
1063       """       """
1064       @brief switches to reduced order for test functions       switches on reduced order for equation representation
1065    
1066         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1067       """       """
1068       new_fs=escript.ReducedSolution(self.getDomain())       if not self.__reduce_equation_order:
1069       if self.getFunctionSpaceForEquation()!=new_fs:           if self.__altered_coefficients:
1070           if self.debug() : print "PDE Debug: Reduced order is used for test functions."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
1071           self.__row_function_space=new_fs           self.trace("Reduced order is used for test functions.")
1072           self.__rebuildSystem(deep=True)           self.__reduce_equation_order=True
1073             self.__resetSystem()
1074    
1075     def setReducedOrderForEquationOff(self):     def setReducedOrderForEquationOff(self):
1076       """       """
1077       @brief switches to full order for test functions       switches off reduced order for equation representation
1078    
1079         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1080       """       """
1081       new_fs=escript.Solution(self.getDomain())       if self.__reduce_equation_order:
1082       if self.getFunctionSpaceForEquation()!=new_fs:           if self.__altered_coefficients:
1083           if self.debug() : print "PDE Debug: Full order is used for test functions."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
1084           self.__row_function_space=new_fs           self.trace("Full order is used for test functions.")
1085           self.__rebuildSystem(deep=True)           self.__reduce_equation_order=False
1086             self.__resetSystem()
1087    
1088     def setReducedOrderForEquationTo(self,flag=False):     def setReducedOrderForEquationTo(self,flag=False):
1089       """       """
1090       @brief sets order for test functions according to flag       sets order for test functions according to flag
1091    
1092       @param flag       @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or
1093                      if flag is not present order reduction is switched off
1094         @type flag: C{bool}
1095         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1096       """       """
1097       if flag:       if flag:
1098          self.setReducedOrderForEquationOn()          self.setReducedOrderForEquationOn()
1099       else:       else:
1100          self.setReducedOrderForEquationOff()          self.setReducedOrderForEquationOff()
1101                                                                                                                                                              
1102     # ==== initialization =====================================================================     # =============================================================================
1103     def __makeNewOperator(self):     # private method:
1104       # =============================================================================
1105       def __checkMatrixType(self):
1106         """
1107         reassess the matrix type and, if a new matrix is needed, resets the system.
1108         """
1109         new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverPackage(),self.isSymmetric())
1110         if not new_matrix_type==self.__matrix_type:
1111             self.trace("Matrix type is now %d."%new_matrix_type)
1112             self.__matrix_type=new_matrix_type
1113             self.__resetSystem()
1114       #
1115       #   rebuild switches :
1116       #
1117       def __invalidateSolution(self):
1118           """
1119           indicates the PDE has to be resolved if the solution is requested
1120           """
1121           if self.__solution_isValid: self.trace("PDE has to be resolved.")
1122           self.__solution_isValid=False
1123    
1124       def __invalidateOperator(self):
1125           """
1126           indicates the operator has to be rebuilt next time it is used
1127           """
1128           if self.__operator_is_Valid: self.trace("Operator has to be rebuilt.")
1129           self.__invalidateSolution()
1130           self.__operator_is_Valid=False
1131    
1132       def __invalidateRightHandSide(self):
1133           """
1134           indicates the right hand side has to be rebuild next time it is used
1135         """         """
1136         @brief         if self.__righthandside_isValid: self.trace("Right hand side has to be rebuilt.")
1137           self.__invalidateSolution()
1138           self.__righthandside_isValid=False
1139    
1140       def __invalidateSystem(self):
1141           """
1142           annonced that everthing has to be rebuild:
1143           """
1144           if self.__righthandside_isValid: self.trace("System has to be rebuilt.")
1145           self.__invalidateSolution()
1146           self.__invalidateOperator()
1147           self.__invalidateRightHandSide()
1148    
1149       def __resetSystem(self):
1150           """
1151           annonced that everthing has to be rebuild:
1152           """
1153           self.trace("New System is built from scratch.")
1154           self.__operator=escript.Operator()
1155           self.__operator_is_Valid=False
1156           self.__righthandside=escript.Data()
1157           self.__righthandside_isValid=False
1158           self.__solution=escript.Data()
1159           self.__solution_isValid=False
1160       #
1161       #    system initialization:
1162       #
1163       def __getNewOperator(self):
1164         """         """
1165           returns an instance of a new operator
1166           """
1167           self.trace("New operator is allocated.")
1168         return self.getDomain().newOperator( \         return self.getDomain().newOperator( \
1169                             self.getNumEquations(), \                             self.getNumEquations(), \
1170                             self.getFunctionSpaceForEquation(), \                             self.getFunctionSpaceForEquation(), \
# Line 545  class LinearPDE: Line 1172  class LinearPDE:
1172                             self.getFunctionSpaceForSolution(), \                             self.getFunctionSpaceForSolution(), \
1173                             self.__matrix_type)                             self.__matrix_type)
1174    
1175     def __makeNewRightHandSide(self):     def __getNewRightHandSide(self):
1176         """         """
1177         @brief         returns an instance of a new right hand side
1178         """         """
1179         return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)         self.trace("New right hand side is allocated.")
1180           if self.getNumEquations()>1:
1181               return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1182           else:
1183               return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1184    
1185     def __makeNewSolution(self):     def __getNewSolution(self):
1186         """         """
1187         @brief         returns an instance of a new solution
1188         """         """
1189         return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)         self.trace("New solution is allocated.")
1190           if self.getNumSolutions()>1:
1191               return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1192           else:
1193               return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1194    
1195     def __getFreshOperator(self):     def __makeFreshSolution(self):
1196         """         """
1197         @brief         makes sure that the solution is instantiated and returns it initialized by zeros
1198         """         """
1199         if self.__operator.isEmpty():         if self.__solution.isEmpty():
1200             self.__operator=self.__makeNewOperator()             self.__solution=self.__getNewSolution()
            if self.debug() : print "PDE Debug: New operator allocated"  
1201         else:         else:
1202             self.__operator.setValue(0.)             self.__solution*=0
1203             self.__operator.resetSolver()             self.trace("Solution is reset to zero.")
1204             if self.debug() : print "PDE Debug: Operator reset to zero"         return self.__solution
        return self.__operator  
1205    
1206     def __getFreshRightHandSide(self):     def __makeFreshRightHandSide(self):
1207         """         """
1208         @brief         makes sure that the right hand side is instantiated and returns it initialized by zeros
1209         """         """
1210         if self.__righthandside.isEmpty():         if self.__righthandside.isEmpty():
1211             self.__righthandside=self.__makeNewRightHandSide()             self.__righthandside=self.__getNewRightHandSide()
            if self.debug() : print "PDE Debug: New right hand side allocated"  
1212         else:         else:
1213             print "fix self.__righthandside*=0"             self.__righthandside*=0
1214             self.__righthandside*=0.             self.trace("Right hand side is reset to zero.")
1215             if self.debug() : print "PDE Debug: Right hand side reset to zero"         return self.__righthandside
        return  self.__righthandside  
1216    
1217     #============ some serivice functions  =====================================================     def __makeFreshOperator(self):
1218     def getDomain(self):         """
1219           makes sure that the operator is instantiated and returns it initialized by zeros
1220           """
1221           if self.__operator.isEmpty():
1222               self.__operator=self.__getNewOperator()
1223           else:
1224               self.__operator.resetValues()
1225               self.trace("Operator reset to zero")
1226           return self.__operator
1227    
1228       def __applyConstraint(self):
1229           """
1230           applies the constraints defined by q and r to the system
1231           """
1232           if not self.isUsingLumping():
1233              q=self.getCoefficientOfGeneralPDE("q")
1234              r=self.getCoefficientOfGeneralPDE("r")
1235              if not q.isEmpty() and not self.__operator.isEmpty():
1236                 # q is the row and column mask to indicate where constraints are set:
1237                 row_q=escript.Data(q,self.getFunctionSpaceForEquation())
1238                 col_q=escript.Data(q,self.getFunctionSpaceForSolution())
1239                 u=self.__getNewSolution()
1240                 if r.isEmpty():
1241                    r_s=self.__getNewSolution()
1242                 else:
1243                    r_s=escript.Data(r,self.getFunctionSpaceForSolution())
1244                 u.copyWithMask(r_s,col_q)
1245                 if not self.__righthandside.isEmpty():
1246                    self.__righthandside-=self.__operator*u
1247                    self.__righthandside=self.copyConstraint(self.__righthandside)
1248                 self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
1249       # =============================================================================
1250       # function giving access to coefficients of the general PDE:
1251       # =============================================================================
1252       def getCoefficientOfGeneralPDE(self,name):
1253         """
1254         return the value of the coefficient name of the general PDE.
1255    
1256         @note: This method is called by the assembling routine it can be overwritten
1257               to map coefficients of a particular PDE to the general PDE.
1258         @param name: name of the coefficient requested.
1259         @type name: C{string}
1260         @return: the value of the coefficient  name
1261         @rtype: L{Data<escript.Data>}
1262         @raise IllegalCoefficient: if name is not one of coefficients
1263                      M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
1264       """       """
1265       @brief returns the domain of the PDE       if self.hasCoefficientOfGeneralPDE(name):
1266            return self.getCoefficient(name)
1267         else:
1268            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1269    
1270       def hasCoefficientOfGeneralPDE(self,name):
1271       """       """
1272       return self.__domain       checks if name is a the name of a coefficient of the general PDE.
1273    
1274         @param name: name of the coefficient enquired.
1275         @type name: C{string}
1276         @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1277         @rtype: C{bool}
1278    
    def getDim(self):  
1279       """       """
1280       @brief returns the spatial dimension of the PDE       return self.__COEFFICIENTS_OF_GENEARL_PDE.has_key(name)
1281    
1282       def createCoefficientOfGeneralPDE(self,name):
1283       """       """
1284       return self.getDomain().getDim()       returns a new instance of a coefficient for coefficient name of the general PDE
1285    
1286     def getNumEquations(self):       @param name: name of the coefficient requested.
1287         @type name: C{string}
1288         @return: a coefficient name initialized to 0.
1289         @rtype: L{Data<escript.Data>}
1290         @raise IllegalCoefficient: if name is not one of coefficients
1291                      M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
1292       """       """
1293       @brief returns the number of equations       if self.hasCoefficientOfGeneralPDE(name):
1294            return escript.Data(0,self.getShapeOfCoefficientOfGeneralPDE(name),self.getFunctionSpaceForCoefficientOfGeneralPDE(name))
1295         else:
1296            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1297    
1298       def getFunctionSpaceForCoefficientOfGeneralPDE(self,name):
1299       """       """
1300       if self.__numEquations>0:       return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name of the general PDE
1301           return self.__numEquations  
1302         @param name: name of the coefficient enquired.
1303         @type name: C{string}
1304         @return: the function space to be used for coefficient name
1305         @rtype: L{FunctionSpace<escript.FunctionSpace>}
1306         @raise IllegalCoefficient: if name is not one of coefficients
1307                      M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
1308         """
1309         if self.hasCoefficientOfGeneralPDE(name):
1310            return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getFunctionSpace(self.getDomain())
1311       else:       else:
1312           raise ValueError,"Number of equations is undefined. Please specify argument numEquations."          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1313    
1314     def getNumSolutions(self):     def getShapeOfCoefficientOfGeneralPDE(self,name):
1315       """       """
1316       @brief returns the number of unknowns       return the shape of the coefficient name of the general PDE
1317    
1318         @param name: name of the coefficient enquired.
1319         @type name: C{string}
1320         @return: the shape of the coefficient name
1321         @rtype: C{tuple} of C{int}
1322         @raise IllegalCoefficient: if name is not one of coefficients
1323                      M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
1324       """       """
1325       if self.__numSolutions>0:       if self.hasCoefficientOfGeneralPDE(name):
1326          return self.__numSolutions          return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1327       else:       else:
1328          raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1329    
1330       # =============================================================================
1331       # functions giving access to coefficients of a particular PDE implementation:
1332       # =============================================================================
1333       def getCoefficient(self,name):
1334         """
1335         returns the value of the coefficient name
1336    
1337     def checkSymmetry(self,verbose=True):       @param name: name of the coefficient requested.
1338        """       @type name: C{string}
1339        @brief returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.       @return: the value of the coefficient name
1340        """       @rtype: L{Data<escript.Data>}
1341        verbose=verbose or self.debug()       @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1342        out=True       """
1343        if self.getNumSolutions()!=self.getNumEquations():       if self.hasCoefficient(name):
1344           if verbose: print "non-symmetric PDE because of different number of equations and solutions"           return self.COEFFICIENTS[name].getValue()
1345           out=False       else:
1346        else:          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
          A=self.getCoefficientOfPDE("A")  
          if not A.isEmpty():  
             tol=util.Lsup(A)*self.TOL  
             if self.getNumSolutions()>1:  
                for i in range(self.getNumEquations()):  
                   for j in range(self.getDim()):  
                      for k in range(self.getNumSolutions()):  
                         for l in range(self.getDim()):  
                             if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:  
                                if verbose: print "non-symmetric PDE because A[%d,%d,%d,%d]!=A[%d,%d,%d,%d]"%(i,j,k,l,k,l,i,j)  
                                out=False  
             else:  
                for j in range(self.getDim()):  
                   for l in range(self.getDim()):  
                      if util.Lsup(A[j,l]-A[l,j])>tol:  
                         if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)  
                         out=False  
          B=self.getCoefficientOfPDE("B")  
          C=self.getCoefficientOfPDE("C")  
          if B.isEmpty() and not C.isEmpty():  
             if verbose: print "non-symmetric PDE because B is not present but C is"  
             out=False  
          elif not B.isEmpty() and C.isEmpty():  
             if verbose: print "non-symmetric PDE because C is not present but B is"  
             out=False  
          elif not B.isEmpty() and not C.isEmpty():  
             tol=(util.Lsup(B)+util.Lsup(C))*self.TOL/2.  
             if self.getNumSolutions()>1:  
                for i in range(self.getNumEquations()):  
                    for j in range(self.getDim()):  
                       for k in range(self.getNumSolutions()):  
                          if util.Lsup(B[i,j,k]-C[i,j,k])>tol:  
                               if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,i,j,k)  
                               out=False  
             else:  
                for j in range(self.getDim()):  
                   if util.Lsup(B[j]-C[j])>tol:  
                      if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)  
                      out=False  
          if self.getNumSolutions()>1:  
            D=self.getCoefficientOfPDE("D")  
            if not D.isEmpty():  
              tol=util.Lsup(D)*self.TOL  
              for i in range(self.getNumEquations()):  
                 for k in range(self.getNumSolutions()):  
                   if util.Lsup(D[i,k]-D[k,i])>tol:  
                       if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)  
                       out=False  
               
       return out  
1347    
1348     def getFlux(self,u):     def hasCoefficient(self,name):
1349         """       """
1350         @brief returns the flux J_ij for a given u       return True if name is the name of a coefficient
1351    
1352         @param name: name of the coefficient enquired.
1353         @type name: C{string}
1354         @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1355         @rtype: C{bool}
1356         """
1357         return self.COEFFICIENTS.has_key(name)
1358    
1359              J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}     def createCoefficient(self, name):
1360         """
1361         create a L{Data<escript.Data>} object corresponding to coefficient name
1362    
1363         @param u argument of the operator       @return: a coefficient name initialized to 0.
1364         @rtype: L{Data<escript.Data>}
1365         @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1366         """
1367         if self.hasCoefficient(name):
1368            return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1369         else:
1370            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1371    
1372         """     def getFunctionSpaceForCoefficient(self,name):
1373         raise SystemError,"getFlux is not implemented yet"       """
1374         return None       return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name
1375    
1376     def applyOperator(self,u):       @param name: name of the coefficient enquired.
1377         """       @type name: C{string}
1378         @brief applies the operator of the PDE to a given solution u in weak from       @return: the function space to be used for coefficient name
1379         @rtype: L{FunctionSpace<escript.FunctionSpace>}
1380         @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1381         """
1382         if self.hasCoefficient(name):
1383            return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1384         else:
1385            raise ValueError,"unknown coefficient %s requested"%name
1386       def getShapeOfCoefficient(self,name):
1387         """
1388         return the shape of the coefficient name
1389    
1390         @param u argument of the operator       @param name: name of the coefficient enquired.
1391         @type name: C{string}
1392         @return: the shape of the coefficient name
1393         @rtype: C{tuple} of C{int}
1394         @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1395         """
1396         if self.hasCoefficient(name):
1397            return self.COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1398         else:
1399            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1400    
1401         """     def resetCoefficients(self):
1402         return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())       """
1403                                                                                                                                                                   resets all coefficients to there default values.
1404     def getResidual(self,u):       """
1405         """       for i in self.COEFFICIENTS.iterkeys():
1406         @brief return the residual of u in the weak from           self.COEFFICIENTS[i].resetValue()
1407    
1408         @param u     def alteredCoefficient(self,name):
1409         """       """
1410         return self.applyOperator(u)-self.getRightHandSide()       announce that coefficient name has been changed
1411    
1412         @param name: name of the coefficient enquired.
1413         @type name: C{string}
1414         @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1415         @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.
1416         """
1417         if self.hasCoefficient(name):
1418            self.trace("Coefficient %s has been altered."%name)
1419            if not ((name=="q" or name=="r") and self.isUsingLumping()):
1420               if self.COEFFICIENTS[name].isAlteringOperator(): self.__invalidateOperator()
1421               if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__invalidateRightHandSide()
1422         else:
1423            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1424    
1425       def copyConstraint(self,u):
1426          """
1427          copies the constraint into u and returns u.
1428    
1429          @param u: a function of rank 0 is a single PDE is solved and of shape (numSolution,) for a system of PDEs
1430          @type u: L{Data<escript.Data>}
1431          @return: the input u modified by the constraints.
1432          @rtype: L{Data<escript.Data>}
1433          @warning: u is altered if it has the appropriate L{FunctionSpace<escript.FunctionSpace>}
1434          """
1435          q=self.getCoefficientOfGeneralPDE("q")
1436          r=self.getCoefficientOfGeneralPDE("r")
1437          if not q.isEmpty():
1438             if u.isEmpty(): u=escript.Data(0.,q.getShape(),q.getFunctionSpace())
1439             if r.isEmpty():
1440                 r=escript.Data(0,u.getShape(),u.getFunctionSpace())
1441             else:
1442                 r=escript.Data(r,u.getFunctionSpace())
1443             u.copyWithMask(r,escript.Data(q,u.getFunctionSpace()))
1444          return u
1445    
1446     def _setValue(self,**coefficients):     def setValue(self,**coefficients):
1447        """        """
1448        @brief sets new values to coefficient        sets new values to coefficients
1449    
1450        @param coefficients        @param coefficients: new values assigned to coefficients
1451          @keyword A: value for coefficient A.
1452          @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1453          @keyword B: value for coefficient B
1454          @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1455          @keyword C: value for coefficient C
1456          @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1457          @keyword D: value for coefficient D
1458          @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1459          @keyword X: value for coefficient X
1460          @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1461          @keyword Y: value for coefficient Y
1462          @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1463          @keyword d: value for coefficient d
1464          @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1465          @keyword y: value for coefficient y
1466          @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1467          @keyword d_contact: value for coefficient d_contact
1468          @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1469                           or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1470          @keyword y_contact: value for coefficient y_contact
1471          @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1472                           or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1473          @keyword r: values prescribed to the solution at the locations of constraints
1474          @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1475                   depending of reduced order is used for the solution.
1476          @keyword q: mask for location of constraints
1477          @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1478                   depending of reduced order is used for the representation of the equation.
1479          @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1480        """        """
1481        # check if the coefficients are  legal:        # check if the coefficients are  legal:
1482        for i in coefficients.iterkeys():        for i in coefficients.iterkeys():
1483           if not self.hasCoefficient(i):           if not self.hasCoefficient(i):
1484              raise ValueError,"Attempt to set unknown coefficient %s"%i              raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
1485        # if the number of unknowns or equations is still unknown we try to estimate them:        # if the number of unknowns or equations is still unknown we try to estimate them:
1486        if self.__numEquations<1 or self.__numSolutions<1:        if self.__numEquations==None or self.__numSolutions==None:
1487           for i,d in coefficients.iteritems():           for i,d in coefficients.iteritems():
1488              if hasattr(d,"shape"):              if hasattr(d,"shape"):
1489                  s=d.shape                  s=d.shape
# Line 726  class LinearPDE: Line 1493  class LinearPDE:
1493                  s=numarray.array(d).shape                  s=numarray.array(d).shape
1494              if s!=None:              if s!=None:
1495                  # get number of equations and number of unknowns:                  # get number of equations and number of unknowns:
1496                  res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(s,self.getDim())                  res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1497                  if res==None:                  if res==None:
1498                      raise ValueError,"Illegal shape %s of coefficient %s"%(s,i)                      raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1499                  else:                  else:
1500                      if self.__numEquations<1: self.__numEquations=res[0]                      if self.__numEquations==None: self.__numEquations=res[0]
1501                      if self.__numSolutions<1: self.__numSolutions=res[1]                      if self.__numSolutions==None: self.__numSolutions=res[1]
1502        if self.__numEquations<1: raise ValueError,"unidententified number of equations"        if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"
1503        if self.__numSolutions<1: raise ValueError,"unidententified number of solutions"        if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions"
1504        # 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:
1505        for i,d in coefficients.iteritems():        for i,d in coefficients.iteritems():
1506          if d==None:          try:
1507               d2=escript.Data()             self.COEFFICIENTS[i].setValue(self.getDomain(),self.getNumEquations(),self.getNumSolutions(),self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1508          elif isinstance(d,escript.Data):          except IllegalCoefficientValue,m:
1509               if d.isEmpty():             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
                 d2=d  
              else:  
                 d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))  
         else:  
               d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))  
         if not d2.isEmpty():  
            if not self.getShapeOfCoefficient(i)==d2.getShape():  
                raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),d2.getShape())  
         # overwrite new values:  
         if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i  
         self.COEFFICIENTS[i].setValue(d2)  
1510          self.alteredCoefficient(i)          self.alteredCoefficient(i)
         
       # reset the HomogeneousConstraintFlag:  
       self.__setHomogeneousConstraintFlag()  
       if len(coefficients)>0 and not self.isUsingLumping() and not self.__homogeneous_constraint: self.__rebuildSystem()  
   
    def __setHomogeneousConstraintFlag(self):  
       """  
       @brief checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.  
       """  
       self.__homogeneous_constraint=True  
       q=self.getCoefficientOfPDE("q")  
       r=self.getCoefficientOfPDE("r")  
       if not q.isEmpty() and not r.isEmpty():  
          if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False  
       if self.debug():  
            if self.__homogeneous_constraint:  
                print "PDE Debug: Constraints are homogeneous."  
            else:  
                print "PDE Debug: Constraints are inhomogeneous."  
   
   
    # ==== rebuild switches =====================================================================  
    def __rebuildSolution(self,deep=False):  
        """  
        @brief indicates the PDE has to be reolved if the solution is requested  
        """  
        if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."  
        self.__solution_isValid=False  
        if deep: self.__solution=escript.Data()  
   
1511    
1512     def __rebuildOperator(self,deep=False):        self.__altered_coefficients=True
1513         """        # check if the systrem is inhomogeneous:
1514         @brief indicates the operator has to be rebuilt next time it is used        if len(coefficients)>0 and not self.isUsingLumping():
1515         """           q=self.getCoefficientOfGeneralPDE("q")
1516         if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."           r=self.getCoefficientOfGeneralPDE("r")
1517         self.__rebuildSolution(deep)           homogeneous_constraint=True
1518         self.__operator_isValid=False           if not q.isEmpty() and not r.isEmpty():
1519         if deep: self.__operator=escript.Operator()               if util.Lsup(q*r)>0.:
1520                   self.trace("Inhomogeneous constraint detected.")
1521     def __rebuildRightHandSide(self,deep=False):                 self.__invalidateSystem()
        """  
        @brief indicates the right hand side has to be rebuild next time it is used  
        """  
        if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."  
        self.__rebuildSolution(deep)  
        self.__righthandside_isValid=False  
        if deep: self.__righthandside=escript.Data()  
   
    def __rebuildSystem(self,deep=False):  
        """  
        @brief annonced that all coefficient name has been changed  
        """  
        self.__rebuildSolution(deep)  
        self.__rebuildOperator(deep)  
        self.__rebuildRightHandSide(deep)  
     
    def __checkMatrixType(self):  
      """  
      @brief reassess the matrix type and, if needed, initiates an operator rebuild  
      """  
      new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())  
      if not new_matrix_type==self.__matrix_type:  
          if self.debug() : print "PDE Debug: Matrix type is now %d."%new_matrix_type  
          self.__matrix_type=new_matrix_type  
          self.__rebuildOperator(deep=True)  
   
    #============ assembling =======================================================  
    def __copyConstraint(self):  
       """  
       @brief copies the constrint condition into u  
       """  
       if not self.__righthandside.isEmpty():  
          q=self.getCoefficientOfPDE("q")  
          r=self.getCoefficientOfPDE("r")  
          if not q.isEmpty():  
              if r.isEmpty():  
                 r2=escript.Data(0,self.__righthandside.getShape(),self.__righthandside.getFunctionSpace())  
              else:  
                 r2=escript.Data(r,self.__righthandside.getFunctionSpace())  
              self.__righthandside.copyWithMask(r2,escript.Data(q,self.__righthandside.getFunctionSpace()))  
   
    def __applyConstraint(self):  
        """  
        @brief applies the constraints defined by q and r to the system  
        """  
        q=self.getCoefficientOfPDE("q")  
        r=self.getCoefficientOfPDE("r")  
        if not q.isEmpty() and not self.__operator.isEmpty():  
           # q is the row and column mask to indicate where constraints are set:  
           row_q=escript.Data(q,self.getFunctionSpaceForEquation())  
           col_q=escript.Data(q,self.getFunctionSpaceForSolution())  
           u=self.__makeNewSolution()  
           if r.isEmpty():  
              r_s=self.__makeNewSolution()  
           else:  
              r_s=escript.Data(r,self.getFunctionSpaceForSolution())  
           u.copyWithMask(r_s,col_q)  
           if self.isUsingLumping():  
              self.__operator.copyWithMask(escript.Data(1,q.getShape(),self.getFunctionSpaceForEquation()),row_q)  
           else:  
              if not self.__righthandside.isEmpty(): self.__righthandside-=self.__operator*u  
              self.__operator.nullifyRowsAndCols(row_q,col_q,1.)  
1522    
1523     def getSystem(self):     def getSystem(self):
1524         """         """
1525         @brief return the operator and right hand side of the PDE         return the operator and right hand side of the PDE
1526    
1527           @return: the discrete version of the PDE
1528           @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1529         """         """
1530         if not self.__operator_isValid or not self.__righthandside_isValid:         if not self.__operator_is_Valid or not self.__righthandside_isValid:
1531            if self.isUsingLumping():            if self.isUsingLumping():
1532                if not self.__operator_isValid:                if not self.__operator_is_Valid:
1533                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1534                         raise TypeError,"Lumped matrix requires same order for equations and unknowns"                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1535                   if not self.getCoefficientOfPDE("A").isEmpty():                   if not self.getCoefficientOfGeneralPDE("A").isEmpty():
1536                            raise Warning,"Lumped matrix does not allow coefficient A"                        raise ValueError,"coefficient A in lumped matrix may not be present."
1537                   if not self.getCoefficientOfPDE("B").isEmpty():                   if not self.getCoefficientOfGeneralPDE("B").isEmpty():
1538                            raise Warning,"Lumped matrix does not allow coefficient B"                        raise ValueError,"coefficient A in lumped matrix may not be present."
1539                   if not self.getCoefficientOfPDE("C").isEmpty():                   if not self.getCoefficientOfGeneralPDE("C").isEmpty():
1540                            raise Warning,"Lumped matrix does not allow coefficient C"                        raise ValueError,"coefficient A in lumped matrix may not be present."
1541                   if self.debug() : print "PDE Debug: New lumped operator is built."                   D=self.getCoefficientOfGeneralPDE("D")
1542                   mat=self.__makeNewOperator()                   if not D.isEmpty():
1543                   self.getDomain().addPDEToSystem(mat,escript.Data(), \                       if self.getNumSolutions()>1:
1544                             self.getCoefficientOfPDE("A"), \                          #D_times_e=util.matrixmult(D,numarray.ones((self.getNumSolutions(),)))
1545                             self.getCoefficientOfPDE("B"), \                          D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1546                             self.getCoefficientOfPDE("C"), \                       else:
1547                             self.getCoefficientOfPDE("D"), \                          D_times_e=D
1548                             escript.Data(), \                   else:
1549                             escript.Data(), \                      D_times_e=escript.Data()
1550                             self.getCoefficientOfPDE("d"), \                   d=self.getCoefficientOfGeneralPDE("d")
1551                             escript.Data(),\                   if not d.isEmpty():
1552                             self.getCoefficientOfPDE("d_contact"), \                       if self.getNumSolutions()>1:
1553                             escript.Data())                          #d_times_e=util.matrixmult(d,numarray.ones((self.getNumSolutions(),)))
1554                   self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)                          d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1555                   self.__applyConstraint()                       else:
1556                   self.__operator_isValid=True                          d_times_e=d
1557                     else:
1558                        d_times_e=escript.Data()
1559                     d_contact=self.getCoefficientOfGeneralPDE("d_contact")
1560                     if not d_contact.isEmpty():
1561                         if self.getNumSolutions()>1:
1562                            d_contact_times_e=util.matrixmult(d_contact,numarray.ones((self.getNumSolutions(),)))
1563                         else:
1564                            d_contact_times_e=d_contact
1565                     else:
1566                        d_contact_times_e=escript.Data()
1567        
1568                     self.__operator=self.__getNewRightHandSide()
1569                     self.getDomain().addPDEToRHS(self.__operator, \
1570                                                  escript.Data(), \
1571                                                  D_times_e, \
1572                                                  d_times_e,\
1573                                                  d_contact_times_e)
1574                     self.__operator=1./self.__operator
1575                     self.trace("New lumped operator has been built.")
1576                     self.__operator_is_Valid=True
1577                if not self.__righthandside_isValid:                if not self.__righthandside_isValid:
1578                   if self.debug() : print "PDE Debug: New right hand side is built."                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1579                   self.getDomain().addPDEToRHS(self.__getFreshRightHandSide(), \                                 self.getCoefficientOfGeneralPDE("X"), \
1580                                 self.getCoefficientOfPDE("X"), \                                 self.getCoefficientOfGeneralPDE("Y"),\
1581                                 self.getCoefficientOfPDE("Y"),\                                 self.getCoefficientOfGeneralPDE("y"),\
1582                                 self.getCoefficientOfPDE("y"),\                                 self.getCoefficientOfGeneralPDE("y_contact"))
1583                                 self.getCoefficientOfPDE("y_contact"))                   self.trace("New right hand side as been built.")
                  self.__copyConstraint()  
1584                   self.__righthandside_isValid=True                   self.__righthandside_isValid=True
1585            else:            else:
1586               if not self.__operator_isValid and not self.__righthandside_isValid:               if not self.__operator_is_Valid and not self.__righthandside_isValid:
1587                   if self.debug() : print "PDE Debug: New system is built."                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
1588                   self.getDomain().addPDEToSystem(self.__getFreshOperator(),self.__getFreshRightHandSide(), \                                 self.getCoefficientOfGeneralPDE("A"), \
1589                                 self.getCoefficientOfPDE("A"), \                                 self.getCoefficientOfGeneralPDE("B"), \
1590                                 self.getCoefficientOfPDE("B"), \                                 self.getCoefficientOfGeneralPDE("C"), \
1591                                 self.getCoefficientOfPDE("C"), \                                 self.getCoefficientOfGeneralPDE("D"), \
1592                                 self.getCoefficientOfPDE("D"), \                                 self.getCoefficientOfGeneralPDE("X"), \
1593                                 self.getCoefficientOfPDE("X"), \                                 self.getCoefficientOfGeneralPDE("Y"), \
1594                                 self.getCoefficientOfPDE("Y"), \                                 self.getCoefficientOfGeneralPDE("d"), \
1595                                 self.getCoefficientOfPDE("d"), \                                 self.getCoefficientOfGeneralPDE("y"), \
1596                                 self.getCoefficientOfPDE("y"), \                                 self.getCoefficientOfGeneralPDE("d_contact"), \
1597                                 self.getCoefficientOfPDE("d_contact"), \                                 self.getCoefficientOfGeneralPDE("y_contact"))
                                self.getCoefficientOfPDE("y_contact"))  
1598                   self.__applyConstraint()                   self.__applyConstraint()
1599                   self.__copyConstraint()                   self.__righthandside=self.copyConstraint(self.__righthandside)
1600                   self.__operator_isValid=True                   self.trace("New system has been built.")
1601                     self.__operator_is_Valid=True
1602                   self.__righthandside_isValid=True                   self.__righthandside_isValid=True
1603               elif not self.__righthandside_isValid:               elif not self.__righthandside_isValid:
1604                   if self.debug() : print "PDE Debug: New right hand side is built."                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1605                   self.getDomain().addPDEToRHS(self.__getFreshRightHandSide(), \                                 self.getCoefficientOfGeneralPDE("X"), \
1606                                 self.getCoefficientOfPDE("X"), \                                 self.getCoefficientOfGeneralPDE("Y"),\
1607                                 self.getCoefficientOfPDE("Y"),\                                 self.getCoefficientOfGeneralPDE("y"),\
1608                                 self.getCoefficientOfPDE("y"),\                                 self.getCoefficientOfGeneralPDE("y_contact"))
1609                                 self.getCoefficientOfPDE("y_contact"))                   self.__righthandside=self.copyConstraint(self.__righthandside)
1610                   self.__copyConstraint()                   self.trace("New right hand side has been built.")
1611                   self.__righthandside_isValid=True                   self.__righthandside_isValid=True
1612               elif not self.__operator_isValid:               elif not self.__operator_is_Valid:
1613                   if self.debug() : print "PDE Debug: New operator is built."                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
1614                   self.getDomain().addPDEToSystem(self.__getFreshOperator(),escript.Data(), \                              self.getCoefficientOfGeneralPDE("A"), \
1615                              self.getCoefficientOfPDE("A"), \                              self.getCoefficientOfGeneralPDE("B"), \
1616                              self.getCoefficientOfPDE("B"), \                              self.getCoefficientOfGeneralPDE("C"), \
1617                              self.getCoefficientOfPDE("C"), \                              self.getCoefficientOfGeneralPDE("D"), \
                             self.getCoefficientOfPDE("D"), \  
1618                              escript.Data(), \                              escript.Data(), \
1619                              escript.Data(), \                              escript.Data(), \
1620                              self.getCoefficientOfPDE("d"), \                              self.getCoefficientOfGeneralPDE("d"), \
1621                              escript.Data(),\                              escript.Data(),\
1622                              self.getCoefficientOfPDE("d_contact"), \                              self.getCoefficientOfGeneralPDE("d_contact"), \
1623                              escript.Data())                              escript.Data())
1624                   self.__applyConstraint()                   self.__applyConstraint()
1625                   self.__operator_isValid=True                   self.trace("New operator has been built.")
1626                     self.__operator_is_Valid=True
1627         return (self.__operator,self.__righthandside)         return (self.__operator,self.__righthandside)
    def getOperator(self):  
        """  
        @brief returns the operator of the PDE  
        """  
        return self.getSystem()[0]  
1628    
    def getRightHandSide(self):  
        """  
        @brief returns the right hand side of the PDE  
        """  
        return self.getSystem()[1]  
1629    
1630     def solve(self,**options):  class Poisson(LinearPDE):
1631        """     """
1632        @brief solve the PDE     Class to define a Poisson equation problem, which is genear L{LinearPDE} of the form
1633    
1634        @param options     M{-grad(grad(u)[j])[j] = f}
       """  
       mat,f=self.getSystem()  
       if self.isUsingLumping():  
          out=f/mat  
       else:  
          options[util.TOLERANCE_KEY]=self.getTolerance()  
          options[util.METHOD_KEY]=self.getSolverMethod()  
          options[util.SYMMETRY_KEY]=self.isSymmetric()  
          if self.debug() : print "PDE Debug: solver options: ",options  
          out=mat.solve(f,options)  
       return out  
1635    
1636     def getSolution(self,**options):     with natural boundary conditons
        """  
        @brief returns the solution of the PDE  
1637    
1638         @param options     M{n[j]*grad(u)[j] = 0 }
1639         """  
1640         if not self.__solution_isValid:     and constraints:
1641             if self.debug() : print "PDE Debug: PDE is resolved."  
1642             self.__solution=self.solve(**options)     M{u=0} where M{q>0}
            self.__solution_isValid=True  
        return self.__solution  
1643    
 class AdvectivePDE(LinearPDE):  
1644     """     """
    @brief Class to handel a linear PDE domineated by advective terms:  
     
    class to define a linear PDE of the form  
1645    
1646       -(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i     def __init__(self,domain,debug=False):
1647         """
1648         initializes a new Poisson equation
1649    
1650       with boundary conditons:       @param domain: domain of the PDE
1651         @type domain: L{Domain<escript.Domain>}
1652         @param debug: if True debug informations are printed.
1653    
1654          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i       """
1655         super(Poisson, self).__init__(domain,1,1,debug)
1656         self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1657                              "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1658         self.setSymmetryOn()
1659    
1660      and contact conditions     def setValue(self,**coefficients):
1661         """
1662         sets new values to coefficients
1663    
1664          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i       @param coefficients: new values assigned to coefficients
1665         @keyword f: value for right hand side M{f}
1666         @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1667         @keyword q: mask for location of constraints
1668         @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>}
1669                   depending of reduced order is used for the representation of the equation.
1670         @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1671         """
1672         super(Poisson, self).setValue(**coefficients)
1673    
1674       def getCoefficientOfGeneralPDE(self,name):
1675         """
1676         return the value of the coefficient name of the general PDE
1677         @param name: name of the coefficient requested.
1678         @type name: C{string}
1679         @return: the value of the coefficient  name
1680         @rtype: L{Data<escript.Data>}
1681         @raise IllegalCoefficient: if name is not one of coefficients
1682                      M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
1683         @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1684         """
1685         if name == "A" :
1686             return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1687         elif name == "B" :
1688             return escript.Data()
1689         elif name == "C" :
1690             return escript.Data()
1691         elif name == "D" :
1692             return escript.Data()
1693         elif name == "X" :
1694             return escript.Data()
1695         elif name == "Y" :
1696             return self.getCoefficient("f")
1697         elif name == "d" :
1698             return escript.Data()
1699         elif name == "y" :
1700             return escript.Data()
1701         elif name == "d_contact" :
1702             return escript.Data()
1703         elif name == "y_contact" :
1704             return escript.Data()
1705         elif name == "r" :
1706             return escript.Data()
1707         elif name == "q" :
1708             return self.getCoefficient("q")
1709         else:
1710            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1711    
1712      and constraints:  class Helmholtz(LinearPDE):
1713       """
1714       Class to define a Helmhotz equation problem, which is genear L{LinearPDE} of the form
1715    
1716           u_i=r_i where q_i>0     M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
1717    
1718      The PDE is solved by stabilizing the advective terms using SUPG approach:     with natural boundary conditons
1719    
1720         A_{ijkl}<-A_{ijkl}+0.5*h*(xi(b_{ik})*B_{ijk}*B_{ilk}/length(B_{i:k})^2)+0.5*h*xi_{c_{ik}}*(C_{ikj}*C_{ikl}/length(C_{ik:})^2)     M{k*n[j]*grad(u)[j] = g- S{alpha}u }
1721    
1722      where       and constraints:
1723    
1724             b_{ik}=length(B_{i:k})*h/2/length(A_{i:k:})     M{u=r} where M{q>0}
            c_{ik}=length(C_{i:k})*h/2/length(A_{i:k:})  
1725    
                       alpha/3        alpha<3  
            xi(alpha)=          for                  approximating cotanh(alpha)-1/alpha  
                        1             alpha>=3  
1726     """     """
    def __getXi(self,alpha):  
          c=alpha-3.  
          return c*c.whereNegative()/3.+1.  
   
    def __getUpdateVector(self,V,hover2,alphaByU):  
      v=util.length(V)  
      v_max=util.Lsup(v)  
      if v_max>0:  
          V/=v+v_max*self.TOL  
          alpha=alphaByU*v  
          A_bar=v*hover2*self.__getXi(alpha)  
          print "-------------"  
          print "@ max alpha ",util.Lsup(alpha)  
          print "-------------"  
      else:  
          A_bar=1.  
      return V,A_bar  
   
    def __getAlphaByU(self,A,hover2):  
       a=util.length(A)  
       a_max=util.Lsup(a)  
       if a_max>0:  
          return hover2/(a+a_max*self.TOL)  
       else:  
          return 1./self.TOL  
   
1727    
1728     def getCoefficientOfPDE(self,name):     def __init__(self,domain,debug=False):
1729       """       """
1730       @brief return the value of the coefficient name of the general PDE       initializes a new Poisson equation
1731       @param name  
1732         @param domain: domain of the PDE
1733         @type domain: L{Domain<escript.Domain>}
1734         @param debug: if True debug informations are printed.
1735    
1736         """
1737         super(Helmholtz, self).__init__(domain,1,1,debug)
1738         self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1739                            "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1740                            "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1741                            "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1742                            "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1743                            "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1744                            "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1745         self.setSymmetryOn()
1746    
1747       def setValue(self,**coefficients):
1748       """       """
1749       if name == "A" :       sets new values to coefficients
1750           A=self.getCoefficient("A")  
1751           B=self.getCoefficient("B")       @param coefficients: new values assigned to coefficients
1752           C=self.getCoefficient("C")       @keyword omega: value for coefficient M{S{omega}}
1753           if not B.isEmpty() or not C.isEmpty():       @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1754               if A.isEmpty():       @keyword k: value for coefficeint M{k}
1755                   A=self.createNewCoefficient("A")       @type k: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1756               else:       @keyword f: value for right hand side M{f}
1757                   A=A[:]       @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1758               hover2=self.getDomain().getSize()/2.       @keyword alpha: value for right hand side M{S{alpha}}
1759               if self.getNumEquations()>1:       @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1760                  if self.getNumSolutions()>1:       @keyword g: value for right hand side M{g}
1761                     for i in range(self.getNumEquations()):       @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1762                        for k in range(self.getNumSolutions()):       @keyword r: prescribed values M{r} for the solution in constraints.
1763                           alphaByU=self.__getAlphaByU(A[i,:,k,:],hover2)       @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1764                           if not B.isEmpty():                 depending of reduced order is used for the representation of the equation.
1765                               b_sub,f=self.__getUpdateVector(B[i,:,k],hover2,alphaByU)       @keyword q: mask for location of constraints
1766                               for j in range(self.getDim()):       @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1767                                  for l in range(self.getDim()):                 depending of reduced order is used for the representation of the equation.
1768                                     A[i,j,k,l]+=f*b_sub[j]*b_sub[l]       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1769                           if not C.isEmpty():       """
1770                               c_sub,f=self.__getUpdateVector(C[i,k,:],hover2,alphaByU)       super(Helmholtz, self).setValue(**coefficients)
1771                               for j in range(self.getDim()):  
1772                                  for l in range(self.getDim()):     def getCoefficientOfGeneralPDE(self,name):
1773                                     A[i,j,k,l]+=f*c_sub[j]*c_sub[l]       """
1774                  else:         return the value of the coefficient name of the general PDE
1775                     for i in range(self.getNumEquations()):  
1776                        alphaByU=self.__getAlphaByU(A[i,:,:],hover2)       @param name: name of the coefficient requested.
1777                        if not B.isEmpty():       @type name: C{string}
1778                            b_sub,f=self.__getUpdateVector(B[i,:],hover2,alphaByU)       @return: the value of the coefficient  name
1779                            for j in range(self.getDim()):       @rtype: L{Data<escript.Data>}
1780                               for l in range(self.getDim()):       @raise IllegalCoefficient: if name is not one of coefficients
1781                                   A[i,j,l]+=f*b_sub[j]*b_sub[l]                    "A", M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
1782                        if not C.isEmpty():       @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1783                             c_sub,f=self.__getUpdateVector(C[i,:],hover2,alphaByU)       """
1784                             for j in range(self.getDim()):       if name == "A" :
1785                                for l in range(self.getDim()):           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
1786                                   A[i,j,l]+=f*c_sub[j]*c_sub[l]       elif name == "B" :
1787               else:           return escript.Data()
1788                  if self.getNumSolutions()>1:       elif name == "C" :
1789                     for k in range(self.getNumSolutions()):           return escript.Data()
1790                        alphaByU=self.__getAlphaByU(A[:,k,:],hover2)       elif name == "D" :
1791                        if not B.isEmpty():           return self.getCoefficient("omega")
1792                           b_sub,f=self.__getUpdateVector(B[:,k],hover2,alphaByU)       elif name == "X" :
1793                           for j in range(self.getDim()):           return escript.Data()
1794                              for l in range(self.getDim()):       elif name == "Y" :
1795                                     A[j,k,l]+=f*b_sub[j]*b_sub[l]           return self.getCoefficient("f")
1796                        if not C.isEmpty():       elif name == "d" :
1797                           c_sub,f=self.__getUpdateVector(C[k,:],hover2,alphaByU)           return self.getCoefficient("alpha")
1798                           for j in range(self.getDim()):       elif name == "y" :
1799                              for l in range(self.getDim()):           return self.getCoefficient("g")
1800                                 A[j,k,l]+=f*c_sub[j]*c_sub[l]       elif name == "d_contact" :
1801                  else:             return escript.Data()
                    alphaByU=self.__getAlphaByU(A[:,:],hover2)  
                    if not B.isEmpty():  
                        b_sub,f=self.__getUpdateVector(B[:],hover2,alphaByU)  
                        for j in range(self.getDim()):  
                           for l in range(self.getDim()):  
                              A[j,l]+=f*b_sub[j]*b_sub[l]  
                    if not C.isEmpty():  
                       c_sub,f=self.__getUpdateVector(C[:],hover2,alphaByU)  
                       for j in range(self.getDim()):  
                           for l in range(self.getDim()):  
                              A[j,l]+=f*c_sub[j]*c_sub[l]  
          return A  
      elif name == "B" :  
          return self.getCoefficient("B")  
      elif name == "C" :  
          return self.getCoefficient("C")  
      elif name == "D" :  
          return self.getCoefficient("D")  
      elif name == "X" :  
          return self.getCoefficient("X")  
      elif name == "Y" :  
          return self.getCoefficient("Y")  
      elif name == "d" :  
          return self.getCoefficient("d")  
      elif name == "y" :  
          return self.getCoefficient("y")  
      elif name == "d_contact" :  
          return self.getCoefficient("d_contact")  
1802       elif name == "y_contact" :       elif name == "y_contact" :
1803           return self.getCoefficient("y_contact")           return escript.Data()
1804       elif name == "r" :       elif name == "r" :
1805           return self.getCoefficient("r")           return self.getCoefficient("r")
1806       elif name == "q" :       elif name == "q" :
1807           return self.getCoefficient("q")           return self.getCoefficient("q")
1808       else:       else:
1809           raise SystemError,"unknown PDE coefficient %s",name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1810    
1811    class LameEquation(LinearPDE):
 class Poisson(LinearPDE):  
1812     """     """
1813     @brief Class to define a Poisson equstion problem:     Class to define a Lame equation problem:
1814                                                                                                                                                                
1815     class to define a linear PDE of the form     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] }
1816                                                                                                                                                                
1817          -u_{,jj} = f     with natural boundary conditons:
1818                                                                                                                                                                
1819       with boundary conditons:     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] }
1820                                                                                                                                                                
1821          n_j*u_{,j} = 0     and constraints:
1822                                                                                                                                                                
1823      and constraints:     M{u[i]=r[i]} where M{q[i]>0}
1824                                                                                                                                                                
          u=0 where q>0  
                                                                                                                                                               
1825     """     """
1826    
1827     def __init__(self,domain,f=escript.Data(),q=escript.Data()):     def __init__(self,domain,debug=False):
1828         LinearPDE.__init__(self,domain,1,1)        super(LameEquation, self).__init__(domain,\
1829         self.COEFFICIENTS={                                           domain.getDim(),domain.getDim(),debug)
1830         "f"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),        self.COEFFICIENTS={ "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1831         "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}                            "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1832         self.setSymmetryOn()                            "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1833         self.setValue(f,q)                            "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
1834                              "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1835     def setValue(self,f=escript.Data(),q=escript.Data()):                            "r"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1836         self._setValue(f=f,q=q)                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1837          self.setSymmetryOn()
1838     def getCoefficientOfPDE(self,name):  
1839       """     def setValues(self,**coefficients):
1840       @brief return the value of the coefficient name of the general PDE       """
1841       @param name       sets new values to coefficients
1842       """  
1843       if name == "A" :       @param coefficients: new values assigned to coefficients
1844           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))       @keyword lame_mu: value for coefficient M{S{mu}}
1845       elif name == "B" :       @type lame_mu: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1846         @keyword lame_lambda: value for coefficient M{S{lambda}}
1847         @type lame_lambda: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1848         @keyword F: value for internal force M{F}
1849         @type F: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}
1850         @keyword sigma: value for initial stress M{S{sigma}}
1851         @type sigma: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}
1852         @keyword f: value for extrenal force M{f}
1853         @type f: any type that can be casted to L{Vector<escript.Vector>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
1854         @keyword r: prescribed values M{r} for the solution in constraints.
1855         @type r: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1856                   depending of reduced order is used for the representation of the equation.
1857         @keyword q: mask for location of constraints
1858         @type q: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1859                   depending of reduced order is used for the representation of the equation.
1860         @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1861         """
1862         super(LameEquation, self).setValues(**coefficients)
1863    
1864       def getCoefficientOfGeneralPDE(self,name):
1865         """
1866         return the value of the coefficient name of the general PDE
1867    
1868         @param name: name of the coefficient requested.
1869         @type name: C{string}
1870         @return: the value of the coefficient  name
1871         @rtype: L{Data<escript.Data>}
1872         @raise IllegalCoefficient: if name is not one of coefficients
1873                      "A", M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
1874         @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1875         """
1876         if name == "A" :
1877             out =self.createCoefficientOfGeneralPDE("A")
1878             for i in range(self.getDim()):
1879               for j in range(self.getDim()):
1880                 out[i,i,j,j] += self.getCoefficient("lame_lambda")
1881                 out[i,j,j,i] += self.getCoefficient("lame_mu")
1882                 out[i,j,i,j] += self.getCoefficient("lame_mu")
1883             return out
1884         elif name == "B" :
1885           return escript.Data()           return escript.Data()
1886       elif name == "C" :       elif name == "C" :
1887           return escript.Data()           return escript.Data()
1888       elif name == "D" :       elif name == "D" :
1889           return escript.Data()           return escript.Data()
1890       elif name == "X" :       elif name == "X" :
1891             return self.getCoefficient("sigma")
1892         elif name == "Y" :
1893             return self.getCoefficient("F")
1894         elif name == "d" :
1895           return escript.Data()           return escript.Data()
1896       elif name == "Y" :       elif name == "y" :
1897           return self.getCoefficient("f")           return self.getCoefficient("f")
1898       elif name == "d" :       elif name == "d_contact" :
          return escript.Data()  
      elif name == "y" :  
          return escript.Data()  
      elif name == "d_contact" :  
1899           return escript.Data()           return escript.Data()
1900       elif name == "y_contact" :       elif name == "y_contact" :
1901           return escript.Data()           return escript.Data()
1902       elif name == "r" :       elif name == "r" :
1903           return escript.Data()           return self.getCoefficient("r")
1904       elif name == "q" :       elif name == "q" :
1905           return self.getCoefficient("q")           return self.getCoefficient("q")
1906       else:       else:
1907           raise SystemError,"unknown PDE coefficient %s",name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1908    
1909    class AdvectivePDE(LinearPDE):
1910       """
1911       In cases of PDEs dominated by the advection terms M{B} and M{C} against the adevctive terms M{A}
1912       up-winding has been used.  The L{AdvectivePDE} class applies SUPG upwinding to the advective terms.
1913    
1914       In the following we set
1915    
1916       M{Z[j]=C[j]-B[j]}
1917    
1918       or
1919    
1920       M{Z[i,k,l]=C[i,k,l]-B[i,l,k]}
1921    
1922       To measure the dominance of the advective terms over the diffusive term M{A} the
1923       X{Pelclet number} M{P} is used. It is defined as
1924    
1925       M{P=h|Z|/(2|A|)}
1926    
1927       where M{|.|} denotes the L{length<util.length>} of the arument and M{h} is the local cell size
1928       from L{getSize<escript.Domain.getSize>}. Where M{|A|==0} M{P} is M{S{infinity}}.
1929    
1930       From the X{Pelclet number} the stabilization parameters M{S{Xi}} and M{S{Xi}} are calculated:
1931    
1932       M{S{Xi}=S{xi}(P) h/|Z|}
1933    
1934       where M{S{xi}} is a suitable function of the Peclet number.
1935    
1936       In the case of a single PDE the coefficient are up-dated in the following way:
1937             - M{A[i,j] S{<-} A[i,j] + S{Xi} * Z[j] * Z[l]}
1938             - M{B[j] S{<-} B[j] + S{Xi} * C[j] * D}
1939             - M{C[j] S{<-} C[j] + S{Xi} * B[j] * D}
1940             - M{X[j] S{<-} X[j] + S{Xi} * Z[j] * Y}
1941    
1942       Similar for the case of a systems of PDEs:
1943             - M{A[i,j,k,l] S{<-} A[i,j,k,l]+ S{delta}[p,m] * S{Xi} * Z[p,i,j] * Z[m,k,l]}
1944             - M{B[i,j,k] S{<-} B[i,j,k] +  S{delta}[p,m] * S{Xi} * D[p,k] * C[m,i,j]}
1945             - M{C[i,k,l] S{<-} C[i,k,l] +  S{delta}[p,m] * S{Xi} * D[p,k] * B[m,l,i]}
1946             - M{X[i,j] S{<-} X[i,j] + S{delta}[p,m] * S{Xi}  * Y[p] * Z[m,i,j]}
1947    
1948       where M{S{delta}} is L{kronecker}.
1949       Using upwinding in this form, introduces an additonal error which is proprtional to the cell size M{h}
1950       but with the intension to stabilize the solution.
1951    
1952       """
1953       def __init__(self,domain,numEquations=None,numSolutions=None,xi=None,debug=False):
1954          """
1955          creates a linear, steady, second order PDE on a L{Domain<escript.Domain>}
1956    
1957          @param domain: domain of the PDE
1958          @type domain: L{Domain<escript.Domain>}
1959          @param numEquations: number of equations. If numEquations==None the number of equations
1960                               is exracted from the PDE coefficients.
1961          @param numSolutions: number of solution components. If  numSolutions==None the number of solution components
1962                               is exracted from the PDE coefficients.
1963          @param xi: defines a function which returns for any given Preclet number as L{Scalar<escript.Scalar>} object the
1964                     M{S{xi}}-value used to define the stabilization parameters. If equal to None, L{ELMAN_RAMAGE} is used.
1965          @type xi: callable object which returns a L{Scalar<escript.Scalar>} object.
1966          @param debug: if True debug informations are printed.
1967          """
1968          super(AdvectivePDE, self).__init__(domain,\
1969                                             numEquations,numSolutions,debug)
1970          if xi==None:
1971             self.__xi=AdvectivePDE.ELMAN_RAMAGE
1972          else:
1973             self.__xi=xi
1974          self.__Xi=escript.Data()
1975    
1976       def setValue(self,**coefficients):
1977          """
1978          sets new values to coefficients
1979    
1980          @param coefficients: new values assigned to coefficients
1981          @keyword A: value for coefficient A.
1982          @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1983          @keyword B: value for coefficient B
1984          @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1985          @keyword C: value for coefficient C
1986          @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1987          @keyword D: value for coefficient D
1988          @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1989          @keyword X: value for coefficient X
1990          @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1991          @keyword Y: value for coefficient Y
1992          @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1993          @keyword d: value for coefficient d
1994          @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1995          @keyword y: value for coefficient y
1996          @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1997          @keyword d_contact: value for coefficient d_contact
1998          @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1999                           or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
2000          @keyword y_contact: value for coefficient y_contact
2001          @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
2002                           or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
2003          @keyword r: values prescribed to the solution at the locations of constraints
2004          @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2005                   depending of reduced order is used for the solution.
2006          @keyword q: mask for location of constraints
2007          @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2008                   depending of reduced order is used for the representation of the equation.
2009          @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2010    
2011          """
2012          if "A" in coefficients.keys()   or "B" in coefficients.keys() or "C" in coefficients.keys(): self.__Xi=escript.Data()
2013          super(AdvectivePDE, self).setValue(**coefficients)
2014    
2015       def ELMAN_RAMAGE(self,P):
2016         """
2017         Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.
2018         This function uses the method suggested by H.C. Elman and A. Ramage, I{SIAM J. Numer. Anal.}, B{40} (2002)
2019              - M{S{xi}(P)=0} for M{P<1}
2020              - M{S{xi}(P)=(1-1/P)/2} otherwise
2021    
2022         @param P: Preclet number
2023         @type P: L{Scalar<escript.Scalar>}
2024         @return: up-wind weightimg factor
2025         @rtype: L{Scalar<escript.Scalar>}
2026         """
2027         return util.wherePositive(P-1.)*0.5*(1.-1./(P+1.e-15))
2028    
2029       def SIMPLIFIED_BROOK_HUGHES(self,P):
2030         """
2031         Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.
2032         The original methods is
2033    
2034         M{S{xi}(P)=coth(P)-1/P}
2035    
2036         As the evaluation of M{coth} is expensive we are using the approximation:
2037    
2038             - M{S{xi}(P)=P/3} where M{P<3}
2039             - M{S{xi}(P)=1/2} otherwise
2040    
2041         @param P: Preclet number
2042         @type P: L{Scalar<escript.Scalar>}
2043         @return: up-wind weightimg factor
2044         @rtype: L{Scalar<escript.Scalar>}
2045         """
2046         c=util.whereNegative(P-3.)
2047         return P/6.*c+1./2.*(1.-c)
2048    
2049       def HALF(self,P):
2050         """
2051         Predefined function to set value M{1/2} for M{S{xi}}
2052    
2053         @param P: Preclet number
2054         @type P: L{Scalar<escript.Scalar>}
2055         @return: up-wind weightimg factor
2056         @rtype: L{Scalar<escript.Scalar>}
2057         """
2058         return escript.Scalar(0.5,P.getFunctionSpace())
2059    
2060       def __getXi(self):
2061          if self.__Xi.isEmpty():
2062             B=self.getCoefficient("B")
2063             C=self.getCoefficient("C")
2064             A=self.getCoefficient("A")
2065             h=self.getDomain().getSize()
2066             self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
2067             if not C.isEmpty() or not B.isEmpty():
2068                if not C.isEmpty() and not B.isEmpty():
2069                    if self.getNumEquations()>1:
2070                       if self.getNumSolutions()>1:
2071                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2072                          for i in range(self.getNumEquations()):
2073                             for k in range(self.getNumSolutions()):
2074                                for l in range(self.getDim()): flux2+=(C[i,k,l]-B[i,l,k])**2
2075                          length_of_flux=util.sqrt(flux2)
2076                          # flux=C-util.reorderComponents(B,[0,2,1])
2077                       else:
2078                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2079                          for i in range(self.getNumEquations()):
2080                             for l in range(self.getDim()): flux2+=(C[i,l]-B[i,l])**2
2081                          length_of_flux=util.sqrt(flux2)
2082                          # flux=C-B
2083                    else:
2084                       if self.getNumSolutions()>1:
2085                          flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2086                          for k in range(self.getNumSolutions()):
2087                             for l in range(self.getDim()): flux2+=(C[k,l]-B[l,k])**2
2088                          # flux=C-util.reorderComponents(B,[1,0])
2089                          length_of_flux=util.sqrt(flux2)
2090                       else:
2091                          length_of_flux=util.length(C-B)
2092                elif C.isEmpty():
2093                  length_of_flux=util.length(B)
2094                else:
2095                  length_of_flux=util.length(C)
2096                flux_max=util.Lsup(length_of_flux)
2097                if flux_max>0.:
2098                  if A.isEmpty():
2099                      inv_A=1./self.SMALL_TOLERANCE
2100                      peclet_number=escript.Scalar(inv_A,length_of_flux.getFunctionSpace())
2101                      xi=self.__xi(self,peclet_number)
2102                  else:
2103                      # length_of_A=util.inner(flux,util.tensormutiply(A,flux))
2104                      length_of_A=util.length(A)
2105                      A_max=util.Lsup(length_of_A)
2106                      if A_max>0:
2107                           inv_A=1./(length_of_A+A_max*self.SMALL_TOLERANCE)
2108                      else:
2109                           inv_A=1./self.SMALL_TOLERANCE
2110                      peclet_number=length_of_flux*h/2*inv_A
2111                      xi=self.__xi(self,peclet_number)
2112                  self.__Xi=h*xi/(length_of_flux+flux_max*self.SMALL_TOLERANCE)
2113                  self.trace("preclet number = %e"%util.Lsup(peclet_number))
2114                else:
2115                  self.__Xi=escript.Scalar(0.,length_of_flux.getFunctionSpace())
2116          return self.__Xi
2117    
2118    
2119       def getCoefficientOfGeneralPDE(self,name):
2120         """
2121         return the value of the coefficient name of the general PDE
2122    
2123         @param name: name of the coefficient requested.
2124         @type name: C{string}
2125         @return: the value of the coefficient name
2126         @rtype: L{Data<escript.Data>}
2127         @raise IllegalCoefficient: if name is not one of coefficients
2128                      M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
2129         @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2130         """
2131         if not self.getNumEquations() == self.getNumSolutions():
2132              raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."
2133    
2134         if name == "A" :
2135             A=self.getCoefficient("A")
2136             B=self.getCoefficient("B")
2137             C=self.getCoefficient("C")
2138             if B.isEmpty() and C.isEmpty():
2139                Aout=A
2140             else:
2141                if A.isEmpty():
2142                   Aout=self.createCoefficientOfGeneralPDE("A")
2143                else:
2144                   Aout=A[:]
2145                Xi=self.__getXi()
2146                if self.getNumEquations()>1:
2147                    for i in range(self.getNumEquations()):
2148                       for j in range(self.getDim()):
2149                          for k in range(self.getNumSolutions()):
2150                             for l in range(self.getDim()):
2151                                if not C.isEmpty() and not B.isEmpty():
2152                                   # tmp=C-util.reorderComponents(B,[0,2,1])
2153                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(tmp,[1,2,0]),tmp,offset=1)
2154                                   for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*(C[p,i,j]-B[p,j,i])*(C[p,k,l]-B[p,l,k])
2155                                elif C.isEmpty():
2156                                   for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]
2157                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(B,[2,1,0]),util.reorder(B,[0,2,1]),offset=1)
2158                                else:
2159                                   for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
2160                                   # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),C,offset=1)
2161                else:
2162                   if not C.isEmpty() and not B.isEmpty():
2163                       delta=(C-B)
2164                       Aout+=util.outer(Xi*delta,delta)
2165                   elif not B.isEmpty():
2166                       Aout+=util.outer(Xi*B,B)
2167                   elif not C.isEmpty():
2168                       Aout+=util.outer(Xi*C,C)
2169             return Aout
2170         elif name == "B" :
2171             # return self.getCoefficient("B")
2172             B=self.getCoefficient("B")
2173             C=self.getCoefficient("C")
2174             D=self.getCoefficient("D")
2175             if C.isEmpty() or D.isEmpty():
2176                Bout=B
2177             else:
2178                Xi=self.__getXi()
2179                if B.isEmpty():
2180                    Bout=self.createCoefficientOfGeneralPDE("B")
2181                else:
2182                    Bout=B[:]
2183                if self.getNumEquations()>1:
2184                   for k in range(self.getNumSolutions()):
2185                      for p in range(self.getNumEquations()):
2186                         tmp=Xi*D[p,k]
2187                         for i in range(self.getNumEquations()):
2188                            for j in range(self.getDim()):
2189                               Bout[i,j,k]+=tmp*C[p,i,j]
2190                               # Bout=Bout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),D,offset=1)
2191                else:
2192                   Bout+=(Xi*D)*C
2193             return Bout
2194         elif name == "C" :
2195             # return self.getCoefficient("C")
2196             B=self.getCoefficient("B")
2197             C=self.getCoefficient("C")
2198             D=self.getCoefficient("D")
2199             if B.isEmpty() or D.isEmpty():
2200                Cout=C
2201             else:
2202                Xi=self.__getXi()
2203                if C.isEmpty():
2204                    Cout=self.createCoefficientOfGeneralPDE("C")
2205                else:
2206                    Cout=C[:]
2207                if self.getNumEquations()>1:
2208                   for k in range(self.getNumSolutions()):
2209                       for p in range(self.getNumEquations()):
2210                          tmp=Xi*D[p,k]
2211                          for i in range(self.getNumEquations()):
2212                            for l in range(self.getDim()):
2213                                     Cout[i,k,l]+=tmp*B[p,l,i]
2214                                     # Cout=Cout+Xi*B[p,l,i]*D[p,k]
2215                else:
2216                   Cout+=(Xi*D)*B
2217             return Cout
2218         elif name == "D" :
2219             return self.getCoefficient("D")
2220         elif name == "X" :
2221             # return self.getCoefficient("X")
2222             X=self.getCoefficient("X")
2223             Y=self.getCoefficient("Y")
2224             B=self.getCoefficient("B")
2225             C=self.getCoefficient("C")
2226             if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):
2227                Xout=X
2228             else:
2229                if X.isEmpty():
2230                    Xout=self.createCoefficientOfGeneralPDE("X")
2231                else:
2232                    Xout=X[:]
2233                Xi=self.__getXi()
2234                if self.getNumEquations()>1:
2235                     for p in range(self.getNumEquations()):
2236                        tmp=Xi*Y[p]
2237                        for i in range(self.getNumEquations()):
2238                           for j in range(self.getDim()):
2239                              if not C.isEmpty() and not B.isEmpty():
2240                                 Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
2241                                 # Xout=X_out+Xi*util.inner(Y,C-util.reorderComponents(B,[0,2,1]),offset=1)
2242                              elif C.isEmpty():
2243                                 Xout[i,j]-=tmp*B[p,j,i]
2244                                 # Xout=X_out-Xi*util.inner(Y,util.reorderComponents(B,[0,2,1]),offset=1)
2245                              else:
2246                                 Xout[i,j]+=tmp*C[p,i,j]
2247                                 # Xout=X_out+Xi*util.inner(Y,C,offset=1)
2248                else:
2249                  if not C.isEmpty() and not B.isEmpty():
2250                    Xout+=(Xi*Y)*(C-B)
2251                  elif C.isEmpty():
2252                    Xout-=(Xi*Y)*B
2253                  else:
2254                    Xout+=(Xi*Y)*C
2255             return Xout
2256         elif name == "Y" :
2257             return self.getCoefficient("Y")
2258         elif name == "d" :
2259             return self.getCoefficient("d")
2260         elif name == "y" :
2261             return self.getCoefficient("y")
2262         elif name == "d_contact" :
2263             return self.getCoefficient("d_contact")
2264         elif name == "y_contact" :
2265             return self.getCoefficient("y_contact")
2266         elif name == "r" :
2267             return self.getCoefficient("r")
2268         elif name == "q" :
2269             return self.getCoefficient("q")
2270         else:
2271            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2272    
2273    # $Log$
2274    # Revision 1.14  2005/09/22 01:54:57  jgs
2275    # Merge of development branch dev-02 back to main trunk on 2005-09-22
2276    #
2277    # Revision 1.13  2005/09/15 03:44:19  jgs
2278    # Merge of development branch dev-02 back to main trunk on 2005-09-15
2279    #
2280    # Revision 1.12  2005/09/01 03:31:28  jgs
2281    # Merge of development branch dev-02 back to main trunk on 2005-09-01
2282    #
2283    # Revision 1.11  2005/08/23 01:24:28  jgs
2284    # Merge of development branch dev-02 back to main trunk on 2005-08-23
2285    #
2286    # Revision 1.10  2005/08/12 01:45:36  jgs
2287    # erge of development branch dev-02 back to main trunk on 2005-08-12
2288    #
2289    # Revision 1.9.2.17  2005/09/21 07:03:33  matt
2290    # PDECoefficient and LinearPDE are now new style classes (introduced in Python
2291    # 2.2). Classes Poisson, Helmholtz, LameEquation and AdvectivePDE have been
2292    # modified to instead use portable/cooperative "super" calls to extend base
2293    # class methods.
2294    #
2295    # Revision 1.9.2.16  2005/09/16 01:54:37  matt
2296    # Removed redundant if-loop.
2297    #
2298    # Revision 1.9.2.15  2005/09/14 08:09:18  matt
2299    # Added "REDUCED" solution PDECoefficient descriptors for LinearPDEs.
2300    #
2301    # Revision 1.9.2.14  2005/09/07 06:26:16  gross
2302    # the solver from finley are put into the standalone package paso now
2303    #
2304    # Revision 1.9.2.13  2005/08/31 08:45:03  gross
2305    # in the case of lumping no new system is allocated if the constraint is changed.
2306    #
2307    # Revision 1.9.2.12  2005/08/31 07:10:23  gross
2308    # test for Lumping added
2309    #
2310    # Revision 1.9.2.11  2005/08/30 01:53:45  gross
2311    # bug in format fixed.
2312    #
2313    # Revision 1.9.2.10  2005/08/26 07:14:17  gross
2314    # a few more bugs in linearPDE fixed. remaining problem are finley problems
2315    #
2316    # Revision 1.9.2.9  2005/08/26 06:30:45  gross
2317    # fix for reported bug  0000004. test_linearPDE passes a few more tests
2318    #
2319    # Revision 1.9.2.8  2005/08/26 04:30:13  gross
2320    # gneric unit testing for linearPDE
2321    #
2322    # Revision 1.9.2.7  2005/08/25 07:06:50  gross
2323    # linearPDE documentation is parsed now by epydoc. there is still a problem with links into escriptcpp.so
2324    #
2325    # Revision 1.9.2.6  2005/08/24 05:01:24  gross
2326    # problem with resetting the matrix in case of resetting its values to 0 fixed.
2327    #
2328    # Revision 1.9.2.5  2005/08/24 02:03:28  gross
2329    # epydoc mark up partially fixed
2330    #
2331    # Revision 1.9.2.4  2005/08/22 07:11:09  gross
2332    # some problems with LinearPDEs fixed.
2333    #
2334    # Revision 1.9.2.3  2005/08/18 04:48:48  gross
2335    # the methods SetLumping*() are removed. Lumping is set trough setSolverMethod(LinearPDE.LUMPING)
2336    #
2337    # Revision 1.9.2.2  2005/08/18 04:39:32  gross
2338    # the constants have been removed from util.py as they not needed anymore. PDE related constants are accessed through LinearPDE attributes now
2339    #
2340    # Revision 1.9.2.1  2005/07/29 07:10:27  gross
2341    # new functions in util and a new pde type in linearPDEs
2342    #
2343    # Revision 1.1.2.25  2005/07/28 04:21:09  gross
2344    # Lame equation: (linear elastic, isotropic) added
2345    #
2346    # Revision 1.1.2.24  2005/07/22 06:37:11  gross
2347    # some extensions to modellib and linearPDEs
2348    #
2349    # Revision 1.1.2.23  2005/05/13 00:55:20  cochrane
2350    # Fixed up some docstrings.  Moved module-level functions to top of file so
2351    # that epydoc and doxygen can pick them up properly.
2352    #
2353    # Revision 1.1.2.22  2005/05/12 11:41:30  gross
2354    # some basic Models have been added
2355    #
2356    # Revision 1.1.2.21  2005/05/12 07:16:12  cochrane
2357    # Moved ELMAN_RAMAGE, SIMPLIFIED_BROOK_HUGHES, and HALF functions to bottom of
2358    # file so that the AdvectivePDE class is picked up by doxygen.  Some
2359    # reformatting of docstrings.  Addition of code to make equations come out
2360    # as proper LaTeX.
2361    #
2362    # Revision 1.1.2.20  2005/04/15 07:09:08  gross
2363    # some problems with functionspace and linearPDEs fixed.
2364    #
2365    # Revision 1.1.2.19  2005/03/04 05:27:07  gross
2366    # bug in SystemPattern fixed.
2367    #
2368    # Revision 1.1.2.18  2005/02/08 06:16:45  gross
2369    # Bugs in AdvectivePDE fixed, AdvectiveTest is stable but more testing is needed
2370    #
2371    # Revision 1.1.2.17  2005/02/08 05:56:19  gross
2372    # Reference Number handling added
2373    #
2374    # Revision 1.1.2.16  2005/02/07 04:41:28  gross
2375    # some function exposed to python to make mesh merging running
2376    #
2377    # Revision 1.1.2.15  2005/02/03 00:14:44  gross
2378    # timeseries add and ESySParameter.py renames esysXML.py for consistence
2379    #
2380    # Revision 1.1.2.14  2005/02/01 06:44:10  gross
2381    # new implementation of AdvectivePDE which now also updates right hand side. systems of PDEs are still not working
2382    #
2383    # Revision 1.1.2.13  2005/01/25 00:47:07  gross
2384    # updates in the documentation
2385    #
2386    # Revision 1.1.2.12  2005/01/12 01:28:04  matt
2387    # Added createCoefficient method for linearPDEs.
2388    #
2389    # Revision 1.1.2.11  2005/01/11 01:55:34  gross
2390    # a problem in linearPDE class fixed
2391    #
2392    # Revision 1.1.2.10  2005/01/07 01:13:29  gross
2393    # some bugs in linearPDE fixed
2394    #
2395    # Revision 1.1.2.9  2005/01/06 06:24:58  gross
2396    # some bugs in slicing fixed
2397    #
2398    # Revision 1.1.2.8  2005/01/05 04:21:40  gross
2399    # FunctionSpace checking/matchig in slicing added
2400    #
2401    # Revision 1.1.2.7  2004/12/29 10:03:41  gross
2402    # bug in setValue fixed
2403    #
2404    # Revision 1.1.2.6  2004/12/29 05:29:59  gross
2405    # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()
2406    #
2407    # Revision 1.1.2.5  2004/12/29 00:18:41  gross
2408    # AdvectivePDE added
2409    #
2410    # Revision 1.1.2.4  2004/12/24 06:05:41  gross
2411    # some changes in linearPDEs to add AdevectivePDE
2412    #
2413    # Revision 1.1.2.3  2004/12/16 00:12:34  gross
2414    # __init__ of LinearPDE does not accept any coefficient anymore
2415    #
2416    # Revision 1.1.2.2  2004/12/14 03:55:01  jgs
2417    # *** empty log message ***
2418    #
2419    # Revision 1.1.2.1  2004/12/12 22:53:47  gross
2420    # linearPDE has been renamed LinearPDE
2421    #
2422    # Revision 1.1.1.1.2.7  2004/12/07 10:13:08  gross
2423    # GMRES added
2424    #
2425    # Revision 1.1.1.1.2.6  2004/12/07 03:19:50  gross
2426    # options for GMRES and PRES20 added
2427    #
2428    # Revision 1.1.1.1.2.5  2004/12/01 06:25:15  gross
2429    # some small changes
2430    #
2431    # Revision 1.1.1.1.2.4  2004/11/24 01:50:21  gross
2432    # Finley solves 4M unknowns now
2433    #
2434    # Revision 1.1.1.1.2.3  2004/11/15 06:05:26  gross
2435    # poisson solver added
2436    #
2437    # Revision 1.1.1.1.2.2  2004/11/12 06:58:15  gross
2438    # a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry
2439    #
2440    # Revision 1.1.1.1.2.1  2004/10/28 22:59:22  gross
2441    # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed
2442    #
2443    # Revision 1.1.1.1  2004/10/26 06:53:56  jgs
2444    # initial import of project esys2
2445    #
2446    # Revision 1.3.2.3  2004/10/26 06:43:48  jgs
2447    # committing Lutz's and Paul's changes to brach jgs
2448    #
2449    # Revision 1.3.4.1  2004/10/20 05:32:51  cochrane
2450    # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.
2451    #
2452    # Revision 1.3  2004/09/23 00:53:23  jgs
2453    # minor fixes
2454    #
2455    # Revision 1.1  2004/08/28 12:58:06  gross
2456    # SimpleSolve is not running yet: problem with == of functionsspace
2457    #

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

  ViewVC Help
Powered by ViewVC 1.1.26