/usr/lib/python2.7/dist-packages/Scientific/Geometry/TensorModule.py is in python-scientific 2.9.4-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 | # This module defines 3d geometrical tensors with the standard
# operations on them. The elements are stored in an array.
#
# Written by Konrad Hinsen <hinsen@cnrs-orleans.fr>
# last revision: 2006-11-23
#
from Scientific import N; Numeric = N
class Tensor:
"""Tensor in 3D space
Tensors support the usual arithmetic operations
('t1', 't2': tensors, 'v': vector, 's': scalar):
- 't1+t2' (addition)
- 't1-t2' (subtraction)
- 't1*t2' (tensorial (outer) product)
- 't1*v' (contraction with a vector, same as t1.dot(v.asTensor()))
- 's*t1', 't1*s' (multiplication with a scalar)
- 't1/s' (division by a scalar)
The coordinates can be extracted by indexing; a tensor of rank N
can be indexed like an array of dimension N.
Tensors are I{immutable}, i.e. their elements cannot be changed.
Tensor elements can be any objects on which the standard
arithmetic operations are defined. However, eigenvalue calculation
is supported only for float elements.
"""
is_tensor = 1
def __init__(self, elements, nocheck = None):
"""
@param elements: 2D array or nested list specifying the nine
tensor components
[[xx, xy, xz], [yx, yy, yz], [zx, zy, zz]]
@type elements: C{Numeric.array} or C{list}
"""
self.array = N.array(elements)
if nocheck is None:
if not N.logical_and.reduce(N.equal(N.array(self.array.shape), 3)):
raise ValueError('Tensor must have length 3 along any axis')
self.rank = len(self.array.shape)
def __repr__(self):
return 'Tensor(' + str(self) + ')'
def __str__(self):
return str(self.array)
def __add__(self, other):
return Tensor(self.array+other.array, 1)
__radd__ = __add__
def __neg__(self):
return Tensor(-self.array, 1)
def __sub__(self, other):
return Tensor(self.array-other.array, 1)
def __rsub__(self, other):
return Tensor(other.array-self.array, 1)
def __mul__(self, other):
from Scientific import Geometry
if isTensor(other):
a = self.array[self.rank*(slice(None),)+(N.NewAxis,)]
b = other.array[other.rank*(slice(None),)+(N.NewAxis,)]
return Tensor(N.innerproduct(a, b), 1)
elif Geometry.isVector(other):
return other.__rmul__(self)
else:
return Tensor(self.array*other, 1)
def __rmul__(self, other):
return Tensor(self.array*other, 1)
def __div__(self, other):
if isTensor(other):
raise TypeError("Can't divide by a tensor")
else:
return Tensor(self.array/(1.*other), 1)
__truediv__ = __div__
def __rdiv__(self, other):
raise TypeError("Can't divide by a tensor")
def __cmp__(self, other):
if not isTensor(other):
return NotImplemented
if self.rank != other.rank:
return 1
else:
return not N.logical_and.reduce(
N.equal(self.array, other.array).flat)
def __len__(self):
return 3
def __getitem__(self, index):
elements = self.array[index]
if type(elements) == type(self.array):
return Tensor(elements)
else:
return elements
def asVector(self):
"""
@returns: an equivalent vector object
@rtype: L{Scientific.Geometry.Vector}
@raises ValueError: if rank > 1
"""
from Scientific import Geometry
if self.rank == 1:
return Geometry.Vector(self.array)
else:
raise ValueError('rank > 1')
def dot(self, other):
"""
@returns: the contraction with other
@rtype: L{Tensor}
"""
if isTensor(other):
a = self.array
b = N.transpose(other.array, range(1, other.rank)+[0])
return Tensor(N.innerproduct(a, b), 1)
else:
return Tensor(self.array*other, 1)
def diagonal(self, axis1=0, axis2=1):
if self.rank == 2:
return Tensor([self.array[0,0], self.array[1,1], self.array[2,2]])
else:
if axis2 < axis1: axis1, axis2 = axis2, axis1
raise ValueError('Not yet implemented')
def trace(self, axis1=0, axis2=1):
"""
@returns: the trace of the tensor
@rtype: type of tensor elements
@raises ValueError: if rank !=2
"""
if self.rank == 2:
return self.array[0,0]+self.array[1,1]+self.array[2,2]
else:
raise ValueError('Not yet implemented')
def transpose(self):
"""
@returns: the transposed (index reversed) tensor
@rtype: L{Tensor}
"""
return Tensor(N.transpose(self.array))
def symmetricalPart(self):
"""
@returns: the symmetrical part of the tensor
@rtype: L{Tensor}
@raises ValueError: if rank !=2
"""
if self.rank == 2:
return Tensor(0.5*(self.array + \
N.transpose(self.array,
N.array([1,0]))),
1)
else:
raise ValueError('Not yet implemented')
def asymmetricalPart(self):
"""
@returns: the asymmetrical part of the tensor
@rtype: L{Tensor}
@raises ValueError: if rank !=2
"""
if self.rank == 2:
return Tensor(0.5*(self.array - \
N.transpose(self.array,
N.array([1,0]))),
1)
else:
raise ValueError('Not yet implemented')
def eigenvalues(self):
"""
@returns: the eigenvalues of the tensor
@rtype: C{Numeric.array}
@raises ValueError: if rank !=2
"""
if self.rank == 2:
from Scientific.LA import eigenvalues
return eigenvalues(self.array)
else:
raise ValueError('Undefined operation')
def diagonalization(self):
"""
@returns: the eigenvalues of the tensor and a tensor
representing the rotation matrix to the diagonal form
@rtype: (C{Numeric.array}, L{Tensor})
@raises ValueError: if rank !=2
"""
if self.rank == 2:
from Scientific.LA import eigenvectors
ev, vectors = eigenvectors(self.array)
return ev, Tensor(vectors)
else:
raise ValueError, 'Undefined operation'
def inverse(self):
"""
@returns: the inverse of the tensor
@rtype: L{Tensor}
@raises ValueError: if rank !=2
"""
if self.rank == 2:
from Scientific.LA import inverse
return Tensor(inverse(self.array))
else:
raise ValueError('Undefined operation')
# Type check
def isTensor(x):
"""
@returns: C{True} if x is a L{Tensor}
"""
return hasattr(x,'is_tensor')
|