.gitignore

 build
+buildDebug
+builds
 .ccls-cache/
 a.out
 .ccls-cache/
 .cache
+.*
+!.gitignore
+!.gitlab-ci.yml
+!.gitattributes
+!.clang-format
+!.dockerignore
+!.gitlab-ci

CMakeLists.txt

@@ -22,7 +22,7 @@ message(STATUS "Found Nektar++: version ${NEKTAR++_VERSION}")
 
 set(CMAKE_INSTALL_RPATH "${NEKTAR++_LIBRARY_DIRS}")
 
-set(SRC Operators/Operator.cpp Operators/OperatorBwdTrans.cpp Operators/BwdTrans/BwdTransImpl.cpp)
+set(SRC Field.cpp Operators/Operator.cpp Operators/OperatorBwdTrans.cpp Operators/BwdTrans/BwdTransImpl.cpp Operators/OperatorIProductWRTBase.cpp Operators/IProductWRTBase/IProductWRTBaseImpl.cpp)
 
 if (NEKTAR_USE_CUDA)
     enable_language(CUDA)

Field.cpp

+#include "Field.hpp"
+#include <MultiRegions/ExpList.h>
+
+using namespace Nektar;
+using namespace LibUtilities;
+
+std::vector<BlockAttributes> GetBlockAttributes(
+    FieldState state, const MultiRegions::ExpListSharedPtr explist)
+{
+    std::vector<BlockAttributes> blockAttr;
+
+    // initialize the basisKeys
+    std::vector<BasisKey> prevbasisKeys(3, NullBasisKey);
+    std::vector<BasisKey> thisbasisKeys(3, NullBasisKey);
+    int prevIsDeformed = -1, thisIsDeformed = -1;
+
+    // loop over elements
+    for (int i = 0; i < explist->GetNumElmts(); i++)
+    {
+        auto expPtr = explist->GetExp(i);
+
+        // fetch basiskeys of current element
+        for (int d = 0; d < expPtr->GetNumBases(); d++)
+        {
+            thisbasisKeys[d] = expPtr->GetBasis(d)->GetBasisKey();
+        }
+
+        thisIsDeformed = expPtr->GetMetricInfo()->GetGtype();
+
+        // if the basis is the same as the previous one,
+        // increment the number of elements
+        if (thisbasisKeys == prevbasisKeys && thisIsDeformed == prevIsDeformed)
+        {
+            blockAttr.back().num_elements++;
+            blockAttr.back().block_size += blockAttr.back().num_pts;
+        }
+        else // if not, create a new block with the number of elements = 1
+        {
+            size_t num_pts = state == FieldState::Phys ? expPtr->GetTotPoints()
+                                                       : expPtr->GetNcoeffs();
+            blockAttr.push_back({1, num_pts});
+            prevbasisKeys  = thisbasisKeys;
+            prevIsDeformed = thisIsDeformed;
+        }
+    }
+
+    return blockAttr;
+}

Field.hpp

 #pragma once
 
-#include <StdRegions/StdExpansion.h>
-#include <StdRegions/StdRegions.hpp>
 #include <array>
 #include <iostream>
 #include <memory>
@@ -14,7 +12,15 @@
 #include <LibUtilities/BasicUtils/ErrorUtil.hpp>
 #include <LibUtilities/BasicUtils/ShapeType.hpp>
 #include <LibUtilities/BasicUtils/SharedArray.hpp>
-#include <StdRegions/StdExpansion.h>
+
+namespace Nektar
+{
+namespace MultiRegions
+{
+class ExpList;
+typedef std::shared_ptr<ExpList> ExpListSharedPtr;
+} // namespace MultiRegions
+} // namespace Nektar
 
 /**
  * @brief Captures the structure of a block of elements of identical shape
@@ -49,6 +55,9 @@ enum class FieldState
 
 static constexpr FieldState DefaultState = FieldState::Phys;
 
+std::vector<BlockAttributes> GetBlockAttributes(
+    FieldState state, const Nektar::MultiRegions::ExpListSharedPtr explist);
+
 /**
  * @brief A Field represents expansion data to be operated on.
  * @tparam TType  The floating-point representation used by the field.
@@ -56,8 +65,6 @@ static constexpr FieldState DefaultState = FieldState::Phys;
  */
 template <typename TType = double, FieldState TState = DefaultState> class Field
 {
-    using ExpansionSharedPtr = Nektar::StdRegions::StdExpansionSharedPtr;
-
 public:
     Field(const Field &) = delete;
     virtual ~Field()     = default;

MemoryRegionCPU.hpp

@@ -13,8 +13,8 @@
 template <typename TData> class MemoryRegionCPU
 {
 public:
-    MemoryRegionCPU()                                      = delete;
-    MemoryRegionCPU(const MemoryRegionCPU &rhs)            = delete;
+    MemoryRegionCPU()                           = delete;
+    MemoryRegionCPU(const MemoryRegionCPU &rhs) = delete;
     MemoryRegionCPU &operator=(const MemoryRegionCPU &rhs) = delete;
 
     MemoryRegionCPU(MemoryRegionCPU &&rhs)

MemoryRegionCUDA.hpp

@@ -42,10 +42,10 @@ public:
     void operator=(MemoryRegionCUDA &&rhs)
     {
         MemoryRegionCPU<TData>::operator=(std::move(rhs));
-        m_device     = rhs.m_device;
-        m_size       = rhs.m_size;
-        rhs.m_device = nullptr;
-        rhs.m_size   = 0;
+        m_device                        = rhs.m_device;
+        m_size                          = rhs.m_size;
+        rhs.m_device                    = nullptr;
+        rhs.m_size                      = 0;
     }
 
     virtual TData *GetCPUPtr() override

Operators/BwdTrans/BwdTransImpl.cpp

-#include "BwdTransStdMat.hpp"
 #include "BwdTransMatFree.hpp"
+#include "BwdTransStdMat.hpp"
 #include "BwdTransSumFac.hpp"
 
 namespace Nektar::Operators::detail

Operators/BwdTrans/BwdTransStdMat.hpp

@@ -20,10 +20,11 @@ public:
         auto const *inptr = in.GetStorage().GetCPUPtr();
         auto *outptr      = out.GetStorage().GetCPUPtr();
 
+        size_t exp_idx = 0;
         for (size_t block_idx = 0; block_idx < in.GetBlocks().size();
              ++block_idx)
         {
-            auto const expPtr = this->m_expansionList->GetExp(block_idx);
+            auto const expPtr = this->m_expansionList->GetExp(exp_idx);
 
             Nektar::StdRegions::StdMatrixKey key(
                 StdRegions::eBwdTrans, expPtr->DetShapeType(), *expPtr);
@@ -38,6 +39,7 @@ public:
 
             inptr += block.block_size;
             outptr += expPtr->GetTotPoints() * block.num_elements;
+            exp_idx += block.num_elements;
         }
     }
 

Operators/IProductWRTBase/IProductWRTBaseImpl.cpp

+#include "IProductWRTBaseStdMat.hpp"
+
+namespace Nektar::Operators::detail
+{
+
+// Add different IProductWRTBase implementations to the factory.
+template <>
+std::string OperatorIProductWRTBaseImpl<double, ImplStdMat>::className =
+    GetOperatorFactory<double>().RegisterCreatorFunction(
+        "IProductWRTBaseStdMat",
+        OperatorIProductWRTBaseImpl<double, ImplStdMat>::instantiate, "...");
+
+} // namespace Nektar::Operators::detail

Operators/IProductWRTBase/IProductWRTBaseStdMat.hpp

+#include "Operators/OperatorIProductWRTBase.hpp"
+#include <StdRegions/StdExpansion.h>
+
+namespace Nektar::Operators::detail
+{
+
+// standard matrix implementation
+template <typename TData>
+class OperatorIProductWRTBaseImpl<TData, ImplStdMat>
+    : public OperatorIProductWRTBase<TData>
+{
+public:
+    OperatorIProductWRTBaseImpl(
+        const MultiRegions::ExpListSharedPtr &expansionList)
+        : OperatorIProductWRTBase<TData>(std::move(expansionList))
+    {
+        size_t jacSize   = 0;
+        size_t nTotElmts = this->m_expansionList->GetNumElmts();
+        // Calculate the jacobian array size
+        for (size_t e = 0; e < nTotElmts; ++e)
+        {
+            // Determine shape and type of the element
+            auto const expPtr = this->m_expansionList->GetExp(e);
+            if (expPtr->GetMetricInfo()->GetGtype() ==
+                SpatialDomains::eDeformed)
+            {
+                jacSize += expPtr->GetTotPoints();
+            }
+            else
+            {
+                jacSize++;
+            }
+        }
+
+        // Allocate memory for the jacobian
+        m_jac = {jacSize, 0.0};
+
+        // Initialise jacobian.
+        size_t index = 0;
+        for (size_t e = 0; e < nTotElmts; ++e)
+        {
+            auto expPtr = this->m_expansionList->GetExp(e);
+            auto &auxJac =
+                expPtr->GetMetricInfo()->GetJac(expPtr->GetPointsKeys());
+            if (expPtr->GetMetricInfo()->GetGtype() ==
+                SpatialDomains::eDeformed)
+            {
+                size_t nqe = expPtr->GetTotPoints();
+                for (size_t i = 0; i < nqe; ++i)
+                {
+                    m_jac[index++] = auxJac[i];
+                }
+            }
+            else
+            {
+                m_jac[index++] = auxJac[0];
+            }
+        }
+    }
+
+    void apply(Field<TData, FieldState::Phys> &in,
+               Field<TData, FieldState::Coeff> &out) override
+    {
+        auto const *inptr = in.GetStorage().GetCPUPtr();
+        auto *outptr      = out.GetStorage().GetCPUPtr();
+
+        size_t expIdx = 0;
+        size_t jacIdx = 0;
+        for (size_t block_idx = 0; block_idx < in.GetBlocks().size();
+             ++block_idx)
+        {
+            auto const expPtr = this->m_expansionList->GetExp(expIdx);
+            auto nElmts       = in.GetBlocks()[block_idx].num_elements;
+            auto nqTot        = expPtr->GetTotPoints();
+            auto nmTot        = expPtr->GetNcoeffs();
+
+            Nektar::StdRegions::StdMatrixKey key(
+                StdRegions::eIProductWRTBase, expPtr->DetShapeType(), *expPtr);
+
+            // This is the B^{T} matrix
+            auto const matPtr = expPtr->GetStdMatrix(key);
+
+            Array<OneD, NekDouble> wsp(nqTot * nElmts, 0.0);
+            if (expPtr->GetMetricInfo()->GetGtype() ==
+                SpatialDomains::eDeformed)
+            {
+                for (size_t i = 0; i < nElmts * nqTot; ++i)
+                {
+                    wsp[i] = m_jac[jacIdx++] * inptr[i];
+                }
+            }
+            else
+            {
+                for (size_t e = 0; e < nElmts; ++e)
+                {
+                    for (size_t i = 0; i < nqTot; ++i)
+                    {
+                        wsp[e * nqTot + i] =
+                            m_jac[jacIdx] * inptr[e * nqTot + i];
+                    }
+                    jacIdx++;
+                }
+            }
+
+            Blas::Dgemm('N', 'N', matPtr->GetRows(), nElmts,
+                        matPtr->GetColumns(), 1.0, matPtr->GetRawPtr(),
+                        matPtr->GetRows(), wsp.get(), nqTot, 0.0, outptr,
+                        nmTot);
+
+            inptr += nqTot * nElmts;
+            outptr += nmTot * nElmts;
+            expIdx += nElmts;
+        }
+    }
+
+    static std::unique_ptr<Operator<TData>> instantiate(
+        const MultiRegions::ExpListSharedPtr &expansionList)
+    {
+        return std::make_unique<OperatorIProductWRTBaseImpl<TData, ImplStdMat>>(
+            std::move(expansionList));
+    }
+
+    static std::string className;
+
+private:
+    Array<OneD, NekDouble> m_jac;
+};
+
+} // namespace Nektar::Operators::detail

Operators/Operator.hpp

 #pragma once
 
+#include "Field.hpp"
+
 #include <string>
 
 #include <LibUtilities/BasicUtils/NekFactory.hpp>

Operators/OperatorIProductWRTBase.cpp

+#include "OperatorIProductWRTBase.hpp"
+
+namespace Nektar::Operators
+{
+
+template <> const std::string IProductWRTBase<>::key = "IProductWRTBase";
+template <> const std::string IProductWRTBase<>::default_impl = "StdMat";
+
+} // namespace Nektar::Operators
\ No newline at end of file

Operators/OperatorIProductWRTBase.hpp

+#pragma once
+
+#include "Field.hpp"
+#include "Operator.hpp"
+
+namespace Nektar::Operators
+{
+
+// IProductWRTBase base class
+// Defines the apply operator to enforce apply parameter types
+template <typename TData> class OperatorIProductWRTBase : public Operator<TData>
+{
+public:
+    virtual ~OperatorIProductWRTBase() = default;
+
+    OperatorIProductWRTBase(const MultiRegions::ExpListSharedPtr &expansionList)
+        : Operator<TData>(std::move(expansionList))
+    {
+    }
+
+    virtual void apply(Field<TData, FieldState::Phys> &in,
+                       Field<TData, FieldState::Coeff> &out) = 0;
+    virtual void operator()(Field<TData, FieldState::Phys> &in,
+                            Field<TData, FieldState::Coeff> &out)
+    {
+        apply(in, out);
+    }
+};
+
+// Descriptor / traits class for IProductWRTBase
+template <typename TData = default_fp_type> struct IProductWRTBase
+{
+    using class_name = OperatorIProductWRTBase<TData>;
+    using FieldIn    = Field<TData, FieldState::Phys>;
+    using FieldOut   = Field<TData, FieldState::Coeff>;
+    static const std::string key;
+    static const std::string default_impl;
+
+    static std::shared_ptr<class_name> create(
+        const MultiRegions::ExpListSharedPtr &expansionList,
+        std::string pKey = "")
+    {
+        return Operator<TData>::template create<IProductWRTBase>(
+            std::move(expansionList), pKey);
+    }
+};
+
+namespace detail
+{
+// Template for IProductWRTBase implementations
+template <typename TData, typename Op> class OperatorIProductWRTBaseImpl;
+} // namespace detail
+
+} // namespace Nektar::Operators

line.xml

+<?xml version="1.0" encoding="utf-8" ?>
+<NEKTAR>
+    <GEOMETRY DIM="1" SPACE="1">
+        <VERTEX COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJxjYMAPGFF4H+zR5QEVbwEx</VERTEX>
+        <ELEMENT>
+            <S COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJxjYEAFjFAaAAAgAAIA</S>
+        </ELEMENT>
+        <COMPOSITE>
+            <C ID="0"> S[0] </C>
+        </COMPOSITE>
+        <DOMAIN>
+            <D ID="0"> C[0] </D>
+        </DOMAIN>
+    </GEOMETRY>
+    <Metadata>
+        <Provenance>
+            <GitBranch>refs/heads/master</GitBranch>
+            <GitSHA1>b9041a54043f57d230bcc1eb5be5c9a048b69a92</GitSHA1>
+            <Hostname>dyn3237-113.wlan.ic.ac.uk</Hostname>
+            <NektarVersion>5.3.0</NektarVersion>
+            <Timestamp>07-Jun-2023 11:25:53</Timestamp>
+        </Provenance>
+        <NekMeshCommandLine>-v line.msh line.xml </NekMeshCommandLine>
+    </Metadata>
+    <EXPANSIONS>
+        <E COMPOSITE="C[0]" NUMMODES="6" TYPE="MODIFIED" FIELDS="u" />
+    </EXPANSIONS>
+</NEKTAR>

main.cpp

@@ -14,11 +14,12 @@
 
 #include "Field.hpp"
 #include "Operators/OperatorBwdTrans.hpp"
+#include "Operators/OperatorIProductWRTBase.hpp"
 
 #include <LibUtilities/BasicUtils/SessionReader.h>
 //#include <LibUtilities/SimdLib/tinysimd.hpp>
-#include <SpatialDomains/MeshGraph.h>
 #include <MultiRegions/ExpList.h>
+#include <SpatialDomains/MeshGraph.h>
 
 #ifdef NEKTAR_USE_CUDA
 #include "MemoryRegionCUDA.hpp"
@@ -41,46 +42,26 @@ std::ostream &operator<<(std::ostream &stream, Field<TType, State> &f)
     return stream;
 }
 
-std::vector<BlockAttributes> GetBlockAttributes(
-        FieldState state,
-        const MultiRegions::ExpListSharedPtr explist)
-{
-    const int n = explist->GetNumElmts();
-    std::map<std::tuple<LibUtilities::ShapeType,unsigned int,unsigned int>,size_t> blockList;
-    for (int i = 0; i < explist->GetNumElmts(); ++i)
-    {
-        auto e = explist->GetExp(i);
-        blockList[{e->DetShapeType(),e->GetNcoeffs(),e->GetTotPoints()}]++;
-    }
-    std::vector<BlockAttributes> blockAttr;
-    for (auto &x : blockList)
-    {
-        auto val = state == FieldState::Phys ? std::get<2>(x.first) : std::get<1>(x.first);
-        blockAttr.push_back( { x.second, val } );
-    }
-    return blockAttr;
-}
-
 int main(int argc, char *argv[])
 {
     // Initialise a session, graph and create an expansion list
     LibUtilities::SessionReaderSharedPtr session;
-    SpatialDomains::MeshGraphSharedPtr   graph;
-    MultiRegions::ExpListSharedPtr       explist;
+    SpatialDomains::MeshGraphSharedPtr graph;
+    MultiRegions::ExpListSharedPtr explist;
 
     session = LibUtilities::SessionReader::CreateInstance(argc, argv);
     graph   = SpatialDomains::MeshGraph::Read(session);
-    explist = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr
-                    (session, graph);
+    explist =
+        MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(session, graph);
 
     // Generate a blocks definition from the expansion list for each state
-    auto blocks_phys  = GetBlockAttributes(FieldState::Phys,  explist);
+    auto blocks_phys  = GetBlockAttributes(FieldState::Phys, explist);
     auto blocks_coeff = GetBlockAttributes(FieldState::Coeff, explist);
 
     // Create two Field objects with a MemoryRegionCPU backend by default
     auto in  = Field<double, FieldState::Coeff>::create(blocks_coeff);
     auto in2 = Field<double, FieldState::Coeff>::create(blocks_coeff);
-    auto out = Field<double, FieldState::Phys >::create(blocks_phys);
+    auto out = Field<double, FieldState::Phys>::create(blocks_phys);
 
     // Populate the field with some data. In this case, we just grab a pointer
     // to the memory on the CPU and populate the array with values. Operators,
@@ -106,15 +87,16 @@ int main(int argc, char *argv[])
         }
     }
 
-    memset(out.GetStorage().GetCPUPtr(), 0, out.GetStorage().size()*sizeof(double));
+    memset(out.GetStorage().GetCPUPtr(), 0,
+           out.GetStorage().size() * sizeof(double));
 
     std::cout << "Initial shape:\n" << in << std::endl << std::endl;
 
+#ifdef NEKTAR_ENABLE_SIMD_AVX2
     // Test out field reshaping
     in.ReshapeStorage<4>();
     std::cout << "Reshaped to 4:\n" << in << std::endl << std::endl;
 
-#ifdef NEKTAR_ENABLE_SIMD_AVX2
     // Test out SIMD instructions
     using vec_t = tinysimd::simd<double>;
 
@@ -195,14 +177,15 @@ int main(int argc, char *argv[])
 
     std::cout << out << std::endl;
 
-
 #ifdef NEKTAR_USE_CUDA
 
     // Test CUDA MemoryRegion
 
     // Create two Fields with memory on the GPU
-    in  = Field<double, FieldState::Coeff>::create<MemoryRegionCUDA>(blocks_coeff);
-    out = Field<double, FieldState::Phys>::create<MemoryRegionCUDA>(blocks_phys);
+    in = Field<double, FieldState::Coeff>::create<MemoryRegionCUDA>(
+        blocks_coeff);
+    out =
+        Field<double, FieldState::Phys>::create<MemoryRegionCUDA>(blocks_phys);
 
     // Perform the BwdTrans on the fields using the CUDA implementation
     // Since this is a CUDA operator, acting on CUDA fields, everything happens
@@ -224,6 +207,40 @@ int main(int argc, char *argv[])
     BwdTrans<>::create(explist, "CUDA")->apply(in, out);
 #endif
 
+    std::cout << "Inner Product of a function WRT the Basis" << std::endl;
+
+    // Create two Field objects with a MemoryRegionCPU backend by default
+    // for the inner product with respect to base
+    auto inPhys   = Field<double, FieldState::Phys>::create(blocks_phys);
+    auto outCoeff = Field<double, FieldState::Coeff>::create(blocks_coeff);
+
+    double *y = inPhys.GetStorage().GetCPUPtr();
+    for (auto const &block : blocks_phys)
+    {
+        for (size_t el = 0; el < block.num_elements; ++el)
+        {
+            for (size_t phys = 0; phys < block.num_pts; ++phys)
+            {
+                // Each element is the index of the quadrature points
+                // this is useful for testing reshapes
+                *(y++) = phys;
+            }
+        }
+    }
+
+    memset(outCoeff.GetStorage().GetCPUPtr(), 0,
+           outCoeff.GetStorage().size() * sizeof(double));
+
+    std::cout << "Initial shape:\n" << inPhys << std::endl << std::endl;
+
+    // IProductWRTBase
+    auto ipwrtb = IProductWRTBase<>::create(explist, "StdMat");
+    // ... and then apply it
+    ipwrtb->apply(inPhys, outCoeff); // out and in is inverted for the IP
+
+    // Let's display the result
+    std::cout << "Out:\n" << outCoeff << std::endl;
+
     std::cout << "END" << std::endl;
     return 0;
 }

square2.xml

+<?xml version="1.0" encoding="utf-8" ?>
+<NEKTAR>
+    <GEOMETRY DIM="2" SPACE="2">
+        <VERTEX COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJx1j8sVgjAABIOIgIrK324o3RYoIaV4yi7guBfey2SySwjbxM/+G0K24+sPPyFfxXN83/yMvlOg71y4f0knJfeLV+hH8Rp98+tx0GHnDbn/847cadhX/4P7xZ/oed8Lufe1f3pTOvZ1r+d+8QF3e9+I3Psm5M7Mvu69uV/8C2/eNsUA</VERTEX>
+        <EDGE COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJx1kEkSglAMBQEVZxxBHEARh/vf0I29oKvyN6l+FUI6STJ8aVCzoI6CyhurH54EnIuZN1U/fTMxfXMx8xbB90sx36+Ceet/5T7MLcTM3Yi5x1Y5/9uJ8d4rZ4+DGL+jcvYrxXhX2p+9T2L2rgOPs3J8LmI8rsrxvInxaJTj34rxu2t/7vIQc5dOjMdTOffqxXi8lHPHtxiPj3Lu+xXj9wOocAa/</EDGE>
+        <ELEMENT>
+            <Q COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJx10EcOwjAABdFQ0kkPvYTO/W/IZrxgJLx50kSyvxJFv2eGc1z86UuMMdH30FPMMNe9oRdY4krvVFhjo12xeoud7knUexy0N1Ufca29Yf8Gt7jT3lx9jwftLdSPeNLeUv2MF+0N/3HCK960t1K/40N7a/UnvrS3UX/jR3u/iOQFQAAA</Q>
+        </ELEMENT>
+        <COMPOSITE>
+            <C ID="0"> Q[0-11] </C>
+            <C ID="1"> E[0,13,22,31,11,21,30,39,3,6,9,12,32,34,36,38] </C>
+            <C ID="2"> Q[12-15] </C>
+        </COMPOSITE>
+        <DOMAIN> C[0,2] </DOMAIN>
+    </GEOMETRY>
+    <Metadata>
+        <Provenance>
+            <GitBranch></GitBranch>
+            <GitSHA1>88b1cb3be29cc9090d8a14af19c26decd88345db</GitSHA1>
+            <Hostname>kingscross</Hostname>
+            <NektarVersion>5.1.0</NektarVersion>
+            <Timestamp>27-Apr-2023 15:32:08</Timestamp>
+        </Provenance>
+        <NekMeshCommandLine>square.msh square.xml </NekMeshCommandLine>
+    </Metadata>
+    <EXPANSIONS>
+        <E COMPOSITE="C[0]" NUMMODES="4" TYPE="MODIFIED" FIELDS="u" />
+        <E COMPOSITE="C[2]" NUMMODES="3" TYPE="MODIFIED" FIELDS="u" />
+    </EXPANSIONS>
+    <CONDITIONS>
+        <SOLVERINFO>
+            <I PROPERTY="EQTYPE" VALUE="Helmholtz" />
+        </SOLVERINFO>
+        <VARIABLES>
+            <V ID="0">u</V>
+        </VARIABLES>
+        <FUNCTION NAME="Forcing">
+            <E VAR="u" VALUE="0" />
+        </FUNCTION>
+    </CONDITIONS>
+</NEKTAR>

squareDeformed.xml

+<?xml version="1.0" encoding="utf-8" ?>
+<NEKTAR>
+    <GEOMETRY DIM="2" SPACE="2">
+        <VERTEX COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJxjYICAx7ZT/D4/OW/vfW1rRBvvBXuoMAMjAyoo9QCRCHkmBvyAGUr/XPwfCM7bo8uzoPAeYMizosmjm8OGVT/CHHas9iPkObC6GiHPiSL+AcN9XGjy6O7jhtJ6YPKFfcsjVHkeFP0vMMznxSqPcB8fVvch5PlR5GHiCHsEoPQntRlz3sx8aS+W9+p69IyXcHlBKK0QYgwEr+195s0Ego9weSGs5iPCSRjmsikQfZ1CEHNg8iJQ2hgMPtvDaJi8KMzn4IRywd4Dzf1iDNgAQl4cSqOmW4T7JLDqR8gDAHwmZuEA</VERTEX>
+        <EDGE COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJx1kkcSwlAMQ0MLoQZC7z3c/4Zs/Ba8GXuj0ZflyJ4UxX91Euwm2EuQ6usdPkh4KU4Nk75KnL6RODXWO/6JOHtP5WPuTD70uXzotfzoC3HusQxkX+Y0yfdXgaX0tXKRc6M+9K3ykm8nTr69OPkO4sw7Ki/zTuLkOisfvovmsvc1kL3JdZMf/R7I/8BdHuLc5SnO//IK5E7c5S2O76P5+NrAWvo3sFH/D7Z4BrgA</EDGE>
+        <ELEMENT>
+            <Q COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJx10ccSwjAMRdEAofeS0CH0//9DNvcueDNkc8a2LEtKUfx+LWxjJ/Zdl9jFXtxz3ccBDiPPCMdRh3G+M8EpznAedSzifPmnTvdXuI56dYNVxFuHfdW4xR3u0X4OeIy8J3Qe5jtHfd5zXpd41/xXdG7mbfCGd7SfBz4j7oXO1b6t7x35/S+fuG+cc/kCmoQFPwAA</Q>
+        </ELEMENT>
+        <COMPOSITE>
+            <C ID="13"> E[36,34,11,1] </C>
+            <C ID="14"> E[2,4,16,15] </C>
+            <C ID="15"> E[12,20,28,31] </C>
+            <C ID="16"> E[30,24,39,35] </C>
+            <C ID="17"> Q[0-15] </C>
+        </COMPOSITE>
+        <DOMAIN>
+            <D ID="0"> C[17] </D>
+        </DOMAIN>
+    </GEOMETRY>
+    <EXPANSIONS>
+        <E COMPOSITE="C[17]" NUMMODES="4" TYPE="MODIFIED" FIELDS="u" />
+    </EXPANSIONS>
+</NEKTAR>

tests/CMakeLists.txt

-set(test_src test_bwdtrans.cpp)
-
 find_package(
   Boost
   REQUIRED
   COMPONENTS unit_test_framework
 )
 
-add_executable(test_bwdtrans ${test_src})
+add_executable(test_bwdtrans test_bwdtrans.cpp)
 target_link_libraries(test_bwdtrans PRIVATE Operators)
 target_link_libraries(test_bwdtrans PRIVATE ${NEKTAR++_LIBRARIES})
 target_link_libraries(test_bwdtrans PRIVATE Boost::unit_test_framework)
-
 target_include_directories(test_bwdtrans PRIVATE "${CMAKE_SOURCE_DIR}")
 target_include_directories(test_bwdtrans PRIVATE ${NEKTAR++_INCLUDE_DIRS})
+target_compile_definitions(test_bwdtrans PRIVATE
+    -DTEST_PATH="${CMAKE_SOURCE_DIR}")
+
+add_executable(test_ipwrtbase test_ipwrtbase.cpp)
+target_link_libraries(test_ipwrtbase PRIVATE Operators)
+target_link_libraries(test_ipwrtbase PRIVATE ${NEKTAR++_LIBRARIES})
+target_link_libraries(test_ipwrtbase PRIVATE Boost::unit_test_framework)
+target_include_directories(test_ipwrtbase PRIVATE "${CMAKE_SOURCE_DIR}")
+target_include_directories(test_ipwrtbase PRIVATE ${NEKTAR++_INCLUDE_DIRS})
+target_compile_definitions(test_ipwrtbase PRIVATE
+    -DTEST_PATH="${CMAKE_SOURCE_DIR}")
 
 add_test(
   NAME BwdTrans
   COMMAND test_bwdtrans
   # Test executable is looking for input session file in the working
   # directory
-  WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/input"
+  WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}"
+)
+
+add_test(
+  NAME IProductWRTBase
+  COMMAND test_ipwrtbase
+  WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}"
 )

tests/init_fields.hpp

+#pragma once
 #include <boost/test/unit_test_log.hpp>
-#include <vector>
 #include <string>
+#include <vector>
 
 #include "Field.hpp"
 #include <MultiRegions/ExpList.h>
 #include <Operators/OperatorBwdTrans.hpp>
+#include <Operators/OperatorIProductWRTBase.hpp>
 
 using namespace Nektar::Operators;
 using namespace Nektar::LibUtilities;
@@ -18,62 +20,57 @@ using namespace Nektar;
  * before (after) each call to BOOST_FIXTURE_TEST_CASE(<test name>,
  * InitFields) macro.
  *
- * See "Single test case fixture" on the Boost.Test documentation for more details:
+ * See "Single test case fixture" on the Boost.Test documentation for more
+ * details:
  * https://www.boost.org/doc/libs/1_82_0/libs/test/doc/html/boost_test/tests_organization/fixtures/case.html
  */
-struct InitFields {
-  Field<double, FieldState::Coeff> *fixt_in;
-  Field<double, FieldState::Phys> *fixt_out;
-  MultiRegions::ExpListSharedPtr fixt_explist {nullptr};
-  ~InitFields() {BOOST_TEST_MESSAGE("teardown fixture");}
 
-
-  static std::vector<BlockAttributes>
-  GetBlockAttributes(FieldState state,
-		     const MultiRegions::ExpListSharedPtr explist) {
-        const int n = explist->GetNumElmts();
-    std::map<std::tuple<LibUtilities::ShapeType,unsigned int,unsigned int>,size_t> blockList;
-    for (int i = 0; i < explist->GetNumElmts(); ++i)
+template <FieldState stateIn  = FieldState::Coeff,
+          FieldState stateOut = FieldState::Phys>
+class InitFields
+{
+public:
+    Field<double, stateIn> *fixt_in;
+    Field<double, stateOut> *fixt_out;
+    MultiRegions::ExpListSharedPtr fixt_explist{nullptr};
+    ~InitFields()
     {
-        auto e = explist->GetExp(i);
-        blockList[{e->DetShapeType(),e->GetNcoeffs(),e->GetTotPoints()}]++;
+        BOOST_TEST_MESSAGE("teardown fixture");
     }
-    std::vector<BlockAttributes> blockAttr;
-    for (auto &x : blockList)
+
+    InitFields()
     {
-        auto val = state == FieldState::Phys ? std::get<2>(x.first) : std::get<1>(x.first);
-        blockAttr.push_back( { x.second, val } );
     }
-    return blockAttr;
-  }
 
-  InitFields() {
-    BOOST_TEST_MESSAGE("Creating input and output fields");
-    // Initialise a session, graph and create an expansion list
-    LibUtilities::SessionReaderSharedPtr session;
-    SpatialDomains::MeshGraphSharedPtr   graph;
+    void Configure()
+    {
+        BOOST_TEST_MESSAGE("Creating input and output fields");
+        // Initialise a session, graph and create an expansion list
+        LibUtilities::SessionReaderSharedPtr session;
+        SpatialDomains::MeshGraphSharedPtr graph;
+
+        // Construct a fake command-line argument array to be fed to
+        // Session::Reader::CreateInstance. The first element stands for
+        // the name of the executable which, in our case, doesn't matter.
+        int argc     = 2;
+        char *argv[] = {(char *)"exe_name", GetMeshName().data()};
 
-    // Construct a fake command-line argument array to be fed to
-    // Session::Reader::CreateInstance. The first element stands for
-    // the name of the executable which, in our case, doesn't matter.
-    int argc = 2;
-    char *argv[] = {
-      (char*)"exe_name", (char*)"square.xml"
-    };
+        session      = LibUtilities::SessionReader::CreateInstance(argc, argv);
+        graph        = SpatialDomains::MeshGraph::Read(session);
+        fixt_explist = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
+            session, graph);
 
-    session = LibUtilities::SessionReader::CreateInstance(argc, argv);
-    graph   = SpatialDomains::MeshGraph::Read(session);
-    fixt_explist = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr
-                    (session, graph);
+        // Generate a blocks definition from the expansion list for each state
+        auto blocks_in  = GetBlockAttributes(stateIn, fixt_explist);
+        auto blocks_out = GetBlockAttributes(stateOut, fixt_explist);
 
-    // Generate a blocks definition from the expansion list for each state
-    auto blocks_phys  = InitFields::GetBlockAttributes(FieldState::Phys,  fixt_explist);
-    auto blocks_coeff = InitFields::GetBlockAttributes(FieldState::Coeff, fixt_explist);
+        // Create two Field objects with a MemoryRegionCPU backend by default
+        auto f_in  = Field<double, stateIn>::create(blocks_in);
+        auto f_out = Field<double, stateOut>::create(blocks_out);
+        fixt_in    = new Field<double, stateIn>(std::move(f_in));
+        fixt_out   = new Field<double, stateOut>(std::move(f_out));
+    }
 
-    // Create two Field objects with a MemoryRegionCPU backend by default
-    auto f_in  = Field<double, FieldState::Coeff>::create(blocks_coeff);
-    auto f_out = Field<double, FieldState::Phys >::create(blocks_phys);
-    fixt_in  = new Field<double, FieldState::Coeff>(std::move(f_in));
-    fixt_out = new Field<double, FieldState::Phys>(std::move(f_out));
-  }
+protected:
+    virtual std::string GetMeshName() = 0;
 };

tests/test_bwdtrans.cpp

@@ -8,8 +8,8 @@
 #include "Operators/OperatorBwdTrans.hpp"
 
 #include <LibUtilities/BasicUtils/SessionReader.h>
-#include <SpatialDomains/MeshGraph.h>
 #include <MultiRegions/ExpList.h>
+#include <SpatialDomains/MeshGraph.h>
 
 #include "init_fields.hpp"
 
@@ -18,8 +18,29 @@ using namespace Nektar::Operators;
 using namespace Nektar::LibUtilities;
 using namespace Nektar;
 
-BOOST_FIXTURE_TEST_CASE( bwdtrans, InitFields )
+class Line : public InitFields<FieldState::Coeff, FieldState::Phys>
+{
+protected:
+    virtual std::string GetMeshName() override
+    {
+        return "line.xml";
+    }
+};
+
+class Square : public InitFields<FieldState::Coeff, FieldState::Phys>
 {
+protected:
+    virtual std::string GetMeshName() override
+    {
+        return "square.xml";
+    }
+};
+
+using init = InitFields<FieldState::Coeff, FieldState::Phys>;
+BOOST_FIXTURE_TEST_CASE(bwdtrans, Line)
+{
+    Configure();
+
     double *x = fixt_in->GetStorage().GetCPUPtr();
 
     // For each element, initialise first coefficient to zero and rest
@@ -30,10 +51,12 @@ BOOST_FIXTURE_TEST_CASE( bwdtrans, InitFields )
         {
             for (size_t coeff = 0; coeff < block.num_pts; ++coeff)
             {
-                if (coeff == 0) {
+                if (coeff == 0)
+                {
                     *(x++) = 1.0;
                 }
-                else {
+                else
+                {
                     *(x++) = 0.0;
                 }
             }
@@ -43,5 +66,5 @@ BOOST_FIXTURE_TEST_CASE( bwdtrans, InitFields )
     BwdTrans<>::create(fixt_explist, "StdMat")->apply(*fixt_in, *fixt_out);
 
     double *y = fixt_out->GetStorage().GetCPUPtr();
-    BOOST_TEST( y[0] == 1.0 );
+    BOOST_TEST(y[0] == 1.0);
 }