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;
- }
- }
+ // [TODO]: when the add_entries_local_to_global below will implement
+ // the version with the dof_mask, this should be uncommented.
+ //
+ // // Construct a dof_mask, used to distribute entries to the sparsity
+ // able< 2, bool > dof_mask(space_fe.dofs_per_cell,
+ // immersed_fe.dofs_per_cell);
+ // of_mask.fill(false);
+ // or (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)
{
if (ocell->is_locally_owned())
{
ocell->get_dof_indices(odofs);
- constraints.add_entries_local_to_global(odofs, dofs, sparsity);
+ // [TODO]: When the following function will be implemented
+ // for the case of non-trivial dof_mask, we should
+ // uncomment the missing part.
+ constraints.add_entries_local_to_global(odofs, dofs, sparsity); //, true, dof_mask);
}
}
}
--- /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.
+//
+// Even when locally refined grids are used.
+
+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(1);
+ space_tria.begin_active()->set_refine_flag();
+ space_tria.execute_coarsening_and_refinement();
+
+ 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);
+
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints(space_dh, constraints);
+
+ constraints.close();
+
+ deallog << "Dofs : " << dh.n_dofs() << std::endl
+ << "Space dofs: " << space_dh.n_dofs() << std::endl
+ << "Constrained dofs: " << constraints.n_constraints() << 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, constraints);
+ sparsity.copy_from(dsp);
+ }
+ SparseMatrix<double> coupling(sparsity);
+ NonMatching::create_coupling_mass_matrix(space_dh, dh, quad,
+ coupling, constraints);
+
+ 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: 4
+DEAL::Constrained dofs: 0
+DEAL::Error on constants: 5.20741e-16
+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: 14
+DEAL::Constrained dofs: 2
+DEAL::Error on constants: 4.00297e-16
+DEAL::dim: 2, spacedim: 2
+DEAL::FE : FE_Q<2>(1)
+DEAL::Space FE: FE_Q<2>(1)
+DEAL::Dofs : 9
+DEAL::Space dofs: 14
+DEAL::Constrained dofs: 2
+DEAL::Error on constants: 1.49365e-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: 46
+DEAL::Constrained dofs: 12
+DEAL::Error on constants: 1.41744e-15
+DEAL::dim: 3, spacedim: 3
+DEAL::FE : FE_Q<3>(1)
+DEAL::Space FE: FE_Q<3>(1)
+DEAL::Dofs : 27
+DEAL::Space dofs: 46
+DEAL::Constrained dofs: 12
+DEAL::Error on constants: 1.84498e-14