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

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

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

revision 123 by jgs, Fri Jul 8 04:08:13 2005 UTC revision 154 by jgs, Mon Nov 7 05:51:17 2005 UTC
# Line 3  Line 3 
3  ## @file util.py  ## @file util.py
4    
5  """  """
6  @brief Utility functions for escript  Utility functions for escript
7    
8  TODO for Data:  @todo:
   
   * binary operations @:               (a@b)[:,*]=a[:]@b[:,*] if rank(a)<rank(b)  
                 @=+,-,*,/,**           (a@b)[:]=a[:]@b[:] if rank(a)=rank(b)  
                                        (a@b)[*,:]=a[*,:]@b[:] if rank(a)>rank(b)  
     
   * implementation of outer outer(a,b)[:,*]=a[:]*b[*]  
   * trace: trace(arg,axis0=a0,axis1=a1)(:,&,*)=sum_i trace(:,i,&,i,*) (i are at index a0 and a1)  
9    
10      - binary operations @ (@=+,-,*,/,**)::
11          (a@b)[:,*]=a[:]@b[:,*] if rank(a)<rank(b)
12          (a@b)[:]=a[:]@b[:] if rank(a)=rank(b)
13          (a@b)[*,:]=a[*,:]@b[:] if rank(a)>rank(b)
14      - implementation of outer::
15          outer(a,b)[:,*]=a[:]*b[*]
16      - trace::
17          trace(arg,axis0=a0,axis1=a1)(:,&,*)=sum_i trace(:,i,&,i,*) (i are at index a0 and a1)
18  """  """
19    
20  import numarray  import numarray
21  import escript  import escript
22  #  import symbols
23  #   escript constants (have to be consistent witj utilC.h  import os
 #  
 UNKNOWN=-1  
 EPSILON=1.e-15  
 Pi=numarray.pi  
 # some solver options:  
 NO_REORDERING=0  
 MINIMUM_FILL_IN=1  
 NESTED_DISSECTION=2  
 # solver methods  
 DEFAULT_METHOD=0  
 DIRECT=1  
 CHOLEVSKY=2  
 PCG=3  
 CR=4  
 CGS=5  
 BICGSTAB=6  
 SSOR=7  
 ILU0=8  
 ILUT=9  
 JACOBI=10  
 GMRES=11  
 PRES20=12  
   
 METHOD_KEY="method"  
 SYMMETRY_KEY="symmetric"  
 TOLERANCE_KEY="tolerance"  
   
 # supported file formats:  
 VRML=1  
 PNG=2  
 JPEG=3  
 JPG=3  
 PS=4  
 OOGL=5  
 BMP=6  
 TIFF=7  
 OPENINVENTOR=8  
 RENDERMAN=9  
 PNM=10  
 #===========================================================  
 #  a simple tool box to deal with _differentials of functions  
 #===========================================================  
   
 class Symbol:  
    """symbol class"""  
    def __init__(self,name="symbol",shape=(),dim=3,args=[]):  
        """creates an instance of a symbol of shape shape and spatial dimension dim  
           The symbol may depending on a list of arguments args which  
           may be symbols or other objects. name gives the name of the symbol."""  
        self.__args=args  
        self.__name=name  
        self.__shape=shape  
        if hasattr(dim,"getDim"):  
            self.__dim=dim.getDim()  
        else:      
            self.__dim=dim  
        #  
        self.__cache_val=None  
        self.__cache_argval=None  
   
    def getArgument(self,i):  
        """returns the i-th argument"""  
        return self.__args[i]  
   
    def getDim(self):  
        """returns the spatial dimension of the symbol"""  
        return self.__dim  
   
    def getRank(self):  
        """returns the rank of the symbol"""  
        return len(self.getShape())  
   
    def getShape(self):  
        """returns the shape of the symbol"""  
        return self.__shape  
   
    def getEvaluatedArguments(self,argval):  
        """returns the list of evaluated arguments by subsituting symbol u by argval[u]."""  
        if argval==self.__cache_argval:  
            print "%s: cached value used"%self  
            return self.__cache_val  
        else:  
            out=[]  
            for a  in self.__args:  
               if isinstance(a,Symbol):  
                 out.append(a.eval(argval))  
               else:  
                 out.append(a)  
            self.__cache_argval=argval  
            self.__cache_val=out  
            return out  
   
    def getDifferentiatedArguments(self,arg):  
        """returns the list of the arguments _differentiated by arg"""  
        out=[]  
        for a in self.__args:  
           if isinstance(a,Symbol):  
             out.append(a.diff(arg))  
           else:  
             out.append(0)  
        return out  
   
    def diff(self,arg):  
        """returns the _differention of self by arg."""  
        if self==arg:  
           out=numarray.zeros(tuple(2*list(self.getShape())),numarray.Float)  
           if self.getRank()==0:  
              out=1.  
           elif self.getRank()==1:  
               for i0 in range(self.getShape()[0]):  
                  out[i0,i0]=1.    
           elif self.getRank()==2:  
               for i0 in range(self.getShape()[0]):  
                 for i1 in range(self.getShape()[1]):  
                      out[i0,i1,i0,i1]=1.    
           elif self.getRank()==3:  
               for i0 in range(self.getShape()[0]):  
                 for i1 in range(self.getShape()[1]):  
                   for i2 in range(self.getShape()[2]):  
                      out[i0,i1,i2,i0,i1,i2]=1.    
           elif self.getRank()==4:  
               for i0 in range(self.getShape()[0]):  
                 for i1 in range(self.getShape()[1]):  
                   for i2 in range(self.getShape()[2]):  
                     for i3 in range(self.getShape()[3]):  
                        out[i0,i1,i2,i3,i0,i1,i2,i3]=1.    
           else:  
              raise ValueError,"differential support rank<5 only."  
           return out  
        else:  
           return self._diff(arg)  
   
    def _diff(self,arg):  
        """return derivate of self with respect to arg (!=self).  
           This method is overwritten by a particular symbol"""  
        return 0  
   
    def eval(self,argval):  
        """subsitutes symbol u in self by argval[u] and returns the result. If  
           self is not a key of argval then self is returned."""  
        if argval.has_key(self):  
          return argval[self]  
        else:  
          return self  
   
   
    def __str__(self):  
        """returns a string representation of the symbol"""  
        return self.__name  
   
    def __add__(self,other):  
        """adds other to symbol self. if _testForZero(other) self is returned."""  
        if _testForZero(other):  
           return self  
        else:  
           a=_matchShape([self,other])  
           return Add_Symbol(a[0],a[1])  
   
   
    def __radd__(self,other):  
        """adds other to symbol self. if _testForZero(other) self is returned."""  
        return self+other  
   
    def __neg__(self):  
        """returns -self."""  
        return self*(-1.)  
   
    def __pos__(self):  
        """returns +self."""  
        return self  
   
    def __abs__(self):  
        """returns absolute value"""  
        return Abs_Symbol(self)  
   
    def __sub__(self,other):  
        """subtracts other from symbol self. if _testForZero(other) self is returned."""  
        if _testForZero(other):  
           return self  
        else:  
           return self+(-other)  
   
    def __rsub__(self,other):  
        """subtracts symbol self from other."""  
        return -self+other  
   
    def __div__(self,other):  
        """divides symbol self by other."""  
        if isinstance(other,Symbol):  
           a=_matchShape([self,other])  
           return Div_Symbol(a[0],a[1])  
        else:  
           return self*(1./other)  
   
    def __rdiv__(self,other):  
        """dived other by symbol self. if _testForZero(other) 0 is returned."""  
        if _testForZero(other):  
           return 0  
        else:  
           a=_matchShape([self,other])  
           return Div_Symbol(a[0],a[1])  
   
    def __pow__(self,other):  
        """raises symbol self to the power of other"""  
        a=_matchShape([self,other])  
        return Power_Symbol(a[0],a[1])  
   
    def __rpow__(self,other):  
        """raises other to the symbol self"""  
        a=_matchShape([self,other])  
        return Power_Symbol(a[1],a[0])  
   
    def __mul__(self,other):  
        """multiplies other by symbol self. if _testForZero(other) 0 is returned."""  
        if _testForZero(other):  
           return 0  
        else:  
           a=_matchShape([self,other])  
           return Mult_Symbol(a[0],a[1])  
   
    def __rmul__(self,other):  
        """multiplies other by symbol self. if _testSForZero(other) 0 is returned."""  
        return self*other  
   
    def __getitem__(self,sl):  
           print sl  
   
 def Float_Symbol(Symbol):  
     def __init__(self,name="symbol",shape=(),args=[]):  
         Symbol.__init__(self,dim=0,name="symbol",shape=(),args=[])  
   
 class ScalarSymbol(Symbol):  
    """a scalar symbol"""  
    def __init__(self,dim=3,name="scalar"):  
       """creates a scalar symbol of spatial dimension dim"""  
       if hasattr(dim,"getDim"):  
            d=dim.getDim()  
       else:      
            d=dim  
       Symbol.__init__(self,shape=(),dim=d,name=name)  
   
 class VectorSymbol(Symbol):  
    """a vector symbol"""  
    def __init__(self,dim=3,name="vector"):  
       """creates a vector symbol of spatial dimension dim"""  
       if hasattr(dim,"getDim"):  
            d=dim.getDim()  
       else:      
            d=dim  
       Symbol.__init__(self,shape=(d,),dim=d,name=name)  
   
 class TensorSymbol(Symbol):  
    """a tensor symbol"""  
    def __init__(self,dim=3,name="tensor"):  
       """creates a tensor symbol of spatial dimension dim"""  
       if hasattr(dim,"getDim"):  
            d=dim.getDim()  
       else:      
            d=dim  
       Symbol.__init__(self,shape=(d,d),dim=d,name=name)  
   
 class Tensor3Symbol(Symbol):  
    """a tensor order 3 symbol"""  
    def __init__(self,dim=3,name="tensor3"):  
       """creates a tensor order 3 symbol of spatial dimension dim"""  
       if hasattr(dim,"getDim"):  
            d=dim.getDim()  
       else:      
            d=dim  
       Symbol.__init__(self,shape=(d,d,d),dim=d,name=name)  
   
 class Tensor4Symbol(Symbol):  
    """a tensor order 4 symbol"""  
    def __init__(self,dim=3,name="tensor4"):  
       """creates a tensor order 4 symbol of spatial dimension dim"""      
       if hasattr(dim,"getDim"):  
            d=dim.getDim()  
       else:      
            d=dim  
       Symbol.__init__(self,shape=(d,d,d,d),dim=d,name=name)  
   
   
 class Add_Symbol(Symbol):  
    """symbol representing the sum of two arguments"""  
    def __init__(self,arg0,arg1):  
        a=[arg0,arg1]  
        Symbol.__init__(self,dim=_extractDim(a),shape=_extractShape(a),args=a)  
    def __str__(self):  
       return "(%s+%s)"%(str(self.getArgument(0)),str(self.getArgument(1)))  
    def eval(self,argval):  
        a=self.getEvaluatedArguments(argval)  
        return a[0]+a[1]  
    def _diff(self,arg):  
        a=self.getDifferentiatedArguments(arg)  
        return a[0]+a[1]  
   
 class Mult_Symbol(Symbol):  
    """symbol representing the product of two arguments"""  
    def __init__(self,arg0,arg1):  
        a=[arg0,arg1]  
        Symbol.__init__(self,dim=_extractDim(a),shape=_extractShape(a),args=a)  
    def __str__(self):  
       return "(%s*%s)"%(str(self.getArgument(0)),str(self.getArgument(1)))  
    def eval(self,argval):  
        a=self.getEvaluatedArguments(argval)  
        return a[0]*a[1]  
    def _diff(self,arg):  
        a=self.getDifferentiatedArguments(arg)  
        return self.getArgument(1)*a[0]+self.getArgument(0)*a[1]  
   
 class Div_Symbol(Symbol):  
    """symbol representing the quotient of two arguments"""  
    def __init__(self,arg0,arg1):  
        a=[arg0,arg1]  
        Symbol.__init__(self,dim=_extractDim(a),shape=_extractShape(a),args=a)  
    def __str__(self):  
       return "(%s/%s)"%(str(self.getArgument(0)),str(self.getArgument(1)))  
    def eval(self,argval):  
        a=self.getEvaluatedArguments(argval)  
        return a[0]/a[1]  
    def _diff(self,arg):  
        a=self.getDifferentiatedArguments(arg)  
        return (a[0]*self.getArgument(1)-self.getArgument(0)*a[1])/ \  
                           (self.getArgument(1)*self.getArgument(1))  
   
 class Power_Symbol(Symbol):  
    """symbol representing the power of the first argument to the power of the second argument"""  
    def __init__(self,arg0,arg1):  
        a=[arg0,arg1]  
        Symbol.__init__(self,dim=_extractDim(a),shape=_extractShape(a),args=a)  
    def __str__(self):  
       return "(%s**%s)"%(str(self.getArgument(0)),str(self.getArgument(1)))  
    def eval(self,argval):  
        a=self.getEvaluatedArguments(argval)  
        return a[0]**a[1]  
    def _diff(self,arg):  
        a=self.getDifferentiatedArguments(arg)  
        return self*(a[1]*log(self.getArgument(0))+self.getArgument(1)/self.getArgument(0)*a[0])  
   
 class Abs_Symbol(Symbol):  
    """symbol representing absolute value of its argument"""  
    def __init__(self,arg):  
        Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])  
    def __str__(self):  
       return "abs(%s)"%str(self.getArgument(0))  
    def eval(self,argval):  
        return abs(self.getEvaluatedArguments(argval)[0])  
    def _diff(self,arg):  
        return sign(self.getArgument(0))*self.getDifferentiatedArguments(arg)[0]  
24    
25  #=========================================================  #=========================================================
26  #   some little helpers  #   some little helpers
27  #=========================================================  #=========================================================
28  def _testForZero(arg):  def _testForZero(arg):
29     """returns True is arg is considered of being zero"""     """
30       Returns True is arg is considered to be zero.
31       """
32     if isinstance(arg,int):     if isinstance(arg,int):
33        return not arg>0        return not arg>0
34     elif isinstance(arg,float):     elif isinstance(arg,float):
35        return not arg>0.        return not arg>0.
36     elif isinstance(arg,numarray.NumArray):     elif isinstance(arg,numarray.NumArray):
37        a=abs(arg)        a=abs(arg)
38        while isinstance(a,numarray.NumArray): a=numarray.sometrue(a)        while isinstance(a,numarray.NumArray): a=numarray.sometrue(a)
39        return not a>0        return not a>0
40     else:     else:
41        return False        return False
42    
43  def _extractDim(args):  #=========================================================
44      dim=None  def saveVTK(filename,domain=None,**data):
45      for a in args:      """
46         if hasattr(a,"getDim"):      writes a L{Data} objects into a files using the the VTK XML file format.
47            d=a.getDim()  
48            if dim==None:      Example:
49               dim=d  
50            else:         tmp=Scalar(..)
51               if dim!=d: raise ValueError,"inconsistent spatial dimension of arguments"         v=Vector(..)
52      if dim==None:         saveVTK("solution.xml",temperature=tmp,velovity=v)
53         raise ValueError,"cannot recover spatial dimension"  
54      return dim      tmp and v are written into "solution.dx" where tmp is named "temperature" and v is named "velovity"
55    
56  def _identifyShape(arg):      @param filename: file name of the output file
57     """identifies the shape of arg."""      @type filename: C(str}
58     if hasattr(arg,"getShape"):      @param domain: domain of the L{Data} object. If not specified, the domain of the given L{Data} objects is used.
59         arg_shape=arg.getShape()      @type domain: L{escript.Domain}
60     elif hasattr(arg,"shape"):      @keyword <name>: writes the assigned value to the VTK file using <name> as identifier.
61       s=arg.shape      @type <name>: L{Data} object.
62       if callable(s):      @note: The data objects have to be defined on the same domain. They may not be in the same L{FunctionSpace} but one cannot expect that all L{FunctionSpace} can be mixed. Typically, data on the boundary and data on the interior cannot be mixed.
63         arg_shape=s()      """
64       else:      if domain==None:
65         arg_shape=s         for i in data.keys():
66     else:            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()
67         arg_shape=()      if domain==None:
68     return arg_shape          raise ValueError,"no domain detected."
69        else:
70            domain.saveVTK(filename,data)
71    
 def _extractShape(args):  
     """extracts the common shape of the list of arguments args"""  
     shape=None  
     for a in args:  
        a_shape=_identifyShape(a)  
        if shape==None: shape=a_shape  
        if shape!=a_shape: raise ValueError,"inconsistent shape"  
     if shape==None:  
        raise ValueError,"cannot recover shape"  
     return shape  
   
 def _matchShape(args,shape=None):  
     """returns the list of arguments args as object which have all the specified shape.  
        if shape is not given the shape "largest" shape of args is used."""  
     # identify the list of shapes:  
     arg_shapes=[]  
     for a in args: arg_shapes.append(_identifyShape(a))  
     # get the largest shape (currently the longest shape):  
     if shape==None: shape=max(arg_shapes)  
       
     out=[]  
     for i in range(len(args)):  
        if shape==arg_shapes[i]:  
           out.append(args[i])  
        else:  
           if len(shape)==0: # then len(arg_shapes[i])>0  
             raise ValueError,"cannot adopt shape of %s to %s"%(str(args[i]),str(shape))  
           else:  
             if len(arg_shapes[i])==0:  
                 out.append(outer(args[i],numarray.ones(shape)))          
             else:    
                 raise ValueError,"cannot adopt shape of %s to %s"%(str(args[i]),str(shape))  
     return out    
72  #=========================================================  #=========================================================
73  #   wrapper for various mathematical functions:  def saveDX(filename,domain=None,**data):
74        """
75        writes a L{Data} objects into a files using the the DX file format.
76    
77        Example:
78    
79           tmp=Scalar(..)
80           v=Vector(..)
81           saveDX("solution.dx",temperature=tmp,velovity=v)
82    
83        tmp and v are written into "solution.dx" where tmp is named "temperature" and v is named "velovity".
84    
85        @param filename: file name of the output file
86        @type filename: C(str}
87        @param domain: domain of the L{Data} object. If not specified, the domain of the given L{Data} objects is used.
88        @type domain: L{escript.Domain}
89        @keyword <name>: writes the assigned value to the DX file using <name> as identifier. The identifier can be used to select the data set when data are imported into DX.
90        @type <name>: L{Data} object.
91        @note: The data objects have to be defined on the same domain. They may not be in the same L{FunctionSpace} but one cannot expect that all L{FunctionSpace} can be mixed. Typically, data on the boundary and data on the interior cannot be mixed.
92        """
93        if domain==None:
94           for i in data.keys():
95              if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()
96        if domain==None:
97            raise ValueError,"no domain detected."
98        else:
99            domain.saveDX(filename,data)
100    
101  #=========================================================  #=========================================================
 def diff(arg,dep):  
     """returns the derivative of arg with respect to dep. If arg is not Symbol object  
        0 is returned"""  
     if isinstance(arg,Symbol):  
        return arg.diff(dep)  
     elif hasattr(arg,"shape"):  
           if callable(arg.shape):  
               return numarray.zeros(arg.shape(),numarray.Float)  
           else:  
               return numarray.zeros(arg.shape,numarray.Float)  
     else:  
        return 0  
102    
103  def exp(arg):  def exp(arg):
104      """      """
105      @brief applies the exponential function to arg      Applies the exponential function to arg.
106      @param arg (input): argument  
107        @param arg: argument
108      """      """
109      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
110         return Exp_Symbol(arg)         return symbols.Exp_Symbol(arg)
111      elif hasattr(arg,"exp"):      elif hasattr(arg,"exp"):
112         return arg.exp()         return arg.exp()
113      else:      else:
114         return numarray.exp(arg)         return numarray.exp(arg)
115    
 class Exp_Symbol(Symbol):  
    """symbol representing the power of the first argument to the power of the second argument"""  
    def __init__(self,arg):  
        Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])  
    def __str__(self):  
       return "exp(%s)"%str(self.getArgument(0))  
    def eval(self,argval):  
        return exp(self.getEvaluatedArguments(argval)[0])  
    def _diff(self,arg):  
        return self*self.getDifferentiatedArguments(arg)[0]  
   
116  def sqrt(arg):  def sqrt(arg):
117      """      """
118      @brief applies the squre root function to arg      Applies the squre root function to arg.
119      @param arg (input): argument  
120        @param arg: argument
121      """      """
122      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
123         return Sqrt_Symbol(arg)         return symbols.Sqrt_Symbol(arg)
124      elif hasattr(arg,"sqrt"):      elif hasattr(arg,"sqrt"):
125         return arg.sqrt()         return arg.sqrt()
126      else:      else:
127         return numarray.sqrt(arg)               return numarray.sqrt(arg)      
128    
 class Sqrt_Symbol(Symbol):  
    """symbol representing square root of argument"""  
    def __init__(self,arg):  
        Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])  
    def __str__(self):  
       return "sqrt(%s)"%str(self.getArgument(0))  
    def eval(self,argval):  
        return sqrt(self.getEvaluatedArguments(argval)[0])  
    def _diff(self,arg):  
        return (-0.5)/self*self.getDifferentiatedArguments(arg)[0]  
   
129  def log(arg):  def log(arg):
130      """      """
131      @brief applies the logarithmic function bases exp(1.) to arg      Applies the logarithmic function base 10 to arg.
132      @param arg (input): argument  
133        @param arg: argument
134      """      """
135      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
136         return Log_Symbol(arg)         return symbols.Log_Symbol(arg)
137      elif hasattr(arg,"log"):      elif hasattr(arg,"log"):
138         return arg.log()         return arg.log()
139      else:      else:
140         return numarray.log(arg)         return numarray.log10(arg)
141    
142    def ln(arg):
143        """
144        Applies the natural logarithmic function to arg.
145    
146  class Log_Symbol(Symbol):      @param arg: argument
147     """symbol representing logarithm of the argument"""      """
148     def __init__(self,arg):      if isinstance(arg,symbols.Symbol):
149         Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])         return symbols.Ln_Symbol(arg)
150     def __str__(self):      elif hasattr(arg,"ln"):
151        return "log(%s)"%str(self.getArgument(0))         return arg.ln()
152     def eval(self,argval):      else:
153         return log(self.getEvaluatedArguments(argval)[0])         return numarray.log(arg)
    def _diff(self,arg):  
        return self.getDifferentiatedArguments(arg)[0]/self.getArgument(0)  
154    
155  def sin(arg):  def sin(arg):
156      """      """
157      @brief applies the sinus function to arg      Applies the sin function to arg.
158      @param arg (input): argument  
159        @param arg: argument
160      """      """
161      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
162         return Sin_Symbol(arg)         return symbols.Sin_Symbol(arg)
163      elif hasattr(arg,"sin"):      elif hasattr(arg,"sin"):
164         return arg.sin()         return arg.sin()
165      else:      else:
166         return numarray.sin(arg)         return numarray.sin(arg)
167    
 class Sin_Symbol(Symbol):  
    """symbol representing logarithm of the argument"""  
    def __init__(self,arg):  
        Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])  
    def __str__(self):  
       return "sin(%s)"%str(self.getArgument(0))  
    def eval(self,argval):  
        return sin(self.getEvaluatedArguments(argval)[0])  
    def _diff(self,arg):  
        return cos(self.getArgument(0))*self.getDifferentiatedArguments(arg)[0]  
   
168  def cos(arg):  def cos(arg):
169      """      """
170      @brief applies the sinus function to arg      Applies the cos function to arg.
171      @param arg (input): argument  
172        @param arg: argument
173      """      """
174      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
175         return Cos_Symbol(arg)         return symbols.Cos_Symbol(arg)
176      elif hasattr(arg,"cos"):      elif hasattr(arg,"cos"):
177         return arg.cos()         return arg.cos()
178      else:      else:
179         return numarray.cos(arg)         return numarray.cos(arg)
180    
 class Cos_Symbol(Symbol):  
    """symbol representing logarithm of the argument"""  
    def __init__(self,arg):  
        Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])  
    def __str__(self):  
       return "cos(%s)"%str(self.getArgument(0))  
    def eval(self,argval):  
        return cos(self.getEvaluatedArguments(argval)[0])  
    def _diff(self,arg):  
        return -sin(self.getArgument(0))*self.getDifferentiatedArguments(arg)[0]  
   
181  def tan(arg):  def tan(arg):
182      """      """
183      @brief applies the sinus function to arg      Applies the tan function to arg.
184      @param arg (input): argument  
185        @param arg: argument
186      """      """
187      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
188         return Tan_Symbol(arg)         return symbols.Tan_Symbol(arg)
189      elif hasattr(arg,"tan"):      elif hasattr(arg,"tan"):
190         return arg.tan()         return arg.tan()
191      else:      else:
192         return numarray.tan(arg)         return numarray.tan(arg)
193    
194  class Tan_Symbol(Symbol):  def asin(arg):
195     """symbol representing logarithm of the argument"""      """
196     def __init__(self,arg):      Applies the asin function to arg.
197         Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])  
198     def __str__(self):      @param arg: argument
199        return "tan(%s)"%str(self.getArgument(0))      """
200     def eval(self,argval):      if isinstance(arg,symbols.Symbol):
201         return tan(self.getEvaluatedArguments(argval)[0])         return symbols.Asin_Symbol(arg)
202     def _diff(self,arg):      elif hasattr(arg,"asin"):
203         s=cos(self.getArgument(0))         return arg.asin()
204         return 1./(s*s)*self.getDifferentiatedArguments(arg)[0]      else:
205           return numarray.asin(arg)
206    
207    def acos(arg):
208        """
209        Applies the acos function to arg.
210    
211        @param arg: argument
212        """
213        if isinstance(arg,symbols.Symbol):
214           return symbols.Acos_Symbol(arg)
215        elif hasattr(arg,"acos"):
216           return arg.acos()
217        else:
218           return numarray.acos(arg)
219    
220    def atan(arg):
221        """
222        Applies the atan function to arg.
223    
224        @param arg: argument
225        """
226        if isinstance(arg,symbols.Symbol):
227           return symbols.Atan_Symbol(arg)
228        elif hasattr(arg,"atan"):
229           return arg.atan()
230        else:
231           return numarray.atan(arg)
232    
233    def sinh(arg):
234        """
235        Applies the sinh function to arg.
236    
237        @param arg: argument
238        """
239        if isinstance(arg,symbols.Symbol):
240           return symbols.Sinh_Symbol(arg)
241        elif hasattr(arg,"sinh"):
242           return arg.sinh()
243        else:
244           return numarray.sinh(arg)
245    
246    def cosh(arg):
247        """
248        Applies the cosh function to arg.
249    
250        @param arg: argument
251        """
252        if isinstance(arg,symbols.Symbol):
253           return symbols.Cosh_Symbol(arg)
254        elif hasattr(arg,"cosh"):
255           return arg.cosh()
256        else:
257           return numarray.cosh(arg)
258    
259    def tanh(arg):
260        """
261        Applies the tanh function to arg.
262    
263        @param arg: argument
264        """
265        if isinstance(arg,symbols.Symbol):
266           return symbols.Tanh_Symbol(arg)
267        elif hasattr(arg,"tanh"):
268           return arg.tanh()
269        else:
270           return numarray.tanh(arg)
271    
272    def asinh(arg):
273        """
274        Applies the asinh function to arg.
275    
276        @param arg: argument
277        """
278        if isinstance(arg,symbols.Symbol):
279           return symbols.Asinh_Symbol(arg)
280        elif hasattr(arg,"asinh"):
281           return arg.asinh()
282        else:
283           return numarray.asinh(arg)
284    
285    def acosh(arg):
286        """
287        Applies the acosh function to arg.
288    
289        @param arg: argument
290        """
291        if isinstance(arg,symbols.Symbol):
292           return symbols.Acosh_Symbol(arg)
293        elif hasattr(arg,"acosh"):
294           return arg.acosh()
295        else:
296           return numarray.acosh(arg)
297    
298    def atanh(arg):
299        """
300        Applies the atanh function to arg.
301    
302        @param arg: argument
303        """
304        if isinstance(arg,symbols.Symbol):
305           return symbols.Atanh_Symbol(arg)
306        elif hasattr(arg,"atanh"):
307           return arg.atanh()
308        else:
309           return numarray.atanh(arg)
310    
311  def sign(arg):  def sign(arg):
312      """      """
313      @brief applies the sign function to arg      Applies the sign function to arg.
314      @param arg (input): argument  
315        @param arg: argument
316      """      """
317      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
318         return Sign_Symbol(arg)         return symbols.Sign_Symbol(arg)
319      elif hasattr(arg,"sign"):      elif hasattr(arg,"sign"):
320         return arg.sign()         return arg.sign()
321      else:      else:
322         return numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))- \         return numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))- \
323                numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))                numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))
324    
 class Sign_Symbol(Symbol):  
    """symbol representing the sign of the argument"""  
    def __init__(self,arg):  
        Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])  
    def __str__(self):  
       return "sign(%s)"%str(self.getArgument(0))  
    def eval(self,argval):  
        return sign(self.getEvaluatedArguments(argval)[0])  
   
325  def maxval(arg):  def maxval(arg):
326      """      """
327      @brief returns the maximum value of argument arg""      Returns the maximum value of argument arg.
328      @param arg (input): argument  
329        @param arg: argument
330      """      """
331      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
332         return Max_Symbol(arg)         return symbols.Max_Symbol(arg)
333      elif hasattr(arg,"maxval"):      elif hasattr(arg,"maxval"):
334         return arg.maxval()         return arg.maxval()
335      elif hasattr(arg,"max"):      elif hasattr(arg,"max"):
# Line 637  def maxval(arg): Line 337  def maxval(arg):
337      else:      else:
338         return arg         return arg
339    
 class Max_Symbol(Symbol):  
    """symbol representing the sign of the argument"""  
    def __init__(self,arg):  
        Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])  
    def __str__(self):  
       return "maxval(%s)"%str(self.getArgument(0))  
    def eval(self,argval):  
        return maxval(self.getEvaluatedArguments(argval)[0])  
   
340  def minval(arg):  def minval(arg):
341      """      """
342      @brief returns the maximum value of argument arg""      Returns the minimum value of argument arg.
343      @param arg (input): argument  
344        @param arg: argument
345      """      """
346      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
347         return Min_Symbol(arg)         return symbols.Min_Symbol(arg)
348      elif hasattr(arg,"maxval"):      elif hasattr(arg,"maxval"):
349         return arg.minval()         return arg.minval()
350      elif hasattr(arg,"min"):      elif hasattr(arg,"min"):
# Line 660  def minval(arg): Line 352  def minval(arg):
352      else:      else:
353         return arg         return arg
354    
 class Min_Symbol(Symbol):  
    """symbol representing the sign of the argument"""  
    def __init__(self,arg):  
        Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])  
    def __str__(self):  
       return "minval(%s)"%str(self.getArgument(0))  
    def eval(self,argval):  
        return minval(self.getEvaluatedArguments(argval)[0])  
   
355  def wherePositive(arg):  def wherePositive(arg):
356      """      """
357      @brief returns the maximum value of argument arg""      Returns the positive values of argument arg.
358      @param arg (input): argument  
359        @param arg: argument
360      """      """
361      if _testForZero(arg):      if _testForZero(arg):
362        return 0        return 0
363      elif isinstance(arg,Symbol):      elif isinstance(arg,symbols.Symbol):
364         return WherePositive_Symbol(arg)         return symbols.WherePositive_Symbol(arg)
365      elif hasattr(arg,"wherePositive"):      elif hasattr(arg,"wherePositive"):
366         return arg.minval()         return arg.minval()
367      elif hasattr(arg,"wherePositive"):      elif hasattr(arg,"wherePositive"):
# Line 688  def wherePositive(arg): Line 372  def wherePositive(arg):
372         else:         else:
373            return 0.            return 0.
374    
 class WherePositive_Symbol(Symbol):  
    """symbol representing the wherePositive function"""  
    def __init__(self,arg):  
        Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])  
    def __str__(self):  
       return "wherePositive(%s)"%str(self.getArgument(0))  
    def eval(self,argval):  
        return wherePositive(self.getEvaluatedArguments(argval)[0])  
   
375  def whereNegative(arg):  def whereNegative(arg):
376      """      """
377      @brief returns the maximum value of argument arg""      Returns the negative values of argument arg.
378      @param arg (input): argument  
379        @param arg: argument
380      """      """
381      if _testForZero(arg):      if _testForZero(arg):
382        return 0        return 0
383      elif isinstance(arg,Symbol):      elif isinstance(arg,symbols.Symbol):
384         return WhereNegative_Symbol(arg)         return symbols.WhereNegative_Symbol(arg)
385      elif hasattr(arg,"whereNegative"):      elif hasattr(arg,"whereNegative"):
386         return arg.whereNegative()         return arg.whereNegative()
387      elif hasattr(arg,"shape"):      elif hasattr(arg,"shape"):
# Line 716  def whereNegative(arg): Line 392  def whereNegative(arg):
392         else:         else:
393            return 0.            return 0.
394    
395  class WhereNegative_Symbol(Symbol):  def maximum(arg0,arg1):
396     """symbol representing the whereNegative function"""      """
397     def __init__(self,arg):      Return arg1 where arg1 is bigger then arg0 otherwise arg0 is returned.
398         Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])      """
399     def __str__(self):      m=whereNegative(arg0-arg1)
400        return "whereNegative(%s)"%str(self.getArgument(0))      return m*arg1+(1.-m)*arg0
401     def eval(self,argval):    
402         return whereNegative(self.getEvaluatedArguments(argval)[0])  def minimum(arg0,arg1):
403        """
404        Return arg0 where arg1 is bigger then arg0 otherwise arg1 is returned.
405        """
406        m=whereNegative(arg0-arg1)
407        return m*arg0+(1.-m)*arg1
408      
409  def outer(arg0,arg1):  def outer(arg0,arg1):
410     if _testForZero(arg0) or _testForZero(arg1):     if _testForZero(arg0) or _testForZero(arg1):
411        return 0        return 0
412     else:     else:
413        if isinstance(arg0,Symbol) or isinstance(arg1,Symbol):        if isinstance(arg0,symbols.Symbol) or isinstance(arg1,symbols.Symbol):
414          return Outer_Symbol(arg0,arg1)          return symbols.Outer_Symbol(arg0,arg1)
415        elif _identifyShape(arg0)==() or _identifyShape(arg1)==():        elif _identifyShape(arg0)==() or _identifyShape(arg1)==():
416          return arg0*arg1          return arg0*arg1
417        elif isinstance(arg0,numarray.NumArray) and isinstance(arg1,numarray.NumArray):        elif isinstance(arg0,numarray.NumArray) and isinstance(arg1,numarray.NumArray):
# Line 745  def outer(arg0,arg1): Line 426  def outer(arg0,arg1):
426          else:          else:
427            raise ValueError,"outer is not fully implemented yet."            raise ValueError,"outer is not fully implemented yet."
428    
 class Outer_Symbol(Symbol):  
    """symbol representing the outer product of its two argument"""  
    def __init__(self,arg0,arg1):  
        a=[arg0,arg1]  
        s=tuple(list(_identifyShape(arg0))+list(_identifyShape(arg1)))  
        Symbol.__init__(self,shape=s,dim=_extractDim(a),args=a)  
    def __str__(self):  
       return "outer(%s,%s)"%(str(self.getArgument(0)),str(self.getArgument(1)))  
    def eval(self,argval):  
        a=self.getEvaluatedArguments(argval)  
        return outer(a[0],a[1])  
    def _diff(self,arg):  
        a=self.getDifferentiatedArguments(arg)  
        return outer(a[0],self.getArgument(1))+outer(self.getArgument(0),a[1])  
   
429  def interpolate(arg,where):  def interpolate(arg,where):
430      """      """
431      @brief interpolates the function into the FunctionSpace where.      Interpolates the function into the FunctionSpace where.
432    
433      @param arg    interpolant      @param arg:    interpolant
434      @param where  FunctionSpace to interpolate to      @param where:  FunctionSpace to interpolate to
435      """      """
436      if _testForZero(arg):      if _testForZero(arg):
437        return 0        return 0
438      elif isinstance(arg,Symbol):      elif isinstance(arg,symbols.Symbol):
439         return Interpolated_Symbol(arg,where)         return symbols.Interpolated_Symbol(arg,where)
440      else:      else:
441         return escript.Data(arg,where)         return escript.Data(arg,where)
442    
443  def Interpolated_Symbol(Symbol):  def div(arg,where=None):
444     """symbol representing the integral of the argument"""      """
445     def __init__(self,arg,where):      Returns the divergence of arg at where.
446          Symbol.__init__(self,shape=_extractShape(arg),dim=_extractDim([arg]),args=[arg,where])  
447     def __str__(self):      @param arg:   Data object representing the function which gradient to
448        return "interpolated(%s)"%(str(self.getArgument(0)))                    be calculated.
449     def eval(self,argval):      @param where: FunctionSpace in which the gradient will be calculated.
450         a=self.getEvaluatedArguments(argval)                    If not present or C{None} an appropriate default is used.
451         return integrate(a[0],where=self.getArgument(1))      """
452     def _diff(self,arg):      return trace(grad(arg,where))
453         a=self.getDifferentiatedArguments(arg)  
454         return integrate(a[0],where=self.getArgument(1))  def jump(arg):
455        """
456        Returns the jump of arg across a continuity.
457    
458        @param arg:   Data object representing the function which gradient
459                      to be calculated.
460        """
461        d=arg.getDomain()
462        return arg.interpolate(escript.FunctionOnContactOne())-arg.interpolate(escript.FunctionOnContactZero())
463      
464    
465  def grad(arg,where=None):  def grad(arg,where=None):
466      """      """
467      @brief returns the spatial gradient of arg at where.      Returns the spatial gradient of arg at where.
468    
469      @param arg:   Data object representing the function which gradient to be calculated.      @param arg:   Data object representing the function which gradient
470      @param where: FunctionSpace in which the gradient will be calculated. If not present or                    to be calculated.
471                    None an appropriate default is used.      @param where: FunctionSpace in which the gradient will be calculated.
472                      If not present or C{None} an appropriate default is used.
473      """      """
474      if _testForZero(arg):      if _testForZero(arg):
475        return 0        return 0
476      elif isinstance(arg,Symbol):      elif isinstance(arg,symbols.Symbol):
477         return Grad_Symbol(arg,where)         return symbols.Grad_Symbol(arg,where)
478      elif hasattr(arg,"grad"):      elif hasattr(arg,"grad"):
479         if where==None:         if where==None:
480            return arg.grad()            return arg.grad()
# Line 807  def grad(arg,where=None): Line 483  def grad(arg,where=None):
483      else:      else:
484         return arg*0.         return arg*0.
485    
 def Grad_Symbol(Symbol):  
    """symbol representing the gradient of the argument"""  
    def __init__(self,arg,where=None):  
        d=_extractDim([arg])  
        s=tuple(list(_identifyShape([arg])).append(d))  
        Symbol.__init__(self,shape=s,dim=_extractDim([arg]),args=[arg,where])  
    def __str__(self):  
       return "grad(%s)"%(str(self.getArgument(0)))  
    def eval(self,argval):  
        a=self.getEvaluatedArguments(argval)  
        return grad(a[0],where=self.getArgument(1))  
    def _diff(self,arg):  
        a=self.getDifferentiatedArguments(arg)  
        return grad(a[0],where=self.getArgument(1))  
   
486  def integrate(arg,where=None):  def integrate(arg,where=None):
487      """      """
488      @brief return the integral if the function represented by Data object arg over its domain.      Return the integral if the function represented by Data object arg over
489        its domain.
490    
491      @param arg:   Data object representing the function which is integrated.      @param arg:   Data object representing the function which is integrated.
492      @param where: FunctionSpace in which the integral is calculated. If not present or      @param where: FunctionSpace in which the integral is calculated.
493                    None an appropriate default is used.                    If not present or C{None} an appropriate default is used.
494      """      """
495      if _testForZero(arg):      if _testForZero(arg):
496        return 0        return 0
497      elif isinstance(arg,Symbol):      elif isinstance(arg,symbols.Symbol):
498         return Integral_Symbol(arg,where)         return symbols.Integral_Symbol(arg,where)
499      else:          else:    
500         if not where==None: arg=escript.Data(arg,where)         if not where==None: arg=escript.Data(arg,where)
501         if arg.getRank()==0:         if arg.getRank()==0:
# Line 841  def integrate(arg,where=None): Line 503  def integrate(arg,where=None):
503         else:         else:
504           return arg.integrate()           return arg.integrate()
505    
 def Integral_Symbol(Float_Symbol):  
    """symbol representing the integral of the argument"""  
    def __init__(self,arg,where=None):  
        Float_Symbol.__init__(self,shape=_identifyShape([arg]),args=[arg,where])  
    def __str__(self):  
       return "integral(%s)"%(str(self.getArgument(0)))  
    def eval(self,argval):  
        a=self.getEvaluatedArguments(argval)  
        return integrate(a[0],where=self.getArgument(1))  
    def _diff(self,arg):  
        a=self.getDifferentiatedArguments(arg)  
        return integrate(a[0],where=self.getArgument(1))  
   
506  #=============================  #=============================
507  #  #
508  # wrapper for various functions: if the argument has attribute the function name  # wrapper for various functions: if the argument has attribute the function name
509  # as an argument it calls the correspong methods. Otherwise the coresponsing numarray  # as an argument it calls the corresponding methods. Otherwise the corresponding
510  # function is called.  # numarray function is called.
511  #  
512  # functions involving the underlying Domain:  # functions involving the underlying Domain:
 #  
513    
514    
515  # functions returning Data objects:  # functions returning Data objects:
516    
517  def transpose(arg,axis=None):  def transpose(arg,axis=None):
518      """      """
519      @brief returns the transpose of the Data object arg.      Returns the transpose of the Data object arg.
520    
521      @param arg      @param arg:
522      """      """
523      if axis==None:      if axis==None:
524         r=0         r=0
525         if hasattr(arg,"getRank"): r=arg.getRank()         if hasattr(arg,"getRank"): r=arg.getRank()
526         if hasattr(arg,"rank"): r=arg.rank         if hasattr(arg,"rank"): r=arg.rank
527         axis=r/2         axis=r/2
528      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
529         return Transpose_Symbol(arg,axis=r)         return symbols.Transpose_Symbol(arg,axis=r)
530      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
531         # hack for transpose         # hack for transpose
532         r=arg.getRank()         r=arg.getRank()
# Line 896  def transpose(arg,axis=None): Line 544  def transpose(arg,axis=None):
544    
545  def trace(arg,axis0=0,axis1=1):  def trace(arg,axis0=0,axis1=1):
546      """      """
547      @brief return      Return
548    
549      @param arg      @param arg:
550      """      """
551      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
552         s=list(arg.getShape())                 s=list(arg.getShape())        
553         s=tuple(s[0:axis0]+s[axis0+1:axis1]+s[axis1+1:])         s=tuple(s[0:axis0]+s[axis0+1:axis1]+s[axis1+1:])
554         return Trace_Symbol(arg,axis0=axis0,axis1=axis1)         return symbols.Trace_Symbol(arg,axis0=axis0,axis1=axis1)
555      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
556         # hack for trace         # hack for trace
557         s=arg.getShape()         s=arg.getShape()
558         if s[axis0]!=s[axis1]:         if s[axis0]!=s[axis1]:
559             raise ValueError,"illegal axis in trace"             raise ValueError,"illegal axis in trace"
560         out=escript.Scalar(0.,arg.getFunctionSpace())         out=escript.Scalar(0.,arg.getFunctionSpace())
561         for i in range(s[0]):         for i in range(s[axis0]):
562            for j in range(s[1]):            out+=arg[i,i]
              out+=arg[i,j]  
563         return out         return out
564         # end hack for transpose         # end hack for trace
        return arg.transpose(axis0=axis0,axis1=axis1)  
565      else:      else:
566         return numarray.trace(arg,axis0=axis0,axis1=axis1)         return numarray.trace(arg,axis0=axis0,axis1=axis1)
567    
   
   
 def Trace_Symbol(Symbol):  
     pass  
   
568  def length(arg):  def length(arg):
569      """      """
570      @brief      @param arg:
   
     @param arg  
571      """      """
572      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
573         if arg.isEmpty(): return escript.Data()         if arg.isEmpty(): return escript.Data()
574         if arg.getRank()==0:         if arg.getRank()==0:
575            return abs(arg)            return abs(arg)
576         elif arg.getRank()==1:         elif arg.getRank()==1:
577            sum=escript.Scalar(0,arg.getFunctionSpace())            out=escript.Scalar(0,arg.getFunctionSpace())
578            for i in range(arg.getShape()[0]):            for i in range(arg.getShape()[0]):
579               sum+=arg[i]**2               out+=arg[i]**2
580            return sqrt(sum)            return sqrt(out)
581         elif arg.getRank()==2:         elif arg.getRank()==2:
582            sum=escript.Scalar(0,arg.getFunctionSpace())            out=escript.Scalar(0,arg.getFunctionSpace())
583            for i in range(arg.getShape()[0]):            for i in range(arg.getShape()[0]):
584               for j in range(arg.getShape()[1]):               for j in range(arg.getShape()[1]):
585                  sum+=arg[i,j]**2                  out+=arg[i,j]**2
586            return sqrt(sum)            return sqrt(out)
587         elif arg.getRank()==3:         elif arg.getRank()==3:
588            sum=escript.Scalar(0,arg.getFunctionSpace())            out=escript.Scalar(0,arg.getFunctionSpace())
589            for i in range(arg.getShape()[0]):            for i in range(arg.getShape()[0]):
590               for j in range(arg.getShape()[1]):               for j in range(arg.getShape()[1]):
591                  for k in range(arg.getShape()[2]):                  for k in range(arg.getShape()[2]):
592                     sum+=arg[i,j,k]**2                     out+=arg[i,j,k]**2
593            return sqrt(sum)            return sqrt(out)
594         elif arg.getRank()==4:         elif arg.getRank()==4:
595            sum=escript.Scalar(0,arg.getFunctionSpace())            out=escript.Scalar(0,arg.getFunctionSpace())
596            for i in range(arg.getShape()[0]):            for i in range(arg.getShape()[0]):
597               for j in range(arg.getShape()[1]):               for j in range(arg.getShape()[1]):
598                  for k in range(arg.getShape()[2]):                  for k in range(arg.getShape()[2]):
599                     for l in range(arg.getShape()[3]):                     for l in range(arg.getShape()[3]):
600                        sum+=arg[i,j,k,l]**2                        out+=arg[i,j,k,l]**2
601            return sqrt(sum)            return sqrt(out)
602         else:         else:
603            raise SystemError,"length is not been implemented yet"            raise SystemError,"length is not been fully implemented yet"
604         # return arg.length()            # return arg.length()
605        elif isinstance(arg,float):
606           return abs(arg)
607      else:      else:
608         return sqrt((arg**2).sum())         return sqrt((arg**2).sum())
609    
610  def deviator(arg):  def deviator(arg):
611      """      """
612      @brief      @param arg:
   
     @param arg0  
613      """      """
614      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
615          shape=arg.getShape()          shape=arg.getShape()
# Line 984  def deviator(arg): Line 623  def deviator(arg):
623    
624  def inner(arg0,arg1):  def inner(arg0,arg1):
625      """      """
626      @brief      @param arg0:
627        @param arg1:
     @param arg0, arg1  
628      """      """
629      sum=escript.Scalar(0,arg0.getFunctionSpace())      if isinstance(arg0,escript.Data):
630           arg=arg0
631        else:
632           arg=arg1
633    
634        out=escript.Scalar(0,arg.getFunctionSpace())
635      if arg.getRank()==0:      if arg.getRank()==0:
636            return arg0*arg1            return arg0*arg1
637      elif arg.getRank()==1:      elif arg.getRank()==1:
638           sum=escript.Scalar(0,arg.getFunctionSpace())           out=escript.Scalar(0,arg.getFunctionSpace())
639           for i in range(arg.getShape()[0]):           for i in range(arg.getShape()[0]):
640              sum+=arg0[i]*arg1[i]              out+=arg0[i]*arg1[i]
641      elif arg.getRank()==2:      elif arg.getRank()==2:
642          sum=escript.Scalar(0,arg.getFunctionSpace())          out=escript.Scalar(0,arg.getFunctionSpace())
643          for i in range(arg.getShape()[0]):          for i in range(arg.getShape()[0]):
644             for j in range(arg.getShape()[1]):             for j in range(arg.getShape()[1]):
645                sum+=arg0[i,j]*arg1[i,j]                out+=arg0[i,j]*arg1[i,j]
646      elif arg.getRank()==3:      elif arg.getRank()==3:
647          sum=escript.Scalar(0,arg.getFunctionSpace())          out=escript.Scalar(0,arg.getFunctionSpace())
648          for i in range(arg.getShape()[0]):          for i in range(arg.getShape()[0]):
649              for j in range(arg.getShape()[1]):              for j in range(arg.getShape()[1]):
650                 for k in range(arg.getShape()[2]):                 for k in range(arg.getShape()[2]):
651                    sum+=arg0[i,j,k]*arg1[i,j,k]                    out+=arg0[i,j,k]*arg1[i,j,k]
652      elif arg.getRank()==4:      elif arg.getRank()==4:
653          sum=escript.Scalar(0,arg.getFunctionSpace())          out=escript.Scalar(0,arg.getFunctionSpace())
654          for i in range(arg.getShape()[0]):          for i in range(arg.getShape()[0]):
655             for j in range(arg.getShape()[1]):             for j in range(arg.getShape()[1]):
656                for k in range(arg.getShape()[2]):                for k in range(arg.getShape()[2]):
657                   for l in range(arg.getShape()[3]):                   for l in range(arg.getShape()[3]):
658                      sum+=arg0[i,j,k,l]*arg1[i,j,k,l]                      out+=arg0[i,j,k,l]*arg1[i,j,k,l]
659      else:      else:
660            raise SystemError,"inner is not been implemented yet"            raise SystemError,"inner is not been implemented yet"
661      return sum      return out
662    
663    def tensormult(arg0,arg1):
664        # check LinearPDE!!!!
665        raise SystemError,"tensormult is not implemented yet!"
666    
667  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
668    
# Line 1031  def matrixmult(arg0,arg1): Line 678  def matrixmult(arg0,arg1):
678            out=escript.Data(0,(arg0.getShape()[0],),arg0.getFunctionSpace())            out=escript.Data(0,(arg0.getShape()[0],),arg0.getFunctionSpace())
679            for i in range(arg0.getShape()[0]):            for i in range(arg0.getShape()[0]):
680               for j in range(arg0.getShape()[1]):               for j in range(arg0.getShape()[1]):
681                   # uses Data object slicing, plus Data * and += operators
682                 out[i]+=arg0[i,j]*arg1[j]                 out[i]+=arg0[i,j]*arg1[j]
683            return out            return out
684          elif arg0.getRank()==1 and arg1.getRank()==1:
685              return inner(arg0,arg1)
686        else:        else:
687            raise SystemError,"matrixmult is not fully implemented yet!"            raise SystemError,"matrixmult is not fully implemented yet!"
688    
689  #=========================================================  #=========================================================
690  # reduction operations:  # reduction operations:
691  #=========================================================  #=========================================================
692  def sum(arg):  def sum(arg):
693      """      """
694      @brief      @param arg:
   
     @param arg  
695      """      """
696      return arg.sum()      return arg.sum()
697    
698  def sup(arg):  def sup(arg):
699      """      """
700      @brief      @param arg:
   
     @param arg  
701      """      """
702      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
703         return arg.sup()         return arg.sup()
# Line 1061  def sup(arg): Line 708  def sup(arg):
708    
709  def inf(arg):  def inf(arg):
710      """      """
711      @brief      @param arg:
   
     @param arg  
712      """      """
713      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
714         return arg.inf()         return arg.inf()
# Line 1074  def inf(arg): Line 719  def inf(arg):
719    
720  def L2(arg):  def L2(arg):
721      """      """
722      @brief returns the L2-norm of the      Returns the L2-norm of the argument
723    
724      @param arg      @param arg:
725      """      """
726      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
727         return arg.L2()         return arg.L2()
# Line 1087  def L2(arg): Line 732  def L2(arg):
732    
733  def Lsup(arg):  def Lsup(arg):
734      """      """
735      @brief      @param arg:
   
     @param arg  
736      """      """
737      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
738         return arg.Lsup()         return arg.Lsup()
739      elif isinstance(arg,float) or isinstance(arg,int):      elif isinstance(arg,float) or isinstance(arg,int):
740         return abs(arg)         return abs(arg)
741      else:      else:
742         return max(numarray.abs(arg))         return numarray.abs(arg).max()
743    
744  def dot(arg0,arg1):  def dot(arg0,arg1):
745      """      """
746      @brief      @param arg0:
747        @param arg1:
     @param arg  
748      """      """
749      if isinstance(arg0,escript.Data):      if isinstance(arg0,escript.Data):
750         return arg0.dot(arg1)         return arg0.dot(arg1)
# Line 1113  def dot(arg0,arg1): Line 755  def dot(arg0,arg1):
755    
756  def kronecker(d):  def kronecker(d):
757     if hasattr(d,"getDim"):     if hasattr(d,"getDim"):
758        return numarray.identity(d.getDim())        return numarray.identity(d.getDim())*1.
759     else:     else:
760        return numarray.identity(d)        return numarray.identity(d)*1.
761    
762  def unit(i,d):  def unit(i,d):
763     """     """
764     @brief return a unit vector of dimension d with nonzero index i     Return a unit vector of dimension d with nonzero index i.
765     @param d dimension  
766     @param i index     @param d: dimension
767       @param i: index
768     """     """
769     e = numarray.zeros((d,),numarray.Float)     e = numarray.zeros((d,),numarray.Float)
770     e[i] = 1.0     e[i] = 1.0
771     return e     return e
   
 #  
 # ============================================  
 #   testing  
 # ============================================  
   
 if __name__=="__main__":  
   u=ScalarSymbol(dim=2,name="u")  
   v=ScalarSymbol(dim=2,name="v")  
   v=VectorSymbol(2,"v")  
   u=VectorSymbol(2,"u")  
   
   
   print u+5,(u+5).diff(u)  
   print 5+u,(5+u).diff(u)  
   print u+v,(u+v).diff(u)  
   print v+u,(v+u).diff(u)  
   
   print u*5,(u*5).diff(u)  
   print 5*u,(5*u).diff(u)  
   print u*v,(u*v).diff(u)  
   print v*u,(v*u).diff(u)  
   
   print u-5,(u-5).diff(u)  
   print 5-u,(5-u).diff(u)  
   print u-v,(u-v).diff(u)  
   print v-u,(v-u).diff(u)  
   
   print u/5,(u/5).diff(u)  
   print 5/u,(5/u).diff(u)  
   print u/v,(u/v).diff(u)  
   print v/u,(v/u).diff(u)  
   
   print u**5,(u**5).diff(u)  
   print 5**u,(5**u).diff(u)  
   print u**v,(u**v).diff(u)  
   print v**u,(v**u).diff(u)  
   
   print exp(u),exp(u).diff(u)  
   print sqrt(u),sqrt(u).diff(u)  
   print log(u),log(u).diff(u)  
   print sin(u),sin(u).diff(u)  
   print cos(u),cos(u).diff(u)  
   print tan(u),tan(u).diff(u)  
   print sign(u),sign(u).diff(u)  
   print abs(u),abs(u).diff(u)  
   print wherePositive(u),wherePositive(u).diff(u)  
   print whereNegative(u),whereNegative(u).diff(u)  
   print maxval(u),maxval(u).diff(u)  
   print minval(u),minval(u).diff(u)  
   
   g=grad(u)  
   print diff(5*g,g)  
   4*(g+transpose(g))/2+6*trace(g)*kronecker(3)  
 #  
 # $Log$  
 # Revision 1.11  2005/07/08 04:07:35  jgs  
 # Merge of development branch back to main trunk on 2005-07-08  
 #  
 # Revision 1.10  2005/06/09 05:37:59  jgs  
 # Merge of development branch back to main trunk on 2005-06-09  
 #  
 # Revision 1.2.2.17  2005/07/07 07:28:58  gross  
 # some stuff added to util.py to improve functionality  
 #  
 # Revision 1.2.2.16  2005/06/30 01:53:55  gross  
 # a bug in coloring fixed  
 #  
 # Revision 1.2.2.15  2005/06/29 02:36:43  gross  
 # Symbols have been introduced and some function clarified. needs much more work  
 #  
 # Revision 1.2.2.14  2005/05/20 04:05:23  gross  
 # some work on a darcy flow started  
 #  
 # Revision 1.2.2.13  2005/03/16 05:17:58  matt  
 # Implemented unit(idx, dim) to create cartesian unit basis vectors to  
 # complement kronecker(dim) function.  
 #  
 # Revision 1.2.2.12  2005/03/10 08:14:37  matt  
 # Added non-member Linf utility function to complement Data::Linf().  
 #  
 # Revision 1.2.2.11  2005/02/17 05:53:25  gross  
 # some bug in saveDX fixed: in fact the bug was in  
 # DataC/getDataPointShape  
 #  
 # Revision 1.2.2.10  2005/01/11 04:59:36  gross  
 # automatic interpolation in integrate switched off  
 #  
 # Revision 1.2.2.9  2005/01/11 03:38:13  gross  
 # Bug in Data.integrate() fixed for the case of rank 0. The problem is not totallly resolved as the method should return a scalar rather than a numarray object in the case of rank 0. This problem is fixed by the util.integrate wrapper.  
 #  
 # Revision 1.2.2.8  2005/01/05 04:21:41  gross  
 # FunctionSpace checking/matchig in slicing added  
 #  
 # Revision 1.2.2.7  2004/12/29 05:29:59  gross  
 # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()  
 #  
 # Revision 1.2.2.6  2004/12/24 06:05:41  gross  
 # some changes in linearPDEs to add AdevectivePDE  
 #  
 # Revision 1.2.2.5  2004/12/17 00:06:53  gross  
 # mk sets ESYS_ROOT is undefined  
 #  
 # Revision 1.2.2.4  2004/12/07 03:19:51  gross  
 # options for GMRES and PRES20 added  
 #  
 # Revision 1.2.2.3  2004/12/06 04:55:18  gross  
 # function wraper extended  
 #  
 # Revision 1.2.2.2  2004/11/22 05:44:07  gross  
 # a few more unitary functions have been added but not implemented in Data yet  
 #  
 # Revision 1.2.2.1  2004/11/12 06:58:15  gross  
 # a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry  
 #  
 # Revision 1.2  2004/10/27 00:23:36  jgs  
 # fixed minor syntax error  
 #  
 # Revision 1.1.1.1  2004/10/26 06:53:56  jgs  
 # initial import of project esys2  
 #  
 # Revision 1.1.2.3  2004/10/26 06:43:48  jgs  
 # committing Lutz's and Paul's changes to brach jgs  
 #  
 # Revision 1.1.4.1  2004/10/20 05:32:51  cochrane  
 # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.  
 #  
 # Revision 1.1  2004/08/05 03:58:27  gross  
 # Bug in Assemble_NodeCoordinates fixed  
 #  
 #  

Legend:
Removed from v.123  
changed lines
  Added in v.154

  ViewVC Help
Powered by ViewVC 1.1.26