* Stored the level of the mesh to be worked on.
*/
unsigned int mg_level;
+
+ /**
+ * Stores the index of the first DoFHandler that is in hp-mode. If no
+ * DoFHandler is in hp-mode, the value is 0.
+ */
+ unsigned int first_hp_dof_handler_index;
};
MatrixFree<dim, Number, VectorizedArrayType>::get_cell_active_fe_index(
const std::pair<unsigned int, unsigned int> range) const
{
- const auto &fe_indices = dof_info[0].cell_active_fe_index;
+ const auto &fe_indices =
+ dof_info[first_hp_dof_handler_index].cell_active_fe_index;
if (fe_indices.empty() == true ||
dof_handlers[0]->get_fe_collection().size() == 1)
const std::pair<unsigned int, unsigned int> range,
const bool is_interior_face) const
{
- const auto &fe_indices = dof_info[0].cell_active_fe_index;
+ const auto &fe_indices =
+ dof_info[first_hp_dof_handler_index].cell_active_fe_index;
if (fe_indices.empty() == true)
return 0;
const unsigned int cell_batch_index) const
{
AssertIndexRange(0, dof_info.size());
- AssertIndexRange(cell_batch_index, dof_info[0].cell_active_fe_index.size());
- if (dof_info[0].cell_active_fe_index.empty())
+ AssertIndexRange(
+ cell_batch_index,
+ dof_info[first_hp_dof_handler_index].cell_active_fe_index.size());
+ if (dof_info[first_hp_dof_handler_index].cell_active_fe_index.empty())
return 0;
else
- return dof_info[0].cell_active_fe_index[cell_batch_index];
+ return dof_info[first_hp_dof_handler_index]
+ .cell_active_fe_index[cell_batch_index];
}
const unsigned int face_batch_index) const
{
AssertIndexRange(face_batch_index, face_info.faces.size());
- if (dof_info[0].cell_active_fe_index.empty())
+ if (dof_info[first_hp_dof_handler_index].cell_active_fe_index.empty())
return std::make_pair(0U, 0U);
std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
face_info.faces[face_batch_index].cells_interior[v] !=
numbers::invalid_unsigned_int;
++v)
- result.first = std::max(
- result.first,
- dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
- .cells_interior[v] /
- VectorizedArrayType::size()]);
+ result.first =
+ std::max(result.first,
+ dof_info[first_hp_dof_handler_index].cell_active_fe_index
+ [face_info.faces[face_batch_index].cells_interior[v] /
+ VectorizedArrayType::size()]);
if (face_info.faces[face_batch_index].cells_exterior[0] !=
numbers::invalid_unsigned_int)
for (unsigned int v = 0;
face_info.faces[face_batch_index].cells_exterior[v] !=
numbers::invalid_unsigned_int;
++v)
- result.second = std::max(
- result.second,
- dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
- .cells_exterior[v] /
- VectorizedArrayType::size()]);
+ result.second =
+ std::max(result.second,
+ dof_info[first_hp_dof_handler_index].cell_active_fe_index
+ [face_info.faces[face_batch_index].cells_exterior[v] /
+ VectorizedArrayType::size()]);
else
result.second = numbers::invalid_unsigned_int;
return result;
, indices_are_initialized(false)
, mapping_is_initialized(false)
, mg_level(numbers::invalid_unsigned_int)
+ , first_hp_dof_handler_index(0)
{}
indices_are_initialized = v.indices_are_initialized;
mapping_is_initialized = v.mapping_is_initialized;
mg_level = v.mg_level;
+ first_hp_dof_handler_index = v.first_hp_dof_handler_index;
}
.reinit(quad[nq][q_no], dof_handler[no]->get_fe(fe_no), b);
}
+ // Vector of DoFHandler indices of those that are in hp-mode
+ const std::vector<unsigned int> hp_dof_handler_index =
+ std::invoke([&dof_handler]() {
+ std::vector<unsigned int> hp_dof_handler_index;
+ hp_dof_handler_index.reserve(dof_handler.size());
+ for (unsigned int i = 0; i < dof_handler.size(); ++i)
+ {
+ if (dof_handler[i]->has_hp_capabilities())
+ hp_dof_handler_index.emplace_back(i);
+ }
+ return hp_dof_handler_index;
+ });
+ if (hp_dof_handler_index.size() > 0)
+ first_hp_dof_handler_index = hp_dof_handler_index[0];
+
// Store pointers to AffineConstraints objects if Number type matches
affine_constraints.resize(constraints.size());
for (unsigned int no = 0; no < constraints.size(); ++no)
}
}
-
// subdivide cell, face and boundary face partitioner data, s.t., all
// ranges have the same active FE indices
if (task_info.scheme != internal::MatrixFreeFunctions::TaskInfo::
for (unsigned int i = 0; i < this->n_active_fe_indices(); ++i)
{
const auto cell_subrange =
- this->create_cell_subrange_hp_by_index(range, i);
+ this->create_cell_subrange_hp_by_index(
+ range, i, first_hp_dof_handler_index);
if (cell_subrange.second <= cell_subrange.first)
continue;
const unsigned int face_type,
const unsigned int fe_index_interior,
const unsigned int fe_index_exterior,
- const unsigned int dof_handler_index =
- 0) -> std::pair<unsigned int, unsigned int> {
+ const unsigned int dof_handler_index)
+ -> std::pair<unsigned int, unsigned int> {
const unsigned int n_face_types =
std::max<unsigned int>(dim - 1, 1);
++j)
{
const auto subrange =
- create_inner_face_subrange_hp_by_index(range,
- t,
- i,
- j);
+ create_inner_face_subrange_hp_by_index(
+ range, t, i, j, first_hp_dof_handler_index);
if (subrange.second <= subrange.first)
continue;
[&](const std::pair<unsigned int, unsigned int> &range,
const unsigned int face_type,
const unsigned int fe_index,
- const unsigned int dof_handler_index =
- 0) -> std::pair<unsigned int, unsigned int> {
+ const unsigned int dof_handler_index)
+ -> std::pair<unsigned int, unsigned int> {
const unsigned int n_face_types =
std::max<unsigned int>(dim - 1, 1);
++i)
{
const auto cell_subrange =
- create_boundary_face_subrange_hp_by_index(range,
- t,
- i);
+ create_boundary_face_subrange_hp_by_index(
+ range, t, i, first_hp_dof_handler_index);
if (cell_subrange.second <= cell_subrange.first)
continue;
// determined in @p extract_local_to_global_indices.
if (additional_data.initialize_mapping == true)
{
- if (dof_handler.size() > 1)
+ // check if active FE indices of all hp-DoFHandlers are the same.
+ for (unsigned int i = 1; i < hp_dof_handler_index.size(); ++i)
{
- // check if all DoHandlers are in the same hp-mode; and if hp-
- // capabilities are enabled: check if active FE indices of all
- // DoFHandlers are the same.
- for (unsigned int i = 1; i < dof_handler.size(); ++i)
- {
- Assert(dof_handler[0]->has_hp_capabilities() ==
- dof_handler[i]->has_hp_capabilities(),
- ExcNotImplemented());
-
- if (dof_handler[0]->has_hp_capabilities())
- {
- Assert(dof_info[0].cell_active_fe_index ==
- dof_info[i].cell_active_fe_index,
- ExcNotImplemented());
- }
- }
+ Assert(dof_info[hp_dof_handler_index[0]].cell_active_fe_index ==
+ dof_info[hp_dof_handler_index[i]].cell_active_fe_index,
+ ExcNotImplemented());
}
// Will the piola transform be used? If so we need to update
dof_handler[0]->get_triangulation(),
cell_level_index,
face_info,
- dof_handler[0]->has_hp_capabilities() ?
- dof_info[0].cell_active_fe_index :
+ hp_dof_handler_index.size() > 0 ?
+ dof_info[hp_dof_handler_index[0]].cell_active_fe_index :
std::vector<unsigned int>(),
mapping,
quad,
const std::shared_ptr<hp::MappingCollection<dim>> &mapping)
{
AssertDimension(shape_info.size(1), mapping_info.cell_data.size());
- mapping_info.update_mapping(dof_handlers[0]->get_triangulation(),
- cell_level_index,
- face_info,
- dof_info[0].cell_active_fe_index,
- mapping);
+ mapping_info.update_mapping(
+ dof_handlers[0]->get_triangulation(),
+ cell_level_index,
+ face_info,
+ dof_info[first_hp_dof_handler_index].cell_active_fe_index,
+ mapping);
}
}
}
- const bool hp_functionality_enabled =
- std::any_of(dof_handlers.begin(), dof_handlers.end(), [](const auto &dh) {
- return dh->has_hp_capabilities();
+ const auto [hp_functionality_enabled, first_hp_dof_handler_index] =
+ std::invoke([&dof_handlers]() -> std::pair<bool, unsigned int> {
+ bool hp_functionality_enabled = false;
+ unsigned int first_hp_dof_handler_index = 0;
+ for (unsigned int i = 0; i < dof_handlers.size(); ++i)
+ {
+ if (dof_handlers[i]->has_hp_capabilities())
+ {
+ hp_functionality_enabled = true;
+ first_hp_dof_handler_index = i;
+ break;
+ }
+ }
+ return {hp_functionality_enabled, first_hp_dof_handler_index};
});
const unsigned int n_lanes = task_info.vectorization_length;
parent_relation[i] = position;
++position;
}
- task_info.create_blocks_serial(subdomain_boundary_cells,
- max_dofs_per_cell,
- hp_functionality_enabled,
- dof_info[0].cell_active_fe_index,
- strict_categories,
- parent_relation,
- renumbering,
- irregular_cells);
+ task_info.create_blocks_serial(
+ subdomain_boundary_cells,
+ max_dofs_per_cell,
+ hp_functionality_enabled,
+ dof_info[first_hp_dof_handler_index].cell_active_fe_index,
+ strict_categories,
+ parent_relation,
+ renumbering,
+ irregular_cells);
}
else
{
task_info.initial_setup_blocks_tasks(subdomain_boundary_cells,
renumbering,
irregular_cells);
- task_info.guess_block_size(dof_info[0].dofs_per_cell[0]);
+ const internal::MatrixFreeFunctions::DoFInfo &dof_info_hp =
+ dof_info[first_hp_dof_handler_index];
+ task_info.guess_block_size(dof_info_hp.dofs_per_cell[0]);
unsigned int n_cell_batches_before =
*(task_info.cell_partition_data.end() - 2);
{
irregular_cells.resize(0);
irregular_cells.resize(task_info.cell_partition_data.back() +
- 2 * dof_info[0].max_fe_index);
+ 2 * dof_info_hp.max_fe_index);
std::vector<std::vector<unsigned int>> renumbering_fe_index;
- renumbering_fe_index.resize(dof_info[0].max_fe_index);
+ renumbering_fe_index.resize(dof_info_hp.max_fe_index);
unsigned int counter;
n_cell_batches_before = 0;
for (counter = 0;
{
AssertIndexRange(counter, renumbering.size());
AssertIndexRange(renumbering[counter],
- dof_info[0].cell_active_fe_index.size());
+ dof_info_hp.cell_active_fe_index.size());
renumbering_fe_index
- [dof_info[0].cell_active_fe_index[renumbering[counter]]]
+ [dof_info_hp.cell_active_fe_index[renumbering[counter]]]
.push_back(renumbering[counter]);
}
counter = 0;
- for (unsigned int j = 0; j < dof_info[0].max_fe_index; ++j)
+ for (unsigned int j = 0; j < dof_info_hp.max_fe_index; ++j)
{
for (const auto jj : renumbering_fe_index[j])
renumbering[counter++] = jj;
counter++)
{
renumbering_fe_index
- [dof_info[0].cell_active_fe_index.empty() ?
+ [dof_info_hp.cell_active_fe_index.empty() ?
0 :
- dof_info[0].cell_active_fe_index[renumbering[counter]]]
+ dof_info_hp.cell_active_fe_index[renumbering[counter]]]
.push_back(renumbering[counter]);
}
counter = start_nonboundary * n_lanes;
- for (unsigned int j = 0; j < dof_info[0].max_fe_index; ++j)
+ for (unsigned int j = 0; j < dof_info_hp.max_fe_index; ++j)
{
for (const auto jj : renumbering_fe_index[j])
renumbering[counter++] = jj;
}
AssertIndexRange(n_cell_batches_before,
task_info.cell_partition_data.back() +
- 2 * dof_info[0].max_fe_index + 1);
+ 2 * dof_info_hp.max_fe_index + 1);
irregular_cells.resize(n_cell_batches_before + n_ghost_slots);
*(task_info.cell_partition_data.end() - 2) =
n_cell_batches_before;
#endif
}
if (task_info.n_active_cells > 0)
- dof_info[0].make_connectivity_graph(task_info,
+ dof_info_hp.make_connectivity_graph(task_info,
renumbering,
connectivity);
- task_info.make_thread_graph(dof_info[0].cell_active_fe_index,
+ task_info.make_thread_graph(dof_info_hp.cell_active_fe_index,
connectivity,
renumbering,
irregular_cells,
true);
}
+ const internal::MatrixFreeFunctions::DoFInfo &dof_info_hp =
+ dof_info[first_hp_dof_handler_index];
+
if (additional_data.mapping_update_flags_inner_faces != update_default)
internal::MatrixFreeFunctions::collect_faces_vectorization(
face_setup.inner_faces,
hard_vectorization_boundary,
task_info.face_partition_data,
face_info.faces,
- dof_info[0].cell_active_fe_index);
+ dof_info_hp.cell_active_fe_index);
// on boundary faces, we must also respect the vectorization boundary of
// the inner faces because we might have dependencies on ghosts of
hard_vectorization_boundary,
task_info.boundary_partition_data,
face_info.faces,
- dof_info[0].cell_active_fe_index);
+ dof_info_hp.cell_active_fe_index);
// for the other ghosted faces, there are no scheduling restrictions
hard_vectorization_boundary.clear();
hard_vectorization_boundary,
task_info.ghost_face_partition_data,
face_info.faces,
- dof_info[0].cell_active_fe_index);
+ dof_info_hp.cell_active_fe_index);
hard_vectorization_boundary.clear();
hard_vectorization_boundary.resize(
task_info.refinement_edge_face_partition_data.size(), false);
hard_vectorization_boundary,
task_info.refinement_edge_face_partition_data,
face_info.faces,
- dof_info[0].cell_active_fe_index);
+ dof_info_hp.cell_active_fe_index);
cell_level_index.resize(
cell_level_index.size() +
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2012 - 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+
+// this function tests the correctness of the implementation of matrix free
+// operations in getting the function values, the function gradients, and the
+// function Laplacians on a hypecube mesh with adaptive refinement with
+// components in two different DoFHandlers, one of which is in hp-mode.
+
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/q_collection.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+#include <utility>
+#include <vector>
+
+#include "../tests.h"
+
+
+template <int dim,
+ int fe_degree,
+ int n_q_points_1d = fe_degree + 1,
+ typename Number = double>
+class MatrixFreeTest
+{
+public:
+ using VectorType = std::vector<Vector<Number>>;
+
+ MatrixFreeTest(const MatrixFree<dim, Number> &data_in,
+ const hp::FECollection<dim> &fe_collection)
+ : data(data_in)
+ , fe_val0(data.get_dof_handler(0).get_fe(),
+ Quadrature<dim>(data.get_quadrature(0)),
+ update_values | update_gradients | update_hessians)
+ , hp_fe_val1(fe_collection,
+ hp::QCollection<dim>(data.get_quadrature(1)),
+ update_values | update_gradients | update_hessians){};
+
+ void
+ operator()(const MatrixFree<dim, Number> &data,
+ VectorType &,
+ const VectorType &src,
+ const std::pair<unsigned int, unsigned int> &cell_range) const
+ {
+ const auto active_fe_index = data.get_cell_range_category(cell_range);
+
+ FEEvaluation<dim, -1, 0, 1, Number> fe_eval0(data, 0, 0);
+ FEEvaluation<dim, -1, 0, 1, Number> fe_eval1(
+ data, 1, 1, 0, active_fe_index);
+ std::vector<double> reference_values0(fe_eval0.n_q_points);
+ std::vector<Tensor<1, dim>> reference_grads0(fe_eval0.n_q_points);
+ std::vector<Tensor<2, dim>> reference_hess0(fe_eval0.n_q_points);
+ std::vector<double> reference_values1(fe_eval1.n_q_points);
+ std::vector<Tensor<1, dim>> reference_grads1(fe_eval1.n_q_points);
+ std::vector<Tensor<2, dim>> reference_hess1(fe_eval1.n_q_points);
+ for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+ {
+ fe_eval0.reinit(cell);
+ fe_eval0.read_dof_values(src[0]);
+ fe_eval0.evaluate(EvaluationFlags::values | EvaluationFlags::gradients |
+ EvaluationFlags::hessians);
+
+ fe_eval1.reinit(cell);
+ fe_eval1.read_dof_values(src[1]);
+ fe_eval1.evaluate(EvaluationFlags::values | EvaluationFlags::gradients |
+ EvaluationFlags::hessians);
+
+ // compare values with the ones the FEValues
+ // gives us. Those are seen as reference
+ for (unsigned int j = 0; j < data.n_active_entries_per_cell_batch(cell);
+ ++j)
+ {
+ // FE 0
+ fe_val0.reinit(data.get_cell_iterator(cell, j, 0));
+ fe_val0.get_function_values(src[0], reference_values0);
+ fe_val0.get_function_gradients(src[0], reference_grads0);
+ fe_val0.get_function_hessians(src[0], reference_hess0);
+
+ for (int q = 0; q < (int)fe_eval0.n_q_points; ++q)
+ {
+ errors[0] +=
+ std::fabs(fe_eval0.get_value(q)[j] - reference_values0[q]);
+ for (unsigned int d = 0; d < dim; ++d)
+ errors[1] += std::fabs(fe_eval0.get_gradient(q)[d][j] -
+ reference_grads0[q][d]);
+ errors[2] += std::fabs(fe_eval0.get_laplacian(q)[j] -
+ trace(reference_hess0[q]));
+ total[0] += std::fabs(reference_values0[q]);
+ for (unsigned int d = 0; d < dim; ++d)
+ total[1] += std::fabs(reference_grads0[q][d]);
+ total[2] += std::fabs(fe_eval0.get_laplacian(q)[j]);
+ }
+
+ // FE 1
+ hp_fe_val1.reinit(data.get_cell_iterator(cell, j, 1));
+ const auto &fe_val1 = hp_fe_val1.get_present_fe_values();
+ fe_val1.get_function_values(src[1], reference_values1);
+ fe_val1.get_function_gradients(src[1], reference_grads1);
+ fe_val1.get_function_hessians(src[1], reference_hess1);
+
+ for (int q = 0; q < (int)fe_eval1.n_q_points; ++q)
+ {
+ errors[3] +=
+ std::fabs(fe_eval1.get_value(q)[j] - reference_values1[q]);
+ for (unsigned int d = 0; d < dim; ++d)
+ errors[4] += std::fabs(fe_eval1.get_gradient(q)[d][j] -
+ reference_grads1[q][d]);
+ errors[5] += std::fabs(fe_eval1.get_laplacian(q)[j] -
+ trace(reference_hess1[q]));
+ total[3] += std::fabs(reference_values1[q]);
+ for (unsigned int d = 0; d < dim; ++d)
+ total[4] += std::fabs(reference_grads1[q][d]);
+ total[5] += std::fabs(fe_eval1.get_laplacian(q)[j]);
+ }
+ }
+ }
+ }
+
+ void
+ test_functions(const VectorType &src) const
+ {
+ for (unsigned int i = 0; i < 3 * 2; ++i)
+ {
+ errors[i] = 0;
+ total[i] = 0;
+ }
+ VectorType dst_dummy;
+ data.cell_loop(
+ &MatrixFreeTest<dim, fe_degree, n_q_points_1d, Number>::operator(),
+ this,
+ dst_dummy,
+ src);
+
+ // for doubles, use a stricter condition then
+ // for floats for the relative error size
+ for (unsigned int i = 0; i < 2; ++i)
+ {
+ if (std::is_same_v<Number, double> == true)
+ {
+ deallog << "Error function values FE " << i << ": "
+ << errors[i * 3 + 0] / total[i * 3 + 0] << std::endl;
+ deallog << "Error function gradients FE " << i << ": "
+ << errors[i * 3 + 1] / total[i * 3 + 1] << std::endl;
+
+ // need to set quite a loose tolerance because
+ // FEValues approximates Hessians with finite
+ // differences, which are not so
+ // accurate. moreover, Hessians are quite
+ // large since we chose random numbers. for
+ // some elements, it might also be zero
+ // (linear elements on quadrilaterals), so
+ // need to check for division by 0, too.
+ const double output2 =
+ total[i * 3 + 2] == 0 ? 0. : errors[i * 3 + 2] / total[i * 3 + 2];
+ deallog << "Error function Laplacians FE " << i << ": " << output2
+ << std::endl;
+ }
+ else if (std::is_same_v<Number, float> == true)
+ {
+ deallog << "Error function values FE " << i << ": "
+ << errors[i * 3 + 0] / total[i * 3 + 0] << std::endl;
+ deallog << "Error function gradients FE " << i << ": "
+ << errors[i * 3 + 1] / total[i * 3 + 1] << std::endl;
+ const double output2 =
+ total[i * 3 + 2] == 0 ? 0. : errors[i * 3 + 2] / total[i * 3 + 2];
+ deallog << "Error function Laplacians FE " << i << ": " << output2
+ << std::endl;
+ }
+ }
+ }
+
+private:
+ const MatrixFree<dim, Number> &data;
+ mutable FEValues<dim> fe_val0;
+ mutable hp::FEValues<dim> hp_fe_val1;
+ mutable double errors[6], total[6];
+};
+
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+ using number = double;
+
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(1);
+ tria.begin(tria.n_levels() - 1)->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ tria.begin(tria.n_levels() - 1)->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ tria.begin(tria.n_levels() - 1)->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ tria.begin(tria.n_levels() - 1)->set_refine_flag();
+ tria.begin_active(tria.n_levels() - 2)->set_refine_flag();
+ tria.begin_active(tria.n_levels() - 3)->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+
+ FE_Q<dim> fe0(fe_degree);
+ FE_Q<dim> fe1(fe_degree + 1);
+
+ DoFHandler<dim> dof0(tria);
+ dof0.distribute_dofs(fe0);
+
+ hp::FECollection<dim> fe_collection;
+ fe_collection.push_back(fe0);
+ fe_collection.push_back(fe1);
+
+ DoFHandler<dim> dof1(tria);
+ for (const auto &cell : dof1.active_cell_iterators())
+ cell->set_active_fe_index(Testing::rand() % fe_collection.size());
+ dof1.distribute_dofs(fe_collection);
+
+ std::vector<const DoFHandler<dim> *> dof(2);
+ dof[0] = &dof0;
+ dof[1] = &dof1;
+
+ deallog << "Testing " << fe0.get_name() << " and " << fe1.get_name()
+ << std::endl;
+ // std::cout << "Number of cells: " << tria.n_active_cells() << std::endl;
+
+ std::vector<Vector<double>> src(dof.size());
+ for (unsigned int no = 0; no < dof.size(); ++no)
+ src[no].reinit(dof[no]->n_dofs());
+
+ std::vector<const AffineConstraints<double> *> constraints(2);
+ AffineConstraints<double> constraint0;
+ DoFTools::make_hanging_node_constraints(*dof[0], constraint0);
+ constraint0.close();
+ constraints[0] = &constraint0;
+ AffineConstraints<double> constraint1;
+ DoFTools::make_hanging_node_constraints(*dof[1], constraint1);
+ constraint1.close();
+ constraints[1] = &constraint1;
+
+ // create vector with random entries
+ for (unsigned int no = 0; no < 2; ++no)
+ for (unsigned int i = 0; i < dof[no]->n_dofs(); ++i)
+ {
+ if (constraints[no]->is_constrained(i))
+ continue;
+ const double entry = random_value<double>();
+ src[no](i) = entry;
+ }
+
+
+ constraints[0]->distribute(src[0]);
+ constraints[1]->distribute(src[1]);
+ MatrixFree<dim, number> mf_data;
+ {
+ std::vector<Quadrature<dim>> quad;
+ for (unsigned int no = 0; no < 2; ++no)
+ quad.push_back(QGauss<dim>(fe_degree + 1 + no));
+ mf_data.reinit(MappingQ1<dim>{},
+ dof,
+ constraints,
+ quad,
+ typename MatrixFree<dim, number>::AdditionalData(
+ MatrixFree<dim, number>::AdditionalData::none));
+ }
+
+ MatrixFreeTest<dim, fe_degree, fe_degree + 1, number> mf(mf_data,
+ fe_collection);
+ mf.test_functions(src);
+ deallog << std::endl;
+}
+
+
+int
+main()
+{
+ initlog();
+ deallog << std::setprecision(7);
+
+ {
+ deallog.push("2d");
+ test<2, 1>();
+ test<2, 2>();
+ test<2, 3>();
+ test<2, 4>();
+ deallog.pop();
+ deallog.push("3d");
+ test<3, 1>();
+ test<3, 2>();
+ test<3, 3>();
+ deallog.pop();
+ }
+}
--- /dev/null
+
+DEAL:2d::Testing FE_Q<2>(1) and FE_Q<2>(2)
+DEAL:2d::Error function values FE 0: 0
+DEAL:2d::Error function gradients FE 0: 0
+DEAL:2d::Error function Laplacians FE 0: 0
+DEAL:2d::Error function values FE 1: 0
+DEAL:2d::Error function gradients FE 1: 0
+DEAL:2d::Error function Laplacians FE 1: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(2) and FE_Q<2>(3)
+DEAL:2d::Error function values FE 0: 0
+DEAL:2d::Error function gradients FE 0: 0
+DEAL:2d::Error function Laplacians FE 0: 0
+DEAL:2d::Error function values FE 1: 0
+DEAL:2d::Error function gradients FE 1: 0
+DEAL:2d::Error function Laplacians FE 1: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(3) and FE_Q<2>(4)
+DEAL:2d::Error function values FE 0: 0
+DEAL:2d::Error function gradients FE 0: 0
+DEAL:2d::Error function Laplacians FE 0: 0
+DEAL:2d::Error function values FE 1: 0
+DEAL:2d::Error function gradients FE 1: 0
+DEAL:2d::Error function Laplacians FE 1: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(4) and FE_Q<2>(5)
+DEAL:2d::Error function values FE 0: 0
+DEAL:2d::Error function gradients FE 0: 0
+DEAL:2d::Error function Laplacians FE 0: 0
+DEAL:2d::Error function values FE 1: 0
+DEAL:2d::Error function gradients FE 1: 0
+DEAL:2d::Error function Laplacians FE 1: 0
+DEAL:2d::
+DEAL:3d::Testing FE_Q<3>(1) and FE_Q<3>(2)
+DEAL:3d::Error function values FE 0: 0
+DEAL:3d::Error function gradients FE 0: 0
+DEAL:3d::Error function Laplacians FE 0: 0
+DEAL:3d::Error function values FE 1: 0
+DEAL:3d::Error function gradients FE 1: 0
+DEAL:3d::Error function Laplacians FE 1: 0
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(2) and FE_Q<3>(3)
+DEAL:3d::Error function values FE 0: 0
+DEAL:3d::Error function gradients FE 0: 0
+DEAL:3d::Error function Laplacians FE 0: 0
+DEAL:3d::Error function values FE 1: 0
+DEAL:3d::Error function gradients FE 1: 0
+DEAL:3d::Error function Laplacians FE 1: 0
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(3) and FE_Q<3>(4)
+DEAL:3d::Error function values FE 0: 0
+DEAL:3d::Error function gradients FE 0: 0
+DEAL:3d::Error function Laplacians FE 0: 0
+DEAL:3d::Error function values FE 1: 0
+DEAL:3d::Error function gradients FE 1: 0
+DEAL:3d::Error function Laplacians FE 1: 0
+DEAL:3d::
projects / dealii.git / commitdiff