From 9af2cd073789816f62f90d98731081c58572e311 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Tue, 6 Mar 2018 12:08:09 +0100 Subject: [PATCH] Split coupling.h into coupling.cc and coupling.inst --- cmake/config/template-arguments.in | 3 - include/deal.II/non_matching/coupling.h | 277 +++--------------------- source/non_matching/CMakeLists.txt | 2 + source/non_matching/coupling.cc | 249 +++++++++++++++++++++ source/non_matching/coupling.inst.in | 48 ++++ 5 files changed, 334 insertions(+), 245 deletions(-) create mode 100644 source/non_matching/coupling.cc create mode 100644 source/non_matching/coupling.inst.in diff --git a/cmake/config/template-arguments.in b/cmake/config/template-arguments.in index 5a239ba1b6..fc9466d78d 100644 --- a/cmake/config/template-arguments.in +++ b/cmake/config/template-arguments.in @@ -174,9 +174,6 @@ VECTORS_WITH_MATRIX := { Vector; SPARSE_MATRICES := { SparseMatrix; SparseMatrix; - SparseMatrixEZ; - SparseMatrixEZ; - BlockSparseMatrix; BlockSparseMatrix; diff --git a/include/deal.II/non_matching/coupling.h b/include/deal.II/non_matching/coupling.h index 3b8be3e477..9bdbd19e5a 100644 --- a/include/deal.II/non_matching/coupling.h +++ b/include/deal.II/non_matching/coupling.h @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 1998 - 2017 by the deal.II authors +// Copyright (C) 2018 by the deal.II authors // // This file is part of the deal.II library. // @@ -43,32 +43,34 @@ namespace NonMatching * * 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 + * $V(\Omega) = \text{span}\{v_i\}_{i=0}^n$ and $Q(B) = + * \text{span}\{w_j\}_{j=0}^m$, 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) + * M_{ij} := \int_{B} v_i(x) w_j(x) dx, \quad i \in [0,n), j \in [0,m), * \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. + * where $V(\Omega)$ is the finite element space associated with the + * `space_dh` passed to this function (or part of it, if specified in + * `space_comps`), while $Q(B)$ is the finite element space associated with + * the `immersed_dh` passed to this function (or part of it, if specified in + * `immersed_comps`). + * + * The `sparsity` is filled 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. + * on. If one of the two masks has more non-zero than 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. + * thrown by the algorithm that computes the quadrature point locations. In + * particular, notice that this function only makes sens 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 * defaults to StaticMappingQ1 for both. @@ -88,20 +90,21 @@ namespace NonMatching /** * 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 - * + * $V(\Omega) = \text{span}\{v_i\}_{i=0}^n$ and $Q(B) = + * \text{span}\{w_j\}_{j=0}^m$, compute the coupling 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) + * M_{ij} := \int_{B} v_i(x) w_j(x) dx, \quad i \in [0,n), j \in [0,m), * \f] + * where $V(\Omega)$ is the finite element space associated with the + * `space_dh` passed to this function (or part of it, if specified in + * `space_comps`), while $Q(B)$ is the finite element space associated with + * the `immersed_dh` passed to this function (or part of it, if specified in + * `immersed_comps`). * - * 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 + * The corresponding sparsity patterns 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$. @@ -109,11 +112,13 @@ namespace NonMatching * 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. + * on. If one of the two masks has more non-zero entries non-zero than 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. + * thrown by the algorithm that computes the quadrature point locations. In + * particular, notice that this function only makes sens 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 * defaults to StaticMappingQ1 for both. @@ -130,218 +135,6 @@ namespace NonMatching const ComponentMask &immersed_comps = ComponentMask(), const Mapping &space_mapping = StaticMappingQ1::mapping, const Mapping &immersed_mapping = StaticMappingQ1::mapping); - - - - -// === inline and template functions === - - - - template - void create_coupling_sparsity_pattern(const DoFHandler &space_dh, - const DoFHandler &immersed_dh, - const Quadrature &quad, - Sparsity &sparsity, - const ComponentMask &space_comps, - const ComponentMask &immersed_comps, - const Mapping &space_mapping, - const Mapping &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::active_cell_iterator - cell = immersed_dh.begin_active(), - endc = immersed_dh.end(); - - // Dof indices - std::vector dofs(immersed_fe.dofs_per_cell); - std::vector odofs(space_fe.dofs_per_cell); - - FEValues fe_v(immersed_mapping, immersed_fe, quad, - update_quadrature_points); - - GridTools::Cache 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 space_gtl(space_fe.n_components(), numbers::invalid_unsigned_int); - std::vector immersed_gtl(immersed_fe.n_components(), numbers::invalid_unsigned_int); - - for (unsigned int i=0, j=0; isubdomain_id() == this_mpi_process) - { - // Reinitialize the cell and the fe_values - fe_v.reinit(cell); - cell->get_dof_indices(dofs); - - const std::vector > &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::cell_iterator - ocell(*cells[c], &space_dh); - ocell->get_dof_indices(odofs); - for (unsigned int i=0; i - void create_coupling_mass_matrix(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 &space_mapping, - const Mapping &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 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 - 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 space_gtl(space_fe.n_components(), numbers::invalid_unsigned_int); - std::vector immersed_gtl(immersed_fe.n_components(), numbers::invalid_unsigned_int); - - for (unsigned int i=0, j=0; i cell_matrix(space_dh.get_fe().dofs_per_cell, - immersed_dh.get_fe().dofs_per_cell); - - FEValues 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::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 > &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::active_cell_iterator - ocell(*cells[c], &space_dh); - 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, - update_values); - o_fe_v.reinit(ocell); - ocell->get_dof_indices(odofs); - - // Reset the matrices. - cell_matrix = 0; - - for (unsigned int i=0; i + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN +namespace NonMatching +{ + template + void create_coupling_sparsity_pattern(const DoFHandler &space_dh, + const DoFHandler &immersed_dh, + const Quadrature &quad, + Sparsity &sparsity, + const ComponentMask &space_comps, + const ComponentMask &immersed_comps, + const Mapping &space_mapping, + const Mapping &immersed_mapping) + { + 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"); + + 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::active_cell_iterator + cell = immersed_dh.begin_active(), + endc = immersed_dh.end(); + + // Dof indices + std::vector dofs(immersed_fe.dofs_per_cell); + std::vector odofs(space_fe.dofs_per_cell); + + 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 ? + ComponentMask(space_fe.n_components(), true) : + space_comps); + + const ComponentMask immersed_c + = (immersed_comps.size() == 0 ? + ComponentMask(immersed_fe.n_components(), true) : + immersed_comps); + + AssertDimension(space_c.size(), space_fe.n_components()); + AssertDimension(immersed_c.size(), immersed_fe.n_components()); + + std::vector space_gtl(space_fe.n_components(), numbers::invalid_unsigned_int); + std::vector immersed_gtl(immersed_fe.n_components(), numbers::invalid_unsigned_int); + + for (unsigned int i=0, j=0; iget_dof_indices(dofs); + + const std::vector > &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::cell_iterator + ocell(*cells[c], &space_dh); + ocell->get_dof_indices(odofs); + for (unsigned int i=0; i + void create_coupling_mass_matrix(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 &space_mapping, + const Mapping &immersed_mapping) + { + 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"); + + const auto &space_fe = space_dh.get_fe(); + const auto &immersed_fe = immersed_dh.get_fe(); + + // Dof indices + 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 ? + ComponentMask(space_fe.n_components(), true) : + space_comps); + + const ComponentMask immersed_c + = (immersed_comps.size() == 0 ? + ComponentMask(immersed_fe.n_components(), true) : + immersed_comps); + + AssertDimension(space_c.size(), space_fe.n_components()); + AssertDimension(immersed_c.size(), immersed_fe.n_components()); + + std::vector space_gtl(space_fe.n_components(), numbers::invalid_unsigned_int); + std::vector immersed_gtl(immersed_fe.n_components(), numbers::invalid_unsigned_int); + + for (unsigned int i=0, j=0; i cell_matrix(space_dh.get_fe().dofs_per_cell, + immersed_dh.get_fe().dofs_per_cell); + + FEValues 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::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 > &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::active_cell_iterator + ocell(*cells[c], &space_dh); + 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, + 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, + const DoFHandler &immersed_dh, + const Quadrature &quad, + Sparsity &sparsity, + const ComponentMask &space_comps, + const ComponentMask &immersed_comps, + const Mapping &space_mapping, + const Mapping &immersed_mapping); +#endif +} + + +for (dim0 : DIMENSIONS; dim1 :DIMENSIONS; spacedim : SPACE_DIMENSIONS; + Matrix : SPARSE_MATRICES) +{ +#if dim1 <= dim0 && dim0 <= spacedim + template + void create_coupling_mass_matrix(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 &space_mapping, + const Mapping &immersed_mapping); +#endif +} -- 2.39.5