/[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 1042 by gross, Mon Mar 19 03:50:34 2007 UTC revision 2100 by gross, Wed Nov 26 08:13:00 2008 UTC
# Line 1  Line 1 
1  # $Id$  
2    ########################################################
3    #
4    # Copyright (c) 2003-2008 by University of Queensland
5    # Earth Systems Science Computational Center (ESSCC)
6    # http://www.uq.edu.au/esscc
7    #
8    # Primary Business: Queensland, Australia
9    # Licensed under the Open Software License version 3.0
10    # http://www.opensource.org/licenses/osl-3.0.php
11    #
12    ########################################################
13    
14    __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
17    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19    http://www.opensource.org/licenses/osl-3.0.php"""
20    __url__="http://www.uq.edu.au/esscc/escript-finley"
21    
22  """  """
23  Utility functions for escript  Utility functions for escript
# Line 9  Utility functions for escript Line 28  Utility functions for escript
28  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
29  @var __version__: version  @var __version__: version
30  @var __date__: date of the version  @var __date__: date of the version
31    @var EPSILON: smallest positive value with 1.<1.+EPSILON
32    @var DBLE_MAX: largest positive float
33  """  """
34                                                                                                                                                                                                        
35  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2006 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys/escript"  
 __version__="$Revision$"  
 __date__="$Date$"  
36    
37    
38  import math  import math
# Line 27  import numarray Line 40  import numarray
40  import escript  import escript
41  import os  import os
42  from esys.escript import C_GeneralTensorProduct  from esys.escript import C_GeneralTensorProduct
43    from esys.escript import getVersion
44    from esys.escript import printParallelThreadCounts
45    
46  #=========================================================  #=========================================================
47  #   some helpers:  #   some helpers:
48  #=========================================================  #=========================================================
49    def getEpsilon():
50         return escript.getMachinePrecision()
51    EPSILON=getEpsilon()
52    
53    def getMaxFloat():
54         return escript.getMaxFloat()
55    DBLE_MAX=getMaxFloat()
56    
57    
58  def getTagNames(domain):  def getTagNames(domain):
59      """      """
60      returns a list of the tag names used by the domain      returns a list of the tag names used by the domain
61    
62            
63      @param domain: a domain object      @param domain: a domain object
64      @type domain: C{escript.Domain}      @type domain: L{escript.Domain}
65      @return: a list of the tag name used by the domain.      @return: a list of the tag name used by the domain.
66      @rtype: C{list} of C{str}      @rtype: C{list} of C{str}
67      """      """
68      return [n.strip() for n in domain.showTagNames().split(",") ]      return [n.strip() for n in domain.showTagNames().split(",") ]
69        
70    def insertTagNames(domain,**kwargs):
71        """
72        inserts tag names into the domain
73    
74        @param domain: a domain object
75        @type domain: C{escript.Domain}
76        @keyword <tag_name>: tag key assigned to <tag_name>
77        @type <tag_name>: C{int}
78        """
79        for  k in kwargs:
80             domain.setTagMap(k,kwargs[k])
81    
82    def insertTaggedValues(target,**kwargs):
83        """
84        inserts tagged values into the tagged using tag names
85    
86        @param target: data to be filled by tagged values
87        @type target: L{escript.Data}
88        @keyword <tag_name>: value to be used for <tag_name>
89        @type <tag_name>: C{float} or {numarray.NumArray}
90        @return: C{target}
91        @rtype: L{escript.Data}
92        """
93        for k in kwargs:
94            target.setTaggedValue(k,kwargs[k])
95        return target
96    
97  def saveVTK(filename,domain=None,**data):  def saveVTK(filename,domain=None,**data):
98      """      """
99      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 51  def saveVTK(filename,domain=None,**data) Line 102  def saveVTK(filename,domain=None,**data)
102    
103         tmp=Scalar(..)         tmp=Scalar(..)
104         v=Vector(..)         v=Vector(..)
105         saveVTK("solution.xml",temperature=tmp,velovity=v)         saveVTK("solution.xml",temperature=tmp,velocity=v)
106    
107      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"
108    
109      @param filename: file name of the output file      @param filename: file name of the output file
110      @type filename: C{str}      @type filename: C{str}
# Line 63  def saveVTK(filename,domain=None,**data) Line 114  def saveVTK(filename,domain=None,**data)
114      @type <name>: L{Data} object.      @type <name>: L{Data} object.
115      @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.
116      """      """
117      if domain==None:      new_data={}
118         for i in data.keys():      for n,d in data.items():
119            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
120                fs=d.getFunctionSpace()
121                domain2=fs.getDomain()
122                if fs == escript.Solution(domain2):
123                   new_data[n]=interpolate(d,escript.ContinuousFunction(domain2))
124                elif fs == escript.ReducedSolution(domain2):
125                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
126                else:
127                   new_data[n]=d
128                if domain==None: domain=domain2
129      if domain==None:      if domain==None:
130          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
131      else:      domain.saveVTK(filename,new_data)
         domain.saveVTK(filename,data)  
132    
133  def saveDX(filename,domain=None,**data):  def saveDX(filename,domain=None,**data):
134      """      """
# Line 79  def saveDX(filename,domain=None,**data): Line 138  def saveDX(filename,domain=None,**data):
138    
139         tmp=Scalar(..)         tmp=Scalar(..)
140         v=Vector(..)         v=Vector(..)
141         saveDX("solution.dx",temperature=tmp,velovity=v)         saveDX("solution.dx",temperature=tmp,velocity=v)
142    
143      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".
144    
145      @param filename: file name of the output file      @param filename: file name of the output file
146      @type filename: C{str}      @type filename: C{str}
# Line 91  def saveDX(filename,domain=None,**data): Line 150  def saveDX(filename,domain=None,**data):
150      @type <name>: L{Data} object.      @type <name>: L{Data} object.
151      @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.
152      """      """
153      if domain==None:      new_data={}
154         for i in data.keys():      for n,d in data.items():
155            if not data[i].isEmpty(): domain=data[i].getFunctionSpace().getDomain()            if not d.isEmpty():
156                fs=d.getFunctionSpace()
157                domain2=fs.getDomain()
158                if fs == escript.Solution(domain2):
159                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
160                elif fs == escript.ReducedSolution(domain2):
161                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
162                elif fs == escript.ContinuousFunction(domain2):
163                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
164                else:
165                   new_data[n]=d
166                if domain==None: domain=domain2
167      if domain==None:      if domain==None:
168          raise ValueError,"no domain detected."          raise ValueError,"no domain detected."
169      else:      domain.saveDX(filename,new_data)
         domain.saveDX(filename,data)  
170    
171  def kronecker(d=3):  def kronecker(d=3):
172     """     """
# Line 3739  class Add_Symbol(DependendSymbol): Line 3808  class Add_Symbol(DependendSymbol):
3808              raise TypeError,"%s: new value is not appropriate."%str(self)              raise TypeError,"%s: new value is not appropriate."%str(self)
3809        else:        else:
3810           args=self.getSubstitutedArguments(argvals)           args=self.getSubstitutedArguments(argvals)
3811           return add(args[0],args[1])           out=add(args[0],args[1])
3812             return out
3813    
3814     def diff(self,arg):     def diff(self,arg):
3815        """        """
# Line 4880  def integrate(arg,where=None): Line 4950  def integrate(arg,where=None):
4950         else:         else:
4951            return arg._integrate()            return arg._integrate()
4952      else:      else:
4953        raise TypeError,"integrate: Unknown argument type."         arg2=escript.Data(arg,where)
4954           if arg2.getRank()==0:
4955              return arg2._integrate()[0]
4956           else:
4957              return arg2._integrate()
4958    
4959  class Integrate_Symbol(DependendSymbol):  class Integrate_Symbol(DependendSymbol):
4960     """     """
# Line 4952  class Integrate_Symbol(DependendSymbol): Line 5026  class Integrate_Symbol(DependendSymbol):
5026    
5027  def interpolate(arg,where):  def interpolate(arg,where):
5028      """      """
5029      interpolates the function into the FunctionSpace where.      interpolates the function into the FunctionSpace where. If the argument C{arg} has the requested function space
5030        C{where} no interpolation is performed and C{arg} is returned.
5031    
5032      @param arg: interpolant      @param arg: interpolant
5033      @type arg: L{escript.Data} or L{Symbol}      @type arg: L{escript.Data} or L{Symbol}
# Line 4963  def interpolate(arg,where): Line 5038  def interpolate(arg,where):
5038      """      """
5039      if isinstance(arg,Symbol):      if isinstance(arg,Symbol):
5040         return Interpolate_Symbol(arg,where)         return Interpolate_Symbol(arg,where)
5041        elif isinstance(arg,escript.Data):
5042           if arg.isEmpty():
5043              return arg
5044           elif where == arg.getFunctionSpace():
5045              return arg
5046           else:
5047              return escript.Data(arg,where)
5048      else:      else:
5049         return escript.Data(arg,where)         return escript.Data(arg,where)
5050    
# Line 5081  def L2(arg): Line 5163  def L2(arg):
5163      @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))      @note: L2(arg) is equivalent to sqrt(integrate(inner(arg,arg)))
5164      """      """
5165      return sqrt(integrate(inner(arg,arg)))      return sqrt(integrate(inner(arg,arg)))
5166    
5167    def getClosestValue(arg,origin=0):
5168        """
5169        returns the value in arg which is closest to origin
5170        
5171        @param arg: function
5172        @type arg: L{escript.Data} or L{Symbol}
5173        @param origin: reference value
5174        @type origin: C{float} or L{escript.Data}
5175        @return: value in arg closest to origin
5176        @rtype: L{numarray.NumArray} or L{Symbol}
5177        """
5178        return arg.getValueOfGlobalDataPoint(*(length(arg-origin).minGlobalDataPoint()))
5179    
5180    def normalize(arg,zerolength=0):
5181        """
5182        returns normalized version of arg (=arg/length(arg))
5183        
5184        @param arg: function
5185        @type arg: L{escript.Data} or L{Symbol}
5186        @param zerolength: realitive tolerance for arg == 0.
5187        @type zerolength: C{float}
5188        @return: normalized arg where arg is non zero and zero elsewhere
5189        @rtype: L{escript.Data} or L{Symbol}
5190        """
5191        l=length(arg)
5192        m=whereZero(l,zerolength*Lsup(l))
5193        mm=1-m
5194        return arg*(mm/(l*mm+m))
5195    
5196    def deviatoric(arg):
5197        """
5198        returns the deviatoric version of arg
5199        """
5200        return arg-(trace(arg)/trace(kronecker(arg.getDomain())))*kronecker(arg.getDomain())
5201    
5202    def diameter(domain):
5203        """
5204        returns the diameter of a domain
5205    
5206        @param domain: a domain
5207        @type domain: L{escript.Domain}
5208        @rtype: C{float}
5209        """
5210        x=domain.getX()
5211        out=0.
5212        for i in xrange(domain.getDim()):
5213           x_i=x[i]
5214           out=max(out,sup(x_i)-inf(x_i))
5215        return out
5216  #=============================  #=============================
5217  #  #
5218    
# Line 5090  def reorderComponents(arg,index): Line 5222  def reorderComponents(arg,index):
5222    
5223      """      """
5224      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.1042  
changed lines
  Added in v.2100

  ViewVC Help
Powered by ViewVC 1.1.26