* 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
// take the absolute value of these factors as the normal could point
// into either positive or negative direction.
const VectorizedArray<number> 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<number> sigma =
inverse_length_normal_to_face * get_penalty_factor();
EvaluationFlags::values |
EvaluationFlags::gradients);
- const VectorizedArray<number> inverse_length_normal_to_face =
- std::abs((phi_inner.get_normal_vector(0) *
- phi_inner.inverse_jacobian(0))[dim - 1]);
+ const VectorizedArray<number> inverse_length_normal_to_face = std::abs((
+ phi_inner.normal_vector(0) * phi_inner.inverse_jacobian(0))[dim - 1]);
const VectorizedArray<number> sigma =
inverse_length_normal_to_face * get_penalty_factor();
{
phi_face.reinit(face);
- const VectorizedArray<double> inverse_length_normal_to_face =
- std::abs((phi_face.get_normal_vector(0) *
- phi_face.inverse_jacobian(0))[dim - 1]);
+ const VectorizedArray<double> inverse_length_normal_to_face = std::abs(
+ (phi_face.normal_vector(0) * phi_face.inverse_jacobian(0))[dim - 1]);
const VectorizedArray<double> sigma =
inverse_length_normal_to_face * system_matrix.get_penalty_factor();
{
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);
}
const auto numerical_flux =
euler_numerical_flux<dim>(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);
}
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)
const auto numerical_flux =
euler_numerical_flux<dim>(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);
}
}
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)
/**
* 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<unsigned int, vectorization_width> cells_interior;
* 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
* 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
* @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;
/** @} */
template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number>
-FEEvaluationData<dim, Number, is_face>::get_normal_vector(
+FEEvaluationData<dim, Number, is_face>::normal_vector(
const unsigned int q_point) const
{
AssertIndexRange(q_point, n_quadrature_points);
+// This function is deprecated.
+template <int dim, typename Number, bool is_face>
+inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number>
+FEEvaluationData<dim, Number, is_face>::get_normal_vector(
+ const unsigned int q_point) const
+{
+ return normal_vector(q_point);
+}
+
+
+
template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE Number
FEEvaluationData<dim, Number, is_face>::JxW(const unsigned int q_point) const
value_type u_minus = phi_m.get_value(q),
u_plus = phi_p.get_value(q);
const VectorizedArray<number> 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<number>(0.5) *
((u_minus + u_plus) * normal_times_advection +
value_type u_minus = fe_eval.get_value(q);
value_type u_plus = -u_minus;
const VectorizedArray<number> 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<number>(0.5) *
((u_minus + u_plus) * normal_times_advection +
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_dgq.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/operators.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+
+template <int dim,
+ int fe_degree,
+ int n_q_points_1d = fe_degree + 1,
+ typename number = double,
+ typename VectorType = Vector<number>,
+ int n_components = 1>
+class MatrixFreeAdvectionBasic
+{
+public:
+ MatrixFreeAdvectionBasic(const MatrixFree<dim, number> &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<dim, number>::DataAccessOnFaces::values,
+ MatrixFree<dim, number>::DataAccessOnFaces::values);
+
+ FEEvaluation<dim, fe_degree, fe_degree + 1, n_components, number> phi(data);
+
+ const unsigned int dofs_per_cell = phi.dofs_per_cell;
+
+ AlignedVector<VectorizedArray<number>> coefficients(phi.dofs_per_cell);
+ MatrixFreeOperators::
+ CellwiseInverseMassMatrix<dim, fe_degree, n_components, number>
+ 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<dim, number> & data,
+ VectorType & dst,
+ const VectorType & src,
+ const std::pair<unsigned int, unsigned int> &cell_range) const
+ {
+ FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, number> 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<dim, number> & data,
+ VectorType & dst,
+ const VectorType & src,
+ const std::pair<unsigned int, unsigned int> &face_range) const
+ {
+ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components, number> phi_m(
+ data, true, 0, 0, start_vector_component);
+ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components, number> phi_p(
+ data, false, 0, 0, start_vector_component);
+ using value_type = typename FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ number>::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<number> normal_times_advection =
+ advection * phi_m.get_normal_vector(q);
+ const value_type flux_times_normal =
+ make_vectorized_array<number>(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<dim, number> & data,
+ VectorType & dst,
+ const VectorType & src,
+ const std::pair<unsigned int, unsigned int> &face_range) const
+ {
+ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components, number>
+ fe_eval(data, true, 0, 0, start_vector_component);
+ using value_type = typename FEFaceEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ number>::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<number> normal_times_advection =
+ advection * fe_eval.get_normal_vector(q);
+ const value_type flux_times_normal =
+ make_vectorized_array<number>(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<dim, number> & data;
+ const bool zero_within_loop;
+ const unsigned int start_vector_component;
+ Tensor<1, dim, VectorizedArray<number>> advection;
+};
+
+
+
+template <int dim>
+class AnalyticFunction : public Function<dim>
+{
+public:
+ static_assert(dim == 1, "Only 1D implemented");
+ AnalyticFunction()
+ : Function<dim>(1)
+ {}
+
+ virtual double
+ value(const Point<dim> &p, const unsigned int) const override
+ {
+ return std::sin(3 * numbers::PI * p[0] / 0.8);
+ }
+};
+
+
+
+template <int dim>
+class AnalyticDerivative : public Function<dim>
+{
+public:
+ static_assert(dim == 1, "Only 1D implemented");
+ AnalyticDerivative()
+ : Function<dim>(1)
+ {}
+
+ virtual double
+ value(const Point<dim> &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 <int dim, int fe_degree>
+void
+test(const unsigned int n_refine)
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria, 0, 0.8);
+ tria.refine_global(n_refine);
+
+ FE_DGQ<dim> fe(fe_degree);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+ AffineConstraints<double> constraints;
+ constraints.close();
+
+ if (n_refine == 3)
+ {
+ deallog << "Testing " << dof.get_fe().get_name();
+ deallog << std::endl;
+ }
+
+ LinearAlgebra::distributed::Vector<double> in, out;
+
+ const QGauss<1> quad(fe_degree + 1);
+ typename MatrixFree<dim, double>::AdditionalData data;
+ data.tasks_parallel_scheme = MatrixFree<dim, double>::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<dim, double> mf_data;
+ mf_data.reinit(MappingQ1<dim>{}, dof, constraints, quad, data);
+
+ mf_data.initialize_dof_vector(in);
+ mf_data.initialize_dof_vector(out);
+
+ VectorTools::interpolate(dof, AnalyticFunction<dim>(), in);
+
+ MatrixFreeAdvectionBasic<dim,
+ fe_degree,
+ fe_degree + 1,
+ double,
+ LinearAlgebra::distributed::Vector<double>>
+ mf2(mf_data);
+ mf2.vmult(out, in);
+
+ VectorTools::interpolate(dof, AnalyticDerivative<dim>(), 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);
+}
--- /dev/null
+
+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
value_type u_minus = phi_m.get_value(q),
u_plus = phi_p.get_value(q);
const VectorizedArray<number> 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<number>(0.5) *
((u_minus + u_plus) * normal_times_advection +
value_type u_minus = fe_eval.get_value(q);
value_type u_plus = -u_minus;
const VectorizedArray<number> 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<number>(0.5) *
((u_minus + u_plus) * normal_times_advection +
value_type u_minus = phi_m.get_value(q),
u_plus = phi_p.get_value(q);
const VectorizedArray<number> 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<number>(0.5) *
((u_minus + u_plus) * normal_times_advection +
value_type u_minus = fe_eval.get_value(q);
value_type u_plus = -u_minus;
const VectorizedArray<number> 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<number>(0.5) *
((u_minus + u_plus) * normal_times_advection +
// penalty parameter because the cell-based method cannot read the
// sigma parameter on the neighbor
VectorizedArray<number> 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
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)
VectorizedArray<number> local_diagonal_vector[phi.static_dofs_per_cell];
VectorizedArray<number> 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)
{
phif.reinit(cell, face);
VectorizedArray<number> 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<types::boundary_id, VectorizedArray<number>::size()>
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)
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)
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;
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));
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;
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);
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);
}
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<number, VectorizedArrayType::size()>(0.5) *
((u_minus + u_plus) * normal_times_advection +
{
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<number, VectorizedArrayType::size()>(0.5) *
((u_minus + u_plus) * normal_times_advection +
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));
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)
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));
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)
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<number, VectorizedArrayType::size()>(0.5) *
((u_minus + u_plus) * normal_times_advection +
{
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<number, VectorizedArrayType::size()>(0.5) *
((u_minus + u_plus) * normal_times_advection +
fe_eval_neighbor.evaluate(EvaluationFlags::values |
EvaluationFlags::gradients);
VectorizedArray<number> 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));
fe_eval.read_dof_values(src);
fe_eval.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
VectorizedArray<number> 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)
phi_outer.reinit(face);
VectorizedArray<number> 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
phi.reinit(face);
VectorizedArray<number> 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)
fe_eval_neighbor.evaluate(EvaluationFlags::values |
EvaluationFlags::gradients);
VectorizedArray<number> 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));
fe_eval.read_dof_values(src);
fe_eval.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
VectorizedArray<number> 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)
phi_outer.reinit(face);
VectorizedArray<number> 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
phi.reinit(face);
VectorizedArray<number> 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)
fe_eval_neighbor.evaluate(EvaluationFlags::values |
EvaluationFlags::gradients);
VectorizedArray<number> 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));
fe_eval.read_dof_values(src);
fe_eval.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
VectorizedArray<number> 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)
{
phif.reinit(cell, face);
VectorizedArray<number> 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<types::boundary_id, VectorizedArray<number>::size()>
boundary_ids = data.get_faces_by_cells_boundary_id(cell, face);
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));
{
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<dim>(w_m, w_p, normal);
const auto grad_w_m = phi_m.get_gradient(q);
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;
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));
{
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<dim>(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);
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);
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;
EvaluationFlags::values |
EvaluationFlags::gradients);
VectorizedArray<number> 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));
fe_eval.evaluate(EvaluationFlags::values |
EvaluationFlags::gradients);
VectorizedArray<number> 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;
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)),
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);
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);
}
const auto numerical_flux =
euler_numerical_flux<dim>(phi_m.get_value(q),
phi_p.get_value(q),
- phi_m.get_normal_vector(q));
+ phi_m.normal_vector(q));
// std::cout<<" phim: "<<phi_m.get_value(q)
// <<" phip: "<<phi_p.get_value(q)<<std::endl;
phi_m.submit_value(-numerical_flux, q);
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)
projects / dealii.git / commitdiff