Support integrated Legendre polynomials in ChangeOfBasis

Implement for L2ElementRestriction (in this case, same as transpose)

fem/CMakeLists.txt

@@ -56,6 +56,7 @@ set(SRCS
   integ/lininteg_domain_vectorfe.cpp
   integ/nonlininteg_vecconvection_pa.cpp
   integ/nonlininteg_vecconvection_mf.cpp
+  change_basis.cpp
   coefficient.cpp
   complex_fem.cpp
   convergence.cpp
@@ -164,6 +165,7 @@ set(HDRS
   integ/bilininteg_hdiv_kernels.hpp
   integ/bilininteg_hcurlhdiv_kernels.hpp
   integ/bilininteg_mass_kernels.hpp
+  change_basis.hpp
   coefficient.hpp
   complex_fem.hpp
   convergence.hpp

fem/change_basis.cpp

@@ -0,0 +1,132 @@
+// Copyright (c) 2010-2024, 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.
+
+#include "change_basis.hpp"
+#include "qinterp/dispatch.hpp"
+
+namespace mfem
+{
+
+ChangeOfBasis::ChangeOfBasis(FiniteElementSpace &fes_, int dest_btype)
+   : Operator(fes_.GetVSize()),
+     fes(fes_)
+{
+   const FiniteElement *fe = fes.GetFE(0);
+   auto *tbe = dynamic_cast<const TensorBasisElement*>(fe);
+   MFEM_VERIFY(tbe != nullptr, "Must be a tensor element.");
+   const int source_btype = tbe->GetBasisType();
+
+   const int order = fes.GetElementOrder(0);
+   IntegrationRules irs(0, Quadrature1D::GaussLobatto);
+   const IntegrationRule &ir = irs.Get(Geometry::SEGMENT, 2*order - 1);
+
+   auto compute_vandermonde = [&](int btype, DenseMatrix &V)
+   {
+      if (btype < BasisType::NumBasisTypes)
+      {
+         Poly_1D::Basis &basis = poly1d.GetBasis(order, btype);
+         for (int i = 0; i < ir.Size(); ++i)
+         {
+            Vector col;
+            V.GetColumnReference(i, col);
+            basis.Eval(ir[i].x, col);
+         }
+      }
+      else if (btype == LEGENDRE)
+      {
+         for (int i = 0; i < ir.Size(); ++i)
+         {
+            Vector col;
+            V.GetColumnReference(i, col);
+            Poly_1D::CalcLegendre(order, ir[i].x, col.HostWrite());
+         }
+      }
+      else if (btype == INTEGRATED_LEGENDRE)
+      {
+         for (int i = 0; i < ir.Size(); ++i)
+         {
+            Vector col;
+            V.GetColumnReference(i, col);
+            Poly_1D::CalcIntegratedLegendre(order, ir[i].x, col.HostWrite());
+         }
+      }
+      else
+      {
+         MFEM_ABORT("");
+      }
+      V.Transpose();
+   };
+
+   DenseMatrix V1(order + 1, order + 1);
+   compute_vandermonde(source_btype, V1);
+
+   DenseMatrix V2(order + 1, order + 1);
+   compute_vandermonde(dest_btype, V2);
+
+   DenseMatrixInverse V2_inv(V2);
+
+   T1D.SetSize(order + 1, order + 1);
+   V2_inv.Mult(V1, T1D);
+
+   {
+      T1D_inv.SetSize(order + 1, order + 1);
+      DenseMatrix A = T1D;
+      DenseMatrixInverse A_inv(A);
+      A_inv.GetInverseMatrix(T1D_inv);
+   }
+
+   dof2quad.FE = fe;
+   dof2quad.mode = DofToQuad::TENSOR;
+   dof2quad.ndof = order + 1;
+   dof2quad.nqpt = order + 1;
+   dof2quad.B.SetSize((order + 1)*(order + 1));
+}
+
+void ChangeOfBasis::Mult_(const DenseMatrix &B1D, const Vector &x,
+                          Vector &y) const
+{
+   using namespace internal::quadrature_interpolator;
+
+   const auto ordering = ElementDofOrdering::LEXICOGRAPHIC;
+   const Operator *restr_op = fes.GetElementRestriction(ordering);
+   auto restr = dynamic_cast<const ElementRestrictionOperator*>(restr_op);
+   MFEM_VERIFY(restr != nullptr, "Unsupported element restriction type.");
+
+   x_e.SetSize(restr->Height());
+   y_e.SetSize(restr->Height());
+
+   restr->Mult(x, x_e);
+
+   dof2quad.B.GetMemory().CopyFrom(B1D.GetMemory(), dof2quad.B.Size());
+
+   if (fes.GetOrdering() == Ordering::byNODES)
+   {
+      TensorValues<QVectorLayout::byNODES>(fes.GetNE(), 1, dof2quad, x_e, y_e);
+   }
+   else
+   {
+      TensorValues<QVectorLayout::byVDIM>(fes.GetNE(), 1, dof2quad, x_e, y_e);
+   }
+
+   restr->MultLeftInverse(y_e, y);
+}
+
+void ChangeOfBasis::Mult(const Vector &x, Vector &y) const
+{
+   Mult_(T1D, x, y);
+}
+
+void ChangeOfBasis::MultInverse(const Vector &x, Vector &y) const
+{
+   Mult_(T1D_inv, x, y);
+}
+
+} // namespace mfem

fem/change_basis.hpp

@@ -0,0 +1,47 @@
+// Copyright (c) 2010-2024, 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_CHANGE_BASIS_HPP
+#define MFEM_CHANGE_BASIS_HPP
+
+#include "../linalg/operator.hpp"
+#include "fespace.hpp"
+
+namespace mfem
+{
+
+class ChangeOfBasis : public Operator
+{
+public:
+   enum
+   {
+      LEGENDRE            = BasisType::NumBasisTypes + 1,
+      INTEGRATED_LEGENDRE = BasisType::NumBasisTypes + 2
+   };
+protected:
+   FiniteElementSpace &fes;
+   mutable Vector x_e, y_e;
+
+   mutable DofToQuad dof2quad;
+   DenseMatrix T1D;
+   DenseMatrix T1D_inv;
+
+   void Mult_(const DenseMatrix &B1D, const Vector &x, Vector &y) const;
+
+public:
+   ChangeOfBasis(FiniteElementSpace &fes_, int dest_btype);
+   void Mult(const Vector &x, Vector &y) const override;
+   void MultInverse(const Vector &x, Vector &y) const;
+};
+
+} // namespace mfem
+
+#endif

fem/fe/fe_base.cpp

@@ -2201,6 +2201,21 @@ void Poly_1D::CalcLegendre(const int p, const real_t x, real_t *u, real_t *d)
    }
 }
 
+void Poly_1D::CalcIntegratedLegendre(const int p, const real_t x, real_t *u)
+{
+   Poly_1D::CalcLegendre(p, x, u);
+   for (int n = p; n > 1; n--)
+   {
+      u[n] = (u[n] - u[n - 2]) / 2;
+   }
+   for (int n = 1; n < p; ++n)
+   {
+      u[n] = u[n+1];
+   }
+   u[p] = x;
+   u[0] = 1 - x;
+}
+
 void Poly_1D::CalcChebyshev(const int p, const real_t x, real_t *u)
 {
    // recursive definition, z in [-1,1]

fem/fe/fe_base.hpp

@@ -1186,6 +1186,8 @@ public:
    static void CalcLegendre(const int p, const real_t x, real_t *u);
    static void CalcLegendre(const int p, const real_t x, real_t *u, real_t *d);
 
+   static void CalcIntegratedLegendre(const int p, const real_t x, real_t *u);
+
    ~Poly_1D();
 };
 

fem/fem.hpp

@@ -49,6 +49,7 @@
 #include "lor/lor.hpp"
 #include "dgmassinv.hpp"
 #include "hyperbolic.hpp"
+#include "change_basis.hpp"
 
 #ifdef MFEM_USE_MPI
 #include "pfespace.hpp"

fem/restriction.cpp

@@ -560,6 +560,13 @@ void L2ElementRestriction::AddMultTranspose(const Vector &x, Vector &y,
    TAddMultTranspose<ADD>(x, y);
 }
 
+void L2ElementRestriction::MultLeftInverse(const Vector &x, Vector &y) const
+{
+   // For L2 spaces, the restriction is a permutation, and so its inverse is
+   // equal to its transpose.
+   MultTranspose(x, y);
+}
+
 void L2ElementRestriction::FillI(SparseMatrix &mat) const
 {
    const int elem_dofs = ndof;

fem/restriction.hpp

@@ -32,6 +32,10 @@ public:
    /// of freedom @a y.
    void AddMultTranspose(const Vector &x, Vector &y,
                          const real_t a = 1.0) const override = 0;
+
+   /// @brief Compute MultTranspose by setting (rather than adding) element
+   /// contributions; this is a left inverse of the Mult() operation
+   virtual void MultLeftInverse(const Vector &x, Vector &y) const = 0;
 };
 
 /// Operator that converts FiniteElementSpace L-vectors to E-vectors.
@@ -62,15 +66,13 @@ public:
    void MultTranspose(const Vector &x, Vector &y) const override;
    void AddMultTranspose(const Vector &x, Vector &y,
                          const real_t a = 1.0) const override;
+   void MultLeftInverse(const Vector &x, Vector &y) const override;
 
    /// Compute Mult without applying signs based on DOF orientations.
    void MultUnsigned(const Vector &x, Vector &y) const;
    /// Compute MultTranspose without applying signs based on DOF orientations.
    void MultTransposeUnsigned(const Vector &x, Vector &y) const;
 
-   /// Compute MultTranspose by setting (rather than adding) element
-   /// contributions; this is a left inverse of the Mult() operation
-   void MultLeftInverse(const Vector &x, Vector &y) const;
 
    /// @brief Fills the E-vector y with `boolean` values 0.0 and 1.0 such that each
    /// each entry of the L-vector is uniquely represented in `y`.
@@ -121,6 +123,7 @@ public:
    void MultTranspose(const Vector &x, Vector &y) const override;
    void AddMultTranspose(const Vector &x, Vector &y,
                          const real_t a = 1.0) const override;
+   void MultLeftInverse(const Vector &x, Vector &y) const override;
    /** Fill the I array of SparseMatrix corresponding to the sparsity pattern
        given by this ElementRestriction. */
    void FillI(SparseMatrix &mat) const;

tests/unit/CMakeLists.txt

@@ -71,6 +71,7 @@ set(UNIT_TESTS_SRCS
   fem/test_calcdshape.cpp
   fem/test_calcshape.cpp
   fem/test_calcvshape.cpp
+  fem/test_change_basis.cpp
   fem/test_coefficient.cpp
   fem/test_datacollection.cpp
   fem/test_derefine.cpp

tests/unit/fem/test_change_basis.cpp

@@ -0,0 +1,102 @@
+// Copyright (c) 2010-2024, 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.
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+using namespace mfem;
+
+TEST_CASE("Change of Basis", "[ChangeOfBasis][CUDA]")
+{
+   const auto mesh_fname = GENERATE(
+                              "../../data/star.mesh",
+                              "../../data/star-q3.mesh",
+                              "../../data/fichera.mesh",
+                              "../../data/fichera-q3.mesh"
+                           );
+   const int order = GENERATE(1, 2, 3);
+   const int btype1 = GENERATE(BasisType::GaussLobatto,
+                               BasisType::ClosedUniform,
+                               BasisType::Positive);
+   const int btype2 = GENERATE(BasisType::GaussLobatto,
+                               BasisType::ClosedUniform,
+                               BasisType::Positive);
+
+   if (btype1 == btype2) { return; }
+
+   CAPTURE(mesh_fname, order, btype1, btype2);
+
+   Mesh mesh = Mesh::LoadFromFile(mesh_fname);
+   H1_FECollection fec1(order, mesh.Dimension(), btype1);
+   H1_FECollection fec2(order, mesh.Dimension(), btype2);
+   FiniteElementSpace fes1(&mesh, &fec1);
+   FiniteElementSpace fes2(&mesh, &fec2);
+
+   DiscreteLinearOperator op1(&fes1, &fes2);
+   op1.AddDomainInterpolator(new IdentityInterpolator);
+   op1.Assemble();
+
+   ChangeOfBasis op2(fes1, btype2);
+
+   GridFunction x1(&fes1), x2(&fes1), y1(&fes2), y2(&fes2);
+   x1.Randomize(1);
+
+   op1.Mult(x1, y1);
+   op2.Mult(x1, y2);
+   op2.MultInverse(y2, x2);
+
+   y2 -= y1;
+   x2 -= x1;
+
+   REQUIRE(y2.Normlinf() == MFEM_Approx(0.0));
+   REQUIRE(x2.Normlinf() == MFEM_Approx(0.0));
+}
+
+TEST_CASE("Change of Basis Legendre", "[ChangeOfBasis][CUDA]")
+{
+   const auto mesh_fname = GENERATE(
+                              "../../data/star.mesh",
+                              "../../data/star-q3.mesh",
+                              "../../data/fichera.mesh",
+                              "../../data/fichera-q3.mesh"
+                           );
+   const int order = GENERATE(1, 2, 3);
+   const int btype = GENERATE(BasisType::GaussLobatto,
+                              BasisType::ClosedUniform,
+                              BasisType::Positive);
+   const int dest_btype = GENERATE(ChangeOfBasis::LEGENDRE,
+                                   ChangeOfBasis::INTEGRATED_LEGENDRE);
+
+   CAPTURE(mesh_fname, order, btype);
+
+   Mesh mesh = Mesh::LoadFromFile(mesh_fname);
+   std::unique_ptr<FiniteElementCollection> fec;
+   if (dest_btype == ChangeOfBasis::LEGENDRE)
+   {
+      fec.reset(new L2_FECollection(order, mesh.Dimension(), btype));
+   }
+   else
+   {
+      fec.reset(new H1_FECollection(order, mesh.Dimension(), btype));
+   }
+   FiniteElementSpace fes(&mesh, fec.get());
+
+   ChangeOfBasis op(fes, dest_btype);
+
+   GridFunction x1(&fes), x2(&fes), y(&fes);
+   x1.Randomize(1);
+   op.Mult(x1, y);
+   op.MultInverse(y, x2);
+
+   x2 -= x1;
+
+   REQUIRE(x2.Normlinf() == MFEM_Approx(0.0));
+}