.gitignore

@@ -5,6 +5,7 @@ builds
 a.out
 .ccls-cache/
 .cache
+*.opt
 .*
 !.gitignore
 !.gitlab-ci.yml

CMakeLists.txt

@@ -27,7 +27,7 @@ set(SRC Field.cpp Operators/Operator.cpp Operators/OperatorBwdTrans.cpp Operator
 if (NEKTAR_USE_CUDA)
     enable_language(CUDA)
     add_definitions(-DNEKTAR_USE_CUDA)
-    set(SRC ${SRC} MemoryRegionCUDA.cu Operators/BwdTrans/BwdTransCUDA.cu)
+    set(SRC ${SRC} MemoryRegionCUDA.cu Operators/BwdTrans/BwdTransCUDA.cu Operators/IProductWRTBase/IProductWRTBaseCUDA.cu)
 endif()
 
 if (NEKTAR_USE_SIMD)

MemoryRegionCUDA.hpp

@@ -76,7 +76,7 @@ public:
     void HostToDevice();
     void DeviceToHost();
 
-    bool GetOnDevice() const
+    bool IsOnDevice() const
     {
         return m_ondevice;
     }

Operators/IProductWRTBase/IProductWRTBaseCUDA.cu

+#include "IProductWRTBaseCUDA.hpp"
+
+namespace Nektar::Operators::detail
+{
+template <>
+std::string OperatorIProductWRTBaseImpl<double, ImplCUDA>::className =
+    GetOperatorFactory<double>().RegisterCreatorFunction(
+        "IProductWRTBaseCUDA",
+        OperatorIProductWRTBaseImpl<double, ImplCUDA>::instantiate, "...");
+}

Operators/IProductWRTBase/IProductWRTBaseCUDA.hpp

+#include "MemoryRegionCUDA.hpp"
+#include "Operators/IProductWRTBase/IProductWRTBaseCUDAKernels.cuh"
+#include "Operators/OperatorHelper.cuh"
+#include "Operators/OperatorIProductWRTBase.hpp"
+
+namespace Nektar::Operators::detail
+{
+
+template <typename TData, bool APPEND = false, bool DEFORMED = false>
+void IProductWRTBase1DKernel(const size_t gridSize, const size_t blockSize,
+                             const size_t nm0, const size_t nq0,
+                             const size_t nElmts, const TData *basis0,
+                             const TData *w0, const TData *jac, const TData *in,
+                             TData *out);
+
+template <typename TData, bool APPEND = false, bool DEFORMED = false>
+void IProductWRTBase2DKernel(const size_t gridSize, const size_t blockSize,
+                             LibUtilities::ShapeType shapetype,
+                             const size_t nm0, const size_t nm1,
+                             const size_t nq0, const size_t nq1,
+                             const size_t nElmts, const bool correct,
+                             const TData *basis0, const TData *basis1,
+                             const TData *w0, const TData *w1, const TData *jac,
+                             const TData *in, TData *out);
+
+template <typename TData, bool APPEND = false, bool DEFORMED = false>
+void IProductWRTBase3DKernel(
+    const size_t gridSize, const size_t blockSize,
+    LibUtilities::ShapeType shapetype, const size_t nm0, const size_t nm1,
+    const size_t nm2, const size_t nq0, const size_t nq1, const size_t nq2,
+    const size_t nElmts, const bool correct, const TData *basis0,
+    const TData *basis1, const TData *basis2, const TData *w0, const TData *w1,
+    const TData *w2, const TData *jac, const TData *in, TData *out);
+
+// IProductWRTBase implementation
+template <typename TData>
+class OperatorIProductWRTBaseImpl<TData, ImplCUDA>
+    : public OperatorIProductWRTBase<TData>
+{
+public:
+    OperatorIProductWRTBaseImpl(
+        const MultiRegions::ExpListSharedPtr &expansionList)
+        : OperatorIProductWRTBase<TData>(expansionList)
+    {
+        size_t nTotElmts = this->m_expansionList->GetNumElmts();
+
+        // Initialise jacobian.
+        size_t jacSize = Operator<TData>::GetGeometricFactorSize();
+        auto jac       = Operator<TData>::SetJacobian(jacSize);
+        cudaMalloc((void **)&m_jac, sizeof(TData) * jacSize);
+        cudaMemcpy(m_jac, jac.get(), sizeof(TData) * jacSize,
+                   cudaMemcpyHostToDevice);
+
+        // Initialize basiskey.
+        m_basis = GetBasisDataCUDA<TData>(expansionList);
+
+        // Initialize weight.
+        m_weight = GetWeightDataCUDA<TData>(expansionList);
+    }
+
+    ~OperatorIProductWRTBaseImpl(void)
+    {
+        for (auto &basis : m_basis)
+        {
+            for (size_t i = 0; i < basis.second.size(); i++)
+            {
+                cudaFree(basis.second[i]);
+            }
+        }
+        for (auto &weight : m_weight)
+        {
+            for (size_t i = 0; i < weight.second.size(); i++)
+            {
+                cudaFree(weight.second[i]);
+            }
+        }
+        cudaFree(m_jac);
+    }
+
+    void apply(Field<TData, FieldState::Phys> &in,
+               Field<TData, FieldState::Coeff> &out) override
+    {
+        // Copy memory to GPU, if necessary and get raw pointers.
+        auto const *inptr =
+            in.template GetStorage<MemoryRegionCUDA>().GetGPUPtr();
+        auto *outptr = out.template GetStorage<MemoryRegionCUDA>().GetGPUPtr();
+
+        // Initialize index.
+        size_t expIdx = 0;
+        auto jacptr   = m_jac;
+
+        // Initialize basiskey.
+        std::vector<LibUtilities::BasisKey> basisKeys(
+            3, LibUtilities::NullBasisKey);
+
+        // Loop over the blocks.
+        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();
+            auto ptsKeys      = expPtr->GetPointsKeys();
+            auto deformed     = expPtr->GetMetricInfo()->GetGtype() ==
+                            SpatialDomains::eDeformed;
+
+            // Deterime CUDA grid parameters.
+            m_gridSize = nElmts / m_blockSize;
+            m_gridSize += (nElmts % m_blockSize == 0) ? 0 : 1;
+
+            // Flag for collapsed coordinate correction.
+            bool correct = expPtr->GetBasis(0)->GetBasisType() ==
+                           LibUtilities::eModified_A;
+
+            // Fetch basis key for the current element type.
+            for (size_t d = 0; d < expPtr->GetShapeDimension(); d++)
+            {
+                basisKeys[d] = expPtr->GetBasis(d)->GetBasisKey();
+            }
+
+            // Function call to kernel functions.
+            if (expPtr->GetShapeDimension() == 1)
+            {
+                auto basis0 = m_basis[basisKeys][0];
+                auto w0     = m_weight[basisKeys][0];
+                auto nm0    = expPtr->GetBasisNumModes(0);
+                auto nq0    = expPtr->GetNumPoints(0);
+                if (deformed)
+                {
+                    IProductWRTBase1DKernel<TData, false, true>(
+                        m_gridSize, m_blockSize, nm0, nq0, nElmts, basis0, w0,
+                        jacptr, inptr, outptr);
+                }
+                else
+                {
+                    IProductWRTBase1DKernel<TData, false, false>(
+                        m_gridSize, m_blockSize, nm0, nq0, nElmts, basis0, w0,
+                        jacptr, inptr, outptr);
+                }
+            }
+            else if (expPtr->GetShapeDimension() == 2)
+            {
+                auto basis0 = m_basis[basisKeys][0];
+                auto basis1 = m_basis[basisKeys][1];
+                auto w0     = m_weight[basisKeys][0];
+                auto w1     = m_weight[basisKeys][1];
+                auto shape  = expPtr->DetShapeType();
+                auto nm0    = expPtr->GetBasisNumModes(0);
+                auto nm1    = expPtr->GetBasisNumModes(1);
+                auto nq0    = expPtr->GetNumPoints(0);
+                auto nq1    = expPtr->GetNumPoints(1);
+                if (deformed)
+                {
+                    IProductWRTBase2DKernel<TData, false, true>(
+                        m_gridSize, m_blockSize, shape, nm0, nm1, nq0, nq1,
+                        nElmts, correct, basis0, basis1, w0, w1, jacptr, inptr,
+                        outptr);
+                }
+                else
+                {
+                    IProductWRTBase2DKernel<TData, false, false>(
+                        m_gridSize, m_blockSize, shape, nm0, nm1, nq0, nq1,
+                        nElmts, correct, basis0, basis1, w0, w1, jacptr, inptr,
+                        outptr);
+                }
+            }
+            else if (expPtr->GetShapeDimension() == 3)
+            {
+                auto basis0 = m_basis[basisKeys][0];
+                auto basis1 = m_basis[basisKeys][1];
+                auto basis2 = m_basis[basisKeys][2];
+                auto w0     = m_weight[basisKeys][0];
+                auto w1     = m_weight[basisKeys][1];
+                auto w2     = m_weight[basisKeys][2];
+                auto shape  = expPtr->DetShapeType();
+                auto nm0    = expPtr->GetBasisNumModes(0);
+                auto nm1    = expPtr->GetBasisNumModes(1);
+                auto nm2    = expPtr->GetBasisNumModes(2);
+                auto nq0    = expPtr->GetNumPoints(0);
+                auto nq1    = expPtr->GetNumPoints(1);
+                auto nq2    = expPtr->GetNumPoints(2);
+                if (deformed)
+                {
+                    IProductWRTBase3DKernel<TData, false, true>(
+                        m_gridSize, m_blockSize, shape, nm0, nm1, nm2, nq0, nq1,
+                        nq2, nElmts, correct, basis0, basis1, basis2, w0, w1,
+                        w2, jacptr, inptr, outptr);
+                }
+                else
+                {
+                    IProductWRTBase3DKernel<TData, false, false>(
+                        m_gridSize, m_blockSize, shape, nm0, nm1, nm2, nq0, nq1,
+                        nq2, nElmts, correct, basis0, basis1, basis2, w0, w1,
+                        w2, jacptr, inptr, outptr);
+                }
+            }
+
+            // Increment pointer and index for next element type.
+            jacptr += deformed ? nqTot * nElmts : nElmts;
+            inptr += nqTot * nElmts;
+            outptr += nmTot * nElmts;
+            expIdx += nElmts;
+        }
+    }
+
+    static std::unique_ptr<Operator<TData>> instantiate(
+        MultiRegions::ExpListSharedPtr expansionList)
+    {
+        return std::make_unique<OperatorIProductWRTBaseImpl<TData, ImplCUDA>>(
+            expansionList);
+    }
+
+    static std::string className;
+
+private:
+    std::map<std::vector<LibUtilities::BasisKey>, std::vector<TData *>> m_basis;
+    std::map<std::vector<LibUtilities::BasisKey>, std::vector<TData *>>
+        m_weight;
+    TData *m_jac;
+    size_t m_blockSize = 32;
+    size_t m_gridSize;
+};
+
+template <typename TData, bool APPEND, bool DEFORMED>
+void IProductWRTBase1DKernel(const size_t gridSize, const size_t blockSize,
+                             const size_t nm0, const size_t nq0,
+                             const size_t nElmts, const TData *basis0,
+                             const TData *w0, const TData *jac, const TData *in,
+                             TData *out)
+{
+    IProductWRTBaseSegKernel<TData, false, APPEND, DEFORMED>
+        <<<gridSize, blockSize>>>(nm0, nq0, nElmts, basis0, w0, jac, in, out);
+}
+
+template <typename TData, bool APPEND, bool DEFORMED>
+void IProductWRTBase2DKernel(const size_t gridSize, const size_t blockSize,
+                             LibUtilities::ShapeType shapetype,
+                             const size_t nm0, const size_t nm1,
+                             const size_t nq0, const size_t nq1,
+                             const size_t nElmts, const bool correct,
+                             const TData *basis0, const TData *basis1,
+                             const TData *w0, const TData *w1, const TData *jac,
+                             const TData *in, TData *out)
+{
+    if (shapetype == LibUtilities::Quad)
+    {
+        IProductWRTBaseQuadKernel<TData, false, APPEND, DEFORMED>
+            <<<gridSize, blockSize>>>(nm0, nm1, nq0, nq1, nElmts, basis0,
+                                      basis1, w0, w1, jac, in, out);
+    }
+    else if (shapetype == LibUtilities::Tri)
+    {
+        size_t nmTot =
+            LibUtilities::StdTriData::getNumberOfCoefficients(nm0, nm1);
+        IProductWRTBaseTriKernel<TData, false, APPEND, DEFORMED>
+            <<<gridSize, blockSize>>>(nm0, nm1, nmTot, nq0, nq1, nElmts,
+                                      correct, basis0, basis1, w0, w1, jac, in,
+                                      out);
+    }
+}
+
+template <typename TData, bool APPEND, bool DEFORMED>
+void IProductWRTBase3DKernel(
+    const size_t gridSize, const size_t blockSize,
+    LibUtilities::ShapeType shapetype, const size_t nm0, const size_t nm1,
+    const size_t nm2, const size_t nq0, const size_t nq1, const size_t nq2,
+    const size_t nElmts, const bool correct, const TData *basis0,
+    const TData *basis1, const TData *basis2, const TData *w0, const TData *w1,
+    const TData *w2, const TData *jac, const TData *in, TData *out)
+{
+    if (shapetype == LibUtilities::Hex)
+    {
+        IProductWRTBaseHexKernel<TData, false, APPEND, DEFORMED>
+            <<<gridSize, blockSize>>>(nm0, nm1, nm2, nq0, nq1, nq2, nElmts,
+                                      basis0, basis1, basis2, w0, w1, w2, jac,
+                                      in, out);
+    }
+    else if (shapetype == LibUtilities::Tet)
+    {
+        size_t nmTot =
+            LibUtilities::StdTetData::getNumberOfCoefficients(nm0, nm1, nm2);
+        IProductWRTBaseTetKernel<TData, false, APPEND, DEFORMED>
+            <<<gridSize, blockSize>>>(nm0, nm1, nm2, nmTot, nq0, nq1, nq2,
+                                      nElmts, correct, basis0, basis1, basis2,
+                                      w0, w1, w2, jac, in, out);
+    }
+    else if (shapetype == LibUtilities::Pyr)
+    {
+        size_t nmTot =
+            LibUtilities::StdPyrData::getNumberOfCoefficients(nm0, nm1, nm2);
+        IProductWRTBasePyrKernel<TData, false, APPEND, DEFORMED>
+            <<<gridSize, blockSize>>>(nm0, nm1, nm2, nmTot, nq0, nq1, nq2,
+                                      nElmts, correct, basis0, basis1, basis2,
+                                      w0, w1, w2, jac, in, out);
+    }
+    else if (shapetype == LibUtilities::Prism)
+    {
+        size_t nmTot =
+            LibUtilities::StdPrismData::getNumberOfCoefficients(nm0, nm1, nm2);
+        IProductWRTBasePrismKernel<TData, false, APPEND, DEFORMED>
+            <<<gridSize, blockSize>>>(nm0, nm1, nm2, nmTot, nq0, nq1, nq2,
+                                      nElmts, correct, basis0, basis1, basis2,
+                                      w0, w1, w2, jac, in, out);
+    }
+}
+
+} // namespace Nektar::Operators::detail

Operators/IProductWRTBase/IProductWRTBaseStdMat.hpp

@@ -14,48 +14,9 @@ public:
         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];
-            }
-        }
+        size_t jacSize = Operator<TData>::GetGeometricFactorSize();
+        m_jac          = Operator<TData>::SetJacobian(jacSize);
     }
 
     void apply(Field<TData, FieldState::Phys> &in,
@@ -80,7 +41,7 @@ public:
             // This is the B^{T} matrix
             auto const matPtr = expPtr->GetStdMatrix(key);
 
-            Array<OneD, NekDouble> wsp(nqTot * nElmts, 0.0);
+            Array<OneD, TData> wsp(nqTot * nElmts, 0.0);
             if (expPtr->GetMetricInfo()->GetGtype() ==
                 SpatialDomains::eDeformed)
             {
@@ -123,7 +84,7 @@ public:
     static std::string className;
 
 private:
-    Array<OneD, NekDouble> m_jac;
+    Array<OneD, TData> m_jac;
 };
 
 } // namespace Nektar::Operators::detail

Operators/Operator.hpp

 #pragma once
 
-#include "Field.hpp"
-
 #include <string>
 
 #include <LibUtilities/BasicUtils/NekFactory.hpp>
 #include <MultiRegions/ExpList.h>
-// #include <StdRegions/StdExpansion.h>
 
 namespace Nektar::Operators
 {
@@ -63,6 +60,60 @@ public:
     }
 
 protected:
+    size_t GetGeometricFactorSize(void)
+    {
+        size_t gfSize    = 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)
+            {
+                gfSize += expPtr->GetTotPoints();
+            }
+            else
+            {
+                gfSize++;
+            }
+        }
+
+        return gfSize;
+    }
+
+    Array<OneD, TData> SetJacobian(size_t jacSize)
+    {
+        // Allocate memory for the jacobian
+        Array<OneD, TData> jac(jacSize, 0.0);
+
+        // Initialise jacobian.
+        size_t index     = 0;
+        size_t nTotElmts = this->m_expansionList->GetNumElmts();
+        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)
+                {
+                    jac[index++] = auxJac[i];
+                }
+            }
+            else
+            {
+                jac[index++] = auxJac[0];
+            }
+        }
+        return jac;
+    }
+
     MultiRegions::ExpListSharedPtr m_expansionList;
 };
 

Operators/OperatorHelper.cuh

+#ifdef NEKTAR_USE_CUDA
+#include "MemoryRegionCUDA.hpp"
+#endif
+#include <MultiRegions/ExpList.h>
+
+namespace Nektar::Operators
+{
+
+template <typename TData>
+using BasisMap =
+    std::map<std::vector<LibUtilities::BasisKey>, std::vector<TData *>>;
+template <typename TData>
+using WeightMap =
+    std::map<std::vector<LibUtilities::BasisKey>, std::vector<TData *>>;
+
+template <typename TData>
+BasisMap<TData> GetBasisDataCUDA(
+    const MultiRegions::ExpListSharedPtr &expansionList)
+{
+    // Initialize data map.
+    BasisMap<TData> basis;
+
+    // Initialize basiskey.
+    std::vector<LibUtilities::BasisKey> basisKeys(3,
+                                                  LibUtilities::NullBasisKey);
+
+    // Loop over the elements of expansionList.
+    size_t nDim = expansionList->GetShapeDimension();
+    for (size_t i = 0; i < expansionList->GetNumElmts(); ++i)
+    {
+        auto const expPtr = expansionList->GetExp(i);
+
+        // Fetch basiskeys of current element.
+        for (size_t d = 0; d < nDim; d++)
+        {
+            basisKeys[d] = expPtr->GetBasis(d)->GetBasisKey();
+        }
+
+        // Copy data to basis, if necessary.
+        if (basis.find(basisKeys) == basis.end())
+        {
+            basis[basisKeys] = std::vector<TData *>(nDim, 0);
+            for (size_t d = 0; d < expPtr->GetShapeDimension(); d++)
+            {
+                auto ndata      = expPtr->GetBasis(d)->GetBdata().size();
+                auto hostPtr    = expPtr->GetBasis(d)->GetBdata().get();
+                auto &devicePtr = basis[basisKeys][d];
+                cudaMalloc((void **)&devicePtr, sizeof(TData) * ndata);
+                cudaMemcpy(devicePtr, hostPtr, sizeof(TData) * ndata,
+                           cudaMemcpyHostToDevice);
+            }
+        }
+    }
+    return basis;
+}
+
+template <typename TData>
+WeightMap<TData> GetWeightDataCUDA(
+    const MultiRegions::ExpListSharedPtr &expansionList)
+{
+    // Initialize data map.
+    WeightMap<TData> weight;
+
+    // Initialize basiskey.
+    std::vector<LibUtilities::BasisKey> basisKeys(3,
+                                                  LibUtilities::NullBasisKey);
+
+    // Loop over the elements of expansionList.
+    size_t nDim = expansionList->GetShapeDimension();
+    for (size_t i = 0; i < expansionList->GetNumElmts(); ++i)
+    {
+        auto const expPtr = expansionList->GetExp(i);
+
+        // Fetch basiskeys of current element.
+        for (size_t d = 0; d < nDim; d++)
+        {
+            basisKeys[d] = expPtr->GetBasis(d)->GetBasisKey();
+        }
+
+        // Copy data to weight, if necessary.
+        if (weight.find(basisKeys) == weight.end())
+        {
+            weight[basisKeys] = std::vector<TData *>(nDim, 0);
+            for (size_t d = 0; d < expPtr->GetShapeDimension(); d++)
+            {
+                auto ndata = expPtr->GetBasis(d)->GetW().size();
+                Array<OneD, TData> w(ndata);
+                if (expPtr->GetBasis(d)->GetPointsType() ==
+                    LibUtilities::eGaussRadauMAlpha1Beta0)
+                {
+                    Vmath::Smul(ndata, 0.5, expPtr->GetBasis(d)->GetW().get(),
+                                1, w.get(), 1);
+                }
+                else if (expPtr->GetBasis(d)->GetPointsType() ==
+                         LibUtilities::eGaussRadauMAlpha2Beta0)
+                {
+                    Vmath::Smul(ndata, 0.25, expPtr->GetBasis(d)->GetW().get(),
+                                1, w.get(), 1);
+                }
+                else
+                {
+                    Vmath::Vcopy(ndata, expPtr->GetBasis(d)->GetW().get(), 1,
+                                 w.get(), 1);
+                }
+                auto hostPtr    = w.get();
+                auto &devicePtr = weight[basisKeys][d];
+                cudaMalloc((void **)&devicePtr, sizeof(TData) * ndata);
+                cudaMemcpy(devicePtr, hostPtr, sizeof(TData) * ndata,
+                           cudaMemcpyHostToDevice);
+            }
+        }
+    }
+    return weight;
+}
+} // namespace Nektar::Operators

cube_all_elements.xml

+<?xml version="1.0" encoding="utf-8" ?>
+<NEKTAR>
+    <EXPANSIONS>
+        <E COMPOSITE="C[0]" NUMMODES="4,5,5" FIELDS="u" BASISTYPE="Modified_A,Modified_B,Modified_C" NUMPOINTS="5,6,6" POINTSTYPE="GaussLobattoLegendre,GaussRadauMAlpha1Beta0,GaussRadauMAlpha2Beta0"/>
+        <E COMPOSITE="C[7]" NUMMODES="4,5,6" FIELDS="u" BASISTYPE="Modified_A,Modified_A,Modified_A" NUMPOINTS="5,6,7" POINTSTYPE="GaussLobattoLegendre,GaussLobattoLegendre,GaussLobattoLegendre"/>
+        <E COMPOSITE="C[8]" NUMMODES="4,5,6" FIELDS="u" BASISTYPE="Modified_A,Modified_A,Modified_B" NUMPOINTS="5,6,7" POINTSTYPE="GaussLobattoLegendre,GaussLobattoLegendre,GaussRadauMAlpha1Beta0"/>
+        <E COMPOSITE="C[9]" NUMMODES="4,5,5" FIELDS="u" BASISTYPE="Modified_A,Modified_A,ModifiedPyr_C" NUMPOINTS="5,6,6" POINTSTYPE="GaussLobattoLegendre,GaussLobattoLegendre,GaussRadauMAlpha2Beta0"/>
+    </EXPANSIONS>
+
+    <GEOMETRY DIM="3" SPACE="3">
+        <VERTEX COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJx9lTlTFFEUhdt9A9cRcQM3BMUVwV2vOyq4gDSiAu6oiEtiomaWEQGRf2ICMquMrfIfQELiaNVLMCExMLN0+syrOfPOe8mr11+fc8+93dOTJMVV8/lPT5pOWnZM6j/+6vpSN2mzkvACny304+0/3ubGv9uc7Dzxqao3Tb8Z9qH/a8rmCn/o59F11EH9+YJDvyDan7OFIt+73D+Hn7Yo2r+zxdF8zpbQdeZVpH/T8CFNU1fqr5ryFU+F0vyWCj36Xyb06G95tL6zFZSf+crgdGZK+VeRnnlO6JF/NeUfnXqfz+enS/3XCH/o10TzOasV/pjP2mh+Z+tobqiA/tYLjnwbynIVKn4fG6P+zuqEHvnrSc/v36agv38+m8mf579F6NHfVqFHvm3R+s4aiPP7t73Mf6Zi/o2CI19TdkYunv8O0vP3Z6fg8G+O5nO2S9THfHZH6zvbQ5y/H3sFR759VL948t+X/VF/Zy1Cj/wHgv3796tVcORrI39+/w5G/Z0dEnrkO1ymL1TkP5LE19Fs5/8XrGPE+b7jdD/zE0KP+Z0kzt8JDsT8lNDD/3Qwl3+PzwSn4vlZURfrHHHu/7zgyHeBONdpFxz6i8F8Pv8l0RfWZeKcv4M45+sM+vs6V7Kdnz9+H1fFfMCvif7Br5M/8y7qF/eNff3d1D82bd08kGyB3xD5wXuEP+aQCn/wXuEPfjPo6/P1CQ79LaqI+UB/m/TM7wg9/PuFHnxAzAX+g2Iu4HeFHvs9ocd+n3Lz83sgOPQPg337+T8SHPrHxHm+Q4JD/yTo7+fzVHDszygP9z8sOM7Pg3P1/Y8IjvMLOnP/LwXH+VXQ3/f/WnDsfwFYCm9Q</VERTEX>
+        <EDGE COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJx1mmWQFkcURRk0sNjiIcgii0MguC9OcHcL7iS4u7sEjbu7u7u7u7u7/nlnq+ZUwZ+ue+b2nf56emem35AnT/pfcpw2r9o80vnE0fmVAy9wnP4F5ad/IfnhJ0iTU1g5nK+INP0yxMkpKs74iqk/vLg0OSXk5/eUjDaveKb8/M5S8sFLKwdeRn7mpazy4eWkySkfbUEdryCO/0T5mPeK8sNPkianknK4TpWjLSReRZqcquKcJytari/Xu5o0460uP+ughvzwmuoPz5afdVNLfnhtaXLqiPO76oqz/urJB6+v88IbyM96bah8eCNpck5WDuu7cbQZ4k3EyTlFPnhTcc7fLFr+bvg7aS7N30kL+eEtxfn7aaUceGv1h7eRn7+3tvLD20mT0z5a1g/j6iDN9ekoPzk58sE7iTOuztLcF7ponPCu0twvuikH3l058B7K4f7SU374qdLk9NJ5WSe9pfH3iba4eF/5GVc/cc7fP9pM8QHi5A+Mlvsl/wZJc98cLD98iDj306HKgQ9Tf/hwcZ4zI6TxjxTnfj1KOfDR8nMfHyMffKxy4OPk5/kwXhr/BHGeE6cpBz5RmpxJ4ty3JkvjnxIt94GsaKdKs66nyc/zabr88BnqD58pP8+zWfLDZ0uTM0ec9TxXnOfiPPngp+u88DPk5zk6X/nwBdLkLIw2U8cXRZstvlicfkukyV8qP3xZtDzP+ftbLs1zfYX88JXiPO9XKQe+Wv3ha8T5u1krjX9dtLxHcP710rxPbJAfvlGc/E3ywTeLM94t0ryvbNX44dukeY/Zrhz4DuXAdyqH+8Euafy7lc970p5o64nvlSZnn/zw/fIzrgMaD8fPlOb+dFCa/EPR8h7HdTsszfvcEfnhR8V5zzumHPhZ6g8/W5z5OEca/7nivEeepxz4+fLzfnmBfPALlQO/SH7m+2Jp/JdEy/zT71Jx/JeJc7+5XP3hV0iTc2W0vC9zfa6S5r35avnh14jzPn2tcuDXqT/8enHm9QZp/DdGm6Xz3yTN8+Fm+cm5RT74reKM6zZp9gm3a5zwO6TZP9ypHPhdyoHfrRyu8z3S+O9VPs+1+6Tx3x9tdfEH5Of8D0abLf6QODkPR5uh449Ey3pmH/Wo/PDH1B/+eLSZ4k9Eyz6MfdeT0uy7npJm3/W0NHnPSFNnelb5/HtO+ezfnlc+dYUX1B//i9L4X1I++8CXpdkHviLN+V+VPyfa1+RnH/m6xoP/DWn8byqffeZb0uwz35bm/O/Iz3vxu/KzD35P48H/frRcL/axH0izj/1Qmvn+SH7eyz6Wn33vJ/LBP5VmXj9TDvvkz6WZ1y/E2T9/qXz4V9I50X6tHPbb30hzHb4Vz4r2O+XDv4+WeWPf/oM0+/YfpZmPn+Tn/fRn+Xnv+EWa+fhVfuoCv0kzH7+L89z7Q5rf/af81Bf+kmb+/hbn+fZPtKwP6hH/RsvvoR7xn8ZBPYIbV05I6hFJks4lJ2+S1uTkS9KanPxJWpNXIDj1DOodBcWp+xZK0r7cOnFw6h8cL6wcjheRn+dMRnDqItRNioqTVyxJ+zhv8eDUUTheQjkcLyk/48rU+Km/lErSmrpFaXFyyui81GvKJmlNv3LijKu85p/6ToUkrXPrXsFHyV9RObn1IPl5nlfSPFMnqpykNflVgo+Rv6pyOJ4lP+evpnUCr67rRb8awalLUYeqmaQ1dcJs+alz1NJ8cry2OHl1glPfmhBt3SSt8deTn/z6mk+ONxAnr6GuFzmNNJ/knBycuhrz2ljrluNN5M/dDwWnDsd1aKp1y/Fm8pPXXOelftciSWv20y3FyWml81Lva52kNf3aiDOutsGpCzKv7bQOOd5eft6LOwSnjsh16Kh1y/Ec+cnrpPmHd9Z80q9LcOqUrIeu4uwnumn9cLx7cOqarJ8e4vh7ar1x/FTND/5e+r34e+v6Uj/tk6Q1+/i+4sxTP11f6q39k7Sm3wBx5nWgxkN9dlCS1tR/BouTP0TjoZ47NElr+g0TZ1zDtR6o/45I0pr96Ejx3HWt9UO9eHSS1vQbI851G6v5xz9O80m/8cGpQ7MeJoiznz1N64fjE4NTt2b9TBLHP1nrjeNTND/4p+r34p+m9zG+k0wX5/8bzAjeXnxmcL7z0m9W8Azx2crne/occb5Hz9V54fP0Hgg/Xe+f8DOC892NfvN1Xr7/LgieKf9CnRf/InG+ey7WeOBLgvOdgn5LlcP3u2V6T6auvFyc72UrgncTX6nfC1+leYavDs7zkvOs0Xn5DrVW+fB1wXmfrxbt+uA9xDdo3uAbNT/wTcGpO3GezcGzxbdoPHxn2ap8+Db9Xuq728X57rBD8w/fKT918l3yU3/eLT98j8ZPnW+vOP59mmf4fuVT1z2gHPiZyqHOeVB++CFx6niHdV7qbEe0fuBHlQM/pusF/x9lNlD9</EDGE>
+        <FACE>
+            <T COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJyFmVm4jlUYhrcQoqJIEUqGRMKuZChzbaESGZvIkLYUTUKZsiOKohBKSlEhQoMGjdIg6WqejuqgDjrqoIOGk+c+2Pd1fdfeJ8/1vu/zPOtd61/ft/9/rZKSyn/Vgkcp5q9GsHoBr2YVvJrKwzs6WKvAj3rtKnQ1xSNfJ3hMgR/1uvK3rrZ4xwjrBY8Vj/xxBfljVa+t/PHB+sFa4lM/IdigQHeiePaj3kDYUH6N5NNQ8Unyg99A9YaqM9/GBfmTC3waq+78KYqbiE8dXVPlzW+qPPM4NdhMPPLNNY751E9Qnj5aFPi1UNxU+dOCpwfrKQ+eEWwpPtiqgFdP9ZbitQ62KfCj3rYKXSvxyJ8ZbKJ8W9XPCrYr0LUv4LVRvZ14HeR3tvgdxO+oPPx2qpPnczwn2Ek86vh3Fh/sIl4nIfXOQvosLfArVdxR+XOD5wXPVx68INhVfLCbeOcLqXcVou8e7CEf8j2V7yZ+T+XRXRi8KMi+7aE6/fSSH7re4qFrrTp6dMy/j/K9lEffTXFf9dtPOuoDgv0LdBeL109Ivb/wEvHK5HOJ+APlV6aYOrpLgzxPzKtM9UEaFx90g8VD10f1AeINEb+v+EPE76/8ZcHLg6XKg/Q7SHzwCvGISxWjHxpkHRmHvq4Uj3iYkPpVweFCdCPEA4epPlwx3y9Hym+k6qPkO0q8EcqPFo4Rf4zqVwfHCuFdIx44UvWxiq/VOPWVB/G7Tgh/rPLXS088Tkh+vPTjxSOmjxuCzRTjN0F1YvwnSjdevEniEfNenqz6jcoTo5uiPHF3xehvEpbLr1z1qfKfKt4U5W8Wlok/TfVblCdGd6vyxOWK0U8XMs4MIXV8bxPCm6b87cI7gncqpo7urgLezCp4M5Unpn/G4X15t2J0zHOWeMSzhdTvCc4RortXPHC26nMUo+c9OVdIfp5854lHjO984YLgQsXU6eu+At6iKniLlCdm/hXyqxDCp//7hcxzsZD6A8ElQnRLxQMXq75E8TLxxigPontQCH+p8g8JlwcrFFNfob5XiPew8sQVitE/ImQc1nOl6vgyz1XBZQX8R5VHz/wfU53/Q/iuVh4/9GvEp844azUOfugeV50Y3Trp1om3Xn7U0W8ILhQfvydUh4//k9KtE2+jfMx/SnVidJuk2yTe0/Kjjv4Z8RiH/bhZdXzRPSu+fZ+Tns8Z3y3K44d+q/jU6eN51RkP3xeUxw+fF8Wnjv82jYMfuu2qE6PbId0O8V6SH3X0O6Vnn6LbpTx8/F8Wnzrj7JYPfL4v7VGdvtDtlW6PeK8U+KF/Ncj87fea8ujp/3WNwzzx3ac6fHzfkG6feG8W+KF/Sz7weB7fVh0dvvvF36f6OxqH73Gsx7uqMx79vyfdfvHelx919B8U9MF6fqg6vszzQJD585zR90fK0y/zPyg+dfw/1vjw+D76iWL09PlpkN83rAfvi89UZ76syyHpDoj3ufyooz8c5PyP78n4faE664X/EekOiFc9Fzk1g3WCtYO1gnWDxwePU75B8ET5wGsYPEk+6BoHmwRPCTYKnhpsEWwePFn9ni6fJuqjpXyaBc8Itg22kd+ZwQ7B9sFW4tP32cq3V98d5cu6dQqeIx3rVhrsIt55wa7ygXdBsLt8ztU8ewZ7yI916xW8SH700Uc+PdVHX/mA/bRuA8S/WOtWFuwvPn0PlK5MfV8q3qDgEM17cPCy4BWa9+XBocFh8rsyODw4Qn5XBUcGR2sdRwXHBK/WuowNXiOfIcpfKx/mcZ2Q9RgXHC9k3SYEJwrHaR6Tg5OkYx5Tgv+UVI5vDN4U5D1XXbxyjcP9Mvqbg0ek4353msZBx73tdI0zlec49Rkan/tUeLfLnzz3lHdo/tPU953qv6V4d8mf+zf0M4Nc2N8avCV4d/Bb+cKbFayhvtDPCX5XUlnH/eA9iadLx73fXK0v68m91jzNj/s6eAuCt2m9uedaqP5niLdI/ZHnnq9C68fnwed2v+bPfQ68xfLvXVJZ/0DwF+m431qqzxkd91IPav583twTPKTxuc+At0L+5LmneSTxv4lna1+sDP6a+nDxVmmfcX6P/lHtf/YL++qx4G/yhbc6WE19oV+j9wc67hfW6vlFx/3C+sT3BtmvnHNv0Py4D4D3pPY3+5nz8Y3qf654T6k/8twTbEr8X2Keh/nBp4N/pM7zAO8Z7W/q9wU3B6srj/9zwT81PufnW/ScoeP8fKv2P+f29PG89ifn7/SxLchzyvPIufF29b9YvB16fqhzjrxT/aPjnHyX+uN8Hp/den54nnnu9wT/Yl7i7Q0uoY+SyvpX9Pyh4/z4Ve1PdPNS36f9vYznMvU3ND7ni/DelD95zkHf0vx4X3H+9bbWhzrvrf0anzznY+9q/z+s99p7wb9TXyve+8HlQc7f0H8Q/F06+vwweJTGpc+P1D/v043xO6j1p8750yfa37xvOYf6VP1zfgfvkPpbqXE+1/rwvuZ86LDWlzrv7S80PnnOkY7o/cC5Dj5f6vni/Aafr+TP/wPOQb7W+lDnd/w38t9fUpn3rdaf/yf8/v4uMe/nNeJ9HyxXnd/fP6g//t/wO/xHfT6rxftJ76fH2VfR/6z321rxfgnWUJ4+/weSxzHH</T>
+            <Q COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJx1mGtsnnMYxt+nHTNaihWdjTqsbDbbqNHNMEwddjQzNptVDTEJc8wiIY1DIlmyOGQRhA/ikPFhIZEQXyYSMea0GRvDaHVrp9gwzOmD6yd5ftH3y5X7vq/7vq7/4e37PF1f+ffzdnBt8J3gu8EN4r1XKX/WBT/WnPeDH6pvo3TgfRT8IPiJfJAfJN1P5Qc95sHfJN31QdaH/82q4xcd/H6mfnjo4P9zzSffGGRdW+QbvU3ifyFd/OIHvS9V36I8fr9SP7xm6W3VfPKTgqzra/lGr1n8b6TPfPq5R53i0TcryL3q0pxO5en7VjqdynO+3dLpEo97tk3zupXnHm1XnbhHPoi3idcb3KgYPz3Kc392aB7xd8HNiqn3Bb+Xjz7xiRcHuT/09SqGx734QfhjcGdwq/LwdgWXBTnfneL/FOxQ3y7x4f0c7FQenQ7xuhQz75fg7mC38uCvwZXS2y2E91twu/LorxSvRzHzfg/uCe5QHvwjuEp6e4TwXgr2KY/+KvG4F38Guc9/BYuiXCfPuv4OVhXlOvlexdUFg8v43+9I4gHSxQeIXrX45PcK7h1syJyB4lMfmvq+iQdpDn1HireP5uwXrFF/rfg10qkLHqC8++DtX5R9HBg8KDi4UvZVo/owreNg6Q0O1mt+jXzWqx9/hyg+VPs4KvqHyRf1MakfnniIdOg7WbwG+R+qdeFjmHzWSqeR8+6nr0G8I5RnX8fK31FFOa5TnnlHy/8xwWODrZk7XPrUp6XOPh2n9dA3Xbwm+Tg+OEI+RsrnCOmwT6P66WsS7wT5AmfEH/s1Wjojlcf3ifLF/VqUeeM4H83h3NrFG6P11Ku/mfsoP3Xqh3eSfJ4SHB9cWCn7Hqt6m/ydWpTj04Itml8vny3qx98E9bEvb0Z3UuLT5Ytzul28ifJVp/6zgmfKV6P64Z2hPOt+TrqTi3LcrDzzztbcc3QPl2buueojvjV11j9FOufpnjG/tZ94itZ/vvLo3qb1XqA5rcrj+0LNvSg4Nfhw5o4Xn/ojqbPv0zRnOt9f5VvEn6GY85mpObOCFweXy99M1R+qlP3O1pxLgnOUbxF/jmL8Xqo5c4OXBZ+J/jjxqT+fOud2uebMC85Xvln8+Yo53ys0Z4H0n47+VPEXau4CrZv6lbqHazJvUfH/+IbW2ya8SveV+e39YJvWe7VwsfRfjz7rv0bYrj7Wf12Q51DeO6uC16Z+Q7Amef6PQN/1qd9YlPsGpl4r3k3BJeLXhcfz2y3JD1E/fTyv3pl4eGL+P8Dz9tLgXUGej5nXVCnzOoqy/gDp8txzb+LR6qeP57cHEk9IzPsp/u7ge6+YeZNDnxhcwd+3SpmPLs8VDybmvYX3ydnBZak/VpTjualzLryXPK59hD8v/Lbgk9qvQn08fzzLOSTmfY33qftSf6Eox9XapyXhv5j4fvHYF363Vye+WTr08bvzsvaP9+B7gpzrq/Qnf3eQ/YH3mvaB/Irw+d1Zkzz7wPvto0HO/63gE/LDeyK8tVof+afSx9/rdfJXpT6ePzYkXi0/7OMryf8DbeJGigAA</Q>
+        </FACE>
+        <ELEMENT>
+            <A COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJxd1edfj1EcxvGMZO89M1JRVpQtMysS2RlZGZGdRNmbUPbIXmUnUvaO0F/kyed64Po9eb9e97l+55z73Od8T0DA/79qWB1r2POaGIi1LC+DsDbWsf+p37rWXg/r23gNsCE2wsbWbxNsavlmNt/m9rwFtrT5t8LW1t4G9X5tbdx22B71fh2sv46W13p0wmDsjF1Q69QVu2GI5bWO3TEUwzActb49sKflI1DrGYm9sDf2QX2fvtZfP8vr+0VhfxyA0ajvGmPtA3EQav0H4xAcisNQ3324jTPC8tonsTgSR+Fo1P4Zg2OtfRxq/8TheJyAE1H7bhJOtrzmoX0Uj1NwKiag9us0m5/y6l/7OBGn4wxMQu3XmTgLZ1te+38OzsV5OB91LpJxgeUXos7LIsstxhTU+VmCS3GZ5YNxOa7AVFyJOk+rcDWusbzOZxquxXWYjjq363GD5TeizvUm3IxbcCvq3GdYf8pvQ9WBTNyOWbgDde53Wns25qDqyC7cjXtwL6rO7LNx9lte9ecAHsRDeBhVT47gUWs/hqpbx/EE5uJJVD07hactr3mo3uXZuPn2P9W3M9Z+Fs+h6uV5vIAX8RKqjl7GK3jV8qqvBXgNr+MNVL29ibcsfxtVj+9Y7i7eQ9Xj+/b8ARai6mwRPsRH+BhV559Y+1N8hqr/z7EYX2AJ6n54aeO8srzujVJ8jWVYjrpX3uBba3+Hukfe4wf8iJ9Q99Jn/GJ5zUP31Vcb9xt+R91jP6w/tVdgLP7EX1iJv1H3zR/8i1WW/wfhWnRd</A>
+            <P COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJxdlUlsjVEUx++LRCQkbCwkIiIkJCIiFkREYiFhQyJCQiIRYip9aKsofU8HqsNr6Tyjg+o8oRSlRWmrWlPRkthY2lhYmBb9/SX3vM0v/3Pu/5zz3ffd+wXd5O9AYJIH4SE4hfxh9DF0FHo6+ojxH4XHyY9Bwv/ryzeN+AkYHfD7TiU+w/n1Y9BBzae+cCb5uTBWfcjPQccY/zzicWaeWc6vL99s4ieh9icWxsEFWoeOR2tj4uEp4z+tOMvPoBeiE4z/rHzkJ/Q8pr5881UXnjN9Ne8i59dPQGs/E2EILia/TD4YJn8eJhn/cs0Bk4kvRacY/wriiWb+kKkv3xLWhaD2MxVegGvJX0SH0WnoVehLxp+u+ch/dX5drV9JfA1MghmmjvqvI78aJovkM2EW3Eh+C0yBEfIb0NnGv5l4qplHdbV+PflNej6YQ/4yvAK3kc/Vc6F/wzzi+cZfANNYp/tqO7ow4PuLtM/kv8Gtzp9HfeRPJ18c8Ofc4fy6mieDeAm6FJbBneT3wExYTr4CVhr/PtZlQd1vu9FXjX8v8YiZf5fz51Ef+bPJ/4XXiF+HVbAa5rCuRvczujbg+2/o/yOfZurV6XtB/qbeE/R358+h+vud789F15u+mld1o1iXBxuIN8Im2Kz7jnX5sIW4vk+txq/vRAGMmHptMJq87vFCM7/mUP2g8/1F6HZ0B+yEug9uoYvROo+6H28b/x1YQv4H1PdA95j6yF9KvEv3Hzrk/Lqap4z4XfQ92G3eR+1LOdT50vm5b/za9wozj+4LvbfqI38lce2n+j1AP4Q9OjfkH6Efw16zH306D+g/zq/3xPif6tyx7qeZQ/9jr/FXEX+m+8D586qu5qkmr/2vQfejn+scOT//gvgAHDTP1wBr4ZCp99L4653fR/P3m30YNP464tr/NvQw+hUc0fnXc5j8qN5j8q/NXL9M3SZ0p/P9eu436Hbn91fdLuKNUPvZg36Lfqd9cH7fFpN/r/eO/ABsNvOobjfxPuf7W9D6Ho3BD/AjHGJdK9T/+wmOG/8w6/T/6L6dgJ+Nf5R12r+w6f/F+EdY1wEDZu5x00f+f9N28hAA</P>
+            <R COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJxdlEdOQ0EQBT8552xyjiYHk3M2JmMwycbGZNiw4CCIKyBxESTEgiOxcL3FbzalUvdjPNMz3+8k/pKTEkyBqbCWehq+hafjGTDT5LNggP5sPAfmQo/jXm8bz8PzYYHJa70d+gvxIlgMB6iX4Lt4KV4Gy02+Au7Rr/1Xwiroddzr7ePavwdWm7zWO6C/RucN68zvrccP8Qa8ETaZfDM8MvkWqPPQPFvxoOPOt5nzUL4dHtOveXbAThil3oWf4JpnN+wxeS8M0d+L98F+GHbc653iA/ggHDJ5rXdGv+Y5DEfgH/VR/BzXPMegz+TH4QX92v8EnIS/jnu9S1z7n4LTJu8z+5/BZ6Hm+UZ9Do/g8/iCmafyi/DK5JegzuOV+rKZt/Ir5jyUX4Ux+tfwdbgBP6hr3tf4Ju7XPTZ5zTtOf0DfF3034LtZ/wbf1fvXuzZ5rX9L/4Hep94d/KKued/hQb0fvQuT17zv6Q/pfuu7AT+p67we8DPdV71rk9d5PdJ/ofsHNc8f6mH8CY/ofph5Kh+Fz/THND/dW/hNXef1gsc1D/P/rs15/QM7bzcM</R>
+            <H COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJxdz8diAQAQRVF+TRefr0ZE7wRREp2FYzOzOZu3uJNKvK7MCqussc4G0/wMuya/2OI3M2yHXYdd9thnlgMOOeI4dE+Y4zTsZpyH7gXzXIbdD1ehe80CN/zllrvQvWeRh7D743/oPvKDp7A78xK6ryzx/e+Ndz6YSIJPE3E3WQAA</H>
+        </ELEMENT>
+        <COMPOSITE>
+            <C ID="0"> A[0-65] </C>
+            <C ID="1"> F[211,217,222,227,232,236,240,245,249] </C>
+            <C ID="2"> F[4,0,17,13,86,82,93,89,126,122,139,135,212,228,241,266,269,302,305,334,337] </C>
+            <C ID="3"> F[133,129,146,142,157,153,164,160,177,173,190,186,250,246,242,336,339,346,348,356,358] </C>
+            <C ID="4"> F[184,180,197,193,106,102,113,109,55,51,68,64,251,238,224,350,353,318,321,284,288] </C>
+            <C ID="5"> F[62,58,75,71,35,31,42,38,10,7,23,20,225,220,215,286,290,272,276,258,262] </C>
+            <C ID="6"> F[263,270,277,282,289,294,299,306,311,316,322,326,332,338,344,349,354,359] </C>
+            <C ID="7"> H[66-75] </C>
+            <C ID="8"> R[76-111] </C>
+            <C ID="9"> P[112-180] </C>
+        </COMPOSITE>
+        <DOMAIN> C[0,7,8,9] </DOMAIN>
+    </GEOMETRY>
+</NEKTAR>

main.cpp

@@ -192,7 +192,7 @@ int main(int argc, char *argv[])
     // on the GPU.
     BwdTrans<>::create(explist, "CUDA")->apply(in, out);
 
-    // Test the GPU-backed fields with a CPU operator
+    /*// Test the GPU-backed fields with a CPU operator
     // This should show a debug warning due to the implicit conversion. The
     // purpose of this is to allow us to transition the code to the new
     // infrastructure and add CUDA operators, without the need to add ALL CUDA
@@ -202,9 +202,9 @@ int main(int argc, char *argv[])
     // Create two CPU backed fields and use them with a GPU operator
     // This call implicitly converts the fields to a GPU backend and warns the
     // user about that fact
-    in  = Field<double, FieldState::Coeff>::create(blocks);
-    out = Field<double, FieldState::Phys>::create(blocks);
-    BwdTrans<>::create(explist, "CUDA")->apply(in, out);
+    in  = Field<double, FieldState::Coeff>::create(blocks_coeff);
+    out = Field<double, FieldState::Phys>::create(blocks_phys);
+    BwdTrans<>::create(explist, "CUDA")->apply(in, out);*/
 #endif
 
     std::cout << "Inner Product of a function WRT the Basis" << std::endl;
@@ -223,7 +223,7 @@ int main(int argc, char *argv[])
             {
                 // Each element is the index of the quadrature points
                 // this is useful for testing reshapes
-                *(y++) = phys;
+                *(y++) = phys + 1;
             }
         }
     }
@@ -238,8 +238,62 @@ int main(int argc, char *argv[])
     // ... 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;
+    // Check output values.
+    std::cout << "Out:" << std::endl;
+    y = outCoeff.template GetStorage().GetCPUPtr();
+    for (auto const &block : outCoeff.GetBlocks())
+    {
+        for (size_t el = 0; el < block.num_elements; ++el)
+        {
+            for (size_t pts = 0; pts < block.num_pts; ++pts)
+            {
+                std::cout << *y << ' ';
+                y++;
+            }
+            std::cout << std::endl;
+        }
+        std::cout << std::endl;
+    }
+
+#ifdef NEKTAR_USE_CUDA
+
+    // Create two Fields with memory on the GPU
+    auto in_cuda =
+        Field<double, FieldState::Phys>::create<MemoryRegionCUDA>(blocks_phys);
+    auto out_cuda = Field<double, FieldState::Coeff>::create<MemoryRegionCUDA>(
+        blocks_coeff);
+
+    // Assign input values from the CPU.
+    y = in_cuda.template GetStorage<MemoryRegionCUDA>().GetCPUPtr();
+    for (size_t i = 0;
+         i < in_cuda.template GetStorage<MemoryRegionCUDA>().size(); ++i)
+    {
+        y[i] = inPhys.GetStorage().GetCPUPtr()[i];
+    }
+
+    // Perform the IProductWRTBase on the fields using the CUDA implementation
+    // Since this is a CUDA operator, acting on CUDA fields, everything happens
+    // on the GPU.
+    IProductWRTBase<>::create(explist, "CUDA")->apply(in_cuda, out_cuda);
+
+    // Check output values.
+    std::cout << "Out:" << std::endl;
+    y = out_cuda.template GetStorage<MemoryRegionCUDA>().GetCPUPtr();
+    for (auto const &block : out_cuda.GetBlocks())
+    {
+        for (size_t el = 0; el < block.num_elements; ++el)
+        {
+            for (size_t pts = 0; pts < block.num_pts; ++pts)
+            {
+                std::cout << *y << ' ';
+                y++;
+            }
+            std::cout << std::endl;
+        }
+        std::cout << std::endl;
+    }
+
+#endif
 
     std::cout << "END" << std::endl;
     return 0;

segment.xml

+<?xml version="1.0" encoding="utf-8" ?>
+<NEKTAR>
+    <GEOMETRY DIM="1" SPACE="1">
+        <VERTEX>
+            <V ID="0"> -1.0  0.0  0.0</V>
+            <V ID="1"> -0.8  0.0  0.0</V>
+            <V ID="2"> -0.6  0.0  0.0</V>
+            <V ID="3"> -0.4  0.0  0.0</V>
+            <V ID="4"> -0.2  0.0  0.0</V>
+            <V ID="5">  0.0  0.0  0.0</V>
+            <V ID="6">  0.2  0.0  0.0</V>
+            <V ID="7">  0.4  0.0  0.0</V>
+            <V ID="8">  0.6  0.0  0.0</V>
+            <V ID="9">  0.8  0.0  0.0</V>
+            <V ID="10"> 1.0  0.0  0.0</V>
+        </VERTEX>
+
+        <ELEMENT>
+            <S ID="0">    0     1 </S>
+            <S ID="1">    1     2 </S>
+            <S ID="2">    2     3 </S>
+            <S ID="3">    3     4 </S>
+            <S ID="4">    4     5 </S>
+            <S ID="5">    5     6 </S>
+            <S ID="6">    6     7 </S>
+            <S ID="7">    7     8 </S>
+            <S ID="8">    8     9 </S>
+            <S ID="9">    9    10 </S>
+        </ELEMENT>
+
+        <COMPOSITE>
+            <C ID="0"> S[0-9] </C>
+            <C ID="1"> V[0]   </C>
+            <C ID="2"> V[10]  </C>
+        </COMPOSITE>
+
+        <DOMAIN> C[0] </DOMAIN>
+    </GEOMETRY>
+
+    <EXPANSIONS>
+        <E COMPOSITE="C[0]" FIELDS="u" TYPE="MODIFIED" NUMMODES="4"/>
+    </EXPANSIONS>
+
+</NEKTAR>

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>

square_all_elements.xml

+<?xml version="1.0" encoding="utf-8"?>
+
+<NEKTAR xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
+    xsi:noNamespaceSchemaLocation="http://www.nektar.info/schema/nektar.xsd">
+
+    <GEOMETRY DIM="2" SPACE="2">
+
+        <VERTEX>
+            <V ID="0"> -1.0 -1.0 0.0 </V>
+            <V ID="1"> 0.0 -1.0 0.0 </V>
+            <V ID="2"> 1.0 -1.0 0.0 </V>
+            <V ID="3"> -1.0 0.0 0.0 </V>
+            <V ID="4"> 0.0 0.0 0.0 </V>
+            <V ID="5"> 1.0 0.0 0.0 </V>
+            <V ID="6"> -1.0 1.0 0.0 </V>
+            <V ID="7"> 0.0 1.0 0.0 </V>
+            <V ID="8"> 1.0 1.0 0.0 </V>
+        </VERTEX>
+
+        <EDGE>
+            <E ID="0"> 0 1  </E>
+            <E ID="1"> 1 2  </E>
+            <E ID="2"> 0 3  </E>
+            <E ID="3"> 1 4  </E>
+            <E ID="4"> 2 5  </E>
+            <E ID="5"> 3 4  </E>
+            <E ID="6"> 4 5  </E>
+            <E ID="7"> 6 3  </E>
+            <E ID="8"> 4 7  </E>
+            <E ID="9"> 5 8  </E>
+            <E ID="10"> 6 7  </E>
+            <E ID="11"> 7 8  </E>
+            <E ID="12"> 0 4  </E>
+            <E ID="13"> 1 5  </E>
+        </EDGE>
+
+        <ELEMENT>
+            <T ID="0"> 0 3 12 </T>
+            <T ID="1"> 2 12 5 </T>
+            <T ID="2"> 1 4 13 </T>
+            <T ID="3"> 3 13 6 </T>
+            <Q ID="4"> 5 8 10 7 </Q>
+            <Q ID="5"> 8 6 9 11 </Q>
+        </ELEMENT>
+
+        <COMPOSITE>
+            <C ID="0"> T[0-3] </C>
+            <C ID="1"> Q[4-5] </C>
+            <C ID="2"> E[2,7,4,9] </C>
+            <C ID="3"> E[0,1,10,11] </C>
+        </COMPOSITE>
+
+        <DOMAIN> C[0-1] </DOMAIN>
+
+    </GEOMETRY>
+
+    <EXPANSIONS>
+        <E COMPOSITE="C[0]" NUMMODES="6,7" FIELDS="u" BASISTYPE="Modified_A,Modified_B" NUMPOINTS="7,8" POINTSTYPE="GaussLobattoLegendre,GaussRadauMAlpha1Beta0"/>
+        <E COMPOSITE="C[1]" NUMMODES="6,7" FIELDS="u" BASISTYPE="Modified_A,Modified_A" NUMPOINTS="7,8" POINTSTYPE="GaussLobattoLegendre,GaussLobattoLegendre"/>
+    </EXPANSIONS>
+
+</NEKTAR>
+

tests/test_bwdtrans.cpp

@@ -77,7 +77,8 @@ BOOST_FIXTURE_TEST_CASE(bwdtrans, Line)
                 if (coeff == 0) {
                     *(x++) = 1.0;
                 }
-                else {
+                else
+                {
                     *(x++) = 0.0;
                 }
             }