IndexSet is_locally_owned_fine(cell_id_translator.size());
IndexSet is_locally_owned_coarse(cell_id_translator.size());
- for (const auto &cell : dof_handler_fine.active_cell_iterators())
- if (cell->is_locally_owned())
- is_locally_owned_fine.add_index(cell_id_translator.translate(cell));
+ for (const auto &cell : dof_handler_fine.active_cell_iterators() |
+ IteratorFilters::LocallyOwnedCell())
+ is_locally_owned_fine.add_index(cell_id_translator.translate(cell));
for (const auto &cell : dof_handler_coarse.active_cell_iterators() |
IteratorFilters::LocallyOwnedCell())
// 4a) set all cells artificial (and set the actual
// (level_)subdomain_ids in the next step)
- for (auto cell = this->begin(); cell != this->end(); ++cell)
+ for (const auto &cell : this->cell_iterators())
{
if (cell->is_active())
cell->set_subdomain_id(
# include <deal.II/distributed/grid_refinement.h>
+# include <deal.II/grid/filtered_iterator.h>
# include <deal.II/grid/grid_refinement.h>
# include <deal.II/grid/tria.h>
# include <deal.II/grid/tria_accessor.h>
ExcInternalError());
unsigned int owned_index = 0;
- for (const auto &cell : tria.active_cell_iterators())
- if (cell->subdomain_id() == tria.locally_owned_subdomain())
- {
- locally_owned_indicators(owned_index) =
- criteria(cell->active_cell_index());
- ++owned_index;
- }
+ for (const auto &cell :
+ tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+ {
+ locally_owned_indicators(owned_index) =
+ criteria(cell->active_cell_index());
+ ++owned_index;
+ }
Assert(owned_index == tria.n_locally_owned_active_cells(),
ExcInternalError());
}
// step 1) check if all processes have enough cells
- unsigned int n_locally_owned_active_cells = 0;
- for (const auto &cell : tria_in.active_cell_iterators())
- if (cell->is_locally_owned())
- ++n_locally_owned_active_cells;
+ const auto locally_owned_cells =
+ tria_in.active_cell_iterators() | IteratorFilters::LocallyOwnedCell();
+ const unsigned int n_locally_owned_active_cells =
+ std::distance(locally_owned_cells.begin(), locally_owned_cells.end());
const auto comm = tria_in.get_communicator();
const auto n_subdomains = Utilities::MPI::n_mpi_processes(mpi_communicator);
// determine weight of each cell
- for (const auto &cell : tria->active_cell_iterators())
- if (cell->is_locally_owned())
- weights[partitioner->global_to_local(
- cell->global_active_cell_index())] =
- weighting_function(
- cell, Triangulation<dim, spacedim>::CellStatus::CELL_PERSIST);
+ for (const auto &cell :
+ tria->active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+ weights[partitioner->global_to_local(cell->global_active_cell_index())] =
+ weighting_function(
+ cell, Triangulation<dim, spacedim>::CellStatus::CELL_PERSIST);
// determine weight of all the cells locally owned by this process
uint64_t process_local_weight = 0;
# ifdef DEBUG
{
// Assert that each cell is owned by a processor
- unsigned int n_my_cells = 0;
- typename parallel::shared::Triangulation<dim,
- spacedim>::active_cell_iterator
- cell = this->begin_active(),
- endc = this->end();
- for (; cell != endc; ++cell)
- if (cell->is_locally_owned())
- n_my_cells += 1;
+ const unsigned int n_my_cells = std::count_if(
+ this->begin_active(),
+ typename Triangulation<dim, spacedim>::active_cell_iterator(
+ this->end()),
+ [](const auto &i) { return (i.is_locally_owned()); });
const unsigned int total_cells =
Utilities::MPI::sum(n_my_cells, this->get_communicator());
// cell is owned by a processor
if (settings & construct_multigrid_hierarchy)
{
- unsigned int n_my_cells = 0;
- typename parallel::shared::Triangulation<dim, spacedim>::cell_iterator
- cell = this->begin(),
- endc = this->end();
- for (; cell != endc; ++cell)
- if (cell->is_locally_owned_on_level())
- n_my_cells += 1;
+ const unsigned int n_my_cells =
+ std::count_if(this->begin(), this->end(), [](const auto &i) {
+ return (i.is_locally_owned_on_level());
+ });
+
const unsigned int total_cells =
Utilities::MPI::sum(n_my_cells, this->get_communicator());
#include <deal.II/distributed/shared_tria.h>
#include <deal.II/distributed/tria_base.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>
}
if (this->n_levels() > 0)
- for (const auto &cell : this->active_cell_iterators())
- if (cell->subdomain_id() == my_subdomain)
- ++number_cache.n_locally_owned_active_cells;
+ number_cache.n_locally_owned_active_cells = std::count_if(
+ this->begin_active(),
+ typename Triangulation<dim, spacedim>::active_cell_iterator(
+ this->end()),
+ [](const auto &i) { return i.is_locally_owned(); });
+ else
+ number_cache.n_locally_owned_active_cells = 0;
- // Potentially cast to a 64 bit type before accumulating to avoid overflow:
+ // Potentially cast to a 64 bit type before accumulating to avoid
+ // overflow:
number_cache.n_global_active_cells =
Utilities::MPI::sum(static_cast<types::global_cell_index>(
number_cache.n_locally_owned_active_cells),
number_cache.n_global_levels =
Utilities::MPI::max(this->n_levels(), this->mpi_communicator);
- // Store MPI ranks of level ghost owners of this processor on all levels.
+ // Store MPI ranks of level ghost owners of this processor on all
+ // levels.
if (this->is_multilevel_hierarchy_constructed() == true)
{
number_cache.level_ghost_owners.clear();
return;
// find level ghost owners
- for (typename Triangulation<dim, spacedim>::cell_iterator cell =
- this->begin();
- cell != this->end();
- ++cell)
+ for (const auto &cell : this->cell_iterators())
if (cell->level_subdomain_id() != numbers::artificial_subdomain_id &&
cell->level_subdomain_id() != this->locally_owned_subdomain())
this->number_cache.level_ghost_owners.insert(
cell->level_subdomain_id());
# ifdef DEBUG
- // Check that level_ghost_owners is symmetric by sending a message to
- // everyone
+ // Check that level_ghost_owners is symmetric by sending a message
+ // to everyone
{
int ierr = MPI_Barrier(this->mpi_communicator);
AssertThrowMPI(ierr);
cell_sizes_variable_cumulative.end(),
cell_sizes_variable_cumulative.begin());
- // Serialize cumulative variable size vector value-by-value.
- // This way we can circumvent the overhead of storing the
- // container object as a whole, since we know its size by
- // the number of registered callback functions.
+ // Serialize cumulative variable size vector
+ // value-by-value. This way we can circumvent the overhead
+ // of storing the container object as a whole, since we
+ // know its size by the number of registered callback
+ // functions.
data_fixed_it->resize(n_callbacks_variable *
sizeof(unsigned int));
for (unsigned int i = 0; i < n_callbacks_variable; ++i)
// functions (i.e. a cell that was not flagged with CELL_INVALID)
// and store the sizes of each buffer.
//
- // To deal with the case that at least one of the processors does not own
- // any cell at all, we will exchange the information about the data sizes
- // among them later. The code in between is still well-defined, since the
- // following loops will be skipped.
+ // To deal with the case that at least one of the processors does not
+ // own any cell at all, we will exchange the information about the data
+ // sizes among them later. The code in between is still well-defined,
+ // since the following loops will be skipped.
std::vector<unsigned int> local_sizes_fixed(
1 + n_callbacks_fixed + (variable_size_data_stored ? 1 : 0));
for (const auto &data_cell : packed_fixed_size_data)
src_sizes_variable.end(),
std::vector<int>::size_type(0));
- // Move every piece of packed fixed size data into the consecutive buffer.
+ // Move every piece of packed fixed size data into the consecutive
+ // buffer.
src_data_fixed.reserve(expected_size_fixed);
for (const auto &data_cell_fixed : packed_fixed_size_data)
{
{
// We decode the handle returned by register_data_attach() back into
// a format we can use. All even handles belong to those callback
- // functions which write/read variable size data, all odd handles interact
- // with fixed size buffers.
+ // functions which write/read variable size data, all odd handles
+ // interact with fixed size buffers.
const bool callback_variable_transfer = (handle % 2 == 0);
const unsigned int callback_index = handle / 2;
spacedim>::CELL_INVALID)
{
// Extract the corresponding values for offset and size from
- // the cumulative sizes array stored in the fixed size buffer.
+ // the cumulative sizes array stored in the fixed size
+ // buffer.
if (callback_index == 0)
offset = 0;
else
// ------------------
// Write cumulative sizes to file.
- // Since each processor owns the same information about the data sizes,
- // it is sufficient to let only the first processor perform this task.
+ // Since each processor owns the same information about the data
+ // sizes, it is sufficient to let only the first processor perform
+ // this task.
if (myrank == 0)
{
ierr = MPI_File_write_at(fh,
const MPI_Offset size_header =
sizes_fixed_cumulative.size() * sizeof(unsigned int);
- // Make sure we do the following computation in 64bit integers to be able
- // to handle 4GB+ files:
+ // Make sure we do the following computation in 64bit integers to be
+ // able to handle 4GB+ files:
const MPI_Offset my_global_file_position =
size_header +
static_cast<MPI_Offset>(global_first_cell) * bytes_per_cell;
const MPI_Offset my_global_file_position =
static_cast<MPI_Offset>(global_first_cell) * sizeof(unsigned int);
- // It is very unlikely that a single process has more than 2 billion
- // cells, but we might as well check.
+ // It is very unlikely that a single process has more than
+ // 2 billion cells, but we might as well check.
AssertThrow(src_sizes_variable.size() <
static_cast<std::size_t>(
std::numeric_limits<int>::max()),
AssertThrowMPI(ierr);
// Read cumulative sizes from file.
- // Since all processors need the same information about the data sizes,
- // let each of them retrieve it by reading from the same location in
- // the file.
+ // Since all processors need the same information about the data
+ // sizes, let each of them retrieve it by reading from the same
+ // location in the file.
sizes_fixed_cumulative.resize(1 + n_attached_deserialize_fixed +
(variable_size_data_stored ? 1 : 0));
ierr = MPI_File_read_at(fh,
const MPI_Offset size_header =
sizes_fixed_cumulative.size() * sizeof(unsigned int);
- // Make sure we do the following computation in 64bit integers to be able
- // to handle 4GB+ files:
+ // Make sure we do the following computation in 64bit integers to be
+ // able to handle 4GB+ files:
const MPI_Offset my_global_file_position =
size_header +
static_cast<MPI_Offset>(global_first_cell) * bytes_per_cell;
AssertThrowMPI(ierr);
- // Compute my data size in bytes and compute prefix sum. We do this in
- // 64 bit to avoid overflow for files larger than 4 GB:
+ // Compute my data size in bytes and compute prefix sum. We do this
+ // in 64 bit to avoid overflow for files larger than 4 GB:
const std::uint64_t size_on_proc =
std::accumulate(dest_sizes_variable.begin(),
dest_sizes_variable.end(),
std::vector<types::global_dof_index> dof_indices;
std::set<types::global_dof_index> global_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 (const types::global_dof_index dof_index : dof_indices)
- if (!dof_set.is_element(dof_index))
- global_dof_indices.insert(dof_index);
- }
+ for (const types::global_dof_index dof_index : dof_indices)
+ if (!dof_set.is_element(dof_index))
+ global_dof_indices.insert(dof_index);
+ }
dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end());
// In a first step, obtain the weights of the locally owned
// cells. For all others, the weight remains at the zero the
// vector was initialized with above.
- for (const auto &cell : triangulation.active_cell_iterators())
- if (cell->is_locally_owned())
- cell_weights[cell->active_cell_index()] =
- triangulation.signals.cell_weight(
- cell, Triangulation<dim, spacedim>::CellStatus::CELL_PERSIST);
+ for (const auto &cell : triangulation.active_cell_iterators() |
+ IteratorFilters::LocallyOwnedCell())
+ cell_weights[cell->active_cell_index()] =
+ triangulation.signals.cell_weight(
+ cell, Triangulation<dim, spacedim>::CellStatus::CELL_PERSIST);
// If this is a parallel triangulation, we then need to also
// get the weights for all other cells. We have asserted above
typename Triangulation<dim, spacedim>::active_cell_iterator>>
boxes;
boxes.reserve(tria->n_active_cells());
- for (const auto &cell : tria->active_cell_iterators())
- if (cell->is_locally_owned())
- boxes.emplace_back(
- std::make_pair(mapping->get_bounding_box(cell), cell));
+ for (const auto &cell : tria->active_cell_iterators() |
+ IteratorFilters::LocallyOwnedCell())
+ boxes.emplace_back(mapping->get_bounding_box(cell), cell);
locally_owned_cell_bounding_boxes_rtree = pack_rtree(boxes);
update_flags =
// deep copy error indicators
predicted_errors = error_indicators;
- for (const auto &cell : dof_handler.active_cell_iterators())
- if (cell->is_locally_owned())
- {
- // current cell will not be adapted
- if (!(cell->future_fe_index_set()) && !(cell->refine_flag_set()) &&
- !(cell->coarsen_flag_set()))
- {
- predicted_errors[cell->active_cell_index()] *= gamma_n;
- continue;
- }
+ for (const auto &cell : dof_handler.active_cell_iterators() |
+ IteratorFilters::LocallyOwnedCell())
+ {
+ // current cell will not be adapted
+ if (!(cell->future_fe_index_set()) && !(cell->refine_flag_set()) &&
+ !(cell->coarsen_flag_set()))
+ {
+ predicted_errors[cell->active_cell_index()] *= gamma_n;
+ continue;
+ }
- // current cell will be adapted
- // determine degree of its future finite element
- if (cell->coarsen_flag_set())
- {
- // cell will be coarsened, thus determine future finite element
- // on parent cell
- const auto &parent = cell->parent();
- if (future_fe_indices_on_coarsened_cells.find(parent) ==
- future_fe_indices_on_coarsened_cells.end())
- {
+ // current cell will be adapted
+ // determine degree of its future finite element
+ if (cell->coarsen_flag_set())
+ {
+ // cell will be coarsened, thus determine future finite element
+ // on parent cell
+ const auto &parent = cell->parent();
+ if (future_fe_indices_on_coarsened_cells.find(parent) ==
+ future_fe_indices_on_coarsened_cells.end())
+ {
#ifdef DEBUG
- for (const auto &child : parent->child_iterators())
- Assert(child->is_active() && child->coarsen_flag_set(),
- typename dealii::Triangulation<
- dim>::ExcInconsistentCoarseningFlags());
+ for (const auto &child : parent->child_iterators())
+ Assert(child->is_active() && child->coarsen_flag_set(),
+ typename dealii::Triangulation<
+ dim>::ExcInconsistentCoarseningFlags());
#endif
- parent_future_fe_index =
- dealii::internal::hp::DoFHandlerImplementation::
- dominated_future_fe_on_children<dim, spacedim>(parent);
+ parent_future_fe_index =
+ dealii::internal::hp::DoFHandlerImplementation::
+ dominated_future_fe_on_children<dim, spacedim>(parent);
- future_fe_indices_on_coarsened_cells.insert(
- {parent, parent_future_fe_index});
- }
- else
- {
- parent_future_fe_index =
- future_fe_indices_on_coarsened_cells[parent];
- }
+ future_fe_indices_on_coarsened_cells.insert(
+ {parent, parent_future_fe_index});
+ }
+ else
+ {
+ parent_future_fe_index =
+ future_fe_indices_on_coarsened_cells[parent];
+ }
- future_fe_degree =
- dof_handler.get_fe_collection()[parent_future_fe_index]
- .degree;
- }
- else
- {
- // future finite element on current cell is already set
- future_fe_degree =
- dof_handler.get_fe_collection()[cell->future_fe_index()]
- .degree;
- }
+ future_fe_degree =
+ dof_handler.get_fe_collection()[parent_future_fe_index].degree;
+ }
+ else
+ {
+ // future finite element on current cell is already set
+ future_fe_degree =
+ dof_handler.get_fe_collection()[cell->future_fe_index()].degree;
+ }
- // step 1: exponential decay with p-adaptation
- if (cell->future_fe_index_set())
- {
- predicted_errors[cell->active_cell_index()] *=
- std::pow(gamma_p,
- int(future_fe_degree) - int(cell->get_fe().degree));
- }
+ // step 1: exponential decay with p-adaptation
+ if (cell->future_fe_index_set())
+ {
+ predicted_errors[cell->active_cell_index()] *=
+ std::pow(gamma_p,
+ int(future_fe_degree) - int(cell->get_fe().degree));
+ }
- // step 2: algebraic decay with h-adaptation
- if (cell->refine_flag_set())
- {
- predicted_errors[cell->active_cell_index()] *=
- (gamma_h * std::pow(.5, future_fe_degree));
+ // step 2: algebraic decay with h-adaptation
+ if (cell->refine_flag_set())
+ {
+ predicted_errors[cell->active_cell_index()] *=
+ (gamma_h * std::pow(.5, future_fe_degree));
- // predicted error will be split on children cells
- // after adaptation via CellDataTransfer
- }
- else if (cell->coarsen_flag_set())
- {
- predicted_errors[cell->active_cell_index()] /=
- (gamma_h * std::pow(.5, future_fe_degree));
+ // predicted error will be split on children cells
+ // after adaptation via CellDataTransfer
+ }
+ else if (cell->coarsen_flag_set())
+ {
+ predicted_errors[cell->active_cell_index()] /=
+ (gamma_h * std::pow(.5, future_fe_degree));
- // predicted error will be summed up on parent cell
- // after adaptation via CellDataTransfer
- }
- }
+ // predicted error will be summed up on parent cell
+ // after adaptation via CellDataTransfer
+ }
+ }
}
partitioner = std::make_shared<Utilities::MPI::Partitioner>(
patch_dof_handler.locally_owned_dofs(), active_dofs, MPI_COMM_WORLD);
- for (const auto &cell : patch_dof_handler.active_cell_iterators())
+ for (const auto &cell : patch_dof_handler.active_cell_iterators() |
+ IteratorFilters::LocallyOwnedCell())
{
- if (cell->is_locally_owned() == false)
- continue;
-
fe_values.reinit(cell);
cell->get_dof_indices(dof_indices);
Vector<number> local_dof_values;
std::vector<double> converted_indices;
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_legendre.get_n_coefficients_per_direction(
- cell->active_fe_index());
- resize(expansion_coefficients, n_modes);
-
- local_dof_values.reinit(cell->get_fe().n_dofs_per_cell());
- cell->get_dof_values(solution, local_dof_values);
-
- fe_legendre.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_sum_less_than_N(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.
- converted_indices.assign(res.first.begin(),
- res.first.end());
-
- for (auto &residual_element : res.second)
- residual_element = std::log(residual_element);
-
- const std::pair<double, double> fit =
- FESeries::linear_regression(converted_indices,
- res.second);
- regularity = static_cast<float>(-fit.first);
- }
-
- 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_legendre.get_n_coefficients_per_direction(
+ cell->active_fe_index());
+ resize(expansion_coefficients, n_modes);
+
+ local_dof_values.reinit(cell->get_fe().n_dofs_per_cell());
+ cell->get_dof_values(solution, local_dof_values);
+
+ fe_legendre.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_sum_less_than_N(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.
+ converted_indices.assign(res.first.begin(), res.first.end());
+
+ for (auto &residual_element : res.second)
+ residual_element = std::log(residual_element);
+
+ const std::pair<double, double> fit =
+ FESeries::linear_regression(converted_indices, res.second);
+ regularity = static_cast<float>(-fit.first);
+ }
+
+ 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_legendre.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_legendre.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
- for (unsigned int i = 0; 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(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
+ 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_legendre.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_legendre.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
+ for (unsigned int i = 0; 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(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()) =
- numbers::signaling_nan<float>();
- }
+ static_cast<float>(k_v);
+ }
+ else
+ smoothness_indicators(cell->active_cell_index()) =
+ numbers::signaling_nan<float>();
+ }
}
parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
GridGenerator::hyper_cube(tria);
tria.refine_global(1);
- for (const auto &cell : tria.active_cell_iterators())
- if (cell->is_locally_owned())
- if (cell->center().norm() < 0.2)
- cell->set_refine_flag();
+ for (const auto &cell :
+ tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+ if (cell->center().norm() < 0.2)
+ cell->set_refine_flag();
tria.execute_coarsening_and_refinement();
if (dim < 3 && fe_degree < 2)
tria.refine_global(2);
// ----- verify ------
deallog << "pre_adaptation" << std::endl;
- for (const auto &cell : dh.active_cell_iterators())
- if (cell->is_locally_owned())
- {
- deallog << " cell:" << cell->id().to_string()
- << " fe_deg:" << cell->get_fe().degree
- << " error:" << error_indicators[cell->active_cell_index()];
+ for (const auto &cell :
+ dh.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+ {
+ deallog << " cell:" << cell->id().to_string()
+ << " fe_deg:" << cell->get_fe().degree
+ << " error:" << error_indicators[cell->active_cell_index()];
- if (cell->coarsen_flag_set())
- deallog << " coarsening";
- else if (cell->refine_flag_set())
- deallog << " refining";
+ if (cell->coarsen_flag_set())
+ deallog << " coarsening";
+ else if (cell->refine_flag_set())
+ deallog << " refining";
- if (cell->future_fe_index_set())
- deallog << " future_fe_deg:" << fes[cell->future_fe_index()].degree;
+ if (cell->future_fe_index_set())
+ deallog << " future_fe_deg:" << fes[cell->future_fe_index()].degree;
- deallog << std::endl;
- }
+ deallog << std::endl;
+ }
// ----- execute adaptation -----
parallel::distributed::CellDataTransfer<dim, dim, Vector<float>>
// ----- verify ------
deallog << "pre_adaptation" << std::endl;
- for (const auto &cell : dh.active_cell_iterators())
- if (cell->is_locally_owned())
- {
- deallog << " cell:" << cell->id().to_string()
- << " fe_deg:" << cell->get_fe().degree
- << " error:" << error_indicators[cell->active_cell_index()];
-
- if (cell->coarsen_flag_set())
- deallog << " coarsening";
- else if (cell->refine_flag_set())
- deallog << " refining";
-
- if (cell->future_fe_index_set())
- deallog << " future_fe_deg:" << fes[cell->future_fe_index()].degree;
-
- deallog << std::endl;
- }
+ for (const auto &cell :
+ dh.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+ {
+ deallog << " cell:" << cell->id().to_string()
+ << " fe_deg:" << cell->get_fe().degree
+ << " error:" << error_indicators[cell->active_cell_index()];
+
+ if (cell->coarsen_flag_set())
+ deallog << " coarsening";
+ else if (cell->refine_flag_set())
+ deallog << " refining";
+
+ if (cell->future_fe_index_set())
+ deallog << " future_fe_deg:" << fes[cell->future_fe_index()].degree;
+
+ deallog << std::endl;
+ }
// ----- execute adaptation -----
parallel::distributed::CellDataTransfer<dim, dim, Vector<float>>
// display number of cells for each FE index
std::vector<unsigned int> count(fes.size(), 0);
- for (const auto &cell : dofh.active_cell_iterators())
- if (cell->is_locally_owned())
- count[cell->active_fe_index()]++;
+ for (const auto &cell :
+ dofh.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+ count[cell->active_fe_index()]++;
Utilities::MPI::sum(count, tria.get_communicator(), count);
deallog << "fe count:" << count << std::endl;
// display number of cells for each FE index
std::vector<unsigned int> count(fes.size(), 0);
- for (const auto &cell : dofh.active_cell_iterators())
- if (cell->is_locally_owned())
- count[cell->active_fe_index()]++;
+ for (const auto &cell :
+ dofh.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+ count[cell->active_fe_index()]++;
Utilities::MPI::sum(count, tria.get_communicator(), count);
deallog << "cycle:" << i << ", fe count:" << count << std::endl;
}
// set future_fe_indices
unsigned int future_feidx = 0;
- for (const auto &cell : dh.active_cell_iterators())
- if (cell->is_locally_owned())
- {
- // check if cell is initialized correctly
- Assert(cell->active_fe_index() == 0, ExcInternalError());
-
- cell->set_future_fe_index(future_feidx);
- future_feidx = ((future_feidx + 1) < fe.size()) ? future_feidx + 1 : 0;
- }
+ for (const auto &cell :
+ dh.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+ {
+ // check if cell is initialized correctly
+ Assert(cell->active_fe_index() == 0, ExcInternalError());
+
+ cell->set_future_fe_index(future_feidx);
+ future_feidx = ((future_feidx + 1) < fe.size()) ? future_feidx + 1 : 0;
+ }
dh.distribute_dofs(fe);
tria.execute_coarsening_and_refinement();
Vector<float> indicator(tria.n_active_cells());
// assign each cell a globally unique cellid
- for (const auto &cell : tria.active_cell_iterators())
- if (cell->is_locally_owned())
- {
- const std::string cellid = cell->id().to_string();
- const unsigned int fine_cellid =
- std::stoul(cellid.substr(cellid.find(':') + 1, std::string::npos));
+ for (const auto &cell :
+ tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+ {
+ const std::string cellid = cell->id().to_string();
+ const unsigned int fine_cellid =
+ std::stoul(cellid.substr(cellid.find(':') + 1, std::string::npos));
- indicator[cell->active_cell_index()] = fine_cellid + 1;
- }
+ indicator[cell->active_cell_index()] = fine_cellid + 1;
+ }
deallog << "l1-norm: ";
parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction(
Vector<float> indicator(tria.n_active_cells());
// assign each cell a globally unique cellid
- for (const auto &cell : tria.active_cell_iterators())
- if (cell->is_locally_owned())
- {
- const std::string cellid = cell->id().to_string();
- const unsigned int fine_cellid =
- std::stoul(cellid.substr(cellid.find(':') + 1, std::string::npos));
+ for (const auto &cell :
+ tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+ {
+ const std::string cellid = cell->id().to_string();
+ const unsigned int fine_cellid =
+ std::stoul(cellid.substr(cellid.find(':') + 1, std::string::npos));
- Testing::srand(fine_cellid);
- indicator[cell->active_cell_index()] = random_value<float>();
- }
+ Testing::srand(fine_cellid);
+ indicator[cell->active_cell_index()] = random_value<float>();
+ }
deallog << "l1-norm: ";
parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction(
dofh.distribute_dofs(fe);
deallog << "Before:" << std::endl;
- for (const auto &cell : dofh.active_cell_iterators())
- if (cell->is_locally_owned())
- {
- deallog << "locally owned cell: " << cell << std::endl;
- deallog << " dof indices: ";
-
- std::vector<types::global_dof_index> cell_dofs(
- cell->get_fe().dofs_per_cell);
- cell->get_dof_indices(cell_dofs);
-
- for (auto i : cell_dofs)
- deallog << i << ' ';
- deallog << std::endl;
- }
+ for (const auto &cell :
+ dofh.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+ {
+ deallog << "locally owned cell: " << cell << std::endl;
+ deallog << " dof indices: ";
+
+ std::vector<types::global_dof_index> cell_dofs(
+ cell->get_fe().dofs_per_cell);
+ cell->get_dof_indices(cell_dofs);
+
+ for (auto i : cell_dofs)
+ deallog << i << ' ';
+ deallog << std::endl;
+ }
std::set<types::global_dof_index> starting_indices;
for (const auto &cell :
// setup dof-handlers
DoFHandler<dim> dof_handler_fine(tria_fine);
- for (const auto &cell : dof_handler_fine.active_cell_iterators())
- if (cell->is_locally_owned())
- cell->set_active_fe_index(0);
+ for (const auto &cell : dof_handler_fine.active_cell_iterators() |
+ IteratorFilters::LocallyOwnedCell())
+ cell->set_active_fe_index(0);
dof_handler_fine.distribute_dofs(fe);
DoFHandler<dim> dof_handler_coarse(tria_coarse);
// ----- gather -----
// store parent id of all locally owned cells
Vector<PetscScalar> cell_ids_pre(tria.n_active_cells());
- for (const auto &cell : tria.active_cell_iterators())
- if (cell->is_locally_owned())
- {
- const std::string parent_cellid = cell->parent()->id().to_string();
- const unsigned int parent_coarse_cell_id =
- static_cast<unsigned int>(std::stoul(parent_cellid));
- cell_ids_pre(cell->active_cell_index()) = parent_coarse_cell_id;
- }
+ for (const auto &cell :
+ tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+ {
+ const std::string parent_cellid = cell->parent()->id().to_string();
+ const unsigned int parent_coarse_cell_id =
+ static_cast<unsigned int>(std::stoul(parent_cellid));
+ cell_ids_pre(cell->active_cell_index()) = parent_coarse_cell_id;
+ }
// distribute local vector (as presented in step-18)
PETScWrappers::MPI::Vector distributed_cell_ids_pre(
// ----- gather -----
// store parent id of all locally owned cells
Vector<PetscScalar> cell_ids_pre(tria.n_active_cells());
- for (const auto &cell : tria.active_cell_iterators())
- if (cell->is_locally_owned())
- {
- const std::string parent_cellid = cell->parent()->id().to_string();
- const unsigned int parent_coarse_cell_id =
- static_cast<unsigned int>(std::stoul(parent_cellid));
- cell_ids_pre(cell->active_cell_index()) = parent_coarse_cell_id;
- }
+ for (const auto &cell :
+ tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+ {
+ const std::string parent_cellid = cell->parent()->id().to_string();
+ const unsigned int parent_coarse_cell_id =
+ static_cast<unsigned int>(std::stoul(parent_cellid));
+ cell_ids_pre(cell->active_cell_index()) = parent_coarse_cell_id;
+ }
// distribute local vector (as presented in step-18)
PETScWrappers::MPI::Vector distributed_cell_ids_pre(
// display number of cells for each FE index
std::vector<unsigned int> count(fes.size(), 0);
- for (const auto &cell : dofh.active_cell_iterators())
- if (cell->is_locally_owned())
- count[cell->active_fe_index()]++;
+ for (const auto &cell :
+ dofh.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+ count[cell->active_fe_index()]++;
Utilities::MPI::sum(count, tria.get_communicator(), count);
deallog << "fe count:" << count << std::endl;
// display number of cells for each FE index
std::vector<unsigned int> count(fes.size(), 0);
- for (const auto &cell : dofh.active_cell_iterators())
- if (cell->is_locally_owned())
- count[cell->active_fe_index()]++;
+ for (const auto &cell :
+ dofh.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+ count[cell->active_fe_index()]++;
Utilities::MPI::sum(count, tria.get_communicator(), count);
deallog << "cycle:" << i << ", fe count:" << count << std::endl;
}
Assert(fe_indices_changed, ExcInternalError());
deallog << "future FE indices before adaptation:" << std::endl;
- for (const auto &cell : dofh.active_cell_iterators())
- if (cell->is_locally_owned())
- deallog << " " << cell->id().to_string() << " " << cell->future_fe_index()
- << std::endl;
+ for (const auto &cell :
+ dofh.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+ deallog << " " << cell->id().to_string() << " " << cell->future_fe_index()
+ << std::endl;
tria.execute_coarsening_and_refinement();
quadrature_formula.size());
unsigned int history_index = 0;
- for (auto &cell : triangulation.active_cell_iterators())
- if (cell->is_locally_owned())
- {
- cell->set_user_pointer(&quadrature_point_history[history_index]);
- history_index += quadrature_formula.size();
- }
+ for (auto &cell : triangulation.active_cell_iterators() |
+ IteratorFilters::LocallyOwnedCell())
+ {
+ cell->set_user_pointer(&quadrature_point_history[history_index]);
+ history_index += quadrature_formula.size();
+ }
Assert(history_index == quadrature_point_history.size(),
ExcInternalError());
projects / dealii.git / commitdiff