--- /dev/null
+New: Two new methods NonMatching::compute_coupling_sparsity_pattern() and
+NonMatching::compute_coupling_mass_matrix() allow to construct the coupling between arbitrary grids,
+using a convolution kernel.
+<br>
+(Luca Heltai, 2019/04/19)
#include <deal.II/base/config.h>
+#include <deal.II/base/function_lib.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/component_mask.h>
#include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/grid/grid_tools_cache.h>
/**
* A namespace for functions offering tools to handle two meshes with no
- * alignment requirements, but where one of the meshes is embedded
- * inside the other in the real-space.
+ * alignment requirements.
*
* Typically these functions allow for computations on the real-space
* intersection between the two meshes e.g. surface integrals and
- * construction of mass matrices.
+ * construction of coupling matrices.
*/
namespace NonMatching
{
const ComponentMask & immersed_comps = ComponentMask(),
const Mapping<dim1, spacedim> &immersed_mapping =
StaticMappingQ1<dim1, spacedim>::mapping);
+
+ /**
+ * Create a coupling sparsity pattern for non-matching independent grids,
+ * using a convolution kernel with compact support of radius epsilon.
+ *
+ * Given two non-matching triangulations, representing the domains $\Omega^0$
+ * and $\Omega^1$, both embedded in $R^d$, and two finite element spaces
+ * $V^0(\Omega^0) = \text{span}\{v_i\}_{i=0}^n$ and $V^1(\Omega^1) =
+ * \text{span}\{w_\alpha\}_{\alpha=0}^m$, compute the sparsity pattern that
+ * would be necessary to assemble the matrix
+ *
+ * \f[
+ * M_{i\alpha} \dealcoloneq \int_{\Omega^0} \int_{\Omega^1}
+ * v_i(x) K^{\epsilon}(x-y) w_\alpha(y) dx \ dy,
+ * \quad i \in [0,n), \alpha \in [0,m),
+ * \f]
+ *
+ * where $V^0(\Omega^0)$ is the finite element space associated with the
+ * `dh0` passed to this function (or part of it, if specified in
+ * `comps0`), while $V^1(\Omega^1)$ is the finite element space associated
+ * with the `dh1` passed to this function (or part of it, if specified
+ * in `comps1`), and $K^\epsilon$ is a function with compact support included
+ * in a ball of radius $\epsilon$, derived from CutOffFunctionBase.
+ *
+ * The `comps0` and `comps1` masks are assumed to be ordered in
+ * the same way: the first component of `comps0` will couple with the
+ * first component of `comps1`, the second with the second, and so
+ * on. If one of the two masks has more non-zero than the other, then the
+ * excess components will be ignored.
+ *
+ * For both spaces, it is possible to specify a custom Mapping, which
+ * defaults to StaticMappingQ1 for both.
+ *
+ * This function will also work in parallel, provided that at least one of the
+ * triangulations is of type parallel::shared::Triangulation<dim1,spacedim>.
+ * An exception is thrown if both triagnulations are of type
+ * parallel::distributed::Triangulation<dim1,spacedim>.
+ *
+ * This function assumes that the convolution has support contained in a box
+ * of radius @p epsilon. If epsilon is set to zero, then we assume that the
+ * kernel is the Dirac delta distribution, and the call is forwarded to the
+ * method in this namespace with the same name, that does not take an epsilon
+ * as input. In this case, more restrictive conditions are required on the
+ * two spaces. See the documentation of the other
+ * create_coupling_sparsity_pattern() function.
+ *
+ * @author Luca Heltai, 2019.
+ */
+ template <int dim0,
+ int dim1,
+ int spacedim,
+ typename Sparsity,
+ typename Number = double>
+ void
+ create_coupling_sparsity_pattern(
+ const double & epsilon,
+ const GridTools::Cache<dim0, spacedim> &cache0,
+ const GridTools::Cache<dim1, spacedim> &cache1,
+ const DoFHandler<dim0, spacedim> & dh0,
+ const DoFHandler<dim1, spacedim> & dh1,
+ Sparsity & sparsity,
+ const AffineConstraints<Number> &constraints0 = AffineConstraints<Number>(),
+ const ComponentMask & comps0 = ComponentMask(),
+ const ComponentMask & comps1 = ComponentMask());
+
+ /**
+ * Create a coupling mass matrix for non-matching independent grids,
+ * using a convolution kernel with compact support.
+ *
+ * Given two non-matching triangulations, representing the domains
+ * $\Omega^0$ and $\Omega^1$, both embedded in $R^d$, and two finite element
+ * spaces $V^0(\Omega^0) = \text{span}\{v_i\}_{i=0}^n$ and $V^1(\Omega^1) =
+ * \text{span}\{w_\alpha\}_{\alpha=0}^m$, compute the matrix
+ *
+ * \f[
+ * M_{i\alpha} \dealcoloneq \int_{\Omega^0} \int_{\Omega^1}
+ * v_i(x) K^{\epsilon}(x-y) w_\alpha(y) dx \ dy,
+ * \quad i \in [0,n), \alpha \in [0,m),
+ * \f]
+ *
+ * where $V^0(\Omega^0)$ is the finite element space associated with the
+ * `dh0` passed to this function (or part of it, if specified in
+ * `comps0`), while $V^1(\Omega^1)$ is the finite element space associated
+ * with the `dh1` passed to this function (or part of it, if specified
+ * in `comps1`), and $K^\epsilon$ is a function with compact support included,
+ * in a ball of radius $\epsilon$, derived from CutOffFunctionBase.
+ *
+ * The corresponding sparsity patterns can be computed by calling the
+ * make_coupling_sparsity_pattern() function.
+ *
+ * The `comps0` and `comps1` masks are assumed to be ordered in
+ * the same way: the first component of `comps0` will couple with the
+ * first component of `comps1`, the second with the second, and so
+ * on. If one of the two masks has more non-zero than the other, then the
+ * excess components will be ignored.
+ *
+ * For both spaces, it is possible to specify a custom Mapping, which
+ * defaults to StaticMappingQ1 for both.
+ *
+ * This function will also work in parallel, provided that one of the two
+ * triangulations is of type parallel::shared::Triangulation<dim1,spacedim>.
+ * An exception is thrown if both triangulations are of type
+ * parallel::distributed::Triangulation<dim1,spacedim>.
+ *
+ * The parameter @p epsilon is used to set the size of the cut-off function
+ * used to compute the convolution. If epsilon is set to zero, then we assume
+ * that the kernel is the Dirac delta distribution, and the call is forwarded
+ * to the method in this namespace with the same name, that does not take an
+ * epsilon as input.
+ *
+ * @author Luca Heltai, 2019.
+ */
+ template <int dim0, int dim1, int spacedim, typename Matrix>
+ void
+ create_coupling_mass_matrix(
+ Functions::CutOffFunctionBase<spacedim> & kernel,
+ const double & epsilon,
+ const GridTools::Cache<dim0, spacedim> & cache0,
+ const GridTools::Cache<dim1, spacedim> & cache1,
+ const DoFHandler<dim0, spacedim> & dh0,
+ const DoFHandler<dim1, spacedim> & dh1,
+ const Quadrature<dim0> & quadrature0,
+ const Quadrature<dim1> & quadrature1,
+ Matrix & matrix,
+ const AffineConstraints<typename Matrix::value_type> &constraints0 =
+ AffineConstraints<typename Matrix::value_type>(),
+ const ComponentMask &comps0 = ComponentMask(),
+ const ComponentMask &comps1 = ComponentMask());
} // namespace NonMatching
DEAL_II_NAMESPACE_CLOSE
}
return {std::move(points_over_local_cells), std::move(used_cells_ids)};
}
+
+
+ /**
+ * Given two ComponentMasks and the corresponding finite element spaces,
+ * compute a pairing between the selected components of the first finite
+ * element space, and the selected components of the second finite element
+ * space.
+ */
+ template <int dim0, int dim1, int spacedim>
+ std::pair<std::vector<unsigned int>, std::vector<unsigned int>>
+ compute_components_coupling(const ComponentMask & comps0,
+ const ComponentMask & comps1,
+ const FiniteElement<dim0, spacedim> &fe0,
+ const FiniteElement<dim1, spacedim> &fe1)
+ {
+ // Take care of components
+ const ComponentMask mask0 =
+ (comps0.size() == 0 ? ComponentMask(fe0.n_components(), true) : comps0);
+
+ const ComponentMask mask1 =
+ (comps1.size() == 0 ? ComponentMask(fe1.n_components(), true) : comps1);
+
+ AssertDimension(mask0.size(), fe0.n_components());
+ AssertDimension(mask1.size(), fe1.n_components());
+
+ // Global to local indices
+ std::vector<unsigned int> gtl0(fe0.n_components(),
+ numbers::invalid_unsigned_int);
+ std::vector<unsigned int> gtl1(fe1.n_components(),
+ numbers::invalid_unsigned_int);
+
+ for (unsigned int i = 0, j = 0; i < gtl0.size(); ++i)
+ if (mask0[i])
+ gtl0[i] = j++;
+
+ for (unsigned int i = 0, j = 0; i < gtl1.size(); ++i)
+ if (mask1[i])
+ gtl1[i] = j++;
+ return {gtl0, gtl1};
+ }
} // namespace internal
template <int dim0,
{
AssertDimension(sparsity.n_rows(), space_dh.n_dofs());
AssertDimension(sparsity.n_cols(), immersed_dh.n_dofs());
- static_assert(dim1 <= dim0, "This function can only work if dim1 <= dim0");
+ Assert(dim1 <= dim0,
+ ExcMessage("This function can only work if dim1 <= dim0"));
Assert((dynamic_cast<
const parallel::distributed::Triangulation<dim1, spacedim> *>(
&immersed_dh.get_triangulation()) == nullptr),
{
AssertDimension(matrix.m(), space_dh.n_dofs());
AssertDimension(matrix.n(), immersed_dh.n_dofs());
- static_assert(dim1 <= dim0, "This function can only work if dim1 <= dim0");
+ Assert(dim1 <= dim0,
+ ExcMessage("This function can only work if dim1 <= dim0"));
Assert((dynamic_cast<
const parallel::distributed::Triangulation<dim1, spacedim> *>(
&immersed_dh.get_triangulation()) == nullptr),
}
}
+ template <int dim0,
+ int dim1,
+ int spacedim,
+ typename Sparsity,
+ typename Number>
+ void
+ create_coupling_sparsity_pattern(
+ const double & epsilon,
+ const GridTools::Cache<dim0, spacedim> &cache0,
+ const GridTools::Cache<dim1, spacedim> &cache1,
+ const DoFHandler<dim0, spacedim> & dh0,
+ const DoFHandler<dim1, spacedim> & dh1,
+ Sparsity & sparsity,
+ const AffineConstraints<Number> & constraints0,
+ const ComponentMask & comps0,
+ const ComponentMask & comps1)
+ {
+ if (epsilon == 0.0)
+ {
+ Assert(dim1 <= dim0,
+ ExcMessage("When epsilon is zero, you can only "
+ "call this function with dim1 <= dim0."));
+ create_coupling_sparsity_pattern(cache0,
+ dh0,
+ dh1,
+ QGauss<dim1>(dh1.get_fe().degree + 1),
+ sparsity,
+ constraints0,
+ comps0,
+ comps1,
+ cache1.get_mapping());
+ return;
+ }
+ AssertDimension(sparsity.n_rows(), dh0.n_dofs());
+ AssertDimension(sparsity.n_cols(), dh1.n_dofs());
+
+ const bool zero_is_distributed =
+ (dynamic_cast<const parallel::distributed::Triangulation<dim1, spacedim>
+ *>(&dh0.get_triangulation()) != nullptr);
+ const bool one_is_distributed =
+ (dynamic_cast<const parallel::distributed::Triangulation<dim1, spacedim>
+ *>(&dh1.get_triangulation()) != nullptr);
+
+ // We bail out if both are distributed triangulations
+ Assert(!zero_is_distributed || !one_is_distributed, ExcNotImplemented());
+
+ // If we can loop on both, we decide where to make the outer loop according
+ // to the size of triangulation. The reasoning is the following:
+ // - Access to the tree: log(N)
+ // - We compute intersection for each of the outer loop cells (M)
+ // Total cost (besides the setup) is: M log(N)
+ // If we can, make sure M is the smallest number
+ const bool outer_loop_on_zero =
+ (zero_is_distributed && !one_is_distributed) ||
+ (dh1.get_triangulation().n_active_cells() >
+ dh0.get_triangulation().n_active_cells());
+
+ const auto &fe0 = dh0.get_fe();
+ const auto &fe1 = dh1.get_fe();
+
+ // Dof indices
+ std::vector<types::global_dof_index> dofs0(fe0.dofs_per_cell);
+ std::vector<types::global_dof_index> dofs1(fe1.dofs_per_cell);
+
+ if (outer_loop_on_zero)
+ {
+ std::cout << "Looping on zero." << std::endl;
+
+ Assert(one_is_distributed == false, ExcInternalError());
+
+ const auto &tree1 = cache1.get_cell_bounding_boxes_rtree();
+
+ std::vector<std::pair<
+ BoundingBox<spacedim>,
+ typename Triangulation<dim1, spacedim>::active_cell_iterator>>
+ intersection;
+
+ for (const auto &cell0 : dh0.active_cell_iterators())
+ if (cell0->is_locally_owned())
+ {
+ intersection.resize(0);
+
+ BoundingBox<spacedim> box0;
+ if (cache0.get_mapping().preserves_vertex_locations())
+ box0 = cell0->bounding_box();
+ else
+ {
+ const auto vertices =
+ cache0.get_mapping().get_vertices(cell0);
+ Point<spacedim> p0 = vertices[0], p1 = vertices[0];
+ for (unsigned int j = 1; j < vertices.size(); ++j)
+ for (unsigned int d = 0; d < spacedim; ++d)
+ {
+ p0[d] = std::min(p0[d], vertices[j][d]);
+ p1[d] = std::max(p1[d], vertices[j][d]);
+ }
+ box0 = BoundingBox<spacedim>({p0, p1});
+ }
+
+ box0.extend(epsilon);
+ boost::geometry::index::query(tree1,
+ boost::geometry::index::intersects(
+ box0),
+ std::back_inserter(intersection));
+ if (!intersection.empty())
+ {
+ cell0->get_dof_indices(dofs0);
+ for (const auto &entry : intersection)
+ {
+ typename DoFHandler<dim1, spacedim>::cell_iterator cell1(
+ *entry.second, &dh1);
+ cell1->get_dof_indices(dofs1);
+ constraints0.add_entries_local_to_global(dofs0,
+ dofs1,
+ sparsity);
+ }
+ }
+ }
+ }
+ else
+ {
+ Assert(zero_is_distributed == false, ExcInternalError());
+ const auto &tree0 = cache0.get_cell_bounding_boxes_rtree();
+
+ std::vector<std::pair<
+ BoundingBox<spacedim>,
+ typename Triangulation<dim0, spacedim>::active_cell_iterator>>
+ intersection;
+
+ for (const auto &cell1 : dh1.active_cell_iterators())
+ if (cell1->is_locally_owned())
+ {
+ intersection.resize(0);
+
+ BoundingBox<spacedim> box1;
+ Assert(&cache1.get_mapping() != nullptr, ExcInternalError());
+ if (cache1.get_mapping().preserves_vertex_locations())
+ box1 = cell1->bounding_box();
+ else
+ {
+ const auto vertices =
+ cache1.get_mapping().get_vertices(cell1);
+ Point<spacedim> p0 = vertices[0], p1 = vertices[0];
+ for (unsigned int j = 1; j < vertices.size(); ++j)
+ for (unsigned int d = 0; d < spacedim; ++d)
+ {
+ p0[d] = std::min(p0[d], vertices[j][d]);
+ p1[d] = std::max(p1[d], vertices[j][d]);
+ }
+ box1 = BoundingBox<spacedim>({p0, p1});
+ }
+
+ box1.extend(epsilon);
+ boost::geometry::index::query(tree0,
+ boost::geometry::index::intersects(
+ box1),
+ std::back_inserter(intersection));
+ if (!intersection.empty())
+ {
+ cell1->get_dof_indices(dofs1);
+ for (const auto &entry : intersection)
+ {
+ typename DoFHandler<dim0, spacedim>::cell_iterator cell0(
+ *entry.second, &dh0);
+ cell0->get_dof_indices(dofs0);
+ constraints0.add_entries_local_to_global(dofs0,
+ dofs1,
+ sparsity);
+ }
+ }
+ }
+ }
+ }
+
+
+
+ template <int dim0, int dim1, int spacedim, typename Matrix>
+ void
+ create_coupling_mass_matrix(
+ Functions::CutOffFunctionBase<spacedim> & kernel,
+ const double & epsilon,
+ const GridTools::Cache<dim0, spacedim> & cache0,
+ const GridTools::Cache<dim1, spacedim> & cache1,
+ const DoFHandler<dim0, spacedim> & dh0,
+ const DoFHandler<dim1, spacedim> & dh1,
+ const Quadrature<dim0> & quadrature0,
+ const Quadrature<dim1> & quadrature1,
+ Matrix & matrix,
+ const AffineConstraints<typename Matrix::value_type> &constraints0,
+ const ComponentMask & comps0,
+ const ComponentMask & comps1)
+ {
+ if (epsilon == 0)
+ {
+ Assert(dim1 <= dim0,
+ ExcMessage("When epsilon is zero, you can only "
+ "call this function with dim1 <= dim0."));
+ create_coupling_mass_matrix(cache0,
+ dh0,
+ dh1,
+ quadrature1,
+ matrix,
+ constraints0,
+ comps0,
+ comps1,
+ cache1.get_mapping());
+ return;
+ }
+
+ AssertDimension(matrix.m(), dh0.n_dofs());
+ AssertDimension(matrix.n(), dh1.n_dofs());
+
+ const bool zero_is_distributed =
+ (dynamic_cast<const parallel::distributed::Triangulation<dim1, spacedim>
+ *>(&dh0.get_triangulation()) != nullptr);
+ const bool one_is_distributed =
+ (dynamic_cast<const parallel::distributed::Triangulation<dim1, spacedim>
+ *>(&dh1.get_triangulation()) != nullptr);
+
+ // We bail out if both are distributed triangulations
+ Assert(!zero_is_distributed || !one_is_distributed, ExcNotImplemented());
+
+ // If we can loop on both, we decide where to make the outer loop according
+ // to the size of triangulation. The reasoning is the following:
+ // - Access to the tree: log(N)
+ // - We compute intersection for each of the outer loop cells (M)
+ // Total cost (besides the setup) is: M log(N)
+ // If we can, make sure M is the smallest number
+ const bool outer_loop_on_zero =
+ (zero_is_distributed && !one_is_distributed) ||
+ (dh1.get_triangulation().n_active_cells() >
+ dh0.get_triangulation().n_active_cells());
+
+ const auto &fe0 = dh0.get_fe();
+ const auto &fe1 = dh1.get_fe();
+
+ FEValues<dim0, spacedim> fev0(cache0.get_mapping(),
+ fe0,
+ quadrature0,
+ update_values | update_JxW_values |
+ update_quadrature_points);
+
+ FEValues<dim1, spacedim> fev1(cache1.get_mapping(),
+ fe1,
+ quadrature1,
+ update_values | update_JxW_values |
+ update_quadrature_points);
+
+ // Dof indices
+ std::vector<types::global_dof_index> dofs0(fe0.dofs_per_cell);
+ std::vector<types::global_dof_index> dofs1(fe1.dofs_per_cell);
+
+ // Local Matrix
+ FullMatrix<typename Matrix::value_type> cell_matrix(fe0.dofs_per_cell,
+ fe1.dofs_per_cell);
+
+ // Global to local indices
+ auto p = internal::compute_components_coupling(comps0, comps1, fe0, fe1);
+ const auto >l0 = p.first;
+ const auto >l1 = p.second;
+
+ kernel.set_radius(epsilon);
+ std::vector<double> kernel_values(quadrature1.size());
+
+ auto assemble_one_pair = [&]() {
+ cell_matrix = 0;
+ for (unsigned int q0 = 0; q0 < quadrature0.size(); ++q0)
+ {
+ kernel.set_center(fev0.quadrature_point(q0));
+ kernel.value_list(fev1.get_quadrature_points(), kernel_values);
+ for (unsigned int q1 = 0; q1 < quadrature1.size(); ++q1)
+ if (kernel_values[q1] != 0.0)
+ {
+ for (unsigned int i = 0; i < fe0.dofs_per_cell; ++i)
+ {
+ const auto comp_i = fe0.system_to_component_index(i).first;
+ if (gtl0[comp_i] != numbers::invalid_unsigned_int)
+ for (unsigned int j = 0; j < fe1.dofs_per_cell; ++j)
+ {
+ const auto comp_j =
+ fe1.system_to_component_index(j).first;
+ if (gtl1[comp_j] == gtl0[comp_i])
+ {
+ cell_matrix(i, j) += fev0.shape_value(i, q0) *
+ fev1.shape_value(j, q1) *
+ kernel_values[q1] *
+ fev0.JxW(q0) * fev1.JxW(q1);
+ }
+ }
+ }
+ }
+ }
+
+ constraints0.distribute_local_to_global(cell_matrix,
+ dofs0,
+ dofs1,
+ matrix);
+ };
+
+ if (outer_loop_on_zero)
+ {
+ Assert(one_is_distributed == false, ExcInternalError());
+
+ const auto &tree1 = cache1.get_cell_bounding_boxes_rtree();
+
+ std::vector<std::pair<
+ BoundingBox<spacedim>,
+ typename Triangulation<dim1, spacedim>::active_cell_iterator>>
+ intersection;
+
+ for (const auto &cell0 : dh0.active_cell_iterators())
+ if (cell0->is_locally_owned())
+ {
+ intersection.resize(0);
+
+ BoundingBox<spacedim> box0;
+ if (cache0.get_mapping().preserves_vertex_locations())
+ box0 = cell0->bounding_box();
+ else
+ {
+ const auto vertices =
+ cache0.get_mapping().get_vertices(cell0);
+ Point<spacedim> p0 = vertices[0], p1 = vertices[0];
+ for (unsigned int j = 1; j < vertices.size(); ++j)
+ for (unsigned int d = 0; d < spacedim; ++d)
+ {
+ p0[d] = std::min(p0[d], vertices[j][d]);
+ p1[d] = std::max(p1[d], vertices[j][d]);
+ }
+ box0 = BoundingBox<spacedim>({p0, p1});
+ }
+
+ box0.extend(epsilon);
+ boost::geometry::index::query(tree1,
+ boost::geometry::index::intersects(
+ box0),
+ std::back_inserter(intersection));
+ if (!intersection.empty())
+ {
+ cell0->get_dof_indices(dofs0);
+ fev0.reinit(cell0);
+ for (const auto &entry : intersection)
+ {
+ typename DoFHandler<dim1, spacedim>::cell_iterator cell1(
+ *entry.second, &dh1);
+ cell1->get_dof_indices(dofs1);
+ fev1.reinit(cell1);
+ assemble_one_pair();
+ }
+ }
+ }
+ }
+ else
+ {
+ Assert(zero_is_distributed == false, ExcInternalError());
+ const auto &tree0 = cache0.get_cell_bounding_boxes_rtree();
+
+ std::vector<std::pair<
+ BoundingBox<spacedim>,
+ typename Triangulation<dim0, spacedim>::active_cell_iterator>>
+ intersection;
+
+ for (const auto &cell1 : dh1.active_cell_iterators())
+ if (cell1->is_locally_owned())
+ {
+ intersection.resize(0);
+ BoundingBox<spacedim> box1;
+ if (cache1.get_mapping().preserves_vertex_locations())
+ box1 = cell1->bounding_box();
+ else
+ {
+ const auto vertices =
+ cache1.get_mapping().get_vertices(cell1);
+ Point<spacedim> p0 = vertices[0], p1 = vertices[0];
+ for (unsigned int j = 1; j < vertices.size(); ++j)
+ for (unsigned int d = 0; d < spacedim; ++d)
+ {
+ p0[d] = std::min(p0[d], vertices[j][d]);
+ p1[d] = std::max(p1[d], vertices[j][d]);
+ }
+ box1 = BoundingBox<spacedim>({p0, p1});
+ }
+ box1.extend(epsilon);
+ boost::geometry::index::query(tree0,
+ boost::geometry::index::intersects(
+ box1),
+ std::back_inserter(intersection));
+ if (!intersection.empty())
+ {
+ cell1->get_dof_indices(dofs1);
+ fev1.reinit(cell1);
+ for (const auto &entry : intersection)
+ {
+ typename DoFHandler<dim0, spacedim>::cell_iterator cell0(
+ *entry.second, &dh0);
+ cell0->get_dof_indices(dofs0);
+ fev0.reinit(cell0);
+ assemble_one_pair();
+ }
+ }
+ }
+ }
+ }
+
#include "coupling.inst"
} // namespace NonMatching
-
DEAL_II_NAMESPACE_CLOSE
Sparsity : SPARSITY_PATTERNS;
S : REAL_AND_COMPLEX_SCALARS)
{
-#if dim1 <= dim0 && dim0 <= spacedim
+#if dim1 <= spacedim && dim0 <= spacedim
template void create_coupling_sparsity_pattern(
const DoFHandler<dim0, spacedim> &space_dh,
const DoFHandler<dim1, spacedim> &immersed_dh,
const ComponentMask & space_comps,
const ComponentMask & immersed_comps,
const Mapping<dim1, spacedim> & immersed_mapping);
+
+ template void create_coupling_sparsity_pattern(
+ const double & epsilon,
+ const GridTools::Cache<dim0, spacedim> &cache0,
+ const GridTools::Cache<dim1, spacedim> &cache1,
+ const DoFHandler<dim0, spacedim> & dh0,
+ const DoFHandler<dim1, spacedim> & dh1,
+ Sparsity & sparsity,
+ const AffineConstraints<S> & constraints0,
+ const ComponentMask & comps0,
+ const ComponentMask & comps1);
#endif
}
for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS;
Matrix : SPARSE_MATRICES)
{
-#if dim1 <= dim0 && dim0 <= spacedim
+#if dim1 <= spacedim && dim0 <= spacedim
template void create_coupling_mass_matrix(
const DoFHandler<dim0, spacedim> & space_dh,
const DoFHandler<dim1, spacedim> & immersed_dh,
const ComponentMask & space_comps,
const ComponentMask & immersed_comps,
const Mapping<dim1, spacedim> & immersed_mapping);
+
+ template void create_coupling_mass_matrix(
+ Functions::CutOffFunctionBase<spacedim> & kernel,
+ const double & epsilon,
+ const GridTools::Cache<dim0, spacedim> & cache0,
+ const GridTools::Cache<dim1, spacedim> & cache1,
+ const DoFHandler<dim0, spacedim> & dh0,
+ const DoFHandler<dim1, spacedim> & dh1,
+ const Quadrature<dim0> & quadrature0,
+ const Quadrature<dim1> & quadrature1,
+ Matrix & matrix,
+ const AffineConstraints<typename Matrix::value_type> &constraints0,
+ const ComponentMask & comps0,
+ const ComponentMask & comps1);
#endif
- }
+ }
\ No newline at end of file
--- /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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/function_parser.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
+
+#include <deal.II/non_matching/coupling.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 dim0, int dim1, int spacedim>
+void
+test()
+{
+ deallog << "dim0: " << dim0 << ", dim1: " << dim1
+ << ", spacedim: " << spacedim << std::endl;
+
+ Triangulation<dim0, spacedim> tria0;
+ Triangulation<dim1, spacedim> tria1;
+
+ GridGenerator::hyper_cube(tria0, -.4, .3);
+ GridGenerator::hyper_cube(tria1, -1, 1);
+
+ tria0.refine_global(3);
+ tria1.refine_global(2);
+
+ FE_Q<dim0, spacedim> fe0(1);
+ FE_Q<dim1, spacedim> fe1(1);
+
+ deallog << "FE0 : " << fe0.get_name() << std::endl
+ << "FE1 : " << fe1.get_name() << std::endl;
+
+ DoFHandler<dim0, spacedim> dh0(tria0);
+ DoFHandler<dim1, spacedim> dh1(tria1);
+
+ GridTools::Cache<dim0, spacedim> cache0(tria0);
+ GridTools::Cache<dim1, spacedim> cache1(tria1);
+
+ dh0.distribute_dofs(fe0);
+ dh1.distribute_dofs(fe1);
+
+ AffineConstraints<double> constraints0;
+ DoFTools::make_hanging_node_constraints(dh0, constraints0);
+
+ constraints0.close();
+
+ deallog << "Dofs 0 : " << dh0.n_dofs() << std::endl
+ << "Dofs 1 : " << dh1.n_dofs() << std::endl
+ << "Constrained dofs: " << constraints0.n_constraints() << std::endl;
+
+ QGauss<dim0> quad0(3); // Quadrature for coupling
+ QGauss<dim1> quad1(3); // Quadrature for coupling
+
+ const double epsilon = 2 * std::max(GridTools::maximal_cell_diameter(tria0),
+ GridTools::maximal_cell_diameter(tria1));
+
+ deallog << "Epsilon: " << epsilon << std::endl;
+
+ SparsityPattern sparsity;
+ {
+ DynamicSparsityPattern dsp(dh0.n_dofs(), dh1.n_dofs());
+ NonMatching::create_coupling_sparsity_pattern(
+ epsilon, cache0, cache1, dh0, dh1, dsp, constraints0);
+ sparsity.copy_from(dsp);
+ }
+ SparseMatrix<double> coupling(sparsity);
+
+ Functions::CutOffFunctionC1<spacedim> dirac(
+ 1,
+ Point<spacedim>(),
+ 1,
+ Functions::CutOffFunctionBase<spacedim>::no_component,
+ true);
+
+ NonMatching::create_coupling_mass_matrix(dirac,
+ epsilon,
+ cache0,
+ cache1,
+ dh0,
+ dh1,
+ quad0,
+ quad1,
+ coupling,
+ constraints0);
+
+ SparsityPattern mass_sparsity1;
+ {
+ DynamicSparsityPattern dsp(dh1.n_dofs(), dh1.n_dofs());
+ DoFTools::make_sparsity_pattern(dh1, dsp);
+ mass_sparsity1.copy_from(dsp);
+ }
+ SparseMatrix<double> mass_matrix1(mass_sparsity1);
+ MatrixTools::create_mass_matrix(dh1, quad1, mass_matrix1);
+
+ SparseDirectUMFPACK mass_matrix1_inv;
+ mass_matrix1_inv.factorize(mass_matrix1);
+
+ // now take ones in dh0, project them onto dh1,
+ // get back ones, and check for the error.
+ //
+ // WARNINGS: Only works if dh1 is immersed in dh0
+
+ Vector<double> ones0(dh0.n_dofs());
+ Vector<double> ones1(dh1.n_dofs());
+
+ ones0 = 1.0;
+ coupling.Tvmult(ones1, ones0);
+ mass_matrix1_inv.solve(ones1);
+
+ Vector<double> exact_ones1(dh1.n_dofs());
+ exact_ones1 = 1.0;
+ ones1 -= exact_ones1;
+
+ deallog << "Error on constants: " << ones1.l2_norm() << std::endl;
+}
+
+
+
+int
+main()
+{
+ initlog(1);
+ test<1, 1, 1>();
+ test<2, 1, 2>();
+ test<2, 2, 2>();
+ test<3, 2, 3>();
+ test<3, 3, 3>();
+}
projects / dealii.git / commitdiff