diff --git a/lib/matplotlib/tests/baseline_images/test_triangulation/tri_smooth_contouring.png b/lib/matplotlib/tests/baseline_images/test_triangulation/tri_smooth_contouring.png index ce81a1f55cb6..f9fec56be0c6 100644 Binary files a/lib/matplotlib/tests/baseline_images/test_triangulation/tri_smooth_contouring.png and b/lib/matplotlib/tests/baseline_images/test_triangulation/tri_smooth_contouring.png differ diff --git a/lib/matplotlib/tests/baseline_images/test_triangulation/tri_smooth_gradient.png b/lib/matplotlib/tests/baseline_images/test_triangulation/tri_smooth_gradient.png index aff9d7ade540..0dd2d6f0f6ba 100644 Binary files a/lib/matplotlib/tests/baseline_images/test_triangulation/tri_smooth_gradient.png and b/lib/matplotlib/tests/baseline_images/test_triangulation/tri_smooth_gradient.png differ diff --git a/lib/matplotlib/tests/test_triangulation.py b/lib/matplotlib/tests/test_triangulation.py index 1847c912dbc9..a5bbe6f65075 100644 --- a/lib/matplotlib/tests/test_triangulation.py +++ b/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): """ diff --git a/lib/matplotlib/tri/triinterpolate.py b/lib/matplotlib/tri/triinterpolate.py index 2026711c51e5..aa2c9c92ff91 100644 --- a/lib/matplotlib/tri/triinterpolate.py +++ b/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