Select Git revision
PugsUtils.cpp
-
Stéphane Del Pino authored
When `--no-preamble` is used, pugs do not print execution info. This at least is useful for documentation generation.
Stéphane Del Pino authoredWhen `--no-preamble` is used, pugs do not print execution info. This at least is useful for documentation generation.
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&);