std::array<unsigned int, VectorizedArrayType::size()>
get_cell_or_face_ids() const;
-
/**
* Return the numbering of local degrees of freedom within the evaluation
* routines of FEEvaluation in terms of the standard numbering on finite
#include <deal.II/base/aligned_vector.h>
#include <deal.II/base/exceptions.h>
-#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/table.h>
#include <deal.II/base/vectorization.h>
-#include <deal.II/fe/fe.h>
-
DEAL_II_NAMESPACE_OPEN
+
+// forward declaration
+template <int dim, int spacedim>
+class FiniteElement;
+
+
namespace internal
{
namespace MatrixFreeFunctions
ElementType element_type;
/**
- * Stores the shape values of the 1D finite element evaluated on all 1D
+ * Stores the shape values of the 1D finite element evaluated at all 1D
* quadrature points. The length of
* this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
* points are the index running fastest.
AlignedVector<Number> shape_values;
/**
- * Stores the shape gradients of the 1D finite element evaluated on all
+ * Stores the shape gradients of the 1D finite element evaluated at all
* 1D quadrature points. The length of
* this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
* points are the index running fastest.
AlignedVector<Number> shape_gradients;
/**
- * Stores the shape Hessians of the 1D finite element evaluated on all
+ * Stores the shape Hessians of the 1D finite element evaluated at all
* 1D quadrature points. The length of
* this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
* points are the index running fastest.
bool nodal_at_cell_boundaries;
/**
- * Stores the shape values of the finite element evaluated on all
+ * Stores the shape values of the finite element evaluated at 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
+ * Stores the shape gradients of the finite element evaluated at all
* quadrature points for all faces, orientations, and directions
- * (no tensor-product structure exploited).
+ * (no tensor-product structure exploited).
*/
Table<4, Number> shape_gradients_face;
};
/**
* Constructor that initializes the data fields using the reinit method.
*/
- template <int dim, int dim_q>
- ShapeInfo(const Quadrature<dim_q> & quad,
- const FiniteElement<dim> &fe,
- const unsigned int base_element = 0);
+ template <int dim, int spacedim, int dim_q>
+ ShapeInfo(const Quadrature<dim_q> & quad,
+ const FiniteElement<dim, spacedim> &fe,
+ const unsigned int base_element = 0);
/**
* Initializes the data fields. Takes a one-dimensional quadrature
* dimensional element by a tensor product and that the zeroth shape
* function in zero evaluates to one.
*/
- template <int dim, int dim_q>
+ template <int dim, int spacedim, int dim_q>
void
- reinit(const Quadrature<dim_q> & quad,
- const FiniteElement<dim> &fe_dim,
- const unsigned int base_element = 0);
+ reinit(const Quadrature<dim_q> & quad,
+ const FiniteElement<dim, spacedim> &fe_dim,
+ const unsigned int base_element = 0);
/**
* Return which kinds of elements are supported by MatrixFree.
// ------------------------------------------ inline functions
- template <typename Number>
- template <int dim, int dim_q>
- inline ShapeInfo<Number>::ShapeInfo(const Quadrature<dim_q> & quad,
- const FiniteElement<dim> &fe_in,
- const unsigned int base_element_number)
- : element_type(tensor_general)
- , n_dimensions(0)
- , n_components(0)
- , n_q_points(0)
- , dofs_per_component_on_cell(0)
- , n_q_points_face(0)
- , dofs_per_component_on_face(0)
- {
- reinit(quad, fe_in, base_element_number);
- }
-
template <typename Number>
inline const UnivariateShapeData<Number> &
ShapeInfo<Number>::get_shape_data(const unsigned int dimension,
#include <deal.II/base/tensor_product_polynomials.h>
#include <deal.II/base/utilities.h>
+#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_dgp.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_poly.h>
- // ----------------- actual ShapeInfo functions --------------------
-
template <typename Number>
Number
get_first_array_element(const Number a)
+ // ----------------- actual ShapeInfo implementation --------------------
+
template <typename Number>
ShapeInfo<Number>::ShapeInfo()
: element_type(tensor_general)
template <typename Number>
- template <int dim, int spacedim>
- bool
- ShapeInfo<Number>::is_supported(const FiniteElement<dim, spacedim> &fe)
+ template <int dim, int spacedim, int dim_q>
+ inline ShapeInfo<Number>::ShapeInfo(
+ const Quadrature<dim_q> & quad,
+ const FiniteElement<dim, spacedim> &fe_in,
+ const unsigned int base_element_number)
+ : element_type(tensor_general)
+ , n_dimensions(0)
+ , n_components(0)
+ , n_q_points(0)
+ , dofs_per_component_on_cell(0)
+ , n_q_points_face(0)
+ , dofs_per_component_on_face(0)
{
- if (dim != spacedim)
- return false;
-
- for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
- {
- const FiniteElement<dim, spacedim> *fe_ptr = &(fe.base_element(base));
- if (fe_ptr->n_components() != 1)
- return false;
-
- // then check if the base element is supported or not
- if (dynamic_cast<const FE_Poly<dim, spacedim> *>(fe_ptr) != nullptr)
- {
- const FE_Poly<dim, spacedim> *fe_poly_ptr =
- dynamic_cast<const FE_Poly<dim, spacedim> *>(fe_ptr);
- // Simplices are a special case since the polynomial family is not
- // indicative of their support
- if (dynamic_cast<const FE_SimplexP<dim> *>(fe_poly_ptr) ||
- dynamic_cast<const FE_SimplexDGP<dim> *>(fe_poly_ptr) ||
- dynamic_cast<const FE_WedgeP<dim> *>(fe_poly_ptr) ||
- dynamic_cast<const FE_PyramidP<dim> *>(fe_poly_ptr))
- return true;
-
- if (dynamic_cast<const TensorProductPolynomials<dim> *>(
- &fe_poly_ptr->get_poly_space()) == nullptr &&
- dynamic_cast<const TensorProductPolynomials<
- dim,
- Polynomials::PiecewisePolynomial<double>> *>(
- &fe_poly_ptr->get_poly_space()) == nullptr &&
- dynamic_cast<const FE_DGP<dim, spacedim> *>(fe_ptr) ==
- nullptr &&
- dynamic_cast<const FE_Q_DG0<dim, spacedim> *>(fe_ptr) ==
- nullptr)
- return false;
- }
- else
- return false;
- }
-
- // if we arrived here, all base elements were supported so we can
- // support the present element
- return true;
+ reinit(quad, fe_in, base_element_number);
}
template <typename Number>
- template <int dim, int dim_q>
+ template <int dim, int spacedim, int dim_q>
void
- ShapeInfo<Number>::reinit(const Quadrature<dim_q> & quad_in,
- const FiniteElement<dim> &fe_in,
- const unsigned int base_element_number)
+ ShapeInfo<Number>::reinit(const Quadrature<dim_q> & quad_in,
+ const FiniteElement<dim, spacedim> &fe_in,
+ const unsigned int base_element_number)
{
+ static_assert(dim == spacedim,
+ "Currently, only the case dim=spacedim is implemented");
if (quad_in.is_tensor_product() == false ||
dynamic_cast<const FE_SimplexP<dim> *>(
&fe_in.base_element(base_element_number)) ||
+ template <typename Number>
+ template <int dim, int spacedim>
+ bool
+ ShapeInfo<Number>::is_supported(const FiniteElement<dim, spacedim> &fe)
+ {
+ if (dim != spacedim)
+ return false;
+
+ for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
+ {
+ const FiniteElement<dim, spacedim> *fe_ptr = &(fe.base_element(base));
+ if (fe_ptr->n_components() != 1)
+ return false;
+
+ // then check if the base element is supported or not
+ if (dynamic_cast<const FE_Poly<dim, spacedim> *>(fe_ptr) != nullptr)
+ {
+ const FE_Poly<dim, spacedim> *fe_poly_ptr =
+ dynamic_cast<const FE_Poly<dim, spacedim> *>(fe_ptr);
+ // Simplices are a special case since the polynomial family is not
+ // indicative of their support
+ if (dynamic_cast<const FE_SimplexP<dim> *>(fe_poly_ptr) ||
+ dynamic_cast<const FE_SimplexDGP<dim> *>(fe_poly_ptr) ||
+ dynamic_cast<const FE_WedgeP<dim> *>(fe_poly_ptr) ||
+ dynamic_cast<const FE_PyramidP<dim> *>(fe_poly_ptr))
+ return true;
+
+ if (dynamic_cast<const TensorProductPolynomials<dim> *>(
+ &fe_poly_ptr->get_poly_space()) == nullptr &&
+ dynamic_cast<const TensorProductPolynomials<
+ dim,
+ Polynomials::PiecewisePolynomial<double>> *>(
+ &fe_poly_ptr->get_poly_space()) == nullptr &&
+ dynamic_cast<const FE_DGP<dim, spacedim> *>(fe_ptr) ==
+ nullptr &&
+ dynamic_cast<const FE_Q_DG0<dim, spacedim> *>(fe_ptr) ==
+ nullptr)
+ return false;
+ }
+ else
+ return false;
+ }
+
+ // if we arrived here, all base elements were supported so we can
+ // support the present element
+ return true;
+ }
+
+
+
template <typename Number>
std::size_t
ShapeInfo<Number>::memory_consumption() const
for (deal_II_dimension : DIMENSIONS; deal_II_scalar : REAL_SCALARS)
{
- template void internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar>::
- reinit<deal_II_dimension>(
- const Quadrature<deal_II_dimension> &,
- const FiniteElement<deal_II_dimension, deal_II_dimension> &,
- const unsigned int);
+ template internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar>::
+ ShapeInfo(const Quadrature<deal_II_dimension> &,
+ const FiniteElement<deal_II_dimension, deal_II_dimension> &,
+ const unsigned int);
+
+ template void
+ internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar>::reinit(
+ const Quadrature<deal_II_dimension> &,
+ const FiniteElement<deal_II_dimension, deal_II_dimension> &,
+ const unsigned int);
#if deal_II_dimension > 1
- template void internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar>::
- reinit<deal_II_dimension>(
- const Quadrature<1> &,
- const FiniteElement<deal_II_dimension, deal_II_dimension> &,
- const unsigned int);
+ template internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar>::
+ ShapeInfo(const Quadrature<1> &,
+ const FiniteElement<deal_II_dimension, deal_II_dimension> &,
+ const unsigned int);
+
+ template void
+ internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar>::reinit(
+ const Quadrature<1> &,
+ const FiniteElement<deal_II_dimension, deal_II_dimension> &,
+ const unsigned int);
#endif
}
for (deal_II_dimension : DIMENSIONS;
deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED)
{
- template void internal::MatrixFreeFunctions::
- ShapeInfo<deal_II_scalar_vectorized>::reinit<deal_II_dimension>(
+ template internal::MatrixFreeFunctions::
+ ShapeInfo<deal_II_scalar_vectorized>::ShapeInfo(
const Quadrature<deal_II_dimension> &,
const FiniteElement<deal_II_dimension, deal_II_dimension> &,
const unsigned int);
+ template void
+ internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar_vectorized>::reinit(
+ const Quadrature<deal_II_dimension> &,
+ const FiniteElement<deal_II_dimension, deal_II_dimension> &,
+ const unsigned int);
+
#if deal_II_dimension > 1
- template void internal::MatrixFreeFunctions::
- ShapeInfo<deal_II_scalar_vectorized>::reinit<deal_II_dimension>(
+ template internal::MatrixFreeFunctions::
+ ShapeInfo<deal_II_scalar_vectorized>::ShapeInfo(
const Quadrature<1> &,
const FiniteElement<deal_II_dimension, deal_II_dimension> &,
const unsigned int);
+
+ template void
+ internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar_vectorized>::reinit(
+ const Quadrature<1> &,
+ const FiniteElement<deal_II_dimension, deal_II_dimension> &,
+ const unsigned int);
#endif
}
projects / dealii.git / commitdiff