Skip to content
Snippets Groups Projects
Select Git revision
  • 7f0b0ce06a57f633ce1f1e3318a6ae380d294cd5
  • develop default protected
  • feature/gmsh-reader
  • origin/stage/bouguettaia
  • feature/kinetic-schemes
  • feature/reconstruction
  • feature/local-dt-fsi
  • feature/composite-scheme-sources
  • feature/composite-scheme-other-fluxes
  • feature/serraille
  • feature/variational-hydro
  • feature/composite-scheme
  • hyperplastic
  • feature/polynomials
  • feature/gks
  • feature/implicit-solver-o2
  • feature/coupling_module
  • feature/implicit-solver
  • feature/merge-local-dt-fsi
  • master protected
  • feature/escobar-smoother
  • v0.5.0 protected
  • v0.4.1 protected
  • v0.4.0 protected
  • v0.3.0 protected
  • v0.2.0 protected
  • v0.1.0 protected
  • Kidder
  • v0.0.4 protected
  • v0.0.3 protected
  • v0.0.2 protected
  • v0 protected
  • v0.0.1 protected
33 results

ScalarNodalScheme.cpp

Blame
  • BoundaryIntegralReconstructionMatrixBuilder.cpp 8.25 KiB
    #include <scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp>
    
    #include <analysis/GaussLegendreQuadratureDescriptor.hpp>
    #include <analysis/QuadratureManager.hpp>
    #include <geometry/SymmetryUtils.hpp>
    #include <mesh/Connectivity.hpp>
    #include <mesh/Mesh.hpp>
    #include <mesh/MeshDataManager.hpp>
    #include <scheme/DiscreteFunctionDPk.hpp>
    
    template <MeshConcept MeshTypeT>
    template <typename ConformTransformationT>
    void
    PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::_computeEjkBoundaryMean(
      const QuadratureFormula<MeshType::Dimension - 1>& quadrature,
      const ConformTransformationT& T,
      const Rd& Xj,
      const double inv_Vi,
      SmallArray<double>& mean_of_ejk)
    {
      const double velocity_perp_e1 = T.velocity()[1] * inv_Vi;
    
      for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
        const double wq          = quadrature.weight(i_q);
        const TinyVector<1> xi_q = quadrature.point(i_q);
    
        const Rd X_Xj = T(xi_q) - Xj;
    
        const double x_xj = X_Xj[0];
        const double y_yj = X_Xj[1];
    
        {
          size_t k                                   = 0;
          m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = x_xj * wq * velocity_perp_e1;
          for (; k <= m_polynomial_degree; ++k) {
            m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k] =
              x_xj * m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k - 1] * m_antiderivative_coef[k];   // ((1. * k) / (k + 1));
          }
    
          for (size_t i_y = 1; i_y <= m_polynomial_degree; ++i_y) {
            const size_t begin_i_y_1 = ((i_y - 1) * (2 * m_polynomial_degree - i_y + 4)) / 2;
            for (size_t l = 0; l <= m_polynomial_degree - i_y; ++l) {
              m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = y_yj * m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[begin_i_y_1 + l];
            }
          }
        }
    
        for (size_t k = 1; k < m_basis_dimension; ++k) {
          mean_of_ejk[k - 1] += m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k];
        }
      }
    }
    
    template <MeshConcept MeshTypeT>
    void
    PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::_computeEjkMeanByBoundary(
      const Rd& Xj,
      const CellId& cell_id,
      SmallArray<double>& mean_of_ejk)
    {
      const auto& quadrature =
        QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
    
      const double inv_Vi = 1. / m_Vj[cell_id];
    
      mean_of_ejk.fill(0);
      auto cell_face_list = m_cell_to_face_matrix[cell_id];
      for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
        bool is_reversed = m_cell_face_is_reversed[cell_id][i_face];
    
        const FaceId face_id = cell_face_list[i_face];
        auto face_node_list  = m_face_to_node_matrix[face_id];
        if (is_reversed) {
          const LineTransformation<2> T{m_xr[face_node_list[1]], m_xr[face_node_list[0]]};
          _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
        } else {
          const LineTransformation<2> T{m_xr[face_node_list[0]], m_xr[face_node_list[1]]};
          _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
        }
      }
    }
    
    template <MeshConcept MeshTypeT>
    void
    PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
      MeshTypeT>::_computeEjkMeanByBoundaryInSymmetricCell(const Rd& origin,
                                                           const Rd& normal,
                                                           const Rd& Xj,
                                                           const CellId& cell_id,
                                                           SmallArray<double>& mean_of_ejk)
    {
      const auto& quadrature =
        QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
    
      const double inv_Vi = 1. / m_Vj[cell_id];
    
      mean_of_ejk.fill(0);
      auto cell_face_list = m_cell_to_face_matrix[cell_id];
      for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
        bool is_reversed = m_cell_face_is_reversed[cell_id][i_face];
    
        const FaceId face_id = cell_face_list[i_face];
        auto face_node_list  = m_face_to_node_matrix[face_id];
    
        const auto x0 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[1]]);
        const auto x1 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[0]]);
    
        if (is_reversed) {
          const LineTransformation<2> T{x1, x0};
          _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
        } else {
          const LineTransformation<2> T{x0, x1};
          _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
        }
      }
    }
    
    template <MeshConcept MeshTypeT>
    void
    PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::build(
      const CellId cell_j_id,
      ShrinkMatrixView<SmallMatrix<double>>& A)
    {
      if constexpr (MeshType::Dimension == 2) {
        const auto& stencil_cell_list = m_stencil_array[cell_j_id];
    
        const Rd& Xj = m_xj[cell_j_id];
    
        _computeEjkMeanByBoundary(Xj, cell_j_id, m_mean_j_of_ejk);
    
        size_t index = 0;
        for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
          const CellId cell_i_id = stencil_cell_list[i];
    
          _computeEjkMeanByBoundary(Xj, cell_i_id, m_mean_i_of_ejk);
    
          for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
            A(index, l) = m_mean_i_of_ejk[l] - m_mean_j_of_ejk[l];
          }
        }
        for (size_t i_symmetry = 0; i_symmetry < m_stencil_array.symmetryBoundaryStencilArrayList().size(); ++i_symmetry) {
          auto& ghost_stencil  = m_stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
          auto ghost_cell_list = ghost_stencil[cell_j_id];
    
          const Rd& origin = m_symmetry_origin_list[i_symmetry];
          const Rd& normal = m_symmetry_normal_list[i_symmetry];
    
          for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
            const CellId cell_i_id = ghost_cell_list[i];
    
            _computeEjkMeanByBoundaryInSymmetricCell(origin, normal, Xj, cell_i_id, m_mean_i_of_ejk);
    
            for (size_t l = 0; l < m_basis_dimension - 1; ++l) {
              A(index, l) = m_mean_i_of_ejk[l] - m_mean_j_of_ejk[l];
            }
          }
        }
      } else {
        throw NotImplementedError("invalid mesh dimension");
      }
    }
    
    template <MeshConcept MeshTypeT>
    PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
      MeshTypeT>::BoundaryIntegralReconstructionMatrixBuilder(const MeshType& mesh,
                                                              const size_t polynomial_degree,
                                                              const SmallArray<const Rd>& symmetry_origin_list,
                                                              const SmallArray<const Rd>& symmetry_normal_list,
                                                              const CellToCellStencilArray& stencil_array)
      : m_mesh{mesh},
        m_basis_dimension{
          DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(polynomial_degree)},
        m_polynomial_degree{polynomial_degree},
        m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek{m_basis_dimension},
        m_mean_j_of_ejk{m_basis_dimension - 1},
        m_mean_i_of_ejk{m_basis_dimension - 1},
        m_cell_to_face_matrix{mesh.connectivity().cellToFaceMatrix()},
        m_face_to_node_matrix{mesh.connectivity().faceToNodeMatrix()},
        m_cell_face_is_reversed{mesh.connectivity().cellFaceIsReversed()},
        m_stencil_array{stencil_array},
        m_symmetry_origin_list{symmetry_origin_list},
        m_symmetry_normal_list{symmetry_normal_list},
        m_Vj{MeshDataManager::instance().getMeshData(mesh).Vj()},
        m_xj{MeshDataManager::instance().getMeshData(mesh).xj()},
        m_xr{mesh.xr()}
    {
      if constexpr (MeshType::Dimension == 2) {
        SmallArray<double> antiderivative_coef(m_polynomial_degree + 1);
        for (size_t k = 0; k < antiderivative_coef.size(); ++k) {
          antiderivative_coef[k] = ((1. * k) / (k + 1));
        }
    
        m_antiderivative_coef = antiderivative_coef;
      }
    }
    
    template void PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<2>>::build(
      const CellId,
      ShrinkMatrixView<SmallMatrix<double>>&);
    
    template PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
      Mesh<2>>::BoundaryIntegralReconstructionMatrixBuilder(const MeshType&,
                                                            const size_t,
                                                            const SmallArray<const Rd>&,
                                                            const SmallArray<const Rd>&,
                                                            const CellToCellStencilArray&);