From: Denis Davydov Date: Fri, 30 Mar 2018 11:37:18 +0000 (+0200) Subject: make VectorTools::get_position_vector() write into vector only on locally owned cells X-Git-Tag: v9.0.0-rc1~251^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F6127%2Fhead;p=dealii.git make VectorTools::get_position_vector() write into vector only on locally owned cells --- diff --git a/doc/news/changes/minor/20180401DenisDavydov b/doc/news/changes/minor/20180401DenisDavydov new file mode 100644 index 0000000000..71838b5f41 --- /dev/null +++ b/doc/news/changes/minor/20180401DenisDavydov @@ -0,0 +1,3 @@ +Improved: VectorTools::get_position_vector() now supports parallel triangulations. +
+(Denis Davydov, 2018/04/01) diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index ca54be3b5a..bcc6399bae 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -8119,7 +8119,6 @@ namespace VectorTools if ( fe.has_support_points() ) { - typename DoFHandlerType::active_cell_iterator cell; const Quadrature quad(fe.get_unit_support_points()); MappingQ map_q(fe.degree); @@ -8129,19 +8128,20 @@ namespace VectorTools AssertDimension(fe.dofs_per_cell, fe.get_unit_support_points().size()); Assert(fe.is_primitive(), ExcMessage("FE is not Primitive! This won't work.")); - for (cell = dh.begin_active(); cell != dh.end(); ++cell) - { - fe_v.reinit(cell); - cell->get_dof_indices(dofs); - const std::vector > &points = fe_v.get_quadrature_points(); - for (unsigned int q = 0; q < points.size(); ++q) - { - unsigned int comp = fe.system_to_component_index(q).first; - if (fe_mask[comp]) - ::dealii::internal::ElementAccess::set( - points[q][fe_to_real[comp]], dofs[q], vector); - } - } + for (const auto &cell : dh.active_cell_iterators()) + if (cell->is_locally_owned()) + { + fe_v.reinit(cell); + cell->get_dof_indices(dofs); + const std::vector > &points = fe_v.get_quadrature_points(); + for (unsigned int q = 0; q < points.size(); ++q) + { + const unsigned int comp = fe.system_to_component_index(q).first; + if (fe_mask[comp]) + ::dealii::internal::ElementAccess::set( + points[q][fe_to_real[comp]], dofs[q], vector); + } + } } else { diff --git a/tests/mappings/mapping_fe_field_03.cc b/tests/mappings/mapping_fe_field_03.cc new file mode 100644 index 0000000000..821ef1bc40 --- /dev/null +++ b/tests/mappings/mapping_fe_field_03.cc @@ -0,0 +1,86 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2001 - 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. +// +// --------------------------------------------------------------------- + +// tests VectorTools::get_position_vector() used in MappingFEField +// with parallel distributed triangulation + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include + +template +void test() +{ + parallel::distributed::Triangulation tria(MPI_COMM_WORLD, + Triangulation::limit_level_difference_at_vertices, + parallel::distributed::Triangulation::construct_multigrid_hierarchy); + + { + Point center; + GridGenerator::hyper_ball (tria, + center, + 12); + + const types::manifold_id sphere_id = 0; + static const SphericalManifold boundary_ball(center); + tria.set_all_manifold_ids_on_boundary(sphere_id); + tria.set_manifold (sphere_id, boundary_ball); + + tria.refine_global (dim==2?1:3); + } + + FESystem fe(FE_Q(4), dim); + DoFHandler dh(tria); + + dh.distribute_dofs(fe); + + deallog << "dim: " << dim << std::endl + << "cells: " << tria.n_global_active_cells() << ", dofs: " + << dh.n_dofs() < map_vector; + + IndexSet locally_relevant_euler; + DoFTools::extract_locally_relevant_dofs (dh, locally_relevant_euler); + map_vector.reinit(dh.locally_owned_dofs(), locally_relevant_euler, MPI_COMM_WORLD); + + VectorTools::get_position_vector(dh, map_vector); + + deallog << "L2=" << map_vector.l2_norm() << std::endl + << "L1=" << map_vector.l1_norm() << std::endl + << "Linfty=" << map_vector.linfty_norm() << std::endl; + + +} + +using namespace dealii; +int main (int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0) + initlog(); + + test<2>(); + test<3>(); +} diff --git a/tests/mappings/mapping_fe_field_03.mpirun=1.output b/tests/mappings/mapping_fe_field_03.mpirun=1.output new file mode 100644 index 0000000000..1bedcc6eb0 --- /dev/null +++ b/tests/mappings/mapping_fe_field_03.mpirun=1.output @@ -0,0 +1,11 @@ + +DEAL::dim: 2 +DEAL::cells: 20, dofs: 674 +DEAL::L2=142.998 +DEAL::L1=3043.68 +DEAL::Linfty=12.0000 +DEAL::dim: 3 +DEAL::cells: 3584, dofs: 697827 +DEAL::L2=3538.82 +DEAL::L1=2.38344e+06 +DEAL::Linfty=12.0000 diff --git a/tests/mappings/mapping_fe_field_03.mpirun=3.output b/tests/mappings/mapping_fe_field_03.mpirun=3.output new file mode 100644 index 0000000000..1bedcc6eb0 --- /dev/null +++ b/tests/mappings/mapping_fe_field_03.mpirun=3.output @@ -0,0 +1,11 @@ + +DEAL::dim: 2 +DEAL::cells: 20, dofs: 674 +DEAL::L2=142.998 +DEAL::L1=3043.68 +DEAL::Linfty=12.0000 +DEAL::dim: 3 +DEAL::cells: 3584, dofs: 697827 +DEAL::L2=3538.82 +DEAL::L1=2.38344e+06 +DEAL::Linfty=12.0000 diff --git a/tests/mpi/step-37.cc b/tests/mpi/step-37.cc new file mode 100644 index 0000000000..1d53080716 --- /dev/null +++ b/tests/mpi/step-37.cc @@ -0,0 +1,568 @@ +/* --------------------------------------------------------------------- + * + * 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. + * + * --------------------------------------------------------------------- + + * + * a light modification of step-37 tutorial to test MappingFEField + * in parallel with non-zero BC and zero body load. + */ + + +#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 + + +namespace Step37 +{ + using namespace dealii; + + + const unsigned int degree_finite_element = 2; + const unsigned int dimension = 3; + + + template + class LaplaceProblem + { + public: + LaplaceProblem (); + void run (); + + ~LaplaceProblem(); + + private: + void setup_system (); + void assemble_rhs (); + void solve (); + void output_results (const unsigned int cycle) const; + +#ifdef DEAL_II_WITH_P4EST + parallel::distributed::Triangulation triangulation; +#else + Triangulation triangulation; +#endif + + FE_Q fe; + DoFHandler dof_handler; + + FE_Q fe_euler; + FESystem fe_system; + DoFHandler dof_euler; + std::shared_ptr>> mapping; + ConstraintMatrix constraints_euler; + LinearAlgebra::distributed::Vector euler_positions; + + IndexSet locally_relevant_dofs; + + ConstraintMatrix constraints; + ConstraintMatrix non_homogeneous_constraints; + typedef MatrixFreeOperators::LaplaceOperator> SystemMatrixType; + SystemMatrixType system_matrix; + + MGConstrainedDoFs mg_constrained_dofs; + typedef MatrixFreeOperators::LaplaceOperator> LevelMatrixType; + MGLevelObject mg_matrices; + + LinearAlgebra::distributed::Vector solution; + LinearAlgebra::distributed::Vector system_rhs; + + ConditionalOStream pcout; + }; + + + + template + LaplaceProblem::LaplaceProblem () + : +#ifdef DEAL_II_WITH_P4EST + triangulation(MPI_COMM_WORLD, + Triangulation::limit_level_difference_at_vertices, + parallel::distributed::Triangulation::construct_multigrid_hierarchy), +#else + triangulation (Triangulation::limit_level_difference_at_vertices), +#endif + fe (degree_finite_element), + dof_handler (triangulation), + fe_euler(degree_finite_element), + fe_system(fe_euler,dim), + dof_euler(triangulation), + pcout (std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + {} + + + + template + LaplaceProblem::~LaplaceProblem() + { + mapping.reset(); + dof_euler.clear(); + } + + + + + + template + void LaplaceProblem::setup_system () + { + system_matrix.clear(); + mg_matrices.clear_elements(); + + dof_handler.distribute_dofs (fe); + dof_handler.distribute_mg_dofs (); + + dof_euler.distribute_dofs(fe_system); + { + IndexSet locally_relevant_euler; + DoFTools::extract_locally_relevant_dofs (dof_euler, + locally_relevant_euler); + euler_positions.reinit(dof_euler.locally_owned_dofs(), locally_relevant_euler, MPI_COMM_WORLD); + } + + // Set up a quadrature formula based on the support points of FE_Q and go + // through the mesh with a MappingQ of your favorite degree; on each cell, + // evaluate the quadrature points and write them into the vector to be + // associated with MappingFEField + VectorTools::get_position_vector(dof_euler, euler_positions); + + // Apply hanging node constraints on that vector + constraints_euler.clear(); + DoFTools::make_hanging_node_constraints(dof_euler, + constraints_euler); + constraints_euler.close(); + constraints_euler.distribute(euler_positions); + + euler_positions.update_ghost_values(); + + mapping = std::make_shared>>(dof_euler, euler_positions); + + pcout << "Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + DoFTools::extract_locally_relevant_dofs (dof_handler, + locally_relevant_dofs); + + constraints.clear(); + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, + constraints); + VectorTools::interpolate_boundary_values (dof_handler, + 0, + Functions::ZeroFunction(), + constraints); + constraints.close(); + + { + typename MatrixFree::AdditionalData additional_data; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + additional_data.mapping_update_flags = (update_gradients | update_JxW_values | + update_quadrature_points); + std::shared_ptr > + system_mf_storage(new MatrixFree()); + system_mf_storage->reinit (*mapping.get(), dof_handler, constraints, QGauss<1>(fe.degree+1), + additional_data); + system_matrix.initialize (system_mf_storage); + } + + system_matrix.compute_diagonal(); + + system_matrix.initialize_dof_vector(solution); + system_matrix.initialize_dof_vector(system_rhs); + + const unsigned int nlevels = triangulation.n_global_levels(); + mg_matrices.resize(0, nlevels-1); + + std::set dirichlet_boundary; + dirichlet_boundary.insert(0); + mg_constrained_dofs.initialize(dof_handler); + mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, dirichlet_boundary); + + for (unsigned int level=0; level::AdditionalData additional_data; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + additional_data.mapping_update_flags = (update_gradients | update_JxW_values | + update_quadrature_points); + additional_data.level_mg_handler = level; + std::shared_ptr > + mg_mf_storage_level(new MatrixFree()); + mg_mf_storage_level->reinit(dof_handler, level_constraints, + QGauss<1>(fe.degree+1), additional_data); + + mg_matrices[level].initialize(mg_mf_storage_level, mg_constrained_dofs, + level); + mg_matrices[level].compute_diagonal(); + } + } + + + + template + class PotentialBCFunction : public Function + { + public: + PotentialBCFunction(const double &charge, const Point &dipole) + : + Function(1), + charge(charge), + x0(std::abs(charge) < 1e-10? Point() : dipole/charge) + { + } + + virtual ~PotentialBCFunction() = default; + + virtual double value(const Point &p, const unsigned int) const + { + const double r = p.distance(x0); + Assert (r > 0, ExcDivideByZero()); + return charge/r; + } + + private: + const double charge; + const Point x0; + }; + + + + + template + void LaplaceProblem::assemble_rhs () + { + Timer time; + + non_homogeneous_constraints.clear(); + { + ConstraintMatrix hanging_nodes_laplace_constraints; + hanging_nodes_laplace_constraints.reinit (locally_relevant_dofs); + non_homogeneous_constraints.reinit (locally_relevant_dofs); + DoFTools::make_hanging_node_constraints (dof_handler, hanging_nodes_laplace_constraints); + + std::set dirichlet_boundary_ids; + typename FunctionMap::type dirichlet_boundary_functions; + PotentialBCFunction bc_func (240, Point()); + dirichlet_boundary_ids.insert(0); + dirichlet_boundary_functions[0] = &bc_func; + VectorTools::interpolate_boundary_values (*mapping.get(), + dof_handler, + dirichlet_boundary_functions, + non_homogeneous_constraints); + // make sure hanging nodes override Dirichlet + non_homogeneous_constraints.merge(hanging_nodes_laplace_constraints, ConstraintMatrix::MergeConflictBehavior::right_object_wins); + non_homogeneous_constraints.close (); + } + + solution = 0; + non_homogeneous_constraints.distribute(solution); + system_rhs = 0; + solution.update_ghost_values(); + FEEvaluation phi(*system_matrix.get_matrix_free()); + for (unsigned int cell=0; celln_macro_cells(); ++cell) + { + phi.reinit(cell); + phi.read_dof_values_plain(solution); + phi.evaluate(false, true); + for (unsigned int q=0; q(1.0), q); + } + //phi.integrate(true, true); + phi.integrate(false, true); + phi.distribute_local_to_global(system_rhs); + } + system_rhs.compress(VectorOperation::add); + } + + + + + template + void LaplaceProblem::solve () + { + MGTransferMatrixFree mg_transfer(mg_constrained_dofs); + mg_transfer.build(dof_handler); + + typedef PreconditionChebyshev > SmootherType; + mg::SmootherRelaxation > + mg_smoother; + MGLevelObject smoother_data; + smoother_data.resize(0, triangulation.n_global_levels()-1); + for (unsigned int level = 0; level 0) + { + smoother_data[level].smoothing_range = 15.; + smoother_data[level].degree = 4; + smoother_data[level].eig_cg_n_iterations = 10; + } + else + { + smoother_data[0].smoothing_range = 1e-3; + smoother_data[0].degree = numbers::invalid_unsigned_int; + smoother_data[0].eig_cg_n_iterations = mg_matrices[0].m(); + } + mg_matrices[level].compute_diagonal(); + smoother_data[level].preconditioner = mg_matrices[level].get_matrix_diagonal_inverse(); + } + mg_smoother.initialize(mg_matrices, smoother_data); + + MGCoarseGridApplySmoother > mg_coarse; + mg_coarse.initialize(mg_smoother); + + mg::Matrix > mg_matrix(mg_matrices); + + MGLevelObject > mg_interface_matrices; + mg_interface_matrices.resize(0, triangulation.n_global_levels()-1); + for (unsigned int level=0; level > mg_interface(mg_interface_matrices); + + Multigrid > mg(mg_matrix, + mg_coarse, + mg_transfer, + mg_smoother, + mg_smoother); + mg.set_edge_matrices(mg_interface, mg_interface); + + PreconditionMG, + MGTransferMatrixFree > + preconditioner(dof_handler, mg, mg_transfer); + + + SolverControl solver_control (500, 1e-12*system_rhs.l2_norm()); + SolverCG > cg (solver_control); + + non_homogeneous_constraints.set_zero(solution); + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + + non_homogeneous_constraints.distribute(solution); + + const double linfty = Utilities::MPI::max(solution.linfty_norm(), MPI_COMM_WORLD); + + pcout << "Solved in " + << solver_control.last_step() + << " iterations" << std::endl + << "Linfty=" << linfty << std::endl; + } + + + + + template + void LaplaceProblem::output_results (const unsigned int cycle) const + { + if (triangulation.n_global_active_cells() > 1000000) + return; + + DataOut data_out; + + solution.update_ghost_values(); + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "solution"); + data_out.build_patches (); + + std::ofstream output ("solution-" + + std::to_string(cycle) + + "." + + std::to_string(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)) + + ".vtu"); + data_out.write_vtu (output); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + { + std::vector filenames; + for (unsigned int i=0; i > support_points; + MappingQ mapping(degree_finite_element); + DoFTools::map_dofs_to_support_points(mapping, dof_handler, support_points); + + const std::string base_filename = "grid" + dealii::Utilities::int_to_string(dim) + "_" + dealii::Utilities::int_to_string(cycle); + const std::string filename = base_filename + ".gp"; + std::ofstream f(filename.c_str()); + + f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" << std::endl + << "set output \"" << base_filename << ".png\"" << std::endl + << "set size square" << std::endl + << "set view equal xy" << std::endl + << "unset xtics" << std::endl + << "unset ytics" << std::endl + << "plot '-' using 1:2 with lines notitle, '-' with labels point pt 2 offset 1,1 notitle" << std::endl; + GridOut().write_gnuplot (triangulation, f); + f << "e" << std::endl; + + DoFTools::write_gnuplot_dof_support_point_info(f, + support_points); + + f << "e" << std::endl; + } + } + + + + + template + void LaplaceProblem::run () + { + for (unsigned int cycle=0; cycle<2; ++cycle) + { + pcout << "Cycle " << cycle << std::endl; + + if (cycle == 0) + { + Point center; + GridGenerator::hyper_ball (triangulation, + center, + 12); + + const types::manifold_id sphere_id = 0; + static const SphericalManifold boundary_ball(center); + triangulation.set_all_manifold_ids_on_boundary(sphere_id); + triangulation.set_manifold (sphere_id, boundary_ball); + + triangulation.refine_global (dim==2?1:3); + } + else + { + for (auto cell: triangulation.active_cell_iterators()) + if (cell->is_locally_owned()) + { + if ((cell->center()[0] < -8. && cell->center()[1] < 0) || + (cell->center()[0] > 9. && cell->center()[1] > 5)) + cell->set_refine_flag (); + } + triangulation.execute_coarsening_and_refinement (); + // triangulation.refine_global(); + } + setup_system (); + assemble_rhs (); + + if (cycle==0 && dim == 2) + { + pcout << "Constraints:" << std::endl; + non_homogeneous_constraints.print(std::cout); + } + + solve (); + output_results (cycle); + pcout << std::endl; + }; + } +} + + + + +int main (int argc, char *argv[]) +{ + try + { + using namespace Step37; + + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + LaplaceProblem laplace_problem; + laplace_problem.run (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/mpi/step-37.with_p4est=true.mpirun=1.output b/tests/mpi/step-37.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..e1e8554fa7 --- /dev/null +++ b/tests/mpi/step-37.with_p4est=true.mpirun=1.output @@ -0,0 +1,10 @@ +Cycle 0 +Number of degrees of freedom: 29521 +Solved in 8 iterations +Linfty=20 + +Cycle 1 +Number of degrees of freedom: 33518 +Solved in 12 iterations +Linfty=20 + diff --git a/tests/mpi/step-37.with_p4est=true.mpirun=7.output b/tests/mpi/step-37.with_p4est=true.mpirun=7.output new file mode 100644 index 0000000000..e1e8554fa7 --- /dev/null +++ b/tests/mpi/step-37.with_p4est=true.mpirun=7.output @@ -0,0 +1,10 @@ +Cycle 0 +Number of degrees of freedom: 29521 +Solved in 8 iterations +Linfty=20 + +Cycle 1 +Number of degrees of freedom: 33518 +Solved in 12 iterations +Linfty=20 +