plane_wave.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_EXC_PLANE_WAVE_H
#define BEM_RWG_EXC_PLANE_WAVE_H
 
#include <stdexcept>
#include <memory>
 
#include "types.hpp"
#include "quadrature/triangle/gauss.hpp"
 
#include "rwg/excitations/base.hpp"
 
 
namespace bem
{
// Forward declarations
template <uint8_t dim> class Triangle;
}
 
 
namespace bem::rwg
{
 
template <uint8_t obs_num_dof>
class PlaneWaveBase: public ExcitationBase<obs_num_dof>
{
public:
 
    PlaneWaveBase(
        ConstEigRef<EigMatNX<Float, 3>> dir,
        ConstEigRef<EigMatNX<Float, 3>> pol,
        ConstEigRef<EigMatNX<Float, 3>> pos,
        ConstEigRef<EigRowVec<Complex>> amp,
        const uint8_t quadrature_order = 4
        ):
            dir_(dir.colwise().normalized()),
            pol_(pol.colwise().normalized()),
            pos_(pos),
            amp_(amp),
            tri_quad_(std::make_shared<GaussTriangleQuadrature<3>> (
                GaussTriangleQuadrature<3> (quadrature_order)
                ))
    {
        if (
            dir_.cols() != pol_.cols() ||
            dir_.cols() != pos_.cols() ||
            dir_.cols() != amp.size()
            )
        {
            throw std::invalid_argument(
                "PlaneWaveBase: `dir`, `pol`, `pos`, and `amp` must have the same number of columns."
                );
        }
 
        return;
    };
 
 
    Index num_excitations() const override { return dir_.cols(); };
 
 
protected:
 
    EigMatNX<Complex, 3> eval(const Complex k, ConstEigRef<EigMatNX<Float, 3>> r_obs, const Index idx = 0)
    {
        EigMatNX<Float, 3> r_diff = -pos_.col(idx).replicate(1, r_obs.cols()) + r_obs;
        EigRowVec<Float> r_dir = dir_.col(idx).transpose() * r_diff;
        return amp_[idx] * pol_.col(idx) * Eigen::exp(-J * k * r_dir.array()).matrix();
    };
 
 
    const EigMatNX<Float, 3> dir_, pol_, pos_;
    const EigRowVec<Complex> amp_;
    std::shared_ptr<GaussTriangleQuadrature<3>> tri_quad_;
 
};
 
 
class RwgPlaneWave: public PlaneWaveBase<3>
{
public:
 
    using PlaneWaveBase<3>::PlaneWaveBase;
 
    EigMatNX<Complex, 3> compute(const Complex k, const Triangle<3>& obs_tri) override;
 
};
 
 
class NxRwgPlaneWave: public PlaneWaveBase<3>
{
public:
 
    using PlaneWaveBase<3>::PlaneWaveBase;
 
    EigMatNX<Complex, 3> compute(const Complex k, const Triangle<3>& obs_tri) override;
 
};
 
}
 
#ifndef BEM_LINKED
#include "rwg/excitations/plane_wave.cpp"
#endif
 
#endif