From d3d3043eb20f1363e28da98bd85d569c5b9d573e Mon Sep 17 00:00:00 2001
From: MARMAJOU ISABELLE <id.counilh@wanadoo.fr>
Date: Wed, 26 Mar 2025 16:22:32 +0100
Subject: [PATCH] Implement variational solver for cubic meshes and plug
polygonal meshes
---
src/scheme/VariationalSolver.cpp | 345 +++++++++++++++----------------
src/scheme/VariationalSolver.hpp | 2 +-
2 files changed, 170 insertions(+), 177 deletions(-)
diff --git a/src/scheme/VariationalSolver.cpp b/src/scheme/VariationalSolver.cpp
index 8ab8f5492..cbd857ddd 100644
--- a/src/scheme/VariationalSolver.cpp
+++ b/src/scheme/VariationalSolver.cpp
@@ -166,7 +166,7 @@ variational_acoustic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c_
c.meshVariant()->variant());
}
-template <MeshConcept MeshTypeT>
+template <MeshConcept MeshTypeT, size_t mesh_edges_degree>
class VariationalSolverHandler::VariationalSolver final : public VariationalSolverHandler::IVariationalSolver
{
private:
@@ -254,7 +254,10 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
mesh.xr(),
mesh_node_boundary.nodeList());
- if (hasFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())) {
+ if constexpr (is_polygonal_mesh_v<MeshType>) {
+ bc_list.emplace_back(VelocityBoundaryCondition{mesh_node_boundary, node_value_list});
+ } else {
+ static_assert(is_polynomial_mesh_v<MeshType>);
MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
Table<const Rd> face_inner_node_value_list =
@@ -264,8 +267,6 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
bc_list.emplace_back(VelocityBoundaryCondition{mesh_node_boundary, node_value_list, //
mesh_face_boundary, face_inner_node_value_list});
- } else {
- bc_list.emplace_back(VelocityBoundaryCondition{mesh_node_boundary, node_value_list});
}
} else if (dirichlet_bc_descriptor.name() == "pressure") {
const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();
@@ -280,7 +281,7 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
Table<TinyVector<Dimension>> quadrature_points(mesh_face_boundary.faceList().size(), qf.numberOfPoints());
- LineTransformationAccessor<MeshType, 2> transformation{mesh};
+ LineTransformationAccessor<MeshType, mesh_edges_degree> transformation{mesh};
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
@@ -340,15 +341,15 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
this->_applyVelocityBC(bc_list, mesh, velocity_bc_treatment, A, b);
}
- void _forcebcSymmetryBC(const BoundaryConditionList& bc_list, const MeshType& mesh, Vector<double>& U) const;
+ void _forcebcSymmetryBC(const BoundaryConditionList& bc_list, Vector<double>& U) const;
- void _forcebcVelocityBC(const BoundaryConditionList& bc_list, const MeshType& mesh, Vector<double>& U) const;
+ void _forcebcVelocityBC(const BoundaryConditionList& bc_list, Vector<double>& U) const;
void
- _forcebcBoundaryConditions(const BoundaryConditionList& bc_list, const MeshType& mesh, Vector<double>& U) const
+ _forcebcBoundaryConditions(const BoundaryConditionList& bc_list, Vector<double>& U) const
{
- this->_forcebcSymmetryBC(bc_list, mesh, U);
- this->_forcebcVelocityBC(bc_list, mesh, U);
+ this->_forcebcSymmetryBC(bc_list, U);
+ this->_forcebcVelocityBC(bc_list, U);
}
std::tuple<NodeValue<const Rd>, FaceArray<const Rd>, CellValue<const Rd>, CellValue<const double>>
@@ -362,9 +363,9 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
std::optional<FunctionSymbolId> function_id) const
{
- std::shared_ptr mesh = getCommonMesh({rho_v, u_v, E_v})->get<MeshType>();
+ std::shared_ptr p_mesh = getCommonMesh({rho_v, u_v, E_v})->get<MeshType>();
- auto xr = mesh->xr();
+ auto xr = p_mesh->xr();
DiscreteScalarFunction rho = rho_v->get<DiscreteScalarFunction>();
DiscreteVectorFunction u = u_v->get<DiscreteVectorFunction>();
@@ -399,21 +400,21 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
return (gamma - 1) * rho_epsilon;
};
- auto node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
- auto face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+ auto node_to_face_matrix = p_mesh->connectivity().nodeToFaceMatrix();
+ auto face_to_node_matrix = p_mesh->connectivity().faceToNodeMatrix();
m_face_indices = [&]() {
- FaceArray<size_t> face_indices{mesh->connectivity(), mesh->degree() + 1};
+ FaceArray<size_t> face_indices{p_mesh->connectivity(), mesh_edges_degree + 1};
{
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ for (FaceId face_id = 0; face_id < p_mesh->numberOfFaces(); ++face_id) {
auto face_nodes = face_to_node_matrix[face_id];
- face_indices[face_id][0] = face_nodes[0];
- face_indices[face_id][mesh->degree()] = face_nodes[1];
+ face_indices[face_id][0] = face_nodes[0];
+ face_indices[face_id][mesh_edges_degree] = face_nodes[1];
}
- size_t cpt = mesh->numberOfNodes();
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- for (size_t i_face_node = 1; i_face_node < mesh->degree(); ++i_face_node) {
+ size_t cpt = p_mesh->numberOfNodes();
+ for (FaceId face_id = 0; face_id < p_mesh->numberOfFaces(); ++face_id) {
+ for (size_t i_face_node = 1; i_face_node < mesh_edges_degree; ++i_face_node) {
face_indices[face_id][i_face_node] = cpt++;
}
}
@@ -422,46 +423,48 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
}();
m_node_index = [&]() {
- NodeValue<size_t> node_index{mesh->connectivity()};
+ NodeValue<size_t> node_index{p_mesh->connectivity()};
{
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ for (FaceId face_id = 0; face_id < p_mesh->numberOfFaces(); ++face_id) {
auto face_nodes = face_to_node_matrix[face_id];
node_index[face_nodes[0]] = m_face_indices[face_id][0];
- node_index[face_nodes[1]] = m_face_indices[face_id][mesh->degree()];
+ node_index[face_nodes[1]] = m_face_indices[face_id][mesh_edges_degree];
}
}
return node_index;
}();
- Array<int> non_zero(Dimension * (mesh->numberOfNodes() + mesh->numberOfEdges() * (mesh->degree() - 1)));
+ Array<int> non_zero(Dimension * (p_mesh->numberOfNodes() + p_mesh->numberOfEdges() * (mesh_edges_degree - 1)));
parallel_for(
- mesh->numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+ p_mesh->numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
const size_t node_idx = m_node_index[node_id];
- const size_t node_non_zero = (node_to_face_matrix[node_id].size() * mesh->degree() + 1) * 2;
+ const size_t node_non_zero = (node_to_face_matrix[node_id].size() * mesh_edges_degree + 1) * 2;
non_zero[2 * node_idx] = node_non_zero;
non_zero[2 * node_idx + 1] = node_non_zero;
});
parallel_for(
- mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
- for (size_t i = 1; i < mesh->degree(); ++i) {
+ p_mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
+ for (size_t i = 1; i < mesh_edges_degree; ++i) {
const size_t face_node_idx = m_face_indices[face_id][i];
- const size_t face_node_non_zero = 2 * (mesh->degree() + 1);
+ const size_t face_node_non_zero = 2 * (mesh_edges_degree + 1);
non_zero[2 * face_node_idx] = face_node_non_zero;
non_zero[2 * face_node_idx + 1] = face_node_non_zero;
}
});
- CRSMatrixDescriptor A_descriptor(Dimension * (mesh->numberOfNodes() + mesh->numberOfEdges() * (mesh->degree() - 1)),
- Dimension * (mesh->numberOfNodes() + mesh->numberOfEdges() * (mesh->degree() - 1)),
+ CRSMatrixDescriptor A_descriptor(Dimension *
+ (p_mesh->numberOfNodes() + p_mesh->numberOfEdges() * (mesh_edges_degree - 1)),
+ Dimension *
+ (p_mesh->numberOfNodes() + p_mesh->numberOfEdges() * (mesh_edges_degree - 1)),
non_zero);
- auto face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
+ auto face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix();
const Rdxd I = identity;
@@ -478,19 +481,29 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
return 0.5 * x[0] * (x[0] + 1);
}};
+ const std::vector<std::function<double(R1)>> w_hat_set_P3 =
+ {[](R1 x) -> double { return -9. / 16 * (x[0] + 1. / 3) * (x[0] - 1. / 3) * (x[0] - 1); }, //
+ [](R1 x) -> double { return 27. / 16 * (x[0] + 1) * (x[0] - 1. / 3) * (x[0] - 1); }, //
+ [](R1 x) -> double { return -27. / 16 * (x[0] + 1) * (x[0] + 1. / 3) * (x[0] - 1); }, //
+ [](R1 x) -> double { return 9. / 16 * (x[0] + 1. / 3) * (x[0] - 1. / 3) * (x[0] + 1); }};
+
std::vector<std::function<double(R1)>> w_hat_set;
- if (mesh->degree() == 1) {
+ if constexpr (mesh_edges_degree == 1) {
w_hat_set = w_hat_set_P1;
- } else if (mesh->degree() == 2) {
+ } else if constexpr (mesh_edges_degree == 2) {
w_hat_set = w_hat_set_P2;
+ } else if constexpr (mesh_edges_degree == 3) {
+ w_hat_set = w_hat_set_P3;
} else {
- throw NotImplementedError("degree > 2");
+ throw NotImplementedError("degree > 3");
}
const QuadratureFormula<1> qf =
QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(m_quadrature_degree));
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ LineTransformationAccessor<MeshType, mesh_edges_degree> transformation{*p_mesh};
+
+ for (FaceId face_id = 0; face_id < p_mesh->numberOfFaces(); ++face_id) {
auto cell_list = face_to_cell_matrix[face_id];
double Z = 0;
@@ -501,18 +514,7 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
Za += a[cell_id];
}
- const Rd& x0 = mesh->xr()[face_to_node_matrix[face_id][0]];
- const Rd& x2 = mesh->xr()[face_to_node_matrix[face_id][1]];
-
- const Rd x1 = [&]() {
- if (mesh->degree() == 2) {
- return mesh->xl()[face_id][0];
- } else {
- return 0.5 * (x0 + x2);
- }
- }();
-
- const LineParabolicTransformation<Dimension> t(x0, x1, x2);
+ auto t = transformation.getTransformation(face_id);
for (size_t n0 = 0; n0 < w_hat_set.size(); ++n0) {
const size_t r = m_face_indices[face_id][n0];
@@ -541,25 +543,15 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
}
}
- const auto& face_local_numbers_in_their_cells = mesh->connectivity().faceLocalNumbersInTheirCells();
- const auto& face_cell_is_reversed = mesh->connectivity().cellFaceIsReversed();
+ const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
+ const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
- Vector<double> b(Dimension * (mesh->numberOfNodes() + mesh->numberOfEdges() * (mesh->degree() - 1)));
+ Vector<double> b(Dimension * (p_mesh->numberOfNodes() + p_mesh->numberOfEdges() * (mesh_edges_degree - 1)));
b.fill(0);
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- const Rd& x0 = mesh->xr()[face_to_node_matrix[face_id][0]];
- const Rd& x2 = mesh->xr()[face_to_node_matrix[face_id][1]];
+ for (FaceId face_id = 0; face_id < p_mesh->numberOfFaces(); ++face_id) {
+ auto t = transformation.getTransformation(face_id);
- const Rd x1 = [&]() {
- if (mesh->degree() == 2) {
- return mesh->xl()[face_id][0];
- } else {
- return 0.5 * (x0 + x2);
- }
- }();
-
- const LineParabolicTransformation<Dimension> t(x0, x1, x2);
auto cell_list = face_to_cell_matrix[face_id];
for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
@@ -596,9 +588,9 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
}
}
- const BoundaryConditionList bc_list = this->_getBCList(*mesh, bc_descriptor_list);
+ const BoundaryConditionList bc_list = this->_getBCList(*p_mesh, bc_descriptor_list);
- this->_applyBoundaryConditions(bc_list, *mesh, velocity_bc_treatment, A_descriptor, b);
+ this->_applyBoundaryConditions(bc_list, *p_mesh, velocity_bc_treatment, A_descriptor, b);
CRSMatrix A = A_descriptor.getCRSMatrix();
@@ -607,14 +599,14 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
LinearSolver solver;
solver.solveLocalSystem(A, U, b);
- this->_forcebcBoundaryConditions(bc_list, *mesh, U);
+ this->_forcebcBoundaryConditions(bc_list, U);
- const auto& cell_to_face_matrix = mesh->connectivity().cellToFaceMatrix();
+ const auto& cell_to_face_matrix = p_mesh->connectivity().cellToFaceMatrix();
- CellValue<Rd> momentum_fluxes{mesh->connectivity()};
- CellValue<double> energy_fluxes{mesh->connectivity()};
+ CellValue<Rd> momentum_fluxes{p_mesh->connectivity()};
+ CellValue<double> energy_fluxes{p_mesh->connectivity()};
- for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
Rd rho_u_j_fluxes = zero;
double rho_E_j_fluxes = 0;
@@ -622,18 +614,7 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
- const Rd& x0 = mesh->xr()[face_to_node_matrix[face_id][0]];
- const Rd& x2 = mesh->xr()[face_to_node_matrix[face_id][1]];
-
- const Rd x1 = [&]() {
- if (mesh->degree() == 2) {
- return mesh->xl()[face_id][0];
- } else {
- return 0.5 * (x0 + x2);
- }
- }();
-
- const LineParabolicTransformation<Dimension> t(x0, x1, x2);
+ auto t = transformation.getTransformation(face_id);
const double sign = face_cell_is_reversed(cell_id, i_face) ? -1 : 1;
@@ -680,33 +661,32 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
}
if (function_id.has_value()) {
- auto Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
+ auto Vj = MeshDataManager::instance().getMeshData(*p_mesh).Vj();
GaussLegendreQuadratureDescriptor quadrature(7);
- CellValue<double> S(mesh->connectivity());
+ CellValue<double> S(p_mesh->connectivity());
- IntegrateOnCells<double(Rd)>::integrateTo(function_id.value(), quadrature, *mesh, S);
+ IntegrateOnCells<double(Rd)>::integrateTo(function_id.value(), quadrature, *p_mesh, S);
- for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
energy_fluxes[cell_id] -= S[cell_id];
}
}
- NodeValue<Rd> ur(mesh->connectivity());
- FaceArray<Rd> ul(mesh->connectivity(), mesh->degree() - 1);
+ NodeValue<Rd> ur(p_mesh->connectivity());
+ FaceArray<Rd> ul(p_mesh->connectivity(), mesh_edges_degree - 1);
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ for (FaceId face_id = 0; face_id < p_mesh->numberOfFaces(); ++face_id) {
for (size_t i = 0; i < Dimension; ++i) {
ur[face_to_node_matrix[face_id][0]][i] = U[Dimension * m_face_indices[face_id][0] + i];
}
for (size_t i = 0; i < Dimension; ++i) {
- for (size_t il = 0; il < mesh->degree() - 1; ++il) {
+ for (size_t il = 0; il < mesh_edges_degree - 1; ++il) {
ul[face_id][il][i] = U[Dimension * m_face_indices[face_id][il + 1] + i];
}
}
for (size_t i = 0; i < Dimension; ++i) {
- ur[face_to_node_matrix[face_id][mesh->degree() - 1]][i] =
- U[Dimension * m_face_indices[face_id][mesh->degree()] + i];
+ ur[face_to_node_matrix[face_id][1]][i] = U[Dimension * m_face_indices[face_id][mesh_edges_degree] + i];
}
}
@@ -747,20 +727,26 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
new_E[cell_id] -= dt_over_Mj * energy_fluxes[cell_id];
}
+ std::shared_ptr<const MeshType> new_mesh;
+
NodeValue<Rd> new_xr = copy(mesh->xr());
- FaceArray<Rd> new_xl = copy(mesh->xl());
parallel_for(
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { new_xr[node_id] += dt * ur[node_id]; });
- parallel_for(
- mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
- for (size_t i = 0; i < mesh->degree() - 1; ++i) {
- new_xl[face_id][i] += dt * ul[face_id][i];
- }
- });
+ if constexpr (is_polygonal_mesh_v<MeshType>) {
+ new_mesh = std::make_shared<MeshType>(mesh->shared_connectivity(), new_xr);
+ } else if constexpr (is_polynomial_mesh_v<MeshType>) {
+ FaceArray<Rd> new_xl = copy(mesh->xl());
+ parallel_for(
+ mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+ for (size_t i = 0; i < mesh_edges_degree - 1; ++i) {
+ new_xl[face_id][i] += dt * ul[face_id][i];
+ }
+ });
- std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh->shared_connectivity(), new_xr, new_xl);
+ new_mesh = std::make_shared<MeshType>(mesh->shared_connectivity(), new_xr, new_xl);
+ }
CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
@@ -893,11 +879,12 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
~VariationalSolver() = default;
};
-template <MeshConcept MeshType>
+template <MeshConcept MeshType, size_t mesh_edges_degree>
void
-VariationalSolverHandler::VariationalSolver<MeshType>::_applyPressureBC(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- Vector<double>& b) const
+VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_applyPressureBC(
+ const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ Vector<double>& b) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
@@ -918,17 +905,24 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applyPressureBC(const Bo
[](R1 x) -> double { return -(x[0] - 1) * (x[0] + 1); },
[](R1 x) -> double { return 0.5 * x[0] * (x[0] + 1); }};
+ const std::vector<std::function<double(R1)>> w_hat_set_P3 =
+ {[](R1 x) -> double { return -9. / 16 * (x[0] + 1. / 3) * (x[0] - 1. / 3) * (x[0] - 1); }, //
+ [](R1 x) -> double { return 27. / 16 * (x[0] + 1) * (x[0] - 1. / 3) * (x[0] - 1); }, //
+ [](R1 x) -> double { return -27. / 16 * (x[0] + 1) * (x[0] + 1. / 3) * (x[0] - 1); }, //
+ [](R1 x) -> double { return 9. / 16 * (x[0] + 1. / 3) * (x[0] - 1. / 3) * (x[0] + 1); }};
+
std::vector<std::function<double(R1)>> w_hat_set;
- if (mesh.degree() == 1) {
+ if constexpr (mesh_edges_degree == 1) {
w_hat_set = w_hat_set_P1;
- } else if (mesh.degree() == 2) {
+ } else if constexpr (mesh_edges_degree == 2) {
w_hat_set = w_hat_set_P2;
+ } else if constexpr (mesh_edges_degree == 3) {
+ w_hat_set = w_hat_set_P3;
} else {
- throw NotImplementedError("degree > 2");
+ throw NotImplementedError("degree > 3");
}
const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
- const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
const Array<const FaceId>& face_list = bc.faceList();
const Table<const double>& p_ext = bc.valueList();
@@ -937,21 +931,13 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applyPressureBC(const Bo
const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
+ LineTransformationAccessor<MeshType, mesh_edges_degree> transformation{mesh};
+
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
FaceId face_id = face_list[i_face];
- const Rd& x0 = mesh.xr()[face_to_node_matrix[face_id][0]];
- const Rd& x2 = mesh.xr()[face_to_node_matrix[face_id][1]];
+ auto t = transformation.getTransformation(face_id);
- const Rd x1 = [&]() {
- if (mesh.degree() == 2) {
- return mesh.xl()[face_id][0];
- } else {
- return 0.5 * (x0 + x2);
- }
- }();
-
- const LineParabolicTransformation<Dimension> t(x0, x1, x2);
auto cell_list = face_to_cell_matrix[face_id];
for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
@@ -988,11 +974,11 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applyPressureBC(const Bo
}
}
-template <MeshConcept MeshType>
+template <MeshConcept MeshType, size_t mesh_edges_degree>
void
-VariationalSolverHandler::VariationalSolver<MeshType>::_forcebcSymmetryBC(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- Vector<double>& U) const
+VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_forcebcSymmetryBC(
+ const BoundaryConditionList& bc_list,
+ Vector<double>& U) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
@@ -1023,15 +1009,13 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_forcebcSymmetryBC(const
}
}
- const size_t degree = mesh.degree();
-
const Array<const FaceId>& face_list = bc.faceList();
// treat inner nodes of the faces
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
- for (size_t r_hat = 1; r_hat < degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_edges_degree; ++r_hat) {
const size_t r = m_face_indices[face_id][r_hat];
TinyVector<Dimension> ur;
@@ -1053,11 +1037,11 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_forcebcSymmetryBC(const
}
}
-template <MeshConcept MeshType>
+template <MeshConcept MeshType, size_t mesh_edges_degree>
void
-VariationalSolverHandler::VariationalSolver<MeshType>::_forcebcVelocityBC(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- Vector<double>& U) const
+VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_forcebcVelocityBC(
+ const BoundaryConditionList& bc_list,
+ Vector<double>& U) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
@@ -1077,8 +1061,6 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_forcebcVelocityBC(const
}
}
- const size_t degree = mesh.degree();
-
const Array<const FaceId>& face_list = bc.faceList();
const Table<const Rd>& face_inner_node_value_list = bc.faceInnerNodeValueList();
@@ -1086,7 +1068,7 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_forcebcVelocityBC(const
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
- for (size_t r_hat = 1; r_hat < degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_edges_degree; ++r_hat) {
const size_t r = m_face_indices[face_id][r_hat];
for (size_t i = 0; i < Dimension; ++i) {
@@ -1107,15 +1089,13 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_forcebcVelocityBC(const
}
}
- const size_t degree = mesh.degree();
-
const Array<const FaceId>& face_list = bc.faceList();
// treat inner nodes of the faces
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
- for (size_t r_hat = 1; r_hat < degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_edges_degree; ++r_hat) {
const size_t r = m_face_indices[face_id][r_hat];
for (size_t i = 0; i < Dimension; ++i) {
@@ -1129,12 +1109,13 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_forcebcVelocityBC(const
}
}
-template <MeshConcept MeshType>
+template <MeshConcept MeshType, size_t mesh_edges_degree>
void
-VariationalSolverHandler::VariationalSolver<MeshType>::_applySymmetryBC(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- CRSMatrixDescriptor<double>& A_descriptor,
- Vector<double>& b) const
+VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_applySymmetryBC(
+ const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ CRSMatrixDescriptor<double>& A_descriptor,
+ Vector<double>& b) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
@@ -1182,13 +1163,11 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applySymmetryBC(const Bo
}
}
- const size_t degree = mesh.degree();
-
const Array<const FaceId>& face_list = bc.faceList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
- for (size_t i_hat = 1; i_hat < degree; ++i_hat) {
+ for (size_t i_hat = 1; i_hat < mesh_edges_degree; ++i_hat) {
const size_t r = m_face_indices[face_id][i_hat];
TinyMatrix<Dimension> Arr;
@@ -1231,7 +1210,7 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applySymmetryBC(const Bo
const FaceId face_id = node_face_list[i_face];
auto face_node_list = face_to_node_matrix[face_id];
- for (size_t i_hat = 0; i_hat <= degree; ++i_hat) {
+ for (size_t i_hat = 0; i_hat <= mesh_edges_degree; ++i_hat) {
const size_t s = m_face_indices[face_id][i_hat];
if (s != r) {
@@ -1256,10 +1235,10 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applySymmetryBC(const Bo
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
- for (size_t i_hat = 1; i_hat < degree; ++i_hat) {
+ for (size_t i_hat = 1; i_hat < mesh_edges_degree; ++i_hat) {
const size_t r = m_face_indices[face_id][i_hat];
- for (size_t j_hat = 0; j_hat <= degree; ++j_hat) {
+ for (size_t j_hat = 0; j_hat <= mesh_edges_degree; ++j_hat) {
const size_t s = m_face_indices[face_id][j_hat];
if (s != r) {
@@ -1286,9 +1265,9 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applySymmetryBC(const Bo
}
}
-template <MeshConcept MeshType>
+template <MeshConcept MeshType, size_t mesh_edges_degree>
void
-VariationalSolverHandler::VariationalSolver<MeshType>::_applyVelocityBC(
+VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_applyVelocityBC(
const BoundaryConditionList& bc_list,
const MeshType& mesh,
const VariationalSolverHandler::VelocityBCTreatment& velocity_bc_treatment,
@@ -1307,7 +1286,6 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applyVelocityBC(
const Array<const NodeId>& node_list = bc.nodeList();
auto node_to_face_matrix = mesh.connectivity().nodeToFaceMatrix();
const Array<const FaceId>& face_list = bc.faceList();
- const size_t degree = mesh.degree();
const Array<const Rd>& node_value_list = bc.nodeValueList();
const Table<const Rd>& face_inner_node_value_list = bc.faceInnerNodeValueList();
@@ -1321,7 +1299,7 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applyVelocityBC(
for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) {
const FaceId face_id = node_face_list[i_face];
- for (size_t s_hat = 0; s_hat <= degree; ++s_hat) {
+ for (size_t s_hat = 0; s_hat <= mesh_edges_degree; ++s_hat) {
const size_t s = m_face_indices[face_id][s_hat];
if (r != s) {
@@ -1345,9 +1323,9 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applyVelocityBC(
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
- for (size_t r_hat = 1; r_hat < degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_edges_degree; ++r_hat) {
const size_t r = m_face_indices[face_id][r_hat];
- for (size_t s_hat = 0; s_hat <= degree; ++s_hat) {
+ for (size_t s_hat = 0; s_hat <= mesh_edges_degree; ++s_hat) {
const size_t s = m_face_indices[face_id][s_hat];
if (r != s) {
@@ -1377,7 +1355,7 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applyVelocityBC(
for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) {
const FaceId face_id = node_face_list[i_face];
- for (size_t s_hat = 0; s_hat <= degree; ++s_hat) {
+ for (size_t s_hat = 0; s_hat <= mesh_edges_degree; ++s_hat) {
const size_t s = m_face_indices[face_id][s_hat];
if (r == s) {
@@ -1403,9 +1381,9 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applyVelocityBC(
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
- for (size_t r_hat = 1; r_hat < degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_edges_degree; ++r_hat) {
const size_t r = m_face_indices[face_id][r_hat];
- for (size_t s_hat = 0; s_hat <= degree; ++s_hat) {
+ for (size_t s_hat = 0; s_hat <= mesh_edges_degree; ++s_hat) {
const size_t s = m_face_indices[face_id][s_hat];
if (r == s) {
@@ -1432,7 +1410,6 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applyVelocityBC(
const Array<const NodeId>& node_list = bc.nodeList();
auto node_to_face_matrix = mesh.connectivity().nodeToFaceMatrix();
const Array<const FaceId>& face_list = bc.faceList();
- const size_t degree = mesh.degree();
// Treats vertices of the faces
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
@@ -1443,7 +1420,7 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applyVelocityBC(
for (size_t i_face = 0; i_face < node_faces.size(); ++i_face) {
const FaceId face_id = node_faces[i_face];
- for (size_t s_hat = 0; s_hat <= degree; ++s_hat) {
+ for (size_t s_hat = 0; s_hat <= mesh_edges_degree; ++s_hat) {
const size_t s = m_face_indices[face_id][s_hat];
if (r == s) {
@@ -1469,9 +1446,9 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applyVelocityBC(
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
- for (size_t r_hat = 1; r_hat < degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_edges_degree; ++r_hat) {
const size_t r = m_face_indices[face_id][r_hat];
- for (size_t s_hat = 0; s_hat <= degree; ++s_hat) {
+ for (size_t s_hat = 0; s_hat <= mesh_edges_degree; ++s_hat) {
const size_t s = m_face_indices[face_id][s_hat];
if (r == s) {
@@ -1513,11 +1490,10 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applyVelocityBC(
A_descriptor(Dimension * r + 1, Dimension * r + 1) += 1.e30;
}
- const size_t degree = mesh.degree();
const Array<const FaceId>& face_list = bc.faceList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
- for (size_t r_hat = 1; r_hat < degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_edges_degree; ++r_hat) {
const size_t r = m_face_indices[face_id][r_hat];
A_descriptor(Dimension * r, Dimension * r) += 1.e30;
A_descriptor(Dimension * r + 1, Dimension * r + 1) += 1.e30;
@@ -1536,7 +1512,7 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applyVelocityBC(
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
- for (size_t r_hat = 1; r_hat < degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_edges_degree; ++r_hat) {
const size_t value_index = r_hat - 1; // nodal values on the face are not stored in this array
const size_t r = m_face_indices[face_id][r_hat];
@@ -1554,11 +1530,10 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applyVelocityBC(
A_descriptor(Dimension * r + 1, Dimension * r + 1) += 1.e30;
}
- const size_t degree = mesh.degree();
const Array<const FaceId>& face_list = bc.faceList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
- for (size_t r_hat = 1; r_hat < degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_edges_degree; ++r_hat) {
const size_t r = m_face_indices[face_id][r_hat];
A_descriptor(Dimension * r, Dimension * r) += 1.e30;
A_descriptor(Dimension * r + 1, Dimension * r + 1) += 1.e30;
@@ -1574,8 +1549,8 @@ VariationalSolverHandler::VariationalSolver<MeshType>::_applyVelocityBC(
}
}
-template <MeshConcept MeshType>
-class VariationalSolverHandler::VariationalSolver<MeshType>::FixedBoundaryCondition
+template <MeshConcept MeshType, size_t mesh_edges_degree>
+class VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::FixedBoundaryCondition
{
private:
const MeshNodeBoundary m_mesh_node_boundary;
@@ -1605,8 +1580,8 @@ class VariationalSolverHandler::VariationalSolver<MeshType>::FixedBoundaryCondit
~FixedBoundaryCondition() = default;
};
-template <MeshConcept MeshType>
-class VariationalSolverHandler::VariationalSolver<MeshType>::PressureBoundaryCondition
+template <MeshConcept MeshType, size_t mesh_edges_degree>
+class VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::PressureBoundaryCondition
{
private:
const MeshFaceBoundary m_mesh_face_boundary;
@@ -1632,8 +1607,8 @@ class VariationalSolverHandler::VariationalSolver<MeshType>::PressureBoundaryCon
~PressureBoundaryCondition() = default;
};
-template <MeshConcept MeshType>
-class VariationalSolverHandler::VariationalSolver<MeshType>::VelocityBoundaryCondition
+template <MeshConcept MeshType, size_t mesh_edges_degree>
+class VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::VelocityBoundaryCondition
{
private:
const MeshNodeBoundary m_mesh_node_boundary;
@@ -1688,8 +1663,8 @@ class VariationalSolverHandler::VariationalSolver<MeshType>::VelocityBoundaryCon
~VelocityBoundaryCondition() = default;
};
-template <MeshConcept MeshType>
-class VariationalSolverHandler::VariationalSolver<MeshType>::SymmetryBoundaryCondition
+template <MeshConcept MeshType, size_t mesh_edges_degree>
+class VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::SymmetryBoundaryCondition
{
public:
using Rd = TinyVector<Dimension, double>;
@@ -1738,8 +1713,26 @@ VariationalSolverHandler::VariationalSolverHandler(const std::shared_ptr<const M
std::visit(
[&](auto&& mesh) {
using MeshType = mesh_type_t<decltype(mesh)>;
- if constexpr ((is_polynomial_mesh_v<MeshType>)and(MeshType::Dimension == 2)) {
- m_acoustic_solver = std::make_unique<VariationalSolver<MeshType>>();
+ if constexpr ((is_polygonal_mesh_v<MeshType>)and(MeshType::Dimension == 2)) {
+ m_acoustic_solver = std::make_unique<VariationalSolver<MeshType, 1>>();
+ } else if constexpr ((is_polynomial_mesh_v<MeshType>)and(MeshType::Dimension == 2)) {
+ switch (mesh->degree()) {
+ case 1: {
+ m_acoustic_solver = std::make_unique<VariationalSolver<MeshType, 1>>();
+ break;
+ }
+ case 2: {
+ m_acoustic_solver = std::make_unique<VariationalSolver<MeshType, 2>>();
+ break;
+ }
+ case 3: {
+ m_acoustic_solver = std::make_unique<VariationalSolver<MeshType, 3>>();
+ break;
+ }
+ default: {
+ throw NotImplementedError("polynomial mesh of degree > 3 are not supported");
+ }
+ }
} else {
throw NormalError("unexpected mesh type");
}
diff --git a/src/scheme/VariationalSolver.hpp b/src/scheme/VariationalSolver.hpp
index 85a2b158b..64bbb6b23 100644
--- a/src/scheme/VariationalSolver.hpp
+++ b/src/scheme/VariationalSolver.hpp
@@ -53,7 +53,7 @@ class VariationalSolverHandler
virtual ~IVariationalSolver() = default;
};
- template <MeshConcept MeshTypeT>
+ template <MeshConcept MeshTypeT, size_t mesh_edges_degree>
class VariationalSolver;
std::unique_ptr<IVariationalSolver> m_acoustic_solver;
--
GitLab