--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/config.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/quadrature.h>
+
+#include <deal.II/grid/grid_tools.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>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+
+
+
+#include <vector>
+
+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<int dim0, int dim1, int spacedim, typename Sparsity>
+ void create_coupling_sparsity_pattern(const DoFHandler<dim0, spacedim> &space_dh,
+ const DoFHandler<dim1, spacedim> &immersed_dh,
+ const Quadrature<dim1> &quad,
+ Sparsity &sparsity,
+ 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.
+ *
+ * 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<int dim0, int dim1, int spacedim, typename Matrix>
+ void create_coupling_mass_matrix(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<dim0, spacedim> &space_mapping = StaticMappingQ1<dim0,spacedim>::mapping,
+ const Mapping<dim1, spacedim> &immersed_mapping = StaticMappingQ1<dim1, spacedim>::mapping);
+
+
+
+
+// === inline and template functions ===
+
+
+
+ template<int dim0, int dim1, int spacedim, typename Sparsity>
+ void create_coupling_sparsity_pattern(const DoFHandler<dim0, spacedim> &space_dh,
+ const DoFHandler<dim1, spacedim> &immersed_dh,
+ const Quadrature<dim1> &quad,
+ Sparsity &sparsity,
+ const ComponentMask &space_comps,
+ const ComponentMask &immersed_comps,
+ const Mapping<dim0, spacedim> &space_mapping,
+ const Mapping<dim1, spacedim> &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<dim1,spacedim>::active_cell_iterator
+ cell = immersed_dh.begin_active(),
+ endc = immersed_dh.end();
+
+ // Dof indices
+ std::vector<types::global_dof_index> dofs(immersed_fe.dofs_per_cell);
+ std::vector<types::global_dof_index> odofs(space_fe.dofs_per_cell);
+
+ 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
+ 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<unsigned int> space_gtl(space_fe.n_components(), numbers::invalid_unsigned_int);
+ std::vector<unsigned int> immersed_gtl(immersed_fe.n_components(), numbers::invalid_unsigned_int);
+
+ for (unsigned int i=0, j=0; i<space_gtl.size(); ++i)
+ if (space_c[i])
+ space_gtl[i] = j++;
+
+ for (unsigned int i=0, j=0; i<immersed_gtl.size(); ++i)
+ if (immersed_c[i])
+ immersed_gtl[i] = j++;
+
+ for (; cell != endc; ++cell)
+ // if(cell->subdomain_id() == this_mpi_process)
+ {
+ // Reinitialize the cell and the fe_values
+ fe_v.reinit(cell);
+ cell->get_dof_indices(dofs);
+
+ const std::vector<Point<spacedim> > &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<cells.size(); ++c)
+ {
+ // Get the ones in the current outer cell
+ typename DoFHandler<dim0,spacedim>::cell_iterator
+ ocell(*cells[c], &space_dh);
+ 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]);
+ }
+ }
+ }
+ }
+ }
+ }
+
+
+
+ template<int dim0, int dim1, int spacedim, typename Matrix>
+ void create_coupling_mass_matrix(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<dim0, spacedim> &space_mapping,
+ const Mapping<dim1, spacedim> &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<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
+ 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<unsigned int> space_gtl(space_fe.n_components(), numbers::invalid_unsigned_int);
+ std::vector<unsigned int> immersed_gtl(immersed_fe.n_components(), numbers::invalid_unsigned_int);
+
+ for (unsigned int i=0, j=0; i<space_gtl.size(); ++i)
+ if (space_c[i])
+ space_gtl[i] = j++;
+
+ for (unsigned int i=0, j=0; i<immersed_gtl.size(); ++i)
+ if (immersed_c[i])
+ immersed_gtl[i] = j++;
+
+ FullMatrix<double> cell_matrix(space_dh.get_fe().dofs_per_cell,
+ immersed_dh.get_fe().dofs_per_cell);
+
+ FEValues<dim1,spacedim> 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<dim1,spacedim>::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<Point<spacedim> > &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<cells.size(); ++c)
+ {
+ // Get the ones in the current outer cell
+ typename DoFHandler<dim0,spacedim>::active_cell_iterator
+ ocell(*cells[c], &space_dh);
+ const std::vector< Point<spacedim> > &qps = qpoints[c];
+ const std::vector< unsigned int > &ids = maps[c];
+
+ FEValues<dim0,spacedim> 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<space_dh.get_fe().dofs_per_cell; ++i)
+ {
+ const auto comp_i = space_dh.get_fe().system_to_component_index(i).first;
+ if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
+ for (unsigned int j=0; j<immersed_dh.get_fe().dofs_per_cell; ++j)
+ {
+ const auto comp_j = immersed_dh.get_fe().system_to_component_index(j).first;
+ if (space_gtl[comp_i] == immersed_gtl[comp_j])
+ for (unsigned int oq=0; oq<o_fe_v.n_quadrature_points; ++oq)
+ {
+ // Get the corrisponding q point
+ const unsigned int q=ids[oq];
+
+ cell_matrix(i,j) += ( fe_v.shape_value(j,q) *
+ o_fe_v.shape_value(i,oq) *
+ fe_v.JxW(q) );
+ }
+ }
+ }
+
+ // Now assemble the matrices
+ constraints.distribute_local_to_global (cell_matrix, odofs, dofs, matrix);
+ }
+ }
+ }
+}
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
--- /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/non_matching/coupling.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/fe_q.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 "../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<int dim, int spacedim>
+void test()
+{
+ deallog << "dim: " << dim << ", spacedim: "
+ << spacedim << std::endl;
+
+ Triangulation<dim,spacedim> tria;
+ Triangulation<spacedim,spacedim> 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<dim,spacedim> fe(1);
+ FE_Q<spacedim,spacedim> space_fe(1);
+
+ 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);
+
+ deallog << "Dofs : " << dh.n_dofs() << std::endl
+ << "Space dofs: " << space_dh.n_dofs() << std::endl;
+
+ QGauss<dim> 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<double> 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<double> 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<double> space_ones(space_dh.n_dofs());
+ Vector<double> ones(dh.n_dofs());
+
+ space_ones = 1.0;
+ coupling.Tvmult(ones, space_ones);
+ mass_matrix_inv.solve(ones);
+
+ Vector<double> 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>();
+}
--- /dev/null
+
+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
--- /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/non_matching/coupling.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.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 "../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<int dim, int spacedim>
+void test()
+{
+ deallog << "dim: " << dim << ", spacedim: "
+ << spacedim << std::endl;
+
+ Triangulation<dim,spacedim> tria;
+ Triangulation<spacedim,spacedim> 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<dim,spacedim> fe(1);
+ FESystem<spacedim,spacedim> space_fe(FE_Q<spacedim,spacedim>(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<dim,spacedim> dh(tria);
+ DoFHandler<spacedim,spacedim> 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<dim> 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<double> 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<double> 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<double> space_ones(space_dh.n_dofs());
+ Vector<double> ones(dh.n_dofs());
+
+ space_ones = 1.0;
+ coupling.Tvmult(ones, space_ones);
+ mass_matrix_inv.solve(ones);
+
+ Vector<double> 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>();
+}
--- /dev/null
+
+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
--- /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/non_matching/coupling.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.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 "../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<int dim, int spacedim>
+void test()
+{
+ deallog << "dim: " << dim << ", spacedim: "
+ << spacedim << std::endl;
+
+ Triangulation<dim,spacedim> tria;
+ Triangulation<spacedim,spacedim> 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<dim,spacedim> fe(FE_Q<dim,spacedim>(1), spacedim+1);
+ FESystem<spacedim,spacedim> space_fe(FE_Q<spacedim,spacedim>(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<dim,spacedim> dh(tria);
+ DoFHandler<spacedim,spacedim> 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<dim> 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<double> 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<double> 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<double> space_ones(space_dh.n_dofs());
+ Vector<double> 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<double> real_ones(dh.n_dofs());
+ for (unsigned int i=0; i<dh.n_dofs()/(spacedim+1); ++i)
+ real_ones((spacedim+1)*i) = 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>();
+}
--- /dev/null
+
+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
--- /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/non_matching/coupling.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.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 "../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;
+
+ Triangulation<dim,spacedim> tria;
+ Triangulation<spacedim,spacedim> 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<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);
+
+ deallog << "Dofs : " << dh.n_dofs() << std::endl
+ << "Space dofs: " << space_dh.n_dofs() << std::endl;
+
+ QGauss<dim> 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<double> 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<double> 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<double> space_square(space_dh.n_dofs());
+ Vector<double> squares(dh.n_dofs());
+ Vector<double> projected_squares(dh.n_dofs());
+
+ VectorTools::interpolate(space_dh, Functions::SquareFunction<spacedim>(),
+ space_square);
+ VectorTools::interpolate(dh, Functions::SquareFunction<spacedim>(),
+ 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>();
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff