Select Git revision
PolynomialReconstruction.cpp
-
Stéphane Del Pino authoredStéphane Del Pino authored
PolynomialReconstruction.cpp 39.95 KiB
#include <scheme/PolynomialReconstruction.hpp>
#include <algebra/Givens.hpp>
#include <algebra/ShrinkMatrixView.hpp>
#include <algebra/ShrinkVectorView.hpp>
#include <algebra/SmallMatrix.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/MeshFlatFaceBoundary.hpp>
#include <mesh/NamedBoundaryDescriptor.hpp>
#include <mesh/StencilDescriptor.hpp>
#include <mesh/StencilManager.hpp>
#include <scheme/DiscreteFunctionDPkVariant.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
#include <scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp>
#include <scheme/reconstruction_utils/CellCenterReconstructionMatrixBuilder.hpp>
#include <scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp>
#include <scheme/reconstruction_utils/MutableDiscreteFunctionDPkVariant.hpp>
#warning put in a file in geometry
template <size_t Dimension>
PUGS_INLINE auto
symmetrize_vector(const TinyVector<Dimension>& normal, const TinyVector<Dimension>& u)
{
return u - 2 * dot(u, normal) * normal;
}
template <size_t Dimension>
PUGS_INLINE auto
symmetrize_matrix(const TinyVector<Dimension>& normal, const TinyMatrix<Dimension>& A)
{
const TinyMatrix S = TinyMatrix<Dimension>{identity} - 2 * tensorProduct(normal, normal);
return S * A * S;
}
template <size_t Dimension>
PUGS_INLINE auto
symmetrize_coordinates(const TinyVector<Dimension>& origin,
const TinyVector<Dimension>& normal,
const TinyVector<Dimension>& u)
{
return u - 2 * dot(u - origin, normal) * normal;
}
size_t
PolynomialReconstruction::_getNumberOfColumns(
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
{
size_t number_of_columns = 0;
for (auto discrete_function_variant_p : discrete_function_variant_list) {
number_of_columns += std::visit(
[](auto&& discrete_function) -> size_t {
using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
using data_type = std::decay_t<typename DiscreteFunctionT::data_type>;
if constexpr (std::is_arithmetic_v<data_type>) {
return 1;
} else if constexpr (is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>) {
return data_type::Dimension;
} else {
// LCOV_EXCL_START
throw UnexpectedError("unexpected data type " + demangle<data_type>());
// LCOV_EXCL_STOP
}
} else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
using data_type = std::decay_t<typename DiscreteFunctionT::data_type>;
if constexpr (std::is_arithmetic_v<data_type>) {
return discrete_function.size();
} else if constexpr (is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>) {
return discrete_function.size() * data_type::Dimension;
} else {
// LCOV_EXCL_START
throw UnexpectedError("unexpected data type " + demangle<data_type>());
// LCOV_EXCL_STOP
}
} else {
// LCOV_EXCL_START
throw UnexpectedError("unexpected discrete function type");
// LCOV_EXCL_STOP
}
},
discrete_function_variant_p->discreteFunction());
}
return number_of_columns;
}
template <MeshConcept MeshType>
std::vector<PolynomialReconstruction::MutableDiscreteFunctionDPkVariant>
PolynomialReconstruction::_createMutableDiscreteFunctionDPKVariantList(
const std::shared_ptr<const MeshType>& p_mesh,
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
{
std::vector<MutableDiscreteFunctionDPkVariant> mutable_discrete_function_dpk_variant_list;
for (size_t i_discrete_function_variant = 0; i_discrete_function_variant < discrete_function_variant_list.size();
++i_discrete_function_variant) {
auto discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
std::visit(
[&](auto&& discrete_function) {
using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
using DataType = std::remove_const_t<std::decay_t<typename DiscreteFunctionT::data_type>>;
mutable_discrete_function_dpk_variant_list.push_back(
DiscreteFunctionDPk<MeshType::Dimension, DataType>(p_mesh, m_descriptor.degree()));
} else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
using DataType = std::remove_const_t<std::decay_t<typename DiscreteFunctionT::data_type>>;
mutable_discrete_function_dpk_variant_list.push_back(
DiscreteFunctionDPkVector<MeshType::Dimension, DataType>(p_mesh, m_descriptor.degree(),
discrete_function.size()));
} else {
// LCOV_EXCL_START
throw UnexpectedError("unexpected discrete function type");
// LCOV_EXCL_STOP
}
},
discrete_function_variant->discreteFunction());
}
return mutable_discrete_function_dpk_variant_list;
}
template <MeshConcept MeshType>
void
PolynomialReconstruction::_checkDataAndSymmetriesCompatibility(
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
{
for (auto&& discrete_function_variant : discrete_function_variant_list) {
std::visit(
[&](auto&& discrete_function) {
using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
if constexpr (is_discrete_function_P0_v<DiscreteFunctionT> or
is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
if constexpr (is_tiny_vector_v<DataType>) {
if constexpr (DataType::Dimension != MeshType::Dimension) {
std::stringstream error_msg;
error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
<< " using a mesh of dimension " << MeshType::Dimension;
throw NormalError(error_msg.str());
}
} else if constexpr (is_tiny_matrix_v<DataType>) {
if constexpr ((DataType::NumberOfRows != MeshType::Dimension) or
(DataType::NumberOfColumns != MeshType::Dimension)) {
std::stringstream error_msg;
error_msg << "cannot symmetrize matrices of dimensions " << DataType::NumberOfRows << 'x'
<< DataType::NumberOfColumns << " using a mesh of dimension " << MeshType::Dimension;
throw NormalError(error_msg.str());
}
}
} else {
// LCOV_EXCL_START
throw UnexpectedError("invalid discrete function type");
// LCOV_EXCL_STOP
}
},
discrete_function_variant->discreteFunction());
}
}
template <MeshConcept MeshType>
[[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
PolynomialReconstruction::_build(
const std::shared_ptr<const MeshType>& p_mesh,
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
{
const MeshType& mesh = *p_mesh;
using Rd = TinyVector<MeshType::Dimension>;
if (m_descriptor.symmetryBoundaryDescriptorList().size() > 0) {
this->_checkDataAndSymmetriesCompatibility<MeshType>(discrete_function_variant_list);
}
const size_t number_of_columns = this->_getNumberOfColumns(discrete_function_variant_list);
const size_t basis_dimension =
DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(m_descriptor.degree());
const auto& stencil_array =
StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), m_descriptor.stencilDescriptor(),
m_descriptor.symmetryBoundaryDescriptorList());
auto xr = mesh.xr();
auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
auto cell_is_owned = mesh.connectivity().cellIsOwned();
auto cell_type = mesh.connectivity().cellType();
auto cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
auto full_stencil_size = [&](const CellId cell_id) {
auto stencil_cell_list = stencil_array[cell_id];
size_t stencil_size = stencil_cell_list.size();
for (size_t i = 0; i < m_descriptor.symmetryBoundaryDescriptorList().size(); ++i) {
auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i].stencilArray();
stencil_size += ghost_stencil[cell_id].size();
}
return stencil_size;
};
const size_t max_stencil_size = [&]() {
size_t max_size = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const size_t stencil_size = full_stencil_size(cell_id);
if (cell_is_owned[cell_id] and stencil_size > max_size) {
max_size = stencil_size;
}
}
return max_size;
}();
SmallArray<const Rd> symmetry_normal_list = [&] {
SmallArray<Rd> normal_list(m_descriptor.symmetryBoundaryDescriptorList().size());
size_t i_symmetry_boundary = 0;
for (auto p_boundary_descriptor : m_descriptor.symmetryBoundaryDescriptorList()) {
const IBoundaryDescriptor& boundary_descriptor = *p_boundary_descriptor;
auto symmetry_boundary = getMeshFlatFaceBoundary(mesh, boundary_descriptor);
normal_list[i_symmetry_boundary++] = symmetry_boundary.outgoingNormal();
}
return normal_list;
}();
SmallArray<const Rd> symmetry_origin_list = [&] {
SmallArray<Rd> origin_list(m_descriptor.symmetryBoundaryDescriptorList().size());
size_t i_symmetry_boundary = 0;
for (auto p_boundary_descriptor : m_descriptor.symmetryBoundaryDescriptorList()) {
const IBoundaryDescriptor& boundary_descriptor = *p_boundary_descriptor;
auto symmetry_boundary = getMeshFlatFaceBoundary(mesh, boundary_descriptor);
origin_list[i_symmetry_boundary++] = symmetry_boundary.origin();
}
return origin_list;
}();
Kokkos::Experimental::UniqueToken<Kokkos::DefaultExecutionSpace::execution_space,
Kokkos::Experimental::UniqueTokenScope::Global>
tokens;
auto mutable_discrete_function_dpk_variant_list =
this->_createMutableDiscreteFunctionDPKVariantList(p_mesh, discrete_function_variant_list);
SmallArray<SmallMatrix<double>> A_pool(Kokkos::DefaultExecutionSpace::concurrency());
SmallArray<SmallMatrix<double>> B_pool(Kokkos::DefaultExecutionSpace::concurrency());
SmallArray<SmallVector<double>> G_pool(Kokkos::DefaultExecutionSpace::concurrency());
SmallArray<SmallMatrix<double>> X_pool(Kokkos::DefaultExecutionSpace::concurrency());
SmallArray<SmallArray<double>> inv_Vj_wq_detJ_ek_pool(Kokkos::DefaultExecutionSpace::concurrency());
SmallArray<SmallArray<double>> mean_j_of_ejk_pool(Kokkos::DefaultExecutionSpace::concurrency());
SmallArray<SmallArray<double>> mean_i_of_ejk_pool(Kokkos::DefaultExecutionSpace::concurrency());
SmallArray<SmallArray<double>> inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool(Kokkos::DefaultExecutionSpace::concurrency());
for (size_t i = 0; i < A_pool.size(); ++i) {
A_pool[i] = SmallMatrix<double>(max_stencil_size, basis_dimension - 1);
B_pool[i] = SmallMatrix<double>(max_stencil_size, number_of_columns);
G_pool[i] = SmallVector<double>(basis_dimension - 1);
X_pool[i] = SmallMatrix<double>(basis_dimension - 1, number_of_columns);
inv_Vj_wq_detJ_ek_pool[i] = SmallArray<double>(basis_dimension);
mean_j_of_ejk_pool[i] = SmallArray<double>(basis_dimension - 1);
mean_i_of_ejk_pool[i] = SmallArray<double>(basis_dimension - 1);
inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool[i] = SmallArray<double>(basis_dimension);
}
std::shared_ptr p_cell_center_reconstruction_matrix_builder =
std::make_shared<CellCenterReconstructionMatrixBuilder<MeshType>>(symmetry_origin_list, symmetry_normal_list,
stencil_array, xj);
parallel_for(
mesh.numberOfCells(), PUGS_CLASS_LAMBDA(const CellId cell_j_id) {
if (cell_is_owned[cell_j_id]) {
const int32_t t = tokens.acquire();
auto stencil_cell_list = stencil_array[cell_j_id];
ShrinkMatrixView B(B_pool[t], full_stencil_size(cell_j_id));
size_t column_begin = 0;
for (size_t i_discrete_function_variant = 0;
i_discrete_function_variant < discrete_function_variant_list.size(); ++i_discrete_function_variant) {
const auto& discrete_function_variant = discrete_function_variant_list[i_discrete_function_variant];
std::visit(
[&](auto&& discrete_function) {
using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
const DataType& qj = discrete_function[cell_j_id];
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];
const DataType& qi_qj = discrete_function[cell_i_id] - qj;
if constexpr (std::is_arithmetic_v<DataType>) {
B(index, column_begin) = qi_qj;
} else if constexpr (is_tiny_vector_v<DataType>) {
for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
B(index, kB) = qi_qj[k];
}
} else if constexpr (is_tiny_matrix_v<DataType>) {
for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
const size_t kB = column_begin + p * DataType::NumberOfColumns;
for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
B(index, kB + q) = qi_qj(p, q);
}
}
}
}
for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
++i_symmetry) {
auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
auto ghost_cell_list = ghost_stencil[cell_j_id];
for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
const CellId cell_i_id = ghost_cell_list[i];
if constexpr (std::is_arithmetic_v<DataType>) {
const DataType& qi_qj = discrete_function[cell_i_id] - qj;
B(index, column_begin) = qi_qj;
} else if constexpr (is_tiny_vector_v<DataType>) {
if constexpr (DataType::Dimension == MeshType::Dimension) {
const Rd& normal = symmetry_normal_list[i_symmetry];
const DataType& qi = discrete_function[cell_i_id];
const DataType& qi_qj = symmetrize_vector(normal, qi) - qj;
for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
B(index, kB) = qi_qj[k];
}
} else {
// LCOV_EXCL_START
std::stringstream error_msg;
error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
<< " using a mesh of dimension " << MeshType::Dimension;
throw UnexpectedError(error_msg.str());
// LCOV_EXCL_STOP
}
} else if constexpr (is_tiny_matrix_v<DataType>) {
if constexpr ((DataType::NumberOfColumns == DataType::NumberOfRows) and
(DataType::NumberOfColumns == MeshType::Dimension)) {
const Rd& normal = symmetry_normal_list[i_symmetry];
const DataType& qi = discrete_function[cell_i_id];
const DataType& qi_qj = symmetrize_matrix(normal, qi) - qj;
for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
B(index, column_begin + p * DataType::NumberOfColumns + q) = qi_qj(p, q);
}
}
} else {
// LCOV_EXCL_START
std::stringstream error_msg;
error_msg << "cannot symmetrize matrices of dimensions " << DataType::NumberOfRows << 'x'
<< DataType::NumberOfColumns << " using a mesh of dimension " << MeshType::Dimension;
throw UnexpectedError(error_msg.str());
// LCOV_EXCL_STOP
}
}
}
}
if constexpr (std::is_arithmetic_v<DataType>) {
++column_begin;
} else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
column_begin += DataType::Dimension;
}
} else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
const auto qj_vector = discrete_function[cell_j_id];
if constexpr (std::is_arithmetic_v<DataType>) {
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];
for (size_t l = 0; l < qj_vector.size(); ++l) {
const DataType& qj = qj_vector[l];
const DataType& qi_qj = discrete_function[cell_i_id][l] - qj;
B(index, column_begin + l) = qi_qj;
}
}
for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
++i_symmetry) {
auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
auto ghost_cell_list = ghost_stencil[cell_j_id];
for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
const CellId cell_i_id = ghost_cell_list[i];
for (size_t l = 0; l < qj_vector.size(); ++l) {
const DataType& qj = qj_vector[l];
const DataType& qi_qj = discrete_function[cell_i_id][l] - qj;
B(index, column_begin + l) = qi_qj;
}
}
}
} else if constexpr (is_tiny_vector_v<DataType>) {
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];
for (size_t l = 0; l < qj_vector.size(); ++l) {
const DataType& qj = qj_vector[l];
const DataType& qi_qj = discrete_function[cell_i_id][l] - qj;
for (size_t kB = column_begin + l * DataType::Dimension, k = 0; k < DataType::Dimension;
++k, ++kB) {
B(index, kB) = qi_qj[k];
}
}
}
for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
++i_symmetry) {
if constexpr (DataType::Dimension == MeshType::Dimension) {
auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
auto ghost_cell_list = ghost_stencil[cell_j_id];
const Rd& normal = 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];
for (size_t l = 0; l < qj_vector.size(); ++l) {
const DataType& qj = qj_vector[l];
const DataType& qi = discrete_function[cell_i_id][l];
const DataType& qi_qj = symmetrize_vector(normal, qi) - qj;
for (size_t kB = column_begin + l * DataType::Dimension, k = 0; k < DataType::Dimension;
++k, ++kB) {
B(index, kB) = qi_qj[k];
}
}
}
} else {
// LCOV_EXCL_START
std::stringstream error_msg;
error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
<< " using a mesh of dimension " << MeshType::Dimension;
throw UnexpectedError(error_msg.str());
// LCOV_EXCL_STOP
}
}
} else if constexpr (is_tiny_matrix_v<DataType>) {
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];
for (size_t l = 0; l < qj_vector.size(); ++l) {
const DataType& qj = qj_vector[l];
const DataType& qi = discrete_function[cell_i_id][l];
const DataType& qi_qj = qi - qj;
for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
const size_t kB = column_begin + l * DataType::Dimension + p * DataType::NumberOfColumns;
for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
B(index, kB + q) = qi_qj(p, q);
}
}
}
}
for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
++i_symmetry) {
if constexpr ((DataType::NumberOfRows == MeshType::Dimension) and
(DataType::NumberOfColumns == MeshType::Dimension)) {
auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
auto ghost_cell_list = ghost_stencil[cell_j_id];
const Rd& normal = 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];
for (size_t l = 0; l < qj_vector.size(); ++l) {
const DataType& qj = qj_vector[l];
const DataType& qi = discrete_function[cell_i_id][l];
const DataType& qi_qj = symmetrize_matrix(normal, qi) - qj;
for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
const size_t kB = column_begin + l * DataType::Dimension + p * DataType::NumberOfColumns;
for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
B(index, kB + q) = qi_qj(p, q);
}
}
}
}
} else {
// LCOV_EXCL_START
std::stringstream error_msg;
error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
<< " using a mesh of dimension " << MeshType::Dimension;
throw UnexpectedError(error_msg.str());
// LCOV_EXCL_STOP
}
}
}
if constexpr (std::is_arithmetic_v<DataType>) {
column_begin += qj_vector.size();
} else if constexpr (is_tiny_vector_v<DataType> or is_tiny_matrix_v<DataType>) {
column_begin += qj_vector.size() * DataType::Dimension;
}
} else {
// LCOV_EXCL_START
throw UnexpectedError("invalid discrete function type");
// LCOV_EXCL_STOP
}
},
discrete_function_variant->discreteFunction());
}
ShrinkMatrixView A(A_pool[t], full_stencil_size(cell_j_id));
if ((m_descriptor.degree() == 1) and
(m_descriptor.integrationMethodType() == IntegrationMethodType::cell_center)) {
p_cell_center_reconstruction_matrix_builder->build(cell_j_id, A);
} else if ((m_descriptor.integrationMethodType() == IntegrationMethodType::element) or
(m_descriptor.integrationMethodType() == IntegrationMethodType::boundary)) {
#warning implement boundary based reconstruction for 1d and 3d
if ((m_descriptor.integrationMethodType() == IntegrationMethodType::boundary) and
(MeshType::Dimension == 2)) {
if constexpr (MeshType::Dimension == 2) {
#warning incorporate in reconstruction_matrix_builder
SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek = inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool[t];
SmallArray<double>& mean_j_of_ejk = mean_j_of_ejk_pool[t];
SmallArray<double>& mean_i_of_ejk = mean_i_of_ejk_pool[t];
BoundaryIntegralReconstructionMatrixBuilder
boundary_integral_reconstruction_matrix_builder(*p_mesh, m_descriptor.degree(),
inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_j_of_ejk,
mean_i_of_ejk, symmetry_origin_list,
symmetry_normal_list, stencil_array);
boundary_integral_reconstruction_matrix_builder.build(cell_j_id, A);
} else {
throw NotImplementedError("invalid mesh dimension");
}
} else {
#warning incorporate in reconstruction_matrix_builder
SmallArray<double>& inv_Vj_wq_detJ_ek = inv_Vj_wq_detJ_ek_pool[t];
SmallArray<double>& mean_j_of_ejk = mean_j_of_ejk_pool[t];
SmallArray<double>& mean_i_of_ejk = mean_i_of_ejk_pool[t];
ElementIntegralReconstructionMatrixBuilder
element_integral_reconstruction_matrix_builder(*p_mesh, m_descriptor.degree(), inv_Vj_wq_detJ_ek,
mean_j_of_ejk, mean_i_of_ejk, symmetry_origin_list,
symmetry_normal_list, stencil_array);
element_integral_reconstruction_matrix_builder.build(cell_j_id, A);
}
} else {
throw UnexpectedError("invalid integration strategy");
}
if (m_descriptor.rowWeighting()) {
// Add row weighting (give more importance to cells that are
// closer to j)
const Rd& Xj = xj[cell_j_id];
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];
const double weight = 1. / l2Norm(xj[cell_i_id] - Xj);
for (size_t l = 0; l < basis_dimension - 1; ++l) {
A(index, l) *= weight;
}
for (size_t l = 0; l < number_of_columns; ++l) {
B(index, l) *= weight;
}
}
for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
++i_symmetry) {
auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
auto ghost_cell_list = ghost_stencil[cell_j_id];
const Rd& origin = symmetry_origin_list[i_symmetry];
const Rd& normal = 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];
const double weight = 1. / l2Norm(symmetrize_coordinates(origin, normal, xj[cell_i_id]) - Xj);
for (size_t l = 0; l < basis_dimension - 1; ++l) {
A(index, l) *= weight;
}
for (size_t l = 0; l < number_of_columns; ++l) {
B(index, l) *= weight;
}
}
}
}
const SmallMatrix<double>& X = X_pool[t];
if (m_descriptor.preconditioning()) {
// Add column weighting preconditioning (increase the presition)
SmallVector<double>& G = G_pool[t];
for (size_t l = 0; l < A.numberOfColumns(); ++l) {
double g = 0;
for (size_t i = 0; i < A.numberOfRows(); ++i) {
const double Ail = A(i, l);
g += Ail * Ail;
}
G[l] = std::sqrt(g);
}
for (size_t l = 0; l < A.numberOfColumns(); ++l) {
const double Gl = G[l];
for (size_t i = 0; i < A.numberOfRows(); ++i) {
A(i, l) *= Gl;
}
}
Givens::solveCollectionInPlace(A, X, B);
for (size_t l = 0; l < X.numberOfRows(); ++l) {
const double Gl = G[l];
for (size_t i = 0; i < X.numberOfColumns(); ++i) {
X(l, i) *= Gl;
}
}
} else {
Givens::solveCollectionInPlace(A, X, B);
}
column_begin = 0;
for (size_t i_dpk_variant = 0; i_dpk_variant < mutable_discrete_function_dpk_variant_list.size();
++i_dpk_variant) {
const auto& dpk_variant = mutable_discrete_function_dpk_variant_list[i_dpk_variant];
const auto& discrete_function_variant = discrete_function_variant_list[i_dpk_variant];
std::visit(
[&](auto&& dpk_function, auto&& p0_function) {
using DPkFunctionT = std::decay_t<decltype(dpk_function)>;
using P0FunctionT = std::decay_t<decltype(p0_function)>;
using DataType = std::remove_const_t<std::decay_t<typename DPkFunctionT::data_type>>;
using P0DataType = std::remove_const_t<std::decay_t<typename P0FunctionT::data_type>>;
if constexpr (std::is_same_v<DataType, P0DataType>) {
if constexpr (is_discrete_function_P0_v<P0FunctionT>) {
if constexpr (is_discrete_function_dpk_scalar_v<DPkFunctionT>) {
auto dpk_j = dpk_function.coefficients(cell_j_id);
dpk_j[0] = p0_function[cell_j_id];
if constexpr (std::is_arithmetic_v<DataType>) {
if (m_descriptor.degree() > 1) {
auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
for (size_t i = 0; i < basis_dimension - 1; ++i) {
dpk_j[0] -= X(i, column_begin) * mean_j_of_ejk[i];
}
}
for (size_t i = 0; i < basis_dimension - 1; ++i) {
auto& dpk_j_ip1 = dpk_j[i + 1];
dpk_j_ip1 = X(i, column_begin);
}
++column_begin;
} else if constexpr (is_tiny_vector_v<DataType>) {
if (m_descriptor.degree() > 1) {
auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
for (size_t i = 0; i < basis_dimension - 1; ++i) {
auto& dpk_j_0 = dpk_j[0];
for (size_t k = 0; k < DataType::Dimension; ++k) {
dpk_j_0[k] -= X(i, column_begin + k) * mean_j_of_ejk[i];
}
}
}
for (size_t i = 0; i < basis_dimension - 1; ++i) {
auto& dpk_j_ip1 = dpk_j[i + 1];
for (size_t k = 0; k < DataType::Dimension; ++k) {
dpk_j_ip1[k] = X(i, column_begin + k);
}
}
column_begin += DataType::Dimension;
} else if constexpr (is_tiny_matrix_v<DataType>) {
if (m_descriptor.degree() > 1) {
auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
for (size_t i = 0; i < basis_dimension - 1; ++i) {
auto& dpk_j_0 = dpk_j[0];
for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
dpk_j_0(k, l) -=
X(i, column_begin + k * DataType::NumberOfColumns + l) * mean_j_of_ejk[i];
}
}
}
}
for (size_t i = 0; i < basis_dimension - 1; ++i) {
auto& dpk_j_ip1 = dpk_j[i + 1];
for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
dpk_j_ip1(k, l) = X(i, column_begin + k * DataType::NumberOfColumns + l);
}
}
}
column_begin += DataType::Dimension;
} else {
// LCOV_EXCL_START
throw UnexpectedError("unexpected data type");
// LCOV_EXCL_STOP
}
} else {
// LCOV_EXCL_START
throw UnexpectedError("unexpected discrete dpk function type");
// LCOV_EXCL_STOP
}
} else if constexpr (is_discrete_function_P0_vector_v<P0FunctionT>) {
if constexpr (is_discrete_function_dpk_vector_v<DPkFunctionT>) {
auto dpk_j = dpk_function.coefficients(cell_j_id);
auto cell_vector = p0_function[cell_j_id];
const size_t size = basis_dimension;
for (size_t l = 0; l < cell_vector.size(); ++l) {
const size_t component_begin = l * size;
dpk_j[component_begin] = cell_vector[l];
if constexpr (std::is_arithmetic_v<DataType>) {
if (m_descriptor.degree() > 1) {
auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
for (size_t i = 0; i < basis_dimension - 1; ++i) {
dpk_j[component_begin] -= X(i, column_begin) * mean_j_of_ejk[i];
}
}
for (size_t i = 0; i < basis_dimension - 1; ++i) {
auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
dpk_j_ip1 = X(i, column_begin);
}
++column_begin;
} else if constexpr (is_tiny_vector_v<DataType>) {
if (m_descriptor.degree() > 1) {
auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
for (size_t i = 0; i < basis_dimension - 1; ++i) {
auto& dpk_j_0 = dpk_j[component_begin];
for (size_t k = 0; k < DataType::Dimension; ++k) {
dpk_j_0[k] -= X(i, column_begin + k) * mean_j_of_ejk[i];
}
}
}
for (size_t i = 0; i < basis_dimension - 1; ++i) {
auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
for (size_t k = 0; k < DataType::Dimension; ++k) {
dpk_j_ip1[k] = X(i, column_begin + k);
}
}
column_begin += DataType::Dimension;
} else if constexpr (is_tiny_matrix_v<DataType>) {
if (m_descriptor.degree() > 1) {
auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
for (size_t i = 0; i < basis_dimension - 1; ++i) {
auto& dpk_j_0 = dpk_j[component_begin];
for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
dpk_j_0(p, q) -=
X(i, column_begin + p * DataType::NumberOfColumns + q) * mean_j_of_ejk[i];
}
}
}
}
for (size_t i = 0; i < basis_dimension - 1; ++i) {
auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
dpk_j_ip1(p, q) = X(i, column_begin + p * DataType::NumberOfColumns + q);
}
}
}
column_begin += DataType::Dimension;
} else {
// LCOV_EXCL_START
throw UnexpectedError("unexpected data type");
// LCOV_EXCL_STOP
}
}
} else {
// LCOV_EXCL_START
throw UnexpectedError("unexpected discrete dpk function type");
// LCOV_EXCL_STOP
}
} else {
// LCOV_EXCL_START
throw UnexpectedError("unexpected discrete function type");
// LCOV_EXCL_STOP
}
} else {
// LCOV_EXCL_START
throw UnexpectedError("incompatible data types");
// LCOV_EXCL_STOP
}
},
dpk_variant.mutableDiscreteFunctionDPk(), discrete_function_variant->discreteFunction());
}
tokens.release(t);
}
});
std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> discrete_function_dpk_variant_list;
for (auto discrete_function_dpk_variant_p : mutable_discrete_function_dpk_variant_list) {
std::visit(
[&](auto&& mutable_function_dpk) {
synchronize(mutable_function_dpk.cellArrays());
discrete_function_dpk_variant_list.push_back(
std::make_shared<DiscreteFunctionDPkVariant>(mutable_function_dpk));
},
discrete_function_dpk_variant_p.mutableDiscreteFunctionDPk());
}
return discrete_function_dpk_variant_list;
}
std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>>
PolynomialReconstruction::build(
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
{
if (not hasSameMesh(discrete_function_variant_list)) {
throw NormalError("cannot reconstruct functions living of different meshes simultaneously");
}
auto mesh_v = getCommonMesh(discrete_function_variant_list);
return std::visit([&](auto&& p_mesh) { return this->_build(p_mesh, discrete_function_variant_list); },
mesh_v->variant());
}