CMakeLists.txt

@@ -22,7 +22,7 @@ message(STATUS "Found Nektar++: version ${NEKTAR++_VERSION}")
 
 set(CMAKE_INSTALL_RPATH "${NEKTAR++_LIBRARY_DIRS}")
 
-set(SRC Field.cpp 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

@@ -12,31 +12,35 @@ std::vector<BlockAttributes> GetBlockAttributes(
     // 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 exp = explist->GetExp(i);
+        auto expPtr = explist->GetExp(i);
 
         // fetch basiskeys of current element
-        for (int d = 0; d < exp->GetNumBases(); d++)
+        for (int d = 0; d < expPtr->GetNumBases(); d++)
         {
-            thisbasisKeys[d] = exp->GetBasis(d)->GetBasisKey();
+            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)
+        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 ? exp->GetTotPoints()
-                                                   : exp->GetNcoeffs();
+            size_t num_pts = state == FieldState::Phys ? expPtr->GetTotPoints()
+                                                       : expPtr->GetNcoeffs();
             blockAttr.push_back({1, num_pts});
-            prevbasisKeys = thisbasisKeys;
+            prevbasisKeys  = thisbasisKeys;
+            prevIsDeformed = thisIsDeformed;
         }
     }
 

Field.hpp

@@ -17,11 +17,10 @@ namespace Nektar
 {
 namespace MultiRegions
 {
-    class ExpList;
-    typedef std::shared_ptr<ExpList> ExpListSharedPtr;
-}
-}
-
+class ExpList;
+typedef std::shared_ptr<ExpList> ExpListSharedPtr;
+} // namespace MultiRegions
+} // namespace Nektar
 
 /**
  * @brief Captures the structure of a block of elements of identical shape
@@ -116,9 +115,8 @@ public:
 
         size_t storage_size = std::accumulate(
             field.block_attributes.begin(), field.block_attributes.end(), 0,
-            [](size_t acc, const BlockAttributes &block) {
-                return acc + block.block_size;
-            });
+            [](size_t acc, const BlockAttributes &block)
+            { return acc + block.block_size; });
 
         // Create new TMemoryRegion and polymorphically store as MemoryRegionCPU
         field.m_storage = std::make_unique<TMemoryRegion<TType>>(

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/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/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,6 +14,7 @@
 
 #include "Field.hpp"
 #include "Operators/OperatorBwdTrans.hpp"
+#include "Operators/OperatorIProductWRTBase.hpp"
 
 #include <LibUtilities/BasicUtils/SessionReader.h>
 //#include <LibUtilities/SimdLib/tinysimd.hpp>
@@ -91,11 +92,11 @@ int main(int argc, char *argv[])
 
     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>;
 
@@ -206,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;
 }

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,42 +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");}
-
-  InitFields() {
-    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", (char*)"square.xml"
-    };
-
-    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_phys  = GetBlockAttributes(FieldState::Phys,  fixt_explist);
-    auto blocks_coeff = GetBlockAttributes(FieldState::Coeff, fixt_explist);
-
-    // 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));
-  }
+
+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()
+    {
+        BOOST_TEST_MESSAGE("teardown fixture");
+    }
+
+    InitFields()
+    {
+    }
+
+    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()};
+
+        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);
+
+        // 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));
+    }
+
+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);
 }

tests/test_ipwrtbase.cpp

+#define BOOST_TEST_MODULE example
+#include <boost/test/unit_test.hpp>
+
+#include <iostream>
+#include <memory>
+
+#include "Field.hpp"
+#include "Operators/OperatorIProductWRTBase.hpp"
+
+#include <LibUtilities/BasicUtils/SessionReader.h>
+#include <MultiRegions/ExpList.h>
+#include <SpatialDomains/MeshGraph.h>
+
+#include "init_fields.hpp"
+
+using namespace std;
+using namespace Nektar::Operators;
+using namespace Nektar::LibUtilities;
+using namespace Nektar;
+
+class Line : public InitFields<FieldState::Phys, FieldState::Coeff>
+{
+protected:
+    virtual std::string GetMeshName() override
+    {
+        return "line.xml";
+    }
+};
+
+class Square : public InitFields<FieldState::Phys, FieldState::Coeff>
+{
+protected:
+    virtual std::string GetMeshName() override
+    {
+        return "square.xml";
+    }
+};
+
+BOOST_FIXTURE_TEST_CASE(ipwrtbase, Line)
+{
+    Configure();
+
+    static double *x = fixt_in->GetStorage().GetCPUPtr();
+
+    for (auto const &block : fixt_in->GetBlocks())
+    {
+        for (size_t el = 0; el < block.num_elements; ++el)
+        {
+            for (size_t phys = 0; phys < block.num_pts; ++phys)
+            {
+                *(x++) = phys;
+            }
+        }
+    }
+
+    IProductWRTBase<>::create(fixt_explist, "StdMat")
+        ->apply(*fixt_in, *fixt_out);
+
+    static double *y = fixt_out->GetStorage().GetCPUPtr();
+    double TOL       = 1e-12;
+    BOOST_CHECK_CLOSE(y[0], 1.09753217840612, TOL);
+    BOOST_CHECK_CLOSE(y[1], 1.90246782159389, TOL);
+    BOOST_CHECK_CLOSE(y[2], 0.5, TOL);
+    BOOST_CHECK_CLOSE(y[3], 0.150377322580158, TOL);
+    BOOST_TEST(std::abs(y[4] - 9.36750677027476e-17) < TOL);
+    BOOST_CHECK_CLOSE(y[5], 0.00746847423694819, TOL);
+}