From a5991d642b3519251949f475e13b422446e7c45d Mon Sep 17 00:00:00 2001 From: kronbichler Date: Sun, 23 Feb 2014 12:50:21 +0000 Subject: [PATCH] Lazy initialization of projection and restriction matrices in FE_DGQ git-svn-id: https://svn.dealii.org/trunk@32531 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/base/aligned_vector.h | 1 + deal.II/include/deal.II/fe/fe_dgq.h | 52 +++++ deal.II/source/fe/fe.cc | 8 +- deal.II/source/fe/fe_dgq.cc | 187 +++++++++++++----- 4 files changed, 202 insertions(+), 46 deletions(-) diff --git a/deal.II/include/deal.II/base/aligned_vector.h b/deal.II/include/deal.II/base/aligned_vector.h index 0876e0f6fe..7c8241a273 100644 --- 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 #include +#include #if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0 #include diff --git a/deal.II/include/deal.II/fe/fe_dgq.h b/deal.II/include/deal.II/fe/fe_dgq.h index c10aaa8ad8..103c135c7e 100644 --- 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 #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -179,6 +180,52 @@ public: const unsigned int subface, FullMatrix &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 & + get_restriction_matrix (const unsigned int child, + const RefinementCase &refinement_case=RefinementCase::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 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; + /** * @name Functions to support hp * @{ @@ -368,6 +415,11 @@ private: void rotate_indices (std::vector &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 165e5da3ab..8237022f8a 100644 --- a/deal.II/source/fe/fe.cc +++ b/deal.II/source/fe/fe.cc @@ -311,10 +311,14 @@ FiniteElement::reinit_restriction_and_prolongation_matrices ( for (unsigned int i=0; i::isotropic_refinement) + if (this->restriction[ref_case-1][i].m() != this->dofs_per_cell + && + (!isotropic_restriction_only || ref_case==RefinementCase::isotropic_refinement)) this->restriction[ref_case-1][i].reinit (this->dofs_per_cell, this->dofs_per_cell); - if (!isotropic_prolongation_only || ref_case==RefinementCase::isotropic_refinement) + if (this->prolongation[ref_case-1][i].m() != this->dofs_per_cell + && + (!isotropic_prolongation_only || ref_case==RefinementCase::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 a61562b881..fd63197417 100644 --- a/deal.II/source/fe/fe_dgq.cc +++ b/deal.II/source/fe/fe_dgq.cc @@ -142,39 +142,7 @@ FE_DGQ::FE_DGQ (const unsigned int degree) std::vector(FiniteElementData( get_dpo_vector(degree),1, degree).dofs_per_cell, std::vector(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 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::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::FE_DGQ (const Quadrature<1> &points) std::vector(FiniteElementData( get_dpo_vector(points.size()-1),1, points.size()-1).dofs_per_cell, std::vector(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 (*this, this->prolongation); - // Fill restriction matrices with L2-projection - FETools::compute_projection_matrices (*this, this->restriction); - // Compute support points, which // are the tensor product of the // Lagrange interpolation points in // the constructor. Quadrature 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 &x_source_f +template +const FullMatrix & +FE_DGQ +::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_DGQ &this_nonconst = const_cast& >(*this); + if (refinement_case == RefinementCase::isotropic_refinement) + { + std::vector > > + isotropic_matrices(RefinementCase::isotropic_refinement); + isotropic_matrices.back(). + resize(GeometryInfo::n_children(RefinementCase(refinement_case)), + FullMatrix(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(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 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 +const FullMatrix & +FE_DGQ +::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(refinement_case), + ExcIndexRange(child,0,GeometryInfo::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 &this_nonconst = const_cast& >(*this); + if (refinement_case == RefinementCase::isotropic_refinement) + { + std::vector > > + isotropic_matrices(RefinementCase::isotropic_refinement); + isotropic_matrices.back(). + resize(GeometryInfo::n_children(RefinementCase(refinement_case)), + FullMatrix(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(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 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 bool FE_DGQ::hp_constraints_are_implemented () const -- 2.39.5