From 0c9ed5bb6382710bd73d1c4a15bd284cd9cdc6f5 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Mon, 5 Mar 2018 15:34:20 +0100 Subject: [PATCH] Non matching sparsity and mass matrix. --- include/deal.II/non_matching/coupling.h | 348 ++++++++++++++++++++++++ tests/non_matching/coupling_01.cc | 118 ++++++++ tests/non_matching/coupling_01.output | 31 +++ tests/non_matching/coupling_02.cc | 123 +++++++++ tests/non_matching/coupling_02.output | 31 +++ tests/non_matching/coupling_03.cc | 134 +++++++++ tests/non_matching/coupling_03.output | 31 +++ tests/non_matching/coupling_04.cc | 126 +++++++++ tests/non_matching/coupling_04.output | 31 +++ 9 files changed, 973 insertions(+) create mode 100644 include/deal.II/non_matching/coupling.h create mode 100644 tests/non_matching/coupling_01.cc create mode 100644 tests/non_matching/coupling_01.output create mode 100644 tests/non_matching/coupling_02.cc create mode 100644 tests/non_matching/coupling_02.output create mode 100644 tests/non_matching/coupling_03.cc create mode 100644 tests/non_matching/coupling_03.output create mode 100644 tests/non_matching/coupling_04.cc create mode 100644 tests/non_matching/coupling_04.output diff --git a/include/deal.II/non_matching/coupling.h b/include/deal.II/non_matching/coupling.h new file mode 100644 index 0000000000..3b8be3e477 --- /dev/null +++ b/include/deal.II/non_matching/coupling.h @@ -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 +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + + + +#include + +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 + void create_coupling_sparsity_pattern(const DoFHandler &space_dh, + const DoFHandler &immersed_dh, + const Quadrature &quad, + Sparsity &sparsity, + const ComponentMask &space_comps = ComponentMask(), + const ComponentMask &immersed_comps = ComponentMask(), + const Mapping &space_mapping = StaticMappingQ1::mapping, + const Mapping &immersed_mapping = StaticMappingQ1::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 + void create_coupling_mass_matrix(const DoFHandler &space_dh, + const DoFHandler &immersed_dh, + const Quadrature &quad, + Matrix &matrix, + const ConstraintMatrix &constraints = ConstraintMatrix(), + const ComponentMask &space_comps = ComponentMask(), + const ComponentMask &immersed_comps = ComponentMask(), + const Mapping &space_mapping = StaticMappingQ1::mapping, + const Mapping &immersed_mapping = StaticMappingQ1::mapping); + + + + +// === inline and template functions === + + + + template + void create_coupling_sparsity_pattern(const DoFHandler &space_dh, + const DoFHandler &immersed_dh, + const Quadrature &quad, + Sparsity &sparsity, + const ComponentMask &space_comps, + const ComponentMask &immersed_comps, + const Mapping &space_mapping, + const Mapping &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::active_cell_iterator + cell = immersed_dh.begin_active(), + endc = immersed_dh.end(); + + // Dof indices + std::vector dofs(immersed_fe.dofs_per_cell); + std::vector odofs(space_fe.dofs_per_cell); + + FEValues fe_v(immersed_mapping, immersed_fe, quad, + update_quadrature_points); + + GridTools::Cache 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 space_gtl(space_fe.n_components(), numbers::invalid_unsigned_int); + std::vector immersed_gtl(immersed_fe.n_components(), numbers::invalid_unsigned_int); + + for (unsigned int i=0, j=0; isubdomain_id() == this_mpi_process) + { + // Reinitialize the cell and the fe_values + fe_v.reinit(cell); + cell->get_dof_indices(dofs); + + const std::vector > &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::cell_iterator + ocell(*cells[c], &space_dh); + ocell->get_dof_indices(odofs); + for (unsigned int i=0; i + void create_coupling_mass_matrix(const DoFHandler &space_dh, + const DoFHandler &immersed_dh, + const Quadrature &quad, + Matrix &matrix, + const ConstraintMatrix &constraints, + const ComponentMask &space_comps, + const ComponentMask &immersed_comps, + const Mapping &space_mapping, + const Mapping &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 dofs(immersed_fe.dofs_per_cell); + std::vector odofs(space_fe.dofs_per_cell); + + GridTools::Cache 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 space_gtl(space_fe.n_components(), numbers::invalid_unsigned_int); + std::vector immersed_gtl(immersed_fe.n_components(), numbers::invalid_unsigned_int); + + for (unsigned int i=0, j=0; i cell_matrix(space_dh.get_fe().dofs_per_cell, + immersed_dh.get_fe().dofs_per_cell); + + FEValues 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::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 > &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::active_cell_iterator + ocell(*cells[c], &space_dh); + const std::vector< Point > &qps = qpoints[c]; + const std::vector< unsigned int > &ids = maps[c]; + + FEValues 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 +#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. + +template +void test() +{ + deallog << "dim: " << dim << ", spacedim: " + << spacedim << std::endl; + + Triangulation tria; + Triangulation 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 fe(1); + FE_Q space_fe(1); + + deallog << "FE : " << fe.get_name() << std::endl + << "Space FE: " << space_fe.get_name() << std::endl; + + DoFHandler dh(tria); + DoFHandler 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 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 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 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 space_ones(space_dh.n_dofs()); + Vector ones(dh.n_dofs()); + + space_ones = 1.0; + coupling.Tvmult(ones, space_ones); + mass_matrix_inv.solve(ones); + + Vector 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 index 0000000000..db8ef199f2 --- /dev/null +++ b/tests/non_matching/coupling_01.output @@ -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 index 0000000000..9e954284e2 --- /dev/null +++ b/tests/non_matching/coupling_02.cc @@ -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 +#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, with a non trivial component selection in the space_dh, +// and check that constants are projected correctly. + +template +void test() +{ + deallog << "dim: " << dim << ", spacedim: " + << spacedim << std::endl; + + Triangulation tria; + Triangulation 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 fe(1); + FESystem space_fe(FE_Q(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 dh(tria); + DoFHandler 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 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 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 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 space_ones(space_dh.n_dofs()); + Vector ones(dh.n_dofs()); + + space_ones = 1.0; + coupling.Tvmult(ones, space_ones); + mass_matrix_inv.solve(ones); + + Vector 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 index 0000000000..62e0b26508 --- /dev/null +++ b/tests/non_matching/coupling_02.output @@ -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 index 0000000000..59c62dc11f --- /dev/null +++ b/tests/non_matching/coupling_03.cc @@ -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 +#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, with a non trivial component selection in both +// space_dh, and immersed_dh, and check that constants are projected correctly. + +template +void test() +{ + deallog << "dim: " << dim << ", spacedim: " + << spacedim << std::endl; + + Triangulation tria; + Triangulation 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 fe(FE_Q(1), spacedim+1); + FESystem space_fe(FE_Q(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 dh(tria); + DoFHandler 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 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 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 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 space_ones(space_dh.n_dofs()); + Vector 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 real_ones(dh.n_dofs()); + for (unsigned int i=0; i(); + 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 index 0000000000..0f51e2bb3a --- /dev/null +++ b/tests/non_matching/coupling_03.output @@ -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 index 0000000000..463efe1102 --- /dev/null +++ b/tests/non_matching/coupling_04.cc @@ -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 +#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 quadratic functions are correctly +// projected. + +template +void test() +{ + deallog << "dim: " << dim << ", spacedim: " + << spacedim << std::endl; + + Triangulation tria; + Triangulation 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 fe(2); + FE_Q space_fe(2); + + deallog << "FE : " << fe.get_name() << std::endl + << "Space FE: " << space_fe.get_name() << std::endl; + + DoFHandler dh(tria); + DoFHandler 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 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 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 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 space_square(space_dh.n_dofs()); + Vector squares(dh.n_dofs()); + Vector projected_squares(dh.n_dofs()); + + VectorTools::interpolate(space_dh, Functions::SquareFunction(), + space_square); + VectorTools::interpolate(dh, Functions::SquareFunction(), + 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 index 0000000000..2ae7ba5607 --- /dev/null +++ b/tests/non_matching/coupling_04.output @@ -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 -- 2.39.5