Method getBoundaryBlurMatrix

tools/Assembler/ref/Assembler.cpp:65–118

};

ZSparseMatrix Assembler::getBoundaryBlurMatrix(double radius) {
    typedef std::priority_queue<Item, std::vector<Item>, ItemComp> Queue;
    typedef Eigen::Triplet<double> T;
    std::vector<T> triplets;

    size_t num_bd_vertices = m_mesh->getNbrBoundaryNodes();
    double sd = radius / 3.0;
    const double A = 2.50662827463; // A = sqrt(2*pi);

    std::vector<int> visited(num_bd_vertices, -1);

    for (size_t i=0; i<num_bd_vertices; i++) {
        size_t seed = i;
        Queue Q;
        Q.push(Item(i, 0.0));
        VectorF seed_n = m_mesh->getBoundaryNodeNormal(seed);

        while (!Q.empty()) {
            Item cur_item = Q.top();
            Q.pop();

            // Check if visited before
            if (visited[cur_item.first] == seed) continue;
            visited[cur_item.first] = seed;

            size_t curr_glob_idx = m_mesh->getBoundaryNode(cur_item.first);
            VectorF curr_n = m_mesh->getBoundaryNodeNormal(cur_item.first);
            double normal_factor = fmax(0.0, curr_n.dot(seed_n));

            // Use gaussian distribution as weights
            triplets.push_back(T(seed, cur_item.first,
                        exp(-0.5 * pow(cur_item.second / sd, 2)) /
                        (sd*A) * normal_factor));

            VectorI neighbors = m_mesh->getBoundaryNodeAdjacentBoundaryNodes(
                    cur_item.first);
            for (size_t j=0; j<neighbors.size(); j++) {
                size_t next = neighbors[j];
                if (visited[next] == seed) continue;

                size_t next_glob_idx = m_mesh->getBoundaryNode(next);
                double next_dist = cur_item.second +
                    (m_mesh->getNode(curr_glob_idx) -
                     m_mesh->getNode(next_glob_idx)).norm();
                if (next_dist > 2*sd) continue;
                Q.push(Item(next, next_dist));
            }
        }
    }

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

VectorF Assembler::getDivergence(double* vector_field, int num_vector_field) {
    size_t num_node = m_mesh->getNbrNodes();