From: Denis Davydov Date: Thu, 23 Mar 2017 16:48:47 +0000 (+0100) Subject: add a function to get DoFs that have support only within the given domain X-Git-Tag: v9.0.0-rc1~1688^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c0c9940bc2916998fd2cbbd47bd9d42a09df1f06;p=dealii.git add a function to get DoFs that have support only within the given domain --- diff --git a/doc/doxygen/images/extract_dofs_with_support_on.png b/doc/doxygen/images/extract_dofs_with_support_on.png new file mode 100644 index 0000000000..9d6c3328e9 Binary files /dev/null and b/doc/doxygen/images/extract_dofs_with_support_on.png differ diff --git a/doc/news/changes/minor/20170324DenisDavydov b/doc/news/changes/minor/20170324DenisDavydov new file mode 100644 index 0000000000..462f1a659c --- /dev/null +++ b/doc/news/changes/minor/20170324DenisDavydov @@ -0,0 +1,4 @@ +New: DoFTools::extract_dofs_with_support_on() returns a set of degrees of +freedom that have support on a subset of locally owned cells. +
+(Denis Davydov, 2017/03/24) diff --git a/include/deal.II/dofs/dof_tools.h b/include/deal.II/dofs/dof_tools.h index aa4d0ababb..d0eead0241 100644 --- a/include/deal.II/dofs/dof_tools.h +++ b/include/deal.II/dofs/dof_tools.h @@ -1436,6 +1436,39 @@ namespace DoFTools std::vector &selected_dofs, const std::set &boundary_ids = std::set()); + /** + * Extract all degrees of freedom (DoFs) that have support only within cells, + * for which @p predicate is true. The result is returned as + * an IndexSet. + * + * Consider the following FE space where predicate returns true + * 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 [9,10], 12, [14,38]. + * 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 + 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 component_mask (see diff --git a/source/dofs/dof_tools.cc b/source/dofs/dof_tools.cc index 19314c16c1..796335ef7f 100644 --- a/source/dofs/dof_tools.cc +++ b/source/dofs/dof_tools.cc @@ -765,6 +765,110 @@ namespace DoFTools + template + 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 local_dof_indices; + local_dof_indices.reserve (max_dofs_per_cell(dof_handler)); + + // Get all the dofs that live on the subdomain: + std::set 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 halo_dofs; + + const std::vector halo_cells = + GridTools::compute_active_cell_halo_layer(dof_handler, predicate_local); + for (typename std::vector::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 > *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 extra_halo; + for (std::set::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 &, const std::set &); + template + IndexSet + DoFTools::extract_dofs_with_support_on > + (const DoFHandler &, + const std::function< bool(const typename DoFHandler::active_cell_iterator &)> &, + const ConstraintMatrix&); + template void DoFTools::extract_hanging_node_dofs diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_01.cc b/tests/dofs/dof_tools_extract_dofs_with_support_on_01.cc new file mode 100644 index 0000000000..741020cd8c --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_01.cc @@ -0,0 +1,137 @@ +// --------------------------------------------------------------------- +// +// 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + + +template +bool +pred_d(const typename Triangulation::active_cell_iterator &cell) +{ + return (cell->center()(0) < 0.49 && + cell->center()(1) < 0.49); +} + +template +bool +pred_r(const typename Triangulation::active_cell_iterator &cell) +{ + return (cell->center()(0) < 0.49 && + cell->center()(0) > 0.25 && + cell->center()(1) < 0.49); +} + + + +template +void test (const unsigned int flag) +{ + + // Setup system + Triangulation triangulation; + + GridGenerator::hyper_rectangle (triangulation, + Point(0,0), + Point(1,1)); + + if (flag == 0) + triangulation.refine_global(2); + else + triangulation.refine_global(1); + + // Extra refinement to generate hanging nodes + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + if ( (flag==1 && pred_d(cell)) || + (flag==2 && !pred_d(cell)) ) + cell->set_refine_flag (); + + triangulation.prepare_coarsening_and_refinement(); + triangulation.execute_coarsening_and_refinement (); + + DoFHandler dh (triangulation); + + FE_Q 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, 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 > support_points; + MappingQ1 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; +} diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_01.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_01.output new file mode 100644 index 0000000000..d326d0283a --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_01.output @@ -0,0 +1,4 @@ + +DEAL::{[0,8], [12,14], [17,18], 20, 24} +DEAL::{[21,42]} +DEAL::{0, 4, 6, 8} diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_02.cc b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.cc new file mode 100644 index 0000000000..d74ba742fb --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.cc @@ -0,0 +1,185 @@ +// --------------------------------------------------------------------- +// +// 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + + +template +bool +pred_d(const typename parallel::distributed::Triangulation::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 +void test (const unsigned int flag) +{ + // Setup system + parallel::distributed::Triangulation triangulation(MPI_COMM_WORLD); + + GridGenerator::hyper_rectangle (triangulation, + Point(0,0), + Point(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::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + if (cell->is_locally_owned() && + ((flag==1 && pred_d(cell)) || + (flag==2 && !pred_d(cell))) ) + cell->set_refine_flag (); + + triangulation.prepare_coarsening_and_refinement(); + triangulation.execute_coarsening_and_refinement (); + + if (flag > 0) + triangulation.refine_global(1); + + DoFHandler dh (triangulation); + + FE_Q 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, 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 data_out; + data_out.attach_dof_handler (dh); + + std::vector > shape_functions(dh.n_dofs()); + for (unsigned int i = 0; i < dh.n_dofs(); ++i) + { + LinearAlgebra::distributed::Vector &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 subdomain (triangulation.n_active_cells()); + for (unsigned int i=0; i filenames; + for (unsigned int i=0; i(0); + test<2>(1); + test<2>(2); + + deallog.pop(); +} diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=1.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=1.output new file mode 100644 index 0000000000..ad14a0a639 --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=1.output @@ -0,0 +1,10 @@ + +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 diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=3.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=3.output new file mode 100644 index 0000000000..24d5d3f17a --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=3.output @@ -0,0 +1,20 @@ + +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 + diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=5.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=5.output new file mode 100644 index 0000000000..4a18c5f509 --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=5.output @@ -0,0 +1,30 @@ + +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 + diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_03.cc b/tests/dofs/dof_tools_extract_dofs_with_support_on_03.cc new file mode 100644 index 0000000000..e19e4ec036 --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_03.cc @@ -0,0 +1,143 @@ +// --------------------------------------------------------------------- +// +// 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + + +template +bool +pred_left(const typename Triangulation::active_cell_iterator &cell) +{ + return (cell->center()(0) < 0.49); +} + +template +bool +pred_right(const typename Triangulation::active_cell_iterator &cell) +{ + return (cell->center()(0) > 0.51); +} + + +template +bool +pred_r(const typename Triangulation::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 +void test (const bool left = true) +{ + + // Setup system + Triangulation triangulation; + + GridGenerator::hyper_rectangle (triangulation, + Point(0,0), + Point(1,1)); + + triangulation.refine_global(1); + + // Extra refinement to generate hanging nodes + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + if (pred_r(cell)) + cell->set_refine_flag (); + + triangulation.prepare_coarsening_and_refinement(); + triangulation.execute_coarsening_and_refinement (); + + DoFHandler dh (triangulation); + + FE_Q 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, cm) : + DoFTools::extract_dofs_with_support_on(dh, pred_right, 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 > support_points; + MappingQ1 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; +} diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_03.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_03.output new file mode 100644 index 0000000000..76ee14d0e0 --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_03.output @@ -0,0 +1,3 @@ + +DEAL::{[9,10], 12, [14,38]} +DEAL::{1, 3, [5,8], [39,60]} diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_04.cc b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.cc new file mode 100644 index 0000000000..7126ab0558 --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.cc @@ -0,0 +1,327 @@ +// --------------------------------------------------------------------- +// +// 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + + +std::string output_name(const unsigned int subdomain) +{ + return "output_" + + Utilities::int_to_string(subdomain) + + ".vtu"; +} + +template +bool +pred_d(const typename Triangulation::active_cell_iterator &cell) +{ + return (cell->center()(0) < 0.49); +} + +template +bool +pred_r(const typename Triangulation::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 +void test () +{ + + // Setup system + parallel::distributed::Triangulation triangulation(MPI_COMM_WORLD); + + GridGenerator::hyper_rectangle (triangulation, + Point(0,0), + Point(1,1)); + + triangulation.refine_global(1); + + // Extra refinement to generate hanging nodes + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + if (cell->is_locally_owned() && pred_r(cell)) + cell->set_refine_flag (); + + triangulation.prepare_coarsening_and_refinement(); + triangulation.execute_coarsening_and_refinement (); + + DoFHandler dh (triangulation); + + FE_Q 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 local_dof_indices(fe.dofs_per_cell); + + // get support on the predicate + IndexSet support = DoFTools::extract_dofs_with_support_on(dh, pred_d, cm); + IndexSet local_support = support & locally_owned_set; + + // rhs vectors: + LinearAlgebra::distributed::Vector sparse_rhs; + LinearAlgebra::distributed::Vector 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 &)> rhs_func + = [=] (const Point &p) -> double + { + return p[0] > 0.5 ? 0. : 0.5 - p[0]; + }; + + Vector local_rhs(fe.dofs_per_cell); + QGauss quadrature(3); + FEValues fe_values(fe, quadrature, + update_values | update_JxW_values | update_quadrature_points); + for (typename DoFHandler::active_cell_iterator + cell = dh.begin_active(); + cell != dh.end(); ++cell) + if (cell->is_locally_owned() && pred_d(cell)) + { + fe_values.reinit(cell); + cell->get_dof_indices (local_dof_indices); + + const std::vector> &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 > *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 > support_points; + MappingQ1 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 data_out; + data_out.attach_dof_handler (dh); + + std::vector > shape_functions(dh.n_dofs()); + for (unsigned int i = 0; i < dh.n_dofs(); ++i) + { + LinearAlgebra::distributed::Vector 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 &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 subdomain (triangulation.n_active_cells()); + for (unsigned int i=0; i filenames; + for (unsigned int i=0; i(); + + return 0; +} diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_04.mpirun=2.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.mpirun=2.output new file mode 100644 index 0000000000..04c93dbf6b --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.mpirun=2.output @@ -0,0 +1 @@ +DEAL::Ok diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_04.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.output new file mode 100644 index 0000000000..04c93dbf6b --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.output @@ -0,0 +1 @@ +DEAL::Ok