#include <deal.II/matrix_free/shape_info.h>
+#include <deal.II/simplex/fe_lib.h>
+
DEAL_II_NAMESPACE_OPEN
const FiniteElement<dim> &fe_in,
const unsigned int base_element_number)
{
+#ifdef DEAL_II_WITH_SIMPLEX_SUPPORT
+ if (quad_in.is_tensor_product() == false ||
+ dynamic_cast<const Simplex::FE_P<dim> *>(&fe_in) ||
+ dynamic_cast<const Simplex::FE_DGP<dim> *>(&fe_in))
+ {
+ // specialization for arbitrary finite elements and quadrature rules
+ // as needed in the context, e.g., of simplices
+
+ AssertDimension(dim, dim_q);
+
+ const auto quad = Quadrature<dim>(quad_in);
+ const auto fe = &fe_in.base_element(base_element_number);
+ n_dimensions = dim;
+ n_components = fe_in.n_components();
+ n_q_points = quad.size();
+ dofs_per_component_on_cell = fe->n_dofs_per_cell();
+ n_q_points_face = 0; // not implemented yet
+ dofs_per_component_on_face = 0; //
+
+ Assert(fe->n_components() == 1,
+ ExcMessage(
+ "FEEvaluation only works for scalar finite elements."));
+
+ data.resize(1);
+ UnivariateShapeData<Number> &univariate_shape_data = data.front();
+ data_access.reinit(n_dimensions, n_components);
+ data_access.fill(&univariate_shape_data);
+
+ // note: we cannot write `univariate_shape_data.quadrature = quad`,
+ // since the quadrature rule within UnivariateShapeData expects
+ // a 1D quadrature rule. However, in this case we are not able to
+ // define that rule anyway so other code cannot use this information.
+
+ univariate_shape_data.fe_degree = fe->degree;
+ univariate_shape_data.n_q_points_1d = quad.size();
+
+ if ((fe->n_dofs_per_cell() == 0) || (quad.size() == 0))
+ return;
+
+ // grant write access to common univariate shape data
+ auto &shape_values = univariate_shape_data.shape_values;
+ auto &shape_gradients = univariate_shape_data.shape_gradients;
+
+ const unsigned int n_dofs = fe->n_dofs_per_cell();
+
+ const unsigned int array_size = n_dofs * n_q_points;
+
+ shape_values.resize_fast(array_size);
+ shape_gradients.resize_fast(array_size * dim);
+
+ for (unsigned int i = 0; i < n_dofs; ++i)
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ {
+ shape_values[i * n_q_points + q] =
+ fe->shape_value(i, quad.point(q));
+
+ const auto grad = fe->shape_grad(i, quad.point(q));
+
+ for (int d = 0; d < dim; ++d)
+ shape_gradients[d * n_dofs * n_q_points + i * n_q_points +
+ q] = grad[d];
+ }
+
+ // TODO: also fill shape_hessians, inverse_shape_values,
+ // shape_data_on_face, quadrature_data_on_face,
+ // values_within_subface, gradients_within_subface,
+ // hessians_within_subface
+
+ // note: shape_gradients_collocation, shape_hessians_collocation,
+ // shape_values_eo, shape_gradients_eo, shape_hessians_eo,
+ // shape_gradients_collocation_eo, shape_hessians_collocation_eo,
+ // inverse_shape_values_eo cannot be filled
+
+ // indicate that no tensor product properties could be exploited
+ element_type = tensor_none;
+ univariate_shape_data.element_type = this->element_type;
+
+ univariate_shape_data.nodal_at_cell_boundaries = true;
+
+ lexicographic_numbering.resize(n_dofs);
+
+ for (unsigned int i = 0; i < n_dofs; ++i)
+ lexicographic_numbering[i] = i;
+
+ // TODO: setup face_to_cell_index_nodal, face_to_cell_index_hermite,
+ // face_orientations
+
+ return;
+ }
+#else
Assert(quad_in.is_tensor_product(), ExcNotImplemented());
+#endif
+
const auto quad = quad_in.get_tensor_basis()[0];
const FiniteElement<dim> *fe = &fe_in.base_element(base_element_number);
--- /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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test ShapeData for Simplex::FE_P and Simplex::QGauss
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+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] + 2 * p[1];
+ }
+};
+
+template <int dim, typename Number>
+class FEEvaluationDummy
+{
+public:
+ FEEvaluationDummy(
+ const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info)
+ : shape_info(shape_info)
+ {
+ values.resize(shape_info.n_q_points);
+ gradients.resize(shape_info.n_q_points * dim);
+ }
+
+ void
+ evaluate(const Vector<Number> &dof_values)
+ {
+ const unsigned int n_dofs = shape_info.dofs_per_component_on_cell;
+ const unsigned int n_q_points = shape_info.n_q_points;
+
+ for (unsigned int q = 0, c = 0; q < shape_info.n_q_points; ++q)
+ {
+ values[q] = 0.0;
+ for (unsigned int i = 0; i < n_dofs; ++i, ++c)
+ values[q] += shape_info.data[0].shape_values[c] * dof_values[i];
+ }
+
+ for (unsigned int d = 0; d < dim; ++d)
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ {
+ gradients[n_q_points * d + q] = 0.0;
+ for (unsigned int i = 0; i < n_dofs; ++i)
+ gradients[n_q_points * d + q] +=
+ shape_info.data[0]
+ .shape_gradients[d * n_dofs * n_q_points + i * n_q_points + q] *
+ dof_values[i];
+ }
+ }
+
+ Number
+ get_value(const unsigned int q_point) const
+ {
+ return values[q_point];
+ }
+
+ Tensor<1, dim, Number>
+ get_gradient(const unsigned int q_point) const
+ {
+ Tensor<1, dim, Number> out;
+
+ for (unsigned int c = 0; c < dim; ++c)
+ out[c] = gradients[shape_info.n_q_points * c + q_point];
+
+ return out;
+ }
+
+private:
+ const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info;
+
+ AlignedVector<Number> values;
+ AlignedVector<Number> gradients;
+};
+
+template <int dim, int spacedim = dim, typename Number = double>
+void
+test(const FiniteElement<dim, spacedim> &fe)
+{
+ Triangulation<dim, spacedim> tria;
+ GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
+
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+
+ MappingFE<dim> mapping(Simplex::FE_P<dim>(1));
+
+ Simplex::QGauss<dim> quadrature(1);
+
+ internal::MatrixFreeFunctions::ShapeInfo<Number> shape_info(quadrature, fe);
+
+ FEValues<dim> fe_values(mapping,
+ fe,
+ quadrature,
+ update_values | update_gradients |
+ update_inverse_jacobians);
+
+ Vector<Number> src(dof_handler.n_dofs());
+ VectorTools::interpolate(mapping, dof_handler, ExactSolution<dim>(), src);
+
+ FEEvaluationDummy<dim, Number> fe_eval(shape_info);
+
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ {
+ fe_values.reinit(cell);
+
+ Vector<Number> dof_values(fe_values.dofs_per_cell);
+ cell->get_dof_values(src, dof_values);
+ fe_eval.evaluate(dof_values);
+
+ // check values
+ {
+ std::vector<Number> values(fe_values.n_quadrature_points);
+ fe_values.get_function_values(src, values);
+
+ for (const auto q : fe_values.quadrature_point_indices())
+ {
+ Assert(std::abs(values[q] - fe_eval.get_value(q)) < 1e-10,
+ ExcMessage("Entries do not match!"));
+ }
+ }
+
+ // check gradients
+ {
+ std::vector<Tensor<1, dim, Number>> gradients(
+ fe_values.n_quadrature_points);
+ fe_values.get_function_gradients(src, gradients);
+
+ for (const auto q : fe_values.quadrature_point_indices())
+ {
+ // apply the transposed inverse Jacobian of the mapping
+ Tensor<1, dim, Number> temp = fe_eval.get_gradient(q);
+ Tensor<1, dim, Number> temp_transformed;
+ for (int d = 0; d < dim; ++d)
+ {
+ Number sum = 0;
+ for (int e = 0; e < dim; ++e)
+ sum += fe_values.inverse_jacobian(q)[e][d] * temp[e];
+ temp_transformed[d] = sum;
+ }
+
+ for (int i = 0; i < dim; ++i)
+ Assert(std::abs(gradients[q][i] - temp_transformed[i]) < 1e-10,
+ ExcMessage("Entries do not match!"));
+ }
+ }
+ }
+}
+
+int
+main()
+{
+ initlog();
+
+ test<2>(Simplex::FE_P<2>(2));
+}
projects / dealii.git / commitdiff