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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26