.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)
 
 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);
+
+    // loop over elements
+    for (int i = 0; i < explist->GetNumElmts(); i++)
+    {
+        auto exp = explist->GetExp(i);
+
+        // fetch basiskeys of current element
+        for (int d = 0; d < exp->GetNumBases(); d++)
+        {
+            thisbasisKeys[d] = exp->GetBasis(d)->GetBasisKey();
+        }
+
+        // if the basis is the same as the previous one,
+        // increment the number of elements
+        if (thisbasisKeys == prevbasisKeys)
+        {
+            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();
+            blockAttr.push_back({1, num_pts});
+            prevbasisKeys = thisbasisKeys;
+        }
+    }
+
+    return blockAttr;
+}

Field.hpp

 #pragma once
 
-#include <StdRegions/StdExpansion.h>
-#include <StdRegions/StdRegions.hpp>
 #include <array>
 #include <iostream>
 #include <memory>
@@ -14,7 +12,16 @@
 #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;
+}
+}
+
 
 /**
  * @brief Captures the structure of a block of elements of identical shape
@@ -49,6 +56,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 +66,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;
@@ -108,8 +116,9 @@ 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>>(

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/Operator.hpp

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

main.cpp

@@ -17,8 +17,8 @@
 
 #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 +41,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,7 +86,8 @@ 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;
 
@@ -195,14 +176,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

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>

tests/init_fields.hpp

@@ -27,26 +27,6 @@ struct InitFields {
   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)
-    {
-        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;
-  }
-
   InitFields() {
     BOOST_TEST_MESSAGE("Creating input and output fields");
     // Initialise a session, graph and create an expansion list
@@ -67,8 +47,8 @@ struct InitFields {
                     (session, graph);
 
     // 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);
+    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);