diff --git a/src/scheme/VariationalSolver.cpp b/src/scheme/VariationalSolver.cpp
index cbd857dddb3ba467ce4675194865fe70ef4f64f0..34f946900ae14ceafb9fdc4720cbd492c3157d98 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, size_t mesh_edges_degree>
+template <MeshConcept MeshTypeT, size_t mesh_degree>
class VariationalSolverHandler::VariationalSolver final : public VariationalSolverHandler::IVariationalSolver
{
private:
@@ -281,7 +281,7 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
Table<TinyVector<Dimension>> quadrature_points(mesh_face_boundary.faceList().size(), qf.numberOfPoints());
- LineTransformationAccessor<MeshType, mesh_edges_degree> transformation{mesh};
+ LineTransformationAccessor<MeshType, mesh_degree> transformation{mesh};
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
@@ -353,7 +353,7 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
}
std::tuple<NodeValue<const Rd>, FaceArray<const Rd>, CellValue<const Rd>, CellValue<const double>>
- _computeFluxes(const size_t& degree,
+ _computeFluxes(const size_t& polynomial_degree,
const VelocityBCTreatment& velocity_bc_treatment,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
@@ -379,7 +379,7 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
}
}
- PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::boundary, degree,
+ PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::boundary, polynomial_degree,
symmetry_boundary_descriptor_list);
auto reconstructions = PolynomialReconstruction{reconstruction_descriptor}.build(rho, rho * u, rho * E);
@@ -404,17 +404,17 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
auto face_to_node_matrix = p_mesh->connectivity().faceToNodeMatrix();
m_face_indices = [&]() {
- FaceArray<size_t> face_indices{p_mesh->connectivity(), mesh_edges_degree + 1};
+ FaceArray<size_t> face_indices{p_mesh->connectivity(), mesh_degree + 1};
{
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_edges_degree] = face_nodes[1];
+ face_indices[face_id][0] = face_nodes[0];
+ face_indices[face_id][mesh_degree] = face_nodes[1];
}
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) {
+ for (size_t i_face_node = 1; i_face_node < mesh_degree; ++i_face_node) {
face_indices[face_id][i_face_node] = cpt++;
}
}
@@ -429,18 +429,18 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
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_edges_degree];
+ node_index[face_nodes[1]] = m_face_indices[face_id][mesh_degree];
}
}
return node_index;
}();
- Array<int> non_zero(Dimension * (p_mesh->numberOfNodes() + p_mesh->numberOfEdges() * (mesh_edges_degree - 1)));
+ Array<int> non_zero(Dimension * (p_mesh->numberOfNodes() + p_mesh->numberOfEdges() * (mesh_degree - 1)));
parallel_for(
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_edges_degree + 1) * 2;
+ const size_t node_non_zero = (node_to_face_matrix[node_id].size() * mesh_degree + 1) * 2;
non_zero[2 * node_idx] = node_non_zero;
non_zero[2 * node_idx + 1] = node_non_zero;
@@ -448,10 +448,10 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
parallel_for(
p_mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
- for (size_t i = 1; i < mesh_edges_degree; ++i) {
+ for (size_t i = 1; i < mesh_degree; ++i) {
const size_t face_node_idx = m_face_indices[face_id][i];
- const size_t face_node_non_zero = 2 * (mesh_edges_degree + 1);
+ const size_t face_node_non_zero = 2 * (mesh_degree + 1);
non_zero[2 * face_node_idx] = face_node_non_zero;
non_zero[2 * face_node_idx + 1] = face_node_non_zero;
@@ -459,9 +459,9 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
});
CRSMatrixDescriptor A_descriptor(Dimension *
- (p_mesh->numberOfNodes() + p_mesh->numberOfEdges() * (mesh_edges_degree - 1)),
+ (p_mesh->numberOfNodes() + p_mesh->numberOfEdges() * (mesh_degree - 1)),
Dimension *
- (p_mesh->numberOfNodes() + p_mesh->numberOfEdges() * (mesh_edges_degree - 1)),
+ (p_mesh->numberOfNodes() + p_mesh->numberOfEdges() * (mesh_degree - 1)),
non_zero);
auto face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix();
@@ -488,11 +488,11 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
[](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 constexpr (mesh_edges_degree == 1) {
+ if constexpr (mesh_degree == 1) {
w_hat_set = w_hat_set_P1;
- } else if constexpr (mesh_edges_degree == 2) {
+ } else if constexpr (mesh_degree == 2) {
w_hat_set = w_hat_set_P2;
- } else if constexpr (mesh_edges_degree == 3) {
+ } else if constexpr (mesh_degree == 3) {
w_hat_set = w_hat_set_P3;
} else {
throw NotImplementedError("degree > 3");
@@ -501,7 +501,7 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
const QuadratureFormula<1> qf =
QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(m_quadrature_degree));
- LineTransformationAccessor<MeshType, mesh_edges_degree> transformation{*p_mesh};
+ LineTransformationAccessor<MeshType, mesh_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];
@@ -546,7 +546,7 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
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 * (p_mesh->numberOfNodes() + p_mesh->numberOfEdges() * (mesh_edges_degree - 1)));
+ Vector<double> b(Dimension * (p_mesh->numberOfNodes() + p_mesh->numberOfEdges() * (mesh_degree - 1)));
b.fill(0);
for (FaceId face_id = 0; face_id < p_mesh->numberOfFaces(); ++face_id) {
@@ -674,19 +674,19 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
}
NodeValue<Rd> ur(p_mesh->connectivity());
- FaceArray<Rd> ul(p_mesh->connectivity(), mesh_edges_degree - 1);
+ FaceArray<Rd> ul(p_mesh->connectivity(), mesh_degree - 1);
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_edges_degree - 1; ++il) {
+ for (size_t il = 0; il < mesh_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][1]][i] = U[Dimension * m_face_indices[face_id][mesh_edges_degree] + i];
+ ur[face_to_node_matrix[face_id][1]][i] = U[Dimension * m_face_indices[face_id][mesh_degree] + i];
}
}
@@ -740,7 +740,7 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
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) {
+ for (size_t i = 0; i < mesh_degree - 1; ++i) {
new_xl[face_id][i] += dt * ul[face_id][i];
}
});
@@ -764,7 +764,7 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
- apply(const size_t& degree,
+ apply(const size_t& polynomial_degree,
const double& dt,
const VelocityBCTreatment& velocity_bc_treatment,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
@@ -781,22 +781,78 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
const CellValue<const double> rhoc =
(rho_v->get<DiscreteFunctionP0<const double>>() * c_v->get<DiscreteFunctionP0<const double>>()).cellValues();
-#if 0
+#if 1
+ // k1
+ auto [ur1, ul1, momentum_fluxes1, energy_fluxes1] =
+ _computeFluxes(polynomial_degree, velocity_bc_treatment, rho_v, u_v, E_v, rhoc, a, bc_descriptor_list,
+ function_id);
+
+ // y1
+ auto [mesh1, rho1, u1, E1] = _applyFluxes(0.5 * dt, rho_v, u_v, E_v, ur1, ul1, momentum_fluxes1, energy_fluxes1);
+
+ // k2
+ auto [ur2, ul2, momentum_fluxes2, energy_fluxes2] =
+ _computeFluxes(polynomial_degree, velocity_bc_treatment, rho1, u1, E1, rhoc, a, bc_descriptor_list, function_id);
+
+ // y2
+ auto [mesh2, rho2, u2, E2] = _applyFluxes(0.5 * dt, rho_v, u_v, E_v, ur2, ul2, momentum_fluxes2, energy_fluxes2);
+
+ // k3
+ auto [ur3, ul3, momentum_fluxes3, energy_fluxes3] =
+ _computeFluxes(polynomial_degree, velocity_bc_treatment, rho2, u2, E2, rhoc, a, bc_descriptor_list, function_id);
+
+ // y3
+ auto [mesh3, rho3, u3, E3] = _applyFluxes(dt, rho_v, u_v, E_v, ur3, ul3, momentum_fluxes3, energy_fluxes3);
+
+ // k4
+ auto [ur4, ul4, momentum_fluxes4, energy_fluxes4] =
+ _computeFluxes(polynomial_degree, velocity_bc_treatment, rho3, u3, E3, rhoc, a, bc_descriptor_list, function_id);
+
+ NodeValue<Rd> ur{mesh->connectivity()};
+ FaceArray<Rd> ul{mesh->connectivity(), ul1.sizeOfArrays()};
+ CellValue<Rd> momentum_fluxes{mesh->connectivity()};
+ CellValue<double> energy_fluxes{mesh->connectivity()};
+
+ for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+ ur[node_id] = (1. / 6.) * ur1[node_id] + (1. / 3.) * ur2[node_id] + //
+ (1. / 3.) * ur3[node_id] + (1. / 6.) * ur4[node_id];
+ }
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ for (size_t i = 0; i < ul.sizeOfArrays(); ++i) {
+ ul[face_id][i] = (1. / 6.) * ul1[face_id][i] + (1. / 3.) * ul2[face_id][i] + //
+ (1. / 3.) * ul3[face_id][i] + (1. / 6.) * ul4[face_id][i];
+ }
+ }
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ momentum_fluxes[cell_id] = (1. / 6.) * momentum_fluxes1[cell_id] + (1. / 3.) * momentum_fluxes2[cell_id] +
+ (1. / 3.) * momentum_fluxes3[cell_id] + (1. / 6.) * momentum_fluxes4[cell_id];
+ }
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ energy_fluxes[cell_id] = (1. / 6.) * energy_fluxes1[cell_id] + (1. / 3.) * energy_fluxes2[cell_id] +
+ (1. / 3.) * energy_fluxes3[cell_id] + (1. / 6.) * energy_fluxes4[cell_id];
+ }
+
+ auto [mesh4, rho4, u4, E4] = _applyFluxes(dt, rho_v, u_v, E_v, ur, ul, momentum_fluxes, energy_fluxes);
+
+ return {mesh4, rho4, u4, E4};
+
+#elif 0
// Heun's RK3 method
auto [ur1, ul1, momentum_fluxes1, energy_fluxes1] =
- _computeFluxes(degree, velocity_bc_treatment, rho_v, u_v, E_v, rhoc, a, bc_descriptor_list, function_id);
+ _computeFluxes(polynomial_degree, velocity_bc_treatment, rho_v, u_v, E_v, rhoc, a, bc_descriptor_list,
+ function_id);
auto [mesh1, rho1, u1, E1] = _applyFluxes(dt / 3., rho_v, u_v, E_v, ur1, ul1, momentum_fluxes1, energy_fluxes1);
auto [ur2, ul2, momentum_fluxes2, energy_fluxes2] =
- _computeFluxes(degree, velocity_bc_treatment, rho1, u1, E1, rhoc, a, bc_descriptor_list, function_id);
+ _computeFluxes(polynomial_degree, velocity_bc_treatment, rho1, u1, E1, rhoc, a, bc_descriptor_list, function_id);
auto [mesh2, rho2, u2, E2] =
_applyFluxes(2. / 3. * dt, rho_v, u_v, E_v, ur2, ul2, momentum_fluxes2, energy_fluxes2);
auto [ur3, ul3, momentum_fluxes3, energy_fluxes3] =
- _computeFluxes(degree, velocity_bc_treatment, rho2, u2, E2, rhoc, a, bc_descriptor_list, function_id);
+ _computeFluxes(polynomial_degree, velocity_bc_treatment, rho2, u2, E2, rhoc, a, bc_descriptor_list, function_id);
NodeValue<Rd> ur{mesh->connectivity()};
FaceArray<Rd> ul{mesh->connectivity(), ul1.sizeOfArrays()};
@@ -826,18 +882,19 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
// Ralston's RK3 method
auto [ur1, ul1, momentum_fluxes1, energy_fluxes1] =
- _computeFluxes(degree, velocity_bc_treatment, rho_v, u_v, E_v, rhoc, a, bc_descriptor_list, function_id);
+ _computeFluxes(polynomial_degree, velocity_bc_treatment, rho_v, u_v, E_v, rhoc, a, bc_descriptor_list,
+ function_id);
auto [mesh1, rho1, u1, E1] = _applyFluxes(dt / 2., rho_v, u_v, E_v, ur1, ul1, momentum_fluxes1, energy_fluxes1);
auto [ur2, ul2, momentum_fluxes2, energy_fluxes2] =
- _computeFluxes(degree, velocity_bc_treatment, rho1, u1, E1, rhoc, a, bc_descriptor_list, function_id);
+ _computeFluxes(polynomial_degree, velocity_bc_treatment, rho1, u1, E1, rhoc, a, bc_descriptor_list, function_id);
auto [mesh2, rho2, u2, E2] =
_applyFluxes(3. / 4. * dt, rho_v, u_v, E_v, ur2, ul2, momentum_fluxes2, energy_fluxes2);
auto [ur3, ul3, momentum_fluxes3, energy_fluxes3] =
- _computeFluxes(degree, velocity_bc_treatment, rho2, u2, E2, rhoc, a, bc_descriptor_list, function_id);
+ _computeFluxes(polynomial_degree, velocity_bc_treatment, rho2, u2, E2, rhoc, a, bc_descriptor_list, function_id);
NodeValue<Rd> ur{mesh->connectivity()};
FaceArray<Rd> ul{mesh->connectivity(), ul1.sizeOfArrays()};
@@ -864,24 +921,39 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
auto [mesh3, rho3, u3, E3] = _applyFluxes(dt, rho_v, u_v, E_v, ur, ul, momentum_fluxes, energy_fluxes);
return {mesh3, rho3, u3, E3};
+#elif 0
+ // RK2
+ // k1
+ auto [ur1, ul1, momentum_fluxes1, energy_fluxes1] =
+ _computeFluxes(polynomial_degree, velocity_bc_treatment, rho_v, u_v, E_v, rhoc, a, bc_descriptor_list,
+ function_id);
+
+ // y1
+ auto [mesh1, rho1, u1, E1] = _applyFluxes(dt / 2., rho_v, u_v, E_v, ur1, ul1, momentum_fluxes1, energy_fluxes1);
+ // k2
+ auto [ur2, ul2, momentum_fluxes2, energy_fluxes2] =
+ _computeFluxes(polynomial_degree, velocity_bc_treatment, rho1, u1, E1, rhoc, a, bc_descriptor_list, function_id);
+
+ return _applyFluxes(dt, rho_v, u_v, E_v, ur2, ul2, momentum_fluxes2, energy_fluxes2);
#else
auto [ur1, ul1, momentum_fluxes1, energy_fluxes1] =
- _computeFluxes(degree, velocity_bc_treatment, rho_v, u_v, E_v, rhoc, a, bc_descriptor_list, function_id);
+ _computeFluxes(polynomial_degree, velocity_bc_treatment, rho_v, u_v, E_v, rhoc, a, bc_descriptor_list,
+ function_id);
return _applyFluxes(dt, rho_v, u_v, E_v, ur1, ul1, momentum_fluxes1, energy_fluxes1);
#endif
}
- VariationalSolver() : m_quadrature_degree(8) {}
+ VariationalSolver() : m_quadrature_degree(15) {}
VariationalSolver(VariationalSolver&&) = default;
~VariationalSolver() = default;
};
-template <MeshConcept MeshType, size_t mesh_edges_degree>
+template <MeshConcept MeshType, size_t mesh_degree>
void
-VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_applyPressureBC(
+VariationalSolverHandler::VariationalSolver<MeshType, mesh_degree>::_applyPressureBC(
const BoundaryConditionList& bc_list,
const MeshType& mesh,
Vector<double>& b) const
@@ -912,11 +984,11 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_apply
[](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 constexpr (mesh_edges_degree == 1) {
+ if constexpr (mesh_degree == 1) {
w_hat_set = w_hat_set_P1;
- } else if constexpr (mesh_edges_degree == 2) {
+ } else if constexpr (mesh_degree == 2) {
w_hat_set = w_hat_set_P2;
- } else if constexpr (mesh_edges_degree == 3) {
+ } else if constexpr (mesh_degree == 3) {
w_hat_set = w_hat_set_P3;
} else {
throw NotImplementedError("degree > 3");
@@ -931,7 +1003,7 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_apply
const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
- LineTransformationAccessor<MeshType, mesh_edges_degree> transformation{mesh};
+ LineTransformationAccessor<MeshType, mesh_degree> transformation{mesh};
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
FaceId face_id = face_list[i_face];
@@ -974,9 +1046,9 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_apply
}
}
-template <MeshConcept MeshType, size_t mesh_edges_degree>
+template <MeshConcept MeshType, size_t mesh_degree>
void
-VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_forcebcSymmetryBC(
+VariationalSolverHandler::VariationalSolver<MeshType, mesh_degree>::_forcebcSymmetryBC(
const BoundaryConditionList& bc_list,
Vector<double>& U) const
{
@@ -1015,7 +1087,7 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_force
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 < mesh_edges_degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_degree; ++r_hat) {
const size_t r = m_face_indices[face_id][r_hat];
TinyVector<Dimension> ur;
@@ -1037,9 +1109,9 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_force
}
}
-template <MeshConcept MeshType, size_t mesh_edges_degree>
+template <MeshConcept MeshType, size_t mesh_degree>
void
-VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_forcebcVelocityBC(
+VariationalSolverHandler::VariationalSolver<MeshType, mesh_degree>::_forcebcVelocityBC(
const BoundaryConditionList& bc_list,
Vector<double>& U) const
{
@@ -1068,7 +1140,7 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_force
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 < mesh_edges_degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_degree; ++r_hat) {
const size_t r = m_face_indices[face_id][r_hat];
for (size_t i = 0; i < Dimension; ++i) {
@@ -1095,7 +1167,7 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_force
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 < mesh_edges_degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_degree; ++r_hat) {
const size_t r = m_face_indices[face_id][r_hat];
for (size_t i = 0; i < Dimension; ++i) {
@@ -1109,9 +1181,9 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_force
}
}
-template <MeshConcept MeshType, size_t mesh_edges_degree>
+template <MeshConcept MeshType, size_t mesh_degree>
void
-VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_applySymmetryBC(
+VariationalSolverHandler::VariationalSolver<MeshType, mesh_degree>::_applySymmetryBC(
const BoundaryConditionList& bc_list,
const MeshType& mesh,
CRSMatrixDescriptor<double>& A_descriptor,
@@ -1167,7 +1239,7 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_apply
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 < mesh_edges_degree; ++i_hat) {
+ for (size_t i_hat = 1; i_hat < mesh_degree; ++i_hat) {
const size_t r = m_face_indices[face_id][i_hat];
TinyMatrix<Dimension> Arr;
@@ -1210,7 +1282,7 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_apply
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 <= mesh_edges_degree; ++i_hat) {
+ for (size_t i_hat = 0; i_hat <= mesh_degree; ++i_hat) {
const size_t s = m_face_indices[face_id][i_hat];
if (s != r) {
@@ -1235,10 +1307,10 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_apply
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 < mesh_edges_degree; ++i_hat) {
+ for (size_t i_hat = 1; i_hat < mesh_degree; ++i_hat) {
const size_t r = m_face_indices[face_id][i_hat];
- for (size_t j_hat = 0; j_hat <= mesh_edges_degree; ++j_hat) {
+ for (size_t j_hat = 0; j_hat <= mesh_degree; ++j_hat) {
const size_t s = m_face_indices[face_id][j_hat];
if (s != r) {
@@ -1265,9 +1337,9 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_apply
}
}
-template <MeshConcept MeshType, size_t mesh_edges_degree>
+template <MeshConcept MeshType, size_t mesh_degree>
void
-VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_applyVelocityBC(
+VariationalSolverHandler::VariationalSolver<MeshType, mesh_degree>::_applyVelocityBC(
const BoundaryConditionList& bc_list,
const MeshType& mesh,
const VariationalSolverHandler::VelocityBCTreatment& velocity_bc_treatment,
@@ -1299,7 +1371,7 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_apply
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 <= mesh_edges_degree; ++s_hat) {
+ for (size_t s_hat = 0; s_hat <= mesh_degree; ++s_hat) {
const size_t s = m_face_indices[face_id][s_hat];
if (r != s) {
@@ -1323,9 +1395,9 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_apply
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 < mesh_edges_degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_degree; ++r_hat) {
const size_t r = m_face_indices[face_id][r_hat];
- for (size_t s_hat = 0; s_hat <= mesh_edges_degree; ++s_hat) {
+ for (size_t s_hat = 0; s_hat <= mesh_degree; ++s_hat) {
const size_t s = m_face_indices[face_id][s_hat];
if (r != s) {
@@ -1355,7 +1427,7 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_apply
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 <= mesh_edges_degree; ++s_hat) {
+ for (size_t s_hat = 0; s_hat <= mesh_degree; ++s_hat) {
const size_t s = m_face_indices[face_id][s_hat];
if (r == s) {
@@ -1381,9 +1453,9 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_apply
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 < mesh_edges_degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_degree; ++r_hat) {
const size_t r = m_face_indices[face_id][r_hat];
- for (size_t s_hat = 0; s_hat <= mesh_edges_degree; ++s_hat) {
+ for (size_t s_hat = 0; s_hat <= mesh_degree; ++s_hat) {
const size_t s = m_face_indices[face_id][s_hat];
if (r == s) {
@@ -1420,7 +1492,7 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_apply
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 <= mesh_edges_degree; ++s_hat) {
+ for (size_t s_hat = 0; s_hat <= mesh_degree; ++s_hat) {
const size_t s = m_face_indices[face_id][s_hat];
if (r == s) {
@@ -1446,9 +1518,9 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_apply
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 < mesh_edges_degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_degree; ++r_hat) {
const size_t r = m_face_indices[face_id][r_hat];
- for (size_t s_hat = 0; s_hat <= mesh_edges_degree; ++s_hat) {
+ for (size_t s_hat = 0; s_hat <= mesh_degree; ++s_hat) {
const size_t s = m_face_indices[face_id][s_hat];
if (r == s) {
@@ -1493,7 +1565,7 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_apply
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 < mesh_edges_degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_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;
@@ -1512,7 +1584,7 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_apply
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 < mesh_edges_degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_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];
@@ -1533,7 +1605,7 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_apply
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 < mesh_edges_degree; ++r_hat) {
+ for (size_t r_hat = 1; r_hat < mesh_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;
@@ -1549,8 +1621,8 @@ VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::_apply
}
}
-template <MeshConcept MeshType, size_t mesh_edges_degree>
-class VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::FixedBoundaryCondition
+template <MeshConcept MeshType, size_t mesh_degree>
+class VariationalSolverHandler::VariationalSolver<MeshType, mesh_degree>::FixedBoundaryCondition
{
private:
const MeshNodeBoundary m_mesh_node_boundary;
@@ -1580,8 +1652,8 @@ class VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::
~FixedBoundaryCondition() = default;
};
-template <MeshConcept MeshType, size_t mesh_edges_degree>
-class VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::PressureBoundaryCondition
+template <MeshConcept MeshType, size_t mesh_degree>
+class VariationalSolverHandler::VariationalSolver<MeshType, mesh_degree>::PressureBoundaryCondition
{
private:
const MeshFaceBoundary m_mesh_face_boundary;
@@ -1607,8 +1679,8 @@ class VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::
~PressureBoundaryCondition() = default;
};
-template <MeshConcept MeshType, size_t mesh_edges_degree>
-class VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::VelocityBoundaryCondition
+template <MeshConcept MeshType, size_t mesh_degree>
+class VariationalSolverHandler::VariationalSolver<MeshType, mesh_degree>::VelocityBoundaryCondition
{
private:
const MeshNodeBoundary m_mesh_node_boundary;
@@ -1663,8 +1735,8 @@ class VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::
~VelocityBoundaryCondition() = default;
};
-template <MeshConcept MeshType, size_t mesh_edges_degree>
-class VariationalSolverHandler::VariationalSolver<MeshType, mesh_edges_degree>::SymmetryBoundaryCondition
+template <MeshConcept MeshType, size_t mesh_degree>
+class VariationalSolverHandler::VariationalSolver<MeshType, mesh_degree>::SymmetryBoundaryCondition
{
public:
using Rd = TinyVector<Dimension, double>;