adaptive.hpp

Go to the documentation of this file.

// OpenBEM - Copyright (C) 2026 Shashwat Sharma
 
// This file is part of OpenBEM.
 
// OpenBEM is free software: you can redistribute it and/or modify it under the terms of the
// GNU General Public License as published by the Free Software Foundation, either version 3
// of the License, or (at your option) any later version.
 
// You should have received a copy of the GNU General Public License along with OpenBEM.
// If not, see <https://www.gnu.org/licenses/>.
 
 
#ifndef ADAPTIVE_TRI_QUAD_H
#define ADAPTIVE_TRI_QUAD_H
 
#include <cassert>
 
#include "types.hpp"
#include "constants.hpp"
 
#include "geometry/primitives/triangle.hpp"
#include "quadrature/utility.hpp"
#include "quadrature/triangle/base.hpp"
 
 
namespace bem
{
 
const uint16_t ADAPTIVE_TRI_DEFAULT_MAX_LEVELS = 200;
const Float ADAPTIVE_TRI_DEFAULT_TOL = 1e-3;
 
template <uint8_t dim>
class AdaptiveTriangleQuadrature: public TriangleQuadratureBase<dim>
{
 
    using base = TriangleQuadratureBase<dim>;
 
public:
 
    void set_max_levels(const uint16_t max_levels)
    {
        max_levels_ = max_levels;
        base::points_weights_computed_ = false;
        converged_ = false;
        return;
    };
 
 
    uint16_t max_levels() const { return max_levels_; };
 
 
    void set_tol(const Float tol)
    {
        tol_ = tol;
        base::points_weights_computed_ = false;
        converged_ = false;
        return;
    };
 
 
    Float tol() const { return tol_; };
 
 
    void compute_points_weights(
        const Triangle<dim>& tri,
        std::function<EigRowVec<Complex> (ConstEigRef<EigMatNX<Float, dim>>)> eval = {}
        ) override;
 
 
    bool converged() const
    {
        if (!base::points_weights_computed_)
            throw std::runtime_error(
                "AdaptiveTriangleQuadrature::converged(): must call `compute_points_weights()` first.");
        return converged_;
    };
 
 
    uint16_t converged_iter() const { return converged_iter_; };
 
 
private:
 
    Complex run_recursion(
        int& level,
        const Triangle<dim>& tri,
        std::function<EigRowVec<Complex> (ConstEigRef<EigMatNX<Float, dim>>)> eval,
        const uint8_t longest_edge
        );
 
 
    static void get_5_points(
        EigMatMN<Float, dim, 5>& p,
        const Triangle<dim>& tri,
        const uint8_t longest_edge
        );
 
 
    void get_subtris(
        EigMatMN<Float, dim, 3>& v_subtri1,
        EigMatMN<Float, dim, 3>& v_subtri2,
        const Triangle<dim>& tri,
        const uint8_t longest_edge
        );
 
 
    bool check_convergence(const Complex val_coarse, const Complex val_fine)
    {
        bool converged = compare_with_tol(val_fine, val_coarse, tol_, 3);
 
        assert((converged || val_fine != zero) &&
               "AdaptiveTriangleQuadrature::check_convergence(): divide by zero.");
        // converged = converged || check_convergence(val_fine);
        // check_convergence(val_fine);
        // std::cout << val_fine << ", " << val_coarse << ", " << std::endl;
        return converged;
    }
 
 
    bool check_convergence(const Complex val)
    {
        // Complex I_new = I_prev_ + val;
        // bool converged = compare_with_tol(I_new, I_prev_, tol, 3);
        // I_prev_ += val;
 
        Complex I_new = base::weights_.dot(vals_) + val;
        bool converged = compare_with_tol(I_new, I_prev_, tol_, 3);
 
        if (converged)
            I_prev_ = I_new;
 
        return converged;
    }
 
 
    uint16_t total_recursions()
    { return total_recursions_; };
 
 
    void reset_recursion_count()
    {
        total_recursions_ = 0;
        return;
    };
 
 
    uint16_t max_levels_ = ADAPTIVE_TRI_DEFAULT_MAX_LEVELS;
    Float tol_ = ADAPTIVE_TRI_DEFAULT_TOL;
    bool converged_ = false;
    uint16_t converged_iter_ = 0;
 
    const EigColVecN<Float, 10> weights10_ = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 2.0};
    const EigColVecN<Float, 5> weights5_ = {1.0, 1.0, 1.0, 1.0, 2.0};
 
    EigMatMN<Float, dim, 5> p_main_;
    EigMatMN<Float, dim, 6> p_sub_;
    EigRowVecN<Complex, 5> vals_main_, vals_subtri1_, vals_subtri2_;
    EigRowVecN<Complex, 6> vals_sub_;
    EigRowVec<Complex> vals_;
 
    EigRowVecN<uint32_t, 2> idx_extra1_, idx_extra2_;
    Complex I_prev_ = 0.0;
    uint16_t total_recursions_ = 0;
 
};
 
}
 
#ifndef BEM_LINKED
#include "quadrature/triangle/adaptive.cpp"
#endif
 
#endif