From: Pasquale Africa Date: Thu, 24 Oct 2019 11:33:25 +0000 (+0000) Subject: Final implementation + update documentation and test case X-Git-Tag: v9.2.0-rc1~947^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4e57dfd388f9a7ee8e08134573f8a0e3983b5636;p=dealii.git Final implementation + update documentation and test case Co-authored-by: Doug Shi-Dong --- diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index 341e3429c5..fafa5d3d7b 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -641,7 +641,7 @@ namespace VectorTools * Interpolate different finite element spaces. The interpolation of vector * @p data_1 (which is assumed to be ghosted, see @ref GlossGhostedVector) * is executed from the FE space represented by @p dof_1 - * to the vector @p data_2 (which must be non-ghosted) on FE space @p dof_2. + * to the vector @p data_2 on FE space @p dof_2. * The interpolation on each cell is represented by the matrix @p transfer. * Curved boundaries are neglected so far. * diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 65316239e9..c95d4b7b7e 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -602,6 +602,9 @@ namespace VectorTools Vector cell_data_1(dof_1.get_fe().dofs_per_cell); Vector cell_data_2(dof_2.get_fe().dofs_per_cell); + // Reset output vector. + data_2 = static_cast(0); + // Store how many cells share each dof (unghosted). OutVector touch_count; touch_count.reinit(data_2); @@ -660,11 +663,12 @@ namespace VectorTools ::dealii::internal::ElementAccess::set(value, i, data_2); } - // Compress data_2 to set the proper values on all the processors. + // Compress data_2 to set the proper values on all the parallel processes. data_2.compress(VectorOperation::insert); } + template #include @@ -16,134 +41,11 @@ template void -run_1d() -{ - using Vector = dealii::Vector; - - // Create grid - dealii::Triangulation triangulation( - typename dealii::Triangulation::MeshSmoothing( - dealii::Triangulation::smoothing_on_refinement | - dealii::Triangulation::smoothing_on_coarsening)); - const unsigned int n_1d_cells = 2; - dealii::GridGenerator::subdivided_hyper_cube(triangulation, n_1d_cells); - - // Create low order system - dealii::DoFHandler dof_handler_coarse(triangulation); - - const unsigned int degree_coarse = 2; - const dealii::FE_Q fe_q_coarse(degree_coarse); - const dealii::FESystem fe_system_coarse(fe_q_coarse, dim); - dof_handler_coarse.initialize(triangulation, fe_system_coarse); - dof_handler_coarse.distribute_dofs(fe_system_coarse); - - Vector solution_coarse; - solution_coarse.reinit(dof_handler_coarse.n_dofs()); - // Initialize dummy coarse solution - // solution_coarse.add(1.0); - for (unsigned int idof = 0; idof < dof_handler_coarse.n_dofs(); ++idof) - { - solution_coarse[idof] = idof; - } - - const unsigned int degree_fine = 3; - - dealii::DoFHandler dof_handler_fine(triangulation); - - const dealii::FE_Q fe_q_fine(degree_fine); - const dealii::FESystem fe_system_fine(fe_q_fine, dim); - dof_handler_fine.initialize(triangulation, fe_system_fine); - dof_handler_fine.distribute_dofs(fe_system_fine); - - Vector solution_fine; - solution_fine.reinit(dof_handler_fine.n_dofs()); - solution_fine.add(1.0); // Initialize to non-zero - - // Interpolate the solution - dealii::FullMatrix interpolation_matrix( - fe_system_fine.dofs_per_cell, fe_system_coarse.dofs_per_cell); - fe_system_fine.get_interpolation_matrix(fe_system_coarse, - interpolation_matrix); - dealii::VectorTools::interpolate(dof_handler_coarse, - dof_handler_fine, - interpolation_matrix, - solution_coarse, - solution_fine); - - - const unsigned int n_dofs_coarse = fe_system_coarse.dofs_per_cell; - const unsigned int n_dofs_fine = fe_system_fine.dofs_per_cell; - const std::vector> &points = - fe_system_fine.get_unit_support_points(); - - // Check that interpolated solution matches at the support points - std::vector dof_indices_coarse( - n_dofs_coarse); - std::vector dof_indices_fine(n_dofs_fine); - - auto cell_coarse = dof_handler_coarse.begin_active(); - auto cell_fine = dof_handler_fine.begin_active(); - auto endcell = dof_handler_fine.end(); - - bool matching = true; - for (; cell_fine != dof_handler_fine.end(); ++cell_fine, ++cell_coarse) - { - if (!cell_fine->is_locally_owned()) - continue; - cell_fine->get_dof_indices(dof_indices_fine); - cell_coarse->get_dof_indices(dof_indices_coarse); - for (unsigned int iquad = 0; iquad < points.size(); ++iquad) - { - dealii::Tensor<1, dim, double> local_solution_coarse; - dealii::Tensor<1, dim, double> local_solution_fine; - for (unsigned int idof = 0; idof < n_dofs_coarse; ++idof) - { - const unsigned int axis = - fe_system_coarse.system_to_component_index(idof).first; - local_solution_coarse[axis] += - solution_coarse[dof_indices_coarse[idof]] * - fe_system_coarse.shape_value(idof, points[iquad]); - } - for (unsigned int idof = 0; idof < n_dofs_fine; ++idof) - { - const unsigned int axis = - fe_system_fine.system_to_component_index(idof).first; - local_solution_fine[axis] += - solution_fine[dof_indices_fine[idof]] * - fe_system_fine.shape_value(idof, points[iquad]); - } - dealii::Tensor<1, dim, double> diff = local_solution_fine; - diff -= local_solution_coarse; - const double diff_norm = diff.norm(); - - std::cout << "Fine solution " << local_solution_fine - << " Coarse solution " << local_solution_coarse - << " Difference " << diff_norm << std::endl; - if (diff_norm > 1e-12) - matching = false; - } - } - if (matching) - dealii::deallog << "OK" << std::endl; - if (!matching) - dealii::deallog << "Non matching" << std::endl; -} - -template -void -run() +run(const dealii::Triangulation &triangulation) { using Vector = dealii::LinearAlgebra::distributed::Vector; - // Create grid - dealii::parallel::distributed::Triangulation triangulation( - MPI_COMM_WORLD, - typename dealii::Triangulation::MeshSmoothing( - dealii::Triangulation::smoothing_on_refinement | - dealii::Triangulation::smoothing_on_coarsening)); - const unsigned int n_1d_cells = 2; - dealii::GridGenerator::subdivided_hyper_cube(triangulation, n_1d_cells); - // Create low order system + // Create low order system. dealii::DoFHandler dof_handler_coarse(triangulation); const unsigned int degree_coarse = 2; @@ -160,10 +62,19 @@ run() ghost_dofs_coarse); dealii::IndexSet locally_relevant_dofs_coarse = ghost_dofs_coarse; ghost_dofs_coarse.subtract_set(locally_owned_dofs_coarse); - solution_coarse.reinit(locally_owned_dofs_coarse, - ghost_dofs_coarse, - MPI_COMM_WORLD); - // Initialize dummy coarse solution + + if (dim == 1) + { + solution_coarse.reinit(locally_owned_dofs_coarse, MPI_COMM_WORLD); + } + else + { + solution_coarse.reinit(locally_owned_dofs_coarse, + ghost_dofs_coarse, + MPI_COMM_WORLD); + } + + // Initialize dummy coarse solution. for (unsigned int idof = 0; idof < dof_handler_coarse.n_dofs(); ++idof) { if (locally_owned_dofs_coarse.is_element(idof)) @@ -190,13 +101,19 @@ run() ghost_dofs_fine); dealii::IndexSet locally_relevant_dofs_fine = ghost_dofs_fine; ghost_dofs_fine.subtract_set(locally_owned_dofs_fine); - solution_fine.reinit(locally_owned_dofs_fine, - ghost_dofs_fine, - MPI_COMM_WORLD); - solution_fine.add(1.0); // Initialize to non-zero - solution_fine.update_ghost_values(); - // Interpolate the solution + if (dim == 1) + { + solution_fine.reinit(locally_owned_dofs_fine, MPI_COMM_WORLD); + } + else + { + solution_fine.reinit(locally_owned_dofs_fine, + ghost_dofs_fine, + MPI_COMM_WORLD); + } + + // Interpolate the solution. dealii::FullMatrix interpolation_matrix( fe_system_fine.dofs_per_cell, fe_system_coarse.dofs_per_cell); fe_system_fine.get_interpolation_matrix(fe_system_coarse, @@ -206,14 +123,14 @@ run() interpolation_matrix, solution_coarse, solution_fine); - + solution_fine.update_ghost_values(); const unsigned int n_dofs_coarse = fe_system_coarse.dofs_per_cell; const unsigned int n_dofs_fine = fe_system_fine.dofs_per_cell; const std::vector> &points = fe_system_fine.get_unit_support_points(); - // Check that interpolated solution matches at the support points + // Check that interpolated solution matches at quadrature points. std::vector dof_indices_coarse( n_dofs_coarse); std::vector dof_indices_fine(n_dofs_fine); @@ -223,16 +140,19 @@ run() auto endcell = dof_handler_fine.end(); bool matching = true; - for (; cell_fine != dof_handler_fine.end(); ++cell_fine, ++cell_coarse) + for (; cell_fine != endcell; ++cell_fine, ++cell_coarse) { if (!cell_fine->is_locally_owned()) continue; + cell_fine->get_dof_indices(dof_indices_fine); cell_coarse->get_dof_indices(dof_indices_coarse); + for (unsigned int iquad = 0; iquad < points.size(); ++iquad) { dealii::Tensor<1, dim, double> local_solution_coarse; dealii::Tensor<1, dim, double> local_solution_fine; + for (unsigned int idof = 0; idof < n_dofs_coarse; ++idof) { const unsigned int axis = @@ -241,6 +161,7 @@ run() solution_coarse[dof_indices_coarse[idof]] * fe_system_coarse.shape_value(idof, points[iquad]); } + for (unsigned int idof = 0; idof < n_dofs_fine; ++idof) { const unsigned int axis = @@ -249,21 +170,31 @@ run() solution_fine[dof_indices_fine[idof]] * fe_system_fine.shape_value(idof, points[iquad]); } + dealii::Tensor<1, dim, double> diff = local_solution_fine; diff -= local_solution_coarse; + const double diff_norm = diff.norm(); - std::cout << "Fine solution " << local_solution_fine - << " Coarse solution " << local_solution_coarse - << " Difference " << diff_norm << std::endl; + std::cout << "Coarse solution = " << local_solution_coarse + << ", fine solution = " << local_solution_fine + << ", difference = " << diff_norm << std::endl; + if (diff_norm > 1e-12) - matching = false; + { + matching = false; + } } } + if (matching) - dealii::deallog << "OK" << std::endl; - if (!matching) - dealii::deallog << "Non matching" << std::endl; + { + dealii::deallog << "OK" << std::endl; + } + else + { + dealii::deallog << "Non matching" << std::endl; + } } int @@ -272,11 +203,34 @@ main(int argc, char *argv[]) dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); initlog(); - // dealii::deallog.push("1d"); - dealii::deallog << "1d" << std::endl; - run_1d<1>(); - dealii::deallog << "2d" << std::endl; - run<2>(); - dealii::deallog << "3d" << std::endl; - run<3>(); + const unsigned int n_1d_cells = 2; + + { + // 1D + dealii::Triangulation<1> triangulation; + dealii::GridGenerator::subdivided_hyper_cube(triangulation, n_1d_cells); + + dealii::deallog << "1D" << std::endl; + run<1>(triangulation); + } + + // 2D. + { + dealii::parallel::distributed::Triangulation<2> triangulation( + MPI_COMM_WORLD); + dealii::GridGenerator::subdivided_hyper_cube(triangulation, n_1d_cells); + + dealii::deallog << "2D" << std::endl; + run<2>(triangulation); + } + + // 3D. + { + dealii::parallel::distributed::Triangulation<3> triangulation( + MPI_COMM_WORLD); + dealii::GridGenerator::subdivided_hyper_cube(triangulation, n_1d_cells); + + dealii::deallog << "3D" << std::endl; + run<3>(triangulation); + } } diff --git a/tests/vector_tools/interpolate.with_p4est=true.mpirun=1.output b/tests/vector_tools/interpolate.with_p4est=true.mpirun=1.output index 4b883267c3..661e4988ab 100644 --- a/tests/vector_tools/interpolate.with_p4est=true.mpirun=1.output +++ b/tests/vector_tools/interpolate.with_p4est=true.mpirun=1.output @@ -1,7 +1,7 @@ -DEAL::1d +DEAL::1D DEAL::OK -DEAL::2d +DEAL::2D DEAL::OK -DEAL::3d +DEAL::3D DEAL::OK diff --git a/tests/vector_tools/interpolate.with_p4est=true.mpirun=2.output b/tests/vector_tools/interpolate.with_p4est=true.mpirun=2.output index 4b883267c3..661e4988ab 100644 --- a/tests/vector_tools/interpolate.with_p4est=true.mpirun=2.output +++ b/tests/vector_tools/interpolate.with_p4est=true.mpirun=2.output @@ -1,7 +1,7 @@ -DEAL::1d +DEAL::1D DEAL::OK -DEAL::2d +DEAL::2D DEAL::OK -DEAL::3d +DEAL::3D DEAL::OK