]> https://gitweb.dealii.org/ - dealii.git/commitdiff
make VectorTools::get_position_vector() write into vector only on locally owned cells 6127/head
authorDenis Davydov <davydden@gmail.com>
Fri, 30 Mar 2018 11:37:18 +0000 (13:37 +0200)
committerDenis Davydov <davydden@gmail.com>
Sun, 1 Apr 2018 10:18:20 +0000 (12:18 +0200)
doc/news/changes/minor/20180401DenisDavydov [new file with mode: 0644]
include/deal.II/numerics/vector_tools.templates.h
tests/mappings/mapping_fe_field_03.cc [new file with mode: 0644]
tests/mappings/mapping_fe_field_03.mpirun=1.output [new file with mode: 0644]
tests/mappings/mapping_fe_field_03.mpirun=3.output [new file with mode: 0644]
tests/mpi/step-37.cc [new file with mode: 0644]
tests/mpi/step-37.with_p4est=true.mpirun=1.output [new file with mode: 0644]
tests/mpi/step-37.with_p4est=true.mpirun=7.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20180401DenisDavydov b/doc/news/changes/minor/20180401DenisDavydov
new file mode 100644 (file)
index 0000000..71838b5
--- /dev/null
@@ -0,0 +1,3 @@
+Improved: VectorTools::get_position_vector() now supports parallel triangulations.
+<br>
+(Denis Davydov, 2018/04/01)
index ca54be3b5af6cc15f23f154fbd76240670ff2893..bcc6399baef717f696d44a6b6e0a9e8c5ef8ae1b 100644 (file)
@@ -8119,7 +8119,6 @@ namespace VectorTools
 
     if ( fe.has_support_points() )
       {
-        typename DoFHandlerType<dim,spacedim>::active_cell_iterator cell;
         const Quadrature<dim> quad(fe.get_unit_support_points());
 
         MappingQ<dim,spacedim> 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<Point<spacedim> > &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<VectorType>::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<Point<spacedim> > &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<VectorType>::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 (file)
index 0000000..821ef1b
--- /dev/null
@@ -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 <deal.II/distributed/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_fe_field.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/dofs/dof_tools.h>
+
+template<int dim>
+void test()
+{
+  parallel::distributed::Triangulation<dim>  tria(MPI_COMM_WORLD,
+                                                  Triangulation<dim>::limit_level_difference_at_vertices,
+                                                  parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
+
+  {
+    Point<dim> center;
+    GridGenerator::hyper_ball (tria,
+                               center,
+                               12);
+
+    const types::manifold_id sphere_id = 0;
+    static const SphericalManifold<dim> 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<dim> fe(FE_Q<dim>(4), dim);
+  DoFHandler<dim> dh(tria);
+
+  dh.distribute_dofs(fe);
+
+  deallog << "dim: " << dim << std::endl
+          << "cells: " << tria.n_global_active_cells() << ", dofs: "
+          << dh.n_dofs() <<std::endl;
+
+  // Create a Mapping
+  LinearAlgebra::distributed::Vector<double> 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 (file)
index 0000000..1bedcc6
--- /dev/null
@@ -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 (file)
index 0000000..1bedcc6
--- /dev/null
@@ -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 (file)
index 0000000..1d53080
--- /dev/null
@@ -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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/timer.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/precondition.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/multigrid/multigrid.h>
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+#include <deal.II/multigrid/mg_tools.h>
+#include <deal.II/multigrid/mg_coarse.h>
+#include <deal.II/multigrid/mg_smoother.h>
+#include <deal.II/multigrid/mg_matrix.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/operators.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/operators.h>
+
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_fe_field.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/grid/grid_out.h>
+
+#include <iostream>
+#include <fstream>
+
+
+namespace Step37
+{
+  using namespace dealii;
+
+
+  const unsigned int degree_finite_element = 2;
+  const unsigned int dimension = 3;
+
+
+  template <int dim>
+  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<dim>  triangulation;
+#else
+    Triangulation<dim>                         triangulation;
+#endif
+
+    FE_Q<dim>                                  fe;
+    DoFHandler<dim>                            dof_handler;
+
+    FE_Q<dim>                                  fe_euler;
+    FESystem<dim>                              fe_system;
+    DoFHandler<dim>                            dof_euler;
+    std::shared_ptr<MappingFEField<dim,dim,LinearAlgebra::distributed::Vector<double>>>       mapping;
+    ConstraintMatrix                           constraints_euler;
+    LinearAlgebra::distributed::Vector<double> euler_positions;
+
+    IndexSet locally_relevant_dofs;
+
+    ConstraintMatrix                           constraints;
+    ConstraintMatrix                           non_homogeneous_constraints;
+    typedef MatrixFreeOperators::LaplaceOperator<dim,degree_finite_element,degree_finite_element+1,1,LinearAlgebra::distributed::Vector<double>> SystemMatrixType;
+    SystemMatrixType                           system_matrix;
+
+    MGConstrainedDoFs                          mg_constrained_dofs;
+    typedef MatrixFreeOperators::LaplaceOperator<dim,degree_finite_element,degree_finite_element+1,1,LinearAlgebra::distributed::Vector<float>>  LevelMatrixType;
+    MGLevelObject<LevelMatrixType>             mg_matrices;
+
+    LinearAlgebra::distributed::Vector<double> solution;
+    LinearAlgebra::distributed::Vector<double> system_rhs;
+
+    ConditionalOStream                         pcout;
+  };
+
+
+
+  template <int dim>
+  LaplaceProblem<dim>::LaplaceProblem ()
+    :
+#ifdef DEAL_II_WITH_P4EST
+    triangulation(MPI_COMM_WORLD,
+                  Triangulation<dim>::limit_level_difference_at_vertices,
+                  parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy),
+#else
+    triangulation (Triangulation<dim>::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 <int dim>
+  LaplaceProblem<dim>::~LaplaceProblem()
+  {
+    mapping.reset();
+    dof_euler.clear();
+  }
+
+
+
+
+
+  template <int dim>
+  void LaplaceProblem<dim>::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<MappingFEField<dim,dim,LinearAlgebra::distributed::Vector<double>>>(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<dim>(),
+                                              constraints);
+    constraints.close();
+
+    {
+      typename MatrixFree<dim,double>::AdditionalData additional_data;
+      additional_data.tasks_parallel_scheme =
+        MatrixFree<dim,double>::AdditionalData::none;
+      additional_data.mapping_update_flags = (update_gradients | update_JxW_values |
+                                              update_quadrature_points);
+      std::shared_ptr<MatrixFree<dim,double> >
+      system_mf_storage(new MatrixFree<dim,double>());
+      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<types::boundary_id> 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<nlevels; ++level)
+      {
+        IndexSet relevant_dofs;
+        DoFTools::extract_locally_relevant_level_dofs(dof_handler, level,
+                                                      relevant_dofs);
+        ConstraintMatrix level_constraints;
+        level_constraints.reinit(relevant_dofs);
+        level_constraints.add_lines(mg_constrained_dofs.get_boundary_indices(level));
+        level_constraints.close();
+
+        typename MatrixFree<dim,float>::AdditionalData additional_data;
+        additional_data.tasks_parallel_scheme =
+          MatrixFree<dim,float>::AdditionalData::none;
+        additional_data.mapping_update_flags = (update_gradients | update_JxW_values |
+                                                update_quadrature_points);
+        additional_data.level_mg_handler = level;
+        std::shared_ptr<MatrixFree<dim,float> >
+        mg_mf_storage_level(new MatrixFree<dim,float>());
+        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 <int dim>
+  class PotentialBCFunction : public Function<dim>
+  {
+  public:
+    PotentialBCFunction(const double &charge, const Point<dim> &dipole)
+      :
+      Function<dim>(1),
+      charge(charge),
+      x0(std::abs(charge) < 1e-10? Point<dim>() : dipole/charge)
+    {
+    }
+
+    virtual ~PotentialBCFunction() = default;
+
+    virtual double value(const Point<dim> &p, const unsigned int) const
+    {
+      const double r = p.distance(x0);
+      Assert (r > 0, ExcDivideByZero());
+      return charge/r;
+    }
+
+  private:
+    const double charge;
+    const Point<dim> x0;
+  };
+
+
+
+
+  template <int dim>
+  void LaplaceProblem<dim>::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<types::boundary_id>       dirichlet_boundary_ids;
+      typename FunctionMap<dim>::type    dirichlet_boundary_functions;
+      PotentialBCFunction<dim>           bc_func (240, Point<dim>());
+      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<dim,degree_finite_element> phi(*system_matrix.get_matrix_free());
+    for (unsigned int cell=0; cell<system_matrix.get_matrix_free()->n_macro_cells(); ++cell)
+      {
+        phi.reinit(cell);
+        phi.read_dof_values_plain(solution);
+        phi.evaluate(false, true);
+        for (unsigned int q=0; q<phi.n_q_points; ++q)
+          {
+            phi.submit_gradient(-phi.get_gradient(q), q);
+            // phi.submit_value(make_vectorized_array<double>(1.0), q);
+          }
+        //phi.integrate(true, true);
+        phi.integrate(false, true);
+        phi.distribute_local_to_global(system_rhs);
+      }
+    system_rhs.compress(VectorOperation::add);
+  }
+
+
+
+
+  template <int dim>
+  void LaplaceProblem<dim>::solve ()
+  {
+    MGTransferMatrixFree<dim,float> mg_transfer(mg_constrained_dofs);
+    mg_transfer.build(dof_handler);
+
+    typedef PreconditionChebyshev<LevelMatrixType,LinearAlgebra::distributed::Vector<float> > SmootherType;
+    mg::SmootherRelaxation<SmootherType, LinearAlgebra::distributed::Vector<float> >
+    mg_smoother;
+    MGLevelObject<typename SmootherType::AdditionalData> smoother_data;
+    smoother_data.resize(0, triangulation.n_global_levels()-1);
+    for (unsigned int level = 0; level<triangulation.n_global_levels(); ++level)
+      {
+        if (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<LinearAlgebra::distributed::Vector<float> > mg_coarse;
+    mg_coarse.initialize(mg_smoother);
+
+    mg::Matrix<LinearAlgebra::distributed::Vector<float> > mg_matrix(mg_matrices);
+
+    MGLevelObject<MatrixFreeOperators::MGInterfaceOperator<LevelMatrixType> > mg_interface_matrices;
+    mg_interface_matrices.resize(0, triangulation.n_global_levels()-1);
+    for (unsigned int level=0; level<triangulation.n_global_levels(); ++level)
+      mg_interface_matrices[level].initialize(mg_matrices[level]);
+    mg::Matrix<LinearAlgebra::distributed::Vector<float> > mg_interface(mg_interface_matrices);
+
+    Multigrid<LinearAlgebra::distributed::Vector<float> > mg(mg_matrix,
+                                                             mg_coarse,
+                                                             mg_transfer,
+                                                             mg_smoother,
+                                                             mg_smoother);
+    mg.set_edge_matrices(mg_interface, mg_interface);
+
+    PreconditionMG<dim, LinearAlgebra::distributed::Vector<float>,
+                   MGTransferMatrixFree<dim,float> >
+                   preconditioner(dof_handler, mg, mg_transfer);
+
+
+    SolverControl solver_control (500, 1e-12*system_rhs.l2_norm());
+    SolverCG<LinearAlgebra::distributed::Vector<double> > 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 <int dim>
+  void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
+  {
+    if (triangulation.n_global_active_cells() > 1000000)
+      return;
+
+    DataOut<dim> 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<std::string> filenames;
+        for (unsigned int i=0; i<Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); ++i)
+          filenames.emplace_back("solution-"
+                                 + std::to_string(cycle)
+                                 + "."
+                                 + std::to_string(i)
+                                 + ".vtu");
+
+        std::string master_name = "solution-" + Utilities::to_string(cycle) + ".pvtu";
+        std::ofstream master_output (master_name.c_str());
+        data_out.write_pvtu_record (master_output, filenames);
+      }
+
+    if (dim == 2 && Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1)
+      {
+        std::map<types::global_dof_index, Point<dim> > support_points;
+        MappingQ<dim> 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 <int dim>
+  void LaplaceProblem<dim>::run ()
+  {
+    for (unsigned int cycle=0; cycle<2; ++cycle)
+      {
+        pcout << "Cycle " << cycle << std::endl;
+
+        if (cycle == 0)
+          {
+            Point<dim> center;
+            GridGenerator::hyper_ball (triangulation,
+                                       center,
+                                       12);
+
+            const types::manifold_id sphere_id = 0;
+            static const SphericalManifold<dim> 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<dimension> 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 (file)
index 0000000..e1e8554
--- /dev/null
@@ -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 (file)
index 0000000..e1e8554
--- /dev/null
@@ -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
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.