/[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 100 by jgs, Wed Dec 15 03:48:48 2004 UTC revision 122 by jgs, Thu Jun 9 05:38:05 2005 UTC
# Line 7  Line 7 
7  """  """
8    
9  import numarray  import numarray
10    import escript
11  #  #
12  #   escript constants:  #   escript constants (have to be consistent witj utilC.h
13  #  #
 FALSE=0  
 TRUE=1  
14  UNKNOWN=-1  UNKNOWN=-1
15  EPSILON=1.e-15  EPSILON=1.e-15
16  Pi=3.1415926535897931  Pi=numarray.pi
 # matrix types  
 CSC=0    
 CSR=1  
 LUMPED=10  
17  # some solver options:  # some solver options:
18  NO_REORDERING=0  NO_REORDERING=0
19  MINIMUM_FILL_IN=1  MINIMUM_FILL_IN=1
20  NESTED_DISSECTION=2  NESTED_DISSECTION=2
21    # solver methods
22  DEFAULT_METHOD=0  DEFAULT_METHOD=0
23  PCG=1  DIRECT=1
24  CR=2  CHOLEVSKY=2
25  CGS=3  PCG=3
26  BICGSTAB=4  CR=4
27  SSOR=5  CGS=5
28  ILU0=6  BICGSTAB=6
29  ILUT=7  SSOR=7
30  JACOBI=8  ILU0=8
31    ILUT=9
32    JACOBI=10
33    GMRES=11
34    PRES20=12
35    
36    METHOD_KEY="method"
37    SYMMETRY_KEY="symmetric"
38    TOLERANCE_KEY="tolerance"
39    
40  # supported file formats:  # supported file formats:
41  VRML=1  VRML=1
42  PNG=2  PNG=2
# Line 44  TIFF=7 Line 49  TIFF=7
49  OPENINVENTOR=8  OPENINVENTOR=8
50  RENDERMAN=9  RENDERMAN=9
51  PNM=10  PNM=10
52    
53  #  #
54  # wrapper for various functions: if the argument has attribute the function name  # wrapper for various functions: if the argument has attribute the function name
55  # as an argument it calls the correspong methods. Otherwise the coresponsing numarray  # as an argument it calls the correspong methods. Otherwise the coresponsing numarray
56  # function is called.  # function is called.
57  #  #
58  def L2(arg):  # functions involving the underlying Domain:
59      """  #
     @brief  
   
     @param arg  
     """  
     return arg.L2()  
   
60  def grad(arg,where=None):  def grad(arg,where=None):
61      """      """
62      @brief      @brief returns the spatial gradient of the Data object arg
63    
64      @param arg      @param arg: Data object representing the function which gradient to be calculated.
65      @param where      @param where: FunctionSpace in which the gradient will be. If None Function(dom) where dom is the
66      """                    domain of the Data object arg.
67      if where==None:      """
68         return arg.grad()      if isinstance(arg,escript.Data):
69           if where==None:
70              return arg.grad()
71           else:
72              return arg.grad(where)
73      else:      else:
74         return arg.grad(where)         return arg*0.
75    
76  def integrate(arg):  def integrate(arg,what=None):
77      """      """
78      @brief      @brief return the integral if the function represented by Data object arg over its domain.
79    
80      @param arg      @param arg
81      """      """
82      return arg.integrate()      if not what==None:
83           arg2=escript.Data(arg,what)
84        else:
85           arg2=arg
86        if arg2.getRank()==0:
87            return arg2.integrate()[0]
88        else:
89            return arg2.integrate()
90    
91  def interpolate(arg,where):  def interpolate(arg,where):
92      """      """
93      @brief      @brief interpolates the function represented by Data object arg into the FunctionSpace where.
94    
95      @param arg      @param arg
96      @param where      @param where
97      """      """
98      return arg.interpolate(where)      if isinstance(arg,escript.Data):
99           return arg.interpolate(where)
100        else:
101           return arg
102    
103    # functions returning Data objects:
104    
105  def transpose(arg):  def transpose(arg,axis=None):
106      """      """
107      @brief      @brief returns the transpose of the Data object arg.
108    
109      @param arg      @param arg
110      """      """
111      if hasattr(arg,"transpose"):      if isinstance(arg,escript.Data):
112         return arg.transpose()         # hack for transpose
113           r=arg.getRank()
114           if r!=2: raise ValueError,"Tranpose only avalaible for rank 2 objects"
115           s=arg.getShape()
116           out=escript.Data(0.,(s[1],s[0]),arg.getFunctionSpace())
117           for i in range(s[0]):
118              for j in range(s[1]):
119                 out[j,i]=arg[i,j]
120           return out
121           # end hack for transpose
122           if axis==None: axis=arg.getRank()/2
123           return arg.transpose(axis)
124      else:      else:
125         return numarray.transpose(arg,axis=None)         if axis==None: axis=arg.rank/2
126           return numarray.transpose(arg,axis=axis)
127    
128  def trace(arg):  def trace(arg):
129      """      """
# Line 103  def trace(arg): Line 131  def trace(arg):
131    
132      @param arg      @param arg
133      """      """
134      if hasattr(arg,"trace"):      if isinstance(arg,escript.Data):
135           # hack for trace
136           r=arg.getRank()
137           if r!=2: raise ValueError,"trace only avalaible for rank 2 objects"
138           s=arg.getShape()
139           out=escript.Scalar(0,arg.getFunctionSpace())
140           for i in range(min(s)):
141                 out+=arg[i,i]
142           return out
143           # end hack for trace
144         return arg.trace()         return arg.trace()
145      else:      else:
146         return numarray.trace(arg,k=0)         return numarray.trace(arg)
147    
148  def exp(arg):  def exp(arg):
149      """      """
# Line 114  def exp(arg): Line 151  def exp(arg):
151    
152      @param arg      @param arg
153      """      """
154      if hasattr(arg,"exp"):      if isinstance(arg,escript.Data):
155         return arg.exp()         return arg.exp()
156      else:      else:
157         return numarray.exp(arg)         return numarray.exp(arg)
# Line 125  def sqrt(arg): Line 162  def sqrt(arg):
162    
163      @param arg      @param arg
164      """      """
165      if hasattr(arg,"sqrt"):      if isinstance(arg,escript.Data):
166         return arg.sqrt()         return arg.sqrt()
167      else:      else:
168         return numarray.sqrt(arg)         return numarray.sqrt(arg)
# Line 136  def sin(arg): Line 173  def sin(arg):
173    
174      @param arg      @param arg
175      """      """
176      if hasattr(arg,"sin"):      if isinstance(arg,escript.Data):
177         return arg.sin()         return arg.sin()
178      else:      else:
179         return numarray.sin(arg)         return numarray.sin(arg)
180    
181    def tan(arg):
182        """
183        @brief
184    
185        @param arg
186        """
187        if isinstance(arg,escript.Data):
188           return arg.tan()
189        else:
190           return numarray.tan(arg)
191    
192  def cos(arg):  def cos(arg):
193      """      """
194      @brief      @brief
195    
196      @param arg      @param arg
197      """      """
198      if hasattr(arg,"cos"):      if isinstance(arg,escript.Data):
199         return arg.cos()         return arg.cos()
200      else:      else:
201         return numarray.cos(arg)         return numarray.cos(arg)
# Line 158  def maxval(arg): Line 206  def maxval(arg):
206    
207      @param arg      @param arg
208      """      """
209      return arg.maxval()      if isinstance(arg,escript.Data):
210           return arg.maxval()
211        elif isinstance(arg,float) or isinstance(arg,int):
212           return arg
213        else:
214           return arg.max()
215    
216  def minval(arg):  def minval(arg):
217      """      """
# Line 166  def minval(arg): Line 219  def minval(arg):
219    
220      @param arg      @param arg
221      """      """
222      return arg.minval()      if isinstance(arg,escript.Data):
223           return arg.minval()
224        elif isinstance(arg,float) or isinstance(arg,int):
225           return arg
226        else:
227           return arg.min()
228    
229    def length(arg):
230        """
231        @brief
232    
233        @param arg
234        """
235        if isinstance(arg,escript.Data):
236           if arg.isEmpty(): return escript.Data()
237           if arg.getRank()==0:
238              return abs(arg)
239           elif arg.getRank()==1:
240              sum=escript.Scalar(0,arg.getFunctionSpace())
241              for i in range(arg.getShape()[0]):
242                 sum+=arg[i]**2
243              return sqrt(sum)
244           elif arg.getRank()==2:
245              sum=escript.Scalar(0,arg.getFunctionSpace())
246              for i in range(arg.getShape()[0]):
247                 for j in range(arg.getShape()[1]):
248                    sum+=arg[i,j]**2
249              return sqrt(sum)
250           elif arg.getRank()==3:
251              sum=escript.Scalar(0,arg.getFunctionSpace())
252              for i in range(arg.getShape()[0]):
253                 for j in range(arg.getShape()[1]):
254                    for k in range(arg.getShape()[2]):
255                       sum+=arg[i,j,k]**2
256              return sqrt(sum)
257           elif arg.getRank()==4:
258              sum=escript.Scalar(0,arg.getFunctionSpace())
259              for i in range(arg.getShape()[0]):
260                 for j in range(arg.getShape()[1]):
261                    for k in range(arg.getShape()[2]):
262                       for l in range(arg.getShape()[3]):
263                          sum+=arg[i,j,k,l]**2
264              return sqrt(sum)
265           else:
266              raise SystemError,"length is not been implemented yet"
267           # return arg.length()
268        else:
269           return sqrt((arg**2).sum())
270    
271    def deviator(arg):
272        """
273        @brief
274    
275        @param arg1
276        """
277        if isinstance(arg,escript.Data):
278            shape=arg.getShape()
279        else:
280            shape=arg.shape
281        if len(shape)!=2:
282              raise ValueError,"Deviator requires rank 2 object"
283        if shape[0]!=shape[1]:
284              raise ValueError,"Deviator requires a square matrix"
285        return arg-1./(shape[0]*1.)*trace(arg)*kronecker(shape[0])
286    
287    def inner(arg1,arg2):
288        """
289        @brief
290    
291        @param arg1, arg2
292        """
293        sum=escript.Scalar(0,arg1.getFunctionSpace())
294        if arg.getRank()==0:
295              return arg1*arg2
296        elif arg.getRank()==1:
297             sum=escript.Scalar(0,arg.getFunctionSpace())
298             for i in range(arg.getShape()[0]):
299                sum+=arg1[i]*arg2[i]
300        elif arg.getRank()==2:
301            sum=escript.Scalar(0,arg.getFunctionSpace())
302            for i in range(arg.getShape()[0]):
303               for j in range(arg.getShape()[1]):
304                  sum+=arg1[i,j]*arg2[i,j]
305        elif arg.getRank()==3:
306            sum=escript.Scalar(0,arg.getFunctionSpace())
307            for i in range(arg.getShape()[0]):
308                for j in range(arg.getShape()[1]):
309                   for k in range(arg.getShape()[2]):
310                      sum+=arg1[i,j,k]*arg2[i,j,k]
311        elif arg.getRank()==4:
312            sum=escript.Scalar(0,arg.getFunctionSpace())
313            for i in range(arg.getShape()[0]):
314               for j in range(arg.getShape()[1]):
315                  for k in range(arg.getShape()[2]):
316                     for l in range(arg.getShape()[3]):
317                        sum+=arg1[i,j,k,l]*arg2[i,j,k,l]
318        else:
319              raise SystemError,"inner is not been implemented yet"
320        return sum
321    
322    def sign(arg):
323        """
324        @brief
325    
326        @param arg
327        """
328        if isinstance(arg,escript.Data):
329           return arg.sign()
330        else:
331           return numarray.greater(arg,numarray.zeros(arg.shape))-numarray.less(arg,numarray.zeros(arg.shape))
332    
333    # reduction operations:
334    
335    def sum(arg):
336        """
337        @brief
338    
339        @param arg
340        """
341        return arg.sum()
342    
343  def sup(arg):  def sup(arg):
344      """      """
# Line 174  def sup(arg): Line 346  def sup(arg):
346    
347      @param arg      @param arg
348      """      """
349      return arg.sup()      if isinstance(arg,escript.Data):
350           return arg.sup()
351        elif isinstance(arg,float) or isinstance(arg,int):
352           return arg
353        else:
354           return arg.max()
355    
356  def inf(arg):  def inf(arg):
357      """      """
# Line 182  def inf(arg): Line 359  def inf(arg):
359    
360      @param arg      @param arg
361      """      """
362      return arg.inf()      if isinstance(arg,escript.Data):
363           return arg.inf()
364        elif isinstance(arg,float) or isinstance(arg,int):
365           return arg
366        else:
367           return arg.min()
368    
369    def L2(arg):
370        """
371        @brief returns the L2-norm of the
372    
373        @param arg
374        """
375        if isinstance(arg,escript.Data):
376           return arg.L2()
377        elif isinstance(arg,float) or isinstance(arg,int):
378           return abs(arg)
379        else:
380           return numarry.sqrt(dot(arg,arg))
381    
382  def Lsup(arg):  def Lsup(arg):
383      """      """
# Line 190  def Lsup(arg): Line 385  def Lsup(arg):
385    
386      @param arg      @param arg
387      """      """
388      return arg.Lsup()      if isinstance(arg,escript.Data):
389           return arg.Lsup()
390        elif isinstance(arg,float) or isinstance(arg,int):
391           return abs(arg)
392        else:
393           return max(numarray.abs(arg))
394    
395  def length(arg):  def Linf(arg):
396      """      """
397      @brief      @brief
398    
399      @param arg      @param arg
400      """      """
401      return arg.length()      if isinstance(arg,escript.Data):
402           return arg.Linf()
403        elif isinstance(arg,float) or isinstance(arg,int):
404           return abs(arg)
405        else:
406           return min(numarray.abs(arg))
407    
408  def sign(arg):  def dot(arg1,arg2):
409      """      """
410      @brief      @brief
411    
412      @param arg      @param arg
413      """      """
414      return arg.sign()      if isinstance(arg1,escript.Data):
415           return arg1.dot(arg2)
416        elif isinstance(arg1,escript.Data):
417           return arg2.dot(arg1)
418        else:
419           return numarray.dot(arg1,arg2)
420    
421    def kronecker(d):
422       if hasattr(d,"getDim"):
423          return numarray.identity(d.getDim())
424       else:
425          return numarray.identity(d)
426    
427    def unit(i,d):
428       """
429       @brief return a unit vector of dimension d with nonzero index i
430       @param d dimension
431       @param i index
432       """
433       e = numarray.zeros((d,))
434       e[i] = 1.0
435       return e
436    
437  #  #
438  # $Log$  # $Log$
439  # Revision 1.4  2004/12/15 03:48:42  jgs  # Revision 1.10  2005/06/09 05:37:59  jgs
440  # *** empty log message ***  # Merge of development branch back to main trunk on 2005-06-09
441    #
442    # Revision 1.2.2.14  2005/05/20 04:05:23  gross
443    # some work on a darcy flow started
444    #
445    # Revision 1.2.2.13  2005/03/16 05:17:58  matt
446    # Implemented unit(idx, dim) to create cartesian unit basis vectors to
447    # complement kronecker(dim) function.
448    #
449    # Revision 1.2.2.12  2005/03/10 08:14:37  matt
450    # Added non-member Linf utility function to complement Data::Linf().
451    #
452    # Revision 1.2.2.11  2005/02/17 05:53:25  gross
453    # some bug in saveDX fixed: in fact the bug was in
454    # DataC/getDataPointShape
455    #
456    # Revision 1.2.2.10  2005/01/11 04:59:36  gross
457    # automatic interpolation in integrate switched off
458    #
459    # Revision 1.2.2.9  2005/01/11 03:38:13  gross
460    # 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.
461    #
462    # Revision 1.2.2.8  2005/01/05 04:21:41  gross
463    # FunctionSpace checking/matchig in slicing added
464    #
465    # Revision 1.2.2.7  2004/12/29 05:29:59  gross
466    # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()
467    #
468    # Revision 1.2.2.6  2004/12/24 06:05:41  gross
469    # some changes in linearPDEs to add AdevectivePDE
470    #
471    # Revision 1.2.2.5  2004/12/17 00:06:53  gross
472    # mk sets ESYS_ROOT is undefined
473    #
474    # Revision 1.2.2.4  2004/12/07 03:19:51  gross
475    # options for GMRES and PRES20 added
476    #
477    # Revision 1.2.2.3  2004/12/06 04:55:18  gross
478    # function wraper extended
479    #
480    # Revision 1.2.2.2  2004/11/22 05:44:07  gross
481    # a few more unitary functions have been added but not implemented in Data yet
482    #
483    # Revision 1.2.2.1  2004/11/12 06:58:15  gross
484    # 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
485  #  #
486  # Revision 1.2  2004/10/27 00:23:36  jgs  # Revision 1.2  2004/10/27 00:23:36  jgs
487  # fixed minor syntax error  # fixed minor syntax error

Legend:
Removed from v.100  
changed lines
  Added in v.122

  ViewVC Help
Powered by ViewVC 1.1.26