Let FETools::compute_node_matrix() return its result, rather than tak…

…e it by reference.

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 <int dim, int spacedim>
-  void compute_node_matrix(FullMatrix<double> &M,
-                           const FiniteElement<dim,spacedim> &fe);
+  FullMatrix<double>
+  compute_node_matrix(const FiniteElement<dim,spacedim> &fe);
 
   /**
    * For all possible (isotropic and anisotropic) refinement cases compute the

source/fe/fe_abf.cc

@@ -62,18 +62,16 @@ FE_ABF<dim>::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<double> 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<double> 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

source/fe/fe_bdm.cc

@@ -57,27 +57,16 @@ FE_BDM<dim>::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<double> 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<double> 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.

source/fe/fe_raviart_thomas.cc

@@ -59,23 +59,16 @@ FE_RaviartThomas<dim>::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<double> 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<double> 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

source/fe/fe_raviart_thomas_nodal.cc

@@ -52,23 +52,16 @@ FE_RaviartThomasNodal<dim>::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<double> 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<double> 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

source/fe/fe_tools.cc

@@ -1612,15 +1612,14 @@ namespace FETools
 
 
   template<int dim, int spacedim>
-  void
-  compute_node_matrix(
-    FullMatrix<double> &N,
-    const FiniteElement<dim,spacedim> &fe)
+  FullMatrix<double>
+  compute_node_matrix(const FiniteElement<dim,spacedim> &fe)
   {
     const unsigned int n_dofs = fe.dofs_per_cell;
+
+    FullMatrix<double> 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<Point<dim> > &points = fe.get_generalized_support_points();
@@ -1662,6 +1661,8 @@ namespace FETools
             Assert (numbers::is_finite(local_dofs[j]), ExcInternalError());
           }
       }
+
+    return N;
   }
 
 

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<deal_II_dimension> (const std::string &);
 
     template
-    void compute_node_matrix(
-        FullMatrix<double> &,
-        const FiniteElement<deal_II_dimension> &);
+    FullMatrix<double>
+    compute_node_matrix(const FiniteElement<deal_II_dimension> &);
 
     template
     void compute_component_wise(