From 9cdebe50c714465b144f0731d12b7d43edd22f47 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 17 Nov 2021 14:09:54 -0700 Subject: [PATCH] Use iterator filters in our code base. --- .../mg_transfer_global_coarsening.templates.h | 7 +- .../vector_tools_interpolate.templates.h | 32 +- .../vector_tools_mean_value.templates.h | 24 +- source/base/mpi_remote_point_evaluation.cc | 7 +- .../repartitioning_policy_tools.cc | 11 +- source/dofs/dof_handler.cc | 11 +- source/dofs/dof_tools.cc | 70 ++--- source/dofs/dof_tools_constraints.cc | 51 ++-- source/grid/grid_tools.cc | 14 +- source/hp/refinement.cc | 37 +-- source/non_matching/coupling.cc | 201 ++++++------ source/numerics/smoothness_estimator.cc | 286 +++++++++--------- 12 files changed, 380 insertions(+), 371 deletions(-) diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h index 97d7cff744..de67003379 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h @@ -33,6 +33,7 @@ #include #include +#include #include #include @@ -2827,9 +2828,9 @@ MGTwoLevelTransfer>::reinit( if (cell->is_locally_owned()) is_locally_owned_fine.add_index(cell_id_translator.translate(cell)); - for (const auto &cell : dof_handler_coarse.active_cell_iterators()) - if (cell->is_locally_owned()) - is_locally_owned_coarse.add_index(cell_id_translator.translate(cell)); + for (const auto &cell : dof_handler_coarse.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + is_locally_owned_coarse.add_index(cell_id_translator.translate(cell)); const MPI_Comm communicator = dof_handler_fine.get_communicator(); diff --git a/include/deal.II/numerics/vector_tools_interpolate.templates.h b/include/deal.II/numerics/vector_tools_interpolate.templates.h index c4900770d6..6636dd18db 100644 --- a/include/deal.II/numerics/vector_tools_interpolate.templates.h +++ b/include/deal.II/numerics/vector_tools_interpolate.templates.h @@ -20,6 +20,7 @@ #include #include +#include #include #include @@ -673,22 +674,21 @@ namespace VectorTools Assert(fe.is_primitive(), ExcMessage("FE is not Primitive! This won't work.")); - for (const auto &cell : dh.active_cell_iterators()) - if (cell->is_locally_owned()) - { - fe_v.reinit(cell); - cell->get_dof_indices(dofs); - const std::vector> &points = - fe_v.get_quadrature_points(); - for (unsigned int q = 0; q < points.size(); ++q) - { - const unsigned int comp = - fe.system_to_component_index(q).first; - if (fe_mask[comp]) - ::dealii::internal::ElementAccess::set( - points[q][fe_to_real[comp]], dofs[q], vector); - } - } + for (const auto &cell : + dh.active_cell_iterators() | IteratorFilters::LocallyOwnedCell()) + { + fe_v.reinit(cell); + cell->get_dof_indices(dofs); + const std::vector> &points = + fe_v.get_quadrature_points(); + for (unsigned int q = 0; q < points.size(); ++q) + { + const unsigned int comp = fe.system_to_component_index(q).first; + if (fe_mask[comp]) + ::dealii::internal::ElementAccess::set( + points[q][fe_to_real[comp]], dofs[q], vector); + } + } } else { diff --git a/include/deal.II/numerics/vector_tools_mean_value.templates.h b/include/deal.II/numerics/vector_tools_mean_value.templates.h index 28636773cc..91804a5d35 100644 --- a/include/deal.II/numerics/vector_tools_mean_value.templates.h +++ b/include/deal.II/numerics/vector_tools_mean_value.templates.h @@ -19,6 +19,8 @@ #include +#include + #include #include #include @@ -143,17 +145,17 @@ namespace VectorTools Number mean = Number(); typename numbers::NumberTraits::real_type area = 0.; // Compute mean value - for (const auto &cell : dof.active_cell_iterators()) - if (cell->is_locally_owned()) - { - fe.reinit(cell); - fe.get_function_values(v, values); - for (unsigned int k = 0; k < quadrature.size(); ++k) - { - mean += fe.JxW(k) * values[k](component); - area += fe.JxW(k); - } - } + for (const auto &cell : + dof.active_cell_iterators() | IteratorFilters::LocallyOwnedCell()) + { + fe.reinit(cell); + fe.get_function_values(v, values); + for (unsigned int k = 0; k < quadrature.size(); ++k) + { + mean += fe.JxW(k) * values[k](component); + area += fe.JxW(k); + } + } #ifdef DEAL_II_WITH_MPI // if this was a distributed DoFHandler, we need to do the reduction diff --git a/source/base/mpi_remote_point_evaluation.cc b/source/base/mpi_remote_point_evaluation.cc index e146270219..5caf30a611 100644 --- a/source/base/mpi_remote_point_evaluation.cc +++ b/source/base/mpi_remote_point_evaluation.cc @@ -24,6 +24,7 @@ #include +#include #include #include #include @@ -80,9 +81,9 @@ namespace Utilities this->mapping = &mapping; std::vector> local_boxes; - for (const auto &cell : tria.active_cell_iterators()) - if (cell->is_locally_owned()) - local_boxes.push_back(mapping.get_bounding_box(cell)); + for (const auto &cell : + tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell()) + local_boxes.push_back(mapping.get_bounding_box(cell)); // create r-tree of bounding boxes const auto local_tree = pack_rtree(local_boxes); diff --git a/source/distributed/repartitioning_policy_tools.cc b/source/distributed/repartitioning_policy_tools.cc index be19473898..ada719d1c2 100644 --- a/source/distributed/repartitioning_policy_tools.cc +++ b/source/distributed/repartitioning_policy_tools.cc @@ -21,6 +21,7 @@ #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -171,11 +172,11 @@ namespace RepartitioningPolicyTools LinearAlgebra::distributed::Vector partition( tria->global_active_cell_index_partitioner().lock()); - for (const auto &cell : tria_coarse_in.active_cell_iterators()) - if (cell->is_locally_owned()) - partition[cell->global_active_cell_index()] = - owning_ranks_of_coarse_cells[is_coarse.index_within_set( - cell_id_translator.translate(cell))]; + for (const auto &cell : tria_coarse_in.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + partition[cell->global_active_cell_index()] = + owning_ranks_of_coarse_cells[is_coarse.index_within_set( + cell_id_translator.translate(cell))]; return partition; } diff --git a/source/dofs/dof_handler.cc b/source/dofs/dof_handler.cc index 236896c10c..33f7ccd18d 100644 --- a/source/dofs/dof_handler.cc +++ b/source/dofs/dof_handler.cc @@ -27,6 +27,7 @@ #include #include +#include #include #include #include @@ -1793,11 +1794,11 @@ namespace internal { std::vector future_fe_indices( tr->n_active_cells(), 0u); - for (const auto &cell : dof_handler.active_cell_iterators()) - if (cell->is_locally_owned()) - future_fe_indices[cell->active_cell_index()] = - dof_handler - .hp_cell_future_fe_indices[cell->level()][cell->index()]; + for (const auto &cell : dof_handler.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + future_fe_indices[cell->active_cell_index()] = + dof_handler + .hp_cell_future_fe_indices[cell->level()][cell->index()]; Utilities::MPI::sum(future_fe_indices, tr->get_communicator(), diff --git a/source/dofs/dof_tools.cc b/source/dofs/dof_tools.cc index 7b6469bc52..3ad2a1025e 100644 --- a/source/dofs/dof_tools.cc +++ b/source/dofs/dof_tools.cc @@ -224,18 +224,18 @@ namespace DoFTools // then loop over all cells and do the work std::vector indices; - for (const auto &c : dof.active_cell_iterators()) - if (c->is_locally_owned()) - { - const unsigned int 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]; - } + for (const auto &c : + dof.active_cell_iterators() | IteratorFilters::LocallyOwnedCell()) + { + const unsigned int 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]; + } } @@ -1272,27 +1272,27 @@ namespace DoFTools // Loop over all owned cells and ask the element for the constant modes std::vector dof_indices; - for (const auto &cell : dof_handler.active_cell_iterators()) - if (cell->is_locally_owned()) - { - dof_indices.resize(cell->get_fe().n_dofs_per_cell()); - cell->get_dof_indices(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); - for (unsigned int i = 0; i < dof_indices.size(); ++i) - if (locally_owned_dofs.is_element(dof_indices[i])) - { - const unsigned int loc_index = - locally_owned_dofs.index_within_set(dof_indices[i]); - const unsigned int comp = dofs_by_component[loc_index]; - if (component_mask[comp]) - for (auto &indices : - constant_mode_to_component_translation[comp]) - constant_modes[indices - .first][component_numbering[loc_index]] = - element_constant_modes[cell->active_fe_index()]( - indices.second, i); - } - } + for (unsigned int i = 0; i < dof_indices.size(); ++i) + if (locally_owned_dofs.is_element(dof_indices[i])) + { + const unsigned int loc_index = + locally_owned_dofs.index_within_set(dof_indices[i]); + const unsigned int comp = dofs_by_component[loc_index]; + if (component_mask[comp]) + for (auto &indices : + constant_mode_to_component_translation[comp]) + constant_modes[indices + .first][component_numbering[loc_index]] = + element_constant_modes[cell->active_fe_index()]( + indices.second, i); + } + } } @@ -1523,9 +1523,9 @@ namespace DoFTools } else { - for (const auto &cell : dof_handler.active_cell_iterators()) - if (cell->is_locally_owned()) - cell_owners[cell->active_cell_index()] = cell->subdomain_id(); + for (const auto &cell : dof_handler.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + cell_owners[cell->active_cell_index()] = cell->subdomain_id(); } // preset all values by an invalid value diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc index d285de7841..874d1e5564 100644 --- a/source/dofs/dof_tools_constraints.cc +++ b/source/dofs/dof_tools_constraints.cc @@ -26,6 +26,7 @@ #include #include +#include #include #include #include @@ -3062,15 +3063,15 @@ namespace DoFTools std::vector local_dof_indices( fine_fe.n_dofs_per_cell()); - for (const auto &cell : fine_grid.active_cell_iterators()) - if (cell->is_locally_owned()) - { - cell->get_dof_indices(local_dof_indices); - for (unsigned int i = 0; i < fine_fe.n_dofs_per_cell(); ++i) - if (fine_fe.system_to_component_index(i).first == - fine_component) - dof_is_interesting[local_dof_indices[i]] = true; - } + for (const auto &cell : fine_grid.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + { + cell->get_dof_indices(local_dof_indices); + for (unsigned int i = 0; i < fine_fe.n_dofs_per_cell(); ++i) + if (fine_fe.system_to_component_index(i).first == + fine_component) + dof_is_interesting[local_dof_indices[i]] = true; + } n_parameters_on_fine_grid = std::count(dof_is_interesting.begin(), dof_is_interesting.end(), @@ -3089,22 +3090,22 @@ namespace DoFTools std::vector local_dof_indices( fine_fe.n_dofs_per_cell()); unsigned int next_free_index = 0; - for (const auto &cell : fine_grid.active_cell_iterators()) - if (cell->is_locally_owned()) - { - cell->get_dof_indices(local_dof_indices); - for (unsigned int i = 0; i < fine_fe.n_dofs_per_cell(); ++i) - // if this DoF is a parameter dof and has not yet been - // numbered, then do so - if ((fine_fe.system_to_component_index(i).first == - fine_component) && - (weight_mapping[local_dof_indices[i]] == - numbers::invalid_dof_index)) - { - weight_mapping[local_dof_indices[i]] = next_free_index; - ++next_free_index; - } - } + for (const auto &cell : fine_grid.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + { + cell->get_dof_indices(local_dof_indices); + for (unsigned int i = 0; i < fine_fe.n_dofs_per_cell(); ++i) + // if this DoF is a parameter dof and has not yet been + // numbered, then do so + if ((fine_fe.system_to_component_index(i).first == + fine_component) && + (weight_mapping[local_dof_indices[i]] == + numbers::invalid_dof_index)) + { + weight_mapping[local_dof_indices[i]] = next_free_index; + ++next_free_index; + } + } Assert(next_free_index == n_parameters_on_fine_grid, ExcInternalError()); diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index ba6e7f2a0b..7989543010 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -6480,10 +6480,10 @@ namespace GridTools // 2) collect vertices belonging to local cells std::vector vertex_of_own_cell(tria.n_vertices(), false); - for (const auto &cell : tria.active_cell_iterators()) - if (cell->is_locally_owned()) - for (const unsigned int v : cell->vertex_indices()) - vertex_of_own_cell[cell->vertex_index(v)] = true; + for (const auto &cell : + tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell()) + for (const unsigned int v : cell->vertex_indices()) + vertex_of_own_cell[cell->vertex_index(v)] = true; // 3) for each vertex belonging to a locally owned cell all ghost // neighbors (including the periodic own) @@ -6574,9 +6574,9 @@ namespace GridTools std::vector> & vertices, std::vector> &cells) const { - for (const auto &cell : background_dof_handler.active_cell_iterators()) - if (cell->is_locally_owned()) - process_cell(cell, ls_vector, iso_level, vertices, cells); + for (const auto &cell : background_dof_handler.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + process_cell(cell, ls_vector, iso_level, vertices, cells); } diff --git a/source/hp/refinement.cc b/source/hp/refinement.cc index 4b28f63aa7..57917f61a2 100644 --- a/source/hp/refinement.cc +++ b/source/hp/refinement.cc @@ -25,6 +25,7 @@ #include #include +#include #include #include @@ -869,10 +870,10 @@ namespace hp dof_handler.get_triangulation().n_active_cells()); } - for (const auto &cell : dof_handler.active_cell_iterators()) - if (cell->is_locally_owned()) - future_levels[cell->global_active_cell_index()] = - hierarchy_level_for_fe_index[cell->future_fe_index()]; + for (const auto &cell : dof_handler.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + future_levels[cell->global_active_cell_index()] = + hierarchy_level_for_fe_index[cell->future_fe_index()]; // @@ -1041,22 +1042,22 @@ namespace hp while (levels_changed_in_cycle); // update future FE indices on locally owned cells - for (const auto &cell : dof_handler.active_cell_iterators()) - if (cell->is_locally_owned()) - { - const level_type cell_level = static_cast( - future_levels[cell->global_active_cell_index()]); + for (const auto &cell : dof_handler.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + { + const level_type cell_level = static_cast( + future_levels[cell->global_active_cell_index()]); - if (cell_level != invalid_level) - { - const unsigned int fe_index = - fe_index_for_hierarchy_level[cell_level]; + if (cell_level != invalid_level) + { + const unsigned int fe_index = + fe_index_for_hierarchy_level[cell_level]; - // only update if necessary - if (fe_index != cell->active_fe_index()) - cell->set_future_fe_index(fe_index); - } - } + // only update if necessary + if (fe_index != cell->active_fe_index()) + cell->set_future_fe_index(fe_index); + } + } return levels_changed; } diff --git a/source/non_matching/coupling.cc b/source/non_matching/coupling.cc index dcaa2b91fa..09aa5fa0d1 100644 --- a/source/non_matching/coupling.cc +++ b/source/non_matching/coupling.cc @@ -22,6 +22,7 @@ #include +#include #include #include @@ -692,31 +693,31 @@ namespace NonMatching typename Triangulation::active_cell_iterator>> intersection; - for (const auto &cell0 : dh0.active_cell_iterators()) - if (cell0->is_locally_owned()) - { - intersection.resize(0); - BoundingBox box0 = - cache0.get_mapping().get_bounding_box(cell0); - box0.extend(epsilon); - boost::geometry::index::query(tree1, - boost::geometry::index::intersects( - box0), - std::back_inserter(intersection)); - if (!intersection.empty()) - { - cell0->get_dof_indices(dofs0); - for (const auto &entry : intersection) - { - typename DoFHandler::cell_iterator cell1( - *entry.second, &dh1); - cell1->get_dof_indices(dofs1); - constraints0.add_entries_local_to_global(dofs0, - dofs1, - sparsity); - } - } - } + for (const auto &cell0 : + dh0.active_cell_iterators() | IteratorFilters::LocallyOwnedCell()) + { + intersection.resize(0); + BoundingBox box0 = + cache0.get_mapping().get_bounding_box(cell0); + box0.extend(epsilon); + boost::geometry::index::query(tree1, + boost::geometry::index::intersects( + box0), + std::back_inserter(intersection)); + if (!intersection.empty()) + { + cell0->get_dof_indices(dofs0); + for (const auto &entry : intersection) + { + typename DoFHandler::cell_iterator cell1( + *entry.second, &dh1); + cell1->get_dof_indices(dofs1); + constraints0.add_entries_local_to_global(dofs0, + dofs1, + sparsity); + } + } + } } else { @@ -728,31 +729,31 @@ namespace NonMatching typename Triangulation::active_cell_iterator>> intersection; - for (const auto &cell1 : dh1.active_cell_iterators()) - if (cell1->is_locally_owned()) - { - intersection.resize(0); - BoundingBox box1 = - cache1.get_mapping().get_bounding_box(cell1); - box1.extend(epsilon); - boost::geometry::index::query(tree0, - boost::geometry::index::intersects( - box1), - std::back_inserter(intersection)); - if (!intersection.empty()) - { - cell1->get_dof_indices(dofs1); - for (const auto &entry : intersection) - { - typename DoFHandler::cell_iterator cell0( - *entry.second, &dh0); - cell0->get_dof_indices(dofs0); - constraints0.add_entries_local_to_global(dofs0, - dofs1, - sparsity); - } - } - } + for (const auto &cell1 : + dh1.active_cell_iterators() | IteratorFilters::LocallyOwnedCell()) + { + intersection.resize(0); + BoundingBox box1 = + cache1.get_mapping().get_bounding_box(cell1); + box1.extend(epsilon); + boost::geometry::index::query(tree0, + boost::geometry::index::intersects( + box1), + std::back_inserter(intersection)); + if (!intersection.empty()) + { + cell1->get_dof_indices(dofs1); + for (const auto &entry : intersection) + { + typename DoFHandler::cell_iterator cell0( + *entry.second, &dh0); + cell0->get_dof_indices(dofs0); + constraints0.add_entries_local_to_global(dofs0, + dofs1, + sparsity); + } + } + } } } @@ -895,31 +896,31 @@ namespace NonMatching typename Triangulation::active_cell_iterator>> intersection; - for (const auto &cell0 : dh0.active_cell_iterators()) - if (cell0->is_locally_owned()) - { - intersection.resize(0); - BoundingBox box0 = - cache0.get_mapping().get_bounding_box(cell0); - box0.extend(epsilon); - boost::geometry::index::query(tree1, - boost::geometry::index::intersects( - box0), - std::back_inserter(intersection)); - if (!intersection.empty()) - { - cell0->get_dof_indices(dofs0); - fev0.reinit(cell0); - for (const auto &entry : intersection) - { - typename DoFHandler::cell_iterator cell1( - *entry.second, &dh1); - cell1->get_dof_indices(dofs1); - fev1.reinit(cell1); - assemble_one_pair(); - } - } - } + for (const auto &cell0 : + dh0.active_cell_iterators() | IteratorFilters::LocallyOwnedCell()) + { + intersection.resize(0); + BoundingBox box0 = + cache0.get_mapping().get_bounding_box(cell0); + box0.extend(epsilon); + boost::geometry::index::query(tree1, + boost::geometry::index::intersects( + box0), + std::back_inserter(intersection)); + if (!intersection.empty()) + { + cell0->get_dof_indices(dofs0); + fev0.reinit(cell0); + for (const auto &entry : intersection) + { + typename DoFHandler::cell_iterator cell1( + *entry.second, &dh1); + cell1->get_dof_indices(dofs1); + fev1.reinit(cell1); + assemble_one_pair(); + } + } + } } else { @@ -931,31 +932,31 @@ namespace NonMatching typename Triangulation::active_cell_iterator>> intersection; - for (const auto &cell1 : dh1.active_cell_iterators()) - if (cell1->is_locally_owned()) - { - intersection.resize(0); - BoundingBox box1 = - cache1.get_mapping().get_bounding_box(cell1); - box1.extend(epsilon); - boost::geometry::index::query(tree0, - boost::geometry::index::intersects( - box1), - std::back_inserter(intersection)); - if (!intersection.empty()) - { - cell1->get_dof_indices(dofs1); - fev1.reinit(cell1); - for (const auto &entry : intersection) - { - typename DoFHandler::cell_iterator cell0( - *entry.second, &dh0); - cell0->get_dof_indices(dofs0); - fev0.reinit(cell0); - assemble_one_pair(); - } - } - } + for (const auto &cell1 : + dh1.active_cell_iterators() | IteratorFilters::LocallyOwnedCell()) + { + intersection.resize(0); + BoundingBox box1 = + cache1.get_mapping().get_bounding_box(cell1); + box1.extend(epsilon); + boost::geometry::index::query(tree0, + boost::geometry::index::intersects( + box1), + std::back_inserter(intersection)); + if (!intersection.empty()) + { + cell1->get_dof_indices(dofs1); + fev1.reinit(cell1); + for (const auto &entry : intersection) + { + typename DoFHandler::cell_iterator cell0( + *entry.second, &dh0); + cell0->get_dof_indices(dofs0); + fev0.reinit(cell0); + assemble_one_pair(); + } + } + } } } diff --git a/source/numerics/smoothness_estimator.cc b/source/numerics/smoothness_estimator.cc index 5ec60c3c57..a92347cf4c 100644 --- a/source/numerics/smoothness_estimator.cc +++ b/source/numerics/smoothness_estimator.cc @@ -18,6 +18,8 @@ #include +#include + #include #include @@ -388,79 +390,77 @@ namespace SmoothnessEstimator Vector local_dof_values; std::vector ln_k; std::pair, std::vector> res; - for (const auto &cell : dof_handler.active_cell_iterators()) - if (cell->is_locally_owned()) - { - if (!only_flagged_cells || cell->refine_flag_set() || - cell->coarsen_flag_set()) - { - n_modes = fe_fourier.get_n_coefficients_per_direction( - cell->active_fe_index()); - resize(expansion_coefficients, n_modes); - - // Inside the loop, we first need to get the values of the local - // degrees of freedom and then need to compute the series - // expansion by multiplying this vector with the matrix ${\cal - // F}$ corresponding to this finite element. - local_dof_values.reinit(cell->get_fe().n_dofs_per_cell()); - cell->get_dof_values(solution, local_dof_values); - - fe_fourier.calculate(local_dof_values, - cell->active_fe_index(), - expansion_coefficients); - - // We fit our exponential decay of expansion coefficients to the - // provided regression_strategy on each possible value of |k|. - // To this end, we use FESeries::process_coefficients() to - // rework coefficients into the desired format. - res = FESeries::process_coefficients( - expansion_coefficients, - [n_modes](const TableIndices &indices) { - return index_norm_greater_than_zero_and_less_than_N_squared( - indices, n_modes); - }, - regression_strategy, - smallest_abs_coefficient); - - Assert(res.first.size() == res.second.size(), - ExcInternalError()); - - // Last, do the linear regression. - float regularity = std::numeric_limits::infinity(); - if (res.first.size() > 1) - { - // Prepare linear equation for the logarithmic least squares - // fit. - // - // First, calculate ln(|k|). - // - // For Fourier expansion, this translates to - // ln(2*pi*sqrt(predicate)) = ln(2*pi) + 0.5*ln(predicate). - // Since we are just interested in the slope of a linear - // regression later, we omit the ln(2*pi) factor. - ln_k.resize(res.first.size()); - for (unsigned int f = 0; f < res.first.size(); ++f) - ln_k[f] = - 0.5 * std::log(static_cast(res.first[f])); - - // Second, calculate ln(U_k). - for (auto &residual_element : res.second) - residual_element = std::log(residual_element); - - const std::pair fit = - FESeries::linear_regression(ln_k, res.second); - // Compute regularity s = mu - dim/2 - regularity = static_cast(-fit.first) - - ((dim > 1) ? (.5 * dim) : 0); - } - - // Store result in the vector of estimated values for each cell. - smoothness_indicators(cell->active_cell_index()) = regularity; - } - else - smoothness_indicators(cell->active_cell_index()) = - numbers::signaling_nan(); - } + for (const auto &cell : dof_handler.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + { + if (!only_flagged_cells || cell->refine_flag_set() || + cell->coarsen_flag_set()) + { + n_modes = fe_fourier.get_n_coefficients_per_direction( + cell->active_fe_index()); + resize(expansion_coefficients, n_modes); + + // Inside the loop, we first need to get the values of the local + // degrees of freedom and then need to compute the series + // expansion by multiplying this vector with the matrix ${\cal + // F}$ corresponding to this finite element. + local_dof_values.reinit(cell->get_fe().n_dofs_per_cell()); + cell->get_dof_values(solution, local_dof_values); + + fe_fourier.calculate(local_dof_values, + cell->active_fe_index(), + expansion_coefficients); + + // We fit our exponential decay of expansion coefficients to the + // provided regression_strategy on each possible value of |k|. + // To this end, we use FESeries::process_coefficients() to + // rework coefficients into the desired format. + res = FESeries::process_coefficients( + expansion_coefficients, + [n_modes](const TableIndices &indices) { + return index_norm_greater_than_zero_and_less_than_N_squared( + indices, n_modes); + }, + regression_strategy, + smallest_abs_coefficient); + + Assert(res.first.size() == res.second.size(), ExcInternalError()); + + // Last, do the linear regression. + float regularity = std::numeric_limits::infinity(); + if (res.first.size() > 1) + { + // Prepare linear equation for the logarithmic least squares + // fit. + // + // First, calculate ln(|k|). + // + // For Fourier expansion, this translates to + // ln(2*pi*sqrt(predicate)) = ln(2*pi) + 0.5*ln(predicate). + // Since we are just interested in the slope of a linear + // regression later, we omit the ln(2*pi) factor. + ln_k.resize(res.first.size()); + for (unsigned int f = 0; f < res.first.size(); ++f) + ln_k[f] = 0.5 * std::log(static_cast(res.first[f])); + + // Second, calculate ln(U_k). + for (auto &residual_element : res.second) + residual_element = std::log(residual_element); + + const std::pair fit = + FESeries::linear_regression(ln_k, res.second); + // Compute regularity s = mu - dim/2 + regularity = static_cast(-fit.first) - + ((dim > 1) ? (.5 * dim) : 0); + } + + // Store result in the vector of estimated values for each cell. + smoothness_indicators(cell->active_cell_index()) = regularity; + } + else + smoothness_indicators(cell->active_cell_index()) = + numbers::signaling_nan(); + } } @@ -498,78 +498,78 @@ namespace SmoothnessEstimator x.reserve(max_degree); y.reserve(max_degree); - for (const auto &cell : dof_handler.active_cell_iterators()) - if (cell->is_locally_owned()) - { - if (!only_flagged_cells || cell->refine_flag_set() || - cell->coarsen_flag_set()) - { - n_modes = fe_fourier.get_n_coefficients_per_direction( - cell->active_fe_index()); - resize(expansion_coefficients, n_modes); - - const unsigned int pe = cell->get_fe().degree; - Assert(pe > 0, ExcInternalError()); - - // since we use coefficients with indices [1,pe] in each - // direction, the number of coefficients we need to calculate is - // at least N=pe+1 - AssertIndexRange(pe, n_modes); - - local_dof_values.reinit(cell->get_fe().n_dofs_per_cell()); - cell->get_dof_values(solution, local_dof_values); - - fe_fourier.calculate(local_dof_values, - cell->active_fe_index(), - expansion_coefficients); - - // choose the smallest decay of coefficients in each direction, - // i.e. the maximum decay slope k_v as in exp(-k_v) - double k_v = std::numeric_limits::infinity(); - for (unsigned int d = 0; d < dim; ++d) - { - x.resize(0); - y.resize(0); - - // will use all non-zero coefficients allowed by the - // predicate function - // - // skip i=0 because of logarithm - for (unsigned int i = 1; i <= pe; ++i) - if (coefficients_predicate[i]) - { - TableIndices ind; - ind[d] = i; - const double coeff_abs = - std::abs(expansion_coefficients(ind)); - - if (coeff_abs > smallest_abs_coefficient) - { - x.push_back(std::log(i)); - y.push_back(std::log(coeff_abs)); - } - } - - // in case we don't have enough non-zero coefficient to fit, - // skip this direction - if (x.size() < 2) - continue; - - const std::pair fit = - FESeries::linear_regression(x, y); - - // decay corresponds to negative slope - // take the lesser negative slope along each direction - k_v = std::min(k_v, -fit.first); - } + for (const auto &cell : dof_handler.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + { + if (!only_flagged_cells || cell->refine_flag_set() || + cell->coarsen_flag_set()) + { + n_modes = fe_fourier.get_n_coefficients_per_direction( + cell->active_fe_index()); + resize(expansion_coefficients, n_modes); + + const unsigned int pe = cell->get_fe().degree; + Assert(pe > 0, ExcInternalError()); + + // since we use coefficients with indices [1,pe] in each + // direction, the number of coefficients we need to calculate is + // at least N=pe+1 + AssertIndexRange(pe, n_modes); + + local_dof_values.reinit(cell->get_fe().n_dofs_per_cell()); + cell->get_dof_values(solution, local_dof_values); + + fe_fourier.calculate(local_dof_values, + cell->active_fe_index(), + expansion_coefficients); + + // choose the smallest decay of coefficients in each direction, + // i.e. the maximum decay slope k_v as in exp(-k_v) + double k_v = std::numeric_limits::infinity(); + for (unsigned int d = 0; d < dim; ++d) + { + x.resize(0); + y.resize(0); + + // will use all non-zero coefficients allowed by the + // predicate function + // + // skip i=0 because of logarithm + for (unsigned int i = 1; i <= pe; ++i) + if (coefficients_predicate[i]) + { + TableIndices ind; + ind[d] = i; + const double coeff_abs = + std::abs(expansion_coefficients(ind)); + + if (coeff_abs > smallest_abs_coefficient) + { + x.push_back(std::log(i)); + y.push_back(std::log(coeff_abs)); + } + } + + // in case we don't have enough non-zero coefficient to fit, + // skip this direction + if (x.size() < 2) + continue; + + const std::pair fit = + FESeries::linear_regression(x, y); + + // decay corresponds to negative slope + // take the lesser negative slope along each direction + k_v = std::min(k_v, -fit.first); + } - smoothness_indicators(cell->active_cell_index()) = - static_cast(k_v); - } - else smoothness_indicators(cell->active_cell_index()) = - numbers::signaling_nan(); - } + static_cast(k_v); + } + else + smoothness_indicators(cell->active_cell_index()) = + numbers::signaling_nan(); + } } -- 2.39.5