Merge pull request matplotlib#2363 from GBillotey/new_test_triinterpo…

…late

[bug correction] trirefine is now independant of triangulation numbering

lib/matplotlib/tests/test_triangulation.py

@@ -829,7 +829,7 @@ def power(x, a):
 
 
 def test_trirefine():
-    # test subdiv=2 refinement
+    # Testing subdiv=2 refinement
     n = 3
     subdiv = 2
     x = np.linspace(-1., 1., n+1)
@@ -854,7 +854,7 @@ def test_trirefine():
                     np.around(x_refi*(2.5+y_refi), 8))
     assert_array_equal(ind1d, True)
 
-    # tests the mask of the refined triangulation
+    # Testing the mask of the refined triangulation
     refi_mask = refi_triang.mask
     refi_tri_barycenter_x = np.sum(refi_triang.x[refi_triang.triangles],
                                    axis=1)/3.
@@ -866,6 +866,23 @@ def test_trirefine():
     refi_tri_mask = triang.mask[refi_tri_indices]
     assert_array_equal(refi_mask, refi_tri_mask)
 
+    # Testing that the numbering of triangles does not change the
+    # interpolation result.
+    x = np.asarray([0.0, 1.0, 0.0, 1.0])
+    y = np.asarray([0.0, 0.0, 1.0, 1.0])
+    triang = [mtri.Triangulation(x, y, [[0, 1, 3], [3, 2, 0]]),
+              mtri.Triangulation(x, y, [[0, 1, 3], [2, 0, 3]])]
+    z = np.sqrt((x-0.3)*(x-0.3) + (y-0.4)*(y-0.4))
+    # Refining the 2 triangulations and reordering the points
+    xyz_data = []
+    for i in range(2):
+        refiner = mtri.UniformTriRefiner(triang[i])
+        refined_triang, refined_z = refiner.refine_field(z, subdiv=1)
+        xyz = np.dstack((refined_triang.x, refined_triang.y, refined_z))[0]
+        xyz = xyz[np.lexsort((xyz[:, 1], xyz[:, 0]))]
+        xyz_data += [xyz]
+    assert_array_almost_equal(xyz_data[0], xyz_data[1])
+
 
 def meshgrid_triangles(n):
     """

lib/matplotlib/tri/triinterpolate.py

@@ -676,11 +676,20 @@ class _ReducedHCT_Element():
                            [1., 1., 1.], [1., 0., 0.], [-2.,  0., -1.]])
 
     # 3) Loads Gauss points & weights on the 3 sub-_triangles for P2
-    #    exact integral - points at the middle of subtriangles apex
-    gauss_pts = np.array([[0.5, 0.5, 0.], [0.5, 0., 0.5], [0., 0.5, 0.5],
-                          [4./6., 1./6., 1./6.], [1./6., 4./6., 1./6.],
-                          [1./6., 1./6., 4./6.]])
-    gauss_w = np.array([1./9., 1./9., 1./9., 2./9., 2./9., 2./9.])
+    #    exact integral - 3 points on each subtriangles.
+    # NOTE: as the 2nd derivative is discontinuous , we really need those 9
+    # points!
+    n_gauss = 9
+    gauss_pts = np.array([[13./18.,  4./18.,  1./18.],
+                          [ 4./18., 13./18.,  1./18.],
+                          [ 7./18.,  7./18.,  4./18.],
+                          [ 1./18., 13./18.,  4./18.],
+                          [ 1./18.,  4./18., 13./18.],
+                          [ 4./18.,  7./18.,  7./18.],
+                          [ 4./18.,  1./18., 13./18.],
+                          [13./18.,  1./18.,  4./18.],
+                          [ 7./18.,  4./18.,  7./18.]], dtype=np.float64)
+    gauss_w = np.ones([9], dtype=np.float64) / 9.
 
     #  4) Stiffness matrix for curvature energy
     E = np.array([[1., 0., 0.], [0., 1., 0.], [0., 0., 2.]])
@@ -755,16 +764,16 @@ def get_function_derivatives(self, alpha, J, ecc, dofs):
         y_sq = y*y
         z_sq = z*z
         dV = _to_matrix_vectorized([
-            [-3.*x_sq, -3.*x_sq],
-            [3.*y_sq, 0.],
-            [0., 3.*z_sq],
-            [-2.*x*z, -2.*x*z+x_sq],
-            [-2.*x*y+x_sq, -2.*x*y],
-            [2.*x*y-y_sq, -y_sq],
-            [2.*y*z, y_sq],
-            [z_sq, 2.*y*z],
-            [-z_sq, 2.*x*z-z_sq],
-            [x*z-y*z, x*y-y*z]])
+            [    -3.*x_sq,     -3.*x_sq],
+            [     3.*y_sq,           0.],
+            [          0.,      3.*z_sq],
+            [     -2.*x*z, -2.*x*z+x_sq],
+            [-2.*x*y+x_sq,      -2.*x*y],
+            [ 2.*x*y-y_sq,        -y_sq],
+            [      2.*y*z,         y_sq],
+            [        z_sq,       2.*y*z],
+            [       -z_sq,  2.*x*z-z_sq],
+            [     x*z-y*z,      x*y-y*z]])
         # Puts back dV in first apex basis
         dV = _prod_vectorized(dV, _extract_submatrices(
             self.rotate_dV, subtri, block_size=2, axis=0))
@@ -831,16 +840,16 @@ def get_d2Sidksij2(self, alpha, ecc):
         y = ksi[:, 1, 0]
         z = ksi[:, 2, 0]
         d2V = _to_matrix_vectorized([
-            [6.*x, 6.*x, 6.*x],
-            [6.*y, 0., 0.],
-            [0., 6.*z, 0.],
-            [2.*z, 2.*z-4.*x, 2.*z-2.*x],
-            [2.*y-4.*x, 2.*y, 2.*y-2.*x],
-            [2.*x-4.*y, 0., -2.*y],
-            [2.*z, 0., 2.*y],
-            [0., 2.*y, 2.*z],
-            [0., 2.*x-4.*z, -2.*z],
-            [-2.*z, -2.*y, x-y-z]])
+            [     6.*x,      6.*x,      6.*x],
+            [     6.*y,        0.,        0.],
+            [       0.,      6.*z,        0.],
+            [     2.*z, 2.*z-4.*x, 2.*z-2.*x],
+            [2.*y-4.*x,      2.*y, 2.*y-2.*x],
+            [2.*x-4.*y,        0.,     -2.*y],
+            [     2.*z,        0.,      2.*y],
+            [       0.,      2.*y,      2.*z],
+            [       0., 2.*x-4.*z,     -2.*z],
+            [    -2.*z,     -2.*y,     x-y-z]])
         # Puts back d2V in first apex basis
         d2V = _prod_vectorized(d2V, _extract_submatrices(
             self.rotate_d2V, subtri, block_size=3, axis=0))
@@ -887,11 +896,11 @@ def get_bending_matrices(self, J, ecc):
         H_rot, area = self.get_Hrot_from_J(J, return_area=True)
 
         # 3) Computes stiffness matrix
-        # Integration according to Gauss P2 rule for each subtri.
+        # Gauss quadrature.
         K = np.zeros([n, 9, 9], dtype=np.float64)
         weights = self.gauss_w
         pts = self.gauss_pts
-        for igauss in range(6):
+        for igauss in range(self.n_gauss):
             alpha = np.tile(pts[igauss, :], n).reshape(n, 3)
             alpha = np.expand_dims(alpha, 3)
             weight = weights[igauss]
@@ -916,7 +925,8 @@ def get_Hrot_from_J(self, J, return_area=False):
 
         Returns
         -------
-        Returns H_rot used to rotate Hessian from local to global coordinates.
+        Returns H_rot used to rotate Hessian from local basis of first apex,
+        to global coordinates.
         if *return_area* is True, returns also the triangle area (0.5*det(J))
         """
         # Here we try to deal with the simpliest colinear cases ; a null