#include <deal.II/fe/fe_tools.h>
#include <deal.II/grid/cell_id_translator.h>
+#include <deal.II/grid/filtered_iterator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria_description.h>
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();
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
+#include <deal.II/grid/filtered_iterator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/intergrid_map.h>
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<Point<spacedim>> &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<VectorType>::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<Point<spacedim>> &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<VectorType>::set(
+ points[q][fe_to_real[comp]], dofs[q], vector);
+ }
+ }
}
else
{
#include <deal.II/fe/fe_values.h>
+#include <deal.II/grid/filtered_iterator.h>
+
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/la_parallel_vector.h>
Number mean = Number();
typename numbers::NumberTraits<Number>::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
#include <deal.II/fe/mapping.h>
+#include <deal.II/grid/filtered_iterator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_tools_cache.h>
#include <deal.II/grid/tria.h>
this->mapping = &mapping;
std::vector<BoundingBox<spacedim>> 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);
#include <deal.II/distributed/tria_base.h>
#include <deal.II/grid/cell_id_translator.h>
+#include <deal.II/grid/filtered_iterator.h>
DEAL_II_NAMESPACE_OPEN
LinearAlgebra::distributed::Vector<double> 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;
}
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_handler_policy.h>
+#include <deal.II/grid/filtered_iterator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
{
std::vector<active_fe_index_type> 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(),
// then loop over all cells and do the work
std::vector<types::global_dof_index> 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];
+ }
}
// 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())
- 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);
+ }
+ }
}
}
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
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/grid/filtered_iterator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/intergrid_map.h>
#include <deal.II/grid/tria.h>
std::vector<types::global_dof_index> 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(),
std::vector<types::global_dof_index> 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());
// 2) collect vertices belonging to local cells
std::vector<bool> 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)
std::vector<Point<dim>> & vertices,
std::vector<CellData<dim - 1>> &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);
}
#include <deal.II/dofs/dof_accessor.templates.h>
#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/filtered_iterator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_tools.h>
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()];
//
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<level_type>(
- 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<level_type>(
+ 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;
}
#include <deal.II/fe/fe_values.h>
+#include <deal.II/grid/filtered_iterator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_tools_cache.h>
typename Triangulation<dim1, spacedim>::active_cell_iterator>>
intersection;
- for (const auto &cell0 : dh0.active_cell_iterators())
- if (cell0->is_locally_owned())
- {
- intersection.resize(0);
- BoundingBox<spacedim> 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<dim1, spacedim>::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<spacedim> 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<dim1, spacedim>::cell_iterator cell1(
+ *entry.second, &dh1);
+ cell1->get_dof_indices(dofs1);
+ constraints0.add_entries_local_to_global(dofs0,
+ dofs1,
+ sparsity);
+ }
+ }
+ }
}
else
{
typename Triangulation<dim0, spacedim>::active_cell_iterator>>
intersection;
- for (const auto &cell1 : dh1.active_cell_iterators())
- if (cell1->is_locally_owned())
- {
- intersection.resize(0);
- BoundingBox<spacedim> 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<dim0, spacedim>::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<spacedim> 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<dim0, spacedim>::cell_iterator cell0(
+ *entry.second, &dh0);
+ cell0->get_dof_indices(dofs0);
+ constraints0.add_entries_local_to_global(dofs0,
+ dofs1,
+ sparsity);
+ }
+ }
+ }
}
}
typename Triangulation<dim1, spacedim>::active_cell_iterator>>
intersection;
- for (const auto &cell0 : dh0.active_cell_iterators())
- if (cell0->is_locally_owned())
- {
- intersection.resize(0);
- BoundingBox<spacedim> 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<dim1, spacedim>::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<spacedim> 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<dim1, spacedim>::cell_iterator cell1(
+ *entry.second, &dh1);
+ cell1->get_dof_indices(dofs1);
+ fev1.reinit(cell1);
+ assemble_one_pair();
+ }
+ }
+ }
}
else
{
typename Triangulation<dim0, spacedim>::active_cell_iterator>>
intersection;
- for (const auto &cell1 : dh1.active_cell_iterators())
- if (cell1->is_locally_owned())
- {
- intersection.resize(0);
- BoundingBox<spacedim> 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<dim0, spacedim>::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<spacedim> 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<dim0, spacedim>::cell_iterator cell0(
+ *entry.second, &dh0);
+ cell0->get_dof_indices(dofs0);
+ fev0.reinit(cell0);
+ assemble_one_pair();
+ }
+ }
+ }
}
}
#include <deal.II/fe/fe_series.h>
+#include <deal.II/grid/filtered_iterator.h>
+
#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/q_collection.h>
Vector<number> local_dof_values;
std::vector<double> ln_k;
std::pair<std::vector<unsigned int>, std::vector<double>> 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<dim>(
- expansion_coefficients,
- [n_modes](const TableIndices<dim> &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<float>::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<double>(res.first[f]));
-
- // Second, calculate ln(U_k).
- for (auto &residual_element : res.second)
- residual_element = std::log(residual_element);
-
- const std::pair<double, double> fit =
- FESeries::linear_regression(ln_k, res.second);
- // Compute regularity s = mu - dim/2
- regularity = static_cast<float>(-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<float>();
- }
+ 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<dim>(
+ expansion_coefficients,
+ [n_modes](const TableIndices<dim> &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<float>::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<double>(res.first[f]));
+
+ // Second, calculate ln(U_k).
+ for (auto &residual_element : res.second)
+ residual_element = std::log(residual_element);
+
+ const std::pair<double, double> fit =
+ FESeries::linear_regression(ln_k, res.second);
+ // Compute regularity s = mu - dim/2
+ regularity = static_cast<float>(-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<float>();
+ }
}
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<double>::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<dim> 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<double, double> 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<double>::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<dim> 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<double, double> 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<float>(k_v);
- }
- else
smoothness_indicators(cell->active_cell_index()) =
- numbers::signaling_nan<float>();
- }
+ static_cast<float>(k_v);
+ }
+ else
+ smoothness_indicators(cell->active_cell_index()) =
+ numbers::signaling_nan<float>();
+ }
}
projects / dealii.git / commitdiff