base.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 BEM_RWG_ASSEMBLER_BASE_H
#define BEM_RWG_ASSEMBLER_BASE_H
 
#include "types.hpp"
#include "geometry/point_cloud.hpp"
#include "geometry/mesh/triangle_mesh.hpp"
#include "geometry/primitives/triangle.hpp"
#include "matrix/base.hpp"
#include "rwg/operators/base.hpp"
#include "rwg/excitations/base.hpp"
#include "rwg/projectors/base.hpp"
 
 
namespace bem::rwg
{
 
template <uint8_t obs_num_dof, uint8_t src_num_dof>
class OperatorAssemblerBase
{
 
    static_assert((obs_num_dof > 0), "OperatorAssemblerBase: `obs_num_dof` must be greater than 0.");
    static_assert((src_num_dof > 0), "OperatorAssemblerBase: `src_num_dof` must be greater than 0.");
 
public:
 
    OperatorAssemblerBase(
        const TriangleMesh<3>& obs_mesh,
        const TriangleMesh<3>& src_mesh,
        ConstEigRef<EigMatNX<Index, 2>> elem_pairs
        ):
            obs_mesh_(obs_mesh),
            src_mesh_(src_mesh),
            elem_pairs_(elem_pairs) {};
 
 
    OperatorAssemblerBase(
        const TriangleMesh<3>& obs_mesh,
        const TriangleMesh<3>& src_mesh
        ):
            obs_mesh_(obs_mesh),
            src_mesh_(src_mesh),
            elem_pairs_(make_pairs(obs_mesh, src_mesh)) {};
 
 
    OperatorAssemblerBase(const TriangleMesh<3>& mesh, ConstEigRef<EigMatNX<Index, 2>> elem_pairs):
        obs_mesh_(mesh),
        src_mesh_(mesh),
        elem_pairs_(elem_pairs) {};
 
 
    OperatorAssemblerBase(const TriangleMesh<3>& mesh):
        obs_mesh_(mesh),
        src_mesh_(mesh),
        elem_pairs_(make_pairs(mesh, mesh)) {};
 
 
    template <typename OperatorType>
    void assemble(
        MatrixBase<Complex>& mat,
        OperatorType op,
        const Complex k
        )
    {
 
        static_assert(
            std::is_base_of<OperatorBase<obs_num_dof, src_num_dof>, OperatorType>::value,
            "OperatorAssemblerBase::assemble(): `OperatorType` must derive from `OperatorBase<obs_num_dof, src_num_dof>`"
            );
 
        prep_matrix(mat);
 
#pragma omp parallel for firstprivate(op)
        for (Index ii = 0; ii < elem_pairs_.cols(); ++ii)
        {
            Triangle<3> obs_tri = obs_mesh_.elem_primitive(elem_pairs_(0, ii));
            Triangle<3> src_tri = src_mesh_.elem_primitive(elem_pairs_(1, ii));
 
            EigMatMN<Complex, obs_num_dof, src_num_dof> values = op.compute(k, obs_tri, src_tri);
 
#pragma omp critical
            fill_matrix(mat, elem_pairs_.col(ii), values);
        }
 
        mat.assemble();
        return;
    };
 
 
    virtual void prep_matrix(MatrixBase<Complex>& mat) = 0;
 
 
    virtual void fill_matrix(
        MatrixBase<Complex>& mat,
        ConstEigRef<EigColVecN<Index, 2>> elem_pair,
        ConstEigRef<EigMatMN<Complex, obs_num_dof, src_num_dof>> values
        ) = 0;
 
 
    static EigMatNX<Index, 2> make_pairs(
        const TriangleMesh<3>& obs_mesh,
        const TriangleMesh<3>& src_mesh
        )
    {
        Index num_pairs = obs_mesh.num_elems() * src_mesh.num_elems();
        EigMatNX<Index, 2> pairs = EigMatNX<Index, 2>::Zero(2, num_pairs);
 
        for (Index ii = 0; ii < num_pairs; ++ii)
        {
            pairs(0, ii) = ii / src_mesh.num_elems();
            pairs(1, ii) = ii % src_mesh.num_elems();
        }
 
        return pairs;
    };
 
 
    static EigMatNX<Index, 2> make_pairs(
        ConstEigRef<EigRowVec<Index>> obs_elems,
        ConstEigRef<EigRowVec<Index>> src_elems
        )
    {
        Index num_pairs = obs_elems.size() * src_elems.size();
        EigMatNX<Index, 2> pairs = EigMatNX<Index, 2>::Zero(2, num_pairs);
 
        for (Index ii = 0; ii < obs_elems.size(); ++ii)
        {
            for (Index jj = 0; jj < src_elems.size(); ++jj)
            {
                pairs(0, jj + ii * src_elems.size()) = obs_elems[ii];
                pairs(1, jj + ii * src_elems.size()) = src_elems[jj];
            }
        }
 
        return pairs;
    };
 
 
protected:
 
    const TriangleMesh<3>& obs_mesh_;
    const TriangleMesh<3>& src_mesh_;
    const EigMatNX<Index, 2> elem_pairs_;
 
};
 
 
template <uint8_t obs_num_dof>
class ExcitationAssemblerBase
{
 
    static_assert((obs_num_dof > 0), "ExcitationAssemblerBase: `obs_num_dof` must be greater than 0.");
 
public:
 
    ExcitationAssemblerBase(const TriangleMesh<3>& mesh):
        obs_mesh_(mesh),
        elems_(EigRowVec<Index>::LinSpaced(obs_mesh_.num_elems(), 0, obs_mesh_.num_elems() - 1)) {};
 
 
    ExcitationAssemblerBase(const TriangleMesh<3>& mesh, ConstEigRef<EigRowVec<Index>> elems):
        obs_mesh_(mesh),
        elems_(elems) {};
 
 
    virtual void assemble(
        MatrixBase<Complex>& mat,
        ExcitationBase<obs_num_dof>& exc,
        const Complex k
        ) = 0;
 
 
protected:
 
    const TriangleMesh<3>& obs_mesh_;
    const EigRowVec<Index> elems_;
 
};
 
 
template <uint8_t obs_dim, uint8_t src_num_dof>
class ProjectorAssemblerBase
{
 
    static_assert((obs_dim > 0), "`obs_dim` must be greater than 0.");
    static_assert((src_num_dof > 0), "`src_num_dof` must be greater than 0.");
 
public:
 
    ProjectorAssemblerBase(
        const PointCloud<3>& obs_cloud,
        const TriangleMesh<3>& src_mesh,
        ConstEigRef<EigRowVec<Index>> elems
        ):
            obs_cloud_(obs_cloud),
            src_mesh_(src_mesh),
            elems_(elems) {};
 
 
    ProjectorAssemblerBase(
        const PointCloud<3>& obs_cloud,
        const TriangleMesh<3>& src_mesh
        ):
            obs_cloud_(obs_cloud),
            src_mesh_(src_mesh),
            elems_(EigRowVec<Index>::LinSpaced(src_mesh.num_elems(), 0, src_mesh.num_elems() - 1)) {};
 
 
    virtual void assemble(
        MatrixBase<Complex>& mat,
        ProjectorBase<src_num_dof>& op,
        const Complex k
        ) = 0;
 
 
protected:
 
    const PointCloud<3>& obs_cloud_;
    const TriangleMesh<3>& src_mesh_;
    const EigRowVec<Index> elems_;
 
};
 
}
 
#endif