EulerKineticSolverMultiD.cpp

#include <scheme/EulerKineticSolverMultiD.hpp>

#include <language/utils/EvaluateAtPoints.hpp>
#include <language/utils/InterpolateItemArray.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshFaceBoundary.hpp>
#include <mesh/MeshFlatFaceBoundary.hpp>
#include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/MeshVariant.hpp>
#include <scheme/DiscreteFunctionDPkVariant.hpp>
#include <scheme/DiscreteFunctionDPkVector.hpp>
#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
#include <scheme/DiscreteFunctionP0Vector.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
#include <scheme/OutflowBoundaryConditionDescriptor.hpp>
#include <scheme/PolynomialReconstruction.hpp>
#include <scheme/PolynomialReconstructionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <scheme/WallBoundaryConditionDescriptor.hpp>

#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/QuadratureManager.hpp>
#include <geometry/LineTransformation.hpp>
#include <tuple>

std::vector<TinyVector<2>>
getLambdaVector2D(const std::shared_ptr<const MeshVariant>& mesh,
                  const double& lambda,
                  const size_t& L,
                  const size_t& M)
{
  return std::visit(
    [&](auto&& p_mesh) -> std::vector<TinyVector<2>> {
      using MeshType = mesh_type_t<decltype(p_mesh)>;
      if constexpr (is_polygonal_mesh_v<MeshType> and (MeshType::Dimension == 2)) {
        const double pi = std::acos(-1);

        std::vector<TinyVector<2>> lambda_vector(4 * L * M);   //!!!!!!!!!!!!!!!!!!!!

        for (size_t l = 1; l < L + 1; ++l) {
          for (size_t m = 1; m < 4 * M + 1; ++m) {
            size_t i_wave         = (l - 1) * 4 * M + m - 1;
            double ll             = l;
            lambda_vector[i_wave] = TinyVector<2>{(ll / L) * lambda * std::cos((m * pi) / (2 * M)),
                                                  (ll / L) * lambda * std::sin((m * pi) / (2 * M))};
          }
        }
        // lambda_vector[4 * L * M] = TinyVector<2>{0, 0};   //!!!!!!!!!!!!!!!!!!!!

        return lambda_vector;
      } else {
        throw NotImplementedError("Invalid mesh type for MultiD Euler Kinetic solver");
      }
    },
    mesh->variant());
}

std::vector<TinyVector<3>>
getLambdaVector3D(const std::shared_ptr<const MeshVariant>& mesh, const double& lambda)
{
  return std::visit(
    [&](auto&& p_mesh) -> std::vector<TinyVector<3>> {
      using MeshType = mesh_type_t<decltype(p_mesh)>;
      if constexpr (is_polygonal_mesh_v<MeshType> and (MeshType::Dimension == 3)) {
        std::vector<TinyVector<3>> lambda_vector{6};