From 5b6f1a6210d2a0ea6da19bf2e0a1accefb45b17e Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Wed, 7 Feb 2018 14:24:37 +0100 Subject: [PATCH] extend VectorTools::project() to take fe_component --- doc/news/changes/minor/20180207DenisDavydov | 4 + include/deal.II/numerics/vector_tools.h | 10 ++- .../deal.II/numerics/vector_tools.templates.h | 63 ++++++------- .../vector_tools_project_qpmf.inst.in | 6 +- tests/numerics/project_parallel_qpmf_common.h | 90 +++++++++++++------ tests/numerics/project_parallel_qpmf_q_03.cc | 43 +++++++++ ..._qpmf_q_03.with_p4est=true.mpirun=3.output | 35 ++++++++ 7 files changed, 191 insertions(+), 60 deletions(-) create mode 100644 doc/news/changes/minor/20180207DenisDavydov create mode 100644 tests/numerics/project_parallel_qpmf_q_03.cc create mode 100644 tests/numerics/project_parallel_qpmf_q_03.with_p4est=true.mpirun=3.output diff --git a/doc/news/changes/minor/20180207DenisDavydov b/doc/news/changes/minor/20180207DenisDavydov new file mode 100644 index 0000000000..4c1769f545 --- /dev/null +++ b/doc/news/changes/minor/20180207DenisDavydov @@ -0,0 +1,4 @@ +Extend VectorTools::project() function with MatrixFree quadrature data to +optionally take the finite element component index. +
+(Denis Davydov, 2018/02/07) diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index 55eb3a4331..d0c121b2cd 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -919,13 +919,18 @@ namespace VectorTools * @endcode * where qp_data is a an object of type Table<2, VectorizedArray >, * which stores quadrature point data. + * + * @p fe_component allow to additionally specify which component of @p data + * to use in case it was constructed with an std::vector*> + * in will be internally used in constructor of FEEvaluation object. */ template void project (std::shared_ptr > data, const ConstraintMatrix &constraints, const unsigned int n_q_points_1d, const std::function< VectorizedArray (const unsigned int, const unsigned int)> &func, - VectorType &vec_result); + VectorType &vec_result, + const unsigned int fe_component = 0); /** * Same as above but for n_q_points_1d = matrix_free.get_dof_handler().get_fe().degree+1. @@ -934,7 +939,8 @@ namespace VectorTools void project (std::shared_ptr > data, const ConstraintMatrix &constraints, const std::function< VectorizedArray (const unsigned int, const unsigned int)> &func, - VectorType &vec_result); + VectorType &vec_result, + const unsigned int fe_component = 0); /** * Compute Dirichlet boundary conditions. This function makes up a map of diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 1e68a740e0..d07d0bbdce 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -1386,9 +1386,10 @@ namespace VectorTools void project_parallel (std::shared_ptr > matrix_free, const ConstraintMatrix &constraints, const std::function< VectorizedArray (const unsigned int, const unsigned int)> &func, - VectorType &vec_result) + VectorType &vec_result, + const unsigned int fe_component) { - const DoFHandler &dof = matrix_free->get_dof_handler(); + const DoFHandler &dof = matrix_free->get_dof_handler(fe_component); typedef typename VectorType::value_type Number; Assert (dof.get_fe(0).n_components() == 1, @@ -1402,20 +1403,20 @@ namespace VectorTools typedef MatrixFreeOperators::MassOperator > MatrixType; MatrixType mass_matrix; - mass_matrix.initialize(matrix_free); + mass_matrix.initialize(matrix_free, {{fe_component}}); mass_matrix.compute_diagonal(); typedef LinearAlgebra::distributed::Vector LocalVectorType; LocalVectorType vec, rhs, inhomogeneities; - matrix_free->initialize_dof_vector(vec); - matrix_free->initialize_dof_vector(rhs); - matrix_free->initialize_dof_vector(inhomogeneities); + matrix_free->initialize_dof_vector(vec,fe_component); + matrix_free->initialize_dof_vector(rhs,fe_component); + matrix_free->initialize_dof_vector(inhomogeneities,fe_component); constraints.distribute(inhomogeneities); inhomogeneities*=-1.; // assemble right hand side: { - FEEvaluation fe_eval(*matrix_free); + FEEvaluation fe_eval(*matrix_free, fe_component); const unsigned int n_cells = matrix_free->n_macro_cells(); const unsigned int n_q_points = fe_eval.n_q_points; @@ -1489,28 +1490,28 @@ namespace VectorTools const ConstraintMatrix &constraints, const unsigned int n_q_points_1d, const std::function< VectorizedArray (const unsigned int, const unsigned int)> &func, - VectorType &vec_result) + VectorType &vec_result, + const unsigned int fe_component) { - (void) n_q_points_1d; - const unsigned int fe_degree = matrix_free->get_dof_handler().get_fe().degree; - - Assert (fe_degree+1 == n_q_points_1d, - ExcNotImplemented()); + const unsigned int fe_degree = matrix_free->get_dof_handler(fe_component).get_fe().degree; - switch (fe_degree) - { - case 1: - project_parallel (matrix_free,constraints,func,vec_result); - break; - case 2: - project_parallel (matrix_free,constraints,func,vec_result); - break; - case 3: - project_parallel (matrix_free,constraints,func,vec_result); - break; - default: - project_parallel (matrix_free,constraints,func,vec_result); - } + if (fe_degree+1 == n_q_points_1d) + switch (fe_degree) + { + case 1: + project_parallel (matrix_free,constraints,func,vec_result,fe_component); + break; + case 2: + project_parallel (matrix_free,constraints,func,vec_result,fe_component); + break; + case 3: + project_parallel (matrix_free,constraints,func,vec_result,fe_component); + break; + default: + project_parallel (matrix_free,constraints,func,vec_result,fe_component); + } + else + project_parallel (matrix_free,constraints,func,vec_result,fe_component); } @@ -1519,13 +1520,15 @@ namespace VectorTools void project (std::shared_ptr > matrix_free, const ConstraintMatrix &constraints, const std::function< VectorizedArray (const unsigned int, const unsigned int)> &func, - VectorType &vec_result) + VectorType &vec_result, + const unsigned int fe_component) { project (matrix_free, constraints, - matrix_free->get_dof_handler().get_fe().degree+1, + matrix_free->get_dof_handler(fe_component).get_fe().degree+1, func, - vec_result); + vec_result, + fe_component); } diff --git a/source/numerics/vector_tools_project_qpmf.inst.in b/source/numerics/vector_tools_project_qpmf.inst.in index 3b9c0968c6..4db6e507cb 100644 --- a/source/numerics/vector_tools_project_qpmf.inst.in +++ b/source/numerics/vector_tools_project_qpmf.inst.in @@ -24,7 +24,8 @@ for (VEC: REAL_NONBLOCK_VECTORS; deal_II_dimension : DIMENSIONS) std::shared_ptr > matrix_free, const ConstraintMatrix &constraints, const std::function< VectorizedArray (const unsigned int, const unsigned int)> &, - VEC &); + VEC &, + const unsigned int); template void project( @@ -32,7 +33,8 @@ for (VEC: REAL_NONBLOCK_VECTORS; deal_II_dimension : DIMENSIONS) const ConstraintMatrix &constraints, const unsigned int, const std::function< VectorizedArray (const unsigned int, const unsigned int)> &, - VEC &); + VEC &, + const unsigned int); \} diff --git a/tests/numerics/project_parallel_qpmf_common.h b/tests/numerics/project_parallel_qpmf_common.h index 79c1ec6143..d83496e977 100644 --- a/tests/numerics/project_parallel_qpmf_common.h +++ b/tests/numerics/project_parallel_qpmf_common.h @@ -48,6 +48,9 @@ #include #include +/** + * A function of polynomial degree @p q with n_components. + */ template class F : public Function { @@ -90,27 +93,44 @@ private: }; +/** + * A general function to perform projection of polynomial function of degrees [0,p] + * onto the FE space given by @p fe and @p triangulation . + */ template void do_project (const parallel::distributed::Triangulation &triangulation, - const FiniteElement &fe, - const unsigned int p) + const std::vector*> &fes, + const unsigned int p, + const unsigned int fe_index = 0) { - AssertThrow(fe.n_components() ==1, + // we only project scalar value function + AssertThrow(fes[fe_index]->n_components() ==1, ExcNotImplemented()); - DoFHandler dof_handler(triangulation); - dof_handler.distribute_dofs (fe); deallog << "n_cells=" << triangulation.n_global_active_cells() << std::endl; - deallog << "n_dofs=" << dof_handler.n_dofs() << std::endl; - IndexSet locally_relevant_dofs; - DoFTools::extract_locally_relevant_dofs (dof_handler, - locally_relevant_dofs); + std::vector>> dof_handlers(fes.size()); + std::vector locally_relevant_dofs(fes.size()); + std::vector constraints(fes.size()); - ConstraintMatrix constraints; - constraints.reinit (locally_relevant_dofs); - DoFTools::make_hanging_node_constraints (dof_handler, constraints); - constraints.close (); + std::vector*> dof_handlers_mf(fes.size()); + std::vector constraints_mf(fes.size()); + for (unsigned int i = 0; i < fes.size(); ++i) + { + dof_handlers[i] = std::make_shared>(triangulation); + dof_handlers[i]->distribute_dofs (*fes[i]); + deallog << "n_dofs=" << dof_handlers[i]->n_dofs() << std::endl; + + DoFTools::extract_locally_relevant_dofs (*dof_handlers[i], + locally_relevant_dofs[i]); + + constraints[i].reinit (locally_relevant_dofs[i]); + DoFTools::make_hanging_node_constraints (*dof_handlers[i], constraints[i]); + constraints[i].close (); + + dof_handlers_mf[i] = dof_handlers[i].get(); + constraints_mf[i] = &constraints[i]; + } QGauss<1> quadrature_formula_1d(n_q_points_1d); @@ -120,16 +140,16 @@ void do_project (const parallel::distributed::Triangulation &triangulation, additional_data.tasks_parallel_scheme = MatrixFree::AdditionalData::partition_color; additional_data.mapping_update_flags = update_values | update_JxW_values | update_quadrature_points; std::shared_ptr > data(new MatrixFree ()); - data->reinit (dof_handler, constraints, quadrature_formula_1d, additional_data); + data->reinit (dof_handlers_mf, constraints_mf, quadrature_formula_1d, additional_data); for (unsigned int q=0; q<=p; ++q) { // setup quadrature data: - F function(q, fe.n_components()); + F function(q, fes[fe_index]->n_components()); // initialize a quadrature data { - FEEvaluation fe_eval(*data); + FEEvaluation fe_eval(*data,fe_index); const unsigned int n_cells = data->n_macro_cells(); const unsigned int n_q_points = fe_eval.n_q_points; @@ -143,13 +163,14 @@ void do_project (const parallel::distributed::Triangulation &triangulation, } LinearAlgebra::distributed::Vector field; - data->initialize_dof_vector(field); + data->initialize_dof_vector(field,fe_index); VectorTools::project > (data, - constraints, + constraints[fe_index], n_q_points_1d, [=] (const unsigned int cell, const unsigned int q) -> VectorizedArray { return qp_data(cell,q); }, - field); + field, + fe_index); field.update_ghost_values(); @@ -159,16 +180,16 @@ void do_project (const parallel::distributed::Triangulation &triangulation, double L2_norm = 0.; { QGauss quadrature_formula_error(std::max(p,q)+1); - FEValues fe_values (fe, quadrature_formula_error, + FEValues fe_values (*fes[fe_index], quadrature_formula_error, update_values | update_quadrature_points | update_JxW_values); - const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int dofs_per_cell = fes[fe_index]->dofs_per_cell; const unsigned int n_q_points = quadrature_formula_error.size(); std::vector values (n_q_points); typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); + cell = dof_handlers[fe_index]->begin_active(), + endc = dof_handlers[fe_index]->end(); for (; cell!=endc; ++cell) if (cell->is_locally_owned()) { @@ -186,7 +207,7 @@ void do_project (const parallel::distributed::Triangulation &triangulation, L2_norm = Utilities::MPI::sum(L2_norm, MPI_COMM_WORLD); L2_norm = std::sqrt(L2_norm); - deallog << fe.get_name() << ", P_" << q + deallog << fes[fe_index]->get_name() << ", P_" << q << ", rel. error=" << L2_norm / field_l2_norm << std::endl; @@ -212,7 +233,7 @@ void test_no_hanging_nodes (const FiniteElement &fe, GridGenerator::hyper_cube (triangulation); triangulation.refine_global (3); - do_project (triangulation, fe, p); + do_project (triangulation, {{&fe}}, p); } @@ -230,6 +251,23 @@ void test_with_hanging_nodes (const FiniteElement &fe, triangulation.execute_coarsening_and_refinement (); triangulation.refine_global (1); - do_project (triangulation, fe, p); + do_project (triangulation, {{&fe}}, p); } + +// same as above but for multiple fes +template +void test_with_hanging_nodes (const std::vector*> &fes, + const unsigned int p, + const unsigned int fe_index) +{ + parallel::distributed::Triangulation triangulation(MPI_COMM_WORLD); + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (2); + if (triangulation.begin_active()->is_locally_owned()) + triangulation.begin_active()->set_refine_flag(); + triangulation.execute_coarsening_and_refinement (); + triangulation.refine_global (1); + + do_project (triangulation, fes, p, fe_index); +} diff --git a/tests/numerics/project_parallel_qpmf_q_03.cc b/tests/numerics/project_parallel_qpmf_q_03.cc new file mode 100644 index 0000000000..fadf83283b --- /dev/null +++ b/tests/numerics/project_parallel_qpmf_q_03.cc @@ -0,0 +1,43 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// check VectorTools::project_parallel() for matrix-free quadrature data for +// FE_Q on a mesh with hanging nodes and specified fe_component + +#include "../tests.h" +#include "project_parallel_qpmf_common.h" + +template +void test () +{ + FE_Q fe1(1); + FE_Q fe2(2); + FE_Q fe4(4); + test_with_hanging_nodes<1, 2,dim> ({{&fe1,&fe2}}, 1, 0); // take first + test_with_hanging_nodes<2, 3,dim> ({{&fe1,&fe2}}, 2, 1); // take second + test_with_hanging_nodes<2, 6,dim> ({{&fe2,&fe4}}, 2, 0); // take first with non-standard n_q_points_1d +} + +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, + numbers::invalid_unsigned_int); + + initlog(); + + test<2>(); + test<3>(); +} diff --git a/tests/numerics/project_parallel_qpmf_q_03.with_p4est=true.mpirun=3.output b/tests/numerics/project_parallel_qpmf_q_03.with_p4est=true.mpirun=3.output new file mode 100644 index 0000000000..9fd643d60d --- /dev/null +++ b/tests/numerics/project_parallel_qpmf_q_03.with_p4est=true.mpirun=3.output @@ -0,0 +1,35 @@ + +DEAL::n_cells=76 +DEAL::n_dofs=97 +DEAL::n_dofs=349 +DEAL::FE_Q<2>(1), P_0, rel. error=0 +DEAL::FE_Q<2>(1), P_1, rel. error=0 +DEAL::n_cells=76 +DEAL::n_dofs=97 +DEAL::n_dofs=349 +DEAL::FE_Q<2>(2), P_0, rel. error=0 +DEAL::FE_Q<2>(2), P_1, rel. error=0 +DEAL::FE_Q<2>(2), P_2, rel. error=0 +DEAL::n_cells=76 +DEAL::n_dofs=349 +DEAL::n_dofs=1309 +DEAL::FE_Q<2>(2), P_0, rel. error=0 +DEAL::FE_Q<2>(2), P_1, rel. error=0 +DEAL::FE_Q<2>(2), P_2, rel. error=0 +DEAL::n_cells=568 +DEAL::n_dofs=827 +DEAL::n_dofs=5559 +DEAL::FE_Q<3>(1), P_0, rel. error=0 +DEAL::FE_Q<3>(1), P_1, rel. error=0 +DEAL::n_cells=568 +DEAL::n_dofs=827 +DEAL::n_dofs=5559 +DEAL::FE_Q<3>(2), P_0, rel. error=0 +DEAL::FE_Q<3>(2), P_1, rel. error=0 +DEAL::FE_Q<3>(2), P_2, rel. error=0 +DEAL::n_cells=568 +DEAL::n_dofs=5559 +DEAL::n_dofs=40319 +DEAL::FE_Q<3>(2), P_0, rel. error=0 +DEAL::FE_Q<3>(2), P_1, rel. error=0 +DEAL::FE_Q<3>(2), P_2, rel. error=0 -- 2.39.5