std::vector<AffineConstraints<double>> boundary_constraints(n_levels);
for (unsigned int level = 0; level < n_levels; ++level)
{
- const IndexSet dofset =
- DoFTools::extract_locally_relevant_level_dofs(dof_handler, level);
- boundary_constraints[level].reinit(dofset);
+ boundary_constraints[level].reinit(
+ dof_handler.locally_owned_mg_dofs(level),
+ DoFTools::extract_locally_relevant_level_dofs(dof_handler, level));
boundary_constraints[level].add_lines(
mg_constrained_dofs.get_refinement_edge_indices(level));
boundary_constraints[level].add_lines(
// of each matrix or vector will be stored where, then call the functions
// that actually set up the matrices, and at the end also resize the
// various vectors we keep around in this program.
- std::vector<IndexSet> stokes_partitioning, stokes_relevant_partitioning;
- IndexSet temperature_partitioning(n_T),
- temperature_relevant_partitioning(n_T);
- IndexSet stokes_relevant_set;
- {
- IndexSet stokes_index_set = stokes_dof_handler.locally_owned_dofs();
- stokes_partitioning.push_back(stokes_index_set.get_view(0, n_u));
- stokes_partitioning.push_back(stokes_index_set.get_view(n_u, n_u + n_p));
-
- stokes_relevant_set =
- DoFTools::extract_locally_relevant_dofs(stokes_dof_handler);
- stokes_relevant_partitioning.push_back(
- stokes_relevant_set.get_view(0, n_u));
- stokes_relevant_partitioning.push_back(
- stokes_relevant_set.get_view(n_u, n_u + n_p));
-
- temperature_partitioning = temperature_dof_handler.locally_owned_dofs();
- temperature_relevant_partitioning =
- DoFTools::extract_locally_relevant_dofs(temperature_dof_handler);
- }
+
+ const IndexSet &stokes_locally_owned_index_set =
+ stokes_dof_handler.locally_owned_dofs();
+ const IndexSet stokes_locally_relevant_set =
+ DoFTools::extract_locally_relevant_dofs(stokes_dof_handler);
+
+ std::vector<IndexSet> stokes_partitioning;
+ stokes_partitioning.push_back(
+ stokes_locally_owned_index_set.get_view(0, n_u));
+ stokes_partitioning.push_back(
+ stokes_locally_owned_index_set.get_view(n_u, n_u + n_p));
+
+ std::vector<IndexSet> stokes_relevant_partitioning;
+ stokes_relevant_partitioning.push_back(
+ stokes_locally_relevant_set.get_view(0, n_u));
+ stokes_relevant_partitioning.push_back(
+ stokes_locally_relevant_set.get_view(n_u, n_u + n_p));
+
+ const IndexSet temperature_partitioning =
+ temperature_dof_handler.locally_owned_dofs();
+ const IndexSet temperature_relevant_partitioning =
+ DoFTools::extract_locally_relevant_dofs(temperature_dof_handler);
// Following this, we can compute constraints for the solution vectors,
// including hanging node constraints and homogeneous and inhomogeneous
// objects.
{
stokes_constraints.clear();
- stokes_constraints.reinit(stokes_relevant_set);
+ stokes_constraints.reinit(stokes_locally_owned_index_set,
+ stokes_locally_relevant_set);
DoFTools::make_hanging_node_constraints(stokes_dof_handler,
stokes_constraints);
}
{
temperature_constraints.clear();
- temperature_constraints.reinit(temperature_relevant_partitioning);
+ temperature_constraints.reinit(temperature_partitioning,
+ temperature_relevant_partitioning);
DoFTools::make_hanging_node_constraints(temperature_dof_handler,
temperature_constraints);
{
Timer time;
setup_time = 0;
+ {
+ system_matrix.clear();
+ mg_matrices.clear_elements();
- system_matrix.clear();
- mg_matrices.clear_elements();
-
- dof_handler.distribute_dofs(fe);
- dof_handler.distribute_mg_dofs();
-
- pcout << "Number of degrees of freedom: " << dof_handler.n_dofs()
- << std::endl;
+ dof_handler.distribute_dofs(fe);
+ dof_handler.distribute_mg_dofs();
- const IndexSet locally_relevant_dofs =
- DoFTools::extract_locally_relevant_dofs(dof_handler);
+ pcout << "Number of degrees of freedom: " << dof_handler.n_dofs()
+ << std::endl;
- constraints.clear();
- constraints.reinit(locally_relevant_dofs);
- DoFTools::make_hanging_node_constraints(dof_handler, constraints);
- VectorTools::interpolate_boundary_values(
- mapping, dof_handler, 0, Functions::ZeroFunction<dim>(), constraints);
- constraints.close();
+ constraints.clear();
+ constraints.reinit(dof_handler.locally_owned_dofs(),
+ DoFTools::extract_locally_relevant_dofs(dof_handler));
+ DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+ VectorTools::interpolate_boundary_values(
+ mapping, dof_handler, 0, Functions::ZeroFunction<dim>(), constraints);
+ constraints.close();
+ }
setup_time += time.wall_time();
time_details << "Distribute DoFs & B.C. (CPU/wall) " << time.cpu_time()
<< "s/" << time.wall_time() << 's' << std::endl;
time.restart();
-
{
- 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,
- dof_handler,
- constraints,
- QGauss<1>(fe.degree + 1),
- additional_data);
- system_matrix.initialize(system_mf_storage);
- }
-
- system_matrix.evaluate_coefficient(Coefficient<dim>());
+ {
+ 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,
+ dof_handler,
+ constraints,
+ QGauss<1>(fe.degree + 1),
+ additional_data);
+ system_matrix.initialize(system_mf_storage);
+ }
- system_matrix.initialize_dof_vector(solution);
- system_matrix.initialize_dof_vector(system_rhs);
+ system_matrix.evaluate_coefficient(Coefficient<dim>());
+ system_matrix.initialize_dof_vector(solution);
+ system_matrix.initialize_dof_vector(system_rhs);
+ }
setup_time += time.wall_time();
time_details << "Setup matrix-free system (CPU/wall) " << time.cpu_time()
<< "s/" << time.wall_time() << 's' << std::endl;
// closely the construction of the system matrix on the original mesh,
// except the slight difference in naming when accessing information on
// the levels rather than the active cells.
- const unsigned int nlevels = triangulation.n_global_levels();
- mg_matrices.resize(0, nlevels - 1);
-
- const std::set<types::boundary_id> dirichlet_boundary_ids = {0};
- mg_constrained_dofs.initialize(dof_handler);
- mg_constrained_dofs.make_zero_boundary_constraints(dof_handler,
- dirichlet_boundary_ids);
-
- for (unsigned int level = 0; level < nlevels; ++level)
- {
- const IndexSet relevant_dofs =
- DoFTools::extract_locally_relevant_level_dofs(dof_handler, level);
- AffineConstraints<double> 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.mg_level = level;
- std::shared_ptr<MatrixFree<dim, float>> mg_mf_storage_level(
- new MatrixFree<dim, float>());
- mg_mf_storage_level->reinit(mapping,
- 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].evaluate_coefficient(Coefficient<dim>());
- }
+ {
+ const unsigned int nlevels = triangulation.n_global_levels();
+ mg_matrices.resize(0, nlevels - 1);
+
+ const std::set<types::boundary_id> dirichlet_boundary_ids = {0};
+ mg_constrained_dofs.initialize(dof_handler);
+ mg_constrained_dofs.make_zero_boundary_constraints(
+ dof_handler, dirichlet_boundary_ids);
+
+ for (unsigned int level = 0; level < nlevels; ++level)
+ {
+ AffineConstraints<double> level_constraints(
+ dof_handler.locally_owned_mg_dofs(level),
+ DoFTools::extract_locally_relevant_level_dofs(dof_handler, level));
+ 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.mg_level = level;
+ std::shared_ptr<MatrixFree<dim, float>> mg_mf_storage_level(
+ new MatrixFree<dim, float>());
+ mg_mf_storage_level->reinit(mapping,
+ 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].evaluate_coefficient(Coefficient<dim>());
+ }
+ }
setup_time += time.wall_time();
time_details << "Setup matrix-free levels (CPU/wall) " << time.cpu_time()
<< "s/" << time.wall_time() << 's' << std::endl;
// expect any information to store. In our case, as explained in the
// @ref distributed module, the degrees of freedom we need to care about on
// each processor are the locally relevant ones, so we pass this to the
- // AffineConstraints::reinit function. As a side note, if you forget to
- // pass this argument, the AffineConstraints class will allocate an array
+ // AffineConstraints::reinit() function as a second argument. A further
+ // optimization, AffineConstraint can avoid certain operations if you also
+ // provide it with the set of locally owned degrees of freedom -- the
+ // first argument to AffineConstraints::reinit().
+ //
+ // (What would happen if we didn't pass this information to
+ // AffineConstraints, for example if we called the argument-less version of
+ // AffineConstraints::reinit() typically used in non-parallel codes? In that
+ // case, the AffineConstraints class will allocate an array
// with length equal to the largest DoF index it has seen so far. For
- // processors with high MPI process number, this may be very large --
+ // processors with large numbers of MPI processes, this may be very large --
// maybe on the order of billions. The program would then allocate more
// memory than for likely all other operations combined for this single
- // array.
+ // array. Fortunately, recent versions of deal.II would trigger an assertion
+ // that tells you that this is considered a bug.)
constraints.clear();
- constraints.reinit(locally_relevant_dofs);
+ constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
VectorTools::interpolate_boundary_values(dof_handler,
0,
/* setup hanging nodes and Dirichlet constraints */
{
TimerOutput::Scope t(computing_timer, "Setup: constraints");
- constraints_hanging_nodes.reinit(locally_relevant_dofs);
+ constraints_hanging_nodes.reinit(locally_owned_dofs,
+ locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler,
constraints_hanging_nodes);
constraints_hanging_nodes.close();
template <int dim>
void PlasticityContactProblem<dim>::compute_dirichlet_constraints()
{
- constraints_dirichlet_and_hanging_nodes.reinit(locally_relevant_dofs);
+ constraints_dirichlet_and_hanging_nodes.reinit(locally_owned_dofs,
+ locally_relevant_dofs);
constraints_dirichlet_and_hanging_nodes.merge(constraints_hanging_nodes);
if (base_mesh == "box")
diag_mass_matrix_vector_relevant = diag_mass_matrix_vector;
- all_constraints.reinit(locally_relevant_dofs);
+ all_constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
active_set.clear();
// The second part is a loop over all cells in which we look at each
{
owned_partitioning.clear();
- IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs();
+ const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
owned_partitioning.push_back(locally_owned_dofs.get_view(0, n_u));
owned_partitioning.push_back(locally_owned_dofs.get_view(n_u, n_u + n_p));
locally_relevant_dofs.get_view(n_u, n_u + n_p));
constraints.clear();
- constraints.reinit(locally_relevant_dofs);
+ constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
const FEValuesExtractors::Vector velocities(0);
// We generate hanging node constraints for ensuring continuity of the
// solution. As in step-40, we need to equip the constraint matrix with
- // the IndexSet of locally relevant degrees of freedom to avoid it to
- // consume too much memory for big problems. Next, the <code> MatrixFree
- // </code> object for the problem is set up. Note that we specify a
+ // the IndexSet of locally active and locally relevant degrees of freedom
+ // to avoid it consuming too much memory for big problems. Next, the
+ // MatrixFree object for the problem is set up. Note that we specify a
// particular scheme for shared-memory parallelization (hence one would
// use multithreading for intra-node parallelism and not MPI; we here
// choose the standard option — if we wanted to disable shared
locally_relevant_dofs =
DoFTools::extract_locally_relevant_dofs(dof_handler);
constraints.clear();
- constraints.reinit(locally_relevant_dofs);
+ constraints.reinit(dof_handler.locally_owned_dofs(), locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
constraints.close();
solution.reinit(locally_owned_dofs, mpi_communicator);
right_hand_side.reinit(locally_owned_dofs, mpi_communicator);
- constraints.reinit(locally_relevant_dofs);
+ constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
VectorTools::interpolate_boundary_values(
for (unsigned int level = 0; level < n_levels; ++level)
{
- const IndexSet relevant_dofs =
+ AffineConstraints<double> level_constraints(
+ dof_handler.locally_owned_mg_dofs(level),
DoFTools::extract_locally_relevant_level_dofs(dof_handler,
- level);
- AffineConstraints<double> level_constraints;
- level_constraints.reinit(relevant_dofs);
+ level));
level_constraints.add_lines(
mg_constrained_dofs.get_boundary_indices(level));
level_constraints.close();
triangulation.n_global_levels());
for (unsigned int level = 0; level < triangulation.n_global_levels(); ++level)
{
- const IndexSet dof_set =
- DoFTools::extract_locally_relevant_level_dofs(dof_handler, level);
- boundary_constraints[level].reinit(dof_set);
+ boundary_constraints[level].reinit(
+ dof_handler.locally_owned_mg_dofs(level),
+ DoFTools::extract_locally_relevant_level_dofs(dof_handler, level));
boundary_constraints[level].add_lines(
mg_constrained_dofs.get_refinement_edge_indices(level));
boundary_constraints[level].add_lines(
// We split up the IndexSet for locally owned and locally relevant DoFs
// into two IndexSets based on how we want to create the block matrices
// and vectors.
+ const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
owned_partitioning.resize(2);
- owned_partitioning[0] = dof_handler.locally_owned_dofs().get_view(0, n_u);
- owned_partitioning[1] =
- dof_handler.locally_owned_dofs().get_view(n_u, n_u + n_p);
+ owned_partitioning[0] = locally_owned_dofs.get_view(0, n_u);
+ owned_partitioning[1] = locally_owned_dofs.get_view(n_u, n_u + n_p);
const IndexSet locally_relevant_dofs =
DoFTools::extract_locally_relevant_dofs(dof_handler);
// to put this function call in, in case adaptive refinement gets
// introduced later.
{
- constraints.reinit(locally_relevant_dofs);
+ constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
const FEValuesExtractors::Vector velocities(0);
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
for (unsigned int level = 0; level < triangulation.n_global_levels();
++level)
{
- const IndexSet locally_owned_level_dof_indices =
- DoFTools::extract_locally_relevant_level_dofs(dof_handler, level);
- boundary_constraints[level].reinit(locally_owned_level_dof_indices);
+ boundary_constraints[level].reinit(
+ dof_handler.locally_owned_mg_dofs(level),
+ DoFTools::extract_locally_relevant_level_dofs(dof_handler, level));
boundary_constraints[level].add_lines(
mg_constrained_dofs.get_refinement_edge_indices(level));
boundary_constraints[level].add_lines(
system_rhs_dev.reinit(locally_owned_dofs, mpi_communicator);
constraints.clear();
- constraints.reinit(locally_relevant_dofs);
+ constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
VectorTools::interpolate_boundary_values(dof_handler,
0,
dof_handler.distribute_dofs(fe);
dof_handler.distribute_mg_dofs();
- const IndexSet locally_relevant_dofs =
- DoFTools::extract_locally_relevant_dofs(dof_handler);
-
constraints.clear();
- constraints.reinit(locally_relevant_dofs);
+ constraints.reinit(dof_handler.locally_owned_dofs(),
+ DoFTools::extract_locally_relevant_dofs(dof_handler));
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
VectorTools::interpolate_boundary_values(dof_handler,
0,
for (unsigned int level = 0; level < nlevels; ++level)
{
- const IndexSet relevant_dofs =
- DoFTools::extract_locally_relevant_level_dofs(dof_handler, level);
-
- AffineConstraints<double> level_constraints;
- level_constraints.reinit(relevant_dofs);
+ AffineConstraints<double> level_constraints(
+ dof_handler.locally_owned_mg_dofs(level),
+ DoFTools::extract_locally_relevant_level_dofs(dof_handler, level));
level_constraints.add_lines(
mg_constrained_dofs.get_boundary_indices(level));
level_constraints.close();
fluid_relevant_dofs[1] = locally_relevant_dofs.get_view(n_u, n_u + n_p);
{
- constraints.reinit(locally_relevant_dofs);
+ constraints.reinit(fluid_dh.locally_owned_dofs(), locally_relevant_dofs);
const FEValuesExtractors::Vector velocities(0);
DoFTools::make_hanging_node_constraints(fluid_dh, constraints);
constraints.close();
}
- auto locally_owned_dofs_per_processor =
+ const auto locally_owned_dofs_per_processor =
Utilities::MPI::all_gather(mpi_communicator,
fluid_dh.locally_owned_dofs());
{
// operator is applied to a vector with only the Dirichlet values set. The
// result is the negative right-hand-side vector.
{
- AffineConstraints<number> constraints_without_dbc;
-
- const IndexSet locally_relevant_dofs =
- DoFTools::extract_locally_relevant_dofs(dof_handler);
- constraints_without_dbc.reinit(locally_relevant_dofs);
+ AffineConstraints<number> constraints_without_dbc(
+ dof_handler.locally_owned_dofs(),
+ DoFTools::extract_locally_relevant_dofs(dof_handler));
DoFTools::make_hanging_node_constraints(dof_handler,
constraints_without_dbc);
const auto &dof_handler = dof_handlers[level];
auto &constraint = constraints[level];
- const IndexSet locally_relevant_dofs =
- DoFTools::extract_locally_relevant_dofs(dof_handler);
- constraint.reinit(locally_relevant_dofs);
+ constraint.reinit(dof_handler.locally_owned_dofs(),
+ DoFTools::extract_locally_relevant_dofs(dof_handler));
DoFTools::make_hanging_node_constraints(dof_handler, constraint);
VectorTools::interpolate_boundary_values(mapping_collection,
system_rhs.reinit(locally_owned_dofs, mpi_communicator);
constraints.clear();
- constraints.reinit(locally_relevant_dofs);
+ constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
VectorTools::interpolate_boundary_values(
mapping_collection, dof_handler, 0, Solution<dim>(), constraints);
projects / dealii.git / commitdiff