// the output array dofs_by_component lists for each dof the
// corresponding vector component. if the DoFHandler is based on a
// parallel distributed triangulation then the output array is index by
- // dof.locally_owned_dofs().index_within_set(indices[i])
+ // dof_handler.locally_owned_dofs().index_within_set(indices[i])
//
// if an element is non-primitive then we assign to each degree of
// freedom the following component:
// in the component_mask that corresponds to this shape function
template <int dim, int spacedim>
void
- get_component_association(const DoFHandler<dim, spacedim> &dof,
- const ComponentMask &component_mask,
- std::vector<unsigned char> &dofs_by_component)
+ get_component_association(
+ const DoFHandler<dim, spacedim> &dof_handler,
+ const ComponentMask &component_mask,
+ std::vector<unsigned char> &dofs_by_component,
+ const unsigned int mg_level = numbers::invalid_unsigned_int)
{
const dealii::hp::FECollection<dim, spacedim> &fe_collection =
- dof.get_fe_collection();
+ dof_handler.get_fe_collection();
Assert(fe_collection.n_components() < 256, ExcNotImplemented());
- Assert(dofs_by_component.size() == dof.n_locally_owned_dofs(),
- ExcDimensionMismatch(dofs_by_component.size(),
- dof.n_locally_owned_dofs()));
+
+ const auto &locally_owned_dofs =
+ (mg_level == numbers::invalid_unsigned_int) ?
+ dof_handler.locally_owned_dofs() :
+ dof_handler.locally_owned_mg_dofs(mg_level);
+
+ AssertDimension(dofs_by_component.size(),
+ locally_owned_dofs.n_elements());
// next set up a table for the degrees of freedom on each of the
// cells (regardless of the fact whether it is listed in the
// then loop over all cells and do the work
std::vector<types::global_dof_index> indices;
- for (const auto &c :
- dof.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
- {
- const types::fe_index fe_index = c->active_fe_index();
- const unsigned int dofs_per_cell = c->get_fe().n_dofs_per_cell();
- indices.resize(dofs_per_cell);
- c->get_dof_indices(indices);
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- if (dof.locally_owned_dofs().is_element(indices[i]))
- dofs_by_component[dof.locally_owned_dofs().index_within_set(
- indices[i])] = local_component_association[fe_index][i];
- }
+
+ const auto runner = [&](const auto &task) {
+ if (mg_level == numbers::invalid_unsigned_int)
+ {
+ for (const auto &cell : dof_handler.active_cell_iterators() |
+ IteratorFilters::LocallyOwnedCell())
+ {
+ indices.resize(cell->get_fe().n_dofs_per_cell());
+ cell->get_dof_indices(indices);
+
+ task(cell);
+ }
+ }
+ else
+ {
+ for (const auto &cell :
+ dof_handler.cell_iterators_on_level(mg_level))
+ if (cell->is_locally_owned_on_level())
+ {
+ indices.resize(cell->get_fe().n_dofs_per_cell());
+ cell->get_mg_dof_indices(indices);
+
+ task(cell);
+ }
+ }
+ };
+
+ runner([&](const auto &cell) {
+ const types::fe_index fe_index = cell->active_fe_index();
+ for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell(); ++i)
+ if (locally_owned_dofs.is_element(indices[i]))
+ dofs_by_component[locally_owned_dofs.index_within_set(indices[i])] =
+ local_component_association[fe_index][i];
+ });
}
}
-
- template <int dim, int spacedim>
- void
- extract_constant_modes(const DoFHandler<dim, spacedim> &dof_handler,
- const ComponentMask &component_mask,
- std::vector<std::vector<bool>> &constant_modes)
+ namespace internal
{
- // If there are no locally owned DoFs, return with an empty
- // constant_modes object:
- if (dof_handler.n_locally_owned_dofs() == 0)
+ template <int dim, int spacedim>
+ void
+ extract_constant_modes(const DoFHandler<dim, spacedim> &dof_handler,
+ const ComponentMask &component_mask,
+ const unsigned int mg_level,
+ std::vector<std::vector<bool>> &constant_modes)
+ {
+ const auto &locally_owned_dofs =
+ (mg_level == numbers::invalid_unsigned_int) ?
+ dof_handler.locally_owned_dofs() :
+ dof_handler.locally_owned_mg_dofs(mg_level);
+
+ // If there are no locally owned DoFs, return with an empty
+ // constant_modes object:
+ if (locally_owned_dofs.n_elements() == 0)
+ {
+ constant_modes = std::vector<std::vector<bool>>(0);
+ return;
+ }
+
+ const unsigned int n_components = dof_handler.get_fe(0).n_components();
+ Assert(component_mask.represents_n_components(n_components),
+ ExcDimensionMismatch(n_components, component_mask.size()));
+
+ std::vector<unsigned char> dofs_by_component(
+ locally_owned_dofs.n_elements());
+ internal::get_component_association(dof_handler,
+ component_mask,
+ dofs_by_component,
+ mg_level);
+ unsigned int n_selected_dofs = 0;
+ for (unsigned int i = 0; i < n_components; ++i)
+ if (component_mask[i] == true)
+ n_selected_dofs +=
+ std::count(dofs_by_component.begin(), dofs_by_component.end(), i);
+
+ // Find local numbering within the selected components
+ std::vector<unsigned int> component_numbering(
+ locally_owned_dofs.n_elements(), numbers::invalid_unsigned_int);
+ for (unsigned int i = 0, count = 0; i < locally_owned_dofs.n_elements();
+ ++i)
+ if (component_mask[dofs_by_component[i]])
+ component_numbering[i] = count++;
+
+ // get the element constant modes and find a translation table between
+ // index in the constant modes and the components.
+ //
+ // TODO: We might be able to extend this also for elements which do not
+ // have the same constant modes, but that is messy...
+ const dealii::hp::FECollection<dim, spacedim> &fe_collection =
+ dof_handler.get_fe_collection();
+ std::vector<Table<2, bool>> element_constant_modes;
+ std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
+ constant_mode_to_component_translation(n_components);
{
- constant_modes = std::vector<std::vector<bool>>(0);
- return;
- }
+ unsigned int n_constant_modes = 0;
+ int first_non_empty_constant_mode = -1;
+ for (unsigned int f = 0; f < fe_collection.size(); ++f)
+ {
+ std::pair<Table<2, bool>, std::vector<unsigned int>> data =
+ fe_collection[f].get_constant_modes();
- const unsigned int n_components = dof_handler.get_fe(0).n_components();
- Assert(component_mask.represents_n_components(n_components),
- ExcDimensionMismatch(n_components, component_mask.size()));
+ // Store the index of the current element if it is the first that
+ // has non-empty constant modes.
+ if (first_non_empty_constant_mode < 0 && data.first.n_rows() > 0)
+ {
+ first_non_empty_constant_mode = f;
+ // This is the first non-empty constant mode, so we figure out
+ // the translation between index in the constant modes and the
+ // components
+ for (unsigned int i = 0; i < data.second.size(); ++i)
+ if (component_mask[data.second[i]])
+ constant_mode_to_component_translation[data.second[i]]
+ .emplace_back(n_constant_modes++, i);
+ }
- std::vector<unsigned char> dofs_by_component(
- dof_handler.n_locally_owned_dofs());
- internal::get_component_association(dof_handler,
- component_mask,
- dofs_by_component);
- unsigned int n_selected_dofs = 0;
- for (unsigned int i = 0; i < n_components; ++i)
- if (component_mask[i] == true)
- n_selected_dofs +=
- std::count(dofs_by_component.begin(), dofs_by_component.end(), i);
-
- // Find local numbering within the selected components
- const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
- std::vector<unsigned int> component_numbering(
- locally_owned_dofs.n_elements(), numbers::invalid_unsigned_int);
- for (unsigned int i = 0, count = 0; i < locally_owned_dofs.n_elements();
- ++i)
- if (component_mask[dofs_by_component[i]])
- component_numbering[i] = count++;
-
- // get the element constant modes and find a translation table between
- // index in the constant modes and the components.
- //
- // TODO: We might be able to extend this also for elements which do not
- // have the same constant modes, but that is messy...
- const dealii::hp::FECollection<dim, spacedim> &fe_collection =
- dof_handler.get_fe_collection();
- std::vector<Table<2, bool>> element_constant_modes;
- std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
- constant_mode_to_component_translation(n_components);
- {
- unsigned int n_constant_modes = 0;
- int first_non_empty_constant_mode = -1;
- for (unsigned int f = 0; f < fe_collection.size(); ++f)
- {
- std::pair<Table<2, bool>, std::vector<unsigned int>> data =
- fe_collection[f].get_constant_modes();
+ // Add the constant modes of this element to the list and assert
+ // that there are as many constant modes as for the other elements
+ // (or zero constant modes).
+ element_constant_modes.push_back(data.first);
+ Assert(element_constant_modes.back().n_rows() == 0 ||
+ element_constant_modes.back().n_rows() ==
+ element_constant_modes[first_non_empty_constant_mode]
+ .n_rows(),
+ ExcInternalError());
+ }
+ AssertIndexRange(first_non_empty_constant_mode, fe_collection.size());
- // Store the index of the current element if it is the first that has
- // non-empty constant modes.
- if (first_non_empty_constant_mode < 0 && data.first.n_rows() > 0)
- {
- first_non_empty_constant_mode = f;
- // This is the first non-empty constant mode, so we figure out the
- // translation between index in the constant modes and the
- // components
- for (unsigned int i = 0; i < data.second.size(); ++i)
- if (component_mask[data.second[i]])
- constant_mode_to_component_translation[data.second[i]]
- .emplace_back(n_constant_modes++, i);
- }
+ // Now we know the number of constant modes and resize the return vector
+ // accordingly
+ constant_modes.clear();
+ constant_modes.resize(n_constant_modes,
+ std::vector<bool>(n_selected_dofs, false));
+ }
- // Add the constant modes of this element to the list and assert that
- // there are as many constant modes as for the other elements (or zero
- // constant modes).
- element_constant_modes.push_back(data.first);
- Assert(
- element_constant_modes.back().n_rows() == 0 ||
- element_constant_modes.back().n_rows() ==
- element_constant_modes[first_non_empty_constant_mode].n_rows(),
- ExcInternalError());
- }
- AssertIndexRange(first_non_empty_constant_mode, fe_collection.size());
+ // Loop over all owned cells and ask the element for the constant modes
+ std::vector<types::global_dof_index> dof_indices;
- // Now we know the number of constant modes and resize the return vector
- // accordingly
- constant_modes.clear();
- constant_modes.resize(n_constant_modes,
- std::vector<bool>(n_selected_dofs, false));
- }
+ const auto runner = [&](const auto &task) {
+ if (mg_level == numbers::invalid_unsigned_int)
+ {
+ for (const auto &cell : dof_handler.active_cell_iterators() |
+ IteratorFilters::LocallyOwnedCell())
+ {
+ dof_indices.resize(cell->get_fe().n_dofs_per_cell());
+ cell->get_dof_indices(dof_indices);
- // Loop over all owned cells and ask the element for the constant modes
- std::vector<types::global_dof_index> dof_indices;
- for (const auto &cell : dof_handler.active_cell_iterators() |
- IteratorFilters::LocallyOwnedCell())
- {
- dof_indices.resize(cell->get_fe().n_dofs_per_cell());
- cell->get_dof_indices(dof_indices);
+ task(cell);
+ }
+ }
+ else
+ {
+ for (const auto &cell :
+ dof_handler.cell_iterators_on_level(mg_level))
+ if (cell->is_locally_owned_on_level())
+ {
+ dof_indices.resize(cell->get_fe().n_dofs_per_cell());
+ cell->get_mg_dof_indices(dof_indices);
+
+ task(cell);
+ }
+ }
+ };
+ runner([&](const auto &cell) {
for (unsigned int i = 0; i < dof_indices.size(); ++i)
if (locally_owned_dofs.is_element(dof_indices[i]))
{
element_constant_modes[cell->active_fe_index()](
indices.second, i);
}
- }
+ });
+ }
+ } // namespace internal
+
+
+
+ template <int dim, int spacedim>
+ void
+ extract_constant_modes(const DoFHandler<dim, spacedim> &dof_handler,
+ const ComponentMask &component_mask,
+ std::vector<std::vector<bool>> &constant_modes)
+ {
+ internal::extract_constant_modes(dof_handler,
+ component_mask,
+ numbers::invalid_unsigned_int,
+ constant_modes);
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ extract_level_constant_modes(const unsigned int level,
+ const DoFHandler<dim, spacedim> &dof_handler,
+ const ComponentMask &component_mask,
+ std::vector<std::vector<bool>> &constant_modes)
+ {
+ internal::extract_constant_modes(dof_handler,
+ component_mask,
+ level,
+ constant_modes);
}
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2010 - 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+
+// test DoFTools::extract_level_constant_modes (see also
+// extract_constant_modes_01.cc).
+
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+test()
+{
+ unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+ parallel::distributed::Triangulation<dim> tr(
+ MPI_COMM_WORLD,
+ Triangulation<dim>::none,
+ parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
+
+ std::vector<unsigned int> sub(2);
+ sub[0] = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+ sub[1] = 1;
+ GridGenerator::subdivided_hyper_rectangle(
+ static_cast<Triangulation<dim> &>(tr), sub, Point<2>(0, 0), Point<2>(1, 1));
+
+ FESystem<dim> fe(FE_Q<dim>(1), 2, FE_DGQ<dim>(0), 1);
+ DoFHandler<dim> dofh(tr);
+ dofh.distribute_dofs(fe);
+ dofh.distribute_mg_dofs();
+
+ if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+ deallog << "Total dofs=" << dofh.n_dofs() << std::endl;
+
+ // extract constant modes and print
+ if (myid == 0)
+ for (unsigned int c = 0; c < fe.n_components(); ++c)
+ {
+ ComponentMask mask(fe.n_components(), false);
+ mask.set(c, true);
+
+ std::vector<std::vector<bool>> constant_modes;
+ DoFTools::extract_level_constant_modes(0, dofh, mask, constant_modes);
+
+ for (unsigned int i = 0; i < constant_modes.size(); ++i)
+ {
+ for (unsigned int j = 0; j < constant_modes[i].size(); ++j)
+ deallog << (constant_modes[i][j] ? '1' : '0') << ' ';
+ deallog << std::endl;
+ }
+ }
+
+ // renumber dofs and do the same again
+ DoFRenumbering::component_wise(dofh, 0);
+ if (myid == 0)
+ for (unsigned int c = 0; c < fe.n_components(); ++c)
+ {
+ ComponentMask mask(fe.n_components(), false);
+ mask.set(c, true);
+
+ std::vector<std::vector<bool>> constant_modes;
+ DoFTools::extract_level_constant_modes(0, dofh, mask, constant_modes);
+
+ for (unsigned int i = 0; i < constant_modes.size(); ++i)
+ {
+ for (unsigned int j = 0; j < constant_modes[i].size(); ++j)
+ deallog << (constant_modes[i][j] ? '1' : '0') << ' ';
+ deallog << std::endl;
+ }
+ }
+}
+
+
+int
+main(int argc, char *argv[])
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+ unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+ deallog.push(Utilities::int_to_string(myid));
+
+ if (myid == 0)
+ {
+ initlog();
+
+ deallog.push("2d");
+ test<2>();
+ deallog.pop();
+ }
+ else
+ {
+ deallog.push("2d");
+ test<2>();
+ deallog.pop();
+ }
+}
projects / dealii.git / commitdiff