FiniteElement: Implement lazy restriction/prolongation via dealii::Lazy<>

[tamiko]

This PR changes the signature of the restriction and prolongation matrices to Lazy<FullMatrix> and introduces lazy initialization for all finite elements.

TODO

Closes #16201 Closes #16169

[tamiko]

@bangerth @kronbichler While we're at it I suggest we clean up initialization of restriction/prolongation matrices for all finite elements.

I am suggesting the following changes to the (primarily internal) interface of the FiniteElement base class:

public:
  // Changed: not virtual
  // Performs lazy initialization via construct_restriction_matrix()
  const FullMatrix<double> &
  get_restriction_matrix(const unsigned int         child,
                         const RefinementCase<dim> &refinement_case) const;

  bool restriction_is_implemented() const;
  bool isotropic_restriction_is_implemented() const;
  bool restriction_is_additive(const unsigned int index) const;

  // Changed: not virtual
  // Performs lazy initialization via construct_prolongation_matrix()
  const FullMatrix<double> &
  get_prolongation_matrix(const unsigned int         child,
                          const RefinementCase<dim> &refinement_case) const;

  bool prolongation_is_implemented() const;
  bool isotropic_prolongation_is_implemented() const;

protected:
  // Removed
  void reinit_restriction_and_prolongation_matrices( bool, bool);

  // New: derived classes override this method
  virtual std::optional<FullMatrix<double>>
  construct_restriction_matrix(
    const unsigned int         child,
    const RefinementCase<dim> &refinement_case) const;

  // New: derived classes override this method
  virtual std::optional<FullMatrix<double>>
  construct_prolongation_matrix(
    const unsigned int         child,
    const RefinementCase<dim> &refinement_case) const;

  // Should not be manipulated directly. Shall we mark them `private:`?
  std::vector<std::vector<Lazy<FullMatrix<double>>>> restriction;
  std::vector<std::vector<Lazy<FullMatrix<double>>>> prolongation;

This will be a breaking change affecting "externally" implement finite elements. But I think that having a standardized interface where most of the difficult logic resides in FiniteElement is worth the price.

[bangerth]

I like the new interface. But did you mean to let the construct_*() functions return a std::optional? I guess you're thinking of elements that aren't implementing these matrices?

[tamiko]

@bangerth I introduced the std::optional as return value for construct_* to entertain the idea how to introduce a failure mode without throwing exceptions and without auxiliary helper functions returning a bool. If we go that route then we can actually keep the *_is_implemented() variants.

If we decide to deprecate the *_is_implemented() booleans then I agree that the construct_* functions should simply return a FullMatrix and otherwise abort with an Assert. I prefer the latter.

[tamiko]

@bangerth Yes, I have not ported anything to the new interface yet :smile:

[bangerth]

On 10/27/23 21:35, Matthias Maier wrote:

@bangerth https://github.com/bangerth I introduced the |std::optional| as return value for |construct_| to entertain the idea how to introduce a failure mode without throwing exceptions and without auxiliary helper functions returning a bool. If we go that route then we can actually keep the |_is_implemented()| variants.

If we decide to deprecate the |is_implemented()| booleans then I agree that the |construct| functions should simply return a FullMatrix and otherwise abort with an |Assert|. I prefer the latter.

Yes, that sounds reasonable. Alternatively, you could have the default implementation of the construct_*() functions return an empty matrix and let that be an indication that a class does not implement these matrices (unless dofs_per_cell==0).

[drwells]

Would this also fix #331?

[bangerth]

Ha, funny. I reported #331 many years ago, and then I found the same problem again in #16169 and reported it with basically the same wording :-)

[bangerth]

@tamiko Are you still on this?