/[escript]/branches/symbolic_from_3470/escript/py_src/util.py
ViewVC logotype

Contents of /branches/symbolic_from_3470/escript/py_src/util.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3496 - (show annotations)
Wed Apr 6 03:58:42 2011 UTC (7 years, 10 months ago) by caltinay
File MIME type: text/x-python
File size: 96654 byte(s)
Added erf(), eigenvectors() and eigenvectors_and_eigenvalues() to the list of
symbolic functions.

1
2 ########################################################
3 #
4 # Copyright (c) 2003-2010 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 """
23 Utility functions for escript
24
25 :var __author__: name of author
26 :var __copyright__: copyrights
27 :var __license__: licence agreement
28 :var __url__: url entry point on documentation
29 :var __version__: version
30 :var __date__: date of the version
31 :var EPSILON: smallest positive value with 1.<1.+EPSILON
32 :var DBLE_MAX: largest positive float
33 """
34
35 __author__="Lutz Gross, l.gross@uq.edu.au"
36
37
38 import math
39 import numpy
40 import escript
41 import os
42 from esys.escript import C_GeneralTensorProduct
43 from esys.escript import getVersion, getMPIRankWorld, getMPIWorldMax
44 from esys.escript import printParallelThreadCounts
45 from esys.escript import listEscriptParams
46 from esys.escript.escriptcpp import Data, _saveDataCSV, _condEval
47 from esys.escript.symbolic import functions as symfn
48 from esys.escript.symbolic import *
49
50 #=========================================================
51 # some helpers:
52 #=========================================================
53 def getEpsilon():
54 return escript.getMachinePrecision()
55 EPSILON=getEpsilon()
56
57 def getMaxFloat():
58 return escript.getMaxFloat()
59 DBLE_MAX=getMaxFloat()
60
61 def getTagNames(domain):
62 """
63 Returns a list of tag names used by the domain.
64
65 :param domain: a domain object
66 :type domain: `escript.Domain`
67 :return: a list of tag names used by the domain
68 :rtype: ``list`` of ``str``
69 """
70 return [n.strip() for n in domain.showTagNames().split(",") ]
71
72 def insertTagNames(domain,**kwargs):
73 """
74 Inserts tag names into the domain.
75
76 :param domain: a domain object
77 :type domain: ``escript.Domain``
78 :keyword <tag_name>: tag key assigned to <tag_name>
79 :type <tag_name>: ``int``
80 """
81 for k in kwargs:
82 domain.setTagMap(k,kwargs[k])
83
84 def insertTaggedValues(target,**kwargs):
85 """
86 Inserts tagged values into the target using tag names.
87
88 :param target: data to be filled by tagged values
89 :type target: `escript.Data`
90 :keyword <tag_name>: value to be used for <tag_name>
91 :type <tag_name>: ``float`` or ``numpy.ndarray``
92 :return: ``target``
93 :rtype: `escript.Data`
94 """
95 for k in kwargs:
96 target.setTaggedValue(k,kwargs[k])
97 return target
98
99
100 def interpolateTable(tab, dat, start, step, undef=1.e50, check_boundaries=False):
101 try:
102 dim=len(start)
103 except TypeError:
104 start=(start,)
105 dim=1
106 try:
107 slen=len(step)
108 except TypeError:
109 step=(step,)
110 slen=1
111 if dim<1 or dim>3:
112 raise ValueError("Length of start list must be between 1 and 3.")
113 if dim!=slen:
114 raise ValueError("Length of start and step must be the same.")
115 dshape=dat.getShape()
116 if len(dshape)==0:
117 datdim=0
118 firstdim=dat
119 else:
120 datdim=dshape[0]
121 firstdim=dat[0]
122 #So now we know firstdim is a scalar
123 if (dim==1 and datdim>1) or (dim>1 and datdim!=dim):
124 print dim, datdim
125 raise ValueError("The dimension of dat must be equal to the length of start.")
126 if dim==3:
127 d1=dat[1]
128 d2=dat[2]
129 return firstdim._interpolateTable3d(tab, start[0], step[0], d1, start[1], step[1], d2, start[2], step[2], undef, check_boundaries)
130 if dim==2:
131 d1=dat[1]
132 return d1.interpolateTable(tab, start[0], step[0], firstdim, start[1], step[1], undef, check_boundaries)
133 else:
134 return firstdim.interpolateTable(tab, start[0], step[0], undef, check_boundaries)
135
136
137 def saveDataCSV(filename, append=False, sep=", ", csep="_", **data):
138 """
139 Writes `Data` objects to a csv file.
140 These objects must have compatible FunctionSpaces. ie it must be possible to
141 interpolate all data to one `FunctionSpace`.
142
143 :param filename: file to save data to.
144 :type filename: ``string``
145 :param append: If ``True``, then open file at end rather than beginning
146 :type append: ``bool``
147 :param sep: separator between fields
148 :type sep: ``string``
149 :param csep: separator for components of rank2 and above eg ('_' -> c0_1)
150
151 The keyword args are Data objects to save.
152 If a scalar `Data` object is passed with the name ``mask``, then only
153 samples which correspond to positive values in ``mask`` will be output.
154 Example::
155
156 s=Scalar(..)
157 v=Vector(..)
158 t=Tensor(..)
159 f=float()
160 saveDataCSV("f.csv",a=s, b=v, c=t, d=f)
161
162 Will result in a file
163
164 a, b0, b1, c0_0, c0_1, .., c1_0 d
165 1.0, 1.5, 2.7, 3.1, 3.4, .., 0.89 0.0
166 0.9, 8.7, 1.9, 3.4, 7.8, .., 1.21 0.0
167
168 The first line is a header, the remaining lines give the values.
169
170 """
171 # find a function space:
172 fs = None
173 for n,d in data.items():
174 if isinstance(d, Data): fs=d.getFunctionSpace()
175 if fs == None:
176 raise ValueError, "saveDataCSV: there must be at least one Data object in the argument list."
177
178 new_data={}
179 for n,d in data.items():
180 if isinstance(d, Data):
181 new_data[n]=d
182 else:
183 try:
184 new_data[n]=Data(d,fs)
185 except:
186 raise ValueError, "saveDataCSV: unknown non-data argument type for %s"%(str(n))
187 _saveDataCSV(filename, new_data,sep, csep, append)
188
189 def saveVTK(filename,domain=None, metadata=None, metadata_schema=None, **data):
190 """
191 Deprecated. See esys.weipa.saveVTK().
192 """
193
194 import warnings
195 msg = "esys.escript.util.saveVTK is deprecated. Use the esys.weipa module."
196 warnings.warn(msg, DeprecationWarning, stacklevel=2)
197
198 from esys import weipa
199 weipa.saveVTK(filename, domain, metadata, metadata_schema, **data)
200
201 def saveDX(filename,domain=None,**data):
202 """
203 Writes `Data` objects into a file using the OpenDX file format.
204
205 Example::
206
207 tmp=Scalar(..)
208 v=Vector(..)
209 saveDX("solution.dx", temperature=tmp, velocity=v)
210
211 ``tmp`` and ``v`` are written into "solution.dx" where ``tmp`` is named
212 "temperature" and ``v`` is named "velocity".
213
214 :param filename: file name of the output file
215 :type filename: ``str``
216 :param domain: domain of the `Data` objects. If not specified, the domain
217 of the given `Data` objects is used.
218 :type domain: `escript.Domain`
219 :keyword <name>: writes the assigned value to the DX file using <name> as
220 identifier. The identifier can be used to select the data
221 set when data are imported into DX.
222 :type <name>: `Data` object
223 :note: The data objects have to be defined on the same domain. They may not
224 be in the same `FunctionSpace` but one cannot expect that all
225 `FunctionSpace` s can be mixed. Typically, data on the boundary and
226 data on the interior cannot be mixed.
227 """
228 new_data={}
229 for n,d in data.items():
230 if not d.isEmpty():
231 fs=d.getFunctionSpace()
232 domain2=fs.getDomain()
233 if fs == escript.Solution(domain2):
234 new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
235 elif fs == escript.ReducedSolution(domain2):
236 new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
237 elif fs == escript.ContinuousFunction(domain2):
238 new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
239 else:
240 new_data[n]=d
241 if domain==None: domain=domain2
242 if domain==None:
243 raise ValueError,"saveDX: no domain detected."
244 domain.saveDX(filename,new_data)
245
246 def saveESD(datasetName, dataDir=".", domain=None, timeStep=0, deltaT=1, dynamicMesh=0, **data):
247 """
248 Saves `Data` objects to files and creates an I{escript dataset} (ESD) file
249 for convenient processing/visualisation.
250
251 Single timestep example::
252
253 tmp = Scalar(..)
254 v = Vector(..)
255 saveESD("solution", "data", temperature=tmp, velocity=v)
256
257 Time series example::
258
259 while t < t_end:
260 tmp = Scalar(..)
261 v = Vector(..)
262 # save every 10 timesteps
263 if t % 10 == 0:
264 saveESD("solution", "data", timeStep=t, deltaT=10, temperature=tmp, velocity=v)
265 t = t + 1
266
267 tmp, v and the domain are saved in native format in the "data"
268 directory and the file "solution.esd" is created that refers to tmp by
269 the name "temperature" and to v by the name "velocity".
270
271 :param datasetName: name of the dataset, used to name the ESD file
272 :type datasetName: ``str``
273 :param dataDir: optional directory where the data files should be saved
274 :type dataDir: ``str``
275 :param domain: domain of the `Data` object(s). If not specified, the
276 domain of the given `Data` objects is used.
277 :type domain: `escript.Domain`
278 :param timeStep: current timestep or sequence number - first one must be 0
279 :type timeStep: `int`
280 :param deltaT: timestep or sequence increment, see example above
281 :type deltaT: `int`
282 :param dynamicMesh: by default the mesh is assumed to be static and thus
283 only saved once at timestep 0 to save disk space.
284 Setting this to 1 changes the behaviour and the mesh
285 is saved at each timestep.
286 :type dynamicMesh: `int`
287 :keyword <name>: writes the assigned value to the file using <name> as
288 identifier
289 :type <name>: `Data` object
290 :note: The ESD concept is experimental and the file format likely to
291 change so use this function with caution.
292 :note: The data objects have to be defined on the same domain (but not
293 necessarily on the same `FunctionSpace`).
294 :note: When saving a time series the first timestep must be 0 and it is
295 assumed that data from all timesteps share the domain. The dataset
296 file is updated in each iteration.
297 """
298 new_data = {}
299 for n,d in data.items():
300 if not d.isEmpty():
301 fs = d.getFunctionSpace()
302 domain2 = fs.getDomain()
303 if fs == escript.Solution(domain2):
304 new_data[n]=interpolate(d,escript.ContinuousFunction(domain2))
305 elif fs == escript.ReducedSolution(domain2):
306 new_data[n]=interpolate(d,escript.ReducedContinuousFunction(domain2))
307 else:
308 new_data[n]=d
309 if domain==None: domain=domain2
310 if domain==None:
311 raise ValueError, "saveESD: no domain detected."
312
313 if domain.onMasterProcessor() and not os.path.isdir(dataDir):
314 os.mkdir(dataDir)
315
316 meshFile = datasetName+"_mesh"
317 fileNumber = timeStep / deltaT
318
319 if dynamicMesh == 0:
320 # later timesteps reuse mesh from t=0
321 if timeStep == 0:
322 domain.dump(os.path.join(dataDir, meshFile + ".nc"))
323 else:
324 meshFile += ".%04d"
325 domain.dump(os.path.join(dataDir, (meshFile + ".nc") % fileNumber))
326
327 outputString = ""
328
329 if domain.onMasterProcessor():
330 outputString += "#escript datafile V1.0\n"
331 # number of timesteps (currently only 1 is supported)
332 outputString += "T=%d\n" % (fileNumber+1)
333 # timestep increment
334 outputString += "DT=1\n"
335 # name of the mesh file
336 outputString += "M=%s\n" % meshFile
337 # number of blocks (MPI size)
338 outputString += "N=%d\n" % domain.getMPISize()
339
340 # now add the variables
341 for varName, d in new_data.items():
342 varFile = datasetName+"_"+varName+".%04d"
343 d.dump(os.path.join(dataDir, (varFile + ".nc") % fileNumber))
344 if domain.onMasterProcessor():
345 outputString += "V=%s:%s\n" % (varFile, varName)
346
347 if domain.onMasterProcessor():
348 esdfile = open(os.path.join(dataDir, datasetName+".esd"), "w")
349 esdfile.write(outputString)
350 esdfile.close()
351
352 def kronecker(d=3):
353 """
354 Returns the kronecker delta-symbol.
355
356 :param d: dimension or an object that has the ``getDim`` method defining the
357 dimension
358 :type d: ``int``, `escript.Domain` or `escript.FunctionSpace`
359 :return: the object u of rank 2 with *u[i,j]=1* for *i=j* and *u[i,j]=0*
360 otherwise
361 :rtype: ``numpy.ndarray`` or `escript.Data` of rank 2
362 """
363 return identityTensor(d)
364
365 def identity(shape=()):
366 """
367 Returns the ``shape`` x ``shape`` identity tensor.
368
369 :param shape: input shape for the identity tensor
370 :type shape: ``tuple`` of ``int``
371 :return: array whose shape is shape x shape where *u[i,k]=1* for *i=k* and
372 *u[i,k]=0* otherwise for len(shape)=1. If len(shape)=2:
373 *u[i,j,k,l]=1* for *i=k and j=l* and *u[i,j,k,l]=0* otherwise.
374 :rtype: ``numpy.ndarray`` of rank 1, rank 2 or rank 4
375 :raise ValueError: if len(shape)>2
376 """
377 if len(shape)>0:
378 out=numpy.zeros(shape+shape,numpy.float64)
379 if len(shape)==1:
380 for i0 in range(shape[0]):
381 out[i0,i0]=1.
382 elif len(shape)==2:
383 for i0 in range(shape[0]):
384 for i1 in range(shape[1]):
385 out[i0,i1,i0,i1]=1.
386 else:
387 raise ValueError,"identity: length of shape is restricted to 2."
388 else:
389 out=1.
390 return out
391
392 def zeros(shape=()):
393 """
394 Returns the ``shape`` zero tensor.
395
396 :param shape: input shape for the identity tensor
397 :type shape: ``tuple`` of ``int``
398 :return: array of shape filled with zeros
399 :rtype: ``numpy.ndarray``
400 """
401 if len(shape)>0:
402 out=numpy.zeros(shape,numpy.float64)
403 else:
404 out=0.
405 return out
406
407 def identityTensor(d=3):
408 """
409 Returns the ``d`` x ``d`` identity matrix.
410
411 :param d: dimension or an object that has the ``getDim`` method defining the
412 dimension
413 :type d: ``int``, `escript.Domain` or `escript.FunctionSpace`
414 :return: the object u of rank 2 with *u[i,j]=1* for *i=j* and *u[i,j]=0*
415 otherwise
416 :rtype: ``numpy.ndarray`` or `escript.Data` of rank 2
417 """
418 if isinstance(d,escript.FunctionSpace):
419 return escript.Data(identity((d.getDim(),)),d)
420 elif isinstance(d,escript.Domain):
421 return identity((d.getDim(),))
422 else:
423 return identity((d,))
424
425 def identityTensor4(d=3):
426 """
427 Returns the ``d`` x ``d`` x ``d`` x ``d`` identity tensor.
428
429 :param d: dimension or an object that has the ``getDim`` method defining the
430 dimension
431 :type d: ``int`` or any object with a ``getDim`` method
432 :return: the object u of rank 4 with *u[i,j,k,l]=1* for *i=k and j=l* and
433 *u[i,j,k,l]=0* otherwise
434 :rtype: ``numpy.ndarray`` or `escript.Data` of rank 4
435 """
436 if isinstance(d,escript.FunctionSpace):
437 return escript.Data(identity((d.getDim(),d.getDim())),d)
438 elif isinstance(d,escript.Domain):
439 return identity((d.getDim(),d.getDim()))
440 else:
441 return identity((d,d))
442
443 def unitVector(i=0,d=3):
444 """
445 Returns a unit vector u of dimension d whose non-zero element is at index i.
446
447 :param i: index for non-zero element
448 :type i: ``int``
449 :param d: dimension or an object that has the ``getDim`` method defining the
450 dimension
451 :type d: ``int``, `escript.Domain` or `escript.FunctionSpace`
452 :return: the object u of rank 1 with *u[j]=1* for *j=index* and *u[j]=0*
453 otherwise
454 :rtype: ``numpy.ndarray`` or `escript.Data` of rank 1
455 """
456 return kronecker(d)[i]
457
458 #=========================================================================
459 # global reduction operations
460 #=========================================================================
461 def Lsup(arg):
462 """
463 Returns the Lsup-norm of argument ``arg``. This is the maximum absolute value
464 over all data points. This function is equivalent to ``sup(abs(arg))``.
465
466 :param arg: argument
467 :type arg: ``float``, ``int``, `escript.Data`, ``numpy.ndarray``
468 :return: maximum value of the absolute value of ``arg`` over all components
469 and all data points
470 :rtype: ``float``
471 :raise TypeError: if type of ``arg`` cannot be processed
472 """
473 if isinstance(arg,numpy.ndarray):
474 return sup(abs(arg))
475 elif isinstance(arg,escript.Data):
476 return arg._Lsup()
477 elif isinstance(arg,float):
478 return abs(arg)
479 elif isinstance(arg,int):
480 return abs(float(arg))
481 else:
482 raise TypeError,"Lsup: Unknown argument type."
483
484 def sup(arg):
485 """
486 Returns the maximum value over all data points.
487
488 :param arg: argument
489 :type arg: ``float``, ``int``, `escript.Data`, ``numpy.ndarray``
490 :return: maximum value of ``arg`` over all components and all data points
491 :rtype: ``float``
492 :raise TypeError: if type of ``arg`` cannot be processed
493 """
494 if isinstance(arg,numpy.ndarray):
495 return arg.max()
496 elif isinstance(arg,escript.Data):
497 return arg._sup()
498 elif isinstance(arg,float):
499 return arg
500 elif isinstance(arg,int):
501 return float(arg)
502 else:
503 raise TypeError,"sup: Unknown argument type."
504
505 def inf(arg):
506 """
507 Returns the minimum value over all data points.
508
509 :param arg: argument
510 :type arg: ``float``, ``int``, `escript.Data`, ``numpy.ndarray``
511 :return: minimum value of ``arg`` over all components and all data points
512 :rtype: ``float``
513 :raise TypeError: if type of ``arg`` cannot be processed
514 """
515 if isinstance(arg,numpy.ndarray):
516 return arg.min()
517 elif isinstance(arg,escript.Data):
518 return arg._inf()
519 elif isinstance(arg,float):
520 return arg
521 elif isinstance(arg,int):
522 return float(arg)
523 else:
524 raise TypeError,"inf: Unknown argument type."
525
526
527 #=========================================================================
528 # some little helpers
529 #=========================================================================
530 def getRank(arg):
531 """
532 Identifies the rank of the argument.
533
534 :param arg: an object whose rank is to be returned
535 :type arg: ``numpy.ndarray``, `escript.Data`, ``float``, ``int``
536 :return: the rank of the argument
537 :rtype: ``int``
538 :raise TypeError: if type of ``arg`` cannot be processed
539 """
540
541 if isinstance(arg,numpy.ndarray):
542 return arg.ndim
543 elif isinstance(arg,escript.Data):
544 return arg.getRank()
545 elif isinstance(arg,float):
546 return 0
547 elif isinstance(arg,int):
548 return 0
549 else:
550 raise TypeError,"getRank: Unknown argument type."
551
552 def getShape(arg):
553 """
554 Identifies the shape of the argument.
555
556 :param arg: an object whose shape is to be returned
557 :type arg: ``numpy.ndarray``, `escript.Data`, ``float``, ``int``
558 :return: the shape of the argument
559 :rtype: ``tuple`` of ``int``
560 :raise TypeError: if type of ``arg`` cannot be processed
561 """
562
563 if isinstance(arg,numpy.ndarray):
564 return arg.shape
565 elif isinstance(arg,list):
566 return numpy.array(arg).shape
567 elif isinstance(arg,escript.Data):
568 return arg.getShape()
569 elif isinstance(arg,float):
570 return ()
571 elif isinstance(arg,int):
572 return ()
573 else:
574 raise TypeError,"getShape: Cannot identify shape"
575
576 def pokeDim(arg):
577 """
578 Identifies the spatial dimension of the argument.
579
580 :param arg: an object whose spatial dimension is to be returned
581 :type arg: any
582 :return: the spatial dimension of the argument, if available, or ``None``
583 :rtype: ``int`` or ``None``
584 """
585
586 if isinstance(arg,escript.Data):
587 return arg.getFunctionSpace().getDim()
588 else:
589 return None
590
591 def commonShape(arg0, arg1):
592 """
593 Returns a shape to which ``arg0`` can be extended from the right and ``arg1``
594 can be extended from the left.
595
596 :param arg0: an object with a shape (see `getShape`)
597 :param arg1: an object with a shape (see `getShape`)
598 :return: the shape of ``arg0`` or ``arg1`` such that the left part equals the
599 shape of ``arg0`` and the right end equals the shape of ``arg1``
600 :rtype: ``tuple`` of ``int``
601 :raise ValueError: if no shape can be found
602 """
603 sh0=getShape(arg0)
604 sh1=getShape(arg1)
605 if len(sh0)<len(sh1):
606 if not sh0==sh1[:len(sh0)]:
607 raise ValueError,"argument 0 cannot be extended to the shape of argument 1"
608 return sh1
609 elif len(sh0)>len(sh1):
610 if not sh1==sh0[:len(sh1)]:
611 raise ValueError,"argument 1 cannot be extended to the shape of argument 0"
612 return sh0
613 else:
614 if not sh0==sh1:
615 raise ValueError,"argument 1 and argument 0 have not the same shape."
616 return sh0
617
618 def commonDim(*args):
619 """
620 Identifies, if possible, the spatial dimension across a set of objects
621 which may or may not have a spatial dimension.
622
623 :param args: given objects
624 :return: the spatial dimension of the objects with identifiable dimension
625 (see `pokeDim`). If none of the objects has a spatial dimension
626 ``None`` is returned.
627 :rtype: ``int`` or ``None``
628 :raise ValueError: if the objects with identifiable dimension don't have
629 the same spatial dimension.
630 """
631 out=None
632 for a in args:
633 d=pokeDim(a)
634 if not out==None:
635 if not (d==None or out==d):
636 raise ValueError,"dimension of arguments don't match"
637 else:
638 out=d
639 return out
640
641 def testForZero(arg):
642 """
643 Tests if the argument is identical to zero.
644
645 :param arg: the object to test for zero
646 :type arg: typically ``numpy.ndarray``, `escript.Data`, ``float``, ``int``
647 :return: True if the argument is identical to zero, False otherwise
648 :rtype: ``bool``
649 """
650 if isinstance(arg,numpy.ndarray):
651 return not Lsup(arg)>0.
652 elif isinstance(arg,escript.Data):
653 return False
654 elif isinstance(arg,float):
655 return not Lsup(arg)>0.
656 elif isinstance(arg,int):
657 return not Lsup(arg)>0.
658 else:
659 return False
660
661 def matchType(arg0=0.,arg1=0.):
662 """
663 Converts ``arg0`` and ``arg1`` both to the same type ``numpy.ndarray`` or
664 `escript.Data`
665
666 :param arg0: first argument
667 :type arg0: ``numpy.ndarray``,`escript.Data`,``float``, ``int``
668 :param arg1: second argument
669 :type arg1: ``numpy.ndarray``,`escript.Data`,``float``, ``int``
670 :return: a tuple representing ``arg0`` and ``arg1`` with the same type
671 :rtype: ``tuple`` of two ``numpy.ndarray`` or two `escript.Data`
672 :raise TypeError: if type of ``arg0`` or ``arg1`` cannot be processed
673 """
674 if isinstance(arg0,numpy.ndarray):
675 if isinstance(arg1,numpy.ndarray):
676 pass
677 elif isinstance(arg1,escript.Data):
678 arg0=escript.Data(arg0,arg1.getFunctionSpace())
679 elif isinstance(arg1,float):
680 arg1=numpy.array(arg1,dtype=numpy.float64)
681 elif isinstance(arg1,int):
682 arg1=numpy.array(float(arg1),dtype=numpy.float64)
683 else:
684 raise TypeError,"function: Unknown type of second argument."
685 elif isinstance(arg0,escript.Data):
686 if isinstance(arg1,numpy.ndarray):
687 arg1=escript.Data(arg1,arg0.getFunctionSpace())
688 elif isinstance(arg1,escript.Data):
689 pass
690 elif isinstance(arg1,float):
691 arg1=escript.Data(arg1,(),arg0.getFunctionSpace())
692 elif isinstance(arg1,int):
693 arg1=escript.Data(float(arg1),(),arg0.getFunctionSpace())
694 else:
695 raise TypeError,"function: Unknown type of second argument."
696 elif isinstance(arg0,float):
697 if isinstance(arg1,numpy.ndarray):
698 arg0=numpy.array(arg0,dtype=numpy.float64)
699 elif isinstance(arg1,escript.Data):
700 arg0=escript.Data(arg0,arg1.getFunctionSpace())
701 elif isinstance(arg1,float):
702 arg0=numpy.array(arg0,dtype=numpy.float64)
703 arg1=numpy.array(arg1,dtype=numpy.float64)
704 elif isinstance(arg1,int):
705 arg0=numpy.array(arg0,dtype=numpy.float64)
706 arg1=numpy.array(float(arg1),dtype=numpy.float64)
707 else:
708 raise TypeError,"function: Unknown type of second argument."
709 elif isinstance(arg0,int):
710 if isinstance(arg1,numpy.ndarray):
711 arg0=numpy.array(float(arg0),dtype=numpy.float64)
712 elif isinstance(arg1,escript.Data):
713 arg0=escript.Data(float(arg0),arg1.getFunctionSpace())
714 elif isinstance(arg1,float):
715 arg0=numpy.array(float(arg0),dtype=numpy.float64)
716 arg1=numpy.array(arg1,dtype=numpy.float64)
717 elif isinstance(arg1,int):
718 arg0=numpy.array(float(arg0),dtype=numpy.float64)
719 arg1=numpy.array(float(arg1),dtype=numpy.float64)
720 else:
721 raise TypeError,"function: Unknown type of second argument."
722 else:
723 raise TypeError,"function: Unknown type of first argument."
724
725 return arg0,arg1
726
727 def matchShape(arg0,arg1):
728 """
729 Returns a representation of ``arg0`` and ``arg1`` which have the same shape.
730
731 :param arg0: first argument
732 :type arg0: ``numpy.ndarray``,`escript.Data`,``float``, ``int``
733 :param arg1: second argument
734 :type arg1: ``numpy.ndarray``,`escript.Data`,``float``, ``int``
735 :return: ``arg0`` and ``arg1`` where copies are returned when the shape has
736 to be changed
737 :rtype: ``tuple``
738 """
739 sh=commonShape(arg0,arg1)
740 sh0=getShape(arg0)
741 sh1=getShape(arg1)
742 if len(sh0)<len(sh):
743 return outer(arg0,numpy.ones(sh[len(sh0):],numpy.float64)),arg1
744 elif len(sh1)<len(sh):
745 return arg0,outer(arg1,numpy.ones(sh[len(sh1):],numpy.float64))
746 else:
747 return arg0,arg1
748
749 def log10(arg):
750 """
751 Returns base-10 logarithm of argument ``arg``.
752
753 :param arg: argument
754 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
755 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
756 on the type of ``arg``
757 :raise TypeError: if the type of the argument is not expected
758 """
759 if isinstance(arg,numpy.ndarray):
760 return numpy.log10(arg)
761 elif isinstance(arg,escript.Data):
762 return arg._log10()
763 elif isinstance(arg,float):
764 return math.log10(arg)
765 elif isinstance(arg,int):
766 return math.log10(float(arg))
767 elif isinstance(arg,Basic):
768 return symfn.log10(arg)
769 else:
770 raise TypeError,"log10: Unknown argument type."
771
772 def wherePositive(arg):
773 """
774 Returns mask of positive values of argument ``arg``.
775
776 :param arg: argument
777 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``.
778 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
779 on the type of ``arg``
780 :raise TypeError: if the type of the argument is not expected
781 """
782 if isinstance(arg,numpy.ndarray):
783 out=numpy.greater(arg,numpy.zeros(arg.shape,numpy.float64))*1.
784 if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
785 return out
786 elif isinstance(arg,escript.Data):
787 return arg._wherePositive()
788 elif isinstance(arg,float):
789 if arg>0:
790 return 1.
791 else:
792 return 0.
793 elif isinstance(arg,int):
794 if arg>0:
795 return 1.
796 else:
797 return 0.
798 elif isinstance(arg,Basic):
799 return symfn.wherePositive(arg)
800 else:
801 raise TypeError,"wherePositive: Unknown argument type."
802
803 def whereNegative(arg):
804 """
805 Returns mask of negative values of argument ``arg``.
806
807 :param arg: argument
808 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
809 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
810 on the type of ``arg``
811 :raise TypeError: if the type of the argument is not expected
812 """
813 if isinstance(arg,numpy.ndarray):
814 out=numpy.less(arg,numpy.zeros(arg.shape,numpy.float64))*1.
815 if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
816 return out
817 elif isinstance(arg,escript.Data):
818 return arg._whereNegative()
819 elif isinstance(arg,float):
820 if arg<0:
821 return 1.
822 else:
823 return 0.
824 elif isinstance(arg,int):
825 if arg<0:
826 return 1.
827 else:
828 return 0.
829 elif isinstance(arg,Basic):
830 return symfn.whereNegative(arg)
831 else:
832 raise TypeError,"whereNegative: Unknown argument type."
833
834 def whereNonNegative(arg):
835 """
836 Returns mask of non-negative values of argument ``arg``.
837
838 :param arg: argument
839 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
840 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
841 on the type of ``arg``
842 :raise TypeError: if the type of the argument is not expected
843 """
844 if isinstance(arg,numpy.ndarray):
845 out=numpy.greater_equal(arg,numpy.zeros(arg.shape,numpy.float64))*1.
846 if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
847 return out
848 elif isinstance(arg,escript.Data):
849 return arg._whereNonNegative()
850 elif isinstance(arg,float):
851 if arg<0:
852 return 0.
853 else:
854 return 1.
855 elif isinstance(arg,int):
856 if arg<0:
857 return 0.
858 else:
859 return 1.
860 elif isinstance(arg,Basic):
861 return symfn.whereNonNegative(arg)
862 else:
863 raise TypeError,"whereNonNegative: Unknown argument type."
864
865 def whereNonPositive(arg):
866 """
867 Returns mask of non-positive values of argument ``arg``.
868
869 :param arg: argument
870 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
871 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
872 on the type of ``arg``
873 :raise TypeError: if the type of the argument is not expected
874 """
875 if isinstance(arg,numpy.ndarray):
876 out=numpy.less_equal(arg,numpy.zeros(arg.shape,numpy.float64))*1.
877 if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
878 return out
879 elif isinstance(arg,escript.Data):
880 return arg._whereNonPositive()
881 elif isinstance(arg,float):
882 if arg>0:
883 return 0.
884 else:
885 return 1.
886 elif isinstance(arg,int):
887 if arg>0:
888 return 0.
889 else:
890 return 1.
891 elif isinstance(arg,Basic):
892 return symfn.whereNonPositive(arg)
893 else:
894 raise TypeError,"whereNonPositive: Unknown argument type."
895
896 def whereZero(arg,tol=None,rtol=math.sqrt(EPSILON)):
897 """
898 Returns mask of zero entries of argument ``arg``.
899
900 :param arg: argument
901 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
902 :param tol: absolute tolerance. Values with absolute value less than tol are accepted
903 as zero. If ``tol`` is not present ``rtol``*```Lsup` (arg)`` is used.
904 :type tol: ``float``
905 :param rtol: relative tolerance used to define the absolute tolerance if ``tol`` is not present.
906 :type rtol: non-negative ``float``
907 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
908 on the type of ``arg``
909 :raise ValueError: if ``rtol`` is non-negative.
910 :raise TypeError: if the type of the argument is not expected
911 """
912 if tol == None and not isinstance(arg, Basic):
913 if rtol<0: raise ValueError,"rtol must be non-negative."
914 tol = Lsup(arg)*rtol
915 if isinstance(arg,numpy.ndarray):
916 out=numpy.less_equal(abs(arg)-tol,numpy.zeros(arg.shape,numpy.float64))*1.
917 if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
918 return out
919 elif isinstance(arg,escript.Data):
920 return arg._whereZero(tol)
921 elif isinstance(arg,float):
922 if abs(arg)<=tol:
923 return 1.
924 else:
925 return 0.
926 elif isinstance(arg,int):
927 if abs(float(arg))<=tol:
928 return 1.
929 else:
930 return 0.
931 elif isinstance(arg,Basic):
932 return symfn.whereZero(arg)
933 else:
934 raise TypeError,"whereZero: Unknown argument type."
935
936 def whereNonZero(arg,tol=0.):
937 """
938 Returns mask of values different from zero of argument ``arg``.
939
940 :param arg: argument
941 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
942 :param tol: absolute tolerance. Values with absolute value less than tol are accepted
943 as zero. If ``tol`` is not present ``rtol``*```Lsup` (arg)`` is used.
944 :type tol: ``float``
945 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
946 on the type of ``arg``
947 :raise ValueError: if ``rtol`` is non-negative.
948 :raise TypeError: if the type of the argument is not expected
949 """
950 if tol == None:
951 if rtol<=0: raise ValueError,"rtol must be non-negative."
952 tol = Lsup(arg)*rtol
953 if isinstance(arg,numpy.ndarray):
954 out=numpy.greater(abs(arg)-tol,numpy.zeros(arg.shape,numpy.float64))*1.
955 if isinstance(out,float): out=numpy.array(out,dtype=numpy.float64)
956 return out
957 elif isinstance(arg,escript.Data):
958 return arg._whereNonZero(tol)
959 elif isinstance(arg,float):
960 if abs(arg)>tol:
961 return 1.
962 else:
963 return 0.
964 elif isinstance(arg,int):
965 if abs(float(arg))>tol:
966 return 1.
967 else:
968 return 0.
969 elif isinstance(arg,Basic):
970 return symfn.whereNonZero(arg)
971 else:
972 raise TypeError,"whereNonZero: Unknown argument type."
973
974 def erf(arg):
975 """
976 Returns the error function *erf* of argument ``arg``.
977
978 :param arg: argument
979 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``.
980 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
981 on the type of ``arg``
982 :raise TypeError: if the type of the argument is not expected
983 """
984 if isinstance(arg,escript.Data):
985 return arg._erf()
986 elif isinstance(arg,Basic):
987 return symfn.erf(arg)
988 else:
989 raise TypeError,"erf: Unknown argument type."
990
991 def sin(arg):
992 """
993 Returns sine of argument ``arg``.
994
995 :param arg: argument
996 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``.
997 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
998 on the type of ``arg``
999 :raise TypeError: if the type of the argument is not expected
1000 """
1001 if isinstance(arg,numpy.ndarray):
1002 return numpy.sin(arg)
1003 elif isinstance(arg,escript.Data):
1004 return arg._sin()
1005 elif isinstance(arg,float):
1006 return math.sin(arg)
1007 elif isinstance(arg,int):
1008 return math.sin(arg)
1009 elif isinstance(arg,Basic):
1010 return symfn.sin(arg)
1011 else:
1012 raise TypeError,"sin: Unknown argument type."
1013
1014 def cos(arg):
1015 """
1016 Returns cosine of argument ``arg``.
1017
1018 :param arg: argument
1019 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
1020 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
1021 on the type of ``arg``
1022 :raise TypeError: if the type of the argument is not expected
1023 """
1024 if isinstance(arg,numpy.ndarray):
1025 return numpy.cos(arg)
1026 elif isinstance(arg,escript.Data):
1027 return arg._cos()
1028 elif isinstance(arg,float):
1029 return math.cos(arg)
1030 elif isinstance(arg,int):
1031 return math.cos(arg)
1032 elif isinstance(arg,Basic):
1033 return symfn.cos(arg)
1034 else:
1035 raise TypeError,"cos: Unknown argument type."
1036
1037 def tan(arg):
1038 """
1039 Returns tangent of argument ``arg``.
1040
1041 :param arg: argument
1042 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
1043 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
1044 on the type of ``arg``
1045 :raise TypeError: if the type of the argument is not expected
1046 """
1047 if isinstance(arg,numpy.ndarray):
1048 return numpy.tan(arg)
1049 elif isinstance(arg,escript.Data):
1050 return arg._tan()
1051 elif isinstance(arg,float):
1052 return math.tan(arg)
1053 elif isinstance(arg,int):
1054 return math.tan(arg)
1055 elif isinstance(arg,Basic):
1056 return symfn.tan(arg)
1057 else:
1058 raise TypeError,"tan: Unknown argument type."
1059
1060 def asin(arg):
1061 """
1062 Returns the inverse sine of argument ``arg``.
1063
1064 :param arg: argument
1065 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
1066 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
1067 on the type of ``arg``
1068 :raise TypeError: if the type of the argument is not expected
1069 """
1070 if isinstance(arg,numpy.ndarray):
1071 return numpy.arcsin(arg)
1072 elif isinstance(arg,escript.Data):
1073 return arg._asin()
1074 elif isinstance(arg,float):
1075 return math.asin(arg)
1076 elif isinstance(arg,int):
1077 return math.asin(arg)
1078 elif isinstance(arg,Basic):
1079 return symfn.asin(arg)
1080 else:
1081 raise TypeError,"asin: Unknown argument type."
1082
1083 def acos(arg):
1084 """
1085 Returns the inverse cosine of argument ``arg``.
1086
1087 :param arg: argument
1088 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
1089 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
1090 on the type of ``arg``
1091 :raise TypeError: if the type of the argument is not expected
1092 """
1093 if isinstance(arg,numpy.ndarray):
1094 return numpy.arccos(arg)
1095 elif isinstance(arg,escript.Data):
1096 return arg._acos()
1097 elif isinstance(arg,float):
1098 return math.acos(arg)
1099 elif isinstance(arg,int):
1100 return math.acos(arg)
1101 elif isinstance(arg,Basic):
1102 return symfn.acos(arg)
1103 else:
1104 raise TypeError,"acos: Unknown argument type."
1105
1106 def atan(arg):
1107 """
1108 Returns inverse tangent of argument ``arg``.
1109
1110 :param arg: argument
1111 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
1112 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
1113 on the type of ``arg``
1114 :raise TypeError: if the type of the argument is not expected
1115 """
1116 if isinstance(arg,numpy.ndarray):
1117 return numpy.arctan(arg)
1118 elif isinstance(arg,escript.Data):
1119 return arg._atan()
1120 elif isinstance(arg,float):
1121 return math.atan(arg)
1122 elif isinstance(arg,int):
1123 return math.atan(arg)
1124 elif isinstance(arg,Basic):
1125 return symfn.atan(arg)
1126 else:
1127 raise TypeError,"atan: Unknown argument type."
1128
1129 def sinh(arg):
1130 """
1131 Returns the hyperbolic sine of argument ``arg``.
1132
1133 :param arg: argument
1134 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
1135 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
1136 on the type of ``arg``
1137 :raise TypeError: if the type of the argument is not expected
1138 """
1139 if isinstance(arg,numpy.ndarray):
1140 return numpy.sinh(arg)
1141 elif isinstance(arg,escript.Data):
1142 return arg._sinh()
1143 elif isinstance(arg,float):
1144 return math.sinh(arg)
1145 elif isinstance(arg,int):
1146 return math.sinh(arg)
1147 elif isinstance(arg,Basic):
1148 return symfn.sinh(arg)
1149 else:
1150 raise TypeError,"sinh: Unknown argument type."
1151
1152 def cosh(arg):
1153 """
1154 Returns the hyperbolic cosine of argument ``arg``.
1155
1156 :param arg: argument
1157 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
1158 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
1159 on the type of ``arg``
1160 :raise TypeError: if the type of the argument is not expected
1161 """
1162 if isinstance(arg,numpy.ndarray):
1163 return numpy.cosh(arg)
1164 elif isinstance(arg,escript.Data):
1165 return arg._cosh()
1166 elif isinstance(arg,float):
1167 return math.cosh(arg)
1168 elif isinstance(arg,int):
1169 return math.cosh(arg)
1170 elif isinstance(arg,Basic):
1171 return symfn.cosh(arg)
1172 else:
1173 raise TypeError,"cosh: Unknown argument type."
1174
1175 def tanh(arg):
1176 """
1177 Returns the hyperbolic tangent of argument ``arg``.
1178
1179 :param arg: argument
1180 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
1181 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
1182 on the type of ``arg``
1183 :raise TypeError: if the type of the argument is not expected
1184 """
1185 if isinstance(arg,numpy.ndarray):
1186 return numpy.tanh(arg)
1187 elif isinstance(arg,escript.Data):
1188 return arg._tanh()
1189 elif isinstance(arg,float):
1190 return math.tanh(arg)
1191 elif isinstance(arg,int):
1192 return math.tanh(arg)
1193 elif isinstance(arg,Basic):
1194 return symfn.tanh(arg)
1195 else:
1196 raise TypeError,"tanh: Unknown argument type."
1197
1198 def asinh(arg):
1199 """
1200 Returns the inverse hyperbolic sine of argument ``arg``.
1201
1202 :param arg: argument
1203 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
1204 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
1205 on the type of ``arg``
1206 :raise TypeError: if the type of the argument is not expected
1207 """
1208 if isinstance(arg,numpy.ndarray):
1209 return numpy.arcsinh(arg)
1210 elif isinstance(arg,escript.Data):
1211 return arg._asinh()
1212 elif isinstance(arg,float):
1213 return numpy.arcsinh(arg)
1214 elif isinstance(arg,int):
1215 return numpy.arcsinh(float(arg))
1216 elif isinstance(arg,Basic):
1217 return symfn.asinh(arg)
1218 else:
1219 raise TypeError,"asinh: Unknown argument type."
1220
1221 def acosh(arg):
1222 """
1223 Returns the inverse hyperbolic cosine of argument ``arg``.
1224
1225 :param arg: argument
1226 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
1227 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
1228 on the type of ``arg``
1229 :raise TypeError: if the type of the argument is not expected
1230 """
1231 if isinstance(arg,numpy.ndarray):
1232 return numpy.arccosh(arg)
1233 elif isinstance(arg,escript.Data):
1234 return arg._acosh()
1235 elif isinstance(arg,float):
1236 return numpy.arccosh(arg)
1237 elif isinstance(arg,int):
1238 return numpy.arccosh(float(arg))
1239 elif isinstance(arg,Basic):
1240 return symfn.acosh(arg)
1241 else:
1242 raise TypeError,"acosh: Unknown argument type."
1243
1244 def atanh(arg):
1245 """
1246 Returns the inverse hyperbolic tangent of argument ``arg``.
1247
1248 :param arg: argument
1249 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
1250 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
1251 on the type of ``arg``
1252 :raise TypeError: if the type of the argument is not expected
1253 """
1254 if isinstance(arg,numpy.ndarray):
1255 return numpy.arctanh(arg)
1256 elif isinstance(arg,escript.Data):
1257 return arg._atanh()
1258 elif isinstance(arg,float):
1259 return numpy.arctanh(arg)
1260 elif isinstance(arg,int):
1261 return numpy.arctanh(float(arg))
1262 elif isinstance(arg,Basic):
1263 return symfn.atanh(arg)
1264 else:
1265 raise TypeError,"atanh: Unknown argument type."
1266
1267 def exp(arg):
1268 """
1269 Returns *e* to the power of argument ``arg``.
1270
1271 :param arg: argument
1272 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``.
1273 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
1274 on the type of arg
1275 :raise TypeError: if the type of the argument is not expected
1276 """
1277 if isinstance(arg,numpy.ndarray):
1278 return numpy.exp(arg)
1279 elif isinstance(arg,escript.Data):
1280 return arg._exp()
1281 elif isinstance(arg,float):
1282 return math.exp(arg)
1283 elif isinstance(arg,int):
1284 return math.exp(arg)
1285 elif isinstance(arg,Basic):
1286 return symfn.exp(arg)
1287 else:
1288 raise TypeError,"exp: Unknown argument type."
1289
1290 def sqrt(arg):
1291 """
1292 Returns the square root of argument ``arg``.
1293
1294 :param arg: argument
1295 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
1296 :rtype: ``float``, `escript.Data`, ``numpy.ndarray``
1297 depending on the type of ``arg``
1298 :raise TypeError: if the type of the argument is not expected
1299 """
1300 if isinstance(arg,numpy.ndarray):
1301 return numpy.sqrt(arg)
1302 elif isinstance(arg,escript.Data):
1303 return arg._sqrt()
1304 elif isinstance(arg,float):
1305 return math.sqrt(arg)
1306 elif isinstance(arg,int):
1307 return math.sqrt(arg)
1308 elif isinstance(arg,Basic):
1309 return symfn.sqrt(arg)
1310 else:
1311 raise TypeError,"sqrt: Unknown argument type."
1312
1313 def log(arg):
1314 """
1315 Returns the natural logarithm of argument ``arg``.
1316
1317 :param arg: argument
1318 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``.
1319 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
1320 on the type of ``arg``
1321 :raise TypeError: if the type of the argument is not expected
1322 """
1323 if isinstance(arg,numpy.ndarray):
1324 return numpy.log(arg)
1325 elif isinstance(arg,escript.Data):
1326 return arg._log()
1327 elif isinstance(arg,float):
1328 return math.log(arg)
1329 elif isinstance(arg,int):
1330 return math.log(arg)
1331 elif isinstance(arg,Basic):
1332 return symfn.log(arg)
1333 else:
1334 raise TypeError,"log: Unknown argument type."
1335
1336 def sign(arg):
1337 """
1338 Returns the sign of argument ``arg``.
1339
1340 :param arg: argument
1341 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
1342 :rtype: ``float``, `escript.Data`, ``numpy.ndarray`` depending
1343 on the type of ``arg``
1344 :raise TypeError: if the type of the argument is not expected
1345 """
1346 if isinstance(arg,numpy.ndarray):
1347 return wherePositive(arg)-whereNegative(arg)
1348 elif isinstance(arg,escript.Data):
1349 return arg._sign()
1350 elif isinstance(arg,float):
1351 if arg>0:
1352 return 1.
1353 elif arg<0:
1354 return -1.
1355 else:
1356 return 0.
1357 elif isinstance(arg,int):
1358 if float(arg)>0:
1359 return 1.
1360 elif float(arg)<0:
1361 return -1.
1362 else:
1363 return 0.
1364 elif isinstance(arg,Basic):
1365 return symfn.sign(arg)
1366 else:
1367 raise TypeError,"sign: Unknown argument type."
1368
1369 def minval(arg):
1370 """
1371 Returns the minimum value over all components of ``arg`` at each data point.
1372
1373 :param arg: argument
1374 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
1375 :rtype: ``float`` or `escript.Data` depending on the type of ``arg``
1376 :raise TypeError: if the type of the argument is not expected
1377 """
1378 if isinstance(arg,numpy.ndarray):
1379 if arg.ndim==0:
1380 return float(arg)
1381 else:
1382 return arg.min()
1383 elif isinstance(arg,escript.Data):
1384 return arg._minval()
1385 elif isinstance(arg,float):
1386 return arg
1387 elif isinstance(arg,int):
1388 return float(arg)
1389 elif isinstance(arg,Basic):
1390 return symfn.minval(arg)
1391 else:
1392 raise TypeError,"minval: Unknown argument type."
1393
1394 def maxval(arg):
1395 """
1396 Returns the maximum value over all components of ``arg`` at each data point.
1397
1398 :param arg: argument
1399 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
1400 :rtype: ``float`` or `escript.Data` depending on the type of ``arg``
1401 :raise TypeError: if the type of the argument is not expected
1402 """
1403 if isinstance(arg,numpy.ndarray):
1404 if arg.ndim==0:
1405 return float(arg)
1406 else:
1407 return arg.max()
1408 elif isinstance(arg,escript.Data):
1409 return arg._maxval()
1410 elif isinstance(arg,float):
1411 return arg
1412 elif isinstance(arg,int):
1413 return float(arg)
1414 elif isinstance(arg,Basic):
1415 return symfn.maxval(arg)
1416 else:
1417 raise TypeError,"maxval: Unknown argument type."
1418
1419 def length(arg):
1420 """
1421 Returns the length (Euclidean norm) of argument ``arg`` at each data point.
1422
1423 :param arg: argument
1424 :type arg: ``float``, `escript.Data`, ``numpy.ndarray``
1425 :rtype: ``float`` or `escript.Data` depending on the type of ``arg``
1426 """
1427 return sqrt(inner(arg,arg))
1428
1429 def trace(arg,axis_offset=0):
1430 """
1431 Returns the trace of ``arg`` which is the sum of ``arg[k,k]`` over k.
1432
1433 :param arg: argument
1434 :type arg: `escript.Data`, ``numpy.ndarray``
1435 :param axis_offset: ``axis_offset`` to components to sum over. ``axis_offset``
1436 must be non-negative and less than the rank of ``arg`` +1.
1437 The dimensions of component ``axis_offset`` and
1438 axis_offset+1 must be equal.
1439 :type axis_offset: ``int``
1440 :return: trace of arg. The rank of the returned object is rank of ``arg``
1441 minus 2.
1442 :rtype: `escript.Data` or ``numpy.ndarray`` depending on the type of ``arg``
1443 """
1444 if isinstance(arg,numpy.ndarray):
1445 sh=arg.shape
1446 if len(sh)<2:
1447 raise ValueError,"rank of argument must be greater than 1"
1448 if axis_offset<0 or axis_offset>len(sh)-2:
1449 raise ValueError,"axis_offset must be between 0 and %d"%(len(sh)-2)
1450 s1=1
1451 for i in range(axis_offset): s1*=sh[i]
1452 s2=1
1453 for i in range(axis_offset+2,len(sh)): s2*=sh[i]
1454 if not sh[axis_offset] == sh[axis_offset+1]:
1455 raise ValueError,"dimensions of component %d and %d must match."%(axis_offset,axis_offset+1)
1456 arg_reshaped=numpy.reshape(arg,(s1,sh[axis_offset],sh[axis_offset],s2))
1457 out=numpy.zeros([s1,s2],numpy.float64)
1458 for i1 in range(s1):
1459 for i2 in range(s2):
1460 for j in range(sh[axis_offset]): out[i1,i2]+=arg_reshaped[i1,j,j,i2]
1461 out.resize(sh[:axis_offset]+sh[axis_offset+2:])
1462 return out
1463 elif isinstance(arg,escript.Data):
1464 if arg.getRank()<2:
1465 raise ValueError,"rank of argument must be greater than 1"
1466 if axis_offset<0 or axis_offset>arg.getRank()-2:
1467 raise ValueError,"axis_offset must be between 0 and %d"%(arg.getRank()-2)
1468 s=list(arg.getShape())
1469 if not s[axis_offset] == s[axis_offset+1]:
1470 raise ValueError,"dimensions of component %d and %d must match."%(axis_offset,axis_offset+1)
1471 return arg._trace(axis_offset)
1472 elif isinstance(arg,Basic):
1473 return symfn.trace(arg,axis_offset)
1474 elif isinstance(arg,float):
1475 raise TypeError,"illegal argument type float."
1476 elif isinstance(arg,int):
1477 raise TypeError,"illegal argument type int."
1478 else:
1479 raise TypeError,"Unknown argument type."
1480
1481 def transpose(arg,axis_offset=None):
1482 """
1483 Returns the transpose of ``arg`` by swapping the first ``axis_offset`` and the
1484 last ``rank-axis_offset`` components.
1485
1486 :param arg: argument
1487 :type arg: `escript.Data`, ``numpy.ndarray``, ``float``, ``int``
1488 :param axis_offset: the first ``axis_offset`` components are swapped with the
1489 rest. ``axis_offset`` must be non-negative and less or
1490 equal to the rank of ``arg``. If ``axis_offset`` is not
1491 present ``int(r/2)`` where r is the rank of ``arg`` is
1492 used.
1493 :type axis_offset: ``int``
1494 :return: transpose of ``arg``
1495 :rtype: `escript.Data`, ``numpy.ndarray``, ``float``, ``int``
1496 depending on the type of ``arg``
1497 """
1498 if isinstance(arg,numpy.ndarray):
1499 if axis_offset==None: axis_offset=int(arg.ndim/2)
1500 return numpy.transpose(arg,axes=range(axis_offset,arg.ndim)+range(0,axis_offset))
1501 elif isinstance(arg,escript.Data):
1502 r=arg.getRank()
1503 if axis_offset==None: axis_offset=int(r/2)
1504 if axis_offset<0 or axis_offset>r:
1505 raise ValueError,"axis_offset must be between 0 and %s"%r
1506 return arg._transpose(axis_offset)
1507 elif isinstance(arg,float):
1508 if not ( axis_offset==0 or axis_offset==None):
1509 raise ValueError,"axis_offset must be 0 for float argument"
1510 return arg
1511 elif isinstance(arg,int):
1512 if not ( axis_offset==0 or axis_offset==None):
1513 raise ValueError,"axis_offset must be 0 for int argument"
1514 return float(arg)
1515 elif isinstance(arg,Basic):
1516 return symfn.transpose(arg,axis_offset)
1517 else:
1518 raise TypeError,"Unknown argument type."
1519
1520 def swap_axes(arg,axis0=0,axis1=1):
1521 """
1522 Returns the swap of ``arg`` by swapping the components ``axis0`` and ``axis1``.
1523
1524 :param arg: argument
1525 :type arg: `escript.Data`, ``numpy.ndarray``
1526 :param axis0: first axis. ``axis0`` must be non-negative and less than the
1527 rank of ``arg``.
1528 :type axis0: ``int``
1529 :param axis1: second axis. ``axis1`` must be non-negative and less than the
1530 rank of ``arg``.
1531 :type axis1: ``int``
1532 :return: ``arg`` with swapped components
1533 :rtype: `escript.Data` or ``numpy.ndarray`` depending on the type of ``arg``
1534 """
1535 if axis0 > axis1:
1536 axis0,axis1=axis1,axis0
1537 if isinstance(arg,numpy.ndarray):
1538 return numpy.swapaxes(arg,axis0,axis1)
1539 elif isinstance(arg,escript.Data):
1540 return arg._swap_axes(axis0,axis1)
1541 elif isinstance(arg,float):
1542 raise TypeError,"float argument is not supported."
1543 elif isinstance(arg,int):
1544 raise TypeError,"int argument is not supported."
1545 elif isinstance(arg,Basic):
1546 return symfn.swap_axes(arg, axis0, axis1)
1547 else:
1548 raise TypeError,"Unknown argument type."
1549
1550 def symmetric(arg):
1551 """
1552 Returns the symmetric part of the square matrix ``arg``. That is,
1553 *(arg+transpose(arg))/2*.
1554
1555 :param arg: input matrix. Must have rank 2 or 4 and be square.
1556 :type arg: ``numpy.ndarray``, `escript.Data`
1557 :return: symmetric part of ``arg``
1558 :rtype: ``numpy.ndarray`` or `escript.Data` depending on the input
1559 """
1560 if isinstance(arg,numpy.ndarray):
1561 if arg.ndim==2:
1562 if not (arg.shape[0]==arg.shape[1]):
1563 raise ValueError,"argument must be square."
1564 elif arg.ndim==4:
1565 if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
1566 raise ValueError,"argument must be square."
1567 else:
1568 raise ValueError,"rank 2 or 4 is required."
1569 return (arg+transpose(arg))/2
1570 elif isinstance(arg,escript.Data):
1571 if arg.getRank()==2:
1572 if not (arg.getShape()[0]==arg.getShape()[1]):
1573 raise ValueError,"argument must be square."
1574 return arg._symmetric()
1575 elif arg.getRank()==4:
1576 if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
1577 raise ValueError,"argument must be square."
1578 return arg._symmetric()
1579 else:
1580 raise ValueError,"rank 2 or 4 is required."
1581 elif isinstance(arg,float):
1582 return arg
1583 elif isinstance(arg,int):
1584 return float(arg)
1585 elif isinstance(arg,Basic):
1586 return symfn.symmetric(arg)
1587 else:
1588 raise TypeError,"symmetric: Unknown argument type."
1589
1590 def nonsymmetric(arg):
1591 """
1592 Returns the non-symmetric part of the square matrix ``arg``. That is,
1593 *(arg-transpose(arg))/2*.
1594
1595 :param arg: input matrix. Must have rank 2 or 4 and be square.
1596 :type arg: ``numpy.ndarray``, `escript.Data`
1597 :return: non-symmetric part of ``arg``
1598 :rtype: ``numpy.ndarray`` or `escript.Data` depending on the input
1599 """
1600 if isinstance(arg,numpy.ndarray):
1601 if arg.ndim==2:
1602 if not (arg.shape[0]==arg.shape[1]):
1603 raise ValueError,"nonsymmetric: argument must be square."
1604 elif arg.ndim==4:
1605 if not (arg.shape[0]==arg.shape[2] and arg.shape[1]==arg.shape[3]):
1606 raise ValueError,"nonsymmetric: argument must be square."
1607 else:
1608 raise ValueError,"nonsymmetric: rank 2 or 4 is required."
1609 return (arg-transpose(arg))/2
1610 elif isinstance(arg,escript.Data):
1611 if arg.getRank()==2:
1612 if not (arg.getShape()[0]==arg.getShape()[1]):
1613 raise ValueError,"argument must be square."
1614 return arg._nonsymmetric()
1615 elif arg.getRank()==4:
1616 if not (arg.getShape()[0]==arg.getShape()[2] and arg.getShape()[1]==arg.getShape()[3]):
1617 raise ValueError,"argument must be square."
1618 return arg._nonsymmetric()
1619 else:
1620 raise ValueError,"rank 2 or 4 is required."
1621 elif isinstance(arg,float):
1622 return arg
1623 elif isinstance(arg,int):
1624 return float(arg)
1625 elif isinstance(arg,Basic):
1626 return symfn.nonsymmetric(arg)
1627 else:
1628 raise TypeError,"nonsymmetric: Unknown argument type."
1629
1630 def inverse(arg):
1631 """
1632 Returns the inverse of the square matrix ``arg``.
1633
1634 :param arg: square matrix. Must have rank 2 and the first and second
1635 dimension must be equal.
1636 :type arg: ``numpy.ndarray``, `escript.Data`
1637 :return: inverse of the argument. ``matrix_mult(inverse(arg),arg)`` will be
1638 almost equal to ``kronecker(arg.getShape()[0])``
1639 :rtype: ``numpy.ndarray`` or `escript.Data` depending on the input
1640 :note: for `escript.Data` objects the dimension is restricted to 3.
1641 """
1642 import numpy.linalg
1643 if isinstance(arg,numpy.ndarray):
1644 return numpy.linalg.tensorinv(arg,ind=1)
1645 elif isinstance(arg,escript.Data):
1646 return escript_inverse(arg)
1647 elif isinstance(arg,float):
1648 return 1./arg
1649 elif isinstance(arg,int):
1650 return 1./float(arg)
1651 elif isinstance(arg,Basic):
1652 return symfn.inverse(arg)
1653 else:
1654 raise TypeError,"inverse: Unknown argument type."
1655
1656 def escript_inverse(arg): # this should be escript._inverse and use LAPACK
1657 "arg is a Data object!"
1658 return arg._inverse()
1659
1660 #if not arg.getRank()==2:
1661 #raise ValueError,"escript_inverse: argument must have rank 2"
1662 #s=arg.getShape()
1663 #if not s[0] == s[1]:
1664 #raise ValueError,"escript_inverse: argument must be a square matrix."
1665 #out=escript.Data(0.,s,arg.getFunctionSpace())
1666 #if s[0]==1:
1667 #if inf(abs(arg[0,0]))==0: # in c this should be done point wise as abs(arg[0,0](i))<=0.
1668 #raise ZeroDivisionError,"escript_inverse: argument not invertible"
1669 #out[0,0]=1./arg[0,0]
1670 #elif s[0]==2:
1671 #A11=arg[0,0]
1672 #A12=arg[0,1]
1673 #A21=arg[1,0]
1674 #A22=arg[1,1]
1675 #D = A11*A22-A12*A21
1676 #if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
1677 #raise ZeroDivisionError,"escript_inverse: argument not invertible"
1678 #D=1./D
1679 #out[0,0]= A22*D
1680 #out[1,0]=-A21*D
1681 #out[0,1]=-A12*D
1682 #out[1,1]= A11*D
1683 #elif s[0]==3:
1684 #A11=arg[0,0]
1685 #A21=arg[1,0]
1686 #A31=arg[2,0]
1687 #A12=arg[0,1]
1688 #A22=arg[1,1]
1689 #A32=arg[2,1]
1690 #A13=arg[0,2]
1691 #A23=arg[1,2]
1692 #A33=arg[2,2]
1693 #D = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22)
1694 #if inf(abs(D))==0: # in c this should be done point wise as abs(D(i))<=0.
1695 #raise ZeroDivisionError,"escript_inverse: argument not invertible"
1696 #D=1./D
1697 #out[0,0]=(A22*A33-A23*A32)*D
1698 #out[1,0]=(A31*A23-A21*A33)*D
1699 #out[2,0]=(A21*A32-A31*A22)*D
1700 #out[0,1]=(A13*A32-A12*A33)*D
1701 #out[1,1]=(A11*A33-A31*A13)*D
1702 #out[2,1]=(A12*A31-A11*A32)*D
1703 #out[0,2]=(A12*A23-A13*A22)*D
1704 #out[1,2]=(A13*A21-A11*A23)*D
1705 #out[2,2]=(A11*A22-A12*A21)*D
1706 #else:
1707 #raise TypeError,"escript_inverse: only matrix dimensions 1,2,3 are supported right now."
1708 #return out
1709
1710 def eigenvalues(arg):
1711 """
1712 Returns the eigenvalues of the square matrix ``arg``.
1713
1714 :param arg: square matrix. Must have rank 2 and the first and second
1715 dimension must be equal. It must also be symmetric, ie.
1716 ``transpose(arg)==arg`` (this is not checked).
1717 :type arg: ``numpy.ndarray``, `escript.Data`
1718 :return: the eigenvalues in increasing order
1719 :rtype: ``numpy.ndarray`` or `escript.Data` depending on the input
1720 :note: for `escript.Data` objects the dimension is restricted to 3.
1721 """
1722 if isinstance(arg,numpy.ndarray):
1723 out=numpy.linalg.eigvals((arg+numpy.transpose(arg))/2.)
1724 out.sort()
1725 return out
1726 elif isinstance(arg,escript.Data):
1727 return arg._eigenvalues()
1728 elif isinstance(arg,float):
1729 return arg
1730 elif isinstance(arg,int):
1731 return float(arg)
1732 elif isinstance(arg,Basic):
1733 return symfn.eigenvalues(arg)
1734 else:
1735 raise TypeError,"eigenvalues: Unknown argument type."
1736
1737 def eigenvalues_and_eigenvectors(arg):
1738 """
1739 Returns the eigenvalues and eigenvectors of the square matrix ``arg``.
1740
1741 :param arg: square matrix. Must have rank 2 and the first and second
1742 dimension must be equal. It must also be symmetric, ie.
1743 ``transpose(arg)==arg`` (this is not checked).
1744 :type arg: `escript.Data`
1745 :return: the eigenvalues and eigenvectors. The eigenvalues are ordered by
1746 increasing value. The eigenvectors are orthogonal and normalized.
1747 If V are the eigenvectors then V[:,i] is the eigenvector
1748 corresponding to the i-th eigenvalue.
1749 :rtype: `tuple` of `escript.Data`
1750 :note: The dimension is restricted to 3.
1751 """
1752 if isinstance(arg,numpy.ndarray):
1753 raise TypeError,"eigenvalues_and_eigenvectors does not support numpy.ndarray arguments"
1754 elif isinstance(arg,escript.Data):
1755 return arg._eigenvalues_and_eigenvectors()
1756 elif isinstance(arg,float):
1757 return (numpy.array([[arg]],numpy.float_),numpy.ones((1,1),numpy.float_))
1758 elif isinstance(arg,int):
1759 return (numpy.array([[arg]],numpy.float_),numpy.ones((1,1),numpy.float_))
1760 elif isinstance(arg,Basic):
1761 return symfn.eigenvalues_and_eigenvectors(arg)
1762 else:
1763 raise TypeError,"eigenvalues: Unknown argument type."
1764
1765 def maximum(*args):
1766 """
1767 The maximum over arguments ``args``.
1768
1769 :param args: arguments
1770 :type args: ``numpy.ndarray``, `escript.Data`, ``int`` or ``float``
1771 :return: an object which in each entry gives the maximum of the
1772 corresponding values in ``args``
1773 :rtype: ``numpy.ndarray``, `escript.Data`, ``int`` or
1774 ``float`` depending on the input
1775 """
1776 if max([isinstance(v,Basic) for v in args]):
1777 return symfn.maximum(*args)
1778 out=None
1779 for a in args:
1780 if out==None:
1781 out=a*1.
1782 else:
1783 if isinstance(out,escript.Data) and isinstance(a,escript.Data):
1784 if out.getRank()==0 and a.getRank()>0:
1785 # We need to consider the case where we have scalars and
1786 # higher ranked objects mixed. If the scalar was first it will
1787 # get picked as the initial out and we have a problem, so we
1788 # swap the objects
1789 res=a.copy() #Deep copy of a
1790 res.copyWithMask(out,wherePositive(out-a))
1791 out=res
1792 else:
1793 out.copyWithMask(a,wherePositive(a-out))
1794 else:
1795 if isinstance(a,numpy.ndarray):
1796 diff=-out+a
1797 else:
1798 diff=a-out
1799 out=wherePositive(diff)*diff+out
1800 return out
1801
1802 def minimum(*args):
1803 """
1804 The minimum over arguments ``args``.
1805
1806 :param args: arguments
1807 :type args: ``numpy.ndarray``, `escript.Data`, ``int`` or ``float``
1808 :return: an object which gives in each entry the minimum of the
1809 corresponding values in ``args``
1810 :rtype: ``numpy.ndarray``, `escript.Data`, ``int`` or ``float``
1811 depending on the input
1812 """
1813 if max([isinstance(v,Basic) for v in args]):
1814 return symfn.minimum(*args)
1815 out=None
1816 for a in args:
1817 if out==None:
1818 out=a*1.
1819 else:
1820 if isinstance(out,escript.Data) and isinstance(a,escript.Data):
1821 if out.getRank()==0 and a.getRank()>0:
1822 # We need to consider the case where we have scalars and
1823 # higher ranked objects mixed. If the scalar was first it will
1824 # get picked as the initial out and we have a problem, so we
1825 # swap the objects
1826 res=a.copy() # Deep copy of a
1827 res.copyWithMask(out,whereNegative(out-a))
1828 out=res
1829 else:
1830 out.copyWithMask(a,whereNegative(a-out))
1831 else:
1832 if isinstance(a,numpy.ndarray):
1833 diff=-out+a
1834 else:
1835 diff=a-out
1836 out=whereNegative(diff)*diff+out
1837 return out
1838
1839 def clip(arg,minval=None,maxval=None):
1840 """
1841 Cuts the values of ``arg`` between ``minval`` and ``maxval``.
1842
1843 :param arg: argument
1844 :type arg: ``numpy.ndarray``, `escript.Data`, ``int`` or ``float``
1845 :param minval: lower range. If None no lower range is applied
1846 :type minval: ``float`` or ``None``
1847 :param maxval: upper range. If None no upper range is applied
1848 :type maxval: ``float`` or ``None``
1849 :return: an object that contains all values from ``arg`` between ``minval``
1850 and ``maxval``
1851 :rtype: ``numpy.ndarray``, `escript.Data`, ``int`` or ``float``
1852 depending on the input
1853 :raise ValueError: if ``minval>maxval``
1854 """
1855 if isinstance(arg, Basic):
1856 return symfn.clip(arg,minval,maxval)
1857 if not minval==None and not maxval==None:
1858 if minval>maxval:
1859 raise ValueError,"minval = %s must be less than maxval %s"%(minval,maxval)
1860 if minval == None:
1861 tmp=arg
1862 else:
1863 tmp=maximum(minval,arg)
1864 if maxval == None:
1865 return tmp
1866 else:
1867 return minimum(tmp,maxval)
1868
1869 def inner(arg0,arg1):
1870 """
1871 Inner product of the two arguments. The inner product is defined as:
1872
1873 C{out=Sigma_s arg0[s]*arg1[s]}
1874
1875 where s runs through ``arg0.Shape``.
1876
1877 ``arg0`` and ``arg1`` must have the same shape.
1878
1879 :param arg0: first argument
1880 :type arg0: ``numpy.ndarray``, `escript.Data`, ``float``, ``int``
1881 :param arg1: second argument
1882 :type arg1: ``numpy.ndarray``, `escript.Data`, ``float``, ``int``
1883 :return: the inner product of ``arg0`` and ``arg1`` at each data point
1884 :rtype: ``numpy.ndarray``, `escript.Data`, ``float``
1885 depending on the input
1886 :raise ValueError: if the shapes of the arguments are not identical
1887 """
1888 if isinstance(arg0,Basic) and isinstance(arg1,Basic):
1889 return symfn.inner(arg0,arg1)
1890
1891 sh0=getShape(arg0)
1892 sh1=getShape(arg1)
1893 if not sh0==sh1:
1894 raise ValueError,"inner: shape of arguments does not match"
1895 return generalTensorProduct(arg0,arg1,axis_offset=len(sh0))
1896
1897 def outer(arg0,arg1):
1898 """
1899 The outer product of the two arguments. The outer product is defined as:
1900
1901 ``out[t,s]=arg0[t]*arg1[s]``
1902
1903 where
1904 - s runs through ``arg0.Shape``
1905 - t runs through ``arg1.Shape``
1906
1907 :param arg0: first argument
1908 :type arg0: ``numpy.ndarray``, `escript.Data`, ``float``, ``int``
1909 :param arg1: second argument
1910 :type arg1: ``numpy.ndarray``, `escript.Data`, ``float``, ``int``
1911 :return: the outer product of ``arg0`` and ``arg1`` at each data point
1912 :rtype: ``numpy.ndarray`` or `escript.Data` depending on the input
1913 """
1914 if isinstance(arg0,Basic) and isinstance(arg1,Basic):
1915 return symfn.outer(arg0,arg1)
1916
1917 return generalTensorProduct(arg0,arg1,axis_offset=0)
1918
1919 def matrixmult(arg0,arg1):
1920 """
1921 See `matrix_mult`.
1922 """
1923 return matrix_mult(arg0,arg1)
1924
1925 def matrix_mult(arg0,arg1):
1926 """
1927 matrix-matrix or matrix-vector product of the two arguments.
1928
1929 C{out[s0]=Sigma_{r0} arg0[s0,r0]*arg1[r0]}
1930
1931 or
1932
1933 C{out[s0,s1]=Sigma_{r0} arg0[s0,r0]*arg1[r0,s1]}
1934
1935 The second dimension of ``arg0`` and the first dimension of ``arg1`` must
1936 match.
1937
1938 :param arg0: first argument of rank 2
1939 :type arg0: ``numpy.ndarray``, `escript.Data`
1940 :param arg1: second argument of at least rank 1
1941 :type arg1: ``numpy.ndarray``, `escript.Data`
1942 :return: the matrix-matrix or matrix-vector product of ``arg0`` and ``arg1``
1943 at each data point
1944 :rtype: ``numpy.ndarray`` or `escript.Data` depending on the input
1945 :raise ValueError: if the shapes of the arguments are not appropriate
1946 """
1947 sh0=getShape(arg0)
1948 sh1=getShape(arg1)
1949 if not len(sh0)==2 :
1950 raise ValueError,"first argument must have rank 2"
1951 if not len(sh1)==2 and not len(sh1)==1:
1952 raise ValueError,"second argument must have rank 1 or 2"
1953 return generalTensorProduct(arg0,arg1,axis_offset=1)
1954
1955 def tensormult(arg0,arg1):
1956 """
1957 See `tensor_mult`.
1958 """
1959 return tensor_mult(arg0,arg1)
1960
1961 def tensor_mult(arg0,arg1):
1962 """
1963 The tensor product of the two arguments.
1964
1965 For ``arg0`` of rank 2 this is
1966
1967 C{out[s0]=Sigma_{r0} arg0[s0,r0]*arg1[r0]}
1968
1969 or
1970
1971 C{out[s0,s1]=Sigma_{r0} arg0[s0,r0]*arg1[r0,s1]}
1972
1973 and for ``arg0`` of rank 4 this is
1974
1975 C{out[s0,s1,s2,s3]=Sigma_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1,s2,s3]}
1976
1977 or
1978
1979 C{out[s0,s1,s2]=Sigma_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1,s2]}
1980
1981 or
1982
1983 C{out[s0,s1]=Sigma_{r0,r1} arg0[s0,s1,r0,r1]*arg1[r0,r1]}
1984
1985 In the first case the second dimension of ``arg0`` and the last dimension of
1986 ``arg1`` must match and in the second case the two last dimensions of ``arg0``
1987 must match the two first dimensions of ``arg1``.
1988
1989 :param arg0: first argument of rank 2 or 4
1990 :type arg0: ``numpy.ndarray``, `escript.Data`
1991 :param arg1: second argument of shape greater than 1 or 2 depending on the
1992 rank of ``arg0``
1993 :type arg1: ``numpy.ndarray``, `escript.Data`
1994 :return: the tensor product of ``arg0`` and ``arg1`` at each data point
1995 :rtype: ``numpy.ndarray`` or `escript.Data` depending on the input
1996 """
1997 sh0=getShape(arg0)
1998 sh1=getShape(arg1)
1999 if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
2000 return generalTensorProduct(arg0,arg1,axis_offset=1)
2001 elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
2002 return generalTensorProduct(arg0,arg1,axis_offset=2)
2003 else:
2004 raise ValueError,"tensor_mult: first argument must have rank 2 or 4"
2005
2006 def generalTensorProduct(arg0,arg1,axis_offset=0):
2007 """
2008 Generalized tensor product.
2009
2010 C{out[s,t]=Sigma_r arg0[s,r]*arg1[r,t]}
2011
2012 where
2013 - s runs through ``arg0.Shape[:arg0.ndim-axis_offset]``
2014 - r runs through ``arg0.Shape[:axis_offset]``
2015 - t runs through ``arg1.Shape[axis_offset:]``
2016
2017 :param arg0: first argument
2018 :type arg0: ``numpy.ndarray``, `escript.Data`, ``float``, ``int``
2019 :param arg1: second argument
2020 :type arg1: ``numpy.ndarray``, `escript.Data`, ``float``, ``int``
2021 :return: the general tensor product of ``arg0`` and ``arg1`` at each data
2022 point
2023 :rtype: ``numpy.ndarray``, `escript.Data` depending on the input
2024 """
2025 if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
2026 arg0,arg1=matchType(arg0,arg1)
2027 # at this stage arg0 and arg1 are both numpy.ndarray or escript.Data
2028 if isinstance(arg0,numpy.ndarray):
2029 if not arg0.shape[arg0.ndim-axis_offset:]==arg1.shape[:axis_offset]:
2030 raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
2031 arg0_c=arg0.copy()
2032 arg1_c=arg1.copy()
2033 sh0,sh1=arg0.shape,arg1.shape
2034 d0,d1,d01=1,1,1
2035 for i in sh0[:arg0.ndim-axis_offset]: d0*=i
2036 for i in sh1[axis_offset:]: d1*=i
2037 for i in sh1[:axis_offset]: d01*=i
2038 arg0_c.resize((d0,d01))
2039 arg1_c.resize((d01,d1))
2040 out=numpy.zeros((d0,d1),numpy.float64)
2041 for i0 in range(d0):
2042 for i1 in range(d1):
2043 out[i0,i1]=numpy.sum(arg0_c[i0,:]*arg1_c[:,i1])
2044 out.resize(sh0[:arg0.ndim-axis_offset]+sh1[axis_offset:])
2045 return out
2046 elif isinstance(arg0,escript.Data):
2047 return escript_generalTensorProduct(arg0,arg1,axis_offset) # this calls has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
2048
2049 def escript_generalTensorProduct(arg0,arg1,axis_offset,transpose=0):
2050 "arg0 and arg1 are both Data objects but not necessarily on the same function space. They could be identical!!!"
2051 return C_GeneralTensorProduct(arg0, arg1, axis_offset, transpose)
2052
2053 def transposed_matrix_mult(arg0,arg1):
2054 """
2055 transposed(matrix)-matrix or transposed(matrix)-vector product of the two
2056 arguments.
2057
2058 C{out[s0]=Sigma_{r0} arg0[r0,s0]*arg1[r0]}
2059
2060 or
2061
2062 C{out[s0,s1]=Sigma_{r0} arg0[r0,s0]*arg1[r0,s1]}
2063
2064 The function call ``transposed_matrix_mult(arg0,arg1)`` is equivalent to
2065 ``matrix_mult(transpose(arg0),arg1)``.
2066
2067 The first dimension of ``arg0`` and ``arg1`` must match.
2068
2069 :param arg0: first argument of rank 2
2070 :type arg0: ``numpy.ndarray``, `escript.Data`
2071 :param arg1: second argument of at least rank 1
2072 :type arg1: ``numpy.ndarray``, `escript.Data`
2073 :return: the product of the transpose of ``arg0`` and ``arg1`` at each data
2074 point
2075 :rtype: ``numpy.ndarray`` or `escript.Data` depending on the input
2076 :raise ValueError: if the shapes of the arguments are not appropriate
2077 """
2078 sh0=getShape(arg0)
2079 sh1=getShape(arg1)
2080 if not len(sh0)==2 :
2081 raise ValueError,"first argument must have rank 2"
2082 if not len(sh1)==2 and not len(sh1)==1:
2083 raise ValueError,"second argument must have rank 1 or 2"
2084 return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
2085
2086 def transposed_tensor_mult(arg0,arg1):
2087 """
2088 The tensor product of the transpose of the first and the second argument.
2089
2090 For ``arg0`` of rank 2 this is
2091
2092 C{out[s0]=Sigma_{r0} arg0[r0,s0]*arg1[r0]}
2093
2094 or
2095
2096 C{out[s0,s1]=Sigma_{r0} arg0[r0,s0]*arg1[r0,s1]}
2097
2098 and for ``arg0`` of rank 4 this is
2099
2100 C{out[s0,s1,s2,s3]=Sigma_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2,s3]}
2101
2102 or
2103
2104 C{out[s0,s1,s2]=Sigma_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1,s2]}
2105
2106 or
2107
2108 C{out[s0,s1]=Sigma_{r0,r1} arg0[r0,r1,s0,s1]*arg1[r0,r1]}
2109
2110 In the first case the first dimension of ``arg0`` and the first dimension of
2111 ``arg1`` must match and in the second case the two first dimensions of
2112 ``arg0`` must match the two first dimensions of ``arg1``.
2113
2114 The function call ``transposed_tensor_mult(arg0,arg1)`` is equivalent to
2115 ``tensor_mult(transpose(arg0),arg1)``.
2116
2117 :param arg0: first argument of rank 2 or 4
2118 :type arg0: ``numpy.ndarray``, `escript.Data`
2119 :param arg1: second argument of shape greater of 1 or 2 depending on the
2120 rank of ``arg0``
2121 :type arg1: ``numpy.ndarray``, `escript.Data`
2122 :return: the tensor product of transpose of arg0 and arg1 at each data point
2123 :rtype: ``numpy.ndarray``, `escript.Data` depending on the input
2124 """
2125 sh0=getShape(arg0)
2126 sh1=getShape(arg1)
2127 if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
2128 return generalTransposedTensorProduct(arg0,arg1,axis_offset=1)
2129 elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
2130 return generalTransposedTensorProduct(arg0,arg1,axis_offset=2)
2131 else:
2132 raise ValueError,"first argument must have rank 2 or 4"
2133
2134 def generalTransposedTensorProduct(arg0,arg1,axis_offset=0):
2135 """
2136 Generalized tensor product of transposed of ``arg0`` and ``arg1``.
2137
2138 C{out[s,t]=Sigma_r arg0[r,s]*arg1[r,t]}
2139
2140 where
2141 - s runs through ``arg0.Shape[axis_offset:]``
2142 - r runs through ``arg0.Shape[:axis_offset]``
2143 - t runs through ``arg1.Shape[axis_offset:]``
2144
2145 The function call ``generalTransposedTensorProduct(arg0,arg1,axis_offset)``
2146 is equivalent to
2147 ``generalTensorProduct(transpose(arg0,arg0.ndim-axis_offset),arg1,axis_offset)``.
2148
2149 :param arg0: first argument
2150 :type arg0: ``numpy.ndarray``, `escript.Data`, ``float``, ``int``
2151 :param arg1: second argument
2152 :type arg1: ``numpy.ndarray``, `escript.Data`, ``float``, ``int``
2153 :return: the general tensor product of ``transpose(arg0)`` and ``arg1`` at
2154 each data point
2155 :rtype: ``numpy.ndarray``, `escript.Data` depending on the input
2156 """
2157 if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
2158 arg0,arg1=matchType(arg0,arg1)
2159 # at this stage arg0 and arg1 are both numpy.ndarray or escript.Data
2160 if isinstance(arg0,numpy.ndarray):
2161 if not arg0.shape[:axis_offset]==arg1.shape[:axis_offset]:
2162 raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
2163 arg0_c=arg0.copy()
2164 arg1_c=arg1.copy()
2165 sh0,sh1=arg0.shape,arg1.shape
2166 d0,d1,d01=1,1,1
2167 for i in sh0[axis_offset:]: d0*=i
2168 for i in sh1[axis_offset:]: d1*=i
2169 for i in sh0[:axis_offset]: d01*=i
2170 arg0_c.resize((d01,d0))
2171 arg1_c.resize((d01,d1))
2172 out=numpy.zeros((d0,d1),numpy.float64)
2173 for i0 in range(d0):
2174 for i1 in range(d1):
2175 out[i0,i1]=numpy.sum(arg0_c[:,i0]*arg1_c[:,i1])
2176 out.resize(sh0[axis_offset:]+sh1[axis_offset:])
2177 return out
2178 elif isinstance(arg0,escript.Data):
2179 # this call has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
2180 return escript_generalTransposedTensorProduct(arg0,arg1,axis_offset)
2181
2182 # this should be escript._generalTransposedTensorProduct
2183 def escript_generalTransposedTensorProduct(arg0,arg1,axis_offset):
2184 "arg0 and arg1 are both Data objects but not necessarily on the same function space. They could be identical!!!"
2185 return C_GeneralTensorProduct(arg0, arg1, axis_offset, 1)
2186
2187 def matrix_transposed_mult(arg0,arg1):
2188 """
2189 matrix-transposed(matrix) product of the two arguments.
2190
2191 C{out[s0,s1]=Sigma_{r0} arg0[s0,r0]*arg1[s1,r0]}
2192
2193 The function call ``matrix_transposed_mult(arg0,arg1)`` is equivalent to
2194 ``matrix_mult(arg0,transpose(arg1))``.
2195
2196 The last dimensions of ``arg0`` and ``arg1`` must match.
2197
2198 :param arg0: first argument of rank 2
2199 :type arg0: ``numpy.ndarray``, `escript.Data`
2200 :param arg1: second argument of rank 1 or 2
2201 :type arg1: ``numpy.ndarray``, `escript.Data`
2202 :return: the product of ``arg0`` and the transposed of ``arg1`` at each data
2203 point
2204 :rtype: ``numpy.ndarray`` or `escript.Data` depending on the input
2205 :raise ValueError: if the shapes of the arguments are not appropriate
2206 """
2207 sh0=getShape(arg0)
2208 sh1=getShape(arg1)
2209 if not len(sh0)==2 :
2210 raise ValueError,"first argument must have rank 2"
2211 if not len(sh1)==2 and not len(sh1)==1:
2212 raise ValueError,"second argument must have rank 1 or 2"
2213 return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
2214
2215 def tensor_transposed_mult(arg0,arg1):
2216 """
2217 The tensor product of the first and the transpose of the second argument.
2218
2219 For ``arg0`` of rank 2 this is
2220
2221 C{out[s0,s1]=Sigma_{r0} arg0[s0,r0]*arg1[s1,r0]}
2222
2223 and for ``arg0`` of rank 4 this is
2224
2225 C{out[s0,s1,s2,s3]=Sigma_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,s3,r0,r1]}
2226
2227 or
2228
2229 C{out[s0,s1,s2]=Sigma_{r0,r1} arg0[s0,s1,r0,r1]*arg1[s2,r0,r1]}
2230
2231 In the first case the the second dimension of ``arg0`` and ``arg1`` must
2232 match and in the second case the two last dimensions of ``arg0`` must match
2233 the two last dimensions of ``arg1``.
2234
2235 The function call ``tensor_transpose_mult(arg0,arg1)`` is equivalent to
2236 ``tensor_mult(arg0,transpose(arg1))``.
2237
2238 :param arg0: first argument of rank 2 or 4
2239 :type arg0: ``numpy.ndarray``, `escript.Data`
2240 :param arg1: second argument of shape greater of 1 or 2 depending on rank
2241 of ``arg0``
2242 :type arg1: ``numpy.ndarray``, `escript.Data`
2243 :return: the tensor product of the transposed of ``arg0`` and ``arg1`` at
2244 each data point
2245 :rtype: ``numpy.ndarray`` or `escript.Data` depending on the input
2246 """
2247 sh0=getShape(arg0)
2248 sh1=getShape(arg1)
2249 if len(sh0)==2 and ( len(sh1)==2 or len(sh1)==1 ):
2250 return generalTensorTransposedProduct(arg0,arg1,axis_offset=1)
2251 elif len(sh0)==4 and (len(sh1)==2 or len(sh1)==3 or len(sh1)==4):
2252 return generalTensorTransposedProduct(arg0,arg1,axis_offset=2)
2253 else:
2254 raise ValueError,"first argument must have rank 2 or 4"
2255
2256 def generalTensorTransposedProduct(arg0,arg1,axis_offset=0):
2257 """
2258 Generalized tensor product of ``arg0`` and transpose of ``arg1``.
2259
2260 C{out[s,t]=Sigma_r arg0[s,r]*arg1[t,r]}
2261
2262 where
2263 - s runs through ``arg0.Shape[:arg0.ndim-axis_offset]``
2264 - r runs through ``arg0.Shape[arg1.ndim-axis_offset:]``
2265 - t runs through ``arg1.Shape[arg1.ndim-axis_offset:]``
2266
2267 The function call ``generalTensorTransposedProduct(arg0,arg1,axis_offset)``
2268 is equivalent to
2269 ``generalTensorProduct(arg0,transpose(arg1,arg1.ndim-axis_offset),axis_offset)``.
2270
2271 :param arg0: first argument
2272 :type arg0: ``numpy.ndarray``, `escript.Data`, ``float``, ``int``
2273 :param arg1: second argument
2274 :type arg1: ``numpy.ndarray``, `escript.Data`, ``float``, ``int``
2275 :return: the general tensor product of ``arg0`` and ``transpose(arg1)`` at
2276 each data point
2277 :rtype: ``numpy.ndarray`` or `escript.Data` depending on the input
2278 """
2279 if isinstance(arg0,float) and isinstance(arg1,float): return arg1*arg0
2280 arg0,arg1=matchType(arg0,arg1)
2281 # at this stage arg0 and arg1 are both numpy.ndarray or escript.Data
2282 if isinstance(arg0,numpy.ndarray):
2283 if not arg0.shape[arg0.ndim-axis_offset:]==arg1.shape[arg1.ndim-axis_offset:]:
2284 raise ValueError,"dimensions of last %s components in left argument don't match the first %s components in the right argument."%(axis_offset,axis_offset)
2285 arg0_c=arg0.copy()
2286 arg1_c=arg1.copy()
2287 sh0,sh1=arg0.shape,arg1.shape
2288 d0,d1,d01=1,1,1
2289 for i in sh0[:arg0.ndim-axis_offset]: d0*=i
2290 for i in sh1[:arg1.ndim-axis_offset]: d1*=i
2291 for i in sh1[arg1.ndim-axis_offset:]: d01*=i
2292 arg0_c.resize((d0,d01))
2293 arg1_c.resize((d1,d01))
2294 out=numpy.zeros((d0,d1),numpy.float64)
2295 for i0 in range(d0):
2296 for i1 in range(d1):
2297 out[i0,i1]=numpy.sum(arg0_c[i0,:]*arg1_c[i1,:])
2298 out.resize(sh0[:arg0.ndim-axis_offset]+sh1[:arg1.ndim-axis_offset])
2299 return out
2300 elif isinstance(arg0,escript.Data):
2301 # this call has to be replaced by escript._generalTensorProduct(arg0,arg1,axis_offset)
2302 return escript_generalTensorTransposedProduct(arg0,arg1,axis_offset)
2303
2304 # this should be escript._generalTensorTransposedProduct
2305 def escript_generalTensorTransposedProduct(arg0,arg1,axis_offset):
2306 "arg0 and arg1 are both Data objects but not necessarily on the same function space. They could be identical!!!"
2307 return C_GeneralTensorProduct(arg0, arg1, axis_offset, 2)
2308
2309 #=========================================================
2310 # functions dealing with spatial dependency
2311 #=========================================================
2312 def grad(arg,where=None):
2313 """
2314 Returns the spatial gradient of ``arg`` at ``where``.
2315
2316 If ``g`` is the returned object, then
2317
2318 - if ``arg`` is rank 0 ``g[s]`` is the derivative of ``arg`` with respect to
2319 the ``s``-th spatial dimension
2320 - if ``arg`` is rank 1 ``g[i,s]`` is the derivative of ``arg[i]`` with
2321 respect to the ``s``-th spatial dimension
2322 - if ``arg`` is rank 2 ``g[i,j,s]`` is the derivative of ``arg[i,j]`` with
2323 respect to the ``s``-th spatial dimension
2324 - if ``arg`` is rank 3 ``g[i,j,k,s]`` is the derivative of ``arg[i,j,k]``
2325 with respect to the ``s``-th spatial dimension.
2326
2327 :param arg: function of which the gradient is to be calculated. Its rank
2328 has to be less than 3.
2329 :type arg: `escript.Data`
2330 :param where: FunctionSpace in which the gradient is calculated.
2331 If not present or ``None`` an appropriate default is used.
2332 :type where: ``None`` or `escript.FunctionSpace`
2333 :return: gradient of ``arg``
2334 :rtype: `escript.Data`
2335 """
2336 if isinstance(arg,escript.Data):
2337 if where==None:
2338 return arg._grad()
2339 else:
2340 return arg._grad(where)
2341 elif isinstance(arg,Basic):
2342 return symfn.grad(arg, where)
2343 else:
2344 raise TypeError,"grad: Unknown argument type."
2345
2346 def integrate(arg,where=None):
2347 """
2348 Returns the integral of the function ``arg`` over its domain. If ``where`` is
2349 present ``arg`` is interpolated to ``where`` before integration.
2350
2351 :param arg: the function which is integrated
2352 :type arg: `escript.Data`
2353 :param where: FunctionSpace in which the integral is calculated.
2354 If not present or ``None`` an appropriate default is used.
2355 :type where: ``None`` or `escript.FunctionSpace`
2356 :return: integral of ``arg``
2357 :rtype: ``float`` or ``numpy.ndarray``
2358 """
2359 if isinstance(arg,escript.Data):
2360 if not where==None: arg=escript.Data(arg,where)
2361 if arg.getRank()==0:
2362 return arg._integrateToTuple()[0]
2363 else:
2364 return numpy.array(arg._integrateToTuple())
2365 elif isinstance(arg,Basic):
2366 return symfn.integrate(arg, where)
2367 else:
2368 arg2=escript.Data(arg,where)
2369 if arg2.getRank()==0:
2370 return arg2._integrateToTuple()[0]
2371 else:
2372 return numpy.array(arg2._integrateToTuple())
2373
2374 def interpolate(arg,where):
2375 """
2376 Interpolates the function into the `FunctionSpace` ``where``. If the
2377 argument ``arg`` has the requested function space ``where`` no interpolation
2378 is performed and ``arg`` is returned.
2379
2380 :param arg: interpolant
2381 :type arg: `escript.Data`
2382 :param where: `FunctionSpace` to be interpolated to
2383 :type where: `escript.FunctionSpace`
2384 :return: interpolated argument
2385 :rtype: ``escript.Data``
2386 """
2387 if isinstance(arg,escript.Data):
2388 if arg.isEmpty():
2389 return arg
2390 elif where == arg.getFunctionSpace():
2391 return arg
2392 else:
2393 return escript.Data(arg,where)
2394 elif isinstance(arg,Basic):
2395 return symfn.interpolate(arg, where)
2396 else:
2397 return escript.Data(arg,where)
2398
2399 def div(arg,where=None):
2400 """
2401 Returns the divergence of ``arg`` at ``where``.
2402
2403 :param arg: function of which the divergence is to be calculated. Its
2404 shape has to be (d,) where d is the spatial dimension.
2405 :type arg: `escript.Data`
2406 :param where: `FunctionSpace` in which the divergence will be calculated.
2407 If not present or ``None`` an appropriate default is used.
2408 :type where: ``None`` or `escript.FunctionSpace`
2409 :return: divergence of ``arg``
2410 :rtype: `escript.Data`
2411 """
2412 if isinstance(arg,escript.Data):
2413 dim=arg.getDomain().getDim()
2414 if not arg.getShape()==(dim,):
2415 raise ValueError,"div: expected shape is (%s,)"%dim
2416 elif not isinstance(arg, Basic):
2417 raise TypeError,"div: argument type not supported"
2418 return trace(grad(arg,where))
2419
2420 def jump(arg,domain=None):
2421 """
2422 Returns the jump of ``arg`` across the continuity of the domain.
2423
2424 :param arg: argument
2425 :type arg: `escript.Data`
2426 :param domain: the domain where the discontinuity is located. If domain is
2427 not present or equal to ``None`` the domain of ``arg`` is used.
2428 :type domain: ``None`` or `escript.Domain`
2429 :return: jump of ``arg``
2430 :rtype: `escript.Data`
2431 """
2432 if domain==None: domain=arg.getDomain()
2433 return interpolate(arg,escript.FunctionOnContactOne(domain))-interpolate(arg,escript.FunctionOnContactZero(domain))
2434
2435 def L2(arg):
2436 """
2437 Returns the L2 norm of ``arg`` at ``where``.
2438
2439 :param arg: function of which the L2 norm is to be calculated
2440 :type arg: `escript.Data`
2441 :return: L2 norm of ``arg``
2442 :rtype: `float`
2443 :note: L2(arg) is equivalent to ``sqrt(integrate(inner(arg,arg)))``
2444 """
2445 if isinstance(arg,Basic):
2446 return symfn.L2(arg)
2447 return sqrt(integrate(inner(arg,arg)))
2448
2449 def getClosestValue(arg,origin=0):
2450 """
2451 Returns the value in ``arg`` which is closest to origin.
2452
2453 :param arg: function
2454 :type arg: `escript.Data`
2455 :param origin: reference value
2456 :type origin: ``float`` or `escript.Data`
2457 :return: value in ``arg`` closest to origin
2458 :rtype: ``numpy.ndarray``
2459 """
2460 return arg.getValueOfGlobalDataPoint(*(length(arg-origin).minGlobalDataPoint()))
2461
2462 def normalize(arg,zerolength=0):
2463 """
2464 Returns the normalized version of ``arg`` (=``arg/length(arg)``).
2465
2466 :param arg: function
2467 :type arg: `escript.Data`
2468 :param zerolength: relative tolerance for arg == 0
2469 :type zerolength: ``float``
2470 :return: normalized ``arg`` where ``arg`` is non-zero, and zero elsewhere
2471 :rtype: `escript.Data`
2472 """
2473 l=length(arg)
2474 m=whereZero(l,zerolength*Lsup(l))
2475 mm=1-m
2476 return arg*(mm/(l*mm+m))
2477
2478 def deviatoric(arg):
2479 """
2480 Returns the deviatoric version of ``arg``.
2481 """
2482 return arg-(trace(arg)/trace(kronecker(arg.getDomain())))*kronecker(arg.getDomain())
2483
2484 def vol(arg):
2485 """
2486 Returns the volume or area of the oject ``arg``
2487
2488 :param arg: a geometrical object
2489 :type arg: `escript.FunctionSpace` or `escript.Domain`
2490 :rtype: ``float``
2491 """
2492 if isinstance(arg,escript.Domain): arg=escript.Function(arg)
2493 return integrate(escript.Scalar(1.,arg))
2494
2495 def meanValue(arg):
2496 """
2497 return the mean value of the argument over its domain
2498
2499 :param arg: function
2500 :type arg: `escript.Data`
2501 :return: mean value
2502 :rtype: ``float`` or {numpy.ndarray}
2503 """
2504 fs=arg.getFunctionSpace()
2505 d=fs.getDomain()
2506 if fs == escript.Solution(d) or fs == escript.ContinuousFunction(d):
2507 fs=escript.Function(d)
2508 if fs == escript.ReducedSolution(d) or fs == escript.ReducedContinuousFunction(d):
2509 fs=escript.ReducedFunction(d)
2510 a=vol(fs)
2511 if a == 0:
2512 raise ValueError,"FunctionSpace %s with zero volume."%str(fs)
2513 return integrate(arg,fs)/a
2514
2515 def diameter(domain):
2516 """
2517 Returns the diameter of a domain.
2518
2519 :param domain: a domain
2520 :type domain: `escript.Domain`
2521 :rtype: ``float``
2522 """
2523 out=0
2524 for v in boundingBox(domain): out+=(v[1]-v[0])**2
2525 return sqrt(out)
2526
2527 def boundingBox(domain):
2528 """
2529 Returns the bounding box of a domain
2530
2531 :param domain: a domain
2532 :type domain: `escript.Domain`
2533 :return: bounding box of the domain
2534 :rtype: ``list`` of pairs of ``float``
2535 """
2536 x=domain.getX()
2537 out=[]
2538 for i in xrange(domain.getDim()):
2539 x_i=x[i]
2540 out.append((inf(x_i),sup(x_i)))
2541 return out
2542
2543 def longestEdge(domain):
2544 """
2545 Returns the length of the longest edge of the domain
2546
2547 :param domain: a domain
2548 :type domain: `escript.Domain`
2549 :return: longest edge of the domain parallel to the Cartisean axis
2550 :rtype: ``float``
2551 """
2552 return max([v[1]-v[0] for v in boundingBox(domain) ])
2553
2554 def mkDir(*pathname):
2555 """
2556 creates a directory of name ``pathname`` if the directory does not exist.
2557
2558 :param pathname: valid path name
2559 :type pathname: ``str`` or ``sequence of strings``
2560 :note: The method is MPI safe.
2561 """
2562 errno=0
2563 p_fail=None
2564 if getMPIRankWorld()==0:
2565 for p in pathname:
2566 if os.path.exists(p):
2567 if not os.path.isdir(p):
2568 errno=2
2569 p_fail=p
2570 else:
2571 try:
2572 os.makedirs(p)
2573 except Exception, e:
2574 errno=1
2575 p_fail=p
2576
2577 errno=getMPIWorldMax(errno)
2578 if errno>0:
2579 if errno==2:
2580 if p_fail == None:
2581 raise IOError,"Unable to create directory."
2582 else:
2583 raise IOError,"Unable to create directory %s. It already exists and is not a directory."%p_fail
2584 elif e==None:
2585 if p_fail == None:
2586 raise IOError,"Unable to create directory."
2587 else:
2588 raise IOError,"Unable to create directory %s."%p_fail
2589 else:
2590 if hasattr(e,"message"):
2591 raise IOError,e.message
2592 else:
2593 if p_fail == None:
2594 raise IOError,"Unable to create directory."
2595 else:
2596 raise IOError,"Unable to create directory %s."%p_fail
2597
2598 class FileWriter(object):
2599 """
2600 Interface to write data to a file. In essence this class wrappes the standard ``file`` object to write data that are global in MPI
2601 to a file. In fact, data are writen on the processor with MPI rank 0 only. It is recommended to use ``FileWriter`` rather than ``open`` in order to write
2602 code that is running with as well as with MPI. It is safe to use ``open`` onder MPI to read data which are global under MPI.
2603 :var name: name of file
2604 :var mode: access mode (='w' or ='a')
2605 :var closed: True to indicate closed file
2606 :var newlines: line seperator
2607 """
2608 def __init__(self,fn,append=False,createLocalFiles=False):
2609 """
2610 Opens a file of name ``fn`` for writing. If running under MPI only the first processor with rank==0
2611 will open the file and write to it. If ``createLocalFiles`` each individual processor will create a file
2612 where for any processor with rank>0 the file name is extended by its rank. This option is normally only used for
2613 debug purposes.
2614
2615 :param fn: filename.
2616 :type fn: ``str``
2617 :param append: switches on the creation of local files.
2618 :type append: ``bool``
2619 :param createLocalFiles: switches on the creation of local files.
2620 :type createLocalFiles: ``bool``
2621 """
2622 errno=0
2623 self.name=fn
2624 if append:
2625 self.mode='a'
2626 else:
2627 self.mode='w'
2628 self.__file=None
2629 self.closed=False
2630 self.newlines=os.linesep
2631 e=None
2632 # if not the master:
2633 if getMPIRankWorld()>0:
2634 if createLocalFiles:
2635 fn2=fn+".%s"%getMPIRankWorld()
2636 try:
2637 self.__file=open(fn2,self.mode)
2638 except Exception, e:
2639 errno=1
2640 else:
2641 try:
2642 self.__file=open(fn,self.mode)
2643 except Exception, e:
2644 errno=1
2645 self.__handelerror(errno,e,"opening")
2646
2647 def __handelerror(self,errno,e,operation):
2648 errno=getMPIWorldMax(errno)
2649 if errno>0:
2650 if e==None:
2651 raise IOError,"Unable to access file %s in mode %s for %s."%(self.name,self.mode,operation)
2652 else:
2653 raise IOError(str(e))
2654
2655 def close(self):
2656 """
2657 Closes the file
2658 """
2659 errno=0
2660 e=None
2661 try:
2662 if not self.__file == None:
2663 self.__file.close()
2664 except Exception, e:
2665 errno=1
2666 self.__handelerror(errno,e,"closing")
2667 self.closed=True
2668
2669 def flush(self):
2670 """
2671 Flush the internal I/O buffer.
2672 """
2673 errno=0
2674 e=None
2675 try:
2676 if not self.__file == None:
2677 self.__file.flush()
2678 except Exception, e:
2679 errno=1
2680 self.__handelerror(errno,e,"flushing")
2681
2682 def write(self,txt):
2683 """
2684 Write string ``txt`` to file.
2685
2686 :param txt: string ``txt`` to be written to file
2687 :type txt: ``str``
2688 """
2689 errno=0
2690 e=None
2691 try:
2692 if not self.__file == None:
2693 self.__file.write(txt)
2694 except Exception, e:
2695 errno=1
2696 self.__handelerror(errno,e,"writing")
2697
2698 def writelines(self, txts):
2699 """
2700 Write the list ``txt`` of strings to the file.
2701
2702 :param txts: sequense of strings to be written to file
2703 :type txts: any iterable object producing strings
2704 :note: Note that newlines are not added. This method is equivalent to call write() for each string.
2705 """
2706 errno=0
2707 e=None
2708 try:
2709 if not self.__file == None:
2710 self.__file.writelines(txts)
2711 except Exception, e:
2712 errno=1
2713 self.__handelerror(errno,e,"writing strings")
2714
2715 def reorderComponents(arg,index):
2716 """
2717 Resorts the components of ``arg`` according to index.
2718
2719 """
2720 raise NotImplementedError
2721
2722 def showEscriptParams():
2723 """
2724 Displays the parameters escript recognises with an explanation.
2725 """
2726 p=listEscriptParams()
2727 for name,desc in p:
2728 print name+':\t'+desc
2729
2730 #Lazy related things
2731 #These are just wrappers
2732 def resolve(arg):
2733 """
2734 Returns the value of arg resolved.
2735 """
2736 if not isinstance(arg,Data):
2737 raise TypeError, "Can only resolve Data."
2738 if arg.isLazy():
2739 arg.resolve()
2740 return arg
2741
2742 def delay(arg):
2743 """
2744 Returns a lazy version of arg
2745 """
2746 if not isinstance(arg,Data):
2747 raise TypeError, "Can only delay Data."
2748 return arg.delay()
2749
2750 def positive(arg):
2751 """
2752 returns the positive part of arg
2753 """
2754 return (abs(arg)+arg)/2.
2755
2756 def negative(arg):
2757 """
2758 returns the negative part of arg
2759 """
2760 return (arg-abs(arg))/2.
2761
2762 def safeDiv(arg0, arg1, rtol=None):
2763 """
2764 returns arg0/arg1 but return 0 where arg1 is (almost) zero
2765 """
2766 if rtol==None:
2767 m1=whereZero(arg1,tol=0)
2768 else:
2769 m1=whereZero(arg1,tol=None, rtol=rtol)
2770 return arg0/(arg1+m1)*(1-m1)
2771
2772 def condEval(f, tval, fval):
2773 """
2774 Wrapper to allow non-data objects to be used.
2775 """
2776 if not isinstance(tval,Data) and not isinstance(fval,Data):
2777 raise TypeError, "At least one of the alternatives must be a Data object."
2778 if isinstance(tval,Data) and isinstance(fval, Data):
2779 return _condEval(f,tval,fval)
2780 if not isinstance(fval, Data):
2781 return _condEval(f, tval, Data(fval, tval.getShape(), tval.getFunctionSpace()))
2782 return _condEval(f, Data(fval, fval.getShape(), fval.getFunctionSpace()), fval )
2783

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26