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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 149 - (hide annotations)
Thu Sep 1 03:31:39 2005 UTC (14 years, 1 month ago) by jgs
File MIME type: text/x-python
File size: 43121 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-01

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 124
23 jgs 123 #===========================================================
24 jgs 124 # a simple tool box to deal with _differentials of functions
25 jgs 123 #===========================================================
26 jgs 113
27 jgs 123 class Symbol:
28 jgs 149 """
29     Symbol class.
30     """
31 jgs 123 def __init__(self,name="symbol",shape=(),dim=3,args=[]):
32 jgs 149 """
33     Creates an instance of a symbol of shape shape and spatial dimension
34     dim.
35    
36     The symbol may depending on a list of arguments args which may be
37     symbols or other objects. name gives the name of the symbol.
38     """
39    
40 jgs 123 self.__args=args
41     self.__name=name
42     self.__shape=shape
43     if hasattr(dim,"getDim"):
44     self.__dim=dim.getDim()
45     else:
46     self.__dim=dim
47     #
48     self.__cache_val=None
49     self.__cache_argval=None
50 jgs 82
51 jgs 123 def getArgument(self,i):
52 jgs 149 """
53     Returns the i-th argument.
54     """
55 jgs 123 return self.__args[i]
56    
57     def getDim(self):
58 jgs 149 """
59     Returns the spatial dimension of the symbol.
60     """
61 jgs 123 return self.__dim
62    
63     def getRank(self):
64 jgs 149 """
65     Returns the rank of the symbol.
66     """
67 jgs 123 return len(self.getShape())
68    
69     def getShape(self):
70 jgs 149 """
71     Returns the shape of the symbol.
72     """
73 jgs 123 return self.__shape
74    
75     def getEvaluatedArguments(self,argval):
76 jgs 149 """
77     Returns the list of evaluated arguments by subsituting symbol u by
78     argval[u].
79     """
80 jgs 123 if argval==self.__cache_argval:
81     print "%s: cached value used"%self
82     return self.__cache_val
83     else:
84     out=[]
85     for a in self.__args:
86     if isinstance(a,Symbol):
87     out.append(a.eval(argval))
88     else:
89     out.append(a)
90     self.__cache_argval=argval
91     self.__cache_val=out
92     return out
93    
94     def getDifferentiatedArguments(self,arg):
95 jgs 149 """
96     Returns the list of the arguments _differentiated by arg.
97     """
98 jgs 123 out=[]
99     for a in self.__args:
100     if isinstance(a,Symbol):
101     out.append(a.diff(arg))
102     else:
103     out.append(0)
104     return out
105    
106     def diff(self,arg):
107 jgs 149 """
108     Returns the _differention of self by arg.
109     """
110 jgs 123 if self==arg:
111     out=numarray.zeros(tuple(2*list(self.getShape())),numarray.Float)
112     if self.getRank()==0:
113     out=1.
114     elif self.getRank()==1:
115     for i0 in range(self.getShape()[0]):
116     out[i0,i0]=1.
117     elif self.getRank()==2:
118     for i0 in range(self.getShape()[0]):
119     for i1 in range(self.getShape()[1]):
120     out[i0,i1,i0,i1]=1.
121     elif self.getRank()==3:
122     for i0 in range(self.getShape()[0]):
123     for i1 in range(self.getShape()[1]):
124     for i2 in range(self.getShape()[2]):
125     out[i0,i1,i2,i0,i1,i2]=1.
126     elif self.getRank()==4:
127     for i0 in range(self.getShape()[0]):
128     for i1 in range(self.getShape()[1]):
129     for i2 in range(self.getShape()[2]):
130     for i3 in range(self.getShape()[3]):
131     out[i0,i1,i2,i3,i0,i1,i2,i3]=1.
132     else:
133     raise ValueError,"differential support rank<5 only."
134     return out
135 jgs 108 else:
136 jgs 123 return self._diff(arg)
137 jgs 82
138 jgs 123 def _diff(self,arg):
139 jgs 149 """
140     Return derivate of self with respect to arg (!=self).
141    
142     This method is overwritten by a particular symbol.
143     """
144 jgs 123 return 0
145 jgs 82
146 jgs 123 def eval(self,argval):
147 jgs 149 """
148     Subsitutes symbol u in self by argval[u] and returns the result. If
149     self is not a key of argval then self is returned.
150     """
151 jgs 123 if argval.has_key(self):
152     return argval[self]
153     else:
154     return self
155    
156     def __str__(self):
157 jgs 149 """
158     Returns a string representation of the symbol.
159     """
160 jgs 123 return self.__name
161    
162     def __add__(self,other):
163 jgs 149 """
164     Adds other to symbol self. if _testForZero(other) self is returned.
165     """
166 jgs 123 if _testForZero(other):
167     return self
168     else:
169     a=_matchShape([self,other])
170     return Add_Symbol(a[0],a[1])
171    
172     def __radd__(self,other):
173 jgs 149 """
174     Adds other to symbol self. if _testForZero(other) self is returned.
175     """
176 jgs 123 return self+other
177    
178     def __neg__(self):
179 jgs 149 """
180     Returns -self.
181     """
182 jgs 123 return self*(-1.)
183    
184     def __pos__(self):
185 jgs 149 """
186     Returns +self.
187     """
188 jgs 123 return self
189    
190     def __abs__(self):
191 jgs 149 """
192     Returns absolute value.
193     """
194 jgs 123 return Abs_Symbol(self)
195    
196     def __sub__(self,other):
197 jgs 149 """
198     Subtracts other from symbol self.
199    
200     If _testForZero(other) self is returned.
201     """
202 jgs 123 if _testForZero(other):
203     return self
204     else:
205     return self+(-other)
206    
207     def __rsub__(self,other):
208 jgs 149 """
209     Subtracts symbol self from other.
210     """
211 jgs 123 return -self+other
212    
213     def __div__(self,other):
214 jgs 149 """
215     Divides symbol self by other.
216     """
217 jgs 123 if isinstance(other,Symbol):
218     a=_matchShape([self,other])
219     return Div_Symbol(a[0],a[1])
220     else:
221     return self*(1./other)
222    
223     def __rdiv__(self,other):
224 jgs 149 """
225     Dived other by symbol self. if _testForZero(other) 0 is returned.
226     """
227 jgs 123 if _testForZero(other):
228     return 0
229     else:
230     a=_matchShape([self,other])
231     return Div_Symbol(a[0],a[1])
232    
233     def __pow__(self,other):
234 jgs 149 """
235     Raises symbol self to the power of other.
236     """
237 jgs 123 a=_matchShape([self,other])
238     return Power_Symbol(a[0],a[1])
239    
240     def __rpow__(self,other):
241 jgs 149 """
242     Raises other to the symbol self.
243     """
244 jgs 123 a=_matchShape([self,other])
245     return Power_Symbol(a[1],a[0])
246    
247     def __mul__(self,other):
248 jgs 149 """
249     Multiplies other by symbol self. if _testForZero(other) 0 is returned.
250     """
251 jgs 123 if _testForZero(other):
252     return 0
253     else:
254     a=_matchShape([self,other])
255     return Mult_Symbol(a[0],a[1])
256    
257     def __rmul__(self,other):
258 jgs 149 """
259     Multiplies other by symbol self. if _testSForZero(other) 0 is returned.
260     """
261 jgs 123 return self*other
262    
263     def __getitem__(self,sl):
264     print sl
265    
266 jgs 149 class Float_Symbol(Symbol):
267 jgs 123 def __init__(self,name="symbol",shape=(),args=[]):
268     Symbol.__init__(self,dim=0,name="symbol",shape=(),args=[])
269    
270     class ScalarSymbol(Symbol):
271 jgs 149 """
272     A scalar symbol.
273     """
274 jgs 123 def __init__(self,dim=3,name="scalar"):
275 jgs 149 """
276     Creates a scalar symbol of spatial dimension dim.
277     """
278 jgs 123 if hasattr(dim,"getDim"):
279     d=dim.getDim()
280     else:
281     d=dim
282     Symbol.__init__(self,shape=(),dim=d,name=name)
283    
284     class VectorSymbol(Symbol):
285 jgs 149 """
286     A vector symbol.
287     """
288 jgs 123 def __init__(self,dim=3,name="vector"):
289 jgs 149 """
290     Creates a vector symbol of spatial dimension dim.
291     """
292 jgs 123 if hasattr(dim,"getDim"):
293     d=dim.getDim()
294     else:
295     d=dim
296     Symbol.__init__(self,shape=(d,),dim=d,name=name)
297    
298     class TensorSymbol(Symbol):
299 jgs 149 """
300     A tensor symbol.
301     """
302 jgs 123 def __init__(self,dim=3,name="tensor"):
303 jgs 149 """
304     Creates a tensor symbol of spatial dimension dim.
305     """
306 jgs 123 if hasattr(dim,"getDim"):
307     d=dim.getDim()
308     else:
309     d=dim
310     Symbol.__init__(self,shape=(d,d),dim=d,name=name)
311    
312     class Tensor3Symbol(Symbol):
313 jgs 149 """
314     A tensor order 3 symbol.
315     """
316 jgs 123 def __init__(self,dim=3,name="tensor3"):
317 jgs 149 """
318     Creates a tensor order 3 symbol of spatial dimension dim.
319     """
320 jgs 123 if hasattr(dim,"getDim"):
321     d=dim.getDim()
322     else:
323     d=dim
324     Symbol.__init__(self,shape=(d,d,d),dim=d,name=name)
325    
326     class Tensor4Symbol(Symbol):
327 jgs 149 """
328     A tensor order 4 symbol.
329     """
330 jgs 123 def __init__(self,dim=3,name="tensor4"):
331 jgs 149 """
332     Creates a tensor order 4 symbol of spatial dimension dim.
333     """
334 jgs 123 if hasattr(dim,"getDim"):
335     d=dim.getDim()
336     else:
337     d=dim
338     Symbol.__init__(self,shape=(d,d,d,d),dim=d,name=name)
339    
340     class Add_Symbol(Symbol):
341 jgs 149 """
342     Symbol representing the sum of two arguments.
343     """
344 jgs 123 def __init__(self,arg0,arg1):
345     a=[arg0,arg1]
346     Symbol.__init__(self,dim=_extractDim(a),shape=_extractShape(a),args=a)
347     def __str__(self):
348     return "(%s+%s)"%(str(self.getArgument(0)),str(self.getArgument(1)))
349     def eval(self,argval):
350     a=self.getEvaluatedArguments(argval)
351     return a[0]+a[1]
352     def _diff(self,arg):
353     a=self.getDifferentiatedArguments(arg)
354     return a[0]+a[1]
355    
356     class Mult_Symbol(Symbol):
357 jgs 149 """
358     Symbol representing the product of two arguments.
359     """
360 jgs 123 def __init__(self,arg0,arg1):
361     a=[arg0,arg1]
362     Symbol.__init__(self,dim=_extractDim(a),shape=_extractShape(a),args=a)
363     def __str__(self):
364     return "(%s*%s)"%(str(self.getArgument(0)),str(self.getArgument(1)))
365     def eval(self,argval):
366     a=self.getEvaluatedArguments(argval)
367     return a[0]*a[1]
368     def _diff(self,arg):
369     a=self.getDifferentiatedArguments(arg)
370     return self.getArgument(1)*a[0]+self.getArgument(0)*a[1]
371    
372     class Div_Symbol(Symbol):
373 jgs 149 """
374     Symbol representing the quotient of two arguments.
375     """
376 jgs 123 def __init__(self,arg0,arg1):
377     a=[arg0,arg1]
378     Symbol.__init__(self,dim=_extractDim(a),shape=_extractShape(a),args=a)
379     def __str__(self):
380     return "(%s/%s)"%(str(self.getArgument(0)),str(self.getArgument(1)))
381     def eval(self,argval):
382     a=self.getEvaluatedArguments(argval)
383     return a[0]/a[1]
384     def _diff(self,arg):
385     a=self.getDifferentiatedArguments(arg)
386     return (a[0]*self.getArgument(1)-self.getArgument(0)*a[1])/ \
387     (self.getArgument(1)*self.getArgument(1))
388    
389     class Power_Symbol(Symbol):
390 jgs 149 """
391     Symbol representing the power of the first argument to the power of the
392     second argument.
393     """
394 jgs 123 def __init__(self,arg0,arg1):
395     a=[arg0,arg1]
396     Symbol.__init__(self,dim=_extractDim(a),shape=_extractShape(a),args=a)
397     def __str__(self):
398     return "(%s**%s)"%(str(self.getArgument(0)),str(self.getArgument(1)))
399     def eval(self,argval):
400     a=self.getEvaluatedArguments(argval)
401     return a[0]**a[1]
402     def _diff(self,arg):
403     a=self.getDifferentiatedArguments(arg)
404     return self*(a[1]*log(self.getArgument(0))+self.getArgument(1)/self.getArgument(0)*a[0])
405    
406     class Abs_Symbol(Symbol):
407 jgs 149 """
408     Symbol representing absolute value of its argument.
409     """
410 jgs 123 def __init__(self,arg):
411     Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])
412     def __str__(self):
413     return "abs(%s)"%str(self.getArgument(0))
414     def eval(self,argval):
415     return abs(self.getEvaluatedArguments(argval)[0])
416     def _diff(self,arg):
417     return sign(self.getArgument(0))*self.getDifferentiatedArguments(arg)[0]
418    
419     #=========================================================
420     # some little helpers
421     #=========================================================
422     def _testForZero(arg):
423 jgs 149 """
424     Returns True is arg is considered to be zero.
425     """
426 jgs 123 if isinstance(arg,int):
427     return not arg>0
428     elif isinstance(arg,float):
429     return not arg>0.
430     elif isinstance(arg,numarray.NumArray):
431     a=abs(arg)
432     while isinstance(a,numarray.NumArray): a=numarray.sometrue(a)
433     return not a>0
434     else:
435     return False
436    
437     def _extractDim(args):
438     dim=None
439     for a in args:
440     if hasattr(a,"getDim"):
441     d=a.getDim()
442     if dim==None:
443     dim=d
444     else:
445     if dim!=d: raise ValueError,"inconsistent spatial dimension of arguments"
446     if dim==None:
447     raise ValueError,"cannot recover spatial dimension"
448     return dim
449    
450     def _identifyShape(arg):
451 jgs 149 """
452     Identifies the shape of arg.
453     """
454 jgs 123 if hasattr(arg,"getShape"):
455     arg_shape=arg.getShape()
456     elif hasattr(arg,"shape"):
457     s=arg.shape
458     if callable(s):
459     arg_shape=s()
460     else:
461     arg_shape=s
462     else:
463     arg_shape=()
464     return arg_shape
465    
466     def _extractShape(args):
467 jgs 149 """
468     Extracts the common shape of the list of arguments args.
469     """
470 jgs 123 shape=None
471     for a in args:
472     a_shape=_identifyShape(a)
473     if shape==None: shape=a_shape
474     if shape!=a_shape: raise ValueError,"inconsistent shape"
475     if shape==None:
476     raise ValueError,"cannot recover shape"
477     return shape
478    
479     def _matchShape(args,shape=None):
480 jgs 149 """
481     Returns the list of arguments args as object which have all the
482     specified shape.
483    
484     If shape is not given the shape "largest" shape of args is used.
485     """
486 jgs 123 # identify the list of shapes:
487     arg_shapes=[]
488     for a in args: arg_shapes.append(_identifyShape(a))
489     # get the largest shape (currently the longest shape):
490     if shape==None: shape=max(arg_shapes)
491    
492     out=[]
493     for i in range(len(args)):
494     if shape==arg_shapes[i]:
495     out.append(args[i])
496     else:
497     if len(shape)==0: # then len(arg_shapes[i])>0
498     raise ValueError,"cannot adopt shape of %s to %s"%(str(args[i]),str(shape))
499     else:
500     if len(arg_shapes[i])==0:
501     out.append(outer(args[i],numarray.ones(shape)))
502     else:
503     raise ValueError,"cannot adopt shape of %s to %s"%(str(args[i]),str(shape))
504     return out
505 jgs 124
506 jgs 123 #=========================================================
507 jgs 124 # wrappers for various mathematical functions:
508 jgs 123 #=========================================================
509     def diff(arg,dep):
510 jgs 149 """
511     Returns the derivative of arg with respect to dep.
512    
513     If arg is not Symbol object 0 is returned.
514     """
515 jgs 123 if isinstance(arg,Symbol):
516     return arg.diff(dep)
517     elif hasattr(arg,"shape"):
518     if callable(arg.shape):
519     return numarray.zeros(arg.shape(),numarray.Float)
520     else:
521     return numarray.zeros(arg.shape,numarray.Float)
522     else:
523     return 0
524    
525     def exp(arg):
526 jgs 82 """
527 jgs 149 Applies the exponential function to arg.
528    
529     @param arg: argument
530 jgs 123 """
531     if isinstance(arg,Symbol):
532     return Exp_Symbol(arg)
533     elif hasattr(arg,"exp"):
534     return arg.exp()
535 jgs 108 else:
536 jgs 123 return numarray.exp(arg)
537 jgs 82
538 jgs 123 class Exp_Symbol(Symbol):
539 jgs 149 """
540     Symbol representing the power of the first argument to the power of the
541     second argument.
542     """
543 jgs 123 def __init__(self,arg):
544     Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])
545     def __str__(self):
546     return "exp(%s)"%str(self.getArgument(0))
547     def eval(self,argval):
548     return exp(self.getEvaluatedArguments(argval)[0])
549     def _diff(self,arg):
550     return self*self.getDifferentiatedArguments(arg)[0]
551    
552     def sqrt(arg):
553 jgs 82 """
554 jgs 149 Applies the squre root function to arg.
555    
556     @param arg: argument
557 jgs 123 """
558     if isinstance(arg,Symbol):
559     return Sqrt_Symbol(arg)
560     elif hasattr(arg,"sqrt"):
561     return arg.sqrt()
562     else:
563     return numarray.sqrt(arg)
564 jgs 82
565 jgs 123 class Sqrt_Symbol(Symbol):
566 jgs 149 """
567     Symbol representing square root of argument.
568     """
569 jgs 123 def __init__(self,arg):
570     Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])
571     def __str__(self):
572     return "sqrt(%s)"%str(self.getArgument(0))
573     def eval(self,argval):
574     return sqrt(self.getEvaluatedArguments(argval)[0])
575     def _diff(self,arg):
576     return (-0.5)/self*self.getDifferentiatedArguments(arg)[0]
577    
578     def log(arg):
579 jgs 82 """
580 jgs 149 Applies the logarithmic function bases exp(1.) to arg
581    
582     @param arg: argument
583 jgs 123 """
584     if isinstance(arg,Symbol):
585     return Log_Symbol(arg)
586     elif hasattr(arg,"log"):
587     return arg.log()
588 jgs 108 else:
589 jgs 123 return numarray.log(arg)
590 jgs 82
591 jgs 123 class Log_Symbol(Symbol):
592 jgs 149 """
593     Symbol representing logarithm of the argument.
594     """
595 jgs 123 def __init__(self,arg):
596     Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])
597     def __str__(self):
598     return "log(%s)"%str(self.getArgument(0))
599     def eval(self,argval):
600     return log(self.getEvaluatedArguments(argval)[0])
601     def _diff(self,arg):
602     return self.getDifferentiatedArguments(arg)[0]/self.getArgument(0)
603 jgs 102
604 jgs 124 def ln(arg):
605     """
606 jgs 149 Applies the natural logarithmic function to arg.
607    
608     @param arg: argument
609 jgs 124 """
610     if isinstance(arg,Symbol):
611     return Ln_Symbol(arg)
612     elif hasattr(arg,"ln"):
613     return arg.log()
614     else:
615     return numarray.log(arg)
616    
617     class Ln_Symbol(Symbol):
618 jgs 149 """
619     Symbol representing natural logarithm of the argument.
620     """
621 jgs 124 def __init__(self,arg):
622     Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])
623     def __str__(self):
624     return "ln(%s)"%str(self.getArgument(0))
625     def eval(self,argval):
626     return ln(self.getEvaluatedArguments(argval)[0])
627     def _diff(self,arg):
628     return self.getDifferentiatedArguments(arg)[0]/self.getArgument(0)
629    
630 jgs 123 def sin(arg):
631 jgs 82 """
632 jgs 149 Applies the sin function to arg.
633    
634     @param arg: argument
635 jgs 123 """
636     if isinstance(arg,Symbol):
637     return Sin_Symbol(arg)
638     elif hasattr(arg,"sin"):
639     return arg.sin()
640     else:
641     return numarray.sin(arg)
642 jgs 82
643 jgs 123 class Sin_Symbol(Symbol):
644 jgs 149 """
645     Symbol representing sin of the argument.
646     """
647 jgs 123 def __init__(self,arg):
648     Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])
649     def __str__(self):
650     return "sin(%s)"%str(self.getArgument(0))
651     def eval(self,argval):
652     return sin(self.getEvaluatedArguments(argval)[0])
653     def _diff(self,arg):
654     return cos(self.getArgument(0))*self.getDifferentiatedArguments(arg)[0]
655    
656     def cos(arg):
657 jgs 82 """
658 jgs 149 Applies the cos function to arg.
659    
660     @param arg: argument
661 jgs 123 """
662     if isinstance(arg,Symbol):
663     return Cos_Symbol(arg)
664     elif hasattr(arg,"cos"):
665     return arg.cos()
666 jgs 82 else:
667 jgs 123 return numarray.cos(arg)
668 jgs 82
669 jgs 123 class Cos_Symbol(Symbol):
670 jgs 149 """
671     Symbol representing cos of the argument.
672     """
673 jgs 123 def __init__(self,arg):
674     Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])
675     def __str__(self):
676     return "cos(%s)"%str(self.getArgument(0))
677     def eval(self,argval):
678     return cos(self.getEvaluatedArguments(argval)[0])
679     def _diff(self,arg):
680     return -sin(self.getArgument(0))*self.getDifferentiatedArguments(arg)[0]
681    
682     def tan(arg):
683 jgs 82 """
684 jgs 149 Applies the tan function to arg.
685    
686     @param arg: argument
687 jgs 123 """
688     if isinstance(arg,Symbol):
689     return Tan_Symbol(arg)
690     elif hasattr(arg,"tan"):
691     return arg.tan()
692     else:
693     return numarray.tan(arg)
694 jgs 82
695 jgs 123 class Tan_Symbol(Symbol):
696 jgs 149 """
697     Symbol representing tan of the argument.
698     """
699 jgs 123 def __init__(self,arg):
700     Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])
701     def __str__(self):
702     return "tan(%s)"%str(self.getArgument(0))
703     def eval(self,argval):
704     return tan(self.getEvaluatedArguments(argval)[0])
705     def _diff(self,arg):
706     s=cos(self.getArgument(0))
707     return 1./(s*s)*self.getDifferentiatedArguments(arg)[0]
708    
709     def sign(arg):
710 jgs 82 """
711 jgs 149 Applies the sign function to arg.
712    
713     @param arg: argument
714 jgs 123 """
715     if isinstance(arg,Symbol):
716     return Sign_Symbol(arg)
717     elif hasattr(arg,"sign"):
718     return arg.sign()
719 jgs 82 else:
720 jgs 123 return numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))- \
721     numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))
722 jgs 82
723 jgs 123 class Sign_Symbol(Symbol):
724 jgs 149 """
725     Symbol representing the sign of the argument.
726     """
727 jgs 123 def __init__(self,arg):
728     Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])
729     def __str__(self):
730     return "sign(%s)"%str(self.getArgument(0))
731     def eval(self,argval):
732     return sign(self.getEvaluatedArguments(argval)[0])
733    
734     def maxval(arg):
735 jgs 82 """
736 jgs 149 Returns the maximum value of argument arg.
737    
738     @param arg: argument
739 jgs 123 """
740     if isinstance(arg,Symbol):
741     return Max_Symbol(arg)
742     elif hasattr(arg,"maxval"):
743     return arg.maxval()
744     elif hasattr(arg,"max"):
745     return arg.max()
746     else:
747     return arg
748 jgs 82
749 jgs 123 class Max_Symbol(Symbol):
750 jgs 149 """
751     Symbol representing the maximum value of the argument.
752     """
753 jgs 123 def __init__(self,arg):
754     Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])
755     def __str__(self):
756     return "maxval(%s)"%str(self.getArgument(0))
757     def eval(self,argval):
758     return maxval(self.getEvaluatedArguments(argval)[0])
759    
760     def minval(arg):
761 jgs 82 """
762 jgs 149 Returns the minimum value of argument arg.
763    
764     @param arg: argument
765 jgs 123 """
766     if isinstance(arg,Symbol):
767     return Min_Symbol(arg)
768     elif hasattr(arg,"maxval"):
769     return arg.minval()
770     elif hasattr(arg,"min"):
771     return arg.min()
772 jgs 82 else:
773 jgs 123 return arg
774 jgs 82
775 jgs 123 class Min_Symbol(Symbol):
776 jgs 149 """
777     Symbol representing the minimum value of the argument.
778     """
779 jgs 123 def __init__(self,arg):
780     Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])
781     def __str__(self):
782     return "minval(%s)"%str(self.getArgument(0))
783     def eval(self,argval):
784     return minval(self.getEvaluatedArguments(argval)[0])
785    
786     def wherePositive(arg):
787 jgs 82 """
788 jgs 149 Returns the positive values of argument arg.
789    
790     @param arg: argument
791 jgs 123 """
792     if _testForZero(arg):
793     return 0
794     elif isinstance(arg,Symbol):
795     return WherePositive_Symbol(arg)
796     elif hasattr(arg,"wherePositive"):
797     return arg.minval()
798     elif hasattr(arg,"wherePositive"):
799     numarray.greater(arg,numarray.zeros(arg.shape,numarray.Float))
800     else:
801     if arg>0:
802     return 1.
803     else:
804     return 0.
805 jgs 82
806 jgs 123 class WherePositive_Symbol(Symbol):
807 jgs 149 """
808     Symbol representing the wherePositive function.
809     """
810 jgs 123 def __init__(self,arg):
811     Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])
812     def __str__(self):
813     return "wherePositive(%s)"%str(self.getArgument(0))
814     def eval(self,argval):
815     return wherePositive(self.getEvaluatedArguments(argval)[0])
816    
817     def whereNegative(arg):
818 jgs 82 """
819 jgs 149 Returns the negative values of argument arg.
820    
821     @param arg: argument
822 jgs 123 """
823     if _testForZero(arg):
824     return 0
825     elif isinstance(arg,Symbol):
826     return WhereNegative_Symbol(arg)
827     elif hasattr(arg,"whereNegative"):
828     return arg.whereNegative()
829     elif hasattr(arg,"shape"):
830     numarray.less(arg,numarray.zeros(arg.shape,numarray.Float))
831 jgs 88 else:
832 jgs 123 if arg<0:
833     return 1.
834     else:
835     return 0.
836 jgs 82
837 jgs 123 class WhereNegative_Symbol(Symbol):
838 jgs 149 """
839     Symbol representing the whereNegative function.
840     """
841 jgs 123 def __init__(self,arg):
842     Symbol.__init__(self,shape=arg.getShape(),dim=arg.getDim(),args=[arg])
843     def __str__(self):
844     return "whereNegative(%s)"%str(self.getArgument(0))
845     def eval(self,argval):
846     return whereNegative(self.getEvaluatedArguments(argval)[0])
847    
848 jgs 147 def maximum(arg0,arg1):
849 jgs 149 """
850     Return arg1 where arg1 is bigger then arg0 otherwise arg0 is returned.
851     """
852 jgs 147 m=whereNegative(arg0-arg1)
853     return m*arg1+(1.-m)*arg0
854    
855     def minimum(arg0,arg1):
856 jgs 149 """
857     Return arg0 where arg1 is bigger then arg0 otherwise arg1 is returned.
858     """
859 jgs 147 m=whereNegative(arg0-arg1)
860     return m*arg0+(1.-m)*arg1
861    
862 jgs 123 def outer(arg0,arg1):
863     if _testForZero(arg0) or _testForZero(arg1):
864     return 0
865     else:
866     if isinstance(arg0,Symbol) or isinstance(arg1,Symbol):
867     return Outer_Symbol(arg0,arg1)
868     elif _identifyShape(arg0)==() or _identifyShape(arg1)==():
869     return arg0*arg1
870     elif isinstance(arg0,numarray.NumArray) and isinstance(arg1,numarray.NumArray):
871     return numarray.outer(arg0,arg1)
872     else:
873     if arg0.getRank()==1 and arg1.getRank()==1:
874     out=escript.Data(0,(arg0.getShape()[0],arg1.getShape()[0]),arg1.getFunctionSpace())
875     for i in range(arg0.getShape()[0]):
876     for j in range(arg1.getShape()[0]):
877     out[i,j]=arg0[i]*arg1[j]
878     return out
879     else:
880     raise ValueError,"outer is not fully implemented yet."
881    
882     class Outer_Symbol(Symbol):
883 jgs 149 """
884     Symbol representing the outer product of its two arguments.
885     """
886 jgs 123 def __init__(self,arg0,arg1):
887     a=[arg0,arg1]
888     s=tuple(list(_identifyShape(arg0))+list(_identifyShape(arg1)))
889     Symbol.__init__(self,shape=s,dim=_extractDim(a),args=a)
890     def __str__(self):
891     return "outer(%s,%s)"%(str(self.getArgument(0)),str(self.getArgument(1)))
892     def eval(self,argval):
893     a=self.getEvaluatedArguments(argval)
894     return outer(a[0],a[1])
895     def _diff(self,arg):
896     a=self.getDifferentiatedArguments(arg)
897     return outer(a[0],self.getArgument(1))+outer(self.getArgument(0),a[1])
898    
899     def interpolate(arg,where):
900 jgs 82 """
901 jgs 149 Interpolates the function into the FunctionSpace where.
902 jgs 82
903 jgs 149 @param arg: interpolant
904     @param where: FunctionSpace to interpolate to
905 jgs 82 """
906 jgs 123 if _testForZero(arg):
907     return 0
908     elif isinstance(arg,Symbol):
909     return Interpolated_Symbol(arg,where)
910 jgs 82 else:
911 jgs 123 return escript.Data(arg,where)
912 jgs 82
913 jgs 149 class Interpolated_Symbol(Symbol):
914     """
915     Symbol representing the integral of the argument.
916     """
917 jgs 123 def __init__(self,arg,where):
918     Symbol.__init__(self,shape=_extractShape(arg),dim=_extractDim([arg]),args=[arg,where])
919     def __str__(self):
920     return "interpolated(%s)"%(str(self.getArgument(0)))
921     def eval(self,argval):
922     a=self.getEvaluatedArguments(argval)
923     return integrate(a[0],where=self.getArgument(1))
924     def _diff(self,arg):
925     a=self.getDifferentiatedArguments(arg)
926     return integrate(a[0],where=self.getArgument(1))
927    
928 jgs 147 def div(arg,where=None):
929     """
930 jgs 149 Returns the divergence of arg at where.
931 jgs 147
932 jgs 149 @param arg: Data object representing the function which gradient to
933     be calculated.
934     @param where: FunctionSpace in which the gradient will be calculated.
935     If not present or C{None} an appropriate default is used.
936 jgs 147 """
937     return trace(grad(arg,where))
938    
939 jgs 149 def jump(arg):
940     """
941     Returns the jump of arg across a continuity.
942    
943     @param arg: Data object representing the function which gradient
944     to be calculated.
945     """
946     d=arg.getDomain()
947     return arg.interpolate(escript.FunctionOnContactOne())-arg.interpolate(escript.FunctionOnContactZero())
948    
949    
950 jgs 123 def grad(arg,where=None):
951 jgs 102 """
952 jgs 149 Returns the spatial gradient of arg at where.
953 jgs 102
954 jgs 149 @param arg: Data object representing the function which gradient
955     to be calculated.
956     @param where: FunctionSpace in which the gradient will be calculated.
957     If not present or C{None} an appropriate default is used.
958 jgs 102 """
959 jgs 123 if _testForZero(arg):
960     return 0
961     elif isinstance(arg,Symbol):
962     return Grad_Symbol(arg,where)
963     elif hasattr(arg,"grad"):
964     if where==None:
965     return arg.grad()
966     else:
967     return arg.grad(where)
968 jgs 102 else:
969 jgs 123 return arg*0.
970 jgs 102
971 jgs 149 class Grad_Symbol(Symbol):
972     """
973     Symbol representing the gradient of the argument.
974     """
975 jgs 123 def __init__(self,arg,where=None):
976     d=_extractDim([arg])
977     s=tuple(list(_identifyShape([arg])).append(d))
978     Symbol.__init__(self,shape=s,dim=_extractDim([arg]),args=[arg,where])
979     def __str__(self):
980     return "grad(%s)"%(str(self.getArgument(0)))
981     def eval(self,argval):
982     a=self.getEvaluatedArguments(argval)
983     return grad(a[0],where=self.getArgument(1))
984     def _diff(self,arg):
985     a=self.getDifferentiatedArguments(arg)
986     return grad(a[0],where=self.getArgument(1))
987    
988     def integrate(arg,where=None):
989 jgs 82 """
990 jgs 149 Return the integral if the function represented by Data object arg over
991     its domain.
992 jgs 82
993 jgs 123 @param arg: Data object representing the function which is integrated.
994 jgs 149 @param where: FunctionSpace in which the integral is calculated.
995     If not present or C{None} an appropriate default is used.
996 jgs 82 """
997 jgs 123 if _testForZero(arg):
998     return 0
999     elif isinstance(arg,Symbol):
1000     return Integral_Symbol(arg,where)
1001     else:
1002     if not where==None: arg=escript.Data(arg,where)
1003     if arg.getRank()==0:
1004     return arg.integrate()[0]
1005     else:
1006     return arg.integrate()
1007 jgs 82
1008 jgs 149 class Integral_Symbol(Float_Symbol):
1009     """
1010     Symbol representing the integral of the argument.
1011     """
1012 jgs 123 def __init__(self,arg,where=None):
1013     Float_Symbol.__init__(self,shape=_identifyShape([arg]),args=[arg,where])
1014     def __str__(self):
1015     return "integral(%s)"%(str(self.getArgument(0)))
1016     def eval(self,argval):
1017     a=self.getEvaluatedArguments(argval)
1018     return integrate(a[0],where=self.getArgument(1))
1019     def _diff(self,arg):
1020     a=self.getDifferentiatedArguments(arg)
1021     return integrate(a[0],where=self.getArgument(1))
1022    
1023     #=============================
1024     #
1025     # wrapper for various functions: if the argument has attribute the function name
1026 jgs 124 # as an argument it calls the corresponding methods. Otherwise the corresponding
1027     # numarray function is called.
1028    
1029 jgs 123 # functions involving the underlying Domain:
1030    
1031    
1032     # functions returning Data objects:
1033    
1034     def transpose(arg,axis=None):
1035 jgs 82 """
1036 jgs 149 Returns the transpose of the Data object arg.
1037 jgs 82
1038 jgs 149 @param arg:
1039 jgs 82 """
1040 jgs 123 if axis==None:
1041     r=0
1042     if hasattr(arg,"getRank"): r=arg.getRank()
1043     if hasattr(arg,"rank"): r=arg.rank
1044     axis=r/2
1045     if isinstance(arg,Symbol):
1046     return Transpose_Symbol(arg,axis=r)
1047 jgs 102 if isinstance(arg,escript.Data):
1048 jgs 123 # hack for transpose
1049     r=arg.getRank()
1050     if r!=2: raise ValueError,"Tranpose only avalaible for rank 2 objects"
1051     s=arg.getShape()
1052     out=escript.Data(0.,(s[1],s[0]),arg.getFunctionSpace())
1053     for i in range(s[0]):
1054     for j in range(s[1]):
1055     out[j,i]=arg[i,j]
1056     return out
1057     # end hack for transpose
1058     return arg.transpose(axis)
1059 jgs 102 else:
1060 jgs 123 return numarray.transpose(arg,axis=axis)
1061 jgs 82
1062 jgs 123 def trace(arg,axis0=0,axis1=1):
1063 jgs 82 """
1064 jgs 149 Return
1065 jgs 82
1066 jgs 149 @param arg:
1067 jgs 82 """
1068 jgs 123 if isinstance(arg,Symbol):
1069     s=list(arg.getShape())
1070     s=tuple(s[0:axis0]+s[axis0+1:axis1]+s[axis1+1:])
1071     return Trace_Symbol(arg,axis0=axis0,axis1=axis1)
1072     elif isinstance(arg,escript.Data):
1073     # hack for trace
1074     s=arg.getShape()
1075     if s[axis0]!=s[axis1]:
1076     raise ValueError,"illegal axis in trace"
1077     out=escript.Scalar(0.,arg.getFunctionSpace())
1078 jgs 147 for i in range(s[axis0]):
1079     out+=arg[i,i]
1080 jgs 123 return out
1081 jgs 126 # end hack for trace
1082 jgs 102 else:
1083 jgs 123 return numarray.trace(arg,axis0=axis0,axis1=axis1)
1084 jgs 82
1085 jgs 123 def Trace_Symbol(Symbol):
1086     pass
1087    
1088 jgs 102 def length(arg):
1089     """
1090    
1091 jgs 149 @param arg:
1092 jgs 102 """
1093     if isinstance(arg,escript.Data):
1094 jgs 108 if arg.isEmpty(): return escript.Data()
1095     if arg.getRank()==0:
1096     return abs(arg)
1097     elif arg.getRank()==1:
1098 jgs 147 out=escript.Scalar(0,arg.getFunctionSpace())
1099 jgs 104 for i in range(arg.getShape()[0]):
1100 jgs 147 out+=arg[i]**2
1101     return sqrt(out)
1102 jgs 108 elif arg.getRank()==2:
1103 jgs 147 out=escript.Scalar(0,arg.getFunctionSpace())
1104 jgs 108 for i in range(arg.getShape()[0]):
1105     for j in range(arg.getShape()[1]):
1106 jgs 147 out+=arg[i,j]**2
1107     return sqrt(out)
1108 jgs 108 elif arg.getRank()==3:
1109 jgs 147 out=escript.Scalar(0,arg.getFunctionSpace())
1110 jgs 108 for i in range(arg.getShape()[0]):
1111     for j in range(arg.getShape()[1]):
1112     for k in range(arg.getShape()[2]):
1113 jgs 147 out+=arg[i,j,k]**2
1114     return sqrt(out)
1115 jgs 108 elif arg.getRank()==4:
1116 jgs 147 out=escript.Scalar(0,arg.getFunctionSpace())
1117 jgs 108 for i in range(arg.getShape()[0]):
1118     for j in range(arg.getShape()[1]):
1119     for k in range(arg.getShape()[2]):
1120     for l in range(arg.getShape()[3]):
1121 jgs 147 out+=arg[i,j,k,l]**2
1122     return sqrt(out)
1123 jgs 104 else:
1124 jgs 126 raise SystemError,"length is not been fully implemented yet"
1125     # return arg.length()
1126 jgs 147 elif isinstance(arg,float):
1127     return abs(arg)
1128 jgs 102 else:
1129     return sqrt((arg**2).sum())
1130    
1131 jgs 113 def deviator(arg):
1132     """
1133 jgs 149 @param arg:
1134 jgs 113 """
1135     if isinstance(arg,escript.Data):
1136     shape=arg.getShape()
1137     else:
1138     shape=arg.shape
1139     if len(shape)!=2:
1140     raise ValueError,"Deviator requires rank 2 object"
1141     if shape[0]!=shape[1]:
1142     raise ValueError,"Deviator requires a square matrix"
1143     return arg-1./(shape[0]*1.)*trace(arg)*kronecker(shape[0])
1144    
1145 jgs 123 def inner(arg0,arg1):
1146 jgs 113 """
1147 jgs 149 @param arg0:
1148     @param arg1:
1149 jgs 113 """
1150 jgs 147 if isinstance(arg0,escript.Data):
1151     arg=arg0
1152     else:
1153     arg=arg1
1154    
1155     out=escript.Scalar(0,arg.getFunctionSpace())
1156     if arg.getRank()==0:
1157 jgs 123 return arg0*arg1
1158 jgs 147 elif arg.getRank()==1:
1159     out=escript.Scalar(0,arg.getFunctionSpace())
1160     for i in range(arg.getShape()[0]):
1161     out+=arg0[i]*arg1[i]
1162     elif arg.getRank()==2:
1163     out=escript.Scalar(0,arg.getFunctionSpace())
1164     for i in range(arg.getShape()[0]):
1165     for j in range(arg.getShape()[1]):
1166     out+=arg0[i,j]*arg1[i,j]
1167     elif arg.getRank()==3:
1168     out=escript.Scalar(0,arg.getFunctionSpace())
1169     for i in range(arg.getShape()[0]):
1170     for j in range(arg.getShape()[1]):
1171     for k in range(arg.getShape()[2]):
1172     out+=arg0[i,j,k]*arg1[i,j,k]
1173     elif arg.getRank()==4:
1174     out=escript.Scalar(0,arg.getFunctionSpace())
1175     for i in range(arg.getShape()[0]):
1176     for j in range(arg.getShape()[1]):
1177     for k in range(arg.getShape()[2]):
1178     for l in range(arg.getShape()[3]):
1179     out+=arg0[i,j,k,l]*arg1[i,j,k,l]
1180 jgs 113 else:
1181     raise SystemError,"inner is not been implemented yet"
1182 jgs 147 return out
1183 jgs 113
1184 jgs 149 def tensormult(arg0,arg1):
1185     # check LinearPDE!!!!
1186     raise SystemError,"tensormult is not implemented yet!"
1187    
1188 jgs 123 def matrixmult(arg0,arg1):
1189 jgs 102
1190 jgs 123 if isinstance(arg1,numarray.NumArray) and isinstance(arg0,numarray.NumArray):
1191     numarray.matrixmult(arg0,arg1)
1192 jgs 102 else:
1193 jgs 123 # escript.matmult(arg0,arg1)
1194     if isinstance(arg1,escript.Data) and not isinstance(arg0,escript.Data):
1195     arg0=escript.Data(arg0,arg1.getFunctionSpace())
1196     elif isinstance(arg0,escript.Data) and not isinstance(arg1,escript.Data):
1197     arg1=escript.Data(arg1,arg0.getFunctionSpace())
1198     if arg0.getRank()==2 and arg1.getRank()==1:
1199     out=escript.Data(0,(arg0.getShape()[0],),arg0.getFunctionSpace())
1200     for i in range(arg0.getShape()[0]):
1201     for j in range(arg0.getShape()[1]):
1202 jgs 126 # uses Data object slicing, plus Data * and += operators
1203 jgs 123 out[i]+=arg0[i,j]*arg1[j]
1204     return out
1205 jgs 147 elif arg0.getRank()==1 and arg1.getRank()==1:
1206     return inner(arg0,arg1)
1207 jgs 123 else:
1208     raise SystemError,"matrixmult is not fully implemented yet!"
1209 jgs 124
1210 jgs 123 #=========================================================
1211 jgs 102 # reduction operations:
1212 jgs 123 #=========================================================
1213 jgs 102 def sum(arg):
1214     """
1215 jgs 149 @param arg:
1216 jgs 102 """
1217     return arg.sum()
1218    
1219 jgs 82 def sup(arg):
1220     """
1221 jgs 149 @param arg:
1222 jgs 82 """
1223 jgs 102 if isinstance(arg,escript.Data):
1224     return arg.sup()
1225 jgs 108 elif isinstance(arg,float) or isinstance(arg,int):
1226     return arg
1227 jgs 102 else:
1228     return arg.max()
1229 jgs 82
1230     def inf(arg):
1231     """
1232 jgs 149 @param arg:
1233 jgs 82 """
1234 jgs 102 if isinstance(arg,escript.Data):
1235     return arg.inf()
1236 jgs 108 elif isinstance(arg,float) or isinstance(arg,int):
1237     return arg
1238 jgs 102 else:
1239     return arg.min()
1240 jgs 82
1241 jgs 102 def L2(arg):
1242 jgs 82 """
1243 jgs 149 Returns the L2-norm of the argument
1244 jgs 82
1245 jgs 149 @param arg:
1246 jgs 82 """
1247 jgs 108 if isinstance(arg,escript.Data):
1248     return arg.L2()
1249     elif isinstance(arg,float) or isinstance(arg,int):
1250     return abs(arg)
1251     else:
1252     return numarry.sqrt(dot(arg,arg))
1253 jgs 82
1254 jgs 102 def Lsup(arg):
1255 jgs 82 """
1256 jgs 149 @param arg:
1257 jgs 82 """
1258 jgs 102 if isinstance(arg,escript.Data):
1259     return arg.Lsup()
1260 jgs 108 elif isinstance(arg,float) or isinstance(arg,int):
1261     return abs(arg)
1262 jgs 102 else:
1263 jgs 149 return numarray.abs(arg).max()
1264 jgs 82
1265 jgs 123 def dot(arg0,arg1):
1266 jgs 117 """
1267 jgs 149 @param arg0:
1268     @param arg1:
1269 jgs 117 """
1270 jgs 123 if isinstance(arg0,escript.Data):
1271     return arg0.dot(arg1)
1272 jgs 102 elif isinstance(arg1,escript.Data):
1273 jgs 123 return arg1.dot(arg0)
1274 jgs 102 else:
1275 jgs 123 return numarray.dot(arg0,arg1)
1276 jgs 113
1277     def kronecker(d):
1278 jgs 122 if hasattr(d,"getDim"):
1279 jgs 147 return numarray.identity(d.getDim())*1.
1280 jgs 122 else:
1281 jgs 147 return numarray.identity(d)*1.
1282 jgs 117
1283     def unit(i,d):
1284     """
1285 jgs 149 Return a unit vector of dimension d with nonzero index i.
1286    
1287     @param d: dimension
1288     @param i: index
1289 jgs 117 """
1290 jgs 123 e = numarray.zeros((d,),numarray.Float)
1291 jgs 117 e[i] = 1.0
1292     return e
1293 jgs 122
1294 jgs 123 # ============================================
1295     # testing
1296     # ============================================
1297    
1298     if __name__=="__main__":
1299     u=ScalarSymbol(dim=2,name="u")
1300     v=ScalarSymbol(dim=2,name="v")
1301     v=VectorSymbol(2,"v")
1302     u=VectorSymbol(2,"u")
1303    
1304     print u+5,(u+5).diff(u)
1305     print 5+u,(5+u).diff(u)
1306     print u+v,(u+v).diff(u)
1307     print v+u,(v+u).diff(u)
1308    
1309     print u*5,(u*5).diff(u)
1310     print 5*u,(5*u).diff(u)
1311     print u*v,(u*v).diff(u)
1312     print v*u,(v*u).diff(u)
1313    
1314     print u-5,(u-5).diff(u)
1315     print 5-u,(5-u).diff(u)
1316     print u-v,(u-v).diff(u)
1317     print v-u,(v-u).diff(u)
1318    
1319     print u/5,(u/5).diff(u)
1320     print 5/u,(5/u).diff(u)
1321     print u/v,(u/v).diff(u)
1322     print v/u,(v/u).diff(u)
1323    
1324     print u**5,(u**5).diff(u)
1325     print 5**u,(5**u).diff(u)
1326     print u**v,(u**v).diff(u)
1327     print v**u,(v**u).diff(u)
1328    
1329     print exp(u),exp(u).diff(u)
1330     print sqrt(u),sqrt(u).diff(u)
1331     print log(u),log(u).diff(u)
1332     print sin(u),sin(u).diff(u)
1333     print cos(u),cos(u).diff(u)
1334     print tan(u),tan(u).diff(u)
1335     print sign(u),sign(u).diff(u)
1336     print abs(u),abs(u).diff(u)
1337     print wherePositive(u),wherePositive(u).diff(u)
1338     print whereNegative(u),whereNegative(u).diff(u)
1339     print maxval(u),maxval(u).diff(u)
1340     print minval(u),minval(u).diff(u)
1341    
1342     g=grad(u)
1343     print diff(5*g,g)
1344     4*(g+transpose(g))/2+6*trace(g)*kronecker(3)
1345 jgs 124
1346 jgs 123 #
1347 jgs 122 # $Log$
1348 jgs 149 # Revision 1.17 2005/09/01 03:31:28 jgs
1349     # Merge of development branch dev-02 back to main trunk on 2005-09-01
1350     #
1351 jgs 148 # Revision 1.16 2005/08/23 01:24:28 jgs
1352     # Merge of development branch dev-02 back to main trunk on 2005-08-23
1353     #
1354 jgs 147 # Revision 1.15 2005/08/12 01:45:36 jgs
1355     # erge of development branch dev-02 back to main trunk on 2005-08-12
1356 jgs 142 #
1357 jgs 149 # Revision 1.14.2.8 2005/08/26 05:06:37 cochrane
1358     # Corrected errors in docstrings. Improved output formatting of docstrings.
1359     # Other minor improvements to code and docs (eg spelling etc).
1360     #
1361     # Revision 1.14.2.7 2005/08/26 04:45:40 cochrane
1362     # Fixed and tidied markup and docstrings. Some *Symbol classes were defined
1363     # as functions, so changed them to classes (hopefully this was the right thing
1364     # to do).
1365     #
1366     # Revision 1.14.2.6 2005/08/26 04:30:13 gross
1367     # gneric unit testing for linearPDE
1368     #
1369     # Revision 1.14.2.5 2005/08/24 02:02:52 gross
1370     # jump function added
1371     #
1372 jgs 148 # Revision 1.14.2.4 2005/08/18 04:39:32 gross
1373     # the constants have been removed from util.py as they not needed anymore. PDE related constants are accessed through LinearPDE attributes now
1374     #
1375 jgs 147 # Revision 1.14.2.3 2005/08/03 09:55:33 gross
1376     # ContactTest is passing now./mk install!
1377 jgs 126 #
1378 jgs 147 # Revision 1.14.2.2 2005/08/02 03:15:14 gross
1379     # bug inb trace fixed!
1380 jgs 124 #
1381 jgs 147 # Revision 1.14.2.1 2005/07/29 07:10:28 gross
1382     # new functions in util and a new pde type in linearPDEs
1383 jgs 123 #
1384 jgs 147 # Revision 1.2.2.21 2005/07/28 04:19:23 gross
1385     # new functions maximum and minimum introduced.
1386 jgs 122 #
1387 jgs 142 # Revision 1.2.2.20 2005/07/25 01:26:27 gross
1388     # bug in inner fixed
1389     #
1390 jgs 126 # Revision 1.2.2.19 2005/07/21 04:01:28 jgs
1391     # minor comment fixes
1392     #
1393     # Revision 1.2.2.18 2005/07/21 01:02:43 jgs
1394     # commit ln() updates to development branch version
1395     #
1396     # Revision 1.12 2005/07/20 06:14:58 jgs
1397     # added ln(data) style wrapper for data.ln() - also added corresponding
1398     # implementation of Ln_Symbol class (not sure if this is right though)
1399     #
1400     # Revision 1.11 2005/07/08 04:07:35 jgs
1401     # Merge of development branch back to main trunk on 2005-07-08
1402     #
1403     # Revision 1.10 2005/06/09 05:37:59 jgs
1404     # Merge of development branch back to main trunk on 2005-06-09
1405     #
1406 jgs 123 # Revision 1.2.2.17 2005/07/07 07:28:58 gross
1407     # some stuff added to util.py to improve functionality
1408     #
1409     # Revision 1.2.2.16 2005/06/30 01:53:55 gross
1410     # a bug in coloring fixed
1411     #
1412     # Revision 1.2.2.15 2005/06/29 02:36:43 gross
1413     # Symbols have been introduced and some function clarified. needs much more work
1414     #
1415 jgs 122 # Revision 1.2.2.14 2005/05/20 04:05:23 gross
1416     # some work on a darcy flow started
1417     #
1418     # Revision 1.2.2.13 2005/03/16 05:17:58 matt
1419     # Implemented unit(idx, dim) to create cartesian unit basis vectors to
1420     # complement kronecker(dim) function.
1421     #
1422     # Revision 1.2.2.12 2005/03/10 08:14:37 matt
1423     # Added non-member Linf utility function to complement Data::Linf().
1424     #
1425     # Revision 1.2.2.11 2005/02/17 05:53:25 gross
1426     # some bug in saveDX fixed: in fact the bug was in
1427     # DataC/getDataPointShape
1428     #
1429     # Revision 1.2.2.10 2005/01/11 04:59:36 gross
1430     # automatic interpolation in integrate switched off
1431     #
1432     # Revision 1.2.2.9 2005/01/11 03:38:13 gross
1433     # 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.
1434     #
1435     # Revision 1.2.2.8 2005/01/05 04:21:41 gross
1436     # FunctionSpace checking/matchig in slicing added
1437     #
1438     # Revision 1.2.2.7 2004/12/29 05:29:59 gross
1439     # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()
1440     #
1441     # Revision 1.2.2.6 2004/12/24 06:05:41 gross
1442     # some changes in linearPDEs to add AdevectivePDE
1443     #
1444     # Revision 1.2.2.5 2004/12/17 00:06:53 gross
1445     # mk sets ESYS_ROOT is undefined
1446     #
1447     # Revision 1.2.2.4 2004/12/07 03:19:51 gross
1448     # options for GMRES and PRES20 added
1449     #
1450     # Revision 1.2.2.3 2004/12/06 04:55:18 gross
1451     # function wraper extended
1452     #
1453     # Revision 1.2.2.2 2004/11/22 05:44:07 gross
1454     # a few more unitary functions have been added but not implemented in Data yet
1455     #
1456     # Revision 1.2.2.1 2004/11/12 06:58:15 gross
1457     # 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
1458     #
1459     # Revision 1.2 2004/10/27 00:23:36 jgs
1460     # fixed minor syntax error
1461     #
1462     # Revision 1.1.1.1 2004/10/26 06:53:56 jgs
1463     # initial import of project esys2
1464     #
1465     # Revision 1.1.2.3 2004/10/26 06:43:48 jgs
1466     # committing Lutz's and Paul's changes to brach jgs
1467     #
1468     # Revision 1.1.4.1 2004/10/20 05:32:51 cochrane
1469     # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.
1470     #
1471     # Revision 1.1 2004/08/05 03:58:27 gross
1472     # Bug in Assemble_NodeCoordinates fixed
1473     #
1474     #
1475 jgs 149
1476     # 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