]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added ConvolutionCoupling
authorLuca Heltai <luca.heltai@sissa.it>
Thu, 11 Apr 2019 13:53:05 +0000 (15:53 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Sat, 11 May 2019 16:15:24 +0000 (18:15 +0200)
doc/news/changes/minor/20190419LucaHeltai [new file with mode: 0644]
include/deal.II/non_matching/coupling.h
source/non_matching/coupling.cc
source/non_matching/coupling.inst.in
tests/non_matching/coupling_08.cc [new file with mode: 0644]
tests/non_matching/coupling_08.with_umfpack=true.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20190419LucaHeltai b/doc/news/changes/minor/20190419LucaHeltai
new file mode 100644 (file)
index 0000000..5b726ac
--- /dev/null
@@ -0,0 +1,5 @@
+New: Two new methods NonMatching::compute_coupling_sparsity_pattern() and 
+NonMatching::compute_coupling_mass_matrix() allow to construct the coupling between arbitrary grids,
+using a convolution kernel. 
+<br>
+(Luca Heltai, 2019/04/19)
index 056c0d67ed82ef212f51607330ba9fdec0ba5378..1894e625fc01e3bec3874670e3a042c6480dc928 100644 (file)
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/base/function_lib.h>
 #include <deal.II/base/quadrature.h>
 
 #include <deal.II/dofs/dof_handler.h>
 
 #include <deal.II/fe/component_mask.h>
 #include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/mapping_q1.h>
 
 #include <deal.II/grid/grid_tools_cache.h>
@@ -34,12 +36,11 @@ DEAL_II_NAMESPACE_OPEN
 
 /**
  * A namespace for functions offering tools to handle two meshes with no
- * alignment requirements, but where one of the meshes is embedded
- * inside the other in the real-space.
+ * alignment requirements.
  *
  * Typically these functions allow for computations on the real-space
  * intersection between the two meshes e.g. surface integrals and
- * construction of mass matrices.
+ * construction of coupling matrices.
  */
 namespace NonMatching
 {
@@ -221,6 +222,134 @@ namespace NonMatching
     const ComponentMask &          immersed_comps = ComponentMask(),
     const Mapping<dim1, spacedim> &immersed_mapping =
       StaticMappingQ1<dim1, spacedim>::mapping);
+
+  /**
+   * Create a coupling sparsity pattern for non-matching independent grids,
+   * using a convolution kernel with compact support of radius epsilon.
+   *
+   * Given two non-matching triangulations, representing the domains $\Omega^0$
+   * and $\Omega^1$, both embedded in $R^d$, and two finite element spaces
+   * $V^0(\Omega^0) = \text{span}\{v_i\}_{i=0}^n$ and $V^1(\Omega^1) =
+   * \text{span}\{w_\alpha\}_{\alpha=0}^m$, compute the sparsity pattern that
+   * would be necessary to assemble the matrix
+   *
+   * \f[
+   * M_{i\alpha} \dealcoloneq \int_{\Omega^0} \int_{\Omega^1}
+   * v_i(x) K^{\epsilon}(x-y) w_\alpha(y) dx \ dy,
+   * \quad i \in [0,n), \alpha \in [0,m),
+   * \f]
+   *
+   * where $V^0(\Omega^0)$ is the finite element space associated with the
+   * `dh0` passed to this function (or part of it, if specified in
+   * `comps0`), while $V^1(\Omega^1)$ is the finite element space associated
+   * with the `dh1` passed to this function (or part of it, if specified
+   * in `comps1`), and $K^\epsilon$ is a function with compact support included
+   * in a ball of radius $\epsilon$, derived from CutOffFunctionBase.
+   *
+   * The `comps0` and `comps1` masks are assumed to be ordered in
+   * the same way: the first component of `comps0` will couple with the
+   * first component of `comps1`, the second with the second, and so
+   * on. If one of the two masks has more non-zero than the other, then the
+   * excess components will be ignored.
+   *
+   * For both spaces, it is possible to specify a custom Mapping, which
+   * defaults to StaticMappingQ1 for both.
+   *
+   * This function will also work in parallel, provided that at least one of the
+   * triangulations is of type parallel::shared::Triangulation<dim1,spacedim>.
+   * An exception is thrown if both triagnulations are of type
+   * parallel::distributed::Triangulation<dim1,spacedim>.
+   *
+   * This function assumes that the convolution has support contained in a box
+   * of radius @p epsilon. If epsilon is set to zero, then we assume that the
+   * kernel is the Dirac delta distribution, and the call is forwarded to the
+   * method in this namespace with the same name, that does not take an epsilon
+   * as input. In this case, more restrictive conditions are required on the
+   * two spaces. See the documentation of the other
+   * create_coupling_sparsity_pattern() function.
+   *
+   * @author Luca Heltai, 2019.
+   */
+  template <int dim0,
+            int dim1,
+            int spacedim,
+            typename Sparsity,
+            typename Number = double>
+  void
+  create_coupling_sparsity_pattern(
+    const double &                          epsilon,
+    const GridTools::Cache<dim0, spacedim> &cache0,
+    const GridTools::Cache<dim1, spacedim> &cache1,
+    const DoFHandler<dim0, spacedim> &      dh0,
+    const DoFHandler<dim1, spacedim> &      dh1,
+    Sparsity &                              sparsity,
+    const AffineConstraints<Number> &constraints0 = AffineConstraints<Number>(),
+    const ComponentMask &            comps0       = ComponentMask(),
+    const ComponentMask &            comps1       = ComponentMask());
+
+  /**
+   * Create a coupling  mass matrix for non-matching independent grids,
+   * using a convolution kernel with compact support.
+   *
+   * Given two non-matching triangulations, representing the domains
+   * $\Omega^0$ and $\Omega^1$, both embedded in $R^d$, and two finite element
+   * spaces $V^0(\Omega^0) = \text{span}\{v_i\}_{i=0}^n$ and $V^1(\Omega^1) =
+   * \text{span}\{w_\alpha\}_{\alpha=0}^m$, compute the matrix
+   *
+   * \f[
+   * M_{i\alpha} \dealcoloneq \int_{\Omega^0} \int_{\Omega^1}
+   * v_i(x) K^{\epsilon}(x-y) w_\alpha(y) dx \ dy,
+   * \quad i \in [0,n), \alpha \in [0,m),
+   * \f]
+   *
+   * where $V^0(\Omega^0)$ is the finite element space associated with the
+   * `dh0` passed to this function (or part of it, if specified in
+   * `comps0`), while $V^1(\Omega^1)$ is the finite element space associated
+   * with the `dh1` passed to this function (or part of it, if specified
+   * in `comps1`), and $K^\epsilon$ is a function with compact support included,
+   * in a ball of radius $\epsilon$, derived from CutOffFunctionBase.
+   *
+   * The corresponding sparsity patterns can be computed by calling the
+   * make_coupling_sparsity_pattern() function.
+   *
+   * The `comps0` and `comps1` masks are assumed to be ordered in
+   * the same way: the first component of `comps0` will couple with the
+   * first component of `comps1`, the second with the second, and so
+   * on. If one of the two masks has more non-zero than the other, then the
+   * excess components will be ignored.
+   *
+   * For both spaces, it is possible to specify a custom Mapping, which
+   * defaults to StaticMappingQ1 for both.
+   *
+   * This function will also work in parallel, provided that one of the two
+   * triangulations is of type parallel::shared::Triangulation<dim1,spacedim>.
+   * An exception is thrown if both triangulations are of type
+   * parallel::distributed::Triangulation<dim1,spacedim>.
+   *
+   * The parameter @p epsilon is used to set the size of the cut-off function
+   * used to compute the convolution. If epsilon is set to zero, then we assume
+   * that the kernel is the Dirac delta distribution, and the call is forwarded
+   * to the method in this namespace with the same name, that does not take an
+   * epsilon as input.
+   *
+   * @author Luca Heltai, 2019.
+   */
+  template <int dim0, int dim1, int spacedim, typename Matrix>
+  void
+  create_coupling_mass_matrix(
+    Functions::CutOffFunctionBase<spacedim> &             kernel,
+    const double &                                        epsilon,
+    const GridTools::Cache<dim0, spacedim> &              cache0,
+    const GridTools::Cache<dim1, spacedim> &              cache1,
+    const DoFHandler<dim0, spacedim> &                    dh0,
+    const DoFHandler<dim1, spacedim> &                    dh1,
+    const Quadrature<dim0> &                              quadrature0,
+    const Quadrature<dim1> &                              quadrature1,
+    Matrix &                                              matrix,
+    const AffineConstraints<typename Matrix::value_type> &constraints0 =
+      AffineConstraints<typename Matrix::value_type>(),
+    const ComponentMask &comps0 = ComponentMask(),
+    const ComponentMask &comps1 = ComponentMask());
 } // namespace NonMatching
 DEAL_II_NAMESPACE_CLOSE
 
index 4df28b9fa2199050543af16b026fe87680100745..16c34f830f534bc67c5fe1bcea000e6632d9cf31 100644 (file)
@@ -124,6 +124,46 @@ namespace NonMatching
       }
       return {std::move(points_over_local_cells), std::move(used_cells_ids)};
     }
+
+
+    /**
+     * Given two ComponentMasks and the corresponding finite element spaces,
+     * compute a pairing between the selected components of the first finite
+     * element space, and the selected components of the second finite element
+     * space.
+     */
+    template <int dim0, int dim1, int spacedim>
+    std::pair<std::vector<unsigned int>, std::vector<unsigned int>>
+    compute_components_coupling(const ComponentMask &                comps0,
+                                const ComponentMask &                comps1,
+                                const FiniteElement<dim0, spacedim> &fe0,
+                                const FiniteElement<dim1, spacedim> &fe1)
+    {
+      // Take care of components
+      const ComponentMask mask0 =
+        (comps0.size() == 0 ? ComponentMask(fe0.n_components(), true) : comps0);
+
+      const ComponentMask mask1 =
+        (comps1.size() == 0 ? ComponentMask(fe1.n_components(), true) : comps1);
+
+      AssertDimension(mask0.size(), fe0.n_components());
+      AssertDimension(mask1.size(), fe1.n_components());
+
+      // Global to local indices
+      std::vector<unsigned int> gtl0(fe0.n_components(),
+                                     numbers::invalid_unsigned_int);
+      std::vector<unsigned int> gtl1(fe1.n_components(),
+                                     numbers::invalid_unsigned_int);
+
+      for (unsigned int i = 0, j = 0; i < gtl0.size(); ++i)
+        if (mask0[i])
+          gtl0[i] = j++;
+
+      for (unsigned int i = 0, j = 0; i < gtl1.size(); ++i)
+        if (mask1[i])
+          gtl1[i] = j++;
+      return {gtl0, gtl1};
+    }
   } // namespace internal
 
   template <int dim0,
@@ -177,7 +217,8 @@ namespace NonMatching
   {
     AssertDimension(sparsity.n_rows(), space_dh.n_dofs());
     AssertDimension(sparsity.n_cols(), immersed_dh.n_dofs());
-    static_assert(dim1 <= dim0, "This function can only work if dim1 <= dim0");
+    Assert(dim1 <= dim0,
+           ExcMessage("This function can only work if dim1 <= dim0"));
     Assert((dynamic_cast<
               const parallel::distributed::Triangulation<dim1, spacedim> *>(
               &immersed_dh.get_triangulation()) == nullptr),
@@ -360,7 +401,8 @@ namespace NonMatching
   {
     AssertDimension(matrix.m(), space_dh.n_dofs());
     AssertDimension(matrix.n(), immersed_dh.n_dofs());
-    static_assert(dim1 <= dim0, "This function can only work if dim1 <= dim0");
+    Assert(dim1 <= dim0,
+           ExcMessage("This function can only work if dim1 <= dim0"));
     Assert((dynamic_cast<
               const parallel::distributed::Triangulation<dim1, spacedim> *>(
               &immersed_dh.get_triangulation()) == nullptr),
@@ -572,8 +614,411 @@ namespace NonMatching
       }
   }
 
+  template <int dim0,
+            int dim1,
+            int spacedim,
+            typename Sparsity,
+            typename Number>
+  void
+  create_coupling_sparsity_pattern(
+    const double &                          epsilon,
+    const GridTools::Cache<dim0, spacedim> &cache0,
+    const GridTools::Cache<dim1, spacedim> &cache1,
+    const DoFHandler<dim0, spacedim> &      dh0,
+    const DoFHandler<dim1, spacedim> &      dh1,
+    Sparsity &                              sparsity,
+    const AffineConstraints<Number> &       constraints0,
+    const ComponentMask &                   comps0,
+    const ComponentMask &                   comps1)
+  {
+    if (epsilon == 0.0)
+      {
+        Assert(dim1 <= dim0,
+               ExcMessage("When epsilon is zero, you can only "
+                          "call this function with dim1 <= dim0."));
+        create_coupling_sparsity_pattern(cache0,
+                                         dh0,
+                                         dh1,
+                                         QGauss<dim1>(dh1.get_fe().degree + 1),
+                                         sparsity,
+                                         constraints0,
+                                         comps0,
+                                         comps1,
+                                         cache1.get_mapping());
+        return;
+      }
+    AssertDimension(sparsity.n_rows(), dh0.n_dofs());
+    AssertDimension(sparsity.n_cols(), dh1.n_dofs());
+
+    const bool zero_is_distributed =
+      (dynamic_cast<const parallel::distributed::Triangulation<dim1, spacedim>
+                      *>(&dh0.get_triangulation()) != nullptr);
+    const bool one_is_distributed =
+      (dynamic_cast<const parallel::distributed::Triangulation<dim1, spacedim>
+                      *>(&dh1.get_triangulation()) != nullptr);
+
+    // We bail out if both are distributed triangulations
+    Assert(!zero_is_distributed || !one_is_distributed, ExcNotImplemented());
+
+    // If we can loop on both, we decide where to make the outer loop according
+    // to the size of triangulation. The reasoning is the following:
+    // - Access to the tree: log(N)
+    // - We compute intersection for each of the outer loop cells (M)
+    // Total cost (besides the setup) is: M log(N)
+    // If we can, make sure M is the smallest number
+    const bool outer_loop_on_zero =
+      (zero_is_distributed && !one_is_distributed) ||
+      (dh1.get_triangulation().n_active_cells() >
+       dh0.get_triangulation().n_active_cells());
+
+    const auto &fe0 = dh0.get_fe();
+    const auto &fe1 = dh1.get_fe();
+
+    // Dof indices
+    std::vector<types::global_dof_index> dofs0(fe0.dofs_per_cell);
+    std::vector<types::global_dof_index> dofs1(fe1.dofs_per_cell);
+
+    if (outer_loop_on_zero)
+      {
+        std::cout << "Looping on zero." << std::endl;
+
+        Assert(one_is_distributed == false, ExcInternalError());
+
+        const auto &tree1 = cache1.get_cell_bounding_boxes_rtree();
+
+        std::vector<std::pair<
+          BoundingBox<spacedim>,
+          typename Triangulation<dim1, spacedim>::active_cell_iterator>>
+          intersection;
+
+        for (const auto &cell0 : dh0.active_cell_iterators())
+          if (cell0->is_locally_owned())
+            {
+              intersection.resize(0);
+
+              BoundingBox<spacedim> box0;
+              if (cache0.get_mapping().preserves_vertex_locations())
+                box0 = cell0->bounding_box();
+              else
+                {
+                  const auto vertices =
+                    cache0.get_mapping().get_vertices(cell0);
+                  Point<spacedim> p0 = vertices[0], p1 = vertices[0];
+                  for (unsigned int j = 1; j < vertices.size(); ++j)
+                    for (unsigned int d = 0; d < spacedim; ++d)
+                      {
+                        p0[d] = std::min(p0[d], vertices[j][d]);
+                        p1[d] = std::max(p1[d], vertices[j][d]);
+                      }
+                  box0 = BoundingBox<spacedim>({p0, p1});
+                }
+
+              box0.extend(epsilon);
+              boost::geometry::index::query(tree1,
+                                            boost::geometry::index::intersects(
+                                              box0),
+                                            std::back_inserter(intersection));
+              if (!intersection.empty())
+                {
+                  cell0->get_dof_indices(dofs0);
+                  for (const auto &entry : intersection)
+                    {
+                      typename DoFHandler<dim1, spacedim>::cell_iterator cell1(
+                        *entry.second, &dh1);
+                      cell1->get_dof_indices(dofs1);
+                      constraints0.add_entries_local_to_global(dofs0,
+                                                               dofs1,
+                                                               sparsity);
+                    }
+                }
+            }
+      }
+    else
+      {
+        Assert(zero_is_distributed == false, ExcInternalError());
+        const auto &tree0 = cache0.get_cell_bounding_boxes_rtree();
+
+        std::vector<std::pair<
+          BoundingBox<spacedim>,
+          typename Triangulation<dim0, spacedim>::active_cell_iterator>>
+          intersection;
+
+        for (const auto &cell1 : dh1.active_cell_iterators())
+          if (cell1->is_locally_owned())
+            {
+              intersection.resize(0);
+
+              BoundingBox<spacedim> box1;
+              Assert(&cache1.get_mapping() != nullptr, ExcInternalError());
+              if (cache1.get_mapping().preserves_vertex_locations())
+                box1 = cell1->bounding_box();
+              else
+                {
+                  const auto vertices =
+                    cache1.get_mapping().get_vertices(cell1);
+                  Point<spacedim> p0 = vertices[0], p1 = vertices[0];
+                  for (unsigned int j = 1; j < vertices.size(); ++j)
+                    for (unsigned int d = 0; d < spacedim; ++d)
+                      {
+                        p0[d] = std::min(p0[d], vertices[j][d]);
+                        p1[d] = std::max(p1[d], vertices[j][d]);
+                      }
+                  box1 = BoundingBox<spacedim>({p0, p1});
+                }
+
+              box1.extend(epsilon);
+              boost::geometry::index::query(tree0,
+                                            boost::geometry::index::intersects(
+                                              box1),
+                                            std::back_inserter(intersection));
+              if (!intersection.empty())
+                {
+                  cell1->get_dof_indices(dofs1);
+                  for (const auto &entry : intersection)
+                    {
+                      typename DoFHandler<dim0, spacedim>::cell_iterator cell0(
+                        *entry.second, &dh0);
+                      cell0->get_dof_indices(dofs0);
+                      constraints0.add_entries_local_to_global(dofs0,
+                                                               dofs1,
+                                                               sparsity);
+                    }
+                }
+            }
+      }
+  }
+
+
+
+  template <int dim0, int dim1, int spacedim, typename Matrix>
+  void
+  create_coupling_mass_matrix(
+    Functions::CutOffFunctionBase<spacedim> &             kernel,
+    const double &                                        epsilon,
+    const GridTools::Cache<dim0, spacedim> &              cache0,
+    const GridTools::Cache<dim1, spacedim> &              cache1,
+    const DoFHandler<dim0, spacedim> &                    dh0,
+    const DoFHandler<dim1, spacedim> &                    dh1,
+    const Quadrature<dim0> &                              quadrature0,
+    const Quadrature<dim1> &                              quadrature1,
+    Matrix &                                              matrix,
+    const AffineConstraints<typename Matrix::value_type> &constraints0,
+    const ComponentMask &                                 comps0,
+    const ComponentMask &                                 comps1)
+  {
+    if (epsilon == 0)
+      {
+        Assert(dim1 <= dim0,
+               ExcMessage("When epsilon is zero, you can only "
+                          "call this function with dim1 <= dim0."));
+        create_coupling_mass_matrix(cache0,
+                                    dh0,
+                                    dh1,
+                                    quadrature1,
+                                    matrix,
+                                    constraints0,
+                                    comps0,
+                                    comps1,
+                                    cache1.get_mapping());
+        return;
+      }
+
+    AssertDimension(matrix.m(), dh0.n_dofs());
+    AssertDimension(matrix.n(), dh1.n_dofs());
+
+    const bool zero_is_distributed =
+      (dynamic_cast<const parallel::distributed::Triangulation<dim1, spacedim>
+                      *>(&dh0.get_triangulation()) != nullptr);
+    const bool one_is_distributed =
+      (dynamic_cast<const parallel::distributed::Triangulation<dim1, spacedim>
+                      *>(&dh1.get_triangulation()) != nullptr);
+
+    // We bail out if both are distributed triangulations
+    Assert(!zero_is_distributed || !one_is_distributed, ExcNotImplemented());
+
+    // If we can loop on both, we decide where to make the outer loop according
+    // to the size of triangulation. The reasoning is the following:
+    // - Access to the tree: log(N)
+    // - We compute intersection for each of the outer loop cells (M)
+    // Total cost (besides the setup) is: M log(N)
+    // If we can, make sure M is the smallest number
+    const bool outer_loop_on_zero =
+      (zero_is_distributed && !one_is_distributed) ||
+      (dh1.get_triangulation().n_active_cells() >
+       dh0.get_triangulation().n_active_cells());
+
+    const auto &fe0 = dh0.get_fe();
+    const auto &fe1 = dh1.get_fe();
+
+    FEValues<dim0, spacedim> fev0(cache0.get_mapping(),
+                                  fe0,
+                                  quadrature0,
+                                  update_values | update_JxW_values |
+                                    update_quadrature_points);
+
+    FEValues<dim1, spacedim> fev1(cache1.get_mapping(),
+                                  fe1,
+                                  quadrature1,
+                                  update_values | update_JxW_values |
+                                    update_quadrature_points);
+
+    // Dof indices
+    std::vector<types::global_dof_index> dofs0(fe0.dofs_per_cell);
+    std::vector<types::global_dof_index> dofs1(fe1.dofs_per_cell);
+
+    // Local Matrix
+    FullMatrix<typename Matrix::value_type> cell_matrix(fe0.dofs_per_cell,
+                                                        fe1.dofs_per_cell);
+
+    // Global to local indices
+    auto p = internal::compute_components_coupling(comps0, comps1, fe0, fe1);
+    const auto &gtl0 = p.first;
+    const auto &gtl1 = p.second;
+
+    kernel.set_radius(epsilon);
+    std::vector<double> kernel_values(quadrature1.size());
+
+    auto assemble_one_pair = [&]() {
+      cell_matrix = 0;
+      for (unsigned int q0 = 0; q0 < quadrature0.size(); ++q0)
+        {
+          kernel.set_center(fev0.quadrature_point(q0));
+          kernel.value_list(fev1.get_quadrature_points(), kernel_values);
+          for (unsigned int q1 = 0; q1 < quadrature1.size(); ++q1)
+            if (kernel_values[q1] != 0.0)
+              {
+                for (unsigned int i = 0; i < fe0.dofs_per_cell; ++i)
+                  {
+                    const auto comp_i = fe0.system_to_component_index(i).first;
+                    if (gtl0[comp_i] != numbers::invalid_unsigned_int)
+                      for (unsigned int j = 0; j < fe1.dofs_per_cell; ++j)
+                        {
+                          const auto comp_j =
+                            fe1.system_to_component_index(j).first;
+                          if (gtl1[comp_j] == gtl0[comp_i])
+                            {
+                              cell_matrix(i, j) += fev0.shape_value(i, q0) *
+                                                   fev1.shape_value(j, q1) *
+                                                   kernel_values[q1] *
+                                                   fev0.JxW(q0) * fev1.JxW(q1);
+                            }
+                        }
+                  }
+              }
+        }
+
+      constraints0.distribute_local_to_global(cell_matrix,
+                                              dofs0,
+                                              dofs1,
+                                              matrix);
+    };
+
+    if (outer_loop_on_zero)
+      {
+        Assert(one_is_distributed == false, ExcInternalError());
+
+        const auto &tree1 = cache1.get_cell_bounding_boxes_rtree();
+
+        std::vector<std::pair<
+          BoundingBox<spacedim>,
+          typename Triangulation<dim1, spacedim>::active_cell_iterator>>
+          intersection;
+
+        for (const auto &cell0 : dh0.active_cell_iterators())
+          if (cell0->is_locally_owned())
+            {
+              intersection.resize(0);
+
+              BoundingBox<spacedim> box0;
+              if (cache0.get_mapping().preserves_vertex_locations())
+                box0 = cell0->bounding_box();
+              else
+                {
+                  const auto vertices =
+                    cache0.get_mapping().get_vertices(cell0);
+                  Point<spacedim> p0 = vertices[0], p1 = vertices[0];
+                  for (unsigned int j = 1; j < vertices.size(); ++j)
+                    for (unsigned int d = 0; d < spacedim; ++d)
+                      {
+                        p0[d] = std::min(p0[d], vertices[j][d]);
+                        p1[d] = std::max(p1[d], vertices[j][d]);
+                      }
+                  box0 = BoundingBox<spacedim>({p0, p1});
+                }
+
+              box0.extend(epsilon);
+              boost::geometry::index::query(tree1,
+                                            boost::geometry::index::intersects(
+                                              box0),
+                                            std::back_inserter(intersection));
+              if (!intersection.empty())
+                {
+                  cell0->get_dof_indices(dofs0);
+                  fev0.reinit(cell0);
+                  for (const auto &entry : intersection)
+                    {
+                      typename DoFHandler<dim1, spacedim>::cell_iterator cell1(
+                        *entry.second, &dh1);
+                      cell1->get_dof_indices(dofs1);
+                      fev1.reinit(cell1);
+                      assemble_one_pair();
+                    }
+                }
+            }
+      }
+    else
+      {
+        Assert(zero_is_distributed == false, ExcInternalError());
+        const auto &tree0 = cache0.get_cell_bounding_boxes_rtree();
+
+        std::vector<std::pair<
+          BoundingBox<spacedim>,
+          typename Triangulation<dim0, spacedim>::active_cell_iterator>>
+          intersection;
+
+        for (const auto &cell1 : dh1.active_cell_iterators())
+          if (cell1->is_locally_owned())
+            {
+              intersection.resize(0);
+              BoundingBox<spacedim> box1;
+              if (cache1.get_mapping().preserves_vertex_locations())
+                box1 = cell1->bounding_box();
+              else
+                {
+                  const auto vertices =
+                    cache1.get_mapping().get_vertices(cell1);
+                  Point<spacedim> p0 = vertices[0], p1 = vertices[0];
+                  for (unsigned int j = 1; j < vertices.size(); ++j)
+                    for (unsigned int d = 0; d < spacedim; ++d)
+                      {
+                        p0[d] = std::min(p0[d], vertices[j][d]);
+                        p1[d] = std::max(p1[d], vertices[j][d]);
+                      }
+                  box1 = BoundingBox<spacedim>({p0, p1});
+                }
+              box1.extend(epsilon);
+              boost::geometry::index::query(tree0,
+                                            boost::geometry::index::intersects(
+                                              box1),
+                                            std::back_inserter(intersection));
+              if (!intersection.empty())
+                {
+                  cell1->get_dof_indices(dofs1);
+                  fev1.reinit(cell1);
+                  for (const auto &entry : intersection)
+                    {
+                      typename DoFHandler<dim0, spacedim>::cell_iterator cell0(
+                        *entry.second, &dh0);
+                      cell0->get_dof_indices(dofs0);
+                      fev0.reinit(cell0);
+                      assemble_one_pair();
+                    }
+                }
+            }
+      }
+  }
+
 #include "coupling.inst"
 } // namespace NonMatching
 
-
 DEAL_II_NAMESPACE_CLOSE
index c1dcd9e5746d6d7557ae1d06aad90bc14270dcc8..cfc2022cb5eb8addf4fb0c02419083b9f089d27e 100644 (file)
@@ -18,7 +18,7 @@ for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS;
      Sparsity : SPARSITY_PATTERNS;
      S : REAL_AND_COMPLEX_SCALARS)
   {
-#if dim1 <= dim0 && dim0 <= spacedim
+#if dim1 <= spacedim && dim0 <= spacedim
     template void create_coupling_sparsity_pattern(
       const DoFHandler<dim0, spacedim> &space_dh,
       const DoFHandler<dim1, spacedim> &immersed_dh,
@@ -40,6 +40,17 @@ for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS;
       const ComponentMask &                   space_comps,
       const ComponentMask &                   immersed_comps,
       const Mapping<dim1, spacedim> &         immersed_mapping);
+
+    template void create_coupling_sparsity_pattern(
+      const double &                          epsilon,
+      const GridTools::Cache<dim0, spacedim> &cache0,
+      const GridTools::Cache<dim1, spacedim> &cache1,
+      const DoFHandler<dim0, spacedim> &      dh0,
+      const DoFHandler<dim1, spacedim> &      dh1,
+      Sparsity &                              sparsity,
+      const AffineConstraints<S> &            constraints0,
+      const ComponentMask &                   comps0,
+      const ComponentMask &                   comps1);
 #endif
   }
 
@@ -47,7 +58,7 @@ for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS;
 for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS;
      Matrix : SPARSE_MATRICES)
   {
-#if dim1 <= dim0 && dim0 <= spacedim
+#if dim1 <= spacedim && dim0 <= spacedim
     template void create_coupling_mass_matrix(
       const DoFHandler<dim0, spacedim> &           space_dh,
       const DoFHandler<dim1, spacedim> &           immersed_dh,
@@ -69,5 +80,19 @@ for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS;
       const ComponentMask &                        space_comps,
       const ComponentMask &                        immersed_comps,
       const Mapping<dim1, spacedim> &              immersed_mapping);
+
+    template void create_coupling_mass_matrix(
+      Functions::CutOffFunctionBase<spacedim> & kernel,
+      const double &                                        epsilon,
+      const GridTools::Cache<dim0, spacedim> &              cache0,
+      const GridTools::Cache<dim1, spacedim> &              cache1,
+      const DoFHandler<dim0, spacedim> &                    dh0,
+      const DoFHandler<dim1, spacedim> &                    dh1,
+      const Quadrature<dim0> &                              quadrature0,
+      const Quadrature<dim1> &                              quadrature1,
+      Matrix &                                              matrix,
+      const AffineConstraints<typename Matrix::value_type> &constraints0,
+      const ComponentMask &                                 comps0,
+      const ComponentMask &                                 comps1);
 #endif
-  }
+  }
\ No newline at end of file
diff --git a/tests/non_matching/coupling_08.cc b/tests/non_matching/coupling_08.cc
new file mode 100644 (file)
index 0000000..eedf9b1
--- /dev/null
@@ -0,0 +1,162 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/function_parser.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
+
+#include <deal.II/non_matching/coupling.h>
+
+#include <deal.II/numerics/matrix_tools.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+// Test that a coupling matrix can be constructed for each pair of dimension and
+// immersed dimension, and check that constants are projected correctly.
+//
+// Even when locally refined grids are used.
+
+template <int dim0, int dim1, int spacedim>
+void
+test()
+{
+  deallog << "dim0: " << dim0 << ", dim1: " << dim1
+          << ", spacedim: " << spacedim << std::endl;
+
+  Triangulation<dim0, spacedim> tria0;
+  Triangulation<dim1, spacedim> tria1;
+
+  GridGenerator::hyper_cube(tria0, -.4, .3);
+  GridGenerator::hyper_cube(tria1, -1, 1);
+
+  tria0.refine_global(3);
+  tria1.refine_global(2);
+
+  FE_Q<dim0, spacedim> fe0(1);
+  FE_Q<dim1, spacedim> fe1(1);
+
+  deallog << "FE0             : " << fe0.get_name() << std::endl
+          << "FE1             : " << fe1.get_name() << std::endl;
+
+  DoFHandler<dim0, spacedim> dh0(tria0);
+  DoFHandler<dim1, spacedim> dh1(tria1);
+
+  GridTools::Cache<dim0, spacedim> cache0(tria0);
+  GridTools::Cache<dim1, spacedim> cache1(tria1);
+
+  dh0.distribute_dofs(fe0);
+  dh1.distribute_dofs(fe1);
+
+  AffineConstraints<double> constraints0;
+  DoFTools::make_hanging_node_constraints(dh0, constraints0);
+
+  constraints0.close();
+
+  deallog << "Dofs 0          : " << dh0.n_dofs() << std::endl
+          << "Dofs 1          : " << dh1.n_dofs() << std::endl
+          << "Constrained dofs: " << constraints0.n_constraints() << std::endl;
+
+  QGauss<dim0> quad0(3); // Quadrature for coupling
+  QGauss<dim1> quad1(3); // Quadrature for coupling
+
+  const double epsilon = 2 * std::max(GridTools::maximal_cell_diameter(tria0),
+                                      GridTools::maximal_cell_diameter(tria1));
+
+  deallog << "Epsilon: " << epsilon << std::endl;
+
+  SparsityPattern sparsity;
+  {
+    DynamicSparsityPattern dsp(dh0.n_dofs(), dh1.n_dofs());
+    NonMatching::create_coupling_sparsity_pattern(
+      epsilon, cache0, cache1, dh0, dh1, dsp, constraints0);
+    sparsity.copy_from(dsp);
+  }
+  SparseMatrix<double> coupling(sparsity);
+
+  Functions::CutOffFunctionC1<spacedim> dirac(
+    1,
+    Point<spacedim>(),
+    1,
+    Functions::CutOffFunctionBase<spacedim>::no_component,
+    true);
+
+  NonMatching::create_coupling_mass_matrix(dirac,
+                                           epsilon,
+                                           cache0,
+                                           cache1,
+                                           dh0,
+                                           dh1,
+                                           quad0,
+                                           quad1,
+                                           coupling,
+                                           constraints0);
+
+  SparsityPattern mass_sparsity1;
+  {
+    DynamicSparsityPattern dsp(dh1.n_dofs(), dh1.n_dofs());
+    DoFTools::make_sparsity_pattern(dh1, dsp);
+    mass_sparsity1.copy_from(dsp);
+  }
+  SparseMatrix<double> mass_matrix1(mass_sparsity1);
+  MatrixTools::create_mass_matrix(dh1, quad1, mass_matrix1);
+
+  SparseDirectUMFPACK mass_matrix1_inv;
+  mass_matrix1_inv.factorize(mass_matrix1);
+
+  // now take ones in dh0, project them onto dh1,
+  // get back ones, and check for the error.
+  //
+  // WARNINGS: Only works if dh1 is immersed in dh0
+
+  Vector<double> ones0(dh0.n_dofs());
+  Vector<double> ones1(dh1.n_dofs());
+
+  ones0 = 1.0;
+  coupling.Tvmult(ones1, ones0);
+  mass_matrix1_inv.solve(ones1);
+
+  Vector<double> exact_ones1(dh1.n_dofs());
+  exact_ones1 = 1.0;
+  ones1 -= exact_ones1;
+
+  deallog << "Error on constants: " << ones1.l2_norm() << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog(1);
+  test<1, 1, 1>();
+  test<2, 1, 2>();
+  test<2, 2, 2>();
+  test<3, 2, 3>();
+  test<3, 3, 3>();
+}
diff --git a/tests/non_matching/coupling_08.with_umfpack=true.output b/tests/non_matching/coupling_08.with_umfpack=true.output
new file mode 100644 (file)
index 0000000..e69de29

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.