Merged in no-numpy-matrix (pull request #50)

Remove reliance on numpy.matrix

Approved-by: Chris Richardson <[email protected]>

Author: wence-

FIAT/hdiv_trace.py

@@ -348,19 +348,13 @@ def barycentric_coordinates(points, vertices):
     """
 
     # Form mapping matrix
-    last = np.asarray(vertices[-1])
-    T = np.matrix([np.array(v) - last for v in vertices[:-1]]).T
+    T = (np.asarray(vertices[:-1]) - vertices[-1]).T
     invT = np.linalg.inv(T)
 
-    # Compute the barycentric coordinates for all points
-    coords = []
-    for p in points:
-        y = np.asarray(p) - last
-        bary = invT.dot(y.T)
-        bary = [bary[(0, i)] for i in range(len(y))]
-        bary += [1.0 - sum(bary)]
-        coords.append(bary)
-    return coords
+    points = np.asarray(points)
+    bary = np.einsum("ij,kj->ki", invT, (points - vertices[-1]))
+    last = (1 - bary.sum(axis=1))
+    return np.concatenate([bary, last[..., np.newaxis]], axis=1)
 
 
 def map_from_reference_facet(point, vertices):

FIAT/nedelec_second_kind.py

@@ -160,11 +160,11 @@ def _generate_face_dofs(self, cell, degree, offset):
 
             # Map Phis -> phis (reference values to physical values)
             J = Q_face.jacobian()
-            scale = 1.0 / numpy.sqrt(numpy.linalg.det(J.transpose() * J))
+            scale = 1.0 / numpy.sqrt(numpy.linalg.det(numpy.dot(J.T, J)))
             phis = numpy.ndarray((d, num_quad_points))
             for i in range(num_rts):
                 for q in range(num_quad_points):
-                    phi_i_q = scale * J * numpy.matrix(Phis[i, :, q]).transpose()
+                    phi_i_q = scale * numpy.dot(J, Phis[numpy.newaxis, i, :, q].T)
                     for j in range(d):
                         phis[j, q] = phi_i_q[j]
 

FIAT/quadrature.py

@@ -211,15 +211,12 @@ def __init__(self, face_number, degree):
         # Use tet to map points and weights on the appropriate face
         vertices = [numpy.array(list(vertex)) for vertex in vertices]
         x0 = vertices[0]
-        J = numpy.matrix([vertices[1] - x0, vertices[2] - x0]).transpose()
-        x0 = numpy.matrix(x0).transpose()
+        J = numpy.vstack([vertices[1] - x0, vertices[2] - x0]).T
         # This is just a very numpyfied way of writing J*p + x0:
-        F = lambda p: \
-            numpy.array(J*numpy.matrix(p).transpose() + x0).flatten()
-        points = numpy.array([F(p) for p in ref_points])
+        points = numpy.einsum("ij,kj->ki", J, ref_points) + x0
 
         # Map weights: multiply reference weights by sqrt(|J^T J|)
-        detJTJ = numpy.linalg.det(J.transpose() * J)
+        detJTJ = numpy.linalg.det(numpy.dot(J.T, J))
         weights = numpy.sqrt(detJTJ) * ref_weights
 
         # Initialize super class with new points and weights