diff --git a/CMakeLists.txt b/CMakeLists.txt index a50a934d..ab49f9ef 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -98,6 +98,7 @@ pybind11_add_module(gpytoolbox_bindings "${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/binding_upper_envelope.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/binding_read_ply.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/binding_write_ply.cpp" + "${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/binding_curved_hessian_intrinsic.cpp" ) pybind11_add_module(gpytoolbox_bindings_copyleft diff --git a/src/cpp/binding_curved_hessian_intrinsic.cpp b/src/cpp/binding_curved_hessian_intrinsic.cpp new file mode 100644 index 00000000..54f3b67d --- /dev/null +++ b/src/cpp/binding_curved_hessian_intrinsic.cpp @@ -0,0 +1,26 @@ +#include +#include +#include +#include +#include +#include +#include + +using namespace Eigen; +namespace py = pybind11; +using EigenDStride = Stride; +template +using EigenDRef = Ref; //allows passing column/row order matrices easily + +void binding_curved_hessian_intrinsic(py::module& m) { + m.def("_curved_hessian_intrinsic_cpp_impl",[](EigenDRef l_sq, + EigenDRef f) + { + Eigen::MatrixXi E, oE; + igl::orient_halfedges(f, E, oE); + Eigen::SparseMatrix Q; + igl::curved_hessian_energy_intrinsic(f, l_sq, E, oE, Q); + return Q; + }); + +} diff --git a/src/cpp/gpytoolbox_bindings_core.cpp b/src/cpp/gpytoolbox_bindings_core.cpp index 09456864..b80cdaa6 100644 --- a/src/cpp/gpytoolbox_bindings_core.cpp +++ b/src/cpp/gpytoolbox_bindings_core.cpp @@ -26,6 +26,7 @@ void binding_remesh_botsch(py::module& m); void binding_upper_envelope(py::module& m); void binding_read_ply(py::module& m); void binding_write_ply(py::module& m); +void binding_curved_hessian_intrinsic(py::module& m); PYBIND11_MODULE(gpytoolbox_bindings, m) { @@ -46,6 +47,7 @@ PYBIND11_MODULE(gpytoolbox_bindings, m) { binding_upper_envelope(m); binding_read_ply(m); binding_write_ply(m); + binding_curved_hessian_intrinsic(m); m.def("help", [&]() {printf("hi"); }); } diff --git a/src/gpytoolbox/__init__.py b/src/gpytoolbox/__init__.py index 6d687a2b..55677f0a 100644 --- a/src/gpytoolbox/__init__.py +++ b/src/gpytoolbox/__init__.py @@ -29,6 +29,7 @@ from .quadtree_laplacian import quadtree_laplacian from .quadtree_boundary import quadtree_boundary from .quadtree_children import quadtree_children +from .grad_intrinsic import grad_intrinsic from .grad import grad from .doublearea import doublearea from .doublearea_intrinsic import doublearea_intrinsic @@ -107,5 +108,7 @@ from .read_dmat import read_dmat from .linear_blend_skinning import linear_blend_skinning from .barycenters import barycenters +from .biharmonic_energy import biharmonic_energy +from .biharmonic_energy_intrinsic import biharmonic_energy_intrinsic from .adjacency_matrix import adjacency_matrix from .connected_components import connected_components diff --git a/src/gpytoolbox/biharmonic_energy.py b/src/gpytoolbox/biharmonic_energy.py new file mode 100644 index 00000000..9b95cbfd --- /dev/null +++ b/src/gpytoolbox/biharmonic_energy.py @@ -0,0 +1,86 @@ +import numpy as np +import scipy as sp +from .halfedge_lengths_squared import halfedge_lengths_squared +from .biharmonic_energy_intrinsic import biharmonic_energy_intrinsic +from .doublearea import doublearea +from .boundary_vertices import boundary_vertices +from .massmatrix import massmatrix +from .grad import grad + + +def biharmonic_energy(V,F, + bc='mixedfem_zero_neumann'): + """Constructs the biharmonic energy matrix Q such that for a per-vertex function u, the discrete biharmonic energy is u'Qu. + + Parameters + ---------- + V : (n,d) numpy array + vertex list of a triangle mesh + F : (m,3) numpy int array + face index list of a triangle mesh + bc : string, optional (default 'mixedfem_zero_neumann') + Which type of discretization and boundary condition to use. + Options are: {'mixedfem_zero_neumann', 'hessian', 'curved_hessian'} + - 'mixedfem_zero_neumann' implements the mixed finite element + discretization from Jacobson et al. 2010. + "Mixed Finite Elements for Variational Surface Modeling". + - 'hessian' implements the Hessian energy from Stein et al. 2018. + "Natural Boundary Conditions for Smoothing in Geometry Processing". + - 'curved_hessian' implements the Hessian energy from Stein et al. + 2020 "A Smoothness Energy without Boundary Distortion for Curved Surfaces" + via libigl C++ binding. + + Returns + ------- + Q : (n,n) scipy csr_matrix + biharmonic energy matrix + + Examples + -------- + ```python + # Mesh in V,F + import gpytoolbox as gpy + Q = gpy.biharmonic_energy(V,F) + ``` + """ + + if bc=='hessian': + return _hessian_energy(V, F) + else: + l_sq = halfedge_lengths_squared(V,F) + return biharmonic_energy_intrinsic(l_sq,F,n=V.shape[0],bc=bc) + + +def _hessian_energy(V, F): + assert F.shape[1]==3, "Only works on triangle meshes." + + # Q = G' * A * D * Mtilde * D' * A * G + n = V.shape[0] + m = F.shape[0] + dim = V.shape[1] + + b = boundary_vertices(F) + i = np.setdiff1d(np.arange(n),b) + ni = len(i) + + a = doublearea(V, F) / 2. + A = sp.sparse.spdiags([np.tile(a, dim)], 0, + m=dim*m, n=dim*m, format='csr') + + M_d = massmatrix(V, F, type='voronoi').diagonal()[i] + M_d_inv = 1. / M_d + M_tilde_inv = sp.sparse.spdiags([np.tile(M_d_inv, dim**2)], 0, + m=ni*dim**2, n=ni*dim**2, format='csr') + + G = grad(V,F) + Gi = G[:,i] + Dlist = [] + for d in range(dim): + Dlist.append(Gi[d*m:(d+1)*m, :]) + Dblock = sp.sparse.bmat([Dlist], format='csr') + D = sp.sparse.block_diag(dim*[Dblock], format='csr') + + Q = G.transpose() * A * D * M_tilde_inv * D.transpose() * A * G + + return Q + diff --git a/src/gpytoolbox/biharmonic_energy_intrinsic.py b/src/gpytoolbox/biharmonic_energy_intrinsic.py new file mode 100644 index 00000000..26c68428 --- /dev/null +++ b/src/gpytoolbox/biharmonic_energy_intrinsic.py @@ -0,0 +1,92 @@ +import numpy as np +import scipy as sp +from .cotangent_laplacian_intrinsic import cotangent_laplacian_intrinsic +from .massmatrix_intrinsic import massmatrix_intrinsic + +def biharmonic_energy_intrinsic(l_sq,F, + n=None, + bc='mixedfem_zero_neumann'): + """Constructs the biharmonic energy matrix Q such that for a per-vertex function u, the discrete biharmonic energy is u'Qu, using only intrinsic information from the mesh. + + Parameters + ---------- + l_sq : (m,3) numpy array + squared halfedge lengths as computed by halfedge_lengths_squared + F : (m,3) numpy int array + face index list of a triangle mesh + n : int, optional (default None) + number of vertices in the mesh. + If absent, will try to infer from F. + bc : string, optional (default 'mixedfem_zero_neumann') + Which type of discretization and boundary condition to use. + Options are: {'mixedfem_zero_neumann', 'hessian', 'curved_hessian'} + - 'mixedfem_zero_neumann' implements the mixed finite element + discretization from Jacobson et al. 2010. + "Mixed Finite Elements for Variational Surface Modeling". + - 'curved_hessian' implements the Hessian energy from Stein et al. + 2020 "A Smoothness Energy without Boundary Distortion for Curved Surfaces" + via libigl C++ binding. + + Returns + ------- + Q : (n,n) scipy csr_matrix + biharmonic energy matrix + + Examples + -------- + ```python + # Mesh in V,F + import gpytoolbox as gpy + l_sq = gpy.halfedge_lengths_squared(V,F) + Q = biharmonic_energy_intrinsic(l_sq,F) + ``` + + """ + + assert F.shape[1] == 3, "Only works on triangle meshes." + assert l_sq.shape == F.shape + assert np.all(l_sq>=0) + + if n==None: + n = np.max(F)+1 + + if bc=='mixedfem_zero_neumann': + Q = _mixedfem_neumann_laplacian_energy(l_sq, F, n) + elif bc=='curved_hessian': + Q = _curved_hessian_energy(l_sq, F, n) + else: + assert False, "Invalid bc" + + return Q + + +def _mixedfem_neumann_laplacian_energy(l_sq, F, n): + # Q = L' * M^{-1} * L + L = cotangent_laplacian_intrinsic(l_sq, F, n=n) + M = massmatrix_intrinsic(l_sq, F, n=n, type='voronoi') + M_inv = M.power(-1) #Sparse matrix inverts componentwise + Q = L.transpose() * M_inv * L + + return Q + + +try: + # Import C++ bindings for curved Hessian + from gpytoolbox_bindings import _curved_hessian_intrinsic_cpp_impl + _CPP_CURVED_HESSIAN_INTRINSIC_AVAILABLE = True +except Exception as e: + _CPP_CURVED_HESSIAN_INTRINSIC_AVAILABLE = False + +def _curved_hessian_energy(l_sq, F, n): + assert _CPP_CURVED_HESSIAN_INTRINSIC_AVAILABLE, \ + "C++ bindings for curved Hessian not available." + + Q = _curved_hessian_intrinsic_cpp_impl(l_sq.astype(np.float64), + F.astype(np.int32)) + if Q.shape[0] != n: + Q = sp.sparse.block_diag([ + Q, sp.sparse.csr_matrix((n-Q.shape[0], n-Q.shape[0]))], + format='csr') + + return Q + diff --git a/src/gpytoolbox/grad.py b/src/gpytoolbox/grad.py index e694fc0f..c10e4b05 100644 --- a/src/gpytoolbox/grad.py +++ b/src/gpytoolbox/grad.py @@ -30,6 +30,7 @@ def grad(V,F=None): Notes ----- + Adapted from https://github.com/alecjacobson/gptoolbox/blob/master/mesh/grad.m Examples -------- diff --git a/src/gpytoolbox/grad_intrinsic.py b/src/gpytoolbox/grad_intrinsic.py new file mode 100644 index 00000000..c9f347cb --- /dev/null +++ b/src/gpytoolbox/grad_intrinsic.py @@ -0,0 +1,68 @@ +import numpy as np +import scipy as sp + +from .grad import grad + +# Lifted from https://github.com/alecjacobson/gptoolbox/blob/master/mesh/grad.m + +def grad_intrinsic(l_sq,F, + n=None,): + """Intrinsic finite element gradient matrix + + Given a triangle mesh, computes the finite element gradient matrix assuming + piecewise linear hat function basis. + + Parameters + ---------- + l_sq : (m,3) numpy array + squared halfedge lengths as computed by halfedge_lengths_squared + F : (m,3) numpy int array + face index list of a triangle mesh + n : int, optional (default None) + number of vertices in the mesh. + If absent, will try to infer from F. + + Returns + ------- + G : (2*m,n) scipy sparse.csr_matrix + Sparse FEM gradient matrix. + The first m rows ar ethe gradient with respect to the edge (1,2). + The m rows after that run in a pi/2 rotation counter-clockwise. + + See Also + -------- + cotangent_laplacian. + + Notes + ----- + Adapted from https://github.com/alecjacobson/gptoolbox/blob/master/mesh/grad_intrinsic.m + + Examples + -------- + ```python + import gpytoolbox as gpy + l = gpy.halfedge_lengths_squared(V,F) + G = gpy.grad_intrinsic(l_sq,F) + ``` + """ + + + assert F.shape[1] == 3, "Only works on triangle meshes." + assert l_sq.shape == F.shape + assert np.all(l_sq>=0) + + if n==None: + n = np.max(F)+1 + m = F.shape[0] + + l0 = np.sqrt(l_sq[:,0])[:,None] + gx = (l_sq[:,1][:,None]-l_sq[:,0][:,None]-l_sq[:,2][:,None]) / (-2.*l0) + gy = np.sqrt(l_sq[:,2][:,None] - gx**2) + V2 = np.block([[gx,gy], [np.zeros((m,2))], [l0, np.zeros((m,1))]]) + F2 = np.reshape(np.arange(3*m), (m,3), order='F') + G2 = grad(V2,F2) + P = sp.sparse.csr_matrix((np.ones(3*m),(F2.ravel(),F.ravel())), + shape=(3*m,n)) + G = G2*P + + return G diff --git a/test/test_biharmonic_energy.py b/test/test_biharmonic_energy.py new file mode 100644 index 00000000..431eea48 --- /dev/null +++ b/test/test_biharmonic_energy.py @@ -0,0 +1,130 @@ +from .context import gpytoolbox as gpy +from .context import numpy as np +from .context import unittest +import scipy as sp + +class TestBiharmonicEnergy(unittest.TestCase): + + def test_single_triangle_2d(self): + v = np.array([[0.0,0.0],[1.0,0.0],[0.0,1.0]]) + f = np.array([[0,1,2]],dtype=int) + + Q = gpy.biharmonic_energy(v, f, bc='mixedfem_zero_neumann') + Q_gt = np.array([[8., -4., -4.], + [-4., 3., 1.], + [-4., 1., 3.]]) + self.assertTrue(np.isclose(Q.toarray(), Q_gt).all()) + + Q = gpy.biharmonic_energy(v, f, bc='hessian') + Q_gt = np.zeros((3,3)) + self.assertTrue(np.isclose(Q.toarray(), Q_gt).all()) + + Q = gpy.biharmonic_energy(v, f, bc='curved_hessian') + Q_gt = np.zeros((3,3)) + self.assertTrue(np.isclose(Q.toarray(), Q_gt).all()) + + + def test_kernel(self): + meshes = ['cube.obj', 'mountain.obj', 'armadillo.obj', 'bunny_oded.obj'] + for mesh in meshes: + V,F = gpy.read_mesh("test/unit_tests_data/" + mesh) + + c = np.ones(V.shape[0]) + + Qmzn = gpy.biharmonic_energy(V, F, bc='mixedfem_zero_neumann') + self.assertTrue(np.isclose(Qmzn*c, 0.).all()) + + Qh = gpy.biharmonic_energy(V, F, bc='hessian') + self.assertTrue(np.isclose(Qh*c, 0.).all()) + + Qch = gpy.biharmonic_energy(V, F, bc='curved_hessian') + self.assertTrue(np.isclose(Qch*c, 0.).all()) + + + def test_regression(self): + V,F = gpy.read_mesh("test/unit_tests_data/cube.obj") + + Qmzn = gpy.biharmonic_energy(V, F, bc='mixedfem_zero_neumann') + Qmzn_gt = np.array([[16. , -8. , -8. , + 2.6666666666666665, 2.6666666666666665, 0. , + -8. , 2.6666666666666665], + [-8. , 16. , 2.6666666666666665, + -8. , 0. , 2.666666666666667 , + 2.6666666666666665, -8. ], + [-8. , 2.6666666666666665, 16. , + -8. , -8. , 2.666666666666667 , + 2.6666666666666665, 0. ], + [ 2.6666666666666665, -8. , -8. , + 16. , 2.666666666666667 , -8. , + 0. , 2.666666666666667 ], + [ 2.6666666666666665, 0. , -8. , + 2.666666666666667 , 16. , -8. , + -8. , 2.666666666666667 ], + [ 0. , 2.666666666666667 , 2.666666666666667 , + -8. , -8. , 16. , + 2.666666666666667 , -8. ], + [-8. , 2.6666666666666665, 2.6666666666666665, + 0. , -8. , 2.666666666666667 , + 16. , -8. ], + [ 2.6666666666666665, -8. , 0. , + 2.666666666666667 , 2.666666666666667 , -8. , + -8. , 16. ]]) + self.assertTrue(np.isclose(Qmzn.toarray(), Qmzn_gt).all()) + + Qh = gpy.biharmonic_energy(V, F, bc='hessian') + Qh_gt = np.array([[11.999999999999998 , -3.999999999999999 , -3.9999999999999996, + 0. , 0. , 0. , + -4. , 0. ], + [-3.999999999999999 , 11.999999999999998 , 0. , + -3.999999999999999 , 0. , 0. , + 0. , -3.999999999999999 ], + [-3.9999999999999996, 0. , 11.999999999999998 , + -3.999999999999999 , -3.999999999999999 , 0. , + 0. , 0. ], + [ 0. , -3.999999999999999 , -3.999999999999999 , + 11.999999999999998 , 0. , -3.9999999999999987, + 0. , 0. ], + [ 0. , 0. , -3.999999999999999 , + 0. , 12. , -3.999999999999999 , + -4. , 0. ], + [ 0. , 0. , 0. , + -3.9999999999999987, -3.999999999999999 , 11.999999999999996 , + 0. , -3.9999999999999987], + [-4. , 0. , 0. , + 0. , -4. , 0. , + 12. , -3.999999999999999 ], + [ 0. , -3.999999999999999 , 0. , + 0. , 0. , -3.9999999999999987, + -3.999999999999999 , 11.999999999999996 ]]) + self.assertTrue(np.isclose(Qh.toarray(), Qh_gt).all()) + + Qch = gpy.biharmonic_energy(V, F, bc='curved_hessian') + Qch_gt = np.array([[ 22.37758040957278 , -10.403392041388944 , -10.403392041388944 , + 1.9999999999999982, 3.0471975511965974, 0.7382006122008503, + -10.403392041388942 , 3.0471975511965974], + [-10.403392041388942 , 21.853981633974485 , 3.5707963267948966, + -10.403392041388942 , 0.73820061220085 , 2.523598775598298 , + 2.523598775598298 , -10.403392041388946 ], + [-10.403392041388944 , 3.5707963267948966, 21.853981633974488 , + -10.403392041388944 , -10.403392041388946 , 2.5235987755982974, + 2.5235987755982974, 0.73820061220085 ], + [ 1.9999999999999982, -10.40339204138894 , -10.40339204138894 , + 22.37758040957278 , 3.047197551196597 , -10.403392041388942 , + 0.7382006122008503, 3.047197551196597 ], + [ 3.0471975511965974, 0.7382006122008502, -10.403392041388944 , + 3.0471975511965974, 22.37758040957278 , -10.403392041388944 , + -10.403392041388944 , 1.9999999999999978], + [ 0.7382006122008499, 2.5235987755982974, 2.5235987755982974, + -10.403392041388942 , -10.403392041388942 , 21.85398163397448 , + 3.5707963267948952, -10.403392041388942 ], + [-10.403392041388942 , 2.5235987755982974, 2.5235987755982974, + 0.7382006122008499, -10.40339204138894 , 3.5707963267948952, + 21.85398163397448 , -10.40339204138894 ], + [ 3.047197551196597 , -10.403392041388944 , 0.7382006122008502, + 3.047197551196597 , 1.9999999999999976, -10.40339204138894 , + -10.40339204138894 , 22.37758040957278 ]]) + self.assertTrue(np.isclose(Qch.toarray(), Qch_gt).all()) + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/test/test_biharmonic_energy_intrinsic.py b/test/test_biharmonic_energy_intrinsic.py new file mode 100644 index 00000000..f5c1656b --- /dev/null +++ b/test/test_biharmonic_energy_intrinsic.py @@ -0,0 +1,84 @@ +from .context import gpytoolbox as gpy +from .context import numpy as np +from .context import unittest + +class TestBiharmonicEnergyIntrinsic(unittest.TestCase): + + def test_uniform_triangle(self): + c = np.random.default_rng().random() + 0.1 + + l_sq = c * np.array([[1., 1., 1.]]) + f = np.array([[0,1,2]],dtype=int) + + Q = gpy.biharmonic_energy_intrinsic(l_sq, f, bc='mixedfem_zero_neumann') + Q_gt = np.sqrt(3)/c * np.array([[2., -1., -1.], + [-1., 2., -1.], + [-1., -1., 2.]]) + self.assertTrue(np.isclose(Q.toarray(), Q_gt).all()) + + Q = gpy.biharmonic_energy_intrinsic(l_sq, f, bc='curved_hessian') + Q_gt = np.zeros((3,3)) + self.assertTrue(np.isclose(Q.toarray(), Q_gt).all()) + + + def test_regression(self): + V,F = gpy.read_mesh("test/unit_tests_data/cube.obj") + l_sq = gpy.halfedge_lengths_squared(V,F) + + Qmzn = gpy.biharmonic_energy_intrinsic(l_sq, F, bc='mixedfem_zero_neumann') + Qmzn_gt = np.array([[16. , -8. , -8. , + 2.6666666666666665, 2.6666666666666665, 0. , + -8. , 2.6666666666666665], + [-8. , 16. , 2.6666666666666665, + -8. , 0. , 2.666666666666667 , + 2.6666666666666665, -8. ], + [-8. , 2.6666666666666665, 16. , + -8. , -8. , 2.666666666666667 , + 2.6666666666666665, 0. ], + [ 2.6666666666666665, -8. , -8. , + 16. , 2.666666666666667 , -8. , + 0. , 2.666666666666667 ], + [ 2.6666666666666665, 0. , -8. , + 2.666666666666667 , 16. , -8. , + -8. , 2.666666666666667 ], + [ 0. , 2.666666666666667 , 2.666666666666667 , + -8. , -8. , 16. , + 2.666666666666667 , -8. ], + [-8. , 2.6666666666666665, 2.6666666666666665, + 0. , -8. , 2.666666666666667 , + 16. , -8. ], + [ 2.6666666666666665, -8. , 0. , + 2.666666666666667 , 2.666666666666667 , -8. , + -8. , 16. ]]) + self.assertTrue(np.isclose(Qmzn.toarray(), Qmzn_gt).all()) + + Qch = gpy.biharmonic_energy_intrinsic(l_sq, F, bc='curved_hessian') + Qch_gt = np.array([[ 22.37758040957278 , -10.403392041388944 , -10.403392041388944 , + 1.9999999999999982, 3.0471975511965974, 0.7382006122008503, + -10.403392041388942 , 3.0471975511965974], + [-10.403392041388942 , 21.853981633974485 , 3.5707963267948966, + -10.403392041388942 , 0.73820061220085 , 2.523598775598298 , + 2.523598775598298 , -10.403392041388946 ], + [-10.403392041388944 , 3.5707963267948966, 21.853981633974488 , + -10.403392041388944 , -10.403392041388946 , 2.5235987755982974, + 2.5235987755982974, 0.73820061220085 ], + [ 1.9999999999999982, -10.40339204138894 , -10.40339204138894 , + 22.37758040957278 , 3.047197551196597 , -10.403392041388942 , + 0.7382006122008503, 3.047197551196597 ], + [ 3.0471975511965974, 0.7382006122008502, -10.403392041388944 , + 3.0471975511965974, 22.37758040957278 , -10.403392041388944 , + -10.403392041388944 , 1.9999999999999978], + [ 0.7382006122008499, 2.5235987755982974, 2.5235987755982974, + -10.403392041388942 , -10.403392041388942 , 21.85398163397448 , + 3.5707963267948952, -10.403392041388942 ], + [-10.403392041388942 , 2.5235987755982974, 2.5235987755982974, + 0.7382006122008499, -10.40339204138894 , 3.5707963267948952, + 21.85398163397448 , -10.40339204138894 ], + [ 3.047197551196597 , -10.403392041388944 , 0.7382006122008502, + 3.047197551196597 , 1.9999999999999976, -10.40339204138894 , + -10.40339204138894 , 22.37758040957278 ]]) + self.assertTrue(np.isclose(Qch.toarray(), Qch_gt).all()) + +if __name__ == '__main__': + unittest.main() + diff --git a/test/test_grad_intrinsic.py b/test/test_grad_intrinsic.py new file mode 100644 index 00000000..1771cb92 --- /dev/null +++ b/test/test_grad_intrinsic.py @@ -0,0 +1,96 @@ +from .context import gpytoolbox +from .context import numpy as np +from .context import scipy as sp +from .context import unittest +from .context import gpytoolbox as gpy + + +class TestGradIntrinsic(unittest.TestCase): + + def test_single_triangle_2d(self): + c = np.random.default_rng().random() + 0.1 + + l_sq = c * np.array([[1., 1., 1.]]) + f = np.array([[0,1,2]],dtype=int) + + G = gpy.grad_intrinsic(l_sq, f) + + G_gt = np.array([[0., -1./np.sqrt(c), 1./np.sqrt(c)], + [2./np.sqrt(3*c), -1./np.sqrt(3*c), -1./np.sqrt(3*c)]]) + + self.assertTrue(np.isclose(G.toarray(), G_gt).all()) + + + def test_compare_with_cotangent_laplacian(self): + meshes = ['cube.obj', 'mountain.obj', 'armadillo.obj'] + for mesh in meshes: + V,F = gpy.read_mesh("test/unit_tests_data/" + mesh) + m = F.shape[0] + l_sq = gpy.halfedge_lengths_squared(V,F) + G = gpy.grad_intrinsic(l_sq, F) + a = gpy.doublearea_intrinsic(l_sq, F) / 2. + A = sp.sparse.spdiags([np.tile(a, 2)], 0, m=2*m, n=2*m, format='csr') + L = G.transpose() * A * G + + L_gt = gpy.cotangent_laplacian(V,F) + + self.assertTrue(np.isclose(L[L_gt.nonzero()], L_gt[L_gt.nonzero()]).all()) + + + def test_example_compared_with_known_gt(self): + l_sq = np.array([[2., 1., 1.], + [1., 1., 2.], + [2., 1., 1.], + [1., 1., 2.], + [2., 1., 1.], + [1., 1., 2.], + [2., 1., 1.], + [1., 1., 2.], + [2., 1., 1.], + [1., 1., 2.], + [2., 1., 1.], + [1., 1., 2.]]) + F = np.array([[1, 2, 3], + [3, 2, 4], + [3, 4, 5], + [5, 4, 6], + [5, 6, 7], + [7, 6, 8], + [7, 8, 1], + [1, 8, 2], + [2, 8, 4], + [4, 8, 6], + [7, 1, 5], + [5, 1, 3]]) - 1 + G_gt = np.array([[0, -0.70711, 0.70711, 0, 0, 0, 0, 0], + [0, -1, 0, 1, 0, 0, 0, 0], + [0, 0, 0, -0.70711, 0.70711, 0, 0, 0], + [0, 0, 0, -1, 0, 1, 0, 0], + [0, 0, 0, 0, 0, -0.70711, 0.70711, 0], + [0, 0, 0, 0, 0, -1, 0, 1], + [0.70711, 0, 0, 0, 0, 0, 0, -0.70711], + [0, 1, 0, 0, 0, 0, 0, -1], + [0, 0, 0, 0.70711, 0, 0, 0, -0.70711], + [0, 0, 0, 0, 0, 1, 0, -1], + [-0.70711, 0, 0, 0, 0.70711, 0, 0, 0], + [-1, 0, 1, 0, 0, 0, 0, 0], + [1.4142, -0.70711, -0.70711, 0, 0, 0, 0, 0], + [0, 2.2204e-16, 1, -1, 0, 0, 0, 0], + [0, 0, 1.4142, -0.70711, -0.70711, 0, 0, 0], + [0, 0, 0, 2.2204e-16, 1, -1, 0, 0], + [0, 0, 0, 0, 1.4142, -0.70711, -0.70711, 0], + [0, 0, 0, 0, 0, 2.2204e-16, 1, -1], + [-0.70711, 0, 0, 0, 0, 0, 1.4142, -0.70711], + [1, -1, 0, 0, 0, 0, 0, 2.2204e-16], + [0, 1.4142, 0, -0.70711, 0, 0, 0, -0.70711], + [0, 0, 0, 1, 0, -1, 0, 2.2204e-16], + [-0.70711, 0, 0, 0, -0.70711, 0, 1.4142, 0], + [2.2204e-16, 0, -1, 0, 1, 0, 0, 0]]) + + G = gpy.grad_intrinsic(l_sq, F) + + self.assertTrue(np.isclose(G.toarray(), G_gt).all()) + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file