From 09af97b677555e5e4942ef5dd79f2bd09c4472cf Mon Sep 17 00:00:00 2001 From: Alexander Grayver Date: Thu, 21 Aug 2014 20:52:10 +0200 Subject: [PATCH] Construct restriction and prolongation matrices on first request. This reduces construction time of high-order FE_Nedelec substantially. Log changes. --- doc/news/changes.h | 8 +++ include/deal.II/fe/fe_nedelec.h | 47 +++++++++++++ source/fe/fe_nedelec.cc | 116 +++++++++++++++++++++++++++----- 3 files changed, 156 insertions(+), 15 deletions(-) diff --git a/doc/news/changes.h b/doc/news/changes.h index a132eace3a..83814bf0ff 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -277,6 +277,14 @@ inconvenience this causes.

Specific improvements

    +
  1. 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. +
    + (Alexander Grayver, 2014/08/22) +
    2. +
  2. 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 35c4f73d07..2017b31531 100644 --- a/include/deal.II/fe/fe_nedelec.h +++ b/include/deal.II/fe/fe_nedelec.h @@ -25,6 +25,7 @@ #include #include #include +#include #include #include #include @@ -271,6 +272,47 @@ public: get_subface_interpolation_matrix (const FiniteElement &source, const unsigned int subface, FullMatrix &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 & + get_restriction_matrix (const unsigned int child, + const RefinementCase &refinement_case=RefinementCase::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 j,k 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 & + get_prolongation_matrix (const unsigned int child, + const RefinementCase &refinement_case=RefinementCase::isotropic_refinement) const; virtual void interpolate (std::vector &local_dofs, const std::vector &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 9fd9e9f7fd..5fc8b868b9 100644 --- a/source/fe/fe_nedelec.cc +++ b/source/fe/fe_nedelec.cc @@ -74,21 +74,9 @@ FE_Nedelec::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::get_subface_interpolation_matrix( } } +template +const FullMatrix & +FE_Nedelec +::get_prolongation_matrix (const unsigned int child, + const RefinementCase &refinement_case) const +{ + Assert (refinement_case::isotropic_refinement+1, + ExcIndexRange(refinement_case,0,RefinementCase::isotropic_refinement+1)); + Assert (refinement_case!=RefinementCase::no_refinement, + ExcMessage("Prolongation matrices are only available for refined cells!")); + Assert (child::n_children(refinement_case), + ExcIndexRange(child,0,GeometryInfo::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 &this_nonconst = const_cast& >(*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::no_refinement (=0) there is no data + // available and so the vector indices are shifted + return this->prolongation[refinement_case-1][child]; +} + +template +const FullMatrix & +FE_Nedelec +::get_restriction_matrix (const unsigned int child, + const RefinementCase &refinement_case) const +{ + Assert (refinement_case::isotropic_refinement+1, + ExcIndexRange(refinement_case,0,RefinementCase::isotropic_refinement+1)); + Assert (refinement_case!=RefinementCase::no_refinement, + ExcMessage("Restriction matrices are only available for refined cells!")); + Assert (child::n_children(RefinementCase(refinement_case)), + ExcIndexRange(child,0,GeometryInfo::n_children(RefinementCase(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 &this_nonconst = const_cast& >(*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::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. -- 2.39.5