diff --git a/CMakeLists.txt b/CMakeLists.txt
index 43526d8534cefb67b8ae7dfcfb41f571845b98a1..a628e8608e034cf7bbde313da690e3b5c27b6bfc 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -47,6 +47,7 @@ set(SRC Field.cpp Operators/Operator.cpp
     Operators/OperatorFwdTrans.cpp Operators/FwdTrans/FwdTransImpl.cpp
     Operators/OperatorHelmSolve.cpp Operators/HelmSolve/HelmSolveImpl.cpp
     Operators/OperatorMatrix.cpp Operators/Matrix/MatrixImpl.cpp
+    Operators/OperatorAddTraceIntegral.cpp Operators/AddTraceIntegral/AddTraceIntegralImpl.cpp
     )
 
 if (NEKTAR_USE_CUDA AND NEKTAR_USE_SIMD)
diff --git a/Operators/AddTraceIntegral/AddTraceIntegralImpl.cpp b/Operators/AddTraceIntegral/AddTraceIntegralImpl.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fd19d79fbd3588e8f0f77a450c54ad48053b61d1
--- /dev/null
+++ b/Operators/AddTraceIntegral/AddTraceIntegralImpl.cpp
@@ -0,0 +1,13 @@
+#include "AddTraceIntegralStdMat.hpp"
+
+namespace Nektar::Operators::detail
+{
+
+// Add different AddTraceIntegral implementations to the factory.
+template <>
+std::string OperatorAddTraceIntegralImpl<double, ImplStdMat>::className =
+    GetOperatorFactory<double>().RegisterCreatorFunction(
+        "AddTraceIntegralStdMat",
+        OperatorAddTraceIntegralImpl<double, ImplStdMat>::instantiate, "...");
+
+} // namespace Nektar::Operators::detail
diff --git a/Operators/AddTraceIntegral/AddTraceIntegralStdMat.hpp b/Operators/AddTraceIntegral/AddTraceIntegralStdMat.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..eddba30a09486ff926dadfb11c7a610763ab387c
--- /dev/null
+++ b/Operators/AddTraceIntegral/AddTraceIntegralStdMat.hpp
@@ -0,0 +1,102 @@
+#pragma once
+
+#include "Operators/IProductWRTBase/IProductWRTBaseStdMat.hpp"
+#include "Operators/OperatorAddTraceIntegral.hpp"
+
+using namespace Nektar::MultiRegions;
+
+namespace Nektar::Operators::detail
+{
+
+// standard matrix implementation
+template <typename TData>
+class OperatorAddTraceIntegralImpl<TData, ImplStdMat>
+    : public OperatorAddTraceIntegral<TData>
+{
+public:
+    OperatorAddTraceIntegralImpl(
+        const MultiRegions::ExpListSharedPtr &expansionList)
+        : OperatorAddTraceIntegral<TData>(std::move(expansionList)),
+          m_wsp_trace(Field<TData, FieldState::Coeff>::create(
+              GetBlockAttributes(FieldState::Coeff, expansionList->GetTrace())))
+    {
+        // Initialise Trace
+        const MultiRegions::ExpListSharedPtr &traceExpansionList =
+            expansionList->GetTrace();
+        nTraceCoeffs = traceExpansionList->GetNcoeffs();
+        nFieldCoeffs = expansionList->GetNcoeffs();
+
+        // Get Trace-to-Element Map
+        m_locTraceToTraceMap = expansionList->GetLocTraceToTraceMap();
+
+        // Initialise IProductWRTBase operator
+        m_IProductWRTBaseOp = IProductWRTBase<>::create(
+            this->m_expansionList->GetTrace(), "StdMat");
+    }
+
+    void apply(Field<TData, FieldState::Phys> &in,
+               Field<TData, FieldState::Coeff> &out) override
+    {
+        // Step 1: Inner product for trace integral
+        m_IProductWRTBaseOp->apply(in, m_wsp_trace);
+
+        // Step 2: LEGACY Map Trace to Element
+        // TODO: Make this Field-only
+        Array<OneD, TData> inarray(nTraceCoeffs);
+        Array<OneD, TData> outarray(nFieldCoeffs, 0.0);
+
+        // Copy data from input field
+        // Field -> Array
+        auto *inarrptr = inarray.data();
+        auto *inptr    = m_wsp_trace.GetStorage().GetCPUPtr();
+        for (auto const &block : m_wsp_trace.GetBlocks())
+        {
+            auto nSize  = block.block_size;
+            auto nElmts = block.num_elements;
+            auto nmTot  = block.num_pts;
+
+            std::copy(inptr, inptr + nElmts * nmTot, inarrptr);
+
+            inarrptr += nElmts * nmTot;
+            inptr += nSize;
+        }
+
+        // The legacy map (that needs to be vectorised)
+        m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(inarray, outarray);
+
+        // Copy data to output field
+        // Array -> Field
+        auto *outarrptr = outarray.data();
+        auto *outptr    = out.GetStorage().GetCPUPtr();
+        for (auto const &block : out.GetBlocks())
+        {
+            auto nSize  = block.block_size;
+            auto nElmts = block.num_elements;
+            auto nmTot  = block.num_pts;
+
+            std::copy(outarrptr, outarrptr + nElmts * nmTot, outptr);
+
+            outarrptr += nElmts * nmTot;
+            outptr += nSize;
+        }
+    }
+
+    static std::unique_ptr<Operator<TData>> instantiate(
+        const MultiRegions::ExpListSharedPtr &expansionList)
+    {
+        return std::make_unique<
+            OperatorAddTraceIntegralImpl<TData, ImplStdMat>>(expansionList);
+    }
+
+    static std::string className;
+
+private:
+    Field<TData, FieldState::Coeff> m_wsp_trace;
+    std::shared_ptr<OperatorIProductWRTBase<TData>> m_IProductWRTBaseOp;
+    Nektar::MultiRegions::LocTraceToTraceMapSharedPtr m_locTraceToTraceMap;
+    Array<OneD, Array<OneD, int>> m_traceCoeffsToElmtMap;
+    Array<OneD, Array<OneD, int>> m_traceCoeffsToElmtSign;
+    int nTraceCoeffs;
+    int nFieldCoeffs;
+};
+} // namespace Nektar::Operators::detail
diff --git a/Operators/OperatorAddTraceIntegral.cpp b/Operators/OperatorAddTraceIntegral.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d32b2a9bf00834c8cc9ecdc763ddaf4cfeb0c50c
--- /dev/null
+++ b/Operators/OperatorAddTraceIntegral.cpp
@@ -0,0 +1,9 @@
+#include "OperatorAddTraceIntegral.hpp"
+
+namespace Nektar::Operators
+{
+
+template <> const std::string AddTraceIntegral<>::key = "AddTraceIntegral";
+template <> const std::string AddTraceIntegral<>::default_impl = "StdMat";
+
+} // namespace Nektar::Operators
diff --git a/Operators/OperatorAddTraceIntegral.hpp b/Operators/OperatorAddTraceIntegral.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..37130a722869c14f35c5f55989ccc4738c027125
--- /dev/null
+++ b/Operators/OperatorAddTraceIntegral.hpp
@@ -0,0 +1,56 @@
+#pragma once
+
+#include "Field.hpp"
+#include "Operator.hpp"
+
+namespace Nektar::Operators
+{
+
+// AddTraceIntegral base class
+// Defines the apply operator to enforce apply parameter types
+template <typename TData>
+class OperatorAddTraceIntegral : public Operator<TData>
+{
+public:
+    virtual ~OperatorAddTraceIntegral() = default;
+
+    OperatorAddTraceIntegral(
+        const MultiRegions::ExpListSharedPtr &expansionList)
+        : Operator<TData>(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 AddTraceIntegral
+template <typename TData = default_fp_type> struct AddTraceIntegral
+{
+    using class_name = OperatorAddTraceIntegral<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<AddTraceIntegral>(expansionList,
+                                                                  pKey);
+    }
+};
+
+namespace detail
+{
+// Template for AddTraceIntegral implementations
+template <typename TData, typename Op> class OperatorAddTraceIntegralImpl;
+} // namespace detail
+
+} // namespace Nektar::Operators
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index be911b26ad232d10c0abd982747ecf4467dd42f9..a37c18c314c20fe7740070e5321314bfd78cd464 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -340,6 +340,15 @@ IF (NEKTAR_USE_CUDA)
         -DTEST_PATH="${CMAKE_SOURCE_DIR}")
 ENDIF()
 
+add_executable(test_addtraceintegral test_addtraceintegral.cpp)
+target_link_libraries(test_addtraceintegral PRIVATE Operators)
+target_link_libraries(test_addtraceintegral PRIVATE ${NEKTAR++_LIBRARIES})
+target_link_libraries(test_addtraceintegral PRIVATE Boost::unit_test_framework)
+target_include_directories(test_addtraceintegral PRIVATE "${CMAKE_SOURCE_DIR}")
+target_include_directories(test_addtraceintegral PRIVATE ${NEKTAR++_INCLUDE_DIRS})
+target_compile_definitions(test_addtraceintegral PRIVATE
+    -DTEST_PATH="${CMAKE_SOURCE_DIR}")
+
 IF (NEKTAR_USE_CUDA)
     add_test(
       NAME CUDAKernels
@@ -573,3 +582,9 @@ IF (NEKTAR_USE_CUDA)
       WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}"
     )
 ENDIF()
+
+add_test(
+    NAME AddTraceIntegral
+  COMMAND test_addtraceintegral
+  WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}"
+)
diff --git a/tests/init_addtraceintegralfields.hpp b/tests/init_addtraceintegralfields.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..fc8943bbe2e67fbc92b55bf0516599e783f3bcb5
--- /dev/null
+++ b/tests/init_addtraceintegralfields.hpp
@@ -0,0 +1,198 @@
+#include "init_fields.hpp"
+
+using namespace std;
+using namespace Nektar::Operators;
+using namespace Nektar::LibUtilities;
+using namespace Nektar;
+
+class AddTraceIntegralField
+    : public InitFields<double, FieldState::Phys, FieldState::Coeff,
+                        MultiRegions::DisContField>
+{
+public:
+    AddTraceIntegralField()
+        : InitFields<double, FieldState::Phys, FieldState::Coeff,
+                     MultiRegions::DisContField>()
+    {
+    }
+
+    /*
+     *  Re-Initialise the input blocks based on the Trace-ExpList for this
+     * operator Delete previouisly defined fixt_in (also for CUDA) and re-define
+     * input based on TraceExpList
+     */
+    void ReConfigure(size_t nin = 1, size_t nout = 1)
+    {
+        if (fixt_in)
+        {
+            delete fixt_in;
+        }
+        const FieldState stateIn = FieldState::Phys;
+        auto blocks_in =
+            GetBlockAttributes(stateIn, fixt_explist->GetTrace(), vec_t::width);
+        auto f_in =
+            Field<double, stateIn>::create(blocks_in, nin, vec_t::alignment);
+        fixt_in = new Field<double, stateIn>(std::move(f_in));
+
+#ifdef NEKTAR_USE_CUDA
+        if (fixtcuda_in)
+        {
+            delete fixtcuda_in;
+        }
+        auto fcuda_in =
+            Field<double, stateIn>::template create<MemoryRegionCUDA>(blocks_in,
+                                                                      nin);
+        fixtcuda_in = new Field<double, stateIn>(std::move(fcuda_in));
+#endif
+    }
+
+    void SetTestCase(const std::vector<BlockAttributes> &blocks, double *inptr,
+                     bool padding = true)
+    {
+        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)
+                {
+                    *(inptr++) = phys;
+                }
+            }
+            if (padding)
+            {
+                for (size_t el = 0; el < block.num_padding_elements; ++el)
+                {
+                    for (size_t phys = 0; phys < block.num_pts; ++phys)
+                    {
+                        inptr++;
+                    }
+                }
+            }
+        }
+    }
+
+    void ExpectedSolution(const std::vector<BlockAttributes> &blocks,
+                          double *inptr)
+    {
+        auto fixt_explist_trace = fixt_explist->GetTrace();
+        Array<OneD, NekDouble> inTracephys(fixt_explist_trace->GetNpoints(),
+                                           0.0);
+        Array<OneD, NekDouble> outFieldcoeffs(fixt_explist->GetNcoeffs(), 0.0);
+
+        // Set test case
+        SetTestCase(blocks, inTracephys.get(), false);
+
+        // Calculate expected result from Nektar++
+        fixt_explist->AddTraceIntegral(inTracephys, outFieldcoeffs);
+
+        // Copy expected result from Array to fixt_expected
+        double *ptr = outFieldcoeffs.get();
+        for (auto const &block : blocks)
+        {
+            for (size_t el = 0; el < block.num_elements; ++el)
+            {
+                for (size_t coeff = 0; coeff < block.num_pts; ++coeff)
+                {
+                    (*inptr++) = (*ptr++);
+                }
+            }
+            for (size_t el = 0; el < block.num_padding_elements; ++el)
+            {
+                for (size_t coeff = 0; coeff < block.num_pts; ++coeff)
+                {
+                    inptr++;
+                }
+            }
+        }
+    }
+};
+
+class Seg : public AddTraceIntegralField
+{
+public:
+    Seg()
+    {
+        meshName = "run/line.xml";
+    }
+};
+
+class Quad : public AddTraceIntegralField
+{
+public:
+    Quad()
+    {
+        meshName = "run/square.xml";
+    }
+};
+
+class Tri : public AddTraceIntegralField
+{
+public:
+    Tri()
+    {
+        meshName = "run/tri.xml";
+    }
+};
+
+class SquareAllElements : public AddTraceIntegralField
+{
+public:
+    SquareAllElements()
+    {
+        meshName = "run/square_all_elements.xml";
+    }
+};
+
+class Hex : public AddTraceIntegralField
+{
+public:
+    Hex()
+    {
+        meshName = "run/hex.xml";
+    }
+};
+
+class Prism : public AddTraceIntegralField
+{
+public:
+    Prism()
+    {
+        meshName = "run/prism.xml";
+    }
+};
+
+class Pyr : public AddTraceIntegralField
+{
+public:
+    Pyr()
+    {
+        meshName = "run/pyr.xml";
+    }
+};
+
+class Tet : public AddTraceIntegralField
+{
+public:
+    Tet()
+    {
+        meshName = "run/tet.xml";
+    }
+};
+
+class CubePrismHex : public AddTraceIntegralField
+{
+public:
+    CubePrismHex()
+    {
+        meshName = "run/cube_prismhex.xml";
+    }
+};
+
+class CubeAllElements : public AddTraceIntegralField
+{
+public:
+    CubeAllElements()
+    {
+        meshName = "run/cube_all_elements.xml";
+    }
+};
diff --git a/tests/init_fields.hpp b/tests/init_fields.hpp
index 381c33380a7acecca65e61061a376b82cbf8f729..97a752a98263b7dbecf0ec5f8d9516abb320dd6b 100644
--- a/tests/init_fields.hpp
+++ b/tests/init_fields.hpp
@@ -6,6 +6,7 @@
 
 #include "Field.hpp"
 #include <MultiRegions/ContField.h>
+#include <MultiRegions/DisContField.h>
 #include <MultiRegions/ExpList.h>
 
 #ifdef NEKTAR_USE_CUDA
@@ -96,6 +97,14 @@ public:
                     Collections::eNoCollection);
         }
 
+        else if constexpr (std::is_same_v<TExpList, MultiRegions::DisContField>)
+        {
+            fixt_explist =
+                MemoryManager<MultiRegions::DisContField>::AllocateSharedPtr(
+                    session, graph, "u", true, true,
+                    Collections::eNoCollection);
+        }
+
         if constexpr (std::is_same_v<TExpList, MultiRegions::ExpList>)
         {
             fixt_explist =
diff --git a/tests/test_addtraceintegral.cpp b/tests/test_addtraceintegral.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8363487ddaea921733029dae6b36540d5a6e0f9d
--- /dev/null
+++ b/tests/test_addtraceintegral.cpp
@@ -0,0 +1,193 @@
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE TestAddTraceIntegral
+#include <boost/test/tools/output_test_stream.hpp>
+#include <boost/test/unit_test.hpp>
+
+#include <iostream>
+#include <memory>
+
+#include "Operators/OperatorAddTraceIntegral.hpp"
+#include "init_addtraceintegralfields.hpp"
+
+BOOST_AUTO_TEST_SUITE(TestAddTraceIntegral)
+
+using namespace std;
+using namespace Nektar::Operators;
+using namespace Nektar::LibUtilities;
+using namespace Nektar;
+
+/*
+ * Currently fails in GetBlockAttributes
+ * GEometry is not initialised for Expansion(0)
+ * Possibly, because trace is not working for 1D expansions
+BOOST_FIXTURE_TEST_CASE(addtraceintegral_seg, Seg)
+{
+    Configure();
+    ReConfigure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    AddTraceIntegral<>::create(fixt_explist, "StdMat")->apply(*fixt_in,
+*fixt_out); ExpectedSolution(fixt_expected->GetBlocks(),
+                     fixt_expected->GetStorage().GetCPUPtr());
+    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+*/
+
+BOOST_FIXTURE_TEST_CASE(addtraceintegral_quad, Quad)
+{
+    Configure();
+    ReConfigure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    AddTraceIntegral<>::create(fixt_explist, "StdMat")
+        ->apply(*fixt_in, *fixt_out);
+    ExpectedSolution(fixt_expected->GetBlocks(),
+                     fixt_expected->GetStorage().GetCPUPtr());
+    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_FIXTURE_TEST_CASE(addtraceintegral_tri, Tri)
+{
+    Configure();
+    ReConfigure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    AddTraceIntegral<>::create(fixt_explist, "StdMat")
+        ->apply(*fixt_in, *fixt_out);
+    ExpectedSolution(fixt_expected->GetBlocks(),
+                     fixt_expected->GetStorage().GetCPUPtr());
+    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_FIXTURE_TEST_CASE(addtraceintegral_square_all_elements, SquareAllElements)
+{
+    Configure();
+    ReConfigure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    AddTraceIntegral<>::create(fixt_explist, "StdMat")
+        ->apply(*fixt_in, *fixt_out);
+    ExpectedSolution(fixt_expected->GetBlocks(),
+                     fixt_expected->GetStorage().GetCPUPtr());
+    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_FIXTURE_TEST_CASE(addtraceintegral_hex, Hex)
+{
+    Configure();
+    ReConfigure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    AddTraceIntegral<>::create(fixt_explist, "StdMat")
+        ->apply(*fixt_in, *fixt_out);
+    ExpectedSolution(fixt_expected->GetBlocks(),
+                     fixt_expected->GetStorage().GetCPUPtr());
+    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_FIXTURE_TEST_CASE(addtraceintegral_prism, Prism)
+{
+    Configure();
+    ReConfigure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    AddTraceIntegral<>::create(fixt_explist, "StdMat")
+        ->apply(*fixt_in, *fixt_out);
+    ExpectedSolution(fixt_expected->GetBlocks(),
+                     fixt_expected->GetStorage().GetCPUPtr());
+    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_FIXTURE_TEST_CASE(addtraceintegral_pyr, Pyr)
+{
+    Configure();
+    ReConfigure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    AddTraceIntegral<>::create(fixt_explist, "StdMat")
+        ->apply(*fixt_in, *fixt_out);
+    ExpectedSolution(fixt_expected->GetBlocks(),
+                     fixt_expected->GetStorage().GetCPUPtr());
+    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_FIXTURE_TEST_CASE(addtraceintegral_tet, Tet)
+{
+    Configure();
+    ReConfigure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    AddTraceIntegral<>::create(fixt_explist, "StdMat")
+        ->apply(*fixt_in, *fixt_out);
+    ExpectedSolution(fixt_expected->GetBlocks(),
+                     fixt_expected->GetStorage().GetCPUPtr());
+    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_FIXTURE_TEST_CASE(addtraceintegral_cube_prism_hex, CubePrismHex)
+{
+    Configure();
+    ReConfigure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    AddTraceIntegral<>::create(fixt_explist, "StdMat")
+        ->apply(*fixt_in, *fixt_out);
+    ExpectedSolution(fixt_expected->GetBlocks(),
+                     fixt_expected->GetStorage().GetCPUPtr());
+    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_FIXTURE_TEST_CASE(addtraceintegral_cube_all_elements, CubeAllElements)
+{
+    Configure();
+    ReConfigure();
+    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
+    AddTraceIntegral<>::create(fixt_explist, "StdMat")
+        ->apply(*fixt_in, *fixt_out);
+    ExpectedSolution(fixt_expected->GetBlocks(),
+                     fixt_expected->GetStorage().GetCPUPtr());
+    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    boost::test_tools::output_test_stream output;
+    {
+        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
+                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+    }
+}
+
+BOOST_AUTO_TEST_SUITE_END()