quadinterpolator.hpp

Go to the documentation of this file.

// Copyright (c) 2010-2026, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
 
#ifndef MFEM_QUADINTERP
#define MFEM_QUADINTERP
 
#include "fespace.hpp"
#include "kernel_dispatch.hpp"
 
namespace mfem
{
 
/** @brief A class that performs interpolation from an E-vector to quadrature
    point values and/or derivatives (Q-vectors). */
/** An E-vector represents the element-wise discontinuous version of the FE
    space and can be obtained, for example, from a GridFunction using the
    Operator returned by FiniteElementSpace::GetElementRestriction().
 
    The target quadrature points in the elements can be described either by an
    IntegrationRule (all mesh elements must be of the same type in this case) or
    by a QuadratureSpace. */
class QuadratureInterpolator
{
protected:
   friend class FiniteElementSpace; // Needs access to qspace and IntRule
 
   const FiniteElementSpace *fespace;  ///< Not owned
   const QuadratureSpace *qspace;      ///< Not owned
   const IntegrationRule *IntRule;     ///< Not owned
   mutable QVectorLayout q_layout;     ///< Output Q-vector layout
 
   mutable bool use_tensor_products;   ///< Tensor product evaluation mode
   mutable Vector d_buffer;            ///< Auxiliary device buffer
 
   /// Auxiliary method called by Mult() when using H(div)-conforming space
   void MultHDiv(const Vector &e_vec, unsigned eval_flags,
                 Vector &q_val, Vector &q_div) const;
 
public:
   static const int MAX_NQ2D = 100;
   static const int MAX_ND2D = 100;
   static const int MAX_VDIM2D = 3;
 
   static const int MAX_NQ3D = 1000;
   static const int MAX_ND3D = 1000;
   static const int MAX_VDIM3D = 3;
 
   enum EvalFlags
   {
      VALUES       = 1 << 0,  ///< Evaluate the values at quadrature points
      DERIVATIVES  = 1 << 1,  ///< Evaluate the derivatives at quadrature points
      /** @brief Assuming the derivative at quadrature points form a matrix,
          this flag can be used to compute and store their determinants. This
          flag can only be used in Mult(). */
      DETERMINANTS = 1 << 2,
      PHYSICAL_DERIVATIVES = 1 << 3, ///< Evaluate the physical derivatives
      /** Evaluate the values in physical space; for fields with
          FiniteElement::MapType other than FiniteElement::MapType::VALUE,
          such as H(div) and H(curl) elements, the physical values are different
          from the reference values. */
      PHYSICAL_VALUES = 1 << 4,
      /** For vector-valued fields, evaluate the magnitudes of the physical
          space vector values at quadrature points. */
      PHYSICAL_MAGNITUDES = 1 << 5
   };
 
   QuadratureInterpolator(const FiniteElementSpace &fes,
                          const IntegrationRule &ir);
 
   QuadratureInterpolator(const FiniteElementSpace &fes,
                          const QuadratureSpace &qs);
 
   /** @brief Disable the use of tensor product evaluations, for tensor-product
       elements, e.g. quads and hexes. By default, tensor product evaluations
       are enabled. */
   /** @sa EnableTensorProducts(), UsesTensorProducts(). */
   void DisableTensorProducts(bool disable = true) const
   { use_tensor_products = !disable; }
 
   /** @brief Enable the use of tensor product evaluations, for tensor-product
       elements, e.g. quads and hexes. By default, this option is enabled. */
   /** @sa DisableTensorProducts(), UsesTensorProducts(). */
   void EnableTensorProducts() const { use_tensor_products = true; }
 
   /** @brief Query the current tensor product evaluation mode. */
   /** @sa DisableTensorProducts(), EnableTensorProducts(). */
   bool UsesTensorProducts() const { return use_tensor_products; }
 
   /** @brief Query the current output Q-vector layout. The default value is
       QVectorLayout::byNODES. */
   /** @sa SetOutputLayout(). */
   QVectorLayout GetOutputLayout() const { return q_layout; }
 
   /** @brief Set the desired output Q-vector layout. The default value is
       QVectorLayout::byNODES. */
   /** @sa GetOutputLayout(). */
   void SetOutputLayout(QVectorLayout layout) const { q_layout = layout; }
 
   /// Interpolate the E-vector @a e_vec to quadrature points.
   /** The @a eval_flags are a bitwise mask of constants from the EvalFlags
       enumeration. When the VALUES flag is set, the values at quadrature points
       are computed and stored in the Vector @a q_val. Similarly, when one of
       the flags DERIVATIVES or PHYSICAL_DERIVATIVES is set, the derivatives
       (with respect to reference or physical coordinates, respectively) are
       computed and stored in @a q_der. Only one of the flags DERIVATIVES or
       PHYSICAL_DERIVATIVES can be set in a call. When the DETERMINANTS flag is
       set, it is assumed that the derivatives (with respect to reference
       coordinates) form a matrix at each quadrature point (i.e. the associated
       FiniteElementSpace is a vector space) and their determinants are computed
       and stored in @a q_det.
 
       For Integral spaces, the flags VALUES requests the computation of the
       scalar field values. The result is stored in @a q_val. Derivative types
       are not supported.
 
       For H(div)-conforming spaces, the flags VALUES / PHYSICAL_VALUES request
       the computation of the vector field values in reference or physical
       space, respectively. The flag PHYSICAL_MAGNITUDES requests the
       computation of the physical space magnitudes. In all 3 cases, the result
       is stored in @a q_val and therefore only one of the 3 cases can be
       requested in a single call.
 
       The layout of the input E-vector, @a e_vec, must be consistent with the
       evaluation mode: if tensor-product evaluations are enabled, then
       tensor-product elements, must use the ElementDofOrdering::LEXICOGRAPHIC
       layout; otherwise -- ElementDofOrdering::NATIVE layout. See
       FiniteElementSpace::GetElementRestriction(). */
   void Mult(const Vector &e_vec, unsigned eval_flags,
             Vector &q_val, Vector &q_der, Vector &q_det) const;
 
   /// Interpolate the values of the E-vector @a e_vec at quadrature points.
   void Values(const Vector &e_vec, Vector &q_val) const;
 
   /// @brief Interpolate the physical values of the E-vector @a e_vec at
   /// quadrature points.
   void PhysValues(const Vector &e_vec, Vector &q_val) const;
 
   /** @brief Interpolate the derivatives (with respect to reference
       coordinates) of the E-vector @a e_vec at quadrature points. */
   void Derivatives(const Vector &e_vec, Vector &q_der) const;
 
   /** @brief Interpolate the derivatives in physical space of the E-vector
       @a e_vec at quadrature points. */
   void PhysDerivatives(const Vector &e_vec, Vector &q_der) const;
 
   /** @brief Compute the determinants of the derivatives (with respect to
       reference coordinates) of the E-vector @a e_vec at quadrature points. */
   void Determinants(const Vector &e_vec, Vector &q_det) const;
 
   /// Perform the transpose operation of Mult(). (TODO)
   void MultTranspose(unsigned eval_flags, const Vector &q_val,
                      const Vector &q_der, Vector &e_vec) const;
 
   /// @brief Returns true if the given finite element space is supported by
   /// QuadratureInterpolator.
   static bool SupportsFESpace(const FiniteElementSpace &fespace);
 
   // value map types
   using TensorEvalKernelType = void (*)(const int ne, const real_t *B,
                                         const real_t *e_vec, real_t *q_val,
                                         const int vdim, const int nd,
                                         const int nq);
   using GradKernelType = void (*)(const int ne, const real_t *B,
                                   const real_t *G, const real_t *J,
                                   const real_t *e_vec, real_t *q_der,
                                   const int s_dim, const int v_dim,
                                   const int nd, const int nq);
   using CollocatedGradKernelType = void (*)(const int ne, const real_t *G,
                                             const real_t *J,
                                             const real_t *e_vec, real_t *q_der,
                                             const int sdim, const int vdim,
                                             const int d1d);
   using DetKernelType = void (*)(const int NE, const real_t *B,
                                  const real_t *G, const real_t *e_vec,
                                  real_t *q_det, const int nd, const int nq,
                                  Vector *d_buffer);
   using EvalKernelType = void (*)(const int NE, const int vdim,
                                   const QVectorLayout q_layout,
                                   const GeometricFactors *geom,
                                   const DofToQuad &maps, const Vector &e_vec,
                                   Vector &q_val, Vector &q_der, Vector &q_det,
                                   const int eval_flags);
 
   // integral map types
   using IntTensorEvalKernelType = void (*)(const int ne, const real_t *B,
                                            const real_t *detJ,
                                            const real_t *e_vec, real_t *q_val,
                                            const int vdim, const int nd,
                                            const int nq);
   using IntEvalKernelType =
      void (*)(const int NE, const int vdim, const QVectorLayout q_layout,
               const real_t *detJ, const GeometricFactors *geom,
               const DofToQuad &maps, const Vector &e_vec, Vector &q_val,
               Vector &q_der, Vector &q_det, const int eval_flags);
 
   using TensorEvalHDivKernelType =
      void(*)(const int, const real_t *, const real_t *, const real_t *,
              const real_t *, real_t *, const int, const int);
 
   // value-type mapping
   MFEM_REGISTER_KERNELS(TensorEvalKernels, TensorEvalKernelType,
                         (int, QVectorLayout, int, int, int), (int));
   MFEM_REGISTER_KERNELS(GradKernels, GradKernelType,
                         (int, QVectorLayout, bool, int, int, int), (int));
   MFEM_REGISTER_KERNELS(DetKernels, DetKernelType, (int, int, int, int));
   MFEM_REGISTER_KERNELS(EvalKernels, EvalKernelType, (int, int, int, int));
   MFEM_REGISTER_KERNELS(CollocatedGradKernels, CollocatedGradKernelType,
                         (int, QVectorLayout, bool, int, int), (int));
 
   // integral-type mapping
   MFEM_REGISTER_KERNELS(IntTensorEvalKernels, IntTensorEvalKernelType,
                         (int, QVectorLayout, int, int, int), (int));
   MFEM_REGISTER_KERNELS(IntEvalKernels, IntEvalKernelType, (int, int, int, int));
 
   MFEM_REGISTER_KERNELS(TensorEvalHDivKernels, TensorEvalHDivKernelType,
                         (int, QVectorLayout, unsigned, int, int));
 
   /// Adds specializations for TensorEvalKernels
   template <int DIM, QVectorLayout Q_LAYOUT, int VDIM, int D1D, int Q1D,
             int NBZ = 0>
   static void AddTensorEvalSpecializations()
   {
      if constexpr (NBZ)
      {
         IntTensorEvalKernels::Specialization<DIM, Q_LAYOUT, VDIM, D1D,
                              Q1D>::template Opt<NBZ>::Add();
         TensorEvalKernels::Specialization<DIM, Q_LAYOUT, VDIM, D1D,
                           Q1D>::template Opt<NBZ>::Add();
      }
      else if constexpr (NBZ == 0)
      {
         IntTensorEvalKernels::Specialization<DIM, Q_LAYOUT, VDIM, D1D,
                              Q1D>::Add();
         TensorEvalKernels::Specialization<DIM, Q_LAYOUT, VDIM, D1D,
                           Q1D>::Add();
      }
   }
 
   /// Adds specializations for EvalKernels
   template <int DIM, int VDIM, int ND, int NQ>
   static void AddEvalSpecializations()
   {
      IntEvalKernels::Specialization<DIM, VDIM, ND, NQ>::Add();
      EvalKernels::Specialization<DIM, VDIM, ND, NQ>::Add();
   }
 
   /// Adds specializations for GradKernels
   template <int DIM, QVectorLayout Q_LAYOUT, bool GRAD_PHYS, int VDIM, int D1D,
             int Q1D, int NBZ = 0>
   static void AddGradSpecializations()
   {
      if constexpr (NBZ)
      {
         GradKernels::Specialization<DIM, Q_LAYOUT, GRAD_PHYS, VDIM, D1D,
                     Q1D>::template Opt<NBZ>::Add();
      }
      else if constexpr (NBZ == 0)
      {
         GradKernels::Specialization<DIM, Q_LAYOUT, GRAD_PHYS, VDIM, D1D,
                     Q1D>::Add();
      }
   }
 
   /// Adds specializations for CollocatedGradKernels
   template <int DIM, QVectorLayout Q_LAYOUT, bool GRAD_PHYS, int VDIM, int D1D,
             int NBZ = 0>
   static void AddCollocatedGradSpecializations()
   {
      if constexpr (NBZ)
      {
         CollocatedGradKernels::Specialization<DIM, Q_LAYOUT, GRAD_PHYS, VDIM,
                               D1D>::template Opt<NBZ>::Add();
      }
      else if constexpr (NBZ == 0)
      {
         CollocatedGradKernels::Specialization<DIM, Q_LAYOUT, GRAD_PHYS, VDIM,
                               D1D>::Add();
      }
   }
 
   /// Adds specializations for DetKernels
   template <int DIM, int SDIM, int D1D, int Q1D>
   static void AddDetSpecializations()
   {
      DetKernels::Specialization<DIM, SDIM, D1D, Q1D>::Add();
   }
};
 
}
 
#endif