From 3600a7c795e7c17656dc2cb101e5767c976ceda4 Mon Sep 17 00:00:00 2001 From: Chris Cantwell Date: Tue, 22 Mar 2022 22:12:12 +0000 Subject: [PATCH 1/6] Initial implementation of TensorRegion and TensorStorage classes. --- library/CMakeLists.txt | 2 +- library/Demos/TensorRegions/CMakeLists.txt | 5 + library/Demos/TensorRegions/Domain.cpp | 44 +++++ library/TensorRegions/CMakeLists.txt | 25 +++ library/TensorRegions/TensorRegion.cpp | 159 ++++++++++++++++ library/TensorRegions/TensorRegion.h | 69 +++++++ library/TensorRegions/TensorRegions.hpp | 49 +++++ library/TensorRegions/TensorRegionsDeclspec.h | 46 +++++ library/TensorRegions/TensorStorage.hpp | 171 ++++++++++++++++++ 9 files changed, 569 insertions(+), 1 deletion(-) create mode 100644 library/Demos/TensorRegions/CMakeLists.txt create mode 100644 library/Demos/TensorRegions/Domain.cpp create mode 100644 library/TensorRegions/CMakeLists.txt create mode 100644 library/TensorRegions/TensorRegion.cpp create mode 100644 library/TensorRegions/TensorRegion.h create mode 100644 library/TensorRegions/TensorRegions.hpp create mode 100644 library/TensorRegions/TensorRegionsDeclspec.h create mode 100644 library/TensorRegions/TensorStorage.hpp diff --git a/library/CMakeLists.txt b/library/CMakeLists.txt index 257b1d5cc..e75979415 100644 --- a/library/CMakeLists.txt +++ b/library/CMakeLists.txt @@ -1,6 +1,6 @@ # Main library sub-directories, required by all of Nektar++. SUBDIRS(LibUtilities LocalRegions SpatialDomains StdRegions Collections - MultiRegions MatrixFreeOps SolverUtils GlobalMapping FieldUtils NekMesh) + MultiRegions TensorRegions MatrixFreeOps SolverUtils GlobalMapping FieldUtils NekMesh) INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR}/library) IF (NEKTAR_BUILD_UNIT_TESTS) diff --git a/library/Demos/TensorRegions/CMakeLists.txt b/library/Demos/TensorRegions/CMakeLists.txt new file mode 100644 index 000000000..f2e5b30cc --- /dev/null +++ b/library/Demos/TensorRegions/CMakeLists.txt @@ -0,0 +1,5 @@ +#ADD_NEKTAR_EXECUTABLE(Project +# COMPONENT demos DEPENDS TensorRegions SOURCES Project.cpp) + +ADD_NEKTAR_EXECUTABLE(Domain + COMPONENT demos DEPENDS TensorRegions SOURCES Domain.cpp) \ No newline at end of file diff --git a/library/Demos/TensorRegions/Domain.cpp b/library/Demos/TensorRegions/Domain.cpp new file mode 100644 index 000000000..ac9e0b34a --- /dev/null +++ b/library/Demos/TensorRegions/Domain.cpp @@ -0,0 +1,44 @@ +#include +#include +#include +#include +#include +#include + +using namespace std; +using namespace Nektar; + +// compile using Builds/Demos/StdRegions -> make DEBUG=1 ProjectLoc1D + +// This routine projects a polynomial which has energy in all mdoes of +// the expansions and report an error. + +int main(int argc, char *argv[]) +{ + LibUtilities::SessionReaderSharedPtr vSession + = LibUtilities::SessionReader::CreateInstance(argc, argv); + + // read in mesh + SpatialDomains::MeshGraphSharedPtr graph1D = SpatialDomains::MeshGraph::Read(vSession); + + // Read the geometry and the expansion information + size_t nDomains = graph1D->GetDomain().size(); + const std::vector domain = graph1D->GetDomain(); + SpatialDomains::BoundaryConditions Allbcs(vSession, graph1D); + bool SetToOneSpaceDimension = true; + + std::vector exps(nDomains); + for (unsigned int i = 0; i < nDomains; ++i) + { + exps[i] = + MemoryManager::AllocateSharedPtr( + vSession, graph1D, domain[i], Allbcs, + vSession->GetVariable(0), SetToOneSpaceDimension); + } + + TensorRegions::TensorRegion tr(exps); + + TensorRegions::TensorRegionStorageSharedPtr x = tr.AllocateStorage(); + + return 0; +} \ No newline at end of file diff --git a/library/TensorRegions/CMakeLists.txt b/library/TensorRegions/CMakeLists.txt new file mode 100644 index 000000000..43461033a --- /dev/null +++ b/library/TensorRegions/CMakeLists.txt @@ -0,0 +1,25 @@ +SET(TENSOR_REGIONS_SOURCES +TensorRegion.cpp +) + +SET(TENSOR_REGIONS_HEADERS +TensorRegion.h +) + +ADD_DEFINITIONS(-DTENSOR_REGIONS_EXPORTS) + +ADD_NEKTAR_LIBRARY(TensorRegions + SOURCES ${TENSOR_REGIONS_SOURCES} + HEADERS ${TENSOR_REGIONS_HEADERS} + DEPENDS MultiRegions + SUMMARY "Nektar++ TensorRegions library" + DESCRIPTION "This library provides global expansions on tensor domains.") + +# IF (NEKTAR_BUILD_PYTHON) +# SUBDIRS(Python) +# ENDIF() + +INSTALL(DIRECTORY ./ DESTINATION ${NEKTAR_INCLUDE_DIR}/TensorRegions COMPONENT dev FILES_MATCHING PATTERN "*.h" PATTERN "*.hpp") + + + diff --git a/library/TensorRegions/TensorRegion.cpp b/library/TensorRegions/TensorRegion.cpp new file mode 100644 index 000000000..058b87d65 --- /dev/null +++ b/library/TensorRegions/TensorRegion.cpp @@ -0,0 +1,159 @@ +#include +#include + +namespace Nektar +{ +namespace TensorRegions +{ + +/** + * @class TensorRegion + * + * A TensorRegion represents a high-dimensional multi-region finite element + * expansion. It is constructed through a tensor-product of two lower- + * dimensional MultiRegions expansions, defined by two independent meshes. + * Operators are constructed by applying the sub-factorisation approach. + */ + +/** + * @brief Construct a new TensorRegion object from one or more MultiRegions + * objects. + * + * @param subspaces Vector of constituent subspaces. + */ +TensorRegion::TensorRegion(std::vector subspaces) +{ + ASSERTL0(subspaces.size() <= 2, "TensorRegion currently only supports 2 subspaces.") + m_exp = subspaces; + m_numPoints = 1; + m_numCoeffs = 1; + for (auto &e : m_exp) + { + m_numPoints *= e->GetNpoints(); + m_numCoeffs *= e->GetNcoeffs(); + } +} + +/** + * @brief Construct a new TensorRegion object from an existing object. + * + * @param pSrc Existing TensorRegion object. + */ +TensorRegion::TensorRegion(const TensorRegion& pSrc) + : m_exp(pSrc.m_exp), + m_numPoints(pSrc.m_numPoints), + m_numCoeffs(pSrc.m_numCoeffs) +{ + // Nothing to do +} + +/** + * @brief Destroy the TensorRegion object + * + */ +TensorRegion::~TensorRegion() +{ + +} + +TensorRegionStorageSharedPtr TensorRegion::AllocateStorage() +{ + std::vector sizes; + for (int i = 0; i < m_exp.size(); ++i) + { + sizes.push_back(m_exp[i]->GetNpoints()); + } + return std::make_shared(sizes); +} + +/** + * @brief Access a constituent MultiRegions object. + * + * @param index The index of the MultiRegion object to access. + * @return MultiRegions::ExpListSharedPtr + */ +MultiRegions::ExpListSharedPtr TensorRegion::GetExpList(int index) +{ + ASSERTL0(index >= m_exp.size(), "Index out of range."); + + return m_exp[index]; +} + + +size_t TensorRegion::GetNumPoints() +{ + return m_numPoints; +} + +size_t TensorRegion::GetNumCoeffs() +{ + return m_numCoeffs; +} + +void TensorRegion::BwdTrans( + const TensorStorage& inarray, + TensorStorage& outarray) +{ + boost::ignore_unused(inarray, outarray); + + // const int nm0 = m_exp[0]->GetNcoeffs(); + // const int nm1 = m_exp[1]->GetNcoeffs(); + // const int np0 = m_exp[0]->GetNpoints(); + // const int np1 = m_exp[1]->GetNpoints(); + + // ASSERTL0(inarray.size() >= nm0*nm1, "Input array is of incorrect size."); + // ASSERTL0(outarray.size() >= np0*np1, "Output array is of incorrect size."); + + // Array x; + + // // Apply first dimension bwdtrans + // for (int i = 0; i < nm1; ++i) { + // m_exp[0]->BwdTrans(inarray + i*nm0, x = outarray + i*np0); + // } + + // // Transpose data + // Array tmp(np0*nm1); + // for (int j = 0; j < nm1; ++j) { + // for (int i = 0; i < np0; ++i) { + // tmp[i*nm1 + j] = outarray[j*np0 + i]; + // } + // } + + // // Apply second dimension bwdtrans + // for (int j = 0; j < np0; ++j) { + // m_exp[1]->BwdTrans(tmp + j*np0, x = outarray + j*np1); + // } +} + +void TensorRegion::IProductWRTBase( + const TensorStorage& inarray, + TensorStorage& outarray) +{ + boost::ignore_unused(inarray, outarray); +} + +void TensorRegion::IProductWRTDerivBase( + unsigned int dir, + const TensorStorage& inarray, + TensorStorage& outarray) +{ + boost::ignore_unused(dir, inarray, outarray); +} + +void TensorRegion::PhysDeriv( + unsigned int dir, + const TensorStorage& inarray, + TensorStorage& outarray) +{ + boost::ignore_unused(dir, inarray, outarray); +} + +void TensorRegion::FwdTrans( + const TensorStorage& inarray, + TensorStorage& outarray) +{ + boost::ignore_unused(inarray, outarray); +} + +} +} \ No newline at end of file diff --git a/library/TensorRegions/TensorRegion.h b/library/TensorRegions/TensorRegion.h new file mode 100644 index 000000000..db3e96f0c --- /dev/null +++ b/library/TensorRegions/TensorRegion.h @@ -0,0 +1,69 @@ +#ifndef NEKTAR_LIBS_TENSORREGIONS_TENSORREGION_H +#define NEKTAR_LIBS_TENSORREGIONS_TENSORREGION_H + +#include +#include + +#include + +namespace Nektar +{ +namespace TensorRegions +{ + +using TensorRegionStorage = TensorStorage; +using TensorRegionStorageSharedPtr = std::shared_ptr; + +/// High-dimensional region, comprising of a tensor-product of two MultiRegions. +class TensorRegion +{ + public: + TensorRegion(std::vector exp); + TensorRegion(const TensorRegion& pSrc); + ~TensorRegion(); + + TensorRegionStorageSharedPtr AllocateStorage(); + + MultiRegions::ExpListSharedPtr GetExpList(int index); + + size_t GetNumPoints(); + size_t GetNumCoeffs(); + + void BwdTrans( + const TensorStorage& inarray, + TensorStorage& outarray); + + void IProductWRTBase( + const TensorStorage& inarray, + TensorStorage& outarray); + + void IProductWRTDerivBase( + unsigned int dir, + const TensorStorage& inarray, + TensorStorage& outarray); + + void PhysDeriv( + unsigned int dir, + const TensorStorage& inarray, + TensorStorage& outarray); + + void FwdTrans( + const TensorStorage& inarray, + TensorStorage& outarray); + + void GetCoords(Array> coords); + + NekDouble PhysEvaluate(Array coords); + + private: + std::vector m_exp; + size_t m_numPoints; + size_t m_numCoeffs; +}; + +using TensorRegionSharedPtr = std::shared_ptr; + +} +} + +#endif \ No newline at end of file diff --git a/library/TensorRegions/TensorRegions.hpp b/library/TensorRegions/TensorRegions.hpp new file mode 100644 index 000000000..4bd573d3d --- /dev/null +++ b/library/TensorRegions/TensorRegions.hpp @@ -0,0 +1,49 @@ +/////////////////////////////////////////////////////////////////////////////// +// +// File TensorRegsions.hpp +// +// For more information, please see: http://www.nektar.info +// +// The MIT License +// +// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA), +// Department of Aeronautics, Imperial College London (UK), and Scientific +// Computing and Imaging Institute, University of Utah (USA). +// +// Permission is hereby granted, free of charge, to any person obtaining a +// copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +// DEALINGS IN THE SOFTWARE. +// +// Description: TensorRegions overall header +// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef TENSORREGIONS_H +#define TENSORREGIONS_H + +#include + +namespace Nektar +{ +namespace TensorRegions +{ + + +}// end of namespace +}// end of namespace + +#endif diff --git a/library/TensorRegions/TensorRegionsDeclspec.h b/library/TensorRegions/TensorRegionsDeclspec.h new file mode 100644 index 000000000..4e2eaae2f --- /dev/null +++ b/library/TensorRegions/TensorRegionsDeclspec.h @@ -0,0 +1,46 @@ +/////////////////////////////////////////////////////////////////////////////// +// +// For more information, please see: http://www.nektar.info +// +// The MIT License +// +// Copyright (c) 2006 Scientific Computing and Imaging Institute, +// University of Utah (USA) and Department of Aeronautics, Imperial +// College London (UK). +// +// Permission is hereby granted, free of charge, to any person obtaining a +// copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +// DEALINGS IN THE SOFTWARE. +// +// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef NEKTAR_TENSOR_REGIONS_DECLSPEC_H +#define NEKTAR_TENSOR_REGIONS_DECLSPEC_H + +#if defined(_MSC_VER) + #ifdef TENSOR_REGIONS_EXPORTS + #define TENSOR_REGIONS_EXPORT _declspec(dllexport) + #else + #define TENSOR_REGIONS_EXPORT _declspec(dllimport) + #endif +#else + #define TENSOR_REGIONS_EXPORT +#endif + +#endif //NEKTAR_TENSOR_REGIONS_DECLSPEC_H + diff --git a/library/TensorRegions/TensorStorage.hpp b/library/TensorRegions/TensorStorage.hpp new file mode 100644 index 000000000..20d7632cb --- /dev/null +++ b/library/TensorRegions/TensorStorage.hpp @@ -0,0 +1,171 @@ +#ifndef NEKTAR_LIBS_TENSORREGIONS_TENSORSTORAGE_HPP +#define NEKTAR_LIBS_TENSORREGIONS_TENSORSTORAGE_HPP + +#include + +namespace Nektar +{ +namespace TensorRegions +{ + +/** + * @brief Storage container for TensorRegions solution. + * + * Data is stored fastest in the lowest dimension. + * + * @tparam DataType The type of data stored by the container. + */ +template +class TensorStorage +{ + public: + /** + * @brief Describes whether storage represents physical values + * or spectral/hp element coefficients. + * + */ + enum State { + ePhys, /// Physical values + eCoeff /// Spectral/hp element coefficients + }; + + /** + * @brief Construct a new TensorStorage object with two given sizes. + * + * @param size1 First component size. + * @param size2 Second component size. + * @param val Initial value for all elements. + */ + TensorStorage(const size_t size1, const size_t size2, const DataType val = 0.0) + : m_data(Array(size1 * size2, val)), + m_sizes({size1, size2}) + { + ASSERTL0(size1 > 0, "Size1 must be greater than zero."); + ASSERTL0(size2 > 0, "Size2 must be greater than zero."); + // Nothing to do + } + + /** + * @brief Construct a new TensorStorage object with two given sizes. + * + * @param sizes Vector of arbitrary number of sizes > 0 + * @param val Initial vlaue for all elements. + */ + TensorStorage(const std::vector sizes, const DataType val = 0.0) + : m_sizes(sizes) + { + ASSERTL1(!sizes.empty(), "Must have at least one size."); + int size = 1; + for (auto &s : sizes) + { + size *= s; + } + m_data = Array(size, val); + } + + /** + * @brief Construct a new TensorStorage object from an existing object. + * + * @param pSrc Existing TensorStorage object. + */ + TensorStorage(const TensorStorage& pSrc) + : m_data(pSrc.m_data), + m_sizes(pSrc.m_sizes) + { + // Nothing to do + } + + /** + * @brief Construct a new TensorStorage object from an existing object + * (move semantics). + * + * @param pSrc Existing temporary TensorStorage object. + */ + TensorStorage(TensorStorage&& pSrc) + : Array(std::move(pSrc)), + m_sizes(std::move(pSrc.m_sizes)) + { + // Nothing to do + } + + /** + * @brief Destroy a TensorRegions object. + * + */ + virtual ~TensorStorage() + { + + } + + /** + * @brief Returns the size of the TensorStorage object. + * + * @param idx Index of sub-component to request size. + * @return size_t + */ + size_t size(const size_t idx) + { + ASSERTL1(idx < m_sizes.size(), + "Requested index out of range."); + return m_sizes[idx]; + } + + /** + * @brief Returns whether the storage represents physical solution + * or spectral/hp element coefficients. + * + * @return enum State + */ + enum State getState() + { + return m_state; + } + + /** + * @brief Sets whether storage is in physical space or coefficient space. + * + */ + void setState(enum State state = ePhys) + { + m_state = state; + } + + /** + * @brief + * + */ + Array slice(unsigned int dim, std::vector coord) + { + ASSERTL0(dim < m_sizes.size(), "dim larger than number of dimensions."); + ASSERTL0(coord.size() == m_sizes.size(), "coord dimension must match dimension of Tensor."); + + Array x(m_sizes[dim]); + size_t stride = 1; + size_t offset = 0; + for (unsigned int i = 0; i < dim; ++i) + { + offset += coord[i] * stride; + stride *= m_sizes[i]; + } + + for (size_t i = 0; i < m_sizes[dim]; ++i) { + x[i] = m_data[offset + i*stride]; + } + return std::move(x); + } + // slice + + private: + Array m_data; /// Storage for actual data. + std::vector m_sizes; /// Sizes for each sub-component of the tensor. + enum TensorStorage::State m_state = ePhys; /// Flag indicating if data is in coeff or physical state. + +}; + +/// A shared pointer to a TensorStorage object. +template using TensorStorageSharedPtr = std::shared_ptr>; + +} +} + +#endif \ No newline at end of file -- GitLab From d22b08d80e3dddcaad8f6029194f201a9a638a49 Mon Sep 17 00:00:00 2001 From: Chris Cantwell Date: Tue, 22 Mar 2022 22:13:09 +0000 Subject: [PATCH 2/6] Initial TensorRegions Projection demo. --- library/Demos/TensorRegions/Project.cpp | 105 ++++++++++++++++++++++++ 1 file changed, 105 insertions(+) create mode 100644 library/Demos/TensorRegions/Project.cpp diff --git a/library/Demos/TensorRegions/Project.cpp b/library/Demos/TensorRegions/Project.cpp new file mode 100644 index 000000000..08f4ede05 --- /dev/null +++ b/library/Demos/TensorRegions/Project.cpp @@ -0,0 +1,105 @@ +#include +#include +#include +#include +#include + +using namespace std; +using namespace Nektar; + +// compile using Builds/Demos/StdRegions -> make DEBUG=1 ProjectLoc1D + +// This routine projects a polynomial which has energy in all mdoes of +// the expansions and report an error. + +int main(int argc, char *argv[]) +{ + LibUtilities::SessionReaderSharedPtr vSession + = LibUtilities::SessionReader::CreateInstance(argc, argv); + + TensorRegions::TensorStorageSharedPtr Exp,Sol; + int i,j; + int nq; + int coordim; + Array sol; + Array xc0,xc1,xc2; + + // read in mesh + SpatialDomains::MeshGraphSharedPtr graph1D = SpatialDomains::MeshGraph::Read(vSession); + + // Define Expansion + const SpatialDomains::ExpansionInfoMap &expansions = graph1D->GetExpansionInfos(); + LibUtilities::BasisKey bkey0 = expansions.begin()->second->m_basisKeyVector[0]; + int nmodes = bkey0.GetNumModes(); + + Exp = MemoryManager::AllocateSharedPtr(vSession, + graph1D); + + //---------------------------------------------- + // Define solution to be projected + coordim = Exp->GetCoordim(0); + nq = Exp->GetTotPoints(); + + // define coordinates and solution + sol = Array(nq); + + xc0 = Array(nq); + xc1 = Array(nq); + xc2 = Array(nq); + + switch(coordim) + { + case 1: + Exp->GetCoords(xc0); + Vmath::Zero(nq,&xc1[0],1); + Vmath::Zero(nq,&xc2[0],1); + break; + case 2: + Exp->GetCoords(xc0,xc1); + Vmath::Zero(nq,&xc2[0],1); + break; + case 3: + Exp->GetCoords(xc0,xc1,xc2); + break; + } + + for(i = 0; i < nq; ++i) + { + sol[i] = 0.0; + for(j = 0; j < nmodes; ++j) + { + sol[i] += pow(xc0[i],j); + sol[i] += pow(xc1[i],j); + sol[i] += pow(xc2[i],j); + } + } + + //--------------------------------------------- + // Set up ExpList containing the solution + Sol = MemoryManager::AllocateSharedPtr(*Exp); + Sol->SetPhys(sol); + //--------------------------------------------- + + //--------------------------------------------- + // Project onto Expansion + Exp->FwdTrans(Sol->GetPhys(), Exp->UpdateCoeffs()); + //--------------------------------------------- + + //------------------------------------------- + // Backward Transform Solution to get projected values + Exp->BwdTrans(Exp->GetCoeffs(), Exp->UpdatePhys()); + //------------------------------------------- + + //-------------------------------------------- + // Calculate L_inf error + if (vSession->GetComm()->GetRank() == 0) + { + cout << "L infinity error: " << Exp->Linf(Sol->GetPhys()) << endl; + cout << "L 2 error: " << Exp->L2 (Sol->GetPhys()) << endl; + } + //-------------------------------------------- + + vSession->Finalise(); + + return 0; +} -- GitLab From 9a3f4eeabbaa9dc31985c28b6139ceb76901575b Mon Sep 17 00:00:00 2001 From: Chris Cantwell Date: Mon, 28 Mar 2022 22:21:37 +0100 Subject: [PATCH 3/6] Extended Demo code, implemented View sub class. --- docs/doxygen/Doxyfile.in | 1 + library/Demos/CMakeLists.txt | 2 +- library/Demos/TensorRegions/Domain.cpp | 20 ++- library/TensorRegions/TensorRegion.cpp | 7 +- library/TensorRegions/TensorStorage.hpp | 226 ++++++++++++++++++++---- 5 files changed, 217 insertions(+), 39 deletions(-) diff --git a/docs/doxygen/Doxyfile.in b/docs/doxygen/Doxyfile.in index 47bcaa55e..90804118a 100644 --- a/docs/doxygen/Doxyfile.in +++ b/docs/doxygen/Doxyfile.in @@ -760,6 +760,7 @@ INPUT = @CMAKE_SOURCE_DIR@/docs/doxygen/ \ @CMAKE_SOURCE_DIR@/library/Collections/ \ @CMAKE_SOURCE_DIR@/library/LocalRegions/ \ @CMAKE_SOURCE_DIR@/library/MultiRegions/ \ + @CMAKE_SOURCE_DIR@/library/TensorRegions/ \ @CMAKE_SOURCE_DIR@/library/GlobalMapping/ \ @CMAKE_SOURCE_DIR@/library/FieldUtils/ \ @CMAKE_SOURCE_DIR@/library/SolverUtils/ \ diff --git a/library/Demos/CMakeLists.txt b/library/Demos/CMakeLists.txt index 7244e3d8c..ead144bfc 100644 --- a/library/Demos/CMakeLists.txt +++ b/library/Demos/CMakeLists.txt @@ -1,4 +1,4 @@ -SUBDIRS(LibUtilities StdRegions SpatialDomains LocalRegions Collections MultiRegions MatrixFreeOps) +SUBDIRS(LibUtilities StdRegions SpatialDomains LocalRegions Collections MultiRegions TensorRegions MatrixFreeOps) IF(NEKTAR_BUILD_PYTHON) SUBDIRS(Python) diff --git a/library/Demos/TensorRegions/Domain.cpp b/library/Demos/TensorRegions/Domain.cpp index ac9e0b34a..139651303 100644 --- a/library/Demos/TensorRegions/Domain.cpp +++ b/library/Demos/TensorRegions/Domain.cpp @@ -38,7 +38,25 @@ int main(int argc, char *argv[]) TensorRegions::TensorRegion tr(exps); - TensorRegions::TensorRegionStorageSharedPtr x = tr.AllocateStorage(); + TensorRegions::TensorRegionStorageSharedPtr st = tr.AllocateStorage(); + cout << "Total tensor storage size: " << st->size() << endl; + for (size_t i = 0; i < st->dims(); ++i) + { + cout << "Dim " << i << " has size " << st->size(i) << " points." << endl; + } + + TensorRegions::TensorRegionStorage::View st0 = st->slice(0, {0,0}); + TensorRegions::TensorRegionStorage::View st1 = st->slice(1, {0,0}); + cout << "Size of slice 0 is " << st0.size() << endl; + for (int i = 0; i < st1.size(); ++i) { + st1[i] = i; + } + Array data0 = st0.extract(); + Array data1 = st1.extract(); + cout << "Size of slice 0 data is " << data0.size() << endl; + for (int i = 0; i < st->size(); ++i) { + cout << (*st)[i] << endl; + } return 0; } \ No newline at end of file diff --git a/library/TensorRegions/TensorRegion.cpp b/library/TensorRegions/TensorRegion.cpp index 058b87d65..2d0cc6c6c 100644 --- a/library/TensorRegions/TensorRegion.cpp +++ b/library/TensorRegions/TensorRegion.cpp @@ -56,6 +56,11 @@ TensorRegion::~TensorRegion() } +/** + * @brief Allocates a compatible TensorStorage object. + * + * @return TensorRegionStorageSharedPtr + */ TensorRegionStorageSharedPtr TensorRegion::AllocateStorage() { std::vector sizes; @@ -63,7 +68,7 @@ TensorRegionStorageSharedPtr TensorRegion::AllocateStorage() { sizes.push_back(m_exp[i]->GetNpoints()); } - return std::make_shared(sizes); + return TensorRegionStorage::create(sizes); } /** diff --git a/library/TensorRegions/TensorStorage.hpp b/library/TensorRegions/TensorStorage.hpp index 20d7632cb..7235f86ba 100644 --- a/library/TensorRegions/TensorStorage.hpp +++ b/library/TensorRegions/TensorStorage.hpp @@ -1,6 +1,8 @@ #ifndef NEKTAR_LIBS_TENSORREGIONS_TENSORSTORAGE_HPP #define NEKTAR_LIBS_TENSORREGIONS_TENSORSTORAGE_HPP +#include + #include namespace Nektar @@ -8,6 +10,7 @@ namespace Nektar namespace TensorRegions { + /** * @brief Storage container for TensorRegions solution. * @@ -16,7 +19,7 @@ namespace TensorRegions * @tparam DataType The type of data stored by the container. */ template -class TensorStorage +class TensorStorage : public std::enable_shared_from_this> { public: /** @@ -30,37 +33,119 @@ class TensorStorage }; /** - * @brief Construct a new TensorStorage object with two given sizes. + * @brief Represents a one-dimensional view onto a TensorStorage object. * - * @param size1 First component size. - * @param size2 Second component size. - * @param val Initial value for all elements. + * @tparam DataType */ - TensorStorage(const size_t size1, const size_t size2, const DataType val = 0.0) - : m_data(Array(size1 * size2, val)), - m_sizes({size1, size2}) + class View + { + public: + /** + * @brief Construct a new TensorStorage::View object + * + * @param ts Base TensorStorage object. + * @param dim The dimension to be viewed. + * @param coord Coordinate of a point in the required view. + */ + View(std::shared_ptr> ts, unsigned int dim, std::vector coord) + : m_ts(ts), m_dim(dim) + { + m_stride = 1; + m_offset = 0; + for (unsigned int i = 0; i < dim; ++i) + { + m_offset += coord[i] * m_stride; + m_stride *= m_ts->m_sizes[i]; + } + } + + /** + * @brief Access element of a TensorStorage object through the view. + * + * This function allows modification of the underlying TensorStorage object. + * + * @param idx Index in the viewed dimension. + * @return DataType& + */ + DataType& operator[](size_t idx) + { + return m_ts->m_data[m_offset + idx * m_stride]; + } + + /** + * @brief Access element of a TensorStorage object through the view. + * + * This function allows read-only access, supporting const TensorStorage objects. + * + * @param idx + * @return DataType + */ + DataType operator[](size_t idx) const + { + return m_ts->m_data[m_offset + idx * m_stride]; + } + + /** + * @brief Creates a copy of the data in the viewed dimension. + * + * This function generates a new Array object and populates + * it with the view data. No link back to the original TensorStorage is + * maintained. + * + * @return Array + */ + Array extract() + { + Array x(m_ts->size(m_dim)); + for (size_t i = 0; i < m_ts->size(m_dim); ++i) { + x[i] = m_ts->m_data[m_offset + i*m_stride]; + } + return x; + } + + /** + * @brief Get the size of the view. + * + * @return size_t + */ + size_t size() + { + return m_ts->size(m_dim); + } + + private: + std::shared_ptr> m_ts; // Underlying TensorStorage object. + unsigned int m_dim; // The dimension being viewed. + size_t m_offset; // Offset of the first entry in the view (based on coordinate). + size_t m_stride; // Stride for accessing subsequent elements of view. + }; + + /** + * @brief Factory for creating a TensorStorage object. + * + * @param size1 Size in first dimension. + * @param size2 Size in second dimension. + * @param val Default value. + * @return std::shared_ptr> + */ + static std::shared_ptr> create(const size_t size1, const size_t size2, const DataType val = 0.0) { ASSERTL0(size1 > 0, "Size1 must be greater than zero."); ASSERTL0(size2 > 0, "Size2 must be greater than zero."); - // Nothing to do + + return std::shared_ptr>(new TensorStorage(size1, size2, val)); } /** - * @brief Construct a new TensorStorage object with two given sizes. + * @brief Factory for creating a TensorStorage object. * - * @param sizes Vector of arbitrary number of sizes > 0 - * @param val Initial vlaue for all elements. + * @param sizes Vector of sizes for all dimensions. + * @param val Default value. + * @return std::shared_ptr> */ - TensorStorage(const std::vector sizes, const DataType val = 0.0) - : m_sizes(sizes) + static std::shared_ptr> create(const std::vector sizes, const DataType val = 0.0) { - ASSERTL1(!sizes.empty(), "Must have at least one size."); - int size = 1; - for (auto &s : sizes) - { - size *= s; - } - m_data = Array(size, val); + return std::shared_ptr>(new TensorStorage(sizes, val)); } /** @@ -97,6 +182,31 @@ class TensorStorage } + /** + * @brief Returns the number of dimensions of the TensorStorage object. + * + */ + size_t dims() + { + return m_sizes.size(); + } + + + /** + * @brief Returns the total size of the TensorStorage object. + * + * @return size_t + */ + size_t size() + { + size_t s = 1; + for (auto &i : m_sizes) + { + s *= i; + } + return s; + } + /** * @brief Returns the size of the TensorStorage object. * @@ -110,6 +220,30 @@ class TensorStorage return m_sizes[idx]; } + /** + * @brief Access an element of the storage container. + * + * @param idx + * @return DataType& + */ + DataType& operator[](size_t idx) + { + return m_data[idx]; + } + + /** + * @brief Access an element of the storage container. + * + * This routine provides read-only access to the storage container, which works with const objects. + * + * @param idx Index to access. + * @return DataType + */ + DataType operator[](size_t idx) const + { + return m_data[idx]; + } + /** * @brief Returns whether the storage represents physical solution * or spectral/hp element coefficients. @@ -134,24 +268,12 @@ class TensorStorage * @brief * */ - Array slice(unsigned int dim, std::vector coord) + View slice(unsigned int dim, std::vector coord) { ASSERTL0(dim < m_sizes.size(), "dim larger than number of dimensions."); ASSERTL0(coord.size() == m_sizes.size(), "coord dimension must match dimension of Tensor."); - - Array x(m_sizes[dim]); - size_t stride = 1; - size_t offset = 0; - for (unsigned int i = 0; i < dim; ++i) - { - offset += coord[i] * stride; - stride *= m_sizes[i]; - } - - for (size_t i = 0; i < m_sizes[dim]; ++i) { - x[i] = m_data[offset + i*stride]; - } - return std::move(x); + std::shared_ptr> ts = this->shared_from_this(); + return std::move(View(ts, dim, coord)); } // slice @@ -160,6 +282,38 @@ class TensorStorage std::vector m_sizes; /// Sizes for each sub-component of the tensor. enum TensorStorage::State m_state = ePhys; /// Flag indicating if data is in coeff or physical state. + /** + * @brief Construct a new TensorStorage object with two given sizes. + * + * @param size1 First component size. + * @param size2 Second component size. + * @param val Initial value for all elements. + */ + TensorStorage(const size_t size1, const size_t size2, const DataType val = 0.0) + : m_data(Array(size1 * size2, val)), + m_sizes({size1, size2}) + { + // Nothing to do + } + + /** + * @brief Construct a new TensorStorage object with two given sizes. + * + * @param sizes Vector of arbitrary number of sizes > 0 + * @param val Initial vlaue for all elements. + */ + TensorStorage(const std::vector sizes, const DataType val = 0.0) + : m_sizes(sizes) + { + ASSERTL1(!sizes.empty(), "Must have at least one size."); + int size = 1; + for (auto &s : sizes) + { + size *= s; + } + m_data = Array(size, val); + } + }; /// A shared pointer to a TensorStorage object. -- GitLab From c386e0840cebb8eaa4f6d671fb2ae2f03be2c501 Mon Sep 17 00:00:00 2001 From: Chris Cantwell Date: Wed, 30 Mar 2022 14:55:29 +0100 Subject: [PATCH 4/6] Updates to Demo code and documentation. --- library/Demos/TensorRegions/Domain.cpp | 21 +++++++++++---- library/TensorRegions/TensorRegion.cpp | 11 +++++--- library/TensorRegions/TensorStorage.hpp | 34 ++++++++++++++++++++++++- 3 files changed, 56 insertions(+), 10 deletions(-) diff --git a/library/Demos/TensorRegions/Domain.cpp b/library/Demos/TensorRegions/Domain.cpp index 139651303..416393f3c 100644 --- a/library/Demos/TensorRegions/Domain.cpp +++ b/library/Demos/TensorRegions/Domain.cpp @@ -8,10 +8,6 @@ using namespace std; using namespace Nektar; -// compile using Builds/Demos/StdRegions -> make DEBUG=1 ProjectLoc1D - -// This routine projects a polynomial which has energy in all mdoes of -// the expansions and report an error. int main(int argc, char *argv[]) { @@ -27,6 +23,7 @@ int main(int argc, char *argv[]) SpatialDomains::BoundaryConditions Allbcs(vSession, graph1D); bool SetToOneSpaceDimension = true; + // Construct the constituent expansions std::vector exps(nDomains); for (unsigned int i = 0; i < nDomains; ++i) { @@ -34,27 +31,41 @@ int main(int argc, char *argv[]) MemoryManager::AllocateSharedPtr( vSession, graph1D, domain[i], Allbcs, vSession->GetVariable(0), SetToOneSpaceDimension); + cout << "Expansion " << i << ": " << exps[i]->GetNumElmts() + << " elements" << endl; } + // Construct a TensorRegion from the expansions TensorRegions::TensorRegion tr(exps); + // Allocate a TensorStorage object to hold a discrete solution TensorRegions::TensorRegionStorageSharedPtr st = tr.AllocateStorage(); - + cout << "Number of dimensions: " << st->dims() << endl; cout << "Total tensor storage size: " << st->size() << endl; for (size_t i = 0; i < st->dims(); ++i) { cout << "Dim " << i << " has size " << st->size(i) << " points." << endl; } + // Construct views on the TensorStorage object to operate on slices of the + // data along different directions. TensorRegions::TensorRegionStorage::View st0 = st->slice(0, {0,0}); TensorRegions::TensorRegionStorage::View st1 = st->slice(1, {0,0}); cout << "Size of slice 0 is " << st0.size() << endl; + cout << "Size of slice 1 is " << st1.size() << endl; + + // Update TensorStorage object through View. for (int i = 0; i < st1.size(); ++i) { st1[i] = i; } + + // Extract Array data from View Array data0 = st0.extract(); Array data1 = st1.extract(); cout << "Size of slice 0 data is " << data0.size() << endl; + cout << "Size of slice 1 data is " << data1.size() << endl; + + // Display all entries of TensorStorage object. for (int i = 0; i < st->size(); ++i) { cout << (*st)[i] << endl; } diff --git a/library/TensorRegions/TensorRegion.cpp b/library/TensorRegions/TensorRegion.cpp index 2d0cc6c6c..37c71dc1b 100644 --- a/library/TensorRegions/TensorRegion.cpp +++ b/library/TensorRegions/TensorRegion.cpp @@ -9,10 +9,13 @@ namespace TensorRegions /** * @class TensorRegion * - * A TensorRegion represents a high-dimensional multi-region finite element - * expansion. It is constructed through a tensor-product of two lower- - * dimensional MultiRegions expansions, defined by two independent meshes. - * Operators are constructed by applying the sub-factorisation approach. + * A TensorRegion represents a high-dimensional finite element + * expansion. It is constructed through a tensor-product of two or more lower- + * dimensional MultiRegions expansions, defined by independent meshes. The + * solution storage for a TensorRegion is provided by the TensorStorage class. + * Finite element operators are constructed by applying the sub-factorisation + * approach, deferring to the lower-dimensional consistuent MultiRegions classes + * for their implementation. */ /** diff --git a/library/TensorRegions/TensorStorage.hpp b/library/TensorRegions/TensorStorage.hpp index 7235f86ba..ab2bb1f4a 100644 --- a/library/TensorRegions/TensorStorage.hpp +++ b/library/TensorRegions/TensorStorage.hpp @@ -10,10 +10,16 @@ namespace Nektar namespace TensorRegions { +/** + * @namespace TensorRegions + * @brief Namespace encapsulating algorithms and data structures + * for constructing high-dimensional high-order finite element spaces. + */ /** * @brief Storage container for TensorRegions solution. * + * * Data is stored fastest in the lowest dimension. * * @tparam DataType The type of data stored by the container. @@ -35,6 +41,12 @@ class TensorStorage : public std::enable_shared_from_this& pData) + { + ASSERTL1(pData.size() == m_ts->size(m_dim), + "Injected array size does not match dimension size."); + for (size_t i = 0; i < m_ts->size(m_dim); ++i) { + m_ts->m_data[m_offset + i*m_stride] = pData[i]; + } + } + /** * @brief Get the size of the view. * @@ -265,8 +292,13 @@ class TensorStorage : public std::enable_shared_from_this coord) { -- GitLab From 09422bdccd09e92742d70ad5dcca7a81b2e96926 Mon Sep 17 00:00:00 2001 From: Chris Cantwell Date: Wed, 30 Mar 2022 14:56:13 +0100 Subject: [PATCH 5/6] Add demo mesh and conditions. --- .../Demos/TensorRegions/tensor-conditions.xml | 16 ++++++ library/Demos/TensorRegions/tensor.geo | 7 +++ library/Demos/TensorRegions/tensor.msh | 50 +++++++++++++++++++ library/Demos/TensorRegions/tensor.xml | 31 ++++++++++++ 4 files changed, 104 insertions(+) create mode 100644 library/Demos/TensorRegions/tensor-conditions.xml create mode 100644 library/Demos/TensorRegions/tensor.geo create mode 100644 library/Demos/TensorRegions/tensor.msh create mode 100644 library/Demos/TensorRegions/tensor.xml diff --git a/library/Demos/TensorRegions/tensor-conditions.xml b/library/Demos/TensorRegions/tensor-conditions.xml new file mode 100644 index 000000000..0f13fd5ce --- /dev/null +++ b/library/Demos/TensorRegions/tensor-conditions.xml @@ -0,0 +1,16 @@ + + + + +

Lambda = 1

+ + + + u + + + + + + + diff --git a/library/Demos/TensorRegions/tensor.geo b/library/Demos/TensorRegions/tensor.geo new file mode 100644 index 000000000..57ff79600 --- /dev/null +++ b/library/Demos/TensorRegions/tensor.geo @@ -0,0 +1,7 @@ +Point(1) = {0.0, 0.0, 0.0, 0.2}; +Point(2) = {1.0, 0.0, 0.0, 0.2}; +Point(3) = {0.0, 1.0, 0.0, 1.0}; +Line(1) = {1, 2}; +Line(2) = {1, 3}; +Physical Line(0) = {1}; +Physical Line(1) = {2}; diff --git a/library/Demos/TensorRegions/tensor.msh b/library/Demos/TensorRegions/tensor.msh new file mode 100644 index 000000000..bde714c6e --- /dev/null +++ b/library/Demos/TensorRegions/tensor.msh @@ -0,0 +1,50 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$Entities +3 2 0 0 +1 0 0 0 0 +2 1 0 0 0 +3 0 1 0 0 +1 0 0 0 1 0 0 1 0 2 1 -2 +2 0 0 0 0 1 0 1 1 2 1 -3 +$EndEntities +$Nodes +5 9 1 9 +0 1 0 1 +1 +0 0 0 +0 2 0 1 +2 +1 0 0 +0 3 0 1 +3 +0 1 0 +1 1 0 4 +4 +5 +6 +7 +0.1999999999995579 0 0 +0.3999999999989749 0 0 +0.5999999999989468 0 0 +0.7999999999994734 0 0 +1 2 0 2 +8 +9 +0 0.2305846405365542 0 +0 0.5519932425989125 0 +$EndNodes +$Elements +2 8 1 8 +1 1 1 5 +1 1 4 +2 4 5 +3 5 6 +4 6 7 +5 7 2 +1 2 1 3 +6 1 8 +7 8 9 +8 9 3 +$EndElements diff --git a/library/Demos/TensorRegions/tensor.xml b/library/Demos/TensorRegions/tensor.xml new file mode 100644 index 000000000..8f13ee231 --- /dev/null +++ b/library/Demos/TensorRegions/tensor.xml @@ -0,0 +1,31 @@ + + + + eJxjYMAPGKF0YvRMIDhpjy7PBKXLA0HyNzHkmaG0Kp8xEDzGkGeB0kLtIP0vMeRZUXgfMOTZ0Pin7Dg0zzSfhatjR5PnyZy0++3Kh3B5DnQD0ewBAN5PGUEA + + eJxjYEAFjDhoJhw0Mw6aBQfNikbDABsOmh0HzQGlARdgAFgA + + + S[0-4] + S[5-7] + + + C[0] + C[1] + + + + + refs/heads/feature/tensorregions + d22b08d80e3dddcaad8f6029194f201a9a638a49 + kingscross + 5.2.0 + 29-Mar-2022 15:27:18 + + tensor.msh tensor.xml + + + + + + -- GitLab From 18a8e982d487fa12a0359b05e945be7d1ca2d5ae Mon Sep 17 00:00:00 2001 From: Chris Cantwell Date: Wed, 30 Mar 2022 19:38:29 +0100 Subject: [PATCH 6/6] Add shift operation. --- library/Demos/TensorRegions/Domain.cpp | 3 ++ library/TensorRegions/TensorStorage.hpp | 46 +++++++++++++++++++------ 2 files changed, 38 insertions(+), 11 deletions(-) diff --git a/library/Demos/TensorRegions/Domain.cpp b/library/Demos/TensorRegions/Domain.cpp index 416393f3c..7d1c1a32a 100644 --- a/library/Demos/TensorRegions/Domain.cpp +++ b/library/Demos/TensorRegions/Domain.cpp @@ -69,5 +69,8 @@ int main(int argc, char *argv[]) for (int i = 0; i < st->size(); ++i) { cout << (*st)[i] << endl; } + + st0.shift({0,1}); + return 0; } \ No newline at end of file diff --git a/library/TensorRegions/TensorStorage.hpp b/library/TensorRegions/TensorStorage.hpp index ab2bb1f4a..5e572da0b 100644 --- a/library/TensorRegions/TensorStorage.hpp +++ b/library/TensorRegions/TensorStorage.hpp @@ -59,16 +59,10 @@ class TensorStorage : public std::enable_shared_from_this> ts, unsigned int dim, std::vector coord) - : m_ts(ts), m_dim(dim) + View(std::shared_ptr> ts, unsigned int dim, const std::vector& coord) + : m_ts(ts), m_dim(dim), m_coord(coord) { - m_stride = 1; - m_offset = 0; - for (unsigned int i = 0; i < dim; ++i) - { - m_offset += coord[i] * m_stride; - m_stride *= m_ts->m_sizes[i]; - } + computeOffsetStride(); } /** @@ -97,6 +91,21 @@ class TensorStorage : public std::enable_shared_from_thism_data[m_offset + idx * m_stride]; } + /** + * @brief Shifts the view reference coordinate by the given + * translation vector. + * + * @param s Translation vector to shift by. + */ + void shift(const std::vector& s) + { + for (unsigned int i = 0; i < m_dim; ++i) { + ASSERTL0(m_coord[i] + s[i] < m_ts->size(i), + "Requested shift takes coordinate out of range."); + } + computeOffsetStride(); + } + /** * @brief Creates a copy of the data in the viewed dimension. * @@ -143,8 +152,24 @@ class TensorStorage : public std::enable_shared_from_this> m_ts; // Underlying TensorStorage object. unsigned int m_dim; // The dimension being viewed. + std::vector m_coord; // The reference coordinate for the view. size_t m_offset; // Offset of the first entry in the view (based on coordinate). size_t m_stride; // Stride for accessing subsequent elements of view. + + /** + * @brief Calculate the member variables m_offset and m_stride + * from the given reference coordinate m_coord. + */ + void computeOffsetStride() + { + m_stride = 1; + m_offset = 0; + for (unsigned int i = 0; i < m_dim; ++i) + { + m_offset += m_coord[i] * m_stride; + m_stride *= m_ts->m_sizes[i]; + } + } }; /** @@ -300,14 +325,13 @@ class TensorStorage : public std::enable_shared_from_this coord) + View slice(unsigned int dim, const std::vector& coord) { ASSERTL0(dim < m_sizes.size(), "dim larger than number of dimensions."); ASSERTL0(coord.size() == m_sizes.size(), "coord dimension must match dimension of Tensor."); std::shared_ptr> ts = this->shared_from_this(); return std::move(View(ts, dim, coord)); } - // slice private: Array m_data; /// Storage for actual data. -- GitLab