From ca6d467ed87ef465c7b281c3a1a0d453125c3aed Mon Sep 17 00:00:00 2001 From: Nils Much Date: Wed, 23 Oct 2024 16:02:17 +0200 Subject: [PATCH] MatrixFree: allow non-hp and hp-DoFHandlers --- include/deal.II/matrix_free/matrix_free.h | 43 ++- .../matrix_free/matrix_free.templates.h | 148 ++++---- tests/matrix_free/get_functions_multife3.cc | 324 ++++++++++++++++++ .../matrix_free/get_functions_multife3.output | 57 +++ 4 files changed, 492 insertions(+), 80 deletions(-) create mode 100644 tests/matrix_free/get_functions_multife3.cc create mode 100644 tests/matrix_free/get_functions_multife3.output diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index 63df0bb8e3..ff35f4a35e 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -2363,6 +2363,12 @@ private: * 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; }; @@ -2695,7 +2701,8 @@ unsigned int MatrixFree::get_cell_active_fe_index( const std::pair 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) @@ -2717,7 +2724,8 @@ MatrixFree::get_face_active_fe_index( const std::pair 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; @@ -2965,11 +2973,14 @@ MatrixFree::get_cell_category( 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]; } @@ -2980,7 +2991,7 @@ MatrixFree::get_face_category( 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 result = std::make_pair(0U, 0U); @@ -2989,11 +3000,11 @@ MatrixFree::get_face_category( 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; @@ -3001,11 +3012,11 @@ MatrixFree::get_face_category( 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; diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index 2174f3cbc4..11c72fd1b6 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -81,6 +81,7 @@ MatrixFree::MatrixFree() , indices_are_initialized(false) , mapping_is_initialized(false) , mg_level(numbers::invalid_unsigned_int) + , first_hp_dof_handler_index(0) {} @@ -356,6 +357,7 @@ MatrixFree::copy_from( 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; } @@ -429,6 +431,21 @@ MatrixFree::internal_reinit( .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 hp_dof_handler_index = + std::invoke([&dof_handler]() { + std::vector 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) @@ -547,7 +564,6 @@ MatrixFree::internal_reinit( } } - // 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:: @@ -576,7 +592,8 @@ MatrixFree::internal_reinit( 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; @@ -605,8 +622,8 @@ MatrixFree::internal_reinit( 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 { + const unsigned int dof_handler_index) + -> std::pair { const unsigned int n_face_types = std::max(dim - 1, 1); @@ -677,10 +694,8 @@ MatrixFree::internal_reinit( ++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; @@ -708,8 +723,8 @@ MatrixFree::internal_reinit( [&](const std::pair &range, const unsigned int face_type, const unsigned int fe_index, - const unsigned int dof_handler_index = - 0) -> std::pair { + const unsigned int dof_handler_index) + -> std::pair { const unsigned int n_face_types = std::max(dim - 1, 1); @@ -774,9 +789,8 @@ MatrixFree::internal_reinit( ++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; @@ -796,24 +810,12 @@ MatrixFree::internal_reinit( // 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 @@ -837,8 +839,8 @@ MatrixFree::internal_reinit( 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(), mapping, quad, @@ -895,11 +897,12 @@ MatrixFree::update_mapping( const std::shared_ptr> &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); } @@ -1491,9 +1494,20 @@ namespace internal } } - 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 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; @@ -1559,14 +1573,15 @@ namespace internal 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 { @@ -1579,7 +1594,9 @@ namespace internal 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); @@ -1622,9 +1639,9 @@ namespace internal { 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> 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; @@ -1634,13 +1651,13 @@ namespace internal { 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; @@ -1657,13 +1674,13 @@ namespace internal 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; @@ -1675,7 +1692,7 @@ namespace internal } 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; @@ -1743,11 +1760,11 @@ namespace internal #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, @@ -1983,13 +2000,16 @@ MatrixFree::initialize_indices( 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 @@ -1999,7 +2019,7 @@ MatrixFree::initialize_indices( 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(); @@ -2010,7 +2030,7 @@ MatrixFree::initialize_indices( 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); @@ -2019,7 +2039,7 @@ MatrixFree::initialize_indices( 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() + diff --git a/tests/matrix_free/get_functions_multife3.cc b/tests/matrix_free/get_functions_multife3.cc new file mode 100644 index 0000000000..226b08d6b3 --- /dev/null +++ b/tests/matrix_free/get_functions_multife3.cc @@ -0,0 +1,324 @@ +// ------------------------------------------------------------------------ +// +// 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 +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include + +#include +#include +#include + +#include "../tests.h" + + +template +class MatrixFreeTest +{ +public: + using VectorType = std::vector>; + + MatrixFreeTest(const MatrixFree &data_in, + const hp::FECollection &fe_collection) + : data(data_in) + , fe_val0(data.get_dof_handler(0).get_fe(), + Quadrature(data.get_quadrature(0)), + update_values | update_gradients | update_hessians) + , hp_fe_val1(fe_collection, + hp::QCollection(data.get_quadrature(1)), + update_values | update_gradients | update_hessians){}; + + void + operator()(const MatrixFree &data, + VectorType &, + const VectorType &src, + const std::pair &cell_range) const + { + const auto active_fe_index = data.get_cell_range_category(cell_range); + + FEEvaluation fe_eval0(data, 0, 0); + FEEvaluation fe_eval1( + data, 1, 1, 0, active_fe_index); + std::vector reference_values0(fe_eval0.n_q_points); + std::vector> reference_grads0(fe_eval0.n_q_points); + std::vector> reference_hess0(fe_eval0.n_q_points); + std::vector reference_values1(fe_eval1.n_q_points); + std::vector> reference_grads1(fe_eval1.n_q_points); + std::vector> 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::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 == 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 == 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 &data; + mutable FEValues fe_val0; + mutable hp::FEValues hp_fe_val1; + mutable double errors[6], total[6]; +}; + + + +template +void +test() +{ + using number = double; + + Triangulation 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 fe0(fe_degree); + FE_Q fe1(fe_degree + 1); + + DoFHandler dof0(tria); + dof0.distribute_dofs(fe0); + + hp::FECollection fe_collection; + fe_collection.push_back(fe0); + fe_collection.push_back(fe1); + + DoFHandler 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 *> 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> src(dof.size()); + for (unsigned int no = 0; no < dof.size(); ++no) + src[no].reinit(dof[no]->n_dofs()); + + std::vector *> constraints(2); + AffineConstraints constraint0; + DoFTools::make_hanging_node_constraints(*dof[0], constraint0); + constraint0.close(); + constraints[0] = &constraint0; + AffineConstraints 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(); + src[no](i) = entry; + } + + + constraints[0]->distribute(src[0]); + constraints[1]->distribute(src[1]); + MatrixFree mf_data; + { + std::vector> quad; + for (unsigned int no = 0; no < 2; ++no) + quad.push_back(QGauss(fe_degree + 1 + no)); + mf_data.reinit(MappingQ1{}, + dof, + constraints, + quad, + typename MatrixFree::AdditionalData( + MatrixFree::AdditionalData::none)); + } + + MatrixFreeTest 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(); + } +} diff --git a/tests/matrix_free/get_functions_multife3.output b/tests/matrix_free/get_functions_multife3.output new file mode 100644 index 0000000000..e205a840fd --- /dev/null +++ b/tests/matrix_free/get_functions_multife3.output @@ -0,0 +1,57 @@ + +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:: -- 2.39.5