Let FETools::compute_node_matrix() return its result, rather than take it by reference.

diff --git a/include/deal.II/fe/fe_tools.h b/include/deal.II/fe/fe_tools.h
index 5eb1dcb976f3246de2ea643e354584ae246f706c..a4cf52596c903d49018962fe36764a92f3911401 100644 (file)
--- 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 <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

diff --git a/source/fe/fe_abf.cc b/source/fe/fe_abf.cc
index 1598d99317389cb2991925f3e7382186e929ee1c..97d04e44e6a0802c467caeeb9983791e0dab1075 100644 (file)
--- a/source/fe/fe_abf.cc
+++ b/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

diff --git a/source/fe/fe_bdm.cc b/source/fe/fe_bdm.cc
index 0312a36583f1f6120a8e590f0f2e18f6053fbe7c..120d2c88080fe97ae1391b8a393963346cf8511f 100644 (file)
--- a/source/fe/fe_bdm.cc
+++ b/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.

diff --git a/source/fe/fe_raviart_thomas.cc b/source/fe/fe_raviart_thomas.cc
index f18cf72bf6284394ad829e1bf805756081e77e30..e21f8db1edb140611c06330f42e5555ce41ec4c6 100644 (file)
--- a/source/fe/fe_raviart_thomas.cc
+++ b/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

diff --git a/source/fe/fe_raviart_thomas_nodal.cc b/source/fe/fe_raviart_thomas_nodal.cc
index 3556c6851ab00d97b457f4b76045e2ed2736ee20..6b9a1e1bc15b14e3e4d3f388943f05c5433ab4d1 100644 (file)
--- a/source/fe/fe_raviart_thomas_nodal.cc
+++ b/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

diff --git a/source/fe/fe_tools.cc b/source/fe/fe_tools.cc
index 7716e50d8619ac006341b412d50148bf2ab29328..65c4a14d128f79fbaf4b29e1330aaadc8e6b87d7 100644 (file)
--- a/source/fe/fe_tools.cc
+++ b/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;
   }
 
 

diff --git a/source/fe/fe_tools.inst.in b/source/fe/fe_tools.inst.in
index c3b04be5d9a03601a617b21508675c7b49ed7459..a1af55ec32a06fb336220c2ce3457ca109220e29 100644 (file)
--- 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<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(