From 267eb68850623c219518f7a9d5ec510dd64d3551 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Thu, 11 Apr 2019 15:53:05 +0200 Subject: [PATCH] Added ConvolutionCoupling --- doc/news/changes/minor/20190419LucaHeltai | 5 + include/deal.II/non_matching/coupling.h | 135 +++++- source/non_matching/coupling.cc | 451 +++++++++++++++++- source/non_matching/coupling.inst.in | 31 +- tests/non_matching/coupling_08.cc | 162 +++++++ .../coupling_08.with_umfpack=true.output | 0 6 files changed, 775 insertions(+), 9 deletions(-) create mode 100644 doc/news/changes/minor/20190419LucaHeltai create mode 100644 tests/non_matching/coupling_08.cc create mode 100644 tests/non_matching/coupling_08.with_umfpack=true.output diff --git a/doc/news/changes/minor/20190419LucaHeltai b/doc/news/changes/minor/20190419LucaHeltai new file mode 100644 index 0000000000..5b726acc2b --- /dev/null +++ b/doc/news/changes/minor/20190419LucaHeltai @@ -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. +
+(Luca Heltai, 2019/04/19) diff --git a/include/deal.II/non_matching/coupling.h b/include/deal.II/non_matching/coupling.h index 056c0d67ed..1894e625fc 100644 --- a/include/deal.II/non_matching/coupling.h +++ b/include/deal.II/non_matching/coupling.h @@ -18,12 +18,14 @@ #include +#include #include #include #include #include +#include #include #include @@ -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 &immersed_mapping = StaticMappingQ1::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. + * An exception is thrown if both triagnulations are of type + * parallel::distributed::Triangulation. + * + * 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 + void + create_coupling_sparsity_pattern( + const double & epsilon, + const GridTools::Cache &cache0, + const GridTools::Cache &cache1, + const DoFHandler & dh0, + const DoFHandler & dh1, + Sparsity & sparsity, + const AffineConstraints &constraints0 = AffineConstraints(), + 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. + * An exception is thrown if both triangulations are of type + * parallel::distributed::Triangulation. + * + * 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 + void + create_coupling_mass_matrix( + Functions::CutOffFunctionBase & kernel, + const double & epsilon, + const GridTools::Cache & cache0, + const GridTools::Cache & cache1, + const DoFHandler & dh0, + const DoFHandler & dh1, + const Quadrature & quadrature0, + const Quadrature & quadrature1, + Matrix & matrix, + const AffineConstraints &constraints0 = + AffineConstraints(), + const ComponentMask &comps0 = ComponentMask(), + const ComponentMask &comps1 = ComponentMask()); } // namespace NonMatching DEAL_II_NAMESPACE_CLOSE diff --git a/source/non_matching/coupling.cc b/source/non_matching/coupling.cc index 4df28b9fa2..16c34f830f 100644 --- a/source/non_matching/coupling.cc +++ b/source/non_matching/coupling.cc @@ -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 + std::pair, std::vector> + compute_components_coupling(const ComponentMask & comps0, + const ComponentMask & comps1, + const FiniteElement &fe0, + const FiniteElement &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 gtl0(fe0.n_components(), + numbers::invalid_unsigned_int); + std::vector 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 *>( &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 *>( &immersed_dh.get_triangulation()) == nullptr), @@ -572,8 +614,411 @@ namespace NonMatching } } + template + void + create_coupling_sparsity_pattern( + const double & epsilon, + const GridTools::Cache &cache0, + const GridTools::Cache &cache1, + const DoFHandler & dh0, + const DoFHandler & dh1, + Sparsity & sparsity, + const AffineConstraints & 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(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 + *>(&dh0.get_triangulation()) != nullptr); + const bool one_is_distributed = + (dynamic_cast + *>(&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 dofs0(fe0.dofs_per_cell); + std::vector 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, + typename Triangulation::active_cell_iterator>> + intersection; + + for (const auto &cell0 : dh0.active_cell_iterators()) + if (cell0->is_locally_owned()) + { + intersection.resize(0); + + BoundingBox box0; + if (cache0.get_mapping().preserves_vertex_locations()) + box0 = cell0->bounding_box(); + else + { + const auto vertices = + cache0.get_mapping().get_vertices(cell0); + Point 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({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::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, + typename Triangulation::active_cell_iterator>> + intersection; + + for (const auto &cell1 : dh1.active_cell_iterators()) + if (cell1->is_locally_owned()) + { + intersection.resize(0); + + BoundingBox 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 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({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::cell_iterator cell0( + *entry.second, &dh0); + cell0->get_dof_indices(dofs0); + constraints0.add_entries_local_to_global(dofs0, + dofs1, + sparsity); + } + } + } + } + } + + + + template + void + create_coupling_mass_matrix( + Functions::CutOffFunctionBase & kernel, + const double & epsilon, + const GridTools::Cache & cache0, + const GridTools::Cache & cache1, + const DoFHandler & dh0, + const DoFHandler & dh1, + const Quadrature & quadrature0, + const Quadrature & quadrature1, + Matrix & matrix, + const AffineConstraints &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 + *>(&dh0.get_triangulation()) != nullptr); + const bool one_is_distributed = + (dynamic_cast + *>(&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 fev0(cache0.get_mapping(), + fe0, + quadrature0, + update_values | update_JxW_values | + update_quadrature_points); + + FEValues fev1(cache1.get_mapping(), + fe1, + quadrature1, + update_values | update_JxW_values | + update_quadrature_points); + + // Dof indices + std::vector dofs0(fe0.dofs_per_cell); + std::vector dofs1(fe1.dofs_per_cell); + + // Local Matrix + FullMatrix 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 >l0 = p.first; + const auto >l1 = p.second; + + kernel.set_radius(epsilon); + std::vector 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, + typename Triangulation::active_cell_iterator>> + intersection; + + for (const auto &cell0 : dh0.active_cell_iterators()) + if (cell0->is_locally_owned()) + { + intersection.resize(0); + + BoundingBox box0; + if (cache0.get_mapping().preserves_vertex_locations()) + box0 = cell0->bounding_box(); + else + { + const auto vertices = + cache0.get_mapping().get_vertices(cell0); + Point 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({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::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, + typename Triangulation::active_cell_iterator>> + intersection; + + for (const auto &cell1 : dh1.active_cell_iterators()) + if (cell1->is_locally_owned()) + { + intersection.resize(0); + BoundingBox box1; + if (cache1.get_mapping().preserves_vertex_locations()) + box1 = cell1->bounding_box(); + else + { + const auto vertices = + cache1.get_mapping().get_vertices(cell1); + Point 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({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::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 diff --git a/source/non_matching/coupling.inst.in b/source/non_matching/coupling.inst.in index c1dcd9e574..cfc2022cb5 100644 --- a/source/non_matching/coupling.inst.in +++ b/source/non_matching/coupling.inst.in @@ -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 &space_dh, const DoFHandler &immersed_dh, @@ -40,6 +40,17 @@ for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS; const ComponentMask & space_comps, const ComponentMask & immersed_comps, const Mapping & immersed_mapping); + + template void create_coupling_sparsity_pattern( + const double & epsilon, + const GridTools::Cache &cache0, + const GridTools::Cache &cache1, + const DoFHandler & dh0, + const DoFHandler & dh1, + Sparsity & sparsity, + const AffineConstraints & 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 & space_dh, const DoFHandler & immersed_dh, @@ -69,5 +80,19 @@ for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS; const ComponentMask & space_comps, const ComponentMask & immersed_comps, const Mapping & immersed_mapping); + + template void create_coupling_mass_matrix( + Functions::CutOffFunctionBase & kernel, + const double & epsilon, + const GridTools::Cache & cache0, + const GridTools::Cache & cache1, + const DoFHandler & dh0, + const DoFHandler & dh1, + const Quadrature & quadrature0, + const Quadrature & quadrature1, + Matrix & matrix, + const AffineConstraints &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 index 0000000000..eedf9b1dab --- /dev/null +++ b/tests/non_matching/coupling_08.cc @@ -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 +#include +#include + +#include + +#include + +#include +#include + +#include +#include +#include +#include + +#include + +#include + +#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 +void +test() +{ + deallog << "dim0: " << dim0 << ", dim1: " << dim1 + << ", spacedim: " << spacedim << std::endl; + + Triangulation tria0; + Triangulation tria1; + + GridGenerator::hyper_cube(tria0, -.4, .3); + GridGenerator::hyper_cube(tria1, -1, 1); + + tria0.refine_global(3); + tria1.refine_global(2); + + FE_Q fe0(1); + FE_Q fe1(1); + + deallog << "FE0 : " << fe0.get_name() << std::endl + << "FE1 : " << fe1.get_name() << std::endl; + + DoFHandler dh0(tria0); + DoFHandler dh1(tria1); + + GridTools::Cache cache0(tria0); + GridTools::Cache cache1(tria1); + + dh0.distribute_dofs(fe0); + dh1.distribute_dofs(fe1); + + AffineConstraints 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 quad0(3); // Quadrature for coupling + QGauss 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 coupling(sparsity); + + Functions::CutOffFunctionC1 dirac( + 1, + Point(), + 1, + Functions::CutOffFunctionBase::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 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 ones0(dh0.n_dofs()); + Vector ones1(dh1.n_dofs()); + + ones0 = 1.0; + coupling.Tvmult(ones1, ones0); + mass_matrix1_inv.solve(ones1); + + Vector 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 index 0000000000..e69de29bb2 -- 2.39.5