21332822416

Committed 25 Jan 2026 12:46PM UTC coverage: 91.017% (-0.01%) from 91.028%

Build # 21332822416

Build Type

push

github

Committed by

dsieger

Commit Message

Disable failing check

Run Details

5238 of 5755 relevant lines covered (91.02%)

643396.87 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

95.79

/src/pmp/algorithms/laplace.cpp

// Copyright 2011-2023 the Polygon Mesh Processing Library developers.
// Copyright 2020 Astrid Bunge, Philipp Herholz, Misha Kazhdan, Mario Botsch.
// SPDX-License-Identifier: MIT

#include "pmp/algorithms/laplace.h"

// #define TFEM

namespace pmp {

namespace {

// compute triangle area (in double)
double triarea(const Eigen::Vector3d& p0, const Eigen::Vector3d& p1,
               const Eigen::Vector3d& p2)
{
    double double_area = (p1 - p0).cross(p2 - p0).norm();

#ifdef TFEM
    double l0 = (p1 - p2).squaredNorm();
    double l1 = (p0 - p2).squaredNorm();
    double l2 = (p0 - p1).squaredNorm();
    double h = pow((sqrt(l0) + sqrt(l1) + sqrt(l2))/3.0, 2);

    const double hmin = 1e-20; // absolute failsafe bound
    const double C = 1e-3; // dynamic TFEM constant

    double_area = fmax(double_area, C*fmax(h,hmin));
#endif

    return 0.5 * double_area;
}

// compute virtual vertex per polygon, represented by affine weights,
// such that the resulting triangle fan minimizes the sum of squared triangle areas
void compute_virtual_vertex(const DenseMatrix& poly, Eigen::VectorXd& weights)
{
    const int n = poly.rows();

    // setup array of positions and edges
    std::vector<dvec3> x(n), d(n);
    for (int i = 0; i < n; ++i)
        x[i] = poly.row(i);
    for (int i = 0; i < n; ++i)
        d[i] = x[(i + 1) % n] - x[i];

    // setup matrix A and rhs b
    // see Equation (38) of "Polygon Laplacian made simple", Eurographics 2020
    DenseMatrix A(n + 1, n);
    Eigen::VectorXd b(n + 1);
    for (int j = 0; j < n; ++j)
    {
        for (int i = j; i < n; ++i)
        {
            double Aij(0.0), bi(0.0);
            for (int k = 0; k < n; ++k)
            {
                Aij += dot(cross(x[j], d[k]), cross(x[i], d[k]));
                bi += dot(cross(x[i], d[k]), cross(x[k], d[k]));
            }
            A(i, j) = A(j, i) = Aij;
            b(i) = bi;
        }
    }
    for (int j = 0; j < n; ++j)
    {
        A(n, j) = 1.0;
    }
    b(n) = 1.0;

    weights = A.completeOrthogonalDecomposition().solve(b).topRows(n);
}

void triangle_mass_matrix(const Eigen::Vector3d& p0, const Eigen::Vector3d& p1,
                          const Eigen::Vector3d& p2, DiagonalMatrix& Mtri)
{
    // three vertex positions
    const std::array<dvec3, 3> p = {p0, p1, p2};

    // edge vectors
    std::array<dvec3, 3> e;
    for (int i = 0; i < 3; ++i)
        e[i] = p[(i + 1) % 3] - p[i];

    // compute and check (twice the) triangle area
#ifndef TFEM
    const auto tri_area = norm(cross(e[0], e[1]));
#else
    const auto tri_area = 2.0 * triarea(p0, p1, p2);
#endif
    if (tri_area <= std::numeric_limits<double>::min())
    {
        Mtri.setZero(3);
        return;
    }

    // dot products for each corner (of its two emanating edge vectors)
    std::array<double, 3> d;
    for (int i = 0; i < 3; ++i)
        d[i] = -dot(e[i], e[(i + 2) % 3]);

    // cotangents for each corner: cot = cos/sin = dot(A,B)/norm(cross(A,B))
    std::array<double, 3> cot;
    for (int i = 0; i < 3; ++i)
        cot[i] = d[i] / tri_area;

    // compute area for each corner
    Eigen::Vector3d area;
    for (int i = 0; i < 3; ++i)
    {
        // angle at corner is obtuse
        if (d[i] < 0.0)
        {
            area[i] = 0.25 * tri_area;
        }
        // angle at some other corner is obtuse
        else if (d[(i + 1) % 3] < 0.0 || d[(i + 2) % 3] < 0.0)
        {
            area[i] = 0.125 * tri_area;
        }
        // no obtuse angles
        else
        {
            area[i] = 0.125 * (sqrnorm(e[i]) * cot[(i + 2) % 3] +
                               sqrnorm(e[(i + 2) % 3]) * cot[(i + 1) % 3]);
        }
    }
    Mtri = area.asDiagonal();
}

void polygon_mass_matrix(const DenseMatrix& polygon, DiagonalMatrix& Mpoly)
{
    const int n = (int)polygon.rows();

    // shortcut for triangles
    if (n == 3)
    {
        triangle_mass_matrix(polygon.row(0), polygon.row(1), polygon.row(2),
                             Mpoly);
        return;
    }

    // compute position of virtual vertex
    Eigen::VectorXd vweights;
    compute_virtual_vertex(polygon, vweights);
    const Eigen::Vector3d vvertex = polygon.transpose() * vweights;

    // laplace matrix of refined triangle fan
    DenseMatrix Mfan = DenseMatrix::Zero(n + 1, n + 1);
    DiagonalMatrix Mtri;
    for (int i = 0; i < n; ++i)
    {
        const int j = (i + 1) % n;

        // build laplace matrix of one triangle
        triangle_mass_matrix(polygon.row(i), polygon.row(j), vvertex, Mtri);

        // assemble to laplace matrix for refined triangle fan
        // (we are dealing with diagonal matrices)
        Mfan.diagonal()[i] += Mtri.diagonal()[0];
        Mfan.diagonal()[j] += Mtri.diagonal()[1];
        Mfan.diagonal()[n] += Mtri.diagonal()[2];
    }

    // build prolongation matrix
    DenseMatrix P(n + 1, n);
    P.setIdentity();
    P.row(n) = vweights;

    // build polygon laplace matrix by sandwiching
    DenseMatrix PMP = P.transpose() * Mfan * P;
    Mpoly = PMP.rowwise().sum().asDiagonal();
}

void triangle_laplace_matrix(const Eigen::Vector3d& p0,
                             const Eigen::Vector3d& p1,
                             const Eigen::Vector3d& p2, DenseMatrix& Ltri)
{
#ifndef TFEM
    std::array<double, 3> l, l2, cot;

    // squared edge lengths
    l2[0] = (p1 - p2).squaredNorm();
    l2[1] = (p0 - p2).squaredNorm();
    l2[2] = (p0 - p1).squaredNorm();

    // edge lengths
    l[0] = sqrt(l2[0]);
    l[1] = sqrt(l2[1]);
    l[2] = sqrt(l2[2]);

    // triangle area
    const double arg = (l[0] + (l[1] + l[2])) * (l[2] - (l[0] - l[1])) *
                       (l[2] + (l[0] - l[1])) * (l[0] + (l[1] - l[2]));
    const double area = 0.5 * sqrt(arg);

    if (area > std::numeric_limits<double>::min())
    {
        // cotangents of angles
        cot[0] = 0.25 * (l2[1] + l2[2] - l2[0]) / area;
        cot[1] = 0.25 * (l2[2] + l2[0] - l2[1]) / area;
        cot[2] = 0.25 * (l2[0] + l2[1] - l2[2]) / area;

        Ltri(0, 0) = cot[1] + cot[2];
        Ltri(1, 1) = cot[0] + cot[2];
        Ltri(2, 2) = cot[0] + cot[1];
        Ltri(1, 0) = Ltri(0, 1) = -cot[2];
        Ltri(2, 0) = Ltri(0, 2) = -cot[1];
        Ltri(2, 1) = Ltri(1, 2) = -cot[0];
    }
#else
    // double triangle area
    double area = 2.0*triarea(p0, p1, p2);

    if (area > std::numeric_limits<double>::min())
    {
        // squared edge lengths
        std::array<double, 3> l2{};
        l2[0] = (p1 - p2).squaredNorm();
        l2[1] = (p0 - p2).squaredNorm();
        l2[2] = (p0 - p1).squaredNorm();

        // cotangent values
        std::array<double, 3> cot{};
        const double inv_area = 0.25 / area;
        cot[0] = inv_area * (l2[1] + l2[2] - l2[0]);
        cot[1] = inv_area * (l2[2] + l2[0] - l2[1]);
        cot[2] = inv_area * (l2[0] + l2[1] - l2[2]);

        // build 3x3 Laplace matrix
        Ltri(0, 0) = cot[1] + cot[2];
        Ltri(1, 1) = cot[0] + cot[2];
        Ltri(2, 2) = cot[0] + cot[1];
        Ltri(1, 0) = Ltri(0, 1) = -cot[2];
        Ltri(2, 0) = Ltri(0, 2) = -cot[1];
        Ltri(2, 1) = Ltri(1, 2) = -cot[0];
    }
    else
    {
        std::cerr << "Should not happen\n";
        Ltri.setZero();
    }
#endif
}

void polygon_laplace_matrix(const DenseMatrix& polygon, DenseMatrix& Lpoly)
{
    const int n = (int)polygon.rows();
    Lpoly = DenseMatrix::Zero(n, n);

    // shortcut for triangles
    if (n == 3)
    {
        triangle_laplace_matrix(polygon.row(0), polygon.row(1), polygon.row(2),
                                Lpoly);
        return;
    }

    // compute position of virtual vertex
    Eigen::VectorXd vweights;
    compute_virtual_vertex(polygon, vweights);
    const Eigen::Vector3d vvertex = polygon.transpose() * vweights;

    // laplace matrix of refined triangle fan
    DenseMatrix Lfan = DenseMatrix::Zero(n + 1, n + 1);
    DenseMatrix Ltri(3, 3);
    for (int i = 0; i < n; ++i)
    {
        const int j = (i + 1) % n;

        // build laplace matrix of one triangle
        triangle_laplace_matrix(polygon.row(i), polygon.row(j), vvertex, Ltri);

        // assemble to laplace matrix for refined triangle fan
        Lfan(i, i) += Ltri(0, 0);
        Lfan(i, j) += Ltri(0, 1);
        Lfan(i, n) += Ltri(0, 2);
        Lfan(j, i) += Ltri(1, 0);
        Lfan(j, j) += Ltri(1, 1);
        Lfan(j, n) += Ltri(1, 2);
        Lfan(n, i) += Ltri(2, 0);
        Lfan(n, j) += Ltri(2, 1);
        Lfan(n, n) += Ltri(2, 2);
    }

    // build prolongation matrix
    DenseMatrix P(n + 1, n);
    P.setIdentity();
    P.row(n) = vweights;

    // build polygon laplace matrix by sandwiching
    Lpoly = P.transpose() * Lfan * P;
}

void triangle_gradient_matrix(const Eigen::Vector3d& p0,
                              const Eigen::Vector3d& p1,
                              const Eigen::Vector3d& p2, DenseMatrix& G)
{
    G.resize(3, 3);
    Eigen::Vector3d n = (p1 - p0).cross(p2 - p0);
#ifndef TFEM
    double double_area = n.norm();
#else
    double double_area = 2.0 * triarea(p0, p1, p2);
#endif
    if (double_area > std::numeric_limits<double>::min())
    {
        // n /= double_area;
        n.normalize();
        G.col(0) = n.cross(p2 - p1) / (double_area);
        G.col(1) = n.cross(p0 - p2) / (double_area);
        G.col(2) = n.cross(p1 - p0) / (double_area);
    }
    else
    {
        G.setZero();
    }
}

void polygon_gradient_matrix(const DenseMatrix& polygon, DenseMatrix& Gpoly)
{
    const int n = (int)polygon.rows();

    // compute position of virtual vertex
    Eigen::VectorXd vweights;
    compute_virtual_vertex(polygon, vweights);
    const Eigen::Vector3d vvertex = polygon.transpose() * vweights;

    DenseMatrix Gfan = DenseMatrix::Zero(3 * n, n + 1);
    DenseMatrix Gtri(3, 3);
    int row = 0;
    for (int i = 0; i < n; ++i)
    {
        const int j = (i + 1) % n;

        // build laplace matrix of one triangle
        triangle_gradient_matrix(polygon.row(i), polygon.row(j), vvertex, Gtri);

        // assemble to matrix for triangle fan
        for (int k = 0; k < 3; ++k)
        {
            Gfan(row + k, i) += Gtri(k, 0);
            Gfan(row + k, j) += Gtri(k, 1);
            Gfan(row + k, n) += Gtri(k, 2);
        }

        row += 3;
    }

    // build prolongation matrix
    DenseMatrix P(n + 1, n);
    P.setIdentity();
    P.row(n) = vweights;

    // build polygon gradient matrix by sandwiching (from left only)
    Gpoly = Gfan * P;
}

void divmass_matrix(const SurfaceMesh& mesh, DiagonalMatrix& M)
{
    // how many virtual triangles will we have after refinement?
    unsigned int nt = 0;
    for (auto f : mesh.faces())
    {
        nt += mesh.valence(f);
    }

    // initialize global matrix
    M.resize(3 * nt);
    auto& diag = M.diagonal();

    std::vector<Vertex> vertices; // polygon vertices
    DenseMatrix polygon;          // positions of polygon vertices

    unsigned int idx = 0;

    for (const auto f : mesh.faces())
    {
        // collect polygon vertices
        vertices.clear();
        for (const auto v : mesh.vertices(f))
        {
            vertices.push_back(v);
        }
        const int n = vertices.size();

        // collect their positions
        polygon.resize(n, 3);
        for (int i = 0; i < n; ++i)
        {
            polygon.row(i) = (Eigen::Vector3d)mesh.position(vertices[i]);
        }

        // compute position of virtual vertex
        Eigen::VectorXd vweights;
        compute_virtual_vertex(polygon, vweights);
        const Eigen::Vector3d vvertex = polygon.transpose() * vweights;

        for (int i = 0; i < n; ++i)
        {
            const double area =
                triarea(polygon.row(i), polygon.row((i + 1) % n), vvertex);

            diag[idx++] = area;
            diag[idx++] = area;
            diag[idx++] = area;
        }
    }

    assert(idx == 3 * nt);
}

} // anonymous namespace

void uniform_mass_matrix(const SurfaceMesh& mesh, DiagonalMatrix& M)
{
    const unsigned int n = mesh.n_vertices();
    Eigen::VectorXd diag(n);
    for (auto v : mesh.vertices())
        diag[v.idx()] = mesh.valence(v);
    M = diag.asDiagonal();
}

void mass_matrix(const SurfaceMesh& mesh, DiagonalMatrix& M)
{
    const int nv = mesh.n_vertices();
    std::vector<Vertex> vertices; // polygon vertices
    DenseMatrix polygon;          // positions of polygon vertices
    DiagonalMatrix Mpoly;         // local mass matrix

    M.setZero(nv);

    for (const auto f : mesh.faces())
    {
        // collect polygon vertices
        vertices.clear();
        for (const auto v : mesh.vertices(f))
        {
            vertices.push_back(v);
        }
        const int n = vertices.size();

        // collect their positions
        polygon.resize(n, 3);
        for (int i = 0; i < n; ++i)
        {
            polygon.row(i) = (Eigen::Vector3d)mesh.position(vertices[i]);
        }

        // setup local mass matrix
        polygon_mass_matrix(polygon, Mpoly);

        // assemble to global mass matrix
        for (int k = 0; k < n; ++k)
        {
            M.diagonal()[vertices[k].idx()] += Mpoly.diagonal()[k];
        }
    }
}

void uniform_laplace_matrix(const SurfaceMesh& mesh, SparseMatrix& L)
{
    const unsigned int n = mesh.n_vertices();

    std::vector<Triplet> triplets;
    triplets.reserve(8 * n); // conservative estimate for triangle meshes

    for (auto vi : mesh.vertices())
    {
        Scalar sum_weights = 0.0;
        for (auto vj : mesh.vertices(vi))
        {
            sum_weights += 1.0;
            triplets.emplace_back(vi.idx(), vj.idx(), 1.0);
        }
        triplets.emplace_back(vi.idx(), vi.idx(), -sum_weights);
    }

    L.resize(n, n);
    L.setFromTriplets(triplets.begin(), triplets.end());
}

void laplace_matrix(const SurfaceMesh& mesh, SparseMatrix& L, bool clamp)
{
    const int nv = mesh.n_vertices();
    const int nf = mesh.n_faces();
    std::vector<Vertex> vertices; // polygon vertices
    DenseMatrix polygon;          // positions of polygon vertices
    DenseMatrix Lpoly;            // local laplace matrix

    std::vector<Triplet> triplets;
    triplets.reserve(9 * nf); // estimate for triangle meshes

    for (const auto f : mesh.faces())
    {
        // collect polygon vertices
        vertices.clear();
        for (const auto v : mesh.vertices(f))
        {
            vertices.push_back(v);
        }
        const int n = vertices.size();

        // collect their positions
        polygon.resize(n, 3);
        for (int i = 0; i < n; ++i)
        {
            polygon.row(i) = (Eigen::Vector3d)mesh.position(vertices[i]);
        }

        // setup local laplace matrix
        polygon_laplace_matrix(polygon, Lpoly);

        // assemble to global laplace matrix
        for (int j = 0; j < n; ++j)
        {
            for (int k = 0; k < n; ++k)
            {
                triplets.emplace_back(vertices[k].idx(), vertices[j].idx(),
                                      -Lpoly(k, j));
            }
        }
    }

    // build sparse matrix from triplets
    L.resize(nv, nv);
    L.setFromTriplets(triplets.begin(), triplets.end());

    // clamp negative off-diagonal entries to zero
    if (clamp)
    {
        for (unsigned int k = 0; k < L.outerSize(); k++)
        {
            double diag_offset(0.0);

            for (SparseMatrix::InnerIterator iter(L, k); iter; ++iter)
            {
                if (iter.row() != iter.col() && iter.value() < 0.0)
                {
                    diag_offset += -iter.value();
                    iter.valueRef() = 0.0;
                }
            }
            for (SparseMatrix::InnerIterator iter(L, k); iter; ++iter)
            {
                if (iter.row() == iter.col() && iter.value() < 0.0)
                    iter.valueRef() -= diag_offset;
            }
        }
    }
}

void gradient_matrix(const SurfaceMesh& mesh, SparseMatrix& G)
{
    const int nv = mesh.n_vertices();

    // how many virtual triangles will we have after refinement?
    // triangles are not refined, other polygons are.
    unsigned int nt = 0;
    for (auto f : mesh.faces())
    {
        nt += mesh.valence(f);
    }

    std::vector<Vertex> vertices; // polygon vertices
    DenseMatrix polygon;          // positions of polygon vertices
    DenseMatrix Gpoly;            // local gradient matrix

    std::vector<Triplet> triplets;
    triplets.reserve(9 * nt);

    unsigned int n_rows = 0;

    for (const auto f : mesh.faces())
    {
        // collect polygon vertices
        vertices.clear();
        for (const auto v : mesh.vertices(f))
        {
            vertices.push_back(v);
        }
        const int n = vertices.size();

        // collect their positions
        polygon.resize(n, 3);
        for (int i = 0; i < n; ++i)
        {
            polygon.row(i) = (Eigen::Vector3d)mesh.position(vertices[i]);
        }

        // setup local element matrix
        polygon_gradient_matrix(polygon, Gpoly);

        // assemble to global matrix
        for (int j = 0; j < Gpoly.cols(); ++j)
        {
            for (int i = 0; i < Gpoly.rows(); ++i)
            {
                triplets.emplace_back(n_rows + i, vertices[j].idx(),
                                      Gpoly(i, j));
            }
        }

        n_rows += Gpoly.rows();
    }
    assert(n_rows == 3 * nt);

    // build sparse matrix from triplets
    G.resize(n_rows, nv);
    G.setFromTriplets(triplets.begin(), triplets.end());
}

void divergence_matrix(const SurfaceMesh& mesh, SparseMatrix& D)
{
    SparseMatrix G;
    gradient_matrix(mesh, G);
    DiagonalMatrix M;
    divmass_matrix(mesh, M);
    D = -G.transpose() * M;
}

} // namespace pmp

1	// Copyright 2011-2023 the Polygon Mesh Processing Library developers.
2	// Copyright 2020 Astrid Bunge, Philipp Herholz, Misha Kazhdan, Mario Botsch.
3	// SPDX-License-Identifier: MIT
4
5	#include "pmp/algorithms/laplace.h"
6
7	// #define TFEM
8
9	namespace pmp {
10
11	namespace {
12
13	// compute triangle area (in double)
14	double triarea(const Eigen::Vector3d& p0, const Eigen::Vector3d& p1,	10,086✔
15	const Eigen::Vector3d& p2)
16	{
17	double double_area = (p1 - p0).cross(p2 - p0).norm();	10,086✔
18
19	#ifdef TFEM
20	double l0 = (p1 - p2).squaredNorm();
21	double l1 = (p0 - p2).squaredNorm();
22	double l2 = (p0 - p1).squaredNorm();
23	double h = pow((sqrt(l0) + sqrt(l1) + sqrt(l2))/3.0, 2);
24
25	const double hmin = 1e-20; // absolute failsafe bound
26	const double C = 1e-3; // dynamic TFEM constant
27
28	double_area = fmax(double_area, C*fmax(h,hmin));
29	#endif
30
31	return 0.5 * double_area;	10,086✔
32	}
33
34	// compute virtual vertex per polygon, represented by affine weights,
35	// such that the resulting triangle fan minimizes the sum of squared triangle areas
36	void compute_virtual_vertex(const DenseMatrix& poly, Eigen::VectorXd& weights)	11,113✔
37	{
38	const int n = poly.rows();	11,113✔
39
40	// setup array of positions and edges
41	std::vector<dvec3> x(n), d(n);	33,339✔
42	for (int i = 0; i < n; ++i)	47,763✔
43	x[i] = poly.row(i);	36,650✔
44	for (int i = 0; i < n; ++i)	47,763✔
45	d[i] = x[(i + 1) % n] - x[i];	36,650✔
46
47	// setup matrix A and rhs b
48	// see Equation (38) of "Polygon Laplacian made simple", Eurographics 2020
49	DenseMatrix A(n + 1, n);	11,113✔
50	Eigen::VectorXd b(n + 1);	11,113✔
51	for (int j = 0; j < n; ++j)	47,763✔
52	{
53	for (int i = j; i < n; ++i)	116,572✔
54	{
55	double Aij(0.0), bi(0.0);	79,922✔
56	for (int k = 0; k < n; ++k)	352,798✔
57	{
58	Aij += dot(cross(x[j], d[k]), cross(x[i], d[k]));	272,876✔
59	bi += dot(cross(x[i], d[k]), cross(x[k], d[k]));	272,876✔
60	}
61	A(i, j) = A(j, i) = Aij;	79,922✔
62	b(i) = bi;	79,922✔
63	}
64	}
65	for (int j = 0; j < n; ++j)	47,763✔
66	{
67	A(n, j) = 1.0;	36,650✔
68	}
69	b(n) = 1.0;	11,113✔
70
71	weights = A.completeOrthogonalDecomposition().solve(b).topRows(n);	11,113✔
72	}	11,113✔
73
74	void triangle_mass_matrix(const Eigen::Vector3d& p0, const Eigen::Vector3d& p1,	110,440✔
75	const Eigen::Vector3d& p2, DiagonalMatrix& Mtri)
76	{
77	// three vertex positions
78	const std::array<dvec3, 3> p = {p0, p1, p2};	110,440✔
79
80	// edge vectors
81	std::array<dvec3, 3> e;
82	for (int i = 0; i < 3; ++i)	441,760✔
83	e[i] = p[(i + 1) % 3] - p[i];	331,320✔
84
85	// compute and check (twice the) triangle area
86	#ifndef TFEM
87	const auto tri_area = norm(cross(e[0], e[1]));	110,440✔
88	#else
89	const auto tri_area = 2.0 * triarea(p0, p1, p2);
90	#endif
91	if (tri_area <= std::numeric_limits<double>::min())	110,440✔
92	{
93	Mtri.setZero(3);	×
94	return;	×
95	}
96
97	// dot products for each corner (of its two emanating edge vectors)
98	std::array<double, 3> d;
99	for (int i = 0; i < 3; ++i)	441,760✔
100	d[i] = -dot(e[i], e[(i + 2) % 3]);	331,320✔
101
102	// cotangents for each corner: cot = cos/sin = dot(A,B)/norm(cross(A,B))
103	std::array<double, 3> cot;
104	for (int i = 0; i < 3; ++i)	441,760✔
105	cot[i] = d[i] / tri_area;	331,320✔
106
107	// compute area for each corner
108	Eigen::Vector3d area;	110,440✔
109	for (int i = 0; i < 3; ++i)	441,760✔
110	{
111	// angle at corner is obtuse
112	if (d[i] < 0.0)	331,320✔
113	{
114	area[i] = 0.25 * tri_area;	1,982✔
115	}
116	// angle at some other corner is obtuse
117	else if (d[(i + 1) % 3] < 0.0 \|\| d[(i + 2) % 3] < 0.0)	329,338✔
118	{
119	area[i] = 0.125 * tri_area;	3,964✔
120	}
121	// no obtuse angles
122	else
123	{
124	area[i] = 0.125 * (sqrnorm(e[i]) * cot[(i + 2) % 3] +	650,748✔
125	sqrnorm(e[(i + 2) % 3]) * cot[(i + 1) % 3]);	325,374✔
126	}
127	}
128	Mtri = area.asDiagonal();	110,440✔
129	}
130
131	void polygon_mass_matrix(const DenseMatrix& polygon, DiagonalMatrix& Mpoly)	107,551✔
132	{
133	const int n = (int)polygon.rows();	107,551✔
134
135	// shortcut for triangles
136	if (n == 3)	107,551✔
137	{
138	triangle_mass_matrix(polygon.row(0), polygon.row(1), polygon.row(2),	106,588✔
139	Mpoly);
140	return;	106,588✔
141	}
142
143	// compute position of virtual vertex
144	Eigen::VectorXd vweights;	963✔
145	compute_virtual_vertex(polygon, vweights);	963✔
146	const Eigen::Vector3d vvertex = polygon.transpose() * vweights;	963✔
147
148	// laplace matrix of refined triangle fan
149	DenseMatrix Mfan = DenseMatrix::Zero(n + 1, n + 1);	963✔
150	DiagonalMatrix Mtri;	963✔
151	for (int i = 0; i < n; ++i)	4,815✔
152	{
153	const int j = (i + 1) % n;	3,852✔
154
155	// build laplace matrix of one triangle
156	triangle_mass_matrix(polygon.row(i), polygon.row(j), vvertex, Mtri);	3,852✔
157
158	// assemble to laplace matrix for refined triangle fan
159	// (we are dealing with diagonal matrices)
160	Mfan.diagonal()[i] += Mtri.diagonal()[0];	3,852✔
161	Mfan.diagonal()[j] += Mtri.diagonal()[1];	3,852✔
162	Mfan.diagonal()[n] += Mtri.diagonal()[2];	3,852✔
163	}
164
165	// build prolongation matrix
166	DenseMatrix P(n + 1, n);	963✔
167	P.setIdentity();	963✔
168	P.row(n) = vweights;	963✔
169
170	// build polygon laplace matrix by sandwiching
171	DenseMatrix PMP = P.transpose() * Mfan * P;	963✔
172	Mpoly = PMP.rowwise().sum().asDiagonal();	963✔
173	}	963✔
174
175	void triangle_laplace_matrix(const Eigen::Vector3d& p0,	209,948✔
176	const Eigen::Vector3d& p1,
177	const Eigen::Vector3d& p2, DenseMatrix& Ltri)
178	{
179	#ifndef TFEM
180	std::array<double, 3> l, l2, cot;
181
182	// squared edge lengths
183	l2[0] = (p1 - p2).squaredNorm();	209,948✔
184	l2[1] = (p0 - p2).squaredNorm();	209,948✔
185	l2[2] = (p0 - p1).squaredNorm();	209,948✔
186
187	// edge lengths
188	l[0] = sqrt(l2[0]);	209,948✔
189	l[1] = sqrt(l2[1]);	209,948✔
190	l[2] = sqrt(l2[2]);	209,948✔
191
192	// triangle area
193	const double arg = (l[0] + (l[1] + l[2])) * (l[2] - (l[0] - l[1])) *	209,948✔
194	(l[2] + (l[0] - l[1])) * (l[0] + (l[1] - l[2]));	209,948✔
195	const double area = 0.5 * sqrt(arg);	209,948✔
196
197	if (area > std::numeric_limits<double>::min())	209,948✔
198	{
199	// cotangents of angles
200	cot[0] = 0.25 * (l2[1] + l2[2] - l2[0]) / area;	209,946✔
201	cot[1] = 0.25 * (l2[2] + l2[0] - l2[1]) / area;	209,946✔
202	cot[2] = 0.25 * (l2[0] + l2[1] - l2[2]) / area;	209,946✔
203
204	Ltri(0, 0) = cot[1] + cot[2];	209,946✔
205	Ltri(1, 1) = cot[0] + cot[2];	209,946✔
206	Ltri(2, 2) = cot[0] + cot[1];	209,946✔
207	Ltri(1, 0) = Ltri(0, 1) = -cot[2];	209,946✔
208	Ltri(2, 0) = Ltri(0, 2) = -cot[1];	209,946✔
209	Ltri(2, 1) = Ltri(1, 2) = -cot[0];	209,946✔
210	}
211	#else
212	// double triangle area
213	double area = 2.0*triarea(p0, p1, p2);
214
215	if (area > std::numeric_limits<double>::min())
216	{
217	// squared edge lengths
218	std::array<double, 3> l2{};
219	l2[0] = (p1 - p2).squaredNorm();
220	l2[1] = (p0 - p2).squaredNorm();
221	l2[2] = (p0 - p1).squaredNorm();
222
223	// cotangent values
224	std::array<double, 3> cot{};
225	const double inv_area = 0.25 / area;
226	cot[0] = inv_area * (l2[1] + l2[2] - l2[0]);
227	cot[1] = inv_area * (l2[2] + l2[0] - l2[1]);
228	cot[2] = inv_area * (l2[0] + l2[1] - l2[2]);
229
230	// build 3x3 Laplace matrix
231	Ltri(0, 0) = cot[1] + cot[2];
232	Ltri(1, 1) = cot[0] + cot[2];
233	Ltri(2, 2) = cot[0] + cot[1];
234	Ltri(1, 0) = Ltri(0, 1) = -cot[2];
235	Ltri(2, 0) = Ltri(0, 2) = -cot[1];
236	Ltri(2, 1) = Ltri(1, 2) = -cot[0];
237	}
238	else
239	{
240	std::cerr << "Should not happen\n";
241	Ltri.setZero();
242	}
243	#endif
244	}	209,948✔
245
246	void polygon_laplace_matrix(const DenseMatrix& polygon, DenseMatrix& Lpoly)	208,163✔
247	{
248	const int n = (int)polygon.rows();	208,163✔
249	Lpoly = DenseMatrix::Zero(n, n);	208,163✔
250
251	// shortcut for triangles
252	if (n == 3)	208,163✔
253	{
254	triangle_laplace_matrix(polygon.row(0), polygon.row(1), polygon.row(2),	207,568✔
255	Lpoly);
256	return;	207,568✔
257	}
258
259	// compute position of virtual vertex
260	Eigen::VectorXd vweights;	595✔
261	compute_virtual_vertex(polygon, vweights);	595✔
262	const Eigen::Vector3d vvertex = polygon.transpose() * vweights;	595✔
263
264	// laplace matrix of refined triangle fan
265	DenseMatrix Lfan = DenseMatrix::Zero(n + 1, n + 1);	595✔
266	DenseMatrix Ltri(3, 3);	595✔
267	for (int i = 0; i < n; ++i)	2,975✔
268	{
269	const int j = (i + 1) % n;	2,380✔
270
271	// build laplace matrix of one triangle
272	triangle_laplace_matrix(polygon.row(i), polygon.row(j), vvertex, Ltri);	2,380✔
273
274	// assemble to laplace matrix for refined triangle fan
275	Lfan(i, i) += Ltri(0, 0);	2,380✔
276	Lfan(i, j) += Ltri(0, 1);	2,380✔
277	Lfan(i, n) += Ltri(0, 2);	2,380✔
278	Lfan(j, i) += Ltri(1, 0);	2,380✔
279	Lfan(j, j) += Ltri(1, 1);	2,380✔
280	Lfan(j, n) += Ltri(1, 2);	2,380✔
281	Lfan(n, i) += Ltri(2, 0);	2,380✔
282	Lfan(n, j) += Ltri(2, 1);	2,380✔
283	Lfan(n, n) += Ltri(2, 2);	2,380✔
284	}
285
286	// build prolongation matrix
287	DenseMatrix P(n + 1, n);	595✔
288	P.setIdentity();	595✔
289	P.row(n) = vweights;	595✔
290
291	// build polygon laplace matrix by sandwiching
292	Lpoly = P.transpose() * Lfan * P;	595✔
293	}	595✔
294
295	void triangle_gradient_matrix(const Eigen::Vector3d& p0,	20,332✔
296	const Eigen::Vector3d& p1,
297	const Eigen::Vector3d& p2, DenseMatrix& G)
298	{
299	G.resize(3, 3);	20,332✔
300	Eigen::Vector3d n = (p1 - p0).cross(p2 - p0);	20,332✔
301	#ifndef TFEM
302	double double_area = n.norm();	20,332✔
303	#else
304	double double_area = 2.0 * triarea(p0, p1, p2);
305	#endif
306	if (double_area > std::numeric_limits<double>::min())	20,332✔
307	{
308	// n /= double_area;
309	n.normalize();	20,332✔
310	G.col(0) = n.cross(p2 - p1) / (double_area);	20,332✔
311	G.col(1) = n.cross(p0 - p2) / (double_area);	20,332✔
312	G.col(2) = n.cross(p1 - p0) / (double_area);	20,332✔
313	}
314	else
315	{
316	G.setZero();	×
317	}
318	}	20,332✔
319
320	void polygon_gradient_matrix(const DenseMatrix& polygon, DenseMatrix& Gpoly)	6,386✔
321	{
322	const int n = (int)polygon.rows();	6,386✔
323
324	// compute position of virtual vertex
325	Eigen::VectorXd vweights;	6,386✔
326	compute_virtual_vertex(polygon, vweights);	6,386✔
327	const Eigen::Vector3d vvertex = polygon.transpose() * vweights;	6,386✔
328
329	DenseMatrix Gfan = DenseMatrix::Zero(3 * n, n + 1);	6,386✔
330	DenseMatrix Gtri(3, 3);	6,386✔
331	int row = 0;	6,386✔
332	for (int i = 0; i < n; ++i)	26,718✔
333	{
334	const int j = (i + 1) % n;	20,332✔
335
336	// build laplace matrix of one triangle
337	triangle_gradient_matrix(polygon.row(i), polygon.row(j), vvertex, Gtri);	20,332✔
338
339	// assemble to matrix for triangle fan
340	for (int k = 0; k < 3; ++k)	81,328✔
341	{
342	Gfan(row + k, i) += Gtri(k, 0);	60,996✔
343	Gfan(row + k, j) += Gtri(k, 1);	60,996✔
344	Gfan(row + k, n) += Gtri(k, 2);	60,996✔
345	}
346
347	row += 3;	20,332✔
348	}
349
350	// build prolongation matrix
351	DenseMatrix P(n + 1, n);	6,386✔
352	P.setIdentity();	6,386✔
353	P.row(n) = vweights;	6,386✔
354
355	// build polygon gradient matrix by sandwiching (from left only)
356	Gpoly = Gfan * P;	6,386✔
357	}	6,386✔
358
359	void divmass_matrix(const SurfaceMesh& mesh, DiagonalMatrix& M)	4✔
360	{
361	// how many virtual triangles will we have after refinement?
362	unsigned int nt = 0;	4✔
363	for (auto f : mesh.faces())	3,173✔
364	{
365	nt += mesh.valence(f);	3,169✔
366	}
367
368	// initialize global matrix
369	M.resize(3 * nt);	4✔
370	auto& diag = M.diagonal();	4✔
371
372	std::vector<Vertex> vertices; // polygon vertices	4✔
373	DenseMatrix polygon; // positions of polygon vertices	4✔
374
375	unsigned int idx = 0;	4✔
376
377	for (const auto f : mesh.faces())	6,342✔
378	{
379	// collect polygon vertices
380	vertices.clear();	3,169✔
381	for (const auto v : mesh.vertices(f))	23,341✔
382	{
383	vertices.push_back(v);	10,086✔
384	}
385	const int n = vertices.size();	3,169✔
386
387	// collect their positions
388	polygon.resize(n, 3);	3,169✔
389	for (int i = 0; i < n; ++i)	13,255✔
390	{
391	polygon.row(i) = (Eigen::Vector3d)mesh.position(vertices[i]);	10,086✔
392	}
393
394	// compute position of virtual vertex
395	Eigen::VectorXd vweights;	3,169✔
396	compute_virtual_vertex(polygon, vweights);	3,169✔
397	const Eigen::Vector3d vvertex = polygon.transpose() * vweights;	3,169✔
398
399	for (int i = 0; i < n; ++i)	13,255✔
400	{
401	const double area =
402	triarea(polygon.row(i), polygon.row((i + 1) % n), vvertex);	10,086✔
403
404	diag[idx++] = area;	10,086✔
405	diag[idx++] = area;	10,086✔
406	diag[idx++] = area;	10,086✔
407	}
408	}	3,169✔
409
410	assert(idx == 3 * nt);	4✔
411	}	4✔
412
413	} // anonymous namespace
414
415	void uniform_mass_matrix(const SurfaceMesh& mesh, DiagonalMatrix& M)	×
416	{
417	const unsigned int n = mesh.n_vertices();	×
418	Eigen::VectorXd diag(n);	×
419	for (auto v : mesh.vertices())	×
420	diag[v.idx()] = mesh.valence(v);	×
421	M = diag.asDiagonal();	×
422	}	×
423
424	void mass_matrix(const SurfaceMesh& mesh, DiagonalMatrix& M)	21✔
425	{
426	const int nv = mesh.n_vertices();	21✔
427	std::vector<Vertex> vertices; // polygon vertices	21✔
428	DenseMatrix polygon; // positions of polygon vertices	21✔
429	DiagonalMatrix Mpoly; // local mass matrix	21✔
430
431	M.setZero(nv);	21✔
432
433	for (const auto f : mesh.faces())	215,123✔
434	{
435	// collect polygon vertices
436	vertices.clear();	107,551✔
437	for (const auto v : mesh.vertices(f))	754,783✔
438	{
439	vertices.push_back(v);	323,616✔
440	}
441	const int n = vertices.size();	107,551✔
442
443	// collect their positions
444	polygon.resize(n, 3);	107,551✔
445	for (int i = 0; i < n; ++i)	431,167✔
446	{
447	polygon.row(i) = (Eigen::Vector3d)mesh.position(vertices[i]);	323,616✔
448	}
449
450	// setup local mass matrix
451	polygon_mass_matrix(polygon, Mpoly);	107,551✔
452
453	// assemble to global mass matrix
454	for (int k = 0; k < n; ++k)	431,167✔
455	{
456	M.diagonal()[vertices[k].idx()] += Mpoly.diagonal()[k];	323,616✔
457	}
458	}
459	}	21✔
460
461	void uniform_laplace_matrix(const SurfaceMesh& mesh, SparseMatrix& L)	2✔
462	{
463	const unsigned int n = mesh.n_vertices();	2✔
464
465	std::vector<Triplet> triplets;	2✔
466	triplets.reserve(8 * n); // conservative estimate for triangle meshes	2✔
467
468	for (auto vi : mesh.vertices())	20✔
469	{
470	Scalar sum_weights = 0.0;	18✔
471	for (auto vj : mesh.vertices(vi))	82✔
472	{
473	sum_weights += 1.0;	64✔
474	triplets.emplace_back(vi.idx(), vj.idx(), 1.0);	64✔
475	}
476	triplets.emplace_back(vi.idx(), vi.idx(), -sum_weights);	18✔
477	}
478
479	L.resize(n, n);	2✔
480	L.setFromTriplets(triplets.begin(), triplets.end());	2✔
481	}	2✔
482
483	void laplace_matrix(const SurfaceMesh& mesh, SparseMatrix& L, bool clamp)	24✔
484	{
485	const int nv = mesh.n_vertices();	24✔
486	const int nf = mesh.n_faces();	24✔
487	std::vector<Vertex> vertices; // polygon vertices	24✔
488	DenseMatrix polygon; // positions of polygon vertices	24✔
489	DenseMatrix Lpoly; // local laplace matrix	24✔
490
491	std::vector<Triplet> triplets;	24✔
492	triplets.reserve(9 * nf); // estimate for triangle meshes	24✔
493
494	for (const auto f : mesh.faces())	416,350✔
495	{
496	// collect polygon vertices
497	vertices.clear();	208,163✔
498	for (const auto v : mesh.vertices(f))	1,458,331✔
499	{
500	vertices.push_back(v);	625,084✔
501	}
502	const int n = vertices.size();	208,163✔
503
504	// collect their positions
505	polygon.resize(n, 3);	208,163✔
506	for (int i = 0; i < n; ++i)	833,247✔
507	{
508	polygon.row(i) = (Eigen::Vector3d)mesh.position(vertices[i]);	625,084✔
509	}
510
511	// setup local laplace matrix
512	polygon_laplace_matrix(polygon, Lpoly);	208,163✔
513
514	// assemble to global laplace matrix
515	for (int j = 0; j < n; ++j)	833,247✔
516	{
517	for (int k = 0; k < n; ++k)	2,502,716✔
518	{
519	triplets.emplace_back(vertices[k].idx(), vertices[j].idx(),	1,877,632✔
520	-Lpoly(k, j));	3,755,264✔
521	}
522	}
523	}
524
525	// build sparse matrix from triplets
526	L.resize(nv, nv);	24✔
527	L.setFromTriplets(triplets.begin(), triplets.end());	24✔
528
529	// clamp negative off-diagonal entries to zero
530	if (clamp)	24✔
531	{
532	for (unsigned int k = 0; k < L.outerSize(); k++)	51,235✔
533	{
534	double diag_offset(0.0);	51,228✔
535
536	for (SparseMatrix::InnerIterator iter(L, k); iter; ++iter)	409,720✔
537	{
538	if (iter.row() != iter.col() && iter.value() < 0.0)	358,492✔
539	{
540	diag_offset += -iter.value();	×
541	iter.valueRef() = 0.0;	×
542	}
543	}
544	for (SparseMatrix::InnerIterator iter(L, k); iter; ++iter)	409,720✔
545	{
546	if (iter.row() == iter.col() && iter.value() < 0.0)	358,492✔
547	iter.valueRef() -= diag_offset;	51,228✔
548	}
549	}
550	}
551	}	24✔
552
553	void gradient_matrix(const SurfaceMesh& mesh, SparseMatrix& G)	10✔
554	{
555	const int nv = mesh.n_vertices();	10✔
556
557	// how many virtual triangles will we have after refinement?
558	// triangles are not refined, other polygons are.
559	unsigned int nt = 0;	10✔
560	for (auto f : mesh.faces())	6,396✔
561	{
562	nt += mesh.valence(f);	6,386✔
563	}
564
565	std::vector<Vertex> vertices; // polygon vertices	10✔
566	DenseMatrix polygon; // positions of polygon vertices	10✔
567	DenseMatrix Gpoly; // local gradient matrix	10✔
568
569	std::vector<Triplet> triplets;	10✔
570	triplets.reserve(9 * nt);	10✔
571
572	unsigned int n_rows = 0;	10✔
573
574	for (const auto f : mesh.faces())	6,396✔
575	{
576	// collect polygon vertices
577	vertices.clear();	6,386✔
578	for (const auto v : mesh.vertices(f))	47,050✔
579	{
580	vertices.push_back(v);	20,332✔
581	}
582	const int n = vertices.size();	6,386✔
583
584	// collect their positions
585	polygon.resize(n, 3);	6,386✔
586	for (int i = 0; i < n; ++i)	26,718✔
587	{
588	polygon.row(i) = (Eigen::Vector3d)mesh.position(vertices[i]);	20,332✔
589	}
590
591	// setup local element matrix
592	polygon_gradient_matrix(polygon, Gpoly);	6,386✔
593
594	// assemble to global matrix
595	for (int j = 0; j < Gpoly.cols(); ++j)	26,718✔
596	{
597	for (int i = 0; i < Gpoly.rows(); ++i)	217,408✔
598	{
599	triplets.emplace_back(n_rows + i, vertices[j].idx(),	197,076✔
600	Gpoly(i, j));	197,076✔
601	}
602	}
603
604	n_rows += Gpoly.rows();	6,386✔
605	}
606	assert(n_rows == 3 * nt);	10✔
607
608	// build sparse matrix from triplets
609	G.resize(n_rows, nv);	10✔
610	G.setFromTriplets(triplets.begin(), triplets.end());	10✔
611	}	10✔
612
613	void divergence_matrix(const SurfaceMesh& mesh, SparseMatrix& D)	4✔
614	{
615	SparseMatrix G;	4✔
616	gradient_matrix(mesh, G);	4✔
617	DiagonalMatrix M;	4✔
618	divmass_matrix(mesh, M);	4✔
619	D = -G.transpose() * M;	4✔
620	}	4✔
621
622	} // namespace pmp

pmp-library / pmp-library / 21332822416

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous