#include <deal.II/base/config.h>
#include <deal.II/base/quadrature.h>
+#include <deal.II/grid/grid_tools_cache.h>
+
#include <deal.II/fe/component_mask.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/fe/fe.h>
const DoFHandler<dim1, spacedim> &immersed_dh,
const Quadrature<dim1> &quad,
Sparsity &sparsity,
+ const ConstraintMatrix &constraints = ConstraintMatrix(),
const ComponentMask &space_comps = ComponentMask(),
const ComponentMask &immersed_comps = ComponentMask(),
const Mapping<dim0, spacedim> &space_mapping = StaticMappingQ1<dim0,spacedim>::mapping,
const Mapping<dim1, spacedim> &immersed_mapping = StaticMappingQ1<dim1, spacedim>::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<int dim0, int dim1, int spacedim, typename Sparsity>
+ void create_coupling_sparsity_pattern(const GridTools::Cache<dim0, spacedim> &cache,
+ const DoFHandler<dim0, spacedim> &space_dh,
+ const DoFHandler<dim1, spacedim> &immersed_dh,
+ const Quadrature<dim1> &quad,
+ Sparsity &sparsity,
+ const ConstraintMatrix &constraints = ConstraintMatrix(),
+ const ComponentMask &space_comps = ComponentMask(),
+ const ComponentMask &immersed_comps = ComponentMask(),
+ const Mapping<dim1, spacedim> &immersed_mapping = StaticMappingQ1<dim1, spacedim>::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
*
* 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
const ComponentMask &immersed_comps = ComponentMask(),
const Mapping<dim0, spacedim> &space_mapping = StaticMappingQ1<dim0,spacedim>::mapping,
const Mapping<dim1, spacedim> &immersed_mapping = StaticMappingQ1<dim1, spacedim>::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<int dim0, int dim1, int spacedim, typename Matrix>
+ void create_coupling_mass_matrix(const GridTools::Cache<dim0, spacedim> &cache,
+ const DoFHandler<dim0, spacedim> &space_dh,
+ const DoFHandler<dim1, spacedim> &immersed_dh,
+ const Quadrature<dim1> &quad,
+ Matrix &matrix,
+ const ConstraintMatrix &constraints = ConstraintMatrix(),
+ const ComponentMask &space_comps = ComponentMask(),
+ const ComponentMask &immersed_comps = ComponentMask(),
+ const Mapping<dim1, spacedim> &immersed_mapping = StaticMappingQ1<dim1, spacedim>::mapping);
}
DEAL_II_NAMESPACE_CLOSE
const DoFHandler<dim1, spacedim> &immersed_dh,
const Quadrature<dim1> &quad,
Sparsity &sparsity,
+ const ConstraintMatrix &constraints,
const ComponentMask &space_comps,
const ComponentMask &immersed_comps,
const Mapping<dim0, spacedim> &space_mapping,
const Mapping<dim1, spacedim> &immersed_mapping)
+ {
+ GridTools::Cache<dim0,spacedim> 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<int dim0, int dim1, int spacedim, typename Sparsity>
+ void create_coupling_sparsity_pattern(const GridTools::Cache<dim0,spacedim> &cache,
+ const DoFHandler<dim0, spacedim> &space_dh,
+ const DoFHandler<dim1, spacedim> &immersed_dh,
+ const Quadrature<dim1> &quad,
+ Sparsity &sparsity,
+ const ConstraintMatrix &constraints,
+ const ComponentMask &space_comps,
+ const ComponentMask &immersed_comps,
+ const Mapping<dim1, spacedim> &immersed_mapping)
{
AssertDimension(sparsity.n_rows(), space_dh.n_dofs());
AssertDimension(sparsity.n_cols(), immersed_dh.n_dofs());
FEValues<dim1,spacedim> fe_v(immersed_mapping, immersed_fe, quad,
update_quadrature_points);
- GridTools::Cache<dim0,spacedim> cache(space_dh.get_triangulation(), space_mapping);
-
// Take care of components
const ComponentMask space_c
= (space_comps.size() == 0 ?
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; i<space_fe.dofs_per_cell; ++i)
+ {
+ const auto comp_i = space_fe.system_to_component_index(i).first;
+ if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
+ for (unsigned int j=0; j<immersed_fe.dofs_per_cell; ++j)
+ {
+ const auto comp_j = immersed_fe.system_to_component_index(j).first;
+ if (immersed_gtl[comp_j] == space_gtl[comp_i])
+ dof_mask(i,j) = true;
+ }
+ }
+
for (; cell != endc; ++cell)
{
// Reinitialize the cell and the fe_values
if (ocell->is_locally_owned())
{
ocell->get_dof_indices(odofs);
- for (unsigned int i=0; i<odofs.size(); ++i)
- {
- const auto comp_i = space_fe.system_to_component_index(i).first;
- if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
- {
- for (unsigned int j=0; j<dofs.size(); ++j)
- {
- const auto comp_j = immersed_fe.system_to_component_index(j).first;
- if (immersed_gtl[comp_j] == space_gtl[comp_i])
- sparsity.add(odofs[i],dofs[j]);
- }
- }
- }
+ constraints.add_entries_local_to_global(odofs, dofs, sparsity);
}
}
}
const ComponentMask &immersed_comps,
const Mapping<dim0, spacedim> &space_mapping,
const Mapping<dim1, spacedim> &immersed_mapping)
+ {
+ GridTools::Cache<dim0,spacedim> 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<int dim0, int dim1, int spacedim, typename Matrix>
+ void create_coupling_mass_matrix(const GridTools::Cache<dim0, spacedim> &cache,
+ const DoFHandler<dim0, spacedim> &space_dh,
+ const DoFHandler<dim1, spacedim> &immersed_dh,
+ const Quadrature<dim1> &quad,
+ Matrix &matrix,
+ const ConstraintMatrix &constraints,
+ const ComponentMask &space_comps,
+ const ComponentMask &immersed_comps,
+ const Mapping<dim1, spacedim> &immersed_mapping)
{
AssertDimension(matrix.m(), space_dh.n_dofs());
AssertDimension(matrix.n(), immersed_dh.n_dofs());
std::vector<types::global_dof_index> dofs(immersed_fe.dofs_per_cell);
std::vector<types::global_dof_index> odofs(space_fe.dofs_per_cell);
- GridTools::Cache<dim0,spacedim> cache(space_dh.get_triangulation(), space_mapping);
-
// Take care of components
const ComponentMask space_c
= (space_comps.size() == 0 ?
const std::vector< Point<dim0> > &qps = qpoints[c];
const std::vector< unsigned int > &ids = maps[c];
- FEValues<dim0,spacedim> o_fe_v(space_mapping, space_dh.get_fe(), qps,
+ FEValues<dim0,spacedim> o_fe_v(cache.get_mapping(), space_dh.get_fe(), qps,
update_values);
o_fe_v.reinit(ocell);
ocell->get_dof_indices(odofs);
Sparsity : SPARSITY_PATTERNS)
{
#if dim1 <= dim0 && dim0 <= spacedim
- template void create_coupling_sparsity_pattern(const DoFHandler<dim0, spacedim> &space_dh,
+ template
+ void create_coupling_sparsity_pattern(const DoFHandler<dim0, spacedim> &space_dh,
const DoFHandler<dim1, spacedim> &immersed_dh,
const Quadrature<dim1> &quad,
Sparsity &sparsity,
+ const ConstraintMatrix &constraints,
const ComponentMask &space_comps,
const ComponentMask &immersed_comps,
const Mapping<dim0, spacedim> &space_mapping,
const Mapping<dim1, spacedim> &immersed_mapping);
+
+ template
+ void create_coupling_sparsity_pattern(const GridTools::Cache<dim0,spacedim> &cache,
+ const DoFHandler<dim0, spacedim> &space_dh,
+ const DoFHandler<dim1, spacedim> &immersed_dh,
+ const Quadrature<dim1> &quad,
+ Sparsity &sparsity,
+ const ConstraintMatrix &constraints,
+ const ComponentMask &space_comps,
+ const ComponentMask &immersed_comps,
+ const Mapping<dim1, spacedim> &immersed_mapping);
#endif
}
const ComponentMask &immersed_comps,
const Mapping<dim0, spacedim> &space_mapping,
const Mapping<dim1, spacedim> &immersed_mapping);
+
+ template
+ void create_coupling_mass_matrix(const GridTools::Cache<dim0,spacedim> &cache,
+ const DoFHandler<dim0, spacedim> &space_dh,
+ const DoFHandler<dim1, spacedim> &immersed_dh,
+ const Quadrature<dim1> &quad,
+ Matrix &matrix,
+ const ConstraintMatrix &constraints,
+ const ComponentMask &space_comps,
+ const ComponentMask &immersed_comps,
+ const Mapping<dim1, spacedim> &immersed_mapping);
#endif
}
QGauss<dim> 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<double> coupling(sparsity);
NonMatching::create_coupling_mass_matrix(space_dh, dh, quad, coupling,
- ConstraintMatrix(), space_mask);
+ constraints, space_mask);
SparsityPattern mass_sparsity;
{
QGauss<dim> 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<double> coupling(sparsity);
NonMatching::create_coupling_mass_matrix(space_dh, dh, quad, coupling,
- ConstraintMatrix(),
+ constraints,
space_mask, immersed_mask);
SparsityPattern mass_sparsity;
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/shared_tria.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools_cache.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparse_direct.h>
+
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/base/function_lib.h>
+
+#include <deal.II/base/mpi.h>
+
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/lac/trilinos_precondition.h>
+
+#include <deal.II/non_matching/coupling.h>
+
+#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<int dim, int spacedim>
+void test()
+{
+ deallog << "dim: " << dim << ", spacedim: "
+ << spacedim << std::endl;
+
+ const auto &comm = MPI_COMM_WORLD;
+
+ parallel::shared::Triangulation<dim,spacedim> tria(comm);
+ parallel::distributed::Triangulation<spacedim,spacedim> 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<dim,spacedim> fe(2);
+ FE_Q<spacedim,spacedim> space_fe(2);
+
+ deallog << "FE : " << fe.get_name() << std::endl
+ << "Space FE: " << space_fe.get_name() << std::endl;
+
+ DoFHandler<dim,spacedim> dh(tria);
+ DoFHandler<spacedim,spacedim> 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<dim> quad(3); // Quadrature for coupling
+
+ GridTools::Cache<spacedim,spacedim> 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<dim> quad(4);
+ FEValues<dim,spacedim> fev(fe, quad, update_values|update_JxW_values);
+ std::vector<types::global_dof_index> dofs(fe.dofs_per_cell);
+ FullMatrix<double> 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<fe.dofs_per_cell; ++i)
+ for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ for (unsigned int q=0; q<quad.size(); ++q)
+ cell_matrix(i,j) +=
+ fev.shape_value(i,q)*
+ fev.shape_value(j,q)*
+ fev.JxW(q);
+ constraints.distribute_local_to_global(cell_matrix, dofs, mass_matrix);
+ }
+ mass_matrix.compress(VectorOperation::add);
+ }
+
+ // now take the square function in space, project them onto the immersed space,
+ // get back ones, and check for the error.
+ TrilinosWrappers::MPI::Vector space_square(space_locally_owned_dofs, comm);
+ TrilinosWrappers::MPI::Vector squares(locally_owned_dofs, comm);
+ TrilinosWrappers::MPI::Vector Mprojected_squares(locally_owned_dofs, comm);
+ TrilinosWrappers::MPI::Vector projected_squares(locally_owned_dofs, comm);
+
+ VectorTools::interpolate(space_dh, Functions::SquareFunction<spacedim>(),
+ space_square);
+ VectorTools::interpolate(dh, Functions::SquareFunction<spacedim>(),
+ 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>();
+}
--- /dev/null
+
+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
+
projects / dealii.git / commitdiff