]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fixed NonMatching::create_coupling_* when constraints are present. 13416/head
authorLuca Heltai <luca.heltai@sissa.it>
Fri, 18 Feb 2022 09:44:01 +0000 (09:44 +0000)
committerLuca Heltai <luca.heltai@sissa.it>
Sat, 19 Feb 2022 07:53:25 +0000 (07:53 +0000)
doc/news/changes/minor/20220218LucaHeltai_b [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_09.cc [new file with mode: 0644]
tests/non_matching/coupling_09.with_umfpack=on.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20220218LucaHeltai_b b/doc/news/changes/minor/20220218LucaHeltai_b
new file mode 100644 (file)
index 0000000..61fd487
--- /dev/null
@@ -0,0 +1,5 @@
+Fixed: NonMatching::create_coupling_sparsity_pattern() and
+NonMatching::create_coupling_mass_matrix() were not working properly with
+hanging node constraints.
+<br>
+(Luca Heltai, 2022/02/18)
index ec345d06c76b18093fd97ab80955090b6f4e6459..cf75588419951356b4f896e8283b8a0f6d5d335e 100644 (file)
@@ -107,7 +107,9 @@ namespace NonMatching
     const Mapping<dim0, spacedim> &   space_mapping =
       StaticMappingQ1<dim0, spacedim>::mapping,
     const Mapping<dim1, spacedim> &immersed_mapping =
-      StaticMappingQ1<dim1, spacedim>::mapping);
+      StaticMappingQ1<dim1, spacedim>::mapping,
+    const AffineConstraints<number> &immersed_constraints =
+      AffineConstraints<number>());
 
   /**
    * Same as above, but takes an additional GridTools::Cache object, instead of
@@ -131,7 +133,9 @@ namespace NonMatching
     const ComponentMask &            space_comps = ComponentMask(),
     const ComponentMask &            immersed_comps = ComponentMask(),
     const Mapping<dim1, spacedim> &  immersed_mapping =
-      StaticMappingQ1<dim1, spacedim>::mapping);
+      StaticMappingQ1<dim1, spacedim>::mapping,
+    const AffineConstraints<number> &immersed_constraints =
+      AffineConstraints<number>());
 
 
   /**
@@ -192,7 +196,9 @@ namespace NonMatching
     const Mapping<dim0, spacedim> &space_mapping =
       StaticMappingQ1<dim0, spacedim>::mapping,
     const Mapping<dim1, spacedim> &immersed_mapping =
-      StaticMappingQ1<dim1, spacedim>::mapping);
+      StaticMappingQ1<dim1, spacedim>::mapping,
+    const AffineConstraints<typename Matrix::value_type> &immersed_constraints =
+      AffineConstraints<typename Matrix::value_type>());
 
   /**
    * Same as above, but takes an additional GridTools::Cache object, instead of
@@ -213,7 +219,9 @@ namespace NonMatching
     const ComponentMask &          space_comps    = ComponentMask(),
     const ComponentMask &          immersed_comps = ComponentMask(),
     const Mapping<dim1, spacedim> &immersed_mapping =
-      StaticMappingQ1<dim1, spacedim>::mapping);
+      StaticMappingQ1<dim1, spacedim>::mapping,
+    const AffineConstraints<typename Matrix::value_type> &immersed_constraints =
+      AffineConstraints<typename Matrix::value_type>());
 
   /**
    * Create a coupling sparsity pattern for non-matching independent grids,
index 44a2b3b05ed001aa5bc14e30e7e0180b410732f6..294dd244eccfd012f576d5591e7937fa03f34148 100644 (file)
@@ -182,7 +182,8 @@ namespace NonMatching
     const ComponentMask &             space_comps,
     const ComponentMask &             immersed_comps,
     const Mapping<dim0, spacedim> &   space_mapping,
-    const Mapping<dim1, spacedim> &   immersed_mapping)
+    const Mapping<dim1, spacedim> &   immersed_mapping,
+    const AffineConstraints<number> & immersed_constraints)
   {
     GridTools::Cache<dim0, spacedim> cache(space_dh.get_triangulation(),
                                            space_mapping);
@@ -194,7 +195,8 @@ namespace NonMatching
                                      constraints,
                                      space_comps,
                                      immersed_comps,
-                                     immersed_mapping);
+                                     immersed_mapping,
+                                     immersed_constraints);
   }
 
 
@@ -214,7 +216,8 @@ namespace NonMatching
     const AffineConstraints<number> &       constraints,
     const ComponentMask &                   space_comps,
     const ComponentMask &                   immersed_comps,
-    const Mapping<dim1, spacedim> &         immersed_mapping)
+    const Mapping<dim1, spacedim> &         immersed_mapping,
+    const AffineConstraints<number> &       immersed_constraints)
   {
     AssertDimension(sparsity.n_rows(), space_dh.n_dofs());
     AssertDimension(sparsity.n_cols(), immersed_dh.n_dofs());
@@ -350,7 +353,10 @@ namespace NonMatching
                 // for the case of non-trivial dof_mask, we should
                 // uncomment the missing part.
                 constraints.add_entries_local_to_global(
-                  odofs, dofs, sparsity); //, true, dof_mask);
+                  odofs,
+                  immersed_constraints,
+                  dofs,
+                  sparsity); //, true, dof_mask);
               }
           }
         ++i;
@@ -370,7 +376,8 @@ namespace NonMatching
     const ComponentMask &                                 space_comps,
     const ComponentMask &                                 immersed_comps,
     const Mapping<dim0, spacedim> &                       space_mapping,
-    const Mapping<dim1, spacedim> &                       immersed_mapping)
+    const Mapping<dim1, spacedim> &                       immersed_mapping,
+    const AffineConstraints<typename Matrix::value_type> &immersed_constraints)
   {
     GridTools::Cache<dim0, spacedim> cache(space_dh.get_triangulation(),
                                            space_mapping);
@@ -382,7 +389,8 @@ namespace NonMatching
                                 constraints,
                                 space_comps,
                                 immersed_comps,
-                                immersed_mapping);
+                                immersed_mapping,
+                                immersed_constraints);
   }
 
 
@@ -398,7 +406,8 @@ namespace NonMatching
     const AffineConstraints<typename Matrix::value_type> &constraints,
     const ComponentMask &                                 space_comps,
     const ComponentMask &                                 immersed_comps,
-    const Mapping<dim1, spacedim> &                       immersed_mapping)
+    const Mapping<dim1, spacedim> &                       immersed_mapping,
+    const AffineConstraints<typename Matrix::value_type> &immersed_constraints)
   {
     AssertDimension(matrix.m(), space_dh.n_dofs());
     AssertDimension(matrix.n(), immersed_dh.n_dofs());
@@ -608,10 +617,8 @@ namespace NonMatching
                   }
 
                 // Now assemble the matrices
-                constraints.distribute_local_to_global(cell_matrix,
-                                                       odofs,
-                                                       dofs,
-                                                       matrix);
+                constraints.distribute_local_to_global(
+                  cell_matrix, odofs, immersed_constraints, dofs, matrix);
               }
           }
       }
index 25daf080c09b0915a3826ea4339aac97f4a3bd7f..ff650ddd3bb00975ee96e0c865b6f6b6fb06d48a 100644 (file)
@@ -28,7 +28,8 @@ for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS;
       const ComponentMask &             space_comps,
       const ComponentMask &             immersed_comps,
       const Mapping<dim0, spacedim> &   space_mapping,
-      const Mapping<dim1, spacedim> &   immersed_mapping);
+      const Mapping<dim1, spacedim> &   immersed_mapping,
+      const AffineConstraints<S> &      immersed_constraints);
 
     template void create_coupling_sparsity_pattern(
       const GridTools::Cache<dim0, spacedim> &cache,
@@ -39,7 +40,8 @@ for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS;
       const AffineConstraints<S> &            constraints,
       const ComponentMask &                   space_comps,
       const ComponentMask &                   immersed_comps,
-      const Mapping<dim1, spacedim> &         immersed_mapping);
+      const Mapping<dim1, spacedim> &         immersed_mapping,
+      const AffineConstraints<S> &            immersed_constraints);
 
     template void create_coupling_sparsity_pattern(
       const double &                          epsilon,
@@ -69,7 +71,8 @@ for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS;
       const ComponentMask &                        space_comps,
       const ComponentMask &                        immersed_comps,
       const Mapping<dim0, spacedim> &              space_mapping,
-      const Mapping<dim1, spacedim> &              immersed_mapping);
+      const Mapping<dim1, spacedim> &              immersed_mapping,
+      const AffineConstraints<Matrix::value_type> &immersed_constraints);
 
     template void create_coupling_mass_matrix(
       const GridTools::Cache<dim0, spacedim> &     cache,
@@ -80,7 +83,8 @@ for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS;
       const AffineConstraints<Matrix::value_type> &constraints,
       const ComponentMask &                        space_comps,
       const ComponentMask &                        immersed_comps,
-      const Mapping<dim1, spacedim> &              immersed_mapping);
+      const Mapping<dim1, spacedim> &              immersed_mapping,
+      const AffineConstraints<Matrix::value_type> &immersed_constraints);
 
     template void create_coupling_mass_matrix(
       Functions::CutOffFunctionBase<spacedim> & kernel,
diff --git a/tests/non_matching/coupling_09.cc b/tests/non_matching/coupling_09.cc
new file mode 100644 (file)
index 0000000..8f8034b
--- /dev/null
@@ -0,0 +1,188 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#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 "/workspace/dealii/dealii-0/source/non_matching/coupling.cc"
+#include <deal.II/non_matching/coupling.h>
+
+#include <deal.II/numerics/matrix_tools.h>
+
+#include "../tests.h"
+
+// 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()
+{
+  static_assert(dim0 >= dim1, "Only dim0>=dim1 is possible");
+
+  deallog << "dim0: " << dim0 << ", dim1: " << dim1
+          << ", spacedim: " << spacedim << std::endl;
+
+  Triangulation<dim0, spacedim> tria0;
+  Triangulation<dim1, spacedim> tria1;
+
+  GridGenerator::hyper_cube(tria0, -1, 1);
+  GridGenerator::hyper_cube(tria1, -.44444, .3333);
+
+  tria0.refine_global((dim0 < 3 ? 3 : 2));
+  tria1.refine_global(2);
+
+  for (unsigned int i = 0; i < 2; ++i)
+    {
+      for (const auto &cell : tria1.active_cell_iterators())
+        if (cell->center()[0] < 0)
+          cell->set_refine_flag();
+      tria1.execute_coarsening_and_refinement();
+
+      for (const auto &cell : tria0.active_cell_iterators())
+        if (cell->center()[0] > 0)
+          cell->set_refine_flag();
+      tria0.execute_coarsening_and_refinement();
+    }
+
+  FE_Q<dim0, spacedim> fe0(1);
+  FE_Q<dim1, spacedim> fe1(1);
+
+  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;
+  AffineConstraints<double> constraints1;
+
+  DoFTools::make_hanging_node_constraints(dh0, constraints0);
+  DoFTools::make_zero_boundary_constraints(dh0, constraints0);
+
+  DoFTools::make_hanging_node_constraints(dh1, constraints1);
+
+  constraints0.close();
+  constraints1.close();
+
+  deallog << "FE0                : " << fe0.get_name() << std::endl
+          << "FE1                : " << fe1.get_name() << std::endl
+          << "Dofs 0             : " << dh0.n_dofs() << std::endl
+          << "Dofs 1             : " << dh1.n_dofs() << std::endl
+          << "Constrained dofs 0 : " << constraints0.n_constraints()
+          << std::endl
+          << "Constrained dofs 1 : " << constraints1.n_constraints()
+          << std::endl;
+
+  QGauss<dim0> quad0(2); // Quadrature for coupling
+  QGauss<dim1> quad1(2); // Quadrature for coupling
+
+  SparsityPattern sparsity;
+  {
+    DynamicSparsityPattern dsp(dh0.n_dofs(), dh1.n_dofs());
+    NonMatching::create_coupling_sparsity_pattern(
+      cache0,
+      dh0,
+      dh1,
+      quad1,
+      dsp,
+      constraints0,
+      ComponentMask(),
+      ComponentMask(),
+      StaticMappingQ1<dim1, spacedim>::mapping,
+      constraints1);
+    sparsity.copy_from(dsp);
+  }
+  SparseMatrix<double> coupling(sparsity);
+
+  NonMatching::create_coupling_mass_matrix(
+    cache0,
+    dh0,
+    dh1,
+    quad1,
+    coupling,
+    constraints0,
+    ComponentMask(),
+    ComponentMask(),
+    StaticMappingQ1<dim1, spacedim>::mapping,
+    constraints1);
+
+  SparsityPattern mass_sparsity1;
+  {
+    DynamicSparsityPattern dsp(dh1.n_dofs(), dh1.n_dofs());
+    DoFTools::make_sparsity_pattern(dh1, dsp, constraints1, false);
+    mass_sparsity1.copy_from(dsp);
+  }
+  SparseMatrix<double> mass_matrix1(mass_sparsity1);
+  MatrixTools::create_mass_matrix(dh1,
+                                  quad1,
+                                  mass_matrix1,
+                                  static_cast<const Function<spacedim> *>(
+                                    nullptr),
+                                  constraints1);
+
+  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;
+  constraints0.distribute(ones0);
+
+  coupling.Tvmult(ones1, ones0);
+  mass_matrix1_inv.solve(ones1);
+  constraints1.distribute(ones1);
+
+  Vector<double> exact_ones1(dh1.n_dofs());
+  exact_ones1 = 1.0;
+  constraints1.distribute(exact_ones1);
+  ones1 -= exact_ones1;
+
+  deallog << "Error on constants: " << ones1.l2_norm() << std::endl;
+}
+
+int
+main()
+{
+  initlog(true);
+  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_09.with_umfpack=on.output b/tests/non_matching/coupling_09.with_umfpack=on.output
new file mode 100644 (file)
index 0000000..b4af943
--- /dev/null
@@ -0,0 +1,41 @@
+
+DEAL::dim0: 1, dim1: 1, spacedim: 1
+DEAL::FE0                : FE_Q<1>(1)
+DEAL::FE1                : FE_Q<1>(1)
+DEAL::Dofs 0             : 21
+DEAL::Dofs 1             : 11
+DEAL::Constrained dofs 0 : 2
+DEAL::Constrained dofs 1 : 0
+DEAL::Error on constants: 4.71028e-16
+DEAL::dim0: 2, dim1: 1, spacedim: 2
+DEAL::FE0                : FE_Q<2>(1)
+DEAL::FE1                : FE_Q<1,2>(1)
+DEAL::Dofs 0             : 622
+DEAL::Dofs 1             : 11
+DEAL::Constrained dofs 0 : 106
+DEAL::Constrained dofs 1 : 0
+DEAL::Error on constants: 3.68219e-16
+DEAL::dim0: 2, dim1: 2, spacedim: 2
+DEAL::FE0                : FE_Q<2>(1)
+DEAL::FE1                : FE_Q<2>(1)
+DEAL::Dofs 0             : 622
+DEAL::Dofs 1             : 176
+DEAL::Constrained dofs 0 : 106
+DEAL::Constrained dofs 1 : 12
+DEAL::Error on constants: 6.70652e-15
+DEAL::dim0: 3, dim1: 2, spacedim: 3
+DEAL::FE0                : FE_Q<3>(1)
+DEAL::FE1                : FE_Q<2,3>(1)
+DEAL::Dofs 0             : 2788
+DEAL::Dofs 1             : 176
+DEAL::Constrained dofs 0 : 1106
+DEAL::Constrained dofs 1 : 12
+DEAL::Error on constants: 6.95996e-15
+DEAL::dim0: 3, dim1: 3, spacedim: 3
+DEAL::FE0                : FE_Q<3>(1)
+DEAL::FE1                : FE_Q<3>(1)
+DEAL::Dofs 0             : 2788
+DEAL::Dofs 1             : 2788
+DEAL::Constrained dofs 0 : 1106
+DEAL::Constrained dofs 1 : 264
+DEAL::Error on constants: 6.31637e-14

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.