diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 0f3982c4f3b6816996906550821308757e251a18..3c257c2d9d5133fe82b3f13ecfdca5b1c3146915 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -185,9 +185,8 @@ PolynomialReconstruction::_build(
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 cell_is_owned = mesh.connectivity().cellIsOwned();
+ auto cell_type = mesh.connectivity().cellType();
auto full_stencil_size = [&](const CellId cell_id) {
auto stencil_cell_list = stencil_array[cell_id];
@@ -270,6 +269,39 @@ PolynomialReconstruction::_build(
std::make_shared<CellCenterReconstructionMatrixBuilder<MeshType>>(symmetry_origin_list, symmetry_normal_list,
stencil_array, xj);
+ SmallArray<std::shared_ptr<ElementIntegralReconstructionMatrixBuilder<MeshType>>>
+ element_integral_reconstruction_matrix_builder_pool{A_pool.size()};
+
+ for (size_t t = 0; t < element_integral_reconstruction_matrix_builder_pool.size(); ++t) {
+ 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];
+
+ element_integral_reconstruction_matrix_builder_pool[t] =
+ std::make_shared<ElementIntegralReconstructionMatrixBuilder<MeshType>>(*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);
+ }
+
+ SmallArray<std::shared_ptr<BoundaryIntegralReconstructionMatrixBuilder<MeshType>>>
+ boundary_integral_reconstruction_matrix_builder_pool{A_pool.size()};
+
+ if constexpr (MeshType::Dimension == 2) {
+ for (size_t t = 0; t < boundary_integral_reconstruction_matrix_builder_pool.size(); ++t) {
+ 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];
+
+ boundary_integral_reconstruction_matrix_builder_pool[t] =
+ std::make_shared<BoundaryIntegralReconstructionMatrixBuilder<MeshType>>(*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);
+ }
+ }
+
parallel_for(
mesh.numberOfCells(), PUGS_CLASS_LAMBDA(const CellId cell_j_id) {
if (cell_is_owned[cell_j_id]) {
@@ -518,35 +550,12 @@ PolynomialReconstruction::_build(
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);
-
+ boundary_integral_reconstruction_matrix_builder_pool[t]->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);
+ element_integral_reconstruction_matrix_builder_pool[t]->build(cell_j_id, A);
}
} else {
throw UnexpectedError("invalid integration strategy");
diff --git a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
index 8e1ce86607ff8ed3f913307bca652701e14a3265..efb92d0db81af6be75026cfb9e926254f382b0d8 100644
--- a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
+++ b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
@@ -30,7 +30,7 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
const size_t m_basis_dimension;
const size_t m_polynomial_degree;
- const SmallArray<double> m_inv_Vj_wq_detJ_ek;
+ const SmallArray<double> m_wq_detJ_ek;
SmallArray<double> m_mean_j_of_ejk;
SmallArray<double> m_mean_i_of_ejk;
@@ -45,13 +45,21 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
const CellValue<const Rd> m_xj;
const NodeValue<const Rd> m_xr;
+ // 2D
+ SmallArray<const size_t> m_y_row_index;
+
+ // 3D
+ SmallArray<const size_t> m_yz_row_index;
+ SmallArray<const size_t> m_z_triangle_index;
+ SmallArray<const size_t> m_yz_row_size;
+
template <typename ConformTransformationT>
void
_computeEjkMean(const QuadratureFormula<1>& quadrature,
const ConformTransformationT& T,
const TinyVector<1>& Xj,
const double Vi,
- SmallArray<double>& mean_of_ejk)
+ SmallArray<double>& mean_of_ejk) noexcept(NO_ASSERT)
{
mean_of_ejk.fill(0);
@@ -65,17 +73,21 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
const double x_xj = X_Xj[0];
{
- size_t k = 0;
- m_inv_Vj_wq_detJ_ek[k++] = wq * detT / Vi;
- for (; k <= m_polynomial_degree; ++k) {
- m_inv_Vj_wq_detJ_ek[k] = x_xj * m_inv_Vj_wq_detJ_ek[k - 1];
+ m_wq_detJ_ek[0] = wq * detT;
+ for (size_t k = 1; k <= m_polynomial_degree; ++k) {
+ m_wq_detJ_ek[k] = x_xj * m_wq_detJ_ek[k - 1];
}
}
for (size_t k = 1; k < m_basis_dimension; ++k) {
- mean_of_ejk[k - 1] += m_inv_Vj_wq_detJ_ek[k];
+ mean_of_ejk[k - 1] += m_wq_detJ_ek[k];
}
}
+
+ const double inv_Vi = 1. / Vi;
+ for (size_t k = 0; k < mean_of_ejk.size(); ++k) {
+ mean_of_ejk[k] *= inv_Vi;
+ }
}
template <typename ConformTransformationT>
@@ -84,7 +96,7 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
const ConformTransformationT& T,
const TinyVector<2>& Xj,
const double Vi,
- SmallArray<double>& mean_of_ejk)
+ SmallArray<double>& mean_of_ejk) noexcept(NO_ASSERT)
{
mean_of_ejk.fill(0);
@@ -105,25 +117,29 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
const double y_yj = X_Xj[1];
{
- size_t k = 0;
- m_inv_Vj_wq_detJ_ek[k++] = wq * detT / Vi;
+ size_t k = 0;
+ m_wq_detJ_ek[k++] = wq * detT;
for (; k <= m_polynomial_degree; ++k) {
- m_inv_Vj_wq_detJ_ek[k] = x_xj * m_inv_Vj_wq_detJ_ek[k - 1];
+ m_wq_detJ_ek[k] = x_xj * m_wq_detJ_ek[k - 1];
}
for (size_t i_y = 1; i_y <= m_polynomial_degree; ++i_y) {
-#warning store y_row index and size
- 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_wq_detJ_ek[k++] = y_yj * m_inv_Vj_wq_detJ_ek[begin_i_y_1 + l];
+ const size_t begin_i_y_1 = m_y_row_index[i_y - 1];
+ for (size_t l = 0; l <= m_polynomial_degree - i_y; ++l, ++k) {
+ m_wq_detJ_ek[k] = y_yj * m_wq_detJ_ek[begin_i_y_1 + l];
}
}
}
for (size_t k = 1; k < m_basis_dimension; ++k) {
- mean_of_ejk[k - 1] += m_inv_Vj_wq_detJ_ek[k];
+ mean_of_ejk[k - 1] += m_wq_detJ_ek[k];
}
}
+
+ const double inv_Vi = 1. / Vi;
+ for (size_t k = 0; k < mean_of_ejk.size(); ++k) {
+ mean_of_ejk[k] *= inv_Vi;
+ }
}
template <typename ConformTransformationT>
@@ -132,31 +148,8 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
const ConformTransformationT& T,
const TinyVector<3>& Xj,
const double Vi,
- SmallArray<double>& mean_of_ejk)
+ SmallArray<double>& mean_of_ejk) noexcept(NO_ASSERT)
{
-#warning Compute this once for all and rework design
- SmallArray<size_t> m_yz_row_index((m_polynomial_degree + 2) * (m_polynomial_degree + 1) / 2 + 1);
- SmallArray<size_t> m_z_triangle_index(m_polynomial_degree + 1);
-
- {
- size_t i_z = 0;
- size_t i_yz = 0;
-
- m_yz_row_index[i_yz++] = 0;
- for (ssize_t n = m_polynomial_degree + 1; n >= 1; --n) {
- m_z_triangle_index[i_z++] = i_yz - 1;
- for (ssize_t i = n; i >= 1; --i) {
- m_yz_row_index[i_yz] = m_yz_row_index[i_yz - 1] + i;
- ++i_yz;
- }
- }
- }
-
- SmallArray<size_t> m_yz_row_size{m_yz_row_index.size() - 1};
- for (size_t i = 0; i < m_yz_row_size.size(); ++i) {
- m_yz_row_size[i] = m_yz_row_index[i + 1] - m_yz_row_index[i];
- }
-
mean_of_ejk.fill(0);
for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
@@ -177,17 +170,17 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
const double z_zj = X_Xj[2];
{
- size_t k = 0;
- m_inv_Vj_wq_detJ_ek[k++] = wq * detT / Vi;
+ size_t k = 0;
+ m_wq_detJ_ek[k++] = wq * detT;
for (; k <= m_polynomial_degree; ++k) {
- m_inv_Vj_wq_detJ_ek[k] = x_xj * m_inv_Vj_wq_detJ_ek[k - 1];
+ m_wq_detJ_ek[k] = x_xj * m_wq_detJ_ek[k - 1];
}
for (size_t i_y = 1; i_y <= m_polynomial_degree; ++i_y) {
const size_t begin_i_y_1 = m_yz_row_index[i_y - 1];
const size_t nb_monoms = m_yz_row_size[i_y];
- for (size_t l = 0; l < nb_monoms; ++l) {
- m_inv_Vj_wq_detJ_ek[k++] = y_yj * m_inv_Vj_wq_detJ_ek[begin_i_y_1 + l];
+ for (size_t l = 0; l < nb_monoms; ++l, ++k) {
+ m_wq_detJ_ek[k] = y_yj * m_wq_detJ_ek[begin_i_y_1 + l];
}
}
@@ -198,17 +191,22 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
for (size_t i_y = 0; i_y < nb_y; ++i_y) {
const size_t begin_i_yz_1 = m_yz_row_index[index_z_1 + i_y];
const size_t nb_monoms = m_yz_row_size[index_z + i_y];
- for (size_t l = 0; l < nb_monoms; ++l) {
- m_inv_Vj_wq_detJ_ek[k++] = z_zj * m_inv_Vj_wq_detJ_ek[begin_i_yz_1 + l];
+ for (size_t l = 0; l < nb_monoms; ++l, ++k) {
+ m_wq_detJ_ek[k] = z_zj * m_wq_detJ_ek[begin_i_yz_1 + l];
}
}
}
}
for (size_t k = 1; k < m_basis_dimension; ++k) {
- mean_of_ejk[k - 1] += m_inv_Vj_wq_detJ_ek[k];
+ mean_of_ejk[k - 1] += m_wq_detJ_ek[k];
}
}
+
+ const double inv_Vi = 1. / Vi;
+ for (size_t k = 0; k < mean_of_ejk.size(); ++k) {
+ mean_of_ejk[k] *= inv_Vi;
+ }
}
void
@@ -514,7 +512,7 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
polynomial_degree)},
m_polynomial_degree{polynomial_degree},
- m_inv_Vj_wq_detJ_ek{inv_Vj_wq_detJ_ek},
+ m_wq_detJ_ek{inv_Vj_wq_detJ_ek},
m_mean_j_of_ejk{mean_j_of_ejk},
m_mean_i_of_ejk{mean_i_of_ejk},
@@ -526,7 +524,48 @@ class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
m_Vj{MeshDataManager::instance().getMeshData(mesh).Vj()},
m_xj{MeshDataManager::instance().getMeshData(mesh).xj()},
m_xr{mesh.xr()}
- {}
+ {
+ if constexpr (MeshType::Dimension == 2) {
+ SmallArray<size_t> y_row_index(m_polynomial_degree + 2);
+
+ size_t i_y = 0;
+
+ y_row_index[i_y++] = 0;
+ for (ssize_t n = m_polynomial_degree + 1; n >= 1; --n) {
+ y_row_index[i_y] = y_row_index[i_y - 1] + n;
+ ++i_y;
+ }
+
+ m_y_row_index = y_row_index;
+
+ } else if constexpr (MeshType::Dimension == 3) {
+ SmallArray<size_t> yz_row_index((m_polynomial_degree + 2) * (m_polynomial_degree + 1) / 2 + 1);
+ SmallArray<size_t> z_triangle_index(m_polynomial_degree + 1);
+
+ {
+ size_t i_z = 0;
+ size_t i_yz = 0;
+
+ yz_row_index[i_yz++] = 0;
+ for (ssize_t n = m_polynomial_degree + 1; n >= 1; --n) {
+ z_triangle_index[i_z++] = i_yz - 1;
+ for (ssize_t i = n; i >= 1; --i) {
+ yz_row_index[i_yz] = yz_row_index[i_yz - 1] + i;
+ ++i_yz;
+ }
+ }
+ }
+
+ SmallArray<size_t> yz_row_size{yz_row_index.size() - 1};
+ for (size_t i = 0; i < yz_row_size.size(); ++i) {
+ yz_row_size[i] = yz_row_index[i + 1] - yz_row_index[i];
+ }
+
+ m_yz_row_index = yz_row_index;
+ m_z_triangle_index = z_triangle_index;
+ m_yz_row_size = yz_row_size;
+ }
+ }
~ElementIntegralReconstructionMatrixBuilder() = default;
};