diff --git a/src/geometry/LineTransformation.hpp b/src/geometry/LineTransformation.hpp index 9e6ffa09e2701220ea40d16ed367711d6fd18e3a..b9e9ada96ca028f07497c910f06476f8b67eb61f 100644 --- a/src/geometry/LineTransformation.hpp +++ b/src/geometry/LineTransformation.hpp @@ -73,12 +73,20 @@ class LineTransformation return m_velocity; } + PUGS_INLINE double velocityNorm() const { return m_velocity_norm; } + PUGS_INLINE + double + velocityNorm(const TinyVector<1>&) const + { + return m_velocity_norm; + } + PUGS_INLINE LineTransformation(const TinyVector<Dimension>& a, const TinyVector<Dimension>& b) : m_velocity{0.5 * (b - a)}, m_velocity_norm{l2Norm(m_velocity)}, m_shift{0.5 * (a + b)} diff --git a/src/scheme/VariationalSolver.cpp b/src/scheme/VariationalSolver.cpp index 043cc4bbd0807333aa1c30df45f369c457ff4c26..8ab8f549217fc30bb0824e312191c815df2033dc 100644 --- a/src/scheme/VariationalSolver.cpp +++ b/src/scheme/VariationalSolver.cpp @@ -24,6 +24,7 @@ #include <analysis/GaussQuadratureDescriptor.hpp> #include <analysis/QuadratureManager.hpp> +#include <geometry/LineCubicTransformation.hpp> #include <geometry/LineParabolicTransformation.hpp> #include <algebra/CRSMatrixDescriptor.hpp> @@ -38,53 +39,126 @@ #include <analysis/GaussLegendreQuadratureDescriptor.hpp> #include <language/utils/IntegrateOnCells.hpp> +template <MeshConcept MeshType, size_t mesh_degree> +struct LineTransformationAccessor +{ +}; + +template <> +struct LineTransformationAccessor<Mesh<2>, 1> +{ + const NodeValue<const TinyVector<2>> m_xr; + + const ItemToItemMatrix<ItemType::face, ItemType::node> m_face_to_node_matrix; + + LineTransformationAccessor(const Mesh<2>& mesh) + : m_xr{mesh.xr()}, m_face_to_node_matrix(mesh.connectivity().faceToNodeMatrix()) + {} + + auto + getTransformation(const FaceId face_id) const + { + auto node_list = m_face_to_node_matrix[face_id]; + return LineTransformation<2>{m_xr[node_list[0]], m_xr[node_list[1]]}; + } +}; + +template <size_t mesh_degree> +struct LineTransformationAccessor<PolynomialMesh<2>, mesh_degree> +{ + const NodeValue<const TinyVector<2>> m_xr; + const FaceArray<const TinyVector<2>> m_xl; + + const ItemToItemMatrix<ItemType::face, ItemType::node> m_face_to_node_matrix; + + LineTransformationAccessor(const PolynomialMesh<2>& mesh) + : m_xr{mesh.xr()}, m_xl{mesh.xl()}, m_face_to_node_matrix(mesh.connectivity().faceToNodeMatrix()) + {} + + auto + getTransformation(const FaceId face_id) const + { + auto node_list = m_face_to_node_matrix[face_id]; + if constexpr (mesh_degree == 1) { + return LineTransformation<2>{m_xr[node_list[0]], m_xr[node_list[1]]}; + } else if constexpr (mesh_degree == 2) { + return LineParabolicTransformation<2>{m_xr[node_list[0]], m_xl[face_id][0], m_xr[node_list[1]]}; + } else if constexpr (mesh_degree == 3) { + return LineCubicTransformation<2>{m_xr[node_list[0]], m_xl[face_id][0], m_xl[face_id][1], m_xr[node_list[1]]}; + } else { + static_assert(mesh_degree == 1, "mesh degree is not supported"); + } + } +}; + +template <size_t mesh_degree, MeshConcept MeshType> double -variational_acoustic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c_v) +variational_acoustic_dt(const DiscreteFunctionP0<const double>& c, const std::shared_ptr<const MeshType>& p_mesh) { - const auto& c = c_v->get<DiscreteFunctionP0<const double>>(); + const auto& mesh = *p_mesh; + const auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj(); - return std::visit( - [&](auto&& p_mesh) -> double { - const auto& mesh = *p_mesh; + auto cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix(); - using MeshType = mesh_type_t<decltype(mesh)>; - if constexpr (is_polynomial_mesh_v<MeshType>) { - const auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj(); - if (mesh.degree() != 2) { - throw NotImplementedError("mesh degree != 2"); - } + auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2)); - auto xr = mesh.xr(); - auto xl = mesh.xl(); + CellValue<double> Sj{mesh.connectivity()}; + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + LineTransformationAccessor<MeshType, mesh_degree> transform(*p_mesh); - auto cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix(); - auto face_to_node_matrix = mesh.connectivity().faceToNodeMatrix(); + double sum = 0; + auto face_list = cell_to_face_matrix[cell_id]; + for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { + const FaceId face_id = face_list[i_face]; - auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2)); + auto t = transform.getTransformation(face_id); - CellValue<double> Sj{mesh.connectivity()}; - parallel_for( - mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { - double sum = 0; - auto face_list = cell_to_face_matrix[cell_id]; - for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { - const FaceId face_id = face_list[i_face]; - auto node_list = face_to_node_matrix[face_id]; + for (size_t iq = 0; iq < qf.numberOfPoints(); ++iq) { + sum += qf.weight(iq) * t.velocityNorm(qf.point(iq)); + } + } + Sj[cell_id] = sum; + }); - LineParabolicTransformation<MeshType::Dimension> t(xr[node_list[0]], xl[face_id][0], xr[node_list[1]]); + CellValue<double> local_dt{mesh.connectivity()}; + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { local_dt[j] = 2 * Vj[j] / (Sj[j] * c[j]); }); - for (size_t iq = 0; iq < qf.numberOfPoints(); ++iq) { - sum += qf.weight(iq) * t.velocityNorm(qf.point(iq)); - } - } - Sj[cell_id] = sum; - }); + return min(local_dt); +} + +double +variational_acoustic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c_v) +{ + const auto& c = c_v->get<DiscreteFunctionP0<const double>>(); - CellValue<double> local_dt{mesh.connectivity()}; - parallel_for( - mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { local_dt[j] = 2 * Vj[j] / (Sj[j] * c[j]); }); + return std::visit( + [&](auto&& p_mesh) -> double { + const auto& mesh = *p_mesh; - return min(local_dt); + using MeshType = mesh_type_t<decltype(mesh)>; + if constexpr (is_polygonal_mesh_v<MeshType>) { + if constexpr (MeshType::Dimension == 2) { + return variational_acoustic_dt<1>(c, p_mesh); + } else { + throw NotImplementedError("not implemented in dimension d != 2"); + } + } else if constexpr (is_polynomial_mesh_v<MeshType>) { + switch (mesh.degree()) { + case 1: { + return variational_acoustic_dt<1>(c, p_mesh); + } + case 2: { + return variational_acoustic_dt<2>(c, p_mesh); + } + case 3: { + return variational_acoustic_dt<3>(c, p_mesh); + } + default: { + throw NotImplementedError("not implemented for mesh degree > 3"); + } + } } else { throw NormalError("unexpected mesh type"); } @@ -206,23 +280,13 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv Table<TinyVector<Dimension>> quadrature_points(mesh_face_boundary.faceList().size(), qf.numberOfPoints()); + LineTransformationAccessor<MeshType, 2> transformation{mesh}; + 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 if (mesh.degree() == 1) { - return 0.5 * (x0 + x2); - } else { - throw NotImplementedError("degree>2"); - } - }(); + auto t = transformation.getTransformation(face_id); - LineParabolicTransformation<Dimension> t(x0, x1, x2); auto face_quadrature_points = quadrature_points[i_face]; for (size_t i_xi = 0; i_xi < qf.numberOfPoints(); ++i_xi) { face_quadrature_points[i_xi] = t(qf.point(i_xi)); @@ -731,7 +795,7 @@ 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 1 +#if 0 // Heun's RK3 method auto [ur1, ul1, momentum_fluxes1, energy_fluxes1] =