From 5902019f9e0936c64048c64e3cc90222691609ff Mon Sep 17 00:00:00 2001 From: Johannes Heinz Date: Mon, 5 Jun 2023 15:20:40 +0200 Subject: [PATCH] use FEEval::normal_vector() to be consistent with FEValues and FEPEval --- doc/doxygen/headers/matrixfree.h | 2 +- examples/step-59/step-59.cc | 16 +- examples/step-67/step-67.cc | 4 +- examples/step-76/step-76.cc | 4 +- include/deal.II/matrix_free/face_info.h | 4 +- include/deal.II/matrix_free/fe_evaluation.h | 2 +- .../deal.II/matrix_free/fe_evaluation_data.h | 21 +- tests/matrix_free/advect_1d.cc | 4 +- tests/matrix_free/advect_1d_deprecated.cc | 327 ++++++++++++++++++ tests/matrix_free/advect_1d_deprecated.output | 22 ++ tests/matrix_free/advect_1d_system.cc | 4 +- .../advect_1d_vectorization_mask.cc | 4 +- tests/matrix_free/compare_faces_by_cells.cc | 12 +- tests/matrix_free/ecl.h | 18 +- .../matrix_free/element_birth_and_death_01.cc | 2 +- .../interpolate_functions_common.h | 4 +- tests/matrix_free/matrix_vector_faces_35.cc | 4 +- .../matrix_free/matrix_vector_faces_common.h | 20 +- tests/matrix_free/multigrid_dg_periodic.cc | 18 +- tests/matrix_free/multigrid_dg_sip_01.cc | 18 +- tests/matrix_free/multigrid_dg_sip_02.cc | 12 +- tests/performance/timing_navier_stokes.cc | 20 +- tests/simplex/matrix_free_03.cc | 6 +- tests/simplex/step-12c.cc | 10 +- tests/simplex/step-67.cc | 4 +- 25 files changed, 461 insertions(+), 101 deletions(-) create mode 100644 tests/matrix_free/advect_1d_deprecated.cc create mode 100644 tests/matrix_free/advect_1d_deprecated.output diff --git a/doc/doxygen/headers/matrixfree.h b/doc/doxygen/headers/matrixfree.h index a153519370..67cdf52fb0 100644 --- a/doc/doxygen/headers/matrixfree.h +++ b/doc/doxygen/headers/matrixfree.h @@ -211,7 +211,7 @@ digraph G * information, respectively. Besides access to the function values with * FEEvaluationAccess::get_value() or gradients with * FEEvaluationAccess::get_gradient(), the face evaluator also enables the - * access to the normal vector by FEEvaluationAccess::get_normal_vector() and + * access to the normal vector by FEEvaluationAccess::normal_vector() and * a specialized field FEEvaluationAccess::get_normal_derivative(), which * returns the derivative of the solution field normal to the face. This * quantity is computed as the gradient (in real space) multiplied by the diff --git a/examples/step-59/step-59.cc b/examples/step-59/step-59.cc index 980185edef..a079ff559f 100644 --- a/examples/step-59/step-59.cc +++ b/examples/step-59/step-59.cc @@ -579,9 +579,9 @@ namespace Step59 // take the absolute value of these factors as the normal could point // into either positive or negative direction. const VectorizedArray inverse_length_normal_to_face = - 0.5 * (std::abs((phi_inner.get_normal_vector(0) * + 0.5 * (std::abs((phi_inner.normal_vector(0) * phi_inner.inverse_jacobian(0))[dim - 1]) + - std::abs((phi_outer.get_normal_vector(0) * + std::abs((phi_outer.normal_vector(0) * phi_outer.inverse_jacobian(0))[dim - 1])); const VectorizedArray sigma = inverse_length_normal_to_face * get_penalty_factor(); @@ -680,9 +680,8 @@ namespace Step59 EvaluationFlags::values | EvaluationFlags::gradients); - const VectorizedArray inverse_length_normal_to_face = - std::abs((phi_inner.get_normal_vector(0) * - phi_inner.inverse_jacobian(0))[dim - 1]); + const VectorizedArray inverse_length_normal_to_face = std::abs(( + phi_inner.normal_vector(0) * phi_inner.inverse_jacobian(0))[dim - 1]); const VectorizedArray sigma = inverse_length_normal_to_face * get_penalty_factor(); @@ -1134,9 +1133,8 @@ namespace Step59 { phi_face.reinit(face); - const VectorizedArray inverse_length_normal_to_face = - std::abs((phi_face.get_normal_vector(0) * - phi_face.inverse_jacobian(0))[dim - 1]); + const VectorizedArray inverse_length_normal_to_face = std::abs( + (phi_face.normal_vector(0) * phi_face.inverse_jacobian(0))[dim - 1]); const VectorizedArray sigma = inverse_length_normal_to_face * system_matrix.get_penalty_factor(); @@ -1169,7 +1167,7 @@ namespace Step59 { Tensor<1, dim> normal; for (unsigned int d = 0; d < dim; ++d) - normal[d] = phi_face.get_normal_vector(q)[d][v]; + normal[d] = phi_face.normal_vector(q)[d][v]; test_normal_derivative[v] = -normal * exact_solution.gradient(single_point); } diff --git a/examples/step-67/step-67.cc b/examples/step-67/step-67.cc index c76c6e09db..95a814c5c3 100644 --- a/examples/step-67/step-67.cc +++ b/examples/step-67/step-67.cc @@ -1125,7 +1125,7 @@ namespace Euler_DG const auto numerical_flux = euler_numerical_flux(phi_m.get_value(q), phi_p.get_value(q), - phi_m.get_normal_vector(q)); + phi_m.normal_vector(q)); phi_m.submit_value(-numerical_flux, q); phi_p.submit_value(numerical_flux, q); } @@ -1205,7 +1205,7 @@ namespace Euler_DG for (unsigned int q = 0; q < phi.n_q_points; ++q) { const auto w_m = phi.get_value(q); - const auto normal = phi.get_normal_vector(q); + const auto normal = phi.normal_vector(q); auto rho_u_dot_n = w_m[1] * normal[0]; for (unsigned int d = 1; d < dim; ++d) diff --git a/examples/step-76/step-76.cc b/examples/step-76/step-76.cc index c66a6ac44b..d720b53409 100644 --- a/examples/step-76/step-76.cc +++ b/examples/step-76/step-76.cc @@ -812,7 +812,7 @@ namespace Euler_DG const auto numerical_flux = euler_numerical_flux(phi_m.get_value(q), phi_p.get_value(q), - phi_m.get_normal_vector(q)); + phi_m.normal_vector(q)); phi_m.submit_value(-numerical_flux, q); } } @@ -825,7 +825,7 @@ namespace Euler_DG for (unsigned int q = 0; q < phi_m.n_q_points; ++q) { const auto w_m = phi_m.get_value(q); - const auto normal = phi_m.get_normal_vector(q); + const auto normal = phi_m.normal_vector(q); auto rho_u_dot_n = w_m[1] * normal[0]; for (unsigned int d = 1; d < dim; ++d) diff --git a/include/deal.II/matrix_free/face_info.h b/include/deal.II/matrix_free/face_info.h index 1fc018fe95..180bfd009b 100644 --- a/include/deal.II/matrix_free/face_info.h +++ b/include/deal.II/matrix_free/face_info.h @@ -56,7 +56,7 @@ namespace internal /** * Indices of the faces in the current face batch as compared to the * numbers of the cells on the logical "interior" side of the face which - * is aligned to the direction of FEEvaluation::get_normal_vector(). + * is aligned to the direction of FEEvaluation::normal_vector(). */ std::array cells_interior; @@ -64,7 +64,7 @@ namespace internal * Indices of the faces in the current face batch as compared to the * numbers of the cells on the logical "exterior" side of the face which * is aligned to the opposite direction of - * FEEvaluation::get_normal_vector(). Note that the distinction into + * FEEvaluation::normal_vector(). Note that the distinction into * interior and exterior faces is purely logical and refers to the * direction of the normal only. In the actual discretization of a * problem, the discretization typically needs to make sure that interior diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 077067d7a6..fb482e5a25 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -365,7 +365,7 @@ public: * to the face: $\boldsymbol \nabla u(\mathbf x_q) \cdot \mathbf n(\mathbf * x_q)$ * - * This call is equivalent to calling get_gradient() * get_normal_vector() + * This call is equivalent to calling get_gradient() * normal_vector() * but will use a more efficient internal representation of data. * * @note The derived class FEEvaluationAccess overloads this operation diff --git a/include/deal.II/matrix_free/fe_evaluation_data.h b/include/deal.II/matrix_free/fe_evaluation_data.h index 20e2fb2ab2..b188a05a6b 100644 --- a/include/deal.II/matrix_free/fe_evaluation_data.h +++ b/include/deal.II/matrix_free/fe_evaluation_data.h @@ -227,6 +227,14 @@ public: * @note Only implemented in case `is_face == true`. */ Tensor<1, dim, Number> + normal_vector(const unsigned int q_point) const; + + /** + * Same as `normal_vector(const unsigned int q_point)`. + * + * @warning This function will be deprecated! + */ + Tensor<1, dim, Number> get_normal_vector(const unsigned int q_point) const; /** @} */ @@ -1260,7 +1268,7 @@ FEEvaluationData::reinit_face( template inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number> -FEEvaluationData::get_normal_vector( +FEEvaluationData::normal_vector( const unsigned int q_point) const { AssertIndexRange(q_point, n_quadrature_points); @@ -1275,6 +1283,17 @@ FEEvaluationData::get_normal_vector( +// This function is deprecated. +template +inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number> +FEEvaluationData::get_normal_vector( + const unsigned int q_point) const +{ + return normal_vector(q_point); +} + + + template inline DEAL_II_ALWAYS_INLINE Number FEEvaluationData::JxW(const unsigned int q_point) const diff --git a/tests/matrix_free/advect_1d.cc b/tests/matrix_free/advect_1d.cc index 063de99da3..af5d943bba 100644 --- a/tests/matrix_free/advect_1d.cc +++ b/tests/matrix_free/advect_1d.cc @@ -151,7 +151,7 @@ private: value_type u_minus = phi_m.get_value(q), u_plus = phi_p.get_value(q); const VectorizedArray normal_times_advection = - advection * phi_m.get_normal_vector(q); + advection * phi_m.normal_vector(q); const value_type flux_times_normal = make_vectorized_array(0.5) * ((u_minus + u_plus) * normal_times_advection + @@ -193,7 +193,7 @@ private: value_type u_minus = fe_eval.get_value(q); value_type u_plus = -u_minus; const VectorizedArray normal_times_advection = - advection * fe_eval.get_normal_vector(q); + advection * fe_eval.normal_vector(q); const value_type flux_times_normal = make_vectorized_array(0.5) * ((u_minus + u_plus) * normal_times_advection + diff --git a/tests/matrix_free/advect_1d_deprecated.cc b/tests/matrix_free/advect_1d_deprecated.cc new file mode 100644 index 0000000000..9374fb65ca --- /dev/null +++ b/tests/matrix_free/advect_1d_deprecated.cc @@ -0,0 +1,327 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 - 2022 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// Same as advect_1d.cc using the deprecated FEEval::get_normal_vector() +// to ensure user code is not broken. This test can be removed once +// get_normal_vector() is removed. + +#include +#include + +#include + +#include +#include +#include + +#include + +#include +#include +#include + +#include + +#include "../tests.h" + + + +template , + int n_components = 1> +class MatrixFreeAdvectionBasic +{ +public: + MatrixFreeAdvectionBasic(const MatrixFree &data, + const bool zero_within_loop = true, + const unsigned int start_vector_component = 0) + : data(data) + , zero_within_loop(zero_within_loop) + , start_vector_component(start_vector_component) + { + for (unsigned int d = 0; d < dim; ++d) + advection[d] = 0.4 + 0.12 * d; + } + + void + vmult(VectorType &dst, const VectorType &src) const + { + if (!zero_within_loop) + dst = 0; + data.loop(&MatrixFreeAdvectionBasic::local_apply, + &MatrixFreeAdvectionBasic::local_apply_face, + &MatrixFreeAdvectionBasic::local_apply_boundary_face, + this, + dst, + src, + zero_within_loop, + MatrixFree::DataAccessOnFaces::values, + MatrixFree::DataAccessOnFaces::values); + + FEEvaluation phi(data); + + const unsigned int dofs_per_cell = phi.dofs_per_cell; + + AlignedVector> coefficients(phi.dofs_per_cell); + MatrixFreeOperators:: + CellwiseInverseMassMatrix + inverse(phi); + + for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell) + { + phi.reinit(cell); + phi.read_dof_values(dst); + + inverse.fill_inverse_JxW_values(coefficients); + inverse.apply(coefficients, + n_components, + phi.begin_dof_values(), + phi.begin_dof_values()); + + phi.set_dof_values(dst); + } + } + +private: + void + local_apply(const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range) const + { + FEEvaluation phi( + data, 0, 0, start_vector_component); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values(src); + phi.evaluate(EvaluationFlags::values); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(advection * phi.get_value(q), q); + phi.integrate(EvaluationFlags::gradients); + phi.distribute_local_to_global(dst); + } + } + + void + local_apply_face( + const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &face_range) const + { + FEFaceEvaluation phi_m( + data, true, 0, 0, start_vector_component); + FEFaceEvaluation phi_p( + data, false, 0, 0, start_vector_component); + using value_type = typename FEFaceEvaluation::value_type; + + for (unsigned int face = face_range.first; face < face_range.second; ++face) + { + phi_m.reinit(face); + phi_m.read_dof_values(src); + phi_m.evaluate(EvaluationFlags::values); + phi_p.reinit(face); + phi_p.read_dof_values(src); + phi_p.evaluate(EvaluationFlags::values); + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + value_type u_minus = phi_m.get_value(q), + u_plus = phi_p.get_value(q); + const VectorizedArray normal_times_advection = + advection * phi_m.get_normal_vector(q); + const value_type flux_times_normal = + make_vectorized_array(0.5) * + ((u_minus + u_plus) * normal_times_advection + + std::abs(normal_times_advection) * (u_minus - u_plus)); + phi_m.submit_value(-flux_times_normal, q); + phi_p.submit_value(flux_times_normal, q); + } + + phi_m.integrate(EvaluationFlags::values); + phi_m.distribute_local_to_global(dst); + phi_p.integrate(EvaluationFlags::values); + phi_p.distribute_local_to_global(dst); + } + } + + void + local_apply_boundary_face( + const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &face_range) const + { + FEFaceEvaluation + fe_eval(data, true, 0, 0, start_vector_component); + using value_type = typename FEFaceEvaluation::value_type; + + for (unsigned int face = face_range.first; face < face_range.second; ++face) + { + fe_eval.reinit(face); + fe_eval.read_dof_values(src); + fe_eval.evaluate(EvaluationFlags::values); + + for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) + { + value_type u_minus = fe_eval.get_value(q); + value_type u_plus = -u_minus; + const VectorizedArray normal_times_advection = + advection * fe_eval.get_normal_vector(q); + const value_type flux_times_normal = + make_vectorized_array(0.5) * + ((u_minus + u_plus) * normal_times_advection + + std::abs(normal_times_advection) * (u_minus - u_plus)); + fe_eval.submit_value(-flux_times_normal, q); + } + + fe_eval.integrate(EvaluationFlags::values); + fe_eval.distribute_local_to_global(dst); + } + } + + const MatrixFree & data; + const bool zero_within_loop; + const unsigned int start_vector_component; + Tensor<1, dim, VectorizedArray> advection; +}; + + + +template +class AnalyticFunction : public Function +{ +public: + static_assert(dim == 1, "Only 1D implemented"); + AnalyticFunction() + : Function(1) + {} + + virtual double + value(const Point &p, const unsigned int) const override + { + return std::sin(3 * numbers::PI * p[0] / 0.8); + } +}; + + + +template +class AnalyticDerivative : public Function +{ +public: + static_assert(dim == 1, "Only 1D implemented"); + AnalyticDerivative() + : Function(1) + {} + + virtual double + value(const Point &p, const unsigned int) const override + { + Tensor<1, dim> advection; + for (unsigned int d = 0; d < dim; ++d) + advection[d] = 0.4 + 0.12 * d; + + return -std::cos(3 * numbers::PI * p[0] / 0.8) * advection[0] * 3 * + numbers::PI / 0.8; + } +}; + + + +template +void +test(const unsigned int n_refine) +{ + Triangulation tria; + GridGenerator::hyper_cube(tria, 0, 0.8); + tria.refine_global(n_refine); + + FE_DGQ fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + constraints.close(); + + if (n_refine == 3) + { + deallog << "Testing " << dof.get_fe().get_name(); + deallog << std::endl; + } + + LinearAlgebra::distributed::Vector in, out; + + const QGauss<1> quad(fe_degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags_inner_faces = + (update_gradients | update_JxW_values); + data.mapping_update_flags_boundary_faces = + (update_gradients | update_JxW_values); + + MatrixFree mf_data; + mf_data.reinit(MappingQ1{}, dof, constraints, quad, data); + + mf_data.initialize_dof_vector(in); + mf_data.initialize_dof_vector(out); + + VectorTools::interpolate(dof, AnalyticFunction(), in); + + MatrixFreeAdvectionBasic> + mf2(mf_data); + mf2.vmult(out, in); + + VectorTools::interpolate(dof, AnalyticDerivative(), in); + out -= in; + + double diff_norm = out.linfty_norm(); + deallog << "Norm of difference: " << diff_norm << ' ' << std::endl; +} + + + +int +main() +{ + initlog(); + + for (unsigned int r = 3; r < 9; ++r) + test<1, 2>(r); + + for (unsigned int r = 3; r < 9; ++r) + test<1, 4>(r); + + for (unsigned int r = 3; r < 9; ++r) + test<1, 5>(r); +} diff --git a/tests/matrix_free/advect_1d_deprecated.output b/tests/matrix_free/advect_1d_deprecated.output new file mode 100644 index 0000000000..b213520047 --- /dev/null +++ b/tests/matrix_free/advect_1d_deprecated.output @@ -0,0 +1,22 @@ + +DEAL::Testing FE_DGQ<1>(2) +DEAL::Norm of difference: 0.531539 +DEAL::Norm of difference: 0.135407 +DEAL::Norm of difference: 0.0340112 +DEAL::Norm of difference: 0.00851279 +DEAL::Norm of difference: 0.00212882 +DEAL::Norm of difference: 0.000532245 +DEAL::Testing FE_DGQ<1>(4) +DEAL::Norm of difference: 0.00529591 +DEAL::Norm of difference: 0.000336014 +DEAL::Norm of difference: 2.10801e-05 +DEAL::Norm of difference: 1.31874e-06 +DEAL::Norm of difference: 8.24413e-08 +DEAL::Norm of difference: 5.15331e-09 +DEAL::Testing FE_DGQ<1>(5) +DEAL::Norm of difference: 0.000346887 +DEAL::Norm of difference: 1.09982e-05 +DEAL::Norm of difference: 3.44938e-07 +DEAL::Norm of difference: 1.07881e-08 +DEAL::Norm of difference: 3.35085e-10 +DEAL::Norm of difference: 9.72689e-12 diff --git a/tests/matrix_free/advect_1d_system.cc b/tests/matrix_free/advect_1d_system.cc index e1d0ddd8ec..313b5d4db3 100644 --- a/tests/matrix_free/advect_1d_system.cc +++ b/tests/matrix_free/advect_1d_system.cc @@ -158,7 +158,7 @@ private: value_type u_minus = phi_m.get_value(q), u_plus = phi_p.get_value(q); const VectorizedArray normal_times_advection = - advection * phi_m.get_normal_vector(q); + advection * phi_m.normal_vector(q); const value_type flux_times_normal = make_vectorized_array(0.5) * ((u_minus + u_plus) * normal_times_advection + @@ -200,7 +200,7 @@ private: value_type u_minus = fe_eval.get_value(q); value_type u_plus = -u_minus; const VectorizedArray normal_times_advection = - advection * fe_eval.get_normal_vector(q); + advection * fe_eval.normal_vector(q); const value_type flux_times_normal = make_vectorized_array(0.5) * ((u_minus + u_plus) * normal_times_advection + diff --git a/tests/matrix_free/advect_1d_vectorization_mask.cc b/tests/matrix_free/advect_1d_vectorization_mask.cc index db7ea98e0b..d83d6827be 100644 --- a/tests/matrix_free/advect_1d_vectorization_mask.cc +++ b/tests/matrix_free/advect_1d_vectorization_mask.cc @@ -159,7 +159,7 @@ private: value_type u_minus = phi_m.get_value(q), u_plus = phi_p.get_value(q); const VectorizedArray normal_times_advection = - advection * phi_m.get_normal_vector(q); + advection * phi_m.normal_vector(q); const value_type flux_times_normal = make_vectorized_array(0.5) * ((u_minus + u_plus) * normal_times_advection + @@ -213,7 +213,7 @@ private: value_type u_minus = fe_eval.get_value(q); value_type u_plus = -u_minus; const VectorizedArray normal_times_advection = - advection * fe_eval.get_normal_vector(q); + advection * fe_eval.normal_vector(q); const value_type flux_times_normal = make_vectorized_array(0.5) * ((u_minus + u_plus) * normal_times_advection + diff --git a/tests/matrix_free/compare_faces_by_cells.cc b/tests/matrix_free/compare_faces_by_cells.cc index 8e23091004..6d1c18b1ae 100644 --- a/tests/matrix_free/compare_faces_by_cells.cc +++ b/tests/matrix_free/compare_faces_by_cells.cc @@ -142,8 +142,7 @@ private: // penalty parameter because the cell-based method cannot read the // sigma parameter on the neighbor VectorizedArray sigmaF = - std::abs( - (phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]) * + std::abs((phi.normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]) * (number)(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0; // Compute phi part @@ -178,7 +177,7 @@ private: phi.distribute_local_to_global(dst); // Compute phi_outer part - sigmaF = std::abs((phi.get_normal_vector(0) * + sigmaF = std::abs((phi.normal_vector(0) * phi_outer.inverse_jacobian(0))[dim - 1]) * (number)(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0; for (unsigned int j = 0; j < phi.dofs_per_cell; ++j) @@ -229,8 +228,7 @@ private: VectorizedArray local_diagonal_vector[phi.static_dofs_per_cell]; VectorizedArray sigmaF = - std::abs( - (phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]) * + std::abs((phi.normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]) * (number)(std::max(1, fe_degree) * (fe_degree + 1.0)) * 2.; for (unsigned int i = 0; i < phi.dofs_per_cell; ++i) @@ -277,8 +275,8 @@ private: { phif.reinit(cell, face); VectorizedArray sigmaF = - std::abs((phif.get_normal_vector(0) * - phif.inverse_jacobian(0))[dim - 1]) * + std::abs( + (phif.normal_vector(0) * phif.inverse_jacobian(0))[dim - 1]) * (number)(std::max(1, fe_degree) * (fe_degree + 1.0)) * 2.; std::array::size()> diff --git a/tests/matrix_free/ecl.h b/tests/matrix_free/ecl.h index c2d0dc9759..b5e253fad2 100644 --- a/tests/matrix_free/ecl.h +++ b/tests/matrix_free/ecl.h @@ -167,10 +167,10 @@ test(const unsigned int geometry = 0, phi_p.read_dof_values(src); phi_p.evaluate(true, true); VectorizedArrayType sigmaF = - (std::abs((phi_m.get_normal_vector(0) * - phi_m.inverse_jacobian(0))[dim - 1]) + - std::abs((phi_m.get_normal_vector(0) * - phi_p.inverse_jacobian(0))[dim - 1])) * + (std::abs( + (phi_m.normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) + + std::abs( + (phi_m.normal_vector(0) * phi_p.inverse_jacobian(0))[dim - 1])) * (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0)); for (unsigned int q = 0; q < phi_m.n_q_points; ++q) @@ -200,8 +200,8 @@ test(const unsigned int geometry = 0, phi_m.read_dof_values(src); phi_m.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); VectorizedArrayType sigmaF = - std::abs((phi_m.get_normal_vector(0) * - phi_m.inverse_jacobian(0))[dim - 1]) * + std::abs( + (phi_m.normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) * Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0; for (unsigned int q = 0; q < phi_m.n_q_points; ++q) @@ -258,7 +258,7 @@ test(const unsigned int geometry = 0, phi_m.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); VectorizedArrayType sigmaF = - std::abs((phi_m.get_normal_vector(0) * + std::abs((phi_m.normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) * Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0; @@ -291,9 +291,9 @@ test(const unsigned int geometry = 0, EvaluationFlags::gradients); VectorizedArrayType sigmaF = - (std::abs((phi_m.get_normal_vector(0) * + (std::abs((phi_m.normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) + - std::abs((phi_m.get_normal_vector(0) * + std::abs((phi_m.normal_vector(0) * phi_p.inverse_jacobian(0))[dim - 1])) * (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0)); diff --git a/tests/matrix_free/element_birth_and_death_01.cc b/tests/matrix_free/element_birth_and_death_01.cc index 53171d0366..ffa5533ff4 100644 --- a/tests/matrix_free/element_birth_and_death_01.cc +++ b/tests/matrix_free/element_birth_and_death_01.cc @@ -150,7 +150,7 @@ test(const unsigned int n_refinements) for (const auto q : phi.quadrature_point_indices()) { auto gradient = phi.get_gradient(q); - auto normal = phi.get_normal_vector(q); + auto normal = phi.normal_vector(q); if (is_interior_face == false) // fix sign! normal *= -1.0; diff --git a/tests/matrix_free/interpolate_functions_common.h b/tests/matrix_free/interpolate_functions_common.h index 05253c65b3..6118067ef3 100644 --- a/tests/matrix_free/interpolate_functions_common.h +++ b/tests/matrix_free/interpolate_functions_common.h @@ -162,7 +162,7 @@ public: double normal_derivative = 0; for (unsigned int d = 0; d < dim; ++d) normal_derivative += function.gradient(p, 0)[d] * - fe_evalm.get_normal_vector(q)[d][j]; + fe_evalm.normal_vector(q)[d][j]; facem_errors[3] += std::abs( fe_evalm.get_normal_derivative(q)[j] - normal_derivative); @@ -228,7 +228,7 @@ public: double normal_derivative = 0; for (unsigned int d = 0; d < dim; ++d) normal_derivative += function.gradient(p, 0)[d] * - fe_evalm.get_normal_vector(q)[d][j]; + fe_evalm.normal_vector(q)[d][j]; boundary_errors[3] += std::abs( fe_evalm.get_normal_derivative(q)[j] - normal_derivative); } diff --git a/tests/matrix_free/matrix_vector_faces_35.cc b/tests/matrix_free/matrix_vector_faces_35.cc index a87ab23c5c..a0b2587573 100644 --- a/tests/matrix_free/matrix_vector_faces_35.cc +++ b/tests/matrix_free/matrix_vector_faces_35.cc @@ -196,7 +196,7 @@ private: value_type u_minus = phi_m.get_value(q), u_plus = phi_p.get_value(q); const VectorizedArrayType normal_times_advection = - advection * phi_m.get_normal_vector(q); + advection * phi_m.normal_vector(q); const value_type flux_times_normal = make_vectorized_array(0.5) * ((u_minus + u_plus) * normal_times_advection + @@ -247,7 +247,7 @@ private: { value_type u_minus = fe_eval.get_value(q); const VectorizedArrayType normal_times_advection = - advection * fe_eval.get_normal_vector(q); + advection * fe_eval.normal_vector(q); const value_type flux_times_normal = make_vectorized_array(0.5) * ((u_minus + u_plus) * normal_times_advection + diff --git a/tests/matrix_free/matrix_vector_faces_common.h b/tests/matrix_free/matrix_vector_faces_common.h index e9c73d46b1..5463777933 100644 --- a/tests/matrix_free/matrix_vector_faces_common.h +++ b/tests/matrix_free/matrix_vector_faces_common.h @@ -156,9 +156,9 @@ private: fe_eval_neighbor.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); VectorizedArrayType sigmaF = - (std::abs((fe_eval.get_normal_vector(0) * + (std::abs((fe_eval.normal_vector(0) * fe_eval.inverse_jacobian(0))[dim - 1]) + - std::abs((fe_eval.get_normal_vector(0) * + std::abs((fe_eval.normal_vector(0) * fe_eval_neighbor.inverse_jacobian(0))[dim - 1])) * (number)(std::max(actual_degree, 1) * (actual_degree + 1.0)); @@ -217,8 +217,8 @@ private: fe_eval.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); VectorizedArrayType sigmaF = 2.0 * - std::abs((fe_eval.get_normal_vector(0) * - fe_eval.inverse_jacobian(0))[dim - 1]) * + std::abs( + (fe_eval.normal_vector(0) * fe_eval.inverse_jacobian(0))[dim - 1]) * number(std::max(actual_degree, 1) * (actual_degree + 1.0)); for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) @@ -362,9 +362,9 @@ private: EvaluationFlags::gradients); VectorizedArrayType sigmaF = - (std::abs((fe_eval.get_normal_vector(0) * + (std::abs((fe_eval.normal_vector(0) * fe_eval.inverse_jacobian(0))[dim - 1]) + - std::abs((fe_eval.get_normal_vector(0) * + std::abs((fe_eval.normal_vector(0) * fe_eval_neighbor.inverse_jacobian(0))[dim - 1])) * (number)(std::max(actual_degree, 1) * (actual_degree + 1.0)); @@ -423,8 +423,8 @@ private: EvaluationFlags::values | EvaluationFlags::gradients); VectorizedArrayType sigmaF = - std::abs((fe_eval.get_normal_vector(0) * - fe_eval.inverse_jacobian(0))[dim - 1]) * + std::abs( + (fe_eval.normal_vector(0) * fe_eval.inverse_jacobian(0))[dim - 1]) * number(std::max(actual_degree, 1) * (actual_degree + 1.0)) * 2.0; for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) @@ -598,7 +598,7 @@ private: value_type u_minus = phi_m.get_value(q), u_plus = phi_p.get_value(q); const VectorizedArrayType normal_times_advection = - advection * phi_m.get_normal_vector(q); + advection * phi_m.normal_vector(q); const value_type flux_times_normal = make_vectorized_array(0.5) * ((u_minus + u_plus) * normal_times_advection + @@ -646,7 +646,7 @@ private: { value_type u_minus = fe_eval.get_value(q); const VectorizedArrayType normal_times_advection = - advection * fe_eval.get_normal_vector(q); + advection * fe_eval.normal_vector(q); const value_type flux_times_normal = make_vectorized_array(0.5) * ((u_minus + u_plus) * normal_times_advection + diff --git a/tests/matrix_free/multigrid_dg_periodic.cc b/tests/matrix_free/multigrid_dg_periodic.cc index 09606de214..bf77ee538f 100644 --- a/tests/matrix_free/multigrid_dg_periodic.cc +++ b/tests/matrix_free/multigrid_dg_periodic.cc @@ -208,9 +208,9 @@ private: fe_eval_neighbor.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); VectorizedArray sigmaF = - (std::abs((fe_eval.get_normal_vector(0) * + (std::abs((fe_eval.normal_vector(0) * fe_eval.inverse_jacobian(0))[dim - 1]) + - std::abs((fe_eval.get_normal_vector(0) * + std::abs((fe_eval.normal_vector(0) * fe_eval_neighbor.inverse_jacobian(0))[dim - 1])) * (number)(std::max(fe_degree, 1) * (fe_degree + 1.0)); @@ -250,8 +250,8 @@ private: fe_eval.read_dof_values(src); fe_eval.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); VectorizedArray sigmaF = - std::abs((fe_eval.get_normal_vector(0) * - fe_eval.inverse_jacobian(0))[dim - 1]) * + std::abs( + (fe_eval.normal_vector(0) * fe_eval.inverse_jacobian(0))[dim - 1]) * (number)(std::max(1, fe_degree) * (fe_degree + 1.0)) * 2.; for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) @@ -339,10 +339,9 @@ private: phi_outer.reinit(face); VectorizedArray sigmaF = - (std::abs( - (phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]) + - std::abs((phi.get_normal_vector(0) * - phi_outer.inverse_jacobian(0))[dim - 1])) * + (std::abs((phi.normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]) + + std::abs( + (phi.normal_vector(0) * phi_outer.inverse_jacobian(0))[dim - 1])) * (number)(std::max(fe_degree, 1) * (fe_degree + 1.0)); // Compute phi part @@ -426,8 +425,7 @@ private: phi.reinit(face); VectorizedArray sigmaF = - std::abs( - (phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]) * + std::abs((phi.normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]) * (number)(std::max(1, fe_degree) * (fe_degree + 1.0)) * 2.; for (unsigned int i = 0; i < phi.dofs_per_cell; ++i) diff --git a/tests/matrix_free/multigrid_dg_sip_01.cc b/tests/matrix_free/multigrid_dg_sip_01.cc index 4b813c95b7..c3d7fdf8d1 100644 --- a/tests/matrix_free/multigrid_dg_sip_01.cc +++ b/tests/matrix_free/multigrid_dg_sip_01.cc @@ -217,9 +217,9 @@ private: fe_eval_neighbor.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); VectorizedArray sigmaF = - (std::abs((fe_eval.get_normal_vector(0) * + (std::abs((fe_eval.normal_vector(0) * fe_eval.inverse_jacobian(0))[dim - 1]) + - std::abs((fe_eval.get_normal_vector(0) * + std::abs((fe_eval.normal_vector(0) * fe_eval_neighbor.inverse_jacobian(0))[dim - 1])) * (number)(std::max(fe_degree, 1) * (fe_degree + 1.0)); @@ -259,8 +259,8 @@ private: fe_eval.read_dof_values(src); fe_eval.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); VectorizedArray sigmaF = - std::abs((fe_eval.get_normal_vector(0) * - fe_eval.inverse_jacobian(0))[dim - 1]) * + std::abs( + (fe_eval.normal_vector(0) * fe_eval.inverse_jacobian(0))[dim - 1]) * number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0; for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) @@ -346,10 +346,9 @@ private: phi_outer.reinit(face); VectorizedArray sigmaF = - (std::abs( - (phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]) + - std::abs((phi.get_normal_vector(0) * - phi_outer.inverse_jacobian(0))[dim - 1])) * + (std::abs((phi.normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]) + + std::abs( + (phi.normal_vector(0) * phi_outer.inverse_jacobian(0))[dim - 1])) * (number)(std::max(fe_degree, 1) * (fe_degree + 1.0)); // Compute phi part @@ -433,8 +432,7 @@ private: phi.reinit(face); VectorizedArray sigmaF = - std::abs( - (phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]) * + std::abs((phi.normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]) * number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0; for (unsigned int i = 0; i < phi.dofs_per_cell; ++i) diff --git a/tests/matrix_free/multigrid_dg_sip_02.cc b/tests/matrix_free/multigrid_dg_sip_02.cc index 609cc72576..e407917846 100644 --- a/tests/matrix_free/multigrid_dg_sip_02.cc +++ b/tests/matrix_free/multigrid_dg_sip_02.cc @@ -220,9 +220,9 @@ private: fe_eval_neighbor.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); VectorizedArray sigmaF = - (std::abs((fe_eval.get_normal_vector(0) * + (std::abs((fe_eval.normal_vector(0) * fe_eval.inverse_jacobian(0))[dim - 1]) + - std::abs((fe_eval.get_normal_vector(0) * + std::abs((fe_eval.normal_vector(0) * fe_eval_neighbor.inverse_jacobian(0))[dim - 1])) * (number)(std::max(fe_degree, 1) * (fe_degree + 1.0)); @@ -262,8 +262,8 @@ private: fe_eval.read_dof_values(src); fe_eval.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); VectorizedArray sigmaF = - std::abs((fe_eval.get_normal_vector(0) * - fe_eval.inverse_jacobian(0))[dim - 1]) * + std::abs( + (fe_eval.normal_vector(0) * fe_eval.inverse_jacobian(0))[dim - 1]) * (number)(std::max(1, fe_degree) * (fe_degree + 1.0)) * 2.; for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) @@ -328,8 +328,8 @@ private: { phif.reinit(cell, face); VectorizedArray sigmaF = - std::abs((phif.get_normal_vector(0) * - phif.inverse_jacobian(0))[dim - 1]) * + std::abs( + (phif.normal_vector(0) * phif.inverse_jacobian(0))[dim - 1]) * (number)(std::max(1, fe_degree) * (fe_degree + 1.0)) * 2.; std::array::size()> boundary_ids = data.get_faces_by_cells_boundary_id(cell, face); diff --git a/tests/performance/timing_navier_stokes.cc b/tests/performance/timing_navier_stokes.cc index 4480b17b2c..85ab590043 100644 --- a/tests/performance/timing_navier_stokes.cc +++ b/tests/performance/timing_navier_stokes.cc @@ -845,9 +845,9 @@ namespace NavierStokes_DG EvaluationFlags::gradients); const auto tau_ip = - (std::abs((phi_m.get_normal_vector(0) * + (std::abs((phi_m.normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) + - std::abs((phi_p.get_normal_vector(0) * + std::abs((phi_p.normal_vector(0) * phi_p.inverse_jacobian(0))[dim - 1])) * Number(viscosity * (degree + 1) * (degree + 1)); @@ -855,7 +855,7 @@ namespace NavierStokes_DG { const auto w_m = phi_m.get_value(q); const auto w_p = phi_p.get_value(q); - const auto normal = phi_m.get_normal_vector(q); + const auto normal = phi_m.normal_vector(q); auto numerical_flux = -euler_numerical_flux(w_m, w_p, normal); const auto grad_w_m = phi_m.get_gradient(q); @@ -882,14 +882,14 @@ namespace NavierStokes_DG else { const auto tau_ip = - std::abs((phi_m.get_normal_vector(0) * + std::abs((phi_m.normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) * Number(2. * viscosity * (degree + 1) * (degree + 1)); for (unsigned int q = 0; q < phi_m.n_q_points; ++q) { const auto w_m = phi_m.get_value(q); - const auto normal = phi_m.get_normal_vector(q); + const auto normal = phi_m.normal_vector(q); const auto grad_w_m = phi_m.get_gradient(q); const auto grad_w_p = grad_w_m; @@ -1225,9 +1225,9 @@ namespace NavierStokes_DG EvaluationFlags::values | EvaluationFlags::gradients); - const auto tau_ip = (std::abs((phi_m.get_normal_vector(0) * + const auto tau_ip = (std::abs((phi_m.normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) + - std::abs((phi_p.get_normal_vector(0) * + std::abs((phi_p.normal_vector(0) * phi_p.inverse_jacobian(0))[dim - 1])) * Number(viscosity * (degree + 1) * (degree + 1)); @@ -1235,7 +1235,7 @@ namespace NavierStokes_DG { const auto w_m = phi_m.get_value(q); const auto w_p = phi_p.get_value(q); - const auto normal = phi_m.get_normal_vector(q); + const auto normal = phi_m.normal_vector(q); auto numerical_flux = -euler_numerical_flux(w_m, w_p, normal); const auto grad_w_m = phi_m.get_gradient(q); const auto grad_w_p = phi_p.get_gradient(q); @@ -1294,7 +1294,7 @@ namespace NavierStokes_DG const auto tau_ip = std::abs( - (phi_m.get_normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) * + (phi_m.normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) * Number(2. * viscosity * (degree + 1) * (degree + 1)); const auto boundary_id = data.get_boundary_id(face); @@ -1302,7 +1302,7 @@ namespace NavierStokes_DG for (unsigned int q = 0; q < phi_m.n_q_points; ++q) { const auto w_m = phi_m.get_value(q); - const auto normal = phi_m.get_normal_vector(q); + const auto normal = phi_m.normal_vector(q); const auto grad_w_m = phi_m.get_gradient(q); const auto grad_w_p = grad_w_m; diff --git a/tests/simplex/matrix_free_03.cc b/tests/simplex/matrix_free_03.cc index 3caf044413..eae9681809 100644 --- a/tests/simplex/matrix_free_03.cc +++ b/tests/simplex/matrix_free_03.cc @@ -187,9 +187,9 @@ public: EvaluationFlags::values | EvaluationFlags::gradients); VectorizedArray sigmaF = PENALTY; - // (std::abs((fe_eval.get_normal_vector(0) * + // (std::abs((fe_eval.normal_vector(0) * // fe_eval.inverse_jacobian(0))[dim - 1]) + - // std::abs((fe_eval.get_normal_vector(0) * + // std::abs((fe_eval.normal_vector(0) * // fe_eval_neighbor.inverse_jacobian(0))[dim - 1])) * // (number)(std::max(fe_degree, 1) * (fe_degree + 1.0)); @@ -225,7 +225,7 @@ public: fe_eval.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); VectorizedArray sigmaF = PENALTY; - // std::abs((fe_eval.get_normal_vector(0) * + // std::abs((fe_eval.normal_vector(0) * // fe_eval.inverse_jacobian(0))[dim - 1]) * // number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0; diff --git a/tests/simplex/step-12c.cc b/tests/simplex/step-12c.cc index 4a8ba815d0..0418da687f 100644 --- a/tests/simplex/step-12c.cc +++ b/tests/simplex/step-12c.cc @@ -763,8 +763,8 @@ public: phi.reinit(face); for (unsigned int q = 0; q < phi.n_q_points; ++q) { - const auto beta_n = this->beta(phi.quadrature_point(q)) * - phi.get_normal_vector(q); + const auto beta_n = + this->beta(phi.quadrature_point(q)) * phi.normal_vector(q); const auto beta_n_m = (-std::abs(beta_n) + beta_n) * 0.5; phi.submit_value(-beta_n_m * this->boundary_values( phi.quadrature_point(q)), @@ -810,7 +810,7 @@ public: for (unsigned int q = 0; q < phi_m.n_q_points; ++q) { const auto beta_n = this->beta(phi_m.quadrature_point(q)) * - phi_m.get_normal_vector(q); + phi_m.normal_vector(q); const auto u_m = phi_m.get_value(q); const auto u_p = phi_p.get_value(q); @@ -834,8 +834,8 @@ public: phi.gather_evaluate(src, true, false); for (unsigned int q = 0; q < phi.n_q_points; ++q) { - const auto beta_n = this->beta(phi.quadrature_point(q)) * - phi.get_normal_vector(q); + const auto beta_n = + this->beta(phi.quadrature_point(q)) * phi.normal_vector(q); const auto beta_n_p = (std::abs(beta_n) + beta_n) * 0.5; phi.submit_value(beta_n_p * phi.get_value(q), q); } diff --git a/tests/simplex/step-67.cc b/tests/simplex/step-67.cc index 4c13a9d894..52bb6b577b 100644 --- a/tests/simplex/step-67.cc +++ b/tests/simplex/step-67.cc @@ -1199,7 +1199,7 @@ namespace Euler_DG const auto numerical_flux = euler_numerical_flux(phi_m.get_value(q), phi_p.get_value(q), - phi_m.get_normal_vector(q)); + phi_m.normal_vector(q)); // std::cout<<" phim: "<