]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add custom constraints 8730/head
authortcclevenger <tcleven@clemson.edu>
Wed, 21 Aug 2019 19:46:02 +0000 (15:46 -0400)
committertcclevenger <tcleven@clemson.edu>
Mon, 16 Sep 2019 19:32:21 +0000 (15:32 -0400)
doc/news/changes/minor/20190913Clevenger [new file with mode: 0644]
include/deal.II/multigrid/mg_constrained_dofs.h
source/lac/affine_constraints.inst.in
source/multigrid/mg_transfer_matrix_free.cc
source/multigrid/mg_transfer_prebuilt.cc
source/multigrid/multigrid.cc
tests/multigrid/transfer_matrix_free_13.cc [new file with mode: 0644]
tests/multigrid/transfer_matrix_free_13.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20190913Clevenger b/doc/news/changes/minor/20190913Clevenger
new file mode 100644 (file)
index 0000000..74c2b41
--- /dev/null
@@ -0,0 +1,5 @@
+New: MGConstrainedDoFs now allows the user to define constraints 
+on each multigrid level to be applied during prolongation in a 
+matrix-free transfer.
+<br>
+(Conrad Clevenger, 2019/09/13)
index eaa5ea9f038dff913c0c222cdef7271b07ccf1b3..614ab15be44bf62afb31c68b90ecbabdb602eee5 100644 (file)
@@ -103,6 +103,24 @@ public:
     const std::set<types::boundary_id> &boundary_ids,
     const ComponentMask &               component_mask = ComponentMask());
 
+  /**
+   * Add user defined constraints to be used on level @p level.
+   *
+   * The user can call this function multiple times and any new,
+   * conflicting constraints will overwrite the previous constraints
+   * for that DoF.
+   *
+   * Before the transfer, the user defined constraints will be distributed
+   * to the source vector, and then any DoF index set using
+   * make_zero_boundary_constraints() will be overwritten with
+   * value zero.
+   *
+   * @note This is currently only implemented for MGTransferMatrixFree.
+   */
+  void
+  add_user_constraints(const unsigned int               level,
+                       const AffineConstraints<double> &constraints_on_level);
+
   /**
    * Fill the internal data structures with information
    * about no normal flux boundary dofs.
@@ -196,6 +214,16 @@ public:
   const AffineConstraints<double> &
   get_level_constraint_matrix(const unsigned int level) const;
 
+  /**
+   * Return the user defined constraint matrix for a given level. These
+   * constraints are set using the function add_user_constraints() and
+   * should not contain constraints for DoF indices set in
+   * make_zero_boundary_constraints() as they will be overwritten during
+   * the transfer.
+   */
+  const AffineConstraints<double> &
+  get_user_constraint_matrix(const unsigned int level) const;
+
 private:
   /**
    * The indices of boundary dofs for each level.
@@ -213,6 +241,11 @@ private:
    * periodic boundary conditions for each level .
    */
   std::vector<AffineConstraints<double>> level_constraints;
+
+  /**
+   * Constraint matrices defined by user.
+   */
+  std::vector<AffineConstraints<double>> user_constraints;
 };
 
 
@@ -223,12 +256,14 @@ MGConstrainedDoFs::initialize(const DoFHandler<dim, spacedim> &dof)
   boundary_indices.clear();
   refinement_edge_indices.clear();
   level_constraints.clear();
+  user_constraints.clear();
 
   const unsigned int nlevels = dof.get_triangulation().n_global_levels();
 
   // At this point level_constraint and refinement_edge_indices are empty.
-  level_constraints.resize(nlevels);
   refinement_edge_indices.resize(nlevels);
+  level_constraints.resize(nlevels);
+  user_constraints.resize(nlevels);
   for (unsigned int l = 0; l < nlevels; ++l)
     {
       IndexSet relevant_dofs;
@@ -399,11 +434,31 @@ MGConstrainedDoFs::make_no_normal_flux_constraints(
 }
 
 
+inline void
+MGConstrainedDoFs::add_user_constraints(
+  const unsigned int               level,
+  const AffineConstraints<double> &constraints_on_level)
+{
+  AssertIndexRange(level, user_constraints.size());
+
+  // Get the relevant DoFs from level_constraints if
+  // the user constraint matrix has not been initialized
+  if (user_constraints[level].get_local_lines().size() == 0)
+    user_constraints[level].reinit(level_constraints[level].get_local_lines());
+
+  user_constraints[level].merge(
+    constraints_on_level,
+    AffineConstraints<double>::MergeConflictBehavior::right_object_wins);
+  user_constraints[level].close();
+}
+
+
 inline void
 MGConstrainedDoFs::clear()
 {
   boundary_indices.clear();
   refinement_edge_indices.clear();
+  user_constraints.clear();
 }
 
 
@@ -487,6 +542,15 @@ MGConstrainedDoFs::get_level_constraint_matrix(const unsigned int level) const
 
 
 
+inline const AffineConstraints<double> &
+MGConstrainedDoFs::get_user_constraint_matrix(const unsigned int level) const
+{
+  AssertIndexRange(level, user_constraints.size());
+  return user_constraints[level];
+}
+
+
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
index 4523cad087ba1cab7335d61d5c5049bc86e105c6..e3877ee611d95fe0f4e23180eb1d74584e47e555 100644 (file)
@@ -344,4 +344,8 @@ for (T : DEAL_II_VEC_TEMPLATES)
       T<float> &) const;
     template void dealii::AffineConstraints<double>::distribute<T<float>>(
       T<float> &) const;
+
+    template void dealii::AffineConstraints<double>::distribute<
+      LinearAlgebra::distributed::T<float>>(
+      LinearAlgebra::distributed::T<float> &) const;
   }
index 49557a66ddc0d4a6fdbb300638ed2c0c126c0928..769c532e8c60f7c179c1ea6c742f83a7a5a94201 100644 (file)
@@ -355,6 +355,30 @@ MGTransferMatrixFree<dim, Number>::do_prolongate_add(
     Utilities::fixed_power<dim>(n_child_dofs_1d);
   constexpr unsigned int three_to_dim = Utilities::pow(3, dim);
 
+  // If we have user defined MG constraints, we must create
+  // a non-const, ghosted version of the source vector to distribute
+  // constraints.
+  const LinearAlgebra::distributed::Vector<Number> *to_use = &src;
+  LinearAlgebra::distributed::Vector<Number>        new_src;
+  if (this->mg_constrained_dofs != nullptr &&
+      this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1)
+          .get_local_lines()
+          .size() > 0)
+    {
+      LinearAlgebra::distributed::Vector<Number> copy_src(src);
+
+      // Distribute any user defined constraints
+      this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1)
+        .distribute(copy_src);
+
+      // Re-initialize new ghosted vector with correct constraints
+      new_src.reinit(copy_src);
+      new_src = copy_src;
+      new_src.update_ghost_values();
+
+      to_use = &new_src;
+    }
+
   for (unsigned int cell = 0; cell < n_owned_level_cells[to_level - 1];
        cell += vec_size)
     {
@@ -382,7 +406,7 @@ MGTransferMatrixFree<dim, Number>::do_prolongate_add(
               for (unsigned int k = 0; k < (dim > 2 ? degree_size : 1); ++k)
                 for (unsigned int j = 0; j < (dim > 1 ? degree_size : 1); ++j)
                   for (unsigned int i = 0; i < degree_size; ++i, ++m)
-                    evaluation_data[m][v] = src.local_element(
+                    evaluation_data[m][v] = to_use->local_element(
                       indices[c * n_scalar_cell_dofs +
                               k * n_child_dofs_1d * n_child_dofs_1d +
                               j * n_child_dofs_1d + i]);
index a65380efffbddaf13adcc63067bd7217f4449012..0eceb6b356e6257405185a1ae2d26b1d5387b2cb 100644 (file)
@@ -106,6 +106,14 @@ MGTransferPrebuilt<VectorType>::prolongate(const unsigned int to_level,
   Assert((to_level >= 1) && (to_level <= prolongation_matrices.size()),
          ExcIndexRange(to_level, 1, prolongation_matrices.size() + 1));
 
+#ifdef DEBUG
+  if (this->mg_constrained_dofs != nullptr)
+    Assert(this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1)
+               .get_local_lines()
+               .size() == 0,
+           ExcNotImplemented());
+#endif
+
   prolongation_matrices[to_level - 1]->vmult(dst, src);
 }
 
index ff37ab9eb503772ecb540d84c696623f15943f25..310fc0da1356755a267171bd8d832f01db152590 100644 (file)
@@ -94,6 +94,14 @@ MGTransferBlock<number>::prolongate(const unsigned int         to_level,
   Assert(dst.n_blocks() == this->n_mg_blocks,
          ExcDimensionMismatch(dst.n_blocks(), this->n_mg_blocks));
 
+#ifdef DEBUG
+  if (this->mg_constrained_dofs != nullptr)
+    Assert(this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1)
+               .get_local_lines()
+               .size() == 0,
+           ExcNotImplemented());
+#endif
+
   // Multiplicate with prolongation
   // matrix, but only those blocks
   // selected.
@@ -280,6 +288,14 @@ MGTransferBlockSelect<number>::prolongate(const unsigned int    to_level,
   Assert((to_level >= 1) && (to_level <= prolongation_matrices.size()),
          ExcIndexRange(to_level, 1, prolongation_matrices.size() + 1));
 
+#ifdef DEBUG
+  if (this->mg_constrained_dofs != nullptr)
+    Assert(this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1)
+               .get_local_lines()
+               .size() == 0,
+           ExcNotImplemented());
+#endif
+
   prolongation_matrices[to_level - 1]
     ->block(selected_block, selected_block)
     .vmult(dst, src);
diff --git a/tests/multigrid/transfer_matrix_free_13.cc b/tests/multigrid/transfer_matrix_free_13.cc
new file mode 100644 (file)
index 0000000..dd1e39b
--- /dev/null
@@ -0,0 +1,98 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 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 MGTransferMatrixFree with custom user constraints
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <deal.II/multigrid/mg_transfer.h>
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+check()
+{
+  Triangulation<dim> tr(Triangulation<dim>::limit_level_difference_at_vertices);
+  GridGenerator::hyper_cube(tr, 0, 1, false);
+
+  typename Triangulation<dim>::active_cell_iterator cell = tr.begin_active();
+  cell->face(0)->set_all_boundary_ids(1);
+
+  tr.refine_global(1);
+
+  FE_Q<dim> fe = FE_Q<dim>(1);
+
+  DoFHandler<dim> mgdof(tr);
+  mgdof.distribute_dofs(fe);
+  mgdof.distribute_mg_dofs();
+
+  MGConstrainedDoFs mg_constrained_dofs;
+  mg_constrained_dofs.initialize(mgdof);
+
+  IndexSet relevant_dofs;
+  DoFTools::extract_locally_relevant_level_dofs(mgdof, 0, relevant_dofs);
+  AffineConstraints<double> user_constraints;
+  user_constraints.reinit(relevant_dofs);
+
+  typename DoFHandler<dim>::level_face_iterator face0 = mgdof.begin(0)->face(0);
+  std::vector<types::global_dof_index>          face_dofs(fe.dofs_per_face);
+  face0->get_mg_dof_indices(0, face_dofs);
+  for (unsigned int i = 0; i < face_dofs.size(); ++i)
+    {
+      if (user_constraints.can_store_line(face_dofs[i]))
+        {
+          user_constraints.add_line(face_dofs[i]);
+          user_constraints.set_inhomogeneity(face_dofs[i], 5.0);
+        }
+    }
+  user_constraints.close();
+  mg_constrained_dofs.add_user_constraints(0, user_constraints);
+
+  // build matrix-free transfer
+  MGTransferMatrixFree<dim, double> transfer(mg_constrained_dofs);
+  transfer.build(mgdof);
+
+  LinearAlgebra::distributed::Vector<double> src_level_0(mgdof.n_dofs(0));
+  for (unsigned int i = 0; i < mgdof.n_dofs(0); ++i)
+    deallog << src_level_0(i) << " ";
+  deallog << std::endl;
+
+  LinearAlgebra::distributed::Vector<double> dst_level_1(mgdof.n_dofs(1));
+  transfer.prolongate(1, dst_level_1, src_level_0);
+  for (unsigned int i = 0; i < mgdof.n_dofs(1); ++i)
+    deallog << dst_level_1(i) << " ";
+  deallog << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+
+  check<2>();
+  deallog << std::endl;
+  check<3>();
+  deallog << std::endl;
+}
diff --git a/tests/multigrid/transfer_matrix_free_13.output b/tests/multigrid/transfer_matrix_free_13.output
new file mode 100644 (file)
index 0000000..f51ed58
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL::0.00000 0.00000 0.00000 0.00000 
+DEAL::5.00000 2.50000 5.00000 2.50000 0.00000 0.00000 5.00000 2.50000 0.00000 
+DEAL::
+DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 
+DEAL::5.00000 2.50000 5.00000 2.50000 5.00000 2.50000 5.00000 2.50000 0.00000 0.00000 0.00000 0.00000 5.00000 2.50000 5.00000 2.50000 0.00000 0.00000 5.00000 2.50000 5.00000 2.50000 0.00000 0.00000 5.00000 2.50000 0.00000 
+DEAL::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.