#include <deal.II/matrix_free/evaluation_template_factory.h>
#include <deal.II/matrix_free/mapping_info.h>
+#include <deal.II/matrix_free/util.h>
DEAL_II_NAMESPACE_OPEN
for (unsigned int i = 0; i < n_q_points; ++i)
quadrature_weights[i] = quadrature.weight(i);
- // note: quadrature_1 and tensor_quadrature_weights are not set up
+ // note: quadrature_1d and tensor_quadrature_weights are not set up
// TODO: set up face_orientations
(void)update_flags_inner_faces;
if (flag == false)
{
#ifdef DEAL_II_WITH_SIMPLEX_SUPPORT
- continue; // TODO: face data is not set up
+ const auto quad_face = get_face_quadrature(quad[my_q][hpq]);
+ face_data[my_q].descriptor[hpq].initialize(quad_face,
+ update_default);
#else
Assert(false, ExcNotImplemented());
#endif
}
-
- face_data[my_q].descriptor[hpq].initialize(
- quad[my_q][hpq].get_tensor_basis()[0],
- update_flags_boundary_faces);
+ else
+ face_data[my_q].descriptor[hpq].initialize(
+ quad[my_q][hpq].get_tensor_basis()[0],
+ update_flags_boundary_faces);
}
face_data_by_cells[my_q].descriptor.resize(n_hp_quads);
if (flag == false)
{
#ifdef DEAL_II_WITH_SIMPLEX_SUPPORT
- continue; // TODO: face data is not set up
+ const auto quad_face = get_face_quadrature(quad[my_q][hpq]);
+ face_data_by_cells[my_q].descriptor[hpq].initialize(
+ quad_face, update_default);
#else
Assert(false, ExcNotImplemented());
#endif
}
-
- face_data_by_cells[my_q].descriptor[hpq].initialize(
- quad[my_q][hpq].get_tensor_basis()[0], update_default);
+ else
+ face_data_by_cells[my_q].descriptor[hpq].initialize(
+ quad[my_q][hpq].get_tensor_basis()[0], update_default);
}
}
#include <deal.II/matrix_free/face_setup_internal.h>
#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
#ifdef DEAL_II_WITH_TBB
# include <deal.II/base/parallel.h>
Table<2, internal::MatrixFreeFunctions::ShapeInfo<double>> shape_info_dummy(
shape_info.size(0), shape_info.size(2));
{
- QGauss<dim> quad(1);
+ Quadrature<dim> quad(QGauss<dim>(1));
+ Quadrature<dim> quad_simplex(Simplex::QGauss<dim>(1));
for (unsigned int no = 0, c = 0; no < dof_handlers.size(); no++)
for (unsigned int b = 0;
b < dof_handlers[no]->get_fe(0).n_base_elements();
for (unsigned int fe_no = 0;
fe_no < dof_handlers[no]->get_fe_collection().size();
++fe_no)
- shape_info_dummy(c, fe_no).reinit(quad,
- dof_handlers[no]->get_fe(fe_no),
- b);
+ shape_info_dummy(c, fe_no).reinit(
+ dof_handlers[no]->get_fe(fe_no).reference_cell_type() ==
+ ReferenceCell::get_hypercube(dim) ?
+ quad :
+ quad_simplex,
+ dof_handlers[no]->get_fe(fe_no),
+ b);
}
const unsigned int n_lanes = VectorizedArrayType::size();
* end points of the unit cell.
*/
bool nodal_at_cell_boundaries;
+
+ /**
+ * Stores the shape values of the finite element evaluated on all
+ * quadrature points for all faces and orientations (no tensor-product
+ * structure exploited).
+ */
+ Table<3, Number> shape_values_face;
+
+ /**
+ * Stores the shape gradients of the finite element evaluated on all
+ * quadrature points for all faces, orientations, and directions
+ * (no tensor-product structure exploited).
+ */
+ Table<4, Number> shape_gradients_face;
};
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/polynomial.h>
#include <deal.II/base/polynomials_piecewise.h>
+#include <deal.II/base/qprojector.h>
#include <deal.II/base/tensor_product_polynomials.h>
#include <deal.II/base/utilities.h>
#include <deal.II/fe/fe_poly.h>
#include <deal.II/fe/fe_q_dg0.h>
+#include <deal.II/grid/reference_cell.h>
+
#include <deal.II/lac/householder.h>
#include <deal.II/matrix_free/shape_info.h>
+#include <deal.II/matrix_free/util.h>
#include <deal.II/simplex/fe_lib.h>
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;
+ auto &shape_values = univariate_shape_data.shape_values;
+ auto &shape_values_face = univariate_shape_data.shape_values_face;
+ auto &shape_gradients = univariate_shape_data.shape_gradients;
+ auto &shape_gradients_face =
+ univariate_shape_data.shape_gradients_face;
const unsigned int n_dofs = fe->n_dofs_per_cell();
q] = grad[d];
}
+ const auto reference_cell_type = ReferenceCell::get_simplex(dim);
+
+ try
+ {
+ const auto quad_face = get_face_quadrature(quad);
+ this->n_q_points_face = quad_face.size();
+
+ const unsigned int n_face_orientations = dim == 2 ? 2 : 6;
+ const unsigned int n_faces =
+ ReferenceCell::internal::Info::get_cell(reference_cell_type)
+ .n_faces();
+
+ const auto projected_quad_face =
+ QProjector<dim>::project_to_all_faces(reference_cell_type,
+ quad_face);
+
+ shape_values_face.reinit(
+ {n_faces, n_face_orientations, n_dofs * n_q_points_face});
+
+ shape_gradients_face.reinit(
+ {n_faces, n_face_orientations, dim, n_dofs * n_q_points_face});
+
+ for (unsigned int f = 0; f < n_faces; ++f)
+ for (unsigned int o = 0; o < n_face_orientations; ++o)
+ {
+ const auto offset =
+ QProjector<dim>::DataSetDescriptor::face(
+ reference_cell_type,
+ f,
+ (o ^ 1) & 1, // face_orientation
+ (o >> 1) & 1, // face_flip
+ (o >> 2) & 1, // face_rotation
+ n_q_points_face);
+
+ for (unsigned int i = 0; i < n_dofs; ++i)
+ for (unsigned int q = 0; q < n_q_points_face; ++q)
+ {
+ const auto point =
+ projected_quad_face.point(q + offset);
+
+ shape_values_face(f, o, i * n_q_points_face + q) =
+ fe->shape_value(i, point);
+
+ const auto grad = fe->shape_grad(i, point);
+
+ for (int d = 0; d < dim; ++d)
+ shape_gradients_face(
+ f, o, d, i * n_q_points_face + q) = grad[d];
+ }
+ }
+ }
+ catch (...)
+ {
+ // TODO: this might happen if the quadrature rule and the
+ // the FE do not match
+ this->n_q_points_face = 0;
+ }
+
// TODO: also fill shape_hessians, inverse_shape_values,
// shape_data_on_face, quadrature_data_on_face,
// values_within_subface, gradients_within_subface,
--- /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.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_matrix_free_util_h
+#define dealii_matrix_free_util_h
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/quadrature.h>
+
+#include <deal.II/simplex/quadrature_lib.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace internal
+{
+ namespace MatrixFreeFunctions
+ {
+ template <int dim>
+ inline Quadrature<dim - 1>
+ get_face_quadrature(const Quadrature<dim> &quad)
+ {
+ if (dim == 2 || dim == 3)
+ for (unsigned int i = 1; i <= 3; ++i)
+ if (quad == Simplex::QGauss<dim>(i))
+ return Simplex::QGauss<dim - 1>(i);
+
+ AssertThrow(false, ExcNotImplemented());
+
+ return Quadrature<dim - 1>();
+ }
+
+ } // end of namespace MatrixFreeFunctions
+} // end of namespace internal
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
projects / dealii.git / commitdiff