Number * output,
const bool do_gradients,
const unsigned int face_no)
+ {
+ Assert(static_cast<unsigned int>(fe_degree) ==
+ data.data.front().fe_degree ||
+ fe_degree == -1,
+ ExcInternalError());
+
+ interpolate_generic<do_evaluate, add_into_output>(
+ input,
+ output,
+ do_gradients,
+ face_no,
+ data.data.front().fe_degree + 1,
+ data.data.front().shape_data_on_face,
+ data.dofs_per_component_on_cell,
+ 2 * data.dofs_per_component_on_face);
+ }
+
+ /**
+ * Interpolate the values on the cell quadrature points onto a face.
+ */
+ template <bool do_evaluate, bool add_into_output>
+ static void
+ interpolate_quadrature(const MatrixFreeFunctions::ShapeInfo<Number> &data,
+ const Number * input,
+ Number * output,
+ const bool do_gradients,
+ const unsigned int face_no)
+ {
+ Assert(static_cast<unsigned int>(fe_degree + 1) ==
+ data.data.front().quadrature.size() ||
+ fe_degree == -1,
+ ExcInternalError());
+
+ interpolate_generic<do_evaluate, add_into_output>(
+ input,
+ output,
+ do_gradients,
+ face_no,
+ data.data.front().quadrature.size(),
+ data.data.front().quadrature_data_on_face,
+ data.n_q_points,
+ data.n_q_points_face);
+ }
+
+ private:
+ template <bool do_evaluate, bool add_into_output>
+ static void
+ interpolate_generic(const Number * input,
+ Number * output,
+ const bool do_gradients,
+ const unsigned int face_no,
+ const unsigned int n_points_1d,
+ const AlignedVector<Number> *shape_data,
+ const unsigned int dofs_per_component_on_cell,
+ const unsigned int dofs_per_component_on_face)
{
internal::EvaluatorTensorProduct<internal::evaluate_general,
dim,
fe_degree + 1,
0,
Number>
- evalf(data.data.front().shape_data_on_face[face_no % 2],
+ evalf(shape_data[face_no % 2],
AlignedVector<Number>(),
AlignedVector<Number>(),
- data.data.front().fe_degree + 1,
+ n_points_1d,
0);
- const unsigned int in_stride = do_evaluate ?
- data.dofs_per_component_on_cell :
- 2 * data.dofs_per_component_on_face;
- const unsigned int out_stride = do_evaluate ?
- 2 * data.dofs_per_component_on_face :
- data.dofs_per_component_on_cell;
+ const unsigned int in_stride =
+ do_evaluate ? dofs_per_component_on_cell : dofs_per_component_on_face;
+ const unsigned int out_stride =
+ do_evaluate ? dofs_per_component_on_face : dofs_per_component_on_cell;
const unsigned int face_direction = face_no / 2;
for (unsigned int c = 0; c < n_components; c++)
{
univariate_shape_data.shape_gradients_collocation;
auto &shape_hessians_collocation =
univariate_shape_data.shape_hessians_collocation;
- auto &inverse_shape_values = univariate_shape_data.inverse_shape_values;
- auto &shape_data_on_face = univariate_shape_data.shape_data_on_face;
+ auto &inverse_shape_values = univariate_shape_data.inverse_shape_values;
+ auto &shape_data_on_face = univariate_shape_data.shape_data_on_face;
+ auto &quadrature_data_on_face =
+ univariate_shape_data.quadrature_data_on_face;
auto &values_within_subface = univariate_shape_data.values_within_subface;
auto &gradients_within_subface =
univariate_shape_data.gradients_within_subface;
fe->shape_grad_grad(my_i, q_point)[0][0];
}
+ if (n_q_points_1d < 200)
+ {
+ quadrature_data_on_face[0].resize(quad.size() * 3);
+ quadrature_data_on_face[1].resize(quad.size() * 3);
+
+ dealii::FE_DGQArbitraryNodes<1> fe_quad(quad);
+
+ for (unsigned int i = 0; i < quad.size(); ++i)
+ {
+ Point<1> q_point;
+ q_point[0] = 0;
+ quadrature_data_on_face[0][i] = fe_quad.shape_value(i, q_point);
+ q_point[0] = 1;
+ quadrature_data_on_face[1][i] = fe_quad.shape_value(i, q_point);
+ }
+ }
+
// get gradient and Hessian transformation matrix for the polynomial
// space associated with the quadrature rule (collocation space). We
// need to avoid the case with more than a few hundreds of quadrature
{
memory +=
MemoryConsumption::memory_consumption(shape_data_on_face[i]);
+ memory +=
+ MemoryConsumption::memory_consumption(quadrature_data_on_face[i]);
memory +=
MemoryConsumption::memory_consumption(values_within_subface[i]);
memory +=
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+
+// Tests internal::FEFaceNormalEvaluationImpl::interpolate_quadrature()
+// by comparing the results of FEFaceEvaluation::gather_evaluate().
+
+template <int dim>
+class ExactSolution : public Function<dim>
+{
+public:
+ ExactSolution()
+ {}
+
+ virtual double
+ value(const Point<dim> &p, const unsigned int /*component*/ = 0) const
+ {
+ return p[0] + p[1];
+ }
+};
+
+template <int dim,
+ int fe_degree,
+ int n_points = fe_degree + 1,
+ typename Number = double,
+ typename VectorizedArrayType = VectorizedArray<Number>>
+void
+test(const unsigned int n_refinements = 1)
+{
+ using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria, 0, 1, true);
+
+ tria.refine_global(n_refinements);
+
+ FE_DGQ<dim> fe(fe_degree);
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+
+ MappingQ<dim> mapping(1);
+ QGauss<1> quad(n_points);
+
+ AffineConstraints<Number> constraint;
+
+ using MF = MatrixFree<dim, Number, VectorizedArrayType>;
+
+ typename MF::AdditionalData additional_data;
+ additional_data.mapping_update_flags = update_values;
+ additional_data.mapping_update_flags_inner_faces = update_values;
+ additional_data.mapping_update_flags_boundary_faces = update_values;
+ additional_data.mapping_update_flags_faces_by_cells = update_values;
+ additional_data.hold_all_faces_to_owned_cells = true;
+
+ MF matrix_free;
+ matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data);
+
+ VectorType src, dst;
+
+ matrix_free.initialize_dof_vector(src);
+ matrix_free.initialize_dof_vector(dst);
+
+ VectorTools::interpolate(dof_handler, ExactSolution<dim>(), src);
+
+ FEEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType> phi(
+ matrix_free);
+ FEFaceEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType>
+ phi_m(matrix_free, true);
+
+ matrix_free.template loop_cell_centric<VectorType, VectorType>(
+ [&](const auto &, auto &, const auto &src, const auto range) {
+ for (unsigned int cell = range.first; cell < range.second; ++cell)
+ {
+ phi.reinit(cell);
+
+ phi.gather_evaluate(src, true, false);
+
+ for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
+ face++)
+ {
+ phi_m.reinit(cell, face);
+ phi_m.gather_evaluate(src, true, false);
+
+ AlignedVector<VectorizedArrayType> temp(
+ phi_m.static_dofs_per_cell);
+
+ internal::FEFaceNormalEvaluationImpl<dim,
+ n_points - 1,
+ 1,
+ VectorizedArrayType>::
+ template interpolate_quadrature<true, false>(
+ matrix_free.get_shape_info(),
+ phi.begin_values(),
+ temp.data(),
+ false,
+ face);
+
+
+ for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+ {
+ const auto u_cell = temp[q];
+ const auto u_face = phi_m.get_value(q);
+
+ for (unsigned int v = 0; v < VectorizedArray<double>::size();
+ ++v)
+ {
+ Assert(std::abs(u_cell[v] - u_face[v]) < 1e-10,
+ ExcMessage("Entries do not match!"));
+ }
+ }
+ }
+ }
+ },
+ dst,
+ src);
+}
+
+int
+main()
+{
+ initlog();
+ test<2, 1, 2, double, VectorizedArray<double>>();
+ test<3, 1, 2, double, VectorizedArray<double>>();
+
+ deallog << "OK!" << std::endl;
+}
projects / dealii.git / commitdiff