Added fixed fraction tests to compare serial/parallel implementations.
--- /dev/null
+New: GridRefinement::refine_and_coarsen_fixed_fraction() and
+parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction()
+now allow to specify a VectorTools::NormType, which determines how
+combined errors on subsets of cells will be calculated.
+<br>
+(Marc Fehling, 2020/10/06)
#include <deal.II/distributed/tria.h>
+#include <deal.II/numerics/vector_tools_common.h>
+
#include <limits>
#include <vector>
* refined at all.
*
* The same is true for the fraction of cells that is coarsened.
+ *
+ * @param[in,out] tria The triangulation whose cells this function is
+ * supposed to mark for coarsening and refinement.
+ *
+ * @param[in] criteria The refinement criterion for each mesh cell active
+ * on the current triangulation. Entries may not be negative.
+ *
+ * @param[in] top_fraction_of_cells The fraction of cells to be refined.
+ * If this number is zero, no cells will be refined. If it equals one,
+ * the result will be flagging for global refinement.
+ *
+ * @param[in] bottom_fraction_of_cells The fraction of cells to be
+ * coarsened. If this number is zero, no cells will be coarsened.
+ *
+ * @param[in] max_n_cells This argument can be used to specify a maximal
+ * number of cells. If this number is going to be exceeded upon
+ * refinement, then refinement and coarsening fractions are going to be
+ * adjusted in an attempt to reach the maximum number of cells. Be aware
+ * though that through proliferation of refinement due to
+ * Triangulation::MeshSmoothing, this number is only an indicator. The
+ * default value of this argument is to impose no limit on the number of
+ * cells.
*/
template <int dim, typename Number, int spacedim>
void
* processors, no cells are refined at all.
*
* The same is true for the fraction of cells that is coarsened.
+ *
+ * @param[in,out] tria The triangulation whose cells this function is
+ * supposed to mark for coarsening and refinement.
+ *
+ * @param[in] criteria The refinement criterion computed on each mesh cell
+ * active on the current triangulation. Entries may not be negative.
+ *
+ * @param[in] top_fraction_of_error The fraction of the total estimate
+ * which should be refined. If this number is zero, no cells will be
+ * refined. If it equals one, the result will be flagging for global
+ * refinement.
+ *
+ * @param[in] bottom_fraction_of_error The fraction of the estimate
+ * coarsened. If this number is zero, no cells will be coarsened.
+ *
+ * @param[in] norm_type To determine thresholds, combined errors on
+ * subsets of cells are calculated as norms of the criteria on these
+ * cells. Different types of norms can be used for this purpose, from
+ * which VectorTools::NormType::L1_norm and
+ * VectorTools::NormType::L2_norm are currently supported.
*/
template <int dim, typename Number, int spacedim>
void
refine_and_coarsen_fixed_fraction(
parallel::distributed::Triangulation<dim, spacedim> &tria,
const dealii::Vector<Number> & criteria,
- const double top_fraction_of_error,
- const double bottom_fraction_of_error);
+ const double top_fraction_of_error,
+ const double bottom_fraction_of_error,
+ const VectorTools::NormType norm_type = VectorTools::NormType::L1_norm);
} // namespace GridRefinement
} // namespace distributed
} // namespace parallel
#include <deal.II/base/exceptions.h>
+#include <deal.II/numerics/vector_tools_common.h>
+
#include <limits>
DEAL_II_NAMESPACE_OPEN
* through proliferation of refinement due to Triangulation::MeshSmoothing,
* this number is only an indicator. The default value of this argument is
* to impose no limit on the number of cells.
+ *
+ * @param[in] norm_type To determine thresholds, combined errors on
+ * subsets of cells are calculated as norms of the criteria on these
+ * cells. Different types of norms can be used for this purpose, from
+ * which VectorTools::NormType::L1_norm and
+ * VectorTools::NormType::L2_norm are currently supported.
*/
template <int dim, typename Number, int spacedim>
void
const Vector<Number> & criteria,
const double top_fraction,
const double bottom_fraction,
- const unsigned int max_n_cells = std::numeric_limits<unsigned int>::max());
+ const unsigned int max_n_cells = std::numeric_limits<unsigned int>::max(),
+ const VectorTools::NormType norm_type = VectorTools::NormType::L1_norm);
cell->clear_coarsen_flag();
}
}
+
+
+
+ /**
+ * Fixed fraction algorithm without a specified vector norm.
+ *
+ * Entries of the criteria vector and fractions are taken as is, so this
+ * function basically evaluates norms on the vector or its subsets as
+ * l1-norms.
+ */
+ template <int dim, int spacedim, typename Number>
+ void
+ refine_and_coarsen_fixed_fraction_via_l1_norm(
+ parallel::distributed::Triangulation<dim, spacedim> &tria,
+ const dealii::Vector<Number> & criteria,
+ const double top_fraction_of_error,
+ const double bottom_fraction_of_error)
+ {
+ // first extract from the vector of indicators the ones that correspond
+ // to cells that we locally own
+ Vector<Number> locally_owned_indicators(
+ tria.n_locally_owned_active_cells());
+ get_locally_owned_indicators(tria, criteria, locally_owned_indicators);
+
+ MPI_Comm mpi_communicator = tria.get_communicator();
+
+ // figure out the global max and min of the indicators. we don't need it
+ // here, but it's a collective communication call
+ const std::pair<double, double> global_min_and_max =
+ dealii::internal::parallel::distributed::GridRefinement::
+ compute_global_min_and_max_at_root(locally_owned_indicators,
+ mpi_communicator);
+
+ const double total_error =
+ compute_global_sum(locally_owned_indicators, mpi_communicator);
+
+ double top_target_error = top_fraction_of_error * total_error,
+ bottom_target_error = (1. - bottom_fraction_of_error) * total_error;
+
+ double top_threshold, bottom_threshold;
+ top_threshold = dealii::internal::parallel::distributed::GridRefinement::
+ RefineAndCoarsenFixedFraction::compute_threshold(locally_owned_indicators,
+ global_min_and_max,
+ top_target_error,
+ mpi_communicator);
+
+ // compute bottom threshold only if necessary. otherwise use the lowest
+ // threshold possible
+ if (bottom_fraction_of_error > 0)
+ bottom_threshold = dealii::internal::parallel::distributed::
+ GridRefinement::RefineAndCoarsenFixedFraction::compute_threshold(
+ locally_owned_indicators,
+ global_min_and_max,
+ bottom_target_error,
+ mpi_communicator);
+ else
+ bottom_threshold = std::numeric_limits<Number>::lowest();
+
+ // now refine the mesh
+ mark_cells(tria, criteria, top_threshold, bottom_threshold);
+ }
} // namespace
}
+
template <int dim, typename Number, int spacedim>
void
refine_and_coarsen_fixed_fraction(
parallel::distributed::Triangulation<dim, spacedim> &tria,
const dealii::Vector<Number> & criteria,
- const double top_fraction_of_error,
- const double bottom_fraction_of_error)
+ const double top_fraction_of_error,
+ const double bottom_fraction_of_error,
+ const VectorTools::NormType norm_type)
{
Assert(criteria.size() == tria.n_active_cells(),
ExcDimensionMismatch(criteria.size(), tria.n_active_cells()));
Assert(criteria.is_non_negative(),
dealii::GridRefinement::ExcNegativeCriteria());
- // first extract from the vector of indicators the ones that correspond
- // to cells that we locally own
- Vector<Number> locally_owned_indicators(
- tria.n_locally_owned_active_cells());
- get_locally_owned_indicators(tria, criteria, locally_owned_indicators);
-
- MPI_Comm mpi_communicator = tria.get_communicator();
-
- // figure out the global max and min of the indicators. we don't need it
- // here, but it's a collective communication call
- const std::pair<double, double> global_min_and_max =
- dealii::internal::parallel::distributed::GridRefinement::
- compute_global_min_and_max_at_root(locally_owned_indicators,
- mpi_communicator);
-
- const double total_error =
- compute_global_sum(locally_owned_indicators, mpi_communicator);
- double top_threshold, bottom_threshold;
- top_threshold = dealii::internal::parallel::distributed::
- GridRefinement::RefineAndCoarsenFixedFraction::compute_threshold(
- locally_owned_indicators,
- global_min_and_max,
- top_fraction_of_error * total_error,
- mpi_communicator);
-
- // compute bottom threshold only if necessary. otherwise use the lowest
- // threshold possible
- if (bottom_fraction_of_error > 0)
- bottom_threshold = dealii::internal::parallel::distributed::
- GridRefinement::RefineAndCoarsenFixedFraction::compute_threshold(
- locally_owned_indicators,
- global_min_and_max,
- (1. - bottom_fraction_of_error) * total_error,
- mpi_communicator);
- else
- bottom_threshold = std::numeric_limits<Number>::lowest();
+ switch (norm_type)
+ {
+ case VectorTools::NormType::L1_norm:
+ // evaluate norms on subsets and compare them as
+ // c_0 + c_1 + ... < fraction * l1-norm(c)
+ refine_and_coarsen_fixed_fraction_via_l1_norm(
+ tria,
+ criteria,
+ top_fraction_of_error,
+ bottom_fraction_of_error);
+ break;
+
+ case VectorTools::NormType::L2_norm:
+ {
+ // we do not want to evaluate norms on subsets as:
+ // sqrt(c_0^2 + c_1^2 + ...) < fraction * l2-norm(c)
+ // instead take the square of both sides of the equation
+ // and evaluate:
+ // c_0^2 + c_1^2 + ... < fraction^2 * l1-norm(c.c)
+ // we adjust all parameters accordingly
+ Vector<Number> criteria_squared(criteria.size());
+ std::transform(criteria.begin(),
+ criteria.end(),
+ criteria_squared.begin(),
+ [](Number c) { return c * c; });
+
+ refine_and_coarsen_fixed_fraction_via_l1_norm(
+ tria,
+ criteria_squared,
+ top_fraction_of_error * top_fraction_of_error,
+ bottom_fraction_of_error * bottom_fraction_of_error);
+ }
+ break;
- // now refine the mesh
- mark_cells(tria, criteria, top_threshold, bottom_threshold);
+ default:
+ Assert(false, ExcNotImplemented());
+ break;
+ }
}
} // namespace GridRefinement
} // namespace distributed
parallel::distributed::Triangulation<deal_II_dimension> &,
const dealii::Vector<S> &,
const double,
- const double);
+ const double,
+ const VectorTools::NormType);
\}
\}
\}
deal_II_dimension> &,
const dealii::Vector<S> &,
const double,
- const double);
+ const double,
+ const VectorTools::NormType);
\}
\}
\}
DEAL_II_NAMESPACE_OPEN
+namespace
+{
+ /**
+ * Fixed fraction algorithm without a specified vector norm.
+ *
+ * Entries of the criteria vector and fractions are taken as is, so this
+ * function basically evaluates norms on the vector or its subsets as
+ * l1-norms.
+ */
+ template <int dim, int spacedim, typename Number>
+ void
+ refine_and_coarsen_fixed_fraction_via_l1_norm(
+ Triangulation<dim, spacedim> &tria,
+ const Vector<Number> & criteria,
+ const double top_fraction,
+ const double bottom_fraction,
+ const unsigned int max_n_cells)
+ {
+ // sort the criteria in descending order in an auxiliary vector, which we
+ // have to sum up and compare with @p{fraction_of_error*total_error}
+ Vector<Number> criteria_sorted = criteria;
+ std::sort(criteria_sorted.begin(),
+ criteria_sorted.end(),
+ std::greater<double>());
+
+ const double total_error = criteria_sorted.l1_norm();
+
+ // compute thresholds
+ typename Vector<Number>::const_iterator pp = criteria_sorted.begin();
+ for (double sum = 0;
+ (sum < top_fraction * total_error) && (pp != criteria_sorted.end());
+ ++pp)
+ sum += *pp;
+ double top_threshold =
+ (pp != criteria_sorted.begin() ? (*pp + *(pp - 1)) / 2 : *pp);
+
+ typename Vector<Number>::const_iterator qq = criteria_sorted.end() - 1;
+ for (double sum = 0; (sum < bottom_fraction * total_error) &&
+ (qq != criteria_sorted.begin() - 1);
+ --qq)
+ sum += *qq;
+ double bottom_threshold =
+ ((qq != criteria_sorted.end() - 1) ? (*qq + *(qq + 1)) / 2 : 0.);
+
+ // we now have an idea how many cells we
+ // are going to refine and coarsen. we use
+ // this information to see whether we are
+ // over the limit and if so use a function
+ // that knows how to deal with this
+ // situation
+
+ // note, that at this point, we have no
+ // information about anisotropically refined
+ // cells, thus use the situation of purely
+ // isotropic refinement as guess for a mixed
+ // refinemnt as well.
+ const unsigned int refine_cells = pp - criteria_sorted.begin(),
+ coarsen_cells = criteria_sorted.end() - 1 - qq;
+
+ if (static_cast<unsigned int>(
+ tria.n_active_cells() +
+ refine_cells * (GeometryInfo<dim>::max_children_per_cell - 1) -
+ (coarsen_cells * (GeometryInfo<dim>::max_children_per_cell - 1) /
+ GeometryInfo<dim>::max_children_per_cell)) > max_n_cells)
+ {
+ GridRefinement::refine_and_coarsen_fixed_number(tria,
+ criteria,
+ 1. * refine_cells /
+ criteria.size(),
+ 1. * coarsen_cells /
+ criteria.size(),
+ max_n_cells);
+ return;
+ }
+
+ // in some rare cases it may happen that
+ // both thresholds are the same (e.g. if
+ // there are many cells with the same
+ // error indicator). That would mean that
+ // all cells will be flagged for
+ // refinement or coarsening, but some will
+ // be flagged for both, namely those for
+ // which the indicator equals the
+ // thresholds. This is forbidden, however.
+ //
+ // In some rare cases with very few cells
+ // we also could get integer round off
+ // errors and get problems with
+ // the top and bottom fractions.
+ //
+ // In these case we arbitrarily reduce the
+ // bottom threshold by one permille below
+ // the top threshold
+ //
+ // Finally, in some cases
+ // (especially involving symmetric
+ // solutions) there are many cells
+ // with the same error indicator
+ // values. if there are many with
+ // indicator equal to the top
+ // threshold, no refinement will
+ // take place below; to avoid this
+ // case, we also lower the top
+ // threshold if it equals the
+ // largest indicator and the
+ // top_fraction!=1
+ const double max_criterion = *(criteria_sorted.begin()),
+ min_criterion = *(criteria_sorted.end() - 1);
+
+ if ((top_threshold == max_criterion) && (top_fraction != 1))
+ top_threshold *= 0.999;
+
+ if (bottom_threshold >= top_threshold)
+ bottom_threshold = 0.999 * top_threshold;
+
+ // actually flag cells
+ if (top_threshold < max_criterion)
+ GridRefinement::refine(tria, criteria, top_threshold, refine_cells);
+
+ if (bottom_threshold > min_criterion)
+ GridRefinement::coarsen(tria, criteria, bottom_threshold);
+ }
+} // namespace
+
+
template <int dim, typename Number, int spacedim>
void
const Vector<Number> & criteria,
const double top_fraction,
const double bottom_fraction,
- const unsigned int max_n_cells)
+ const unsigned int max_n_cells,
+ const VectorTools::NormType norm_type)
{
- // correct number of cells is
- // checked in @p{refine}
+ // correct number of cells is checked in @p{refine}
Assert((top_fraction >= 0) && (top_fraction <= 1),
ExcInvalidParameterValue());
Assert((bottom_fraction >= 0) && (bottom_fraction <= 1),
ExcInvalidParameterValue());
Assert(criteria.is_non_negative(), ExcNegativeCriteria());
- // let tmp be the cellwise square of the
- // error, which is what we have to sum
- // up and compare with
- // @p{fraction_of_error*total_error}.
- Vector<Number> tmp;
- tmp = criteria;
- const double total_error = tmp.l1_norm();
-
- // sort the largest criteria to the
- // beginning of the vector
- std::sort(tmp.begin(), tmp.end(), std::greater<double>());
-
- // compute thresholds
- typename Vector<Number>::const_iterator pp = tmp.begin();
- for (double sum = 0; (sum < top_fraction * total_error) && (pp != tmp.end());
- ++pp)
- sum += *pp;
- double top_threshold = (pp != tmp.begin() ? (*pp + *(pp - 1)) / 2 : *pp);
- typename Vector<Number>::const_iterator qq = (tmp.end() - 1);
- for (double sum = 0;
- (sum < bottom_fraction * total_error) && (qq != tmp.begin() - 1);
- --qq)
- sum += *qq;
- double bottom_threshold = (qq != (tmp.end() - 1) ? (*qq + *(qq + 1)) / 2 : 0);
-
- // we now have an idea how many cells we
- // are going to refine and coarsen. we use
- // this information to see whether we are
- // over the limit and if so use a function
- // that knows how to deal with this
- // situation
-
- // note, that at this point, we have no
- // information about anisotropically refined
- // cells, thus use the situation of purely
- // isotropic refinement as guess for a mixed
- // refinemnt as well.
- {
- const unsigned int refine_cells = pp - tmp.begin(),
- coarsen_cells = tmp.end() - 1 - qq;
-
- if (static_cast<unsigned int>(
- tria.n_active_cells() +
- refine_cells * (GeometryInfo<dim>::max_children_per_cell - 1) -
- (coarsen_cells * (GeometryInfo<dim>::max_children_per_cell - 1) /
- GeometryInfo<dim>::max_children_per_cell)) > max_n_cells)
- {
- refine_and_coarsen_fixed_number(tria,
- criteria,
- 1. * refine_cells / criteria.size(),
- 1. * coarsen_cells / criteria.size(),
- max_n_cells);
- return;
- }
- }
-
+ switch (norm_type)
+ {
+ case VectorTools::NormType::L1_norm:
+ // evaluate norms on subsets and compare them as
+ // c_0 + c_1 + ... < fraction * l1-norm(c)
+ refine_and_coarsen_fixed_fraction_via_l1_norm(
+ tria, criteria, top_fraction, bottom_fraction, max_n_cells);
+ break;
+
+ case VectorTools::NormType::L2_norm:
+ {
+ // we do not want to evaluate norms on subsets as:
+ // sqrt(c_0^2 + c_1^2 + ...) < fraction * l2-norm(c)
+ // instead take the square of both sides of the equation
+ // and evaluate:
+ // c_0^2 + c_1^2 + ... < fraction^2 * l1-norm(c.c)
+ // we adjust all parameters accordingly
+ Vector<Number> criteria_squared(criteria.size());
+ std::transform(criteria.begin(),
+ criteria.end(),
+ criteria_squared.begin(),
+ [](Number c) { return c * c; });
+
+ refine_and_coarsen_fixed_fraction_via_l1_norm(tria,
+ criteria_squared,
+ top_fraction *
+ top_fraction,
+ bottom_fraction *
+ bottom_fraction,
+ max_n_cells);
+ }
+ break;
- // in some rare cases it may happen that
- // both thresholds are the same (e.g. if
- // there are many cells with the same
- // error indicator). That would mean that
- // all cells will be flagged for
- // refinement or coarsening, but some will
- // be flagged for both, namely those for
- // which the indicator equals the
- // thresholds. This is forbidden, however.
- //
- // In some rare cases with very few cells
- // we also could get integer round off
- // errors and get problems with
- // the top and bottom fractions.
- //
- // In these case we arbitrarily reduce the
- // bottom threshold by one permille below
- // the top threshold
- //
- // Finally, in some cases
- // (especially involving symmetric
- // solutions) there are many cells
- // with the same error indicator
- // values. if there are many with
- // indicator equal to the top
- // threshold, no refinement will
- // take place below; to avoid this
- // case, we also lower the top
- // threshold if it equals the
- // largest indicator and the
- // top_fraction!=1
- const auto minmax_element =
- std::minmax_element(criteria.begin(), criteria.end());
- if ((top_threshold == *minmax_element.second) && (top_fraction != 1))
- top_threshold *= 0.999;
-
- if (bottom_threshold >= top_threshold)
- bottom_threshold = 0.999 * top_threshold;
-
- // actually flag cells
- if (top_threshold < *minmax_element.second)
- refine(tria, criteria, top_threshold, pp - tmp.begin());
-
- if (bottom_threshold > *minmax_element.first)
- coarsen(tria, criteria, bottom_threshold);
+ default:
+ Assert(false, ExcNotImplemented());
+ break;
+ }
}
const dealii::Vector<S> &,
const double,
const double,
- const unsigned int);
+ const unsigned int,
+ const VectorTools::NormType);
template void GridRefinement::
refine_and_coarsen_optimize<deal_II_dimension, S, deal_II_dimension>(
const dealii::Vector<S> &,
const double,
const double,
- const unsigned int);
+ const unsigned int,
+ const VectorTools::NormType);
template void GridRefinement::
refine_and_coarsen_optimize<deal_II_dimension, S, deal_II_dimension + 1>(
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Verify fixed fraction algorithm with l1-norm and l2-norm
+// Equidistant indicators
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <limits>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+count_flags(const Triangulation<dim> &tria)
+{
+ unsigned int n_refine = 0, n_coarsen = 0;
+ for (const auto &cell : tria.active_cell_iterators())
+ {
+ if (cell->refine_flag_set())
+ ++n_refine;
+ else if (cell->coarsen_flag_set())
+ ++n_coarsen;
+ }
+
+ deallog << "n_refine_flags: " << n_refine
+ << ", n_coarsen_flags: " << n_coarsen;
+}
+
+
+
+template <int dim>
+void
+test()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(2);
+
+ Vector<float> indicator(tria.n_active_cells());
+ // assign each cell a globally unique cellid
+ for (const auto &cell : tria.active_cell_iterators())
+ {
+ 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;
+ }
+
+ deallog << "l1-norm: ";
+ GridRefinement::refine_and_coarsen_fixed_fraction(
+ tria,
+ indicator,
+ 0.3,
+ 0.3,
+ std::numeric_limits<unsigned int>::max(),
+ VectorTools::NormType::L1_norm);
+ count_flags(tria);
+ deallog << std::endl;
+
+ // reset refinement flags
+ for (const auto &cell : tria.active_cell_iterators())
+ {
+ cell->clear_refine_flag();
+ cell->clear_coarsen_flag();
+ }
+
+ deallog << "l2-norm: ";
+ GridRefinement::refine_and_coarsen_fixed_fraction(
+ tria,
+ indicator,
+ 0.3,
+ 0.3,
+ std::numeric_limits<unsigned int>::max(),
+ VectorTools::NormType::L2_norm);
+ count_flags(tria);
+ deallog << std::endl;
+}
+
+
+
+int
+main(int argc, const char *argv[])
+{
+ initlog();
+
+ test<2>();
+}
--- /dev/null
+
+DEAL::l1-norm: n_refine_flags: 3, n_coarsen_flags: 10
+DEAL::l2-norm: n_refine_flags: 1, n_coarsen_flags: 8
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Verify fixed fraction algorithm with l1-norm and l2-norm
+// Random indicators
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <limits>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+count_flags(const Triangulation<dim> &tria)
+{
+ unsigned int n_refine = 0, n_coarsen = 0;
+ for (const auto &cell : tria.active_cell_iterators())
+ {
+ if (cell->refine_flag_set())
+ ++n_refine;
+ else if (cell->coarsen_flag_set())
+ ++n_coarsen;
+ }
+
+ deallog << "n_refine_flags: " << n_refine
+ << ", n_coarsen_flags: " << n_coarsen;
+}
+
+
+
+template <int dim>
+void
+test()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(2);
+
+ Vector<float> indicator(tria.n_active_cells());
+ // assign each cell a globally unique cellid
+ for (const auto &cell : tria.active_cell_iterators())
+ {
+ 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>();
+ }
+
+ deallog << "l1-norm: ";
+ GridRefinement::refine_and_coarsen_fixed_fraction(
+ tria,
+ indicator,
+ 0.3,
+ 0.3,
+ std::numeric_limits<unsigned int>::max(),
+ VectorTools::NormType::L1_norm);
+ count_flags(tria);
+ deallog << std::endl;
+
+ // reset refinement flags
+ for (const auto &cell : tria.active_cell_iterators())
+ {
+ cell->clear_refine_flag();
+ cell->clear_coarsen_flag();
+ }
+
+ deallog << "l2-norm: ";
+ GridRefinement::refine_and_coarsen_fixed_fraction(
+ tria,
+ indicator,
+ 0.3,
+ 0.3,
+ std::numeric_limits<unsigned int>::max(),
+ VectorTools::NormType::L2_norm);
+ count_flags(tria);
+ deallog << std::endl;
+}
+
+
+
+int
+main(int argc, const char *argv[])
+{
+ initlog();
+
+ test<2>();
+}
--- /dev/null
+
+DEAL::l1-norm: n_refine_flags: 4, n_coarsen_flags: 8
+DEAL::l2-norm: n_refine_flags: 1, n_coarsen_flags: 5
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Verify fixed fraction algorithm with l1-norm and l2-norm
+// Equidistant indicators
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/vector.h>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+count_flags(const parallel::distributed::Triangulation<dim> &tria)
+{
+ unsigned int n_refine = 0, n_coarsen = 0;
+ for (const auto &cell : tria.active_cell_iterators())
+ if (cell->is_locally_owned())
+ {
+ if (cell->refine_flag_set())
+ ++n_refine;
+ else if (cell->coarsen_flag_set())
+ ++n_coarsen;
+ }
+
+ const unsigned int n_refine_global =
+ Utilities::MPI::sum(n_refine, MPI_COMM_WORLD),
+ n_coarsen_global =
+ Utilities::MPI::sum(n_coarsen, MPI_COMM_WORLD);
+
+ deallog << "n_refine_flags: " << n_refine_global
+ << ", n_coarsen_flags: " << n_coarsen_global;
+}
+
+
+
+template <int dim>
+void
+test()
+{
+ parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(2);
+
+ 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));
+
+ indicator[cell->active_cell_index()] = fine_cellid + 1;
+ }
+
+ deallog << "l1-norm: ";
+ parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction(
+ tria, indicator, 0.3, 0.3, VectorTools::NormType::L1_norm);
+ count_flags(tria);
+ deallog << std::endl;
+
+ // reset refinement flags
+ for (const auto &cell : tria.active_cell_iterators())
+ if (cell->is_locally_owned())
+ {
+ cell->clear_refine_flag();
+ cell->clear_coarsen_flag();
+ }
+
+ deallog << "l2-norm: ";
+ parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction(
+ tria, indicator, 0.3, 0.3, VectorTools::NormType::L2_norm);
+ count_flags(tria);
+ deallog << std::endl;
+}
+
+
+
+int
+main(int argc, char *argv[])
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+ const unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+ if (myid == 0)
+ {
+ initlog();
+
+ test<2>();
+ }
+ else
+ test<2>();
+}
--- /dev/null
+
+DEAL::l1-norm: n_refine_flags: 2, n_coarsen_flags: 9
+DEAL::l2-norm: n_refine_flags: 1, n_coarsen_flags: 7
--- /dev/null
+
+DEAL::l1-norm: n_refine_flags: 2, n_coarsen_flags: 9
+DEAL::l2-norm: n_refine_flags: 1, n_coarsen_flags: 7
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Verify fixed fraction algorithm with l1-norm and l2-norm
+// Random indicators
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/vector.h>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+count_flags(const parallel::distributed::Triangulation<dim> &tria)
+{
+ unsigned int n_refine = 0, n_coarsen = 0;
+ for (const auto &cell : tria.active_cell_iterators())
+ if (cell->is_locally_owned())
+ {
+ if (cell->refine_flag_set())
+ ++n_refine;
+ else if (cell->coarsen_flag_set())
+ ++n_coarsen;
+ }
+
+ const unsigned int n_refine_global =
+ Utilities::MPI::sum(n_refine, MPI_COMM_WORLD),
+ n_coarsen_global =
+ Utilities::MPI::sum(n_coarsen, MPI_COMM_WORLD);
+
+ deallog << "n_refine_flags: " << n_refine_global
+ << ", n_coarsen_flags: " << n_coarsen_global;
+}
+
+
+
+template <int dim>
+void
+test()
+{
+ parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(2);
+
+ 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));
+
+ Testing::srand(fine_cellid);
+ indicator[cell->active_cell_index()] = random_value<float>();
+ }
+
+ deallog << "l1-norm: ";
+ parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction(
+ tria, indicator, 0.3, 0.3, VectorTools::NormType::L1_norm);
+ count_flags(tria);
+ deallog << std::endl;
+
+ // reset refinement flags
+ for (const auto &cell : tria.active_cell_iterators())
+ if (cell->is_locally_owned())
+ {
+ cell->clear_refine_flag();
+ cell->clear_coarsen_flag();
+ }
+
+ deallog << "l2-norm: ";
+ parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction(
+ tria, indicator, 0.3, 0.3, VectorTools::NormType::L2_norm);
+ count_flags(tria);
+ deallog << std::endl;
+}
+
+
+
+int
+main(int argc, char *argv[])
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+ const unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+ if (myid == 0)
+ {
+ initlog();
+
+ test<2>();
+ }
+ else
+ test<2>();
+}
--- /dev/null
+
+DEAL::l1-norm: n_refine_flags: 4, n_coarsen_flags: 8
+DEAL::l2-norm: n_refine_flags: 1, n_coarsen_flags: 5
--- /dev/null
+
+DEAL::l1-norm: n_refine_flags: 4, n_coarsen_flags: 8
+DEAL::l2-norm: n_refine_flags: 1, n_coarsen_flags: 5
projects / dealii.git / commitdiff