Skip to content

Commit 300b393

Browse files
GBilloteyGBillotey
committed
[bug correction] triinterpolate is independant of triangulation numbering.
1 parent c556204 commit 300b393

File tree

4 files changed

+57
-30
lines changed

4 files changed

+57
-30
lines changed

‎lib/matplotlib/tests/baseline_images/test_triangulation/tri_smooth_contouring.png

-1.79 KB
Loading

‎lib/matplotlib/tests/baseline_images/test_triangulation/tri_smooth_gradient.png

-1.73 KB
Loading

‎lib/matplotlib/tests/test_triangulation.py

+19-2
Original file line numberDiff line numberDiff line change
@@ -829,7 +829,7 @@ def power(x, a):
829829

830830

831831
def test_trirefine():
832-
# test subdiv=2 refinement
832+
# Testing subdiv=2 refinement
833833
n = 3
834834
subdiv = 2
835835
x = np.linspace(-1., 1., n+1)
@@ -854,7 +854,7 @@ def test_trirefine():
854854
np.around(x_refi*(2.5+y_refi), 8))
855855
assert_array_equal(ind1d, True)
856856

857-
# tests the mask of the refined triangulation
857+
# Testing the mask of the refined triangulation
858858
refi_mask = refi_triang.mask
859859
refi_tri_barycenter_x = np.sum(refi_triang.x[refi_triang.triangles],
860860
axis=1)/3.
@@ -866,6 +866,23 @@ def test_trirefine():
866866
refi_tri_mask = triang.mask[refi_tri_indices]
867867
assert_array_equal(refi_mask, refi_tri_mask)
868868

869+
# Testing that the numbering of triangles does not change the
870+
# interpolation result.
871+
x = np.asarray([0.0, 1.0, 0.0, 1.0])
872+
y = np.asarray([0.0, 0.0, 1.0, 1.0])
873+
triang = [mtri.Triangulation(x, y, [[0, 1, 3], [3, 2, 0]]),
874+
mtri.Triangulation(x, y, [[0, 1, 3], [2, 0, 3]])]
875+
z = np.sqrt((x-0.3)*(x-0.3) + (y-0.4)*(y-0.4))
876+
# Refining the 2 triangulations and reordering the points
877+
xyz_data = []
878+
for i in range(2):
879+
refiner = mtri.UniformTriRefiner(triang[i])
880+
refined_triang, refined_z = refiner.refine_field(z, subdiv=1)
881+
xyz = np.dstack((refined_triang.x, refined_triang.y, refined_z))[0]
882+
xyz = xyz[np.lexsort((xyz[:, 1], xyz[:, 0]))]
883+
xyz_data += [xyz]
884+
assert_array_almost_equal(xyz_data[0], xyz_data[1])
885+
869886

870887
def meshgrid_triangles(n):
871888
"""

‎lib/matplotlib/tri/triinterpolate.py

+38-28
Original file line numberDiff line numberDiff line change
@@ -676,11 +676,20 @@ class _ReducedHCT_Element():
676676
[1., 1., 1.], [1., 0., 0.], [-2., 0., -1.]])
677677

678678
# 3) Loads Gauss points & weights on the 3 sub-_triangles for P2
679-
# exact integral - points at the middle of subtriangles apex
680-
gauss_pts = np.array([[0.5, 0.5, 0.], [0.5, 0., 0.5], [0., 0.5, 0.5],
681-
[4./6., 1./6., 1./6.], [1./6., 4./6., 1./6.],
682-
[1./6., 1./6., 4./6.]])
683-
gauss_w = np.array([1./9., 1./9., 1./9., 2./9., 2./9., 2./9.])
679+
# exact integral - 3 points on each subtriangles.
680+
# NOTE: as the 2nd derivative is discontinuous , we really need those 9
681+
# points!
682+
n_gauss = 9
683+
gauss_pts = np.array([[13./18., 4./18., 1./18.],
684+
[ 4./18., 13./18., 1./18.],
685+
[ 7./18., 7./18., 4./18.],
686+
[ 1./18., 13./18., 4./18.],
687+
[ 1./18., 4./18., 13./18.],
688+
[ 4./18., 7./18., 7./18.],
689+
[ 4./18., 1./18., 13./18.],
690+
[13./18., 1./18., 4./18.],
691+
[ 7./18., 4./18., 7./18.]], dtype=np.float64)
692+
gauss_w = np.ones([9], dtype=np.float64) / 9.
684693

685694
# 4) Stiffness matrix for curvature energy
686695
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):
755764
y_sq = y*y
756765
z_sq = z*z
757766
dV = _to_matrix_vectorized([
758-
[-3.*x_sq, -3.*x_sq],
759-
[3.*y_sq, 0.],
760-
[0., 3.*z_sq],
761-
[-2.*x*z, -2.*x*z+x_sq],
762-
[-2.*x*y+x_sq, -2.*x*y],
763-
[2.*x*y-y_sq, -y_sq],
764-
[2.*y*z, y_sq],
765-
[z_sq, 2.*y*z],
766-
[-z_sq, 2.*x*z-z_sq],
767-
[x*z-y*z, x*y-y*z]])
767+
[ -3.*x_sq, -3.*x_sq],
768+
[ 3.*y_sq, 0.],
769+
[ 0., 3.*z_sq],
770+
[ -2.*x*z, -2.*x*z+x_sq],
771+
[-2.*x*y+x_sq, -2.*x*y],
772+
[ 2.*x*y-y_sq, -y_sq],
773+
[ 2.*y*z, y_sq],
774+
[ z_sq, 2.*y*z],
775+
[ -z_sq, 2.*x*z-z_sq],
776+
[ x*z-y*z, x*y-y*z]])
768777
# Puts back dV in first apex basis
769778
dV = _prod_vectorized(dV, _extract_submatrices(
770779
self.rotate_dV, subtri, block_size=2, axis=0))
@@ -831,16 +840,16 @@ def get_d2Sidksij2(self, alpha, ecc):
831840
y = ksi[:, 1, 0]
832841
z = ksi[:, 2, 0]
833842
d2V = _to_matrix_vectorized([
834-
[6.*x, 6.*x, 6.*x],
835-
[6.*y, 0., 0.],
836-
[0., 6.*z, 0.],
837-
[2.*z, 2.*z-4.*x, 2.*z-2.*x],
838-
[2.*y-4.*x, 2.*y, 2.*y-2.*x],
839-
[2.*x-4.*y, 0., -2.*y],
840-
[2.*z, 0., 2.*y],
841-
[0., 2.*y, 2.*z],
842-
[0., 2.*x-4.*z, -2.*z],
843-
[-2.*z, -2.*y, x-y-z]])
843+
[ 6.*x, 6.*x, 6.*x],
844+
[ 6.*y, 0., 0.],
845+
[ 0., 6.*z, 0.],
846+
[ 2.*z, 2.*z-4.*x, 2.*z-2.*x],
847+
[2.*y-4.*x, 2.*y, 2.*y-2.*x],
848+
[2.*x-4.*y, 0., -2.*y],
849+
[ 2.*z, 0., 2.*y],
850+
[ 0., 2.*y, 2.*z],
851+
[ 0., 2.*x-4.*z, -2.*z],
852+
[ -2.*z, -2.*y, x-y-z]])
844853
# Puts back d2V in first apex basis
845854
d2V = _prod_vectorized(d2V, _extract_submatrices(
846855
self.rotate_d2V, subtri, block_size=3, axis=0))
@@ -887,11 +896,11 @@ def get_bending_matrices(self, J, ecc):
887896
H_rot, area = self.get_Hrot_from_J(J, return_area=True)
888897

889898
# 3) Computes stiffness matrix
890-
# Integration according to Gauss P2 rule for each subtri.
899+
# Gauss quadrature.
891900
K = np.zeros([n, 9, 9], dtype=np.float64)
892901
weights = self.gauss_w
893902
pts = self.gauss_pts
894-
for igauss in range(6):
903+
for igauss in range(self.n_gauss):
895904
alpha = np.tile(pts[igauss, :], n).reshape(n, 3)
896905
alpha = np.expand_dims(alpha, 3)
897906
weight = weights[igauss]
@@ -916,7 +925,8 @@ def get_Hrot_from_J(self, J, return_area=False):
916925
917926
Returns
918927
-------
919-
Returns H_rot used to rotate Hessian from local to global coordinates.
928+
Returns H_rot used to rotate Hessian from local basis of first apex,
929+
to global coordinates.
920930
if *return_area* is True, returns also the triangle area (0.5*det(J))
921931
"""
922932
# Here we try to deal with the simpliest colinear cases ; a null

0 commit comments

Comments
 (0)