Lazy initialization of projection and restriction matrices in FE_DGQ

git-svn-id: https://svn.dealii.org/trunk@32531 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/include/deal.II/base/aligned_vector.h b/deal.II/include/deal.II/base/aligned_vector.h
index 0876e0f6fe10036dbc875aa355a3bee24cadabc3..7c8241a2734813858309a13cdc82ef8bf37708c1 100644 (file)
--- a/deal.II/include/deal.II/base/aligned_vector.h
+++ b/deal.II/include/deal.II/base/aligned_vector.h
@@ -24,6 +24,7 @@
 #include <deal.II/base/memory_consumption.h>
 #include <deal.II/base/parallel.h>
 
+#include <cstring>
 
 #if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0
 #include <mm_malloc.h>

diff --git a/deal.II/include/deal.II/fe/fe_dgq.h b/deal.II/include/deal.II/fe/fe_dgq.h
index c10aaa8ad8c09a9b8d1b529612a1a98d85459512..103c135c7e6c2257da58a277bcd5b85d21c21040 100644 (file)
--- a/deal.II/include/deal.II/fe/fe_dgq.h
+++ b/deal.II/include/deal.II/fe/fe_dgq.h
@@ -19,6 +19,7 @@
 
 #include <deal.II/base/config.h>
 #include <deal.II/base/tensor_product_polynomials.h>
+#include <deal.II/base/thread_management.h>
 #include <deal.II/fe/fe_poly.h>
 
 DEAL_II_NAMESPACE_OPEN
@@ -179,6 +180,52 @@ public:
                                     const unsigned int        subface,
                                     FullMatrix<double>       &matrix) const;
 
+  /**
+   * Projection from a fine grid space onto a coarse grid space. Overrides the
+   * respective method in FiniteElement, implementing lazy evaluation
+   * (initialize when requested).
+   *
+   * 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. Overrides the respective method in
+   * FiniteElement, implementing lazy evaluation (initialize when queried).
+   *
+   * 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;
+
   /**
    * @name Functions to support hp
    * @{
@@ -368,6 +415,11 @@ private:
   void rotate_indices (std::vector<unsigned int> &indices,
                        const char                 direction) const;
 
+  /*
+   * Mutex for protecting initialization of restriction and embedding matrix.
+   */
+  mutable Threads::Mutex mutex;
+
   /**
    * Allow access from other dimensions.
    */

diff --git a/deal.II/source/fe/fe.cc b/deal.II/source/fe/fe.cc
index 165e5da3ab647cd637bcebe10c26a9388e4ce7d2..8237022f8ab76c8d9395ce0ad03635e11c2ba39b 100644 (file)
--- a/deal.II/source/fe/fe.cc
+++ b/deal.II/source/fe/fe.cc
@@ -311,10 +311,14 @@ FiniteElement<dim,spacedim>::reinit_restriction_and_prolongation_matrices (
 
       for (unsigned int i=0; i<nc; ++i)
         {
-          if (!isotropic_restriction_only || ref_case==RefinementCase<dim>::isotropic_refinement)
+          if (this->restriction[ref_case-1][i].m() != this->dofs_per_cell
+              &&
+              (!isotropic_restriction_only || ref_case==RefinementCase<dim>::isotropic_refinement))
             this->restriction[ref_case-1][i].reinit (this->dofs_per_cell,
                                                      this->dofs_per_cell);
-          if (!isotropic_prolongation_only || ref_case==RefinementCase<dim>::isotropic_refinement)
+          if (this->prolongation[ref_case-1][i].m() != this->dofs_per_cell
+              &&
+              (!isotropic_prolongation_only || ref_case==RefinementCase<dim>::isotropic_refinement))
             this->prolongation[ref_case-1][i].reinit (this->dofs_per_cell,
                                                       this->dofs_per_cell);
         }

diff --git a/deal.II/source/fe/fe_dgq.cc b/deal.II/source/fe/fe_dgq.cc
index a61562b881b4aa7034c5f859991a4591bc952e96..fd6319741748e7818067745cc75f77a27a08d80b 100644 (file)
--- a/deal.II/source/fe/fe_dgq.cc
+++ b/deal.II/source/fe/fe_dgq.cc
@@ -142,39 +142,7 @@ FE_DGQ<dim, spacedim>::FE_DGQ (const unsigned int degree)
     std::vector<ComponentMask>(FiniteElementData<dim>(
                                  get_dpo_vector(degree),1, degree).dofs_per_cell, std::vector<bool>(1,true)))
 {
-  // Reinit the vectors of
-  // restriction and prolongation
-  // matrices to the right sizes
-  this->reinit_restriction_and_prolongation_matrices();
-
-  // Fill prolongation matrices with
-  // embedding operators and
-  // restriction with L2 projection
-  //
-  // we can call the respective
-  // functions in the case
-  // dim==spacedim, but they are not
-  // currently implemented for the
-  // codim-1 case. there's no harm
-  // here, though, since these
-  // matrices are independent of the
-  // space dimension and so we can
-  // copy them from the respective
-  // elements (not cheap, but works)
-  if (dim == spacedim)
-    {
-      FETools::compute_embedding_matrices (*this, this->prolongation);
-      FETools::compute_projection_matrices (*this, this->restriction);
-    }
-  else
-    {
-      FE_DGQ<dim> tmp (degree);
-      this->prolongation = tmp.prolongation;
-      this->restriction  = tmp.restriction;
-    }
-
-
-  // finally fill in support points
+  // fill in support points
   if (degree == 0)
     {
       // constant elements, take
@@ -210,8 +178,11 @@ FE_DGQ<dim, spacedim>::FE_DGQ (const unsigned int degree)
             };
     };
 
-  // note: no face support points for
-  // DG elements
+  // do not initialize embedding and restriction here. these matrices are
+  // initialized on demand in get_restriction_matrix and
+  // get_prolongation_matrix
+
+  // note: no face support points for DG elements
 }
 
 
@@ -226,21 +197,17 @@ FE_DGQ<dim, spacedim>::FE_DGQ (const Quadrature<1> &points)
     std::vector<ComponentMask>(FiniteElementData<dim>(
                                  get_dpo_vector(points.size()-1),1, points.size()-1).dofs_per_cell, std::vector<bool>(1,true)))
 {
-  // Reinit the vectors of
-  // restriction and prolongation
-  // matrices to the right sizes
-  this->reinit_restriction_and_prolongation_matrices();
-  // Fill prolongation matrices with embedding operators
-  FETools::compute_embedding_matrices<dim, double, spacedim> (*this, this->prolongation);
-  // Fill restriction matrices with L2-projection
-  FETools::compute_projection_matrices<dim, double, spacedim> (*this, this->restriction);
-
   // Compute support points, which
   // are the tensor product of the
   // Lagrange interpolation points in
   // the constructor.
   Quadrature<dim> support_quadrature(points);
   this->unit_support_points = support_quadrature.get_points();
+
+
+  // do not initialize embedding and restriction here. these matrices are
+  // initialized on demand in get_restriction_matrix and
+  // get_prolongation_matrix
 }
 
 
@@ -509,6 +476,138 @@ get_subface_interpolation_matrix (const FiniteElement<dim, spacedim> &x_source_f
 
 
 
+template <int dim, int spacedim>
+const FullMatrix<double> &
+FE_DGQ<dim,spacedim>
+::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_DGQ<dim,spacedim> &this_nonconst = const_cast<FE_DGQ<dim,spacedim>& >(*this);
+      if (refinement_case == RefinementCase<dim>::isotropic_refinement)
+        {
+          std::vector<std::vector<FullMatrix<double> > >
+            isotropic_matrices(RefinementCase<dim>::isotropic_refinement);
+          isotropic_matrices.back().
+            resize(GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)),
+                   FullMatrix<double>(this->dofs_per_cell, this->dofs_per_cell));
+          if (dim == spacedim)
+            FETools::compute_embedding_matrices (*this, isotropic_matrices, true);
+          else
+            FETools::compute_embedding_matrices (FE_DGQ<dim>(this->degree),
+                                                 isotropic_matrices, true);
+          this_nonconst.prolongation[refinement_case-1].swap(isotropic_matrices.back());
+        }
+      else
+        {
+          // must compute both restriction and prolongation matrices because
+          // we only check for their size and the reinit call initializes them
+          // all
+          this_nonconst.reinit_restriction_and_prolongation_matrices();
+          if (dim == spacedim)
+            {
+              FETools::compute_embedding_matrices (*this, this_nonconst.prolongation);
+              FETools::compute_projection_matrices (*this, this_nonconst.restriction);
+            }
+          else
+            {
+              FE_DGQ<dim> tmp(this->degree);
+              FETools::compute_embedding_matrices (tmp, this_nonconst.prolongation);
+              FETools::compute_projection_matrices (tmp, this_nonconst.restriction);
+            }
+        }
+    }
+
+  // finally return the matrix
+  return this->prolongation[refinement_case-1][child];
+}
+
+
+
+template <int dim, int spacedim>
+const FullMatrix<double> &
+FE_DGQ<dim,spacedim>
+::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(refinement_case),
+          ExcIndexRange(child,0,GeometryInfo<dim>::n_children(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_DGQ<dim,spacedim> &this_nonconst = const_cast<FE_DGQ<dim,spacedim>& >(*this);
+      if (refinement_case == RefinementCase<dim>::isotropic_refinement)
+        {
+          std::vector<std::vector<FullMatrix<double> > >
+            isotropic_matrices(RefinementCase<dim>::isotropic_refinement);
+          isotropic_matrices.back().
+            resize(GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)),
+                   FullMatrix<double>(this->dofs_per_cell, this->dofs_per_cell));
+          if (dim == spacedim)
+            FETools::compute_projection_matrices (*this, isotropic_matrices, true);
+          else
+            FETools::compute_projection_matrices (FE_DGQ<dim>(this->degree),
+                                                 isotropic_matrices, true);
+          this_nonconst.restriction[refinement_case-1].swap(isotropic_matrices.back());
+        }
+      else
+        {
+          // must compute both restriction and prolongation matrices because
+          // we only check for their size and the reinit call initializes them
+          // all
+          this_nonconst.reinit_restriction_and_prolongation_matrices();
+          if (dim == spacedim)
+            {
+              FETools::compute_embedding_matrices (*this, this_nonconst.prolongation);
+              FETools::compute_projection_matrices (*this, this_nonconst.restriction);
+            }
+          else
+            {
+              FE_DGQ<dim> tmp(this->degree);
+              FETools::compute_embedding_matrices (tmp, this_nonconst.prolongation);
+              FETools::compute_projection_matrices (tmp, this_nonconst.restriction);
+            }
+        }
+    }
+
+  // finally return the matrix
+  return this->restriction[refinement_case-1][child];
+}
+
+
+
 template <int dim, int spacedim>
 bool
 FE_DGQ<dim, spacedim>::hp_constraints_are_implemented () const