]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Non matching sparsity and mass matrix.
authorLuca Heltai <luca.heltai@sissa.it>
Mon, 5 Mar 2018 14:34:20 +0000 (15:34 +0100)
committerLuca Heltai <luca.heltai@sissa.it>
Mon, 5 Mar 2018 14:38:47 +0000 (15:38 +0100)
include/deal.II/non_matching/coupling.h [new file with mode: 0644]
tests/non_matching/coupling_01.cc [new file with mode: 0644]
tests/non_matching/coupling_01.output [new file with mode: 0644]
tests/non_matching/coupling_02.cc [new file with mode: 0644]
tests/non_matching/coupling_02.output [new file with mode: 0644]
tests/non_matching/coupling_03.cc [new file with mode: 0644]
tests/non_matching/coupling_03.output [new file with mode: 0644]
tests/non_matching/coupling_04.cc [new file with mode: 0644]
tests/non_matching/coupling_04.output [new file with mode: 0644]

diff --git a/include/deal.II/non_matching/coupling.h b/include/deal.II/non_matching/coupling.h
new file mode 100644 (file)
index 0000000..3b8be3e
--- /dev/null
@@ -0,0 +1,348 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_non_matching_coupling
+#define dealii_non_matching_coupling
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/quadrature.h>
+
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_cache.h>
+
+#include <deal.II/fe/component_mask.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+
+
+
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+namespace NonMatching
+{
+  /**
+   * Create a coupling sparsity pattern for non-matching, overlapping grids.
+   *
+   * Given two non-matching triangulations, representing the domains $\Omega$
+   * and $B$, with $B \subseteq \Omega$, and two finite element spaces
+   * $V(\Omega)$ and $Q(B)$, compute the sparsity pattern that would be
+   * necessary to assemble the matrix
+   *
+   * \f[
+   * M_{ij} := \int_{B} v_i(x) w_j(x) dx, \quad \foreach v_i \in V(\Omega), w_j \in Q(B)
+   * \f]
+   *
+   * where $V(\Omega)$ is the finite element space associated with `space_dh`
+   * (or part of it, if specified in `space_comps`), while $Q(B)$ is the finite
+   * element space associated with `immersed_dh` (or part of it, if specified
+   * in a `immersed_comps`).
+   *
+   * The `sparsity` is constructed by locating the position of quadrature
+   * points (obtained by the reference quadrature `quad`) defined on elements
+   * of $B$ with respect to the embedding triangulation $\Omega$. For each
+   * overlapping cell, the entries corresponding to `space_comps` in `space_dh`
+   * and `immersed_comps` in `immersed_dh` are added to the sparsity pattern.
+   *
+   * The `space_comps` and `immersed_comps` masks are assumed to be ordered in
+   * the same way: the first component of `space_comps` will couple with the
+   * first component of `immersed_comps`, the second with the second, and so
+   * on. If one of the two masks has more non-zero entries w.r.t. the other,
+   * then the excess components will be ignored.
+   *
+   * If the domain $B$ does not fall withing $\Omega$, an exception will be
+   * thrown by the algorithm that computes the quadrature point locations.
+   *
+   * For both spaces, it is possible to specify a custom Mapping, which
+   * defaults to StaticMappingQ1 for both.
+   *
+   * @author Luca Heltai, 2018
+   */
+  template<int dim0, int dim1, int spacedim, typename Sparsity>
+  void create_coupling_sparsity_pattern(const DoFHandler<dim0, spacedim> &space_dh,
+                                        const DoFHandler<dim1, spacedim> &immersed_dh,
+                                        const Quadrature<dim1>           &quad,
+                                        Sparsity                         &sparsity,
+                                        const ComponentMask              &space_comps = ComponentMask(),
+                                        const ComponentMask              &immersed_comps = ComponentMask(),
+                                        const Mapping<dim0, spacedim>    &space_mapping = StaticMappingQ1<dim0,spacedim>::mapping,
+                                        const Mapping<dim1, spacedim>    &immersed_mapping = StaticMappingQ1<dim1, spacedim>::mapping);
+
+  /**
+   * Create a coupling mass matrix for non-matching, overlapping grids.
+   *
+   * Given two non-matching triangulations, representing the domains $\Omega$
+   * and $B$, with $B \subseteq \Omega$, and two finite element spaces
+   * $V(\Omega)$ and $Q(B)$, compute the matrix
+   *
+   * \f[
+   * M_{ij} := \int_{B} v_i(x) w_j(x) dx, \quad \foreach v_i \in V(\Omega), w_j \in Q(B)
+   * \f]
+   *
+   * where $V(\Omega)$ is the finite element space associated with `space_dh`
+   * (or part of it, if specified in `space_comps`), while $Q(B)$ is the finite
+   * element space associated with `immersed_dh` (or part of it, if specified
+   * in a `immersed_comps`).
+   *
+   * The corresponding `sparsity` can be computed by calling the
+   * make_coupling_sparsity_pattern function. The elements of the matrix are
+   * computed by locating the position of quadrature points defined on elements
+   * of $B$ with respect to the embedding triangulation $\Omega$.
+   *
+   * The `space_comps` and `immersed_comps` masks are assumed to be ordered in
+   * the same way: the first component of `space_comps` will couple with the
+   * first component of `immersed_comps`, the second with the second, and so
+   * on. If one of the two masks has more non-zero entries w.r.t. the other,
+   * then the excess components will be ignored.
+   *
+   * If the domain $B$ does not fall withing $\Omega$, an exception will be
+   * thrown by the algorithm that computes the quadrature point locations.
+   *
+   * For both spaces, it is possible to specify a custom Mapping, which
+   * defaults to StaticMappingQ1 for both.
+   *
+   * @author Luca Heltai, 2018
+   */
+  template<int dim0, int dim1, int spacedim, typename Matrix>
+  void create_coupling_mass_matrix(const DoFHandler<dim0, spacedim> &space_dh,
+                                   const DoFHandler<dim1, spacedim> &immersed_dh,
+                                   const Quadrature<dim1>           &quad,
+                                   Matrix                           &matrix,
+                                   const ConstraintMatrix           &constraints = ConstraintMatrix(),
+                                   const ComponentMask              &space_comps = ComponentMask(),
+                                   const ComponentMask              &immersed_comps = ComponentMask(),
+                                   const Mapping<dim0, spacedim>    &space_mapping = StaticMappingQ1<dim0,spacedim>::mapping,
+                                   const Mapping<dim1, spacedim>    &immersed_mapping = StaticMappingQ1<dim1, spacedim>::mapping);
+
+
+
+
+// === inline and template functions ===
+
+
+
+  template<int dim0, int dim1, int spacedim, typename Sparsity>
+  void create_coupling_sparsity_pattern(const DoFHandler<dim0, spacedim> &space_dh,
+                                        const DoFHandler<dim1, spacedim> &immersed_dh,
+                                        const Quadrature<dim1>           &quad,
+                                        Sparsity                         &sparsity,
+                                        const ComponentMask              &space_comps,
+                                        const ComponentMask              &immersed_comps,
+                                        const Mapping<dim0, spacedim>    &space_mapping,
+                                        const Mapping<dim1, spacedim>    &immersed_mapping)
+  {
+    AssertDimension(sparsity.n_rows(), space_dh.n_dofs());
+    AssertDimension(sparsity.n_cols(), immersed_dh.n_dofs());
+
+    const auto &space_fe = space_dh.get_fe();
+    const auto &immersed_fe = immersed_dh.get_fe();
+
+    // Now we run on ech cell, get a quadrature formula
+    typename DoFHandler<dim1,spacedim>::active_cell_iterator
+    cell = immersed_dh.begin_active(),
+    endc = immersed_dh.end();
+
+    // Dof indices
+    std::vector<types::global_dof_index> dofs(immersed_fe.dofs_per_cell);
+    std::vector<types::global_dof_index> odofs(space_fe.dofs_per_cell);
+
+    FEValues<dim1,spacedim> fe_v(immersed_mapping, immersed_fe, quad,
+                                 update_quadrature_points);
+
+    GridTools::Cache<dim0,spacedim> cache(space_dh.get_triangulation(), space_mapping);
+
+    // Take care of components
+    ComponentMask space_c = space_comps;
+    ComponentMask immersed_c = immersed_comps;
+    if (space_c.size() == 0)
+      space_c = ComponentMask(space_fe.n_components(), true);
+    if (immersed_c.size() == 0)
+      immersed_c = ComponentMask(immersed_fe.n_components(), true);
+
+    AssertDimension(space_c.size(), space_fe.n_components());
+    AssertDimension(immersed_c.size(), immersed_fe.n_components());
+
+
+    std::vector<unsigned int> space_gtl(space_fe.n_components(), numbers::invalid_unsigned_int);
+    std::vector<unsigned int> immersed_gtl(immersed_fe.n_components(), numbers::invalid_unsigned_int);
+
+    for (unsigned int i=0, j=0; i<space_gtl.size(); ++i)
+      if (space_c[i])
+        space_gtl[i] = j++;
+
+    for (unsigned int i=0, j=0; i<immersed_gtl.size(); ++i)
+      if (immersed_c[i])
+        immersed_gtl[i] = j++;
+
+    for (; cell != endc; ++cell)
+      //    if(cell->subdomain_id() == this_mpi_process)
+      {
+        // Reinitialize the cell and the fe_values
+        fe_v.reinit(cell);
+        cell->get_dof_indices(dofs);
+
+        const std::vector<Point<spacedim> > &Xpoints = fe_v.get_quadrature_points();
+
+        // Get a list of outer cells, qpoints and maps.
+        const auto cpm = GridTools::compute_point_locations(cache, Xpoints);
+        const auto &cells = std::get<0>(cpm);
+
+        for (unsigned int c=0; c<cells.size(); ++c)
+          {
+            // Get the ones in the current outer cell
+            typename DoFHandler<dim0,spacedim>::cell_iterator
+            ocell(*cells[c], &space_dh);
+            ocell->get_dof_indices(odofs);
+            for (unsigned int i=0; i<odofs.size(); ++i)
+              {
+                const auto comp_i = space_fe.system_to_component_index(i).first;
+                if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
+                  {
+                    for (unsigned int j=0; j<dofs.size(); ++j)
+                      {
+                        const auto comp_j = immersed_fe.system_to_component_index(j).first;
+                        if (immersed_gtl[comp_j] == space_gtl[comp_i])
+                          sparsity.add(odofs[i],dofs[j]);
+                      }
+                  }
+              }
+          }
+      }
+  }
+
+
+
+  template<int dim0, int dim1, int spacedim, typename Matrix>
+  void create_coupling_mass_matrix(const DoFHandler<dim0, spacedim> &space_dh,
+                                   const DoFHandler<dim1, spacedim> &immersed_dh,
+                                   const Quadrature<dim1>           &quad,
+                                   Matrix                           &matrix,
+                                   const ConstraintMatrix           &constraints,
+                                   const ComponentMask              &space_comps,
+                                   const ComponentMask              &immersed_comps,
+                                   const Mapping<dim0, spacedim>    &space_mapping,
+                                   const Mapping<dim1, spacedim>    &immersed_mapping)
+  {
+    AssertDimension(matrix.m(), space_dh.n_dofs());
+    AssertDimension(matrix.n(), immersed_dh.n_dofs());
+
+    const auto &space_fe = space_dh.get_fe();
+    const auto &immersed_fe = immersed_dh.get_fe();
+
+    // Dof indices
+    std::vector<types::global_dof_index> dofs(immersed_fe.dofs_per_cell);
+    std::vector<types::global_dof_index> odofs(space_fe.dofs_per_cell);
+
+    GridTools::Cache<dim0,spacedim> cache(space_dh.get_triangulation(), space_mapping);
+
+    // Take care of components
+    ComponentMask space_c = space_comps;
+    ComponentMask immersed_c = immersed_comps;
+    if (space_c.size() == 0)
+      space_c = ComponentMask(space_fe.n_components(), true);
+    if (immersed_c.size() == 0)
+      immersed_c = ComponentMask(immersed_fe.n_components(), true);
+
+    AssertDimension(space_c.size(), space_fe.n_components());
+    AssertDimension(immersed_c.size(), immersed_fe.n_components());
+
+    std::vector<unsigned int> space_gtl(space_fe.n_components(), numbers::invalid_unsigned_int);
+    std::vector<unsigned int> immersed_gtl(immersed_fe.n_components(), numbers::invalid_unsigned_int);
+
+    for (unsigned int i=0, j=0; i<space_gtl.size(); ++i)
+      if (space_c[i])
+        space_gtl[i] = j++;
+
+    for (unsigned int i=0, j=0; i<immersed_gtl.size(); ++i)
+      if (immersed_c[i])
+        immersed_gtl[i] = j++;
+
+    FullMatrix<double> cell_matrix(space_dh.get_fe().dofs_per_cell,
+                                   immersed_dh.get_fe().dofs_per_cell);
+
+    FEValues<dim1,spacedim> fe_v(immersed_mapping, immersed_dh.get_fe(), quad,
+                                 update_JxW_values |
+                                 update_quadrature_points |
+                                 update_values);
+
+    // Now we run on ech cell, get a quadrature formula
+    typename DoFHandler<dim1,spacedim>::active_cell_iterator
+    cell = immersed_dh.begin_active(),
+    endc = immersed_dh.end();
+
+    for (; cell != endc; ++cell)
+      {
+        // Reinitialize the cell and the fe_values
+        fe_v.reinit(cell);
+        cell->get_dof_indices(dofs);
+
+        const std::vector<Point<spacedim> > &Xpoints = fe_v.get_quadrature_points();
+
+        // Get a list of outer cells, qpoints and maps.
+        const auto cpm      = GridTools::compute_point_locations(cache, Xpoints);
+        const auto &cells   = std::get<0>(cpm);
+        const auto &qpoints = std::get<1>(cpm);
+        const auto &maps    = std::get<2>(cpm);
+
+        for (unsigned int c=0; c<cells.size(); ++c)
+          {
+            // Get the ones in the current outer cell
+            typename DoFHandler<dim0,spacedim>::active_cell_iterator
+            ocell(*cells[c], &space_dh);
+            const std::vector< Point<spacedim> > &qps = qpoints[c];
+            const std::vector< unsigned int > &ids = maps[c];
+
+            FEValues<dim0,spacedim> o_fe_v(space_mapping, space_dh.get_fe(), qps,
+                                           update_values);
+            o_fe_v.reinit(ocell);
+            ocell->get_dof_indices(odofs);
+
+            // Reset the matrices.
+            cell_matrix = 0;
+
+            for (unsigned int i=0; i<space_dh.get_fe().dofs_per_cell; ++i)
+              {
+                const auto comp_i = space_dh.get_fe().system_to_component_index(i).first;
+                if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
+                  for (unsigned int j=0; j<immersed_dh.get_fe().dofs_per_cell; ++j)
+                    {
+                      const auto comp_j = immersed_dh.get_fe().system_to_component_index(j).first;
+                      if (space_gtl[comp_i] == immersed_gtl[comp_j])
+                        for (unsigned int oq=0; oq<o_fe_v.n_quadrature_points; ++oq)
+                          {
+                            // Get the corrisponding q point
+                            const unsigned int q=ids[oq];
+
+                            cell_matrix(i,j) += ( fe_v.shape_value(j,q) *
+                                                  o_fe_v.shape_value(i,oq) *
+                                                  fe_v.JxW(q) );
+                          }
+                    }
+              }
+
+            // Now assemble the matrices
+            constraints.distribute_local_to_global (cell_matrix, odofs, dofs, matrix);
+          }
+      }
+  }
+}
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/tests/non_matching/coupling_01.cc b/tests/non_matching/coupling_01.cc
new file mode 100644 (file)
index 0000000..655b17f
--- /dev/null
@@ -0,0 +1,118 @@
+// ---------------------------------------------------------------------
+//
+// 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/non_matching/coupling.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparse_direct.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.
+
+template<int dim, int spacedim>
+void test()
+{
+  deallog << "dim: " << dim << ", spacedim: "
+          << spacedim << std::endl;
+
+  Triangulation<dim,spacedim> tria;
+  Triangulation<spacedim,spacedim> space_tria;
+
+  GridGenerator::hyper_cube(tria,-.4,.3);
+  GridGenerator::hyper_cube(space_tria,-1,1);
+
+  tria.refine_global(1);
+  space_tria.refine_global(2);
+
+  FE_Q<dim,spacedim> fe(1);
+  FE_Q<spacedim,spacedim> space_fe(1);
+
+  deallog << "FE      : " << fe.get_name() << std::endl
+          << "Space FE: " << space_fe.get_name() << std::endl;
+
+  DoFHandler<dim,spacedim> dh(tria);
+  DoFHandler<spacedim,spacedim> space_dh(space_tria);
+
+  dh.distribute_dofs(fe);
+  space_dh.distribute_dofs(space_fe);
+
+  deallog << "Dofs      : " << dh.n_dofs() << std::endl
+          << "Space dofs: " << space_dh.n_dofs() << std::endl;
+
+  QGauss<dim> quad(3); // Quadrature for coupling
+
+
+  SparsityPattern sparsity;
+  {
+    DynamicSparsityPattern dsp(space_dh.n_dofs(), dh.n_dofs());
+    NonMatching::create_coupling_sparsity_pattern(space_dh, dh,
+                                                  quad, dsp);
+    sparsity.copy_from(dsp);
+  }
+  SparseMatrix<double> coupling(sparsity);
+  NonMatching::create_coupling_mass_matrix(space_dh, dh, quad, coupling);
+
+  SparsityPattern mass_sparsity;
+  {
+    DynamicSparsityPattern dsp(dh.n_dofs(), dh.n_dofs());
+    DoFTools::make_sparsity_pattern(dh,dsp);
+    mass_sparsity.copy_from(dsp);
+  }
+  SparseMatrix<double> mass_matrix(mass_sparsity);
+  MatrixTools::create_mass_matrix(dh, quad, mass_matrix);
+
+  SparseDirectUMFPACK mass_matrix_inv;
+  mass_matrix_inv.factorize(mass_matrix);
+
+  // now take ones in space, project them onto the immersed space,
+  // get back ones, and check for the error.
+  Vector<double> space_ones(space_dh.n_dofs());
+  Vector<double> ones(dh.n_dofs());
+
+  space_ones = 1.0;
+  coupling.Tvmult(ones, space_ones);
+  mass_matrix_inv.solve(ones);
+
+  Vector<double> real_ones(dh.n_dofs());
+  real_ones = 1.0;
+  ones -= real_ones;
+
+  deallog << "Error on constants: " << ones.l2_norm() << std::endl;
+}
+
+
+
+int main()
+{
+  initlog();
+  test<1,1>();
+  test<1,2>();
+  test<2,2>();
+  test<2,3>();
+  test<3,3>();
+}
diff --git a/tests/non_matching/coupling_01.output b/tests/non_matching/coupling_01.output
new file mode 100644 (file)
index 0000000..db8ef19
--- /dev/null
@@ -0,0 +1,31 @@
+
+DEAL::dim: 1, spacedim: 1
+DEAL::FE      : FE_Q<1>(1)
+DEAL::Space FE: FE_Q<1>(1)
+DEAL::Dofs      : 3
+DEAL::Space dofs: 5
+DEAL::Error on constants: 0.00000
+DEAL::dim: 1, spacedim: 2
+DEAL::FE      : FE_Q<1,2>(1)
+DEAL::Space FE: FE_Q<2>(1)
+DEAL::Dofs      : 3
+DEAL::Space dofs: 25
+DEAL::Error on constants: 0.00000
+DEAL::dim: 2, spacedim: 2
+DEAL::FE      : FE_Q<2>(1)
+DEAL::Space FE: FE_Q<2>(1)
+DEAL::Dofs      : 9
+DEAL::Space dofs: 25
+DEAL::Error on constants: 1.14304e-15
+DEAL::dim: 2, spacedim: 3
+DEAL::FE      : FE_Q<2,3>(1)
+DEAL::Space FE: FE_Q<3>(1)
+DEAL::Dofs      : 9
+DEAL::Space dofs: 125
+DEAL::Error on constants: 7.19507e-16
+DEAL::dim: 3, spacedim: 3
+DEAL::FE      : FE_Q<3>(1)
+DEAL::Space FE: FE_Q<3>(1)
+DEAL::Dofs      : 27
+DEAL::Space dofs: 125
+DEAL::Error on constants: 1.07909e-14
diff --git a/tests/non_matching/coupling_02.cc b/tests/non_matching/coupling_02.cc
new file mode 100644 (file)
index 0000000..9e95428
--- /dev/null
@@ -0,0 +1,123 @@
+// ---------------------------------------------------------------------
+//
+// 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/non_matching/coupling.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparse_direct.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, with a non trivial component selection in the space_dh,
+// and check that constants are projected correctly.
+
+template<int dim, int spacedim>
+void test()
+{
+  deallog << "dim: " << dim << ", spacedim: "
+          << spacedim << std::endl;
+
+  Triangulation<dim,spacedim> tria;
+  Triangulation<spacedim,spacedim> space_tria;
+
+  GridGenerator::hyper_cube(tria,-.4,.3);
+  GridGenerator::hyper_cube(space_tria,-1,1);
+
+  tria.refine_global(1);
+  space_tria.refine_global(2);
+
+  FE_Q<dim,spacedim> fe(1);
+  FESystem<spacedim,spacedim> space_fe(FE_Q<spacedim,spacedim>(1), spacedim+1);
+  ComponentMask space_mask(spacedim+1, false);
+  space_mask.set(spacedim, true);
+
+  deallog << "FE      : " << fe.get_name() << std::endl
+          << "Space FE: " << space_fe.get_name() << std::endl;
+
+  DoFHandler<dim,spacedim> dh(tria);
+  DoFHandler<spacedim,spacedim> space_dh(space_tria);
+
+  dh.distribute_dofs(fe);
+  space_dh.distribute_dofs(space_fe);
+
+  deallog << "Dofs      : " << dh.n_dofs() << std::endl
+          << "Space dofs: " << space_dh.n_dofs() << std::endl;
+
+  QGauss<dim> quad(3); // Quadrature for coupling
+
+
+  SparsityPattern sparsity;
+  {
+    DynamicSparsityPattern dsp(space_dh.n_dofs(), dh.n_dofs());
+    NonMatching::create_coupling_sparsity_pattern(space_dh, dh,
+                                                  quad, dsp, space_mask);
+    sparsity.copy_from(dsp);
+  }
+  SparseMatrix<double> coupling(sparsity);
+  NonMatching::create_coupling_mass_matrix(space_dh, dh, quad, coupling,
+                                           ConstraintMatrix(), space_mask);
+
+  SparsityPattern mass_sparsity;
+  {
+    DynamicSparsityPattern dsp(dh.n_dofs(), dh.n_dofs());
+    DoFTools::make_sparsity_pattern(dh,dsp);
+    mass_sparsity.copy_from(dsp);
+  }
+  SparseMatrix<double> mass_matrix(mass_sparsity);
+  MatrixTools::create_mass_matrix(dh, quad, mass_matrix);
+
+  SparseDirectUMFPACK mass_matrix_inv;
+  mass_matrix_inv.factorize(mass_matrix);
+
+  // now take ones in space, project them onto the immersed space,
+  // get back ones, and check for the error.
+  Vector<double> space_ones(space_dh.n_dofs());
+  Vector<double> ones(dh.n_dofs());
+
+  space_ones = 1.0;
+  coupling.Tvmult(ones, space_ones);
+  mass_matrix_inv.solve(ones);
+
+  Vector<double> real_ones(dh.n_dofs());
+  real_ones = 1.0;
+  ones -= real_ones;
+
+  deallog << "Error on constants: " << ones.l2_norm() << std::endl;
+}
+
+
+
+int main()
+{
+  initlog();
+  test<1,1>();
+  test<1,2>();
+  test<2,2>();
+  test<2,3>();
+  test<3,3>();
+}
diff --git a/tests/non_matching/coupling_02.output b/tests/non_matching/coupling_02.output
new file mode 100644 (file)
index 0000000..62e0b26
--- /dev/null
@@ -0,0 +1,31 @@
+
+DEAL::dim: 1, spacedim: 1
+DEAL::FE      : FE_Q<1>(1)
+DEAL::Space FE: FESystem<1>[FE_Q<1>(1)^2]
+DEAL::Dofs      : 3
+DEAL::Space dofs: 10
+DEAL::Error on constants: 0.00000
+DEAL::dim: 1, spacedim: 2
+DEAL::FE      : FE_Q<1,2>(1)
+DEAL::Space FE: FESystem<2>[FE_Q<2>(1)^3]
+DEAL::Dofs      : 3
+DEAL::Space dofs: 75
+DEAL::Error on constants: 0.00000
+DEAL::dim: 2, spacedim: 2
+DEAL::FE      : FE_Q<2>(1)
+DEAL::Space FE: FESystem<2>[FE_Q<2>(1)^3]
+DEAL::Dofs      : 9
+DEAL::Space dofs: 75
+DEAL::Error on constants: 1.14304e-15
+DEAL::dim: 2, spacedim: 3
+DEAL::FE      : FE_Q<2,3>(1)
+DEAL::Space FE: FESystem<3>[FE_Q<3>(1)^4]
+DEAL::Dofs      : 9
+DEAL::Space dofs: 500
+DEAL::Error on constants: 7.19507e-16
+DEAL::dim: 3, spacedim: 3
+DEAL::FE      : FE_Q<3>(1)
+DEAL::Space FE: FESystem<3>[FE_Q<3>(1)^4]
+DEAL::Dofs      : 27
+DEAL::Space dofs: 500
+DEAL::Error on constants: 1.07909e-14
diff --git a/tests/non_matching/coupling_03.cc b/tests/non_matching/coupling_03.cc
new file mode 100644 (file)
index 0000000..59c62dc
--- /dev/null
@@ -0,0 +1,134 @@
+// ---------------------------------------------------------------------
+//
+// 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/non_matching/coupling.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparse_direct.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, with a non trivial component selection in both
+// space_dh, and immersed_dh, and check that constants are projected correctly.
+
+template<int dim, int spacedim>
+void test()
+{
+  deallog << "dim: " << dim << ", spacedim: "
+          << spacedim << std::endl;
+
+  Triangulation<dim,spacedim> tria;
+  Triangulation<spacedim,spacedim> space_tria;
+
+  GridGenerator::hyper_cube(tria,-.4,.3);
+  GridGenerator::hyper_cube(space_tria,-1,1);
+
+  tria.refine_global(1);
+  space_tria.refine_global(2);
+
+  FESystem<dim,spacedim> fe(FE_Q<dim,spacedim>(1), spacedim+1);
+  FESystem<spacedim,spacedim> space_fe(FE_Q<spacedim,spacedim>(1), spacedim+1);
+
+  ComponentMask immersed_mask(spacedim+1, false);
+  ComponentMask space_mask(spacedim+1, false);
+
+  immersed_mask.set(0, true);
+  space_mask.set(spacedim, true);
+
+  deallog << "FE      : " << fe.get_name() << std::endl
+          << "Space FE: " << space_fe.get_name() << std::endl;
+
+  DoFHandler<dim,spacedim> dh(tria);
+  DoFHandler<spacedim,spacedim> space_dh(space_tria);
+
+  dh.distribute_dofs(fe);
+  space_dh.distribute_dofs(space_fe);
+
+  deallog << "Dofs      : " << dh.n_dofs() << std::endl
+          << "Space dofs: " << space_dh.n_dofs() << std::endl;
+
+  QGauss<dim> quad(3); // Quadrature for coupling
+
+
+  SparsityPattern sparsity;
+  {
+    DynamicSparsityPattern dsp(space_dh.n_dofs(), dh.n_dofs());
+    NonMatching::create_coupling_sparsity_pattern(space_dh, dh,
+                                                  quad, dsp,
+                                                  space_mask, immersed_mask);
+    sparsity.copy_from(dsp);
+  }
+  SparseMatrix<double> coupling(sparsity);
+  NonMatching::create_coupling_mass_matrix(space_dh, dh, quad, coupling,
+                                           ConstraintMatrix(),
+                                           space_mask, immersed_mask);
+
+  SparsityPattern mass_sparsity;
+  {
+    DynamicSparsityPattern dsp(dh.n_dofs(), dh.n_dofs());
+    DoFTools::make_sparsity_pattern(dh,dsp);
+    mass_sparsity.copy_from(dsp);
+  }
+  SparseMatrix<double> mass_matrix(mass_sparsity);
+  MatrixTools::create_mass_matrix(dh, quad, mass_matrix);
+
+  SparseDirectUMFPACK mass_matrix_inv;
+  mass_matrix_inv.factorize(mass_matrix);
+
+  // now take ones in space, project them onto the immersed space,
+  // get back ones, and check for the error.
+  Vector<double> space_ones(space_dh.n_dofs());
+  Vector<double> ones(dh.n_dofs());
+
+  space_ones = 1.0;
+  coupling.Tvmult(ones, space_ones);
+  mass_matrix_inv.solve(ones);
+
+  // Only the first component is different from zero in the result, since the
+  // matrix couples only the first component of the immersed space, with the
+  // last component of the embedding space.
+  Vector<double> real_ones(dh.n_dofs());
+  for (unsigned int i=0; i<dh.n_dofs()/(spacedim+1); ++i)
+    real_ones((spacedim+1)*i) = 1.0;
+
+  ones -= real_ones;
+
+  deallog << "Error on constants: " << ones.l2_norm() << std::endl;
+}
+
+
+
+int main()
+{
+  initlog();
+  test<1,1>();
+  test<1,2>();
+  test<2,2>();
+  test<2,3>();
+  test<3,3>();
+}
diff --git a/tests/non_matching/coupling_03.output b/tests/non_matching/coupling_03.output
new file mode 100644 (file)
index 0000000..0f51e2b
--- /dev/null
@@ -0,0 +1,31 @@
+
+DEAL::dim: 1, spacedim: 1
+DEAL::FE      : FESystem<1>[FE_Q<1>(1)^2]
+DEAL::Space FE: FESystem<1>[FE_Q<1>(1)^2]
+DEAL::Dofs      : 6
+DEAL::Space dofs: 10
+DEAL::Error on constants: 0.00000
+DEAL::dim: 1, spacedim: 2
+DEAL::FE      : FESystem<1,2>[FE_Q<1,2>(1)^3]
+DEAL::Space FE: FESystem<2>[FE_Q<2>(1)^3]
+DEAL::Dofs      : 9
+DEAL::Space dofs: 75
+DEAL::Error on constants: 0.00000
+DEAL::dim: 2, spacedim: 2
+DEAL::FE      : FESystem<2>[FE_Q<2>(1)^3]
+DEAL::Space FE: FESystem<2>[FE_Q<2>(1)^3]
+DEAL::Dofs      : 27
+DEAL::Space dofs: 75
+DEAL::Error on constants: 1.14304e-15
+DEAL::dim: 2, spacedim: 3
+DEAL::FE      : FESystem<2,3>[FE_Q<2,3>(1)^4]
+DEAL::Space FE: FESystem<3>[FE_Q<3>(1)^4]
+DEAL::Dofs      : 36
+DEAL::Space dofs: 500
+DEAL::Error on constants: 7.19507e-16
+DEAL::dim: 3, spacedim: 3
+DEAL::FE      : FESystem<3>[FE_Q<3>(1)^4]
+DEAL::Space FE: FESystem<3>[FE_Q<3>(1)^4]
+DEAL::Dofs      : 108
+DEAL::Space dofs: 500
+DEAL::Error on constants: 1.07909e-14
diff --git a/tests/non_matching/coupling_04.cc b/tests/non_matching/coupling_04.cc
new file mode 100644 (file)
index 0000000..463efe1
--- /dev/null
@@ -0,0 +1,126 @@
+// ---------------------------------------------------------------------
+//
+// 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/non_matching/coupling.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparse_direct.h>
+
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/base/function_lib.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 quadratic functions are correctly
+// projected.
+
+template<int dim, int spacedim>
+void test()
+{
+  deallog << "dim: " << dim << ", spacedim: "
+          << spacedim << std::endl;
+
+  Triangulation<dim,spacedim> tria;
+  Triangulation<spacedim,spacedim> space_tria;
+
+  GridGenerator::hyper_cube(tria,-.4,.3);
+  GridGenerator::hyper_cube(space_tria,-1,1);
+
+  tria.refine_global(1);
+  space_tria.refine_global(2);
+
+  FE_Q<dim,spacedim> fe(2);
+  FE_Q<spacedim,spacedim> space_fe(2);
+
+  deallog << "FE      : " << fe.get_name() << std::endl
+          << "Space FE: " << space_fe.get_name() << std::endl;
+
+  DoFHandler<dim,spacedim> dh(tria);
+  DoFHandler<spacedim,spacedim> space_dh(space_tria);
+
+  dh.distribute_dofs(fe);
+  space_dh.distribute_dofs(space_fe);
+
+  deallog << "Dofs      : " << dh.n_dofs() << std::endl
+          << "Space dofs: " << space_dh.n_dofs() << std::endl;
+
+  QGauss<dim> quad(3); // Quadrature for coupling
+
+
+  SparsityPattern sparsity;
+  {
+    DynamicSparsityPattern dsp(space_dh.n_dofs(), dh.n_dofs());
+    NonMatching::create_coupling_sparsity_pattern(space_dh, dh,
+                                                  quad, dsp);
+    sparsity.copy_from(dsp);
+  }
+  SparseMatrix<double> coupling(sparsity);
+  NonMatching::create_coupling_mass_matrix(space_dh, dh, quad, coupling,
+                                           ConstraintMatrix());
+
+  SparsityPattern mass_sparsity;
+  {
+    DynamicSparsityPattern dsp(dh.n_dofs(), dh.n_dofs());
+    DoFTools::make_sparsity_pattern(dh,dsp);
+    mass_sparsity.copy_from(dsp);
+  }
+  SparseMatrix<double> mass_matrix(mass_sparsity);
+  MatrixTools::create_mass_matrix(dh, quad, mass_matrix);
+
+  SparseDirectUMFPACK mass_matrix_inv;
+  mass_matrix_inv.factorize(mass_matrix);
+
+  // now take the square function in space, project them onto the immersed space,
+  // get back ones, and check for the error.
+  Vector<double> space_square(space_dh.n_dofs());
+  Vector<double> squares(dh.n_dofs());
+  Vector<double> projected_squares(dh.n_dofs());
+
+  VectorTools::interpolate(space_dh, Functions::SquareFunction<spacedim>(),
+                           space_square);
+  VectorTools::interpolate(dh, Functions::SquareFunction<spacedim>(),
+                           squares);
+
+  coupling.Tvmult(projected_squares, space_square);
+  mass_matrix_inv.solve(projected_squares);
+
+  projected_squares -= squares;
+
+  deallog << "Error on squares: " << projected_squares.l2_norm() << std::endl;
+}
+
+
+
+int main()
+{
+  initlog();
+  test<1,1>();
+  test<1,2>();
+  test<2,2>();
+  test<2,3>();
+  test<3,3>();
+}
diff --git a/tests/non_matching/coupling_04.output b/tests/non_matching/coupling_04.output
new file mode 100644 (file)
index 0000000..2ae7ba5
--- /dev/null
@@ -0,0 +1,31 @@
+
+DEAL::dim: 1, spacedim: 1
+DEAL::FE      : FE_Q<1>(2)
+DEAL::Space FE: FE_Q<1>(2)
+DEAL::Dofs      : 5
+DEAL::Space dofs: 9
+DEAL::Error on squares: 6.24154e-17
+DEAL::dim: 1, spacedim: 2
+DEAL::FE      : FE_Q<1,2>(2)
+DEAL::Space FE: FE_Q<2>(2)
+DEAL::Dofs      : 5
+DEAL::Space dofs: 81
+DEAL::Error on squares: 6.24154e-17
+DEAL::dim: 2, spacedim: 2
+DEAL::FE      : FE_Q<2>(2)
+DEAL::Space FE: FE_Q<2>(2)
+DEAL::Dofs      : 25
+DEAL::Space dofs: 81
+DEAL::Error on squares: 3.00905e-16
+DEAL::dim: 2, spacedim: 3
+DEAL::FE      : FE_Q<2,3>(2)
+DEAL::Space FE: FE_Q<3>(2)
+DEAL::Dofs      : 25
+DEAL::Space dofs: 729
+DEAL::Error on squares: 3.54580e-16
+DEAL::dim: 3, spacedim: 3
+DEAL::FE      : FE_Q<3>(2)
+DEAL::Space FE: FE_Q<3>(2)
+DEAL::Dofs      : 125
+DEAL::Space dofs: 729
+DEAL::Error on squares: 5.36692e-15

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.