Construct restriction and prolongation matrices on first request. This reduces constr...

Log changes.

diff --git a/doc/news/changes.h b/doc/news/changes.h
index a132eace3a3324e8210c8324077131ec36adc278..83814bf0ff497e72ca3cd9be381a0947626f3256 100644 (file)
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -277,6 +277,14 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> Optimize construction of high-order FE_Nedelec by moving out some 
+  non-essential computations. Namely, construct restriction and prolongation 
+  matrices on first request. This reduces time spent in FE_Nedelec constructor
+  substantially.
+  <br>
+  (Alexander Grayver, 2014/08/22)
+  </li>
+
   <li> Changed: The functions GridTools::extract_boundary_mesh() and
   GridTools::create_union_triangulation() have been moved to
   GridGenerator::extract_boundary_mesh() and

diff --git a/include/deal.II/fe/fe_nedelec.h b/include/deal.II/fe/fe_nedelec.h
index 35c4f73d079f02a75ff21bb04af587e898630595..2017b31531adfd0cbf710e6002e1df12c3162734 100644 (file)
--- a/include/deal.II/fe/fe_nedelec.h
+++ b/include/deal.II/fe/fe_nedelec.h
@@ -25,6 +25,7 @@
 #include <deal.II/base/polynomial.h>
 #include <deal.II/base/tensor_product_polynomials.h>
 #include <deal.II/base/geometry_info.h>
+#include <deal.II/base/thread_management.h>
 #include <deal.II/fe/fe.h>
 #include <deal.II/fe/fe_poly_tensor.h>
 #include <vector>
@@ -271,6 +272,47 @@ public:
   get_subface_interpolation_matrix (const FiniteElement<dim> &source,
                                     const unsigned int subface,
                                     FullMatrix<double> &matrix) const;
+  /**
+   * Projection from a fine grid space onto a coarse grid space. If this
+   * projection operator is associated with a matrix @p P, then the
+   * restriction of this matrix @p P_i to a single child cell is returned
+   * here.
+   *
+   * The matrix @p P is the concatenation or the sum of the cell matrices @p
+   * P_i, depending on the #restriction_is_additive_flags. This distinguishes
+   * interpolation (concatenation) and projection with respect to scalar
+   * products (summation).
+   *
+   * Row and column indices are related to coarse grid and fine grid spaces,
+   * respectively, consistent with the definition of the associated operator.
+   */
+  virtual const FullMatrix<double> &
+  get_restriction_matrix (const unsigned int child,
+                          const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;
+
+  /**
+   * Embedding matrix between grids.
+   *
+   * The identity operator from a coarse grid space into a fine grid space is
+   * associated with a matrix @p P. The restriction of this matrix @p P_i to a
+   * single child cell is returned here.
+   *
+   * The matrix @p P is the concatenation, not the sum of the cell matrices @p
+   * P_i. That is, if the same non-zero entry <tt>j,k</tt> exists in in two
+   * different child matrices @p P_i, the value should be the same in both
+   * matrices and it is copied into the matrix @p P only once.
+   *
+   * Row and column indices are related to fine grid and coarse grid spaces,
+   * respectively, consistent with the definition of the associated operator.
+   *
+   * These matrices are used by routines assembling the prolongation matrix
+   * for multi-level methods.  Upon assembling the transfer matrix between
+   * cells using this matrix array, zero elements in the prolongation matrix
+   * are discarded and will not fill up the transfer matrix.
+   */
+  virtual const FullMatrix<double> &
+  get_prolongation_matrix (const unsigned int child,
+                           const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;
 
   virtual void interpolate (std::vector<double> &local_dofs,
                             const std::vector<double> &values) const;
@@ -407,6 +449,11 @@ private:
    */
   Table<2, double> boundary_weights;
 
+  /*
+   * Mutex for protecting initialization of restriction and embedding matrix.
+   */
+  mutable Threads::Mutex mutex;
+
   /**
    * Allow access from other
    * dimensions.

diff --git a/source/fe/fe_nedelec.cc b/source/fe/fe_nedelec.cc
index 9fd9e9f7fd2332a60a234a48dcf7746f5e9c48ba..5fc8b868b934c84628d49cc179143f4d4d09154d 100644 (file)
--- a/source/fe/fe_nedelec.cc
+++ b/source/fe/fe_nedelec.cc
@@ -74,21 +74,9 @@ FE_Nedelec<dim>::FE_Nedelec (const unsigned int p) :
   // will be the correct ones, not
   // the raw shape functions anymore.
 
-  // Reinit the vectors of
-  // restriction and prolongation
-  // matrices to the right sizes.
-  // Restriction only for isotropic
-  // refinement
-#ifdef DEBUG_NEDELEC
-  deallog << "Embedding" << std::endl;
-#endif
-  this->reinit_restriction_and_prolongation_matrices ();
-  // Fill prolongation matrices with embedding operators
-  FETools::compute_embedding_matrices (*this, this->prolongation, true);
-#ifdef DEBUG_NEDELEC
-  deallog << "Restriction" << std::endl;
-#endif
-  initialize_restriction ();
+  // do not initialize embedding and restriction here. these matrices are
+  // initialized on demand in get_restriction_matrix and
+  // get_prolongation_matrix
 
 #ifdef DEBUG_NEDELEC
   deallog << "Face Embedding" << std::endl;
@@ -2967,7 +2955,105 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix(
     }
 }
 
+template <int dim>
+const FullMatrix<double> &
+FE_Nedelec<dim>
+::get_prolongation_matrix (const unsigned int child,
+                           const RefinementCase<dim> &refinement_case) const
+{
+  Assert (refinement_case<RefinementCase<dim>::isotropic_refinement+1,
+          ExcIndexRange(refinement_case,0,RefinementCase<dim>::isotropic_refinement+1));
+  Assert (refinement_case!=RefinementCase<dim>::no_refinement,
+          ExcMessage("Prolongation matrices are only available for refined cells!"));
+  Assert (child<GeometryInfo<dim>::n_children(refinement_case),
+          ExcIndexRange(child,0,GeometryInfo<dim>::n_children(refinement_case)));
+
+  // initialization upon first request
+  if (this->prolongation[refinement_case-1][child].n() == 0)
+    {
+      Threads::Mutex::ScopedLock lock(this->mutex);
+
+      // if matrix got updated while waiting for the lock
+      if (this->prolongation[refinement_case-1][child].n() ==
+          this->dofs_per_cell)
+        return this->prolongation[refinement_case-1][child];
+
+      // now do the work. need to get a non-const version of data in order to
+      // be able to modify them inside a const function
+      FE_Nedelec<dim> &this_nonconst = const_cast<FE_Nedelec<dim>& >(*this);
+
+      // Reinit the vectors of
+      // restriction and prolongation
+      // matrices to the right sizes.
+      // Restriction only for isotropic
+      // refinement
+#ifdef DEBUG_NEDELEC
+  deallog << "Embedding" << std::endl;
+#endif
+      this_nonconst.reinit_restriction_and_prolongation_matrices ();
+      // Fill prolongation matrices with embedding operators
+      FETools::compute_embedding_matrices (this_nonconst, this_nonconst.prolongation, true);
+#ifdef DEBUG_NEDELEC
+  deallog << "Restriction" << std::endl;
+#endif
+      this_nonconst.initialize_restriction ();
+    }
+
+  // we use refinement_case-1 here. the -1 takes care of the origin of the
+  // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
+  // available and so the vector indices are shifted
+  return this->prolongation[refinement_case-1][child];
+}
+
+template <int dim>
+const FullMatrix<double> &
+FE_Nedelec<dim>
+::get_restriction_matrix (const unsigned int child,
+                          const RefinementCase<dim> &refinement_case) const
+{
+  Assert (refinement_case<RefinementCase<dim>::isotropic_refinement+1,
+          ExcIndexRange(refinement_case,0,RefinementCase<dim>::isotropic_refinement+1));
+  Assert (refinement_case!=RefinementCase<dim>::no_refinement,
+          ExcMessage("Restriction matrices are only available for refined cells!"));
+  Assert (child<GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)),
+          ExcIndexRange(child,0,GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case))));
+
+  // initialization upon first request
+  if (this->restriction[refinement_case-1][child].n() == 0)
+    {
+      Threads::Mutex::ScopedLock lock(this->mutex);
+
+      // if matrix got updated while waiting for the lock...
+      if (this->restriction[refinement_case-1][child].n() ==
+          this->dofs_per_cell)
+        return this->restriction[refinement_case-1][child];
+
+      // now do the work. need to get a non-const version of data in order to
+      // be able to modify them inside a const function
+      FE_Nedelec<dim> &this_nonconst = const_cast<FE_Nedelec<dim>& >(*this);
+
+      // Reinit the vectors of
+      // restriction and prolongation
+      // matrices to the right sizes.
+      // Restriction only for isotropic
+      // refinement
+#ifdef DEBUG_NEDELEC
+  deallog << "Embedding" << std::endl;
+#endif
+      this_nonconst.reinit_restriction_and_prolongation_matrices ();
+      // Fill prolongation matrices with embedding operators
+      FETools::compute_embedding_matrices (this_nonconst, this_nonconst.prolongation, true);
+#ifdef DEBUG_NEDELEC
+  deallog << "Restriction" << std::endl;
+#endif
+      this_nonconst.initialize_restriction ();
+    }
 
+  // we use refinement_case-1 here. the -1 takes care of the origin of the
+  // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
+  // available and so the vector indices are shifted
+  return this->restriction[refinement_case-1][child];
+}
 
 // Since this is a vector valued element,
 // we cannot interpolate a scalar function.