From c3512c64e7f5d7079939fedb61463024fd240e8d Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Wed, 2 May 2018 20:45:39 +0200 Subject: [PATCH] Moved Cache outside compute_coupling_mass/sparsity. --- include/deal.II/non_matching/coupling.h | 45 ++++- source/non_matching/coupling.cc | 76 ++++++-- source/non_matching/coupling.inst.in | 26 ++- tests/non_matching/coupling_02.cc | 5 +- tests/non_matching/coupling_03.cc | 5 +- tests/non_matching/coupling_06.cc | 177 ++++++++++++++++++ ...ling_06.with_trilinos=true,mpirun=2.output | 83 ++++++++ 7 files changed, 392 insertions(+), 25 deletions(-) create mode 100644 tests/non_matching/coupling_06.cc create mode 100644 tests/non_matching/coupling_06.with_trilinos=true,mpirun=2.output diff --git a/include/deal.II/non_matching/coupling.h b/include/deal.II/non_matching/coupling.h index 51590f97af..a7c08685a1 100644 --- a/include/deal.II/non_matching/coupling.h +++ b/include/deal.II/non_matching/coupling.h @@ -19,6 +19,8 @@ #include #include +#include + #include #include #include @@ -89,14 +91,34 @@ namespace NonMatching const DoFHandler &immersed_dh, const Quadrature &quad, Sparsity &sparsity, + 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); /** - * Create a coupling mass matrix for non-matching, overlapping grids. + * Same as above, but takes an additional GridTools::Cache object, instead of + * creating one internally. In this version of the function, the parameter @p + * space_mapping cannot be specified, since it is taken from the @p cache + * parameter. * + * @author Luca Heltai, 2018 + */ + template + void create_coupling_sparsity_pattern(const GridTools::Cache &cache, + const DoFHandler &space_dh, + const DoFHandler &immersed_dh, + const Quadrature &quad, + Sparsity &sparsity, + const ConstraintMatrix &constraints = ConstraintMatrix(), + const ComponentMask &space_comps = ComponentMask(), + const ComponentMask &immersed_comps = ComponentMask(), + 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 @@ -124,7 +146,7 @@ namespace NonMatching * * If the domain $B$ does not fall within $\Omega$, an exception will be * thrown by the algorithm that computes the quadrature point locations. In - * particular, notice that this function only makes sens for `dim1` lower or + * particular, notice that this function only makes sense for `dim1` lower or * equal than `dim0`. A static assert guards that this is actually the case. * * For both spaces, it is possible to specify a custom Mapping, which @@ -147,6 +169,25 @@ namespace NonMatching const ComponentMask &immersed_comps = ComponentMask(), const Mapping &space_mapping = StaticMappingQ1::mapping, const Mapping &immersed_mapping = StaticMappingQ1::mapping); + + /** + * Same as above, but takes an additional GridTools::Cache object, instead of + * creating one internally. In this version of the function, the parameter @p + * space_mapping cannot specified, since it is taken from the @p cache + * parameter. + * + * @author Luca Heltai, 2018 + */ + template + void create_coupling_mass_matrix(const GridTools::Cache &cache, + 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 &immersed_mapping = StaticMappingQ1::mapping); } DEAL_II_NAMESPACE_CLOSE diff --git a/source/non_matching/coupling.cc b/source/non_matching/coupling.cc index bd4ae01cc6..2ed5670d25 100644 --- a/source/non_matching/coupling.cc +++ b/source/non_matching/coupling.cc @@ -49,10 +49,31 @@ namespace NonMatching const DoFHandler &immersed_dh, const Quadrature &quad, Sparsity &sparsity, + const ConstraintMatrix &constraints, const ComponentMask &space_comps, const ComponentMask &immersed_comps, const Mapping &space_mapping, const Mapping &immersed_mapping) + { + GridTools::Cache cache(space_dh.get_triangulation(), space_mapping); + create_coupling_sparsity_pattern(cache, space_dh, immersed_dh, + quad, sparsity, constraints, + space_comps, immersed_comps, + immersed_mapping); + } + + + + template + void create_coupling_sparsity_pattern(const GridTools::Cache &cache, + const DoFHandler &space_dh, + const DoFHandler &immersed_dh, + const Quadrature &quad, + Sparsity &sparsity, + const ConstraintMatrix &constraints, + const ComponentMask &space_comps, + const ComponentMask &immersed_comps, + const Mapping &immersed_mapping) { AssertDimension(sparsity.n_rows(), space_dh.n_dofs()); AssertDimension(sparsity.n_cols(), immersed_dh.n_dofs()); @@ -75,8 +96,6 @@ namespace NonMatching FEValues fe_v(immersed_mapping, immersed_fe, quad, update_quadrature_points); - GridTools::Cache cache(space_dh.get_triangulation(), space_mapping); - // Take care of components const ComponentMask space_c = (space_comps.size() == 0 ? @@ -102,6 +121,22 @@ namespace NonMatching if (immersed_c[i]) immersed_gtl[i] = j++; + // Construct a dof_mask, used to distribute entries to the sparsity + Table< 2, bool > dof_mask(space_fe.dofs_per_cell, + immersed_fe.dofs_per_cell); + dof_mask.fill(false); + for (unsigned int i=0; iis_locally_owned()) { ocell->get_dof_indices(odofs); - for (unsigned int i=0; i &space_mapping, const Mapping &immersed_mapping) + { + GridTools::Cache cache(space_dh.get_triangulation(), space_mapping); + create_coupling_mass_matrix(cache, space_dh, immersed_dh, + quad, matrix, constraints, space_comps, immersed_comps, + immersed_mapping); + } + + + + template + void create_coupling_mass_matrix(const GridTools::Cache &cache, + 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 &immersed_mapping) { AssertDimension(matrix.m(), space_dh.n_dofs()); AssertDimension(matrix.n(), immersed_dh.n_dofs()); @@ -167,8 +209,6 @@ namespace NonMatching 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 const ComponentMask space_c = (space_comps.size() == 0 ? @@ -232,7 +272,7 @@ namespace NonMatching 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, + FEValues o_fe_v(cache.get_mapping(), space_dh.get_fe(), qps, update_values); o_fe_v.reinit(ocell); ocell->get_dof_indices(odofs); diff --git a/source/non_matching/coupling.inst.in b/source/non_matching/coupling.inst.in index be535179ac..06eb4ad1b3 100644 --- a/source/non_matching/coupling.inst.in +++ b/source/non_matching/coupling.inst.in @@ -18,14 +18,27 @@ for (dim0 : DIMENSIONS; dim1 :DIMENSIONS; spacedim : SPACE_DIMENSIONS; Sparsity : SPARSITY_PATTERNS) { #if dim1 <= dim0 && dim0 <= spacedim - template void create_coupling_sparsity_pattern(const DoFHandler &space_dh, + template + void create_coupling_sparsity_pattern(const DoFHandler &space_dh, const DoFHandler &immersed_dh, const Quadrature &quad, Sparsity &sparsity, + const ConstraintMatrix &constraints, const ComponentMask &space_comps, const ComponentMask &immersed_comps, const Mapping &space_mapping, const Mapping &immersed_mapping); + + template + void create_coupling_sparsity_pattern(const GridTools::Cache &cache, + const DoFHandler &space_dh, + const DoFHandler &immersed_dh, + const Quadrature &quad, + Sparsity &sparsity, + const ConstraintMatrix &constraints, + const ComponentMask &space_comps, + const ComponentMask &immersed_comps, + const Mapping &immersed_mapping); #endif } @@ -44,5 +57,16 @@ for (dim0 : DIMENSIONS; dim1 :DIMENSIONS; spacedim : SPACE_DIMENSIONS; const ComponentMask &immersed_comps, const Mapping &space_mapping, const Mapping &immersed_mapping); + + template + void create_coupling_mass_matrix(const GridTools::Cache &cache, + 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 &immersed_mapping); #endif } diff --git a/tests/non_matching/coupling_02.cc b/tests/non_matching/coupling_02.cc index 9e954284e2..02e4b79182 100644 --- a/tests/non_matching/coupling_02.cc +++ b/tests/non_matching/coupling_02.cc @@ -70,17 +70,18 @@ void test() QGauss quad(3); // Quadrature for coupling + ConstraintMatrix constraints; SparsityPattern sparsity; { DynamicSparsityPattern dsp(space_dh.n_dofs(), dh.n_dofs()); NonMatching::create_coupling_sparsity_pattern(space_dh, dh, - quad, dsp, space_mask); + quad, dsp, constraints, space_mask); sparsity.copy_from(dsp); } SparseMatrix coupling(sparsity); NonMatching::create_coupling_mass_matrix(space_dh, dh, quad, coupling, - ConstraintMatrix(), space_mask); + constraints, space_mask); SparsityPattern mass_sparsity; { diff --git a/tests/non_matching/coupling_03.cc b/tests/non_matching/coupling_03.cc index 59c62dc11f..df61e50414 100644 --- a/tests/non_matching/coupling_03.cc +++ b/tests/non_matching/coupling_03.cc @@ -74,18 +74,19 @@ void test() QGauss quad(3); // Quadrature for coupling + ConstraintMatrix constraints; SparsityPattern sparsity; { DynamicSparsityPattern dsp(space_dh.n_dofs(), dh.n_dofs()); NonMatching::create_coupling_sparsity_pattern(space_dh, dh, - quad, dsp, + quad, dsp, constraints, space_mask, immersed_mask); sparsity.copy_from(dsp); } SparseMatrix coupling(sparsity); NonMatching::create_coupling_mass_matrix(space_dh, dh, quad, coupling, - ConstraintMatrix(), + constraints, space_mask, immersed_mask); SparsityPattern mass_sparsity; diff --git a/tests/non_matching/coupling_06.cc b/tests/non_matching/coupling_06.cc new file mode 100644 index 0000000000..fb88bdf582 --- /dev/null +++ b/tests/non_matching/coupling_06.cc @@ -0,0 +1,177 @@ +// --------------------------------------------------------------------- +// +// 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 +#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; + + const auto &comm = MPI_COMM_WORLD; + + parallel::shared::Triangulation tria(comm); + parallel::distributed::Triangulation space_tria(comm); + + GridGenerator::hyper_cube(tria,-.4,.3); + GridGenerator::hyper_cube(space_tria,-1,1); + + tria.refine_global(3); + space_tria.refine_global(3); + + 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); + + auto space_locally_owned_dofs = space_dh.locally_owned_dofs (); + auto locally_owned_dofs = dh.locally_owned_dofs (); + + + deallog << "Dofs : " << dh.n_dofs() << std::endl + << "Space dofs: " << space_dh.n_dofs() << std::endl; + + deallog << "Local dofs : " << locally_owned_dofs.n_elements() + << std::endl; + deallog << "Local space dofs: " << space_locally_owned_dofs.n_elements() + << std::endl; + QGauss quad(3); // Quadrature for coupling + + GridTools::Cache cache(space_tria); + + TrilinosWrappers::SparsityPattern sparsity(space_locally_owned_dofs, + locally_owned_dofs, comm); + NonMatching::create_coupling_sparsity_pattern(cache, space_dh, dh, + quad, sparsity); + sparsity.compress(); + + TrilinosWrappers::SparseMatrix coupling(sparsity); + NonMatching::create_coupling_mass_matrix(cache, space_dh, dh, quad, coupling); + coupling.compress(VectorOperation::add); + + TrilinosWrappers::SparsityPattern mass_sparsity(locally_owned_dofs, comm); + DoFTools::make_sparsity_pattern(dh,mass_sparsity); + mass_sparsity.compress(); + + TrilinosWrappers::SparseMatrix mass_matrix(mass_sparsity); + + { + QGauss quad(4); + FEValues fev(fe, quad, update_values|update_JxW_values); + std::vector dofs(fe.dofs_per_cell); + FullMatrix cell_matrix(fe.dofs_per_cell, fe.dofs_per_cell); + ConstraintMatrix constraints; + + for (auto cell : dh.active_cell_iterators()) + if (cell->is_locally_owned()) + { + cell_matrix = 0; + fev.reinit(cell); + cell->get_dof_indices(dofs); + for (unsigned int i=0; i(), + space_square); + VectorTools::interpolate(dh, Functions::SquareFunction(), + squares); + + coupling.Tvmult(Mprojected_squares, space_square); + + SolverControl cn(100, 1e-12); + TrilinosWrappers::SolverCG solver(cn); + TrilinosWrappers::PreconditionILU prec; + prec.initialize(mass_matrix); + + solver.solve(mass_matrix, projected_squares, Mprojected_squares, prec); + + deallog << "Squares norm : " << projected_squares.l2_norm() << std::endl; + + projected_squares -= squares; + + deallog << "Error on squares: " << projected_squares.l2_norm() << std::endl; +} + + + +int main(int argc, char **argv) +{ + auto init = Utilities::MPI::MPI_InitFinalize(argc, argv, 1); + MPILogInitAll log(true); + test<1,2>(); + test<2,2>(); + test<2,3>(); + test<3,3>(); +} diff --git a/tests/non_matching/coupling_06.with_trilinos=true,mpirun=2.output b/tests/non_matching/coupling_06.with_trilinos=true,mpirun=2.output new file mode 100644 index 0000000000..0f3760084c --- /dev/null +++ b/tests/non_matching/coupling_06.with_trilinos=true,mpirun=2.output @@ -0,0 +1,83 @@ + +DEAL:0::dim: 1, spacedim: 2 +DEAL:0::FE : FE_Q<1,2>(2) +DEAL:0::Space FE: FE_Q<2>(2) +DEAL:0::Dofs : 17 +DEAL:0::Space dofs: 289 +DEAL:0::Local dofs : 9 +DEAL:0::Local space dofs: 153 +DEAL:0::Convergence step 9 value 4.51115e-14 +DEAL:0::Squares norm : 0.275853 +DEAL:0::Error on squares: 1.07673e-12 +DEAL:0::dim: 2, spacedim: 2 +DEAL:0::FE : FE_Q<2>(2) +DEAL:0::Space FE: FE_Q<2>(2) +DEAL:0::Dofs : 289 +DEAL:0::Space dofs: 289 +DEAL:0::Local dofs : 153 +DEAL:0::Local space dofs: 153 +DEAL:0::Convergence step 12 value 2.66570e-13 +DEAL:0::Squares norm : 1.98578 +DEAL:0::Error on squares: 4.16199e-10 +DEAL:0::dim: 2, spacedim: 3 +DEAL:0::FE : FE_Q<2,3>(2) +DEAL:0::Space FE: FE_Q<3>(2) +DEAL:0::Dofs : 289 +DEAL:0::Space dofs: 4913 +DEAL:0::Local dofs : 153 +DEAL:0::Local space dofs: 2601 +DEAL:0::Convergence step 12 value 2.66570e-13 +DEAL:0::Squares norm : 1.98578 +DEAL:0::Error on squares: 4.16199e-10 +DEAL:0::dim: 3, spacedim: 3 +DEAL:0::FE : FE_Q<3>(2) +DEAL:0::Space FE: FE_Q<3>(2) +DEAL:0::Dofs : 4913 +DEAL:0::Space dofs: 4913 +DEAL:0::Local dofs : 2601 +DEAL:0::Local space dofs: 2601 +DEAL:0::Convergence step 12 value 6.65914e-13 +DEAL:0::Squares norm : 11.6248 +DEAL:0::Error on squares: 4.61425e-08 + +DEAL:1::dim: 1, spacedim: 2 +DEAL:1::FE : FE_Q<1,2>(2) +DEAL:1::Space FE: FE_Q<2>(2) +DEAL:1::Dofs : 17 +DEAL:1::Space dofs: 289 +DEAL:1::Local dofs : 8 +DEAL:1::Local space dofs: 136 +DEAL:1::Convergence step 9 value 4.51115e-14 +DEAL:1::Squares norm : 0.275853 +DEAL:1::Error on squares: 1.07673e-12 +DEAL:1::dim: 2, spacedim: 2 +DEAL:1::FE : FE_Q<2>(2) +DEAL:1::Space FE: FE_Q<2>(2) +DEAL:1::Dofs : 289 +DEAL:1::Space dofs: 289 +DEAL:1::Local dofs : 136 +DEAL:1::Local space dofs: 136 +DEAL:1::Convergence step 12 value 2.66570e-13 +DEAL:1::Squares norm : 1.98578 +DEAL:1::Error on squares: 4.16199e-10 +DEAL:1::dim: 2, spacedim: 3 +DEAL:1::FE : FE_Q<2,3>(2) +DEAL:1::Space FE: FE_Q<3>(2) +DEAL:1::Dofs : 289 +DEAL:1::Space dofs: 4913 +DEAL:1::Local dofs : 136 +DEAL:1::Local space dofs: 2312 +DEAL:1::Convergence step 12 value 2.66570e-13 +DEAL:1::Squares norm : 1.98578 +DEAL:1::Error on squares: 4.16199e-10 +DEAL:1::dim: 3, spacedim: 3 +DEAL:1::FE : FE_Q<3>(2) +DEAL:1::Space FE: FE_Q<3>(2) +DEAL:1::Dofs : 4913 +DEAL:1::Space dofs: 4913 +DEAL:1::Local dofs : 2312 +DEAL:1::Local space dofs: 2312 +DEAL:1::Convergence step 12 value 6.65914e-13 +DEAL:1::Squares norm : 11.6248 +DEAL:1::Error on squares: 4.61425e-08 + -- 2.39.5