diff --git a/external/flat/LICENSE_1_0.txt b/external/flat/LICENSE_1_0.txt new file mode 100644 index 000000000..36b7cd93c --- /dev/null +++ b/external/flat/LICENSE_1_0.txt @@ -0,0 +1,23 @@ +Boost Software License - Version 1.0 - August 17th, 2003 + +Permission is hereby granted, free of charge, to any person or organization +obtaining a copy of the software and accompanying documentation covered by +this license (the "Software") to use, reproduce, display, distribute, +execute, and transmit the Software, and to prepare derivative works of the +Software, and to permit third-parties to whom the Software is furnished to +do so, all subject to the following: + +The copyright notices in the Software and this entire statement, including +the above license grant, this restriction and the following disclaimer, +must be included in all copies of the Software, in whole or in part, and +all derivative works of the Software, unless such copies or derivative +works are solely in the form of machine-executable object code generated by +a source language processor. + +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, TITLE AND NON-INFRINGEMENT. IN NO EVENT +SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE +FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, +ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +DEALINGS IN THE SOFTWARE. diff --git a/external/flat/README.md b/external/flat/README.md new file mode 100644 index 000000000..e3b89eab0 --- /dev/null +++ b/external/flat/README.md @@ -0,0 +1,65 @@ +# C++ Flat Containers Library + +Fast+efficient associative containers using sorted arrays, with an interface based on standard containers. + +This was part of a C++ standard proposal and I recommend that you read it for more details: http://pubby.github.io/proposal.html + +#### Container Adaptors +- `fc::flat_map` +- `fc::flat_multimap` +- `fc::flat_set` +- `fc::flat_multiset` + +#### Class Aliases +- `fc::vector_map` +- `fc::vector_multimap` +- `fc::vector_set` +- `fc::vector_multiset` + +#### New Member Functions +- `has` + - `map.has(key)` returns a pointer to key's mapped value if it exists, otherwise returns null. +- `insert` (delayed sort) + - `map.insert(first, last, fc::delay_sort)` performs insertion with a delayed sort optimization. +- Constructors (delayed sort overload) + - Constructs flat container using a delayed sort optimization. +- Constructors (container construct overload) + - Constructs flat container by forwarding arguments to the underlying container member. + +#### What's an adaptor? + +Container adaptors allow you to use any vector-like sequence container for the underlying storage. +`std::vector` is the natural choice, but classes like `boost::small_vector` and `boost::static_vector` are useful too. Because `std::vector` is the most common in usage, aliases are provided for convenience. + +For basic use, just use the `std::vector` aliases. + +#### The public Members: `container` and `underlying` + +The public member `container` allows access to the underlying storage container. Note that it is the user's responsibility to keep the container in a sorted, unique state. + +Likewise, the public member of iterators: `underlying`, allows access to the underlying storage container's iterators. + +*Example: Delayed sort optimization using `.container`* + + std::flat_multiset set; + while(cpu_temperature() < 100) + set.container.push_back(cpu_temperature()); + std::sort(set.begin(), set.end()); + +#### Const Iteration by Default + +For safety reasons, flat container iterators are const by default. To bypass this safety and get non-const iterators, one can either iterate `.container` or take the `.underlying` of the iterator. + +*Example: Modify values in a way that preserves sortedness* + + for(auto& v : set.container) + v *= 2; + +*Example: Same thing but with `.underlying`* + + for(auto it = v.begin(); it != v.end(); ++it) + (*v.underlying) *= 2; + +#### Helper Types + +The directory `include_extra` contains convenience typedefs for use with Boost.Container. diff --git a/external/flat/flat_map.hpp b/external/flat/flat_map.hpp new file mode 100644 index 000000000..7e6ce507c --- /dev/null +++ b/external/flat/flat_map.hpp @@ -0,0 +1,251 @@ +// Copyright Pubby 2016 +// Distributed under the Boost Software License, Version 1.0. +// http://www.boost.org/LICENSE_1_0.txt + +#ifndef LIB_FLAT_FLAT_MAP_HPP +#define LIB_FLAT_FLAT_MAP_HPP + +#include "impl/flat_impl.hpp" + +namespace fc { +namespace impl { + +template +class flat_map_base +: public flat_container_base +{ +#include "impl/container_traits.hpp" + using mapped_type = typename value_type::second_type; + using B = flat_container_base; + D const* self() const { return static_cast(this); } + D* self() { return static_cast(this); } +public: + using value_compare = first_compare; + value_compare value_comp() const { return value_compare(B::key_comp()); } + + using B::insert; + using B::erase; + + // Element access + + mapped_type const* has(key_type const& key) const + { + const_iterator it = self()->find(key); + return it == self()->end() ? nullptr : &it.underlying->second; + } + + mapped_type* has(key_type const& key) + { + iterator it = self()->find(key); + return it == self()->end() ? nullptr : &it.underlying->second; + } + + mapped_type const& at(key_type const& key) const + { + if(mapped_type const* ptr = has(key)) + return *ptr; + throw std::out_of_range("flat_map::at"); + } + + mapped_type& at(key_type const& key) + { + if(mapped_type* ptr = has(key)) + return *ptr; + throw std::out_of_range("flat_map::at"); + } + + mapped_type& operator[](key_type const& key) + { return self()->try_emplace(key).first.underlying->second; } + + mapped_type& operator[](key_type&& key) + { + return self()->try_emplace(std::move(key)).first.underlying->second; + } + + // Modifiers + + std::pair insert(value_type const& value) + { return insert_(value); } + + std::pair insert(value_type&& value) + { return insert_(std::move(value)); } + + template + void insert(InputIt first, InputIt last, delay_sort_t) + { + this->ds_insert_(first, last); + auto it = std::unique( + self()->container.begin(), self()->container.end(), + impl::eq_comp{value_comp()}); + self()->container.erase(it, self()->container.end()); + } + + template + std::pair insert_or_assign(key_type const& key, M&& obj) + { return insert_or_assign_(key, std::forward(obj)); } + + template + std::pair insert_or_assign(key_type&& key, M&& obj) + { return insert_or_assign_(std::move(key), std::forward(obj)); } + + template + std::pair insert_or_assign(const_iterator hint, + key_type const& key, M&& obj) + { return insert_or_assign(key, std::forward(obj)); } + + template + std::pair insert_or_assign(const_iterator hint, + key_type&& key, M&& obj) + { return insert_or_assign(std::move(key), std::forward(obj)); } + + template + std::pair try_emplace(key_type const& key, Args&&... args) + { return try_emplace_(key, std::forward(args)...); } + + template + std::pair try_emplace(key_type&& key, Args&&... args) + { return try_emplace_(std::move(key), std::forward(args)...); } + + template + iterator try_emplace(const_iterator hint, + key_type const& key, Args&&... args) + { return try_emplace_(key, std::forward(args)...).first; } + + template + iterator try_emplace(const_iterator hint, key_type&& key, Args&&... args) + { + return try_emplace_(std::move(key), + std::forward(args)...).first; + } + + size_type erase(key_type const& key) + { + const_iterator it = self()->find(key); + if(it == self()->end()) + return 0; + self()->container.erase(it.underlying); + return 1; + } + + // Lookup + + size_type count(key_type const& key) const + { + return self()->find(key) != self()->end(); + } + +private: + template + std::pair insert_(V&& value) + { + iterator it = self()->lower_bound(value.first); + if(it == self()->end() || self()->value_comp()(value, *it)) + { + it = self()->container.insert(it.underlying, + std::forward(value)); + return std::make_pair(it, true); + } + return std::make_pair(it, false); + } + + template + std::pair insert_or_assign_(K&& key, M&& obj) + { + iterator it = self()->lower_bound(key); + if(it == self()->end() || self()->key_comp()(key, it->first)) + { + it = self()->container.insert(it.underlying, + value_type(std::forward(key), std::forward(obj))); + return std::make_pair(it, true); + } + it.underlying->second = std::forward(obj); + return std::make_pair(it, false); + } + + template + std::pair try_emplace_(K&& key, Args&&... args) + { + iterator it = self()->lower_bound(key); + if(it == self()->end() || self()->key_comp()(key, it->first)) + { + it = self()->container.emplace(it.underlying, + value_type(std::piecewise_construct, + std::forward_as_tuple(std::forward(key)), + std::forward_as_tuple(std::forward(args)...))); + return std::make_pair(it, true); + } + return std::make_pair(it, false); + } +}; + +template +class flat_map_base> +: public flat_map_base +{ +#include "impl/container_traits.hpp" + using mapped_type = typename value_type::second_type; + using B = flat_map_base; + D const* self() const { return static_cast(this); } + D* self() { return static_cast(this); } +public: + + using B::insert; + using B::count; + + // Modifiers + + template + std::pair insert(P&& value) + { + iterator it = self()->lower_bound(value.first); + if(it == self()->end() || self()->value_comp()(value, *it)) + { + it = self()->container.insert( + it.underlying, std::forward

(value)); + return std::make_pair(it, true); + } + return std::make_pair(it, false); + } + + // Lookup + + template + size_type count(K const& key) const + { + return self()->find(key) != self()->end(); + } +}; + +} // namespace impl + +template> +class flat_map +: public impl::flat_map_base, + typename Container::value_type::first_type, Container, Compare> +{ +#define FLATNAME flat_map +#define FLATKEY typename Container::value_type::first_type +#include "impl/class_def.hpp" +#undef FLATNAME +#undef FLATKEY +}; + +template> +using vector_map = flat_map>, Compare>; + +template +inline bool operator==(const flat_map& lhs, const flat_map& rhs) +{ + return lhs.container == rhs.container; +} +template +inline bool operator!=(const flat_map& lhs, const flat_map& rhs) +{ + return lhs.container != rhs.container; +} + +} // namespace fc + +#endif diff --git a/external/flat/flat_multimap.hpp b/external/flat/flat_multimap.hpp new file mode 100644 index 000000000..169be155e --- /dev/null +++ b/external/flat/flat_multimap.hpp @@ -0,0 +1,133 @@ +// Copyright Pubby 2016 +// Distributed under the Boost Software License, Version 1.0. +// http://www.boost.org/LICENSE_1_0.txt + +#ifndef LIB_FLAT_FLAT_MULTIMAP_HPP +#define LIB_FLAT_FLAT_MULTIMAP_HPP + +#include "impl/flat_impl.hpp" + +namespace fc { +namespace impl { + +template +class flat_multimap_base +: public flat_container_base +{ +#include "impl/container_traits.hpp" + using mapped_type = typename value_type::second_type; + using B = flat_container_base; + D const* self() const { return static_cast(this); } + D* self() { return static_cast(this); } +public: + using value_compare = first_compare; + value_compare value_comp() const { return value_compare(B::key_comp()); } + + using B::insert; + using B::erase; + + // Modifiers + + iterator insert(value_type const& value) + { + iterator it = self()->upper_bound(value.first); + return self()->container.insert(it.underlying, value); + } + + iterator insert(value_type&& value) + { + iterator it = self()->upper_bound(value.first); + return self()->container.insert(it.underlying, std::move(value)); + } + + template + void insert(InputIt first, InputIt last, delay_sort_t) + { this->ds_insert_(first, last); } + + size_type erase(key_type const& key) + { + auto it_pair = self()->equal_range(key); + std::size_t ret = std::distance(it_pair.first, it_pair.second); + self()->container.erase(it_pair.first.underlying, + it_pair.second.underlying); + return ret; + } + + // Lookup + + size_type count(key_type const& key) const + { + auto it_pair = self()->equal_range(key); + return std::distance(it_pair.first, it_pair.second); + } + +}; + +template +class flat_multimap_base> +: public flat_multimap_base +{ +#include "impl/container_traits.hpp" + using mapped_type = typename value_type::second_type; + using B = flat_multimap_base; + D const* self() const { return static_cast(this); } + D* self() { return static_cast(this); } +public: + + using B::insert; + using B::count; + + // Modifiers + + template + iterator insert(P&& value) + { + iterator it = self()->upper_bound(value.first); + return self()->container.insert(it.underlying, + std::forward

(value)); + } + + // Lookup + + template + size_type count(K const& key) const + { + auto it_pair = self()->equal_range(key); + return std::distance(it_pair.first, it_pair.second); + } +}; + +} // namespace impl + +template> +class flat_multimap +: public impl::flat_multimap_base, + typename Container::value_type::first_type, Container, Compare> +{ +#define FLATNAME flat_multimap +#define FLATKEY typename Container::value_type::first_type +#include "impl/class_def.hpp" +#undef FLATNAME +#undef FLATKEY +}; + +template> +using vector_multimap + = flat_multimap>, Compare>; + +template +inline bool operator==(const flat_multimap& lhs, const flat_multimap& rhs) +{ + return lhs.container == rhs.container; +} +template +inline bool operator!=(const flat_multimap& lhs, const flat_multimap& rhs) +{ + return lhs.container != rhs.container; +} + +} // namespace fc + +#endif diff --git a/external/flat/flat_multiset.hpp b/external/flat/flat_multiset.hpp new file mode 100644 index 000000000..beca937ff --- /dev/null +++ b/external/flat/flat_multiset.hpp @@ -0,0 +1,117 @@ +// Copyright Pubby 2016 +// Distributed under the Boost Software License, Version 1.0. +// http://www.boost.org/LICENSE_1_0.txt + +#ifndef LIB_FLAT_FLAT_MULTISET_HPP +#define LIB_FLAT_FLAT_MULTISET_HPP + +#include "impl/flat_impl.hpp" + +namespace fc { +namespace impl { + +template +class flat_multiset_base +: public flat_container_base +{ +#include "impl/container_traits.hpp" + using B = flat_container_base; + D const* self() const { return static_cast(this); } + D* self() { return static_cast(this); } +public: + using value_compare = Compare; + value_compare value_comp() const { return value_compare(B::key_comp()); } + + using B::insert; + using B::erase; + + // Modifiers + + iterator insert(value_type const& value) + { + iterator it = self()->upper_bound(value); + return self()->container.insert(it.underlying, value); + } + + iterator insert(value_type&& value) + { + iterator it = self()->upper_bound(value); + return self()->container.insert(it.underlying, std::move(value)); + } + + template + void insert(InputIt first, InputIt last, delay_sort_t) + { this->ds_insert_(first, last); } + + size_type erase(key_type const& key) + { + auto it_pair = self()->equal_range(key); + std::size_t ret = std::distance(it_pair.first, it_pair.second); + self()->container.erase(it_pair.first.underlying, + it_pair.second.underlying); + return ret; + } + + // Lookup + + size_type count(key_type const& key) const + { + auto it_pair = self()->equal_range(key); + return std::distance(it_pair.first, it_pair.second); + } +}; + +template +class flat_multiset_base> +: public flat_multiset_base +{ +#include "impl/container_traits.hpp" + using B = flat_multiset_base; + D const* self() const { return static_cast(this); } + D* self() { return static_cast(this); } +public: + using B::count; + + // Lookup + + template + size_type count(K const& key) const + { + auto it_pair = self()->equal_range(key); + return std::distance(it_pair.first, it_pair.second); + } +}; + +} // namespace impl + +template> +class flat_multiset +: public impl::flat_multiset_base, + typename Container::value_type, Container, Compare> +{ +#define FLATNAME flat_multiset +#define FLATKEY typename Container::value_type +#include "impl/class_def.hpp" +#undef FLATNAME +#undef FLATKEY +}; + +template> +using vector_multiset = flat_multiset, Compare>; + +template +inline bool operator==(const flat_multiset& lhs, const flat_multiset& rhs) +{ + return lhs.container == rhs.container; +} +template +inline bool operator!=(const flat_multiset& lhs, const flat_multiset& rhs) +{ + return lhs.container != rhs.container; +} + +} // namespace fc + +#endif diff --git a/external/flat/flat_set.hpp b/external/flat/flat_set.hpp new file mode 100644 index 000000000..7f95049cb --- /dev/null +++ b/external/flat/flat_set.hpp @@ -0,0 +1,129 @@ +// Copyright Pubby 2016 +// Distributed under the Boost Software License, Version 1.0. +// http://www.boost.org/LICENSE_1_0.txt + +#ifndef LIB_FLAT_FLAT_SET_HPP +#define LIB_FLAT_FLAT_SET_HPP + +#include "impl/flat_impl.hpp" + +namespace fc { +namespace impl { + +template +class flat_set_base +: public flat_container_base +{ +#include "impl/container_traits.hpp" + using B = flat_container_base; + D const* self() const { return static_cast(this); } + D* self() { return static_cast(this); } +public: + using value_compare = Compare; + value_compare value_comp() const { return value_compare(B::key_comp()); } + + using B::insert; + using B::erase; + + // Modifiers + + std::pair insert(value_type const& value) + { return insert_(value); } + + std::pair insert(value_type&& value) + { return insert_(std::move(value)); } + + template + void insert(InputIt first, InputIt last, delay_sort_t) + { + this->ds_insert_(first, last); + auto it = std::unique( + self()->container.begin(), self()->container.end(), + impl::eq_comp{value_comp()}); + self()->container.erase(it, self()->container.end()); + } + + size_type erase(key_type const& key) + { + const_iterator it = self()->find(key); + if(it == self()->end()) + return 0; + self()->container.erase(it.underlying); + return 1; + } + + // Lookup + + size_type count(key_type const& key) const + { + return self()->find(key) != self()->end(); + } + +private: + template + std::pair insert_(V&& value) + { + iterator it = self()->lower_bound(value); + if(it == self()->end() || self()->value_comp()(value, *it)) + { + it = self()->container.insert(it.underlying, + std::forward(value)); + return std::make_pair(it, true); + } + return std::make_pair(it, false); + } +}; + +template +class flat_set_base> +: public flat_set_base +{ +#include "impl/container_traits.hpp" + using B = flat_set_base; + D const* self() const { return static_cast(this); } + D* self() { return static_cast(this); } +public: + using B::count; + + // Lookup + + template + size_type count(K const& key) const + { + return self()->find(key) != self()->end(); + } +}; + +} // namespace impl + +template> +class flat_set +: public impl::flat_set_base, + typename Container::value_type, Container, Compare> +{ +#define FLATNAME flat_set +#define FLATKEY typename Container::value_type +#include "impl/class_def.hpp" +#undef FLATNAME +#undef FLATKEY +}; + +template> +using vector_set = flat_set, Compare>; + +template +inline bool operator==(const flat_set& lhs, const flat_set& rhs) +{ + return lhs.container == rhs.container; +} +template +inline bool operator!=(const flat_set& lhs, const flat_set& rhs) +{ + return lhs.container != rhs.container; +} + +} // namespace fc + +#endif diff --git a/external/flat/impl/class_def.hpp b/external/flat/impl/class_def.hpp new file mode 100644 index 000000000..d29e24b2f --- /dev/null +++ b/external/flat/impl/class_def.hpp @@ -0,0 +1,60 @@ +// Copyright Pubby 2016 +// Distributed under the Boost Software License, Version 1.0. +// http://www.boost.org/LICENSE_1_0.txt + +private: + using D = FLATNAME; + using Key = FLATKEY; +public: +#include "container_traits.hpp" + + FLATNAME() = default; + explicit FLATNAME(Compare const& comp) : comp(comp), container() {} + + template + FLATNAME(InputIt first, InputIt last) + : FLATNAME() { this->insert(first, last); } + + template + FLATNAME(InputIt first, InputIt last, Compare const& comp) + : FLATNAME(comp) { this->insert(first, last); } + + template + FLATNAME(InputIt first, InputIt last, delay_sort_t d) + : FLATNAME() { this->insert(first, last, d); } + + template + FLATNAME(InputIt first, InputIt last, Compare const& comp, delay_sort_t d) + : FLATNAME(comp) { this->insert(first, last, d); } + + FLATNAME(FLATNAME const&) = default; + FLATNAME(FLATNAME&&) = default; + + FLATNAME(std::initializer_list ilist) + : FLATNAME() { this->insert(ilist); } + + FLATNAME(std::initializer_list ilist, delay_sort_t d) + : FLATNAME() { this->insert(ilist, d); } + + FLATNAME(std::initializer_list ilist, + Compare const& comp, delay_sort_t d) + : FLATNAME(comp) { this->insert(ilist, d); } + + template + explicit FLATNAME(container_construct_t, Args&&... args) + : container(std::forward(args)...), comp() {} + + template + FLATNAME(Compare const& comp, container_construct_t, Args&&... args) + : container(std::forward(args)...), comp(comp) {} + + FLATNAME& operator=(FLATNAME const&) = default; + FLATNAME& operator=(FLATNAME&&) = default; + FLATNAME& operator=(std::initializer_list ilist) + { this->clear(); this->insert(ilist); return *this; } + + Container container; + Compare comp; + +#undef FLATNAME +#undef FLATKEY diff --git a/external/flat/impl/container_traits.hpp b/external/flat/impl/container_traits.hpp new file mode 100644 index 000000000..31e16423d --- /dev/null +++ b/external/flat/impl/container_traits.hpp @@ -0,0 +1,25 @@ +// Copyright Pubby 2016 +// Distributed under the Boost Software License, Version 1.0. +// http://www.boost.org/LICENSE_1_0.txt + +using container_type = Container; +using key_type = Key; +using size_type = typename container_type::size_type; +using difference_type = typename container_type::difference_type; +using value_type = typename container_type::value_type; +using iterator + = impl::flat_iterator< + typename container_type::iterator, + impl::dummy_iterator>; +using const_iterator + = impl::flat_iterator< + typename container_type::const_iterator, + iterator>; +using reverse_iterator + = impl::flat_iterator< + typename container_type::reverse_iterator, + impl::dummy_iterator>; +using const_reverse_iterator + = impl::flat_iterator< + typename container_type::const_reverse_iterator, + reverse_iterator>; diff --git a/external/flat/impl/flat_impl.hpp b/external/flat/impl/flat_impl.hpp new file mode 100644 index 000000000..2d639620e --- /dev/null +++ b/external/flat/impl/flat_impl.hpp @@ -0,0 +1,470 @@ +// Copyright Pubby 2016 +// Distributed under the Boost Software License, Version 1.0. +// http://www.boost.org/LICENSE_1_0.txt + +#ifndef LIB_FLAT_FLAT_IMPL_HPP +#define LIB_FLAT_FLAT_IMPL_HPP + +#include +#include +#include +#include +#include + +namespace fc { + +struct delay_sort_t {}; +constexpr delay_sort_t delay_sort = {}; + +struct container_construct_t {}; +constexpr container_construct_t container_construct = {}; + +namespace impl { + +template +using transparent_key_t = std::reference_wrapper; + +template +struct eq_comp +{ + template + bool operator()(A const& lhs, B const& rhs) + { return !comp(lhs, rhs) && !comp(rhs, lhs); } + Comp comp; +}; + + +template +struct dummy_iterator +{ + dummy_iterator() = delete; + dummy_iterator(dummy_iterator const&) = delete; + It underlying; +}; + +template::iterator_category> +class flat_iterator; + +template +class flat_iterator +{ + using traits = std::iterator_traits; +public: + using difference_type = typename traits::difference_type; + using value_type = typename traits::value_type const; + using pointer = value_type*; + using reference = value_type&; + using iterator_category = std::random_access_iterator_tag; + + flat_iterator() = default; + flat_iterator(flat_iterator const&) = default; + flat_iterator(flat_iterator&&) = default; + flat_iterator(Convert const& c) : underlying(c.underlying) {} + flat_iterator(Convert&& c) : underlying(std::move(c.underlying)) {} + flat_iterator(It const& underlying) : underlying(underlying) {} + flat_iterator(It&& underlying) : underlying(std::move(underlying)) {} + + flat_iterator& operator=(flat_iterator const& u) = default; + flat_iterator& operator=(flat_iterator&& u) = default; + flat_iterator& operator=(It const& u) + { this->underlying = u; return *this; } + flat_iterator& operator=(It&& u) + { this->underlying = std::move(u); return *this; } + + reference operator*() const { return *underlying; } + pointer operator->() const { return std::addressof(*underlying); } + + flat_iterator& operator++() { ++this->underlying; return *this; } + flat_iterator operator++(int) + { flat_iterator it = *this; ++this->underlying; return it; } + + flat_iterator& operator--() { --this->underlying; return *this; } + flat_iterator operator--(int) + { flat_iterator it = *this; --this->underlying; return it; } + + flat_iterator& operator+=(difference_type d) + { this->underlying += d; return *this; } + flat_iterator& operator-=(difference_type d) + { this->underlying -= d; return *this; } + + flat_iterator operator+(difference_type d) const + { return this->underlying + d; } + flat_iterator operator-(difference_type d) const + { return this->underlying - d; } + + difference_type operator-(flat_iterator const& o) const + { return this->underlying - o.underlying; } + + reference operator[](difference_type d) const { return *(*this + d); } + + auto operator==(flat_iterator const& o) const + { + using namespace std::rel_ops; + return this->underlying == o.underlying; + } + auto operator!=(flat_iterator const& o) const + { + using namespace std::rel_ops; + return this->underlying != o.underlying; + } + auto operator<(flat_iterator const& o) const + { return this->underlying < o.underlying; } + auto operator<=(flat_iterator const& o) const + { return this->underlying <= o.underlying; } + auto operator>(flat_iterator const& o) const + { return this->underlying > o.underlying; } + auto operator>=(flat_iterator const& o) const + { return this->underlying >= o.underlying; } + + It underlying; +}; + +template, typename TransparentKey = void> +struct first_compare +{ + first_compare(const Compare& comp) : + compare(comp) { + } + + bool operator()(Pair const& lhs, Pair const& rhs) const + { + return compare.get()(lhs.first, rhs.first); + } + + bool operator()(typename Pair::first_type const& lhs, Pair const& rhs) const + { + return compare.get()(lhs, rhs.first); + } + + bool operator()(Pair const& lhs, typename Pair::first_type const& rhs) const + { + return compare.get()(lhs.first, rhs); + } + + template + bool operator()(transparent_key_t const& lhs, Pair const& rhs) const + { + return compare.get()(lhs.get(), rhs.first); + } + + template + bool operator()(transparent_key_t const& lhs, typename Pair::first_type const& rhs) const + { + return compare.get()(lhs.get(), rhs); + } + + template + bool operator()(Pair const& lhs, transparent_key_t const& rhs) const + { + return compare.get()(lhs.first, rhs.get()); + } + + template + bool operator()(typename Pair::first_type const& lhs, transparent_key_t const& rhs) const + { + return compare.get()(lhs, rhs.get()); + } + + std::reference_wrapper compare; +}; + +template +class flat_container_base +{ +#include "container_traits.hpp" + D const* self() const { return static_cast(this); } + D* self() { return static_cast(this); } +public: + using key_compare = Compare; + const key_compare &key_comp() const { return self()->comp; } + + // Iterators + + const_iterator cbegin() const + noexcept(noexcept(std::declval().container.cbegin())) + { return self()->container.cbegin(); } + const_iterator begin() const + noexcept(noexcept(std::declval().container.begin())) + { return self()->container.begin(); } + iterator begin() + noexcept(noexcept(std::declval().container.begin())) + { return self()->container.begin(); } + + const_iterator cend() const + noexcept(noexcept(std::declval().container.cend())) + { return self()->container.cend(); } + const_iterator end() const + noexcept(noexcept(std::declval().container.end())) + { return self()->container.end(); } + iterator end() + noexcept(noexcept(std::declval().container.end())) + { return self()->container.end(); } + + const_reverse_iterator crbegin() const + noexcept(noexcept(std::declval().container.crbegin())) + { return self()->container.crbegin(); } + const_reverse_iterator rbegin() const + noexcept(noexcept(std::declval().container.rbegin())) + { return self()->container.rbegin(); } + reverse_iterator rbegin() + noexcept(noexcept(std::declval().container.rbegin())) + { return self()->container.rbegin(); } + + const_reverse_iterator crend() const + noexcept(noexcept(std::declval().container.crend())) + { return self()->container.crend(); } + const_reverse_iterator rend() const + noexcept(noexcept(std::declval().container.rend())) + { return self()->container.rend(); } + reverse_iterator rend() + noexcept(noexcept(std::declval().container.rend())) + { return self()->container.rend(); } + + // Capacity + + bool empty() const + noexcept(noexcept(std::declval().container.empty())) + { return self()->container.empty(); } + + size_type size() const + noexcept(noexcept(std::declval().container.size())) + { return self()->container.size(); } + + // Modifiers + + iterator insert(const_iterator hint, value_type const& value) + { return self()->insert(value).first; } + + iterator insert(const_iterator hint, value_type&& value) + { return self()->insert(std::move(value)).first; } + + template + void insert(InputIt first, InputIt last) + { + for(InputIt it = first; it != last; ++it) + self()->insert(*it); + } + + void insert(std::initializer_list ilist) + { self()->insert(ilist.begin(), ilist.end()); } + + void insert(std::initializer_list ilist, delay_sort_t d) + { self()->insert(ilist.begin(), ilist.end(), d); } + + template + auto emplace(Args&&... args) + { return self()->insert(value_type(std::forward(args)...)); } + + template + auto emplace_hint(const_iterator hint, Args&&... args) + { return self()->insert(value_type(std::forward(args)...)); } + + iterator erase(const_iterator pos) + noexcept(noexcept(std::declval().container.erase(pos.underlying))) + { return self()->container.erase(pos.underlying); } + + iterator erase(const_iterator first, const_iterator last) + noexcept(noexcept(std::declval().container.erase(first.underlying, + last.underlying))) + { return self()->container.erase(first.underlying, last.underlying); } + + void clear() + noexcept(noexcept(std::declval().container.clear())) + { self()->container.clear(); } + + void swap(D& other) + noexcept(D::has_noexcept_swap()) + { + using std::swap; + swap(self()->container, other.container); + } + + // Lookup + + const_iterator find(key_type const& key) const + { + const_iterator it = self()->lower_bound(key); + if(it == self()->end() || self()->value_comp()(key, *it)) + return self()->end(); + return it; + } + + iterator find(key_type const& key) + { + iterator it = self()->lower_bound(key); + if(it == self()->end() || self()->value_comp()(key, *it)) + return self()->end(); + return it; + } + + const_iterator lower_bound(key_type const& key) const + { + return std::lower_bound( + self()->begin(), self()->end(), + key, self()->value_comp()); + } + + iterator lower_bound(key_type const& key) + { + return std::lower_bound( + self()->begin(), self()->end(), + key, self()->value_comp()); + } + + const_iterator upper_bound(key_type const& key) const + { + return std::upper_bound( + self()->begin(), self()->end(), + key, self()->value_comp()); + } + + iterator upper_bound(key_type const& key) + { + return std::upper_bound( + self()->begin(), self()->end(), + key, self()->value_comp()); + } + + std::pair + equal_range(key_type const& key) const + { + return std::equal_range( + self()->begin(), self()->end(), + key, self()->value_comp()); + } + + std::pair equal_range(key_type const& key) + { + return std::equal_range( + self()->begin(), self()->end(), + key, self()->value_comp()); + } + +private: + static constexpr bool has_noexcept_swap() + { + using std::swap; + return noexcept(swap(*static_cast(nullptr), + *static_cast(nullptr))); + } + +protected: + template + void ds_insert_(InputIt first, InputIt last) + { + size_type const i = self()->size(); + for(InputIt it = first; it != last; ++it) + self()->container.push_back(*it); + std::sort( + self()->container.begin()+i, + self()->container.end(), + self()->value_comp()); + std::inplace_merge( + self()->container.begin(), + self()->container.begin()+i, + self()->container.end()); + // Note: Not calling unique here. Do it in the caller. + } +}; + +template +class flat_container_base> +: public flat_container_base +{ +#include "container_traits.hpp" + D const* self() const { return static_cast(this); } + D* self() { return static_cast(this); } + using B = flat_container_base; +public: + + using B::insert; + using B::find; + using B::lower_bound; + using B::upper_bound; + using B::equal_range; + + // Modifiers + + template + iterator insert(const_iterator hint, P&& value) + { return insert(std::forward

(value)).first; } + + // Lookup + + template + const_iterator find(K const& key) const + { + const_iterator it = self()->lower_bound(key); + if (it == self()->end() || self()->value_comp()(std::ref(key), *it)) + return self()->end(); + return it; + } + + template + iterator find(K const& key) + { + iterator it = self()->lower_bound(key); + if (it == self()->end() || self()->value_comp()(std::ref(key), *it)) + return self()->end(); + return it; + } + + template + const_iterator lower_bound(K const& key) const + { + return std::lower_bound( + self()->begin(), self()->end(), + std::ref(key), self()->value_comp()); + } + + template + iterator lower_bound(K const& key) + { + return std::lower_bound( + self()->begin(), self()->end(), + std::ref(key), self()->value_comp()); + } + + template + const_iterator upper_bound(K const& key) const + { + return std::upper_bound( + self()->begin(), self()->end(), + std::ref(key), self()->value_comp()); + } + + template + iterator upper_bound(K const& key) + { + return std::upper_bound( + self()->begin(), self()->end(), + std::ref(key), self()->value_comp()); + } + + template + std::pair + equal_range(K const& key) const + { + return std::equal_range( + self()->begin(), self()->end(), + std::ref(key), self()->value_comp()); + } + + template + std::pair equal_range(K const& key) + { + return std::equal_range( + self()->begin(), self()->end(), + std::ref(key), self()->value_comp()); + } +}; + +} // namespace fc +} // namespace impl + +#endif diff --git a/external/tsl/ordered_hash.h b/external/tsl/ordered_hash.h new file mode 100644 index 000000000..46d9b4877 --- /dev/null +++ b/external/tsl/ordered_hash.h @@ -0,0 +1,1617 @@ +/** + * MIT License + * + * Copyright (c) 2017 Thibaut Goetghebuer-Planchon + * + * 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 TSL_ORDERED_HASH_H +#define TSL_ORDERED_HASH_H + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +/** + * Macros for compatibility with GCC 4.8 + */ +#if (defined(__GNUC__) && (__GNUC__ == 4) && (__GNUC_MINOR__ < 9)) +# define TSL_OH_NO_CONTAINER_ERASE_CONST_ITERATOR +# define TSL_OH_NO_CONTAINER_EMPLACE_CONST_ITERATOR +#endif + +/** + * Only activate tsl_oh_assert if TSL_DEBUG is defined. + * This way we avoid the performance hit when NDEBUG is not defined with assert as tsl_oh_assert is used a lot + * (people usually compile with "-O3" and not "-O3 -DNDEBUG"). + */ +#ifdef TSL_DEBUG +# define tsl_oh_assert(expr) assert(expr) +#else +# define tsl_oh_assert(expr) (static_cast(0)) +#endif + +/** + * If exceptions are enabled, throw the exception passed in parameter, otherwise call std::terminate. + */ +#if (defined(__cpp_exceptions) || defined(__EXCEPTIONS) || (defined (_MSC_VER) && defined (_CPPUNWIND))) && !defined(TSL_NO_EXCEPTIONS) +# define TSL_OH_THROW_OR_TERMINATE(ex, msg) throw ex(msg) +#else +# define TSL_OH_NO_EXCEPTIONS +# ifdef NDEBUG +# define TSL_OH_THROW_OR_TERMINATE(ex, msg) std::terminate() +# else +# include +# define TSL_OH_THROW_OR_TERMINATE(ex, msg) do { std::cerr << msg << std::endl; std::terminate(); } while(0) +# endif +#endif + + +namespace tsl { + +namespace detail_ordered_hash { + +template +struct make_void { + using type = void; +}; + +template +struct has_is_transparent: std::false_type { +}; + +template +struct has_is_transparent::type>: std::true_type { +}; + + +template +struct is_vector: std::false_type { +}; + +template +struct is_vector>::value + >::type>: std::true_type { +}; + +// Only available in C++17, we need to be compatible with C++11 +template +const T& clamp( const T& v, const T& lo, const T& hi) { + return std::min(hi, std::max(lo, v)); +} + +template +static T numeric_cast(U value, const char* error_message = "numeric_cast() failed.") { + T ret = static_cast(value); + if(static_cast(ret) != value) { + TSL_OH_THROW_OR_TERMINATE(std::runtime_error, error_message); + } + + const bool is_same_signedness = (std::is_unsigned::value && std::is_unsigned::value) || + (std::is_signed::value && std::is_signed::value); + if(!is_same_signedness && (ret < T{}) != (value < U{})) { + TSL_OH_THROW_OR_TERMINATE(std::runtime_error, error_message); + } + + return ret; +} + + +/** + * Fixed size type used to represent size_type values on serialization. Need to be big enough + * to represent a std::size_t on 32 and 64 bits platforms, and must be the same size on both platforms. + */ +using slz_size_type = std::uint64_t; +static_assert(std::numeric_limits::max() >= std::numeric_limits::max(), + "slz_size_type must be >= std::size_t"); + +template +static T deserialize_value(Deserializer& deserializer) { + // MSVC < 2017 is not conformant, circumvent the problem by removing the template keyword +#if defined (_MSC_VER) && _MSC_VER < 1910 + return deserializer.Deserializer::operator()(); +#else + return deserializer.Deserializer::template operator()(); +#endif +} + + +/** + * Each bucket entry stores an index which is the index in m_values corresponding to the bucket's value + * and a hash (which may be truncated to 32 bits depending on IndexType) corresponding to the hash of the value. + * + * The size of IndexType limits the size of the hash table to std::numeric_limits::max() - 1 elements (-1 due to + * a reserved value used to mark a bucket as empty). + */ +template +class bucket_entry { + static_assert(std::is_unsigned::value, "IndexType must be an unsigned value."); + static_assert(std::numeric_limits::max() <= std::numeric_limits::max(), + "std::numeric_limits::max() must be <= std::numeric_limits::max()."); + +public: + using index_type = IndexType; + using truncated_hash_type = typename std::conditional::max() <= + std::numeric_limits::max(), + std::uint_least32_t, + std::size_t>::type; + + bucket_entry() noexcept: m_index(EMPTY_MARKER_INDEX), m_hash(0) { + } + + bool empty() const noexcept { + return m_index == EMPTY_MARKER_INDEX; + } + + void clear() noexcept { + m_index = EMPTY_MARKER_INDEX; + } + + index_type index() const noexcept { + tsl_oh_assert(!empty()); + return m_index; + } + + index_type& index_ref() noexcept { + tsl_oh_assert(!empty()); + return m_index; + } + + void set_index(index_type index) noexcept { + tsl_oh_assert(index <= max_size()); + + m_index = index; + } + + truncated_hash_type truncated_hash() const noexcept { + tsl_oh_assert(!empty()); + return m_hash; + } + + truncated_hash_type& truncated_hash_ref() noexcept { + tsl_oh_assert(!empty()); + return m_hash; + } + + void set_hash(std::size_t hash) noexcept { + m_hash = truncate_hash(hash); + } + + template + void serialize(Serializer& serializer) const { + const slz_size_type index = m_index; + serializer(index); + + const slz_size_type hash = m_hash; + serializer(hash); + } + + template + static bucket_entry deserialize(Deserializer& deserializer) { + const slz_size_type index = deserialize_value(deserializer); + const slz_size_type hash = deserialize_value(deserializer); + + bucket_entry bentry; + bentry.m_index = numeric_cast(index, "Deserialized index is too big."); + bentry.m_hash = numeric_cast(hash, "Deserialized hash is too big."); + + return bentry; + } + + + + static truncated_hash_type truncate_hash(std::size_t hash) noexcept { + return truncated_hash_type(hash); + } + + static std::size_t max_size() noexcept { + return static_cast(std::numeric_limits::max()) - NB_RESERVED_INDEXES; + } + +private: + static const index_type EMPTY_MARKER_INDEX = std::numeric_limits::max(); + static const std::size_t NB_RESERVED_INDEXES = 1; + + index_type m_index; + truncated_hash_type m_hash; +}; + + + +/** + * Internal common class used by ordered_map and ordered_set. + * + * ValueType is what will be stored by ordered_hash (usually std::pair for map and Key for set). + * + * KeySelect should be a FunctionObject which takes a ValueType in parameter and return a reference to the key. + * + * ValueSelect should be a FunctionObject which takes a ValueType in parameter and return a reference to the value. + * ValueSelect should be void if there is no value (in set for example). + * + * ValueTypeContainer is the container which will be used to store ValueType values. + * Usually a std::deque or std::vector. + * + * + * + * The ordered_hash structure is a hash table which preserves the order of insertion of the elements. + * To do so, it stores the values in the ValueTypeContainer (m_values) using emplace_back at each + * insertion of a new element. Another structure (m_buckets of type std::vector) will + * serve as buckets array for the hash table part. Each bucket stores an index which corresponds to + * the index in m_values where the bucket's value is and the (truncated) hash of this value. An index + * is used instead of a pointer to the value to reduce the size of each bucket entry. + * + * To resolve collisions in the buckets array, the structures use robin hood linear probing with + * backward shift deletion. + */ +template +class ordered_hash: private Hash, private KeyEqual { +private: + template + using has_mapped_type = typename std::integral_constant::value>; + + static_assert(std::is_same::value, + "ValueTypeContainer::value_type != ValueType. " + "Check that the ValueTypeContainer has 'Key' as type for a set or 'std::pair' as type for a map."); + + static_assert(std::is_same::value, + "ValueTypeContainer::allocator_type != Allocator. " + "Check that the allocator for ValueTypeContainer is the same as Allocator."); + + static_assert(std::is_same::value, + "Allocator::value_type != ValueType. " + "Check that the allocator has 'Key' as type for a set or 'std::pair' as type for a map."); + + +public: + template + class ordered_iterator; + + using key_type = typename KeySelect::key_type; + using value_type = ValueType; + using size_type = std::size_t; + using difference_type = std::ptrdiff_t; + using hasher = Hash; + using key_equal = KeyEqual; + using allocator_type = Allocator; + using reference = value_type&; + using const_reference = const value_type&; + using pointer = value_type*; + using const_pointer = const value_type*; + using iterator = ordered_iterator; + using const_iterator = ordered_iterator; + using reverse_iterator = std::reverse_iterator; + using const_reverse_iterator = std::reverse_iterator; + + using values_container_type = ValueTypeContainer; + +public: + template + class ordered_iterator { + friend class ordered_hash; + + private: + using iterator = typename std::conditional::type; + + + ordered_iterator(iterator it) noexcept: m_iterator(it) { + } + + public: + using iterator_category = std::random_access_iterator_tag; + using value_type = const typename ordered_hash::value_type; + using difference_type = typename iterator::difference_type; + using reference = value_type&; + using pointer = value_type*; + + + ordered_iterator() noexcept { + } + + // Copy constructor from iterator to const_iterator. + template::type* = nullptr> + ordered_iterator(const ordered_iterator& other) noexcept: m_iterator(other.m_iterator) { + } + + ordered_iterator(const ordered_iterator& other) = default; + ordered_iterator(ordered_iterator&& other) = default; + ordered_iterator& operator=(const ordered_iterator& other) = default; + ordered_iterator& operator=(ordered_iterator&& other) = default; + + const typename ordered_hash::key_type& key() const { + return KeySelect()(*m_iterator); + } + + template::value && IsConst>::type* = nullptr> + const typename U::value_type& value() const { + return U()(*m_iterator); + } + + template::value && !IsConst>::type* = nullptr> + typename U::value_type& value() { + return U()(*m_iterator); + } + + reference operator*() const { return *m_iterator; } + pointer operator->() const { return m_iterator.operator->(); } + + ordered_iterator& operator++() { ++m_iterator; return *this; } + ordered_iterator& operator--() { --m_iterator; return *this; } + + ordered_iterator operator++(int) { ordered_iterator tmp(*this); ++(*this); return tmp; } + ordered_iterator operator--(int) { ordered_iterator tmp(*this); --(*this); return tmp; } + + reference operator[](difference_type n) const { return m_iterator[n]; } + + ordered_iterator& operator+=(difference_type n) { m_iterator += n; return *this; } + ordered_iterator& operator-=(difference_type n) { m_iterator -= n; return *this; } + + ordered_iterator operator+(difference_type n) { ordered_iterator tmp(*this); tmp += n; return tmp; } + ordered_iterator operator-(difference_type n) { ordered_iterator tmp(*this); tmp -= n; return tmp; } + + friend bool operator==(const ordered_iterator& lhs, const ordered_iterator& rhs) { + return lhs.m_iterator == rhs.m_iterator; + } + + friend bool operator!=(const ordered_iterator& lhs, const ordered_iterator& rhs) { + return lhs.m_iterator != rhs.m_iterator; + } + + friend bool operator<(const ordered_iterator& lhs, const ordered_iterator& rhs) { + return lhs.m_iterator < rhs.m_iterator; + } + + friend bool operator>(const ordered_iterator& lhs, const ordered_iterator& rhs) { + return lhs.m_iterator > rhs.m_iterator; + } + + friend bool operator<=(const ordered_iterator& lhs, const ordered_iterator& rhs) { + return lhs.m_iterator <= rhs.m_iterator; + } + + friend bool operator>=(const ordered_iterator& lhs, const ordered_iterator& rhs) { + return lhs.m_iterator >= rhs.m_iterator; + } + + friend ordered_iterator operator+(difference_type n, const ordered_iterator& it) { + return n + it.m_iterator; + } + + friend difference_type operator-(const ordered_iterator& lhs, const ordered_iterator& rhs) { + return lhs.m_iterator - rhs.m_iterator; + } + + private: + iterator m_iterator; + }; + + +private: + using bucket_entry = tsl::detail_ordered_hash::bucket_entry; + + using buckets_container_allocator = typename + std::allocator_traits::template rebind_alloc; + + using buckets_container_type = std::vector; + + + using truncated_hash_type = typename bucket_entry::truncated_hash_type; + using index_type = typename bucket_entry::index_type; + +public: + ordered_hash(size_type bucket_count, + const Hash& hash, + const KeyEqual& equal, + const Allocator& alloc, + float max_load_factor): Hash(hash), + KeyEqual(equal), + m_buckets_data(alloc), + m_buckets(static_empty_bucket_ptr()), + m_hash_mask(0), + m_values(alloc), + m_grow_on_next_insert(false) + { + if(bucket_count > max_bucket_count()) { + TSL_OH_THROW_OR_TERMINATE(std::length_error, "The map exceeds its maximum size."); + } + + if(bucket_count > 0) { + bucket_count = round_up_to_power_of_two(bucket_count); + + m_buckets_data.resize(bucket_count); + m_buckets = m_buckets_data.data(), + m_hash_mask = bucket_count - 1; + } + + this->max_load_factor(max_load_factor); + } + + ordered_hash(const ordered_hash& other): Hash(other), + KeyEqual(other), + m_buckets_data(other.m_buckets_data), + m_buckets(m_buckets_data.empty()?static_empty_bucket_ptr(): + m_buckets_data.data()), + m_hash_mask(other.m_hash_mask), + m_values(other.m_values), + m_load_threshold(other.m_load_threshold), + m_max_load_factor(other.m_max_load_factor), + m_grow_on_next_insert(other.m_grow_on_next_insert) + { + } + + ordered_hash(ordered_hash&& other) noexcept(std::is_nothrow_move_constructible::value && + std::is_nothrow_move_constructible::value && + std::is_nothrow_move_constructible::value && + std::is_nothrow_move_constructible::value) + : Hash(std::move(static_cast(other))), + KeyEqual(std::move(static_cast(other))), + m_buckets_data(std::move(other.m_buckets_data)), + m_buckets(m_buckets_data.empty()?static_empty_bucket_ptr(): + m_buckets_data.data()), + m_hash_mask(other.m_hash_mask), + m_values(std::move(other.m_values)), + m_load_threshold(other.m_load_threshold), + m_max_load_factor(other.m_max_load_factor), + m_grow_on_next_insert(other.m_grow_on_next_insert) + { + other.m_buckets_data.clear(); + other.m_buckets = static_empty_bucket_ptr(); + other.m_hash_mask = 0; + other.m_values.clear(); + other.m_load_threshold = 0; + other.m_grow_on_next_insert = false; + } + + ordered_hash& operator=(const ordered_hash& other) { + if(&other != this) { + Hash::operator=(other); + KeyEqual::operator=(other); + + m_buckets_data = other.m_buckets_data; + m_buckets = m_buckets_data.empty()?static_empty_bucket_ptr(): + m_buckets_data.data(); + + m_hash_mask = other.m_hash_mask; + m_values = other.m_values; + m_load_threshold = other.m_load_threshold; + m_max_load_factor = other.m_max_load_factor; + m_grow_on_next_insert = other.m_grow_on_next_insert; + } + + return *this; + } + + ordered_hash& operator=(ordered_hash&& other) { + other.swap(*this); + other.clear(); + + return *this; + } + + allocator_type get_allocator() const { + return m_values.get_allocator(); + } + + + /* + * Iterators + */ + iterator begin() noexcept { + return iterator(m_values.begin()); + } + + const_iterator begin() const noexcept { + return cbegin(); + } + + const_iterator cbegin() const noexcept { + return const_iterator(m_values.cbegin()); + } + + iterator end() noexcept { + return iterator(m_values.end()); + } + + const_iterator end() const noexcept { + return cend(); + } + + const_iterator cend() const noexcept { + return const_iterator(m_values.cend()); + } + + + reverse_iterator rbegin() noexcept { + return reverse_iterator(m_values.end()); + } + + const_reverse_iterator rbegin() const noexcept { + return rcbegin(); + } + + const_reverse_iterator rcbegin() const noexcept { + return const_reverse_iterator(m_values.cend()); + } + + reverse_iterator rend() noexcept { + return reverse_iterator(m_values.begin()); + } + + const_reverse_iterator rend() const noexcept { + return rcend(); + } + + const_reverse_iterator rcend() const noexcept { + return const_reverse_iterator(m_values.cbegin()); + } + + + /* + * Capacity + */ + bool empty() const noexcept { + return m_values.empty(); + } + + size_type size() const noexcept { + return m_values.size(); + } + + size_type max_size() const noexcept { + return std::min(bucket_entry::max_size(), m_values.max_size()); + } + + + /* + * Modifiers + */ + void clear() noexcept { + for(auto& bucket: m_buckets_data) { + bucket.clear(); + } + + m_values.clear(); + m_grow_on_next_insert = false; + } + + template + std::pair insert(P&& value) { + return insert_impl(KeySelect()(value), std::forward

(value)); + } + + template + iterator insert_hint(const_iterator hint, P&& value) { + if(hint != cend() && compare_keys(KeySelect()(*hint), KeySelect()(value))) { + return mutable_iterator(hint); + } + + return insert(std::forward

(value)).first; + } + + template + void insert(InputIt first, InputIt last) { + if(std::is_base_of::iterator_category>::value) + { + const auto nb_elements_insert = std::distance(first, last); + const size_type nb_free_buckets = m_load_threshold - size(); + tsl_oh_assert(m_load_threshold >= size()); + + if(nb_elements_insert > 0 && nb_free_buckets < size_type(nb_elements_insert)) { + reserve(size() + size_type(nb_elements_insert)); + } + } + + for(; first != last; ++first) { + insert(*first); + } + } + + + + template + std::pair insert_or_assign(K&& key, M&& value) { + auto it = try_emplace(std::forward(key), std::forward(value)); + if(!it.second) { + it.first.value() = std::forward(value); + } + + return it; + } + + template + iterator insert_or_assign(const_iterator hint, K&& key, M&& obj) { + if(hint != cend() && compare_keys(KeySelect()(*hint), key)) { + auto it = mutable_iterator(hint); + it.value() = std::forward(obj); + + return it; + } + + return insert_or_assign(std::forward(key), std::forward(obj)).first; + } + + + + template + std::pair emplace(Args&&... args) { + return insert(value_type(std::forward(args)...)); + } + + template + iterator emplace_hint(const_iterator hint, Args&&... args) { + return insert_hint(hint, value_type(std::forward(args)...)); + } + + + + template + std::pair try_emplace(K&& key, Args&&... value_args) { + return insert_impl(key, std::piecewise_construct, + std::forward_as_tuple(std::forward(key)), + std::forward_as_tuple(std::forward(value_args)...)); + } + + template + iterator try_emplace_hint(const_iterator hint, K&& key, Args&&... args) { + if(hint != cend() && compare_keys(KeySelect()(*hint), key)) { + return mutable_iterator(hint); + } + + return try_emplace(std::forward(key), std::forward(args)...).first; + } + + + + /** + * Here to avoid `template size_type erase(const K& key)` being used when + * we use an `iterator` instead of a `const_iterator`. + */ + iterator erase(iterator pos) { + return erase(const_iterator(pos)); + } + + iterator erase(const_iterator pos) { + tsl_oh_assert(pos != cend()); + + const std::size_t index_erase = iterator_to_index(pos); + + auto it_bucket = find_key(pos.key(), hash_key(pos.key())); + tsl_oh_assert(it_bucket != m_buckets_data.end()); + + erase_value_from_bucket(it_bucket); + + /* + * One element was removed from m_values, due to the left shift the next element + * is now at the position of the previous element (or end if none). + */ + return begin() + index_erase; + } + + iterator erase(const_iterator first, const_iterator last) { + if(first == last) { + return mutable_iterator(first); + } + + tsl_oh_assert(std::distance(first, last) > 0); + const std::size_t start_index = iterator_to_index(first); + const std::size_t nb_values = std::size_t(std::distance(first, last)); + const std::size_t end_index = start_index + nb_values; + + // Delete all values +#ifdef TSL_OH_NO_CONTAINER_ERASE_CONST_ITERATOR + auto next_it = m_values.erase(mutable_iterator(first).m_iterator, mutable_iterator(last).m_iterator); +#else + auto next_it = m_values.erase(first.m_iterator, last.m_iterator); +#endif + + /* + * Mark the buckets corresponding to the values as empty and do a backward shift. + * + * Also, the erase operation on m_values has shifted all the values on the right of last.m_iterator. + * Adapt the indexes for these values. + */ + std::size_t ibucket = 0; + while(ibucket < m_buckets_data.size()) { + if(m_buckets[ibucket].empty()) { + ibucket++; + } + else if(m_buckets[ibucket].index() >= start_index && m_buckets[ibucket].index() < end_index) { + m_buckets[ibucket].clear(); + backward_shift(ibucket); + // Don't increment ibucket, backward_shift may have replaced current bucket. + } + else if(m_buckets[ibucket].index() >= end_index) { + m_buckets[ibucket].set_index(index_type(m_buckets[ibucket].index() - nb_values)); + ibucket++; + } + else { + ibucket++; + } + } + + return iterator(next_it); + } + + + template + size_type erase(const K& key) { + return erase(key, hash_key(key)); + } + + template + size_type erase(const K& key, std::size_t hash) { + return erase_impl(key, hash); + } + + void swap(ordered_hash& other) { + using std::swap; + + swap(static_cast(*this), static_cast(other)); + swap(static_cast(*this), static_cast(other)); + swap(m_buckets_data, other.m_buckets_data); + swap(m_buckets, other.m_buckets); + swap(m_hash_mask, other.m_hash_mask); + swap(m_values, other.m_values); + swap(m_load_threshold, other.m_load_threshold); + swap(m_max_load_factor, other.m_max_load_factor); + swap(m_grow_on_next_insert, other.m_grow_on_next_insert); + } + + + + + /* + * Lookup + */ + template::value>::type* = nullptr> + typename U::value_type& at(const K& key) { + return at(key, hash_key(key)); + } + + template::value>::type* = nullptr> + typename U::value_type& at(const K& key, std::size_t hash) { + return const_cast(static_cast(this)->at(key, hash)); + } + + template::value>::type* = nullptr> + const typename U::value_type& at(const K& key) const { + return at(key, hash_key(key)); + } + + template::value>::type* = nullptr> + const typename U::value_type& at(const K& key, std::size_t hash) const { + auto it = find(key, hash); + if(it != end()) { + return it.value(); + } + else { + TSL_OH_THROW_OR_TERMINATE(std::out_of_range, "Couldn't find the key."); + } + } + + + template::value>::type* = nullptr> + typename U::value_type& operator[](K&& key) { + return try_emplace(std::forward(key)).first.value(); + } + + + template + size_type count(const K& key) const { + return count(key, hash_key(key)); + } + + template + size_type count(const K& key, std::size_t hash) const { + if(find(key, hash) == cend()) { + return 0; + } + else { + return 1; + } + } + + template + iterator find(const K& key) { + return find(key, hash_key(key)); + } + + template + iterator find(const K& key, std::size_t hash) { + auto it_bucket = find_key(key, hash); + return (it_bucket != m_buckets_data.end())?iterator(m_values.begin() + it_bucket->index()):end(); + } + + template + const_iterator find(const K& key) const { + return find(key, hash_key(key)); + } + + template + const_iterator find(const K& key, std::size_t hash) const { + auto it_bucket = find_key(key, hash); + return (it_bucket != m_buckets_data.cend())?const_iterator(m_values.begin() + it_bucket->index()):end(); + } + + + template + std::pair equal_range(const K& key) { + return equal_range(key, hash_key(key)); + } + + template + std::pair equal_range(const K& key, std::size_t hash) { + iterator it = find(key, hash); + return std::make_pair(it, (it == end())?it:std::next(it)); + } + + template + std::pair equal_range(const K& key) const { + return equal_range(key, hash_key(key)); + } + + template + std::pair equal_range(const K& key, std::size_t hash) const { + const_iterator it = find(key, hash); + return std::make_pair(it, (it == cend())?it:std::next(it)); + } + + + /* + * Bucket interface + */ + size_type bucket_count() const { + return m_buckets_data.size(); + } + + size_type max_bucket_count() const { + return m_buckets_data.max_size(); + } + + /* + * Hash policy + */ + float load_factor() const { + if(bucket_count() == 0) { + return 0; + } + + return float(size())/float(bucket_count()); + } + + float max_load_factor() const { + return m_max_load_factor; + } + + void max_load_factor(float ml) { + m_max_load_factor = clamp(ml, float(MAX_LOAD_FACTOR__MINIMUM), + float(MAX_LOAD_FACTOR__MAXIMUM)); + + m_max_load_factor = ml; + m_load_threshold = size_type(float(bucket_count())*m_max_load_factor); + } + + void rehash(size_type count) { + count = std::max(count, size_type(std::ceil(float(size())/max_load_factor()))); + rehash_impl(count); + } + + void reserve(size_type count) { + reserve_space_for_values(count); + + count = size_type(std::ceil(float(count)/max_load_factor())); + rehash(count); + } + + + /* + * Observers + */ + hasher hash_function() const { + return static_cast(*this); + } + + key_equal key_eq() const { + return static_cast(*this); + } + + + /* + * Other + */ + iterator mutable_iterator(const_iterator pos) { + return iterator(m_values.begin() + iterator_to_index(pos)); + } + + iterator nth(size_type index) { + tsl_oh_assert(index <= size()); + return iterator(m_values.begin() + index); + } + + const_iterator nth(size_type index) const { + tsl_oh_assert(index <= size()); + return const_iterator(m_values.cbegin() + index); + } + + const_reference front() const { + tsl_oh_assert(!empty()); + return m_values.front(); + } + + const_reference back() const { + tsl_oh_assert(!empty()); + return m_values.back(); + } + + const values_container_type& values_container() const noexcept { + return m_values; + } + + template::value>::type* = nullptr> + const typename values_container_type::value_type* data() const noexcept { + return m_values.data(); + } + + template::value>::type* = nullptr> + size_type capacity() const noexcept { + return m_values.capacity(); + } + + void shrink_to_fit() { + m_values.shrink_to_fit(); + } + + + template + std::pair insert_at_position(const_iterator pos, P&& value) { + return insert_at_position_impl(pos.m_iterator, KeySelect()(value), std::forward

(value)); + } + + template + std::pair emplace_at_position(const_iterator pos, Args&&... args) { + return insert_at_position(pos, value_type(std::forward(args)...)); + } + + template + std::pair try_emplace_at_position(const_iterator pos, K&& key, Args&&... value_args) { + return insert_at_position_impl(pos.m_iterator, key, + std::piecewise_construct, + std::forward_as_tuple(std::forward(key)), + std::forward_as_tuple(std::forward(value_args)...)); + } + + + void pop_back() { + tsl_oh_assert(!empty()); + erase(std::prev(end())); + } + + + /** + * Here to avoid `template size_type unordered_erase(const K& key)` being used when + * we use a iterator instead of a const_iterator. + */ + iterator unordered_erase(iterator pos) { + return unordered_erase(const_iterator(pos)); + } + + iterator unordered_erase(const_iterator pos) { + const std::size_t index_erase = iterator_to_index(pos); + unordered_erase(pos.key()); + + /* + * One element was deleted, index_erase now points to the next element as the elements after + * the deleted value were shifted to the left in m_values (will be end() if we deleted the last element). + */ + return begin() + index_erase; + } + + template + size_type unordered_erase(const K& key) { + return unordered_erase(key, hash_key(key)); + } + + template + size_type unordered_erase(const K& key, std::size_t hash) { + auto it_bucket_key = find_key(key, hash); + if(it_bucket_key == m_buckets_data.end()) { + return 0; + } + + /** + * If we are not erasing the last element in m_values, we swap + * the element we are erasing with the last element. We then would + * just have to do a pop_back() in m_values. + */ + if(!compare_keys(key, KeySelect()(back()))) { + auto it_bucket_last_elem = find_key(KeySelect()(back()), hash_key(KeySelect()(back()))); + tsl_oh_assert(it_bucket_last_elem != m_buckets_data.end()); + tsl_oh_assert(it_bucket_last_elem->index() == m_values.size() - 1); + + using std::swap; + swap(m_values[it_bucket_key->index()], m_values[it_bucket_last_elem->index()]); + swap(it_bucket_key->index_ref(), it_bucket_last_elem->index_ref()); + } + + erase_value_from_bucket(it_bucket_key); + + return 1; + } + + template + void serialize(Serializer& serializer) const { + serialize_impl(serializer); + } + + template + void deserialize(Deserializer& deserializer, bool hash_compatible) { + deserialize_impl(deserializer, hash_compatible); + } + + friend bool operator==(const ordered_hash& lhs, const ordered_hash& rhs) { + return lhs.m_values == rhs.m_values; + } + + friend bool operator!=(const ordered_hash& lhs, const ordered_hash& rhs) { + return lhs.m_values != rhs.m_values; + } + + friend bool operator<(const ordered_hash& lhs, const ordered_hash& rhs) { + return lhs.m_values < rhs.m_values; + } + + friend bool operator<=(const ordered_hash& lhs, const ordered_hash& rhs) { + return lhs.m_values <= rhs.m_values; + } + + friend bool operator>(const ordered_hash& lhs, const ordered_hash& rhs) { + return lhs.m_values > rhs.m_values; + } + + friend bool operator>=(const ordered_hash& lhs, const ordered_hash& rhs) { + return lhs.m_values >= rhs.m_values; + } + + +private: + template + std::size_t hash_key(const K& key) const { + return Hash::operator()(key); + } + + template + bool compare_keys(const K1& key1, const K2& key2) const { + return KeyEqual::operator()(key1, key2); + } + + template + typename buckets_container_type::iterator find_key(const K& key, std::size_t hash) { + auto it = static_cast(this)->find_key(key, hash); + return m_buckets_data.begin() + std::distance(m_buckets_data.cbegin(), it); + } + + /** + * Return bucket which has the key 'key' or m_buckets_data.end() if none. + * + * From the bucket_for_hash, search for the value until we either find an empty bucket + * or a bucket which has a value with a distance from its ideal bucket longer + * than the probe length for the value we are looking for. + */ + template + typename buckets_container_type::const_iterator find_key(const K& key, std::size_t hash) const { + for(std::size_t ibucket = bucket_for_hash(hash), dist_from_ideal_bucket = 0; ; + ibucket = next_bucket(ibucket), dist_from_ideal_bucket++) + { + if(m_buckets[ibucket].empty()) { + return m_buckets_data.end(); + } + else if(m_buckets[ibucket].truncated_hash() == bucket_entry::truncate_hash(hash) && + compare_keys(key, KeySelect()(m_values[m_buckets[ibucket].index()]))) + { + return m_buckets_data.begin() + ibucket; + } + else if(dist_from_ideal_bucket > distance_from_ideal_bucket(ibucket)) { + return m_buckets_data.end(); + } + } + } + + void rehash_impl(size_type bucket_count) { + tsl_oh_assert(bucket_count >= size_type(std::ceil(float(size())/max_load_factor()))); + + if(bucket_count > max_bucket_count()) { + TSL_OH_THROW_OR_TERMINATE(std::length_error, "The map exceeds its maximum size."); + } + + if(bucket_count > 0) { + bucket_count = round_up_to_power_of_two(bucket_count); + } + + if(bucket_count == this->bucket_count()) { + return; + } + + + buckets_container_type old_buckets(bucket_count); + m_buckets_data.swap(old_buckets); + m_buckets = m_buckets_data.empty()?static_empty_bucket_ptr(): + m_buckets_data.data(); + // Everything should be noexcept from here. + + m_hash_mask = (bucket_count > 0)?(bucket_count - 1):0; + this->max_load_factor(m_max_load_factor); + m_grow_on_next_insert = false; + + + + for(const bucket_entry& old_bucket: old_buckets) { + if(old_bucket.empty()) { + continue; + } + + truncated_hash_type insert_hash = old_bucket.truncated_hash(); + index_type insert_index = old_bucket.index(); + + for(std::size_t ibucket = bucket_for_hash(insert_hash), dist_from_ideal_bucket = 0; ; + ibucket = next_bucket(ibucket), dist_from_ideal_bucket++) + { + if(m_buckets[ibucket].empty()) { + m_buckets[ibucket].set_index(insert_index); + m_buckets[ibucket].set_hash(insert_hash); + break; + } + + const std::size_t distance = distance_from_ideal_bucket(ibucket); + if(dist_from_ideal_bucket > distance) { + std::swap(insert_index, m_buckets[ibucket].index_ref()); + std::swap(insert_hash, m_buckets[ibucket].truncated_hash_ref()); + dist_from_ideal_bucket = distance; + } + } + } + } + + template::value>::type* = nullptr> + void reserve_space_for_values(size_type count) { + m_values.reserve(count); + } + + template::value>::type* = nullptr> + void reserve_space_for_values(size_type /*count*/) { + } + + /** + * Swap the empty bucket with the values on its right until we cross another empty bucket + * or if the other bucket has a distance_from_ideal_bucket == 0. + */ + void backward_shift(std::size_t empty_ibucket) noexcept { + tsl_oh_assert(m_buckets[empty_ibucket].empty()); + + std::size_t previous_ibucket = empty_ibucket; + for(std::size_t current_ibucket = next_bucket(previous_ibucket); + !m_buckets[current_ibucket].empty() && distance_from_ideal_bucket(current_ibucket) > 0; + previous_ibucket = current_ibucket, current_ibucket = next_bucket(current_ibucket)) + { + std::swap(m_buckets[current_ibucket], m_buckets[previous_ibucket]); + } + } + + void erase_value_from_bucket(typename buckets_container_type::iterator it_bucket) { + tsl_oh_assert(it_bucket != m_buckets_data.end() && !it_bucket->empty()); + + m_values.erase(m_values.begin() + it_bucket->index()); + + /* + * m_values.erase shifted all the values on the right of the erased value, + * shift the indexes by -1 in the buckets array for these values. + */ + if(it_bucket->index() != m_values.size()) { + shift_indexes_in_buckets(it_bucket->index(), -1); + } + + // Mark the bucket as empty and do a backward shift of the values on the right + it_bucket->clear(); + backward_shift(std::size_t(std::distance(m_buckets_data.begin(), it_bucket))); + } + + /** + * Go through each value from [from_ivalue, m_values.size()) in m_values and for each + * bucket corresponding to the value, shift the index by delta. + * + * delta must be equal to 1 or -1. + */ + void shift_indexes_in_buckets(index_type from_ivalue, int delta) noexcept { + tsl_oh_assert(delta == 1 || delta == -1); + + for(std::size_t ivalue = from_ivalue; ivalue < m_values.size(); ivalue++) { + // All the values in m_values have been shifted by delta. Find the bucket corresponding + // to the value m_values[ivalue] + const index_type old_index = static_cast(ivalue - delta); + + std::size_t ibucket = bucket_for_hash(hash_key(KeySelect()(m_values[ivalue]))); + while(m_buckets[ibucket].index() != old_index) { + ibucket = next_bucket(ibucket); + } + + m_buckets[ibucket].set_index(index_type(ivalue)); + } + } + + template + size_type erase_impl(const K& key, std::size_t hash) { + auto it_bucket = find_key(key, hash); + if(it_bucket != m_buckets_data.end()) { + erase_value_from_bucket(it_bucket); + + return 1; + } + else { + return 0; + } + } + + /** + * Insert the element at the end. + */ + template + std::pair insert_impl(const K& key, Args&&... value_type_args) { + const std::size_t hash = hash_key(key); + + std::size_t ibucket = bucket_for_hash(hash); + std::size_t dist_from_ideal_bucket = 0; + + while(!m_buckets[ibucket].empty() && dist_from_ideal_bucket <= distance_from_ideal_bucket(ibucket)) { + if(m_buckets[ibucket].truncated_hash() == bucket_entry::truncate_hash(hash) && + compare_keys(key, KeySelect()(m_values[m_buckets[ibucket].index()]))) + { + return std::make_pair(begin() + m_buckets[ibucket].index(), false); + } + + ibucket = next_bucket(ibucket); + dist_from_ideal_bucket++; + } + + if(size() >= max_size()) { + TSL_OH_THROW_OR_TERMINATE(std::length_error, "We reached the maximum size for the hash table."); + } + + + if(grow_on_high_load()) { + ibucket = bucket_for_hash(hash); + dist_from_ideal_bucket = 0; + } + + + m_values.emplace_back(std::forward(value_type_args)...); + insert_index(ibucket, dist_from_ideal_bucket, + index_type(m_values.size() - 1), bucket_entry::truncate_hash(hash)); + + + return std::make_pair(std::prev(end()), true); + } + + /** + * Insert the element before insert_position. + */ + template + std::pair insert_at_position_impl(typename values_container_type::const_iterator insert_position, + const K& key, Args&&... value_type_args) + { + const std::size_t hash = hash_key(key); + + std::size_t ibucket = bucket_for_hash(hash); + std::size_t dist_from_ideal_bucket = 0; + + while(!m_buckets[ibucket].empty() && dist_from_ideal_bucket <= distance_from_ideal_bucket(ibucket)) { + if(m_buckets[ibucket].truncated_hash() == bucket_entry::truncate_hash(hash) && + compare_keys(key, KeySelect()(m_values[m_buckets[ibucket].index()]))) + { + return std::make_pair(begin() + m_buckets[ibucket].index(), false); + } + + ibucket = next_bucket(ibucket); + dist_from_ideal_bucket++; + } + + if(size() >= max_size()) { + TSL_OH_THROW_OR_TERMINATE(std::length_error, "We reached the maximum size for the hash table."); + } + + + if(grow_on_high_load()) { + ibucket = bucket_for_hash(hash); + dist_from_ideal_bucket = 0; + } + + + const index_type index_insert_position = index_type(std::distance(m_values.cbegin(), insert_position)); + +#ifdef TSL_OH_NO_CONTAINER_EMPLACE_CONST_ITERATOR + m_values.emplace(m_values.begin() + std::distance(m_values.cbegin(), insert_position), std::forward(value_type_args)...); +#else + m_values.emplace(insert_position, std::forward(value_type_args)...); +#endif + + insert_index(ibucket, dist_from_ideal_bucket, + index_insert_position, bucket_entry::truncate_hash(hash)); + + /* + * The insertion didn't happend at the end of the m_values container, + * we need to shift the indexes in m_buckets_data. + */ + if(index_insert_position != m_values.size() - 1) { + shift_indexes_in_buckets(index_insert_position + 1, 1); + } + + return std::make_pair(iterator(m_values.begin() + index_insert_position), true); + } + + void insert_index(std::size_t ibucket, std::size_t dist_from_ideal_bucket, + index_type index_insert, truncated_hash_type hash_insert) noexcept + { + while(!m_buckets[ibucket].empty()) { + const std::size_t distance = distance_from_ideal_bucket(ibucket); + if(dist_from_ideal_bucket > distance) { + std::swap(index_insert, m_buckets[ibucket].index_ref()); + std::swap(hash_insert, m_buckets[ibucket].truncated_hash_ref()); + + dist_from_ideal_bucket = distance; + } + + + ibucket = next_bucket(ibucket); + dist_from_ideal_bucket++; + + + if(dist_from_ideal_bucket > REHASH_ON_HIGH_NB_PROBES__NPROBES && !m_grow_on_next_insert && + load_factor() >= REHASH_ON_HIGH_NB_PROBES__MIN_LOAD_FACTOR) + { + // We don't want to grow the map now as we need this method to be noexcept. + // Do it on next insert. + m_grow_on_next_insert = true; + } + } + + + m_buckets[ibucket].set_index(index_insert); + m_buckets[ibucket].set_hash(hash_insert); + } + + std::size_t distance_from_ideal_bucket(std::size_t ibucket) const noexcept { + const std::size_t ideal_bucket = bucket_for_hash(m_buckets[ibucket].truncated_hash()); + + if(ibucket >= ideal_bucket) { + return ibucket - ideal_bucket; + } + // If the bucket is smaller than the ideal bucket for the value, there was a wrapping at the end of the + // bucket array due to the modulo. + else { + return (bucket_count() + ibucket) - ideal_bucket; + } + } + + std::size_t next_bucket(std::size_t index) const noexcept { + tsl_oh_assert(index < m_buckets_data.size()); + + index++; + return (index < m_buckets_data.size())?index:0; + } + + std::size_t bucket_for_hash(std::size_t hash) const noexcept { + return hash & m_hash_mask; + } + + std::size_t iterator_to_index(const_iterator it) const noexcept { + const auto dist = std::distance(cbegin(), it); + tsl_oh_assert(dist >= 0); + + return std::size_t(dist); + } + + /** + * Return true if the map has been rehashed. + */ + bool grow_on_high_load() { + if(m_grow_on_next_insert || size() >= m_load_threshold) { + rehash_impl(std::max(size_type(1), bucket_count() * 2)); + m_grow_on_next_insert = false; + + return true; + } + else { + return false; + } + } + + template + void serialize_impl(Serializer& serializer) const { + const slz_size_type version = SERIALIZATION_PROTOCOL_VERSION; + serializer(version); + + const slz_size_type nb_elements = m_values.size(); + serializer(nb_elements); + + const slz_size_type bucket_count = m_buckets_data.size(); + serializer(bucket_count); + + const float max_load_factor = m_max_load_factor; + serializer(max_load_factor); + + + for(const value_type& value: m_values) { + serializer(value); + } + + for(const bucket_entry& bucket: m_buckets_data) { + bucket.serialize(serializer); + } + } + + template + void deserialize_impl(Deserializer& deserializer, bool hash_compatible) { + tsl_oh_assert(m_buckets_data.empty()); // Current hash table must be empty + + const slz_size_type version = deserialize_value(deserializer); + // For now we only have one version of the serialization protocol. + // If it doesn't match there is a problem with the file. + if(version != SERIALIZATION_PROTOCOL_VERSION) { + TSL_OH_THROW_OR_TERMINATE(std::runtime_error, "Can't deserialize the ordered_map/set. " + "The protocol version header is invalid."); + } + + const slz_size_type nb_elements = deserialize_value(deserializer); + const slz_size_type bucket_count_ds = deserialize_value(deserializer); + const float max_load_factor = deserialize_value(deserializer); + + if(max_load_factor < MAX_LOAD_FACTOR__MINIMUM || max_load_factor > MAX_LOAD_FACTOR__MAXIMUM) { + TSL_OH_THROW_OR_TERMINATE(std::runtime_error, "Invalid max_load_factor. Check that the serializer " + "and deserializer support floats correctly as they " + "can be converted implicitly to ints."); + } + + + this->max_load_factor(max_load_factor); + + if(bucket_count_ds == 0) { + tsl_oh_assert(nb_elements == 0); + return; + } + + + if(!hash_compatible) { + reserve(numeric_cast(nb_elements, "Deserialized nb_elements is too big.")); + for(slz_size_type el = 0; el < nb_elements; el++) { + insert(deserialize_value(deserializer)); + } + } + else { + m_buckets_data.reserve(numeric_cast(bucket_count_ds, "Deserialized bucket_count is too big.")); + m_buckets = m_buckets_data.data(), + m_hash_mask = m_buckets_data.capacity() - 1; + + reserve_space_for_values(numeric_cast(nb_elements, "Deserialized nb_elements is too big.")); + for(slz_size_type el = 0; el < nb_elements; el++) { + m_values.push_back(deserialize_value(deserializer)); + } + + for(slz_size_type b = 0; b < bucket_count_ds; b++) { + m_buckets_data.push_back(bucket_entry::deserialize(deserializer)); + } + } + } + + static std::size_t round_up_to_power_of_two(std::size_t value) { + if(is_power_of_two(value)) { + return value; + } + + if(value == 0) { + return 1; + } + + --value; + for(std::size_t i = 1; i < sizeof(std::size_t) * CHAR_BIT; i *= 2) { + value |= value >> i; + } + + return value + 1; + } + + static constexpr bool is_power_of_two(std::size_t value) { + return value != 0 && (value & (value - 1)) == 0; + } + + +public: + static const size_type DEFAULT_INIT_BUCKETS_SIZE = 0; + static constexpr float DEFAULT_MAX_LOAD_FACTOR = 0.75f; + +private: + static constexpr float MAX_LOAD_FACTOR__MINIMUM = 0.1f; + static constexpr float MAX_LOAD_FACTOR__MAXIMUM = 0.95f; + + static const size_type REHASH_ON_HIGH_NB_PROBES__NPROBES = 128; + static constexpr float REHASH_ON_HIGH_NB_PROBES__MIN_LOAD_FACTOR = 0.15f; + + /** + * Protocol version currenlty used for serialization. + */ + static const slz_size_type SERIALIZATION_PROTOCOL_VERSION = 1; + + /** + * Return an always valid pointer to an static empty bucket_entry with last_bucket() == true. + */ + bucket_entry* static_empty_bucket_ptr() { + static bucket_entry empty_bucket; + return &empty_bucket; + } + +private: + buckets_container_type m_buckets_data; + + /** + * Points to m_buckets_data.data() if !m_buckets_data.empty() otherwise points to static_empty_bucket_ptr. + * This variable is useful to avoid the cost of checking if m_buckets_data is empty when trying + * to find an element. + * + * TODO Remove m_buckets_data and only use a pointer+size instead of a pointer+vector to save some space in the ordered_hash object. + */ + bucket_entry* m_buckets; + + size_type m_hash_mask; + + values_container_type m_values; + + size_type m_load_threshold; + float m_max_load_factor; + + bool m_grow_on_next_insert; +}; + + +} // end namespace detail_ordered_hash + +} // end namespace tsl + +#endif diff --git a/external/tsl/ordered_map.h b/external/tsl/ordered_map.h new file mode 100644 index 000000000..c8b8b6bb2 --- /dev/null +++ b/external/tsl/ordered_map.h @@ -0,0 +1,833 @@ +/** + * MIT License + * + * Copyright (c) 2017 Thibaut Goetghebuer-Planchon + * + * 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 TSL_ORDERED_MAP_H +#define TSL_ORDERED_MAP_H + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "ordered_hash.h" + + +namespace tsl { + + +/** + * Implementation of an hash map using open addressing with robin hood with backshift delete to resolve collisions. + * + * The particularity of this hash map is that it remembers the order in which the elements were added and + * provide a way to access the structure which stores these values through the 'values_container()' method. + * The used container is defined by ValueTypeContainer, by default a std::deque is used (grows faster) but + * a std::vector may be used. In this case the map provides a 'data()' method which give a direct access + * to the memory used to store the values (which can be useful to communicate with C API's). + * + * The Key and T must be copy constructible and/or move constructible. To use `unordered_erase` they both + * must be swappable. + * + * The behaviour of the hash map is undefined if the destructor of Key or T throws an exception. + * + * By default the maximum size of a map is limited to 2^32 - 1 values, if needed this can be changed through + * the IndexType template parameter. Using an `uint64_t` will raise this limit to 2^64 - 1 values but each + * bucket will use 16 bytes instead of 8 bytes in addition to the space needed to store the values. + * + * Iterators invalidation: + * - clear, operator=, reserve, rehash: always invalidate the iterators (also invalidate end()). + * - insert, emplace, emplace_hint, operator[]: when a std::vector is used as ValueTypeContainer + * and if size() < capacity(), only end(). + * Otherwise all the iterators are invalidated if an insert occurs. + * - erase, unordered_erase: when a std::vector is used as ValueTypeContainer invalidate the iterator of + * the erased element and all the ones after the erased element (including end()). + * Otherwise all the iterators are invalidated if an erase occurs. + */ +template, + class KeyEqual = std::equal_to, + class Allocator = std::allocator>, + class ValueTypeContainer = std::deque, Allocator>, + class IndexType = std::uint_least32_t> +class ordered_map { +private: + template + using has_is_transparent = tsl::detail_ordered_hash::has_is_transparent; + + class KeySelect { + public: + using key_type = Key; + + const key_type& operator()(const std::pair& key_value) const noexcept { + return key_value.first; + } + + key_type& operator()(std::pair& key_value) noexcept { + return key_value.first; + } + }; + + class ValueSelect { + public: + using value_type = T; + + const value_type& operator()(const std::pair& key_value) const noexcept { + return key_value.second; + } + + value_type& operator()(std::pair& key_value) noexcept { + return key_value.second; + } + }; + + using ht = detail_ordered_hash::ordered_hash, KeySelect, ValueSelect, + Hash, KeyEqual, Allocator, ValueTypeContainer, IndexType>; + +public: + using key_type = typename ht::key_type; + using mapped_type = T; + using value_type = typename ht::value_type; + using size_type = typename ht::size_type; + using difference_type = typename ht::difference_type; + using hasher = typename ht::hasher; + using key_equal = typename ht::key_equal; + using allocator_type = typename ht::allocator_type; + using reference = typename ht::reference; + using const_reference = typename ht::const_reference; + using pointer = typename ht::pointer; + using const_pointer = typename ht::const_pointer; + using iterator = typename ht::iterator; + using const_iterator = typename ht::const_iterator; + using reverse_iterator = typename ht::reverse_iterator; + using const_reverse_iterator = typename ht::const_reverse_iterator; + + using values_container_type = typename ht::values_container_type; + + + /* + * Constructors + */ + ordered_map(): ordered_map(ht::DEFAULT_INIT_BUCKETS_SIZE) { + } + + explicit ordered_map(size_type bucket_count, + const Hash& hash = Hash(), + const KeyEqual& equal = KeyEqual(), + const Allocator& alloc = Allocator()): + m_ht(bucket_count, hash, equal, alloc, ht::DEFAULT_MAX_LOAD_FACTOR) + { + } + + ordered_map(size_type bucket_count, + const Allocator& alloc): ordered_map(bucket_count, Hash(), KeyEqual(), alloc) + { + } + + ordered_map(size_type bucket_count, + const Hash& hash, + const Allocator& alloc): ordered_map(bucket_count, hash, KeyEqual(), alloc) + { + } + + explicit ordered_map(const Allocator& alloc): ordered_map(ht::DEFAULT_INIT_BUCKETS_SIZE, alloc) { + } + + template + ordered_map(InputIt first, InputIt last, + size_type bucket_count = ht::DEFAULT_INIT_BUCKETS_SIZE, + const Hash& hash = Hash(), + const KeyEqual& equal = KeyEqual(), + const Allocator& alloc = Allocator()): ordered_map(bucket_count, hash, equal, alloc) + { + insert(first, last); + } + + template + ordered_map(InputIt first, InputIt last, + size_type bucket_count, + const Allocator& alloc): ordered_map(first, last, bucket_count, Hash(), KeyEqual(), alloc) + { + } + + template + ordered_map(InputIt first, InputIt last, + size_type bucket_count, + const Hash& hash, + const Allocator& alloc): ordered_map(first, last, bucket_count, hash, KeyEqual(), alloc) + { + } + + ordered_map(std::initializer_list init, + size_type bucket_count = ht::DEFAULT_INIT_BUCKETS_SIZE, + const Hash& hash = Hash(), + const KeyEqual& equal = KeyEqual(), + const Allocator& alloc = Allocator()): + ordered_map(init.begin(), init.end(), bucket_count, hash, equal, alloc) + { + } + + ordered_map(std::initializer_list init, + size_type bucket_count, + const Allocator& alloc): + ordered_map(init.begin(), init.end(), bucket_count, Hash(), KeyEqual(), alloc) + { + } + + ordered_map(std::initializer_list init, + size_type bucket_count, + const Hash& hash, + const Allocator& alloc): + ordered_map(init.begin(), init.end(), bucket_count, hash, KeyEqual(), alloc) + { + } + + + ordered_map& operator=(std::initializer_list ilist) { + m_ht.clear(); + + m_ht.reserve(ilist.size()); + m_ht.insert(ilist.begin(), ilist.end()); + + return *this; + } + + allocator_type get_allocator() const { return m_ht.get_allocator(); } + + + + /* + * Iterators + */ + iterator begin() noexcept { return m_ht.begin(); } + const_iterator begin() const noexcept { return m_ht.begin(); } + const_iterator cbegin() const noexcept { return m_ht.cbegin(); } + + iterator end() noexcept { return m_ht.end(); } + const_iterator end() const noexcept { return m_ht.end(); } + const_iterator cend() const noexcept { return m_ht.cend(); } + + reverse_iterator rbegin() noexcept { return m_ht.rbegin(); } + const_reverse_iterator rbegin() const noexcept { return m_ht.rbegin(); } + const_reverse_iterator rcbegin() const noexcept { return m_ht.rcbegin(); } + + reverse_iterator rend() noexcept { return m_ht.rend(); } + const_reverse_iterator rend() const noexcept { return m_ht.rend(); } + const_reverse_iterator rcend() const noexcept { return m_ht.rcend(); } + + + /* + * Capacity + */ + bool empty() const noexcept { return m_ht.empty(); } + size_type size() const noexcept { return m_ht.size(); } + size_type max_size() const noexcept { return m_ht.max_size(); } + + /* + * Modifiers + */ + void clear() noexcept { m_ht.clear(); } + + + + std::pair insert(const value_type& value) { return m_ht.insert(value); } + + template::value>::type* = nullptr> + std::pair insert(P&& value) { return m_ht.emplace(std::forward

(value)); } + + std::pair insert(value_type&& value) { return m_ht.insert(std::move(value)); } + + + iterator insert(const_iterator hint, const value_type& value) { + return m_ht.insert_hint(hint, value); + } + + template::value>::type* = nullptr> + iterator insert(const_iterator hint, P&& value) { + return m_ht.emplace_hint(hint, std::forward

(value)); + } + + iterator insert(const_iterator hint, value_type&& value) { + return m_ht.insert_hint(hint, std::move(value)); + } + + + template + void insert(InputIt first, InputIt last) { m_ht.insert(first, last); } + void insert(std::initializer_list ilist) { m_ht.insert(ilist.begin(), ilist.end()); } + + + + + template + std::pair insert_or_assign(const key_type& k, M&& obj) { + return m_ht.insert_or_assign(k, std::forward(obj)); + } + + template + std::pair insert_or_assign(key_type&& k, M&& obj) { + return m_ht.insert_or_assign(std::move(k), std::forward(obj)); + } + + + template + iterator insert_or_assign(const_iterator hint, const key_type& k, M&& obj) { + return m_ht.insert_or_assign(hint, k, std::forward(obj)); + } + + template + iterator insert_or_assign(const_iterator hint, key_type&& k, M&& obj) { + return m_ht.insert_or_assign(hint, std::move(k), std::forward(obj)); + } + + /** + * Due to the way elements are stored, emplace will need to move or copy the key-value once. + * The method is equivalent to insert(value_type(std::forward(args)...)); + * + * Mainly here for compatibility with the std::unordered_map interface. + */ + template + std::pair emplace(Args&&... args) { return m_ht.emplace(std::forward(args)...); } + + /** + * Due to the way elements are stored, emplace_hint will need to move or copy the key-value once. + * The method is equivalent to insert(hint, value_type(std::forward(args)...)); + * + * Mainly here for compatibility with the std::unordered_map interface. + */ + template + iterator emplace_hint(const_iterator hint, Args&&... args) { + return m_ht.emplace_hint(hint, std::forward(args)...); + } + + + + + template + std::pair try_emplace(const key_type& k, Args&&... args) { + return m_ht.try_emplace(k, std::forward(args)...); + } + + template + std::pair try_emplace(key_type&& k, Args&&... args) { + return m_ht.try_emplace(std::move(k), std::forward(args)...); + } + + template + iterator try_emplace(const_iterator hint, const key_type& k, Args&&... args) { + return m_ht.try_emplace_hint(hint, k, std::forward(args)...); + } + + template + iterator try_emplace(const_iterator hint, key_type&& k, Args&&... args) { + return m_ht.try_emplace_hint(hint, std::move(k), std::forward(args)...); + } + + + + + /** + * When erasing an element, the insert order will be preserved and no holes will be present in the container + * returned by 'values_container()'. + * + * The method is in O(n), if the order is not important 'unordered_erase(...)' method is faster with an O(1) + * average complexity. + */ + iterator erase(iterator pos) { return m_ht.erase(pos); } + + /** + * @copydoc erase(iterator pos) + */ + iterator erase(const_iterator pos) { return m_ht.erase(pos); } + + /** + * @copydoc erase(iterator pos) + */ + iterator erase(const_iterator first, const_iterator last) { return m_ht.erase(first, last); } + + /** + * @copydoc erase(iterator pos) + */ + size_type erase(const key_type& key) { return m_ht.erase(key); } + + /** + * @copydoc erase(iterator pos) + * + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup to the value if you already have the hash. + */ + size_type erase(const key_type& key, std::size_t precalculated_hash) { + return m_ht.erase(key, precalculated_hash); + } + + /** + * @copydoc erase(iterator pos) + * + * This overload only participates in the overload resolution if the typedef KeyEqual::is_transparent exists. + * If so, K must be hashable and comparable to Key. + */ + template::value>::type* = nullptr> + size_type erase(const K& key) { return m_ht.erase(key); } + + /** + * @copydoc erase(const key_type& key, std::size_t precalculated_hash) + * + * This overload only participates in the overload resolution if the typedef KeyEqual::is_transparent exists. + * If so, K must be hashable and comparable to Key. + */ + template::value>::type* = nullptr> + size_type erase(const K& key, std::size_t precalculated_hash) { + return m_ht.erase(key, precalculated_hash); + } + + + + void swap(ordered_map& other) { other.m_ht.swap(m_ht); } + + /* + * Lookup + */ + T& at(const Key& key) { return m_ht.at(key); } + + /** + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + T& at(const Key& key, std::size_t precalculated_hash) { return m_ht.at(key, precalculated_hash); } + + + const T& at(const Key& key) const { return m_ht.at(key); } + + /** + * @copydoc at(const Key& key, std::size_t precalculated_hash) + */ + const T& at(const Key& key, std::size_t precalculated_hash) const { return m_ht.at(key, precalculated_hash); } + + + /** + * This overload only participates in the overload resolution if the typedef KeyEqual::is_transparent exists. + * If so, K must be hashable and comparable to Key. + */ + template::value>::type* = nullptr> + T& at(const K& key) { return m_ht.at(key); } + + /** + * @copydoc at(const K& key) + * + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + template::value>::type* = nullptr> + T& at(const K& key, std::size_t precalculated_hash) { return m_ht.at(key, precalculated_hash); } + + /** + * @copydoc at(const K& key) + */ + template::value>::type* = nullptr> + const T& at(const K& key) const { return m_ht.at(key); } + + /** + * @copydoc at(const K& key, std::size_t precalculated_hash) + */ + template::value>::type* = nullptr> + const T& at(const K& key, std::size_t precalculated_hash) const { return m_ht.at(key, precalculated_hash); } + + + + T& operator[](const Key& key) { return m_ht[key]; } + T& operator[](Key&& key) { return m_ht[std::move(key)]; } + + + + size_type count(const Key& key) const { return m_ht.count(key); } + + /** + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + size_type count(const Key& key, std::size_t precalculated_hash) const { + return m_ht.count(key, precalculated_hash); + } + + /** + * This overload only participates in the overload resolution if the typedef KeyEqual::is_transparent exists. + * If so, K must be hashable and comparable to Key. + */ + template::value>::type* = nullptr> + size_type count(const K& key) const { return m_ht.count(key); } + + /** + * @copydoc count(const K& key) const + * + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + template::value>::type* = nullptr> + size_type count(const K& key, std::size_t precalculated_hash) const { + return m_ht.count(key, precalculated_hash); + } + + + + iterator find(const Key& key) { return m_ht.find(key); } + + /** + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + iterator find(const Key& key, std::size_t precalculated_hash) { return m_ht.find(key, precalculated_hash); } + + const_iterator find(const Key& key) const { return m_ht.find(key); } + + /** + * @copydoc find(const Key& key, std::size_t precalculated_hash) + */ + const_iterator find(const Key& key, std::size_t precalculated_hash) const { + return m_ht.find(key, precalculated_hash); + } + + /** + * This overload only participates in the overload resolution if the typedef KeyEqual::is_transparent exists. + * If so, K must be hashable and comparable to Key. + */ + template::value>::type* = nullptr> + iterator find(const K& key) { return m_ht.find(key); } + + /** + * @copydoc find(const K& key) + * + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + template::value>::type* = nullptr> + iterator find(const K& key, std::size_t precalculated_hash) { return m_ht.find(key, precalculated_hash); } + + /** + * @copydoc find(const K& key) + */ + template::value>::type* = nullptr> + const_iterator find(const K& key) const { return m_ht.find(key); } + + /** + * @copydoc find(const K& key) + * + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + template::value>::type* = nullptr> + const_iterator find(const K& key, std::size_t precalculated_hash) const { + return m_ht.find(key, precalculated_hash); + } + + + + std::pair equal_range(const Key& key) { return m_ht.equal_range(key); } + + /** + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + std::pair equal_range(const Key& key, std::size_t precalculated_hash) { + return m_ht.equal_range(key, precalculated_hash); + } + + std::pair equal_range(const Key& key) const { return m_ht.equal_range(key); } + + /** + * @copydoc equal_range(const Key& key, std::size_t precalculated_hash) + */ + std::pair equal_range(const Key& key, std::size_t precalculated_hash) const { + return m_ht.equal_range(key, precalculated_hash); + } + + /** + * This overload only participates in the overload resolution if the typedef KeyEqual::is_transparent exists. + * If so, K must be hashable and comparable to Key. + */ + template::value>::type* = nullptr> + std::pair equal_range(const K& key) { return m_ht.equal_range(key); } + + /** + * @copydoc equal_range(const K& key) + * + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + template::value>::type* = nullptr> + std::pair equal_range(const K& key, std::size_t precalculated_hash) { + return m_ht.equal_range(key, precalculated_hash); + } + + /** + * @copydoc equal_range(const K& key) + */ + template::value>::type* = nullptr> + std::pair equal_range(const K& key) const { return m_ht.equal_range(key); } + + /** + * @copydoc equal_range(const K& key, std::size_t precalculated_hash) + */ + template::value>::type* = nullptr> + std::pair equal_range(const K& key, std::size_t precalculated_hash) const { + return m_ht.equal_range(key, precalculated_hash); + } + + + + /* + * Bucket interface + */ + size_type bucket_count() const { return m_ht.bucket_count(); } + size_type max_bucket_count() const { return m_ht.max_bucket_count(); } + + + /* + * Hash policy + */ + float load_factor() const { return m_ht.load_factor(); } + float max_load_factor() const { return m_ht.max_load_factor(); } + void max_load_factor(float ml) { m_ht.max_load_factor(ml); } + + void rehash(size_type count) { m_ht.rehash(count); } + void reserve(size_type count) { m_ht.reserve(count); } + + + /* + * Observers + */ + hasher hash_function() const { return m_ht.hash_function(); } + key_equal key_eq() const { return m_ht.key_eq(); } + + + + /* + * Other + */ + + /** + * Convert a const_iterator to an iterator. + */ + iterator mutable_iterator(const_iterator pos) { + return m_ht.mutable_iterator(pos); + } + + /** + * Requires index <= size(). + * + * Return an iterator to the element at index. Return end() if index == size(). + */ + iterator nth(size_type index) { return m_ht.nth(index); } + + /** + * @copydoc nth(size_type index) + */ + const_iterator nth(size_type index) const { return m_ht.nth(index); } + + + /** + * Return const_reference to the first element. Requires the container to not be empty. + */ + const_reference front() const { return m_ht.front(); } + + /** + * Return const_reference to the last element. Requires the container to not be empty. + */ + const_reference back() const { return m_ht.back(); } + + + /** + * Only available if ValueTypeContainer is a std::vector. Same as calling 'values_container().data()'. + */ + template::value>::type* = nullptr> + const typename values_container_type::value_type* data() const noexcept { return m_ht.data(); } + + /** + * Return the container in which the values are stored. The values are in the same order as the insertion order + * and are contiguous in the structure, no holes (size() == values_container().size()). + */ + const values_container_type& values_container() const noexcept { return m_ht.values_container(); } + + template::value>::type* = nullptr> + size_type capacity() const noexcept { return m_ht.capacity(); } + + void shrink_to_fit() { m_ht.shrink_to_fit(); } + + + + /** + * Insert the value before pos shifting all the elements on the right of pos (including pos) one position + * to the right. + * + * Amortized linear time-complexity in the distance between pos and end(). + */ + std::pair insert_at_position(const_iterator pos, const value_type& value) { + return m_ht.insert_at_position(pos, value); + } + + /** + * @copydoc insert_at_position(const_iterator pos, const value_type& value) + */ + std::pair insert_at_position(const_iterator pos, value_type&& value) { + return m_ht.insert_at_position(pos, std::move(value)); + } + + /** + * @copydoc insert_at_position(const_iterator pos, const value_type& value) + * + * Same as insert_at_position(pos, value_type(std::forward(args)...), mainly + * here for coherence. + */ + template + std::pair emplace_at_position(const_iterator pos, Args&&... args) { + return m_ht.emplace_at_position(pos, std::forward(args)...); + } + + /** + * @copydoc insert_at_position(const_iterator pos, const value_type& value) + */ + template + std::pair try_emplace_at_position(const_iterator pos, const key_type& k, Args&&... args) { + return m_ht.try_emplace_at_position(pos, k, std::forward(args)...); + } + + /** + * @copydoc insert_at_position(const_iterator pos, const value_type& value) + */ + template + std::pair try_emplace_at_position(const_iterator pos, key_type&& k, Args&&... args) { + return m_ht.try_emplace_at_position(pos, std::move(k), std::forward(args)...); + } + + + + void pop_back() { m_ht.pop_back(); } + + /** + * Faster erase operation with an O(1) average complexity but it doesn't preserve the insertion order. + * + * If an erasure occurs, the last element of the map will take the place of the erased element. + */ + iterator unordered_erase(iterator pos) { return m_ht.unordered_erase(pos); } + + /** + * @copydoc unordered_erase(iterator pos) + */ + iterator unordered_erase(const_iterator pos) { return m_ht.unordered_erase(pos); } + + /** + * @copydoc unordered_erase(iterator pos) + */ + size_type unordered_erase(const key_type& key) { return m_ht.unordered_erase(key); } + + /** + * @copydoc unordered_erase(iterator pos) + * + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + size_type unordered_erase(const key_type& key, std::size_t precalculated_hash) { + return m_ht.unordered_erase(key, precalculated_hash); + } + + /** + * @copydoc unordered_erase(iterator pos) + * + * This overload only participates in the overload resolution if the typedef KeyEqual::is_transparent exists. + * If so, K must be hashable and comparable to Key. + */ + template::value>::type* = nullptr> + size_type unordered_erase(const K& key) { return m_ht.unordered_erase(key); } + + /** + * @copydoc unordered_erase(const K& key) + * + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + template::value>::type* = nullptr> + size_type unordered_erase(const K& key, std::size_t precalculated_hash) { + return m_ht.unordered_erase(key, precalculated_hash); + } + + /** + * Serialize the map through the `serializer` parameter. + * + * The `serializer` parameter must be a function object that supports the following call: + * - `template void operator()(const U& value);` where the types `std::uint64_t`, `float` and `std::pair` must be supported for U. + * + * The implementation leaves binary compatibility (endianness, IEEE 754 for floats, ...) of the types it serializes + * in the hands of the `Serializer` function object if compatibility is required. + */ + template + void serialize(Serializer& serializer) const { + m_ht.serialize(serializer); + } + + /** + * Deserialize a previously serialized map through the `deserializer` parameter. + * + * The `deserializer` parameter must be a function object that supports the following calls: + * - `template U operator()();` where the types `std::uint64_t`, `float` and `std::pair` must be supported for U. + * + * If the deserialized hash map type is hash compatible with the serialized map, the deserialization process can be + * sped up by setting `hash_compatible` to true. To be hash compatible, the Hash and KeyEqual must behave the same way + * than the ones used on the serialized map. The `std::size_t` must also be of the same size as the one on the platform used + * to serialize the map, the same apply for `IndexType`. If these criteria are not met, the behaviour is undefined with + * `hash_compatible` sets to true. + * + * The behaviour is undefined if the type `Key` and `T` of the `ordered_map` are not the same as the + * types used during serialization. + * + * The implementation leaves binary compatibility (endianness, IEEE 754 for floats, size of int, ...) of the types it + * deserializes in the hands of the `Deserializer` function object if compatibility is required. + */ + template + static ordered_map deserialize(Deserializer& deserializer, bool hash_compatible = false) { + ordered_map map(0); + map.m_ht.deserialize(deserializer, hash_compatible); + + return map; + } + + + + friend bool operator==(const ordered_map& lhs, const ordered_map& rhs) { return lhs.m_ht == rhs.m_ht; } + friend bool operator!=(const ordered_map& lhs, const ordered_map& rhs) { return lhs.m_ht != rhs.m_ht; } + friend bool operator<(const ordered_map& lhs, const ordered_map& rhs) { return lhs.m_ht < rhs.m_ht; } + friend bool operator<=(const ordered_map& lhs, const ordered_map& rhs) { return lhs.m_ht <= rhs.m_ht; } + friend bool operator>(const ordered_map& lhs, const ordered_map& rhs) { return lhs.m_ht > rhs.m_ht; } + friend bool operator>=(const ordered_map& lhs, const ordered_map& rhs) { return lhs.m_ht >= rhs.m_ht; } + + friend void swap(ordered_map& lhs, ordered_map& rhs) { lhs.swap(rhs); } + +private: + ht m_ht; +}; + +} // end namespace tsl + +#endif diff --git a/external/tsl/ordered_set.h b/external/tsl/ordered_set.h new file mode 100644 index 000000000..0932f54b9 --- /dev/null +++ b/external/tsl/ordered_set.h @@ -0,0 +1,688 @@ +/** + * MIT License + * + * Copyright (c) 2017 Thibaut Goetghebuer-Planchon + * + * 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 TSL_ORDERED_SET_H +#define TSL_ORDERED_SET_H + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "ordered_hash.h" + + +namespace tsl { + + +/** + * Implementation of an hash set using open addressing with robin hood with backshift delete to resolve collisions. + * + * The particularity of this hash set is that it remembers the order in which the elements were added and + * provide a way to access the structure which stores these values through the 'values_container()' method. + * The used container is defined by ValueTypeContainer, by default a std::deque is used (grows faster) but + * a std::vector may be used. In this case the set provides a 'data()' method which give a direct access + * to the memory used to store the values (which can be useful to communicate with C API's). + * + * The Key must be copy constructible and/or move constructible. To use `unordered_erase` it also must be swappable. + * + * The behaviour of the hash set is undefined if the destructor of Key throws an exception. + * + * By default the maximum size of a set is limited to 2^32 - 1 values, if needed this can be changed through + * the IndexType template parameter. Using an `uint64_t` will raise this limit to 2^64 - 1 values but each + * bucket will use 16 bytes instead of 8 bytes in addition to the space needed to store the values. + * + * Iterators invalidation: + * - clear, operator=, reserve, rehash: always invalidate the iterators (also invalidate end()). + * - insert, emplace, emplace_hint, operator[]: when a std::vector is used as ValueTypeContainer + * and if size() < capacity(), only end(). + * Otherwise all the iterators are invalidated if an insert occurs. + * - erase, unordered_erase: when a std::vector is used as ValueTypeContainer invalidate the iterator of + * the erased element and all the ones after the erased element (including end()). + * Otherwise all the iterators are invalidated if an erase occurs. + */ +template, + class KeyEqual = std::equal_to, + class Allocator = std::allocator, + class ValueTypeContainer = std::deque, + class IndexType = std::uint_least32_t> +class ordered_set { +private: + template + using has_is_transparent = tsl::detail_ordered_hash::has_is_transparent; + + class KeySelect { + public: + using key_type = Key; + + const key_type& operator()(const Key& key) const noexcept { + return key; + } + + key_type& operator()(Key& key) noexcept { + return key; + } + }; + + using ht = detail_ordered_hash::ordered_hash; + +public: + using key_type = typename ht::key_type; + using value_type = typename ht::value_type; + using size_type = typename ht::size_type; + using difference_type = typename ht::difference_type; + using hasher = typename ht::hasher; + using key_equal = typename ht::key_equal; + using allocator_type = typename ht::allocator_type; + using reference = typename ht::reference; + using const_reference = typename ht::const_reference; + using pointer = typename ht::pointer; + using const_pointer = typename ht::const_pointer; + using iterator = typename ht::iterator; + using const_iterator = typename ht::const_iterator; + using reverse_iterator = typename ht::reverse_iterator; + using const_reverse_iterator = typename ht::const_reverse_iterator; + + using values_container_type = typename ht::values_container_type; + + + /* + * Constructors + */ + ordered_set(): ordered_set(ht::DEFAULT_INIT_BUCKETS_SIZE) { + } + + explicit ordered_set(size_type bucket_count, + const Hash& hash = Hash(), + const KeyEqual& equal = KeyEqual(), + const Allocator& alloc = Allocator()): + m_ht(bucket_count, hash, equal, alloc, ht::DEFAULT_MAX_LOAD_FACTOR) + { + } + + ordered_set(size_type bucket_count, + const Allocator& alloc): ordered_set(bucket_count, Hash(), KeyEqual(), alloc) + { + } + + ordered_set(size_type bucket_count, + const Hash& hash, + const Allocator& alloc): ordered_set(bucket_count, hash, KeyEqual(), alloc) + { + } + + explicit ordered_set(const Allocator& alloc): ordered_set(ht::DEFAULT_INIT_BUCKETS_SIZE, alloc) { + } + + template + ordered_set(InputIt first, InputIt last, + size_type bucket_count = ht::DEFAULT_INIT_BUCKETS_SIZE, + const Hash& hash = Hash(), + const KeyEqual& equal = KeyEqual(), + const Allocator& alloc = Allocator()): ordered_set(bucket_count, hash, equal, alloc) + { + insert(first, last); + } + + template + ordered_set(InputIt first, InputIt last, + size_type bucket_count, + const Allocator& alloc): ordered_set(first, last, bucket_count, Hash(), KeyEqual(), alloc) + { + } + + template + ordered_set(InputIt first, InputIt last, + size_type bucket_count, + const Hash& hash, + const Allocator& alloc): ordered_set(first, last, bucket_count, hash, KeyEqual(), alloc) + { + } + + ordered_set(std::initializer_list init, + size_type bucket_count = ht::DEFAULT_INIT_BUCKETS_SIZE, + const Hash& hash = Hash(), + const KeyEqual& equal = KeyEqual(), + const Allocator& alloc = Allocator()): + ordered_set(init.begin(), init.end(), bucket_count, hash, equal, alloc) + { + } + + ordered_set(std::initializer_list init, + size_type bucket_count, + const Allocator& alloc): + ordered_set(init.begin(), init.end(), bucket_count, Hash(), KeyEqual(), alloc) + { + } + + ordered_set(std::initializer_list init, + size_type bucket_count, + const Hash& hash, + const Allocator& alloc): + ordered_set(init.begin(), init.end(), bucket_count, hash, KeyEqual(), alloc) + { + } + + + ordered_set& operator=(std::initializer_list ilist) { + m_ht.clear(); + + m_ht.reserve(ilist.size()); + m_ht.insert(ilist.begin(), ilist.end()); + + return *this; + } + + allocator_type get_allocator() const { return m_ht.get_allocator(); } + + + /* + * Iterators + */ + iterator begin() noexcept { return m_ht.begin(); } + const_iterator begin() const noexcept { return m_ht.begin(); } + const_iterator cbegin() const noexcept { return m_ht.cbegin(); } + + iterator end() noexcept { return m_ht.end(); } + const_iterator end() const noexcept { return m_ht.end(); } + const_iterator cend() const noexcept { return m_ht.cend(); } + + reverse_iterator rbegin() noexcept { return m_ht.rbegin(); } + const_reverse_iterator rbegin() const noexcept { return m_ht.rbegin(); } + const_reverse_iterator rcbegin() const noexcept { return m_ht.rcbegin(); } + + reverse_iterator rend() noexcept { return m_ht.rend(); } + const_reverse_iterator rend() const noexcept { return m_ht.rend(); } + const_reverse_iterator rcend() const noexcept { return m_ht.rcend(); } + + + /* + * Capacity + */ + bool empty() const noexcept { return m_ht.empty(); } + size_type size() const noexcept { return m_ht.size(); } + size_type max_size() const noexcept { return m_ht.max_size(); } + + /* + * Modifiers + */ + void clear() noexcept { m_ht.clear(); } + + + + std::pair insert(const value_type& value) { return m_ht.insert(value); } + std::pair insert(value_type&& value) { return m_ht.insert(std::move(value)); } + + iterator insert(const_iterator hint, const value_type& value) { + return m_ht.insert_hint(hint, value); + } + + iterator insert(const_iterator hint, value_type&& value) { + return m_ht.insert_hint(hint, std::move(value)); + } + + template + void insert(InputIt first, InputIt last) { m_ht.insert(first, last); } + void insert(std::initializer_list ilist) { m_ht.insert(ilist.begin(), ilist.end()); } + + + + /** + * Due to the way elements are stored, emplace will need to move or copy the key-value once. + * The method is equivalent to insert(value_type(std::forward(args)...)); + * + * Mainly here for compatibility with the std::unordered_map interface. + */ + template + std::pair emplace(Args&&... args) { return m_ht.emplace(std::forward(args)...); } + + /** + * Due to the way elements are stored, emplace_hint will need to move or copy the key-value once. + * The method is equivalent to insert(hint, value_type(std::forward(args)...)); + * + * Mainly here for compatibility with the std::unordered_map interface. + */ + template + iterator emplace_hint(const_iterator hint, Args&&... args) { + return m_ht.emplace_hint(hint, std::forward(args)...); + } + + /** + * When erasing an element, the insert order will be preserved and no holes will be present in the container + * returned by 'values_container()'. + * + * The method is in O(n), if the order is not important 'unordered_erase(...)' method is faster with an O(1) + * average complexity. + */ + iterator erase(iterator pos) { return m_ht.erase(pos); } + + /** + * @copydoc erase(iterator pos) + */ + iterator erase(const_iterator pos) { return m_ht.erase(pos); } + + /** + * @copydoc erase(iterator pos) + */ + iterator erase(const_iterator first, const_iterator last) { return m_ht.erase(first, last); } + + /** + * @copydoc erase(iterator pos) + */ + size_type erase(const key_type& key) { return m_ht.erase(key); } + + /** + * @copydoc erase(iterator pos) + * + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup to the value if you already have the hash. + */ + size_type erase(const key_type& key, std::size_t precalculated_hash) { + return m_ht.erase(key, precalculated_hash); + } + + /** + * @copydoc erase(iterator pos) + * + * This overload only participates in the overload resolution if the typedef KeyEqual::is_transparent exists. + * If so, K must be hashable and comparable to Key. + */ + template::value>::type* = nullptr> + size_type erase(const K& key) { return m_ht.erase(key); } + + /** + * @copydoc erase(const key_type& key, std::size_t precalculated_hash) + * + * This overload only participates in the overload resolution if the typedef KeyEqual::is_transparent exists. + * If so, K must be hashable and comparable to Key. + */ + template::value>::type* = nullptr> + size_type erase(const K& key, std::size_t precalculated_hash) { + return m_ht.erase(key, precalculated_hash); + } + + + + void swap(ordered_set& other) { other.m_ht.swap(m_ht); } + + /* + * Lookup + */ + size_type count(const Key& key) const { return m_ht.count(key); } + + /** + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + size_type count(const Key& key, std::size_t precalculated_hash) const { + return m_ht.count(key, precalculated_hash); + } + + /** + * This overload only participates in the overload resolution if the typedef KeyEqual::is_transparent exists. + * If so, K must be hashable and comparable to Key. + */ + template::value>::type* = nullptr> + size_type count(const K& key) const { return m_ht.count(key); } + + /** + * @copydoc count(const K& key) const + * + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + template::value>::type* = nullptr> + size_type count(const K& key, std::size_t precalculated_hash) const { + return m_ht.count(key, precalculated_hash); + } + + + + + iterator find(const Key& key) { return m_ht.find(key); } + + /** + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + iterator find(const Key& key, std::size_t precalculated_hash) { return m_ht.find(key, precalculated_hash); } + + const_iterator find(const Key& key) const { return m_ht.find(key); } + + /** + * @copydoc find(const Key& key, std::size_t precalculated_hash) + */ + const_iterator find(const Key& key, std::size_t precalculated_hash) const { + return m_ht.find(key, precalculated_hash); + } + + /** + * This overload only participates in the overload resolution if the typedef KeyEqual::is_transparent exists. + * If so, K must be hashable and comparable to Key. + */ + template::value>::type* = nullptr> + iterator find(const K& key) { return m_ht.find(key); } + + /** + * @copydoc find(const K& key) + * + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + template::value>::type* = nullptr> + iterator find(const K& key, std::size_t precalculated_hash) { return m_ht.find(key, precalculated_hash); } + + /** + * @copydoc find(const K& key) + */ + template::value>::type* = nullptr> + const_iterator find(const K& key) const { return m_ht.find(key); } + + /** + * @copydoc find(const K& key) + * + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + template::value>::type* = nullptr> + const_iterator find(const K& key, std::size_t precalculated_hash) const { + return m_ht.find(key, precalculated_hash); + } + + + + std::pair equal_range(const Key& key) { return m_ht.equal_range(key); } + + /** + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + std::pair equal_range(const Key& key, std::size_t precalculated_hash) { + return m_ht.equal_range(key, precalculated_hash); + } + + std::pair equal_range(const Key& key) const { return m_ht.equal_range(key); } + + /** + * @copydoc equal_range(const Key& key, std::size_t precalculated_hash) + */ + std::pair equal_range(const Key& key, std::size_t precalculated_hash) const { + return m_ht.equal_range(key, precalculated_hash); + } + + /** + * This overload only participates in the overload resolution if the typedef KeyEqual::is_transparent exists. + * If so, K must be hashable and comparable to Key. + */ + template::value>::type* = nullptr> + std::pair equal_range(const K& key) { return m_ht.equal_range(key); } + + /** + * @copydoc equal_range(const K& key) + * + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + template::value>::type* = nullptr> + std::pair equal_range(const K& key, std::size_t precalculated_hash) { + return m_ht.equal_range(key, precalculated_hash); + } + + /** + * @copydoc equal_range(const K& key) + */ + template::value>::type* = nullptr> + std::pair equal_range(const K& key) const { return m_ht.equal_range(key); } + + /** + * @copydoc equal_range(const K& key, std::size_t precalculated_hash) + */ + template::value>::type* = nullptr> + std::pair equal_range(const K& key, std::size_t precalculated_hash) const { + return m_ht.equal_range(key, precalculated_hash); + } + + + /* + * Bucket interface + */ + size_type bucket_count() const { return m_ht.bucket_count(); } + size_type max_bucket_count() const { return m_ht.max_bucket_count(); } + + + /* + * Hash policy + */ + float load_factor() const { return m_ht.load_factor(); } + float max_load_factor() const { return m_ht.max_load_factor(); } + void max_load_factor(float ml) { m_ht.max_load_factor(ml); } + + void rehash(size_type count) { m_ht.rehash(count); } + void reserve(size_type count) { m_ht.reserve(count); } + + + /* + * Observers + */ + hasher hash_function() const { return m_ht.hash_function(); } + key_equal key_eq() const { return m_ht.key_eq(); } + + + /* + * Other + */ + + /** + * Convert a const_iterator to an iterator. + */ + iterator mutable_iterator(const_iterator pos) { + return m_ht.mutable_iterator(pos); + } + + /** + * Requires index <= size(). + * + * Return an iterator to the element at index. Return end() if index == size(). + */ + iterator nth(size_type index) { return m_ht.nth(index); } + + /** + * @copydoc nth(size_type index) + */ + const_iterator nth(size_type index) const { return m_ht.nth(index); } + + + /** + * Return const_reference to the first element. Requires the container to not be empty. + */ + const_reference front() const { return m_ht.front(); } + + /** + * Return const_reference to the last element. Requires the container to not be empty. + */ + const_reference back() const { return m_ht.back(); } + + + /** + * Only available if ValueTypeContainer is a std::vector. Same as calling 'values_container().data()'. + */ + template::value>::type* = nullptr> + const typename values_container_type::value_type* data() const noexcept { return m_ht.data(); } + + /** + * Return the container in which the values are stored. The values are in the same order as the insertion order + * and are contiguous in the structure, no holes (size() == values_container().size()). + */ + const values_container_type& values_container() const noexcept { return m_ht.values_container(); } + + template::value>::type* = nullptr> + size_type capacity() const noexcept { return m_ht.capacity(); } + + void shrink_to_fit() { m_ht.shrink_to_fit(); } + + + + /** + * Insert the value before pos shifting all the elements on the right of pos (including pos) one position + * to the right. + * + * Amortized linear time-complexity in the distance between pos and end(). + */ + std::pair insert_at_position(const_iterator pos, const value_type& value) { + return m_ht.insert_at_position(pos, value); + } + + /** + * @copydoc insert_at_position(const_iterator pos, const value_type& value) + */ + std::pair insert_at_position(const_iterator pos, value_type&& value) { + return m_ht.insert_at_position(pos, std::move(value)); + } + + /** + * @copydoc insert_at_position(const_iterator pos, const value_type& value) + * + * Same as insert_at_position(pos, value_type(std::forward(args)...), mainly + * here for coherence. + */ + template + std::pair emplace_at_position(const_iterator pos, Args&&... args) { + return m_ht.emplace_at_position(pos, std::forward(args)...); + } + + + + void pop_back() { m_ht.pop_back(); } + + /** + * Faster erase operation with an O(1) average complexity but it doesn't preserve the insertion order. + * + * If an erasure occurs, the last element of the map will take the place of the erased element. + */ + iterator unordered_erase(iterator pos) { return m_ht.unordered_erase(pos); } + + /** + * @copydoc unordered_erase(iterator pos) + */ + iterator unordered_erase(const_iterator pos) { return m_ht.unordered_erase(pos); } + + /** + * @copydoc unordered_erase(iterator pos) + */ + size_type unordered_erase(const key_type& key) { return m_ht.unordered_erase(key); } + + /** + * @copydoc unordered_erase(iterator pos) + * + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + size_type unordered_erase(const key_type& key, std::size_t precalculated_hash) { + return m_ht.unordered_erase(key, precalculated_hash); + } + + /** + * @copydoc unordered_erase(iterator pos) + * + * This overload only participates in the overload resolution if the typedef KeyEqual::is_transparent exists. + * If so, K must be hashable and comparable to Key. + */ + template::value>::type* = nullptr> + size_type unordered_erase(const K& key) { return m_ht.unordered_erase(key); } + + /** + * @copydoc unordered_erase(const K& key) + * + * Use the hash value 'precalculated_hash' instead of hashing the key. The hash value should be the same + * as hash_function()(key). Useful to speed-up the lookup if you already have the hash. + */ + template::value>::type* = nullptr> + size_type unordered_erase(const K& key, std::size_t precalculated_hash) { + return m_ht.unordered_erase(key, precalculated_hash); + } + + /** + * Serialize the set through the `serializer` parameter. + * + * The `serializer` parameter must be a function object that supports the following call: + * - `void operator()(const U& value);` where the types `std::uint64_t`, `float` and `Key` must be supported for U. + * + * The implementation leaves binary compatibility (endianness, IEEE 754 for floats, ...) of the types it serializes + * in the hands of the `Serializer` function object if compatibility is required. + */ + template + void serialize(Serializer& serializer) const { + m_ht.serialize(serializer); + } + + /** + * Deserialize a previously serialized set through the `deserializer` parameter. + * + * The `deserializer` parameter must be a function object that supports the following calls: + * - `template U operator()();` where the types `std::uint64_t`, `float` and `Key` must be supported for U. + * + * If the deserialized hash set type is hash compatible with the serialized set, the deserialization process can be + * sped up by setting `hash_compatible` to true. To be hash compatible, the Hash and KeyEqual must behave the same way + * than the ones used on the serialized map. The `std::size_t` must also be of the same size as the one on the platform used + * to serialize the map, the same apply for `IndexType`. If these criteria are not met, the behaviour is undefined with + * `hash_compatible` sets to true. + * + * The behaviour is undefined if the type `Key` of the `ordered_set` is not the same as the + * type used during serialization. + * + * The implementation leaves binary compatibility (endianness, IEEE 754 for floats, size of int, ...) of the types it + * deserializes in the hands of the `Deserializer` function object if compatibility is required. + */ + template + static ordered_set deserialize(Deserializer& deserializer, bool hash_compatible = false) { + ordered_set set(0); + set.m_ht.deserialize(deserializer, hash_compatible); + + return set; + } + + + + friend bool operator==(const ordered_set& lhs, const ordered_set& rhs) { return lhs.m_ht == rhs.m_ht; } + friend bool operator!=(const ordered_set& lhs, const ordered_set& rhs) { return lhs.m_ht != rhs.m_ht; } + friend bool operator<(const ordered_set& lhs, const ordered_set& rhs) { return lhs.m_ht < rhs.m_ht; } + friend bool operator<=(const ordered_set& lhs, const ordered_set& rhs) { return lhs.m_ht <= rhs.m_ht; } + friend bool operator>(const ordered_set& lhs, const ordered_set& rhs) { return lhs.m_ht > rhs.m_ht; } + friend bool operator>=(const ordered_set& lhs, const ordered_set& rhs) { return lhs.m_ht >= rhs.m_ht; } + + friend void swap(ordered_set& lhs, ordered_set& rhs) { lhs.swap(rhs); } + +private: + ht m_ht; +}; + +} // end namespace tsl + +#endif diff --git a/include/mesh.tcc b/include/mesh.tcc index c603a2fea..26f1b5820 100644 --- a/include/mesh.tcc +++ b/include/mesh.tcc @@ -1207,8 +1207,10 @@ void mpm::Mesh::apply_traction_on_particles(double current_time) { std::placeholders::_1, dir, traction)); } if (!particle_tractions_.empty()) { - this->iterate_over_particles(std::bind( - &mpm::ParticleBase::map_traction_force, std::placeholders::_1)); + this->iterate_over_particles( + [](std::shared_ptr> ptr) { + return mpm::particle::map_traction_force(ptr); + }); } } @@ -1771,8 +1773,8 @@ void mpm::Mesh::inject_particles(double current_time) { } } for (auto particle : injected_particles) { - particle->compute_volume(); - particle->compute_mass(); + mpm::particle::compute_volume(particle); + mpm::particle::compute_mass(particle); } } } diff --git a/include/mpm_properties.h b/include/mpm_properties.h new file mode 100644 index 000000000..170c2bec6 --- /dev/null +++ b/include/mpm_properties.h @@ -0,0 +1,28 @@ +#ifndef MPM_PROPERTIES_H_ +#define MPM_PROPERTIES_H_ + +namespace mpm { +namespace properties { +//! Boolean Properties +enum Boolean : unsigned int { SetTraction, Friction, GenericBC }; +//! Scalar Properties +enum Scalar : unsigned int { + Mass, + Volume, + MassDensity, + MassPressure, + Pressure +}; +//! Vector Properties +enum Vector : unsigned int { + Displacement, + Velocity, + Acceleration, + Momentum, + ExternalForce, + InternalForce +}; +} // namespace properties +} // namespace mpm + +#endif // MPM_PROPERTIES_H_ diff --git a/include/node.h b/include/node.h index 31a6b9a03..993abca0b 100644 --- a/include/node.h +++ b/include/node.h @@ -66,6 +66,48 @@ class Node : public NodeBase { //! Return status bool status() const override { return status_; } + //! Assign boolean property at the nodes + //! \param[in] property Name of the property to assign + //! \param[in] boolean Property boolean (true/false) of the node + void assign_boolean_property(mpm::properties::Boolean property, + bool boolean) noexcept override; + + //! Return boolean property + //! \param[in] property Name of the property to update + //! \retval boolean property at node + bool boolean_property(mpm::properties::Boolean property) const override; + + //! Update scalar property at the nodes + //! \param[in] property Name of the property to update + //! \param[in] update A boolean to update (true) or assign (false) + //! \param[in] phase Index corresponding to the phase + //! \param[in] value Property value from the particles in a cell + void update_scalar_property(mpm::properties::Scalar property, bool update, + unsigned phase, double value) noexcept override; + + //! Return property at a given node for a given phase + //! \param[in] property Name of the property to return + //! \param[in] phase Index corresponding to the phase + //! \retval scalar property at the designated phase + double scalar_property(mpm::properties::Scalar property, + unsigned phase) const override; + + //! Update vector property at the nodes + //! \param[in] property Name of the property to update + //! \param[in] update A boolean to update (true) or assign (false) + //! \param[in] phase Index corresponding to the phase + //! \param[in] value Property value from the particles in a cell + virtual void update_vector_property( + mpm::properties::Vector property, bool update, unsigned phase, + const Eigen::Matrix& value) noexcept override; + + //! Return property at a given node for a given phase + //! \param[in] property Name of the property to return + //! \param[in] phase Index corresponding to the phase + //! \retval vector property at the designated phase + virtual Eigen::Matrix vector_property( + mpm::properties::Vector property, unsigned phase) const override; + //! Update mass at the nodes from particle //! \param[in] update A boolean to update (true) or assign (false) //! \param[in] phase Index corresponding to the phase @@ -74,7 +116,9 @@ class Node : public NodeBase { //! Return mass at a given node for a given phase //! \param[in] phase Index corresponding to the phase - double mass(unsigned phase) const override { return mass_(phase); } + double mass(unsigned phase) const override { + return this->scalar_property(mpm::properties::Scalar::Mass, phase); + } //! Update volume at the nodes from particle //! \param[in] update A boolean to update (true) or assign (false) @@ -85,7 +129,9 @@ class Node : public NodeBase { //! Return volume at a given node for a given phase //! \param[in] phase Index corresponding to the phase - double volume(unsigned phase) const override { return volume_(phase); } + double volume(unsigned phase) const override { + return this->scalar_property(mpm::properties::Scalar::Volume, phase); + } //! Assign concentrated force to the node //! \param[in] phase Index corresponding to the phase @@ -112,7 +158,7 @@ class Node : public NodeBase { //! Return external force at a given node for a given phase //! \param[in] phase Index corresponding to the phase VectorDim external_force(unsigned phase) const override { - return external_force_.col(phase); + return this->vector_property(mpm::properties::Vector::ExternalForce, phase); } //! Update internal force (body force / traction force) @@ -125,24 +171,31 @@ class Node : public NodeBase { //! Return internal force at a given node for a given phase //! \param[in] phase Index corresponding to the phase VectorDim internal_force(unsigned phase) const override { - return internal_force_.col(phase); + return this->vector_property(mpm::properties::Vector::InternalForce, phase); } //! Update pressure at the nodes from particle + //! \param[in] update A boolean to update (true) or assign (false) //! \param[in] phase Index corresponding to the phase //! \param[in] mass_pressure Product of mass x pressure of a particle - void update_mass_pressure(unsigned phase, + void update_mass_pressure(bool update, unsigned phase, double mass_pressure) noexcept override; + //! Compute pressure from the mass pressure + virtual void compute_pressure() override; + //! Assign pressure at the nodes from particle //! \param[in] update A boolean to update (true) or assign (false) //! \param[in] phase Index corresponding to the phase - //! \param[in] mass_pressure Product of mass x pressure of a particle - void assign_pressure(unsigned phase, double mass_pressure) override; + //! \param[in] pressure Pressure of a particle + virtual void update_pressure(bool update, unsigned phase, + double pressure) override; //! Return pressure at a given node for a given phase //! \param[in] phase Index corresponding to the phase - double pressure(unsigned phase) const override { return pressure_(phase); } + double pressure(unsigned phase) const override { + return this->scalar_property(mpm::properties::Scalar::Pressure, phase); + } //! Update momentum at the nodes //! \param[in] update A boolean to update (true) or assign (false) @@ -154,7 +207,7 @@ class Node : public NodeBase { //! Return momentum at a given node for a given phase //! \param[in] phase Index corresponding to the phase VectorDim momentum(unsigned phase) const override { - return momentum_.col(phase); + return this->vector_property(mpm::properties::Vector::Momentum, phase); } //! Compute velocity from the momentum @@ -163,7 +216,7 @@ class Node : public NodeBase { //! Return velocity at a given node for a given phase //! \param[in] phase Index corresponding to the phase VectorDim velocity(unsigned phase) const override { - return velocity_.col(phase); + return this->vector_property(mpm::properties::Vector::Velocity, phase); } //! Update nodal acceleration @@ -176,7 +229,7 @@ class Node : public NodeBase { //! Return acceleration at a given node for a given phase //! \param[in] phase Index corresponding to the phase VectorDim acceleration(unsigned phase) const override { - return acceleration_.col(phase); + return this->vector_property(mpm::properties::Vector::Acceleration, phase); } //! Compute acceleration and velocity @@ -218,7 +271,7 @@ class Node : public NodeBase { void assign_rotation_matrix( const Eigen::Matrix& rotation_matrix) override { rotation_matrix_ = rotation_matrix; - generic_boundary_constraints_ = true; + this->assign_boolean_property(mpm::properties::Boolean::GenericBC, true); } //! Add material id from material points to list of materials in materials_ @@ -279,35 +332,23 @@ class Node : public NodeBase { unsigned dof_{std::numeric_limits::max()}; //! Status bool status_{false}; - //! Mass - Eigen::Matrix mass_; - //! Volume - Eigen::Matrix volume_; - //! External force - Eigen::Matrix external_force_; - //! Internal force - Eigen::Matrix internal_force_; - //! Pressure - Eigen::Matrix pressure_; + //! Boolean properties + fc::vector_map boolean_properties_; + //! Scalar properties + fc::vector_map> + scalar_properties_; + //! Vector properties + fc::vector_map> + vector_properties_; //! Displacement Eigen::Matrix contact_displacement_; - //! Velocity - Eigen::Matrix velocity_; - //! Momentum - Eigen::Matrix momentum_; - //! Acceleration - Eigen::Matrix acceleration_; //! Velocity constraints std::map velocity_constraints_; //! Rotation matrix for general velocity constraints Eigen::Matrix rotation_matrix_; //! Material ids whose information was passed to this node std::set material_ids_; - //! A general velocity (non-Cartesian/inclined) constraint is specified at the - //! node - bool generic_boundary_constraints_{false}; //! Frictional constraints - bool friction_{false}; std::tuple friction_constraint_; //! Concentrated force Eigen::Matrix concentrated_force_; diff --git a/include/node.tcc b/include/node.tcc index 7d3a9cfcb..6fb6f5bad 100644 --- a/include/node.tcc +++ b/include/node.tcc @@ -17,23 +17,79 @@ mpm::Node::Node( // Clear any velocity constraints velocity_constraints_.clear(); concentrated_force_.setZero(); + + // Initialize boolean properties + // Friction + boolean_properties_.emplace( + std::make_pair(mpm::properties::Boolean::Friction, false)); + // GenericBC + boolean_properties_.emplace( + std::make_pair(mpm::properties::Boolean::GenericBC, false)); + + // Initialize scalar properties + // Mass + scalar_properties_.emplace( + std::make_pair(mpm::properties::Scalar::Mass, + Eigen::Matrix::Zero())); + // Volume + scalar_properties_.emplace( + std::make_pair(mpm::properties::Scalar::Volume, + Eigen::Matrix::Zero())); + // MassPressure + scalar_properties_.emplace( + std::make_pair(mpm::properties::Scalar::MassPressure, + Eigen::Matrix::Zero())); + // Pressure + scalar_properties_.emplace( + std::make_pair(mpm::properties::Scalar::Pressure, + Eigen::Matrix::Zero())); + + // Initialize vector properties + // Velocity + vector_properties_.emplace( + std::make_pair(mpm::properties::Vector::Velocity, + Eigen::Matrix::Zero())); + // Acceleration + vector_properties_.emplace( + std::make_pair(mpm::properties::Vector::Acceleration, + Eigen::Matrix::Zero())); + // Momentum + vector_properties_.emplace( + std::make_pair(mpm::properties::Vector::Momentum, + Eigen::Matrix::Zero())); + // ExternalForce + vector_properties_.emplace( + std::make_pair(mpm::properties::Vector::ExternalForce, + Eigen::Matrix::Zero())); + // InternalForce + vector_properties_.emplace( + std::make_pair(mpm::properties::Vector::InternalForce, + Eigen::Matrix::Zero())); + this->initialise(); } //! Initialise nodal properties template void mpm::Node::initialise() noexcept { - mass_.setZero(); - volume_.setZero(); - external_force_.setZero(); - internal_force_.setZero(); - pressure_.setZero(); - contact_displacement_.setZero(); - velocity_.setZero(); - momentum_.setZero(); - acceleration_.setZero(); status_ = false; + + // Initialise nodal scalar properties + scalar_properties_.at(mpm::properties::Scalar::Mass).setZero(); + scalar_properties_.at(mpm::properties::Scalar::Volume).setZero(); + scalar_properties_.at(mpm::properties::Scalar::MassPressure).setZero(); + scalar_properties_.at(mpm::properties::Scalar::Pressure).setZero(); + + // Initialise nodal vector properties + vector_properties_.at(mpm::properties::Vector::Velocity).setZero(); + vector_properties_.at(mpm::properties::Vector::Acceleration).setZero(); + vector_properties_.at(mpm::properties::Vector::Momentum).setZero(); + vector_properties_.at(mpm::properties::Vector::ExternalForce).setZero(); + vector_properties_.at(mpm::properties::Vector::InternalForce).setZero(); + + // Initialise variables for contact material_ids_.clear(); + contact_displacement_.setZero(); } //! Initialise shared pointer to nodal properties pool @@ -46,32 +102,94 @@ void mpm::Node::initialise_property_handle( this->prop_id_ = prop_id; } -//! Update mass at the nodes from particle +//! Assign boolean property at the nodes template -void mpm::Node::update_mass(bool update, unsigned phase, - double mass) noexcept { +void mpm::Node::assign_boolean_property( + mpm::properties::Boolean property, bool boolean) noexcept { + // Update/assign value + node_mutex_.lock(); + boolean_properties_.at(property) = boolean; + node_mutex_.unlock(); +} + +//! Return boolean property +template +bool mpm::Node::boolean_property( + mpm::properties::Boolean property) const { + return boolean_properties_.at(property); +} + +//! Update scalar property at the nodes from particle +template +void mpm::Node::update_scalar_property( + mpm::properties::Scalar property, bool update, unsigned phase, + double value) noexcept { + // Assert phase + assert(phase < Tnphases); + // Decide to update or assign const double factor = (update == true) ? 1. : 0.; - // Update/assign mass + // Update/assign value node_mutex_.lock(); - mass_(phase) = (mass_(phase) * factor) + mass; + scalar_properties_.at(property)[phase] = + (scalar_properties_.at(property)[phase] * factor) + value; node_mutex_.unlock(); } -//! Update volume at the nodes from particle +//! Update scalar property at the nodes from particle template -void mpm::Node::update_volume(bool update, unsigned phase, - double volume) noexcept { +double mpm::Node::scalar_property( + mpm::properties::Scalar property, unsigned phase) const { + // Assert phase + assert(phase < Tnphases); + return scalar_properties_.at(property)[phase]; +} + +//! Update vector property at the nodes from particle +template +void mpm::Node::update_vector_property( + mpm::properties::Vector property, bool update, unsigned phase, + const Eigen::Matrix& value) noexcept { + // Assert phase + assert(phase < Tnphases); + // Decide to update or assign const double factor = (update == true) ? 1. : 0.; - // Update/assign volume + // Update/assign value node_mutex_.lock(); - volume_(phase) = volume_(phase) * factor + volume; + Eigen::Matrix vecvalue = + vector_properties_.at(property).col(phase); + vector_properties_.at(property).col(phase) = (vecvalue * factor) + value; node_mutex_.unlock(); } +//! Update vector property at the nodes from particle +template +Eigen::Matrix mpm::Node::vector_property( + mpm::properties::Vector property, unsigned phase) const { + // Assert phase + assert(phase < Tnphases); + return vector_properties_.at(property).col(phase); +} + +//! Update mass at the nodes from particle +template +void mpm::Node::update_mass(bool update, unsigned phase, + double mass) noexcept { + this->update_scalar_property(mpm::properties::Scalar::Mass, update, phase, + mass); +} + +//! Update volume at the nodes from particle +template +void mpm::Node::update_volume(bool update, unsigned phase, + double volume) noexcept { + this->update_scalar_property(mpm::properties::Scalar::Volume, update, phase, + volume); +} + // Assign concentrated force to the node template bool mpm::Node::assign_concentrated_force( @@ -110,16 +228,8 @@ template void mpm::Node::update_external_force( bool update, unsigned phase, const Eigen::Matrix& force) noexcept { - // Assert - assert(phase < Tnphases); - - // Decide to update or assign - const double factor = (update == true) ? 1. : 0.; - - // Update/assign external force - node_mutex_.lock(); - external_force_.col(phase) = external_force_.col(phase) * factor + force; - node_mutex_.unlock(); + this->update_vector_property(mpm::properties::Vector::ExternalForce, update, + phase, force); } //! Update internal force (body force / traction force) @@ -127,16 +237,8 @@ template void mpm::Node::update_internal_force( bool update, unsigned phase, const Eigen::Matrix& force) noexcept { - // Assert - assert(phase < Tnphases); - - // Decide to update or assign - const double factor = (update == true) ? 1. : 0.; - - // Update/assign internal force - node_mutex_.lock(); - internal_force_.col(phase) = internal_force_.col(phase) * factor + force; - node_mutex_.unlock(); + this->update_vector_property(mpm::properties::Vector::InternalForce, update, + phase, force); } //! Assign nodal momentum @@ -144,42 +246,43 @@ template void mpm::Node::update_momentum( bool update, unsigned phase, const Eigen::Matrix& momentum) noexcept { - // Assert - assert(phase < Tnphases); - - // Decide to update or assign - const double factor = (update == true) ? 1. : 0.; - - // Update/assign momentum - node_mutex_.lock(); - momentum_.col(phase) = momentum_.col(phase) * factor + momentum; - node_mutex_.unlock(); + this->update_vector_property(mpm::properties::Vector::Momentum, update, phase, + momentum); } //! Update pressure at the nodes from particle template void mpm::Node::update_mass_pressure( - unsigned phase, double mass_pressure) noexcept { - // Assert - assert(phase < Tnphases); + bool update, unsigned phase, double mass_pressure) noexcept { + this->update_scalar_property(mpm::properties::Scalar::MassPressure, update, + phase, mass_pressure); +} +//! Compute pressure from mass pressure +//! pressure = mass pressure / mass +template +void mpm::Node::compute_pressure() { const double tolerance = 1.E-16; - // Compute pressure from mass*pressure - if (mass_(phase) > tolerance) { - node_mutex_.lock(); - pressure_(phase) += mass_pressure / mass_(phase); - node_mutex_.unlock(); + for (unsigned phase = 0; phase < Tnphases; ++phase) { + if (this->mass(phase) > tolerance) { + scalar_properties_.at(mpm::properties::Scalar::Pressure)(phase) = + scalar_properties_.at(mpm::properties::Scalar::MassPressure)(phase) / + this->mass(phase); + + // Check to see if value is below threshold + if (std::abs(this->pressure(phase)) < 1.E-15) + scalar_properties_.at(mpm::properties::Scalar::Pressure)(phase) = 0.; + } } } //! Assign pressure at the nodes from particle template -void mpm::Node::assign_pressure(unsigned phase, +void mpm::Node::update_pressure(bool update, + unsigned phase, double pressure) { - // Compute pressure from mass*pressure - node_mutex_.lock(); - pressure_(phase) = pressure; - node_mutex_.unlock(); + this->update_scalar_property(mpm::properties::Scalar::Pressure, update, phase, + pressure); } //! Compute velocity from momentum @@ -188,13 +291,15 @@ template void mpm::Node::compute_velocity() { const double tolerance = 1.E-16; for (unsigned phase = 0; phase < Tnphases; ++phase) { - if (mass_(phase) > tolerance) { - velocity_.col(phase) = momentum_.col(phase) / mass_(phase); + if (this->mass(phase) > tolerance) { + vector_properties_.at(mpm::properties::Vector::Velocity).col(phase) = + this->momentum(phase) / this->mass(phase); // Check to see if value is below threshold - for (unsigned i = 0; i < velocity_.rows(); ++i) - if (std::abs(velocity_.col(phase)(i)) < 1.E-15) - velocity_.col(phase)(i) = 0.; + for (unsigned i = 0; i < Tdim; ++i) + if (std::abs(this->velocity(phase)(i)) < 1.E-15) + vector_properties_.at(mpm::properties::Vector::Velocity) + .col(phase)(i) = 0.; } } @@ -208,15 +313,8 @@ template void mpm::Node::update_acceleration( bool update, unsigned phase, const Eigen::Matrix& acceleration) noexcept { - assert(phase < Tnphases); - - // Decide to update or assign - const double factor = (update == true) ? 1. : 0.; - - //! Update/assign acceleration - node_mutex_.lock(); - acceleration_.col(phase) = acceleration_.col(phase) * factor + acceleration; - node_mutex_.unlock(); + this->update_vector_property(mpm::properties::Vector::Acceleration, update, + phase, acceleration); } //! Compute acceleration and velocity @@ -225,28 +323,31 @@ bool mpm::Node::compute_acceleration_velocity( unsigned phase, double dt) noexcept { bool status = false; const double tolerance = 1.0E-15; - if (mass_(phase) > tolerance) { + if (this->mass(phase) > tolerance) { // acceleration = (unbalaced force / mass) - this->acceleration_.col(phase) = - (this->external_force_.col(phase) + this->internal_force_.col(phase)) / - this->mass_(phase); + vector_properties_.at(mpm::properties::Vector::Acceleration).col(phase) = + (this->external_force(phase) + this->internal_force(phase)) / + this->mass(phase); // Apply friction constraints this->apply_friction_constraints(dt); // Velocity += acceleration * dt - this->velocity_.col(phase) += this->acceleration_.col(phase) * dt; + vector_properties_.at(mpm::properties::Vector::Velocity).col(phase) += + this->acceleration(phase) * dt; // Apply velocity constraints, which also sets acceleration to 0, // when velocity is set. this->apply_velocity_constraints(); // Set a threshold for (unsigned i = 0; i < Tdim; ++i) - if (std::abs(velocity_.col(phase)(i)) < tolerance) - velocity_.col(phase)(i) = 0.; + if (std::abs(this->velocity(phase)(i)) < tolerance) + vector_properties_.at(mpm::properties::Vector::Velocity).col(phase)(i) = + 0.; for (unsigned i = 0; i < Tdim; ++i) - if (std::abs(acceleration_.col(phase)(i)) < tolerance) - acceleration_.col(phase)(i) = 0.; + if (std::abs(this->acceleration(phase)(i)) < tolerance) + vector_properties_.at(mpm::properties::Vector::Acceleration) + .col(phase)(i) = 0.; status = true; } return status; @@ -258,31 +359,34 @@ bool mpm::Node::compute_acceleration_velocity_cundall( unsigned phase, double dt, double damping_factor) noexcept { bool status = false; const double tolerance = 1.0E-15; - if (mass_(phase) > tolerance) { + if (this->mass(phase) > tolerance) { // acceleration = (unbalaced force / mass) auto unbalanced_force = - this->external_force_.col(phase) + this->internal_force_.col(phase); - this->acceleration_.col(phase) = + this->external_force(phase) + this->internal_force(phase); + vector_properties_.at(mpm::properties::Vector::Acceleration).col(phase) = (unbalanced_force - damping_factor * unbalanced_force.norm() * - this->velocity_.col(phase).cwiseSign()) / - this->mass_(phase); + this->velocity(phase).cwiseSign()) / + this->mass(phase); // Apply friction constraints this->apply_friction_constraints(dt); // Velocity += acceleration * dt - this->velocity_.col(phase) += this->acceleration_.col(phase) * dt; + vector_properties_.at(mpm::properties::Vector::Velocity).col(phase) += + this->acceleration(phase) * dt; // Apply velocity constraints, which also sets acceleration to 0, // when velocity is set. this->apply_velocity_constraints(); // Set a threshold for (unsigned i = 0; i < Tdim; ++i) - if (std::abs(velocity_.col(phase)(i)) < tolerance) - velocity_.col(phase)(i) = 0.; + if (std::abs(this->velocity(phase)(i)) < tolerance) + vector_properties_.at(mpm::properties::Vector::Velocity).col(phase)(i) = + 0.; for (unsigned i = 0; i < Tdim; ++i) - if (std::abs(acceleration_.col(phase)(i)) < tolerance) - acceleration_.col(phase)(i) = 0.; + if (std::abs(this->acceleration(phase)(i)) < tolerance) + vector_properties_.at(mpm::properties::Vector::Acceleration) + .col(phase)(i) = 0.; status = true; } return status; @@ -321,11 +425,13 @@ void mpm::Node::apply_velocity_constraints() { // Phase: Integer value of division (dir / Tdim) const auto phase = static_cast(dir / Tdim); - if (!generic_boundary_constraints_) { + if (!this->boolean_property(mpm::properties::Boolean::GenericBC)) { // Velocity constraints are applied on Cartesian boundaries - this->velocity_(direction, phase) = constraint.second; + vector_properties_.at(mpm::properties::Vector::Velocity)( + direction, phase) = constraint.second; // Set acceleration to 0 in direction of velocity constraint - this->acceleration_(direction, phase) = 0.; + vector_properties_.at(mpm::properties::Vector::Acceleration)(direction, + phase) = 0.; } else { // Velocity constraints on general boundaries // Compute inverse rotation matrix @@ -333,15 +439,19 @@ void mpm::Node::apply_velocity_constraints() { rotation_matrix_.inverse(); // Transform to local coordinate Eigen::Matrix local_velocity = - inverse_rotation_matrix * this->velocity_; + inverse_rotation_matrix * + vector_properties_.at(mpm::properties::Vector::Velocity); Eigen::Matrix local_acceleration = - inverse_rotation_matrix * this->acceleration_; + inverse_rotation_matrix * + vector_properties_.at(mpm::properties::Vector::Acceleration); // Apply boundary condition in local coordinate local_velocity(direction, phase) = constraint.second; local_acceleration(direction, phase) = 0.; // Transform back to global coordinate - this->velocity_ = rotation_matrix_ * local_velocity; - this->acceleration_ = rotation_matrix_ * local_acceleration; + vector_properties_.at(mpm::properties::Vector::Velocity) = + rotation_matrix_ * local_velocity; + vector_properties_.at(mpm::properties::Vector::Acceleration) = + rotation_matrix_ * local_acceleration; } } } @@ -358,7 +468,7 @@ bool mpm::Node::assign_friction_constraint( this->friction_constraint_ = std::make_tuple(static_cast(dir), static_cast(sign_n), static_cast(friction)); - this->friction_ = true; + this->assign_boolean_property(mpm::properties::Boolean::Friction, true); } else throw std::runtime_error("Constraint direction is out of bounds"); @@ -372,7 +482,7 @@ bool mpm::Node::assign_friction_constraint( //! Apply friction constraints template void mpm::Node::apply_friction_constraints(double dt) { - if (friction_) { + if (this->boolean_property(mpm::properties::Boolean::Friction)) { auto sign = [](double value) { return (value > 0.) ? 1. : -1.; }; // Set friction constraint @@ -394,13 +504,16 @@ void mpm::Node::apply_friction_constraints(double dt) { // tangential direction to boundary const unsigned dir_t = (Tdim - 1) - dir_n; - if (!generic_boundary_constraints_) { + if (!this->boolean_property(mpm::properties::Boolean::GenericBC)) { // Cartesian case // Normal and tangential acceleration - acc_n = this->acceleration_(dir_n, phase); - acc_t = this->acceleration_(dir_t, phase); + acc_n = vector_properties_.at(mpm::properties::Vector::Acceleration)( + dir_n, phase); + acc_t = vector_properties_.at(mpm::properties::Vector::Acceleration)( + dir_t, phase); // Velocity tangential - vel_t = this->velocity_(dir_t, phase); + vel_t = vector_properties_.at(mpm::properties::Vector::Velocity)(dir_t, + phase); } else { // General case, transform to local coordinate // Compute inverse rotation matrix @@ -408,9 +521,11 @@ void mpm::Node::apply_friction_constraints(double dt) { rotation_matrix_.inverse(); // Transform to local coordinate Eigen::Matrix local_acceleration = - inverse_rotation_matrix * this->acceleration_; + inverse_rotation_matrix * + vector_properties_.at(mpm::properties::Vector::Acceleration); Eigen::Matrix local_velocity = - inverse_rotation_matrix * this->velocity_; + inverse_rotation_matrix * + vector_properties_.at(mpm::properties::Vector::Velocity); // Normal and tangential acceleration acc_n = local_acceleration(dir_n, phase); acc_t = local_acceleration(dir_t, phase); @@ -433,9 +548,10 @@ void mpm::Node::apply_friction_constraints(double dt) { acc_t -= sign(acc_t) * mu * std::abs(acc_n); } - if (!generic_boundary_constraints_) { + if (!this->boolean_property(mpm::properties::Boolean::GenericBC)) { // Cartesian case - this->acceleration_(dir_t, phase) = acc_t; + vector_properties_.at(mpm::properties::Vector::Acceleration)( + dir_t, phase) = acc_t; } else { // Local acceleration in terms of tangential and normal Eigen::Matrix acc; @@ -443,7 +559,8 @@ void mpm::Node::apply_friction_constraints(double dt) { acc(dir_n, phase) = acc_n; // General case, transform to global coordinate - this->acceleration_.col(phase) = rotation_matrix_ * acc.col(phase); + vector_properties_.at(mpm::properties::Vector::Acceleration) + .col(phase) = rotation_matrix_ * acc.col(phase); } } } else if (Tdim == 3) { @@ -459,18 +576,18 @@ void mpm::Node::apply_friction_constraints(double dt) { const unsigned dir_t1 = dir(dir_n, 1); Eigen::Matrix acc, vel; - if (!generic_boundary_constraints_) { + if (!this->boolean_property(mpm::properties::Boolean::GenericBC)) { // Cartesian case - acc = this->acceleration_.col(phase); - vel = this->velocity_.col(phase); + acc = this->acceleration(phase); + vel = this->velocity(phase); } else { // General case, transform to local coordinate // Compute inverse rotation matrix const Eigen::Matrix inverse_rotation_matrix = rotation_matrix_.inverse(); // Transform to local coordinate - acc = inverse_rotation_matrix * this->acceleration_.col(phase); - vel = inverse_rotation_matrix * this->velocity_.col(phase); + acc = inverse_rotation_matrix * this->acceleration(phase); + vel = inverse_rotation_matrix * this->velocity(phase); } const auto acc_n = acc(dir_n); @@ -508,12 +625,14 @@ void mpm::Node::apply_friction_constraints(double dt) { } } - if (!generic_boundary_constraints_) { + if (!this->boolean_property(mpm::properties::Boolean::GenericBC)) { // Cartesian case - this->acceleration_.col(phase) = acc; + vector_properties_.at(mpm::properties::Vector::Acceleration) + .col(phase) = acc; } else { // General case, transform to global coordinate - this->acceleration_.col(phase) = rotation_matrix_ * acc; + vector_properties_.at(mpm::properties::Vector::Acceleration) + .col(phase) = rotation_matrix_ * acc; } } } @@ -563,7 +682,8 @@ void mpm::Node momentum = property_handle_->property("momenta", prop_id_, *mitr, Tdim); const Eigen::Matrix change_in_momenta = - velocity_ * mass - momentum; + vector_properties_.at(mpm::properties::Vector::Velocity) * mass - + momentum; property_handle_->update_property("change_in_momenta", prop_id_, *mitr, change_in_momenta, Tdim); } @@ -585,7 +705,9 @@ void mpm::Nodeproperty("masses", prop_id_, *mitr); // displacement of the center of mass - contact_displacement_ += material_displacement / mass_(0, 0); + contact_displacement_ += + material_displacement / + scalar_properties_.at(mpm::properties::Scalar::Mass)(0, 0); // assign nodal-multimaterial displacement by dividing it by this material's // mass property_handle_->assign_property( @@ -603,8 +725,10 @@ void mpm::Nodeupdate_property("separation_vectors", prop_id_, *mitr, separation_vector, Tdim); } diff --git a/include/node_base.h b/include/node_base.h index 17eddfdca..e305a193f 100644 --- a/include/node_base.h +++ b/include/node_base.h @@ -10,9 +10,11 @@ #include #include +#include #include "data_types.h" #include "function_base.h" +#include "mpm_properties.h" #include "nodal_properties.h" namespace mpm { @@ -70,6 +72,49 @@ class NodeBase { //! Return status virtual bool status() const = 0; + //! Assign boolean property at the nodes + //! \param[in] property Name of the property to assign + //! \param[in] boolean Property boolean (true/false) of the node + virtual void assign_boolean_property(mpm::properties::Boolean property, + bool boolean) noexcept = 0; + + //! Return boolean property + //! \param[in] property Name of the property to update + //! \retval boolean property at node + virtual bool boolean_property(mpm::properties::Boolean property) const = 0; + + //! Update scalar property at the nodes + //! \param[in] property Name of the property to update + //! \param[in] update A boolean to update (true) or assign (false) + //! \param[in] phase Index corresponding to the phase + //! \param[in] value Property value from the particles in a cell + virtual void update_scalar_property(mpm::properties::Scalar property, + bool update, unsigned phase, + double value) noexcept = 0; + + //! Return property at a given node for a given phase + //! \param[in] property Name of the property to return + //! \param[in] phase Index corresponding to the phase + //! \retval scalar property at the designated phase + virtual double scalar_property(mpm::properties::Scalar property, + unsigned phase) const = 0; + + //! Update vector property at the nodes + //! \param[in] property Name of the property to update + //! \param[in] update A boolean to update (true) or assign (false) + //! \param[in] phase Index corresponding to the phase + //! \param[in] value Property value from the particles in a cell + virtual void update_vector_property( + mpm::properties::Vector property, bool update, unsigned phase, + const Eigen::Matrix& value) noexcept = 0; + + //! Return property at a given node for a given phase + //! \param[in] property Name of the property to return + //! \param[in] phase Index corresponding to the phase + //! \retval vector property at the designated phase + virtual Eigen::Matrix vector_property( + mpm::properties::Vector property, unsigned phase) const = 0; + //! Update mass at the nodes from particle //! \param[in] update A boolean to update (true) or assign (false) //! \param[in] phase Index corresponding to the phase @@ -128,16 +173,21 @@ class NodeBase { virtual VectorDim internal_force(unsigned phase) const = 0; //! Update pressure at the nodes from particle + //! \param[in] update A boolean to update (true) or assign (false) //! \param[in] phase Index corresponding to the phase //! \param[in] mass_pressure Product of mass x pressure of a particle - virtual void update_mass_pressure(unsigned phase, + virtual void update_mass_pressure(bool update, unsigned phase, double mass_pressure) noexcept = 0; + //! Compute pressure from the mass pressure + virtual void compute_pressure() = 0; + //! Assign pressure at the nodes from particle //! \param[in] update A boolean to update (true) or assign (false) //! \param[in] phase Index corresponding to the phase - //! \param[in] mass_pressure Product of mass x pressure of a particle - virtual void assign_pressure(unsigned phase, double mass_pressure) = 0; + //! \param[in] pressure Pressure of a particle + virtual void update_pressure(bool update, unsigned phase, + double pressure) = 0; //! Return pressure at a given node for a given phase //! \param[in] phase Index corresponding to the phase diff --git a/include/particles/particle.h b/include/particles/particle.h index 849b72f56..1488d7761 100644 --- a/include/particles/particle.h +++ b/include/particles/particle.h @@ -55,7 +55,7 @@ class Particle : public ParticleBase { //! \param[in] particle HDF5 data of particle //! \param[in] material Material associated with the particle //! \retval status Status of reading HDF5 particle - virtual bool initialise_particle( + bool initialise_particle( const HDF5Particle& particle, const std::shared_ptr>& material) override; @@ -105,6 +105,9 @@ class Particle : public ParticleBase { //! Return cell id Index cell_id() const override { return cell_id_; } + //! Return cell ptr + std::shared_ptr> cell() const override { return cell_; } + //! Return cell ptr status bool cell_ptr() const override { return cell_ != nullptr; } @@ -119,26 +122,17 @@ class Particle : public ParticleBase { bool assign_volume(double volume) override; //! Return volume - double volume() const override { return volume_; } + double volume() const override { + return this->scalar_property(mpm::properties::Scalar::Volume); + } //! Return size of particle in natural coordinates VectorDim natural_size() const override { return natural_size_; } - //! Compute volume as cell volume / nparticles - void compute_volume() noexcept override; - - //! Update volume based on centre volumetric strain rate - void update_volume() noexcept override; - //! Return mass density - //! \param[in] phase Index corresponding to the phase - double mass_density() const override { return mass_density_; } - - //! Compute mass as volume * density - void compute_mass() noexcept override; - - //! Map particle mass and momentum to nodes - void map_mass_momentum_to_nodes() noexcept override; + double mass_density() const override { + return this->scalar_property(mpm::properties::Scalar::MassDensity); + } //! Map multimaterial properties to nodes void map_multimaterial_mass_momentum_to_nodes() noexcept override; @@ -151,11 +145,14 @@ class Particle : public ParticleBase { //! Assign nodal mass to particles //! \param[in] mass Mass from the particles in a cell - //! \retval status Assignment status - void assign_mass(double mass) override { mass_ = mass; } + void assign_mass(double mass) override { + scalar_properties_.at(mpm::properties::Scalar::Mass) = mass; + } //! Return mass of the particles - double mass() const override { return mass_; } + double mass() const override { + return this->scalar_property(mpm::properties::Scalar::Mass); + } //! Assign material //! \param[in] material Pointer to a material @@ -197,23 +194,24 @@ class Particle : public ParticleBase { //! Return stress of the particle Eigen::Matrix stress() const override { return stress_; } - //! Map body force - //! \param[in] pgravity Gravity of a particle - void map_body_force(const VectorDim& pgravity) noexcept override; - //! Map internal force inline void map_internal_force() noexcept override; //! Assign velocity to the particle //! \param[in] velocity A vector of particle velocity - //! \retval status Assignment status - bool assign_velocity(const VectorDim& velocity) override; + void assign_velocity(const VectorDim& velocity) override { + vector_properties_.at(mpm::properties::Vector::Velocity) = velocity; + }; //! Return velocity of the particle - VectorDim velocity() const override { return velocity_; } + VectorDim velocity() const override { + return this->vector_property(mpm::properties::Vector::Velocity); + } //! Return displacement of the particle - VectorDim displacement() const override { return displacement_; } + VectorDim displacement() const override { + return this->vector_property(mpm::properties::Vector::Displacement); + } //! Assign traction to the particle //! \param[in] direction Index corresponding to the direction of traction @@ -225,15 +223,23 @@ class Particle : public ParticleBase { //! \param[in] phase Index corresponding to the phase VectorDim traction() const override { return traction_; } - //! Map traction force - void map_traction_force() noexcept override; - //! Compute updated position of the particle //! \param[in] dt Analysis time step //! \param[in] velocity_update Update particle velocity from nodal vel void compute_updated_position(double dt, bool velocity_update = false) noexcept override; + //! Assign a state variable + //! \param[in] var State variable + //! \param[in] value State variable to be assigned + //! \param[in] phase Index to indicate phase + void assign_state_variable( + const std::string& var, double value, + unsigned phase = mpm::ParticlePhase::Solid) override { + if (state_variables_[phase].find(var) != state_variables_[phase].end()) + state_variables_[phase].at(var) = value; + } + //! Return a state variable //! \param[in] var State variable //! \param[in] phase Index to indicate phase @@ -246,14 +252,13 @@ class Particle : public ParticleBase { : std::numeric_limits::quiet_NaN(); } - //! Map particle pressure to nodes - bool map_pressure_to_nodes( - unsigned phase = mpm::ParticlePhase::Solid) noexcept override; - - //! Compute pressure smoothing of the particle based on nodal pressure - //! $$\hat{p}_p = \sum_{i = 1}^{n_n} N_i(x_p) p_i$$ - bool compute_pressure_smoothing( - unsigned phase = mpm::ParticlePhase::Solid) noexcept override; + //! Assign a state variable + //! \param[in] value Particle pressure to be assigned + //! \param[in] phase Index to indicate phase + void assign_pressure(double pressure, + unsigned phase = mpm::ParticlePhase::Solid) override { + this->assign_state_variable("pressure", pressure, phase); + } //! Return pressure of the particles //! \param[in] phase Index to indicate phase @@ -329,12 +334,14 @@ class Particle : public ParticleBase { using ParticleBase::state_variables_; //! Neighbour particles using ParticleBase::neighbours_; - //! Volumetric mass density (mass / volume) - double mass_density_{0.}; - //! Mass - double mass_{0.}; - //! Volume - double volume_{0.}; + //! Scalar properties + using ParticleBase::boolean_properties_; + //! Scalar properties + using ParticleBase::scalar_properties_; + //! Vector properties + using ParticleBase::vector_properties_; + //! Shape functions + using ParticleBase::shapefn_; //! Size of particle Eigen::Matrix size_; //! Size of particle in natural coordinates @@ -351,18 +358,10 @@ class Particle : public ParticleBase { Eigen::Matrix strain_rate_; //! dstrains Eigen::Matrix dstrain_; - //! Velocity - Eigen::Matrix velocity_; - //! Displacement - Eigen::Matrix displacement_; //! Particle velocity constraints std::map particle_velocity_constraints_; - //! Set traction - bool set_traction_{false}; //! Surface Traction (given as a stress; force/area) Eigen::Matrix traction_; - //! Shape functions - Eigen::VectorXd shapefn_; //! dN/dX Eigen::MatrixXd dn_dx_; //! dN/dX at cell centroid @@ -376,5 +375,6 @@ class Particle : public ParticleBase { } // namespace mpm #include "particle.tcc" +#include "particle_functions.tcc" #endif // MPM_PARTICLE_H__ diff --git a/include/particles/particle.tcc b/include/particles/particle.tcc index c3d945b69..b74a5629b 100644 --- a/include/particles/particle.tcc +++ b/include/particles/particle.tcc @@ -37,11 +37,12 @@ bool mpm::Particle::initialise_particle(const HDF5Particle& particle) { // Assign id this->id_ = particle.id; // Mass - this->mass_ = particle.mass; + scalar_properties_.at(mpm::properties::Scalar::Mass) = particle.mass; // Volume - this->volume_ = particle.volume; + scalar_properties_.at(mpm::properties::Scalar::Volume) = particle.volume; // Mass Density - this->mass_density_ = particle.mass / particle.volume; + scalar_properties_.at(mpm::properties::Scalar::MassDensity) = + particle.mass / particle.volume; // Set local size of particle Eigen::Vector3d psize; psize << particle.nsize_x, particle.nsize_y, particle.nsize_z; @@ -59,13 +60,16 @@ bool mpm::Particle::initialise_particle(const HDF5Particle& particle) { displacement << particle.displacement_x, particle.displacement_y, particle.displacement_z; // Initialise displacement - for (unsigned i = 0; i < Tdim; ++i) this->displacement_(i) = displacement(i); + for (unsigned i = 0; i < Tdim; ++i) + vector_properties_.at(mpm::properties::Vector::Displacement)(i) = + displacement(i); // Velocity Eigen::Vector3d velocity; velocity << particle.velocity_x, particle.velocity_y, particle.velocity_z; // Initialise velocity - for (unsigned i = 0; i < Tdim; ++i) this->velocity_(i) = velocity(i); + for (unsigned i = 0; i < Tdim; ++i) + vector_properties_.at(mpm::properties::Vector::Velocity)(i) = velocity(i); // Stress this->stress_[0] = particle.stress_xx; @@ -145,11 +149,11 @@ mpm::HDF5Particle mpm::Particle::hdf5() const { Eigen::Vector3d displacement; displacement.setZero(); - for (unsigned j = 0; j < Tdim; ++j) displacement[j] = this->displacement_[j]; + for (unsigned j = 0; j < Tdim; ++j) displacement[j] = this->displacement()[j]; Eigen::Vector3d velocity; velocity.setZero(); - for (unsigned j = 0; j < Tdim; ++j) velocity[j] = this->velocity_[j]; + for (unsigned j = 0; j < Tdim; ++j) velocity[j] = this->velocity()[j]; // Particle local size Eigen::Vector3d nsize; @@ -229,20 +233,33 @@ mpm::HDF5Particle mpm::Particle::hdf5() const { // Initialise particle properties template void mpm::Particle::initialise() { - displacement_.setZero(); dstrain_.setZero(); - mass_ = 0.; natural_size_.setZero(); - set_traction_ = false; size_.setZero(); strain_rate_.setZero(); strain_.setZero(); stress_.setZero(); traction_.setZero(); - velocity_.setZero(); - volume_ = std::numeric_limits::max(); volumetric_strain_centroid_ = 0.; + // Initialize boolean properties + boolean_properties_.emplace( + std::make_pair(mpm::properties::Boolean::SetTraction, false)); + + // Initialize scalar properties + scalar_properties_.emplace( + std::make_pair(mpm::properties::Scalar::Mass, double(0.))); + scalar_properties_.emplace( + std::make_pair(mpm::properties::Scalar::MassDensity, double(0.))); + scalar_properties_.emplace(std::make_pair( + mpm::properties::Scalar::Volume, std::numeric_limits::max())); + + // Initialize vector properties + vector_properties_.emplace( + std::make_pair(mpm::properties::Vector::Displacement, VectorDim::Zero())); + vector_properties_.emplace( + std::make_pair(mpm::properties::Vector::Velocity, VectorDim::Zero())); + // Initialize vector data properties this->properties_["stresses"] = [&]() { return stress(); }; this->properties_["strains"] = [&]() { return strain(); }; @@ -452,10 +469,10 @@ bool mpm::Particle::assign_volume(double volume) { if (volume <= 0.) throw std::runtime_error("Particle volume cannot be negative"); - this->volume_ = volume; + scalar_properties_.at(mpm::properties::Scalar::Volume) = volume; // Compute size of particle in each direction const double length = - std::pow(this->volume_, static_cast(1. / Tdim)); + std::pow(this->volume(), static_cast(1. / Tdim)); // Set particle size as length on each side this->size_.fill(length); @@ -476,99 +493,57 @@ bool mpm::Particle::assign_volume(double volume) { return status; } -// Compute volume of the particle -template -void mpm::Particle::compute_volume() noexcept { - // Check if particle has a valid cell ptr - assert(cell_ != nullptr); - // Volume of the cell / # of particles - this->assign_volume(cell_->volume() / cell_->nparticles()); -} - -// Update volume based on the central strain rate -template -void mpm::Particle::update_volume() noexcept { - // Check if particle has a valid cell ptr and a valid volume - assert(cell_ != nullptr && volume_ != std::numeric_limits::max()); - // Compute at centroid - // Strain rate for reduced integration - this->volume_ *= (1. + dvolumetric_strain_); - this->mass_density_ = this->mass_density_ / (1. + dvolumetric_strain_); -} - -// Compute mass of particle -template -void mpm::Particle::compute_mass() noexcept { - // Check if particle volume is set and material ptr is valid - assert(volume_ != std::numeric_limits::max() && - this->material() != nullptr); - // Mass = volume of particle * mass_density - this->mass_density_ = - (this->material())->template property(std::string("density")); - this->mass_ = volume_ * mass_density_; -} - -//! Map particle mass and momentum to nodes -template -void mpm::Particle::map_mass_momentum_to_nodes() noexcept { - // Check if particle mass is set - assert(mass_ != std::numeric_limits::max()); - - // Map mass and momentum to nodes - for (unsigned i = 0; i < nodes_.size(); ++i) { - nodes_[i]->update_mass(true, mpm::ParticlePhase::Solid, - mass_ * shapefn_[i]); - nodes_[i]->update_momentum(true, mpm::ParticlePhase::Solid, - mass_ * shapefn_[i] * velocity_); - } -} - //! Map multimaterial properties to nodes +// TODO: Contact function to be refactored template void mpm::Particle::map_multimaterial_mass_momentum_to_nodes() noexcept { // Check if particle mass is set - assert(mass_ != std::numeric_limits::max()); + assert(this->mass() != std::numeric_limits::max()); // Unit 1x1 Eigen matrix to be used with scalar quantities Eigen::Matrix nodal_mass; // Map mass and momentum to nodal property taking into account the material id for (unsigned i = 0; i < nodes_.size(); ++i) { - nodal_mass(0, 0) = mass_ * shapefn_[i]; + nodal_mass(0, 0) = this->mass() * shapefn_[i]; nodes_[i]->update_property(true, "masses", nodal_mass, this->material_id(), 1); - nodes_[i]->update_property(true, "momenta", velocity_ * nodal_mass, + nodes_[i]->update_property(true, "momenta", this->velocity() * nodal_mass, this->material_id(), Tdim); } } //! Map multimaterial displacements to nodes +// TODO: Contact function to be refactored template void mpm::Particle::map_multimaterial_displacements_to_nodes() noexcept { // Check if particle mass is set - assert(mass_ != std::numeric_limits::max()); + assert(this->mass() != std::numeric_limits::max()); // Map displacements to nodal property and divide it by the respective // nodal-material mass for (unsigned i = 0; i < nodes_.size(); ++i) { - const auto& displacement = mass_ * shapefn_[i] * displacement_; + const auto& displacement = + this->mass() * shapefn_[i] * this->displacement(); nodes_[i]->update_property(true, "displacements", displacement, this->material_id(), Tdim); } } //! Map multimaterial domain gradients to nodes +// TODO: Contact function to be refactored template void mpm::Particle< Tdim>::map_multimaterial_domain_gradients_to_nodes() noexcept { // Check if particle volume is set - assert(volume_ != std::numeric_limits::max()); + assert(this->volume() != std::numeric_limits::max()); // Map domain gradients to nodal property. The domain gradients is defined as // the gradient of the particle volume for (unsigned i = 0; i < nodes_.size(); ++i) { Eigen::Matrix gradient; - for (unsigned j = 0; j < Tdim; ++j) gradient[j] = volume_ * dn_dx_(i, j); + for (unsigned j = 0; j < Tdim; ++j) + gradient[j] = this->volume() * dn_dx_(i, j); nodes_[i]->update_property(true, "domain_gradients", gradient, this->material_id(), Tdim); } @@ -664,15 +639,6 @@ void mpm::Particle::compute_stress() noexcept { &state_variables_[mpm::ParticlePhase::Solid]); } -//! Map body force -template -void mpm::Particle::map_body_force(const VectorDim& pgravity) noexcept { - // Compute nodal body forces - for (unsigned i = 0; i < nodes_.size(); ++i) - nodes_[i]->update_external_force(true, mpm::ParticlePhase::Solid, - (pgravity * mass_ * shapefn_(i))); -} - //! Map internal force template <> inline void mpm::Particle<1>::map_internal_force() noexcept { @@ -680,7 +646,7 @@ inline void mpm::Particle<1>::map_internal_force() noexcept { for (unsigned i = 0; i < nodes_.size(); ++i) { // Compute force: -pstress * volume Eigen::Matrix force; - force[0] = -1. * dn_dx_(i, 0) * volume_ * stress_[0]; + force[0] = -1. * dn_dx_(i, 0) * this->volume() * stress_[0]; nodes_[i]->update_internal_force(true, mpm::ParticlePhase::Solid, force); } @@ -696,7 +662,7 @@ inline void mpm::Particle<2>::map_internal_force() noexcept { force[0] = dn_dx_(i, 0) * stress_[0] + dn_dx_(i, 1) * stress_[3]; force[1] = dn_dx_(i, 1) * stress_[1] + dn_dx_(i, 0) * stress_[3]; - force *= -1. * this->volume_; + force *= -1. * this->volume(); nodes_[i]->update_internal_force(true, mpm::ParticlePhase::Solid, force); } @@ -718,35 +684,26 @@ inline void mpm::Particle<3>::map_internal_force() noexcept { force[2] = dn_dx_(i, 2) * stress_[2] + dn_dx_(i, 1) * stress_[4] + dn_dx_(i, 0) * stress_[5]; - force *= -1. * this->volume_; + force *= -1. * this->volume(); nodes_[i]->update_internal_force(true, mpm::ParticlePhase::Solid, force); } } -// Assign velocity to the particle -template -bool mpm::Particle::assign_velocity( - const Eigen::Matrix& velocity) { - // Assign velocity - velocity_ = velocity; - return true; -} - // Assign traction to the particle template bool mpm::Particle::assign_traction(unsigned direction, double traction) { bool status = false; try { if (direction >= Tdim || - this->volume_ == std::numeric_limits::max()) { + this->volume() == std::numeric_limits::max()) { throw std::runtime_error( "Particle traction property: volume / direction is invalid"); } // Assign traction - traction_(direction) = traction * this->volume_ / this->size_(direction); + traction_(direction) = traction * this->volume() / this->size_(direction); status = true; - this->set_traction_ = true; + this->assign_boolean_property(mpm::properties::Boolean::SetTraction, true); } catch (std::exception& exception) { console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); status = false; @@ -754,17 +711,6 @@ bool mpm::Particle::assign_traction(unsigned direction, double traction) { return status; } -//! Map traction force -template -void mpm::Particle::map_traction_force() noexcept { - if (this->set_traction_) { - // Map particle traction forces to nodes - for (unsigned i = 0; i < nodes_.size(); ++i) - nodes_[i]->update_external_force(true, mpm::ParticlePhase::Solid, - (shapefn_[i] * traction_)); - } -} - // Compute updated position of the particle template void mpm::Particle::compute_updated_position( @@ -772,84 +718,40 @@ void mpm::Particle::compute_updated_position( // Check if particle has a valid cell ptr assert(cell_ != nullptr); // Get interpolated nodal velocity - Eigen::Matrix nodal_velocity = - Eigen::Matrix::Zero(); - - for (unsigned i = 0; i < nodes_.size(); ++i) - nodal_velocity += - shapefn_[i] * nodes_[i]->velocity(mpm::ParticlePhase::Solid); + const auto& nodal_velocity = this->interpolate_vector_property_from_nodes( + mpm::properties::Vector::Velocity, mpm::ParticlePhase::Solid); // Acceleration update if (!velocity_update) { // Get interpolated nodal acceleration - Eigen::Matrix nodal_acceleration = - Eigen::Matrix::Zero(); - for (unsigned i = 0; i < nodes_.size(); ++i) - nodal_acceleration += - shapefn_[i] * nodes_[i]->acceleration(mpm::ParticlePhase::Solid); + const auto& nodal_acceleration = + this->interpolate_vector_property_from_nodes( + mpm::properties::Vector::Acceleration, mpm::ParticlePhase::Solid); // Update particle velocity from interpolated nodal acceleration - this->velocity_ += nodal_acceleration * dt; + vector_properties_.at(mpm::properties::Vector::Velocity) += + nodal_acceleration * dt; } // Update particle velocity using interpolated nodal velocity else - this->velocity_ = nodal_velocity; + vector_properties_.at(mpm::properties::Vector::Velocity) = nodal_velocity; // New position current position + velocity * dt this->coordinates_ += nodal_velocity * dt; - // Update displacement (displacement is initialized from zero) - this->displacement_ += nodal_velocity * dt; -} - -//! Map particle pressure to nodes -template -bool mpm::Particle::map_pressure_to_nodes(unsigned phase) noexcept { - // Mass is initialized - assert(mass_ != std::numeric_limits::max()); - bool status = false; - // Check if particle mass is set and state variable pressure is found - if (mass_ != std::numeric_limits::max() && - (state_variables_[phase].find("pressure") != - state_variables_[phase].end())) { - // Map particle pressure to nodes - for (unsigned i = 0; i < nodes_.size(); ++i) - nodes_[i]->update_mass_pressure( - phase, shapefn_[i] * mass_ * state_variables_[phase]["pressure"]); - - status = true; - } - return status; -} - -// Compute pressure smoothing of the particle based on nodal pressure -template -bool mpm::Particle::compute_pressure_smoothing(unsigned phase) noexcept { - // Assert - assert(cell_ != nullptr); - - bool status = false; - // Check if particle has a valid cell ptr - if (cell_ != nullptr && (state_variables_[phase].find("pressure") != - state_variables_[phase].end())) { - - double pressure = 0.; - // Update particle pressure to interpolated nodal pressure - for (unsigned i = 0; i < this->nodes_.size(); ++i) - pressure += shapefn_[i] * nodes_[i]->pressure(phase); - - state_variables_[phase]["pressure"] = pressure; - status = true; - } - return status; + // Update displacement (displacement is initialized from zero) + vector_properties_.at(mpm::properties::Vector::Displacement) += + nodal_velocity * dt; } //! Apply particle velocity constraints +// TODO: revisit constraint application in vector properties, should assign +// direction template void mpm::Particle::apply_particle_velocity_constraints(unsigned dir, double velocity) { // Set particle velocity constraint - this->velocity_(dir) = velocity; + vector_properties_.at(mpm::properties::Vector::Velocity)(dir) = velocity; } //! Return particle tensor data @@ -859,6 +761,7 @@ Eigen::VectorXd mpm::Particle::tensor_data(const std::string& property) { } //! Assign material id of this particle to nodes +// TODO: Contact function to be refactored template void mpm::Particle::append_material_id_to_nodes() const { for (unsigned i = 0; i < nodes_.size(); ++i) diff --git a/include/particles/particle_base.h b/include/particles/particle_base.h index da0d22a99..b4b719890 100644 --- a/include/particles/particle_base.h +++ b/include/particles/particle_base.h @@ -108,6 +108,9 @@ class ParticleBase { //! Return cell id virtual Index cell_id() const = 0; + //! Return cell ptr + virtual std::shared_ptr> cell() const = 0; + //! Return cell ptr status virtual bool cell_ptr() const = 0; @@ -126,21 +129,93 @@ class ParticleBase { //! Return size of particle in natural coordinates virtual VectorDim natural_size() const = 0; - //! Compute volume of particle - virtual void compute_volume() noexcept = 0; - - //! Update volume based on centre volumetric strain rate - virtual void update_volume() noexcept = 0; + //! Assign boolean property at the particle + //! \param[in] property Name of the property to assign + //! \param[in] boolean Property boolean (true/false) of the particles in a + //! cell + void assign_boolean_property(mpm::properties::Boolean property, + bool boolean) noexcept; + + //! Return boolean property + //! \param[in] property Name of the property to update + //! \retval boolean property at particle + bool boolean_property(mpm::properties::Boolean property) const; + + //! Update scalar property at the particle + //! \param[in] property Name of the property to update + //! \param[in] update A boolean to update (true) or assign (false) + //! \param[in] value Property value from the particles in a cell + void update_scalar_property(mpm::properties::Scalar property, bool update, + double value) noexcept; + + //! Return scalar property + //! \param[in] property Name of the property to return + //! \retval scalar property at particle + double scalar_property(mpm::properties::Scalar property) const; + + //! Map scalar property to the nodes + //! \param[in] property Name of the property to update + //! \param[in] update A boolean to update (true) or assign (false) + //! \param[in] phase Index corresponding to the phase + void map_scalar_property_to_nodes(mpm::properties::Scalar property, + bool update, unsigned phase) noexcept; + + //! Map an arbitrary scalar value to nodal scalar property + //! \param[in] property Name of the property to update + //! \param[in] update A boolean to update (true) or assign (false) + //! \param[in] phase Index corresponding to the phase + //! \param[in] value Scalar value to be mapped from particle to node + void map_scalar_property_to_nodes(mpm::properties::Scalar property, + bool update, unsigned phase, + double value) noexcept; + + //! Return an interpolation of scalar property in particle from nodes + //! \param[in] property Name of the property to update + //! \param[in] phase Index corresponding to the phase + //! \retval interpolated scalar property at particle + double interpolate_scalar_property_from_nodes( + mpm::properties::Scalar property, unsigned phase) const; + + //! Update vector property at the particle + //! \param[in] property Name of the property to update + //! \param[in] update A boolean to update (true) or assign (false) + //! \param[in] value Property value from the particles in a cell + void update_vector_property( + mpm::properties::Vector property, bool update, + const Eigen::Matrix& value) noexcept; + + //! Return vector property + //! \param[in] property Name of the property to return + //! \retval vector property at particle + Eigen::Matrix vector_property( + mpm::properties::Vector property) const; + + //! Map vector property to the nodes + //! \param[in] property Name of the property to update + //! \param[in] update A boolean to update (true) or assign (false) + //! \param[in] phase Index corresponding to the phase + void map_vector_property_to_nodes(mpm::properties::Vector property, + bool update, unsigned phase) noexcept; + + //! Map an arbitrary vector value to nodal vector property + //! \param[in] property Name of the property to update + //! \param[in] update A boolean to update (true) or assign (false) + //! \param[in] phase Index corresponding to the phase + //! \param[in] value Vector value to be mapped from particle to node + void map_vector_property_to_nodes( + mpm::properties::Vector property, bool update, unsigned phase, + const Eigen::Matrix& value) noexcept; + + //! Return an interpolation of vector property in particle from nodes + //! \param[in] property Name of the property to update + //! \param[in] phase Index corresponding to the phase + //! \retval interpolated vector property at particle + Eigen::Matrix interpolate_vector_property_from_nodes( + mpm::properties::Vector property, unsigned phase) const; //! Return mass density virtual double mass_density() const = 0; - //! Compute mass of particle - virtual void compute_mass() noexcept = 0; - - //! Map particle mass and momentum to nodes - virtual void map_mass_momentum_to_nodes() noexcept = 0; - //! Map multimaterial properties to nodes virtual void map_multimaterial_mass_momentum_to_nodes() noexcept = 0; @@ -174,6 +249,16 @@ class ParticleBase { return state_variables_[phase]; } + //! Assign a state variable + virtual void assign_state_variable( + const std::string& var, double value, + unsigned phase = mpm::ParticlePhase::Solid) = 0; + + //! Return a state variable + virtual double state_variable( + const std::string& var, + unsigned phase = mpm::ParticlePhase::Solid) const = 0; + //! Assign status void assign_status(bool status) { status_ = status; } @@ -189,6 +274,10 @@ class ParticleBase { //! Return mass virtual double mass() const = 0; + //! Assign pressure + virtual void assign_pressure(double pressure, + unsigned phase = mpm::ParticlePhase::Solid) = 0; + //! Return pressure virtual double pressure(unsigned phase = mpm::ParticlePhase::Solid) const = 0; @@ -216,22 +305,11 @@ class ParticleBase { //! Return stress virtual Eigen::Matrix stress() const = 0; - //! Map body force - virtual void map_body_force(const VectorDim& pgravity) noexcept = 0; - //! Map internal force virtual void map_internal_force() noexcept = 0; - //! Map particle pressure to nodes - virtual bool map_pressure_to_nodes( - unsigned phase = mpm::ParticlePhase::Solid) noexcept = 0; - - //! Compute pressure smoothing of the particle based on nodal pressure - virtual bool compute_pressure_smoothing( - unsigned phase = mpm::ParticlePhase::Solid) noexcept = 0; - //! Assign velocity - virtual bool assign_velocity(const VectorDim& velocity) = 0; + virtual void assign_velocity(const VectorDim& velocity) = 0; //! Return velocity virtual VectorDim velocity() const = 0; @@ -245,18 +323,10 @@ class ParticleBase { //! Return traction virtual VectorDim traction() const = 0; - //! Map traction force - virtual void map_traction_force() noexcept = 0; - //! Compute updated position virtual void compute_updated_position( double dt, bool velocity_update = false) noexcept = 0; - //! Return a state variable - virtual double state_variable( - const std::string& var, - unsigned phase = mpm::ParticlePhase::Solid) const = 0; - //! Return tensor data of particles //! \param[in] property Property string //! \retval vecdata Tensor data of particle property @@ -305,6 +375,15 @@ class ParticleBase { std::vector state_variables_; //! Vector of particle neighbour ids std::vector neighbours_; + //! Shape functions + Eigen::VectorXd shapefn_; + //! Boolean properties + fc::vector_map boolean_properties_; + //! Scalar properties + fc::vector_map scalar_properties_; + //! Vector properties + fc::vector_map> + vector_properties_; }; // ParticleBase class } // namespace mpm diff --git a/include/particles/particle_base.tcc b/include/particles/particle_base.tcc index a75e8ca0b..092340f25 100644 --- a/include/particles/particle_base.tcc +++ b/include/particles/particle_base.tcc @@ -15,3 +15,120 @@ mpm::ParticleBase::ParticleBase(Index id, const VectorDim& coord, : mpm::ParticleBase::ParticleBase(id, coord) { status_ = status; } + +//! Assign boolean property at the particle +template +void mpm::ParticleBase::assign_boolean_property( + mpm::properties::Boolean property, bool boolean) noexcept { + boolean_properties_.at(property) = boolean; +} + +//! Return boolean property +template +bool mpm::ParticleBase::boolean_property( + mpm::properties::Boolean property) const { + return boolean_properties_.at(property); +} + +//! Update scalar property at particle +template +void mpm::ParticleBase::update_scalar_property( + mpm::properties::Scalar property, bool update, double value) noexcept { + // Decide to update or assign + const double factor = (update == true) ? 1. : 0.; + scalar_properties_.at(property) = + scalar_properties_.at(property) * factor + value; +} + +//! Update scalar property at particle +template +double mpm::ParticleBase::scalar_property( + mpm::properties::Scalar property) const { + return scalar_properties_.at(property); +} + +//! Map scalar property to nodes +template +void mpm::ParticleBase::map_scalar_property_to_nodes( + mpm::properties::Scalar property, bool update, unsigned phase) noexcept { + // Check if particle property is set + assert(scalar_properties_.at(property) != std::numeric_limits::max()); + + // Map scalar property to nodes + for (unsigned i = 0; i < nodes_.size(); ++i) + nodes_[i]->update_scalar_property( + property, update, phase, scalar_properties_.at(property) * shapefn_[i]); +} + +//! Map an arbitrary scalar value to nodal scalar property +template +void mpm::ParticleBase::map_scalar_property_to_nodes( + mpm::properties::Scalar property, bool update, unsigned phase, + double value) noexcept { + // Map scalar value to nodes + for (unsigned i = 0; i < nodes_.size(); ++i) + nodes_[i]->update_scalar_property(property, update, phase, + value * shapefn_[i]); +} + +//! Interpolate scalar property from nodes +template +double mpm::ParticleBase::interpolate_scalar_property_from_nodes( + mpm::properties::Scalar property, unsigned phase) const { + double value = 0.; + // Interpolate scalar property from nodes + for (unsigned i = 0; i < nodes_.size(); ++i) + value += nodes_[i]->scalar_property(property, phase) * shapefn_[i]; + return value; +} + +//! Update vector property at particle +template +void mpm::ParticleBase::update_vector_property( + mpm::properties::Vector property, bool update, + const Eigen::Matrix& value) noexcept { + // Decide to update or assign + const double factor = (update == true) ? 1. : 0.; + vector_properties_.at(property) = + vector_properties_.at(property) * factor + value; +} + +//! Update vector property at particle +template +Eigen::Matrix mpm::ParticleBase::vector_property( + mpm::properties::Vector property) const { + return vector_properties_.at(property); +} + +//! Map vector property to nodes +template +void mpm::ParticleBase::map_vector_property_to_nodes( + mpm::properties::Vector property, bool update, unsigned phase) noexcept { + // Map vector property to nodes + for (unsigned i = 0; i < nodes_.size(); ++i) + nodes_[i]->update_vector_property( + property, update, phase, vector_properties_.at(property) * shapefn_[i]); +} + +//! Map an arbitrary vector value to nodal vector property +template +void mpm::ParticleBase::map_vector_property_to_nodes( + mpm::properties::Vector property, bool update, unsigned phase, + const Eigen::Matrix& value) noexcept { + // Map vector property to nodes + for (unsigned i = 0; i < nodes_.size(); ++i) + nodes_[i]->update_vector_property(property, update, phase, + value * shapefn_[i]); +} + +//! Interpolate vector property from nodes +template +Eigen::Matrix + mpm::ParticleBase::interpolate_vector_property_from_nodes( + mpm::properties::Vector property, unsigned phase) const { + Eigen::Matrix value = Eigen::Matrix::Zero(); + // Interpolate vector property from nodes + for (unsigned i = 0; i < nodes_.size(); ++i) + value += nodes_[i]->vector_property(property, phase) * shapefn_[i]; + return value; +} diff --git a/include/particles/particle_functions.tcc b/include/particles/particle_functions.tcc new file mode 100644 index 000000000..fb53fc1ea --- /dev/null +++ b/include/particles/particle_functions.tcc @@ -0,0 +1,123 @@ +// Compute mass of particle +namespace mpm { +namespace particle { + +// Compute particle mass +template +void compute_mass(std::shared_ptr> particle) noexcept { + // Check if particle volume is set and material ptr is valid + assert(particle->volume() != std::numeric_limits::max() && + particle->material() != nullptr); + + // Mass = volume of particle * mass_density + particle->update_scalar_property( + mpm::properties::Scalar::MassDensity, false, + particle->material()->template property(std::string("density"))); + + // Update particle mass + particle->update_scalar_property( + mpm::properties::Scalar::Mass, false, + particle->volume() * particle->mass_density()); +} + +// Compute volume of the particle +template +void compute_volume( + std::shared_ptr> particle) noexcept { + // Check if particle has a valid cell ptr + assert(particle->cell_ptr()); + + // Volume of the cell / # of particles + particle->assign_volume(particle->cell()->volume() / + particle->cell()->nparticles()); +} + +// Update volume based on the central strain rate +template +void update_volume(std::shared_ptr> particle) noexcept { + // Check if particle has a valid cell ptr and a valid volume + assert(particle->cell_ptr() && + particle->volume() != std::numeric_limits::max()); + + // Compute at centroid + // Strain rate for reduced integration + particle->update_scalar_property( + mpm::properties::Scalar::Volume, false, + (particle->volume() * (1. + particle->dvolumetric_strain()))); + particle->update_scalar_property( + mpm::properties::Scalar::MassDensity, false, + (particle->mass_density() / (1. + particle->dvolumetric_strain()))); +} + +//! Map particle mass and momentum to nodes +template +void map_mass_momentum_to_nodes( + std::shared_ptr> particle) noexcept { + // Check if particle mass is set + assert(particle->mass() != std::numeric_limits::max()); + + // Map mass and momentum to nodes + particle->map_scalar_property_to_nodes(mpm::properties::Scalar::Mass, true, + mpm::ParticlePhase::Solid); + particle->map_vector_property_to_nodes( + mpm::properties::Vector::Momentum, true, mpm::ParticlePhase::Solid, + particle->mass() * particle->velocity()); +} + +//! Map particle pressure to nodes +template +void map_mass_pressure_to_nodes( + std::shared_ptr> particle, + unsigned phase = mpm::ParticlePhase::Solid) noexcept { + // Mass is initialized + assert(particle->mass() != std::numeric_limits::max()); + + // Check if state variable pressure is found + if (particle->pressure(phase) != std::numeric_limits::quiet_NaN()) { + // Map particle pressure to nodes + particle->map_scalar_property_to_nodes( + mpm::properties::Scalar::MassPressure, true, phase, + particle->mass() * particle->pressure(phase)); + } +} + +//! Map body force to nodes +template +void map_body_force(std::shared_ptr> particle, + const Eigen::Matrix& pgravity) noexcept { + // Compute nodal body forces + particle->map_vector_property_to_nodes(mpm::properties::Vector::ExternalForce, + true, mpm::ParticlePhase::Solid, + pgravity * particle->mass()); +} + +//! Map traction force +template +void map_traction_force( + std::shared_ptr> particle) noexcept { + if (particle->boolean_property(mpm::properties::Boolean::SetTraction)) { + // Map particle traction forces to nodes + particle->map_vector_property_to_nodes( + mpm::properties::Vector::ExternalForce, true, mpm::ParticlePhase::Solid, + particle->traction()); + } +} + +// Compute pressure smoothing of the particle based on nodal pressure +template +void compute_pressure_smoothing( + std::shared_ptr> particle, + unsigned phase = mpm::ParticlePhase::Solid) noexcept { + // Assert + assert(particle->cell_ptr()); + + // Check if particle has pressure + if (particle->pressure(phase) != std::numeric_limits::quiet_NaN()) { + double pressure = particle->interpolate_scalar_property_from_nodes( + mpm::properties::Scalar::Pressure, phase); + particle->assign_state_variable("pressure", pressure, phase); + } +} + +} // namespace particle +} // namespace mpm diff --git a/include/solvers/mpm_base.tcc b/include/solvers/mpm_base.tcc index 8426ccf4f..65d5031da 100644 --- a/include/solvers/mpm_base.tcc +++ b/include/solvers/mpm_base.tcc @@ -300,8 +300,10 @@ bool mpm::MPMBase::initialise_particles() { auto particles_volume_begin = std::chrono::steady_clock::now(); // Compute volume - mesh_->iterate_over_particles(std::bind( - &mpm::ParticleBase::compute_volume, std::placeholders::_1)); + mesh_->iterate_over_particles( + [](std::shared_ptr> ptr) { + return mpm::particle::compute_volume(ptr); + }); // Read and assign particles volumes this->particles_volumes(mesh_props, particle_io); @@ -1164,10 +1166,16 @@ void mpm::MPMBase::mpi_domain_decompose(bool initial_step) { //! MPM pressure smoothing template void mpm::MPMBase::pressure_smoothing(unsigned phase) { - // Assign pressure to nodes + // Assign mass pressure to nodes mesh_->iterate_over_particles( - std::bind(&mpm::ParticleBase::map_pressure_to_nodes, - std::placeholders::_1, phase)); + [&phase = phase](std::shared_ptr> ptr) { + return mpm::particle::map_mass_pressure_to_nodes(ptr, phase); + }); + + // Compute nodal pressure + mesh_->iterate_over_nodes_predicate( + std::bind(&mpm::NodeBase::compute_pressure, std::placeholders::_1), + std::bind(&mpm::NodeBase::status, std::placeholders::_1)); #ifdef USE_MPI int mpi_size = 1; @@ -1180,13 +1188,14 @@ void mpm::MPMBase::pressure_smoothing(unsigned phase) { // MPI all reduce nodal pressure mesh_->template nodal_halo_exchange( std::bind(&mpm::NodeBase::pressure, std::placeholders::_1, phase), - std::bind(&mpm::NodeBase::assign_pressure, std::placeholders::_1, - phase, std::placeholders::_2)); + std::bind(&mpm::NodeBase::update_pressure, std::placeholders::_1, + false, phase, std::placeholders::_2)); } #endif // Smooth pressure over particles mesh_->iterate_over_particles( - std::bind(&mpm::ParticleBase::compute_pressure_smoothing, - std::placeholders::_1, phase)); + [&phase = phase](std::shared_ptr> ptr) { + return mpm::particle::compute_pressure_smoothing(ptr, phase); + }); } diff --git a/include/solvers/mpm_explicit.tcc b/include/solvers/mpm_explicit.tcc index 67261d041..d72a1cede 100644 --- a/include/solvers/mpm_explicit.tcc +++ b/include/solvers/mpm_explicit.tcc @@ -14,8 +14,10 @@ void mpm::MPMExplicit::compute_stress_strain(unsigned phase) { &mpm::ParticleBase::compute_strain, std::placeholders::_1, dt_)); // Iterate over each particle to update particle volume - mesh_->iterate_over_particles(std::bind( - &mpm::ParticleBase::update_volume, std::placeholders::_1)); + mesh_->iterate_over_particles( + [](std::shared_ptr> ptr) { + return mpm::particle::update_volume(ptr); + }); // Pressure smoothing if (pressure_smoothing_) this->pressure_smoothing(phase); @@ -93,7 +95,9 @@ bool mpm::MPMExplicit::solve() { // Compute mass mesh_->iterate_over_particles( - std::bind(&mpm::ParticleBase::compute_mass, std::placeholders::_1)); + [](std::shared_ptr> ptr) { + return mpm::particle::compute_mass(ptr); + }); // Check point resume if (resume) this->checkpoint_resume(); @@ -153,8 +157,9 @@ bool mpm::MPMExplicit::solve() { // Assign mass and momentum to nodes mesh_->iterate_over_particles( - std::bind(&mpm::ParticleBase::map_mass_momentum_to_nodes, - std::placeholders::_1)); + [](std::shared_ptr> ptr) { + return mpm::particle::map_mass_momentum_to_nodes(ptr); + }); #ifdef USE_MPI // Run if there is more than a single MPI task @@ -223,8 +228,9 @@ bool mpm::MPMExplicit::solve() { { // Iterate over each particle to compute nodal body force mesh_->iterate_over_particles( - std::bind(&mpm::ParticleBase::map_body_force, - std::placeholders::_1, this->gravity_)); + [&value = gravity_](std::shared_ptr> ptr) { + return mpm::particle::map_body_force(ptr, value); + }); // Apply particle traction and map to nodes mesh_->apply_traction_on_particles(this->step_ * this->dt_); diff --git a/tests/interface_test.cc b/tests/interface_test.cc index cb3b151d7..e9104e799 100644 --- a/tests/interface_test.cc +++ b/tests/interface_test.cc @@ -152,9 +152,9 @@ TEST_CASE("Interface functions are checked", "[interface]") { particle3->map_multimaterial_mass_momentum_to_nodes(); // Map masses and momenta from particles to nodes - particle1->map_mass_momentum_to_nodes(); - particle2->map_mass_momentum_to_nodes(); - particle3->map_mass_momentum_to_nodes(); + mpm::particle::map_mass_momentum_to_nodes(particle1); + mpm::particle::map_mass_momentum_to_nodes(particle2); + mpm::particle::map_mass_momentum_to_nodes(particle3); // Compute velocities at nodes node0->compute_velocity(); diff --git a/tests/mesh_test_2d.cc b/tests/mesh_test_2d.cc index c82b04f2e..eba42e693 100644 --- a/tests/mesh_test_2d.cc +++ b/tests/mesh_test_2d.cc @@ -967,8 +967,9 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { // Compute volume mesh->iterate_over_particles( - std::bind(&mpm::ParticleBase::compute_volume, - std::placeholders::_1)); + [](std::shared_ptr> ptr) { + return mpm::particle::compute_volume(ptr); + }); mesh->apply_traction_on_particles(10); } diff --git a/tests/mesh_test_3d.cc b/tests/mesh_test_3d.cc index a6c1db5bb..85d5d2bca 100644 --- a/tests/mesh_test_3d.cc +++ b/tests/mesh_test_3d.cc @@ -1062,8 +1062,9 @@ TEST_CASE("Mesh is checked for 3D case", "[mesh][3D]") { // Compute volume mesh->iterate_over_particles( - std::bind(&mpm::ParticleBase::compute_volume, - std::placeholders::_1)); + [](std::shared_ptr> ptr) { + return mpm::particle::compute_volume(ptr); + }); mesh->apply_traction_on_particles(10); } diff --git a/tests/node_test.cc b/tests/node_test.cc index 3cc92d2a5..3a610136a 100644 --- a/tests/node_test.cc +++ b/tests/node_test.cc @@ -134,30 +134,47 @@ TEST_CASE("Node is checked for 1D case", "[node][1D]") { mass = 100.; REQUIRE_NOTHROW(node->update_mass(false, Nphase, mass)); REQUIRE(node->mass(Nphase) == Approx(100.0).epsilon(Tolerance)); + // Assign mass to 200 using scalar property update true + REQUIRE_NOTHROW(node->update_scalar_property(mpm::properties::Scalar::Mass, + true, Nphase, mass)); + REQUIRE(node->scalar_property(mpm::properties::Scalar::Mass, Nphase) == + Approx(200.0).epsilon(Tolerance)); + // Assign mass to 100 using scalar property update false + REQUIRE_NOTHROW(node->update_scalar_property(mpm::properties::Scalar::Mass, + false, Nphase, mass)); + REQUIRE(node->scalar_property(mpm::properties::Scalar::Mass, Nphase) == + Approx(100.0).epsilon(Tolerance)); SECTION("Check nodal pressure") { // Check pressure REQUIRE(node->pressure(Nphase) == Approx(0.0).epsilon(Tolerance)); double pressure = 1000.7; // Update pressure to 1000.7 - REQUIRE_NOTHROW(node->update_mass_pressure(Nphase, mass * pressure)); + REQUIRE_NOTHROW( + node->update_mass_pressure(false, Nphase, mass * pressure)); + node->compute_pressure(); REQUIRE(node->pressure(Nphase) == Approx(1000.7).epsilon(Tolerance)); // Update pressure to 2001.4 - REQUIRE_NOTHROW(node->update_mass_pressure(Nphase, mass * pressure)); + REQUIRE_NOTHROW( + node->update_mass_pressure(true, Nphase, mass * pressure)); + node->compute_pressure(); REQUIRE(node->pressure(Nphase) == Approx(2001.4).epsilon(Tolerance)); // Assign pressure to 1000 pressure = 1000.; - node->assign_pressure(Nphase, pressure); + REQUIRE_NOTHROW(node->update_pressure(false, Nphase, pressure)); REQUIRE(node->pressure(Nphase) == Approx(1000.0).epsilon(Tolerance)); // Assign mass to 0 mass = 0.; REQUIRE_NOTHROW(node->update_mass(false, Nphase, mass)); REQUIRE(node->mass(Nphase) == Approx(0.0).epsilon(Tolerance)); - // Try to update pressure to 2000, should throw and keep to 1000. + // Try to update pressure, should throw and keep to 1000. pressure = 1000.; - const double pmass = 1.5; - node->assign_pressure(Nphase, pressure); + node->update_mass_pressure(false, Nphase, pressure); + node->compute_pressure(); REQUIRE(node->pressure(Nphase) == Approx(1000.0).epsilon(Tolerance)); + // Update pressure to 2000. + REQUIRE_NOTHROW(node->update_pressure(true, Nphase, pressure)); + REQUIRE(node->pressure(Nphase) == Approx(2000.0).epsilon(Tolerance)); } SECTION("Check external force") { @@ -261,6 +278,22 @@ TEST_CASE("Node is checked for 1D case", "[node][1D]") { for (unsigned i = 0; i < force.size(); ++i) REQUIRE(node->internal_force(Nphase)(i) == Approx(10.).epsilon(Tolerance)); + + // Assign force to 20 using vector property update true + REQUIRE_NOTHROW(node->update_vector_property( + mpm::properties::Vector::InternalForce, true, Nphase, force)); + for (unsigned i = 0; i < force.size(); ++i) + REQUIRE(node->vector_property(mpm::properties::Vector::InternalForce, + Nphase)(i) == + Approx(20.).epsilon(Tolerance)); + + // Assign force to 10 using vector property update false + REQUIRE_NOTHROW(node->update_vector_property( + mpm::properties::Vector::InternalForce, false, Nphase, force)); + for (unsigned i = 0; i < force.size(); ++i) + REQUIRE(node->vector_property(mpm::properties::Vector::InternalForce, + Nphase)(i) == + Approx(10.).epsilon(Tolerance)); } SECTION("Check compute acceleration and velocity") { @@ -308,12 +341,20 @@ TEST_CASE("Node is checked for 1D case", "[node][1D]") { REQUIRE(node->velocity(Nphase)(i) == Approx(velocity(i)).epsilon(Tolerance)); + // Check boolean properties Friction, should be false + REQUIRE(node->boolean_property(mpm::properties::Boolean::Friction) == + false); + // Apply friction constraints REQUIRE(node->assign_friction_constraint(0, 1., 0.5) == true); // Apply friction constraints REQUIRE(node->assign_friction_constraint(-1, 1., 0.5) == false); REQUIRE(node->assign_friction_constraint(3, 1., 0.5) == false); + // Check boolean properties Friction, should be true + REQUIRE(node->boolean_property(mpm::properties::Boolean::Friction) == + true); + // Test acceleration with constraints acceleration[0] = 0.5 * acceleration[0]; for (unsigned i = 0; i < acceleration.size(); ++i) @@ -607,24 +648,31 @@ TEST_CASE("Node is checked for 2D case", "[node][2D]") { REQUIRE(node->pressure(Nphase) == Approx(0.0).epsilon(Tolerance)); double pressure = 1000.7; // Update pressure to 1000.7 - REQUIRE_NOTHROW(node->update_mass_pressure(Nphase, mass * pressure)); + REQUIRE_NOTHROW( + node->update_mass_pressure(false, Nphase, mass * pressure)); + node->compute_pressure(); REQUIRE(node->pressure(Nphase) == Approx(1000.7).epsilon(Tolerance)); // Update pressure to 2001.4 - REQUIRE_NOTHROW(node->update_mass_pressure(Nphase, mass * pressure)); + REQUIRE_NOTHROW( + node->update_mass_pressure(true, Nphase, mass * pressure)); + node->compute_pressure(); REQUIRE(node->pressure(Nphase) == Approx(2001.4).epsilon(Tolerance)); // Assign pressure to 1000 pressure = 1000.; - node->assign_pressure(Nphase, pressure); + REQUIRE_NOTHROW(node->update_pressure(false, Nphase, pressure)); REQUIRE(node->pressure(Nphase) == Approx(1000.0).epsilon(Tolerance)); // Assign mass to 0 mass = 0.; REQUIRE_NOTHROW(node->update_mass(false, Nphase, mass)); REQUIRE(node->mass(Nphase) == Approx(0.0).epsilon(Tolerance)); - // Try to update pressure to 2000, should throw and keep to 1000. + // Try to update pressure, should throw and keep to 1000. pressure = 1000.; - const double pmass = 1.5; - node->assign_pressure(Nphase, pressure); + node->update_mass_pressure(false, Nphase, pressure); + node->compute_pressure(); REQUIRE(node->pressure(Nphase) == Approx(1000.0).epsilon(Tolerance)); + // Update pressure to 2000. + REQUIRE_NOTHROW(node->update_pressure(true, Nphase, pressure)); + REQUIRE(node->pressure(Nphase) == Approx(2000.0).epsilon(Tolerance)); } SECTION("Check volume") { @@ -939,6 +987,10 @@ TEST_CASE("Node is checked for 2D case", "[node][2D]") { // Apply velocity constraints REQUIRE(node->assign_velocity_constraint(0, -12.5) == true); + // Check boolean properties GenericBC, should be false + REQUIRE(node->boolean_property(mpm::properties::Boolean::GenericBC) == + false); + // Apply rotation matrix with Euler angles alpha = 10 deg, beta = 30 deg Eigen::Matrix euler_angles; euler_angles << 10. * M_PI / 180, 30. * M_PI / 180; @@ -947,6 +999,10 @@ TEST_CASE("Node is checked for 2D case", "[node][2D]") { node->assign_rotation_matrix(rotation_matrix); const auto inverse_rotation_matrix = rotation_matrix.inverse(); + // Check boolean properties GenericBC, should be true + REQUIRE(node->boolean_property(mpm::properties::Boolean::GenericBC) == + true); + // Apply inclined velocity constraints node->apply_velocity_constraints(); @@ -976,6 +1032,10 @@ TEST_CASE("Node is checked for 2D case", "[node][2D]") { REQUIRE(node->assign_velocity_constraint(0, -12.5) == true); REQUIRE(node->assign_velocity_constraint(1, 7.5) == true); + // Check boolean properties GenericBC, should be false + REQUIRE(node->boolean_property(mpm::properties::Boolean::GenericBC) == + false); + // Apply rotation matrix with Euler angles alpha = -10 deg, beta = 30 // deg Eigen::Matrix euler_angles; @@ -985,6 +1045,10 @@ TEST_CASE("Node is checked for 2D case", "[node][2D]") { node->assign_rotation_matrix(rotation_matrix); const auto inverse_rotation_matrix = rotation_matrix.inverse(); + // Check boolean properties GenericBC, should be true + REQUIRE(node->boolean_property(mpm::properties::Boolean::GenericBC) == + true); + // Apply inclined velocity constraints node->apply_velocity_constraints(); @@ -1014,11 +1078,19 @@ TEST_CASE("Node is checked for 2D case", "[node][2D]") { } SECTION("Check Cartesian friction constraints") { + // Check boolean properties Friction, should be false + REQUIRE(node->boolean_property(mpm::properties::Boolean::Friction) == + false); + // Apply friction constraints REQUIRE(node->assign_friction_constraint(1, 1, 0.2) == true); // Check out of bounds condition REQUIRE(node->assign_friction_constraint(2, 1, 0.2) == false); + // Check boolean properties Friction, should be true + REQUIRE(node->boolean_property(mpm::properties::Boolean::Friction) == + true); + // Apply friction constraints node->apply_friction_constraints(dt); @@ -1030,9 +1102,21 @@ TEST_CASE("Node is checked for 2D case", "[node][2D]") { } SECTION("Check general friction constraints in 1 direction") { + // Check boolean properties Friction, should be false + REQUIRE(node->boolean_property(mpm::properties::Boolean::Friction) == + false); + // Apply friction constraints REQUIRE(node->assign_friction_constraint(1, 1, 0.2) == true); + // Check boolean properties Friction, should be true + REQUIRE(node->boolean_property(mpm::properties::Boolean::Friction) == + true); + + // Check boolean properties GenericBC, should be false + REQUIRE(node->boolean_property(mpm::properties::Boolean::GenericBC) == + false); + // Apply rotation matrix with Euler angles alpha = 10 deg, beta = 30 deg Eigen::Matrix euler_angles; euler_angles << 10. * M_PI / 180, 30. * M_PI / 180; @@ -1044,6 +1128,10 @@ TEST_CASE("Node is checked for 2D case", "[node][2D]") { // Apply general friction constraints node->apply_friction_constraints(dt); + // Check boolean properties GenericBC, should be true + REQUIRE(node->boolean_property(mpm::properties::Boolean::GenericBC) == + true); + // Check applied constraints on acceleration in the global coordinates acceleration << 4.905579787672637, 4.920772034660430; for (unsigned i = 0; i < Dim; ++i) @@ -1211,24 +1299,31 @@ TEST_CASE("Node is checked for 3D case", "[node][3D]") { REQUIRE(node->pressure(Nphase) == Approx(0.0).epsilon(Tolerance)); double pressure = 1000.7; // Update pressure to 1000.7 - REQUIRE_NOTHROW(node->update_mass_pressure(Nphase, mass * pressure)); + REQUIRE_NOTHROW( + node->update_mass_pressure(false, Nphase, mass * pressure)); + node->compute_pressure(); REQUIRE(node->pressure(Nphase) == Approx(1000.7).epsilon(Tolerance)); // Update pressure to 2001.4 - REQUIRE_NOTHROW(node->update_mass_pressure(Nphase, mass * pressure)); + REQUIRE_NOTHROW( + node->update_mass_pressure(true, Nphase, mass * pressure)); + node->compute_pressure(); REQUIRE(node->pressure(Nphase) == Approx(2001.4).epsilon(Tolerance)); // Assign pressure to 1000 pressure = 1000.; - node->assign_pressure(Nphase, pressure); + REQUIRE_NOTHROW(node->update_pressure(false, Nphase, pressure)); REQUIRE(node->pressure(Nphase) == Approx(1000.0).epsilon(Tolerance)); // Assign mass to 0 mass = 0.; REQUIRE_NOTHROW(node->update_mass(false, Nphase, mass)); REQUIRE(node->mass(Nphase) == Approx(0.0).epsilon(Tolerance)); - // Try to update pressure to 2000, should throw and keep to 1000. + // Try to update pressure, should throw and keep to 1000. pressure = 1000.; - const double pmass = 1.5; - node->assign_pressure(Nphase, pressure); + node->update_mass_pressure(false, Nphase, pressure); + node->compute_pressure(); REQUIRE(node->pressure(Nphase) == Approx(1000.0).epsilon(Tolerance)); + // Update pressure to 2000. + REQUIRE_NOTHROW(node->update_pressure(true, Nphase, pressure)); + REQUIRE(node->pressure(Nphase) == Approx(2000.0).epsilon(Tolerance)); } SECTION("Check external force") { @@ -1516,6 +1611,10 @@ TEST_CASE("Node is checked for 3D case", "[node][3D]") { REQUIRE(node->assign_velocity_constraint(0, 10.5) == true); REQUIRE(node->assign_velocity_constraint(2, -12.5) == true); + // Check boolean properties GenericBC, should be false + REQUIRE(node->boolean_property(mpm::properties::Boolean::GenericBC) == + false); + // Apply rotation matrix with Euler angles alpha = 10 deg, beta = 20 deg // and gamma = 30 deg Eigen::Matrix euler_angles; @@ -1525,6 +1624,10 @@ TEST_CASE("Node is checked for 3D case", "[node][3D]") { node->assign_rotation_matrix(rotation_matrix); const auto inverse_rotation_matrix = rotation_matrix.inverse(); + // Check boolean properties GenericBC, should be true + REQUIRE(node->boolean_property(mpm::properties::Boolean::GenericBC) == + true); + // Apply constraints node->apply_velocity_constraints(); @@ -1560,6 +1663,10 @@ TEST_CASE("Node is checked for 3D case", "[node][3D]") { REQUIRE(node->assign_velocity_constraint(1, -12.5) == true); REQUIRE(node->assign_velocity_constraint(2, 7.5) == true); + // Check boolean properties GenericBC, should be false + REQUIRE(node->boolean_property(mpm::properties::Boolean::GenericBC) == + false); + // Apply rotation matrix with Euler angles alpha = -10 deg, beta = 20 // deg and gamma = -30 deg Eigen::Matrix euler_angles; @@ -1569,6 +1676,10 @@ TEST_CASE("Node is checked for 3D case", "[node][3D]") { node->assign_rotation_matrix(rotation_matrix); const auto inverse_rotation_matrix = rotation_matrix.inverse(); + // Check boolean properties GenericBC, should be true + REQUIRE(node->boolean_property(mpm::properties::Boolean::GenericBC) == + true); + // Apply constraints node->apply_velocity_constraints(); @@ -1602,11 +1713,19 @@ TEST_CASE("Node is checked for 3D case", "[node][3D]") { } SECTION("Check Cartesian friction constraints") { + // Check boolean properties Friction, should be false + REQUIRE(node->boolean_property(mpm::properties::Boolean::Friction) == + false); + // Apply friction constraints REQUIRE(node->assign_friction_constraint(2, 2, 0.3) == true); // Check out of bounds condition REQUIRE(node->assign_friction_constraint(4, 1, 0.2) == false); + // Check boolean properties Friction, should be true + REQUIRE(node->boolean_property(mpm::properties::Boolean::Friction) == + true); + // Apply constraints node->apply_friction_constraints(dt); @@ -1618,9 +1737,21 @@ TEST_CASE("Node is checked for 3D case", "[node][3D]") { } SECTION("Check general friction constraints in 1 direction") { + // Check boolean properties Friction, should be false + REQUIRE(node->boolean_property(mpm::properties::Boolean::Friction) == + false); + // Apply friction constraints REQUIRE(node->assign_friction_constraint(2, 2, 0.3) == true); + // Check boolean properties Friction, should be true + REQUIRE(node->boolean_property(mpm::properties::Boolean::Friction) == + true); + + // Check boolean properties GenericBC, should be false + REQUIRE(node->boolean_property(mpm::properties::Boolean::GenericBC) == + false); + // Apply rotation matrix with Euler angles alpha = 10 deg, beta = 20 deg // and gamma = 30 deg Eigen::Matrix euler_angles; @@ -1633,6 +1764,10 @@ TEST_CASE("Node is checked for 3D case", "[node][3D]") { // Apply inclined velocity constraints node->apply_friction_constraints(dt); + // Check boolean properties GenericBC, should be true + REQUIRE(node->boolean_property(mpm::properties::Boolean::GenericBC) == + true); + // Check applied constraints on acceleration in the global coordinates acceleration << 4.602895052828914, 4.492575657560740, 4.751301246937935; for (unsigned i = 0; i < Dim; ++i) diff --git a/tests/particle_cell_crossing_test.cc b/tests/particle_cell_crossing_test.cc index ea3a1eb35..c8f616e74 100644 --- a/tests/particle_cell_crossing_test.cc +++ b/tests/particle_cell_crossing_test.cc @@ -165,12 +165,14 @@ TEST_CASE("Particle cell crossing is checked for 2D case", material, mpm::ParticlePhase::Solid)); // Compute volume - mesh->iterate_over_particles(std::bind( - &mpm::ParticleBase::compute_volume, std::placeholders::_1)); + mesh->iterate_over_particles([](std::shared_ptr> ptr) { + return mpm::particle::compute_volume(ptr); + }); // Compute mass - mesh->iterate_over_particles( - std::bind(&mpm::ParticleBase::compute_mass, std::placeholders::_1)); + mesh->iterate_over_particles([](std::shared_ptr> ptr) { + return mpm::particle::compute_mass(ptr); + }); // Initialise nodes mesh->iterate_over_nodes( @@ -190,9 +192,9 @@ TEST_CASE("Particle cell crossing is checked for 2D case", &mpm::ParticleBase::compute_shapefn, std::placeholders::_1)); // Assign mass and momentum to nodes - mesh->iterate_over_particles( - std::bind(&mpm::ParticleBase::map_mass_momentum_to_nodes, - std::placeholders::_1)); + mesh->iterate_over_particles([](std::shared_ptr> ptr) { + return mpm::particle::map_mass_momentum_to_nodes(ptr); + }); // Iterate over active nodes to compute acceleratation and velocity mesh->iterate_over_nodes_predicate( @@ -421,12 +423,14 @@ TEST_CASE("Particle cell crossing is checked for 3D case", material, mpm::ParticlePhase::Solid)); // Compute volume - mesh->iterate_over_particles(std::bind( - &mpm::ParticleBase::compute_volume, std::placeholders::_1)); + mesh->iterate_over_particles([](std::shared_ptr> ptr) { + return mpm::particle::compute_volume(ptr); + }); // Compute mass - mesh->iterate_over_particles( - std::bind(&mpm::ParticleBase::compute_mass, std::placeholders::_1)); + mesh->iterate_over_particles([](std::shared_ptr> ptr) { + return mpm::particle::compute_mass(ptr); + }); // Initialise nodes mesh->iterate_over_nodes( @@ -450,9 +454,9 @@ TEST_CASE("Particle cell crossing is checked for 3D case", &mpm::ParticleBase::compute_shapefn, std::placeholders::_1)); // Assign mass and momentum to nodes - mesh->iterate_over_particles( - std::bind(&mpm::ParticleBase::map_mass_momentum_to_nodes, - std::placeholders::_1)); + mesh->iterate_over_particles([](std::shared_ptr> ptr) { + return mpm::particle::map_mass_momentum_to_nodes(ptr); + }); // Iterate over active nodes to compute acceleratation and velocity mesh->iterate_over_nodes_predicate( diff --git a/tests/particle_test.cc b/tests/particle_test.cc index 0d2656f54..b8e25107a 100644 --- a/tests/particle_test.cc +++ b/tests/particle_test.cc @@ -150,6 +150,19 @@ TEST_CASE("Particle is checked for 1D case", "[particle][1D]") { particle->assign_mass(mass); REQUIRE(particle->mass() == Approx(100.5).epsilon(Tolerance)); + // Check mass using scalar properties + REQUIRE(particle->scalar_property(mpm::properties::Scalar::Mass) == + Approx(100.5).epsilon(Tolerance)); + mass = 200.5; + particle->update_scalar_property(mpm::properties::Scalar::Mass, false, + mass); + REQUIRE(particle->scalar_property(mpm::properties::Scalar::Mass) == + Approx(200.5).epsilon(Tolerance)); + + particle->update_scalar_property(mpm::properties::Scalar::Mass, true, mass); + REQUIRE(particle->scalar_property(mpm::properties::Scalar::Mass) == + Approx(401.).epsilon(Tolerance)); + // Check stress Eigen::Matrix stress; for (unsigned i = 0; i < stress.size(); ++i) stress(i) = 17.51; @@ -165,10 +178,29 @@ TEST_CASE("Particle is checked for 1D case", "[particle][1D]") { for (unsigned i = 0; i < velocity.size(); ++i) REQUIRE(particle->velocity()(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(particle->assign_velocity(velocity) == true); + particle->assign_velocity(velocity); for (unsigned i = 0; i < velocity.size(); ++i) REQUIRE(particle->velocity()(i) == Approx(17.51).epsilon(Tolerance)); + // Check velocity using vector properties + for (unsigned i = 0; i < velocity.size(); ++i) + REQUIRE(particle->vector_property(mpm::properties::Vector::Velocity)(i) == + Approx(17.51).epsilon(Tolerance)); + + for (unsigned i = 0; i < velocity.size(); ++i) velocity(i) = 13.88; + + particle->update_vector_property(mpm::properties::Vector::Velocity, false, + velocity); + for (unsigned i = 0; i < velocity.size(); ++i) + REQUIRE(particle->vector_property(mpm::properties::Vector::Velocity)(i) == + Approx(13.88).epsilon(Tolerance)); + + particle->update_vector_property(mpm::properties::Vector::Velocity, true, + velocity); + for (unsigned i = 0; i < velocity.size(); ++i) + REQUIRE(particle->vector_property(mpm::properties::Vector::Velocity)(i) == + Approx(27.76).epsilon(Tolerance)); + // Assign volume REQUIRE(particle->assign_volume(0.0) == false); REQUIRE(particle->assign_volume(-5.0) == false); @@ -202,6 +234,14 @@ TEST_CASE("Particle is checked for 1D case", "[particle][1D]") { else REQUIRE(particle->traction()(i) == Approx(0.).epsilon(Tolerance)); } + + // Check for boolean property assignment and return + REQUIRE(particle->boolean_property(mpm::properties::Boolean::SetTraction) == + true); + particle->assign_boolean_property(mpm::properties::Boolean::SetTraction, + false); + REQUIRE(particle->boolean_property(mpm::properties::Boolean::SetTraction) == + false); } SECTION("Check initialise particle HDF5") { @@ -618,7 +658,8 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { // Add particle mpm::Index id = 0; coords << 0.75, 0.75; - auto particle = std::make_shared>(id, coords); + std::shared_ptr> particle = + std::make_shared>(id, coords); // Time-step const double dt = 0.1; @@ -676,10 +717,7 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { // REQUIRE_NOTHROW(particle->compute_updated_position(dt) == false); // Compute updated particle location from nodal velocity should fail // TODO Assert: REQUIRE_NOTHROW(particle->compute_updated_position(dt, - // true)); Compute volume - // TODO Assert: REQUIRE(particle->compute_volume() == false); - // Update volume should fail - // TODO Assert: REQUIRE(particle->update_volume() == false); + // true)); REQUIRE(particle->assign_cell(cell) == true); REQUIRE(cell->status() == true); @@ -698,7 +736,7 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { // Check volume REQUIRE(particle->volume() == Approx(2.0).epsilon(Tolerance)); // Compute volume - REQUIRE_NOTHROW(particle->compute_volume()); + REQUIRE_NOTHROW(mpm::particle::compute_volume(particle)); // Check volume REQUIRE(particle->volume() == Approx(1.0).epsilon(Tolerance)); @@ -721,9 +759,6 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { Factory, unsigned, const Json&>::instance()->create( "LinearElastic2D", std::move(mid), jmaterial); - // Check compute mass before material and volume - // TODO Assert: REQUIRE(particle->compute_mass() == false); - // Test compute stress before material assignment // TODO Assert: REQUIRE(particle->compute_stress() == false); @@ -734,19 +769,19 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { REQUIRE(particle->material_id() == 1); // Compute volume - REQUIRE_NOTHROW(particle->compute_volume()); + REQUIRE_NOTHROW(mpm::particle::compute_volume(particle)); // Compute mass - REQUIRE_NOTHROW(particle->compute_mass()); + REQUIRE_NOTHROW(mpm::particle::compute_mass(particle)); // Mass REQUIRE(particle->mass() == Approx(1000.).epsilon(Tolerance)); + // Compute mass function + REQUIRE_NOTHROW(mpm::particle::compute_mass(particle)); + REQUIRE(particle->scalar_property(mpm::properties::Scalar::Mass) == + Approx(1000.).epsilon(Tolerance)); // Map particle mass to nodes particle->assign_mass(std::numeric_limits::max()); - // TODO Assert: REQUIRE_NOTHROW(particle->map_mass_momentum_to_nodes()); - - // Map particle pressure to nodes - // TODO Assert: REQUIRE(particle->map_pressure_to_nodes() == false); // Assign mass to nodes REQUIRE(particle->compute_reference_location() == true); @@ -756,15 +791,12 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { Eigen::VectorXd velocity; velocity.resize(Dim); for (unsigned i = 0; i < velocity.size(); ++i) velocity(i) = i; - REQUIRE(particle->assign_velocity(velocity) == true); + particle->assign_velocity(velocity); for (unsigned i = 0; i < velocity.size(); ++i) REQUIRE(particle->velocity()(i) == Approx(i).epsilon(Tolerance)); - REQUIRE_NOTHROW(particle->compute_mass()); - REQUIRE_NOTHROW(particle->map_mass_momentum_to_nodes()); - - // TODO Assert: REQUIRE(particle->map_pressure_to_nodes() == false); - REQUIRE(particle->compute_pressure_smoothing() == false); + REQUIRE_NOTHROW(mpm::particle::compute_mass(particle)); + REQUIRE_NOTHROW(mpm::particle::map_mass_momentum_to_nodes(particle)); // Values of nodal mass std::array nodal_mass{562.5, 187.5, 62.5, 187.5}; @@ -854,7 +886,7 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { // Update volume strain rate REQUIRE(particle->volume() == Approx(1.0).epsilon(Tolerance)); particle->compute_strain(dt); - REQUIRE_NOTHROW(particle->update_volume()); + REQUIRE_NOTHROW(mpm::particle::update_volume(particle)); REQUIRE(particle->volume() == Approx(1.2).epsilon(Tolerance)); // Compute stress @@ -877,7 +909,7 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { Eigen::Matrix gravity; gravity << 0., -9.81; - particle->map_body_force(gravity); + mpm::particle::map_body_force(particle, gravity); // Body force Eigen::Matrix body_force; @@ -906,7 +938,7 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { // Assign traction to particle particle->assign_traction(direction, mfunction->value(current_time) * traction); - particle->map_traction_force(); + mpm::particle::map_traction_force(particle); // Traction force Eigen::Matrix traction_force; @@ -929,7 +961,7 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { particle->assign_traction(direction, -traction * mfunction->value(current_time)); // Map traction force - particle->map_traction_force(); + mpm::particle::map_traction_force(particle); // Check nodal external force for (unsigned i = 0; i < traction_force.rows(); ++i) for (unsigned j = 0; j < traction_force.cols(); ++j) @@ -1049,23 +1081,22 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { Factory, unsigned, const Json&>::instance() ->create("Newtonian2D", std::move(mid1), jmaterial1); + // Reset nodal properties + for (const auto& node : nodes) node->initialise(); + // Assign material properties REQUIRE(particle->assign_material(material1) == true); // Compute volume - REQUIRE_NOTHROW(particle->compute_volume()); + REQUIRE_NOTHROW(mpm::particle::compute_volume(particle)); // Compute mass - REQUIRE_NOTHROW(particle->compute_mass()); + REQUIRE_NOTHROW(mpm::particle::compute_mass(particle)); // Mass REQUIRE(particle->mass() == Approx(1000.).epsilon(Tolerance)); // Map particle mass to nodes particle->assign_mass(std::numeric_limits::max()); - // TODO Assert: REQUIRE(particle->map_mass_momentum_to_nodes() == false); - - // Map particle pressure to nodes - // TODO Assert: REQUIRE(particle->map_pressure_to_nodes() == false); // Assign mass to nodes REQUIRE(particle->compute_reference_location() == true); @@ -1074,12 +1105,12 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { // Check velocity velocity.resize(Dim); for (unsigned i = 0; i < velocity.size(); ++i) velocity(i) = i; - REQUIRE(particle->assign_velocity(velocity) == true); + particle->assign_velocity(velocity); for (unsigned i = 0; i < velocity.size(); ++i) REQUIRE(particle->velocity()(i) == Approx(i).epsilon(Tolerance)); - REQUIRE_NOTHROW(particle->compute_mass()); - REQUIRE_NOTHROW(particle->map_mass_momentum_to_nodes()); + REQUIRE_NOTHROW(mpm::particle::compute_mass(particle)); + REQUIRE_NOTHROW(mpm::particle::map_mass_momentum_to_nodes(particle)); // Check volumetric strain at centroid volumetric_strain = 0.2; @@ -1093,8 +1124,38 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { particle->pressure() == Approx(-8333333.333333333 * volumetric_strain).epsilon(Tolerance)); - REQUIRE_NOTHROW(particle->map_pressure_to_nodes()); - REQUIRE(particle->compute_pressure_smoothing() == true); + // Check return and assign state variable + REQUIRE( + particle->state_variable("pressure") == + Approx(-8333333.333333333 * volumetric_strain).epsilon(Tolerance)); + + REQUIRE_NOTHROW( + particle->assign_pressure(-8333333.333333333 * volumetric_strain)); + + // Check pressure smoothing + REQUIRE_NOTHROW(mpm::particle::map_mass_pressure_to_nodes(particle)); + for (const auto& node : nodes) node->compute_pressure(); + REQUIRE_NOTHROW(mpm::particle::compute_pressure_smoothing(particle)); + REQUIRE( + particle->pressure() == + Approx(-8333333.333333333 * volumetric_strain).epsilon(Tolerance)); + + // Assign pressure equal to 1.0, we expect its smoothed value is also 1.0 + REQUIRE_NOTHROW(particle->assign_state_variable("pressure", 1.0)); + REQUIRE(particle->pressure() == Approx(1.0).epsilon(Tolerance)); + + for (const auto& node : nodes) + node->update_scalar_property(mpm::properties::Scalar::MassPressure, + false, 0, 0.0); + REQUIRE_NOTHROW(mpm::particle::map_mass_pressure_to_nodes(particle)); + for (const auto& node : nodes) { + REQUIRE( + node->scalar_property(mpm::properties::Scalar::MassPressure, 0) == + node->scalar_property(mpm::properties::Scalar::Mass, 0)); + node->compute_pressure(); + } + REQUIRE_NOTHROW(mpm::particle::compute_pressure_smoothing(particle)); + REQUIRE(particle->pressure() == Approx(1.0).epsilon(Tolerance)); } SECTION("Particle assign state variables") { @@ -1228,6 +1289,19 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { particle->assign_mass(mass); REQUIRE(particle->mass() == Approx(100.5).epsilon(Tolerance)); + // Check mass using scalar properties + REQUIRE(particle->scalar_property(mpm::properties::Scalar::Mass) == + Approx(100.5).epsilon(Tolerance)); + mass = 111.11; + particle->update_scalar_property(mpm::properties::Scalar::Mass, false, + mass); + REQUIRE(particle->scalar_property(mpm::properties::Scalar::Mass) == + Approx(111.11).epsilon(Tolerance)); + + particle->update_scalar_property(mpm::properties::Scalar::Mass, true, mass); + REQUIRE(particle->scalar_property(mpm::properties::Scalar::Mass) == + Approx(222.22).epsilon(Tolerance)); + // Check stress Eigen::Matrix stress; for (unsigned i = 0; i < stress.size(); ++i) stress(i) = 17.52; @@ -1243,10 +1317,29 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { for (unsigned i = 0; i < velocity.size(); ++i) REQUIRE(particle->velocity()(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(particle->assign_velocity(velocity) == true); + particle->assign_velocity(velocity); for (unsigned i = 0; i < velocity.size(); ++i) REQUIRE(particle->velocity()(i) == Approx(19.745).epsilon(Tolerance)); + // Check velocity using vector properties + for (unsigned i = 0; i < velocity.size(); ++i) + REQUIRE(particle->vector_property(mpm::properties::Vector::Velocity)(i) == + Approx(19.745).epsilon(Tolerance)); + + for (unsigned i = 0; i < velocity.size(); ++i) velocity(i) = 11.22; + + particle->update_vector_property(mpm::properties::Vector::Velocity, false, + velocity); + for (unsigned i = 0; i < velocity.size(); ++i) + REQUIRE(particle->vector_property(mpm::properties::Vector::Velocity)(i) == + Approx(11.22).epsilon(Tolerance)); + + particle->update_vector_property(mpm::properties::Vector::Velocity, true, + velocity); + for (unsigned i = 0; i < velocity.size(); ++i) + REQUIRE(particle->vector_property(mpm::properties::Vector::Velocity)(i) == + Approx(22.44).epsilon(Tolerance)); + // Assign volume REQUIRE(particle->assign_volume(0.0) == false); REQUIRE(particle->assign_volume(-5.0) == false); @@ -1283,6 +1376,14 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { else REQUIRE(particle->traction()(i) == Approx(0.).epsilon(Tolerance)); } + + // Check for boolean property assignment and return + REQUIRE(particle->boolean_property(mpm::properties::Boolean::SetTraction) == + true); + particle->assign_boolean_property(mpm::properties::Boolean::SetTraction, + false); + REQUIRE(particle->boolean_property(mpm::properties::Boolean::SetTraction) == + false); } // Check initialise particle from HDF5 file @@ -1934,11 +2035,6 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") { // TODO Assert: REQUIRE_NOTHROW(particle->compute_updated_position(dt, // true)); - // Compute volume - // TODO Assert: REQUIRE(particle->compute_volume() == false); - // Update volume should fail - // TODO Assert: REQUIRE(particle->update_volume() == false); - REQUIRE(particle->assign_cell(cell) == true); REQUIRE(cell->status() == true); REQUIRE(particle->cell_id() == 10); @@ -1956,7 +2052,7 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") { // Check volume REQUIRE(particle->volume() == Approx(2.0).epsilon(Tolerance)); // Compute volume - REQUIRE_NOTHROW(particle->compute_volume()); + REQUIRE_NOTHROW(mpm::particle::compute_volume(particle)); // Check volume REQUIRE(particle->volume() == Approx(8.0).epsilon(Tolerance)); @@ -1979,9 +2075,6 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") { Factory, unsigned, const Json&>::instance()->create( "LinearElastic3D", std::move(mid), jmaterial); - // Check compute mass before material and volume - // TODO Assert: REQUIRE(particle->compute_mass() == false); - // Test compute stress before material assignment // TODO Assert: REQUIRE(particle->compute_stress() == false); @@ -1992,19 +2085,15 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") { REQUIRE(particle->material_id() == 0); // Compute volume - REQUIRE_NOTHROW(particle->compute_volume()); + REQUIRE_NOTHROW(mpm::particle::compute_volume(particle)); // Compute mass - REQUIRE_NOTHROW(particle->compute_mass()); + REQUIRE_NOTHROW(mpm::particle::compute_mass(particle)); // Mass REQUIRE(particle->mass() == Approx(8000.).epsilon(Tolerance)); // Map particle mass to nodes particle->assign_mass(std::numeric_limits::max()); - // TODO Assert: REQUIRE(particle->map_mass_momentum_to_nodes() == false); - - // Map particle pressure to nodes - // TODO Assert: REQUIRE(particle->map_pressure_to_nodes() == false); // Assign mass to nodes REQUIRE(particle->compute_reference_location() == true); @@ -2014,15 +2103,12 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") { Eigen::VectorXd velocity; velocity.resize(Dim); for (unsigned i = 0; i < velocity.size(); ++i) velocity(i) = i; - REQUIRE(particle->assign_velocity(velocity) == true); + particle->assign_velocity(velocity); for (unsigned i = 0; i < velocity.size(); ++i) REQUIRE(particle->velocity()(i) == Approx(i).epsilon(Tolerance)); - REQUIRE_NOTHROW(particle->compute_mass()); - REQUIRE_NOTHROW(particle->map_mass_momentum_to_nodes()); - - // TODO Assert: REQUIRE(particle->map_pressure_to_nodes() == false); - REQUIRE(particle->compute_pressure_smoothing() == false); + REQUIRE_NOTHROW(mpm::particle::compute_mass(particle)); + REQUIRE_NOTHROW(mpm::particle::map_mass_momentum_to_nodes(particle)); // Values of nodal mass std::array nodal_mass{125., 375., 1125., 375., @@ -2130,7 +2216,7 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") { // Update volume strain rate REQUIRE(particle->volume() == Approx(8.0).epsilon(Tolerance)); particle->compute_strain(dt); - REQUIRE_NOTHROW(particle->update_volume()); + REQUIRE_NOTHROW(mpm::particle::update_volume(particle)); REQUIRE(particle->volume() == Approx(12.0).epsilon(Tolerance)); // Compute stress @@ -2153,7 +2239,7 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") { Eigen::Matrix gravity; gravity << 0., 0., -9.81; - particle->map_body_force(gravity); + mpm::particle::map_body_force(particle, gravity); // Body force Eigen::Matrix body_force; @@ -2185,7 +2271,7 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") { particle->assign_traction(direction, mfunction->value(current_time) * traction); // Map traction force - particle->map_traction_force(); + mpm::particle::map_traction_force(particle); // Traction force Eigen::Matrix traction_force; @@ -2212,7 +2298,7 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") { particle->assign_traction(direction, mfunction->value(current_time) * -traction); // Map traction force - particle->map_traction_force(); + mpm::particle::map_traction_force(particle); // Check nodal external force for (unsigned i = 0; i < traction_force.rows(); ++i) for (unsigned j = 0; j < traction_force.cols(); ++j) @@ -2305,19 +2391,15 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") { REQUIRE(particle->assign_material(material1) == true); // Compute volume - REQUIRE_NOTHROW(particle->compute_volume()); + REQUIRE_NOTHROW(mpm::particle::compute_volume(particle)); // Compute mass - REQUIRE_NOTHROW(particle->compute_mass()); + REQUIRE_NOTHROW(mpm::particle::compute_mass(particle)); // Mass REQUIRE(particle->mass() == Approx(8000.).epsilon(Tolerance)); // Map particle mass to nodes particle->assign_mass(std::numeric_limits::max()); - // TODO Assert: REQUIRE(particle->map_mass_momentum_to_nodes() == false); - - // Map particle pressure to nodes - // TODO Assert: REQUIRE(particle->map_pressure_to_nodes() == false); // Assign mass to nodes REQUIRE(particle->compute_reference_location() == true); @@ -2326,12 +2408,12 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") { // Check velocity velocity.resize(Dim); for (unsigned i = 0; i < velocity.size(); ++i) velocity(i) = i; - REQUIRE(particle->assign_velocity(velocity) == true); + particle->assign_velocity(velocity); for (unsigned i = 0; i < velocity.size(); ++i) REQUIRE(particle->velocity()(i) == Approx(i).epsilon(Tolerance)); - REQUIRE_NOTHROW(particle->compute_mass()); - REQUIRE_NOTHROW(particle->map_mass_momentum_to_nodes()); + REQUIRE_NOTHROW(mpm::particle::compute_mass(particle)); + REQUIRE_NOTHROW(mpm::particle::map_mass_momentum_to_nodes(particle)); // Check volumetric strain at centroid volumetric_strain = 0.5; @@ -2345,8 +2427,33 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") { particle->pressure() == Approx(-8333333.333333333 * volumetric_strain).epsilon(Tolerance)); - REQUIRE_NOTHROW(particle->map_pressure_to_nodes()); - REQUIRE(particle->compute_pressure_smoothing() == true); + // Check return and assign state variable + REQUIRE( + particle->state_variable("pressure") == + Approx(-8333333.333333333 * volumetric_strain).epsilon(Tolerance)); + + REQUIRE_NOTHROW( + particle->assign_pressure(-8333333.333333333 * volumetric_strain)); + + // Check pressure smoothing + REQUIRE_NOTHROW(mpm::particle::map_mass_pressure_to_nodes(particle)); + for (const auto& node : nodes) node->compute_pressure(); + REQUIRE_NOTHROW(mpm::particle::compute_pressure_smoothing(particle)); + REQUIRE(particle->pressure() == + Approx(-2083333.3333333333).epsilon(Tolerance)); + + // Assign pressure equal to 1.0, we expect its smoothed value is 0.5 as + // the nodal mass is mapped twice in this case + REQUIRE_NOTHROW(particle->assign_state_variable("pressure", 1.0)); + REQUIRE(particle->pressure() == Approx(1.0).epsilon(Tolerance)); + + for (const auto& node : nodes) + node->update_scalar_property(mpm::properties::Scalar::MassPressure, + false, 0, 0.0); + REQUIRE_NOTHROW(mpm::particle::map_mass_pressure_to_nodes(particle)); + for (const auto& node : nodes) node->compute_pressure(); + REQUIRE_NOTHROW(mpm::particle::compute_pressure_smoothing(particle)); + REQUIRE(particle->pressure() == Approx(0.5).epsilon(Tolerance)); } SECTION("Particle assign state variables") { @@ -2521,6 +2628,19 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") { particle->assign_mass(mass); REQUIRE(particle->mass() == Approx(100.5).epsilon(Tolerance)); + // Check mass using scalar properties + REQUIRE(particle->scalar_property(mpm::properties::Scalar::Mass) == + Approx(100.5).epsilon(Tolerance)); + mass = 111.11; + particle->update_scalar_property(mpm::properties::Scalar::Mass, false, + mass); + REQUIRE(particle->scalar_property(mpm::properties::Scalar::Mass) == + Approx(111.11).epsilon(Tolerance)); + + particle->update_scalar_property(mpm::properties::Scalar::Mass, true, mass); + REQUIRE(particle->scalar_property(mpm::properties::Scalar::Mass) == + Approx(222.22).epsilon(Tolerance)); + // Check stress Eigen::Matrix stress; for (unsigned i = 0; i < stress.size(); ++i) stress(i) = 1.; @@ -2536,10 +2656,29 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") { for (unsigned i = 0; i < velocity.size(); ++i) REQUIRE(particle->velocity()(i) == Approx(0.).epsilon(Tolerance)); - REQUIRE(particle->assign_velocity(velocity) == true); + particle->assign_velocity(velocity); for (unsigned i = 0; i < velocity.size(); ++i) REQUIRE(particle->velocity()(i) == Approx(17.51).epsilon(Tolerance)); + // Check velocity using vector properties + for (unsigned i = 0; i < velocity.size(); ++i) + REQUIRE(particle->vector_property(mpm::properties::Vector::Velocity)(i) == + Approx(17.51).epsilon(Tolerance)); + + for (unsigned i = 0; i < velocity.size(); ++i) velocity(i) = 11.22; + + particle->update_vector_property(mpm::properties::Vector::Velocity, false, + velocity); + for (unsigned i = 0; i < velocity.size(); ++i) + REQUIRE(particle->vector_property(mpm::properties::Vector::Velocity)(i) == + Approx(11.22).epsilon(Tolerance)); + + particle->update_vector_property(mpm::properties::Vector::Velocity, true, + velocity); + for (unsigned i = 0; i < velocity.size(); ++i) + REQUIRE(particle->vector_property(mpm::properties::Vector::Velocity)(i) == + Approx(22.44).epsilon(Tolerance)); + // Assign volume REQUIRE(particle->assign_volume(0.0) == false); REQUIRE(particle->assign_volume(-5.0) == false);