--- /dev/null
+New: DoFTools::extract_dofs_with_support_on() returns a set of degrees of
+freedom that have support on a subset of locally owned cells.
+<br>
+(Denis Davydov, 2017/03/24)
std::vector<bool> &selected_dofs,
const std::set<types::boundary_id> &boundary_ids = std::set<types::boundary_id>());
+ /**
+ * Extract all degrees of freedom (DoFs) that have support only within cells,
+ * for which @p predicate is <code>true</code>. The result is returned as
+ * an IndexSet.
+ *
+ * Consider the following FE space where predicate returns <code>true</code>
+ * for all cells on the left side:
+ *
+ * @image html extract_dofs_with_support_on.png
+ *
+ * This functions will return the union of all DoFs on those cells minus
+ * DoF 11, 13, 2 and 0; the result will be <code>[9,10], 12, [14,38]</code>.
+ * On the image above the returned DoFs are separated from the rest by the
+ * red line
+ *
+ * Essentially, the question this functions answers is the following:
+ * Given a subdomain with associated DoFs, what is the largest subset of
+ * these DoFs that are allowed to be non-zero such that after calling
+ * ConstraintMatrix::distribute() the resulting solution vector will have
+ * support only within the given domain. Here @p cm is the ConstraintMatrix
+ * containing hanging nodes constraints.
+ *
+ * In case of parallel::distributed::Triangulation @p predicate will be called
+ * only for locally owned and ghost cells. The resulting index set may contain
+ * DoFs that are associated to the locally owned or ghost cells, but are not
+ * owned by the current MPI core.
+ */
+ template <typename DoFHandlerType>
+ IndexSet
+ extract_dofs_with_support_on(const DoFHandlerType &dof_handler,
+ const std::function< bool(const typename DoFHandlerType::active_cell_iterator &)> &predicate,
+ const ConstraintMatrix &cm = ConstraintMatrix());
+
/**
* Extract a vector that represents the constant modes of the DoFHandler for
* the components chosen by <tt>component_mask</tt> (see
+ template <typename DoFHandlerType>
+ IndexSet
+ extract_dofs_with_support_on(const DoFHandlerType &dof_handler,
+ const std::function< bool(const typename DoFHandlerType::active_cell_iterator &)> &predicate,
+ const ConstraintMatrix &cm)
+ {
+ const std::function< bool(const typename DoFHandlerType::active_cell_iterator &)> predicate_local
+ = [=] (const typename DoFHandlerType::active_cell_iterator & cell) -> bool
+ {
+ return cell->is_locally_owned() && predicate(cell);
+ };
+
+ std::vector<types::global_dof_index> local_dof_indices;
+ local_dof_indices.reserve (max_dofs_per_cell(dof_handler));
+
+ // Get all the dofs that live on the subdomain:
+ std::set<types::global_dof_index> predicate_dofs;
+
+ typename DoFHandlerType::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ if (!cell->is_artificial() && predicate(cell))
+ {
+ local_dof_indices.resize (cell->get_fe().dofs_per_cell);
+ cell->get_dof_indices (local_dof_indices);
+ predicate_dofs.insert(local_dof_indices.begin(),
+ local_dof_indices.end());
+ }
+
+ // Get halo layer and accumulate its DoFs
+ std::set<types::global_dof_index> halo_dofs;
+
+ const std::vector<typename DoFHandlerType::active_cell_iterator> halo_cells =
+ GridTools::compute_active_cell_halo_layer(dof_handler, predicate_local);
+ for (typename std::vector<typename DoFHandlerType::active_cell_iterator>::const_iterator it = halo_cells.begin();
+ it != halo_cells.end(); ++it)
+ {
+ // skip halo cells that satisfy the given predicate and thereby
+ // shall be a part of the index set on another MPI core.
+ // Those could only be ghost cells with p::d::Tria.
+ if (predicate(*it))
+ {
+ Assert ((*it)->is_ghost(),
+ ExcInternalError());
+ continue;
+ }
+
+ const unsigned int dofs_per_cell = (*it)->get_fe().dofs_per_cell;
+ local_dof_indices.resize (dofs_per_cell);
+ (*it)->get_dof_indices (local_dof_indices);
+ halo_dofs.insert(local_dof_indices.begin(),
+ local_dof_indices.end());
+ }
+
+ // A complication coming from the fact that DoFs living on locally
+ // owned cells which satisfy predicate may participate in constrains for
+ // DoFs outside of this region.
+ if (cm.n_constraints() > 0)
+ {
+ const std::vector<std::pair<types::global_dof_index,double> > *line_ptr;
+
+ // Remove DoFs that are part of constrains for DoFs outside
+ // of predicate. Since we will subtract halo_dofs from predicate_dofs,
+ // achieve this by extending halo_dofs with DoFs to which
+ // halo_dofs are constrained.
+ std::set<types::global_dof_index> extra_halo;
+ for (std::set<types::global_dof_index>::iterator it = halo_dofs.begin();
+ it != halo_dofs.end(); ++it)
+ {
+ line_ptr = cm.get_constraint_entries(*it);
+ // if halo DoF is constrained, add all DoFs to which it's constrained.
+ if (line_ptr!=NULL)
+ {
+ const unsigned int line_size = line_ptr->size();
+ for (unsigned int j=0; j<line_size; ++j)
+ extra_halo.insert((*line_ptr)[j].first);
+ }
+ }
+
+ halo_dofs.insert(extra_halo.begin(),
+ extra_halo.end());
+ }
+
+ // Rework our std::set's into IndexSet and subtract halo DoFs from the
+ // predicate_dofs:
+ IndexSet support_set(dof_handler.n_dofs());
+ support_set.add_indices(predicate_dofs.begin(), predicate_dofs.end());
+ support_set.compress();
+
+ IndexSet halo_set(dof_handler.n_dofs());
+ halo_set.add_indices(halo_dofs.begin(), halo_dofs.end());
+ halo_set.compress();
+
+ support_set.subtract_set(halo_set);
+
+ // we intentionally do not want to limit the output to locally owned DoFs.
+ // In other words, the output may contain DoFs that are on the locally
+ // owned cells, but are not owned by this MPI core.
+
+ return support_set;
+ }
+
+
namespace internal
{
std::vector<bool> &,
const std::set<types::boundary_id> &);
+ template
+ IndexSet
+ DoFTools::extract_dofs_with_support_on<DoFHandler<deal_II_dimension> >
+ (const DoFHandler<deal_II_dimension> &,
+ const std::function< bool(const typename DoFHandler<deal_II_dimension>::active_cell_iterator &)> &,
+ const ConstraintMatrix&);
+
template
void
DoFTools::extract_hanging_node_dofs
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 DoFTools::extract_dofs_with_support_on()
+
+#include "../tests.h"
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/point.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+
+#include <fstream>
+#include <sstream>
+#include <iomanip>
+#include <list>
+#include <set>
+#include <cstdio>
+
+
+template<int dim>
+bool
+pred_d(const typename Triangulation<dim>::active_cell_iterator &cell)
+{
+ return (cell->center()(0) < 0.49 &&
+ cell->center()(1) < 0.49);
+}
+
+template<int dim>
+bool
+pred_r(const typename Triangulation<dim>::active_cell_iterator &cell)
+{
+ return (cell->center()(0) < 0.49 &&
+ cell->center()(0) > 0.25 &&
+ cell->center()(1) < 0.49);
+}
+
+
+
+template <int dim>
+void test (const unsigned int flag)
+{
+
+ // Setup system
+ Triangulation<dim> triangulation;
+
+ GridGenerator::hyper_rectangle (triangulation,
+ Point<dim>(0,0),
+ Point<dim>(1,1));
+
+ if (flag == 0)
+ triangulation.refine_global(2);
+ else
+ triangulation.refine_global(1);
+
+ // Extra refinement to generate hanging nodes
+ for (typename Triangulation<dim>::active_cell_iterator
+ cell = triangulation.begin_active();
+ cell != triangulation.end(); ++cell)
+ if ( (flag==1 && pred_d<dim>(cell)) ||
+ (flag==2 && !pred_d<dim>(cell)) )
+ cell->set_refine_flag ();
+
+ triangulation.prepare_coarsening_and_refinement();
+ triangulation.execute_coarsening_and_refinement ();
+
+ DoFHandler<dim> dh (triangulation);
+
+ FE_Q<dim> fe(2);
+ dh.distribute_dofs (fe);
+
+ ConstraintMatrix cm;
+ DoFTools::make_hanging_node_constraints(dh,cm);
+
+ IndexSet support = DoFTools::extract_dofs_with_support_on(dh, pred_d<dim>, cm);
+ support.print(deallog);
+
+ // print grid and DoFs for visual inspection
+ if (false)
+ {
+ std::cout << "-------------------- " << Utilities::int_to_string(dim) << Utilities::int_to_string(flag) << std::endl;
+ cm.print(std::cout);
+
+ std::map<types::global_dof_index, Point<dim> > support_points;
+ MappingQ1<dim> mapping;
+ DoFTools::map_dofs_to_support_points(mapping, dh, support_points);
+
+ const std::string filename =
+ "grid" + Utilities::int_to_string(dim) + Utilities::int_to_string(flag) + ".gp";
+ std::ofstream f(filename.c_str());
+
+ f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" << std::endl
+ << "set output \"grid" << Utilities::int_to_string(dim) << Utilities::int_to_string(flag) << ".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;
+ }
+
+ dh.clear();
+}
+
+int main(int argc, char **argv)
+{
+ initlog();
+
+ test<2>(0);
+ test<2>(1);
+ test<2>(2);
+
+ return 0;
+}
--- /dev/null
+
+DEAL::{[0,8], [12,14], [17,18], 20, 24}
+DEAL::{[21,42]}
+DEAL::{0, 4, 6, 8}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+//
+
+// similar to dof_tools_23, but for the p::d::Tria. As an MPI test, make sure
+// that the total number of shape functions that are non-zero within the domain
+// is the same.
+
+#include "../tests.h"
+#include <deal.II/distributed/tria.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/point.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <fstream>
+#include <sstream>
+#include <iomanip>
+#include <list>
+#include <set>
+#include <cstdio>
+
+
+template<int dim>
+bool
+pred_d(const typename parallel::distributed::Triangulation<dim>::active_cell_iterator &cell)
+{
+ return (cell->center()(0) < 0.5 &&
+ cell->center()(1) < 0.5);
+}
+
+
+std::string output_name(const unsigned int flag, const unsigned int subdomain)
+{
+ return "output" + Utilities::int_to_string(flag) + "_" +
+ Utilities::int_to_string(subdomain) +
+ ".vtu";
+}
+
+
+template <int dim>
+void test (const unsigned int flag)
+{
+ // Setup system
+ parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+
+ GridGenerator::hyper_rectangle (triangulation,
+ Point<dim>(0,0),
+ Point<dim>(1,1));
+
+ if (flag == 0)
+ triangulation.refine_global(2);
+ else
+ triangulation.refine_global(1);
+
+ // Extra refinement to generate hanging nodes
+ for (typename parallel::distributed::Triangulation<dim>::active_cell_iterator
+ cell = triangulation.begin_active();
+ cell != triangulation.end(); ++cell)
+ if (cell->is_locally_owned() &&
+ ((flag==1 && pred_d<dim>(cell)) ||
+ (flag==2 && !pred_d<dim>(cell))) )
+ cell->set_refine_flag ();
+
+ triangulation.prepare_coarsening_and_refinement();
+ triangulation.execute_coarsening_and_refinement ();
+
+ if (flag > 0)
+ triangulation.refine_global(1);
+
+ DoFHandler<dim> dh (triangulation);
+
+ FE_Q<dim> fe(2);
+ dh.distribute_dofs (fe);
+
+ IndexSet locally_relevant_set;
+ DoFTools::extract_locally_relevant_dofs (dh,
+ locally_relevant_set);
+
+ ConstraintMatrix cm;
+ cm.reinit (locally_relevant_set);
+ DoFTools::make_hanging_node_constraints(dh,cm);
+ cm.close ();
+
+ const IndexSet support = DoFTools::extract_dofs_with_support_on(dh, pred_d<dim>, cm);
+ const IndexSet support_local = support & dh.locally_owned_dofs();
+
+ deallog << support.n_elements() << std::endl;
+
+ // now accumulate the number of indices and make sure it's the same for
+ // various runs with different number of MPI cores
+ const unsigned int dofs_support = Utilities::MPI::sum(support_local.n_elements(), MPI_COMM_WORLD);
+ if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+ {
+ deallog << "accumulated: " << std::endl;
+ deallog << dofs_support << std::endl;
+ }
+
+ // print grid and DoFs for visual inspection
+ if (false)
+ {
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler (dh);
+
+ std::vector<LinearAlgebra::distributed::Vector<double> > shape_functions(dh.n_dofs());
+ for (unsigned int i = 0; i < dh.n_dofs(); ++i)
+ {
+ LinearAlgebra::distributed::Vector<double> &s = shape_functions[i];
+ s.reinit(dh.locally_owned_dofs(),
+ locally_relevant_set,
+ MPI_COMM_WORLD);
+ s = 0.;
+ if (dh.locally_owned_dofs().is_element(i))
+ s[i] = 1.0;
+ s.compress(VectorOperation::insert);
+ cm.distribute(s);
+ s.update_ghost_values();
+
+ data_out.add_data_vector (s,
+ std::string("N_") +
+ Utilities::int_to_string(i));
+ }
+
+ Vector<float> subdomain (triangulation.n_active_cells());
+ for (unsigned int i=0; i<subdomain.size(); ++i)
+ subdomain(i) = triangulation.locally_owned_subdomain();
+ data_out.add_data_vector (subdomain, "subdomain");
+ data_out.build_patches ();
+
+ const std::string filename = output_name(flag, triangulation.locally_owned_subdomain());
+
+ std::ofstream output (filename.c_str());
+ 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.push_back (output_name(flag, i));
+
+ const std::string master_name = "output" + Utilities::int_to_string(flag) + ".pvtu";
+ std::ofstream pvtu_master (master_name.c_str());
+ data_out.write_pvtu_record (pvtu_master, filenames);
+ }
+
+ }
+
+ dh.clear();
+}
+
+int main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, testing_max_num_threads());
+ const unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+
+ MPILogInitAll log;
+
+ deallog.push(Utilities::int_to_string(myid));
+
+ test<2>(0);
+ test<2>(1);
+ test<2>(2);
+
+ deallog.pop();
+}
--- /dev/null
+
+DEAL:0:0::16
+DEAL:0:0::accumulated:
+DEAL:0:0::16
+DEAL:0:0::76
+DEAL:0:0::accumulated:
+DEAL:0:0::76
+DEAL:0:0::16
+DEAL:0:0::accumulated:
+DEAL:0:0::16
--- /dev/null
+
+DEAL:0:0::16
+DEAL:0:0::accumulated:
+DEAL:0:0::16
+DEAL:0:0::61
+DEAL:0:0::accumulated:
+DEAL:0:0::76
+DEAL:0:0::16
+DEAL:0:0::accumulated:
+DEAL:0:0::16
+
+DEAL:1:1::21
+DEAL:1:1::64
+DEAL:1:1::15
+
+
+DEAL:2:2::9
+DEAL:2:2::27
+DEAL:2:2::9
+
--- /dev/null
+
+DEAL:0:0::16
+DEAL:0:0::accumulated:
+DEAL:0:0::16
+DEAL:0:0::49
+DEAL:0:0::accumulated:
+DEAL:0:0::76
+DEAL:0:0::16
+DEAL:0:0::accumulated:
+DEAL:0:0::16
+
+DEAL:1:1::15
+DEAL:1:1::69
+DEAL:1:1::15
+
+
+DEAL:2:2::0
+DEAL:2:2::46
+DEAL:2:2::15
+
+
+DEAL:3:3::15
+DEAL:3:3::45
+DEAL:3:3::9
+
+
+DEAL:4:4::9
+DEAL:4:4::9
+DEAL:4:4::0
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+//
+
+// same as dof_tools_32 test DoFTools::extract_dofs_with_support_on() but
+// for a slightly different configuration
+
+#include "../tests.h"
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/point.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+
+#include <fstream>
+#include <sstream>
+#include <iomanip>
+#include <list>
+#include <set>
+#include <cstdio>
+
+
+template<int dim>
+bool
+pred_left(const typename Triangulation<dim>::active_cell_iterator &cell)
+{
+ return (cell->center()(0) < 0.49);
+}
+
+template<int dim>
+bool
+pred_right(const typename Triangulation<dim>::active_cell_iterator &cell)
+{
+ return (cell->center()(0) > 0.51);
+}
+
+
+template<int dim>
+bool
+pred_r(const typename Triangulation<dim>::active_cell_iterator &cell)
+{
+ return (cell->center()(0) < 0.49 &&
+ cell->center()(1) < 0.49) ||
+ (cell->center()(0) > 0.49 &&
+ cell->center()(1) > 0.49);
+}
+
+
+
+template <int dim>
+void test (const bool left = true)
+{
+
+ // Setup system
+ Triangulation<dim> triangulation;
+
+ GridGenerator::hyper_rectangle (triangulation,
+ Point<dim>(0,0),
+ Point<dim>(1,1));
+
+ triangulation.refine_global(1);
+
+ // Extra refinement to generate hanging nodes
+ for (typename Triangulation<dim>::active_cell_iterator
+ cell = triangulation.begin_active();
+ cell != triangulation.end(); ++cell)
+ if (pred_r<dim>(cell))
+ cell->set_refine_flag ();
+
+ triangulation.prepare_coarsening_and_refinement();
+ triangulation.execute_coarsening_and_refinement ();
+
+ DoFHandler<dim> dh (triangulation);
+
+ FE_Q<dim> fe(2);
+ dh.distribute_dofs (fe);
+
+ ConstraintMatrix cm;
+ DoFTools::make_hanging_node_constraints(dh,cm);
+
+ IndexSet support = left ?
+ DoFTools::extract_dofs_with_support_on(dh, pred_left<dim>, cm) :
+ DoFTools::extract_dofs_with_support_on(dh, pred_right<dim>, cm);
+ support.print(deallog);
+
+ // print grid and DoFs for visual inspection
+ if (false)
+ {
+ std::cout << "-------------------- " << Utilities::int_to_string(dim) << std::endl;
+ cm.print(std::cout);
+
+ std::map<types::global_dof_index, Point<dim> > support_points;
+ MappingQ1<dim> mapping;
+ DoFTools::map_dofs_to_support_points(mapping, dh, support_points);
+
+ const std::string filename =
+ "grid" + Utilities::int_to_string(dim) + ".gp";
+ std::ofstream f(filename.c_str());
+
+ f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" << std::endl
+ << "set output \"grid" << Utilities::int_to_string(dim) << ".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;
+ }
+
+ dh.clear();
+}
+
+int main(int argc, char **argv)
+{
+ initlog();
+
+ test<2>(true);
+ test<2>(false);
+
+ return 0;
+}
--- /dev/null
+
+DEAL::{[9,10], 12, [14,38]}
+DEAL::{1, 3, [5,8], [39,60]}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 DoFTools::extract_dofs_with_support_on() for calculation of the RHS
+// function in parallel
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/point.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <fstream>
+#include <sstream>
+#include <iomanip>
+#include <list>
+#include <set>
+#include <cstdio>
+
+
+std::string output_name(const unsigned int subdomain)
+{
+ return "output_" +
+ Utilities::int_to_string(subdomain) +
+ ".vtu";
+}
+
+template<int dim>
+bool
+pred_d(const typename Triangulation<dim>::active_cell_iterator &cell)
+{
+ return (cell->center()(0) < 0.49);
+}
+
+template<int dim>
+bool
+pred_r(const typename Triangulation<dim>::active_cell_iterator &cell)
+{
+ return (cell->center()(0) < 0.49 &&
+ cell->center()(1) < 0.49) ||
+ (cell->center()(0) > 0.49 &&
+ cell->center()(1) > 0.49);
+}
+
+
+
+template <int dim>
+void test ()
+{
+
+ // Setup system
+ parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+
+ GridGenerator::hyper_rectangle (triangulation,
+ Point<dim>(0,0),
+ Point<dim>(1,1));
+
+ triangulation.refine_global(1);
+
+ // Extra refinement to generate hanging nodes
+ for (typename Triangulation<dim>::active_cell_iterator
+ cell = triangulation.begin_active();
+ cell != triangulation.end(); ++cell)
+ if (cell->is_locally_owned() && pred_r<dim>(cell))
+ cell->set_refine_flag ();
+
+ triangulation.prepare_coarsening_and_refinement();
+ triangulation.execute_coarsening_and_refinement ();
+
+ DoFHandler<dim> dh (triangulation);
+
+ FE_Q<dim> fe(2);
+ dh.distribute_dofs (fe);
+
+ const IndexSet &locally_owned_set = dh.locally_owned_dofs();
+ IndexSet locally_relevant_set;
+ DoFTools::extract_locally_relevant_dofs (dh,
+ locally_relevant_set);
+
+ ConstraintMatrix cm;
+ cm.reinit (locally_relevant_set);
+ DoFTools::make_hanging_node_constraints(dh,cm);
+ cm.close ();
+
+ std::vector<types::global_dof_index> local_dof_indices(fe.dofs_per_cell);
+
+ // get support on the predicate
+ IndexSet support = DoFTools::extract_dofs_with_support_on(dh, pred_d<dim>, cm);
+ IndexSet local_support = support & locally_owned_set;
+
+ // rhs vectors:
+ LinearAlgebra::distributed::Vector<double> sparse_rhs;
+ LinearAlgebra::distributed::Vector<double> rhs;
+ sparse_rhs.reinit(locally_owned_set, locally_relevant_set, MPI_COMM_WORLD);
+ rhs.reinit( locally_owned_set, locally_relevant_set, MPI_COMM_WORLD);
+
+ sparse_rhs.zero_out_ghosts();
+ rhs.zero_out_ghosts();
+
+ // assemble RHS which has a local support:
+ const std::function< double(const Point<dim> &)> rhs_func
+ = [=] (const Point<dim> &p) -> double
+ {
+ return p[0] > 0.5 ? 0. : 0.5 - p[0];
+ };
+
+ Vector<double> local_rhs(fe.dofs_per_cell);
+ QGauss<dim> quadrature(3);
+ FEValues<dim> fe_values(fe, quadrature,
+ update_values | update_JxW_values | update_quadrature_points);
+ for (typename DoFHandler<dim>::active_cell_iterator
+ cell = dh.begin_active();
+ cell != dh.end(); ++cell)
+ if (cell->is_locally_owned() && pred_d<dim>(cell))
+ {
+ fe_values.reinit(cell);
+ cell->get_dof_indices (local_dof_indices);
+
+ const std::vector<Point<dim>> &q_points = fe_values.get_quadrature_points();
+
+ local_rhs = 0.;
+ for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
+ for (unsigned int q = 0; q < quadrature.size(); ++q)
+ local_rhs[i] += fe_values.shape_value (i, q) *
+ rhs_func(q_points[q]) *
+ fe_values.JxW (q);
+
+ cm.distribute_local_to_global (local_rhs,
+ local_dof_indices,
+ rhs);
+
+ // copy-paste of CM distribute_local_to_global and
+ // add is_element() checks:
+ auto local_vector_begin = local_rhs.begin();
+ const auto local_vector_end = local_rhs.end();
+ auto local_indices_begin = local_dof_indices.begin();
+ const std::vector<std::pair<types::global_dof_index,double> > *line_ptr;
+ for ( ; local_vector_begin != local_vector_end;
+ ++local_vector_begin, ++local_indices_begin)
+ {
+ line_ptr = cm.get_constraint_entries(*local_indices_begin);
+ if (line_ptr==NULL) // unconstrained
+ {
+ if (support.is_element(*local_indices_begin))
+ sparse_rhs(*local_indices_begin) += *local_vector_begin;
+ }
+ else
+ {
+ const unsigned int line_size = line_ptr->size();
+ for (unsigned int j=0; j<line_size; ++j)
+ if (support.is_element((*line_ptr)[j].first))
+ sparse_rhs((*line_ptr)[j].first)
+ += *local_vector_begin * (*line_ptr)[j].second;
+ }
+ }
+ /*
+ cm.distribute_local_to_global (local_rhs,
+ local_dof_indices,
+ sparse_rhs);
+ */
+ }
+
+ rhs.compress(VectorOperation::add);
+ sparse_rhs.compress(VectorOperation::add);
+
+ rhs.update_ghost_values();
+ sparse_rhs.update_ghost_values();
+
+ rhs.add(-1., sparse_rhs);
+
+ // print grid and DoFs for visual inspection
+ if (false)
+ {
+ const unsigned int this_mpi_process = dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+ const unsigned int n_mpi_processes = dealii::Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+
+ for (unsigned int i = 0; i < n_mpi_processes; ++i)
+ {
+ MPI_Barrier(MPI_COMM_WORLD);
+ if (i == this_mpi_process)
+ {
+ std::cout << "-------------------- " << this_mpi_process << std::endl;
+ cm.print(std::cout);
+
+ std::cout << "local support:" << std::endl;
+ local_support.print(std::cout);
+ std::cout << "support:" << std::endl;
+ support.print(std::cout);
+ std::cout << "locally owned:" << std::endl;
+ locally_owned_set.print(std::cout);
+ std::cout << "locally relevant set:" << std::endl;
+ locally_relevant_set.print(std::cout);
+ }
+ }
+
+ {
+ std::map<types::global_dof_index, Point<dim> > support_points;
+ MappingQ1<dim> mapping;
+ DoFTools::map_dofs_to_support_points(mapping, dh, support_points);
+
+ const std::string filename =
+ "grid" + dealii::Utilities::int_to_string(n_mpi_processes) + dealii::Utilities::int_to_string(this_mpi_process);
+ std::ofstream f((filename+".gp").c_str());
+
+ f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" << std::endl
+ << "set output \"" << 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;
+
+ }
+
+ {
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler (dh);
+
+ std::vector<LinearAlgebra::distributed::Vector<double> > shape_functions(dh.n_dofs());
+ for (unsigned int i = 0; i < dh.n_dofs(); ++i)
+ {
+ LinearAlgebra::distributed::Vector<double> sl(locally_owned_set, MPI_COMM_WORLD);
+ sl = 0.;
+ if (locally_owned_set.is_element(i))
+ sl[i] = 1.0;
+ cm.distribute(sl);
+
+ LinearAlgebra::distributed::Vector<double> &s = shape_functions[i];
+ s.reinit(locally_owned_set,
+ locally_relevant_set,
+ MPI_COMM_WORLD);
+ s = 0.;
+ s = sl;
+
+ data_out.add_data_vector (s,
+ std::string("N_") +
+ dealii::Utilities::int_to_string(i));
+ }
+
+ Vector<float> subdomain (triangulation.n_active_cells());
+ for (unsigned int i=0; i<subdomain.size(); ++i)
+ subdomain(i) = triangulation.locally_owned_subdomain();
+ data_out.add_data_vector (subdomain, "subdomain");
+ data_out.build_patches ();
+
+ const std::string filename = output_name(triangulation.locally_owned_subdomain());
+
+ std::ofstream output (filename.c_str());
+ data_out.write_vtu (output);
+
+ if (dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+ {
+ std::vector<std::string> filenames;
+ for (unsigned int i=0; i<dealii::Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); ++i)
+ filenames.push_back (output_name(i));
+
+ const std::string master_name = "output.pvtu";
+ std::ofstream pvtu_master (master_name.c_str());
+ data_out.write_pvtu_record (pvtu_master, filenames);
+ }
+ }
+ }
+
+ for (unsigned int i = 0; i < locally_owned_set.n_elements(); ++i)
+ {
+ const unsigned int ind = locally_owned_set.nth_index_in_set(i);
+ const double v = rhs[ind];
+ AssertThrow(std::abs(v) < 1e-12,
+ ExcMessage(
+ "Element " + std::to_string(ind) +
+ " has an error " + std::to_string(v) +
+ " on the process " + std::to_string(dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD))
+ ));
+ }
+
+ if (dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+ deallog << "Ok" << std::endl;
+
+
+ dh.clear();
+}
+
+int main(int argc, char **argv)
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile,/*do not print job id*/false);
+ deallog.depth_console(0);
+ deallog.precision(4);
+
+ dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+ test<2>();
+
+ return 0;
+}
projects / dealii.git / commitdiff
