Skip to content

[bug correction] trirefine is now independant of triangulation numbering #2363

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Sep 3, 2013
Merged

[bug correction] trirefine is now independant of triangulation numbering #2363

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Binary file modified lib/matplotlib/tests/baseline_images/test_triangulation/tri_smooth_contouring.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified lib/matplotlib/tests/baseline_images/test_triangulation/tri_smooth_gradient.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
21 changes: 19 additions & 2 deletions lib/matplotlib/tests/test_triangulation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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.
Expand All @@ -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):
"""
Expand Down
66 changes: 38 additions & 28 deletions lib/matplotlib/tri/triinterpolate.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.]])
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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]
Expand All @@ -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
Expand Down