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

revision 104 by jgs, Fri Dec 17 07:43:12 2004 UTC revision 122 by jgs, Thu Jun 9 05:38:05 2005 UTC
# Line 3  Line 3
3  ## @file linearPDEs.py  ## @file linearPDEs.py
4
5  """  """
6  @brief Functions and classes for linear PDEs  Functions and classes for linear PDEs
7  """  """
8
9  import escript  import escript
10  import util  import util
11  import numarray  import numarray
12
def identifyDomain(domain=None,data={}):
"""
@brief Return the Domain which is equal to the input domain (if not None)
and is the domain of all Data objects in the dictionary data.
An exception is raised if this is not possible

@param domain
@param data
"""
# get the domain used by any Data object in the list data:
data_domain=None
for d in data.itervalues():
if isinstance(d,escript.Data):
if not d.isEmpty(): data_domain=d.getDomain()
# check if domain and data_domain are identical?
if domain == None:
if data_domain == None:
raise ValueError,"Undefined PDE domain. Specify a domain or use a Data class object as coefficient"
else:
if data_domain == None:
data_domain=domain
else:
if not data_domain == domain:
raise ValueError,"Domain of coefficients doesnot match specified domain"
# now we check if all Data class object coefficients are defined on data_domain:
for i,d in data.iteritems():
if isinstance(d,escript.Data):
if not d.isEmpty():
if not data_domain==d.getDomain():
raise ValueError,"Illegal domain for coefficient %s."%i
# done:
return data_domain

def identifyNumEquationsAndSolutions(dim,coef={}):
# get number of equations and number of unknowns:
numEquations=0
numSolutions=0
for i in coef.iterkeys():
if not coef[i].isEmpty():
res=_PDECoefficientTypes[i].estimateNumEquationsAndNumSolutions(coef[i].getShape(),dim)
if res==None:
raise ValueError,"Illegal shape %s of coefficient %s"%(coef[i].getShape().__str__(),i)
else:
numEquations=max(numEquations,res[0])
numSolutions=max(numSolutions,res[1])
return numEquations,numSolutions

13
14  def _CompTuple2(t1,t2):  def _CompTuple2(t1,t2):
15     """     """
16     @brief     Compare two tuples
17
18     @param t1     \param t1 The first tuple
19     @param t2     \param t2 The second tuple
20     """     """
21
22     dif=t1[0]+t1[1]-(t2[0]+t2[1])     dif=t1[0]+t1[1]-(t2[0]+t2[1])
23     if dif<0: return 1     if dif<0: return 1
24     elif dif>0: return -1     elif dif>0: return -1
25     else: return 0     else: return 0
26
27  class PDECoefficientType:  def ELMAN_RAMAGE(P):
28        return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
29
30    def SIMPLIFIED_BROOK_HUGHES(P):
31        c=(P-3.).whereNegative()
32        return P/6.*c+1./2.*(1.-c)
33
34    def HALF(P):
35        return escript.Scalar(0.5,P.getFunctionSpace())
36
37    class PDECoefficient:
38      """      """
39      @brief      A class for PDE coefficients
40      """      """
41      # identifier for location of Data objects defining coefficients      # identifier for location of Data objects defining COEFFICIENTS
42      INTERIOR=0      INTERIOR=0
43      BOUNDARY=1      BOUNDARY=1
44      CONTACT=2      CONTACT=2
45      CONTINUOUS=3      CONTINUOUS=3
46      # identifier in the pattern of coefficients:      # identifier in the pattern of COEFFICIENTS:
47      # the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the      # the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the
48      # number of unknowns.      # number of unknowns.
49      EQUATION=3      EQUATION=3
# Line 91  class PDECoefficientType: Line 55  class PDECoefficientType:
55      BOTH=7      BOTH=7
56      def __init__(self,where,pattern,altering):      def __init__(self,where,pattern,altering):
57         """         """
58         @brief Initialise a PDE Coefficient type         Initialise a PDE Coefficient type
59         """         """
60         self.what=where         self.what=where
61         self.pattern=pattern         self.pattern=pattern
62         self.altering=altering         self.altering=altering
63           self.resetValue()
64
65        def resetValue(self):
66           """
67           resets coefficient value to default
68           """
69           self.value=escript.Data()
70
71      def getFunctionSpace(self,domain):      def getFunctionSpace(self,domain):
72         """         """
73         @brief defines the FunctionSpace of the coefficient on the domain         defines the FunctionSpace of the coefficient on the domain
74
75         @param domain         @param domain:
76         """         """
77         if self.what==self.INTERIOR: return escript.Function(domain)         if self.what==self.INTERIOR: return escript.Function(domain)
78         elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)         elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)
79         elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)         elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)
80         elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)         elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)
81
82        def getValue(self):
83           """
84           returns the value of the coefficient:
85           """
86           return self.value
87
88        def setValue(self,newValue):
89           """
90           set the value of the coefficient to new value
91           """
92           self.value=newValue
93
94      def isAlteringOperator(self):      def isAlteringOperator(self):
95          """          """
96      @brief return true if the operator of the PDE is changed when the coefficient is changed      return true if the operator of the PDE is changed when the coefficient is changed
97      """      """
98          if self.altering==self.OPERATOR or self.altering==self.BOTH:          if self.altering==self.OPERATOR or self.altering==self.BOTH:
99              return not None              return not None
# Line 119  class PDECoefficientType: Line 102  class PDECoefficientType:
102
103      def isAlteringRightHandSide(self):      def isAlteringRightHandSide(self):
104          """          """
105      @brief return true if the right hand side of the PDE is changed when the coefficient is changed      return true if the right hand side of the PDE is changed when the coefficient is changed
106      """      """
107          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
108              return not None              return not None
# Line 128  class PDECoefficientType: Line 111  class PDECoefficientType:
111
112      def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3):      def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3):
113         """         """
114         @brief tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim         tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim
115
116         @param shape         @param shape:
117         @param dim         @param dim:
118         """         """
119         if len(shape)>0:         if len(shape)>0:
120             num=max(shape)+1             num=max(shape)+1
# Line 152  class PDECoefficientType: Line 135  class PDECoefficientType:
135
136      def buildShape(self,e=1,u=1,dim=3):      def buildShape(self,e=1,u=1,dim=3):
137          """          """
138      @brief builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim      builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim
139
140      @param e      @param e:
141      @param u      @param u:
142      @param dim      @param dim:
143      """      """
144          s=()          s=()
145          for i in self.pattern:          for i in self.pattern:
# Line 168  class PDECoefficientType: Line 151  class PDECoefficientType:
151                  s=s+(dim,)                  s=s+(dim,)
152          return s          return s
153
_PDECoefficientTypes={
"A"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM,PDECoefficientType.SOLUTION,PDECoefficientType.DIM),PDECoefficientType.OPERATOR),
"B"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),
"C"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION,PDECoefficientType.DIM),PDECoefficientType.OPERATOR),
"D"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),
"X"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM),PDECoefficientType.RIGHTHANDSIDE),
"Y"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),
"d"         : PDECoefficientType(PDECoefficientType.BOUNDARY,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),
"y"         : PDECoefficientType(PDECoefficientType.BOUNDARY,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),
"d_contact" : PDECoefficientType(PDECoefficientType.CONTACT,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),
"y_contact" : PDECoefficientType(PDECoefficientType.CONTACT,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),
"r"         : PDECoefficientType(PDECoefficientType.CONTINUOUS,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),
"q"         : PDECoefficientType(PDECoefficientType.CONTINUOUS,(PDECoefficientType.SOLUTION,),PDECoefficientType.BOTH),
}

154  class LinearPDE:  class LinearPDE:
155     """     """
156     @brief Class to define a linear PDE     Class to handle a linear PDE
157
158     class to define a linear PDE of the form     class to define a linear PDE of the form
159
160       \f[
161       -(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i       -(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i
162       \f]
163
164       with boundary conditons:     with boundary conditons:
165
166          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i     \f[
167       n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
168       \f]
169
170      and contact conditions     and contact conditions
171
172          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i     \f[
173       n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i
174       \f]
175
176      and constraints:     and constraints:
177
178           u_i=r_i where q_i>0     \f[
180       \f]
181
182     """     """
183       TOL=1.e-13
184     DEFAULT_METHOD=util.DEFAULT_METHOD     DEFAULT_METHOD=util.DEFAULT_METHOD
185     DIRECT=util.DIRECT     DIRECT=util.DIRECT
186     CHOLEVSKY=util.CHOLEVSKY     CHOLEVSKY=util.CHOLEVSKY
# Line 214  class LinearPDE: Line 191  class LinearPDE:
191     SSOR=util.SSOR     SSOR=util.SSOR
192     GMRES=util.GMRES     GMRES=util.GMRES
193     PRES20=util.PRES20     PRES20=util.PRES20
194       ILU0=util.ILU0
195       JACOBI=util.JACOBI
196
197     def __init__(self,domain,numEquations=0,numSolutions=0):     def __init__(self,domain,numEquations=0,numSolutions=0):
198       """       """
199       @brief initializes a new linear PDE.       initializes a new linear PDE.
200
201       @param args       @param args:
202       """       """
203         # COEFFICIENTS can be overwritten by subclasses:
204         self.COEFFICIENTS={
205           "A"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
206           "B"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
207           "C"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
208           "D"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
209           "X"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),
210           "Y"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
211           "d"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
212           "y"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
213           "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
214           "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
215           "r"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
216           "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)}
217
218       # initialize attributes       # initialize attributes
219       self.__debug=None       self.__debug=None
# Line 246  class LinearPDE: Line 239  class LinearPDE:
239       self.__sym=False       self.__sym=False
240       self.__lumping=False       self.__lumping=False
241
242       def createCoefficient(self, name):
243         """
244         create a data object corresponding to coefficient name
245         @param name:
246         """
247         return escript.Data(shape = getShapeOfCoefficient(name), \
248                             what = getFunctionSpaceForCoefficient(name))
249
250       def __del__(self):
251         pass
252
253     def getCoefficient(self,name):     def getCoefficient(self,name):
254       """       """
255       @brief return the value of the coefficient name       return the value of the parameter name
256
257       @param name       @param name:
258       """       """
259       return self.__coefficient[name]       return self.COEFFICIENTS[name].getValue()
260
261     def setValue(self,**coefficients):     def getCoefficientOfPDE(self,name):
262         """
263         return the value of the coefficient name of the general PDE.
264         This method is called by the assembling routine it can be
265         overwritten to map coefficients of a particualr PDE to the general PDE.
266
267         @param name:
268         """
269         return self.getCoefficient(name)
270
271       def hasCoefficient(self,name):
272        """        """
273        @brief sets new values to coefficients        return true if name is the name of a coefficient
274
275        @param coefficients        @param name:
276        """        """
277        self._setValue(**coefficients)        return self.COEFFICIENTS.has_key(name)

278
279     def _setValue(self,**coefficients):     def getFunctionSpaceForEquation(self):
280         """
281         return true if the test functions should use reduced order
282         """
283         return self.__row_function_space
284
285       def getFunctionSpaceForSolution(self):
286         """
287         return true if the interpolation of the solution should use reduced order
288         """
289         return self.__column_function_space
290
291       def setValue(self,**coefficients):
292        """        """
293        @brief sets new values to coefficients        sets new values to coefficients
294
295        @param coefficients        @param coefficients:
296        """        """
297          self.__setValue(**coefficients)
298
# get the dictionary of the coefficinets been altered:
alteredCoefficients={}
for i,d in coefficients.iteritems():
if self.hasCoefficient(i):
if d == None:
alteredCoefficients[i]=escript.Data()
elif isinstance(d,escript.Data):
if d.isEmpty():
alteredCoefficients[i]=escript.Data()
else:
alteredCoefficients[i]=escript.Data(d,self.getFunctionSpaceOfCoefficient(i))
else:
alteredCoefficients[i]=escript.Data(d,self.getFunctionSpaceOfCoefficient(i))
else:
raise ValueError,"Attempt to set undefined coefficient %s"%i
# if numEquations and numSolutions is undefined we try identify their values based on the coefficients:
if self.__numEquations<1 or self.__numSolutions<1:
numEquations,numSolutions=identifyNumEquationsAndSolutions(self.getDomain().getDim(),alteredCoefficients)
if self.__numEquations<1 and numEquations>0: self.__numEquations=numEquations
if self.__numSolutions<1 and numSolutions>0: self.__numSolutions=numSolutions
if self.debug() and self.__numEquations>0: print "PDE Debug: identified number of equations is ",self.__numEquations
if self.debug() and self.__numSolutions>0: print "PDE Debug: identified number of solutions is ",self.__numSolutions

# now we check the shape of the coefficient if numEquations and numSolutions are set:
if  self.__numEquations>0 and  self.__numSolutions>0:
for i in self.__coefficient.iterkeys():
if alteredCoefficients.has_key(i) and not alteredCoefficients[i].isEmpty():
if not self.getShapeOfCoefficient(i)==alteredCoefficients[i].getShape():
raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),alteredCoefficients[i].getShape())
else:
if not self.__coefficient[i].isEmpty():
if not self.getShapeOfCoefficient(i)==self.__coefficient[i].getShape():
raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),self.__coefficient[i].getShape())
# overwrite new values:
for i,d in alteredCoefficients.iteritems():
if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i
self.__coefficient[i]=d
self.alteredCoefficient(i)

# reset the HomogeneousConstraintFlag:
self.__setHomogeneousConstraintFlag()
if not "q" in alteredCoefficients and not self.__homogeneous_constraint: self.__rebuildSystem()
299
300     def cleanCoefficients(self):     def cleanCoefficients(self):
301       """       """
302       @brief resets all coefficients to default values.       resets all coefficients to default values.
303       """       """
304       self.__coefficient={}       for i in self.COEFFICIENTS.iterkeys():
305       for i in _PDECoefficientTypes.iterkeys():           self.COEFFICIENTS[i].resetValue()
306           self.__coefficient[i]=escript.Data()
307       def createNewCoefficient(self,name):
308         """
309         returns a new coefficient appropriate for coefficient name:
310         """
311         return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
312
313
314     def getShapeOfCoefficient(self,name):     def getShapeOfCoefficient(self,name):
315       """       """
316       @brief return the shape of the coefficient name       return the shape of the coefficient name
317
318       @param name       @param name:
319       """       """
320       if self.hasCoefficient(name):       if self.hasCoefficient(name):
321          return _PDECoefficientTypes[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())          return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())
322       else:       else:
323          raise ValueError,"Solution coefficient %s requested"%name          raise ValueError,"Solution coefficient %s requested"%name
324
325     def getFunctionSpaceOfCoefficient(self,name):     def getFunctionSpaceForCoefficient(self,name):
326       """       """
327       @brief return the atoms of the coefficient name       return the atoms of the coefficient name
328
329       @param name       @param name:
330       """       """
331       if self.hasCoefficient(name):       if self.hasCoefficient(name):
332          return _PDECoefficientTypes[name].getFunctionSpace(self.getDomain())          return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
333       else:       else:
334          raise ValueError,"Solution coefficient %s requested"%name          raise ValueError,"Solution coefficient %s requested"%name
335
336     def alteredCoefficient(self,name):     def alteredCoefficient(self,name):
337       """       """
338       @brief annonced that coefficient name has been changed       announce that coefficient name has been changed
339
340       @param name       @param name:
341       """       """
342       if self.hasCoefficient(name):       if self.hasCoefficient(name):
343          if _PDECoefficientTypes[name].isAlteringOperator(): self.__rebuildOperator()          if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator()
344          if _PDECoefficientTypes[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()          if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()
345       else:       else:
346          raise ValueError,"Solution coefficient %s requested"%name          raise ValueError,"unknown coefficient %s requested"%name

def __setHomogeneousConstraintFlag(self):
"""
@brief checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.
"""
self.__homogeneous_constraint=True
q=self.getCoefficient("q")
r=self.getCoefficient("r")
if not q.isEmpty() and not r.isEmpty():
print (q*r).Lsup(), 1.e-13*r.Lsup()
if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False
if self.debug():
if self.__homogeneous_constraint:
print "PDE Debug: Constraints are homogeneous."
else:
print "PDE Debug: Constraints are inhomogeneous."

def hasCoefficient(self,name):
"""
@brief return true if name is the name of a coefficient

@param name
"""
return self.__coefficient.has_key(name)

def getFunctionSpaceForEquation(self):
"""
@brief return true if the test functions should use reduced order
"""
return self.__row_function_space

def getFunctionSpaceForSolution(self):
"""
@brief return true if the interpolation of the solution should use reduced order
"""
return self.__column_function_space
347
348     # ===== debug ==============================================================     # ===== debug ==============================================================
349     def setDebugOn(self):     def setDebugOn(self):
350         """         """
@brief
351         """         """
352         self.__debug=not None         self.__debug=not None
353
354     def setDebugOff(self):     def setDebugOff(self):
355         """         """
@brief
356         """         """
357         self.__debug=None         self.__debug=None
358
359     def debug(self):     def debug(self):
360         """         """
361         @brief returns true if the PDE is in the debug mode         returns true if the PDE is in the debug mode
362         """         """
363         return self.__debug         return self.__debug
364
365     #===== Lumping ===========================     #===== Lumping ===========================
366     def setLumpingOn(self):     def setLumpingOn(self):
367        """        """
368        @brief indicates to use matrix lumping        indicates to use matrix lumping
369        """        """
370        if not self.isUsingLumping():        if not self.isUsingLumping():
raise SystemError,"Lumping is not working yet! Talk to the experts"
371           if self.debug() : print "PDE Debug: lumping is set on"           if self.debug() : print "PDE Debug: lumping is set on"
372           self.__rebuildOperator()           self.__rebuildOperator()
373           self.__lumping=True           self.__lumping=True
374
375     def setLumpingOff(self):     def setLumpingOff(self):
376        """        """
377        @brief switches off matrix lumping        switches off matrix lumping
378        """        """
379        if self.isUsingLumping():        if self.isUsingLumping():
380           if self.debug() : print "PDE Debug: lumping is set off"           if self.debug() : print "PDE Debug: lumping is set off"
# Line 433  class LinearPDE: Line 383  class LinearPDE:
383
384     def setLumping(self,flag=False):     def setLumping(self,flag=False):
385        """        """
386        @brief set the matrix lumping flag to flag        set the matrix lumping flag to flag
387        """        """
388        if flag:        if flag:
389           self.setLumpingOn()           self.setLumpingOn()
# Line 442  class LinearPDE: Line 392  class LinearPDE:
392
393     def isUsingLumping(self):     def isUsingLumping(self):
394        """        """
395        @brief
396        """        """
397        return self.__lumping        return self.__lumping
398
400     def setSolverMethod(self,solver=util.DEFAULT_METHOD):     def setSolverMethod(self,solver=util.DEFAULT_METHOD):
401         """         """
402         @brief sets a new solver         sets a new solver
403         """         """
404         if not solver==self.getSolverMethod():         if not solver==self.getSolverMethod():
405             self.__solver_method=solver             self.__solver_method=solver
# Line 458  class LinearPDE: Line 408  class LinearPDE:
408
409     def getSolverMethod(self):     def getSolverMethod(self):
410         """         """
411         @brief returns the solver method         returns the solver method
412         """         """
413         return self.__solver_method         return self.__solver_method
414
416     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
417         """         """
418         @brief resets the tolerance to tol.         resets the tolerance to tol.
419         """         """
420         if not tol>0:         if not tol>0:
421             raise ValueException,"Tolerance as to be positive"             raise ValueException,"Tolerance as to be positive"
# Line 475  class LinearPDE: Line 425  class LinearPDE:
425         return         return
426     def getTolerance(self):     def getTolerance(self):
427         """         """
428         @brief returns the tolerance set for the solution         returns the tolerance set for the solution
429         """         """
430         return self.__tolerance         return self.__tolerance
431
432     #===== symmetry  flag ==========================     #===== symmetry  flag ==========================
433     def isSymmetric(self):     def isSymmetric(self):
434        """        """
435        @brief returns true is the operator is considered to be symmetric        returns true is the operator is considered to be symmetric
436        """        """
437        return self.__sym        return self.__sym
438
439     def setSymmetryOn(self):     def setSymmetryOn(self):
440        """        """
441        @brief sets the symmetry flag to true        sets the symmetry flag to true
442        """        """
443        if not self.isSymmetric():        if not self.isSymmetric():
444           if self.debug() : print "PDE Debug: Operator is set to be symmetric"           if self.debug() : print "PDE Debug: Operator is set to be symmetric"
# Line 497  class LinearPDE: Line 447  class LinearPDE:
447
448     def setSymmetryOff(self):     def setSymmetryOff(self):
449        """        """
450        @brief sets the symmetry flag to false        sets the symmetry flag to false
451        """        """
452        if self.isSymmetric():        if self.isSymmetric():
453           if self.debug() : print "PDE Debug: Operator is set to be unsymmetric"           if self.debug() : print "PDE Debug: Operator is set to be unsymmetric"
# Line 506  class LinearPDE: Line 456  class LinearPDE:
456
457     def setSymmetryTo(self,flag=False):     def setSymmetryTo(self,flag=False):
458       """       """
459       @brief sets the symmetry flag to flag       sets the symmetry flag to flag
460
461       @param flag       @param flag:
462       """       """
463       if flag:       if flag:
464          self.setSymmetryOn()          self.setSymmetryOn()
# Line 518  class LinearPDE: Line 468  class LinearPDE:
468     #===== order reduction ==========================     #===== order reduction ==========================
469     def setReducedOrderOn(self):     def setReducedOrderOn(self):
470       """       """
471       @brief switches to on reduced order       switches to on reduced order
472       """       """
473       self.setReducedOrderForSolutionOn()       self.setReducedOrderForSolutionOn()
474       self.setReducedOrderForEquationOn()       self.setReducedOrderForEquationOn()
475
476     def setReducedOrderOff(self):     def setReducedOrderOff(self):
477       """       """
478       @brief switches to full order       switches to full order
479       """       """
480       self.setReducedOrderForSolutionOff()       self.setReducedOrderForSolutionOff()
481       self.setReducedOrderForEquationOff()       self.setReducedOrderForEquationOff()
482
483     def setReducedOrderTo(self,flag=False):     def setReducedOrderTo(self,flag=False):
484       """       """
485       @brief sets order according to flag       sets order according to flag
486
487       @param flag       @param flag:
488       """       """
489       self.setReducedOrderForSolutionTo(flag)       self.setReducedOrderForSolutionTo(flag)
490       self.setReducedOrderForEquationTo(flag)       self.setReducedOrderForEquationTo(flag)
# Line 543  class LinearPDE: Line 493  class LinearPDE:
493     #===== order reduction solution ==========================     #===== order reduction solution ==========================
494     def setReducedOrderForSolutionOn(self):     def setReducedOrderForSolutionOn(self):
495       """       """
496       @brief switches to reduced order to interpolate solution       switches to reduced order to interpolate solution
497       """       """
498       new_fs=escript.ReducedSolution(self.getDomain())       new_fs=escript.ReducedSolution(self.getDomain())
499       if self.getFunctionSpaceForSolution()!=new_fs:       if self.getFunctionSpaceForSolution()!=new_fs:
# Line 553  class LinearPDE: Line 503  class LinearPDE:
503
504     def setReducedOrderForSolutionOff(self):     def setReducedOrderForSolutionOff(self):
505       """       """
506       @brief switches to full order to interpolate solution       switches to full order to interpolate solution
507       """       """
508       new_fs=escript.Solution(self.getDomain())       new_fs=escript.Solution(self.getDomain())
509       if self.getFunctionSpaceForSolution()!=new_fs:       if self.getFunctionSpaceForSolution()!=new_fs:
# Line 563  class LinearPDE: Line 513  class LinearPDE:
513
514     def setReducedOrderForSolutionTo(self,flag=False):     def setReducedOrderForSolutionTo(self,flag=False):
515       """       """
516       @brief sets order for test functions according to flag       sets order for test functions according to flag
517
518       @param flag       @param flag:
519       """       """
520       if flag:       if flag:
521          self.setReducedOrderForSolutionOn()          self.setReducedOrderForSolutionOn()
# Line 575  class LinearPDE: Line 525  class LinearPDE:
525     #===== order reduction equation ==========================     #===== order reduction equation ==========================
526     def setReducedOrderForEquationOn(self):     def setReducedOrderForEquationOn(self):
527       """       """
528       @brief switches to reduced order for test functions       switches to reduced order for test functions
529       """       """
530       new_fs=escript.ReducedSolution(self.getDomain())       new_fs=escript.ReducedSolution(self.getDomain())
531       if self.getFunctionSpaceForEquation()!=new_fs:       if self.getFunctionSpaceForEquation()!=new_fs:
# Line 585  class LinearPDE: Line 535  class LinearPDE:
535
536     def setReducedOrderForEquationOff(self):     def setReducedOrderForEquationOff(self):
537       """       """
538       @brief switches to full order for test functions       switches to full order for test functions
539       """       """
540       new_fs=escript.Solution(self.getDomain())       new_fs=escript.Solution(self.getDomain())
541       if self.getFunctionSpaceForEquation()!=new_fs:       if self.getFunctionSpaceForEquation()!=new_fs:
# Line 595  class LinearPDE: Line 545  class LinearPDE:
545
546     def setReducedOrderForEquationTo(self,flag=False):     def setReducedOrderForEquationTo(self,flag=False):
547       """       """
548       @brief sets order for test functions according to flag       sets order for test functions according to flag
549
550       @param flag       @param flag:
551       """       """
552       if flag:       if flag:
553          self.setReducedOrderForEquationOn()          self.setReducedOrderForEquationOn()
554       else:       else:
555          self.setReducedOrderForEquationOff()          self.setReducedOrderForEquationOff()
556

557     # ==== initialization =====================================================================     # ==== initialization =====================================================================
558     def __makeNewOperator(self):     def __getNewOperator(self):
559         """         """
@brief
560         """         """
561         return self.getDomain().newOperator( \         return self.getDomain().newOperator( \
562                             self.getNumEquations(), \                             self.getNumEquations(), \
# Line 617  class LinearPDE: Line 565  class LinearPDE:
565                             self.getFunctionSpaceForSolution(), \                             self.getFunctionSpaceForSolution(), \
566                             self.__matrix_type)                             self.__matrix_type)
567
568     def __makeNewRightHandSide(self):     def __makeFreshRightHandSide(self):
569         """         """
@brief
570         """         """
571         return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)         if self.debug() : print "PDE Debug: New right hand side allocated"
572           if self.getNumEquations()>1:
573               self.__righthandside=escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
574           else:
575               self.__righthandside=escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
576           return self.__righthandside
577
578     def __makeNewSolution(self):     def __getNewSolution(self):
579         """         """
@brief
580         """         """
581         return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)         if self.debug() : print "PDE Debug: New right hand side allocated"
582           if self.getNumSolutions()>1:
583               return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
584           else:
585               return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
586
587     def __getFreshOperator(self):     def __makeFreshOperator(self):
588         """         """
@brief
589         """         """
590         if self.__operator.isEmpty():         if self.__operator.isEmpty():
591             self.__operator=self.__makeNewOperator()             self.__operator=self.__getNewOperator()
592             if self.debug() : print "PDE Debug: New operator allocated"             if self.debug() : print "PDE Debug: New operator allocated"
593         else:         else:
594             self.__operator.setValue(0.)             self.__operator.setValue(0.)
595               self.__operator.resetSolver()
596             if self.debug() : print "PDE Debug: Operator reset to zero"             if self.debug() : print "PDE Debug: Operator reset to zero"
597         return self.__operator         return self.__operator
598
599     def __getFreshRightHandSide(self):     #============ some serivice functions  =====================================================
600       def getDomain(self):
601         """
602         returns the domain of the PDE
603         """
604         return self.__domain
605
606       def getDim(self):
607         """
608         returns the spatial dimension of the PDE
609         """
610         return self.getDomain().getDim()
611
612       def getNumEquations(self):
613         """
614         returns the number of equations
615         """
616         if self.__numEquations>0:
617             return self.__numEquations
618         else:
619             raise ValueError,"Number of equations is undefined. Please specify argument numEquations."
620
621       def getNumSolutions(self):
622         """
623         returns the number of unknowns
624         """
625         if self.__numSolutions>0:
626            return self.__numSolutions
627         else:
628            raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."
629
630
631       def checkSymmetry(self,verbose=True):
632          """
633          returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.
634          """
635          verbose=verbose or self.debug()
636          out=True
637          if self.getNumSolutions()!=self.getNumEquations():
638             if verbose: print "non-symmetric PDE because of different number of equations and solutions"
639             out=False
640          else:
641             A=self.getCoefficientOfPDE("A")
642             if not A.isEmpty():
643                tol=util.Lsup(A)*self.TOL
644                if self.getNumSolutions()>1:
645                   for i in range(self.getNumEquations()):
646                      for j in range(self.getDim()):
647                         for k in range(self.getNumSolutions()):
648                            for l in range(self.getDim()):
649                                if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
650                                   if verbose: print "non-symmetric PDE because A[%d,%d,%d,%d]!=A[%d,%d,%d,%d]"%(i,j,k,l,k,l,i,j)
651                                   out=False
652                else:
653                   for j in range(self.getDim()):
654                      for l in range(self.getDim()):
655                         if util.Lsup(A[j,l]-A[l,j])>tol:
656                            if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)
657                            out=False
658             B=self.getCoefficientOfPDE("B")
659             C=self.getCoefficientOfPDE("C")
660             if B.isEmpty() and not C.isEmpty():
661                if verbose: print "non-symmetric PDE because B is not present but C is"
662                out=False
663             elif not B.isEmpty() and C.isEmpty():
664                if verbose: print "non-symmetric PDE because C is not present but B is"
665                out=False
666             elif not B.isEmpty() and not C.isEmpty():
667                tol=(util.Lsup(B)+util.Lsup(C))*self.TOL/2.
668                if self.getNumSolutions()>1:
669                   for i in range(self.getNumEquations()):
670                       for j in range(self.getDim()):
671                          for k in range(self.getNumSolutions()):
672                             if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
673                                  if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
674                                  out=False
675                else:
676                   for j in range(self.getDim()):
677                      if util.Lsup(B[j]-C[j])>tol:
678                         if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)
679                         out=False
680             if self.getNumSolutions()>1:
681               D=self.getCoefficientOfPDE("D")
682               if not D.isEmpty():
683                 tol=util.Lsup(D)*self.TOL
684                 for i in range(self.getNumEquations()):
685                    for k in range(self.getNumSolutions()):
686                      if util.Lsup(D[i,k]-D[k,i])>tol:
687                          if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)
688                          out=False
689
690          return out
691
692       def getFlux(self,u):
693         """         """
694         @brief         returns the flux J_ij for a given u
695
696           \f[
697           J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}
698           \f]
699
700           @param u: argument of the operator
701         """         """
702         if self.__righthandside.isEmpty():         raise SystemError,"getFlux is not implemented yet"
703             self.__righthandside=self.__makeNewRightHandSide()         return None
704             if self.debug() : print "PDE Debug: New right hand side allocated"
705         else:     def applyOperator(self,u):
706             print "fix self.__righthandside*=0"         """
707             self.__righthandside*=0.         applies the operator of the PDE to a given solution u in weak from
708             if self.debug() : print "PDE Debug: Right hand side reset to zero"
709         return  self.__righthandside         @param u: argument of the operator
710           """
711           return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
712
713       def getResidual(self,u):
714           """
715           return the residual of u in the weak from
716
717           @param u:
718           """
719           return self.applyOperator(u)-self.getRightHandSide()
720
721       def __setValue(self,**coefficients):
722          """
723          sets new values to coefficient
724
725          @param coefficients:
726          """
727          # check if the coefficients are  legal:
728          for i in coefficients.iterkeys():
729             if not self.hasCoefficient(i):
730                raise ValueError,"Attempt to set unknown coefficient %s"%i
731          # if the number of unknowns or equations is still unknown we try to estimate them:
732          if self.__numEquations<1 or self.__numSolutions<1:
733             for i,d in coefficients.iteritems():
734                if hasattr(d,"shape"):
735                    s=d.shape
736                elif hasattr(d,"getShape"):
737                    s=d.getShape()
738                else:
739                    s=numarray.array(d).shape
740                if s!=None:
741                    # get number of equations and number of unknowns:
742                    res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(s,self.getDim())
743                    if res==None:
744                        raise ValueError,"Illegal shape %s of coefficient %s"%(s,i)
745                    else:
746                        if self.__numEquations<1: self.__numEquations=res[0]
747                        if self.__numSolutions<1: self.__numSolutions=res[1]
748          if self.__numEquations<1: raise ValueError,"unidententified number of equations"
749          if self.__numSolutions<1: raise ValueError,"unidententified number of solutions"
750          # now we check the shape of the coefficient if numEquations and numSolutions are set:
751          for i,d in coefficients.iteritems():
752            if d==None:
753                 d2=escript.Data()
754            elif isinstance(d,escript.Data):
755                 if d.isEmpty():
756                    d2=d
757                 else:
758                    d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))
759            else:
760                  d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))
761            if not d2.isEmpty():
762               if not self.getShapeOfCoefficient(i)==d2.getShape():
763                   raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),d2.getShape())
764            # overwrite new values:
765            if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i
766            self.COEFFICIENTS[i].setValue(d2)
767            self.alteredCoefficient(i)
768
769          # reset the HomogeneousConstraintFlag:
770          self.__setHomogeneousConstraintFlag()
771          if len(coefficients)>0 and not self.isUsingLumping() and not self.__homogeneous_constraint: self.__rebuildSystem()
772
773       def __setHomogeneousConstraintFlag(self):
774          """
775          checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.
776          """
777          self.__homogeneous_constraint=True
778          q=self.getCoefficientOfPDE("q")
779          r=self.getCoefficientOfPDE("r")
780          if not q.isEmpty() and not r.isEmpty():
781             if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False
782          if self.debug():
783               if self.__homogeneous_constraint:
784                   print "PDE Debug: Constraints are homogeneous."
785               else:
786                   print "PDE Debug: Constraints are inhomogeneous."
787
788
789     # ==== rebuild switches =====================================================================     # ==== rebuild switches =====================================================================
790     def __rebuildSolution(self,deep=False):     def __rebuildSolution(self,deep=False):
791         """         """
792         @brief indicates the PDE has to be reolved if the solution is requested         indicates the PDE has to be reolved if the solution is requested
793         """         """
794         if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."         if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."
795         self.__solution_isValid=False         self.__solution_isValid=False
796         if deep: self.__solution=escript.Data(deep)         if deep: self.__solution=escript.Data()
797
798
799     def __rebuildOperator(self,deep=False):     def __rebuildOperator(self,deep=False):
800         """         """
801         @brief indicates the operator has to be rebuilt next time it is used         indicates the operator has to be rebuilt next time it is used
802         """         """
803         if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."         if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."
804         self.__rebuildSolution(deep)         self.__rebuildSolution(deep)
# Line 675  class LinearPDE: Line 807  class LinearPDE:
807
808     def __rebuildRightHandSide(self,deep=False):     def __rebuildRightHandSide(self,deep=False):
809         """         """
810         @brief indicates the right hand side has to be rebuild next time it is used         indicates the right hand side has to be rebuild next time it is used
811         """         """
812         if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."         if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."
813         self.__rebuildSolution(deep)         self.__rebuildSolution(deep)
814         self.__righthandside_isValid=False         self.__righthandside_isValid=False
if not self.__homogeneous_constraint: self.__rebuildOperator()
815         if deep: self.__righthandside=escript.Data()         if deep: self.__righthandside=escript.Data()
816
817     def __rebuildSystem(self,deep=False):     def __rebuildSystem(self,deep=False):
818         """         """
819         @brief annonced that all coefficient name has been changed         annonced that all coefficient name has been changed
820         """         """
821         self.__rebuildSolution(deep)         self.__rebuildSolution(deep)
822         self.__rebuildOperator(deep)         self.__rebuildOperator(deep)
# Line 693  class LinearPDE: Line 824  class LinearPDE:
824
825     def __checkMatrixType(self):     def __checkMatrixType(self):
826       """       """
827       @brief reassess the matrix type and, if needed, initiates an operator rebuild       reassess the matrix type and, if needed, initiates an operator rebuild
828       """       """
829       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())
830       if not new_matrix_type==self.__matrix_type:       if not new_matrix_type==self.__matrix_type:
# Line 702  class LinearPDE: Line 833  class LinearPDE:
833           self.__rebuildOperator(deep=True)           self.__rebuildOperator(deep=True)
834
835     #============ assembling =======================================================     #============ assembling =======================================================
836     def __copyConstraint(self,u):     def __copyConstraint(self):
837        """        """
838        @brief copies the constrint condition into u        copies the constrint condition into u
839        """        """
840        q=self.getCoefficient("q")        if not self.__righthandside.isEmpty():
841        r=self.getCoefficient("r")           q=self.getCoefficientOfPDE("q")
842        if not q.isEmpty():           r=self.getCoefficientOfPDE("r")
843            if r.isEmpty():           if not q.isEmpty():
844               r2=escript.Data(0,u.getShape(),u.getFunctionSpace())               if r.isEmpty():
845            else:                  r2=escript.Data(0,self.__righthandside.getShape(),self.__righthandside.getFunctionSpace())
846               r2=escript.Data(r,u.getFunctionSpace())               else:
849
850     def __applyConstraint(self,rhs_update=True):     def __applyConstraint(self):
851         """         """
852         @brief applies the constraints  defined by q and r to the system         applies the constraints defined by q and r to the system
853         """         """
854         q=self.getCoefficient("q")         q=self.getCoefficientOfPDE("q")
855         r=self.getCoefficient("r")         r=self.getCoefficientOfPDE("r")
856         if not q.isEmpty() and not self.__operator.isEmpty():         if not q.isEmpty() and not self.__operator.isEmpty():
857            # q is the row and column mask to indicate where constraints are set:            # q is the row and column mask to indicate where constraints are set:
858            row_q=escript.Data(q,self.getFunctionSpaceForEquation())            row_q=escript.Data(q,self.getFunctionSpaceForEquation())
859            col_q=escript.Data(q,self.getFunctionSpaceForSolution())            col_q=escript.Data(q,self.getFunctionSpaceForSolution())
860            u=self.__makeNewSolution()            u=self.__getNewSolution()
861            if r.isEmpty():            if r.isEmpty():
862               r_s=self.__makeNewSolution()               r_s=self.__getNewSolution()
863            else:            else:
864               r_s=escript.Data(r,self.getFunctionSpaceForSolution())               r_s=escript.Data(r,self.getFunctionSpaceForSolution())
866            if not self.__righthandside.isEmpty() and rhs_update: self.__righthandside-=self.__operator*u            if self.isUsingLumping():
868              else:
869                 if not self.__righthandside.isEmpty(): self.__righthandside-=self.__operator*u
870                 self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
871
872     def getOperator(self):     def getSystem(self):
873         """         """
874         @brief returns the operator of the PDE         return the operator and right hand side of the PDE
875         """         """
876         if not self.__operator_isValid:         if not self.__operator_isValid or not self.__righthandside_isValid:
877             # some Constraints are applying for a lumpled stifness matrix:            if self.isUsingLumping():
878             if self.isUsingLumping():                if not self.__operator_isValid:
879                if self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
880                         raise TypeError,"Lumped matrix requires same order for equations and unknowns"                         raise TypeError,"Lumped matrix requires same order for equations and unknowns"
881                if not self.getCoefficient("A").isEmpty():                   if not self.getCoefficientOfPDE("A").isEmpty():
882                         raise Warning,"Lumped matrix does not allow coefficient A"                            raise Warning,"Lumped matrix does not allow coefficient A"
883                if not self.getCoefficient("B").isEmpty():                   if not self.getCoefficientOfPDE("B").isEmpty():
884                         raise Warning,"Lumped matrix does not allow coefficient B"                            raise Warning,"Lumped matrix does not allow coefficient B"
885                if not self.getCoefficient("C").isEmpty():                   if not self.getCoefficientOfPDE("C").isEmpty():
886                         raise Warning,"Lumped matrix does not allow coefficient C"                            raise Warning,"Lumped matrix does not allow coefficient C"
887                if not self.getCoefficient("D").isEmpty():                   if self.debug() : print "PDE Debug: New lumped operator is built."
888                         raise Warning,"Lumped matrix does not allow coefficient D"                   mat=self.__getNewOperator()
889                if self.debug() : print "PDE Debug: New lumped operator is built."                   self.getDomain().addPDEToSystem(mat,escript.Data(), \
890                mat=self.__makeNewOperator(self)                             self.getCoefficientOfPDE("A"), \
891             else:                             self.getCoefficientOfPDE("B"), \
892                if self.debug() : print "PDE Debug: New operator is built."                             self.getCoefficientOfPDE("C"), \
893                mat=self.__getFreshOperator()                             self.getCoefficientOfPDE("D"), \
894                               escript.Data(), \
896                          self.getCoefficient("A"), \                             self.getCoefficientOfPDE("d"), \
897                          self.getCoefficient("B"), \                             escript.Data(),\
898                          self.getCoefficient("C"), \                             self.getCoefficientOfPDE("d_contact"), \
899                          self.getCoefficient("D"), \                             escript.Data())
900                          escript.Data(), \                   self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
901                          escript.Data(), \                   self.__applyConstraint()
902                          self.getCoefficient("d"), \                   self.__operator_isValid=True
903                          escript.Data(),\                if not self.__righthandside_isValid:
904                          self.getCoefficient("d_contact"), \                   if self.debug() : print "PDE Debug: New right hand side is built."
906             if self.isUsingLumping():                                 self.getCoefficientOfPDE("X"), \
907                self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceOfSolution(),True)                                 self.getCoefficientOfPDE("Y"),\
908             else:                                 self.getCoefficientOfPDE("y"),\
909                self.__applyConstraint(rhs_update=False)                                 self.getCoefficientOfPDE("y_contact"))
910             self.__operator_isValid=True                   self.__copyConstraint()
911         return self.__operator                   self.__righthandside_isValid=True
912              else:
913     def getRightHandSide(self,ignoreConstraint=False):               if not self.__operator_isValid and not self.__righthandside_isValid:
914                     if self.debug() : print "PDE Debug: New system is built."
916                                   self.getCoefficientOfPDE("A"), \
917                                   self.getCoefficientOfPDE("B"), \
918                                   self.getCoefficientOfPDE("C"), \
919                                   self.getCoefficientOfPDE("D"), \
920                                   self.getCoefficientOfPDE("X"), \
921                                   self.getCoefficientOfPDE("Y"), \
922                                   self.getCoefficientOfPDE("d"), \
923                                   self.getCoefficientOfPDE("y"), \
924                                   self.getCoefficientOfPDE("d_contact"), \
925                                   self.getCoefficientOfPDE("y_contact"))
926                     self.__applyConstraint()
927                     self.__copyConstraint()
928                     self.__operator_isValid=True
929                     self.__righthandside_isValid=True
930                 elif not self.__righthandside_isValid:
931                     if self.debug() : print "PDE Debug: New right hand side is built."
933                                   self.getCoefficientOfPDE("X"), \
934                                   self.getCoefficientOfPDE("Y"),\
935                                   self.getCoefficientOfPDE("y"),\
936                                   self.getCoefficientOfPDE("y_contact"))
937                     self.__copyConstraint()
938                     self.__righthandside_isValid=True
939                 elif not self.__operator_isValid:
940                     if self.debug() : print "PDE Debug: New operator is built."
942                                self.getCoefficientOfPDE("A"), \
943                                self.getCoefficientOfPDE("B"), \
944                                self.getCoefficientOfPDE("C"), \
945                                self.getCoefficientOfPDE("D"), \
946                                escript.Data(), \
947                                escript.Data(), \
948                                self.getCoefficientOfPDE("d"), \
949                                escript.Data(),\
950                                self.getCoefficientOfPDE("d_contact"), \
951                                escript.Data())
952                     self.__applyConstraint()
953                     self.__operator_isValid=True
954           return (self.__operator,self.__righthandside)
955       def getOperator(self):
956         """         """
957         @brief returns the right hand side of the PDE         returns the operator of the PDE

@param ignoreConstraint
958         """         """
959         if not self.__righthandside_isValid:         return self.getSystem()[0]
if self.debug() : print "PDE Debug: New right hand side is built."
self.getCoefficient("X"), \
self.getCoefficient("Y"),\
self.getCoefficient("y"),\
self.getCoefficient("y_contact"))
self.__righthandside_isValid=True
if ignoreConstraint: self.__copyConstraint(self.__righthandside)
return self.__righthandside
960
961     def getSystem(self):     def getRightHandSide(self):
962         """         """
963         @brief         returns the right hand side of the PDE
964         """         """
965         if not self.__operator_isValid and not self.__righthandside_isValid:         return self.getSystem()[1]
if self.debug() : print "PDE Debug: New PDE is built."
if self.isUsingLumping():
self.getRightHandSide(ignoreConstraint=True)
self.getOperator()
else:
self.getCoefficient("A"), \
self.getCoefficient("B"), \
self.getCoefficient("C"), \
self.getCoefficient("D"), \
self.getCoefficient("X"), \
self.getCoefficient("Y"), \
self.getCoefficient("d"), \
self.getCoefficient("y"), \
self.getCoefficient("d_contact"), \
self.getCoefficient("y_contact"))
self.__operator_isValid=True
self.__righthandside_isValid=True
self.__applyConstraint()
self.__copyConstraint(self.__righthandside)
elif not self.__operator_isValid:
self.getOperator()
elif not self.__righthandside_isValid:
self.getRightHandSide()
return (self.__operator,self.__righthandside)
966
967     def solve(self,**options):     def solve(self,**options):
968        """        """
969        @brief solve the PDE        solve the PDE
970
971        @param options        @param options:
972        """        """
973        mat,f=self.getSystem()        mat,f=self.getSystem()
974        if self.isUsingLumping():        if self.isUsingLumping():
975           out=f/mat           out=f/mat
self.__copyConstraint(out)
976        else:        else:
977           options[util.TOLERANCE_KEY]=self.getTolerance()           options[util.TOLERANCE_KEY]=self.getTolerance()
978           options[util.METHOD_KEY]=self.getSolverMethod()           options[util.METHOD_KEY]=self.getSolverMethod()
# Line 843  class LinearPDE: Line 983  class LinearPDE:
983
984     def getSolution(self,**options):     def getSolution(self,**options):
985         """         """
986         @brief returns the solution of the PDE         returns the solution of the PDE
987
988         @param options         @param options:
989         """         """
990         if not self.__solution_isValid:         if not self.__solution_isValid:
991             if self.debug() : print "PDE Debug: PDE is resolved."             if self.debug() : print "PDE Debug: PDE is resolved."
992             self.__solution=self.solve(**options)             self.__solution=self.solve(**options)
993             self.__solution_isValid=True             self.__solution_isValid=True
994         return self.__solution         return self.__solution
#============ some serivice functions  =====================================================
def getDomain(self):
"""
@brief returns the domain of the PDE
"""
return self.__domain
995
def getDim(self):
"""
@brief returns the spatial dimension of the PDE
"""
return self.getDomain().getDim()
996
def getNumEquations(self):
"""
@brief returns the number of equations
"""
if self.__numEquations>0:
return self.__numEquations
else:
raise ValueError,"Number of equations is undefined. Please specify argument numEquations."
997
998     def getNumSolutions(self):  def ELMAN_RAMAGE(P):
999       """       """   """
1000       @brief returns the number of unknowns       return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
1001       """  def SIMPLIFIED_BROOK_HUGHES(P):
1002       if self.__numSolutions>0:       """   """
1003          return self.__numSolutions       c=(P-3.).whereNegative()
1004       else:       return P/6.*c+1./2.*(1.-c)
1005          raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."
1006    def HALF(P):
1007        """ """
1008        return escript.Scalar(0.5,P.getFunctionSpace())
1009
1011       """
1012       Class to handle a linear PDE dominated by advective terms:
1013
1014       class to define a linear PDE of the form
1015
1016     def checkSymmetry(self):     \f[
1017        """     -(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i
1018        @brief returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.     \f]
"""
raise SystemError,"checkSymmetry is not implemented yet"
1019
1020        return None     with boundary conditons:
1021
1022     def getFlux(self,u):     \f[
1023         """     n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
1024         @brief returns the flux J_ij for a given u     \f]
1025
1026              J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}     and contact conditions
1027
1028         @param u argument of the operator     \f[
1029       n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d^{contact}_{ik}[u_k] = - n_j*X_{ij} + y^{contact}_{i}
1030       \f]
1031
1032         """     and constraints:
raise SystemError,"getFlux is not implemented yet"
return None
1033
1034     def applyOperator(self,u):     \f[
1036         @brief applies the operator of the PDE to a given solution u in weak from     \f]
1037       """
1038       def __init__(self,domain,numEquations=0,numSolutions=0,xi=ELMAN_RAMAGE):
1039          LinearPDE.__init__(self,domain,numEquations,numSolutions)
1040          self.__xi=xi
1041          self.__Xi=escript.Data()
1042
1043       def __calculateXi(self,peclet_factor,Z,h):
1044           Z_max=util.Lsup(Z)
1045           if Z_max>0.:
1046              return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.TOL)
1047           else:
1048              return 0.
1049
1050         @param u argument of the operator     def setValue(self,**args):
1051           if "A" in args.keys()   or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data()
1052           self._LinearPDE__setValue(**args)
1053
1054       def getXi(self):
1055          if self.__Xi.isEmpty():
1056             B=self.getCoefficient("B")
1057             C=self.getCoefficient("C")
1058             A=self.getCoefficient("A")
1059             h=self.getDomain().getSize()
1060             self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
1061             if not C.isEmpty() or not B.isEmpty():
1062                if not C.isEmpty() and not B.isEmpty():
1063                    Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
1064                    if self.getNumEquations()>1:
1065                       if self.getNumSolutions()>1:
1066                          for i in range(self.getNumEquations()):
1067                             for k in range(self.getNumSolutions()):
1068                                for l in range(self.getDim()): Z2+=(C[i,k,l]-B[i,l,k])**2
1069                       else:
1070                          for i in range(self.getNumEquations()):
1071                             for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2
1072                    else:
1073                       if self.getNumSolutions()>1:
1074                          for k in range(self.getNumSolutions()):
1075                             for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2
1076                       else:
1077                          for l in range(self.getDim()): Z2+=(C[l]-B[l])**2
1078                    length_of_Z=util.sqrt(Z2)
1079                elif C.isEmpty():
1080                  length_of_Z=util.length(B)
1081                else:
1082                  length_of_Z=util.length(C)
1083
1084         """              Z_max=util.Lsup(length_of_Z)
1085         return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())              if Z_max>0.:
1086                                                                                                                                                                             length_of_A=util.length(A)
1087     def getResidual(self,u):                 A_max=util.Lsup(length_of_A)
1088         """                 if A_max>0:
1089         @brief return the residual of u in the weak from                      inv_A=1./(length_of_A+A_max*self.TOL)
1090                   else:
1091                        inv_A=1./self.TOL
1092                   peclet_number=length_of_Z*h/2*inv_A
1093                   xi=self.__xi(peclet_number)
1094                   self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL)
1095                   print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z)
1096          return self.__Xi
1097
1098
1099       def getCoefficientOfPDE(self,name):
1100         """
1101         return the value of the coefficient name of the general PDE
1102
1103         @param name:
1104         """
1105         if not self.getNumEquations() == self.getNumSolutions():
1106              raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."
1107
1108         if name == "A" :
1109             A=self.getCoefficient("A")
1110             B=self.getCoefficient("B")
1111             C=self.getCoefficient("C")
1112             if B.isEmpty() and C.isEmpty():
1113                Aout=A
1114             else:
1115                if A.isEmpty():
1116                   Aout=self.createNewCoefficient("A")
1117                else:
1118                   Aout=A[:]
1119                Xi=self.getXi()
1120                if self.getNumEquations()>1:
1121                    for i in range(self.getNumEquations()):
1122                       for j in range(self.getDim()):
1123                          for k in range(self.getNumSolutions()):
1124                             for l in range(self.getDim()):
1125                                if not C.isEmpty() and not B.isEmpty():
1126                                   for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*(C[p,i,j]-B[p,j,i])*(C[p,k,l]-B[p,l,k])
1127                                elif C.isEmpty():
1128                                   for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]
1129                                else:
1130                                   for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
1131                else:
1132                    for j in range(self.getDim()):
1133                       for l in range(self.getDim()):
1134                          if not C.isEmpty() and not B.isEmpty():
1135                              Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])
1136                          elif C.isEmpty():
1137                              Aout[j,l]+=Xi*B[j]*B[l]
1138                          else:
1139                              Aout[j,l]+=Xi*C[j]*C[l]
1140             return Aout
1141         elif name == "B" :
1142             B=self.getCoefficient("B")
1143             C=self.getCoefficient("C")
1144             D=self.getCoefficient("D")
1145             if C.isEmpty() or D.isEmpty():
1146                Bout=B
1147             else:
1148                Xi=self.getXi()
1149                if B.isEmpty():
1150                    Bout=self.createNewCoefficient("B")
1151                else:
1152                    Bout=B[:]
1153                if self.getNumEquations()>1:
1154                   for k in range(self.getNumSolutions()):
1155                      for p in range(self.getNumEquations()):
1156                         tmp=Xi*D[p,k]
1157                         for i in range(self.getNumEquations()):
1158                            for j in range(self.getDim()):
1159                               Bout[i,j,k]+=tmp*C[p,i,j]
1160                else:
1161                   tmp=Xi*D
1162                   for j in range(self.getDim()): Bout[j]+=tmp*C[j]
1163             return Bout
1164         elif name == "C" :
1165             B=self.getCoefficient("B")
1166             C=self.getCoefficient("C")
1167             D=self.getCoefficient("D")
1168             if B.isEmpty() or D.isEmpty():
1169                Cout=C
1170             else:
1171                Xi=self.getXi()
1172                if C.isEmpty():
1173                    Cout=self.createNewCoefficient("C")
1174                else:
1175                    Cout=C[:]
1176                if self.getNumEquations()>1:
1177                   for k in range(self.getNumSolutions()):
1178                       for p in range(self.getNumEquations()):
1179                          tmp=Xi*D[p,k]
1180                          for i in range(self.getNumEquations()):
1181                            for l in range(self.getDim()):
1182                                     Cout[i,k,l]+=tmp*B[p,l,i]
1183                else:
1184                   tmp=Xi*D
1185                   for j in range(self.getDim()): Cout[j]+=tmp*B[j]
1186             return Cout
1187         elif name == "D" :
1188             return self.getCoefficient("D")
1189         elif name == "X" :
1190             X=self.getCoefficient("X")
1191             Y=self.getCoefficient("Y")
1192             B=self.getCoefficient("B")
1193             C=self.getCoefficient("C")
1194             if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):
1195                Xout=X
1196             else:
1197                if X.isEmpty():
1198                    Xout=self.createNewCoefficient("X")
1199                else:
1200                    Xout=X[:]
1201                Xi=self.getXi()
1202                if self.getNumEquations()>1:
1203                     for p in range(self.getNumEquations()):
1204                        tmp=Xi*Y[p]
1205                        for i in range(self.getNumEquations()):
1206                           for j in range(self.getDim()):
1207                              if not C.isEmpty() and not B.isEmpty():
1208                                 Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
1209                              elif C.isEmpty():
1210                                 Xout[i,j]-=tmp*B[p,j,i]
1211                              else:
1212                                 Xout[i,j]+=tmp*C[p,i,j]
1213                else:
1214                     tmp=Xi*Y
1215                     for j in range(self.getDim()):
1216                        if not C.isEmpty() and not B.isEmpty():
1217                           Xout[j]+=tmp*(C[j]-B[j])
1218                        elif C.isEmpty():
1219                           Xout[j]-=tmp*B[j]
1220                        else:
1221                           Xout[j]+=tmp*C[j]
1222             return Xout
1223         elif name == "Y" :
1224             return self.getCoefficient("Y")
1225         elif name == "d" :
1226             return self.getCoefficient("d")
1227         elif name == "y" :
1228             return self.getCoefficient("y")
1229         elif name == "d_contact" :
1230             return self.getCoefficient("d_contact")
1231         elif name == "y_contact" :
1232             return self.getCoefficient("y_contact")
1233         elif name == "r" :
1234             return self.getCoefficient("r")
1235         elif name == "q" :
1236             return self.getCoefficient("q")
1237         else:
1238             raise SystemError,"unknown PDE coefficient %s",name
1239
@param u
"""
return self.applyOperator(u)-self.getRightHandSide()
1240
1241  class Poisson(LinearPDE):  class Poisson(LinearPDE):
1242     """     """
1243     @brief Class to define a Poisson equstion problem:     Class to define a Poisson equstion problem:
1244
1245     class to define a linear PDE of the form     class to define a linear PDE of the form
1246                                                                                                                                                                   \f[
1247          -u_{,jj} = f     -u_{,jj} = f
1248                                                                                                                                                                   \f]
1249       with boundary conditons:
1250                                                                                                                                                                   with boundary conditons:
1251          n_j*u_{,j} = 0
1252                                                                                                                                                                   \f[
1253      and constraints:     n_j*u_{,j} = 0
1254                                                                                                                                                                   \f]
1255           u=0 where q>0
1256                                                                                                                                                                   and constraints:
1257
1258       \f[
1260       \f]
1261     """     """
1262
1263     def __init__(self,domain=None,f=escript.Data(),q=escript.Data()):     def __init__(self,domain,f=escript.Data(),q=escript.Data()):
1264         LinearPDE.__init__(self,domain=identifyDomain(domain,{"f":f, "q":q}))         LinearPDE.__init__(self,domain,1,1)
1265         self._setValue(A=numarray.identity(self.getDomain().getDim()))         self.COEFFICIENTS={
1266           "f"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1267           "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}
1268         self.setSymmetryOn()         self.setSymmetryOn()
1269         self.setValue(f,q)         self.setValue(f,q)
1270
1271     def setValue(self,f=escript.Data(),q=escript.Data()):     def setValue(self,f=escript.Data(),q=escript.Data()):
1272         self._setValue(Y=f,q=q)         self._LinearPDE__setValue(f=f,q=q)
1273
1274                                                                                                                                                                 def getCoefficientOfPDE(self,name):
1275         """
1276         return the value of the coefficient name of the general PDE
1277
1278         @param name:
1279         """
1280         if name == "A" :
1281             return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))
1282         elif name == "B" :
1283             return escript.Data()
1284         elif name == "C" :
1285             return escript.Data()
1286         elif name == "D" :
1287             return escript.Data()
1288         elif name == "X" :
1289             return escript.Data()
1290         elif name == "Y" :
1291             return self.getCoefficient("f")
1292         elif name == "d" :
1293             return escript.Data()
1294         elif name == "y" :
1295             return escript.Data()
1296         elif name == "d_contact" :
1297             return escript.Data()
1298         elif name == "y_contact" :
1299             return escript.Data()
1300         elif name == "r" :
1301             return escript.Data()
1302         elif name == "q" :
1303             return self.getCoefficient("q")
1304         else:
1305             raise SystemError,"unknown PDE coefficient %s",name
1306
1307  # $Log$  # $Log$
1308  # Revision 1.3  2004/12/17 07:43:10  jgs  # Revision 1.8  2005/06/09 05:37:59  jgs
1309  # *** empty log message ***  # Merge of development branch back to main trunk on 2005-06-09
1310    #
1311    # Revision 1.7  2005/05/06 04:26:10  jgs
1312    # Merge of development branch back to main trunk on 2005-05-06
1313    #
1314    # Revision 1.1.2.23  2005/05/13 00:55:20  cochrane
1315    # Fixed up some docstrings.  Moved module-level functions to top of file so
1316    # that epydoc and doxygen can pick them up properly.
1317    #
1318    # Revision 1.1.2.22  2005/05/12 11:41:30  gross
1319    # some basic Models have been added
1320    #
1321    # Revision 1.1.2.21  2005/05/12 07:16:12  cochrane
1322    # Moved ELMAN_RAMAGE, SIMPLIFIED_BROOK_HUGHES, and HALF functions to bottom of
1323    # file so that the AdvectivePDE class is picked up by doxygen.  Some
1324    # reformatting of docstrings.  Addition of code to make equations come out
1325    # as proper LaTeX.
1326    #
1327    # Revision 1.1.2.20  2005/04/15 07:09:08  gross
1328    # some problems with functionspace and linearPDEs fixed.
1329    #
1330    # Revision 1.1.2.19  2005/03/04 05:27:07  gross
1331    # bug in SystemPattern fixed.
1332    #
1333    # Revision 1.1.2.18  2005/02/08 06:16:45  gross
1334    # Bugs in AdvectivePDE fixed, AdvectiveTest is stable but more testing is needed
1335    #
1336    # Revision 1.1.2.17  2005/02/08 05:56:19  gross
1337    # Reference Number handling added
1338    #
1339    # Revision 1.1.2.16  2005/02/07 04:41:28  gross
1340    # some function exposed to python to make mesh merging running
1341    #
1342    # Revision 1.1.2.15  2005/02/03 00:14:44  gross
1343    # timeseries add and ESySParameter.py renames esysXML.py for consistence
1344    #
1345    # Revision 1.1.2.14  2005/02/01 06:44:10  gross
1346    # new implementation of AdvectivePDE which now also updates right hand side. systems of PDEs are still not working
1347    #
1348    # Revision 1.1.2.13  2005/01/25 00:47:07  gross
1349    # updates in the documentation
1350    #
1351    # Revision 1.1.2.12  2005/01/12 01:28:04  matt
1352    # Added createCoefficient method for linearPDEs.
1353    #
1354    # Revision 1.1.2.11  2005/01/11 01:55:34  gross
1355    # a problem in linearPDE class fixed
1356    #
1357    # Revision 1.1.2.10  2005/01/07 01:13:29  gross
1358    # some bugs in linearPDE fixed
1359    #
1360    # Revision 1.1.2.9  2005/01/06 06:24:58  gross
1361    # some bugs in slicing fixed
1362    #
1363    # Revision 1.1.2.8  2005/01/05 04:21:40  gross
1364    # FunctionSpace checking/matchig in slicing added
1365    #
1366    # Revision 1.1.2.7  2004/12/29 10:03:41  gross
1367    # bug in setValue fixed
1368    #
1369    # Revision 1.1.2.6  2004/12/29 05:29:59  gross
1370    # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()
1371    #
1372    # Revision 1.1.2.5  2004/12/29 00:18:41  gross
1374    #
1375    # Revision 1.1.2.4  2004/12/24 06:05:41  gross