# Diff of /trunk/escript/py_src/util.py

revision 122 by jgs, Thu Jun 9 05:38:05 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:
9
10      - binary operations @ (@=+,-,*,/,**)::
11          (a@b)[:,*]=a[:]@b[:,*] if rank(a)<rank(b)
12          (a@b)[:]=a[:]@b[:] if rank(a)=rank(b)
13          (a@b)[*,:]=a[*,:]@b[:] if rank(a)>rank(b)
14      - implementation of outer::
15          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  #   escript constants (have to be consistent witj utilC.h  import os
#
UNKNOWN=-1
EPSILON=1.e-15
Pi=numarray.pi
# some solver options:
NO_REORDERING=0
MINIMUM_FILL_IN=1
NESTED_DISSECTION=2
# solver methods
DEFAULT_METHOD=0
DIRECT=1
CHOLEVSKY=2
PCG=3
CR=4
CGS=5
BICGSTAB=6
SSOR=7
ILU0=8
ILUT=9
JACOBI=10
GMRES=11
PRES20=12

METHOD_KEY="method"
SYMMETRY_KEY="symmetric"
TOLERANCE_KEY="tolerance"

# supported file formats:
VRML=1
PNG=2
JPEG=3
JPG=3
PS=4
OOGL=5
BMP=6
TIFF=7
OPENINVENTOR=8
RENDERMAN=9
PNM=10
24
25  #  #=========================================================
26  # wrapper for various functions: if the argument has attribute the function name  #   some little helpers
27  # as an argument it calls the correspong methods. Otherwise the coresponsing numarray  #=========================================================
28  # function is called.  def _testForZero(arg):
29  #     """
30  # functions involving the underlying Domain:     Returns True is arg is considered to be zero.
31  #     """
32  def grad(arg,where=None):     if isinstance(arg,int):
33      """        return not arg>0
34      @brief returns the spatial gradient of the Data object arg     elif isinstance(arg,float):
35          return not arg>0.
36       elif isinstance(arg,numarray.NumArray):
37          a=abs(arg)
38          while isinstance(a,numarray.NumArray): a=numarray.sometrue(a)
39          return not a>0
40       else:
41          return False
42
43      @param arg: Data object representing the function which gradient to be calculated.  #=========================================================
44      @param where: FunctionSpace in which the gradient will be. If None Function(dom) where dom is the  def saveVTK(filename,domain=None,**data):
domain of the Data object arg.
45      """      """
46      if isinstance(arg,escript.Data):      writes a L{Data} objects into a files using the the VTK XML file format.
if where==None:
else:
else:
return arg*0.
47
48  def integrate(arg,what=None):      Example:
"""
@brief return the integral if the function represented by Data object arg over its domain.
49
50      @param arg         tmp=Scalar(..)
51           v=Vector(..)
52           saveVTK("solution.xml",temperature=tmp,velovity=v)
53
54        tmp and v are written into "solution.dx" where tmp is named "temperature" and v is named "velovity"
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 not what==None:      if domain==None:
65         arg2=escript.Data(arg,what)         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:      else:
70         arg2=arg          domain.saveVTK(filename,data)
71      if arg2.getRank()==0:  #=========================================================
72          return arg2.integrate()[0]  def saveDX(filename,domain=None,**data):
else:
return arg2.integrate()

def interpolate(arg,where):
73      """      """
74      @brief interpolates the function represented by Data object arg into the FunctionSpace where.      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 arg      @param filename: file name of the output file
85      @param where      @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 isinstance(arg,escript.Data):      if domain==None:
93         return arg.interpolate(where)         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:      else:
98         return arg          domain.saveDX(filename,data)
99    #=========================================================
# functions returning Data objects:
100
101  def transpose(arg,axis=None):  def exp(arg):
102      """      """
103      @brief returns the transpose of the Data object arg.      Applies the exponential function to arg.
104
105      @param arg      @param arg: argument
106      """      """
107      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
108         # hack for transpose         return symbols.Exp_Symbol(arg)
109         r=arg.getRank()      elif hasattr(arg,"exp"):
110         if r!=2: raise ValueError,"Tranpose only avalaible for rank 2 objects"         return arg.exp()
s=arg.getShape()
out=escript.Data(0.,(s[1],s[0]),arg.getFunctionSpace())
for i in range(s[0]):
for j in range(s[1]):
out[j,i]=arg[i,j]
return out
# end hack for transpose
if axis==None: axis=arg.getRank()/2
return arg.transpose(axis)
111      else:      else:
112         if axis==None: axis=arg.rank/2         return numarray.exp(arg)
return numarray.transpose(arg,axis=axis)
113
114  def trace(arg):  def sqrt(arg):
115      """      """
116      @brief      Applies the squre root function to arg.
117
118      @param arg      @param arg: argument
119      """      """
120      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
121         # hack for trace         return symbols.Sqrt_Symbol(arg)
122         r=arg.getRank()      elif hasattr(arg,"sqrt"):
123         if r!=2: raise ValueError,"trace only avalaible for rank 2 objects"         return arg.sqrt()
s=arg.getShape()
out=escript.Scalar(0,arg.getFunctionSpace())
for i in range(min(s)):
out+=arg[i,i]
return out
# end hack for trace
return arg.trace()
124      else:      else:
125         return numarray.trace(arg)         return numarray.sqrt(arg)
126
127  def exp(arg):  def log(arg):
128      """      """
129      @brief      Applies the logarithmic function base 10 to arg.
130
131      @param arg      @param arg: argument
132      """      """
133      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
134         return arg.exp()         return symbols.Log_Symbol(arg)
135        elif hasattr(arg,"log"):
136           return arg.log()
137      else:      else:
138         return numarray.exp(arg)         return numarray.log10(arg)
139
140  def sqrt(arg):  def ln(arg):
141      """      """
142      @brief      Applies the natural logarithmic function to arg.
143
144      @param arg      @param arg: argument
145      """      """
146      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
147         return arg.sqrt()         return symbols.Ln_Symbol(arg)
148        elif hasattr(arg,"ln"):
149           return arg.ln()
150      else:      else:
151         return numarray.sqrt(arg)         return numarray.log(arg)
152
153  def sin(arg):  def sin(arg):
154      """      """
155      @brief      Applies the sin function to arg.
156
157      @param arg      @param arg: argument
158      """      """
159      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
160           return symbols.Sin_Symbol(arg)
161        elif hasattr(arg,"sin"):
162         return arg.sin()         return arg.sin()
163      else:      else:
164         return numarray.sin(arg)         return numarray.sin(arg)
165
166    def cos(arg):
167        """
168        Applies the cos function to arg.
169
170        @param arg: argument
171        """
172        if isinstance(arg,symbols.Symbol):
173           return symbols.Cos_Symbol(arg)
174        elif hasattr(arg,"cos"):
175           return arg.cos()
176        else:
177           return numarray.cos(arg)
178
179  def tan(arg):  def tan(arg):
180      """      """
181      @brief      Applies the tan function to arg.
182
183      @param arg      @param arg: argument
184      """      """
185      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
186           return symbols.Tan_Symbol(arg)
187        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  def cos(arg):  def asin(arg):
193      """      """
194      @brief      Applies the asin function to arg.
195
196      @param arg      @param arg: argument
197      """      """
198      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
199         return arg.cos()         return symbols.Asin_Symbol(arg)
200        elif hasattr(arg,"asin"):
201           return arg.asin()
202      else:      else:
203         return numarray.cos(arg)         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):
310        """
311        Applies the sign function to arg.
312
313        @param arg: argument
314        """
315        if isinstance(arg,symbols.Symbol):
316           return symbols.Sign_Symbol(arg)
317        elif hasattr(arg,"sign"):
318           return arg.sign()
319        else:
320           return numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))- \
321                  numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))
322
323  def maxval(arg):  def maxval(arg):
324      """      """
325      @brief      Returns the maximum value of argument arg.
326
327      @param arg      @param arg: argument
328      """      """
329      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
330           return symbols.Max_Symbol(arg)
331        elif hasattr(arg,"maxval"):
332         return arg.maxval()         return arg.maxval()
333      elif isinstance(arg,float) or isinstance(arg,int):      elif hasattr(arg,"max"):
return arg
else:
334         return arg.max()         return arg.max()
335        else:
336           return arg
337
338  def minval(arg):  def minval(arg):
339      """      """
340      @brief      Returns the minimum value of argument arg.
341
342      @param arg      @param arg: argument
343      """      """
344      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
345           return symbols.Min_Symbol(arg)
346        elif hasattr(arg,"maxval"):
347         return arg.minval()         return arg.minval()
348      elif isinstance(arg,float) or isinstance(arg,int):      elif hasattr(arg,"min"):
349           return arg.min()
350        else:
351         return arg         return arg
352
353    def wherePositive(arg):
354        """
355        Returns the positive values of argument arg.
356
357        @param arg: argument
358        """
359        if _testForZero(arg):
360          return 0
361        elif isinstance(arg,symbols.Symbol):
362           return symbols.WherePositive_Symbol(arg)
363        elif hasattr(arg,"wherePositive"):
364           return arg.minval()
365        elif hasattr(arg,"wherePositive"):
366           numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))
367      else:      else:
368         return arg.min()         if arg>0:
369              return 1.
370           else:
371              return 0.
372
373    def whereNegative(arg):
374        """
375        Returns the negative values of argument arg.
376
377        @param arg: argument
378        """
379        if _testForZero(arg):
380          return 0
381        elif isinstance(arg,symbols.Symbol):
382           return symbols.WhereNegative_Symbol(arg)
383        elif hasattr(arg,"whereNegative"):
384           return arg.whereNegative()
385        elif hasattr(arg,"shape"):
386           numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))
387        else:
388           if arg<0:
389              return 1.
390           else:
391              return 0.
392
393    def maximum(arg0,arg1):
394        """
395        Return arg1 where arg1 is bigger then arg0 otherwise arg0 is returned.
396        """
397        m=whereNegative(arg0-arg1)
398        return m*arg1+(1.-m)*arg0
399
400    def minimum(arg0,arg1):
401        """
402        Return arg0 where arg1 is bigger then arg0 otherwise arg1 is returned.
403        """
404        m=whereNegative(arg0-arg1)
405        return m*arg0+(1.-m)*arg1
406
407    def outer(arg0,arg1):
408       if _testForZero(arg0) or _testForZero(arg1):
409          return 0
410       else:
411          if isinstance(arg0,symbols.Symbol) or isinstance(arg1,symbols.Symbol):
412            return symbols.Outer_Symbol(arg0,arg1)
413          elif _identifyShape(arg0)==() or _identifyShape(arg1)==():
414            return arg0*arg1
415          elif isinstance(arg0,numarray.NumArray) and isinstance(arg1,numarray.NumArray):
416            return numarray.outer(arg0,arg1)
417          else:
418            if arg0.getRank()==1 and arg1.getRank()==1:
419              out=escript.Data(0,(arg0.getShape()[0],arg1.getShape()[0]),arg1.getFunctionSpace())
420              for i in range(arg0.getShape()[0]):
421                for j in range(arg1.getShape()[0]):
422                    out[i,j]=arg0[i]*arg1[j]
423              return out
424            else:
425              raise ValueError,"outer is not fully implemented yet."
426
427    def interpolate(arg,where):
428        """
429        Interpolates the function into the FunctionSpace where.
430
431        @param arg:    interpolant
432        @param where:  FunctionSpace to interpolate to
433        """
434        if _testForZero(arg):
435          return 0
436        elif isinstance(arg,symbols.Symbol):
437           return symbols.Interpolated_Symbol(arg,where)
438        else:
439           return escript.Data(arg,where)
440
441    def div(arg,where=None):
442        """
443        Returns the divergence of arg at where.
444
445        @param arg:   Data object representing the function which gradient to
446                      be calculated.
447        @param where: FunctionSpace in which the gradient will be calculated.
448                      If not present or C{None} an appropriate default is used.
449        """
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
464        """
465        Returns the spatial gradient of arg at where.
466
467        @param arg:   Data object representing the function which gradient
468                      to be calculated.
469        @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):
473          return 0
474        elif isinstance(arg,symbols.Symbol):
477           if where==None:
479           else:
481        else:
482           return arg*0.
483
484    def integrate(arg,where=None):
485        """
486        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.
490        @param where: FunctionSpace in which the integral is calculated.
491                      If not present or C{None} an appropriate default is used.
492        """
493        if _testForZero(arg):
494          return 0
495        elif isinstance(arg,symbols.Symbol):
496           return symbols.Integral_Symbol(arg,where)
497        else:
498           if not where==None: arg=escript.Data(arg,where)
499           if arg.getRank()==0:
500             return arg.integrate()[0]
501           else:
502             return arg.integrate()
503
504    #=============================
505    #
506    # wrapper for various functions: if the argument has attribute the function name
507    # as an argument it calls the corresponding methods. Otherwise the corresponding
508    # numarray function is called.
509
510    # functions involving the underlying Domain:
511
512
513    # functions returning Data objects:
514
515    def transpose(arg,axis=None):
516        """
517        Returns the transpose of the Data object arg.
518
519        @param arg:
520        """
521        if axis==None:
522           r=0
523           if hasattr(arg,"getRank"): r=arg.getRank()
524           if hasattr(arg,"rank"): r=arg.rank
525           axis=r/2
526        if isinstance(arg,symbols.Symbol):
527           return symbols.Transpose_Symbol(arg,axis=r)
528        if isinstance(arg,escript.Data):
529           # hack for transpose
530           r=arg.getRank()
531           if r!=2: raise ValueError,"Tranpose only avalaible for rank 2 objects"
532           s=arg.getShape()
533           out=escript.Data(0.,(s[1],s[0]),arg.getFunctionSpace())
534           for i in range(s[0]):
535              for j in range(s[1]):
536                 out[j,i]=arg[i,j]
537           return out
538           # end hack for transpose
539           return arg.transpose(axis)
540        else:
541           return numarray.transpose(arg,axis=axis)
542
543    def trace(arg,axis0=0,axis1=1):
544        """
545        Return
546
547        @param arg:
548        """
549        if isinstance(arg,symbols.Symbol):
550           s=list(arg.getShape())
551           s=tuple(s[0:axis0]+s[axis0+1:axis1]+s[axis1+1:])
552           return symbols.Trace_Symbol(arg,axis0=axis0,axis1=axis1)
553        elif isinstance(arg,escript.Data):
554           # hack for trace
555           s=arg.getShape()
556           if s[axis0]!=s[axis1]:
557               raise ValueError,"illegal axis in trace"
558           out=escript.Scalar(0.,arg.getFunctionSpace())
559           for i in range(s[axis0]):
560              out+=arg[i,i]
561           return out
562           # end hack for trace
563        else:
564           return numarray.trace(arg,axis0=axis0,axis1=axis1)
565
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()
573         if arg.getRank()==0:         if arg.getRank()==0:
574            return abs(arg)            return abs(arg)
575         elif arg.getRank()==1:         elif arg.getRank()==1:
576            sum=escript.Scalar(0,arg.getFunctionSpace())            out=escript.Scalar(0,arg.getFunctionSpace())
577            for i in range(arg.getShape()[0]):            for i in range(arg.getShape()[0]):
578               sum+=arg[i]**2               out+=arg[i]**2
579            return sqrt(sum)            return sqrt(out)
580         elif arg.getRank()==2:         elif arg.getRank()==2:
581            sum=escript.Scalar(0,arg.getFunctionSpace())            out=escript.Scalar(0,arg.getFunctionSpace())
582            for i in range(arg.getShape()[0]):            for i in range(arg.getShape()[0]):
583               for j in range(arg.getShape()[1]):               for j in range(arg.getShape()[1]):
584                  sum+=arg[i,j]**2                  out+=arg[i,j]**2
585            return sqrt(sum)            return sqrt(out)
586         elif arg.getRank()==3:         elif arg.getRank()==3:
587            sum=escript.Scalar(0,arg.getFunctionSpace())            out=escript.Scalar(0,arg.getFunctionSpace())
588            for i in range(arg.getShape()[0]):            for i in range(arg.getShape()[0]):
589               for j in range(arg.getShape()[1]):               for j in range(arg.getShape()[1]):
590                  for k in range(arg.getShape()[2]):                  for k in range(arg.getShape()[2]):
591                     sum+=arg[i,j,k]**2                     out+=arg[i,j,k]**2
592            return sqrt(sum)            return sqrt(out)
593         elif arg.getRank()==4:         elif arg.getRank()==4:
594            sum=escript.Scalar(0,arg.getFunctionSpace())            out=escript.Scalar(0,arg.getFunctionSpace())
595            for i in range(arg.getShape()[0]):            for i in range(arg.getShape()[0]):
596               for j in range(arg.getShape()[1]):               for j in range(arg.getShape()[1]):
597                  for k in range(arg.getShape()[2]):                  for k in range(arg.getShape()[2]):
598                     for l in range(arg.getShape()[3]):                     for l in range(arg.getShape()[3]):
599                        sum+=arg[i,j,k,l]**2                        out+=arg[i,j,k,l]**2
600            return sqrt(sum)            return sqrt(out)
601         else:         else:
602            raise SystemError,"length is not been implemented yet"            raise SystemError,"length is not been fully implemented yet"
603         # return arg.length()            # return arg.length()
604        elif isinstance(arg,float):
605           return abs(arg)
606      else:      else:
607         return sqrt((arg**2).sum())         return sqrt((arg**2).sum())
608
609  def deviator(arg):  def deviator(arg):
610      """      """
611      @brief      @param arg:

@param arg1
612      """      """
613      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
614          shape=arg.getShape()          shape=arg.getShape()
# Line 284  def deviator(arg): Line 620  def deviator(arg):
620            raise ValueError,"Deviator requires a square matrix"            raise ValueError,"Deviator requires a square matrix"
621      return arg-1./(shape[0]*1.)*trace(arg)*kronecker(shape[0])      return arg-1./(shape[0]*1.)*trace(arg)*kronecker(shape[0])
622
623  def inner(arg1,arg2):  def inner(arg0,arg1):
624      """      """
625      @brief      @param arg0:
626        @param arg1:
@param arg1, arg2
627      """      """
628      sum=escript.Scalar(0,arg1.getFunctionSpace())      if isinstance(arg0,escript.Data):
629           arg=arg0
630        else:
631           arg=arg1
632
633        out=escript.Scalar(0,arg.getFunctionSpace())
634      if arg.getRank()==0:      if arg.getRank()==0:
635            return arg1*arg2            return arg0*arg1
636      elif arg.getRank()==1:      elif arg.getRank()==1:
637           sum=escript.Scalar(0,arg.getFunctionSpace())           out=escript.Scalar(0,arg.getFunctionSpace())
638           for i in range(arg.getShape()[0]):           for i in range(arg.getShape()[0]):
639              sum+=arg1[i]*arg2[i]              out+=arg0[i]*arg1[i]
640      elif arg.getRank()==2:      elif arg.getRank()==2:
641          sum=escript.Scalar(0,arg.getFunctionSpace())          out=escript.Scalar(0,arg.getFunctionSpace())
642          for i in range(arg.getShape()[0]):          for i in range(arg.getShape()[0]):
643             for j in range(arg.getShape()[1]):             for j in range(arg.getShape()[1]):
644                sum+=arg1[i,j]*arg2[i,j]                out+=arg0[i,j]*arg1[i,j]
645      elif arg.getRank()==3:      elif arg.getRank()==3:
646          sum=escript.Scalar(0,arg.getFunctionSpace())          out=escript.Scalar(0,arg.getFunctionSpace())
647          for i in range(arg.getShape()[0]):          for i in range(arg.getShape()[0]):
648              for j in range(arg.getShape()[1]):              for j in range(arg.getShape()[1]):
649                 for k in range(arg.getShape()[2]):                 for k in range(arg.getShape()[2]):
650                    sum+=arg1[i,j,k]*arg2[i,j,k]                    out+=arg0[i,j,k]*arg1[i,j,k]
651      elif arg.getRank()==4:      elif arg.getRank()==4:
652          sum=escript.Scalar(0,arg.getFunctionSpace())          out=escript.Scalar(0,arg.getFunctionSpace())
653          for i in range(arg.getShape()[0]):          for i in range(arg.getShape()[0]):
654             for j in range(arg.getShape()[1]):             for j in range(arg.getShape()[1]):
655                for k in range(arg.getShape()[2]):                for k in range(arg.getShape()[2]):
656                   for l in range(arg.getShape()[3]):                   for l in range(arg.getShape()[3]):
657                      sum+=arg1[i,j,k,l]*arg2[i,j,k,l]                      out+=arg0[i,j,k,l]*arg1[i,j,k,l]
658      else:      else:
659            raise SystemError,"inner is not been implemented yet"            raise SystemError,"inner is not been implemented yet"
660      return sum      return out
661
662  def sign(arg):  def tensormult(arg0,arg1):
663      """      # check LinearPDE!!!!
664      @brief      raise SystemError,"tensormult is not implemented yet!"
665
666      @param arg  def matrixmult(arg0,arg1):
667      """
668      if isinstance(arg,escript.Data):      if isinstance(arg1,numarray.NumArray) and isinstance(arg0,numarray.NumArray):
669         return arg.sign()          numarray.matrixmult(arg0,arg1)
670      else:      else:
671         return numarray.greater(arg,numarray.zeros(arg.shape))-numarray.less(arg,numarray.zeros(arg.shape))        # escript.matmult(arg0,arg1)
672          if isinstance(arg1,escript.Data) and not isinstance(arg0,escript.Data):
673            arg0=escript.Data(arg0,arg1.getFunctionSpace())
674          elif isinstance(arg0,escript.Data) and not isinstance(arg1,escript.Data):
675            arg1=escript.Data(arg1,arg0.getFunctionSpace())
676          if arg0.getRank()==2 and arg1.getRank()==1:
677              out=escript.Data(0,(arg0.getShape()[0],),arg0.getFunctionSpace())
678              for i in range(arg0.getShape()[0]):
679                 for j in range(arg0.getShape()[1]):
680                   # uses Data object slicing, plus Data * and += operators
681                   out[i]+=arg0[i,j]*arg1[j]
682              return out
683          elif arg0.getRank()==1 and arg1.getRank()==1:
684              return inner(arg0,arg1)
685          else:
686              raise SystemError,"matrixmult is not fully implemented yet!"
687
688    #=========================================================
689  # reduction operations:  # reduction operations:
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 355  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 368  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 381  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 Linf(arg):  def dot(arg0,arg1):
744      """      """
745      @brief      @param arg0:
746        @param arg1:
@param arg
"""
if isinstance(arg,escript.Data):
return arg.Linf()
elif isinstance(arg,float) or isinstance(arg,int):
return abs(arg)
else:
return min(numarray.abs(arg))

def dot(arg1,arg2):
"""
@brief

@param arg
747      """      """
748      if isinstance(arg1,escript.Data):      if isinstance(arg0,escript.Data):
749         return arg1.dot(arg2)         return arg0.dot(arg1)
750      elif isinstance(arg1,escript.Data):      elif isinstance(arg1,escript.Data):
751         return arg2.dot(arg1)         return arg1.dot(arg0)
752      else:      else:
753         return numarray.dot(arg1,arg2)         return numarray.dot(arg0,arg1)
754
755  def kronecker(d):  def kronecker(d):
756     if hasattr(d,"getDim"):     if hasattr(d,"getDim"):
757        return numarray.identity(d.getDim())        return numarray.identity(d.getDim())*1.
758     else:     else:
759        return numarray.identity(d)        return numarray.identity(d)*1.
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,))     e = numarray.zeros((d,),numarray.Float)
769     e[i] = 1.0     e[i] = 1.0
770     return e     return e

#
# \$Log\$
# 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.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.122 changed lines Added in v.153