laplacian_2d_grid.h

Go to the documentation of this file.

/*
 * Copyright 2023 MusicScience37 (Kenta Kabashima)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#pragma once
 
#include <algorithm>
#include <cstddef>
#include <vector>
 
#include <Eigen/SparseCore>
 
#include "num_collect/base/concepts/sparse_matrix.h"
#include "num_collect/base/index_type.h"
#include "num_collect/util/assert.h"
 
namespace num_prob_collect::linear {
 
template <num_collect::base::concepts::sparse_matrix Matrix>
class laplacian_2d_grid {
public:
    using matrix_type = Matrix;
 
    using scalar_type = typename Matrix::Scalar;
 
    laplacian_2d_grid(num_collect::index_type grid_rows,
        num_collect::index_type grid_cols,
        typename Matrix::value_type grid_width)
        : grid_rows_(grid_rows),
          grid_cols_(grid_cols),
          mat_size_(grid_rows_ * grid_cols_),
          grid_width_(grid_width),
          mat_(mat_size_, mat_size_) {
        NUM_COLLECT_ASSERT(grid_rows > 0);
        NUM_COLLECT_ASSERT(grid_cols > 0);
        NUM_COLLECT_ASSERT(grid_width > static_cast<scalar_type>(0));
        calc();
    }
 
    [[nodiscard]] auto grid_rows() const noexcept -> num_collect::index_type {
        return grid_rows_;
    }
 
    [[nodiscard]] auto grid_cols() const noexcept -> num_collect::index_type {
        return grid_cols_;
    }
 
    [[nodiscard]] auto mat_size() const noexcept -> num_collect::index_type {
        return mat_size_;
    }
 
    [[nodiscard]] auto diag_coeff() const noexcept -> scalar_type {
        return diag_coeff_;
    }
 
    [[nodiscard]] auto off_diag_coeff() const noexcept -> scalar_type {
        return off_diag_coeff_;
    }
 
    [[nodiscard]] auto mat() const noexcept -> const Matrix& { return mat_; }
 
    [[nodiscard]] auto index(num_collect::index_type x,
        num_collect::index_type y) const noexcept -> num_collect::index_type {
        return x + grid_cols_ * y;
    }
 
private:
    void calc() {
        const scalar_type inv_area =
            static_cast<scalar_type>(1) / (grid_width_ * grid_width_);
        // NOLINTNEXTLINE
        diag_coeff_ = static_cast<scalar_type>(8) /
            static_cast<scalar_type>(3) * inv_area;
        // NOLINTNEXTLINE
        off_diag_coeff_ = -inv_area / static_cast<scalar_type>(3);
 
        std::vector<Eigen::Triplet<scalar_type>> triplets;
        const std::size_t approx_elements =
            9U * static_cast<std::size_t>(mat_size_);
        triplets.reserve(approx_elements);
        for (num_collect::index_type xi = 0; xi < grid_cols_; ++xi) {
            for (num_collect::index_type yi = 0; yi < grid_rows_; ++yi) {
                const num_collect::index_type i = index(xi, yi);
                for (num_collect::index_type
                         xj = std::max<num_collect::index_type>(xi - 1, 0),
                         xj_end = std::min<num_collect::index_type>(
                             xi + 2, grid_cols_);
                    xj < xj_end; ++xj) {
                    for (num_collect::index_type
                             yj = std::max<num_collect::index_type>(yi - 1, 0),
                             yj_end = std::min<num_collect::index_type>(
                                 yi + 2, grid_rows_);
                        yj < yj_end; ++yj) {
                        const num_collect::index_type j = index(xj, yj);
                        const scalar_type coeff =
                            (i == j) ? diag_coeff_ : off_diag_coeff_;
                        triplets.emplace_back(
                            static_cast<int>(i), static_cast<int>(j), coeff);
                    }
                }
            }
        }
 
        mat_.setFromTriplets(triplets.begin(), triplets.end());
    }
 
    num_collect::index_type grid_rows_;
 
    num_collect::index_type grid_cols_;
 
    num_collect::index_type mat_size_;
 
    scalar_type grid_width_;
 
    scalar_type diag_coeff_{};
 
    scalar_type off_diag_coeff_{};
 
    Matrix mat_;
};
 
}  // namespace num_prob_collect::linear