diff --git a/AssemblerLib/CMakeLists.txt b/AssemblerLib/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..27cdff53542f05b3562d1329544306e07b2d6d09 --- /dev/null +++ b/AssemblerLib/CMakeLists.txt @@ -0,0 +1,20 @@ +#Source files grouped by a directory +GET_SOURCE_FILES(SOURCES_ASSEMBLERLIB) +SET ( SOURCES ${SOURCES_ASSEMBLERLIB}) + +# Create the library +ADD_LIBRARY( AssemblerLib STATIC ${SOURCES} ) + +INCLUDE_DIRECTORIES ( + . + ../BaseLib + ../GeoLib + ../MathLib + ../MeshLib +) + +TARGET_LINK_LIBRARIES (AssemblerLib + ${Boost_LIBRARIES} + MeshLib +) + diff --git a/AssemblerLib/ComponentGlobalIndexDict.h b/AssemblerLib/ComponentGlobalIndexDict.h new file mode 100644 index 0000000000000000000000000000000000000000..d812225a1bca04cc8ed31128593003b482750dc6 --- /dev/null +++ b/AssemblerLib/ComponentGlobalIndexDict.h @@ -0,0 +1,123 @@ +/** + * \author Norihiro Watanabe + * \date 2013-04-16 + * + * \copyright + * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef ASSEMBLERLIB_COMPONENTGLOBALINDEXDICT_H_ +#define ASSEMBLERLIB_COMPONENTGLOBALINDEXDICT_H_ + +#include <iostream> +#include <limits> + +#include <boost/multi_index/composite_key.hpp> +#include <boost/multi_index/member.hpp> +#include <boost/multi_index/ordered_index.hpp> +#include <boost/multi_index_container.hpp> + +#include "MeshLib/Location.h" + +namespace AssemblerLib +{ + +/// \internal +namespace detail +{ + +struct Line +{ + MeshLib::Location location; + + // Physical component + std::size_t comp_id; + + // Position in global matrix or vector + std::size_t global_index; + + Line(MeshLib::Location const& l, std::size_t c, std::size_t i) + : location(l), comp_id(c), global_index(i) + {} + + Line(MeshLib::Location const& l, std::size_t c) + : location(l), comp_id(c), + global_index(std::numeric_limits<std::size_t>::max()) + {} + + explicit Line(MeshLib::Location const& l) + : location(l), + comp_id(std::numeric_limits<std::size_t>::max()), + global_index(std::numeric_limits<std::size_t>::max()) + {} + + friend std::ostream& operator<<(std::ostream& os, Line const& l) + { + return os << l.location << ", " << l.comp_id << ", " << l.global_index; + } +}; + +struct LineByLocationComparator +{ + bool operator()(Line const& a, Line const& b) const + { + return a.location < b.location; + } +}; + +struct LineByLocationAndComponentComparator +{ + bool operator()(Line const& a, Line const& b) const + { + if (a.location < b.location) + return true; + if (b.location < a.location) + return false; + + // a.loc == b.loc + return a.comp_id < b.comp_id; + } +}; + +struct ByLocation {}; +struct ByLocationAndComponent {}; +struct ByComponent {}; +struct ByGlobalIndex {}; + +typedef boost::multi_index::multi_index_container< + Line, + boost::multi_index::indexed_by + < + boost::multi_index::ordered_unique + < + boost::multi_index::tag<ByLocationAndComponent>, + boost::multi_index::identity<Line>, + LineByLocationAndComponentComparator + >, + boost::multi_index::ordered_non_unique + < + boost::multi_index::tag<ByLocation>, + boost::multi_index::identity<Line>, + LineByLocationComparator + >, + boost::multi_index::ordered_non_unique + < + boost::multi_index::tag<ByComponent>, + boost::multi_index::member<Line, std::size_t, &Line::comp_id> + >, + boost::multi_index::ordered_non_unique + < + boost::multi_index::tag<ByGlobalIndex>, + boost::multi_index::member<Line, std::size_t, &Line::global_index> + > + > + > ComponentGlobalIndexDict; + +} // namespace detail +} // namespace AssemblerLib + +#endif // ASSEMBLERLIB_COMPONENTGLOBALINDEXDICT_H_ diff --git a/AssemblerLib/MeshComponentMap.cpp b/AssemblerLib/MeshComponentMap.cpp new file mode 100644 index 0000000000000000000000000000000000000000..8c550d51005c44b84aadf7e533b71d4ffcf38fdb --- /dev/null +++ b/AssemblerLib/MeshComponentMap.cpp @@ -0,0 +1,152 @@ +/** + * \author Norihiro Watanabe + * \date 2013-04-16 + * + * \copyright + * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#include <iostream> + +#include <boost/range/algorithm/for_each.hpp> + +#include "MeshLib/MeshSubsets.h" + +#include "MeshComponentMap.h" + +namespace AssemblerLib +{ + +using namespace detail; + +std::size_t const MeshComponentMap::nop = std::numeric_limits<std::size_t>::max(); + +MeshComponentMap::MeshComponentMap( + const std::vector<MeshLib::MeshSubsets*> &components, ComponentOrder order) +{ + // construct dict (and here we number global_index by component type) + std::size_t global_index = 0; + std::size_t comp_id = 0; + for (auto c = components.cbegin(); c != components.cend(); ++c) + { + for (unsigned mesh_subset_index = 0; mesh_subset_index < (*c)->size(); mesh_subset_index++) + { + MeshLib::MeshSubset const& mesh_subset = (*c)->getMeshSubset(mesh_subset_index); + std::size_t const mesh_id = mesh_subset.getMeshID(); + // mesh items are ordered first by node, cell, .... + for (std::size_t j=0; j<mesh_subset.getNNodes(); j++) + _dict.insert(Line(Location(mesh_id, MeshLib::MeshItemType::Node, j), comp_id, global_index++)); + for (std::size_t j=0; j<mesh_subset.getNElements(); j++) + _dict.insert(Line(Location(mesh_id, MeshLib::MeshItemType::Cell, j), comp_id, global_index++)); + } + comp_id++; + } + + if (order == ComponentOrder::BY_LOCATION) + renumberByLocation(); +} + +void MeshComponentMap::renumberByLocation(std::size_t offset) +{ + std::size_t global_index = offset; + + auto &m = _dict.get<ByLocation>(); // view as sorted by mesh item + for (auto itr_mesh_item=m.begin(); itr_mesh_item!=m.end(); ++itr_mesh_item) + { + Line pos = *itr_mesh_item; + pos.global_index = global_index++; + m.replace(itr_mesh_item, pos); + } +} + +std::vector<std::size_t> MeshComponentMap::getComponentIDs(const Location &l) const +{ + auto const &m = _dict.get<ByLocation>(); + auto const p = m.equal_range(Line(l)); + std::vector<std::size_t> vec_compID; + for (auto itr=p.first; itr!=p.second; ++itr) + vec_compID.push_back(itr->comp_id); + return vec_compID; +} + +std::size_t MeshComponentMap::getGlobalIndex(Location const& l, + std::size_t const c) const +{ + auto const &m = _dict.get<ByLocationAndComponent>(); + auto const itr = m.find(Line(l, c)); + return itr!=m.end() ? itr->global_index : nop; +} + +std::vector<std::size_t> MeshComponentMap::getGlobalIndices(const Location &l) const +{ + auto const &m = _dict.get<ByLocation>(); + auto const p = m.equal_range(Line(l)); + std::vector<std::size_t> global_indices; + for (auto itr=p.first; itr!=p.second; ++itr) + global_indices.push_back(itr->global_index); + return global_indices; +} + +template <> +std::vector<std::size_t> +MeshComponentMap::getGlobalIndices<ComponentOrder::BY_LOCATION>( + std::vector<Location> const &ls) const +{ + // Create vector of global indices sorted by location containing all + // locations given in ls parameter. + + std::vector<std::size_t> global_indices; + global_indices.reserve(ls.size()); + + auto const &m = _dict.get<ByLocation>(); + for (auto l = ls.cbegin(); l != ls.cend(); ++l) + { + auto const p = m.equal_range(Line(*l)); + for (auto itr = p.first; itr != p.second; ++itr) + global_indices.push_back(itr->global_index); + } + + return global_indices; +} + +template <> +std::vector<std::size_t> +MeshComponentMap::getGlobalIndices<ComponentOrder::BY_COMPONENT>( + std::vector<Location> const &ls) const +{ + // vector of (Component, global Index) pairs. + typedef std::pair<std::size_t, std::size_t> CIPair; + std::vector<CIPair> pairs; + pairs.reserve(ls.size()); + + // Create a sub dictionary containing all lines with location from ls. + auto const &m = _dict.get<ByLocation>(); + for (auto l = ls.cbegin(); l != ls.cend(); ++l) + { + auto const p = m.equal_range(Line(*l)); + for (auto itr = p.first; itr != p.second; ++itr) + pairs.emplace_back(itr->comp_id, itr->global_index); + } + + auto CIPairLess = [](CIPair const& a, CIPair const& b) + { + return a.first < b.first; + }; + + // Create vector of global indices from sub dictionary sorting by component. + if (!std::is_sorted(pairs.begin(), pairs.end(), CIPairLess)) + std::stable_sort(pairs.begin(), pairs.end(), CIPairLess); + + std::vector<std::size_t> global_indices; + global_indices.reserve(pairs.size()); + for (auto p = pairs.cbegin(); p != pairs.cend(); ++p) + global_indices.push_back(p->second); + + return global_indices; +} + +} // namespace AssemblerLib diff --git a/AssemblerLib/MeshComponentMap.h b/AssemblerLib/MeshComponentMap.h new file mode 100644 index 0000000000000000000000000000000000000000..b4543e1c7432936d868acc633c323fb6e00465c9 --- /dev/null +++ b/AssemblerLib/MeshComponentMap.h @@ -0,0 +1,140 @@ +/** + * \author Norihiro Watanabe + * \date 2013-04-16 + * + * \copyright + * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef ASSEMBLERLIB_MESHCOMPONENTMAP_H_ +#define ASSEMBLERLIB_MESHCOMPONENTMAP_H_ + +#include <vector> + +#include "MeshLib/Location.h" + +#include "ComponentGlobalIndexDict.h" + +namespace MeshLib +{ + class MeshSubsets; +} + +namespace AssemblerLib +{ + +/// Ordering of components in global matrix/vector. +enum class ComponentOrder +{ + BY_COMPONENT, ///< Ordering data by component type + BY_LOCATION ///< Ordering data by spatial location +}; + +/// Multidirectional mapping between mesh entities and degrees of freedom. +class MeshComponentMap +{ +public: + typedef MeshLib::Location Location; +public: + /// \param components a vector of components + /// \param order type of ordering values in a vector + MeshComponentMap(std::vector<MeshLib::MeshSubsets*> const& components, + ComponentOrder order); + + /// The number of components in the map. + std::size_t size() const + { + return _dict.size(); + } + + /// Component ids at given location \c l. + /// + /// | Location | ComponentID | + /// | -------- | ----------- | + /// | l | comp_id_1 | + /// | l | ... | + /// | l | comp_id_n | + std::vector<std::size_t> getComponentIDs(const Location &l) const; + + /// Global index of the given component \c c at given location \c l. + /// + /// | Location | ComponentID | GlobalIndex | + /// | -------- | ----------- | ----------- | + /// | l | c | gi | + std::size_t getGlobalIndex(Location const &l, std::size_t const c) const; + + /// Global indices for all components at the given location \c l. + /// + /// If there is more than one component at the given location, the + /// function returns a vector containing global indices for each component. + /// + /// | Location | ComponentID | GlobalIndex | + /// | -------- | ----------- | ----------- | + /// | l | comp_id_1 | gi23 | + /// | ... | ... | ... | + /// | l | comp_id_k | gi45 | + std::vector<std::size_t> getGlobalIndices(const Location &l) const; + + /// Global indices for all components at all given locations \c ls ordered + /// as required by the template parameter ORDER. + /// + /// In case ORDER is ComponentOrder::BY_LOCATION the return list is sorted + /// first by location. + /// | Location | ComponentID | GlobalIndex | + /// | -------- | ----------- | ----------- | + /// | l_1 | comp_id_1 | gi23 | + /// | ... | ... | ... | + /// | l_1 | comp_id_k | gi45 | + /// | l_2 | comp_id_1 | gi46 | + /// | ... | ... | ... | + /// | l_2 | comp_id_m | gi67 | + /// | ... | ... | ... | + /// | l_n | comp_id_n | gi78 | + /// + /// In case ORDER is ComponentOrder::BY_COMPONENT the return list is sorted + /// first by component ids. + /// | Location | ComponentID | GlobalIndex | + /// | -------- | ----------- | ----------- | + /// | l_1 | comp_id_1 | gi23 | + /// | ... | ... | ... | + /// | l_k | comp_id_1 | gi45 | + /// | l_1 | comp_id_2 | gi46 | + /// | ... | ... | ... | + /// | l_m | comp_id_2 | gi78 | + /// | ... | ... | ... | + /// | l_n | comp_id_n | gi89 | + template <ComponentOrder ORDER> + std::vector<std::size_t> getGlobalIndices(const std::vector<Location> &ls) const; + + /// A value returned if no global index was found for the requested + /// location/component. The value is implementation dependent. + static std::size_t const nop; + +#ifndef NDEBUG + const detail::ComponentGlobalIndexDict& getDictionary() const + { + return _dict; + } + + friend std::ostream& operator<<(std::ostream& os, MeshComponentMap const& m) + { + for (auto l = m._dict.begin(); l != m._dict.end(); ++l) + os << *l << "\n"; + return os; + } +#endif // NDEBUG + +private: + void renumberByLocation(std::size_t offset=0); + +private: + detail::ComponentGlobalIndexDict _dict; +}; + +} // namespace AssemblerLib + +#endif // ASSEMBLERLIB_MESHCOMPONENTMAP_H_ diff --git a/CMakeLists.txt b/CMakeLists.txt index 7b56c3a8696ac592782112382b6794e2bfff4951..5c4666fdb7e2f20fafeb1d136ac128af0d09e6a1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -105,6 +105,7 @@ IF(OGS_USE_LIS) ADD_DEFINITIONS(-DUSE_LIS) ENDIF() +ADD_SUBDIRECTORY( AssemblerLib ) ADD_SUBDIRECTORY( BaseLib ) # TODO This is a hack but we have to make sure that Boost is built first ADD_DEPENDENCIES(BaseLib Boost) diff --git a/MeshLib/MeshSubsets.h b/MeshLib/MeshSubsets.h index 94dad6efbbae05468dcf368410de853742683569..5842ae846de80b729f59aea91aa9077b523a7de8 100644 --- a/MeshLib/MeshSubsets.h +++ b/MeshLib/MeshSubsets.h @@ -58,8 +58,8 @@ public: return _n_total_items; } - /// return the number of related meshes - std::size_t getNMeshes() const + /// Number of saved mesh subsets. + std::size_t size() const { return _mesh_subsets.size(); } diff --git a/Tests/AssemblerLib/TestMeshComponentMap.cpp b/Tests/AssemblerLib/TestMeshComponentMap.cpp new file mode 100644 index 0000000000000000000000000000000000000000..6b42d5f3da8149a763d4fcde51f5ee90a05e60c4 --- /dev/null +++ b/Tests/AssemblerLib/TestMeshComponentMap.cpp @@ -0,0 +1,138 @@ +/** + * \author Norihiro Watanabe + * \date 2013-04-16 + * + * \copyright + * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#include <gtest/gtest.h> +#include <vector> + +#include "AssemblerLib/MeshComponentMap.h" + +#include "MeshLib/MeshGenerators/MeshGenerator.h" +#include "MeshLib/MeshSubsets.h" + +class AssemblerLibMeshComponentMapTest : public ::testing::Test +{ + public: + typedef MeshLib::MeshItemType MeshItemType; + typedef MeshLib::Location Location; + typedef AssemblerLib::MeshComponentMap MeshComponentMap; + + public: + AssemblerLibMeshComponentMapTest() + : mesh(nullptr), nodesSubset(nullptr), cmap(nullptr) + { + mesh = MeshLib::MeshGenerator::generateLineMesh(1.0, mesh_size); + nodesSubset = new MeshLib::MeshSubset(*mesh, mesh->getNodes()); + + // Add two components both based on the same nodesSubset. + components.emplace_back(new MeshLib::MeshSubsets(nodesSubset)); + components.emplace_back(new MeshLib::MeshSubsets(nodesSubset)); + } + + ~AssemblerLibMeshComponentMapTest() + { + delete cmap; + std::remove_if(components.begin(), components.end(), + [](MeshLib::MeshSubsets* p) { delete p; return true; }); + delete nodesSubset; + delete mesh; + } + + static std::size_t const mesh_size = 9; + MeshLib::Mesh const* mesh; + MeshLib::MeshSubset const* nodesSubset; + + //data component 0 and 1 are assigned to all nodes in the mesh + static std::size_t const comp0_id = 0; + static std::size_t const comp1_id = 1; + std::vector<MeshLib::MeshSubsets*> components; + MeshComponentMap const* cmap; + + // + // Functions used for checking. + // + + // Returns global index of a node location and a component. + std::size_t giAtNodeForComponent(std::size_t const n, std::size_t const c) const + { + return cmap->getGlobalIndex(Location(mesh->getID(), MeshItemType::Node, n), c); + } + +}; + +TEST_F(AssemblerLibMeshComponentMapTest, CheckOrderByComponent) +{ + // - Entries in the vector are arranged in the order of a component type and then node ID + // - For example, x=[(node 0, comp 0) (node 1, comp 0) ... (node n, comp0), (node 0, comp1) ... ] + + cmap = new MeshComponentMap(components, + AssemblerLib::ComponentOrder::BY_COMPONENT); + + //std::cout << "# database \n" << *cmap << std::endl; + + ASSERT_EQ(2 * mesh->getNNodes(), cmap->size()); + for (std::size_t i = 0; i < mesh_size; i++) + { + // Test global indices for the different components of the node. + ASSERT_EQ(i , giAtNodeForComponent(i, comp0_id)); + ASSERT_EQ(mesh_size + 1 + i, giAtNodeForComponent(i, comp1_id)); + + // Test component ids of the node. + std::vector<std::size_t> const vecCompIDs = cmap->getComponentIDs( + Location(mesh->getID(), MeshItemType::Node, i)); + ASSERT_EQ(2u, vecCompIDs.size()); + ASSERT_EQ(0u, vecCompIDs[0]); + ASSERT_EQ(1u, vecCompIDs[1]); + } +} + +TEST_F(AssemblerLibMeshComponentMapTest, CheckOrderByLocation) +{ + // - Entries in the vector are arranged in the order of node ID and then a component type + // - For example, x=[(node 0, comp 0) (node 0, comp 1) ... (node n, comp0), (node n, comp1) ] + + cmap = new MeshComponentMap(components, + AssemblerLib::ComponentOrder::BY_LOCATION); + + //std::cout << "# database \n" << *cmap << std::endl; + + ASSERT_EQ(2 * mesh->getNNodes(), cmap->size()); + for (std::size_t i = 0; i < mesh_size; i++) + { + // Test global indices for the different components of the node. + ASSERT_EQ(2 * i , giAtNodeForComponent(i, comp0_id)); + ASSERT_EQ(2 * i + 1, giAtNodeForComponent(i, comp1_id)); + + // Test component ids of the node. + std::vector<std::size_t> const vecCompIDs = cmap->getComponentIDs( + Location(mesh->getID(), MeshItemType::Node, i)); + ASSERT_EQ(2u, vecCompIDs.size()); + ASSERT_EQ(0u, vecCompIDs[0]); + ASSERT_EQ(1u, vecCompIDs[1]); + } +} + +TEST_F(AssemblerLibMeshComponentMapTest, OutOfRangeAccess) +{ + cmap = new MeshComponentMap(components, + AssemblerLib::ComponentOrder::BY_COMPONENT); + + //std::cout << "# database \n" << *cmap << std::endl; + + ASSERT_EQ(MeshComponentMap::nop, cmap->getGlobalIndex( + Location(mesh->getID(), MeshItemType::Node, mesh_size + 1), comp0_id)); + ASSERT_EQ(MeshComponentMap::nop, cmap->getGlobalIndex( + Location(mesh->getID() + 1, MeshItemType::Node, 0), comp0_id)); + ASSERT_EQ(MeshComponentMap::nop, cmap->getGlobalIndex( + Location(mesh->getID(), MeshItemType::Cell, 0), comp0_id)); + ASSERT_EQ(MeshComponentMap::nop, cmap->getGlobalIndex( + Location(mesh->getID(), MeshItemType::Node, 0), 10)); +} diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt index f4a48fb52ec898af4cf323f45bf0b1bd8929bc0d..14335f8df9e4c58f0bba31f7e90cd1e0cef0d46d 100644 --- a/Tests/CMakeLists.txt +++ b/Tests/CMakeLists.txt @@ -5,6 +5,7 @@ IF( MSVC ) ENDIF() APPEND_SOURCE_FILES(TEST_SOURCES) +APPEND_SOURCE_FILES(TEST_SOURCES AssemblerLib) APPEND_SOURCE_FILES(TEST_SOURCES BaseLib) APPEND_SOURCE_FILES(TEST_SOURCES GeoLib) APPEND_SOURCE_FILES(TEST_SOURCES MathLib) @@ -15,6 +16,7 @@ IF (QT4_FOUND) ENDIF() INCLUDE_DIRECTORIES( + ${CMAKE_SOURCE_DIR}/AssemblerLib ${CMAKE_SOURCE_DIR}/BaseLib ${CMAKE_SOURCE_DIR}/FemLib ${CMAKE_SOURCE_DIR}/FileIO @@ -33,6 +35,7 @@ SET_TARGET_PROPERTIES(testrunner PROPERTIES FOLDER Testing) TARGET_LINK_LIBRARIES(testrunner GTest + AssemblerLib BaseLib FemLib FileIO diff --git a/scripts/docs/Doxyfile.in b/scripts/docs/Doxyfile.in index 47a6b57c92f65a2aca9df8d970853cecafafee54..02f34fa368a0a4271303fc36ae81798e833b4c9d 100644 --- a/scripts/docs/Doxyfile.in +++ b/scripts/docs/Doxyfile.in @@ -664,7 +664,8 @@ WARN_LOGFILE = # directories like "/usr/src/myproject". Separate the files or directories # with spaces. -INPUT = ${CMAKE_SOURCE_DIR}/BaseLib \ +INPUT = ${CMAKE_SOURCE_DIR}/AssemblerLib \ + ${CMAKE_SOURCE_DIR}/BaseLib \ ${CMAKE_SOURCE_DIR}/FemLib \ ${CMAKE_SOURCE_DIR}/FileIO \ ${CMAKE_SOURCE_DIR}/GeoLib \