/[escript]/trunk/escript/py_src/util.py
ViewVC logotype

Annotation of /trunk/escript/py_src/util.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (14 years ago) by jgs
Original Path: trunk/esys2/escript/py_src/util.py
File MIME type: text/x-python
File size: 25359 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

1 jgs 82 # $Id$
2    
3     ## @file util.py
4    
5     """
6 jgs 149 Utility functions for escript
7 jgs 123
8 jgs 149 @todo:
9 jgs 123
10 jgs 149 - 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 jgs 82 """
19    
20     import numarray
21 jgs 102 import escript
22 jgs 150 import symbols
23     import os
24 jgs 124
25 jgs 123 #=========================================================
26     # some little helpers
27     #=========================================================
28     def _testForZero(arg):
29 jgs 149 """
30     Returns True is arg is considered to be zero.
31     """
32 jgs 123 if isinstance(arg,int):
33     return not arg>0
34     elif isinstance(arg,float):
35     return not arg>0.
36 jgs 150 elif isinstance(arg,numarray.NumArray):
37 jgs 123 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 jgs 150 #=========================================================
44     def saveVTK(filename,**data):
45     """
46     writes arg into files in the vtk file format
47 jgs 123
48 jgs 150 saveVTK(<filename>,<data name 1>=<data object 1>,...,<data name n>=<data object n>)
49 jgs 123
50 jgs 150 This will create VTK files of the name <dir name>+<data name i>+"."+<extension> where <filename>=<dir name>+<extension>
51    
52 jgs 149 """
53 jgs 150 ex=os.path.split(filename)
54     for i in data.keys():
55     data[i].saveVTK(os.path.join(ex[0],i+"."+ex[1]))
56 jgs 123
57 jgs 150 #=========================================================
58     def saveDX(filename,**data):
59 jgs 149 """
60 jgs 150 writes arg into file in the openDX file format
61 jgs 149
62 jgs 150 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    
66 jgs 149 """
67 jgs 150 ex=os.path.split(filename)
68     for i in data.keys():
69     data[i].saveDX(os.path.join(ex[0],i+"."+ex[1]))
70 jgs 124
71 jgs 123 #=========================================================
72    
73     def exp(arg):
74 jgs 82 """
75 jgs 149 Applies the exponential function to arg.
76    
77     @param arg: argument
78 jgs 123 """
79 jgs 150 if isinstance(arg,symbols.Symbol):
80     return symbols.Exp_Symbol(arg)
81 jgs 123 elif hasattr(arg,"exp"):
82     return arg.exp()
83 jgs 108 else:
84 jgs 123 return numarray.exp(arg)
85 jgs 82
86 jgs 123 def sqrt(arg):
87 jgs 82 """
88 jgs 149 Applies the squre root function to arg.
89    
90     @param arg: argument
91 jgs 123 """
92 jgs 150 if isinstance(arg,symbols.Symbol):
93     return symbols.Sqrt_Symbol(arg)
94 jgs 123 elif hasattr(arg,"sqrt"):
95     return arg.sqrt()
96     else:
97     return numarray.sqrt(arg)
98 jgs 82
99 jgs 123 def log(arg):
100 jgs 82 """
101 jgs 149 Applies the logarithmic function bases exp(1.) to arg
102    
103     @param arg: argument
104 jgs 123 """
105 jgs 150 if isinstance(arg,symbols.Symbol):
106     return symbols.Log_Symbol(arg)
107 jgs 123 elif hasattr(arg,"log"):
108     return arg.log()
109 jgs 108 else:
110 jgs 123 return numarray.log(arg)
111 jgs 82
112 jgs 124 def ln(arg):
113     """
114 jgs 149 Applies the natural logarithmic function to arg.
115    
116     @param arg: argument
117 jgs 124 """
118 jgs 150 if isinstance(arg,symbols.Symbol):
119     return symbols.Ln_Symbol(arg)
120 jgs 124 elif hasattr(arg,"ln"):
121     return arg.log()
122     else:
123     return numarray.log(arg)
124    
125 jgs 123 def sin(arg):
126 jgs 82 """
127 jgs 149 Applies the sin function to arg.
128    
129     @param arg: argument
130 jgs 123 """
131 jgs 150 if isinstance(arg,symbols.Symbol):
132     return symbols.Sin_Symbol(arg)
133 jgs 123 elif hasattr(arg,"sin"):
134     return arg.sin()
135     else:
136     return numarray.sin(arg)
137 jgs 82
138 jgs 123 def cos(arg):
139 jgs 82 """
140 jgs 149 Applies the cos function to arg.
141    
142     @param arg: argument
143 jgs 123 """
144 jgs 150 if isinstance(arg,symbols.Symbol):
145     return symbols.Cos_Symbol(arg)
146 jgs 123 elif hasattr(arg,"cos"):
147     return arg.cos()
148 jgs 82 else:
149 jgs 123 return numarray.cos(arg)
150 jgs 82
151 jgs 123 def tan(arg):
152 jgs 82 """
153 jgs 149 Applies the tan function to arg.
154    
155     @param arg: argument
156 jgs 123 """
157 jgs 150 if isinstance(arg,symbols.Symbol):
158     return symbols.Tan_Symbol(arg)
159 jgs 123 elif hasattr(arg,"tan"):
160     return arg.tan()
161     else:
162     return numarray.tan(arg)
163 jgs 82
164 jgs 150 def asin(arg):
165     """
166     Applies the asin function to arg.
167 jgs 123
168 jgs 150 @param arg: argument
169     """
170     if isinstance(arg,symbols.Symbol):
171     return symbols.Asin_Symbol(arg)
172     elif hasattr(arg,"asin"):
173     return arg.asin()
174     else:
175     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 jgs 123 def sign(arg):
282 jgs 82 """
283 jgs 149 Applies the sign function to arg.
284    
285     @param arg: argument
286 jgs 123 """
287 jgs 150 if isinstance(arg,symbols.Symbol):
288     return symbols.Sign_Symbol(arg)
289 jgs 123 elif hasattr(arg,"sign"):
290     return arg.sign()
291 jgs 82 else:
292 jgs 123 return numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))- \
293     numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))
294 jgs 82
295 jgs 123 def maxval(arg):
296 jgs 82 """
297 jgs 149 Returns the maximum value of argument arg.
298    
299     @param arg: argument
300 jgs 123 """
301 jgs 150 if isinstance(arg,symbols.Symbol):
302     return symbols.Max_Symbol(arg)
303 jgs 123 elif hasattr(arg,"maxval"):
304     return arg.maxval()
305     elif hasattr(arg,"max"):
306     return arg.max()
307     else:
308     return arg
309 jgs 82
310 jgs 123 def minval(arg):
311 jgs 82 """
312 jgs 149 Returns the minimum value of argument arg.
313    
314     @param arg: argument
315 jgs 123 """
316 jgs 150 if isinstance(arg,symbols.Symbol):
317     return symbols.Min_Symbol(arg)
318 jgs 123 elif hasattr(arg,"maxval"):
319     return arg.minval()
320     elif hasattr(arg,"min"):
321     return arg.min()
322 jgs 82 else:
323 jgs 123 return arg
324 jgs 82
325 jgs 123 def wherePositive(arg):
326 jgs 82 """
327 jgs 149 Returns the positive values of argument arg.
328    
329     @param arg: argument
330 jgs 123 """
331     if _testForZero(arg):
332     return 0
333 jgs 150 elif isinstance(arg,symbols.Symbol):
334     return symbols.WherePositive_Symbol(arg)
335 jgs 123 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:
340     if arg>0:
341     return 1.
342     else:
343     return 0.
344 jgs 82
345 jgs 123 def whereNegative(arg):
346 jgs 82 """
347 jgs 149 Returns the negative values of argument arg.
348    
349     @param arg: argument
350 jgs 123 """
351     if _testForZero(arg):
352     return 0
353 jgs 150 elif isinstance(arg,symbols.Symbol):
354     return symbols.WhereNegative_Symbol(arg)
355 jgs 123 elif hasattr(arg,"whereNegative"):
356     return arg.whereNegative()
357     elif hasattr(arg,"shape"):
358     numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))
359 jgs 88 else:
360 jgs 123 if arg<0:
361     return 1.
362     else:
363     return 0.
364 jgs 82
365 jgs 147 def maximum(arg0,arg1):
366 jgs 149 """
367     Return arg1 where arg1 is bigger then arg0 otherwise arg0 is returned.
368     """
369 jgs 147 m=whereNegative(arg0-arg1)
370     return m*arg1+(1.-m)*arg0
371    
372     def minimum(arg0,arg1):
373 jgs 149 """
374     Return arg0 where arg1 is bigger then arg0 otherwise arg1 is returned.
375     """
376 jgs 147 m=whereNegative(arg0-arg1)
377     return m*arg0+(1.-m)*arg1
378    
379 jgs 123 def outer(arg0,arg1):
380     if _testForZero(arg0) or _testForZero(arg1):
381     return 0
382     else:
383 jgs 150 if isinstance(arg0,symbols.Symbol) or isinstance(arg1,symbols.Symbol):
384     return symbols.Outer_Symbol(arg0,arg1)
385 jgs 123 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 jgs 82 """
401 jgs 149 Interpolates the function into the FunctionSpace where.
402 jgs 82
403 jgs 149 @param arg: interpolant
404     @param where: FunctionSpace to interpolate to
405 jgs 82 """
406 jgs 123 if _testForZero(arg):
407     return 0
408 jgs 150 elif isinstance(arg,symbols.Symbol):
409     return symbols.Interpolated_Symbol(arg,where)
410 jgs 82 else:
411 jgs 123 return escript.Data(arg,where)
412 jgs 82
413 jgs 147 def div(arg,where=None):
414     """
415 jgs 149 Returns the divergence of arg at where.
416 jgs 147
417 jgs 149 @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 jgs 147 """
422     return trace(grad(arg,where))
423    
424 jgs 149 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 jgs 123 def grad(arg,where=None):
436 jgs 102 """
437 jgs 149 Returns the spatial gradient of arg at where.
438 jgs 102
439 jgs 149 @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 jgs 102 """
444 jgs 123 if _testForZero(arg):
445     return 0
446 jgs 150 elif isinstance(arg,symbols.Symbol):
447     return symbols.Grad_Symbol(arg,where)
448 jgs 123 elif hasattr(arg,"grad"):
449     if where==None:
450     return arg.grad()
451     else:
452     return arg.grad(where)
453 jgs 102 else:
454 jgs 123 return arg*0.
455 jgs 102
456 jgs 123 def integrate(arg,where=None):
457 jgs 82 """
458 jgs 149 Return the integral if the function represented by Data object arg over
459     its domain.
460 jgs 82
461 jgs 123 @param arg: Data object representing the function which is integrated.
462 jgs 149 @param where: FunctionSpace in which the integral is calculated.
463     If not present or C{None} an appropriate default is used.
464 jgs 82 """
465 jgs 123 if _testForZero(arg):
466     return 0
467 jgs 150 elif isinstance(arg,symbols.Symbol):
468     return symbols.Integral_Symbol(arg,where)
469 jgs 123 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 jgs 82
476 jgs 123 #=============================
477     #
478     # wrapper for various functions: if the argument has attribute the function name
479 jgs 124 # as an argument it calls the corresponding methods. Otherwise the corresponding
480     # numarray function is called.
481    
482 jgs 123 # functions involving the underlying Domain:
483    
484    
485     # functions returning Data objects:
486    
487     def transpose(arg,axis=None):
488 jgs 82 """
489 jgs 149 Returns the transpose of the Data object arg.
490 jgs 82
491 jgs 149 @param arg:
492 jgs 82 """
493 jgs 123 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 jgs 150 if isinstance(arg,symbols.Symbol):
499     return symbols.Transpose_Symbol(arg,axis=r)
500 jgs 102 if isinstance(arg,escript.Data):
501 jgs 123 # 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 jgs 102 else:
513 jgs 123 return numarray.transpose(arg,axis=axis)
514 jgs 82
515 jgs 123 def trace(arg,axis0=0,axis1=1):
516 jgs 82 """
517 jgs 149 Return
518 jgs 82
519 jgs 149 @param arg:
520 jgs 82 """
521 jgs 150 if isinstance(arg,symbols.Symbol):
522 jgs 123 s=list(arg.getShape())
523     s=tuple(s[0:axis0]+s[axis0+1:axis1]+s[axis1+1:])
524 jgs 150 return symbols.Trace_Symbol(arg,axis0=axis0,axis1=axis1)
525 jgs 123 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 jgs 147 for i in range(s[axis0]):
532     out+=arg[i,i]
533 jgs 123 return out
534 jgs 126 # end hack for trace
535 jgs 102 else:
536 jgs 123 return numarray.trace(arg,axis0=axis0,axis1=axis1)
537 jgs 82
538 jgs 102 def length(arg):
539     """
540    
541 jgs 149 @param arg:
542 jgs 102 """
543     if isinstance(arg,escript.Data):
544 jgs 108 if arg.isEmpty(): return escript.Data()
545     if arg.getRank()==0:
546     return abs(arg)
547     elif arg.getRank()==1:
548 jgs 147 out=escript.Scalar(0,arg.getFunctionSpace())
549 jgs 104 for i in range(arg.getShape()[0]):
550 jgs 147 out+=arg[i]**2
551     return sqrt(out)
552 jgs 108 elif arg.getRank()==2:
553 jgs 147 out=escript.Scalar(0,arg.getFunctionSpace())
554 jgs 108 for i in range(arg.getShape()[0]):
555     for j in range(arg.getShape()[1]):
556 jgs 147 out+=arg[i,j]**2
557     return sqrt(out)
558 jgs 108 elif arg.getRank()==3:
559 jgs 147 out=escript.Scalar(0,arg.getFunctionSpace())
560 jgs 108 for i in range(arg.getShape()[0]):
561     for j in range(arg.getShape()[1]):
562     for k in range(arg.getShape()[2]):
563 jgs 147 out+=arg[i,j,k]**2
564     return sqrt(out)
565 jgs 108 elif arg.getRank()==4:
566 jgs 147 out=escript.Scalar(0,arg.getFunctionSpace())
567 jgs 108 for i in range(arg.getShape()[0]):
568     for j in range(arg.getShape()[1]):
569     for k in range(arg.getShape()[2]):
570     for l in range(arg.getShape()[3]):
571 jgs 147 out+=arg[i,j,k,l]**2
572     return sqrt(out)
573 jgs 104 else:
574 jgs 126 raise SystemError,"length is not been fully implemented yet"
575     # return arg.length()
576 jgs 147 elif isinstance(arg,float):
577     return abs(arg)
578 jgs 102 else:
579     return sqrt((arg**2).sum())
580    
581 jgs 113 def deviator(arg):
582     """
583 jgs 149 @param arg:
584 jgs 113 """
585     if isinstance(arg,escript.Data):
586     shape=arg.getShape()
587     else:
588     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 jgs 123 def inner(arg0,arg1):
596 jgs 113 """
597 jgs 149 @param arg0:
598     @param arg1:
599 jgs 113 """
600 jgs 147 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 jgs 123 return arg0*arg1
608 jgs 147 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 jgs 113 else:
631     raise SystemError,"inner is not been implemented yet"
632 jgs 147 return out
633 jgs 113
634 jgs 149 def tensormult(arg0,arg1):
635     # check LinearPDE!!!!
636     raise SystemError,"tensormult is not implemented yet!"
637    
638 jgs 123 def matrixmult(arg0,arg1):
639 jgs 102
640 jgs 123 if isinstance(arg1,numarray.NumArray) and isinstance(arg0,numarray.NumArray):
641     numarray.matrixmult(arg0,arg1)
642 jgs 102 else:
643 jgs 123 # 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 jgs 126 # uses Data object slicing, plus Data * and += operators
653 jgs 123 out[i]+=arg0[i,j]*arg1[j]
654     return out
655 jgs 147 elif arg0.getRank()==1 and arg1.getRank()==1:
656     return inner(arg0,arg1)
657 jgs 123 else:
658     raise SystemError,"matrixmult is not fully implemented yet!"
659 jgs 124
660 jgs 123 #=========================================================
661 jgs 102 # reduction operations:
662 jgs 123 #=========================================================
663 jgs 102 def sum(arg):
664     """
665 jgs 149 @param arg:
666 jgs 102 """
667     return arg.sum()
668    
669 jgs 82 def sup(arg):
670     """
671 jgs 149 @param arg:
672 jgs 82 """
673 jgs 102 if isinstance(arg,escript.Data):
674     return arg.sup()
675 jgs 108 elif isinstance(arg,float) or isinstance(arg,int):
676     return arg
677 jgs 102 else:
678     return arg.max()
679 jgs 82
680     def inf(arg):
681     """
682 jgs 149 @param arg:
683 jgs 82 """
684 jgs 102 if isinstance(arg,escript.Data):
685     return arg.inf()
686 jgs 108 elif isinstance(arg,float) or isinstance(arg,int):
687     return arg
688 jgs 102 else:
689     return arg.min()
690 jgs 82
691 jgs 102 def L2(arg):
692 jgs 82 """
693 jgs 149 Returns the L2-norm of the argument
694 jgs 82
695 jgs 149 @param arg:
696 jgs 82 """
697 jgs 108 if isinstance(arg,escript.Data):
698     return arg.L2()
699     elif isinstance(arg,float) or isinstance(arg,int):
700     return abs(arg)
701     else:
702     return numarry.sqrt(dot(arg,arg))
703 jgs 82
704 jgs 102 def Lsup(arg):
705 jgs 82 """
706 jgs 149 @param arg:
707 jgs 82 """
708 jgs 102 if isinstance(arg,escript.Data):
709     return arg.Lsup()
710 jgs 108 elif isinstance(arg,float) or isinstance(arg,int):
711     return abs(arg)
712 jgs 102 else:
713 jgs 149 return numarray.abs(arg).max()
714 jgs 82
715 jgs 123 def dot(arg0,arg1):
716 jgs 117 """
717 jgs 149 @param arg0:
718     @param arg1:
719 jgs 117 """
720 jgs 123 if isinstance(arg0,escript.Data):
721     return arg0.dot(arg1)
722 jgs 102 elif isinstance(arg1,escript.Data):
723 jgs 123 return arg1.dot(arg0)
724 jgs 102 else:
725 jgs 123 return numarray.dot(arg0,arg1)
726 jgs 113
727     def kronecker(d):
728 jgs 122 if hasattr(d,"getDim"):
729 jgs 147 return numarray.identity(d.getDim())*1.
730 jgs 122 else:
731 jgs 147 return numarray.identity(d)*1.
732 jgs 117
733     def unit(i,d):
734     """
735 jgs 149 Return a unit vector of dimension d with nonzero index i.
736    
737     @param d: dimension
738     @param i: index
739 jgs 117 """
740 jgs 123 e = numarray.zeros((d,),numarray.Float)
741 jgs 117 e[i] = 1.0
742     return e
743 jgs 123 #
744 jgs 122 # $Log$
745 jgs 150 # 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 jgs 149 # 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 jgs 148 # 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 jgs 147 # 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 jgs 142 #
757 jgs 150 # 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 jgs 149 # 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 jgs 148 # 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 jgs 147 # Revision 1.14.2.3 2005/08/03 09:55:33 gross
792     # ContactTest is passing now./mk install!
793 jgs 126 #
794 jgs 147 # Revision 1.14.2.2 2005/08/02 03:15:14 gross
795     # bug inb trace fixed!
796 jgs 124 #
797 jgs 147 # 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 jgs 123 #
800 jgs 147 # Revision 1.2.2.21 2005/07/28 04:19:23 gross
801     # new functions maximum and minimum introduced.
802 jgs 122 #
803 jgs 142 # Revision 1.2.2.20 2005/07/25 01:26:27 gross
804     # bug in inner fixed
805     #
806 jgs 126 # 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 jgs 123 # 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 jgs 122 # 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 jgs 149
892     # vim: expandtab shiftwidth=4:

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26