diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8707cbcd382044632773d355cde24d26b61dea47..ada084f53d15d441b3deac7a722522f444e28d18 100644
--- a/CMakeLists.txt
+++ b/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)
diff --git a/Field.cpp b/Field.cpp
index d5d8d47e5a0824c1aa45c313996ce665cea339ff..0a161abd59e2a9842793ae3a3479691423b0080d 100644
--- a/Field.cpp
+++ b/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;
}
}
diff --git a/Field.hpp b/Field.hpp
index 843175d32dd43a93e7a2418107eb6f09437e6c6c..5e97e51bf20f02fb6c858133c0879cf4ed29f8ba 100644
--- a/Field.hpp
+++ b/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>>(
diff --git a/MemoryRegionCPU.hpp b/MemoryRegionCPU.hpp
index 877f27cf47d1a9afa7a70332954e374a089f66af..1c1e169c8c337a42a896054b7f6b2ee424a44b79 100644
--- a/MemoryRegionCPU.hpp
+++ b/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)
diff --git a/MemoryRegionCUDA.hpp b/MemoryRegionCUDA.hpp
index 2187341bee8d1e5033656e0746b0c576cd9fd40f..ff61995b92af25660b9e3533555dccf0c0fc3d76 100644
--- a/MemoryRegionCUDA.hpp
+++ b/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
diff --git a/Operators/BwdTrans/BwdTransImpl.cpp b/Operators/BwdTrans/BwdTransImpl.cpp
index d2290d59c497bbf205b527775b964f8bb516c4eb..9a92e3cde029bae54d33af694277a276f9e3a758 100644
--- a/Operators/BwdTrans/BwdTransImpl.cpp
+++ b/Operators/BwdTrans/BwdTransImpl.cpp
@@ -1,5 +1,5 @@
-#include "BwdTransStdMat.hpp"
#include "BwdTransMatFree.hpp"
+#include "BwdTransStdMat.hpp"
#include "BwdTransSumFac.hpp"
namespace Nektar::Operators::detail
diff --git a/Operators/IProductWRTBase/IProductWRTBaseImpl.cpp b/Operators/IProductWRTBase/IProductWRTBaseImpl.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..44364c06ddf91312ad7ae04d9fb192cd7c561b84
--- /dev/null
+++ b/Operators/IProductWRTBase/IProductWRTBaseImpl.cpp
@@ -0,0 +1,13 @@
+#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
diff --git a/Operators/IProductWRTBase/IProductWRTBaseStdMat.hpp b/Operators/IProductWRTBase/IProductWRTBaseStdMat.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6d8264455995dfb0067caf3267cfe3e86da291aa
--- /dev/null
+++ b/Operators/IProductWRTBase/IProductWRTBaseStdMat.hpp
@@ -0,0 +1,129 @@
+#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
diff --git a/Operators/OperatorIProductWRTBase.cpp b/Operators/OperatorIProductWRTBase.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4dcdbf1411b1911788743bc3355eed8f7d74c5aa
--- /dev/null
+++ b/Operators/OperatorIProductWRTBase.cpp
@@ -0,0 +1,9 @@
+#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
diff --git a/Operators/OperatorIProductWRTBase.hpp b/Operators/OperatorIProductWRTBase.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..fe6adfdd22a9d4c76dec75b2e9b6a085b7ef7abe
--- /dev/null
+++ b/Operators/OperatorIProductWRTBase.hpp
@@ -0,0 +1,54 @@
+#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
diff --git a/line.xml b/line.xml
new file mode 100644
index 0000000000000000000000000000000000000000..45b9d23682019010d27a7fe34b6db30275419824
--- /dev/null
+++ b/line.xml
@@ -0,0 +1,28 @@
+<?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>
diff --git a/main.cpp b/main.cpp
index 00ce6f150fc262bfb91dc0bdb750df8d5a055baf..94521fba49a392368c2ca4c15eb1a9225626fb11 100644
--- a/main.cpp
+++ b/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;
}
diff --git a/squareDeformed.xml b/squareDeformed.xml
new file mode 100644
index 0000000000000000000000000000000000000000..a5a0f6cad71ad4fc0bbc48261dca9a2d64ba5556
--- /dev/null
+++ b/squareDeformed.xml
@@ -0,0 +1,23 @@
+<?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>
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index b86f61ac7c1251219954dc086a262b38f7bd1f20..1a4f15028a80f74ef5c49515f64334c88b5d3599 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,23 +1,37 @@
-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}"
)
diff --git a/tests/init_fields.hpp b/tests/init_fields.hpp
index 0d44eb0c78b9d0a5abae54d507277c03c77bcee5..9aeb7f82f39bc72bb5219feff10ef56e3139e456 100644
--- a/tests/init_fields.hpp
+++ b/tests/init_fields.hpp
@@ -1,10 +1,12 @@
+#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;
};
diff --git a/tests/test_bwdtrans.cpp b/tests/test_bwdtrans.cpp
index 33f338d29f04b9c0225064903568dcfe20c3e665..5116addecf30f5e78ce15aa7b77a1f4725423c6e 100644
--- a/tests/test_bwdtrans.cpp
+++ b/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);
}
diff --git a/tests/test_ipwrtbase.cpp b/tests/test_ipwrtbase.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..459f0e24a1a651ba4c016bf3d496b17725048967
--- /dev/null
+++ b/tests/test_ipwrtbase.cpp
@@ -0,0 +1,67 @@
+#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);
+}