From 089ae8d06030d2c31fd0e867388c41bc6ea34a5c Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 16 Jan 2017 15:38:07 -0700 Subject: [PATCH] Let FETools::compute_node_matrix() return its result, rather than take it by reference. --- include/deal.II/fe/fe_tools.h | 8 ++++++-- source/fe/fe_abf.cc | 16 +++++++--------- source/fe/fe_bdm.cc | 25 +++++++------------------ source/fe/fe_raviart_thomas.cc | 23 ++++++++--------------- source/fe/fe_raviart_thomas_nodal.cc | 23 ++++++++--------------- source/fe/fe_tools.cc | 13 +++++++------ source/fe/fe_tools.inst.in | 5 ++--- 7 files changed, 45 insertions(+), 68 deletions(-) diff --git a/include/deal.II/fe/fe_tools.h b/include/deal.II/fe/fe_tools.h index 5eb1dcb976..a4cf52596c 100644 --- a/include/deal.II/fe/fe_tools.h +++ b/include/deal.II/fe/fe_tools.h @@ -291,10 +291,14 @@ namespace FETools * - That the finite element has exactly @p dim vector components. * - That the function $f_j$ is given by whatever the element implements * through the FiniteElement::interpolate() function. + * + * @param fe The finite element for which the operations above are to be + * performed. + * @return The matrix $X$ as discussed above. */ template - void compute_node_matrix(FullMatrix &M, - const FiniteElement &fe); + FullMatrix + compute_node_matrix(const FiniteElement &fe); /** * For all possible (isotropic and anisotropic) refinement cases compute the diff --git a/source/fe/fe_abf.cc b/source/fe/fe_abf.cc index 1598d99317..97d04e44e6 100644 --- a/source/fe/fe_abf.cc +++ b/source/fe/fe_abf.cc @@ -62,18 +62,16 @@ FE_ABF::FE_ABF (const unsigned int deg) // quadrature weights, since they // are required for interpolation. initialize_support_points(deg); - // Now compute the inverse node - //matrix, generating the correct - //basis functions from the raw - //ones. - FullMatrix M(n_dofs, n_dofs); - FETools::compute_node_matrix(M, *this); + // Now compute the inverse node matrix, generating the correct + // basis functions from the raw ones. For a discussion of what + // exactly happens here, see FETools::compute_node_matrix. + const FullMatrix M = FETools::compute_node_matrix(*this); this->inverse_node_matrix.reinit(n_dofs, n_dofs); this->inverse_node_matrix.invert(M); - // From now on, the shape functions - // will be the correct ones, not - // the raw shape functions anymore. + // From now on, the shape functions provided by FiniteElement::shape_value + // and similar functions will be the correct ones, not + // the raw shape functions from the polynomial space anymore. // Reinit the vectors of // restriction and prolongation diff --git a/source/fe/fe_bdm.cc b/source/fe/fe_bdm.cc index 0312a36583..120d2c8808 100644 --- a/source/fe/fe_bdm.cc +++ b/source/fe/fe_bdm.cc @@ -57,27 +57,16 @@ FE_BDM::FE_BDM (const unsigned int deg) // Set up the generalized support // points initialize_support_points (deg); - //Now compute the inverse node - //matrix, generating the correct - //basis functions from the raw - //ones. - - // We use an auxiliary matrix in - // this function. Therefore, - // inverse_node_matrix is still - // empty and shape_value_component - // returns the 'raw' shape values. - FullMatrix M(n_dofs, n_dofs); - FETools::compute_node_matrix(M, *this); - -// std::cout << std::endl; -// M.print_formatted(std::cout, 2, true); + // Now compute the inverse node matrix, generating the correct + // basis functions from the raw ones. For a discussion of what + // exactly happens here, see FETools::compute_node_matrix. + const FullMatrix M = FETools::compute_node_matrix(*this); this->inverse_node_matrix.reinit(n_dofs, n_dofs); this->inverse_node_matrix.invert(M); - // From now on, the shape functions - // will be the correct ones, not - // the raw shape functions anymore. + // From now on, the shape functions provided by FiniteElement::shape_value + // and similar functions will be the correct ones, not + // the raw shape functions from the polynomial space anymore. // Embedding errors become pretty large, so we just replace the // regular threshold in both "computing_..." functions by 1. diff --git a/source/fe/fe_raviart_thomas.cc b/source/fe/fe_raviart_thomas.cc index f18cf72bf6..e21f8db1ed 100644 --- a/source/fe/fe_raviart_thomas.cc +++ b/source/fe/fe_raviart_thomas.cc @@ -59,23 +59,16 @@ FE_RaviartThomas::FE_RaviartThomas (const unsigned int deg) // quadrature weights, since they // are required for interpolation. initialize_support_points(deg); - // Now compute the inverse node - //matrix, generating the correct - //basis functions from the raw - //ones. - - // We use an auxiliary matrix in - // this function. Therefore, - // inverse_node_matrix is still - // empty and shape_value_component - // returns the 'raw' shape values. - FullMatrix M(n_dofs, n_dofs); - FETools::compute_node_matrix(M, *this); + + // Now compute the inverse node matrix, generating the correct + // basis functions from the raw ones. For a discussion of what + // exactly happens here, see FETools::compute_node_matrix. + const FullMatrix M = FETools::compute_node_matrix(*this); this->inverse_node_matrix.reinit(n_dofs, n_dofs); this->inverse_node_matrix.invert(M); - // From now on, the shape functions - // will be the correct ones, not - // the raw shape functions anymore. + // From now on, the shape functions provided by FiniteElement::shape_value + // and similar functions will be the correct ones, not + // the raw shape functions from the polynomial space anymore. // Reinit the vectors of // restriction and prolongation diff --git a/source/fe/fe_raviart_thomas_nodal.cc b/source/fe/fe_raviart_thomas_nodal.cc index 3556c6851a..6b9a1e1bc1 100644 --- a/source/fe/fe_raviart_thomas_nodal.cc +++ b/source/fe/fe_raviart_thomas_nodal.cc @@ -52,23 +52,16 @@ FE_RaviartThomasNodal::FE_RaviartThomasNodal (const unsigned int deg) // quadrature weights, since they // are required for interpolation. initialize_support_points(deg); - // Now compute the inverse node - //matrix, generating the correct - //basis functions from the raw - //ones. - - // We use an auxiliary matrix in - // this function. Therefore, - // inverse_node_matrix is still - // empty and shape_value_component - // returns the 'raw' shape values. - FullMatrix M(n_dofs, n_dofs); - FETools::compute_node_matrix(M, *this); + + // Now compute the inverse node matrix, generating the correct + // basis functions from the raw ones. For a discussion of what + // exactly happens here, see FETools::compute_node_matrix. + const FullMatrix M = FETools::compute_node_matrix(*this); this->inverse_node_matrix.reinit(n_dofs, n_dofs); this->inverse_node_matrix.invert(M); - // From now on, the shape functions - // will be the correct ones, not - // the raw shape functions anymore. + // From now on, the shape functions provided by FiniteElement::shape_value + // and similar functions will be the correct ones, not + // the raw shape functions from the polynomial space anymore. // Reinit the vectors of // prolongation matrices to the diff --git a/source/fe/fe_tools.cc b/source/fe/fe_tools.cc index 7716e50d86..65c4a14d12 100644 --- a/source/fe/fe_tools.cc +++ b/source/fe/fe_tools.cc @@ -1612,15 +1612,14 @@ namespace FETools template - void - compute_node_matrix( - FullMatrix &N, - const FiniteElement &fe) + FullMatrix + compute_node_matrix(const FiniteElement &fe) { const unsigned int n_dofs = fe.dofs_per_cell; + + FullMatrix N (n_dofs, n_dofs); + Assert (fe.has_generalized_support_points(), ExcNotInitialized()); - Assert (N.n()==n_dofs, ExcDimensionMismatch(N.n(), n_dofs)); - Assert (N.m()==n_dofs, ExcDimensionMismatch(N.m(), n_dofs)); Assert (fe.n_components() == dim, ExcNotImplemented()); const std::vector > &points = fe.get_generalized_support_points(); @@ -1662,6 +1661,8 @@ namespace FETools Assert (numbers::is_finite(local_dofs[j]), ExcInternalError()); } } + + return N; } diff --git a/source/fe/fe_tools.inst.in b/source/fe/fe_tools.inst.in index c3b04be5d9..a1af55ec32 100644 --- a/source/fe/fe_tools.inst.in +++ b/source/fe/fe_tools.inst.in @@ -143,9 +143,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS get_fe_from_name (const std::string &); template - void compute_node_matrix( - FullMatrix &, - const FiniteElement &); + FullMatrix + compute_node_matrix(const FiniteElement &); template void compute_component_wise( -- 2.39.5