/[escript]/trunk-mpi-branch/escript/py_src/linearPDEs.py
ViewVC logotype

Diff of /trunk-mpi-branch/escript/py_src/linearPDEs.py

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

trunk/esys2/escript/py_src/linearPDEs.py revision 148 by jgs, Tue Aug 23 01:24:31 2005 UTC trunk-mpi-branch/escript/py_src/linearPDEs.py revision 1302 by ksteube, Mon Sep 17 06:57:51 2007 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
   
 ## @file linearPDEs.py  
   
2  """  """
3  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):  class IllegalCoefficient(ValueError):
38     """     """
39     raised if an illegal coefficient of the general ar particular PDE is requested.     raised if an illegal coefficient of the general ar particular PDE is requested.
40     """     """
41       pass
42    
43  class IllegalCoefficientValue(ValueError):  class IllegalCoefficientValue(ValueError):
44     """     """
45     raised if an incorrect value for a coefficient is used.     raised if an incorrect value for a coefficient is used.
46     """     """
47       pass
48    
49    class IllegalCoefficientFunctionSpace(ValueError):
50       """
51       raised if an incorrect function space for a coefficient is used.
52       """
53    
54  class UndefinedPDEError(ValueError):  class UndefinedPDEError(ValueError):
55     """     """
56     raised if a PDE is not fully defined yet.     raised if a PDE is not fully defined yet.
57     """     """
58       pass
59    
60  def _CompTuple2(t1,t2):  class PDECoefficient(object):
       """  
       Compare two tuples  
     
       @param t1 The first tuple  
       @param t2 The second tuple  
     
       """  
     
       dif=t1[0]+t1[1]-(t2[0]+t2[1])  
       if dif<0: return 1  
       elif dif>0: return -1  
       else: return 0  
     
 class PDECoefficient:  
61      """      """
62      A class for PDE coefficients      A class for describing a PDE coefficient
63    
64        @cvar INTERIOR: indicator that coefficient is defined on the interior of the PDE domain
65        @cvar BOUNDARY: indicator that coefficient is defined on the boundary of the PDE domain
66        @cvar CONTACT: indicator that coefficient is defined on the contact region within the PDE domain
67        @cvar INTERIOR_REDUCED: indicator that coefficient is defined on the interior of the PDE domain using a reduced integration order
68        @cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the boundary of the PDE domain using a reduced integration order
69        @cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact region within the PDE domain using a reduced integration order
70        @cvar SOLUTION: indicator that coefficient is defined trough a solution of the PDE
71        @cvar REDUCED: indicator that coefficient is defined trough a reduced solution of the PDE
72        @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is defined by the number PDE equations
73        @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is defined by the number PDE solutions
74        @cvar BY_DIM: indicator that the dimension of the coefficient shape is defined by the spatial dimension
75        @cvar OPERATOR: indicator that the the coefficient alters the operator of the PDE
76        @cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right hand side of the PDE
77        @cvar BOTH: indicator that the the coefficient alters the operator as well as the right hand side of the PDE
78    
79      """      """
     # identifier for location of Data objects defining COEFFICIENTS  
80      INTERIOR=0      INTERIOR=0
81      BOUNDARY=1      BOUNDARY=1
82      CONTACT=2      CONTACT=2
83      CONTINUOUS=3      SOLUTION=3
84      # identifier in the pattern of COEFFICIENTS:      REDUCED=4
85      # 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
86      # number of unknowns.      BY_SOLUTION=6
87      EQUATION=3      BY_DIM=7
88      SOLUTION=4      OPERATOR=10
89      DIM=5      RIGHTHANDSIDE=11
90      # indicator for what is altered if the coefficient is altered:      BOTH=12
91      OPERATOR=5      INTERIOR_REDUCED=13
92      RIGHTHANDSIDE=6      BOUNDARY_REDUCED=14
93      BOTH=7      CONTACT_REDUCED=15
94      def __init__(self,where,pattern,altering):  
95        def __init__(self, where, pattern, altering):
96         """         """
97         Initialise a PDE Coefficient type         Initialise a PDE Coefficient type
98    
99           @param where: describes where the coefficient lives
100           @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION}, L{REDUCED},
101                               L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED}, L{CONTACT_REDUCED}.
102           @param pattern: describes the shape of the coefficient and how the shape is build for a given
103                  spatial dimension and numbers of equation and solution in then PDE. For instance,
104                  (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) descrbes a rank 3 coefficient which
105                  is instanciated as shape (3,2,2) in case of a three equations and two solution components
106                  on a 2-dimensional domain. In the case of single equation and a single solution component
107                  the shape compoments marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped. In this case
108                  the example would be read as (2,).
109           @type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM}
110           @param altering: indicates what part of the PDE is altered if the coefficiennt is altered
111           @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}
112           @param reduced: indicates if reduced
113           @type reduced: C{bool}
114         """         """
115           super(PDECoefficient, self).__init__()
116         self.what=where         self.what=where
117         self.pattern=pattern         self.pattern=pattern
118         self.altering=altering         self.altering=altering
# Line 74  class PDECoefficient: Line 124  class PDECoefficient:
124         """         """
125         self.value=escript.Data()         self.value=escript.Data()
126    
127      def getFunctionSpace(self,domain):      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
128         """         """
129         defines the FunctionSpace of the coefficient on the domain         defines the L{FunctionSpace<escript.FunctionSpace>} of the coefficient
130    
131         @param domain:         @param domain: domain on which the PDE uses the coefficient
132         """         @type domain: L{Domain<escript.Domain>}
133         if self.what==self.INTERIOR: return escript.Function(domain)         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
134         elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)         @type reducedEquationOrder: C{bool}
135         elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
136         elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)         @type reducedSolutionOrder: C{bool}
137           @return:  L{FunctionSpace<escript.FunctionSpace>} of the coefficient
138           @rtype:  L{FunctionSpace<escript.FunctionSpace>}
139           """
140           if self.what==self.INTERIOR:
141                return escript.Function(domain)
142           elif self.what==self.INTERIOR_REDUCED:
143                return escript.ReducedFunction(domain)
144           elif self.what==self.BOUNDARY:
145                return escript.FunctionOnBoundary(domain)
146           elif self.what==self.BOUNDARY_REDUCED:
147                return escript.ReducedFunctionOnBoundary(domain)
148           elif self.what==self.CONTACT:
149                return escript.FunctionOnContactZero(domain)
150           elif self.what==self.CONTACT_REDUCED:
151                return escript.ReducedFunctionOnContactZero(domain)
152           elif self.what==self.SOLUTION:
153                if reducedEquationOrder and reducedSolutionOrder:
154                    return escript.ReducedSolution(domain)
155                else:
156                    return escript.Solution(domain)
157           elif self.what==self.REDUCED:
158                return escript.ReducedSolution(domain)
159    
160      def getValue(self):      def getValue(self):
161         """         """
162         returns the value of the coefficient:         returns the value of the coefficient
163    
164           @return:  value of the coefficient
165           @rtype:  L{Data<escript.Data>}
166         """         """
167         return self.value         return self.value
168    
169      def setValue(self,domain,numEquations=1,numSolutions=1,newValue=None):      def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
170         """         """
171         set the value of the coefficient to new value         set the value of the coefficient to a new value
172    
173           @param domain: domain on which the PDE uses the coefficient
174           @type domain: L{Domain<escript.Domain>}
175           @param numEquations: number of equations of the PDE
176           @type numEquations: C{int}
177           @param numSolutions: number of components of the PDE solution
178           @type numSolutions: C{int}
179           @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
180           @type reducedEquationOrder: C{bool}
181           @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
182           @type reducedSolutionOrder: C{bool}
183           @param newValue: number of components of the PDE solution
184           @type newValue: any object that can be converted into a L{Data<escript.Data>} object with the appropriate shape and L{FunctionSpace<escript.FunctionSpace>}
185           @raise IllegalCoefficientValue: if the shape of the assigned value does not match the shape of the coefficient
186           @raise IllegalCoefficientFunctionSpace: if unable to interploate value to appropriate function space
187         """         """
188         if newValue==None:         if newValue==None:
189             newValue=escript.Data()             newValue=escript.Data()
190         elif isinstance(newValue,escript.Data):         elif isinstance(newValue,escript.Data):
191             if not newValue.isEmpty():             if not newValue.isEmpty():
192                newValue=escript.Data(newValue,self.getFunctionSpace(domain))                if not newValue.getFunctionSpace() == self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder):
193                    try:
194                      newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
195                    except:
196                      raise IllegalCoefficientFunctionSpace,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)
197         else:         else:
198             newValue=escript.Data(newValue,self.getFunctionSpace(domain))             newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
199         if not newValue.isEmpty():         if not newValue.isEmpty():
200             if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():             if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
201                 raise IllegalCoefficientValue,"Expected shape for coefficient %s is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())                 raise IllegalCoefficientValue,"Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())
202         self.value=newValue         self.value=newValue
203    
204      def isAlteringOperator(self):      def isAlteringOperator(self):
205          """          """
206      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
207    
208            @return:  True if the operator of the PDE is changed when the coefficient is changed
209            @rtype:  C{bool}
210      """      """
211          if self.altering==self.OPERATOR or self.altering==self.BOTH:          if self.altering==self.OPERATOR or self.altering==self.BOTH:
212              return not None              return not None
# Line 118  class PDECoefficient: Line 215  class PDECoefficient:
215    
216      def isAlteringRightHandSide(self):      def isAlteringRightHandSide(self):
217          """          """
218      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
219    
220        @rtype:  C{bool}
221            @return:  True if the right hand side of the PDE is changed when the coefficient is changed
222      """      """
223          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
224              return not None              return not None
# Line 127  class PDECoefficient: Line 227  class PDECoefficient:
227    
228      def estimateNumEquationsAndNumSolutions(self,domain,shape=()):      def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
229         """         """
230         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
231    
232         @param shape:         @param domain: domain on which the PDE uses the coefficient
233         @param dim:         @type domain: L{Domain<escript.Domain>}
234           @param shape: suggested shape of the coefficient
235           @type shape: C{tuple} of C{int} values
236           @return: the number of equations and number of solutions of the PDE is the coefficient has shape s.
237                     If no appropriate numbers could be identified, C{None} is returned
238           @rtype: C{tuple} of two C{int} values or C{None}
239         """         """
240         dim=domain.getDim()         dim=domain.getDim()
241         if len(shape)>0:         if len(shape)>0:
# Line 138  class PDECoefficient: Line 243  class PDECoefficient:
243         else:         else:
244             num=1             num=1
245         search=[]         search=[]
246         for u in range(num):         if self.definesNumEquation() and self.definesNumSolutions():
247            for e in range(num):            for u in range(num):
248               search.append((e,u))               for e in range(num):
249         search.sort(_CompTuple2)                  search.append((e,u))
250         for item in search:            search.sort(self.__CompTuple2)
251              for item in search:
252               s=self.getShape(domain,item[0],item[1])               s=self.getShape(domain,item[0],item[1])
253               if len(s)==0 and len(shape)==0:               if len(s)==0 and len(shape)==0:
254                   return (1,1)                   return (1,1)
255               else:               else:
256                   if s==shape: return item                   if s==shape: return item
257           elif self.definesNumEquation():
258              for e in range(num,0,-1):
259                 s=self.getShape(domain,e,0)
260                 if len(s)==0 and len(shape)==0:
261                     return (1,None)
262                 else:
263                     if s==shape: return (e,None)
264    
265           elif self.definesNumSolutions():
266              for u in range(num,0,-1):
267                 s=self.getShape(domain,0,u)
268                 if len(s)==0 and len(shape)==0:
269                     return (None,1)
270                 else:
271                     if s==shape: return (None,u)
272         return None         return None
273        def definesNumSolutions(self):
274           """
275           checks if the coefficient allows to estimate the number of solution components
276    
277           @return: True if the coefficient allows an estimate of the number of solution components
278           @rtype: C{bool}
279           """
280           for i in self.pattern:
281                 if i==self.BY_SOLUTION: return True
282           return False
283    
284        def definesNumEquation(self):
285           """
286           checks if the coefficient allows to estimate the number of equations
287    
288           @return: True if the coefficient allows an estimate of the number of equations
289           @rtype: C{bool}
290           """
291           for i in self.pattern:
292                 if i==self.BY_EQUATION: return True
293           return False
294    
295        def __CompTuple2(self,t1,t2):
296          """
297          Compare two tuples of possible number of equations and number of solutions
298    
299          @param t1: The first tuple
300          @param t2: The second tuple
301    
302          """
303    
304          dif=t1[0]+t1[1]-(t2[0]+t2[1])
305          if dif<0: return 1
306          elif dif>0: return -1
307          else: return 0
308    
309      def getShape(self,domain,numEquations=1,numSolutions=1):      def getShape(self,domain,numEquations=1,numSolutions=1):
310          """         """
311      builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim         builds the required shape of the coefficient
312    
313      @param e:         @param domain: domain on which the PDE uses the coefficient
314      @param u:         @type domain: L{Domain<escript.Domain>}
315      @param dim:         @param numEquations: number of equations of the PDE
316      """         @type numEquations: C{int}
317          dim=domain.getDim()         @param numSolutions: number of components of the PDE solution
318          s=()         @type numSolutions: C{int}
319          for i in self.pattern:         @return: shape of the coefficient
320               if i==self.EQUATION:         @rtype: C{tuple} of C{int} values
321           """
322           dim=domain.getDim()
323           s=()
324           for i in self.pattern:
325                 if i==self.BY_EQUATION:
326                  if numEquations>1: s=s+(numEquations,)                  if numEquations>1: s=s+(numEquations,)
327               elif i==self.SOLUTION:               elif i==self.BY_SOLUTION:
328                  if numSolutions>1: s=s+(numSolutions,)                  if numSolutions>1: s=s+(numSolutions,)
329               else:               else:
330                  s=s+(dim,)                  s=s+(dim,)
331          return s         return s
332    
333  class LinearPDE:  class LinearPDE(object):
334     """     """
335     Class to define a linear PDE of the form     This class is used to define a general linear, steady, second order PDE
336       for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
337    
338     \f[     For a single PDE with a solution with a single component the linear PDE is defined in the following form:
      -(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i  
    \f]  
339    
340     with boundary conditons:     M{-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)=-grad(X+X_reduced)[j,j]+(Y+Y_reduced)}
341    
    \f[  
    n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i  
    \f]  
342    
343     and contact conditions     where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,
344       ie. summation over indexes appearing twice in a term of a sum is performed, is used.
345       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
346       L{Function<escript.Function>} and the coefficients M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through L{Data<escript.Data>} objects in the
347       L{ReducedFunction<escript.ReducedFunction>}. It is also allowd to use objects that can be converted into
348       such L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B_reduced}, M{C_reduced}, M{X_reduced}
349       M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and M{D}, M{D_reduced} and M{Y_reduced} are scalar.
350    
351     \f[     The following natural boundary conditions are considered:
    n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i  
    \f]  
352    
353     and constraints:     M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u)+(d+d_reduced)*u=n[j]*(X[j]+X_reduced[j])+y}
354    
355       where M{n} is the outer normal field. Notice that the coefficients M{A}, M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the PDE. The coefficients M{d} and M{y} and are each a scalar in the L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients M{d_reduced} and M{y_reduced} and are each a scalar in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
356    
357    
358       Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
359    
360       M{u=r}  where M{q>0}
361    
362       M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
363       The constraints override any other condition set by the PDE or the boundary condition.
364    
365       The PDE is symmetrical if
366    
367       M{A[i,j]=A[j,i]}  and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]}  and M{B_reduced[j]=C_reduced[j]}
368    
369       For a system of PDEs and a solution with several components the PDE has the form
370    
371       M{-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k] =-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] }
372    
373       M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and M{C_reduced} are each of rank three, M{D}, M{D_reduced}, M{X_reduced} and M{X} are each a rank two and M{Y} and M{Y_reduced} are of rank one.
374       The natural boundary conditions take the form:
375    
376       M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])+(d[i,k]+d_reduced[i,k])*u[k]=n[j]*(X[i,j]+X_reduced[i,j])+y[i]+y_reduced[i]}
377    
378    
379       The coefficient M{d} is a rank two and M{y} is a rank one both in the L{FunctionOnBoundary<escript.FunctionOnBoundary>}. Constraints take the form and the coefficients M{d_reduced} is a rank two and M{y_reduced} is a rank one both in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
380    
381       Constraints take the form
382    
383       M{u[i]=r[i]}  where  M{q[i]>0}
384    
385       M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
386    
387       The system of PDEs is symmetrical if
388    
389            - M{A[i,j,k,l]=A[k,l,i,j]}
390            - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
391            - M{B[i,j,k]=C[k,i,j]}
392            - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
393            - M{D[i,k]=D[i,k]}
394            - M{D_reduced[i,k]=D_reduced[i,k]}
395            - M{d[i,k]=d[k,i]}
396            - M{d_reduced[i,k]=d_reduced[k,i]}
397    
398       L{LinearPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
399       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
400       defined as
401    
402       M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]-X[i,j]-X_reduced[i,j]}
403    
404       For the case of single solution component and single PDE M{J} is defined
405    
406       M{J_{j}=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}
407    
408       In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
409       calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of PDEs
410       the contact condition takes the form
411    
412       M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]}
413    
414       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
415       of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calcualted by
416       L{jump<util.jump>}.
417       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>}.
418       The coefficient M{d_contact_reduced} is a rank two and M{y_contact_reduced} is a rank one both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
419       In case of a single PDE and a single component solution the contact condition takes the form
420    
421       M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
422    
423     \f[     In this case the coefficient M{d_contact} and M{y_contact} are each scalar both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}
424     u_i=r_i \quad \mathrm{where} \quad q_i>0  
425     \f]     @cvar DEFAULT: The default method used to solve the system of linear equations
426       @cvar DIRECT: The direct solver based on LDU factorization
427       @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
428       @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
429       @cvar CR: The conjugate residual method
430       @cvar CGS: The conjugate gardient square method
431       @cvar BICGSTAB: The stabilized BiConjugate Gradient method.
432       @cvar SSOR: The symmetric overrealaxtion method
433       @cvar ILU0: The incomplete LU factorization preconditioner  with no fill in
434       @cvar ILUT: The incomplete LU factorization preconditioner with will in
435       @cvar JACOBI: The Jacobi preconditioner
436       @cvar GMRES: The Gram-Schmidt minimum residual method
437       @cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals
438       @cvar LUMPING: Matrix lumping.
439       @cvar NO_REORDERING: No matrix reordering allowed
440       @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
441       @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
442       @cvar PASO: PASO solver package
443       @cvar SCSL: SGI SCSL solver library
444       @cvar MKL: Intel's MKL solver library
445       @cvar UMFPACK: the UMFPACK library
446       @cvar TRILINOS: the TRILINOS parallel solver class library from Sandia Natl Labs
447       @cvar ITERATIVE: The default iterative solver
448       @cvar AMG: algebraic multi grid
449       @cvar RILU: recursive ILU
450    
451     """     """
452     TOL=1.e-13     DEFAULT= 0
453     # solver methods     DIRECT= 1
454     UNKNOWN=-1     CHOLEVSKY= 2
455     DEFAULT_METHOD=0     PCG= 3
456     DIRECT=1     CR= 4
457     CHOLEVSKY=2     CGS= 5
458     PCG=3     BICGSTAB= 6
459     CR=4     SSOR= 7
460     CGS=5     ILU0= 8
461     BICGSTAB=6     ILUT= 9
462     SSOR=7     JACOBI= 10
463     ILU0=8     GMRES= 11
464     ILUT=9     PRES20= 12
465     JACOBI=10     LUMPING= 13
466     GMRES=11     NO_REORDERING= 17
467     PRES20=12     MINIMUM_FILL_IN= 18
468     LUMPING=13     NESTED_DISSECTION= 19
469     # matrix reordering:     SCSL= 14
470     NO_REORDERING=0     MKL= 15
471     MINIMUM_FILL_IN=1     UMFPACK= 16
472     NESTED_DISSECTION=2     ITERATIVE= 20
473     # important keys in the dictonary used to hand over options to the solver library:     PASO= 21
474     METHOD_KEY="method"     AMG= 22
475     SYMMETRY_KEY="symmetric"     RILU = 23
476     TOLERANCE_KEY="tolerance"     TRILINOS = 24
477    
478       SMALL_TOLERANCE=1.e-13
479       __PACKAGE_KEY="package"
480       __METHOD_KEY="method"
481       __SYMMETRY_KEY="symmetric"
482       __TOLERANCE_KEY="tolerance"
483       __PRECONDITIONER_KEY="preconditioner"
484    
485    
486     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
# Line 228  class LinearPDE: Line 488  class LinearPDE:
488       initializes a new linear PDE       initializes a new linear PDE
489    
490       @param domain: domain of the PDE       @param domain: domain of the PDE
491       @type domain: L{Domain}       @type domain: L{Domain<escript.Domain>}
492       @param numEquations: number of equations. If numEquations==None the number of equations       @param numEquations: number of equations. If numEquations==None the number of equations
493                            is exracted from the PDE coefficients.                            is exracted from the PDE coefficients.
494       @param numSolutions: number of solution components. If  numSolutions==None the number of solution components       @param numSolutions: number of solution components. If  numSolutions==None the number of solution components
495                            is exracted from the PDE coefficients.                            is exracted from the PDE coefficients.
496       @param debug: if True debug informations are printed.       @param debug: if True debug informations are printed.
497    
   
498       """       """
499         super(LinearPDE, self).__init__()
500       #       #
501       #   the coefficients of the general PDE:       #   the coefficients of the general PDE:
502       #       #
503       self.__COEFFICIENTS_OF_GENEARL_PDE={       self.__COEFFICIENTS_OF_GENEARL_PDE={
504         "A"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),         "A"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
505         "B"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),         "B"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
506         "C"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),         "C"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
507         "D"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),         "D"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
508         "X"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),         "X"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
509         "Y"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),         "Y"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
510         "d"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),         "d"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
511         "y"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),         "y"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
512         "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),         "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
513         "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),         "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
514         "r"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),         "A_reduced"         : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
515         "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)}         "B_reduced"         : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
516           "C_reduced"         : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
517           "D_reduced"         : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
518           "X_reduced"         : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
519           "Y_reduced"         : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
520           "d_reduced"         : PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
521           "y_reduced"         : PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
522           "d_contact_reduced" : PDECoefficient(PDECoefficient.CONTACT_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
523           "y_contact_reduced" : PDECoefficient(PDECoefficient.CONTACT_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
524           "r"         : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.RIGHTHANDSIDE),
525           "q"         : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.BOTH)}
526    
527       # COEFFICIENTS can be overwritten by subclasses:       # COEFFICIENTS can be overwritten by subclasses:
528       self.COEFFICIENTS=self.__COEFFICIENTS_OF_GENEARL_PDE       self.COEFFICIENTS=self.__COEFFICIENTS_OF_GENEARL_PDE
529         self.__altered_coefficients=False
530       # initialize attributes       # initialize attributes
531       self.__debug=debug       self.__debug=debug
532       self.__domain=domain       self.__domain=domain
# Line 264  class LinearPDE: Line 535  class LinearPDE:
535       self.__resetSystem()       self.__resetSystem()
536    
537       # set some default values:       # set some default values:
538       self.__homogeneous_constraint=True       self.__reduce_equation_order=False
539       self.__row_function_space=escript.Solution(self.__domain)       self.__reduce_solution_order=False
      self.__column_function_space=escript.Solution(self.__domain)  
540       self.__tolerance=1.e-8       self.__tolerance=1.e-8
541       self.__solver_method=self.DEFAULT_METHOD       self.__solver_method=self.DEFAULT
542       self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT_METHOD,False)       self.__solver_package=self.DEFAULT
543         self.__preconditioner=self.DEFAULT
544         self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)
545       self.__sym=False       self.__sym=False
546    
547       self.resetCoefficients()       self.resetCoefficients()
# Line 278  class LinearPDE: Line 550  class LinearPDE:
550     #    general stuff:     #    general stuff:
551     # =============================================================================     # =============================================================================
552     def __str__(self):     def __str__(self):
553         return "<LinearPDE %d>"%id(self)       """
554         returns string representation of the PDE
555    
556         @return: a simple representation of the PDE
557         @rtype: C{str}
558         """
559         return "<LinearPDE %d>"%id(self)
560     # =============================================================================     # =============================================================================
561     #    debug :     #    debug :
562     # =============================================================================     # =============================================================================
563     def setDebugOn(self):     def setDebugOn(self):
564       """       """
565       switches on debugging       switches on debugging
566       """       """
567       self.__debug=not None       self.__debug=not None
# Line 296  class LinearPDE: Line 574  class LinearPDE:
574    
575     def trace(self,text):     def trace(self,text):
576       """       """
577       print the text message if debugging is swiched on.       print the text message if debugging is swiched on.
578         @param text: message
579       @param name: name of the coefficient enquired.       @type text: C{string}
      @type name: C{string}  
580       """       """
581       if self.__debug: print "%s: %s"%(str(self),text)       if self.__debug: print "%s: %s"%(str(self),text)
582    
# Line 309  class LinearPDE: Line 586  class LinearPDE:
586     def getDomain(self):     def getDomain(self):
587       """       """
588       returns the domain of the PDE       returns the domain of the PDE
       
      @return : the domain of the PDE  
      @rtype : C{Domain}  
589    
590         @return: the domain of the PDE
591         @rtype: L{Domain<escript.Domain>}
592       """       """
593       return self.__domain       return self.__domain
594    
# Line 320  class LinearPDE: Line 596  class LinearPDE:
596       """       """
597       returns the spatial dimension of the PDE       returns the spatial dimension of the PDE
598    
599       @return : the spatial dimension of the PDE domain       @return: the spatial dimension of the PDE domain
600       @rtype : C{int}       @rtype: C{int}
601       """       """
602       return self.getDomain().getDim()       return self.getDomain().getDim()
603    
# Line 329  class LinearPDE: Line 605  class LinearPDE:
605       """       """
606       returns the number of equations       returns the number of equations
607    
608       @return : the number of equations       @return: the number of equations
609       @rtype : C{int}       @rtype: C{int}
610       @raise UndefinedPDEError: if the number of equations is not be specified yet.       @raise UndefinedPDEError: if the number of equations is not be specified yet.
611       """       """
612       if self.__numEquations==None:       if self.__numEquations==None:
# Line 342  class LinearPDE: Line 618  class LinearPDE:
618       """       """
619       returns the number of unknowns       returns the number of unknowns
620    
621       @return : the number of unknowns       @return: the number of unknowns
622       @rtype : C{int}       @rtype: C{int}
623       @raise UndefinedPDEError: if the number of unknowns is not be specified yet.       @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
624       """       """
625       if self.__numSolutions==None:       if self.__numSolutions==None:
# Line 351  class LinearPDE: Line 627  class LinearPDE:
627       else:       else:
628          return self.__numSolutions          return self.__numSolutions
629    
630       def reduceEquationOrder(self):
631         """
632         return status for order reduction for equation
633    
634         @return: return True is reduced interpolation order is used for the represenation of the equation
635         @rtype: L{bool}
636         """
637         return self.__reduce_equation_order
638    
639       def reduceSolutionOrder(self):
640         """
641         return status for order reduction for the solution
642    
643         @return: return True is reduced interpolation order is used for the represenation of the solution
644         @rtype: L{bool}
645         """
646         return self.__reduce_solution_order
647    
648     def getFunctionSpaceForEquation(self):     def getFunctionSpaceForEquation(self):
649       """       """
650       returns the L{escript.FunctionSpace} used to discretize the equation       returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation
       
      @return : representation space of equation  
      @rtype : L{escript.FunctionSpace}  
651    
652         @return: representation space of equation
653         @rtype: L{FunctionSpace<escript.FunctionSpace>}
654       """       """
655       return self.__row_function_space       if self.reduceEquationOrder():
656             return escript.ReducedSolution(self.getDomain())
657         else:
658             return escript.Solution(self.getDomain())
659    
660     def getFunctionSpaceForSolution(self):     def getFunctionSpaceForSolution(self):
661       """       """
662       returns the L{escript.FunctionSpace} used to represent the solution       returns the L{FunctionSpace<escript.FunctionSpace>} used to represent the solution
       
      @return : representation space of solution  
      @rtype : L{escript.FunctionSpace}  
663    
664         @return: representation space of solution
665         @rtype: L{FunctionSpace<escript.FunctionSpace>}
666       """       """
667       return self.__column_function_space       if self.reduceSolutionOrder():
668             return escript.ReducedSolution(self.getDomain())
669         else:
670             return escript.Solution(self.getDomain())
671    
672    
673     def getOperator(self):     def getOperator(self):
674       """       """
675       provides access to the operator of the PDE       provides access to the operator of the PDE
676    
677       @return : the operator of the PDE       @return: the operator of the PDE
678       @rtype : L{Operator}       @rtype: L{Operator<escript.Operator>}
679       """       """
680       m=self.getSystem()[0]       m=self.getSystem()[0]
681       if self.isUsingLumping():       if self.isUsingLumping():
# Line 388  class LinearPDE: Line 686  class LinearPDE:
686     def getRightHandSide(self):     def getRightHandSide(self):
687       """       """
688       provides access to the right hand side of the PDE       provides access to the right hand side of the PDE
689         @return: the right hand side of the PDE
690       @return : the right hand side of the PDE       @rtype: L{Data<escript.Data>}
      @rtype : L{escript.Data}  
691       """       """
692       r=self.getSystem()[1]       r=self.getSystem()[1]
693       if self.isUsingLumping():       if self.isUsingLumping():
# Line 404  class LinearPDE: Line 701  class LinearPDE:
701    
702       @param u: argument of the operator. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}       @param u: argument of the operator. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
703                 the current solution is used.                 the current solution is used.
704       @type u: L{escript.Data} or None       @type u: L{Data<escript.Data>} or None
705       @return : image of u       @return: image of u
706       @rtype : L{escript.Data}       @rtype: L{Data<escript.Data>}
707       """       """
708       if u==None:       if u==None:
709            return self.getOperator()*self.getSolution()          return self.getOperator()*self.getSolution()
710       else:       else:
711          self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())          return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
712    
713     def getResidual(self,u=None):     def getResidual(self,u=None):
714       """       """
# Line 419  class LinearPDE: Line 716  class LinearPDE:
716    
717       @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}       @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
718                 the current solution is used.                 the current solution is used.
719       @type u: L{escript.Data} or None       @type u: L{Data<escript.Data>} or None
720       @return : residual of u       @return: residual of u
721       @rtype : L{escript.Data}       @rtype: L{Data<escript.Data>}
722       """       """
723       return self.applyOperator(u)-self.getRightHandSide()       return self.applyOperator(u)-self.getRightHandSide()
724    
# Line 429  class LinearPDE: Line 726  class LinearPDE:
726        """        """
727        test the PDE for symmetry.        test the PDE for symmetry.
728    
729          @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
730       @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.        @type verbose: C{bool}
731       @type verbose: C{bool}        @return:  True if the PDE is symmetric.
732       @return:  True if the PDE is symmetric.        @rtype: L{Data<escript.Data>}
      @rtype : C{escript.Data}  
   
733        @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.        @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
734        """        """
735        verbose=verbose or self.debug()        verbose=verbose or self.__debug
736        out=True        out=True
737        if self.getNumSolutions()!=self.getNumEquations():        if self.getNumSolutions()!=self.getNumEquations():
738           if verbose: print "non-symmetric PDE because of different number of equations and solutions"           if verbose: print "non-symmetric PDE because of different number of equations and solutions"
# Line 445  class LinearPDE: Line 740  class LinearPDE:
740        else:        else:
741           A=self.getCoefficientOfGeneralPDE("A")           A=self.getCoefficientOfGeneralPDE("A")
742           if not A.isEmpty():           if not A.isEmpty():
743              tol=util.Lsup(A)*self.TOL              tol=util.Lsup(A)*self.SMALL_TOLERANCE
744              if self.getNumSolutions()>1:              if self.getNumSolutions()>1:
745                 for i in range(self.getNumEquations()):                 for i in range(self.getNumEquations()):
746                    for j in range(self.getDim()):                    for j in range(self.getDim()):
# Line 469  class LinearPDE: Line 764  class LinearPDE:
764              if verbose: print "non-symmetric PDE because C is not present but B is"              if verbose: print "non-symmetric PDE because C is not present but B is"
765              out=False              out=False
766           elif not B.isEmpty() and not C.isEmpty():           elif not B.isEmpty() and not C.isEmpty():
767              tol=(util.Lsup(B)+util.Lsup(C))*self.TOL/2.              tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
768              if self.getNumSolutions()>1:              if self.getNumSolutions()>1:
769                 for i in range(self.getNumEquations()):                 for i in range(self.getNumEquations()):
770                     for j in range(self.getDim()):                     for j in range(self.getDim()):
# Line 485  class LinearPDE: Line 780  class LinearPDE:
780           if self.getNumSolutions()>1:           if self.getNumSolutions()>1:
781             D=self.getCoefficientOfGeneralPDE("D")             D=self.getCoefficientOfGeneralPDE("D")
782             if not D.isEmpty():             if not D.isEmpty():
783               tol=util.Lsup(D)*self.TOL               tol=util.Lsup(D)*self.SMALL_TOLERANCE
784               for i in range(self.getNumEquations()):               for i in range(self.getNumEquations()):
785                  for k in range(self.getNumSolutions()):                  for k in range(self.getNumSolutions()):
786                    if util.Lsup(D[i,k]-D[k,i])>tol:                    if util.Lsup(D[i,k]-D[k,i])>tol:
787                        if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)                        if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)
788                        out=False                        out=False
789               d=self.getCoefficientOfGeneralPDE("d")
790               if not d.isEmpty():
791                 tol=util.Lsup(d)*self.SMALL_TOLERANCE
792                 for i in range(self.getNumEquations()):
793                    for k in range(self.getNumSolutions()):
794                      if util.Lsup(d[i,k]-d[k,i])>tol:
795                          if verbose: print "non-symmetric PDE because d[%d,%d]!=d[%d,%d]"%(i,k,k,i)
796                          out=False
797               d_contact=self.getCoefficientOfGeneralPDE("d_contact")
798               if not d_contact.isEmpty():
799                 tol=util.Lsup(d_contact)*self.SMALL_TOLERANCE
800                 for i in range(self.getNumEquations()):
801                    for k in range(self.getNumSolutions()):
802                      if util.Lsup(d_contact[i,k]-d_contact[k,i])>tol:
803                          if verbose: print "non-symmetric PDE because d_contact[%d,%d]!=d_contact[%d,%d]"%(i,k,k,i)
804                          out=False
805             # and now the reduced coefficients
806             A_reduced=self.getCoefficientOfGeneralPDE("A_reduced")
807             if not A_reduced.isEmpty():
808                tol=util.Lsup(A_reduced)*self.SMALL_TOLERANCE
809                if self.getNumSolutions()>1:
810                   for i in range(self.getNumEquations()):
811                      for j in range(self.getDim()):
812                         for k in range(self.getNumSolutions()):
813                            for l in range(self.getDim()):
814                                if util.Lsup(A_reduced[i,j,k,l]-A_reduced[k,l,i,j])>tol:
815                                   if verbose: print "non-symmetric PDE because A_reduced[%d,%d,%d,%d]!=A_reduced[%d,%d,%d,%d]"%(i,j,k,l,k,l,i,j)
816                                   out=False
817                else:
818                   for j in range(self.getDim()):
819                      for l in range(self.getDim()):
820                         if util.Lsup(A_reduced[j,l]-A_reduced[l,j])>tol:
821                            if verbose: print "non-symmetric PDE because A_reduced[%d,%d]!=A_reduced[%d,%d]"%(j,l,l,j)
822                            out=False
823             B_reduced=self.getCoefficientOfGeneralPDE("B_reduced")
824             C_reduced=self.getCoefficientOfGeneralPDE("C_reduced")
825             if B_reduced.isEmpty() and not C_reduced.isEmpty():
826                if verbose: print "non-symmetric PDE because B_reduced is not present but C_reduced is"
827                out=False
828             elif not B_reduced.isEmpty() and C_reduced.isEmpty():
829                if verbose: print "non-symmetric PDE because C_reduced is not present but B_reduced is"
830                out=False
831             elif not B_reduced.isEmpty() and not C_reduced.isEmpty():
832                tol=(util.Lsup(B_reduced)+util.Lsup(C_reduced))*self.SMALL_TOLERANCE/2.
833                if self.getNumSolutions()>1:
834                   for i in range(self.getNumEquations()):
835                       for j in range(self.getDim()):
836                          for k in range(self.getNumSolutions()):
837                             if util.Lsup(B_reduced[i,j,k]-C_reduced[k,i,j])>tol:
838                                  if verbose: print "non-symmetric PDE because B_reduced[%d,%d,%d]!=C_reduced[%d,%d,%d]"%(i,j,k,k,i,j)
839                                  out=False
840                else:
841                   for j in range(self.getDim()):
842                      if util.Lsup(B_reduced[j]-C_reduced[j])>tol:
843                         if verbose: print "non-symmetric PDE because B_reduced[%d]!=C_reduced[%d]"%(j,j)
844                         out=False
845             if self.getNumSolutions()>1:
846               D_reduced=self.getCoefficientOfGeneralPDE("D_reduced")
847               if not D_reduced.isEmpty():
848                 tol=util.Lsup(D_reduced)*self.SMALL_TOLERANCE
849                 for i in range(self.getNumEquations()):
850                    for k in range(self.getNumSolutions()):
851                      if util.Lsup(D_reduced[i,k]-D_reduced[k,i])>tol:
852                          if verbose: print "non-symmetric PDE because D_reduced[%d,%d]!=D_reduced[%d,%d]"%(i,k,k,i)
853                          out=False
854               d_reduced=self.getCoefficientOfGeneralPDE("d_reduced")
855               if not d_reduced.isEmpty():
856                 tol=util.Lsup(d_reduced)*self.SMALL_TOLERANCE
857                 for i in range(self.getNumEquations()):
858                    for k in range(self.getNumSolutions()):
859                      if util.Lsup(d_reduced[i,k]-d_reduced[k,i])>tol:
860                          if verbose: print "non-symmetric PDE because d_reduced[%d,%d]!=d_reduced[%d,%d]"%(i,k,k,i)
861                          out=False
862               d_contact_reduced=self.getCoefficientOfGeneralPDE("d_contact_reduced")
863               if not d_contact_reduced.isEmpty():
864                 tol=util.Lsup(d_contact_reduced)*self.SMALL_TOLERANCE
865                 for i in range(self.getNumEquations()):
866                    for k in range(self.getNumSolutions()):
867                      if util.Lsup(d_contact_reduced[i,k]-d_contact_reduced[k,i])>tol:
868                          if verbose: print "non-symmetric PDE because d_contact_reduced[%d,%d]!=d_contact_reduced[%d,%d]"%(i,k,k,i)
869                          out=False
870        return out        return out
871    
872     def getSolution(self,**options):     def getSolution(self,**options):
873         """         """
874         returns the solution of the PDE. If the solution is not valid the PDE is solved.         returns the solution of the PDE. If the solution is not valid the PDE is solved.
875    
876         @return: the solution         @return: the solution
877         @rtype: L{escript.Data}         @rtype: L{Data<escript.Data>}
878         @param options: solver options         @param options: solver options
879         @keyword verbose:         @keyword verbose: True to get some information during PDE solution
880         @keyword reordering: reordering scheme to be used during elimination         @type verbose: C{bool}
881         @keyword preconditioner: preconditioner method to be used         @keyword reordering: reordering scheme to be used during elimination. Allowed values are
882                                L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}
883         @keyword iter_max: maximum number of iteration steps allowed.         @keyword iter_max: maximum number of iteration steps allowed.
884         @keyword drop_tolerance:         @keyword drop_tolerance: threshold for drupping in L{ILUT}
885         @keyword drop_storage:         @keyword drop_storage: maximum of allowed memory in L{ILUT}
886         @keyword truncation:         @keyword truncation: maximum number of residuals in L{GMRES}
887         @keyword restart:         @keyword restart: restart cycle length in L{GMRES}
888         """         """
889         if not self.__solution_isValid:         if not self.__solution_isValid:
890            mat,f=self.getSystem()            mat,f=self.getSystem()
891            if self.isUsingLumping():            if self.isUsingLumping():
892               self.__solution=self.copyConstraint(f*mat)               self.__solution=self.copyConstraint(f*mat)
893            else:            else:
894               options[self.TOLERANCE_KEY]=self.getTolerance()               options[self.__TOLERANCE_KEY]=self.getTolerance()
895               options[self.METHOD_KEY]=self.getSolverMethod()               options[self.__METHOD_KEY]=self.getSolverMethod()[0]
896               options[self.SYMMETRY_KEY]=self.isSymmetric()               options[self.__PRECONDITIONER_KEY]=self.getSolverMethod()[1]
897                 options[self.__PACKAGE_KEY]=self.getSolverPackage()
898                 options[self.__SYMMETRY_KEY]=self.isSymmetric()
899               self.trace("PDE is resolved.")               self.trace("PDE is resolved.")
900               self.trace("solver options: %s"%str(options))               self.trace("solver options: %s"%str(options))
901               self.__solution=mat.solve(f,options)               self.__solution=mat.solve(f,options)
# Line 526  class LinearPDE: Line 904  class LinearPDE:
904    
905     def getFlux(self,u=None):     def getFlux(self,u=None):
906       """       """
907       returns the flux J_ij for a given u       returns the flux M{J} for a given M{u}
908    
909         \f[       M{J[i,j]=(A[i,j,k,l]+A_reduced[A[i,j,k,l]]*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])u[k]-X[i,j]-X_reduced[i,j]}
        J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}  
        \f]  
910    
911       @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.       or
      @type u: L{escript.Data} or None  
      @return : flux  
      @rtype : L{escript.Data}  
912    
913         M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}
914    
915         @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.
916         @type u: L{Data<escript.Data>} or None
917         @return: flux
918         @rtype: L{Data<escript.Data>}
919       """       """
920       if u==None: u=self.getSolution()       if u==None: u=self.getSolution()
921       return util.tensormult(self.getCoefficientOfGeneralPDE("A"),util.grad(u))+util.matrixmult(self.getCoefficientOfGeneralPDE("B"),u)-util.self.getCoefficientOfGeneralPDE("X")       return util.tensormult(self.getCoefficientOfGeneralPDE("A"),util.grad(u,Funtion(self.getDomain))) \
922               +util.matrixmult(self.getCoefficientOfGeneralPDE("B"),u) \
923               -util.self.getCoefficientOfGeneralPDE("X") \
924               +util.tensormult(self.getCoefficientOfGeneralPDE("A_reduced"),util.grad(u,ReducedFuntion(self.getDomain))) \
925               +util.matrixmult(self.getCoefficientOfGeneralPDE("B_reduced"),u) \
926               -util.self.getCoefficientOfGeneralPDE("X_reduced")
927     # =============================================================================     # =============================================================================
928     #   solver settings:     #   solver settings:
929     # =============================================================================     # =============================================================================
930     def setSolverMethod(self,solver=None):     def setSolverMethod(self,solver=None,preconditioner=None):
931         """         """
932         sets a new solver         sets a new solver
933    
934         @param solver: sets a new solver method.         @param solver: sets a new solver method.
935         @type solver: C{int}         @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}
936           @param preconditioner: sets a new solver method.
937         """         @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}
938         if solver==None: solve=self.DEFAULT_METHOD         """
939         if not solver==self.getSolverMethod():         if solver==None: solver=self.__solver_method
940           if preconditioner==None: preconditioner=self.__preconditioner
941           if solver==None: solver=self.DEFAULT
942           if preconditioner==None: preconditioner=self.DEFAULT
943           if not (solver,preconditioner)==self.getSolverMethod():
944             self.__solver_method=solver             self.__solver_method=solver
945               self.__preconditioner=preconditioner
946             self.__checkMatrixType()             self.__checkMatrixType()
947             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
948    
# Line 562  class LinearPDE: Line 950  class LinearPDE:
950         """         """
951         returns the name of the solver currently used         returns the name of the solver currently used
952    
953         @return : the name of the solver currently used.         @return: the name of the solver currently used.
954         @rtype: C{string}         @rtype: C{string}
955         """         """
956    
957         m=self.getSolverMethod()         m=self.getSolverMethod()
958         if m==self.DEFAULT_METHOD: return "DEFAULT_METHOD"         p=self.getSolverPackage()
959         elif m==self.DIRECT: return "DIRECT"         method=""
960         elif m==self.CHOLEVSKY: return "CHOLEVSKY"         if m[0]==self.DEFAULT: method="DEFAULT"
961         elif m==self.PCG: return "PCG"         elif m[0]==self.DIRECT: method= "DIRECT"
962         elif m==self.CR: return "CR"         elif m[0]==self.ITERATIVE: method= "ITERATIVE"
963         elif m==self.CGS: return "CGS"         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
964         elif m==self.BICGSTAB: return "BICGSTAB"         elif m[0]==self.PCG: method= "PCG"
965         elif m==self.SSOR: return "SSOR"         elif m[0]==self.CR: method= "CR"
966         elif m==self.GMRES: return "GMRES"         elif m[0]==self.CGS: method= "CGS"
967         elif m==self.PRES20: return "PRES20"         elif m[0]==self.BICGSTAB: method= "BICGSTAB"
968         elif m==self.LUMPING: return "LUMPING"         elif m[0]==self.SSOR: method= "SSOR"
969         return None         elif m[0]==self.GMRES: method= "GMRES"
970                 elif m[0]==self.PRES20: method= "PRES20"
971           elif m[0]==self.LUMPING: method= "LUMPING"
972           elif m[0]==self.AMG: method= "AMG"
973           if m[1]==self.DEFAULT: method+="+DEFAULT"
974           elif m[1]==self.JACOBI: method+= "+JACOBI"
975           elif m[1]==self.ILU0: method+= "+ILU0"
976           elif m[1]==self.ILUT: method+= "+ILUT"
977           elif m[1]==self.SSOR: method+= "+SSOR"
978           elif m[1]==self.AMG: method+= "+AMG"
979           elif m[1]==self.RILU: method+= "+RILU"
980           if p==self.DEFAULT: package="DEFAULT"
981           elif p==self.PASO: package= "PASO"
982           elif p==self.MKL: package= "MKL"
983           elif p==self.SCSL: package= "SCSL"
984           elif p==self.UMFPACK: package= "UMFPACK"
985           elif p==self.TRILINOS: package= "TRILINOS"
986           else : method="unknown"
987           return "%s solver of %s package"%(method,package)
988    
989    
990     def getSolverMethod(self):     def getSolverMethod(self):
991         """         """
992         returns the solver method         returns the solver method
993      
994         @return : the solver method currently be used.         @return: the solver method currently be used.
995         @rtype : C{int}         @rtype: C{int}
996           """
997           return self.__solver_method,self.__preconditioner
998    
999       def setSolverPackage(self,package=None):
1000           """
1001           sets a new solver package
1002    
1003           @param package: sets a new solver method.
1004           @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}, L{TRILINOS}
1005           """
1006           if package==None: package=self.DEFAULT
1007           if not package==self.getSolverPackage():
1008               self.__solver_package=package
1009               self.__checkMatrixType()
1010               self.trace("New solver is %s"%self.getSolverMethodName())
1011    
1012       def getSolverPackage(self):
1013           """
1014           returns the package of the solver
1015    
1016           @return: the solver package currently being used.
1017           @rtype: C{int}
1018         """         """
1019         return self.__solver_method         return self.__solver_package
1020    
1021     def isUsingLumping(self):     def isUsingLumping(self):
1022        """        """
1023        checks if matrix lumping is used a solver method        checks if matrix lumping is used a solver method
1024    
1025        @return : True is lumping is currently used a solver method.        @return: True is lumping is currently used a solver method.
1026        @rtype: C{bool}        @rtype: C{bool}
1027        """        """
1028        return self.getSolverMethod()==self.LUMPING        return self.getSolverMethod()[0]==self.LUMPING
1029    
1030     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
1031         """         """
1032         resets the tolerance for the solver method to tol where for an appropriate norm |.|         resets the tolerance for the solver method to tol where for an appropriate norm M{|.|}
1033    
1034                 |self.getResidual()|<tol*|self.getRightHandSide()|         M{|L{getResidual}()|<tol*|L{getRightHandSide}()|}
1035    
1036         defines the stopping criterion.         defines the stopping criterion.
1037    
1038         @param tol: new tolerance for the solver. If the tol is lower then the current tolerence         @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
1039                     the system will be resolved.                     the system will be resolved.
1040         @type solver: C{float}         @type tol: positive C{float}
1041         @raise ValueException: if tolerance is not positive.         @raise ValueError: if tolerance is not positive.
1042         """         """
1043         if not tol>0:         if not tol>0:
1044             raise ValueException,"Tolerance as to be positive"             raise ValueError,"Tolerance as to be positive"
1045         if tol<self.getTolerance(): self.__invalidateSolution()         if tol<self.getTolerance(): self.__invalidateSolution()
1046         self.trace("New tolerance %e"%tol)         self.trace("New tolerance %e"%tol)
1047         self.__tolerance=tol         self.__tolerance=tol
# Line 634  class LinearPDE: Line 1062  class LinearPDE:
1062     def isSymmetric(self):     def isSymmetric(self):
1063        """        """
1064        checks if symmetry is indicated.        checks if symmetry is indicated.
1065        
1066        @return : True is a symmetric PDE is indicated, otherwise False is returned        @return: True is a symmetric PDE is indicated, otherwise False is returned
1067        @rtype : C{bool}        @rtype: C{bool}
1068        """        """
1069        return self.__sym        return self.__sym
1070    
# Line 661  class LinearPDE: Line 1089  class LinearPDE:
1089     def setSymmetryTo(self,flag=False):     def setSymmetryTo(self,flag=False):
1090        """        """
1091        sets the symmetry flag to flag        sets the symmetry flag to flag
1092    
1093        @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.        @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.
1094        @type flag: C{bool}        @type flag: C{bool}
1095        """        """
# Line 670  class LinearPDE: Line 1098  class LinearPDE:
1098        else:        else:
1099           self.setSymmetryOff()           self.setSymmetryOff()
1100    
     
1101     # =============================================================================     # =============================================================================
1102     # function space handling for the equation as well as the solution     # function space handling for the equation as well as the solution
1103     # =============================================================================     # =============================================================================
1104     def setReducedOrderOn(self):     def setReducedOrderOn(self):
1105       """       """
1106       switches on reduced order for solution and equation representation       switches on reduced order for solution and equation representation
1107    
1108         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1109       """       """
1110       self.setReducedOrderForSolutionOn()       self.setReducedOrderForSolutionOn()
1111       self.setReducedOrderForEquationOn()       self.setReducedOrderForEquationOn()
# Line 684  class LinearPDE: Line 1113  class LinearPDE:
1113     def setReducedOrderOff(self):     def setReducedOrderOff(self):
1114       """       """
1115       switches off reduced order for solution and equation representation       switches off reduced order for solution and equation representation
1116    
1117         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1118       """       """
1119       self.setReducedOrderForSolutionOff()       self.setReducedOrderForSolutionOff()
1120       self.setReducedOrderForEquationOff()       self.setReducedOrderForEquationOff()
1121    
1122     def setReducedOrderTo(self,flag=False):     def setReducedOrderTo(self,flag=False):
1123       """       """
1124       sets order reduction for both solution and equation representation according to flag.       sets order reduction for both solution and equation representation according to flag.
1125         @param flag: if flag is True, the order reduction is switched on for both  solution and equation representation, otherwise or
      @param flag: if flag is True, the order reduction is switched on for both  solution and equation representation, otherwise or  
1126                    if flag is not present order reduction is switched off                    if flag is not present order reduction is switched off
1127       @type flag: C{bool}       @type flag: C{bool}
1128         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1129       """       """
1130       self.setReducedOrderForSolutionTo(flag)       self.setReducedOrderForSolutionTo(flag)
1131       self.setReducedOrderForEquationTo(flag)       self.setReducedOrderForEquationTo(flag)
# Line 703  class LinearPDE: Line 1134  class LinearPDE:
1134     def setReducedOrderForSolutionOn(self):     def setReducedOrderForSolutionOn(self):
1135       """       """
1136       switches on reduced order for solution representation       switches on reduced order for solution representation
1137    
1138         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1139       """       """
1140       new_fs=escript.ReducedSolution(self.getDomain())       if not self.__reduce_solution_order:
1141       if self.getFunctionSpaceForSolution()!=new_fs:           if self.__altered_coefficients:
1142                  raise RuntimeError,"order cannot be altered after coefficients have been defined."
1143           self.trace("Reduced order is used to solution representation.")           self.trace("Reduced order is used to solution representation.")
1144           self.__column_function_space=new_fs           self.__reduce_solution_order=True
1145           self.__resetSystem()           self.__resetSystem()
1146    
1147     def setReducedOrderForSolutionOff(self):     def setReducedOrderForSolutionOff(self):
1148       """       """
1149       switches off reduced order for solution representation       switches off reduced order for solution representation
1150    
1151         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1152       """       """
1153       new_fs=escript.Solution(self.getDomain())       if self.__reduce_solution_order:
1154       if self.getFunctionSpaceForSolution()!=new_fs:           if self.__altered_coefficients:
1155                  raise RuntimeError,"order cannot be altered after coefficients have been defined."
1156           self.trace("Full order is used to interpolate solution.")           self.trace("Full order is used to interpolate solution.")
1157           self.__column_function_space=new_fs           self.__reduce_solution_order=False
1158           self.__resetSystem()           self.__resetSystem()
1159    
1160     def setReducedOrderForSolutionTo(self,flag=False):     def setReducedOrderForSolutionTo(self,flag=False):
1161       """       """
1162       sets order for test functions according to flag       sets order for test functions according to flag
1163    
1164       @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or       @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or
1165                    if flag is not present order reduction is switched off                    if flag is not present order reduction is switched off
1166       @type flag: C{bool}       @type flag: C{bool}
1167         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1168       """       """
1169       if flag:       if flag:
1170          self.setReducedOrderForSolutionOn()          self.setReducedOrderForSolutionOn()
# Line 736  class LinearPDE: Line 1174  class LinearPDE:
1174     def setReducedOrderForEquationOn(self):     def setReducedOrderForEquationOn(self):
1175       """       """
1176       switches on reduced order for equation representation       switches on reduced order for equation representation
1177    
1178         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1179       """       """
1180       new_fs=escript.ReducedSolution(self.getDomain())       if not self.__reduce_equation_order:
1181       if self.getFunctionSpaceForEquation()!=new_fs:           if self.__altered_coefficients:
1182                  raise RuntimeError,"order cannot be altered after coefficients have been defined."
1183           self.trace("Reduced order is used for test functions.")           self.trace("Reduced order is used for test functions.")
1184           self.__row_function_space=new_fs           self.__reduce_equation_order=True
1185           self.__resetSystem()           self.__resetSystem()
1186    
1187     def setReducedOrderForEquationOff(self):     def setReducedOrderForEquationOff(self):
1188       """       """
1189       switches off reduced order for equation representation       switches off reduced order for equation representation
1190    
1191         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1192       """       """
1193       new_fs=escript.Solution(self.getDomain())       if self.__reduce_equation_order:
1194       if self.getFunctionSpaceForEquation()!=new_fs:           if self.__altered_coefficients:
1195                  raise RuntimeError,"order cannot be altered after coefficients have been defined."
1196           self.trace("Full order is used for test functions.")           self.trace("Full order is used for test functions.")
1197           self.__row_function_space=new_fs           self.__reduce_equation_order=False
1198           self.__resetSystem()           self.__resetSystem()
1199    
1200     def setReducedOrderForEquationTo(self,flag=False):     def setReducedOrderForEquationTo(self,flag=False):
1201       """       """
1202       sets order for test functions according to flag       sets order for test functions according to flag
1203    
1204       @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or       @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or
1205                    if flag is not present order reduction is switched off                    if flag is not present order reduction is switched off
1206       @type flag: C{bool}       @type flag: C{bool}
1207         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1208       """       """
1209       if flag:       if flag:
1210          self.setReducedOrderForEquationOn()          self.setReducedOrderForEquationOn()
# Line 773  class LinearPDE: Line 1218  class LinearPDE:
1218       """       """
1219       reassess the matrix type and, if a new matrix is needed, resets the system.       reassess the matrix type and, if a new matrix is needed, resets the system.
1220       """       """
1221       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverPackage(),self.isSymmetric())
1222       if not new_matrix_type==self.__matrix_type:       if not new_matrix_type==self.__matrix_type:
1223           self.trace("Matrix type is now %d."%new_matrix_type)           self.trace("Matrix type is now %d."%new_matrix_type)
1224           self.__matrix_type=new_matrix_type           self.__matrix_type=new_matrix_type
1225           self.__resetSystem()           self.__resetSystem()
1226     #     #
1227     #   rebuild switches :     #   rebuild switches :
1228     #     #
1229     def __invalidateSolution(self):     def __invalidateSolution(self):
1230         """         """
1231         indicates the PDE has to be resolved if the solution is requested         indicates the PDE has to be resolved if the solution is requested
# Line 792  class LinearPDE: Line 1237  class LinearPDE:
1237         """         """
1238         indicates the operator has to be rebuilt next time it is used         indicates the operator has to be rebuilt next time it is used
1239         """         """
1240         if self.__operator_isValid: self.trace("Operator has to be rebuilt.")         if self.__operator_is_Valid: self.trace("Operator has to be rebuilt.")
1241         self.__invalidateSolution()         self.__invalidateSolution()
1242         self.__operator_isValid=False         self.__operator_is_Valid=False
1243    
1244     def __invalidateRightHandSide(self):     def __invalidateRightHandSide(self):
1245         """         """
# Line 819  class LinearPDE: Line 1264  class LinearPDE:
1264         """         """
1265         self.trace("New System is built from scratch.")         self.trace("New System is built from scratch.")
1266         self.__operator=escript.Operator()         self.__operator=escript.Operator()
1267         self.__operator_isValid=False         self.__operator_is_Valid=False
1268         self.__righthandside=escript.Data()         self.__righthandside=escript.Data()
1269         self.__righthandside_isValid=False         self.__righthandside_isValid=False
1270         self.__solution=escript.Data()         self.__solution=escript.Data()
1271         self.__solution_isValid=False         self.__solution_isValid=False
1272     #     #
1273     #    system initialization:     #    system initialization:
1274     #     #
1275     def __getNewOperator(self):     def __getNewOperator(self):
1276         """         """
1277         returns an instance of a new operator         returns an instance of a new operator
# Line 877  class LinearPDE: Line 1322  class LinearPDE:
1322         if self.__righthandside.isEmpty():         if self.__righthandside.isEmpty():
1323             self.__righthandside=self.__getNewRightHandSide()             self.__righthandside=self.__getNewRightHandSide()
1324         else:         else:
1325             self.__righthandside*=0             self.__righthandside.setToZero()
1326             self.trace("Right hand side is reset to zero.")             self.trace("Right hand side is reset to zero.")
1327         return self.__righthandside         return self.__righthandside
1328    
# Line 888  class LinearPDE: Line 1333  class LinearPDE:
1333         if self.__operator.isEmpty():         if self.__operator.isEmpty():
1334             self.__operator=self.__getNewOperator()             self.__operator=self.__getNewOperator()
1335         else:         else:
1336             self.__operator.setValue(0.)             self.__operator.resetValues()
1337             self.trace("Operator reset to zero")             self.trace("Operator reset to zero")
1338         return self.__operator         return self.__operator
1339    
# Line 909  class LinearPDE: Line 1354  class LinearPDE:
1354               else:               else:
1355                  r_s=escript.Data(r,self.getFunctionSpaceForSolution())                  r_s=escript.Data(r,self.getFunctionSpaceForSolution())
1356               u.copyWithMask(r_s,col_q)               u.copyWithMask(r_s,col_q)
1357               if not self.__righthandside.isEmpty():               if not self.__righthandside.isEmpty():
1358                  self.__righthandside-=self.__operator*u                  self.__righthandside-=self.__operator*u
1359                  self.__righthandside=self.copyConstraint(self.__righthandside)                  self.__righthandside=self.copyConstraint(self.__righthandside)
1360               self.__operator.nullifyRowsAndCols(row_q,col_q,1.)               self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
# Line 920  class LinearPDE: Line 1365  class LinearPDE:
1365       """       """
1366       return the value of the coefficient name of the general PDE.       return the value of the coefficient name of the general PDE.
1367    
1368       @note This method is called by the assembling routine it can be overwritten       @note: This method is called by the assembling routine it can be overwritten
1369             to map coefficients of a particular PDE to the general PDE.             to map coefficients of a particular PDE to the general PDE.
1370         @param name: name of the coefficient requested.
      @param name: name of the coefficient requested.  
1371       @type name: C{string}       @type name: C{string}
1372       @return : the value of the coefficient  name       @return: the value of the coefficient  name
1373       @rtype : L{escript.Data}       @rtype: L{Data<escript.Data>}
1374       @raise IllegalCoefficient: if name is not one of coefficients       @raise IllegalCoefficient: if name is not one of coefficients
1375                    "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q".                    M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1376                      M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced}, M{Y_reduced}, M{d_reduced}, M{y_reduced}, M{d_contact_reduced}, M{y_contact_reduced}, M{r} or M{q}.
1377       """       """
1378       if self.hasCoefficientOfGeneralPDE(name):       if self.hasCoefficientOfGeneralPDE(name):
1379          return self.getCoefficient(name)          return self.getCoefficient(name)
# Line 938  class LinearPDE: Line 1383  class LinearPDE:
1383     def hasCoefficientOfGeneralPDE(self,name):     def hasCoefficientOfGeneralPDE(self,name):
1384       """       """
1385       checks if name is a the name of a coefficient of the general PDE.       checks if name is a the name of a coefficient of the general PDE.
1386        
1387       @param name: name of the coefficient enquired.       @param name: name of the coefficient enquired.
1388       @type name: C{string}       @type name: C{string}
1389       @return : True if name is the name of a coefficient of the general PDE. Otherwise False.       @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1390       @rtype : C{bool}       @rtype: C{bool}
1391        
1392       """       """
1393       return self.__COEFFICIENTS_OF_GENEARL_PDE.has_key(name)       return self.__COEFFICIENTS_OF_GENEARL_PDE.has_key(name)
1394    
# Line 953  class LinearPDE: Line 1398  class LinearPDE:
1398    
1399       @param name: name of the coefficient requested.       @param name: name of the coefficient requested.
1400       @type name: C{string}       @type name: C{string}
1401       @return : a coefficient name initialized to 0.       @return: a coefficient name initialized to 0.
1402       @rtype : L{escript.Data}       @rtype: L{Data<escript.Data>}
1403       @raise IllegalCoefficient: if name is not one of coefficients       @raise IllegalCoefficient: if name is not one of coefficients
1404                    "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q".                    M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1405                      M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced}, M{Y_reduced}, M{d_reduced}, M{y_reduced}, M{d_contact_reduced}, M{y_contact_reduced}, M{r} or M{q}.
1406       """       """
1407       if self.hasCoefficientOfGeneralPDE(name):       if self.hasCoefficientOfGeneralPDE(name):
1408          return escript.Data(0,self.getShapeOfCoefficientOfGeneralPDE(name),self.getFunctionSpaceForCoefficientOfGeneralPDE(name))          return escript.Data(0,self.getShapeOfCoefficientOfGeneralPDE(name),self.getFunctionSpaceForCoefficientOfGeneralPDE(name))
# Line 965  class LinearPDE: Line 1411  class LinearPDE:
1411    
1412     def getFunctionSpaceForCoefficientOfGeneralPDE(self,name):     def getFunctionSpaceForCoefficientOfGeneralPDE(self,name):
1413       """       """
1414       return the L{escript.FunctionSpace} to be used for coefficient name of the general PDE       return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name of the general PDE
1415    
1416       @param name: name of the coefficient enquired.       @param name: name of the coefficient enquired.
1417       @type name: C{string}       @type name: C{string}
1418       @return : the function space to be used for coefficient name       @return: the function space to be used for coefficient name
1419       @rtype : L{escript.FunctionSpace}       @rtype: L{FunctionSpace<escript.FunctionSpace>}
1420       @raise IllegalCoefficient: if name is not one of coefficients       @raise IllegalCoefficient: if name is not one of coefficients
1421                    "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q".                    M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1422                      M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced}, M{Y_reduced}, M{d_reduced}, M{y_reduced}, M{d_contact_reduced}, M{y_contact_reduced}, M{r} or M{q}.
1423       """       """
1424       if self.hasCoefficientOfGeneralPDE(name):       if self.hasCoefficientOfGeneralPDE(name):
1425          return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getFunctionSpace(self.getDomain())          return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getFunctionSpace(self.getDomain())
# Line 985  class LinearPDE: Line 1432  class LinearPDE:
1432    
1433       @param name: name of the coefficient enquired.       @param name: name of the coefficient enquired.
1434       @type name: C{string}       @type name: C{string}
1435       @return : the shape of the coefficient name       @return: the shape of the coefficient name
1436       @rtype : C{tuple} of C{int}       @rtype: C{tuple} of C{int}
1437       @raise IllegalCoefficient: if name is not one of coefficients       @raise IllegalCoefficient: if name is not one of coefficients
1438                    "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q".                    M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1439                      M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced}, M{Y_reduced}, M{d_reduced}, M{y_reduced}, M{d_contact_reduced}, M{y_contact_reduced}, M{r} or M{q}.
1440       """       """
1441       if self.hasCoefficientOfGeneralPDE(name):       if self.hasCoefficientOfGeneralPDE(name):
1442          return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())          return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
# Line 1003  class LinearPDE: Line 1450  class LinearPDE:
1450       """       """
1451       returns the value of the coefficient name       returns the value of the coefficient name
1452    
1453       @param name: name of the coefficient requested.       @param name: name of the coefficient requested.
1454       @type name: C{string}       @type name: C{string}
1455       @return : the value of the coefficient name       @return: the value of the coefficient name
1456       @rtype : L{escript.Data}       @rtype: L{Data<escript.Data>}
1457       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1458       """       """
1459       if self.hasCoefficient(name):       if self.hasCoefficient(name):
# Line 1020  class LinearPDE: Line 1467  class LinearPDE:
1467    
1468       @param name: name of the coefficient enquired.       @param name: name of the coefficient enquired.
1469       @type name: C{string}       @type name: C{string}
1470       @return : True if name is the name of a coefficient of the general PDE. Otherwise False.       @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1471       @rtype : C{bool}       @rtype: C{bool}
1472       """       """
1473       return self.COEFFICIENTS.has_key(name)       return self.COEFFICIENTS.has_key(name)
1474    
1475     def createCoefficient(self, name):     def createCoefficient(self, name):
1476       """       """
1477       create a L{escript.Data} object corresponding to coefficient name       create a L{Data<escript.Data>} object corresponding to coefficient name
1478    
1479       @return : a coefficient name initialized to 0.       @return: a coefficient name initialized to 0.
1480       @rtype : L{escript.Data}       @rtype: L{Data<escript.Data>}
1481       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1482       """       """
1483       if self.hasCoefficient(name):       if self.hasCoefficient(name):
# Line 1040  class LinearPDE: Line 1487  class LinearPDE:
1487    
1488     def getFunctionSpaceForCoefficient(self,name):     def getFunctionSpaceForCoefficient(self,name):
1489       """       """
1490       return the L{escript.FunctionSpace} to be used for coefficient name       return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name
1491    
1492       @param name: name of the coefficient enquired.       @param name: name of the coefficient enquired.
1493       @type name: C{string}       @type name: C{string}
1494       @return : the function space to be used for coefficient name       @return: the function space to be used for coefficient name
1495       @rtype : L{escript.FunctionSpace}       @rtype: L{FunctionSpace<escript.FunctionSpace>}
1496       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1497       """       """
1498       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1499          return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())          return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1500       else:       else:
1501          raise ValueError,"unknown coefficient %s requested"%name          raise ValueError,"unknown coefficient %s requested"%name
   
1502     def getShapeOfCoefficient(self,name):     def getShapeOfCoefficient(self,name):
1503       """       """
1504       return the shape of the coefficient name       return the shape of the coefficient name
1505    
1506       @param name: name of the coefficient enquired.       @param name: name of the coefficient enquired.
1507       @type name: C{string}       @type name: C{string}
1508       @return : the shape of the coefficient name       @return: the shape of the coefficient name
1509       @rtype : C{tuple} of C{int}       @rtype: C{tuple} of C{int}
1510       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1511       """       """
1512       if self.hasCoefficient(name):       if self.hasCoefficient(name):
# Line 1082  class LinearPDE: Line 1528  class LinearPDE:
1528       @param name: name of the coefficient enquired.       @param name: name of the coefficient enquired.
1529       @type name: C{string}       @type name: C{string}
1530       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1531         @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.
1532       """       """
1533       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1534          self.trace("Coefficient %s has been altered."%name)          self.trace("Coefficient %s has been altered."%name)
1535          if self.COEFFICIENTS[name].isAlteringOperator(): self.__invalidateOperator()          if not ((name=="q" or name=="r") and self.isUsingLumping()):
1536          if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__invalidateRightHandSide()             if self.COEFFICIENTS[name].isAlteringOperator(): self.__invalidateOperator()
1537               if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__invalidateRightHandSide()
1538       else:       else:
1539          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1540    
1541     def copyConstraint(self,u):     def copyConstraint(self,u):
1542        """        """
1543        copies the constraint into u and returns u.        copies the constraint into u and returns u.
   
       @param u: a function of rank 0 is a single PDE is solved and of shape (numSolution,) for a system of PDEs  
       @type u: L{escript.Data}  
       @return : the input u modified by the constraints.  
       @rtype : L{escript.Data}  
       @warning: u is altered if it has the appropriate L{escript.FunctionSpace}  
1544    
1545          @param u: a function of rank 0 is a single PDE is solved and of shape (numSolution,) for a system of PDEs
1546          @type u: L{Data<escript.Data>}
1547          @return: the input u modified by the constraints.
1548          @rtype: L{Data<escript.Data>}
1549          @warning: u is altered if it has the appropriate L{FunctionSpace<escript.FunctionSpace>}
1550        """        """
1551        q=self.getCoefficientOfGeneralPDE("q")        q=self.getCoefficientOfGeneralPDE("q")
1552        r=self.getCoefficientOfGeneralPDE("r")        r=self.getCoefficientOfGeneralPDE("r")
# Line 1116  class LinearPDE: Line 1563  class LinearPDE:
1563        """        """
1564        sets new values to coefficients        sets new values to coefficients
1565    
1566        @note This method is called by the assembling routine it can be overwritten        @param coefficients: new values assigned to coefficients
1567             to map coefficients of a particular PDE to the general PDE.        @keyword A: value for coefficient A.
1568          @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1569        @param name: name of the coefficient requested.        @keyword A_reduced: value for coefficient A_reduced.
1570        @type name: C{string}        @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
       @keyword A: value for coefficient A.  
       @type A: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.  
1571        @keyword B: value for coefficient B        @keyword B: value for coefficient B
1572        @type B: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.        @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1573          @keyword B_reduced: value for coefficient B_reduced
1574          @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1575        @keyword C: value for coefficient C        @keyword C: value for coefficient C
1576        @type C: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.        @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1577          @keyword C_reduced: value for coefficient C_reduced
1578          @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1579        @keyword D: value for coefficient D        @keyword D: value for coefficient D
1580        @type D: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.        @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1581          @keyword D_reduced: value for coefficient D_reduced
1582          @type D_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1583        @keyword X: value for coefficient X        @keyword X: value for coefficient X
1584        @type X: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.        @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1585          @keyword X_reduced: value for coefficient X_reduced
1586          @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1587        @keyword Y: value for coefficient Y        @keyword Y: value for coefficient Y
1588        @type Y: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.        @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1589          @keyword Y_reduced: value for coefficient Y_reduced
1590          @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.
1591        @keyword d: value for coefficient d        @keyword d: value for coefficient d
1592        @type d: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnBoundary}.        @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1593          @keyword d_reduced: value for coefficient d_reduced
1594          @type d_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1595        @keyword y: value for coefficient y        @keyword y: value for coefficient y
1596        @type y: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnBoundary}.        @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1597        @keyword d_contact: value for coefficient d_contact        @keyword d_contact: value for coefficient d_contact
1598        @type d_contact: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnContactOne}.        @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1599                         or  L{escript.FunctionOnContactZero}.        @keyword d_contact_reduced: value for coefficient d_contact_reduced
1600          @type d_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>} or  L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}.
1601        @keyword y_contact: value for coefficient y_contact        @keyword y_contact: value for coefficient y_contact
1602        @type y_contact: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnContactOne}.        @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1603                         or  L{escript.FunctionOnContactZero}.        @keyword y_contact_reduced: value for coefficient y_contact_reduced
1604          @type y_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}.
1605        @keyword r: values prescribed to the solution at the locations of constraints        @keyword r: values prescribed to the solution at the locations of constraints
1606        @type r: any type that can be interpreted as L{escript.Data} object on L{escript.Solution} or L{escript.ReducedSolution}        @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1607                 depending of reduced order is used for the solution.                 depending of reduced order is used for the solution.
1608        @keyword q: mask for location of constraints        @keyword q: mask for location of constraints
1609        @type q: any type that can be interpreted as L{escript.Data} object on L{escript.Solution} or L{escript.ReducedSolution}        @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1610                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
1611        @raise IllegalCoefficient: if an unknown coefficient keyword is used.        @raise IllegalCoefficient: if an unknown coefficient keyword is used.
   
1612        """        """
1613        # check if the coefficients are  legal:        # check if the coefficients are  legal:
1614        for i in coefficients.iterkeys():        for i in coefficients.iterkeys():
# Line 1178  class LinearPDE: Line 1636  class LinearPDE:
1636        # 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:
1637        for i,d in coefficients.iteritems():        for i,d in coefficients.iteritems():
1638          try:          try:
1639             self.COEFFICIENTS[i].setValue(self.getDomain(),self.getNumEquations(),self.getNumSolutions(),d)             self.COEFFICIENTS[i].setValue(self.getDomain(),
1640                                             self.getNumEquations(),self.getNumSolutions(),
1641                                             self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1642               self.alteredCoefficient(i)
1643            except IllegalCoefficientFunctionSpace,m:
1644                # if the function space is wrong then we try the reduced version:
1645                i_red=i+"_reduced"
1646                if (not i_red in coefficients.keys()) and i_red in self.COEFFICIENTS.keys():
1647                    try:
1648                        self.COEFFICIENTS[i_red].setValue(self.getDomain(),
1649                                                          self.getNumEquations(),self.getNumSolutions(),
1650                                                          self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1651                        self.alteredCoefficient(i_red)
1652                    except IllegalCoefficientValue,m:
1653                        raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1654                    except IllegalCoefficientFunctionSpace,m:
1655                        raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1656                else:
1657                    raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1658          except IllegalCoefficientValue,m:          except IllegalCoefficientValue,m:
1659             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1660          self.alteredCoefficient(i)        self.__altered_coefficients=True
   
1661        # check if the systrem is inhomogeneous:        # check if the systrem is inhomogeneous:
1662        if len(coefficients)>0 and not self.isUsingLumping():        if len(coefficients)>0 and not self.isUsingLumping():
1663           q=self.getCoefficientOfGeneralPDE("q")           q=self.getCoefficientOfGeneralPDE("q")
1664           r=self.getCoefficientOfGeneralPDE("r")           r=self.getCoefficientOfGeneralPDE("r")
1665           homogeneous_constraint=True           homogeneous_constraint=True
1666           if not q.isEmpty() and not r.isEmpty():           if not q.isEmpty() and not r.isEmpty():
1667               if util.Lsup(q*r)>=1.e-13*util.Lsup(r):               if util.Lsup(q*r)>0.:
1668                 self.trace("Inhomogeneous constraint detected.")                 self.trace("Inhomogeneous constraint detected.")
1669                 self.__invalidateSystem()                 self.__invalidateSystem()
1670    
   
1671     def getSystem(self):     def getSystem(self):
1672         """         """
1673         return the operator and right hand side of the PDE         return the operator and right hand side of the PDE
1674    
1675           @return: the discrete version of the PDE
1676           @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1677         """         """
1678         if not self.__operator_isValid or not self.__righthandside_isValid:         if not self.__operator_is_Valid or not self.__righthandside_isValid:
1679            if self.isUsingLumping():            if self.isUsingLumping():
1680                if not self.__operator_isValid:                if not self.__operator_is_Valid:
1681                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution(): raise TypeError,"Lumped matrix requires same order for equations and unknowns"                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1682                   if not self.getCoefficientOfGeneralPDE("A").isEmpty(): raise Warning,"Lumped matrix does not allow coefficient A"                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1683                   if not self.getCoefficientOfGeneralPDE("B").isEmpty(): raise Warning,"Lumped matrix does not allow coefficient B"                   if not self.getCoefficientOfGeneralPDE("A").isEmpty():
1684                   if not self.getCoefficientOfGeneralPDE("C").isEmpty(): raise Warning,"Lumped matrix does not allow coefficient C"                        raise ValueError,"coefficient A in lumped matrix may not be present."
1685                   mat=self.__getNewOperator()                   if not self.getCoefficientOfGeneralPDE("B").isEmpty():
1686                   self.getDomain().addPDEToSystem(mat,escript.Data(), \                        raise ValueError,"coefficient B in lumped matrix may not be present."
1687                             self.getCoefficientOfGeneralPDE("A"), \                   if not self.getCoefficientOfGeneralPDE("C").isEmpty():
1688                             self.getCoefficientOfGeneralPDE("B"), \                        raise ValueError,"coefficient C in lumped matrix may not be present."
1689                             self.getCoefficientOfGeneralPDE("C"), \                   if not self.getCoefficientOfGeneralPDE("d_contact").isEmpty():
1690                             self.getCoefficientOfGeneralPDE("D"), \                        raise ValueError,"coefficient d_contact in lumped matrix may not be present."
1691                             escript.Data(), \                   if not self.getCoefficientOfGeneralPDE("A_reduced").isEmpty():
1692                             escript.Data(), \                        raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
1693                             self.getCoefficientOfGeneralPDE("d"), \                   if not self.getCoefficientOfGeneralPDE("B_reduced").isEmpty():
1694                             escript.Data(),\                        raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
1695                             self.getCoefficientOfGeneralPDE("d_contact"), \                   if not self.getCoefficientOfGeneralPDE("C_reduced").isEmpty():
1696                             escript.Data())                        raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
1697                   self.__operator=1./(mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True))                   if not self.getCoefficientOfGeneralPDE("d_contact_reduced").isEmpty():
1698                   del mat                        raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
1699                     D=self.getCoefficientOfGeneralPDE("D")
1700                     d=self.getCoefficientOfGeneralPDE("d")
1701                     D_reduced=self.getCoefficientOfGeneralPDE("D_reduced")
1702                     d_reduced=self.getCoefficientOfGeneralPDE("d_reduced")
1703                     if not D.isEmpty():
1704                         if self.getNumSolutions()>1:
1705                            D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1706                         else:
1707                            D_times_e=D
1708                     else:
1709                        D_times_e=escript.Data()
1710                     if not d.isEmpty():
1711                         if self.getNumSolutions()>1:
1712                            d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1713                         else:
1714                            d_times_e=d
1715                     else:
1716                        d_times_e=escript.Data()
1717          
1718                     if not D_reduced.isEmpty():
1719                         if self.getNumSolutions()>1:
1720                            D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))
1721                         else:
1722                            D_reduced_times_e=D_reduced
1723                     else:
1724                        D_reduced_times_e=escript.Data()
1725                     if not d_reduced.isEmpty():
1726                         if self.getNumSolutions()>1:
1727                            d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))
1728                         else:
1729                            d_reduced_times_e=d_reduced
1730                     else:
1731                        d_reduced_times_e=escript.Data()
1732    
1733                     self.__operator=self.__getNewRightHandSide()
1734                     if hasattr(self.getDomain(), "addPDEToLumpedSystem") :
1735                        self.getDomain().addPDEToLumpedSystem(self.__operator, D_times_e, d_times_e)
1736                        self.getDomain().addPDEToLumpedSystem(self.__operator, D_reduced_times_e, d_reduced_times_e)
1737                     else:
1738                        self.getDomain().addPDEToRHS(self.__operator, \
1739                                                     escript.Data(), \
1740                                                     D_times_e, \
1741                                                     d_times_e,\
1742                                                     escript.Data())
1743                        self.getDomain().addPDEToRHS(self.__operator, \
1744                                                     escript.Data(), \
1745                                                     D_reduced_times_e, \
1746                                                     d_reduced_times_e,\
1747                                                     escript.Data())
1748                     self.__operator=1./self.__operator
1749                   self.trace("New lumped operator has been built.")                   self.trace("New lumped operator has been built.")
1750                   self.__operator_isValid=True                   self.__operator_is_Valid=True
1751                if not self.__righthandside_isValid:                if not self.__righthandside_isValid:
1752                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1753                                 self.getCoefficientOfGeneralPDE("X"), \                                 self.getCoefficientOfGeneralPDE("X"), \
1754                                 self.getCoefficientOfGeneralPDE("Y"),\                                 self.getCoefficientOfGeneralPDE("Y"),\
1755                                 self.getCoefficientOfGeneralPDE("y"),\                                 self.getCoefficientOfGeneralPDE("y"),\
1756                                 self.getCoefficientOfGeneralPDE("y_contact"))                                 self.getCoefficientOfGeneralPDE("y_contact"))
1757                     self.getDomain().addPDEToRHS(self.__righthandside, \
1758                                   self.getCoefficientOfGeneralPDE("X_reduced"), \
1759                                   self.getCoefficientOfGeneralPDE("Y_reduced"),\
1760                                   self.getCoefficientOfGeneralPDE("y_reduced"),\
1761                                   self.getCoefficientOfGeneralPDE("y_contact_reduced"))
1762                   self.trace("New right hand side as been built.")                   self.trace("New right hand side as been built.")
1763                   self.__righthandside_isValid=True                   self.__righthandside_isValid=True
1764            else:            else:
1765               if not self.__operator_isValid and not self.__righthandside_isValid:               if not self.__operator_is_Valid and not self.__righthandside_isValid:
1766                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
1767                                 self.getCoefficientOfGeneralPDE("A"), \                                 self.getCoefficientOfGeneralPDE("A"), \
1768                                 self.getCoefficientOfGeneralPDE("B"), \                                 self.getCoefficientOfGeneralPDE("B"), \
# Line 1242  class LinearPDE: Line 1774  class LinearPDE:
1774                                 self.getCoefficientOfGeneralPDE("y"), \                                 self.getCoefficientOfGeneralPDE("y"), \
1775                                 self.getCoefficientOfGeneralPDE("d_contact"), \                                 self.getCoefficientOfGeneralPDE("d_contact"), \
1776                                 self.getCoefficientOfGeneralPDE("y_contact"))                                 self.getCoefficientOfGeneralPDE("y_contact"))
1777                     self.getDomain().addPDEToSystem(self.__operator,self.__righthandside, \
1778                                   self.getCoefficientOfGeneralPDE("A_reduced"), \
1779                                   self.getCoefficientOfGeneralPDE("B_reduced"), \
1780                                   self.getCoefficientOfGeneralPDE("C_reduced"), \
1781                                   self.getCoefficientOfGeneralPDE("D_reduced"), \
1782                                   self.getCoefficientOfGeneralPDE("X_reduced"), \
1783                                   self.getCoefficientOfGeneralPDE("Y_reduced"), \
1784                                   self.getCoefficientOfGeneralPDE("d_reduced"), \
1785                                   self.getCoefficientOfGeneralPDE("y_reduced"), \
1786                                   self.getCoefficientOfGeneralPDE("d_contact_reduced"), \
1787                                   self.getCoefficientOfGeneralPDE("y_contact_reduced"))
1788                   self.__applyConstraint()                   self.__applyConstraint()
1789                   self.__righthandside=self.copyConstraint(self.__righthandside)                   self.__righthandside=self.copyConstraint(self.__righthandside)
1790                   self.trace("New system has been built.")                   self.trace("New system has been built.")
1791                   self.__operator_isValid=True                   self.__operator_is_Valid=True
1792                   self.__righthandside_isValid=True                   self.__righthandside_isValid=True
1793               elif not self.__righthandside_isValid:               elif not self.__righthandside_isValid:
1794                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
# Line 1253  class LinearPDE: Line 1796  class LinearPDE:
1796                                 self.getCoefficientOfGeneralPDE("Y"),\                                 self.getCoefficientOfGeneralPDE("Y"),\
1797                                 self.getCoefficientOfGeneralPDE("y"),\                                 self.getCoefficientOfGeneralPDE("y"),\
1798                                 self.getCoefficientOfGeneralPDE("y_contact"))                                 self.getCoefficientOfGeneralPDE("y_contact"))
1799                     self.getDomain().addPDEToRHS(self.__righthandside, \
1800                                   self.getCoefficientOfGeneralPDE("X_reduced"), \
1801                                   self.getCoefficientOfGeneralPDE("Y_reduced"),\
1802                                   self.getCoefficientOfGeneralPDE("y_reduced"),\
1803                                   self.getCoefficientOfGeneralPDE("y_contact_reduced"))
1804                   self.__righthandside=self.copyConstraint(self.__righthandside)                   self.__righthandside=self.copyConstraint(self.__righthandside)
1805                   self.trace("New right hand side has been built.")                   self.trace("New right hand side has been built.")
1806                   self.__righthandside_isValid=True                   self.__righthandside_isValid=True
1807               elif not self.__operator_isValid:               elif not self.__operator_is_Valid:
1808                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
1809                              self.getCoefficientOfGeneralPDE("A"), \                              self.getCoefficientOfGeneralPDE("A"), \
1810                              self.getCoefficientOfGeneralPDE("B"), \                              self.getCoefficientOfGeneralPDE("B"), \
# Line 1268  class LinearPDE: Line 1816  class LinearPDE:
1816                              escript.Data(),\                              escript.Data(),\
1817                              self.getCoefficientOfGeneralPDE("d_contact"), \                              self.getCoefficientOfGeneralPDE("d_contact"), \
1818                              escript.Data())                              escript.Data())
1819                     self.getDomain().addPDEToSystem(self.__operator,escript.Data(), \
1820                                self.getCoefficientOfGeneralPDE("A_reduced"), \
1821                                self.getCoefficientOfGeneralPDE("B_reduced"), \
1822                                self.getCoefficientOfGeneralPDE("C_reduced"), \
1823                                self.getCoefficientOfGeneralPDE("D_reduced"), \
1824                                escript.Data(), \
1825                                escript.Data(), \
1826                                self.getCoefficientOfGeneralPDE("d_reduced"), \
1827                                escript.Data(),\
1828                                self.getCoefficientOfGeneralPDE("d_contact_reduced"), \
1829                                escript.Data())
1830                   self.__applyConstraint()                   self.__applyConstraint()
1831                   self.trace("New operator has been built.")                   self.trace("New operator has been built.")
1832                   self.__operator_isValid=True                   self.__operator_is_Valid=True
1833         return (self.__operator,self.__righthandside)         return (self.__operator,self.__righthandside)
1834    
1835    
1836    class Poisson(LinearPDE):
   
 class AdvectivePDE(LinearPDE):  
1837     """     """
1838     Class to handle a linear PDE dominated by advective terms:     Class to define a Poisson equation problem, which is genear L{LinearPDE} of the form
   
    class to define a linear PDE of the form  
   
    \f[  
    -(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i  
    \f]  
1839    
1840     with boundary conditons:     M{-grad(grad(u)[j])[j] = f}
1841    
1842     \f[     with natural boundary conditons
    n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i  
    \f]  
1843    
1844     and contact conditions     M{n[j]*grad(u)[j] = 0 }
   
    \f[  
    n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d^{contact}_{ik}[u_k] = - n_j*X_{ij} + y^{contact}_{i}  
    \f]  
1845    
1846     and constraints:     and constraints:
1847    
1848     \f[     M{u=0} where M{q>0}
1849     u_i=r_i \quad \mathrm{where} \quad q_i>0  
    \f]  
1850     """     """
    def __init__(self,domain,numEquations=0,numSolutions=0,xi=None,debug=False):  
       LinearPDE.__init__(self,domain,numEquations,numSolutions,debug)  
       if xi==None:  
          self.__xi=AdvectivePDE.ELMAN_RAMAGE  
       else:  
          self.__xi=xi  
       self.__Xi=escript.Data()  
1851    
1852     def __calculateXi(self,peclet_factor,Z,h):     def __init__(self,domain,debug=False):
1853         Z_max=util.Lsup(Z)       """
1854         if Z_max>0.:       initializes a new Poisson equation
           return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.TOL)  
        else:  
           return 0.  
1855    
1856     def setValue(self,**args):       @param domain: domain of the PDE
1857         if "A" in args.keys()   or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data()       @type domain: L{Domain<escript.Domain>}
1858         LinearPDE.setValue(**args)       @param debug: if True debug informations are printed.
   
    def ELMAN_RAMAGE(P):  
      """   """  
      return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))  
    def SIMPLIFIED_BROOK_HUGHES(P):  
      """   """  
      c=(P-3.).whereNegative()  
      return P/6.*c+1./2.*(1.-c)  
   
    def HALF(P):  
     """ """  
     return escript.Scalar(0.5,P.getFunctionSpace())  
   
    def getXi(self):  
       if self.__Xi.isEmpty():  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          A=self.getCoefficient("A")  
          h=self.getDomain().getSize()  
          self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))  
          if not C.isEmpty() or not B.isEmpty():  
             if not C.isEmpty() and not B.isEmpty():  
                 Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))  
                 if self.getNumEquations()>1:  
                    if self.getNumSolutions()>1:  
                       for i in range(self.getNumEquations()):  
                          for k in range(self.getNumSolutions()):  
                             for l in range(self.getDim()): Z2+=(C[i,k,l]-B[i,l,k])**2  
                    else:  
                       for i in range(self.getNumEquations()):  
                          for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2  
                 else:  
                    if self.getNumSolutions()>1:  
                       for k in range(self.getNumSolutions()):  
                          for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2  
                    else:  
                       for l in range(self.getDim()): Z2+=(C[l]-B[l])**2  
                 length_of_Z=util.sqrt(Z2)  
             elif C.isEmpty():  
               length_of_Z=util.length(B)  
             else:  
               length_of_Z=util.length(C)  
1859    
1860              Z_max=util.Lsup(length_of_Z)       """
1861              if Z_max>0.:       super(Poisson, self).__init__(domain,1,1,debug)
1862                 length_of_A=util.length(A)       self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1863                 A_max=util.Lsup(length_of_A)                          "f_reduced": PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1864                 if A_max>0:                          "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1865                      inv_A=1./(length_of_A+A_max*self.TOL)       self.setSymmetryOn()
                else:  
                     inv_A=1./self.TOL  
                peclet_number=length_of_Z*h/2*inv_A  
                xi=self.__xi(peclet_number)  
                self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL)  
                print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z)  
       return self.__Xi  
1866    
1867       def setValue(self,**coefficients):
1868         """
1869         sets new values to coefficients
1870    
1871         @param coefficients: new values assigned to coefficients
1872         @keyword f: value for right hand side M{f}
1873         @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1874         @keyword q: mask for location of constraints
1875         @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>}
1876                   depending of reduced order is used for the representation of the equation.
1877         @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1878         """
1879         super(Poisson, self).setValue(**coefficients)
1880    
1881     def getCoefficientOfGeneralPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
1882       """       """
1883       return the value of the coefficient name of the general PDE       return the value of the coefficient name of the general PDE
1884         @param name: name of the coefficient requested.
1885       @param name:       @type name: C{string}
1886         @return: the value of the coefficient  name
1887         @rtype: L{Data<escript.Data>}
1888         @raise IllegalCoefficient: if name is not one of coefficients
1889                      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}.
1890         @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1891       """       """
      if not self.getNumEquations() == self.getNumSolutions():  
           raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."  
   
1892       if name == "A" :       if name == "A" :
1893           A=self.getCoefficient("A")           return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          if B.isEmpty() and C.isEmpty():  
             Aout=A  
          else:  
             if A.isEmpty():  
                Aout=self.createNewCoefficient("A")  
             else:  
                Aout=A[:]  
             Xi=self.getXi()  
             if self.getNumEquations()>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 not C.isEmpty() and not B.isEmpty():  
                                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])  
                             elif C.isEmpty():  
                                for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]  
                             else:  
                                for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]  
             else:  
                 for j in range(self.getDim()):  
                    for l in range(self.getDim()):  
                       if not C.isEmpty() and not B.isEmpty():  
                           Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])  
                       elif C.isEmpty():  
                           Aout[j,l]+=Xi*B[j]*B[l]  
                       else:  
                           Aout[j,l]+=Xi*C[j]*C[l]  
          return Aout  
1894       elif name == "B" :       elif name == "B" :
1895           B=self.getCoefficient("B")           return escript.Data()
          C=self.getCoefficient("C")  
          D=self.getCoefficient("D")  
          if C.isEmpty() or D.isEmpty():  
             Bout=B  
          else:  
             Xi=self.getXi()  
             if B.isEmpty():  
                 Bout=self.createNewCoefficient("B")  
             else:  
                 Bout=B[:]  
             if self.getNumEquations()>1:  
                for k in range(self.getNumSolutions()):  
                   for p in range(self.getNumEquations()):  
                      tmp=Xi*D[p,k]  
                      for i in range(self.getNumEquations()):  
                         for j in range(self.getDim()):  
                            Bout[i,j,k]+=tmp*C[p,i,j]  
             else:  
                tmp=Xi*D  
                for j in range(self.getDim()): Bout[j]+=tmp*C[j]  
          return Bout  
1896       elif name == "C" :       elif name == "C" :
1897           B=self.getCoefficient("B")           return escript.Data()
          C=self.getCoefficient("C")  
          D=self.getCoefficient("D")  
          if B.isEmpty() or D.isEmpty():  
             Cout=C  
          else:  
             Xi=self.getXi()  
             if C.isEmpty():  
                 Cout=self.createNewCoefficient("C")  
             else:  
                 Cout=C[:]  
             if self.getNumEquations()>1:  
                for k in range(self.getNumSolutions()):  
                    for p in range(self.getNumEquations()):  
                       tmp=Xi*D[p,k]  
                       for i in range(self.getNumEquations()):  
                         for l in range(self.getDim()):  
                                  Cout[i,k,l]+=tmp*B[p,l,i]  
             else:  
                tmp=Xi*D  
                for j in range(self.getDim()): Cout[j]+=tmp*B[j]  
          return Cout  
1898       elif name == "D" :       elif name == "D" :
1899           return self.getCoefficient("D")           return escript.Data()
1900       elif name == "X" :       elif name == "X" :
1901           X=self.getCoefficient("X")           return escript.Data()
          Y=self.getCoefficient("Y")  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):  
             Xout=X  
          else:  
             if X.isEmpty():  
                 Xout=self.createNewCoefficient("X")  
             else:  
                 Xout=X[:]  
             Xi=self.getXi()  
             if self.getNumEquations()>1:  
                  for p in range(self.getNumEquations()):  
                     tmp=Xi*Y[p]  
                     for i in range(self.getNumEquations()):  
                        for j in range(self.getDim()):  
                           if not C.isEmpty() and not B.isEmpty():  
                              Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])  
                           elif C.isEmpty():  
                              Xout[i,j]-=tmp*B[p,j,i]  
                           else:  
                              Xout[i,j]+=tmp*C[p,i,j]  
             else:  
                  tmp=Xi*Y  
                  for j in range(self.getDim()):  
                     if not C.isEmpty() and not B.isEmpty():  
                        Xout[j]+=tmp*(C[j]-B[j])  
                     elif C.isEmpty():  
                        Xout[j]-=tmp*B[j]  
                     else:  
                        Xout[j]+=tmp*C[j]  
          return Xout  
1902       elif name == "Y" :       elif name == "Y" :
1903           return self.getCoefficient("Y")           return self.getCoefficient("f")
1904       elif name == "d" :       elif name == "d" :
1905           return self.getCoefficient("d")           return escript.Data()
1906       elif name == "y" :       elif name == "y" :
1907           return self.getCoefficient("y")           return escript.Data()
1908       elif name == "d_contact" :       elif name == "d_contact" :
1909           return self.getCoefficient("d_contact")           return escript.Data()
1910       elif name == "y_contact" :       elif name == "y_contact" :
1911           return self.getCoefficient("y_contact")           return escript.Data()
1912         elif name == "A_reduced" :
1913             return escript.Data()
1914         elif name == "B_reduced" :
1915             return escript.Data()
1916         elif name == "C_reduced" :
1917             return escript.Data()
1918         elif name == "D_reduced" :
1919             return escript.Data()
1920         elif name == "X_reduced" :
1921             return escript.Data()
1922         elif name == "Y_reduced" :
1923             return self.getCoefficient("f_reduced")
1924         elif name == "d_reduced" :
1925             return escript.Data()
1926         elif name == "y_reduced" :
1927             return escript.Data()
1928         elif name == "d_contact_reduced" :
1929             return escript.Data()
1930         elif name == "y_contact_reduced" :
1931             return escript.Data()
1932       elif name == "r" :       elif name == "r" :
1933           return self.getCoefficient("r")           return escript.Data()
1934       elif name == "q" :       elif name == "q" :
1935           return self.getCoefficient("q")           return self.getCoefficient("q")
1936       else:       else:
1937           raise SystemError,"unknown PDE coefficient %s",name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
   
1938    
1939  class Poisson(LinearPDE):  class Helmholtz(LinearPDE):
1940     """     """
1941     Class to define a Poisson equation problem:     Class to define a Helmhotz equation problem, which is genear L{LinearPDE} of the form
1942    
1943       M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
1944    
1945     class to define a linear PDE of the form     with natural boundary conditons
1946     \f[  
1947     -u_{,jj} = f     M{k*n[j]*grad(u)[j] = g- S{alpha}u }
    \f]  
   
    with boundary conditons:  
   
    \f[  
    n_j*u_{,j} = 0  
    \f]  
1948    
1949     and constraints:     and constraints:
1950    
1951     \f[     M{u=r} where M{q>0}
1952     u=0 \quad \mathrm{where} \quad q>0  
    \f]  
1953     """     """
1954    
1955     def __init__(self,domain,f=escript.Data(),q=escript.Data(),debug=False):     def __init__(self,domain,debug=False):
1956         LinearPDE.__init__(self,domain,1,1,debug)       """
1957         self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),       initializes a new Poisson equation
1958                            "q": PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}  
1959         self.setSymmetryOn()       @param domain: domain of the PDE
1960         self.setValue(f,q)       @type domain: L{Domain<escript.Domain>}
1961         @param debug: if True debug informations are printed.
1962     def setValue(self,f=escript.Data(),q=escript.Data()):  
1963         """set value of PDE parameters f and q"""       """
1964         self._LinearPDE__setValue(f=f,q=q)       super(Helmholtz, self).__init__(domain,1,1,debug)
1965         self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1966                            "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1967                            "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1968                            "f_reduced": PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1969                            "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1970                            "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1971                            "g_reduced": PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1972                            "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1973                            "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1974         self.setSymmetryOn()
1975    
1976       def setValue(self,**coefficients):
1977         """
1978         sets new values to coefficients
1979    
1980         @param coefficients: new values assigned to coefficients
1981         @keyword omega: value for coefficient M{S{omega}}
1982         @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1983         @keyword k: value for coefficeint M{k}
1984         @type k: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1985         @keyword f: value for right hand side M{f}
1986         @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1987         @keyword alpha: value for right hand side M{S{alpha}}
1988         @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1989         @keyword g: value for right hand side M{g}
1990         @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1991         @keyword r: prescribed values M{r} for the solution in constraints.
1992         @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1993                   depending of reduced order is used for the representation of the equation.
1994         @keyword q: mask for location of constraints
1995         @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1996                   depending of reduced order is used for the representation of the equation.
1997         @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1998         """
1999         super(Helmholtz, self).setValue(**coefficients)
2000    
2001     def getCoefficientOfGeneralPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
2002       """       """
2003       return the value of the coefficient name of the general PDE       return the value of the coefficient name of the general PDE
2004    
2005       @param name:       @param name: name of the coefficient requested.
2006         @type name: C{string}
2007         @return: the value of the coefficient  name
2008         @rtype: L{Data<escript.Data>}
2009         @raise IllegalCoefficient: if name is not one of coefficients
2010                      "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}.
2011         @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2012       """       """
2013       if name == "A" :       if name == "A" :
2014           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
2015       elif name == "B" :       elif name == "B" :
2016           return escript.Data()           return escript.Data()
2017       elif name == "C" :       elif name == "C" :
2018           return escript.Data()           return escript.Data()
2019       elif name == "D" :       elif name == "D" :
2020           return escript.Data()           return self.getCoefficient("omega")
2021       elif name == "X" :       elif name == "X" :
2022           return escript.Data()           return escript.Data()
2023       elif name == "Y" :       elif name == "Y" :
2024           return self.getCoefficient("f")           return self.getCoefficient("f")
2025       elif name == "d" :       elif name == "d" :
2026           return escript.Data()           return self.getCoefficient("alpha")
2027       elif name == "y" :       elif name == "y" :
2028           return escript.Data()           return self.getCoefficient("g")
2029       elif name == "d_contact" :       elif name == "d_contact" :
2030           return escript.Data()           return escript.Data()
2031       elif name == "y_contact" :       elif name == "y_contact" :
2032           return escript.Data()           return escript.Data()
2033       elif name == "r" :       elif name == "A_reduced" :
2034             return escript.Data()
2035         elif name == "B_reduced" :
2036             return escript.Data()
2037         elif name == "C_reduced" :
2038             return escript.Data()
2039         elif name == "D_reduced" :
2040             return escript.Data()
2041         elif name == "X_reduced" :
2042           return escript.Data()           return escript.Data()
2043         elif name == "Y_reduced" :
2044             return self.getCoefficient("f_reduced")
2045         elif name == "d_reduced" :
2046             return escript.Data()
2047         elif name == "y_reduced" :
2048            return self.getCoefficient("g_reduced")
2049         elif name == "d_contact_reduced" :
2050             return escript.Data()
2051         elif name == "y_contact_reduced" :
2052             return escript.Data()
2053         elif name == "r" :
2054             return self.getCoefficient("r")
2055       elif name == "q" :       elif name == "q" :
2056           return self.getCoefficient("q")           return self.getCoefficient("q")
2057       else:       else:
2058           raise SystemError,"unknown PDE coefficient %s",name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2059    
2060  class LameEquation(LinearPDE):  class LameEquation(LinearPDE):
2061     """     """
2062     Class to define a Lame equation problem:     Class to define a Lame equation problem:
2063    
2064     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] }
2065     \f[  
2066     -(\mu (u_{i,j}+u_{j,i}))_{,j} - \lambda u_{j,ji}} = F_i -\sigma_{ij,j}     with natural boundary conditons:
2067     \f]  
2068       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] }
    with boundary conditons:  
   
    \f[  
    n_j(\mu (u_{i,j}+u_{j,i})-sigma_{ij}) + n_i\lambda u_{j,j} = f_i  
    \f]  
2069    
2070     and constraints:     and constraints:
2071    
2072     \f[     M{u[i]=r[i]} where M{q[i]>0}
2073     u_i=r_i \quad \mathrm{where} \quad q_i>0  
    \f]  
2074     """     """
2075    
2076     def __init__(self,domain,f=escript.Data(),q=escript.Data(),debug=False):     def __init__(self,domain,debug=False):
2077         LinearPDE.__init__(self,domain,domain.getDim(),domain.getDim(),debug)        super(LameEquation, self).__init__(domain,\
2078         self.COEFFICIENTS={ "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),                                           domain.getDim(),domain.getDim(),debug)
2079          self.COEFFICIENTS={ "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
2080                            "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),                            "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
2081                            "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),                            "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2082                            "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),                            "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
2083                            "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),                            "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2084                            "r"            : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH),                            "r"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
2085                            "q"            : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
2086         self.setSymmetryOn()        self.setSymmetryOn()
2087    
2088     def setValue(self,lame_lambda=escript.Data(),lame_mu=escript.Data(),F=escript.Data(),sigma=escript.Data(),f=escript.Data(),r=escript.Data(),q=escript.Data()):     def setValues(self,**coefficients):
2089         """set value of PDE parameters"""       """
2090         self._LinearPDE__setValue(lame_lambda=lame_lambda, \       sets new values to coefficients
2091                                   lame_mu=lame_mu, \  
2092                                   F=F, \       @param coefficients: new values assigned to coefficients
2093                                   sigma=sigma, \       @keyword lame_mu: value for coefficient M{S{mu}}
2094                                   f=f, \       @type lame_mu: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2095                                   r=r, \       @keyword lame_lambda: value for coefficient M{S{lambda}}
2096                                   q=q)       @type lame_lambda: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2097         @keyword F: value for internal force M{F}
2098         @type F: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}
2099         @keyword sigma: value for initial stress M{S{sigma}}
2100         @type sigma: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}
2101         @keyword f: value for extrenal force M{f}
2102         @type f: any type that can be casted to L{Vector<escript.Vector>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2103         @keyword r: prescribed values M{r} for the solution in constraints.
2104         @type r: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2105                   depending of reduced order is used for the representation of the equation.
2106         @keyword q: mask for location of constraints
2107         @type q: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2108                   depending of reduced order is used for the representation of the equation.
2109         @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2110         """
2111         super(LameEquation, self).setValues(**coefficients)
2112    
2113     def getCoefficientOfGeneralPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
2114       """       """
2115       return the value of the coefficient name of the general PDE       return the value of the coefficient name of the general PDE
2116    
2117       @param name:       @param name: name of the coefficient requested.
2118         @type name: C{string}
2119         @return: the value of the coefficient  name
2120         @rtype: L{Data<escript.Data>}
2121         @raise IllegalCoefficient: if name is not one of coefficients
2122                      "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}.
2123         @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2124       """       """
2125       if name == "A" :       if name == "A" :
2126           out =self.createCoefficientOfGeneralPDE("A")           out =self.createCoefficientOfGeneralPDE("A")
# Line 1662  class LameEquation(LinearPDE): Line 2148  class LameEquation(LinearPDE):
2148           return escript.Data()           return escript.Data()
2149       elif name == "y_contact" :       elif name == "y_contact" :
2150           return escript.Data()           return escript.Data()
2151         elif name == "A_reduced" :
2152             return escript.Data()
2153         elif name == "B_reduced" :
2154             return escript.Data()
2155         elif name == "C_reduced" :
2156             return escript.Data()
2157         elif name == "D_reduced" :
2158             return escript.Data()
2159         elif name == "X_reduced" :
2160             return escript.Data()
2161         elif name == "Y_reduced" :
2162             return escript.Data()
2163         elif name == "d_reduced" :
2164             return escript.Data()
2165         elif name == "y_reduced" :
2166             return escript.Data()
2167         elif name == "d_contact_reduced" :
2168             return escript.Data()
2169         elif name == "y_contact_reduced" :
2170             return escript.Data()
2171       elif name == "r" :       elif name == "r" :
2172           return self.getCoefficient("r")           return self.getCoefficient("r")
2173       elif name == "q" :       elif name == "q" :
2174           return self.getCoefficient("q")           return self.getCoefficient("q")
2175       else:       else:
2176           raise SystemError,"unknown PDE coefficient %s",name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2177    
 # $Log$  
 # Revision 1.11  2005/08/23 01:24:28  jgs  
 # Merge of development branch dev-02 back to main trunk on 2005-08-23  
 #  
 # Revision 1.10  2005/08/12 01:45:36  jgs  
 # erge of development branch dev-02 back to main trunk on 2005-08-12  
 #  
 # Revision 1.9.2.4  2005/08/22 07:11:09  gross  
 # some problems with LinearPDEs fixed.  
 #  
 # Revision 1.9.2.3  2005/08/18 04:48:48  gross  
 # the methods SetLumping*() are removed. Lumping is set trough setSolverMethod(LinearPDE.LUMPING)  
 #  
 # Revision 1.9.2.2  2005/08/18 04:39:32  gross  
 # the constants have been removed from util.py as they not needed anymore. PDE related constants are accessed through LinearPDE attributes now  
 #  
 # Revision 1.9.2.1  2005/07/29 07:10:27  gross  
 # new functions in util and a new pde type in linearPDEs  
 #  
 # Revision 1.1.2.25  2005/07/28 04:21:09  gross  
 # Lame equation: (linear elastic, isotropic) added  
 #  
 # Revision 1.1.2.24  2005/07/22 06:37:11  gross  
 # some extensions to modellib and linearPDEs  
 #  
 # Revision 1.1.2.23  2005/05/13 00:55:20  cochrane  
 # Fixed up some docstrings.  Moved module-level functions to top of file so  
 # that epydoc and doxygen can pick them up properly.  
 #  
 # Revision 1.1.2.22  2005/05/12 11:41:30  gross  
 # some basic Models have been added  
 #  
 # Revision 1.1.2.21  2005/05/12 07:16:12  cochrane  
 # Moved ELMAN_RAMAGE, SIMPLIFIED_BROOK_HUGHES, and HALF functions to bottom of  
 # file so that the AdvectivePDE class is picked up by doxygen.  Some  
 # reformatting of docstrings.  Addition of code to make equations come out  
 # as proper LaTeX.  
 #  
 # Revision 1.1.2.20  2005/04/15 07:09:08  gross  
 # some problems with functionspace and linearPDEs fixed.  
 #  
 # Revision 1.1.2.19  2005/03/04 05:27:07  gross  
 # bug in SystemPattern fixed.  
 #  
 # Revision 1.1.2.18  2005/02/08 06:16:45  gross  
 # Bugs in AdvectivePDE fixed, AdvectiveTest is stable but more testing is needed  
 #  
 # Revision 1.1.2.17  2005/02/08 05:56:19  gross  
 # Reference Number handling added  
 #  
 # Revision 1.1.2.16  2005/02/07 04:41:28  gross  
 # some function exposed to python to make mesh merging running  
 #  
 # Revision 1.1.2.15  2005/02/03 00:14:44  gross  
 # timeseries add and ESySParameter.py renames esysXML.py for consistence  
 #  
 # Revision 1.1.2.14  2005/02/01 06:44:10  gross  
 # new implementation of AdvectivePDE which now also updates right hand side. systems of PDEs are still not working  
 #  
 # Revision 1.1.2.13  2005/01/25 00:47:07  gross  
 # updates in the documentation  
 #  
 # Revision 1.1.2.12  2005/01/12 01:28:04  matt  
 # Added createCoefficient method for linearPDEs.  
 #  
 # Revision 1.1.2.11  2005/01/11 01:55:34  gross  
 # a problem in linearPDE class fixed  
 #  
 # Revision 1.1.2.10  2005/01/07 01:13:29  gross  
 # some bugs in linearPDE fixed  
 #  
 # Revision 1.1.2.9  2005/01/06 06:24:58  gross  
 # some bugs in slicing fixed  
 #  
 # Revision 1.1.2.8  2005/01/05 04:21:40  gross  
 # FunctionSpace checking/matchig in slicing added  
 #  
 # Revision 1.1.2.7  2004/12/29 10:03:41  gross  
 # bug in setValue fixed  
 #  
 # Revision 1.1.2.6  2004/12/29 05:29:59  gross  
 # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()  
 #  
 # Revision 1.1.2.5  2004/12/29 00:18:41  gross  
 # AdvectivePDE added  
 #  
 # Revision 1.1.2.4  2004/12/24 06:05:41  gross  
 # some changes in linearPDEs to add AdevectivePDE  
 #  
 # Revision 1.1.2.3  2004/12/16 00:12:34  gross  
 # __init__ of LinearPDE does not accept any coefficient anymore  
 #  
 # Revision 1.1.2.2  2004/12/14 03:55:01  jgs  
 # *** empty log message ***  
 #  
 # Revision 1.1.2.1  2004/12/12 22:53:47  gross  
 # linearPDE has been renamed LinearPDE  
 #  
 # Revision 1.1.1.1.2.7  2004/12/07 10:13:08  gross  
 # GMRES added  
 #  
 # Revision 1.1.1.1.2.6  2004/12/07 03:19:50  gross  
 # options for GMRES and PRES20 added  
 #  
 # Revision 1.1.1.1.2.5  2004/12/01 06:25:15  gross  
 # some small changes  
 #  
 # Revision 1.1.1.1.2.4  2004/11/24 01:50:21  gross  
 # Finley solves 4M unknowns now  
 #  
 # Revision 1.1.1.1.2.3  2004/11/15 06:05:26  gross  
 # poisson solver added  
 #  
 # Revision 1.1.1.1.2.2  2004/11/12 06:58:15  gross  
 # 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  
 #  
 # Revision 1.1.1.1.2.1  2004/10/28 22:59:22  gross  
 # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed  
 #  
 # Revision 1.1.1.1  2004/10/26 06:53:56  jgs  
 # initial import of project esys2  
 #  
 # Revision 1.3.2.3  2004/10/26 06:43:48  jgs  
 # committing Lutz's and Paul's changes to brach jgs  
 #  
 # Revision 1.3.4.1  2004/10/20 05:32:51  cochrane  
 # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.  
 #  
 # Revision 1.3  2004/09/23 00:53:23  jgs  
 # minor fixes  
 #  
 # Revision 1.1  2004/08/28 12:58:06  gross  
 # SimpleSolve is not running yet: problem with == of functionsspace  
 #  
 #  

Legend:
Removed from v.148  
changed lines
  Added in v.1302

  ViewVC Help
Powered by ViewVC 1.1.26