Select Git revision
ParallelChecker.hpp
-
Stéphane Del Pino authored
This tool checks if chosen data are synchronized and have same values comparing two runs By now only ItemValue are checked
Stéphane Del Pino authoredThis tool checks if chosen data are synchronized and have same values comparing two runs By now only ItemValue are checked
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());
}