Method getBdLaplacianMatrix

tools/Assembler/ref/Assembler2D.cpp:239–276

}

ZSparseMatrix Assembler2D::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() == 2);

        double total_weight =
            m_mesh->getBoundaryFaceArea(neighbor_faces[0]) +
            m_mesh->getBoundaryFaceArea(neighbor_faces[1]);
        lumped_mass[i] = 0.5 * total_weight;
    }

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

        double l = m_mesh->getBoundaryFaceArea(i);
        size_t v1 = m_mesh->getBoundaryIndex(face[0]);
        size_t v2 = m_mesh->getBoundaryIndex(face[1]);
        double weight = 1.0 / l;
        triplets.push_back(T(v1, v1, -weight / lumped_mass[v1]));
        triplets.push_back(T(v1, v2,  weight / lumped_mass[v1]));
        triplets.push_back(T(v2, v1,  weight / lumped_mass[v2]));
        triplets.push_back(T(v2, v2, -weight / lumped_mass[v2]));
    }

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