/[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 2169 by caltinay, Wed Dec 17 03:08:58 2008 UTC revision 2288 by gross, Tue Feb 24 06:11:48 2009 UTC
# Line 182  def saveDX(filename,domain=None,**data): Line 182  def saveDX(filename,domain=None,**data):
182          raise ValueError,"saveDX: no domain detected."          raise ValueError,"saveDX: no domain detected."
183      domain.saveDX(filename,new_data)      domain.saveDX(filename,new_data)
184    
185    def saveESD(datasetName, dataDir=".", domain=None, **data):
186        """
187        Saves L{Data} objects to files and creates an I{escript dataset} (ESD) file
188        for convenient processing/visualisation.
189    
190        Example::
191    
192            tmp = Scalar(..)
193            v = Vector(..)
194            saveESD("solution", "data", temperature=tmp, velocity=v)
195    
196        tmp, v and the domain are saved in native format in the "data"
197        directory and the file "solution.esd" is created that refers to tmp by
198        the name "temperature" and to v by the name "velocity".
199    
200        @param datasetName: name of the dataset, used to name the ESD file
201        @type datasetName: C{str}
202        @param dataDir: optional directory where the data files should be saved
203        @type dataDir: C{str}
204        @param domain: domain of the L{Data} object(s). If not specified, the
205                       domain of the given L{Data} objects is used.
206        @type domain: L{escript.Domain}
207        @keyword <name>: writes the assigned value to the file using <name> as
208                         identifier
209        @type <name>: L{Data} object.
210        @note: The data objects have to be defined on the same domain. They may not
211               be in the same L{FunctionSpace} but one cannot expect that all
212               L{FunctionSpace}s can be mixed. Typically, data on the boundary and
213               data on the interior cannot be mixed.
214        """
215        new_data = {}
216        for n,d in data.items():
217              if not d.isEmpty():
218                fs = d.getFunctionSpace()
219                domain2 = fs.getDomain()
220                if fs == escript.Solution(domain2):
221                   new_data[n]=interpolate(d,escript.ContinuousFunction(domain2))
222                elif fs == escript.ReducedSolution(domain2):
223                   new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
224                else:
225                   new_data[n]=d
226                if domain==None: domain=domain2
227        if domain==None:
228            raise ValueError, "saveESD: no domain detected."
229    
230        if domain.onMasterProcessor() and not os.path.isdir(dataDir):
231            os.mkdir(dataDir)
232    
233        meshFile = os.path.join(dataDir, datasetName+"_mesh")
234        domain.dump(meshFile + ".nc")
235        outputString = ""
236        if domain.onMasterProcessor():
237            outputString += "#escript datafile V1.0\n"
238            # number of timesteps (currently only 1 is supported)
239            outputString += "T=1\n"
240            # name of the mesh file
241            outputString += "M=%s\n" % meshFile
242            # number of blocks (MPI size)
243            outputString += "N=%d\n" % domain.getMPISize()
244    
245        # now add the variables
246        for varName, d in new_data.items():
247            varFile = os.path.join(dataDir, datasetName+"_"+varName)
248            d.dump(varFile + ".nc")
249            if domain.onMasterProcessor():
250                outputString += "V=%s:%s\n" % (varFile, varName)
251    
252        if domain.onMasterProcessor():
253            esdfile = open(datasetName+".esd", "w")
254            esdfile.write(outputString)
255            esdfile.close()
256    
257  def kronecker(d=3):  def kronecker(d=3):
258     """     """
259     Returns the kronecker S{delta}-symbol.     Returns the kronecker S{delta}-symbol.
# Line 3368  def trace(arg,axis_offset=0): Line 3440  def trace(arg,axis_offset=0):
3440        if len(sh)<2:        if len(sh)<2:
3441          raise ValueError,"rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
3442        if axis_offset<0 or axis_offset>len(sh)-2:        if axis_offset<0 or axis_offset>len(sh)-2:
3443          raise ValueError,"axis_offset must be between 0 and %s"%len(sh)-2          raise ValueError,"axis_offset must be between 0 and %d"%(len(sh)-2)
3444        s1=1        s1=1
3445        for i in range(axis_offset): s1*=sh[i]        for i in range(axis_offset): s1*=sh[i]
3446        s2=1        s2=1
3447        for i in range(axis_offset+2,len(sh)): s2*=sh[i]        for i in range(axis_offset+2,len(sh)): s2*=sh[i]
3448        if not sh[axis_offset] == sh[axis_offset+1]:        if not sh[axis_offset] == sh[axis_offset+1]:
3449          raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)          raise ValueError,"dimensions of component %d and %d must match."%(axis_offset,axis_offset+1)
3450        arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))        arg_reshaped=numarray.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
3451        out=numarray.zeros([s1,s2],numarray.Float64)        out=numarray.zeros([s1,s2],numarray.Float64)
3452        for i1 in range(s1):        for i1 in range(s1):
# Line 3386  def trace(arg,axis_offset=0): Line 3458  def trace(arg,axis_offset=0):
3458        if arg.getRank()<2:        if arg.getRank()<2:
3459          raise ValueError,"rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
3460        if axis_offset<0 or axis_offset>arg.getRank()-2:        if axis_offset<0 or axis_offset>arg.getRank()-2:
3461          raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2          raise ValueError,"axis_offset must be between 0 and %d"%(arg.getRank()-2)
3462        s=list(arg.getShape())        s=list(arg.getShape())
3463        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
3464          raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)          raise ValueError,"dimensions of component %d and %d must match."%(axis_offset,axis_offset+1)
3465        return arg._trace(axis_offset)        return arg._trace(axis_offset)
3466     elif isinstance(arg,float):     elif isinstance(arg,float):
3467        raise TypeError,"illegal argument type float."        raise TypeError,"illegal argument type float."
# Line 3419  class Trace_Symbol(DependendSymbol): Line 3491  class Trace_Symbol(DependendSymbol):
3491        if arg.getRank()<2:        if arg.getRank()<2:
3492          raise ValueError,"rank of argument must be greater than 1"          raise ValueError,"rank of argument must be greater than 1"
3493        if axis_offset<0 or axis_offset>arg.getRank()-2:        if axis_offset<0 or axis_offset>arg.getRank()-2:
3494          raise ValueError,"axis_offset must be between 0 and %s"%arg.getRank()-2          raise ValueError,"axis_offset must be between 0 and %d"%(arg.getRank()-2)
3495        s=list(arg.getShape())        s=list(arg.getShape())
3496        if not s[axis_offset] == s[axis_offset+1]:        if not s[axis_offset] == s[axis_offset+1]:
3497          raise ValueError,"dimensions of component %s and %s must match."%(axis_offset.axis_offset+1)          raise ValueError,"dimensions of component %d and %d must match."%(axis_offset,axis_offset+1)
3498        super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())        super(Trace_Symbol,self).__init__(args=[arg,axis_offset],shape=tuple(s[0:axis_offset]+s[axis_offset+2:]),dim=arg.getDim())
3499    
3500     def getMyCode(self,argstrs,format="escript"):     def getMyCode(self,argstrs,format="escript"):
# Line 4578  def maximum(*args): Line 4650  def maximum(*args):
4650      out=None      out=None
4651      for a in args:      for a in args:
4652         if out==None:         if out==None:
4653            out=a            out=a*1.
4654         else:         else:
4655            diff=add(a,-out)            if isinstance(out,escript.Data) and isinstance(a,escript.Data):
4656            out=add(out,mult(wherePositive(diff),diff))           if out.getRank()==0 and a.getRank()>0:
4657            #We need to consider the case where we have scalars and higher
4658            #ranked objects mixed. If the scalar was first it will get
4659            #picked as the initial out and we have a problem,
4660            #so we swap the objects
4661            res=a.copy()    #Deep copy of a
4662            res.copyWithMask(out,wherePositive(out-a))
4663            out=res
4664             else:
4665                    out.copyWithMask(a,wherePositive(a-out))
4666              else:
4667                 diff=add(a,-out)
4668                 out=add(out,mult(wherePositive(diff),diff))
4669      return out      return out
4670    
4671  def minimum(*args):  def minimum(*args):
# Line 4599  def minimum(*args): Line 4683  def minimum(*args):
4683      out=None      out=None
4684      for a in args:      for a in args:
4685         if out==None:         if out==None:
4686            out=a            out=a*1.
4687         else:         else:
4688            diff=add(a,-out)            if isinstance(out,escript.Data) and isinstance(a,escript.Data):
4689            out=add(out,mult(whereNegative(diff),diff))           if out.getRank()==0 and a.getRank()>0:
4690            #We need to consider the case where we have scalars and higher
4691            #ranked objects mixed. If the scalar was first it will get
4692            #picked as the initial out and we have a problem,
4693            #so we swap the objects
4694            res=a.copy()    #Deep copy of a
4695            res.copyWithMask(out,whereNegative(out-a))
4696            out=res
4697             else:
4698            out.copyWithMask(a,whereNegative(a-out))
4699              else:
4700                 diff=add(a,-out)
4701                 out=add(out,mult(whereNegative(diff),diff))
4702      return out      return out
4703    
4704  def clip(arg,minval=None,maxval=None):  def clip(arg,minval=None,maxval=None):
# Line 5753  def diameter(domain): Line 5849  def diameter(domain):
5849      @type domain: L{escript.Domain}      @type domain: L{escript.Domain}
5850      @rtype: C{float}      @rtype: C{float}
5851      """      """
5852        out=0
5853        for v in boundingBox(domain): out+=(v[1]-v[0])**2
5854        return sqrt(out)
5855    
5856    def boundingBox(domain):
5857        """
5858        Returns the bounding box of a domain
5859    
5860        @param domain: a domain
5861        @type domain: L{escript.Domain}
5862        @return: bounding box of the domain
5863        @rtype: C{list} of pairs of C{float}
5864        """
5865      x=domain.getX()      x=domain.getX()
5866      out=0.      out=[]
5867      for i in xrange(domain.getDim()):      for i in xrange(domain.getDim()):
5868         x_i=x[i]         x_i=x[i]
5869         out=max(out,sup(x_i)-inf(x_i))         out.append((inf(x_i),sup(x_i)))
5870      return out      return out
5871    
5872    def longestEdge(domain):
5873        """
5874        Returns the length of the longest edge of the domain
5875    
5876        @param domain: a domain
5877        @type domain: L{escript.Domain}
5878        @return: longest edge of the domain parallel to the Cartisean axis
5879        @rtype: C{float}
5880        """
5881        return max([v[1]-v[0] for v in boundingBox(domain) ])
5882    
5883  #=============================  #=============================
5884  #  #
5885    

Legend:
Removed from v.2169  
changed lines
  Added in v.2288

  ViewVC Help
Powered by ViewVC 1.1.26