/[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 1247 by ksteube, Tue Aug 14 01:29:20 2007 UTC revision 1719 by gross, Thu Aug 21 06:24:29 2008 UTC
# Line 1  Line 1 
1    #
2  # $Id$  # $Id$
3    #
4    #######################################################
5    #
6    #           Copyright 2003-2007 by ACceSS MNRF
7    #       Copyright 2007 by University of Queensland
8    #
9    #                http://esscc.uq.edu.au
10    #        Primary Business: Queensland, Australia
11    #  Licensed under the Open Software License version 3.0
12    #     http://www.opensource.org/licenses/osl-3.0.php
13    #
14    #######################################################
15    #
16    
17  """  """
18  Utility functions for escript  Utility functions for escript
# Line 9  Utility functions for escript Line 23  Utility functions for escript
23  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
24  @var __version__: version  @var __version__: version
25  @var __date__: date of the version  @var __date__: date of the version
26    @var EPSILON: smallest positive value with 1.<1.+EPSILON
27  """  """
28                                                                                                                                                                                                                                                                                                                                                                                                            
29  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
# Line 28  import escript Line 43  import escript
43  import os  import os
44  from esys.escript import C_GeneralTensorProduct  from esys.escript import C_GeneralTensorProduct
45  from esys.escript import getVersion  from esys.escript import getVersion
46    from esys.escript import printParallelThreadCounts
47    
48  #=========================================================  #=========================================================
49  #   some helpers:  #   some helpers:
50  #=========================================================  #=========================================================
51    def getEpsilon():
52         #     ------------------------------------------------------------------
53         #     Compute EPSILON, the machine precision.  The call to daxpp is
54         #     inTENded to fool compilers that use extra-length registers.
55         #     31 Map 1999: Hardwire EPSILON so the debugger can step thru easily.
56         #     ------------------------------------------------------------------
57         eps    = 2.**(-12)
58         p=2.
59         while p>1.:
60                eps/=2.
61                p=1.+eps
62         return eps*2.
63    
64    EPSILON=getEpsilon()
65    
66  def getTagNames(domain):  def getTagNames(domain):
67      """      """
68      returns a list of the tag names used by the domain      returns a list of the tag names used by the domain
# Line 71  def insertTaggedValues(target,**kwargs): Line 102  def insertTaggedValues(target,**kwargs):
102          target.setTaggedValue(k,kwargs[k])          target.setTaggedValue(k,kwargs[k])
103      return target      return target
104    
       
105  def saveVTK(filename,domain=None,**data):  def saveVTK(filename,domain=None,**data):
106      """      """
107      writes a L{Data} objects into a files using the the VTK XML file format.      writes a L{Data} objects into a files using the the VTK XML file format.
# Line 80  def saveVTK(filename,domain=None,**data) Line 110  def saveVTK(filename,domain=None,**data)
110    
111         tmp=Scalar(..)         tmp=Scalar(..)
112         v=Vector(..)         v=Vector(..)
113         saveVTK("solution.xml",temperature=tmp,velovity=v)         saveVTK("solution.xml",temperature=tmp,velocity=v)
114    
115      tmp and v are written into "solution.xml" where tmp is named "temperature" and v is named "velovity"      tmp and v are written into "solution.xml" where tmp is named "temperature" and v is named "velocity"
116    
117      @param filename: file name of the output file      @param filename: file name of the output file
118      @type filename: C{str}      @type filename: C{str}
# Line 92  def saveVTK(filename,domain=None,**data) Line 122  def saveVTK(filename,domain=None,**data)
122      @type <name>: L{Data} object.      @type <name>: L{Data} object.
123      @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.      @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.
124      """      """
125      if domain==None:      new_data={}
126         for i in data.keys():      for n,d in data.items():
127            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
128                fs=d.getFunctionSpace()
129                domain2=fs.getDomain()
130                if fs == escript.Solution(domain2):
131                   new_data[n]=interpolate(d,escript.ContinuousFunction(domain2))
132                elif fs == escript.ReducedSolution(domain2):
133                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
134                else:
135                   new_data[n]=d
136                if domain==None: domain=domain2
137      if domain==None:      if domain==None:
138          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
139      else:      domain.saveVTK(filename,new_data)
         domain.saveVTK(filename,data)  
140    
141  def saveDX(filename,domain=None,**data):  def saveDX(filename,domain=None,**data):
142      """      """
# Line 108  def saveDX(filename,domain=None,**data): Line 146  def saveDX(filename,domain=None,**data):
146    
147         tmp=Scalar(..)         tmp=Scalar(..)
148         v=Vector(..)         v=Vector(..)
149         saveDX("solution.dx",temperature=tmp,velovity=v)         saveDX("solution.dx",temperature=tmp,velocity=v)
150    
151      tmp and v are written into "solution.dx" where tmp is named "temperature" and v is named "velovity".      tmp and v are written into "solution.dx" where tmp is named "temperature" and v is named "velocity".
152    
153      @param filename: file name of the output file      @param filename: file name of the output file
154      @type filename: C{str}      @type filename: C{str}
# Line 120  def saveDX(filename,domain=None,**data): Line 158  def saveDX(filename,domain=None,**data):
158      @type <name>: L{Data} object.      @type <name>: L{Data} object.
159      @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.      @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.
160      """      """
161      if domain==None:      new_data={}
162         for i in data.keys():      for n,d in data.items():
163            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
164                fs=d.getFunctionSpace()
165                domain2=fs.getDomain()
166                if fs == escript.Solution(domain2):
167                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
168                elif fs == escript.ReducedSolution(domain2):
169                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
170                elif fs == escript.ContinuousFunction(domain2):
171                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
172                else:
173                   new_data[n]=d
174                if domain==None: domain=domain2
175      if domain==None:      if domain==None:
176          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
177      else:      domain.saveDX(filename,new_data)
         domain.saveDX(filename,data)  
178    
179  def kronecker(d=3):  def kronecker(d=3):
180     """     """
# Line 3768  class Add_Symbol(DependendSymbol): Line 3816  class Add_Symbol(DependendSymbol):
3816              raise TypeError,"%s: new value is not appropriate."%str(self)              raise TypeError,"%s: new value is not appropriate."%str(self)
3817        else:        else:
3818           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
3819           return add(args[0],args[1])           out=add(args[0],args[1])
3820             return out
3821    
3822     def diff(self,arg):     def diff(self,arg):
3823        """        """
# Line 4909  def integrate(arg,where=None): Line 4958  def integrate(arg,where=None):
4958         else:         else:
4959            return arg._integrate()            return arg._integrate()
4960      else:      else:
4961        raise TypeError,"integrate: Unknown argument type."         arg2=escript.Data(arg,where)
4962           if arg2.getRank()==0:
4963              return arg2._integrate()[0]
4964           else:
4965              return arg2._integrate()
4966    
4967  class Integrate_Symbol(DependendSymbol):  class Integrate_Symbol(DependendSymbol):
4968     """     """
# Line 4981  class Integrate_Symbol(DependendSymbol): Line 5034  class Integrate_Symbol(DependendSymbol):
5034    
5035  def interpolate(arg,where):  def interpolate(arg,where):
5036      """      """
5037      interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where. If the argument C{arg} has the requested function space
5038        C{where} no interpolation is performed and C{arg} is returned.
5039    
5040      @param arg: interpolant      @param arg: interpolant
5041      @type arg: L{escript.Data} or L{Symbol}      @type arg: L{escript.Data} or L{Symbol}
# Line 4992  def interpolate(arg,where): Line 5046  def interpolate(arg,where):
5046      """      """
5047      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5048         return Interpolate_Symbol(arg,where)         return Interpolate_Symbol(arg,where)
5049        elif isinstance(arg,escript.Data):
5050           if where == arg.getFunctionSpace():
5051              return arg
5052           else:
5053              return escript.Data(arg,where)
5054      else:      else:
5055         return escript.Data(arg,where)         return escript.Data(arg,where)
5056    
# Line 5110  def L2(arg): Line 5169  def L2(arg):
5169      @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))      @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5170      """      """
5171      return sqrt(integrate(inner(arg,arg)))      return sqrt(integrate(inner(arg,arg)))
5172    
5173    def getClosestValue(arg,origin=0):
5174        """
5175        returns the value in arg which is closest to origin
5176        
5177        @param arg: function
5178        @type arg: L{escript.Data} or L{Symbol}
5179        @param origin: reference value
5180        @type origin: C{float} or L{escript.Data}
5181        @return: value in arg closest to origin
5182        @rtype: L{numarray.NumArray} or L{Symbol}
5183        """
5184        return arg.getValueOfGlobalDataPoint(*(length(arg-origin).minGlobalDataPoint()))
5185    
5186    def normalize(arg,zerolength=0):
5187        """
5188        returns normalized version of arg (=arg/length(arg))
5189        
5190        @param arg: function
5191        @type arg: L{escript.Data} or L{Symbol}
5192        @param zerolength: realitive tolerance for arg == 0.
5193        @type zerolength: C{float}
5194        @return: normalized arg where arg is non zero and zero elsewhere
5195        @rtype: L{escript.Data} or L{Symbol}
5196        """
5197        l=length(arg)
5198        m=whereZero(l,zerolength*Lsup(l))
5199        mm=1-m
5200        return arg*(mm/(l*mm+m))
5201    
5202  #=============================  #=============================
5203  #  #
5204    
# Line 5119  def reorderComponents(arg,index): Line 5208  def reorderComponents(arg,index):
5208    
5209      """      """
5210      raise NotImplementedError      raise NotImplementedError
 #  
 # $Log: util.py,v $  
 # Revision 1.14.2.16  2005/10/19 06:09:57  gross  
 # new versions of saveVTK and saveDX which allow to write several Data objects into a single file  
 #  
 # Revision 1.14.2.15  2005/10/19 03:25:27  jgs  
 # numarray uses log10 for log, and log for ln - log()/ln() updated to reflect this  
 #  
 # Revision 1.14.2.14  2005/10/19 02:10:23  jgs  
 # fixed ln() to call Data::ln  
 #  
 # Revision 1.14.2.13  2005/09/12 03:32:14  gross  
 # test_visualiztion has been aded to mk  
 #  
 # Revision 1.14.2.12  2005/09/09 01:56:24  jgs  
 # added implementations of acos asin atan sinh cosh tanh asinh acosh atanh  
 # and some associated testing  
 #  
 # Revision 1.14.2.11  2005/09/08 08:28:39  gross  
 # some cleanup in savevtk  
 #  
 # Revision 1.14.2.10  2005/09/08 00:25:32  gross  
 # test for finley mesh generators added  
 #  
 # Revision 1.14.2.9  2005/09/07 10:32:05  gross  
 # Symbols removed from util and put into symmbols.py.  
 #  
 # Revision 1.14.2.8  2005/08/26 05:06:37  cochrane  
 # Corrected errors in docstrings.  Improved output formatting of docstrings.  
 # Other minor improvements to code and docs (eg spelling etc).  
 #  
 # Revision 1.14.2.7  2005/08/26 04:45:40  cochrane  
 # Fixed and tidied markup and docstrings.  Some *Symbol classes were defined  
 # as functions, so changed them to classes (hopefully this was the right thing  
 # to do).  
 #  
 # Revision 1.14.2.6  2005/08/26 04:30:13  gross  
 # gneric unit testing for linearPDE  
 #  
 # Revision 1.14.2.5  2005/08/24 02:02:52  gross  
 # jump function added  
 #  
 # 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  
 #  
 #  
   
 # vim: expandtab shiftwidth=4:  

Legend:
Removed from v.1247  
changed lines
  Added in v.1719

  ViewVC Help
Powered by ViewVC 1.1.26