*/
QCollection() = default;
+ /**
+ * Copy constructor.
+ */
+ template <int dim_in>
+ QCollection(const QCollection<dim_in> &other);
+
/**
* Conversion constructor. This constructor creates a QCollection from a
* single quadrature rule. More quadrature formulas can be added with
* push_back(), if desired, though it would probably be clearer to add all
* mappings the same way.
*/
- explicit QCollection(const Quadrature<dim> &quadrature);
+ template <int dim_in>
+ explicit QCollection(const Quadrature<dim_in> &quadrature);
/**
* Constructor. This constructor creates a QCollection from one or
* is later destroyed by this object upon destruction of the entire
* collection.
*/
+ template <int dim_in>
void
- push_back(const Quadrature<dim> &new_quadrature);
+ push_back(const Quadrature<dim_in> &new_quadrature);
/**
* Return a reference to the quadrature rule specified by the argument.
/* --------------- inline functions ------------------- */
+ template <int dim>
+ template <int dim_in>
+ QCollection<dim>::QCollection(const QCollection<dim_in> &other)
+ {
+ for (unsigned int i = 0; i < other.size(); ++i)
+ push_back(other[i]);
+ }
+
+
+
template <int dim>
template <class... QTypes>
QCollection<dim>::QCollection(const QTypes &... quadrature_objects)
{
- static_assert(is_base_of_all<Quadrature<dim>, QTypes...>::value,
- "Not all of the input arguments of this function "
- "are derived from Quadrature<dim>!");
-
// loop over all of the given arguments and add the quadrature objects to
// this collection. Inlining the definition of q_pointers causes internal
// compiler errors on GCC 7.1.1 so we define it separately:
- const auto q_pointers = {
- (static_cast<const Quadrature<dim> *>(&quadrature_objects))...};
- for (const auto p : q_pointers)
- push_back(*p);
+ if (is_base_of_all<Quadrature<dim>, QTypes...>::value)
+ {
+ const auto q_pointers = {
+ (reinterpret_cast<const Quadrature<dim> *>(&quadrature_objects))...};
+ for (const auto p : q_pointers)
+ push_back(*p);
+ }
+ else if (is_base_of_all<Quadrature<1>, QTypes...>::value)
+ {
+ const auto q_pointers = {
+ (reinterpret_cast<const Quadrature<1> *>(&quadrature_objects))...};
+ for (const auto p : q_pointers)
+ push_back(*p);
+ }
+ else
+ {
+ Assert(false, ExcNotImplemented());
+ }
}
template <int dim>
- inline QCollection<dim>::QCollection(const Quadrature<dim> &quadrature)
+ template <int dim_in>
+ inline QCollection<dim>::QCollection(const Quadrature<dim_in> &quadrature)
{
quadratures.push_back(std::make_shared<const Quadrature<dim>>(quadrature));
}
template <int dim>
+ template <int dim_in>
inline void
- QCollection<dim>::push_back(const Quadrature<dim> &new_quadrature)
+ QCollection<dim>::push_back(const Quadrature<dim_in> &new_quadrature)
{
quadratures.push_back(
std::make_shared<const Quadrature<dim>>(new_quadrature));
,
// select the correct base element from the given FE component
data(new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
- quadrature,
+ Quadrature<dim - is_face>(quadrature),
fe,
fe.component_to_base_index(first_selected_component).first))
, jacobian(nullptr)
* CellIterator::level() and CellIterator::index(), in order to allow
* for different kinds of iterators, e.g. standard DoFHandler,
* multigrid, etc.) on a fixed Triangulation. In addition, a mapping
- * and several quadrature formulas are given.
+ * and several 1D quadrature formulas are given.
*/
void
initialize(
const dealii::Triangulation<dim> & tria,
const std::vector<std::pair<unsigned int, unsigned int>> &cells,
const FaceInfo<VectorizedArrayType::size()> & faces,
- const std::vector<unsigned int> & active_fe_index,
- const Mapping<dim> & mapping,
- const std::vector<dealii::hp::QCollection<1>> &quad,
- const UpdateFlags update_flags_cells,
+ const std::vector<unsigned int> & active_fe_index,
+ const Mapping<dim> & mapping,
+ const std::vector<dealii::hp::QCollection<dim>> &quad,
+ const UpdateFlags update_flags_cells,
const UpdateFlags update_flags_boundary_faces,
const UpdateFlags update_flags_inner_faces,
const UpdateFlags update_flags_faces_by_cells);
* internal functions to initialize all data as requested by the user.
*/
static UpdateFlags
- compute_update_flags(const UpdateFlags update_flags,
- const std::vector<dealii::hp::QCollection<1>> &quad =
- std::vector<dealii::hp::QCollection<1>>());
+ compute_update_flags(
+ const UpdateFlags update_flags,
+ const std::vector<dealii::hp::QCollection<dim>> &quad =
+ std::vector<dealii::hp::QCollection<dim>>());
};
template <int dim, typename Number, typename VectorizedArrayType>
UpdateFlags
MappingInfo<dim, Number, VectorizedArrayType>::compute_update_flags(
- const UpdateFlags update_flags,
- const std::vector<dealii::hp::QCollection<1>> &quad)
+ const UpdateFlags update_flags,
+ const std::vector<dealii::hp::QCollection<dim>> &quad)
{
// this class is build around the evaluation of jacobians, so compute
// them in any case. The Jacobians will be inverted manually. Since we
const FaceInfo<VectorizedArrayType::size()> & face_info,
const std::vector<unsigned int> & active_fe_index,
const Mapping<dim> & mapping,
- const std::vector<dealii::hp::QCollection<1>> & quad,
+ const std::vector<dealii::hp::QCollection<dim>> & quad,
const UpdateFlags update_flags_cells,
const UpdateFlags update_flags_boundary_faces,
const UpdateFlags update_flags_inner_faces,
AssertIndexRange(0, n_hp_quads);
cell_data[my_q].descriptor.resize(n_hp_quads);
for (unsigned int q = 0; q < n_hp_quads; ++q)
- cell_data[my_q].descriptor[q].initialize(quad[my_q][q],
- update_default);
+ {
+ Assert(quad[my_q][q].is_tensor_product(), ExcNotImplemented());
+ for (unsigned int i = 1; i < dim; ++i)
+ {
+ Assert(quad[my_q][q].get_tensor_basis()[0] ==
+ quad[my_q][q].get_tensor_basis()[i],
+ ExcNotImplemented());
+ }
+
+ cell_data[my_q].descriptor[q].initialize(
+ quad[my_q][q].get_tensor_basis()[0], update_default);
+ }
face_data[my_q].descriptor.resize(n_hp_quads);
for (unsigned int hpq = 0; hpq < n_hp_quads; ++hpq)
- face_data[my_q].descriptor[hpq].initialize(
- quad[my_q][hpq], update_flags_boundary_faces);
+ {
+ Assert(quad[my_q][hpq].is_tensor_product(), ExcNotImplemented());
+ for (unsigned int i = 1; i < dim; ++i)
+ {
+ Assert(quad[my_q][hpq].get_tensor_basis()[0] ==
+ quad[my_q][hpq].get_tensor_basis()[i],
+ ExcNotImplemented());
+ }
+
+ 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);
for (unsigned int hpq = 0; hpq < n_hp_quads; ++hpq)
- face_data_by_cells[my_q].descriptor[hpq].initialize(quad[my_q][hpq],
- update_default);
+ {
+ Assert(quad[my_q][hpq].is_tensor_product(), ExcNotImplemented());
+ for (unsigned int i = 1; i < dim; ++i)
+ {
+ Assert(quad[my_q][hpq].get_tensor_basis()[0] ==
+ quad[my_q][hpq].get_tensor_basis()[i],
+ ExcNotImplemented());
+ }
+
+ face_data_by_cells[my_q].descriptor[hpq].initialize(
+ quad[my_q][hpq].get_tensor_basis()[0], update_default);
+ }
}
// In case we have no hp adaptivity (active_fe_index is empty), we have
{
FE_DGQ<dim> fe_geometry(mapping_degree);
for (unsigned int my_q = 0; my_q < cell_data.size(); ++my_q)
- shape_infos[my_q].reinit(cell_data[my_q].descriptor[0].quadrature_1d,
+ shape_infos[my_q].reinit(cell_data[my_q].descriptor[0].quadrature,
fe_geometry);
}
* This is the actual reinit function that sets up the indices for the
* DoFHandler case.
*/
- template <typename number2>
+ template <typename number2, int q_dim>
void
internal_reinit(
const Mapping<dim> & mapping,
const std::vector<const DoFHandler<dim, dim> *> & dof_handlers,
const std::vector<const AffineConstraints<number2> *> &constraint,
const std::vector<IndexSet> & locally_owned_set,
- const std::vector<hp::QCollection<1>> & quad,
+ const std::vector<hp::QCollection<q_dim>> & quad,
const AdditionalData & additional_data);
/**
internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
dof_handlers, additional_data.mg_level);
- std::vector<hp::QCollection<1>> quad_hp;
+ std::vector<hp::QCollection<dim>> quad_hp;
quad_hp.emplace_back(quad);
internal_reinit(StaticMappingQ1<dim>::mapping,
internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
dof_handlers, additional_data.mg_level);
- std::vector<hp::QCollection<1>> quad_hp;
+ std::vector<hp::QCollection<dim>> quad_hp;
quad_hp.emplace_back(quad);
internal_reinit(mapping,
std::vector<IndexSet> locally_owned_set =
internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
dof_handler, additional_data.mg_level);
- std::vector<hp::QCollection<1>> quad_hp;
+ std::vector<hp::QCollection<dim>> quad_hp;
for (unsigned int q = 0; q < quad.size(); ++q)
quad_hp.emplace_back(quad[q]);
internal_reinit(StaticMappingQ1<dim>::mapping,
std::vector<IndexSet> locally_owned_set =
internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
dof_handler, additional_data.mg_level);
- std::vector<hp::QCollection<1>> quad_hp;
+ std::vector<hp::QCollection<dim>> quad_hp;
quad_hp.emplace_back(quad);
internal_reinit(StaticMappingQ1<dim>::mapping,
dof_handler,
std::vector<IndexSet> locally_owned_set =
internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
dof_handler, additional_data.mg_level);
- std::vector<hp::QCollection<1>> quad_hp;
+ std::vector<hp::QCollection<dim>> quad_hp;
quad_hp.emplace_back(quad);
internal_reinit(mapping,
dof_handler,
std::vector<IndexSet> locally_owned_set =
internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
dof_handler, additional_data.mg_level);
- std::vector<hp::QCollection<1>> quad_hp;
+ std::vector<hp::QCollection<dim>> quad_hp;
for (unsigned int q = 0; q < quad.size(); ++q)
quad_hp.emplace_back(quad[q]);
internal_reinit(mapping,
template <int dim, typename Number, typename VectorizedArrayType>
-template <typename number2>
+template <typename number2, int q_dim>
void
MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
const Mapping<dim> & mapping,
const std::vector<const DoFHandler<dim, dim> *> & dof_handler,
const std::vector<const AffineConstraints<number2> *> &constraint,
const std::vector<IndexSet> & locally_owned_dofs,
- const std::vector<hp::QCollection<1>> & quad,
+ const std::vector<hp::QCollection<q_dim>> & quad,
const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
&additional_data)
{
Table<2, internal::MatrixFreeFunctions::ShapeInfo<double>> shape_info_dummy(
shape_info.size(0), shape_info.size(2));
{
- QGauss<1> quad(1);
+ QGauss<dim> quad(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();
/**
* Constructor that initializes the data fields using the reinit method.
*/
- template <int dim>
- ShapeInfo(const Quadrature<1> & quad,
+ template <int dim, int dim_q>
+ ShapeInfo(const Quadrature<dim_q> & quad,
const FiniteElement<dim> &fe,
const unsigned int base_element = 0);
* dimensional element by a tensor product and that the zeroth shape
* function in zero evaluates to one.
*/
- template <int dim>
+ template <int dim, int dim_q>
void
- reinit(const Quadrature<1> & quad,
+ reinit(const Quadrature<dim_q> & quad,
const FiniteElement<dim> &fe_dim,
const unsigned int base_element = 0);
// ------------------------------------------ inline functions
template <typename Number>
- template <int dim>
- inline ShapeInfo<Number>::ShapeInfo(const Quadrature<1> & quad,
+ 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)
template <typename Number>
- template <int dim>
+ template <int dim, int dim_q>
void
- ShapeInfo<Number>::reinit(const Quadrature<1> & quad,
+ ShapeInfo<Number>::reinit(const Quadrature<dim_q> & quad_in,
const FiniteElement<dim> &fe_in,
const unsigned int base_element_number)
{
+ Assert(quad_in.is_tensor_product(), ExcNotImplemented());
+ const auto quad = quad_in.get_tensor_basis()[0];
+
const FiniteElement<dim> *fe = &fe_in.base_element(base_element_number);
n_dimensions = dim;
n_components = fe_in.n_components();
const std::vector<
const AffineConstraints<deal_II_scalar_vectorized::value_type> *> &,
const std::vector<IndexSet> &,
- const std::vector<hp::QCollection<1>> &,
+ const std::vector<hp::QCollection<deal_II_dimension>> &,
const AdditionalData &);
template const DoFHandler<deal_II_dimension> &
const std::vector<const DoFHandler<deal_II_dimension> *> &,
const std::vector<const AffineConstraints<double> *> &,
const std::vector<IndexSet> &,
- const std::vector<hp::QCollection<1>> &,
+ const std::vector<hp::QCollection<deal_II_dimension>> &,
const AdditionalData &);
}
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);
+
+#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);
+#endif
}
for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
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>(
+ 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>(
const Quadrature<1> &,
const FiniteElement<deal_II_dimension, deal_II_dimension> &,
const unsigned int);
+#endif
}
for (deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED)
projects / dealii.git / commitdiff