Method getBdLaplacianMatrix

tools/Assembler/ref/Assembler3D.cpp:323–378

}

ZSparseMatrix Assembler3D::getBdLaplacianMatrix() {
    typedef Eigen::Triplet<double> T;
    std::vector<T> triplets;

    size_t num_bdv = m_mesh->getNbrBoundaryNodes();
    size_t num_bdf = m_mesh->getNbrBoundaryFaces();

    // Compute lumped mass
    VectorF lumped_mass(num_bdv);
    for (size_t i=0; i<num_bdv; i++) {
        VectorI neighbor_faces = m_mesh->getBoundaryNodeAdjacentBoundaryFaces(i);
        assert(neighbor_faces.size() == 3);

        double total_weight =
            m_mesh->getBoundaryFaceArea(neighbor_faces[0]) +
            m_mesh->getBoundaryFaceArea(neighbor_faces[1]) +
            m_mesh->getBoundaryFaceArea(neighbor_faces[2]);
        lumped_mass[i] = total_weight / 3.0;
    }

    // Compute laplacian matrix.
    for (size_t i=0; i<num_bdf; i++) {
        VectorI face = m_mesh->getBoundaryFace(i);
        assert(face.size() == 3);

        size_t v[3] = {
            m_mesh->getBoundaryIndex(face[0]),
            m_mesh->getBoundaryIndex(face[1]),
            m_mesh->getBoundaryIndex(face[2])
        };
        Eigen::Vector3d p[3] = {
            Eigen::Vector3d(m_mesh->getNode(face[0])),
            Eigen::Vector3d(m_mesh->getNode(face[1])),
            Eigen::Vector3d(m_mesh->getNode(face[2]))
        };

        for (size_t j=0; j<3; j++) {
            size_t i0 = j;
            size_t i1 = (j+1)%3;
            size_t i2 = (j+2)%3;

            double cos_angle = (p[i1] - p[i0]).dot(p[i2] - p[i0]);
            double sin_angle = ((p[i1] - p[i0]).cross(p[i2] - p[i0])).norm();
            double weight = 0.5 * cos_angle / sin_angle;

            triplets.push_back(T(v[i1], v[i1], -weight / lumped_mass[v[i1]]));
            triplets.push_back(T(v[i1], v[i2],  weight / lumped_mass[v[i1]]));
            triplets.push_back(T(v[i2], v[i1],  weight / lumped_mass[v[i2]]));
            triplets.push_back(T(v[i2], v[i2], -weight / lumped_mass[v[i2]]));
        }
    }

    Eigen::SparseMatrix<double> Lb = Eigen::SparseMatrix<double>(num_bdv, num_bdv);
    Lb.setFromTriplets(triplets.begin(), triplets.end());
    return ZSparseMatrix(Lb);
}