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

revision 108 by jgs, Thu Jan 27 06:21:59 2005 UTC revision 150 by jgs, Thu Sep 15 03:44:45 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
#
# wrapper for various functions: if the argument has attribute the function name
# as an argument it calls the correspong methods. Otherwise the coresponsing numarray
# function is called.
#
# functions involving the underlying Domain:
#
def grad(arg,where=None):
"""
@brief returns the spatial gradient of the Data object arg
24
25      @param arg: Data object representing the function which gradient to be calculated.  #=========================================================
26      @param where: FunctionSpace in which the gradient will be. If None Function(dom) where dom is the  #   some little helpers
27                    domain of the Data object arg.  #=========================================================
28      """  def _testForZero(arg):
29      if isinstance(arg,escript.Data):     """
30         if where==None:     Returns True is arg is considered to be zero.
31            return arg.grad()     """
32         else:     if isinstance(arg,int):
33            return arg.grad(where)        return not arg>0
34      else:     elif isinstance(arg,float):
35         return arg*0.        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  def integrate(arg,what=None):  #=========================================================
44    def saveVTK(filename,**data):
45      """      """
46      @brief return the integral if the function represented by Data object arg over its domain.      writes arg into files in the vtk file format
47
48               saveVTK(<filename>,<data name 1>=<data object 1>,...,<data name n>=<data object n>)
49
50          This will create VTK files of the name <dir name>+<data name i>+"."+<extension> where <filename>=<dir name>+<extension>
51
@param arg
52      """      """
53      if not what==None:      ex=os.path.split(filename)
54         arg2=escript.Data(arg,what)      for i in data.keys():
55      else:         data[i].saveVTK(os.path.join(ex[0],i+"."+ex[1]))
arg2=arg
if arg2.getRank()==0:
return arg2.integrate()[0]
else:
return arg2.integrate()
56
57  def interpolate(arg,where):  #=========================================================
58    def saveDX(filename,**data):
59      """      """
60      @brief interpolates the function represented by Data object arg into the FunctionSpace where.      writes arg into file in the openDX file format
61
62               saveDX(<filename>,<data name 1>=<data object 1>,...,<data name n>=<data object n>)
63
64          This will create DX files of the name <dir name>+<data name i>+"."+<extension> where <filename>=<dir name>+<extension>
65
@param arg
@param where
66      """      """
67      if isinstance(arg,escript.Data):      ex=os.path.split(filename)
68         return arg.interpolate(where)      for i in data.keys():
69      else:         data[i].saveDX(os.path.join(ex[0],i+"."+ex[1]))
return arg
70
71  # functions returning Data objects:  #=========================================================
72
73  def transpose(arg,axis=None):  def exp(arg):
74      """      """
75      @brief returns the transpose of the Data object arg.      Applies the exponential function to arg.
76
77      @param arg      @param arg: argument
78      """      """
79      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
80         # hack for transpose         return symbols.Exp_Symbol(arg)
81         r=arg.getRank()      elif hasattr(arg,"exp"):
82         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)
83      else:      else:
84         if axis==None: axis=arg.rank/2         return numarray.exp(arg)
return numarray.transpose(arg,axis=axis)
85
86  def trace(arg):  def sqrt(arg):
87      """      """
88      @brief      Applies the squre root function to arg.
89
90      @param arg      @param arg: argument
91      """      """
92      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
93         # hack for trace         return symbols.Sqrt_Symbol(arg)
94         r=arg.getRank()      elif hasattr(arg,"sqrt"):
95         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()
96      else:      else:
97         return numarray.trace(arg)         return numarray.sqrt(arg)
98
99  def exp(arg):  def log(arg):
100      """      """
101      @brief      Applies the logarithmic function bases exp(1.) to arg
102
103      @param arg      @param arg: argument
104      """      """
105      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
106         return arg.exp()         return symbols.Log_Symbol(arg)
107        elif hasattr(arg,"log"):
108           return arg.log()
109      else:      else:
110         return numarray.exp(arg)         return numarray.log(arg)
111
112  def sqrt(arg):  def ln(arg):
113      """      """
114      @brief      Applies the natural logarithmic function to arg.
115
116      @param arg      @param arg: argument
117      """      """
118      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
119         return arg.sqrt()         return symbols.Ln_Symbol(arg)
120        elif hasattr(arg,"ln"):
121           return arg.log()
122      else:      else:
123         return numarray.sqrt(arg)         return numarray.log(arg)
124
125  def sin(arg):  def sin(arg):
126      """      """
127      @brief      Applies the sin function to arg.
128
129      @param arg      @param arg: argument
130      """      """
131      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
132           return symbols.Sin_Symbol(arg)
133        elif hasattr(arg,"sin"):
134         return arg.sin()         return arg.sin()
135      else:      else:
136         return numarray.sin(arg)         return numarray.sin(arg)
137
138    def cos(arg):
139        """
140        Applies the cos function to arg.
141
142        @param arg: argument
143        """
144        if isinstance(arg,symbols.Symbol):
145           return symbols.Cos_Symbol(arg)
146        elif hasattr(arg,"cos"):
147           return arg.cos()
148        else:
149           return numarray.cos(arg)
150
151  def tan(arg):  def tan(arg):
152      """      """
153      @brief      Applies the tan function to arg.
154
155      @param arg      @param arg: argument
156      """      """
157      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
158           return symbols.Tan_Symbol(arg)
159        elif hasattr(arg,"tan"):
160         return arg.tan()         return arg.tan()
161      else:      else:
162         return numarray.tan(arg)         return numarray.tan(arg)
163
164  def cos(arg):  def asin(arg):
165      """      """
166      @brief      Applies the asin function to arg.
167
168      @param arg      @param arg: argument
169      """      """
170      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
171         return arg.cos()         return symbols.Asin_Symbol(arg)
172        elif hasattr(arg,"asin"):
173           return arg.asin()
174      else:      else:
175         return numarray.cos(arg)         return numarray.asin(arg)
176
177    def acos(arg):
178        """
179        Applies the acos function to arg.
180
181        @param arg: argument
182        """
183        if isinstance(arg,symbols.Symbol):
184           return symbols.Acos_Symbol(arg)
185        elif hasattr(arg,"acos"):
186           return arg.acos()
187        else:
188           return numarray.acos(arg)
189
190    def atan(arg):
191        """
192        Applies the atan function to arg.
193
194        @param arg: argument
195        """
196        if isinstance(arg,symbols.Symbol):
197           return symbols.Atan_Symbol(arg)
198        elif hasattr(arg,"atan"):
199           return arg.atan()
200        else:
201           return numarray.atan(arg)
202
203    def sinh(arg):
204        """
205        Applies the sinh function to arg.
206
207        @param arg: argument
208        """
209        if isinstance(arg,symbols.Symbol):
210           return symbols.Sinh_Symbol(arg)
211        elif hasattr(arg,"sinh"):
212           return arg.sinh()
213        else:
214           return numarray.sinh(arg)
215
216    def cosh(arg):
217        """
218        Applies the cosh function to arg.
219
220        @param arg: argument
221        """
222        if isinstance(arg,symbols.Symbol):
223           return symbols.Cosh_Symbol(arg)
224        elif hasattr(arg,"cosh"):
225           return arg.cosh()
226        else:
227           return numarray.cosh(arg)
228
229    def tanh(arg):
230        """
231        Applies the tanh function to arg.
232
233        @param arg: argument
234        """
235        if isinstance(arg,symbols.Symbol):
236           return symbols.Tanh_Symbol(arg)
237        elif hasattr(arg,"tanh"):
238           return arg.tanh()
239        else:
240           return numarray.tanh(arg)
241
242    def asinh(arg):
243        """
244        Applies the asinh function to arg.
245
246        @param arg: argument
247        """
248        if isinstance(arg,symbols.Symbol):
249           return symbols.Asinh_Symbol(arg)
250        elif hasattr(arg,"asinh"):
251           return arg.asinh()
252        else:
253           return numarray.asinh(arg)
254
255    def acosh(arg):
256        """
257        Applies the acosh function to arg.
258
259        @param arg: argument
260        """
261        if isinstance(arg,symbols.Symbol):
262           return symbols.Acosh_Symbol(arg)
263        elif hasattr(arg,"acosh"):
264           return arg.acosh()
265        else:
266           return numarray.acosh(arg)
267
268    def atanh(arg):
269        """
270        Applies the atanh function to arg.
271
272        @param arg: argument
273        """
274        if isinstance(arg,symbols.Symbol):
275           return symbols.Atanh_Symbol(arg)
276        elif hasattr(arg,"atanh"):
277           return arg.atanh()
278        else:
279           return numarray.atanh(arg)
280
281    def sign(arg):
282        """
283        Applies the sign function to arg.
284
285        @param arg: argument
286        """
287        if isinstance(arg,symbols.Symbol):
288           return symbols.Sign_Symbol(arg)
289        elif hasattr(arg,"sign"):
290           return arg.sign()
291        else:
292           return numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))- \
293                  numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))
294
295  def maxval(arg):  def maxval(arg):
296      """      """
297      @brief      Returns the maximum value of argument arg.
298
299      @param arg      @param arg: argument
300      """      """
301      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
302           return symbols.Max_Symbol(arg)
303        elif hasattr(arg,"maxval"):
304         return arg.maxval()         return arg.maxval()
305      elif isinstance(arg,float) or isinstance(arg,int):      elif hasattr(arg,"max"):
return arg
else:
306         return arg.max()         return arg.max()
307        else:
308           return arg
309
310  def minval(arg):  def minval(arg):
311      """      """
312      @brief      Returns the minimum value of argument arg.
313
314      @param arg      @param arg: argument
315      """      """
316      if isinstance(arg,escript.Data):      if isinstance(arg,symbols.Symbol):
317           return symbols.Min_Symbol(arg)
318        elif hasattr(arg,"maxval"):
319         return arg.minval()         return arg.minval()
320      elif isinstance(arg,float) or isinstance(arg,int):      elif hasattr(arg,"min"):
321           return arg.min()
322        else:
323         return arg         return arg
324
325    def wherePositive(arg):
326        """
327        Returns the positive values of argument arg.
328
329        @param arg: argument
330        """
331        if _testForZero(arg):
332          return 0
333        elif isinstance(arg,symbols.Symbol):
334           return symbols.WherePositive_Symbol(arg)
335        elif hasattr(arg,"wherePositive"):
336           return arg.minval()
337        elif hasattr(arg,"wherePositive"):
338           numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))
339      else:      else:
340         return arg.min()         if arg>0:
341              return 1.
342           else:
343              return 0.
344
345    def whereNegative(arg):
346        """
347        Returns the negative values of argument arg.
348
349        @param arg: argument
350        """
351        if _testForZero(arg):
352          return 0
353        elif isinstance(arg,symbols.Symbol):
354           return symbols.WhereNegative_Symbol(arg)
355        elif hasattr(arg,"whereNegative"):
356           return arg.whereNegative()
357        elif hasattr(arg,"shape"):
358           numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))
359        else:
360           if arg<0:
361              return 1.
362           else:
363              return 0.
364
365    def maximum(arg0,arg1):
366        """
367        Return arg1 where arg1 is bigger then arg0 otherwise arg0 is returned.
368        """
369        m=whereNegative(arg0-arg1)
370        return m*arg1+(1.-m)*arg0
371
372    def minimum(arg0,arg1):
373        """
374        Return arg0 where arg1 is bigger then arg0 otherwise arg1 is returned.
375        """
376        m=whereNegative(arg0-arg1)
377        return m*arg0+(1.-m)*arg1
378
379    def outer(arg0,arg1):
380       if _testForZero(arg0) or _testForZero(arg1):
381          return 0
382       else:
383          if isinstance(arg0,symbols.Symbol) or isinstance(arg1,symbols.Symbol):
384            return symbols.Outer_Symbol(arg0,arg1)
385          elif _identifyShape(arg0)==() or _identifyShape(arg1)==():
386            return arg0*arg1
387          elif isinstance(arg0,numarray.NumArray) and isinstance(arg1,numarray.NumArray):
388            return numarray.outer(arg0,arg1)
389          else:
390            if arg0.getRank()==1 and arg1.getRank()==1:
391              out=escript.Data(0,(arg0.getShape()[0],arg1.getShape()[0]),arg1.getFunctionSpace())
392              for i in range(arg0.getShape()[0]):
393                for j in range(arg1.getShape()[0]):
394                    out[i,j]=arg0[i]*arg1[j]
395              return out
396            else:
397              raise ValueError,"outer is not fully implemented yet."
398
399    def interpolate(arg,where):
400        """
401        Interpolates the function into the FunctionSpace where.
402
403        @param arg:    interpolant
404        @param where:  FunctionSpace to interpolate to
405        """
406        if _testForZero(arg):
407          return 0
408        elif isinstance(arg,symbols.Symbol):
409           return symbols.Interpolated_Symbol(arg,where)
410        else:
411           return escript.Data(arg,where)
412
413    def div(arg,where=None):
414        """
415        Returns the divergence of arg at where.
416
417        @param arg:   Data object representing the function which gradient to
418                      be calculated.
419        @param where: FunctionSpace in which the gradient will be calculated.
420                      If not present or C{None} an appropriate default is used.
421        """
422        return trace(grad(arg,where))
423
424    def jump(arg):
425        """
426        Returns the jump of arg across a continuity.
427
428        @param arg:   Data object representing the function which gradient
429                      to be calculated.
430        """
431        d=arg.getDomain()
432        return arg.interpolate(escript.FunctionOnContactOne())-arg.interpolate(escript.FunctionOnContactZero())
433
434
435    def grad(arg,where=None):
436        """
437        Returns the spatial gradient of arg at where.
438
439        @param arg:   Data object representing the function which gradient
440                      to be calculated.
441        @param where: FunctionSpace in which the gradient will be calculated.
442                      If not present or C{None} an appropriate default is used.
443        """
444        if _testForZero(arg):
445          return 0
446        elif isinstance(arg,symbols.Symbol):
447           return symbols.Grad_Symbol(arg,where)
448        elif hasattr(arg,"grad"):
449           if where==None:
450              return arg.grad()
451           else:
452              return arg.grad(where)
453        else:
454           return arg*0.
455
456    def integrate(arg,where=None):
457        """
458        Return the integral if the function represented by Data object arg over
459        its domain.
460
461        @param arg:   Data object representing the function which is integrated.
462        @param where: FunctionSpace in which the integral is calculated.
463                      If not present or C{None} an appropriate default is used.
464        """
465        if _testForZero(arg):
466          return 0
467        elif isinstance(arg,symbols.Symbol):
468           return symbols.Integral_Symbol(arg,where)
469        else:
470           if not where==None: arg=escript.Data(arg,where)
471           if arg.getRank()==0:
472             return arg.integrate()[0]
473           else:
474             return arg.integrate()
475
476    #=============================
477    #
478    # wrapper for various functions: if the argument has attribute the function name
479    # as an argument it calls the corresponding methods. Otherwise the corresponding
480    # numarray function is called.
481
482    # functions involving the underlying Domain:
483
484
485    # functions returning Data objects:
486
487    def transpose(arg,axis=None):
488        """
489        Returns the transpose of the Data object arg.
490
491        @param arg:
492        """
493        if axis==None:
494           r=0
495           if hasattr(arg,"getRank"): r=arg.getRank()
496           if hasattr(arg,"rank"): r=arg.rank
497           axis=r/2
498        if isinstance(arg,symbols.Symbol):
499           return symbols.Transpose_Symbol(arg,axis=r)
500        if isinstance(arg,escript.Data):
501           # hack for transpose
502           r=arg.getRank()
503           if r!=2: raise ValueError,"Tranpose only avalaible for rank 2 objects"
504           s=arg.getShape()
505           out=escript.Data(0.,(s[1],s[0]),arg.getFunctionSpace())
506           for i in range(s[0]):
507              for j in range(s[1]):
508                 out[j,i]=arg[i,j]
509           return out
510           # end hack for transpose
511           return arg.transpose(axis)
512        else:
513           return numarray.transpose(arg,axis=axis)
514
515    def trace(arg,axis0=0,axis1=1):
516        """
517        Return
518
519        @param arg:
520        """
521        if isinstance(arg,symbols.Symbol):
522           s=list(arg.getShape())
523           s=tuple(s[0:axis0]+s[axis0+1:axis1]+s[axis1+1:])
524           return symbols.Trace_Symbol(arg,axis0=axis0,axis1=axis1)
525        elif isinstance(arg,escript.Data):
526           # hack for trace
527           s=arg.getShape()
528           if s[axis0]!=s[axis1]:
529               raise ValueError,"illegal axis in trace"
530           out=escript.Scalar(0.,arg.getFunctionSpace())
531           for i in range(s[axis0]):
532              out+=arg[i,i]
533           return out
534           # end hack for trace
535        else:
536           return numarray.trace(arg,axis0=axis0,axis1=axis1)
537
538  def length(arg):  def length(arg):
539      """      """
@brief
540
541      @param arg      @param arg:
542      """      """
543      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
544         if arg.isEmpty(): return escript.Data()         if arg.isEmpty(): return escript.Data()
545         if arg.getRank()==0:         if arg.getRank()==0:
546            return abs(arg)            return abs(arg)
547         elif arg.getRank()==1:         elif arg.getRank()==1:
548            sum=escript.Scalar(0,arg.getFunctionSpace())            out=escript.Scalar(0,arg.getFunctionSpace())
549            for i in range(arg.getShape()[0]):            for i in range(arg.getShape()[0]):
550               sum+=arg[i]**2               out+=arg[i]**2
551            return sqrt(sum)            return sqrt(out)
552         elif arg.getRank()==2:         elif arg.getRank()==2:
553            sum=escript.Scalar(0,arg.getFunctionSpace())            out=escript.Scalar(0,arg.getFunctionSpace())
554            for i in range(arg.getShape()[0]):            for i in range(arg.getShape()[0]):
555               for j in range(arg.getShape()[1]):               for j in range(arg.getShape()[1]):
556                  sum+=arg[i,j]**2                  out+=arg[i,j]**2
557            return sqrt(sum)            return sqrt(out)
558         elif arg.getRank()==3:         elif arg.getRank()==3:
559            sum=escript.Scalar(0,arg.getFunctionSpace())            out=escript.Scalar(0,arg.getFunctionSpace())
560            for i in range(arg.getShape()[0]):            for i in range(arg.getShape()[0]):
561               for j in range(arg.getShape()[1]):               for j in range(arg.getShape()[1]):
562                  for k in range(arg.getShape()[2]):                  for k in range(arg.getShape()[2]):
563                     sum+=arg[i,j,k]**2                     out+=arg[i,j,k]**2
564            return sqrt(sum)            return sqrt(out)
565         elif arg.getRank()==4:         elif arg.getRank()==4:
566            sum=escript.Scalar(0,arg.getFunctionSpace())            out=escript.Scalar(0,arg.getFunctionSpace())
567            for i in range(arg.getShape()[0]):            for i in range(arg.getShape()[0]):
568               for j in range(arg.getShape()[1]):               for j in range(arg.getShape()[1]):
569                  for k in range(arg.getShape()[2]):                  for k in range(arg.getShape()[2]):
570                     for l in range(arg.getShape()[3]):                     for l in range(arg.getShape()[3]):
571                        sum+=arg[i,j,k,l]**2                        out+=arg[i,j,k,l]**2
572            return sqrt(sum)            return sqrt(out)
573         else:         else:
574            raise SystemError,"length is not been implemented yet"            raise SystemError,"length is not been fully implemented yet"
575         # return arg.length()            # return arg.length()
576        elif isinstance(arg,float):
577           return abs(arg)
578      else:      else:
579         return sqrt((arg**2).sum())         return sqrt((arg**2).sum())
580
581  def sign(arg):  def deviator(arg):
582      """      """
583      @brief      @param arg:

@param arg
584      """      """
585      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
586         return arg.sign()          shape=arg.getShape()
587      else:      else:
588         return numarray.greater(arg,numarray.zeros(arg.shape))-numarray.less(arg,numarray.zeros(arg.shape))          shape=arg.shape
589        if len(shape)!=2:
590              raise ValueError,"Deviator requires rank 2 object"
591        if shape[0]!=shape[1]:
592              raise ValueError,"Deviator requires a square matrix"
593        return arg-1./(shape[0]*1.)*trace(arg)*kronecker(shape[0])
594
595    def inner(arg0,arg1):
596        """
597        @param arg0:
598        @param arg1:
599        """
600        if isinstance(arg0,escript.Data):
601           arg=arg0
602        else:
603           arg=arg1
604
605        out=escript.Scalar(0,arg.getFunctionSpace())
606        if arg.getRank()==0:
607              return arg0*arg1
608        elif arg.getRank()==1:
609             out=escript.Scalar(0,arg.getFunctionSpace())
610             for i in range(arg.getShape()[0]):
611                out+=arg0[i]*arg1[i]
612        elif arg.getRank()==2:
613            out=escript.Scalar(0,arg.getFunctionSpace())
614            for i in range(arg.getShape()[0]):
615               for j in range(arg.getShape()[1]):
616                  out+=arg0[i,j]*arg1[i,j]
617        elif arg.getRank()==3:
618            out=escript.Scalar(0,arg.getFunctionSpace())
619            for i in range(arg.getShape()[0]):
620                for j in range(arg.getShape()[1]):
621                   for k in range(arg.getShape()[2]):
622                      out+=arg0[i,j,k]*arg1[i,j,k]
623        elif arg.getRank()==4:
624            out=escript.Scalar(0,arg.getFunctionSpace())
625            for i in range(arg.getShape()[0]):
626               for j in range(arg.getShape()[1]):
627                  for k in range(arg.getShape()[2]):
628                     for l in range(arg.getShape()[3]):
629                        out+=arg0[i,j,k,l]*arg1[i,j,k,l]
630        else:
631              raise SystemError,"inner is not been implemented yet"
632        return out
633
634    def tensormult(arg0,arg1):
635        # check LinearPDE!!!!
636        raise SystemError,"tensormult is not implemented yet!"
637
638    def matrixmult(arg0,arg1):
639
640        if isinstance(arg1,numarray.NumArray) and isinstance(arg0,numarray.NumArray):
641            numarray.matrixmult(arg0,arg1)
642        else:
643          # escript.matmult(arg0,arg1)
644          if isinstance(arg1,escript.Data) and not isinstance(arg0,escript.Data):
645            arg0=escript.Data(arg0,arg1.getFunctionSpace())
646          elif isinstance(arg0,escript.Data) and not isinstance(arg1,escript.Data):
647            arg1=escript.Data(arg1,arg0.getFunctionSpace())
648          if arg0.getRank()==2 and arg1.getRank()==1:
649              out=escript.Data(0,(arg0.getShape()[0],),arg0.getFunctionSpace())
650              for i in range(arg0.getShape()[0]):
651                 for j in range(arg0.getShape()[1]):
652                   # uses Data object slicing, plus Data * and += operators
653                   out[i]+=arg0[i,j]*arg1[j]
654              return out
655          elif arg0.getRank()==1 and arg1.getRank()==1:
656              return inner(arg0,arg1)
657          else:
658              raise SystemError,"matrixmult is not fully implemented yet!"
659
660    #=========================================================
661  # reduction operations:  # reduction operations:
662    #=========================================================
663  def sum(arg):  def sum(arg):
664      """      """
665      @brief      @param arg:

@param arg
666      """      """
667      return arg.sum()      return arg.sum()
668
669  def sup(arg):  def sup(arg):
670      """      """
671      @brief      @param arg:

@param arg
672      """      """
673      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
674         return arg.sup()         return arg.sup()
# Line 303  def sup(arg): Line 679  def sup(arg):
679
680  def inf(arg):  def inf(arg):
681      """      """
682      @brief      @param arg:

@param arg
683      """      """
684      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
685         return arg.inf()         return arg.inf()
# Line 316  def inf(arg): Line 690  def inf(arg):
690
691  def L2(arg):  def L2(arg):
692      """      """
693      @brief returns the L2-norm of the      Returns the L2-norm of the argument
694
695      @param arg      @param arg:
696      """      """
697      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
698         return arg.L2()         return arg.L2()
# Line 329  def L2(arg): Line 703  def L2(arg):
703
704  def Lsup(arg):  def Lsup(arg):
705      """      """
706      @brief      @param arg:

@param arg
707      """      """
708      if isinstance(arg,escript.Data):      if isinstance(arg,escript.Data):
709         return arg.Lsup()         return arg.Lsup()
710      elif isinstance(arg,float) or isinstance(arg,int):      elif isinstance(arg,float) or isinstance(arg,int):
711         return abs(arg)         return abs(arg)
712      else:      else:
713         return max(numarray.abs(arg))         return numarray.abs(arg).max()
714
715  def dot(arg1,arg2):  def dot(arg0,arg1):
716      """      """
717      @brief      @param arg0:
718        @param arg1:
@param arg
719      """      """
720      if isinstance(arg1,escript.Data):      if isinstance(arg0,escript.Data):
721         return arg1.dot(arg2)         return arg0.dot(arg1)
722      elif isinstance(arg1,escript.Data):      elif isinstance(arg1,escript.Data):
723         return arg2.dot(arg1)         return arg1.dot(arg0)
724      else:      else:
725         return numarray.dot(arg1,arg2)         return numarray.dot(arg0,arg1)
726
727    def kronecker(d):
728       if hasattr(d,"getDim"):
729          return numarray.identity(d.getDim())*1.
730       else:
731          return numarray.identity(d)*1.
732
733    def unit(i,d):
734       """
735       Return a unit vector of dimension d with nonzero index i.
736
737       @param d: dimension
738       @param i: index
739       """
740       e = numarray.zeros((d,),numarray.Float)
741       e[i] = 1.0
742       return e
743    #
744    # \$Log\$
745    # Revision 1.18  2005/09/15 03:44:19  jgs
746    # Merge of development branch dev-02 back to main trunk on 2005-09-15
747    #
748    # Revision 1.17  2005/09/01 03:31:28  jgs
749    # Merge of development branch dev-02 back to main trunk on 2005-09-01
750    #
751    # Revision 1.16  2005/08/23 01:24:28  jgs
752    # Merge of development branch dev-02 back to main trunk on 2005-08-23
753    #
754    # Revision 1.15  2005/08/12 01:45:36  jgs
755    # erge of development branch dev-02 back to main trunk on 2005-08-12
756    #
757    # Revision 1.14.2.13  2005/09/12 03:32:14  gross
758    # test_visualiztion has been aded to mk
759    #
760    # Revision 1.14.2.12  2005/09/09 01:56:24  jgs
761    # added implementations of acos asin atan sinh cosh tanh asinh acosh atanh
762    # and some associated testing
763    #
764    # Revision 1.14.2.11  2005/09/08 08:28:39  gross
765    # some cleanup in savevtk
766    #
767    # Revision 1.14.2.10  2005/09/08 00:25:32  gross
768    # test for finley mesh generators added
769    #
770    # Revision 1.14.2.9  2005/09/07 10:32:05  gross
771    # Symbols removed from util and put into symmbols.py.
772    #
773    # Revision 1.14.2.8  2005/08/26 05:06:37  cochrane
774    # Corrected errors in docstrings.  Improved output formatting of docstrings.
775    # Other minor improvements to code and docs (eg spelling etc).
776    #
777    # Revision 1.14.2.7  2005/08/26 04:45:40  cochrane
778    # Fixed and tidied markup and docstrings.  Some *Symbol classes were defined
779    # as functions, so changed them to classes (hopefully this was the right thing
780    # to do).
781    #
782    # Revision 1.14.2.6  2005/08/26 04:30:13  gross
783    # gneric unit testing for linearPDE
784    #
785    # Revision 1.14.2.5  2005/08/24 02:02:52  gross
786    # jump function added
787    #
788    # Revision 1.14.2.4  2005/08/18 04:39:32  gross
789    # the constants have been removed from util.py as they not needed anymore. PDE related constants are accessed through LinearPDE attributes now
790    #
791    # Revision 1.14.2.3  2005/08/03 09:55:33  gross
792    # ContactTest is passing now./mk install!
793    #
794    # Revision 1.14.2.2  2005/08/02 03:15:14  gross
795    # bug inb trace fixed!
796    #
797    # Revision 1.14.2.1  2005/07/29 07:10:28  gross
798    # new functions in util and a new pde type in linearPDEs
799    #
800    # Revision 1.2.2.21  2005/07/28 04:19:23  gross
801    # new functions maximum and minimum introduced.
802    #
803    # Revision 1.2.2.20  2005/07/25 01:26:27  gross
804    # bug in inner fixed
805    #
806    # Revision 1.2.2.19  2005/07/21 04:01:28  jgs
807    # minor comment fixes
808    #
809    # Revision 1.2.2.18  2005/07/21 01:02:43  jgs
810    # commit ln() updates to development branch version
811    #
812    # Revision 1.12  2005/07/20 06:14:58  jgs
813    # added ln(data) style wrapper for data.ln() - also added corresponding
814    # implementation of Ln_Symbol class (not sure if this is right though)
815    #
816    # Revision 1.11  2005/07/08 04:07:35  jgs
817    # Merge of development branch back to main trunk on 2005-07-08
818    #
819    # Revision 1.10  2005/06/09 05:37:59  jgs
820    # Merge of development branch back to main trunk on 2005-06-09
821    #
822    # Revision 1.2.2.17  2005/07/07 07:28:58  gross
823    # some stuff added to util.py to improve functionality
824    #
825    # Revision 1.2.2.16  2005/06/30 01:53:55  gross
826    # a bug in coloring fixed
827    #
828    # Revision 1.2.2.15  2005/06/29 02:36:43  gross
829    # Symbols have been introduced and some function clarified. needs much more work
830    #
831    # Revision 1.2.2.14  2005/05/20 04:05:23  gross
832    # some work on a darcy flow started
833    #
834    # Revision 1.2.2.13  2005/03/16 05:17:58  matt
835    # Implemented unit(idx, dim) to create cartesian unit basis vectors to
836    # complement kronecker(dim) function.
837    #
838    # Revision 1.2.2.12  2005/03/10 08:14:37  matt
839    # Added non-member Linf utility function to complement Data::Linf().
840    #
841    # Revision 1.2.2.11  2005/02/17 05:53:25  gross
842    # some bug in saveDX fixed: in fact the bug was in
843    # DataC/getDataPointShape
844    #
845    # Revision 1.2.2.10  2005/01/11 04:59:36  gross
846    # automatic interpolation in integrate switched off
847    #
848    # Revision 1.2.2.9  2005/01/11 03:38:13  gross
849    # 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.
850    #
851    # Revision 1.2.2.8  2005/01/05 04:21:41  gross
852    # FunctionSpace checking/matchig in slicing added
853    #
854    # Revision 1.2.2.7  2004/12/29 05:29:59  gross
855    # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()
856    #
857    # Revision 1.2.2.6  2004/12/24 06:05:41  gross
858    # some changes in linearPDEs to add AdevectivePDE
859    #
860    # Revision 1.2.2.5  2004/12/17 00:06:53  gross
861    # mk sets ESYS_ROOT is undefined
862    #
863    # Revision 1.2.2.4  2004/12/07 03:19:51  gross
864    # options for GMRES and PRES20 added
865    #
866    # Revision 1.2.2.3  2004/12/06 04:55:18  gross
867    # function wraper extended
868    #
869    # Revision 1.2.2.2  2004/11/22 05:44:07  gross
870    # a few more unitary functions have been added but not implemented in Data yet
871    #
872    # Revision 1.2.2.1  2004/11/12 06:58:15  gross
873    # 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
874    #
875    # Revision 1.2  2004/10/27 00:23:36  jgs
876    # fixed minor syntax error
877    #
878    # Revision 1.1.1.1  2004/10/26 06:53:56  jgs
879    # initial import of project esys2
880    #
881    # Revision 1.1.2.3  2004/10/26 06:43:48  jgs
882    # committing Lutz's and Paul's changes to brach jgs
883    #
884    # Revision 1.1.4.1  2004/10/20 05:32:51  cochrane
885    # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.
886    #
887    # Revision 1.1  2004/08/05 03:58:27  gross
888    # Bug in Assemble_NodeCoordinates fixed
889    #
890    #
891
892    # vim: expandtab shiftwidth=4:

Legend:
 Removed from v.108 changed lines Added in v.150

 ViewVC Help Powered by ViewVC 1.1.26