From c776a9a9a829e67ee892ef5a0bbea754ebce9512 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Fri, 18 Feb 2022 09:44:01 +0000 Subject: [PATCH] Fixed NonMatching::create_coupling_* when constraints are present. --- doc/news/changes/minor/20220218LucaHeltai_b | 5 + include/deal.II/non_matching/coupling.h | 16 +- source/non_matching/coupling.cc | 29 ++- source/non_matching/coupling.inst.in | 12 +- tests/non_matching/coupling_09.cc | 188 ++++++++++++++++++ .../coupling_09.with_umfpack=on.output | 41 ++++ 6 files changed, 272 insertions(+), 19 deletions(-) create mode 100644 doc/news/changes/minor/20220218LucaHeltai_b create mode 100644 tests/non_matching/coupling_09.cc create mode 100644 tests/non_matching/coupling_09.with_umfpack=on.output diff --git a/doc/news/changes/minor/20220218LucaHeltai_b b/doc/news/changes/minor/20220218LucaHeltai_b new file mode 100644 index 0000000000..61fd4875b8 --- /dev/null +++ b/doc/news/changes/minor/20220218LucaHeltai_b @@ -0,0 +1,5 @@ +Fixed: NonMatching::create_coupling_sparsity_pattern() and +NonMatching::create_coupling_mass_matrix() were not working properly with +hanging node constraints. +
+(Luca Heltai, 2022/02/18) diff --git a/include/deal.II/non_matching/coupling.h b/include/deal.II/non_matching/coupling.h index ec345d06c7..cf75588419 100644 --- a/include/deal.II/non_matching/coupling.h +++ b/include/deal.II/non_matching/coupling.h @@ -107,7 +107,9 @@ namespace NonMatching const Mapping & space_mapping = StaticMappingQ1::mapping, const Mapping &immersed_mapping = - StaticMappingQ1::mapping); + StaticMappingQ1::mapping, + const AffineConstraints &immersed_constraints = + AffineConstraints()); /** * 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 & immersed_mapping = - StaticMappingQ1::mapping); + StaticMappingQ1::mapping, + const AffineConstraints &immersed_constraints = + AffineConstraints()); /** @@ -192,7 +196,9 @@ namespace NonMatching const Mapping &space_mapping = StaticMappingQ1::mapping, const Mapping &immersed_mapping = - StaticMappingQ1::mapping); + StaticMappingQ1::mapping, + const AffineConstraints &immersed_constraints = + AffineConstraints()); /** * 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 &immersed_mapping = - StaticMappingQ1::mapping); + StaticMappingQ1::mapping, + const AffineConstraints &immersed_constraints = + AffineConstraints()); /** * Create a coupling sparsity pattern for non-matching independent grids, diff --git a/source/non_matching/coupling.cc b/source/non_matching/coupling.cc index 44a2b3b05e..294dd244ec 100644 --- a/source/non_matching/coupling.cc +++ b/source/non_matching/coupling.cc @@ -182,7 +182,8 @@ namespace NonMatching const ComponentMask & space_comps, const ComponentMask & immersed_comps, const Mapping & space_mapping, - const Mapping & immersed_mapping) + const Mapping & immersed_mapping, + const AffineConstraints & immersed_constraints) { GridTools::Cache 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 & constraints, const ComponentMask & space_comps, const ComponentMask & immersed_comps, - const Mapping & immersed_mapping) + const Mapping & immersed_mapping, + const AffineConstraints & 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 & space_mapping, - const Mapping & immersed_mapping) + const Mapping & immersed_mapping, + const AffineConstraints &immersed_constraints) { GridTools::Cache 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 &constraints, const ComponentMask & space_comps, const ComponentMask & immersed_comps, - const Mapping & immersed_mapping) + const Mapping & immersed_mapping, + const AffineConstraints &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); } } } diff --git a/source/non_matching/coupling.inst.in b/source/non_matching/coupling.inst.in index 25daf080c0..ff650ddd3b 100644 --- a/source/non_matching/coupling.inst.in +++ b/source/non_matching/coupling.inst.in @@ -28,7 +28,8 @@ for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS; const ComponentMask & space_comps, const ComponentMask & immersed_comps, const Mapping & space_mapping, - const Mapping & immersed_mapping); + const Mapping & immersed_mapping, + const AffineConstraints & immersed_constraints); template void create_coupling_sparsity_pattern( const GridTools::Cache &cache, @@ -39,7 +40,8 @@ for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS; const AffineConstraints & constraints, const ComponentMask & space_comps, const ComponentMask & immersed_comps, - const Mapping & immersed_mapping); + const Mapping & immersed_mapping, + const AffineConstraints & 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 & space_mapping, - const Mapping & immersed_mapping); + const Mapping & immersed_mapping, + const AffineConstraints &immersed_constraints); template void create_coupling_mass_matrix( const GridTools::Cache & cache, @@ -80,7 +83,8 @@ for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS; const AffineConstraints &constraints, const ComponentMask & space_comps, const ComponentMask & immersed_comps, - const Mapping & immersed_mapping); + const Mapping & immersed_mapping, + const AffineConstraints &immersed_constraints); template void create_coupling_mass_matrix( Functions::CutOffFunctionBase & kernel, diff --git a/tests/non_matching/coupling_09.cc b/tests/non_matching/coupling_09.cc new file mode 100644 index 0000000000..8f8034b45f --- /dev/null +++ b/tests/non_matching/coupling_09.cc @@ -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 +#include +#include + +#include + +#include + +#include +#include + +#include +#include +#include +#include + +// #include "/workspace/dealii/dealii-0/source/non_matching/coupling.cc" +#include + +#include + +#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 +void +test() +{ + static_assert(dim0 >= dim1, "Only dim0>=dim1 is possible"); + + deallog << "dim0: " << dim0 << ", dim1: " << dim1 + << ", spacedim: " << spacedim << std::endl; + + Triangulation tria0; + Triangulation 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 fe0(1); + FE_Q fe1(1); + + DoFHandler dh0(tria0); + DoFHandler dh1(tria1); + + GridTools::Cache cache0(tria0); + GridTools::Cache cache1(tria1); + + dh0.distribute_dofs(fe0); + dh1.distribute_dofs(fe1); + + AffineConstraints constraints0; + AffineConstraints 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 quad0(2); // Quadrature for coupling + QGauss 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::mapping, + constraints1); + sparsity.copy_from(dsp); + } + SparseMatrix coupling(sparsity); + + NonMatching::create_coupling_mass_matrix( + cache0, + dh0, + dh1, + quad1, + coupling, + constraints0, + ComponentMask(), + ComponentMask(), + StaticMappingQ1::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 mass_matrix1(mass_sparsity1); + MatrixTools::create_mass_matrix(dh1, + quad1, + mass_matrix1, + static_cast *>( + 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 ones0(dh0.n_dofs()); + Vector ones1(dh1.n_dofs()); + + ones0 = 1.0; + constraints0.distribute(ones0); + + coupling.Tvmult(ones1, ones0); + mass_matrix1_inv.solve(ones1); + constraints1.distribute(ones1); + + Vector 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 index 0000000000..b4af94323d --- /dev/null +++ b/tests/non_matching/coupling_09.with_umfpack=on.output @@ -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 -- 2.39.5