/[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 148 by jgs, Tue Aug 23 01:24:31 2005 UTC revision 153 by jgs, Tue Oct 25 01:51:20 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:
9    
10    * binary operations @:               (a@b)[:,*]=a[:]@b[:,*] if rank(a)<rank(b)    - binary operations @ (@=+,-,*,/,**)::
11                  @=+,-,*,/,**           (a@b)[:]=a[:]@b[:] if rank(a)=rank(b)        (a@b)[:,*]=a[:]@b[:,*] if rank(a)<rank(b)
12                                         (a@b)[*,:]=a[*,:]@b[:] if rank(a)>rank(b)        (a@b)[:]=a[:]@b[:] if rank(a)=rank(b)
13            (a@b)[*,:]=a[*,:]@b[:] if rank(a)>rank(b)
14    * implementation of outer outer(a,b)[:,*]=a[:]*b[*]    - implementation of outer::
15    * trace: trace(arg,axis0=a0,axis1=a1)(:,&,*)=sum_i trace(:,i,&,i,*) (i are at index a0 and a1)        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  #===========================================================  import os
 # 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.
           d=a.getDim()  
           if dim==None:  
              dim=d  
           else:  
              if dim!=d: raise ValueError,"inconsistent spatial dimension of arguments"  
     if dim==None:  
        raise ValueError,"cannot recover spatial dimension"  
     return dim  
   
 def _identifyShape(arg):  
    """identifies the shape of arg."""  
    if hasattr(arg,"getShape"):  
        arg_shape=arg.getShape()  
    elif hasattr(arg,"shape"):  
      s=arg.shape  
      if callable(s):  
        arg_shape=s()  
      else:  
        arg_shape=s  
    else:  
        arg_shape=()  
    return arg_shape  
47    
48  def _extractShape(args):      Example:
49      """extracts the common shape of the list of arguments args"""  
50      shape=None         tmp=Scalar(..)
51      for a in args:         v=Vector(..)
52         a_shape=_identifyShape(a)         saveVTK("solution.xml",temperature=tmp,velovity=v)
53         if shape==None: shape=a_shape  
54         if shape!=a_shape: raise ValueError,"inconsistent shape"      tmp and v are written into "solution.dx" where tmp is named "temperature" and v is named "velovity"
     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    
55    
56        @param filename: file name of the output file
57        @type filename: C(str}
58        @param domain: domain of the L{Data} object. If not specified, the domain of the given L{Data} objects is used.
59        @type domain: L{escript.Domain}
60        @keyword <name>: writes the assigned value to the VTK file using <name> as identifier.
61        @type <name>: L{Data} object.
62        @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        """
64        if domain==None:
65           for i in data.keys():
66              if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()
67        if domain==None:
68            raise ValueError,"no domain detected."
69        else:
70            domain.saveVTK(filename,data)
71  #=========================================================  #=========================================================
72  #   wrappers for various mathematical functions:  def saveDX(filename,domain=None,**data):
73        """
74        writes a L{Data} objects into a files using the the DX file format.
75    
76        Example:
77    
78           tmp=Scalar(..)
79           v=Vector(..)
80           saveDX("solution.dx",temperature=tmp,velovity=v)
81    
82        tmp and v are written into "solution.dx" where tmp is named "temperature" and v is named "velovity".
83    
84        @param filename: file name of the output file
85        @type filename: C(str}
86        @param domain: domain of the L{Data} object. If not specified, the domain of the given L{Data} objects is used.
87        @type domain: L{escript.Domain}
88        @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.
89        @type <name>: L{Data} object.
90        @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.
91        """
92        if domain==None:
93           for i in data.keys():
94              if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()
95        if domain==None:
96            raise ValueError,"no domain detected."
97        else:
98            domain.saveDX(filename,data)
99  #=========================================================  #=========================================================
 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  
100    
101  def exp(arg):  def exp(arg):
102      """      """
103      @brief applies the exponential function to arg      Applies the exponential function to arg.
104      @param arg (input): argument  
105        @param arg: argument
106      """      """
107      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
108         return Exp_Symbol(arg)         return symbols.Exp_Symbol(arg)
109      elif hasattr(arg,"exp"):      elif hasattr(arg,"exp"):
110         return arg.exp()         return arg.exp()
111      else:      else:
112         return numarray.exp(arg)         return numarray.exp(arg)
113    
 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]  
   
114  def sqrt(arg):  def sqrt(arg):
115      """      """
116      @brief applies the squre root function to arg      Applies the squre root function to arg.
117      @param arg (input): argument  
118        @param arg: argument
119      """      """
120      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
121         return Sqrt_Symbol(arg)         return symbols.Sqrt_Symbol(arg)
122      elif hasattr(arg,"sqrt"):      elif hasattr(arg,"sqrt"):
123         return arg.sqrt()         return arg.sqrt()
124      else:      else:
125         return numarray.sqrt(arg)               return numarray.sqrt(arg)      
126    
 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]  
   
127  def log(arg):  def log(arg):
128      """      """
129      @brief applies the logarithmic function bases exp(1.) to arg      Applies the logarithmic function base 10 to arg.
130      @param arg (input): argument  
131        @param arg: argument
132      """      """
133      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
134         return Log_Symbol(arg)         return symbols.Log_Symbol(arg)
135      elif hasattr(arg,"log"):      elif hasattr(arg,"log"):
136         return arg.log()         return arg.log()
137      else:      else:
138         return numarray.log(arg)         return numarray.log10(arg)
   
 class Log_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 "log(%s)"%str(self.getArgument(0))  
    def eval(self,argval):  
        return log(self.getEvaluatedArguments(argval)[0])  
    def _diff(self,arg):  
        return self.getDifferentiatedArguments(arg)[0]/self.getArgument(0)  
139    
140  def ln(arg):  def ln(arg):
141      """      """
142      @brief applies the natural logarithmic function to arg      Applies the natural logarithmic function to arg.
143      @param arg (input): argument  
144        @param arg: argument
145      """      """
146      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
147         return Ln_Symbol(arg)         return symbols.Ln_Symbol(arg)
148      elif hasattr(arg,"ln"):      elif hasattr(arg,"ln"):
149         return arg.log()         return arg.ln()
150      else:      else:
151         return numarray.log(arg)         return numarray.log(arg)
152    
 class Ln_Symbol(Symbol):  
    """symbol representing natural logarithm of the argument"""  
    def __init__(self,arg):  
        Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])  
    def __str__(self):  
       return "ln(%s)"%str(self.getArgument(0))  
    def eval(self,argval):  
        return ln(self.getEvaluatedArguments(argval)[0])  
    def _diff(self,arg):  
        return self.getDifferentiatedArguments(arg)[0]/self.getArgument(0)  
   
153  def sin(arg):  def sin(arg):
154      """      """
155      @brief applies the sin function to arg      Applies the sin function to arg.
156      @param arg (input): argument  
157        @param arg: argument
158      """      """
159      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
160         return Sin_Symbol(arg)         return symbols.Sin_Symbol(arg)
161      elif hasattr(arg,"sin"):      elif hasattr(arg,"sin"):
162         return arg.sin()         return arg.sin()
163      else:      else:
164         return numarray.sin(arg)         return numarray.sin(arg)
165    
 class Sin_Symbol(Symbol):  
    """symbol representing sin 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]  
   
166  def cos(arg):  def cos(arg):
167      """      """
168      @brief applies the cos function to arg      Applies the cos function to arg.
169      @param arg (input): argument  
170        @param arg: argument
171      """      """
172      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
173         return Cos_Symbol(arg)         return symbols.Cos_Symbol(arg)
174      elif hasattr(arg,"cos"):      elif hasattr(arg,"cos"):
175         return arg.cos()         return arg.cos()
176      else:      else:
177         return numarray.cos(arg)         return numarray.cos(arg)
178    
 class Cos_Symbol(Symbol):  
    """symbol representing cos 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]  
   
179  def tan(arg):  def tan(arg):
180      """      """
181      @brief applies the tan function to arg      Applies the tan function to arg.
182      @param arg (input): argument  
183        @param arg: argument
184      """      """
185      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
186         return Tan_Symbol(arg)         return symbols.Tan_Symbol(arg)
187      elif hasattr(arg,"tan"):      elif hasattr(arg,"tan"):
188         return arg.tan()         return arg.tan()
189      else:      else:
190         return numarray.tan(arg)         return numarray.tan(arg)
191    
192  class Tan_Symbol(Symbol):  def asin(arg):
193     """symbol representing tan of the argument"""      """
194     def __init__(self,arg):      Applies the asin function to arg.
195         Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])  
196     def __str__(self):      @param arg: argument
197        return "tan(%s)"%str(self.getArgument(0))      """
198     def eval(self,argval):      if isinstance(arg,symbols.Symbol):
199         return tan(self.getEvaluatedArguments(argval)[0])         return symbols.Asin_Symbol(arg)
200     def _diff(self,arg):      elif hasattr(arg,"asin"):
201         s=cos(self.getArgument(0))         return arg.asin()
202         return 1./(s*s)*self.getDifferentiatedArguments(arg)[0]      else:
203           return numarray.asin(arg)
204    
205    def acos(arg):
206        """
207        Applies the acos function to arg.
208    
209        @param arg: argument
210        """
211        if isinstance(arg,symbols.Symbol):
212           return symbols.Acos_Symbol(arg)
213        elif hasattr(arg,"acos"):
214           return arg.acos()
215        else:
216           return numarray.acos(arg)
217    
218    def atan(arg):
219        """
220        Applies the atan function to arg.
221    
222        @param arg: argument
223        """
224        if isinstance(arg,symbols.Symbol):
225           return symbols.Atan_Symbol(arg)
226        elif hasattr(arg,"atan"):
227           return arg.atan()
228        else:
229           return numarray.atan(arg)
230    
231    def sinh(arg):
232        """
233        Applies the sinh function to arg.
234    
235        @param arg: argument
236        """
237        if isinstance(arg,symbols.Symbol):
238           return symbols.Sinh_Symbol(arg)
239        elif hasattr(arg,"sinh"):
240           return arg.sinh()
241        else:
242           return numarray.sinh(arg)
243    
244    def cosh(arg):
245        """
246        Applies the cosh function to arg.
247    
248        @param arg: argument
249        """
250        if isinstance(arg,symbols.Symbol):
251           return symbols.Cosh_Symbol(arg)
252        elif hasattr(arg,"cosh"):
253           return arg.cosh()
254        else:
255           return numarray.cosh(arg)
256    
257    def tanh(arg):
258        """
259        Applies the tanh function to arg.
260    
261        @param arg: argument
262        """
263        if isinstance(arg,symbols.Symbol):
264           return symbols.Tanh_Symbol(arg)
265        elif hasattr(arg,"tanh"):
266           return arg.tanh()
267        else:
268           return numarray.tanh(arg)
269    
270    def asinh(arg):
271        """
272        Applies the asinh function to arg.
273    
274        @param arg: argument
275        """
276        if isinstance(arg,symbols.Symbol):
277           return symbols.Asinh_Symbol(arg)
278        elif hasattr(arg,"asinh"):
279           return arg.asinh()
280        else:
281           return numarray.asinh(arg)
282    
283    def acosh(arg):
284        """
285        Applies the acosh function to arg.
286    
287        @param arg: argument
288        """
289        if isinstance(arg,symbols.Symbol):
290           return symbols.Acosh_Symbol(arg)
291        elif hasattr(arg,"acosh"):
292           return arg.acosh()
293        else:
294           return numarray.acosh(arg)
295    
296    def atanh(arg):
297        """
298        Applies the atanh function to arg.
299    
300        @param arg: argument
301        """
302        if isinstance(arg,symbols.Symbol):
303           return symbols.Atanh_Symbol(arg)
304        elif hasattr(arg,"atanh"):
305           return arg.atanh()
306        else:
307           return numarray.atanh(arg)
308    
309  def sign(arg):  def sign(arg):
310      """      """
311      @brief applies the sign function to arg      Applies the sign function to arg.
312      @param arg (input): argument  
313        @param arg: argument
314      """      """
315      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
316         return Sign_Symbol(arg)         return symbols.Sign_Symbol(arg)
317      elif hasattr(arg,"sign"):      elif hasattr(arg,"sign"):
318         return arg.sign()         return arg.sign()
319      else:      else:
320         return numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))- \         return numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))- \
321                numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))                numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))
322    
 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])  
   
323  def maxval(arg):  def maxval(arg):
324      """      """
325      @brief returns the maximum value of argument arg""      Returns the maximum value of argument arg.
326      @param arg (input): argument  
327        @param arg: argument
328      """      """
329      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
330         return Max_Symbol(arg)         return symbols.Max_Symbol(arg)
331      elif hasattr(arg,"maxval"):      elif hasattr(arg,"maxval"):
332         return arg.maxval()         return arg.maxval()
333      elif hasattr(arg,"max"):      elif hasattr(arg,"max"):
# Line 617  def maxval(arg): Line 335  def maxval(arg):
335      else:      else:
336         return arg         return arg
337    
 class Max_Symbol(Symbol):  
    """symbol representing the maximum value 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])  
   
338  def minval(arg):  def minval(arg):
339      """      """
340      @brief returns the minimum value of argument arg""      Returns the minimum value of argument arg.
341      @param arg (input): argument  
342        @param arg: argument
343      """      """
344      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
345         return Min_Symbol(arg)         return symbols.Min_Symbol(arg)
346      elif hasattr(arg,"maxval"):      elif hasattr(arg,"maxval"):
347         return arg.minval()         return arg.minval()
348      elif hasattr(arg,"min"):      elif hasattr(arg,"min"):
# Line 640  def minval(arg): Line 350  def minval(arg):
350      else:      else:
351         return arg         return arg
352    
 class Min_Symbol(Symbol):  
    """symbol representing the minimum value 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])  
   
353  def wherePositive(arg):  def wherePositive(arg):
354      """      """
355      @brief returns the positive values of argument arg""      Returns the positive values of argument arg.
356      @param arg (input): argument  
357        @param arg: argument
358      """      """
359      if _testForZero(arg):      if _testForZero(arg):
360        return 0        return 0
361      elif isinstance(arg,Symbol):      elif isinstance(arg,symbols.Symbol):
362         return WherePositive_Symbol(arg)         return symbols.WherePositive_Symbol(arg)
363      elif hasattr(arg,"wherePositive"):      elif hasattr(arg,"wherePositive"):
364         return arg.minval()         return arg.minval()
365      elif hasattr(arg,"wherePositive"):      elif hasattr(arg,"wherePositive"):
# Line 668  def wherePositive(arg): Line 370  def wherePositive(arg):
370         else:         else:
371            return 0.            return 0.
372    
 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])  
   
373  def whereNegative(arg):  def whereNegative(arg):
374      """      """
375      @brief returns the negative values of argument arg""      Returns the negative values of argument arg.
376      @param arg (input): argument  
377        @param arg: argument
378      """      """
379      if _testForZero(arg):      if _testForZero(arg):
380        return 0        return 0
381      elif isinstance(arg,Symbol):      elif isinstance(arg,symbols.Symbol):
382         return WhereNegative_Symbol(arg)         return symbols.WhereNegative_Symbol(arg)
383      elif hasattr(arg,"whereNegative"):      elif hasattr(arg,"whereNegative"):
384         return arg.whereNegative()         return arg.whereNegative()
385      elif hasattr(arg,"shape"):      elif hasattr(arg,"shape"):
# Line 696  def whereNegative(arg): Line 390  def whereNegative(arg):
390         else:         else:
391            return 0.            return 0.
392    
 class WhereNegative_Symbol(Symbol):  
    """symbol representing the whereNegative function"""  
    def __init__(self,arg):  
        Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])  
    def __str__(self):  
       return "whereNegative(%s)"%str(self.getArgument(0))  
    def eval(self,argval):  
        return whereNegative(self.getEvaluatedArguments(argval)[0])  
   
393  def maximum(arg0,arg1):  def maximum(arg0,arg1):
394      """return arg1 where arg1 is bigger then arg0 otherwise arg0 is returned"""      """
395        Return arg1 where arg1 is bigger then arg0 otherwise arg0 is returned.
396        """
397      m=whereNegative(arg0-arg1)      m=whereNegative(arg0-arg1)
398      return m*arg1+(1.-m)*arg0      return m*arg1+(1.-m)*arg0
399        
400  def minimum(arg0,arg1):  def minimum(arg0,arg1):
401      """return arg0 where arg1 is bigger then arg0 otherwise arg1 is returned"""      """
402        Return arg0 where arg1 is bigger then arg0 otherwise arg1 is returned.
403        """
404      m=whereNegative(arg0-arg1)      m=whereNegative(arg0-arg1)
405      return m*arg0+(1.-m)*arg1      return m*arg0+(1.-m)*arg1
406        
# Line 719  def outer(arg0,arg1): Line 408  def outer(arg0,arg1):
408     if _testForZero(arg0) or _testForZero(arg1):     if _testForZero(arg0) or _testForZero(arg1):
409        return 0        return 0
410     else:     else:
411        if isinstance(arg0,Symbol) or isinstance(arg1,Symbol):        if isinstance(arg0,symbols.Symbol) or isinstance(arg1,symbols.Symbol):
412          return Outer_Symbol(arg0,arg1)          return symbols.Outer_Symbol(arg0,arg1)
413        elif _identifyShape(arg0)==() or _identifyShape(arg1)==():        elif _identifyShape(arg0)==() or _identifyShape(arg1)==():
414          return arg0*arg1          return arg0*arg1
415        elif isinstance(arg0,numarray.NumArray) and isinstance(arg1,numarray.NumArray):        elif isinstance(arg0,numarray.NumArray) and isinstance(arg1,numarray.NumArray):
# Line 735  def outer(arg0,arg1): Line 424  def outer(arg0,arg1):
424          else:          else:
425            raise ValueError,"outer is not fully implemented yet."            raise ValueError,"outer is not fully implemented yet."
426    
 class Outer_Symbol(Symbol):  
    """symbol representing the outer product of its two arguments"""  
    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])  
   
427  def interpolate(arg,where):  def interpolate(arg,where):
428      """      """
429      @brief interpolates the function into the FunctionSpace where.      Interpolates the function into the FunctionSpace where.
430    
431      @param arg    interpolant      @param arg:    interpolant
432      @param where  FunctionSpace to interpolate to      @param where:  FunctionSpace to interpolate to
433      """      """
434      if _testForZero(arg):      if _testForZero(arg):
435        return 0        return 0
436      elif isinstance(arg,Symbol):      elif isinstance(arg,symbols.Symbol):
437         return Interpolated_Symbol(arg,where)         return symbols.Interpolated_Symbol(arg,where)
438      else:      else:
439         return escript.Data(arg,where)         return escript.Data(arg,where)
440    
 def Interpolated_Symbol(Symbol):  
    """symbol representing the integral of the argument"""  
    def __init__(self,arg,where):  
         Symbol.__init__(self,shape=_extractShape(arg),dim=_extractDim([arg]),args=[arg,where])  
    def __str__(self):  
       return "interpolated(%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))  
   
441  def div(arg,where=None):  def div(arg,where=None):
442      """      """
443      @brief returns the divergence of arg at where.      Returns the divergence of arg at where.
444    
445      @param arg:   Data object representing the function which gradient to be calculated.      @param arg:   Data object representing the function which gradient to
446      @param where: FunctionSpace in which the gradient will be calculated. If not present or                    be calculated.
447                    None an appropriate default is used.      @param where: FunctionSpace in which the gradient will be calculated.
448                      If not present or C{None} an appropriate default is used.
449      """      """
450      return trace(grad(arg,where))      return trace(grad(arg,where))
451    
452    def jump(arg):
453        """
454        Returns the jump of arg across a continuity.
455    
456        @param arg:   Data object representing the function which gradient
457                      to be calculated.
458        """
459        d=arg.getDomain()
460        return arg.interpolate(escript.FunctionOnContactOne())-arg.interpolate(escript.FunctionOnContactZero())
461      
462    
463  def grad(arg,where=None):  def grad(arg,where=None):
464      """      """
465      @brief returns the spatial gradient of arg at where.      Returns the spatial gradient of arg at where.
466    
467      @param arg:   Data object representing the function which gradient to be calculated.      @param arg:   Data object representing the function which gradient
468      @param where: FunctionSpace in which the gradient will be calculated. If not present or                    to be calculated.
469                    None an appropriate default is used.      @param where: FunctionSpace in which the gradient will be calculated.
470                      If not present or C{None} an appropriate default is used.
471      """      """
472      if _testForZero(arg):      if _testForZero(arg):
473        return 0        return 0
474      elif isinstance(arg,Symbol):      elif isinstance(arg,symbols.Symbol):
475         return Grad_Symbol(arg,where)         return symbols.Grad_Symbol(arg,where)
476      elif hasattr(arg,"grad"):      elif hasattr(arg,"grad"):
477         if where==None:         if where==None:
478            return arg.grad()            return arg.grad()
# Line 807  def grad(arg,where=None): Line 481  def grad(arg,where=None):
481      else:      else:
482         return arg*0.         return arg*0.
483    
 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))  
   
484  def integrate(arg,where=None):  def integrate(arg,where=None):
485      """      """
486      @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
487        its domain.
488    
489      @param arg:   Data object representing the function which is integrated.      @param arg:   Data object representing the function which is integrated.
490      @param where: FunctionSpace in which the integral is calculated. If not present or      @param where: FunctionSpace in which the integral is calculated.
491                    None an appropriate default is used.                    If not present or C{None} an appropriate default is used.
492      """      """
493      if _testForZero(arg):      if _testForZero(arg):
494        return 0        return 0
495      elif isinstance(arg,Symbol):      elif isinstance(arg,symbols.Symbol):
496         return Integral_Symbol(arg,where)         return symbols.Integral_Symbol(arg,where)
497      else:          else:    
498         if not where==None: arg=escript.Data(arg,where)         if not where==None: arg=escript.Data(arg,where)
499         if arg.getRank()==0:         if arg.getRank()==0:
# Line 841  def integrate(arg,where=None): Line 501  def integrate(arg,where=None):
501         else:         else:
502           return arg.integrate()           return arg.integrate()
503    
 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))  
   
504  #=============================  #=============================
505  #  #
506  # wrapper for various functions: if the argument has attribute the function name  # wrapper for various functions: if the argument has attribute the function name
# Line 867  def Integral_Symbol(Float_Symbol): Line 514  def Integral_Symbol(Float_Symbol):
514    
515  def transpose(arg,axis=None):  def transpose(arg,axis=None):
516      """      """
517      @brief returns the transpose of the Data object arg.      Returns the transpose of the Data object arg.
518    
519      @param arg      @param arg:
520      """      """
521      if axis==None:      if axis==None:
522         r=0         r=0
523         if hasattr(arg,"getRank"): r=arg.getRank()         if hasattr(arg,"getRank"): r=arg.getRank()
524         if hasattr(arg,"rank"): r=arg.rank         if hasattr(arg,"rank"): r=arg.rank
525         axis=r/2         axis=r/2
526      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
527         return Transpose_Symbol(arg,axis=r)         return symbols.Transpose_Symbol(arg,axis=r)
528      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
529         # hack for transpose         # hack for transpose
530         r=arg.getRank()         r=arg.getRank()
# Line 895  def transpose(arg,axis=None): Line 542  def transpose(arg,axis=None):
542    
543  def trace(arg,axis0=0,axis1=1):  def trace(arg,axis0=0,axis1=1):
544      """      """
545      @brief return      Return
546    
547      @param arg      @param arg:
548      """      """
549      if isinstance(arg,Symbol):      if isinstance(arg,symbols.Symbol):
550         s=list(arg.getShape())                 s=list(arg.getShape())        
551         s=tuple(s[0:axis0]+s[axis0+1:axis1]+s[axis1+1:])         s=tuple(s[0:axis0]+s[axis0+1:axis1]+s[axis1+1:])
552         return Trace_Symbol(arg,axis0=axis0,axis1=axis1)         return symbols.Trace_Symbol(arg,axis0=axis0,axis1=axis1)
553      elif isinstance(arg,escript.Data):      elif isinstance(arg,escript.Data):
554         # hack for trace         # hack for trace
555         s=arg.getShape()         s=arg.getShape()
# Line 916  def trace(arg,axis0=0,axis1=1): Line 563  def trace(arg,axis0=0,axis1=1):
563      else:      else:
564         return numarray.trace(arg,axis0=axis0,axis1=axis1)         return numarray.trace(arg,axis0=axis0,axis1=axis1)
565    
 def Trace_Symbol(Symbol):  
     pass  
   
566  def length(arg):  def length(arg):
567      """      """
     @brief  
568    
569      @param arg      @param arg:
570      """      """
571      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
572         if arg.isEmpty(): return escript.Data()         if arg.isEmpty(): return escript.Data()
# Line 965  def length(arg): Line 608  def length(arg):
608    
609  def deviator(arg):  def deviator(arg):
610      """      """
611      @brief      @param arg:
   
     @param arg0  
612      """      """
613      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
614          shape=arg.getShape()          shape=arg.getShape()
# Line 981  def deviator(arg): Line 622  def deviator(arg):
622    
623  def inner(arg0,arg1):  def inner(arg0,arg1):
624      """      """
625      @brief      @param arg0:
626        @param arg1:
     @param arg0, arg1  
627      """      """
628      if isinstance(arg0,escript.Data):      if isinstance(arg0,escript.Data):
629         arg=arg0         arg=arg0
# Line 1019  def inner(arg0,arg1): Line 659  def inner(arg0,arg1):
659            raise SystemError,"inner is not been implemented yet"            raise SystemError,"inner is not been implemented yet"
660      return out      return out
661    
662    def tensormult(arg0,arg1):
663        # check LinearPDE!!!!
664        raise SystemError,"tensormult is not implemented yet!"
665    
666  def matrixmult(arg0,arg1):  def matrixmult(arg0,arg1):
667    
668      if isinstance(arg1,numarray.NumArray) and isinstance(arg0,numarray.NumArray):      if isinstance(arg1,numarray.NumArray) and isinstance(arg0,numarray.NumArray):
# Line 1046  def matrixmult(arg0,arg1): Line 690  def matrixmult(arg0,arg1):
690  #=========================================================  #=========================================================
691  def sum(arg):  def sum(arg):
692      """      """
693      @brief      @param arg:
   
     @param arg  
694      """      """
695      return arg.sum()      return arg.sum()
696    
697  def sup(arg):  def sup(arg):
698      """      """
699      @brief      @param arg:
   
     @param arg  
700      """      """
701      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
702         return arg.sup()         return arg.sup()
# Line 1067  def sup(arg): Line 707  def sup(arg):
707    
708  def inf(arg):  def inf(arg):
709      """      """
710      @brief      @param arg:
   
     @param arg  
711      """      """
712      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
713         return arg.inf()         return arg.inf()
# Line 1080  def inf(arg): Line 718  def inf(arg):
718    
719  def L2(arg):  def L2(arg):
720      """      """
721      @brief returns the L2-norm of the      Returns the L2-norm of the argument
722    
723      @param arg      @param arg:
724      """      """
725      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
726         return arg.L2()         return arg.L2()
# Line 1093  def L2(arg): Line 731  def L2(arg):
731    
732  def Lsup(arg):  def Lsup(arg):
733      """      """
734      @brief      @param arg:
   
     @param arg  
735      """      """
736      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
737         return arg.Lsup()         return arg.Lsup()
738      elif isinstance(arg,float) or isinstance(arg,int):      elif isinstance(arg,float) or isinstance(arg,int):
739         return abs(arg)         return abs(arg)
740      else:      else:
741         return max(numarray.abs(arg))         return numarray.abs(arg).max()
742    
743  def dot(arg0,arg1):  def dot(arg0,arg1):
744      """      """
745      @brief      @param arg0:
746        @param arg1:
     @param arg  
747      """      """
748      if isinstance(arg0,escript.Data):      if isinstance(arg0,escript.Data):
749         return arg0.dot(arg1)         return arg0.dot(arg1)
# Line 1125  def kronecker(d): Line 760  def kronecker(d):
760    
761  def unit(i,d):  def unit(i,d):
762     """     """
763     @brief return a unit vector of dimension d with nonzero index i     Return a unit vector of dimension d with nonzero index i.
764     @param d dimension  
765     @param i index     @param d: dimension
766       @param i: index
767     """     """
768     e = numarray.zeros((d,),numarray.Float)     e = numarray.zeros((d,),numarray.Float)
769     e[i] = 1.0     e[i] = 1.0
770     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.16  2005/08/23 01:24:28  jgs  
 # Merge of development branch dev-02 back to main trunk on 2005-08-23  
 #  
 # Revision 1.15  2005/08/12 01:45:36  jgs  
 # erge of development branch dev-02 back to main trunk on 2005-08-12  
 #  
 # Revision 1.14.2.4  2005/08/18 04:39:32  gross  
 # the constants have been removed from util.py as they not needed anymore. PDE related constants are accessed through LinearPDE attributes now  
 #  
 # Revision 1.14.2.3  2005/08/03 09:55:33  gross  
 # ContactTest is passing now./mk install!  
 #  
 # Revision 1.14.2.2  2005/08/02 03:15:14  gross  
 # bug inb trace fixed!  
 #  
 # Revision 1.14.2.1  2005/07/29 07:10:28  gross  
 # new functions in util and a new pde type in linearPDEs  
 #  
 # Revision 1.2.2.21  2005/07/28 04:19:23  gross  
 # new functions maximum and minimum introduced.  
 #  
 # Revision 1.2.2.20  2005/07/25 01:26:27  gross  
 # bug in inner fixed  
 #  
 # Revision 1.2.2.19  2005/07/21 04:01:28  jgs  
 # minor comment fixes  
 #  
 # Revision 1.2.2.18  2005/07/21 01:02:43  jgs  
 # commit ln() updates to development branch version  
 #  
 # Revision 1.12  2005/07/20 06:14:58  jgs  
 # added ln(data) style wrapper for data.ln() - also added corresponding  
 # implementation of Ln_Symbol class (not sure if this is right though)  
 #  
 # 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.148  
changed lines
  Added in v.153

  ViewVC Help
Powered by ViewVC 1.1.26