From: Denis Davydov Date: Wed, 27 Sep 2017 14:24:19 +0000 (+0200) Subject: add MGTransferMatrixFree::interpolate_to_mg() to transfer fine-level solution to... X-Git-Tag: v9.0.0-rc1~1012^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c7b37bfa7e0b7c7eb89dd5a249c322cb61b25864;p=dealii.git add MGTransferMatrixFree::interpolate_to_mg() to transfer fine-level solution to a multi-grid vector --- diff --git a/include/deal.II/multigrid/mg_transfer.h b/include/deal.II/multigrid/mg_transfer.h index 04c4629928..d1da88be38 100644 --- a/include/deal.II/multigrid/mg_transfer.h +++ b/include/deal.II/multigrid/mg_transfer.h @@ -368,6 +368,18 @@ public: protected: + /** + * Internal function to perform transfer of residuals or solutions + * basesd on the flag @p solution_transfer. + */ + template + void + copy_to_mg (const DoFHandler &mg_dof, + MGLevelObject > &dst, + const LinearAlgebra::distributed::Vector &src, + const bool solution_transfer) const; + + /** * Internal function to @p fill copy_indices*. Called by derived classes. */ @@ -389,6 +401,13 @@ protected: std::vector > > copy_indices; + + /** + * Same as above, but used to transfer solution vectors. + */ + std::vector > > + solution_copy_indices; + /** * Additional degrees of freedom for the copy_to_mg() function. These are * the ones where the global degree of freedom is locally owned and the @@ -399,6 +418,12 @@ protected: std::vector > > copy_indices_global_mine; + /** + * Same as above, but used to transfer solution vectors. + */ + std::vector > > + solution_copy_indices_global_mine; + /** * Additional degrees of freedom for the copy_from_mg() function. These are * the ones where the level degree of freedom is locally owned and the @@ -409,6 +434,12 @@ protected: std::vector > > copy_indices_level_mine; + /** + * Same as above, but used to transfer solution vectors. + */ + std::vector > > + solution_copy_indices_level_mine; + /** * Stores whether the copy operation from the global to the level vector is * actually a plain copy to the finest level. This means that the grid has @@ -435,11 +466,22 @@ protected: */ mutable LinearAlgebra::distributed::Vector ghosted_global_vector; + /** + * Same as above but used when working with solution vectors. + */ + mutable LinearAlgebra::distributed::Vector solution_ghosted_global_vector; + /** * In the function copy_from_mg, we access all level vectors with certain * ghost entries for inserting the result into a global vector. */ mutable MGLevelObject > ghosted_level_vector; + + /** + * Same as above but used when working with solution vectors. + */ + mutable MGLevelObject > solution_ghosted_level_vector; + }; diff --git a/include/deal.II/multigrid/mg_transfer.templates.h b/include/deal.II/multigrid/mg_transfer.templates.h index 3454f37cc7..d0858185ce 100644 --- a/include/deal.II/multigrid/mg_transfer.templates.h +++ b/include/deal.II/multigrid/mg_transfer.templates.h @@ -369,6 +369,29 @@ MGLevelGlobalTransfer >::copy_to_mg MGLevelObject > &dst, const LinearAlgebra::distributed::Vector &src) const { + copy_to_mg (mg_dof_handler, + dst, + src, + false); +} + + +template +template +void +MGLevelGlobalTransfer >::copy_to_mg +(const DoFHandler &mg_dof_handler, + MGLevelObject > &dst, + const LinearAlgebra::distributed::Vector &src, + const bool solution_transfer) const +{ + LinearAlgebra::distributed::Vector &this_ghosted_global_vector = + solution_transfer ? solution_ghosted_global_vector : ghosted_global_vector; + const std::vector > > &this_copy_indices = + solution_transfer ? solution_copy_indices : copy_indices; + const std::vector > > &this_copy_indices_level_mine = + solution_transfer ? solution_copy_indices_level_mine : copy_indices_level_mine; + AssertIndexRange(dst.max_level(), mg_dof_handler.get_triangulation().n_global_levels()); AssertIndexRange(dst.min_level(), dst.max_level()+1); reinit_vector(mg_dof_handler, component_to_block_map, dst); @@ -390,8 +413,8 @@ MGLevelGlobalTransfer >::copy_to_mg // copy the source vector to the temporary vector that we hold for the // purpose of data exchange - ghosted_global_vector = src; - ghosted_global_vector.update_ghost_values(); + this_ghosted_global_vector = src; + this_ghosted_global_vector.update_ghost_values(); for (unsigned int level=dst.max_level()+1; level != dst.min_level();) { @@ -401,15 +424,15 @@ MGLevelGlobalTransfer >::copy_to_mg LinearAlgebra::distributed::Vector &dst_level = dst[level]; // first copy local unknowns - for (dof_pair_iterator i = copy_indices[level].begin(); - i != copy_indices[level].end(); ++i) - dst_level.local_element(i->second) = ghosted_global_vector.local_element(i->first); + for (dof_pair_iterator i = this_copy_indices[level].begin(); + i != this_copy_indices[level].end(); ++i) + dst_level.local_element(i->second) = this_ghosted_global_vector.local_element(i->first); // Do the same for the indices where the level index is local, but the // global index is not - for (dof_pair_iterator i = copy_indices_level_mine[level].begin(); - i != copy_indices_level_mine[level].end(); ++i) - dst_level.local_element(i->second) = ghosted_global_vector.local_element(i->first); + for (dof_pair_iterator i = this_copy_indices_level_mine[level].begin(); + i != this_copy_indices_level_mine[level].end(); ++i) + dst_level.local_element(i->second) = this_ghosted_global_vector.local_element(i->first); dst_level.compress(VectorOperation::insert); } diff --git a/include/deal.II/multigrid/mg_transfer_internal.h b/include/deal.II/multigrid/mg_transfer_internal.h index 8fa1c86388..7eea670cdc 100644 --- a/include/deal.II/multigrid/mg_transfer_internal.h +++ b/include/deal.II/multigrid/mg_transfer_internal.h @@ -32,13 +32,18 @@ namespace internal /** * Internal function for filling the copy indices from global to level * indices + * + * If @p skip_interface_dofs is false, the mapping will also contain + * DoFs at the interface between levels. This is desirable when + * transfering solution vectors instead of residuals. */ template void fill_copy_indices(const dealii::DoFHandler &mg_dof, const MGConstrainedDoFs *mg_constrained_dofs, std::vector > > ©_indices, std::vector > > ©_indices_global_mine, - std::vector > > ©_indices_level_mine); + std::vector > > ©_indices_level_mine, + const bool skip_interface_dofs = true); diff --git a/include/deal.II/multigrid/mg_transfer_matrix_free.h b/include/deal.II/multigrid/mg_transfer_matrix_free.h index ddb94109c9..7b0328dc64 100644 --- a/include/deal.II/multigrid/mg_transfer_matrix_free.h +++ b/include/deal.II/multigrid/mg_transfer_matrix_free.h @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -126,6 +127,18 @@ public: LinearAlgebra::distributed::Vector &dst, const LinearAlgebra::distributed::Vector &src) const; + /** + * Restrict fine-mesh field @p src to each multigrid level in @p mg_dof and + * store the result in @p dst. + * + * If @p dst is empty or has incorrect locally owned size, it will be + * resized to locally relevant vectors on each level. + */ + template + void interpolate_to_mg(const DoFHandler< dim, spacedim > &mg_dof, + MGLevelObject> &dst, + const LinearAlgebra::distributed::Vector &src) const; + /** * Finite element does not provide prolongation matrices. */ @@ -372,6 +385,98 @@ private: //------------------------ templated functions ------------------------- #ifndef DOXYGEN + +template +template +void +MGTransferMatrixFree:: +interpolate_to_mg(const DoFHandler &dof_handler, + MGLevelObject> &dst, + const LinearAlgebra::distributed::Vector &src) const +{ + const unsigned int min_level = dst.min_level(); + const unsigned int max_level = dst.max_level(); + + Assert (max_level == dof_handler.get_triangulation().n_global_levels()-1, + ExcDimensionMismatch(max_level,dof_handler.get_triangulation().n_global_levels()-1)); + + const parallel::Triangulation *p_tria = + (dynamic_cast*> + (&dof_handler.get_triangulation())); + MPI_Comm mpi_communicator = p_tria != nullptr ? p_tria->get_communicator() : MPI_COMM_SELF; + + // resize the dst vector if it's empty or has incorrect size + MGLevelObject relevant_dofs(min_level,max_level); + for (unsigned int level = min_level; level <= max_level; ++level) + { + DoFTools::extract_locally_relevant_level_dofs(dof_handler, level, + relevant_dofs[level]); + if (dst[level].size() != dof_handler.locally_owned_mg_dofs(level).size() || + dst[level].local_size() != dof_handler.locally_owned_mg_dofs(level).n_elements()) + dst[level].reinit(dof_handler.locally_owned_mg_dofs(level), + relevant_dofs[level], + mpi_communicator); + } + + const FiniteElement &fe = dof_handler.get_fe(); + + // copy fine level vector to active cells in MG hierarchy + this->copy_to_mg (dof_handler, dst, src, true); + + // FIXME: maybe need to store hanging nodes constraints per level? + // MGConstrainedDoFs does NOT keep this info right now, only periodicity constraints... + dst[max_level].update_ghost_values(); + // do the transfer from level to level-1: + for (unsigned int level=max_level; level > min_level; --level) + { + // auxiliary vector which always has ghost elements + LinearAlgebra::distributed::Vector + ghosted_vector(dof_handler.locally_owned_mg_dofs(level), + relevant_dofs[level], + mpi_communicator); + ghosted_vector = dst[level]; + ghosted_vector.update_ghost_values(); + + std::vector dof_values_coarse(fe.dofs_per_cell); + Vector dof_values_fine(fe.dofs_per_cell); + Vector tmp(fe.dofs_per_cell); + std::vector dof_indices(fe.dofs_per_cell); + typename DoFHandler::cell_iterator cell=dof_handler.begin(level-1); + typename DoFHandler::cell_iterator endc=dof_handler.end(level-1); + for ( ; cell != endc; ++cell) + if (cell->is_locally_owned_on_level()) + { + // if we get to a cell without children (== active), we can + // skip it as there values should be already set by the + // equivalent of copy_to_mg() + if (!cell->has_children()) + continue; + + std::fill(dof_values_coarse.begin(), dof_values_coarse.end(), 0.); + for (unsigned int child=0; childn_children(); ++child) + { + cell->child(child)->get_mg_dof_indices(dof_indices); + for (unsigned int i=0; irefinement_case()).vmult (tmp, dof_values_fine); + for (unsigned int i=0; iget_mg_dof_indices(dof_indices); + for (unsigned int i=0; i template void diff --git a/source/multigrid/mg_level_global_transfer.cc b/source/multigrid/mg_level_global_transfer.cc index 4676b5a27e..49603c7b09 100644 --- a/source/multigrid/mg_level_global_transfer.cc +++ b/source/multigrid/mg_level_global_transfer.cc @@ -143,6 +143,110 @@ MGLevelGlobalTransfer::memory_consumption () const /* ------------------ MGLevelGlobalTransfer ----------------- */ +namespace +{ + template + void fill_internal(const DoFHandler &mg_dof, + SmartPointer > > mg_constrained_dofs, + const MPI_Comm mpi_communicator, + const bool transfer_solution_vectors, + std::vector > > ©_indices, + std::vector > > ©_indices_global_mine, + std::vector > > ©_indices_level_mine, + LinearAlgebra::distributed::Vector &ghosted_global_vector, + MGLevelObject > &ghosted_level_vector) + { + // first go to the usual routine... + std::vector > > + my_copy_indices; + std::vector > > + my_copy_indices_global_mine; + std::vector > > + my_copy_indices_level_mine; + + internal::MGTransfer::fill_copy_indices(mg_dof, mg_constrained_dofs, my_copy_indices, + my_copy_indices_global_mine, my_copy_indices_level_mine, + !transfer_solution_vectors); + + // get all degrees of freedom that we need read access to in copy_to_mg + // and copy_from_mg, respectively. We fill an IndexSet once on each level + // (for the global_mine indices accessing remote level indices) and once + // globally (for the level_mine indices accessing remote global indices). + + // the variables index_set and level_index_set are going to define the + // ghost indices of the respective vectors (due to construction, these are + // precisely the indices that we need) + + IndexSet index_set(mg_dof.locally_owned_dofs().size()); + std::vector accessed_indices; + ghosted_level_vector.resize(0, mg_dof.get_triangulation().n_global_levels()-1); + std::vector level_index_set(mg_dof.get_triangulation().n_global_levels()); + for (unsigned int l=0; l accessed_level_indices; + for (unsigned int i=0; i + (global_partitioner.global_to_local(my_copy_indices[level][i].first), + level_partitioner.global_to_local(my_copy_indices[level][i].second)); + + // remote-owned case: the locally owned indices for the level and the + // ghost dofs for the global indices set the local index + copy_indices_level_mine[level]. + resize(my_copy_indices_level_mine[level].size()); + for (unsigned int i=0; i + (global_partitioner.global_to_local(my_copy_indices_level_mine[level][i].first), + level_partitioner.global_to_local(my_copy_indices_level_mine[level][i].second)); + + // owned-remote case: the locally owned indices for the global dofs + // and the ghost dofs for the level indices set the local index + copy_indices_global_mine[level]. + resize(my_copy_indices_global_mine[level].size()); + for (unsigned int i=0; i + (global_partitioner.global_to_local(my_copy_indices_global_mine[level][i].first), + level_partitioner.global_to_local(my_copy_indices_global_mine[level][i].second)); + } + } +} template template @@ -150,99 +254,32 @@ void MGLevelGlobalTransfer >::fill_and_communicate_copy_indices (const DoFHandler &mg_dof) { - // first go to the usual routine... - std::vector > > - my_copy_indices; - std::vector > > - my_copy_indices_global_mine; - std::vector > > - my_copy_indices_level_mine; - - internal::MGTransfer::fill_copy_indices(mg_dof, mg_constrained_dofs, my_copy_indices, - my_copy_indices_global_mine, my_copy_indices_level_mine); - - // get all degrees of freedom that we need read access to in copy_to_mg - // and copy_from_mg, respectively. We fill an IndexSet once on each level - // (for the global_mine indices accessing remote level indices) and once - // globally (for the level_mine indices accessing remote global indices). - - // the variables index_set and level_index_set are going to define the - // ghost indices of the respective vectors (due to construction, these are - // precisely the indices that we need) + const parallel::Triangulation *ptria = dynamic_cast *> (&mg_dof.get_triangulation()); const MPI_Comm mpi_communicator = ptria != nullptr ? ptria->get_communicator() : MPI_COMM_SELF; - IndexSet index_set(mg_dof.locally_owned_dofs().size()); - std::vector accessed_indices; - ghosted_level_vector.resize(0, mg_dof.get_triangulation().n_global_levels()-1); - std::vector level_index_set(mg_dof.get_triangulation().n_global_levels()); - for (unsigned int l=0; l accessed_level_indices; - for (unsigned int i=0; icopy_indices.resize(mg_dof.get_triangulation().n_global_levels()); - this->copy_indices_level_mine.resize(mg_dof.get_triangulation().n_global_levels()); - this->copy_indices_global_mine.resize(mg_dof.get_triangulation().n_global_levels()); - for (unsigned int level=0; levelcopy_indices[level].resize(my_copy_indices[level].size()); - for (unsigned int i=0; icopy_indices[level][i] = - std::pair - (global_partitioner.global_to_local(my_copy_indices[level][i].first), - level_partitioner.global_to_local(my_copy_indices[level][i].second)); - - // remote-owned case: the locally owned indices for the level and the - // ghost dofs for the global indices set the local index - this->copy_indices_level_mine[level]. - resize(my_copy_indices_level_mine[level].size()); - for (unsigned int i=0; icopy_indices_level_mine[level][i] = - std::pair - (global_partitioner.global_to_local(my_copy_indices_level_mine[level][i].first), - level_partitioner.global_to_local(my_copy_indices_level_mine[level][i].second)); - - // owned-remote case: the locally owned indices for the global dofs - // and the ghost dofs for the level indices set the local index - this->copy_indices_global_mine[level]. - resize(my_copy_indices_global_mine[level].size()); - for (unsigned int i=0; icopy_indices_global_mine[level][i] = - std::pair - (global_partitioner.global_to_local(my_copy_indices_global_mine[level][i].first), - level_partitioner.global_to_local(my_copy_indices_global_mine[level][i].second)); - } + fill_internal(mg_dof, + mg_constrained_dofs, + mpi_communicator, + false, + this->copy_indices, + this->copy_indices_global_mine, + this->copy_indices_level_mine, + ghosted_global_vector, + ghosted_level_vector); + + fill_internal(mg_dof, + mg_constrained_dofs, + mpi_communicator, + true, + this->solution_copy_indices, + this->solution_copy_indices_global_mine, + this->solution_copy_indices_level_mine, + solution_ghosted_global_vector, + solution_ghosted_level_vector); perform_plain_copy = this->copy_indices.back().size() == mg_dof.locally_owned_dofs().n_elements(); @@ -271,6 +308,8 @@ MGLevelGlobalTransfer >::fill_and_com { ghosted_global_vector.reinit(0); ghosted_level_vector.resize(0, 0); + solution_ghosted_global_vector.reinit(0); + solution_ghosted_level_vector.resize(0, 0); } } diff --git a/source/multigrid/mg_level_global_transfer.inst.in b/source/multigrid/mg_level_global_transfer.inst.in index 64a1f9bd9a..f6b02513be 100644 --- a/source/multigrid/mg_level_global_transfer.inst.in +++ b/source/multigrid/mg_level_global_transfer.inst.in @@ -45,6 +45,13 @@ for (deal_II_dimension : DIMENSIONS) template void MGLevelGlobalTransfer< LinearAlgebra::distributed::Vector >::fill_and_communicate_copy_indices( const DoFHandler &mg_dof); + template + void MGLevelGlobalTransfer< LinearAlgebra::distributed::Vector >::copy_to_mg ( + const DoFHandler&, MGLevelObject>&, const LinearAlgebra::distributed::Vector&, const bool) const; + template + void MGLevelGlobalTransfer< LinearAlgebra::distributed::Vector >::copy_to_mg ( + const DoFHandler&, MGLevelObject>&, const LinearAlgebra::distributed::Vector&, const bool) const; + } for (deal_II_dimension : DIMENSIONS; S1, S2 : REAL_SCALARS) diff --git a/source/multigrid/mg_transfer_internal.cc b/source/multigrid/mg_transfer_internal.cc index 9d157ddf49..b9cb4fb7d1 100644 --- a/source/multigrid/mg_transfer_internal.cc +++ b/source/multigrid/mg_transfer_internal.cc @@ -53,19 +53,21 @@ namespace internal {} }; - // Internal function for filling the copy indices from global to level - // indices + + + template void fill_copy_indices(const dealii::DoFHandler &mg_dof, const MGConstrainedDoFs *mg_constrained_dofs, std::vector > > ©_indices, std::vector > > ©_indices_global_mine, - std::vector > > ©_indices_level_mine) + std::vector > > ©_indices_level_mine, + const bool skip_interface_dofs) { // Now we are filling the variables copy_indices*, which are essentially // maps from global to mgdof for each level stored as a std::vector of // pairs. We need to split this map on each level depending on the - // ownership of the global and mgdof, so that we later not access + // ownership of the global and mgdof, so that we later do not access // non-local elements in copy_to/from_mg. // We keep track in the bitfield dof_touched which global dof has been // processed already (on the current level). This is the same as the @@ -114,9 +116,11 @@ namespace internal for (unsigned int i=0; iat_refinement_edge(level, level_dof_indices[i])) continue; + types::global_dof_index global_idx = globally_relevant.index_within_set(global_dof_indices[i]); //skip if we did this global dof already (on this or a coarser level) if (dof_touched[global_idx]) diff --git a/source/multigrid/mg_transfer_internal.inst.in b/source/multigrid/mg_transfer_internal.inst.in index d1b848dd9b..81ac91a26e 100644 --- a/source/multigrid/mg_transfer_internal.inst.in +++ b/source/multigrid/mg_transfer_internal.inst.in @@ -30,7 +30,8 @@ for (deal_II_dimension : DIMENSIONS; const MGConstrainedDoFs*, std::vector > > &, std::vector > > &, - std::vector > > &); + std::vector > > &, + const bool); #endif \} \} diff --git a/tests/matrix_free/interpolate_to_mg.cc b/tests/matrix_free/interpolate_to_mg.cc new file mode 100644 index 0000000000..ce2cca5958 --- /dev/null +++ b/tests/matrix_free/interpolate_to_mg.cc @@ -0,0 +1,250 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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. +// +// --------------------------------------------------------------------- + + +// test MGTransferMatrixFree::interpolate_to_mg() +// for a scalar field + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + + + + +using namespace dealii; + + +template +class SimpleField : public Function +{ +public: + SimpleField() : + Function(1) + {} + + double value (const Point &p, + const unsigned int component = 0) const + { + (void)component; + return p[0]*2. + p[1] - 10.; + } +}; + + +template +void test (const unsigned int n_glob_ref=2, const unsigned int n_ref = 0) +{ + SimpleField function; + + deallog << "dim=" << dim << std::endl; + MPI_Comm mpi_communicator(MPI_COMM_WORLD); + unsigned int myid = Utilities::MPI::this_mpi_process (mpi_communicator); + unsigned int numproc = Utilities::MPI::n_mpi_processes (mpi_communicator); + + deallog << "numproc=" << numproc << std::endl; + + parallel::distributed::Triangulation + triangulation (mpi_communicator, + Triangulation::limit_level_difference_at_vertices, + parallel::distributed::Triangulation::construct_multigrid_hierarchy); + GridGenerator::subdivided_hyper_cube (triangulation, 5, -1, 1); + // now do some refinement + triangulation.refine_global(n_glob_ref); + for (unsigned int ref=0; ref::active_cell_iterator cell=triangulation.begin_active(); + cell != triangulation.end(); ++cell) + if (cell->is_locally_owned() && + (cell->center().norm() < 0.5 && (cell->level() < 5 || + cell->center().norm() > 0.45) + || + (dim == 2 && cell->center().norm() > 1.2))) + cell->set_refine_flag(); + triangulation.execute_coarsening_and_refinement(); + } + + FE_Q fe (fe_degree); + DoFHandler dof_handler (triangulation); + + dof_handler.distribute_dofs(fe); + dof_handler.distribute_mg_dofs(); + + IndexSet locally_owned_dofs, locally_relevant_dofs; + locally_owned_dofs = dof_handler.locally_owned_dofs (); + DoFTools::extract_locally_relevant_dofs (dof_handler,locally_relevant_dofs); + + // constraints: + ConstraintMatrix constraints; + constraints.reinit (locally_relevant_dofs); + DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + + // interpolate: + LinearAlgebra::distributed::Vector fine_projection; + fine_projection.reinit(locally_owned_dofs, + locally_relevant_dofs, + mpi_communicator); + VectorTools::project (dof_handler, + constraints, + QGauss(n_q_points+2), + function, + fine_projection); + fine_projection.update_ghost_values(); + + // output for debug purposes: + if (true) + { + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (fine_projection, "projection"); + + Vector subdomain (triangulation.n_active_cells()); + for (unsigned int i=0; i mg_transfer(mg_constrained_dofs); + mg_transfer.build(dof_handler); + + // now the core of the test: + const unsigned int max_level = dof_handler.get_triangulation().n_global_levels()-1; + const unsigned int min_level = 0; + MGLevelObject> level_projection(min_level, max_level); + mg_transfer.interpolate_to_mg(dof_handler,level_projection, fine_projection); + + // now go through all GMG levels and make sure FE field can represent + // analytic function exactly: + QGauss quadrature(n_q_points); + std::vector q_values(quadrature.size()); + + FEValues fe_values(fe, quadrature, update_values | update_quadrature_points); + for (unsigned int level = max_level+1; level != min_level;) + { + --level; + + std::vector dof_indices(fe.dofs_per_cell); + typename DoFHandler::cell_iterator cell=dof_handler.begin(level); + typename DoFHandler::cell_iterator endc=dof_handler.end(level); + for ( ; cell != endc; ++cell) + if (cell->is_locally_owned_on_level()) + { + fe_values.reinit(cell); + cell->get_mg_dof_indices(dof_indices); + fe_values.get_function_values(level_projection[level], + dof_indices, + q_values); + + const std::vector> &q_points = fe_values.get_quadrature_points(); + + for (unsigned int q = 0; q < q_points.size(); ++q) + { + const double diff = q_values[q] - function.value(q_points[q]); + if (std::abs(diff) > 1e-10) + { + std::cout<<"dofs: "; + for (const auto i : dof_indices) + std::cout << i << " "; + std::cout << std::endl << "values: "; + std::vector local_values(dof_indices.size()); + level_projection[level].extract_subvector_to(dof_indices, local_values); + for (const auto v : local_values) + std::cout << v << " "; + std::cout << std::endl + << "val(q)=" << q_values[q] << std::endl; + std::cout << "MGTransfer indices:" << std::endl; + mg_transfer.print_indices(std::cout); + AssertThrow (false, + ExcMessage("Level " + std::to_string(level) + + " Diff " + std::to_string(diff) + + " center_x " + std::to_string(cell->center()[0]) + + " center_y " + std::to_string(cell->center()[1]) + )); + } + } + } + } + + deallog << "Ok" << std::endl; +} + + + +int main (int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + MPILogInitAll log; + + test<2> (); + test<2,1> (0,1); +} + + diff --git a/tests/matrix_free/interpolate_to_mg.with_p4est=true.mpirun=1.output b/tests/matrix_free/interpolate_to_mg.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..dad0ed2a3a --- /dev/null +++ b/tests/matrix_free/interpolate_to_mg.with_p4est=true.mpirun=1.output @@ -0,0 +1,7 @@ + +DEAL:0::dim=2 +DEAL:0::numproc=1 +DEAL:0::Ok +DEAL:0::dim=2 +DEAL:0::numproc=1 +DEAL:0::Ok diff --git a/tests/matrix_free/interpolate_to_mg.with_p4est=true.mpirun=3.output b/tests/matrix_free/interpolate_to_mg.with_p4est=true.mpirun=3.output new file mode 100644 index 0000000000..a4a680237e --- /dev/null +++ b/tests/matrix_free/interpolate_to_mg.with_p4est=true.mpirun=3.output @@ -0,0 +1,23 @@ + +DEAL:0::dim=2 +DEAL:0::numproc=3 +DEAL:0::Ok +DEAL:0::dim=2 +DEAL:0::numproc=3 +DEAL:0::Ok + +DEAL:1::dim=2 +DEAL:1::numproc=3 +DEAL:1::Ok +DEAL:1::dim=2 +DEAL:1::numproc=3 +DEAL:1::Ok + + +DEAL:2::dim=2 +DEAL:2::numproc=3 +DEAL:2::Ok +DEAL:2::dim=2 +DEAL:2::numproc=3 +DEAL:2::Ok +