* Constructor.
*
* @param[in] dof_handler The DoFHandler on which all the operations
- * will happen. At the time when this constructor is call, the
- * DoFHandler still points to the Triangulation before the refinement in
- * question happens.
+ * will happen. This constructor must be called before the underlying
+ * Triangulation is coarsened/refined.
*/
FieldTransfer(const DoFHandler<dim, spacedim> &dof_handler);
- /**
- * Destructor.
- */
- ~FieldTransfer() = default;
-
/**
* Prepare the current object for coarsening and refinement.
* @param[in] in The vector that will be interpolated
/**
* Interpolate the data previously stored in this object before the mesh
- * was refined or coarsenend onto the current set of cells. @p new_value
+ * was refined or coarsened onto the current set of cells. @p new_value
* is the value associated to the new degrees of freedom that where
* created during the element activation. @p affine_constraints is the
- * AffinieConstraints after refinement.
+ * AffineConstraints after refinement.
*/
void
interpolate(const Number & new_value,
const DoFHandler<dim, spacedim> &dof)
: dof_handler(dof)
{
+ // When coarsening, we want to mimic the behavior of SolutionTransfer
+ // and interpolate from child cells to parent. Define this strategy here
+ // since it is not readily available
+ const auto coarsening_strategy =
+ [this](
+ const typename dealii::Triangulation<dim, spacedim>::cell_iterator
+ & parent,
+ const std::vector<Vector<Number>> &children_values) {
+ // get the equivalent DoFCellAccessor
+ typename DoFHandler<dim, spacedim>::cell_iterator dof_cell_iterator(
+ &dof_handler.get_triangulation(),
+ parent->level(),
+ parent->index(),
+ &dof_handler);
+
+ int fe_index = 0;
+ if (dof_handler.has_hp_capabilities())
+ fe_index = dealii::internal::hp::DoFHandlerImplementation::
+ dominated_future_fe_on_children<dim, spacedim>(
+ dof_cell_iterator);
+
+ const auto &fe = dof_handler.get_fe(fe_index);
+ Assert(fe.n_dofs_per_cell() > 0,
+ ExcMessage(
+ "Cannot coarsen onto a FiniteElement with no DoFs."));
+ AssertDimension(dof_cell_iterator->n_children(),
+ children_values.size());
+
+ const auto child_iterators = dof_cell_iterator->child_iterators();
+ const unsigned int n_children_with_fe_nothing =
+ std::count_if(child_iterators.begin(),
+ child_iterators.end(),
+ [](const auto &child_cell) {
+ return child_cell->get_fe().n_dofs_per_cell() ==
+ 0;
+ });
+
+ Assert(
+ n_children_with_fe_nothing == 0 ||
+ n_children_with_fe_nothing == dof_cell_iterator->n_children(),
+ ExcMessage(
+ "Coarsening is only supported for parent cells where either all"
+ " or none of the child cells are FE_Nothing."));
+
+ // in case all children are FE_Nothing there is nothing to
+ // interpolate and we just return the first entry from the children
+ // values (containing invalid entries)
+ if (n_children_with_fe_nothing == dof_cell_iterator->n_children())
+ {
+ return children_values[0];
+ }
+
+ const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
+ Vector<Number> tmp(dofs_per_cell);
+ Vector<Number> interpolated_values(dofs_per_cell);
+
+ // Otherwise, perform the actual interpolation here. Due to the
+ // assert above, we know that all child cells have data to
+ // interpolate.
+ for (unsigned int child = 0;
+ child < dof_cell_iterator->n_children();
+ ++child)
+ {
+ // interpolate the previously stored values on a child to the
+ // mother cell
+ fe.get_restriction_matrix(child,
+ dof_cell_iterator->refinement_case())
+ .vmult(tmp, children_values[child]);
+
+ // and add up or set them in the output vector
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ if (fe.restriction_is_additive(i))
+ interpolated_values(i) += tmp(i);
+ else if (tmp(i) != Number())
+ interpolated_values(i) = tmp(i);
+ }
+
+ return interpolated_values;
+ };
+
cell_data_transfer = std::make_unique<
CellDataTransfer<dim, spacedim, std::vector<Vector<Number>>>>(
dynamic_cast<
dealii::parallel::distributed::Triangulation<dim, spacedim> &>(
const_cast<dealii::Triangulation<dim, spacedim> &>(
- dof_handler.get_triangulation())));
+ dof_handler.get_triangulation())),
+ false,
+ &dealii::AdaptationStrategies::Refinement::
+ preserve<dim, spacedim, Vector<Number>>,
+ coarsening_strategy);
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 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.
+//
+// ---------------------------------------------------------------------
+
+// Check that we can use FieldTransfer when coarsening cells with non-uniform
+// data
+
+#include <deal.II/distributed/field_transfer.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+class TestFunction : public Function<2>
+{
+public:
+ double
+ value(const dealii::Point<2> &p, const unsigned int component) const override
+ {
+ return std::exp(10 - 10 * p[1]);
+ }
+};
+
+void
+output(const DoFHandler<2> & dof_handler,
+ const LinearAlgebra::distributed::Vector<double> &vector,
+ const std::string & file_name)
+{
+ DataOut<2> data_out;
+ data_out.attach_dof_handler(dof_handler);
+ data_out.add_data_vector(vector, "solution");
+ data_out.build_patches();
+ data_out.write_vtu_in_parallel(file_name, MPI_COMM_WORLD);
+}
+
+int
+main(int argc, char *argv[])
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+ MPILogInitAll mpi_log;
+
+ parallel::distributed::Triangulation<2> triangulation(MPI_COMM_WORLD);
+ GridGenerator::hyper_cube(triangulation);
+ triangulation.refine_global(3);
+
+ hp::FECollection<2> fe_collection;
+ fe_collection.push_back(FE_Q<2>(1));
+ fe_collection.push_back(FE_Nothing<2>());
+
+ DoFHandler<2> dof_handler(triangulation);
+
+ // Assign FE_Nothing to half of the domain
+ for (auto cell : dof_handler.active_cell_iterators())
+ {
+ if (cell->is_locally_owned())
+ {
+ if (cell->center()[0] < 0.5)
+ {
+ cell->set_active_fe_index(0);
+ }
+ else
+ {
+ cell->set_active_fe_index(1);
+ }
+ }
+ }
+
+ dof_handler.distribute_dofs(fe_collection);
+
+ // Initialize solution
+ auto locally_relevant_dofs =
+ DoFTools::extract_locally_relevant_dofs(dof_handler);
+ LinearAlgebra::distributed::Vector<double> solution(
+ dof_handler.locally_owned_dofs(), locally_relevant_dofs, MPI_COMM_WORLD);
+
+ // Initialize with a non-constant function
+ VectorTools::interpolate(dof_handler, TestFunction(), solution);
+
+ {
+ if (false)
+ output(dof_handler, solution, "result_0.vtu");
+
+ deallog << "solution before: " << std::endl;
+ std::stringstream ss;
+ solution.print(ss);
+ deallog << ss.str() << std::endl;
+ }
+
+ parallel::distributed::experimental::
+ FieldTransfer<2, LinearAlgebra::distributed::Vector<double>>
+ field_transfer(dof_handler);
+ // Assign FE_Q to all the cells
+ for (auto cell : dof_handler.active_cell_iterators())
+ {
+ if (cell->is_locally_owned())
+ {
+ cell->set_future_fe_index(0);
+ }
+ }
+
+ // coarsen some cells
+ for (auto cell : dof_handler.active_cell_iterators())
+ {
+ if (cell->is_locally_owned())
+ {
+ if (cell->center()[0] < 0.25 && cell->center()[1] < 0.25)
+ {
+ cell->set_coarsen_flag();
+ }
+ }
+ }
+
+ triangulation.prepare_coarsening_and_refinement();
+
+ field_transfer.prepare_for_coarsening_and_refinement(solution, 1);
+
+ triangulation.execute_coarsening_and_refinement();
+
+ dof_handler.distribute_dofs(fe_collection);
+
+ locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler);
+ LinearAlgebra::distributed::Vector<double> new_solution(
+ dof_handler.locally_owned_dofs(), locally_relevant_dofs, MPI_COMM_WORLD);
+
+ AffineConstraints<double> affine_constraints;
+ DoFTools::make_hanging_node_constraints(dof_handler, affine_constraints);
+ affine_constraints.close();
+
+ const double new_value = 11.;
+
+ field_transfer.interpolate(new_value, affine_constraints, new_solution);
+
+ affine_constraints.distribute(new_solution);
+
+ {
+ if (false)
+ output(dof_handler, new_solution, "result_1.vtu");
+
+ deallog << "solution after coarsening and activation: " << std::endl;
+ std::stringstream ss;
+ new_solution.print(ss);
+ deallog << ss.str() << std::endl;
+ }
+}
--- /dev/null
+
+DEAL:0::solution before:
+DEAL:0::Process #0
+Local range: [0, 25), global size: 45
+Vector data:
+2.203e+04 2.203e+04 6.311e+03 6.311e+03 2.203e+04 6.311e+03 1.808e+03 1.808e+03 1.808e+03 2.203e+04 6.311e+03 2.203e+04 6.311e+03 1.808e+03 1.808e+03 5.180e+02 5.180e+02 5.180e+02 1.484e+02 1.484e+02 1.484e+02 5.180e+02 5.180e+02 1.484e+02 1.484e+02
+
+DEAL:0::solution after coarsening and activation:
+DEAL:0::Process #0
+Local range: [0, 42), global size: 78
+Vector data:
+2.203e+04 2.203e+04 1.808e+03 1.808e+03 2.203e+04 1.192e+04 6.311e+03 2.203e+04 6.311e+03 1.808e+03 1.808e+03 1.808e+03 5.180e+02 5.180e+02 5.180e+02 1.484e+02 1.484e+02 1.484e+02 5.180e+02 5.180e+02 1.484e+02 1.484e+02 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01
+
+
+DEAL:1::solution before:
+DEAL:1::Process #1
+Local range: [25, 45), global size: 45
+Vector data:
+4.252e+01 4.252e+01 4.252e+01 1.218e+01 1.218e+01 1.218e+01 4.252e+01 4.252e+01 1.218e+01 1.218e+01 3.490e+00 3.490e+00 3.490e+00 1.000e+00 1.000e+00 1.000e+00 3.490e+00 3.490e+00 1.000e+00 1.000e+00
+
+DEAL:1::solution after coarsening and activation:
+DEAL:1::Process #1
+Local range: [42, 78), global size: 78
+Vector data:
+4.252e+01 4.252e+01 4.252e+01 1.218e+01 1.218e+01 1.218e+01 4.252e+01 4.252e+01 1.218e+01 1.218e+01 3.490e+00 3.490e+00 3.490e+00 1.000e+00 1.000e+00 1.000e+00 3.490e+00 3.490e+00 1.000e+00 1.000e+00 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01
+
+
projects / dealii.git / commitdiff