diff --git a/src/common/field.hh b/src/common/field.hh index d0833d0..97b9b11 100644 --- a/src/common/field.hh +++ b/src/common/field.hh @@ -1,389 +1,408 @@ /** * file field.hh * * @author Till Junge * * @date 07 Sep 2017 * * @brief header-only implementation of a field for field collections * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef FIELD_H #define FIELD_H #include #include #include #include #include #include #include #include #include namespace muSpectre { namespace internal { /* ---------------------------------------------------------------------- */ template class FieldBase { protected: //! constructor //! unique name (whithin Collection) //! number of components //! collection to which this field belongs (eg, material, system) FieldBase(std::string unique_name, size_t nb_components, FieldCollection & collection); public: //! Copy constructor FieldBase(const FieldBase &other) = delete; //! Move constructor FieldBase(FieldBase &&other) noexcept = delete; //! Destructor virtual ~FieldBase() noexcept = default; //! Copy assignment operator FieldBase& operator=(const FieldBase &other) = delete; //! Move assignment operator FieldBase& operator=(FieldBase &&other) noexcept = delete; /* ---------------------------------------------------------------------- */ //!Identifying accessors //! return field name inline const std::string & get_name() const; //! return field type //inline const Field_t & get_type() const; //! return my collection (for iterating) inline const FieldCollection & get_collection() const; //! return my collection (for iterating) inline const size_t & get_nb_components() const; //! return type_id of stored type virtual const std::type_info & get_stored_typeid() const = 0; - inline size_t size() const; + virtual size_t size() const = 0; //! initialise field to zero (do more complicated initialisations through //! fully typed maps) virtual void set_zero() = 0; //! give access to collections friend FieldCollection; protected: /* ---------------------------------------------------------------------- */ //! allocate memory etc - virtual void initialise(size_t size) = 0; + virtual void resize(size_t size) = 0; const std::string name; const size_t nb_components; const FieldCollection & collection; private: }; /* ---------------------------------------------------------------------- */ //! declaraton for friending template class FieldMap; /* ---------------------------------------------------------------------- */ template class TypedFieldBase: public FieldBase { friend class FieldMap; public: - using parent = FieldBase; - using base = parent; + using Parent = FieldBase; + using Base = Parent; //using storage_type = Eigen::Array; - using stored_type = Eigen::Array; - using storage_type = std::vector>; + using StoredType = Eigen::Array; + using StorageType = std::vector>; TypedFieldBase(std::string unique_name, FieldCollection& collection); virtual ~TypedFieldBase() = default; //! return type_id of stored type virtual const std::type_info & get_stored_typeid() const; //! initialise field to zero (do more complicated initialisations through //! fully typed maps) inline void set_zero() override final; + inline void push_back(const StoredType & value); + template< class... Args> + inline void emplace_back(Args&&... args); + size_t size() const override final; protected: inline T* get_ptr_to_entry(const size_t&& index); inline T& get_ref_to_entry(const size_t&& index); - inline virtual void initialise(size_t size) override final; - storage_type array; + inline virtual void resize(size_t size) override final; + StorageType array; }; } // internal /* ---------------------------------------------------------------------- */ template class TensorField: public internal::TypedFieldBase { public: - using parent = internal::TypedFieldBase; - using base = typename parent::base; + using Base = typename Parent::Base; using Field_p = std::unique_ptr>; using component_type = T; //! Copy constructor TensorField(const TensorField &other) = delete; //! Move constructor TensorField(TensorField &&other) noexcept = delete; //! Destructor virtual ~TensorField() noexcept = default; //! Copy assignment operator TensorField& operator=(const TensorField &other) = delete; //! Move assignment operator TensorField& operator=(TensorField &&other) noexcept = delete; //! accessors inline Dim_t get_order() const; inline Dim_t get_dim() const; //! factory function template - friend typename FieldType::base& make_field(std::string unique_name, + friend typename FieldType::Base& make_field(std::string unique_name, CollectionType & collection, Args&&... args); protected: //! constructor protected! TensorField(std::string unique_name, FieldCollection & collection); private: }; /* ---------------------------------------------------------------------- */ template class MatrixField: public internal::TypedFieldBase { public: - using parent = internal::TypedFieldBase; - using base = typename parent::base; + using Base = typename Parent::Base; using Field_p = std::unique_ptr>; using component_type = T; //! Copy constructor MatrixField(const MatrixField &other) = delete; //! Move constructor MatrixField(MatrixField &&other) noexcept = delete; //! Destructor virtual ~MatrixField() noexcept = default; //! Copy assignment operator MatrixField& operator=(const MatrixField &other) = delete; //! Move assignment operator MatrixField& operator=(MatrixField &&other) noexcept = delete; //! accessors inline Dim_t get_nb_row() const; inline Dim_t get_nb_col() const; //! factory function template - friend typename FieldType::base& make_field(std::string unique_name, + friend typename FieldType::Base& make_field(std::string unique_name, CollectionType & collection, Args&&... args); protected: //! constructor protected! MatrixField(std::string unique_name, FieldCollection & collection); private: }; /* ---------------------------------------------------------------------- */ //! convenience alias template using ScalarField = TensorField; /* ---------------------------------------------------------------------- */ // Implementations namespace internal { /* ---------------------------------------------------------------------- */ template FieldBase::FieldBase(std::string unique_name, size_t nb_components_, FieldCollection & collection_) :name(unique_name), nb_components(nb_components_), collection(collection_) {} /* ---------------------------------------------------------------------- */ template inline const std::string & FieldBase::get_name() const { return this->name; } /* ---------------------------------------------------------------------- */ template inline const FieldCollection & FieldBase:: get_collection() const { return this->collection; } /* ---------------------------------------------------------------------- */ template inline const size_t & FieldBase:: get_nb_components() const { return this->nb_components; } - /* ---------------------------------------------------------------------- */ - template - inline size_t FieldBase:: - size() const { - return this->collection.size(); - } - /* ---------------------------------------------------------------------- */ template TypedFieldBase:: TypedFieldBase(std::string unique_name, FieldCollection & collection) :FieldBase(unique_name, NbComponents, collection){ static_assert ((std::is_arithmetic::value || std::is_same::value), "Use TypedFieldBase for integer, real or complex scalars for T"); static_assert(NbComponents > 0, "Only fields with more than 0 components"); } /* ---------------------------------------------------------------------- */ //! return type_id of stored type template const std::type_info & TypedFieldBase:: get_stored_typeid() const { return typeid(T); } /* ---------------------------------------------------------------------- */ template void TypedFieldBase:: set_zero() { std::fill(this->array.begin(), this->array.end(), T{}); } + /* ---------------------------------------------------------------------- */ + template + size_t TypedFieldBase:: + size() const { + return this->array.size(); + } + /* ---------------------------------------------------------------------- */ template T* TypedFieldBase:: get_ptr_to_entry(const size_t&& index) { return &this->array[std::move(index)](0, 0); } /* ---------------------------------------------------------------------- */ template T& TypedFieldBase:: get_ref_to_entry(const size_t && index) { return this->array[std::move(index)](0, 0); } /* ---------------------------------------------------------------------- */ template void TypedFieldBase:: - initialise(size_t size) { + resize(size_t size) { this->array.resize(size); } + /* ---------------------------------------------------------------------- */ + template + void TypedFieldBase:: + push_back(const StoredType & value) { + this->array.push_back(value); + } + + /* ---------------------------------------------------------------------- */ + template + template + void TypedFieldBase:: + emplace_back(Args&&... args) { + this->array.emplace_back(std::move(args...)); + } + } // internal /* ---------------------------------------------------------------------- */ //! Factory function, guarantees that only fields get created that are //! properly registered and linked to a collection. template - typename FieldType::base & + typename FieldType::Base & make_field(std::string unique_name, FieldCollection & collection, Args&&... args) { auto && ptr = std::unique_ptr{ new FieldType(unique_name, collection, args...)}; auto& retref = *ptr; collection.register_field(std::move(ptr)); return retref; } /* ---------------------------------------------------------------------- */ template TensorField:: TensorField(std::string unique_name, FieldCollection & collection) - :parent(unique_name, collection) {} + :Parent(unique_name, collection) {} /* ---------------------------------------------------------------------- */ template Dim_t TensorField:: get_order() const { return order; } /* ---------------------------------------------------------------------- */ template Dim_t TensorField:: get_dim() const { return dim; } /* ---------------------------------------------------------------------- */ template MatrixField:: MatrixField(std::string unique_name, FieldCollection & collection) - :parent(unique_name, collection) {} + :Parent(unique_name, collection) {} /* ---------------------------------------------------------------------- */ template Dim_t MatrixField:: get_nb_col() const { return NbCol; } /* ---------------------------------------------------------------------- */ template Dim_t MatrixField:: get_nb_row() const { return NbRow; } } // muSpectre #endif /* FIELD_H */ diff --git a/src/common/field_collection.hh b/src/common/field_collection.hh index 18f243f..c002656 100644 --- a/src/common/field_collection.hh +++ b/src/common/field_collection.hh @@ -1,277 +1,307 @@ /** * file field_collection.hh * * @author Till Junge * * @date 07 Sep 2017 * * @brief Provides pixel-iterable containers for scalar and tensorial fields, * addressable by field name * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef FIELD_COLLECTION_H #define FIELD_COLLECTION_H #include "common/common.hh" #include "common/ccoord_operations.hh" #include "common/field.hh" #include #include #include #include #include #include #include #include #include namespace muSpectre { class FieldCollectionError: public std::runtime_error { public: explicit FieldCollectionError(const std::string& what) :std::runtime_error(what){} explicit FieldCollectionError(const char * what) :std::runtime_error(what){} }; class FieldError: public FieldCollectionError { using Parent = FieldCollectionError; public: explicit FieldError(const std::string& what) :Parent(what){} explicit FieldError(const char * what) :Parent(what){} }; class FieldInterpretationError: public FieldError { public: explicit FieldInterpretationError(const std::string & what) :FieldError(what){} explicit FieldInterpretationError(const char * what) :FieldError(what){} }; /* ---------------------------------------------------------------------- */ //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) //! Global determines whether this field is present everywhere or per material template class FieldCollection { public: using Field = internal::FieldBase; using Field_p = std::unique_ptr; using Ccoord = Ccoord_t; //! Default constructor FieldCollection(); //! Copy constructor FieldCollection(const FieldCollection &other) = delete; //! Move constructor FieldCollection(FieldCollection &&other) noexcept = delete; //! Destructor virtual ~FieldCollection() noexcept = default; //! Copy assignment operator FieldCollection& operator=(const FieldCollection &other) = delete; //! Move assignment operator FieldCollection& operator=(FieldCollection &&other) noexcept = delete; //! add a pixel/voxel to the field collection template inline std::enable_if_t add_pixel(const Ccoord & local_ccoord); - //! allocate memory, etc + //! allocate memory, etc. At this point, the field knows how many + //! entries it should have (global fields know from the 'sizes' + //! parameter and non-global fields know from the size of the + //! pixel container, which was grown using 'add_pixel(...)'. The + //! job of initialise is to check that all fields are either of + //! size 0, in which case they need to be allocated, or are of the + //! same size as the pixel_container. Any field of a different + //! size is wrong. TODO: check whether it makes sense to put a + //! runtime check here template inline std::enable_if_t initialise(); template inline std::enable_if_t initialise(Ccoord sizes); //! Register a new field (fields need to be in heap, so I want to keep them //! yas shared pointers void register_field(Field_p&& field); //! for return values of iterators constexpr inline static Dim_t spatial_dim(); //! for return values of iterators inline Dim_t get_spatial_dim() const; //! retrieve field by unique_name inline Field& operator[](std::string unique_name); //! retrieve field by unique_name with bounds checking inline Field& at(std::string unique_name); //! returns size of collection, this refers to the number of pixels handled //! by the collection, not the number of fields inline size_t size() const {return this->size_;} //! return the pixel sizes template inline const std::enable_if_t & get_sizes() const; //! returns the linear index corresponding to cell coordinates inline size_t get_index(Ccoord && ccoord) const; //! returns the cell coordinates corresponding to a linear index inline Ccoord get_ccoord(size_t index) const; protected: std::map fields; bool is_initialised = false; const Uint id; static Uint counter; size_t size_ = 0; Ccoord sizes; //! container of pixel coords for non-global collections std::vector ccoords; //! container of indices for non-global collections (slow!) std::map indices; private: }; /* ---------------------------------------------------------------------- */ template Uint FieldCollection::counter{0}; /* ---------------------------------------------------------------------- */ template FieldCollection::FieldCollection() :id(counter++){} /* ---------------------------------------------------------------------- */ template template std::enable_if_t FieldCollection:: add_pixel(const Ccoord & local_ccoord) { this->indices[local_ccoord] = this->ccoords.size(); this->ccoords.push_back(local_ccoord); } /* ---------------------------------------------------------------------- */ template template std::enable_if_t FieldCollection:: initialise() { this->size_ = this->ccoords.size(); std::for_each(std::begin(this->fields), std::end(this->fields), [this](auto && item) { - item.second->initialise(this->size()); + auto && field = *item.second; + const auto field_size = field.size(); + if (field_size == 0) { + field.resize(this->size()); + } else if (field_size != this->size()) { + std::stringstream err_stream; + err_stream << "Field '" << field.get_name() + << "' contains " << field_size + << " entries, but the field collection " + << "has " << this->size() << " pixels"; + throw FieldCollectionError(err_stream.str()); + } }); this->is_initialised = true; } /* ---------------------------------------------------------------------- */ template template std::enable_if_t FieldCollection:: initialise(Ccoord sizes) { this->size_ = std::accumulate(sizes.begin(), sizes.end(), 1, std::multiplies()); this->sizes = sizes; std::for_each(std::begin(this->fields), std::end(this->fields), [this](auto && item) { - item.second->initialise(this->size()); + auto && field = *item.second; + const auto field_size = field.size(); + if (field_size == 0) { + field.resize(this->size()); + } else if (field_size != this->size()) { + std::stringstream err_stream; + err_stream << "Field '" << field.get_name() + << "' contains " << field_size + << " entries, but the field collection " + << "has " << this->size() << " pixels"; + throw FieldCollectionError(err_stream.str()); + } }); this->is_initialised = true; } /* ---------------------------------------------------------------------- */ template void FieldCollection::register_field(Field_p &&field) { auto&& search_it = this->fields.find(field->get_name()); auto&& does_exist = search_it != this->fields.end(); if (does_exist) { std::stringstream err_str; err_str << "a field named " << field->get_name() << "is already registered in this field collection. " << "Currently registered fields: "; for (const auto& name_field_pair: this->fields) { err_str << ", " << name_field_pair.first; } throw FieldCollectionError(err_str.str()); } this->fields[field->get_name()] = std::move(field); } /* ---------------------------------------------------------------------- */ template constexpr Dim_t FieldCollection::spatial_dim() { return DimS; } /* ---------------------------------------------------------------------- */ template Dim_t FieldCollection::get_spatial_dim() const { return DimS; } /* ---------------------------------------------------------------------- */ template typename FieldCollection::Field& FieldCollection::operator[](std::string unique_name) { return *(this->fields[unique_name]); } /* ---------------------------------------------------------------------- */ template typename FieldCollection::Field& FieldCollection::at(std::string unique_name) { return *(this->fields.at(unique_name)); } //----------------------------------------------------------------------------// //! return the pixel sizes template template const std::enable_if_t>& FieldCollection::get_sizes() const { return this->sizes; } //----------------------------------------------------------------------------// //! returns the linear index corresponding to cell coordinates template size_t FieldCollection::get_index(Ccoord && ccoord) const { return CcoordOps::get_index(this->sizes, std::move(ccoord)); } //----------------------------------------------------------------------------// //! returns the cell coordinates corresponding to a linear index template typename FieldCollection::Ccoord FieldCollection::get_ccoord(size_t index) const { return CcoordOps::get_ccoord(this->sizes, std::move(index)); } } // muSpectre #endif /* FIELD_COLLECTION_H */ diff --git a/src/common/field_map_base.hh b/src/common/field_map_base.hh index f83a579..391ac26 100644 --- a/src/common/field_map_base.hh +++ b/src/common/field_map_base.hh @@ -1,499 +1,499 @@ /** * file field_map.hh * * @author Till Junge * * @date 12 Sep 2017 * * @brief Defined a strongly defines proxy that iterates efficiently over a field * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef FIELD_MAP_H #define FIELD_MAP_H #include #include #include #include #include #include "common/field.hh" #include "common/common.hh" namespace muSpectre { namespace internal { //----------------------------------------------------------------------------// //! little helper to automate creation of const maps without duplication template struct const_corrector { using type = typename T::reference; }; template struct const_corrector { using type = typename T::const_reference; }; template using const_corrector_t = typename const_corrector::type; //----------------------------------------------------------------------------// template class FieldMap { public: using TypedField = ::muSpectre::internal::TypedFieldBase ; - using Field = typename TypedField::parent; + using Field = typename TypedField::Parent; using size_type = std::size_t; //! Default constructor FieldMap() = delete; FieldMap(Field & field); //! Copy constructor FieldMap(const FieldMap &other) = delete; //! Move constructor FieldMap(FieldMap &&other) noexcept = delete; //! Destructor virtual ~FieldMap() noexcept = default; //! Copy assignment operator FieldMap& operator=(const FieldMap &other) = delete; //! Move assignment operator FieldMap& operator=(FieldMap &&other) noexcept = delete; //! give static access to nb_components inline constexpr static size_t nb_components(); //! give human-readable field map type virtual std::string info_string() const = 0; //! return field name inline const std::string & get_name() const; //! return my collection (for iterating) inline const FieldCollection & get_collection() const; //! member access needs to be implemented by inheriting classes //inline value_type operator[](size_t index); //inline value_type operator[](Ccoord ccord); //! check compatibility (must be called by all inheriting classes at the //! end of their constructors inline void check_compatibility(); //! convenience call to collection's size method inline size_t size() const; template class iterator { public: using value_type = const_corrector_t; using pointer = typename FullyTypedFieldMap::pointer; using difference_type = std::ptrdiff_t; using iterator_category = std::random_access_iterator_tag; using Ccoord = typename FieldCollection::Ccoord; //! Default constructor iterator() = delete; //! constructor inline iterator(FullyTypedFieldMap & fieldmap, bool begin=true); //! constructor for random access inline iterator(FullyTypedFieldMap & fieldmap, size_t index); //! Copy constructor iterator(const iterator &other)= default; //! Move constructor iterator(iterator &&other) noexcept = default; //! Destructor virtual ~iterator() noexcept = default; //! Copy assignment operator iterator& operator=(const iterator &other) = default; //! Move assignment operator iterator& operator=(iterator &&other) noexcept = default; //! pre-increment inline iterator & operator++(); //! post-increment inline iterator operator++(int); //! dereference inline value_type operator*(); //! member of pointer inline pointer operator->(); //! pre-decrement inline iterator & operator--(); //! post-decrement inline iterator operator--(int); //! access subscripting inline value_type operator[](difference_type diff); //! equality inline bool operator==(const iterator & other) const; //! inequality inline bool operator!=(const iterator & other) const; //! div. comparisons inline bool operator<(const iterator & other) const; inline bool operator<=(const iterator & other) const; inline bool operator>(const iterator & other) const; inline bool operator>=(const iterator & other) const; //! additions, subtractions and corresponding assignments inline iterator operator+(difference_type diff) const; inline iterator operator-(difference_type diff) const; inline iterator& operator+=(difference_type diff); inline iterator& operator-=(difference_type diff); //! get pixel coordinates inline Ccoord get_ccoord() const; //! ostream operator (mainly for debug friend std::ostream & operator<<(std::ostream & os, const iterator& it) { if (isConst) { os << "const "; } os << "iterator on field '" << it.fieldmap.get_name() << "', entry " << it.index; return os; } protected: const FieldCollection & collection; FullyTypedFieldMap & fieldmap; size_t index; private: }; protected: inline T* get_ptr_to_entry(size_t index); inline T& get_ref_to_entry(size_t index); const FieldCollection & collection; TypedField & field; private: }; } // internal namespace internal { /* ---------------------------------------------------------------------- */ template FieldMap:: FieldMap(Field & field) :collection(field.get_collection()), field(static_cast(field)) { static_assert(NbComponents>0, "Only fields with more than 0 components allowed"); } /* ---------------------------------------------------------------------- */ template void FieldMap:: check_compatibility() { if (typeid(T).hash_code() != this->field.get_stored_typeid().hash_code()) { throw FieldInterpretationError ("Cannot create a Map of type '" + this->info_string() + "' for field '" + this->field.get_name() + "' of type '" + this->field.get_stored_typeid().name() + "'"); } //check size compatibility if (NbComponents != this->field.get_nb_components()) { throw FieldInterpretationError ("Cannot create a Map of type '" + this->info_string() + "' for field '" + this->field.get_name() + "' with " + std::to_string(this->field.get_nb_components()) + " components"); } } /* ---------------------------------------------------------------------- */ template size_t FieldMap:: size() const { return this->collection.size(); } /* ---------------------------------------------------------------------- */ template constexpr size_t FieldMap:: nb_components() { return NbComponents; } /* ---------------------------------------------------------------------- */ template const std::string & FieldMap:: get_name() const { return this->field.get_name(); } /* ---------------------------------------------------------------------- */ template const FieldCollection & FieldMap:: get_collection() const { return this->collection; } /* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */ // Iterator implementations //! constructor template template FieldMap::iterator:: iterator(FullyTypedFieldMap & fieldmap, bool begin) :collection(fieldmap.get_collection()), fieldmap(fieldmap), index(begin ? 0 : fieldmap.field.size()) {} /* ---------------------------------------------------------------------- */ //! constructor for random access template template FieldMap::iterator:: iterator(FullyTypedFieldMap & fieldmap, size_t index) :collection(fieldmap.collection), fieldmap(fieldmap), index(index) {} /* ---------------------------------------------------------------------- */ //! pre-increment template template FieldMap::iterator & FieldMap::iterator:: operator++() { this->index++; return *this; } /* ---------------------------------------------------------------------- */ //! post-increment template template FieldMap::iterator FieldMap::iterator:: operator++(int) { iterator current = *this; this->index++; return current; } /* ---------------------------------------------------------------------- */ //! dereference template template typename FieldMap::template iterator::value_type FieldMap::iterator:: operator*() { return this->fieldmap.template operator[](this->index); } /* ---------------------------------------------------------------------- */ //! member of pointer template template typename FullyTypedFieldMap::pointer FieldMap::iterator:: operator->() { return this->fieldmap.ptr_to_val_t(this->index); } /* ---------------------------------------------------------------------- */ //! pre-decrement template template FieldMap::iterator & FieldMap::iterator:: operator--() { this->index--; return *this; } /* ---------------------------------------------------------------------- */ //! post-decrement template template FieldMap::iterator FieldMap::iterator:: operator--(int) { iterator current = *this; this->index--; return current; } /* ---------------------------------------------------------------------- */ //! Access subscripting template template typename FieldMap::template iterator::value_type FieldMap::iterator:: operator[](difference_type diff) { return this->fieldmap[this->index+diff]; } /* ---------------------------------------------------------------------- */ //! equality template template bool FieldMap::iterator:: operator==(const iterator & other) const { return (this->index == other.index && &this->fieldmap == &other.fieldmap); } /* ---------------------------------------------------------------------- */ //! inquality template template bool FieldMap::iterator:: operator!=(const iterator & other) const { return !(*this == other); } /* ---------------------------------------------------------------------- */ //! div. comparisons template template bool FieldMap::iterator:: operator<(const iterator & other) const { return (this->index < other.index); } template template bool FieldMap::iterator:: operator<=(const iterator & other) const { return (this->index <= other.index); } template template bool FieldMap::iterator:: operator>(const iterator & other) const { return (this->index > other.index); } template template bool FieldMap::iterator:: operator>=(const iterator & other) const { return (this->index >= other.index); } /* ---------------------------------------------------------------------- */ //! additions, subtractions and corresponding assignments template template FieldMap::iterator FieldMap::iterator:: operator+(difference_type diff) const { return iterator(this->fieldmap, this->index + diff); } template template FieldMap::iterator FieldMap::iterator:: operator-(difference_type diff) const { return iterator(this->fieldmap, this->index - diff); } template template FieldMap::iterator & FieldMap::iterator:: operator+=(difference_type diff) { this->index += diff; return *this; } template template FieldMap::iterator & FieldMap::iterator:: operator-=(difference_type diff) { this->index -= diff; return *this; } /* ---------------------------------------------------------------------- */ //! get pixel coordinates template template typename FieldCollection::Ccoord FieldMap::iterator:: get_ccoord() const { return this->collection.get_ccoord(this->index); } ////----------------------------------------------------------------------------// //template //std::ostream & operator << //(std::ostream &os, // const typename FieldMap:: // template iterator & it) { // os << "iterator on field '" // << it.field.get_name() // << "', entry " << it.index; // return os; //} /* ---------------------------------------------------------------------- */ template T* FieldMap::get_ptr_to_entry(size_t index) { return this->field.get_ptr_to_entry(std::move(index)); } /* ---------------------------------------------------------------------- */ template T& FieldMap::get_ref_to_entry(size_t index) { return this->field.get_ref_to_entry(std::move(index)); } } // internal } // muSpectre #endif /* FIELD_MAP_H */ diff --git a/src/common/field_map_matrixlike.hh b/src/common/field_map_matrixlike.hh index d6968ef..686081a 100644 --- a/src/common/field_map_matrixlike.hh +++ b/src/common/field_map_matrixlike.hh @@ -1,194 +1,194 @@ /** * file field_map_matrixlike.hh * * @author Till Junge * * @date 26 Sep 2017 * * @brief Eigen-Matrix and -Array maps over strongly typed fields * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef FIELD_MAP_MATRIXLIKE_H #define FIELD_MAP_MATRIXLIKE_H #include #include #include "common/field_map_base.hh" namespace muSpectre { namespace internal { /* ---------------------------------------------------------------------- */ enum class Map_t{Matrix, Array}; template struct NameWrapper { const static std::string field_info_root; }; template<> const std::string NameWrapper::field_info_root{"Array"}; template<> const std::string NameWrapper::field_info_root{"Matrix"}; /* ---------------------------------------------------------------------- */ template class MatrixLikeFieldMap: public FieldMap { public: - using parent = FieldMap; using Ccoord = Ccoord_t; using value_type = EigenArray; using reference = value_type; // since it's a resource handle using const_reference = const value_type; - using size_type = typename parent::size_type; + using size_type = typename Parent::size_type; using pointer = std::unique_ptr; - using TypedField = typename parent::TypedField; - using Field = typename TypedField::parent; - using iterator = typename parent::template iterator; - using const_iterator= typename parent::template iterator; + using TypedField = typename Parent::TypedField; + using Field = typename TypedField::Parent; + using iterator = typename Parent::template iterator; + using const_iterator= typename Parent::template iterator; using reverse_iterator = std::reverse_iterator; using const_reverse_iterator = std::reverse_iterator; friend iterator; //! Default constructor MatrixLikeFieldMap() = delete; MatrixLikeFieldMap(Field & field); //! Copy constructor MatrixLikeFieldMap(const MatrixLikeFieldMap &other) = delete; //! Move constructor MatrixLikeFieldMap(MatrixLikeFieldMap &&other) noexcept = delete; //! Destructor virtual ~MatrixLikeFieldMap() noexcept = default; //! Copy assignment operator MatrixLikeFieldMap& operator=(const MatrixLikeFieldMap &other) = delete; //! Move assignment operator MatrixLikeFieldMap& operator=(MatrixLikeFieldMap &&other) noexcept = delete; //! give human-readable field map type inline std::string info_string() const override final; //! member access template inline ref operator[](size_type index); inline reference operator[](const Ccoord& ccoord); //! return an iterator to head of field for ranges inline iterator begin(){return iterator(*this);} inline const_iterator cbegin(){return const_iterator(*this);} //! return an iterator to tail of field for ranges inline iterator end(){return iterator(*this, false);}; inline const_iterator cend(){return const_iterator(*this, false);} protected: //! for sad, legacy iterator use inline pointer ptr_to_val_t(size_type index); const static std::string field_info_root; private: }; template const std::string MatrixLikeFieldMap:: field_info_root{NameWrapper::field_info_root}; /* ---------------------------------------------------------------------- */ template MatrixLikeFieldMap:: MatrixLikeFieldMap(Field & field) - :parent(field) { + :Parent(field) { this->check_compatibility(); } /* ---------------------------------------------------------------------- */ //! human-readable field map type template std::string MatrixLikeFieldMap:: info_string() const { std::stringstream info; info << field_info_root << "(" << typeid(typename EigenArray::value_type).name() << ", " << EigenArray::RowsAtCompileTime << "x" << EigenArray::ColsAtCompileTime << ")"; return info.str(); } /* ---------------------------------------------------------------------- */ //! member access template template ref MatrixLikeFieldMap:: operator[](size_type index) { return ref(this->get_ptr_to_entry(index)); } template typename MatrixLikeFieldMap::reference MatrixLikeFieldMap:: operator[](const Ccoord & ccoord) { auto && index = this->collection.get_index(ccoord); return reference(this->get_ptr_to_entry(std::move(index))); } /* ---------------------------------------------------------------------- */ template typename MatrixLikeFieldMap::pointer MatrixLikeFieldMap:: ptr_to_val_t(size_type index) { return std::make_unique (this->get_ptr_to_entry(std::move(index))); } } // internal /* ---------------------------------------------------------------------- */ //! short-hand for an Eigen matrix map as iterate template using MatrixFieldMap = internal::MatrixLikeFieldMap >, internal::Map_t::Matrix>; /* ---------------------------------------------------------------------- */ //! short-hand for an Eigen array map as iterate template using ArrayFieldMap = internal::MatrixLikeFieldMap >, internal::Map_t::Array>; } // muSpectre #endif /* FIELD_MAP_MATRIXLIKE_H */ diff --git a/src/common/field_map_tensor.hh b/src/common/field_map_tensor.hh index 88bf4c8..ee27a3e 100644 --- a/src/common/field_map_tensor.hh +++ b/src/common/field_map_tensor.hh @@ -1,169 +1,169 @@ /** * file field_map_tensor.hh * * @author Till Junge * * @date 26 Sep 2017 * * @brief Defines an Eigen-Tensor map over strongly typed fields * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef FIELD_MAP_TENSOR_H #define FIELD_MAP_TENSOR_H #include #include #include "common/eigen_tools.hh" #include "common/field_map_base.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ template class TensorFieldMap: public internal::FieldMap::Sizes::total_size > { public: - using parent = internal::FieldMap; using Ccoord = Ccoord_t; using Sizes = typename SizesByOrder::Sizes; using T_t = Eigen::TensorFixedSize; using value_type = Eigen::TensorMap; using reference = value_type; // since it's a resource handle using const_reference = Eigen::TensorMap; - using size_type = typename parent::size_type; + using size_type = typename Parent::size_type; using pointer = std::unique_ptr; - using TypedField = typename parent::TypedField; - using Field = typename TypedField::parent; - using iterator = typename parent::template iterator; - using const_iterator = typename parent::template iterator; + using TypedField = typename Parent::TypedField; + using Field = typename TypedField::Parent; + using iterator = typename Parent::template iterator; + using const_iterator = typename Parent::template iterator; using reverse_iterator = std::reverse_iterator; using const_reverse_iterator = std::reverse_iterator; friend iterator; //! Default constructor TensorFieldMap() = delete; TensorFieldMap(Field & field); //! Copy constructor TensorFieldMap(const TensorFieldMap &other) = delete; //! Move constructor TensorFieldMap(TensorFieldMap &&other) noexcept = delete; //! Destructor virtual ~TensorFieldMap() noexcept = default; //! Copy assignment operator TensorFieldMap& operator=(const TensorFieldMap &other) = delete; //! Move assignment operator TensorFieldMap& operator=(TensorFieldMap &&other) noexcept = delete; //! give human-readable field map type inline std::string info_string() const override final; //! member access template inline ref_t operator[](size_type index); inline reference operator[](const Ccoord & ccoord); //! return an iterator to head of field for ranges inline iterator begin(){return iterator(*this);} inline const_iterator cbegin(){return const_iterator(*this);} //! return an iterator to tail of field for ranges inline iterator end(){return iterator(*this, false);} inline const_iterator cend(){return const_iterator(*this, false);} protected: //! for sad, legacy iterator use inline pointer ptr_to_val_t(size_type index); private: }; /* ---------------------------------------------------------------------- */ template TensorFieldMap:: TensorFieldMap(Field & field) - :parent(field) { + :Parent(field) { this->check_compatibility(); } /* ---------------------------------------------------------------------- */ //! human-readable field map type template std::string TensorFieldMap::info_string() const { std::stringstream info; info << "Tensor(" << typeid(T).name() << ", " << order << "_o, " << dim << "_d)"; return info.str(); } /* ---------------------------------------------------------------------- */ //! member access template template ref_t TensorFieldMap::operator[](size_type index) { auto && lambda = [this, &index](auto&&...tens_sizes) { return ref_t(this->get_ptr_to_entry(index), tens_sizes...); }; return call_sizes(lambda); } template typename TensorFieldMap::reference TensorFieldMap:: operator[](const Ccoord & ccoord) { auto && index = this->collection.get_index(ccoord); auto && lambda = [this, &index](auto&&...sizes) { return ref_t(this->get_ptr_to_entry(index), sizes...); }; return call_sizes(lambda); } /* ---------------------------------------------------------------------- */ //! for sad, legacy iterator use. Don't use unless you have to. template typename TensorFieldMap::pointer TensorFieldMap:: ptr_to_val_t(size_type index) { auto && lambda = [this, &index](auto&&... tens_sizes) { return std::make_unique(this->get_ptr_to_entry(index), tens_sizes...); }; return call_sizes(lambda); } } // muSpectre #endif /* FIELD_MAP_TENSOR_H */ diff --git a/src/materials/material.hh b/src/materials/material.hh index 099a8bb..d8e1bd5 100644 --- a/src/materials/material.hh +++ b/src/materials/material.hh @@ -1,135 +1,136 @@ /** * file material.hh * * @author Till Junge * * @date 25 Oct 2017 * * @brief Base class for materials (constitutive models) * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include #include "common/common.hh" #include "common/field_map_tensor.hh" #include "common/field_collection.hh" #ifndef MATERIAL_H #define MATERIAL_H namespace muSpectre { //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) template class MaterialBase { public: //! typedefs for data handled by this interface //! global field collection for system-wide fields, like stress, strain, etc using GFieldCollection_t = FieldCollection; //! field collection for internal variables, such as eigen-strains, //! plastic strains, damage variables, etc, but also for managing which //! pixels the material is responsible for using MFieldCollection_t = FieldCollection; using StressMap_t = TensorFieldMap; using StrainMap_t = StressMap_t; using StiffnessMap_t = TensorFieldMap; using Ccoord = Ccoord_t; //! Default constructor MaterialBase() = delete; //! Construct by name MaterialBase(std::string name); //! Copy constructor MaterialBase(const MaterialBase &other) = delete; //! Move constructor MaterialBase(MaterialBase &&other) noexcept = delete; //! Destructor virtual ~MaterialBase() noexcept = default; //! Copy assignment operator MaterialBase& operator=(const MaterialBase &other) = delete; //! Move assignment operator MaterialBase& operator=(MaterialBase &&other) noexcept = delete; //! take responsibility for a pixel identified by its cell coordinates //! TODO: this won't work. for materials with additional info per pixel (as, e.g. for eigenstrain), we need to pass more parameters, so this needs to be a variadic function. void add_pixel(const Ccoord & ccord); //! allocate memory, etc //! TODO: this won't work. for materials with additional info per pixel (see above TODO), we neet to allocate before we know for sure how many pixels the material is responsible for. - virtual void initialize() = 0; + virtual void initialize(bool stiffness = false) = 0; //! computes the first Piola-Kirchhoff stress for finite strain problems virtual void compute_stress(const StrainMap_t & F, StressMap_t & P); //! computes the first Piola-Kirchhoff stress and the tangent stiffness //! for finite strain problems virtual void compute_stress_stiffness(const StrainMap_t & F, StressMap_t & P, StiffnessMap_t & K); //! computes Cauchy stress for small strain problems virtual void compute_cauchy(const StrainMap_t & eps, StressMap_t &sig); //! computes Cauchy stress and stiffness for small strain problems virtual void compute_cauchy_stiffness(const StrainMap_t & eps, StressMap_t &sig, StiffnessMap_t & K); //! return the materil's name const std::string & get_name() const; //! for static inheritance stuff constexpr static Dim_t sdim() {return DimS;} constexpr static Dim_t mdim() {return DimM;} protected: //! computes stress virtual void compute_stress(const StrainMap_t & F, StressMap_t & P, Formulation form); //! computes stress and tangent stiffness //! for finite strain problems virtual void compute_stress_stiffness(const StrainMap_t & F, StressMap_t & P, StiffnessMap_t & K, Formulation form); //! members const std::string name; MFieldCollection_t internal_fields; + private: }; } // muSpectre #endif /* MATERIAL_H */ diff --git a/src/materials/material_muSpectre_base.hh b/src/materials/material_muSpectre_base.hh index ab65a76..89dd4ec 100644 --- a/src/materials/material_muSpectre_base.hh +++ b/src/materials/material_muSpectre_base.hh @@ -1,78 +1,80 @@ /** * file material_muSpectre_base.hh * * @author Till Junge * * @date 25 Oct 2017 * * @brief Base class for materials written for µSpectre specifically. These * can take full advantage of the configuration-change utilities of * µSpectre. The user can inherit from them to define new constitutive * laws and is merely required to provide the methods for computing the * second Piola-Kirchhoff stress and Stiffness. This class uses the * "curiously recurring template parameter" to avoid virtual calls. * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "materials/material.hh" #ifndef MATERIAL_MUSPECTRE_BASE_H #define MATERIAL_MUSPECTRE_BASE_H namespace muSpectre { + //! 'Material' is a CRTP template class MaterialMuSpectre: public> { public: + //! Default constructor MaterialMuSpectre() = delete; //! Construct by name MaterialMuSpectre(std::string name); //! Copy constructor MaterialMuSpectre(const MaterialMuSpectre &other) = delete; //! Move constructor MaterialMuSpectre(MaterialMuSpectre &&other) noexcept = delete; //! Destructor virtual ~MaterialMuSpectre() noexcept = default; //! Copy assignment operator MaterialMuSpectre& operator=(const MaterialMuSpectre &other) = delete; //! Move assignment operator MaterialMuSpectre& operator=(MaterialMuSpectre &&other) noexcept = delete; protected: private: }; } // muSpectre #endif /* MATERIAL_MUSPECTRE_BASE_H */ diff --git a/src/materials/materials_toolbox.hh b/src/materials/materials_toolbox.hh new file mode 100644 index 0000000..026980c --- /dev/null +++ b/src/materials/materials_toolbox.hh @@ -0,0 +1,154 @@ +/** + * file materials_toolbox.hh + * + * @author Till Junge + * + * @date 02 Nov 2017 + * + * @brief collection of common continuum mechanics tools + * + * @section LICENCE + * + * Copyright (C) 2017 Till Junge + * + * µSpectre is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 3, or (at + * your option) any later version. + * + * µSpectre is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Emacs; see the file COPYING. If not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ + +#ifndef MATERIALS_TOOLBOX_H +#define MATERIALS_TOOLBOX_H + +#include +#include +#include +#include "common/common.hh" +#include "materials/tensor_algebra.hh" + +namespace muSpectre { + + namespace MatTB { + class MaterialsToolboxError:public std::runtime_error{ + public: + explicit MaterialsToolboxError(const std::string& what) + :std::runtime_error(what){} + explicit MaterialsToolboxError(const char * what) + :std::runtime_error(what){} + }; + + /* ---------------------------------------------------------------------- */ + //! Material laws can declare which type of stress measure they provide, + //! and µSpectre will handle conversions + enum class StressMeasure { + Cauchy, PK1, PK2, Kirchhoff, Biot, Mandel}; + std::ostream & operator<<(std::ostream & os, StressMeasure s) { + switch (s) { + case StressMeasure::Cauchy: {os << "Cauchy"; break;} + case StressMeasure::PK1: {os << "PK1"; break;} + case StressMeasure::PK1: {os << "PK1"; break;} + case StressMeasure::PK2: {os << "PK2"; break;} + case StressMeasure::Kirchhoff: {os << "Kirchhoff"; break;} + case StressMeasure::Biot: {os << "Biot"; break;} + case StressMeasure::Mandel: {os << "Mandel"; break;} + default: + throw MaterialsToolboxError + ("a stress measure must be missing"); + break; + } + return os; + } + + /* ---------------------------------------------------------------------- */ + //! Material laws can declare which type of strain measure they require and + //! µSpectre will provide it + enum class StrainMeasure { + GreenLagrange, Biot, Log, Almansi, RCauchyGreen, LCauchyGreen}; + std::ostream & operator<<(std::ostream & os, StrainMeasure s) { + switch (s) { + case StrainMeasure::GreenLagrange: {os << "Green-Lagrange"; break;} + case StrainMeasure::Biot: {os << "Biot"; break;} + case StrainMeasure::Log: {os << "Logarithmic"; break;} + case StrainMeasure::Almansi: {os << "Almansi"; break;} + case StrainMeasure::RCauchyGreen: {os << "Right Cauchy-Green"; break;} + case StrainMeasure::LCauchyGreen: {os << "Left Cauchy-Green"; break;} + default: + throw MaterialsToolboxError + ("a strain measure must be missing"); + break; + } + } + + /* ---------------------------------------------------------------------- */ + //! set of functions returning an expression for PK2 stress based on + template + decltype(auto) PK2_stress(Tens1 && stress, Tens2 && Strain) { + // the following test always fails to generate a compile-time error + static_assert((StressM == StressMeasure::Cauchy) && + (StressM == StressMeasure::PK1), + "The requested Stress conversion is not implemented. " + "You either made a programming mistake or need to " + "implement it as a specialisation of this function. " + "See PK2stress for an example."); + } + + /* ---------------------------------------------------------------------- */ + template + decltype(auto) + PK2_stress(Tens1 && S, Tens2 && F) { + return F*S; + } + + /* ---------------------------------------------------------------------- */ + //! set of functions returning ∂strain/∂F (where F is the transformation + //! gradient. (e.g. for a material law expressed as PK2 stress S as + //! function of Green-Lagrange strain E, this would return F^T ⊗̲ I) (⊗_ + //! refers to Curnier's notation) WARNING: this function assumes that your + //! stress measure has the two minor symmetries due to stress equilibrium + //! and the major symmetry due to the existance of an elastic potential + //! W(E). Do not rely on this function if you are using exotic stress + //! measures + template + decltype(auto) StrainDerivative(Tens && F) { + // the following test always fails to generate a compile-time error + static_assert((StrainM == StrainMeasure::Almansi) && + (StrainM == StrainMeasure::Biot), + "The requested StrainDerivative calculation is not " + "implemented. You either made a programming mistake or " + "need to implement it as a specialisation of this " + "function. See " + "StrainDerivative for an " + "example."); + } + + /* ---------------------------------------------------------------------- */ + template + decltype(auto) + StrainDerivative (Tens && F) { + constexpr size_t order{2}; + constexpr Dim_t dim{F.Dimensions[0]}; + return Tensors::outer_under + (F.shuffle(std::array{1,0}), Tensors::I2()); + + /* ---------------------------------------------------------------------- */ + //! function returning expressions for PK2 stress and stiffness + template + decltype(auto) PK2_stress_stiffness(Tens1 && stress, Tens2 && F) { + } + + } // MatTB + +} // muSpectre + +#endif /* MATERIALS_TOOLBOX_H */ diff --git a/src/materials/tensor_algebra.hh b/src/materials/tensor_algebra.hh new file mode 100644 index 0000000..2ec2345 --- /dev/null +++ b/src/materials/tensor_algebra.hh @@ -0,0 +1,95 @@ +/** + * file tensor_algebra.hh + * + * @author Till Junge + * + * @date 05 Nov 2017 + * + * @brief collection of compile-time quantities and algrebraic functions for + * tensor operations + * + * @section LICENCE + * + * Copyright (C) 2017 Till Junge + * + * µSpectre is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 3, or (at + * your option) any later version. + * + * µSpectre is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Emacs; see the file COPYING. If not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ + +#ifndef TENSOR_ALGEBRA_H +#define TENSOR_ALGEBRA_H + +#include "common/common.hh" +#include +#include + + +namespace muSpectre { + + namespace Tensors { + + template + using Tens2_t = Eigen::TensorFixedSize>; + template + using Tens4_t = Eigen::TensorFixedSize>; + + //----------------------------------------------------------------------------// + //! compile-time second-order identity + template + constexpr inline Tens2_t I2() { + Tens2_t T; + using Mat_t = Eigen::Matrix; + Eigen::Map(&T(0, 0)) = Mat_t::Identity(); + return T; + } + + /* ---------------------------------------------------------------------- */ + /** compile-time outer tensor product as defined by Curnier + * R_ijkl = A_ij.B_klxx + * 0123 01 23 + */ + template + constexpr inline decltype(auto) outer(T1 && A, T2 && B) { + std::array, 0> dims{}; + return A.contract(B, dims); + } + + /* ---------------------------------------------------------------------- */ + /** compile-time underlined outer tensor product as defined by Curnier + * R_ijkl = A_ik.B_jlxx + * 0123 02 13 + */ + template + constexpr inline decltype(auto) outer_under(T1 && A, T2 && B) { + constexpr size_t order{4}; + return outer(A, B).shuffle(std::array{0, 2, 1, 3}); + } + + /* ---------------------------------------------------------------------- */ + /** compile-time overlined outer tensor product as defined by Curnier + * R_ijkl = A_il.B_jkxx + * 0123 03 12 + */ + template + constexpr inline decltype(auto) outer_over(T1 && A, T2 && B) { + constexpr size_t order{4}; + return outer(A, B).shuffle(std::array{0, 3, 1, 2}); + } + + +} // muSpectre + + +#endif /* TENSOR_ALGEBRA_H */ diff --git a/tests/test_field_collections.cc b/tests/test_field_collections.cc index 21b81e3..43fdd61 100644 --- a/tests/test_field_collections.cc +++ b/tests/test_field_collections.cc @@ -1,523 +1,547 @@ /** * file test_field_collections.cc * * @author Till Junge * * @date 20 Sep 2017 * * @brief Test the FieldCollection classes which provide fast optimized iterators * over run-time typed fields * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include #include #include #include #include #include #include "common/common.hh" #include "common/ccoord_operations.hh" #include "common/test_goodies.hh" #include "tests.hh" #include "common/field_collection.hh" #include "common/field.hh" #include "common/field_map.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(field_collection_tests); //! Test fixture for simple tests on single field in collection template struct FC_fixture: public FieldCollection { FC_fixture() :fc() {} inline static constexpr Dim_t sdim(){return DimS;} inline static constexpr Dim_t mdim(){return DimM;} inline static constexpr bool global(){return Global;} using FC_t = FieldCollection; FC_t fc; }; using test_collections = boost::mpl::list, FC_fixture<2, 3, true>, FC_fixture<3, 3, true>, FC_fixture<2, 2, false>, FC_fixture<2, 3, false>, FC_fixture<3, 3, false>>; BOOST_AUTO_TEST_CASE(simple) { const Dim_t sdim = 2; const Dim_t mdim = 2; using FC_t = FieldCollection; FC_t fc; BOOST_CHECK_EQUAL(FC_t::spatial_dim(), sdim); BOOST_CHECK_EQUAL(fc.get_spatial_dim(), sdim); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(Simple_construction_test, F, test_collections, F) { BOOST_CHECK_EQUAL(F::FC_t::spatial_dim(), F::sdim()); BOOST_CHECK_EQUAL(F::fc.get_spatial_dim(), F::sdim()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(get_field2_test, F, test_collections, F) { const auto order{2}; using FC_t = typename F::FC_t; using TF_t = TensorField; auto && myfield = make_field("TensorField real 2", F::fc); using TensorMap = TensorFieldMap; using MatrixMap = MatrixFieldMap; using ArrayMap = ArrayFieldMap; TensorMap TFM(myfield); MatrixMap MFM(myfield); ArrayMap AFM(myfield); BOOST_CHECK_EQUAL(TFM.info_string(), "Tensor(d, "+ std::to_string(order) + "_o, " + std::to_string(F::mdim()) + "_d)"); BOOST_CHECK_EQUAL(MFM.info_string(), "Matrix(d, "+ std::to_string(F::mdim()) + "x" + std::to_string(F::mdim()) + ")"); BOOST_CHECK_EQUAL(AFM.info_string(), "Array(d, "+ std::to_string(F::mdim()) + "x" + std::to_string(F::mdim()) + ")"); } constexpr Dim_t order{4}, matrix_order{2}; //! Test fixture for multiple fields in the collection template struct FC_multi_fixture{ FC_multi_fixture() :fc() { //add Real tensor field make_field> ("Tensorfield Real o4", fc); //add Real tensor field make_field> ("Tensorfield Real o2", fc); //add integer scalar field make_field> ("integer Scalar", fc); //add complex matrix field make_field> ("Matrixfield Complex sdim x mdim", fc); } inline static constexpr Dim_t sdim(){return DimS;} inline static constexpr Dim_t mdim(){return DimM;} inline static constexpr bool global(){return Global;} using FC_t = FieldCollection; FC_t fc; }; using mult_collections = boost::mpl::list, FC_multi_fixture<2, 3, true>, FC_multi_fixture<3, 3, true>, FC_multi_fixture<2, 2, false>, FC_multi_fixture<2, 3, false>, FC_multi_fixture<3, 3, false>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(multi_field_test, F, mult_collections, F) { using FC_t = typename F::FC_t; // possible maptypes for Real tensor fields using T_type = Real; using T_TFM1_t = TensorFieldMap; using T_TFM2_t = TensorFieldMap; //! dangerous // impossible maptypes for Real tensor fields using T_SFM_t = ScalarFieldMap; using T_MFM_t = MatrixFieldMap; using T_AFM_t = ArrayFieldMap; using T_MFMw1_t = MatrixFieldMap; using T_MFMw2_t = MatrixFieldMap; using T_MFMw3_t = MatrixFieldMap; const std::string T_name{"Tensorfield Real o4"}; const std::string T_name_w{"TensorField Real o4 wrongname"}; BOOST_CHECK_THROW(T_SFM_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_NO_THROW(T_TFM1_t(F::fc.at(T_name))); BOOST_CHECK_NO_THROW(T_TFM2_t(F::fc.at(T_name))); BOOST_CHECK_THROW(T_MFM_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_AFM_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw1_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw2_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw2_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw3_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_SFM_t(F::fc.at(T_name_w)), std::out_of_range); // possible maptypes for integer scalar fields using S_type = Int; using S_SFM_t = ScalarFieldMap; using S_TFM1_t = TensorFieldMap; using S_TFM2_t = TensorFieldMap; using S_MFM_t = MatrixFieldMap; using S_AFM_t = ArrayFieldMap; // impossible maptypes for integer scalar fields using S_MFMw1_t = MatrixFieldMap; using S_MFMw2_t = MatrixFieldMap; using S_MFMw3_t = MatrixFieldMap; const std::string S_name{"integer Scalar"}; const std::string S_name_w{"integer Scalar wrongname"}; BOOST_CHECK_NO_THROW(S_SFM_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_TFM1_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_TFM2_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_MFM_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_AFM_t(F::fc.at(S_name))); BOOST_CHECK_THROW(S_MFMw1_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_MFMw2_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_MFMw2_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_MFMw3_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_SFM_t(F::fc.at(S_name_w)), std::out_of_range); // possible maptypes for complex matrix fields using M_type = Complex; using M_MFM_t = MatrixFieldMap; using M_AFM_t = ArrayFieldMap; // impossible maptypes for complex matrix fields using M_SFM_t = ScalarFieldMap; using M_MFMw1_t = MatrixFieldMap; using M_MFMw2_t = MatrixFieldMap; using M_MFMw3_t = MatrixFieldMap; const std::string M_name{"Matrixfield Complex sdim x mdim"}; const std::string M_name_w{"Matrixfield Complex sdim x mdim wrongname"}; BOOST_CHECK_THROW(M_SFM_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_NO_THROW(M_MFM_t(F::fc.at(M_name))); BOOST_CHECK_NO_THROW(M_AFM_t(F::fc.at(M_name))); BOOST_CHECK_THROW(M_MFMw1_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_MFMw2_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_MFMw2_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_MFMw3_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_SFM_t(F::fc.at(M_name_w)), std::out_of_range); } /* ---------------------------------------------------------------------- */ //! Check whether fields can be initialized using mult_collections_t = boost::mpl::list, FC_multi_fixture<2, 3, true>, FC_multi_fixture<3, 3, true>>; using mult_collections_f = boost::mpl::list, FC_multi_fixture<2, 3, false>, FC_multi_fixture<3, 3, false>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_glob, F, mult_collections_t, F) { Ccoord_t size; for (auto && s: size) { s = 3; } - F::fc.initialise(size); + BOOST_CHECK_NO_THROW(F::fc.initialise(size)); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_loca, F, mult_collections_f, F) { testGoodies::RandRange rng; for (int i = 0; i < 7; ++i) { Ccoord_t pixel; for (auto && s: pixel) { s = rng.randval(0, 7); } F::fc.add_pixel(pixel); } - F::fc.initialise(); + BOOST_CHECK_NO_THROW(F::fc.initialise()); + } + + BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_loca_with_push_back, F, mult_collections_f, F) { + constexpr auto mdim = F::mdim(); + constexpr int nb_pix = 7; + testGoodies::RandRange rng; + using ftype = internal::TypedFieldBase; + using stype = Eigen::Array; + auto & field = reinterpret_cast(F::fc["Tensorfield Real o4"]); + field.push_back(stype()); + for (int i = 0; i < nb_pix; ++i) { + Ccoord_t pixel; + for (auto && s: pixel) { + s = rng.randval(0, 7); + } + F::fc.add_pixel(pixel); + } + + BOOST_CHECK_THROW(F::fc.initialise(), FieldCollectionError); + for (int i = 0; i < nb_pix-1; ++i) { + field.push_back(stype()); + } + BOOST_CHECK_NO_THROW(F::fc.initialise()); + } //! Test fixture for iterators over multiple fields template struct FC_iterator_fixture : public FC_multi_fixture { - using parent = FC_multi_fixture; + using Parent = FC_multi_fixture; FC_iterator_fixture() - :parent() { + :Parent() { this-> fill(); } template std::enable_if_t fill() { static_assert(Global==isGlobal, "You're breaking my SFINAE plan"); - Ccoord_t size; + Ccoord_t size; for (auto && s: size) { s = cube_size(); } this->fc.initialise(size); } template std::enable_if_t fill (int dummy = 0) { static_assert(notGlobal != Global, "You're breaking my SFINAE plan"); testGoodies::RandRange rng; for (int i = 0*dummy; i < sele_size(); ++i) { - Ccoord_t pixel; + Ccoord_t pixel; for (auto && s: pixel) { s = rng.randval(0, 7); } this->fc.add_pixel(pixel); } this->fc.initialise(); } constexpr static Dim_t cube_size() {return 3;} constexpr static Dim_t sele_size() {return 7;} }; using iter_collections = boost::mpl::list, FC_iterator_fixture<2, 3, true>, FC_iterator_fixture<3, 3, true>, FC_iterator_fixture<2, 2, false>, FC_iterator_fixture<2, 3, false>, FC_iterator_fixture<3, 3, false>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(iter_field_test, F, iter_collections, F) { - using FC_t = typename F::parent::FC_t; - using Tensor4Map = TensorFieldMap; + using FC_t = typename F::Parent::FC_t; + using Tensor4Map = TensorFieldMap; Tensor4Map T4map{F::fc["Tensorfield Real o4"]}; F::fc["Tensorfield Real o4"].set_zero(); for (auto && tens:T4map) { BOOST_CHECK_EQUAL(Real(Eigen::Tensor(tens.abs().sum().eval())()), 0); } - using Tensor2Map = TensorFieldMap; - using MSqMap = MatrixFieldMap; - using ASqMap = ArrayFieldMap; + using Tensor2Map = TensorFieldMap; + using MSqMap = MatrixFieldMap; + using ASqMap = ArrayFieldMap; Tensor2Map T2map{F::fc["Tensorfield Real o2"]}; MSqMap Mmap{F::fc["Tensorfield Real o2"]}; ASqMap Amap{F::fc["Tensorfield Real o2"]}; auto t2_it = T2map.begin(); auto t2_it_end = T2map.end(); auto m_it = Mmap.begin(); auto a_it = Amap.begin(); for (; t2_it != t2_it_end; ++t2_it, ++m_it, ++a_it) { t2_it->setRandom(); auto && m = *m_it; bool comp = (m == a_it->matrix()); BOOST_CHECK(comp); } using ScalarMap = ScalarFieldMap; ScalarMap s_map{F::fc["integer Scalar"]}; for (Uint i = 0; i < s_map.size(); ++i) { s_map[i] = i; } Uint counter{0}; for (const auto& val: s_map) { BOOST_CHECK_EQUAL(counter++, val); } } using glob_iter_colls = boost::mpl::list, FC_iterator_fixture<2, 3, true>, FC_iterator_fixture<3, 3, true>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(ccoord_indexing_test, F, glob_iter_colls, F) { - using FC_t = typename F::parent::FC_t; + using FC_t = typename F::Parent::FC_t; using ScalarMap = ScalarFieldMap; ScalarMap s_map{F::fc["integer Scalar"]}; for (Uint i = 0; i < s_map.size(); ++i) { s_map[i] = i; } for (size_t i = 0; i < CcoordOps::get_size(F::fc.get_sizes()); ++i) { Ccoord_t sizes(F::fc.get_sizes()); Ccoord_t ccoords(CcoordOps::get_ccoord(sizes, i)); BOOST_CHECK_EQUAL(CcoordOps::get_index(F::fc.get_sizes(), CcoordOps::get_ccoord(F::fc.get_sizes(), i)), i); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(iterator_methods_test, F, iter_collections, F) { - using FC_t = typename F::parent::FC_t; - using Tensor4Map = TensorFieldMap; + using FC_t = typename F::Parent::FC_t; + using Tensor4Map = TensorFieldMap; Tensor4Map T4map{F::fc["Tensorfield Real o4"]}; using it_t = typename Tensor4Map::iterator; std::ptrdiff_t diff{3}; // arbitrary, as long as it is smaller than the container size // check constructors auto itstart = T4map.begin(); // standard way of obtaining iterator auto itend = T4map.end(); // ditto it_t it1{T4map}; it_t it2{T4map, false}; it_t it3{T4map, size_t(diff)}; BOOST_CHECK(itstart == itstart); BOOST_CHECK(itstart != itend); BOOST_CHECK_EQUAL(itstart, it1); BOOST_CHECK_EQUAL(itend, it2); // check ostream operator std::stringstream response; response << it3; BOOST_CHECK_EQUAL (response.str(), std::string ("iterator on field 'Tensorfield Real o4', entry ") + std::to_string(diff)); // check copy, move, and assigment constructor (and operator+) it_t itcopy = it1; it_t itmove = std::move(T4map.begin()); it_t it4 = it1+diff; it_t it5(it2); it_t tmp(it2); it_t it6(std::move(tmp)); it_t it7 = it4 -diff; BOOST_CHECK_EQUAL(itcopy, it1); BOOST_CHECK_EQUAL(itmove, it1); BOOST_CHECK_EQUAL(it4, it3); BOOST_CHECK_EQUAL(it5, it2); BOOST_CHECK_EQUAL(it6, it5); BOOST_CHECK_EQUAL(it7, it1); // check increments/decrements BOOST_CHECK_EQUAL(it1++, itcopy); // post-increment BOOST_CHECK_EQUAL(it1, itcopy+1); BOOST_CHECK_EQUAL(--it1, itcopy); // pre -decrement BOOST_CHECK_EQUAL(++it1, itcopy+1); // pre -increment BOOST_CHECK_EQUAL(it1--, itcopy+1); // post-decrement BOOST_CHECK_EQUAL(it1, itcopy); // dereference and member-of-pointer check Eigen::Tensor Tens = *it1; Eigen::Tensor Tens2 = *itstart; Eigen::Tensor check = (Tens==Tens2).all(); BOOST_CHECK_EQUAL(bool(check()), true); BOOST_CHECK_NO_THROW(itstart->setZero()); //check access subscripting auto T3a = *it3; auto T3b = itstart[diff]; BOOST_CHECK(bool(Eigen::Tensor((T3a==T3b).all())())); // div. comparisons BOOST_CHECK_LT(itstart, itend); BOOST_CHECK(!(itend < itstart)); BOOST_CHECK(!(itstart < itstart)); BOOST_CHECK_LE(itstart, itend); BOOST_CHECK_LE(itstart, itstart); BOOST_CHECK(!(itend <= itstart)); BOOST_CHECK_GT(itend, itstart); BOOST_CHECK(!(itend>itend)); BOOST_CHECK(!(itstart>itend)); BOOST_CHECK_GE(itend, itstart); BOOST_CHECK_GE(itend, itend); BOOST_CHECK(!(itstart >= itend)); // check assignment increment/decrement BOOST_CHECK_EQUAL(it1+=diff, it3); BOOST_CHECK_EQUAL(it1-=diff, itstart); // check cell coordinates using Ccoord = Ccoord_t; Ccoord a{itstart.get_ccoord()}; Ccoord b{}; // Weirdly, boost::has_left_shift::value is falso for Ccoords, even though the operator is implemented :( //BOOST_CHECK_EQUAL(a, b); bool check2 = (a==b); BOOST_CHECK(check2); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_tensor_iter_test, F, iter_collections, F) { - using FC_t = typename F::parent::FC_t; - using Tensor4Map = TensorFieldMap; + using FC_t = typename F::Parent::FC_t; + using Tensor4Map = TensorFieldMap; Tensor4Map T4map{F::fc["Tensorfield Real o4"]}; using T_t = typename Tensor4Map::T_t; - Eigen::TensorMap Tens2(T4map[0].data(), F::parent::sdim(), F::parent::sdim(), F::parent::sdim(), F::parent::sdim()); + Eigen::TensorMap Tens2(T4map[0].data(), F::Parent::sdim(), F::Parent::sdim(), F::Parent::sdim(), F::Parent::sdim()); for (auto it = T4map.cbegin(); it != T4map.cend(); ++it) { // maps to const tensors can't be initialised with a const pointer this sucks auto&& tens = *it; auto&& ptr = tens.data(); static_assert(std::is_pointer>::value, "should be getting a pointer"); //static_assert(std::is_const>::value, "should be const"); // If Tensor were written well, above static_assert should pass, and the // following check shouldn't. If it get triggered, it means that a newer // version of Eigen now does have const-correct // TensorMap. This means that const-correct field maps // are then also possible for tensors BOOST_CHECK(!std::is_const>::value); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_matrix_iter_test, F, iter_collections, F) { - using FC_t = typename F::parent::FC_t; + using FC_t = typename F::Parent::FC_t; using MatrixMap = MatrixFieldMap; MatrixMap Mmap{F::fc["Matrixfield Complex sdim x mdim"]}; for (auto it = Mmap.cbegin(); it != Mmap.cend(); ++it) { // maps to const tensors can't be initialised with a const pointer this sucks auto&& mat = *it; auto&& ptr = mat.data(); static_assert(std::is_pointer>::value, "should be getting a pointer"); //static_assert(std::is_const>::value, "should be const"); // If Matrix were written well, above static_assert should pass, and the // following check shouldn't. If it get triggered, it means that a newer // version of Eigen now does have const-correct // MatrixMap. This means that const-correct field maps // are then also possible for matrices BOOST_CHECK(!std::is_const>::value); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(const_scalar_iter_test, F, iter_collections, F) { - using FC_t = typename F::parent::FC_t; + using FC_t = typename F::Parent::FC_t; using ScalarMap = ScalarFieldMap; ScalarMap Smap{F::fc["integer Scalar"]}; for (auto it = Smap.cbegin(); it != Smap.cend(); ++it) { auto&& scal = *it; static_assert(std::is_const>::value, "referred type should be const"); static_assert(std::is_lvalue_reference::value, "Should have returned an lvalue ref"); } } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/test_tensor_algebra.cc b/tests/test_tensor_algebra.cc new file mode 100644 index 0000000..366341e --- /dev/null +++ b/tests/test_tensor_algebra.cc @@ -0,0 +1,46 @@ +/** + * file test_tensor_algebra.cc + * + * @author Till Junge + * + * @date 05 Nov 2017 + * + * @brief Tests for the tensor algebra functions + * + * @section LICENCE + * + * Copyright (C) 2017 Till Junge + * + * µSpectre is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 3, or (at + * your option) any later version. + * + * µSpectre is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Emacs; see the file COPYING. If not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ + +#include "tests.hh" +#include + + +namespace muSpectre { + + BOOST_AUTO_TEST_SUITE(tensor_algebra) + + BOOST_AUTO_TEST_CASE(outer_product_test) { + Eigen::TensorFixedSize> A, B; + A.setRandom(); + B.setRandom(); + } + + BOOST_AUTO_TEST_SUITE_END() + +} // muSpectre