diff --git a/src/common/field.hh b/src/common/field.hh
index 3aaa29f..1643470 100644
--- a/src/common/field.hh
+++ b/src/common/field.hh
@@ -1,971 +1,1028 @@
/**
* @file field.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 07 Sep 2017
*
* @brief header-only implementation of a field for field collections
*
* Copyright © 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 "common/T4_map_proxy.hh"
#include <Eigen/Dense>
#include <string>
#include <sstream>
#include <utility>
#include <typeinfo>
#include <vector>
#include <algorithm>
#include <cmath>
#include <memory>
#include <type_traits>
namespace muSpectre {
/* ---------------------------------------------------------------------- */
/**
* base class for field collection-related exceptions
*/
class FieldCollectionError: public std::runtime_error {
public:
//! constructor
explicit FieldCollectionError(const std::string& what)
:std::runtime_error(what){}
//! constructor
explicit FieldCollectionError(const char * what)
:std::runtime_error(what){}
};
/// base class for field-related exceptions
class FieldError: public FieldCollectionError {
using Parent = FieldCollectionError;
public:
//! constructor
explicit FieldError(const std::string& what)
:Parent(what){}
//! constructor
explicit FieldError(const char * what)
:Parent(what){}
};
/**
* Thrown when a associating a field map to and incompatible field
* is attempted
*/
class FieldInterpretationError: public FieldError
{
public:
//! constructor
explicit FieldInterpretationError(const std::string & what)
:FieldError(what){}
//! constructor
explicit FieldInterpretationError(const char * what)
:FieldError(what){}
};
namespace internal {
/* ---------------------------------------------------------------------- */
/**
- * Virtual base class for all fields. A field represents the per-voxel
- * storage for a scalar, vector or tensor quantity. It is used for type and
- * size checking at runtime and for storage of polymorphic pointers to
- * fully typed and sized fields. FieldBase (and its children) are
- * templated with a specific FieldCollection. A FieldCollection stores
- * multiple fields that all apply to the same set of voxels. Allocating
- * and managing the data for all voxels is handled by the FieldCollection.
- * Note that FieldBase does not know anything about the data type of the
- * field.
+ * Virtual base class for all fields. A field represents meta-information
+ * for the per-pixel storage for a scalar, vector or tensor quantity and
+ * is therefore the abstract class defining the field. It is used for type
+ * and size checking at runtime and for storage of polymorphic pointers to
+ * fully typed and sized fields. `FieldBase` (and its children) are
+ * templated with a specific `FieldCollection`. A `FieldCollection` stores
+ * multiple fields that all apply to the same set of pixels. Allocating
+ * and managing the data for all pixels is handled by the `FieldCollection`.
+ * Note that `FieldBase` does not know anything about about mathematical
+ * operations on the data or how to iterate over all pixels. Mapping the
+ * raw data onto Eigen classes and iterating over those is handled by
+ * the `FieldMap`.
+ * `FieldBase` has the specialisation `TypedFieldBase`.
*/
template <class FieldCollection>
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:
using collection_t = FieldCollection; //!< for type checks
//! Copy constructor
FieldBase(const FieldBase &other) = delete;
//! Move constructor
FieldBase(FieldBase &&other) = delete;
//! Destructor
virtual ~FieldBase() = default;
//! Copy assignment operator
FieldBase& operator=(const FieldBase &other) = delete;
//! Move assignment operator
FieldBase& operator=(FieldBase &&other) = 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;
//! number of pixels in the field
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;
//! give access to collection's base class
friend typename FieldCollection::Parent;
protected:
/* ---------------------------------------------------------------------- */
//! allocate memory etc
virtual void resize(size_t size) = 0;
const std::string name; //!< the field's unique name
const size_t nb_components; //!< number of components per entry
//! reference to the collection this field belongs to
const FieldCollection & collection;
private:
};
/**
* Dummy intermediate class to provide a run-time polymorphic
- * typed field. Mainly for binding python. TypedFieldBase specifies methods
+ * typed field. Mainly for binding Python. TypedFieldBase specifies methods
* that return typed Eigen maps and vectors in addition to pointers to the
* raw data.
+ * `TypedFieldBase` has the specialisation `TypedSizedFieldBase`.
*/
template <class FieldCollection, typename T>
class TypedFieldBase: public FieldBase<FieldCollection>
{
public:
using Parent = FieldBase<FieldCollection>; //!< base class
//! for type checks when mapping this field
using collection_t = typename Parent::collection_t;
using Scalar = T; //!< for type checks
using Base = Parent; //!< for uniformity of interface
//! Plain Eigen type to map
using EigenRep = Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic>;
//! map returned when iterating over field
using EigenMap = Eigen::Map<EigenRep>;
//! Plain eigen vector to map
using EigenVec = Eigen::Map<Eigen::VectorXd>;
//! vector map returned when iterating over field
using EigenVecConst = Eigen::Map<const Eigen::VectorXd>;
//! Default constructor
TypedFieldBase() = delete;
//! constructor
TypedFieldBase(std::string unique_name,
size_t nb_components,
FieldCollection& collection);
//! Copy constructor
TypedFieldBase(const TypedFieldBase &other) = delete;
//! Move constructor
TypedFieldBase(TypedFieldBase &&other) = delete;
//! Destructor
virtual ~TypedFieldBase() = default;
//! Copy assignment operator
TypedFieldBase& operator=(const TypedFieldBase &other) = delete;
//! Move assignment operator
TypedFieldBase& operator=(TypedFieldBase &&other) = delete;
//! return type_id of stored type
virtual const std::type_info & get_stored_typeid() const override final;
virtual size_t size() const override = 0;
//! initialise field to zero (do more complicated initialisations through
//! fully typed maps)
virtual void set_zero() override = 0;
//! raw pointer to content (e.g., for Eigen maps)
virtual T* data() = 0;
//! raw pointer to content (e.g., for Eigen maps)
virtual const T* data() const = 0;
//! return a map representing the entire field as a single `Eigen::Array`
EigenMap eigen();
//! return a map representing the entire field as a single Eigen vector
EigenVec eigenvec();
//! return a map representing the entire field as a single Eigen vector
EigenVecConst eigenvec() const;
protected:
private:
};
/* ---------------------------------------------------------------------- */
//! declaraton for friending
template <class FieldCollection, typename T, Dim_t NbComponents, bool isConst>
class FieldMap;
/* ---------------------------------------------------------------------- */
/**
- * implements Fields that contain a statically known number of
- * scalars of a statically known type per pixel in a
- * `muSpectre::FieldCollection`
+ * A `TypedSizeFieldBase` is the base class for fields that contain
+ * statically known number of scalars of a statically known type per pixel
+ * in a `muSpectre::FieldCollection`. The actual data for all pixels is
+ * stored in `TypedSizeFieldBase::values`.
+ * `TypedSizedFieldBase` has the specialisations `MatrixField` and
+ * `TensorField`.
*/
template <class FieldCollection, typename T, Dim_t NbComponents,
bool ArrayStore=false>
class TypedSizedFieldBase: public TypedFieldBase<FieldCollection, T>
{
friend class FieldMap<FieldCollection, T, NbComponents, true>;
friend class FieldMap<FieldCollection, T, NbComponents, false>;
public:
//! for compatibility checks
constexpr static auto nb_components{NbComponents};
using Parent = TypedFieldBase<FieldCollection, T>; //!< base class
using Scalar = T; //!< for type checking
using Base = typename Parent::Base; //!< root base class
//! type stored if ArrayStore is true
using StoredType = Eigen::Array<T, NbComponents, 1>;
//! storage container
using StorageType = std::conditional_t
<ArrayStore,
std::vector<StoredType,
Eigen::aligned_allocator<StoredType>>,
std::vector<T,Eigen::aligned_allocator<T>>>;
//! Plain type that is being mapped (Eigen lingo)
using EigenRep = Eigen::Array<T, NbComponents, Eigen::Dynamic>;
//! maps returned when iterating over field
using EigenMap = std::conditional_t<
ArrayStore,
Eigen::Map<EigenRep, Eigen::Aligned,
Eigen::OuterStride<sizeof(StoredType)/sizeof(T)>>,
Eigen::Map<EigenRep>>;
//! maps returned when iterating over field
using ConstEigenMap = std::conditional_t<
ArrayStore,
Eigen::Map<const EigenRep, Eigen::Aligned,
Eigen::OuterStride<sizeof(StoredType)/sizeof(T)>>,
Eigen::Map<const EigenRep>>;
//! constructor
TypedSizedFieldBase(std::string unique_name,
FieldCollection& collection);
virtual ~TypedSizedFieldBase() = default;
//! initialise field to zero (do more complicated initialisations through
//! fully typed maps)
inline void set_zero() override final;
//! add a new value at the end of the field
template <bool isArrayStore = ArrayStore>
inline void push_back(const
std::enable_if_t<isArrayStore, StoredType> &
value);
//! add a new value at the end of the field
template <bool componentStore = !ArrayStore>
inline std::enable_if_t<componentStore> push_back(const StoredType & value);
//! Number of stored arrays (i.e. total number of stored
//! scalars/NbComponents)
size_t size() const override final;
/**
* returns an upcasted reference to a field, or throws an
* exception if the field is incompatible
*/
static TypedSizedFieldBase & check_ref(Base & other);
/**
* returns an upcasted reference to a field, or throws an
* exception if the field is incompatible
*/
static const TypedSizedFieldBase & check_ref(const Base & other);
//! return raw pointer to stored data (necessary for Eigen maps)
inline T* data() override final {return this->get_ptr_to_entry(0);}
//! return raw pointer to stored data (necessary for Eigen maps)
inline const T* data() const override final {return this->get_ptr_to_entry(0);}
//! return a map representing the entire field as a single `Eigen::Array`
inline EigenMap eigen();
//! return a map representing the entire field as a single `Eigen::Array`
inline ConstEigenMap eigen() const;
/**
* return a map representing the entire field as a single
* dynamically sized `Eigen::Array` (for python bindings)
*/
inline typename Parent::EigenMap dyn_eigen() {return Parent::eigen();}
//! inner product between compatible fields
template<typename otherT>
inline Real inner_product(const TypedSizedFieldBase<FieldCollection, otherT, NbComponents,
ArrayStore> & other) const;
protected:
//! returns a raw pointer to the entry, for `Eigen::Map`
template <bool isArray=ArrayStore>
inline std::enable_if_t<isArray, T*> get_ptr_to_entry(const size_t&& index);
//! returns a raw pointer to the entry, for `Eigen::Map`
template <bool isArray=ArrayStore>
inline std::enable_if_t<isArray, const T*>
get_ptr_to_entry(const size_t&& index) const;
//! returns a raw pointer to the entry, for `Eigen::Map`
template <bool noArray = !ArrayStore>
inline T* get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index);
//! returns a raw pointer to the entry, for `Eigen::Map`
template <bool noArray = !ArrayStore>
inline const T*
get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index) const;
//! set the storage size of this field
inline virtual void resize(size_t size) override final;
//! The actual storage container
StorageType values{};
};
} // internal
/* ---------------------------------------------------------------------- */
/**
- * subclass of `muSpectre::internal::TypedSizedFieldBase` to
- * represent tensorial fields defined by a stored scalar type,
- * spatial dimension, and tensorial order
+ * The `TensorField` is a subclass of `muSpectre::internal::TypedSizedFieldBase`
+ * that represents tensorial fields, i.e. arbitrary-dimensional arrays with
+ * identical number of rows/columns (that typically correspond to the spatial
+ * cartesian dimensions). It is defined by the stored scalar type @a T, the
+ * tensorial order @a order (often also called degree or rank) and the
+ * number of spatial dimensions @a dim.
*/
template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
class TensorField: public internal::TypedSizedFieldBase<FieldCollection,
T,
ipow(dim,order)>
{
public:
//! base class
using Parent = internal::TypedSizedFieldBase<FieldCollection,
T,
ipow(dim,order)>;
using Base = typename Parent::Base; //!< root base class
//! polymorphic base class
using Field_p = typename FieldCollection::Field_p;
using Scalar = typename Parent::Scalar; //!< for type checking
//! Copy constructor
TensorField(const TensorField &other) = delete;
//! Move constructor
TensorField(TensorField &&other) = delete;
//! Destructor
virtual ~TensorField() = default;
//! Copy assignment operator
TensorField& operator=(const TensorField &other) = delete;
//! Move assignment operator
TensorField& operator=(TensorField &&other) = delete;
//! return the order of the stored tensor
inline Dim_t get_order() const;
//! return the dimension of the stored tensor
inline Dim_t get_dim() const;
//! factory function
template<class FieldType, class CollectionType, typename... Args>
friend FieldType& make_field(std::string unique_name,
CollectionType & collection,
Args&&... args);
//! return a reference or throw an exception if `other` is incompatible
static TensorField & check_ref(Base & other) {
return static_cast<TensorField &>(Parent::check_ref(other));}
//! return a reference or throw an exception if `other` is incompatible
static const TensorField & check_ref(const Base & other) {
return static_cast<const TensorField &>(Parent::check_ref(other));}
/**
- * Pure convenience functions to get a MatrixFieldMap of
- * appropriate dimensions mapped to this field. You can also
- * create other types of maps, as long as they have the right
- * fundamental type (T) and the correct size (nbComponents).
+ * Convenience functions to return a map onto this field. A map allows
+ * iteration over all pixels. The map's iterator returns an object that
+ * represents the underlying mathematical structure of the field and
+ * implements common linear algebra operations on it.
+ * Specifically, this function returns
+ * - A `MatrixFieldMap` with @a dim rows and one column if the tensorial
+ * order @a order is unity.
+ * - A `MatrixFieldMap` with @a dim rows and @a dim columns if the tensorial
+ * order @a order is 2.
+ * - A `T4MatrixFieldMap` if the tensorial order is 4.
*/
decltype(auto) get_map();
/**
- * Pure convenience functions to get a MatrixFieldMap of
- * appropriate dimensions mapped to this field. You can also
- * create other types of maps, as long as they have the right
- * fundamental type (T) and the correct size (nbComponents).
+ * Convenience functions to return a map onto this field. A map allows
+ * iteration over all pixels. The map's iterator returns an object that
+ * represents the underlying mathematical structure of the field and
+ * implements common linear algebra operations on it.
+ * Specifically, this function returns
+ * - A `MatrixFieldMap` with @a dim rows and one column if the tensorial
+ * order @a order is unity.
+ * - A `MatrixFieldMap` with @a dim rows and @a dim columns if the tensorial
+ * order @a order is 2.
+ * - A `T4MatrixFieldMap` if the tensorial order is 4.
*/
decltype(auto) get_const_map();
/**
- * Pure convenience functions to get a MatrixFieldMap of
- * appropriate dimensions mapped to this field. You can also
- * create other types of maps, as long as they have the right
- * fundamental type (T) and the correct size (nbComponents).
+ * Convenience functions to return a map onto this field. A map allows
+ * iteration over all pixels. The map's iterator returns an object that
+ * represents the underlying mathematical structure of the field and
+ * implements common linear algebra operations on it.
+ * Specifically, this function returns
+ * - A `MatrixFieldMap` with @a dim rows and one column if the tensorial
+ * order @a order is unity.
+ * - A `MatrixFieldMap` with @a dim rows and @a dim columns if the tensorial
+ * order @a order is 2.
+ * - A `T4MatrixFieldMap` if the tensorial order is 4.
*/
decltype(auto) get_map() const;
protected:
//! constructor protected!
TensorField(std::string unique_name, FieldCollection & collection);
private:
};
/* ---------------------------------------------------------------------- */
/**
- * subclass of `muSpectre::internal::TypedSizedFieldBase` to
- * represent matrix fields defined by a stored scalar type and the
- * number of rows and columns of the field
+ * The `MatrixField` is subclass of `muSpectre::internal::TypedSizedFieldBase`
+ * that represents matrix fields, i.e. a two dimensional arrays, defined by
+ * the stored scalar type @a T and the number of rows @a NbRow and columns
+ * @a NbCol of the matrix.
*/
template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol=NbRow>
class MatrixField: public internal::TypedSizedFieldBase<FieldCollection,
T,
NbRow*NbCol>
{
public:
//! base class
using Parent = internal::TypedSizedFieldBase<FieldCollection,
T,
NbRow*NbCol>;
using Base = typename Parent::Base; //!< root base class
//! polymorphic base field ptr to store
using Field_p = std::unique_ptr<internal::FieldBase<FieldCollection>>;
//! Copy constructor
MatrixField(const MatrixField &other) = delete;
//! Move constructor
MatrixField(MatrixField &&other) = delete;
//! Destructor
virtual ~MatrixField() = default;
//! Copy assignment operator
MatrixField& operator=(const MatrixField &other) = delete;
//! Move assignment operator
MatrixField& operator=(MatrixField &&other) = delete;
//! returns the number of rows
inline Dim_t get_nb_row() const;
//! returns the number of columns
inline Dim_t get_nb_col() const;
//! factory function
template<class FieldType, class CollectionType, typename... Args>
friend FieldType& make_field(std::string unique_name,
CollectionType & collection,
Args&&... args);
//! returns a `MatrixField` reference if `other` is a compatible field
static MatrixField & check_ref(Base & other) {
return static_cast<MatrixField &>(Parent::check_ref(other));}
//! returns a `MatrixField` reference if `other` is a compatible field
static const MatrixField & check_ref(const Base & other) {
return static_cast<const MatrixField &>(Parent::check_ref(other));}
- //! returns the default map type
+ /**
+ * Convenience functions to return a map onto this field. A map allows
+ * iteration over all pixels. The map's iterator returns an object that
+ * represents the underlying mathematical structure of the field and
+ * implements common linear algebra operations on it.
+ * Specifically, this function returns
+ * - A `ScalarFieldMap` if @a NbRows and @a NbCols are unity.
+ * - A `MatrixFieldMap` with @a NbRows rows and @a NbCols columns
+ * otherwise.
+ */
decltype(auto) get_map();
- //! returns the default map type
+ /**
+ * Convenience functions to return a map onto this field. A map allows
+ * iteration over all pixels. The map's iterator returns an object that
+ * represents the underlying mathematical structure of the field and
+ * implements common linear algebra operations on it.
+ * Specifically, this function returns
+ * - A `ScalarFieldMap` if @a NbRows and @a NbCols are unity.
+ * - A `MatrixFieldMap` with @a NbRows rows and @a NbCols columns
+ * otherwise.
+ */
decltype(auto) get_const_map();
- //! returns the default map type
+ /**
+ * Convenience functions to return a map onto this field. A map allows
+ * iteration over all pixels. The map's iterator returns an object that
+ * represents the underlying mathematical structure of the field and
+ * implements common linear algebra operations on it.
+ * Specifically, this function returns
+ * - A `ScalarFieldMap` if @a NbRows and @a NbCols are unity.
+ * - A `MatrixFieldMap` with @a NbRows rows and @a NbCols columns
+ * otherwise.
+ */
decltype(auto) get_map() const;
protected:
//! constructor protected!
MatrixField(std::string unique_name, FieldCollection & collection);
private:
};
/* ---------------------------------------------------------------------- */
//! convenience alias (
template <class FieldCollection, typename T>
using ScalarField = MatrixField<FieldCollection, T, 1, 1>;
/* ---------------------------------------------------------------------- */
// Implementations
namespace internal {
/* ---------------------------------------------------------------------- */
template <class FieldCollection>
FieldBase<FieldCollection>::FieldBase(std::string unique_name,
size_t nb_components_,
FieldCollection & collection_)
:name(unique_name), nb_components(nb_components_),
collection(collection_) {}
/* ---------------------------------------------------------------------- */
template <class FieldCollection>
inline const std::string & FieldBase<FieldCollection>::get_name() const {
return this->name;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection>
inline const FieldCollection & FieldBase<FieldCollection>::
get_collection() const {
return this->collection;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection>
inline const size_t & FieldBase<FieldCollection>::
get_nb_components() const {
return this->nb_components;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T>
TypedFieldBase<FieldCollection, T>::
TypedFieldBase(std::string unique_name, size_t nb_components,
FieldCollection & collection)
:Parent(unique_name, nb_components, collection)
{}
/* ---------------------------------------------------------------------- */
//! return type_id of stored type
template <class FieldCollection, typename T>
const std::type_info & TypedFieldBase<FieldCollection, T>::
get_stored_typeid() const {
return typeid(T);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T>
typename TypedFieldBase<FieldCollection, T>::EigenMap
TypedFieldBase<FieldCollection, T>::
eigen() {
return EigenMap(this->data(), this->get_nb_components(), this->size());
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T>
typename TypedFieldBase<FieldCollection, T>::EigenVec
TypedFieldBase<FieldCollection, T>::
eigenvec() {
return EigenVec(this->data(), this->get_nb_components() * this->size());
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T>
typename TypedFieldBase<FieldCollection, T>::EigenVecConst
TypedFieldBase<FieldCollection, T>::
eigenvec() const{
return EigenVecConst(this->data(), this->get_nb_components() * this->size());
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
TypedSizedFieldBase(std::string unique_name, FieldCollection & collection)
:Parent(unique_name, NbComponents, collection){
static_assert
((std::is_arithmetic<T>::value ||
std::is_same<T, Complex>::value),
"Use TypedSizedFieldBase for integer, real or complex scalars for T");
static_assert(NbComponents > 0,
"Only fields with more than 0 components");
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
void
TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
set_zero() {
std::fill(this->values.begin(), this->values.end(), T{});
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
size_t TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
size() const {
if (ArrayStore) {
return this->values.size();
} else {
return this->values.size()/NbComponents;
}
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore> &
TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
check_ref(Base & other) {
if (typeid(T).hash_code() != other.get_stored_typeid().hash_code()) {
std::string err ="Cannot create a Reference of requested type " +(
"for field '" + other.get_name() + "' of type '" +
other.get_stored_typeid().name() + "'");
throw std::runtime_error
(err);
}
//check size compatibility
if (NbComponents != other.get_nb_components()) {
throw std::runtime_error
("Cannot create a Reference to a field with " +
std::to_string(NbComponents) + " components " +
"for field '" + other.get_name() + "' with " +
std::to_string(other.get_nb_components()) + " components");
}
return static_cast<TypedSizedFieldBase&>(other);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents,
bool ArrayStore>
const TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore> &
TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
check_ref(const Base & other) {
if (typeid(T).hash_code() != other.get_stored_typeid().hash_code()) {
std::stringstream err_str{};
err_str << "Cannot create a Reference of requested type "
<< "for field '" << other.get_name() << "' of type '"
<< other.get_stored_typeid().name() << "'";
throw std::runtime_error
(err_str.str());
}
//check size compatibility
if (NbComponents != other.get_nb_components()) {
throw std::runtime_error
("Cannot create a Reference to a field with " +
std::to_string(NbComponents) + " components " +
"for field '" + other.get_name() + "' with " +
std::to_string(other.get_nb_components()) + " components");
}
return static_cast<const TypedSizedFieldBase&>(other);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents,
bool ArrayStore>
typename TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::EigenMap
TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
eigen() {
return EigenMap(this->data(), NbComponents, this->size());
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents,
bool ArrayStore>
typename TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::ConstEigenMap
TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
eigen() const {
return ConstEigenMap(this->data(), NbComponents, this->size());
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents,
bool ArrayStore>
template<typename otherT>
Real
TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
inner_product(const TypedSizedFieldBase<FieldCollection, otherT, NbComponents,
ArrayStore> & other) const {
return (this->eigen() * other.eigen()).sum();
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool isArray>
std::enable_if_t<isArray, T*>
TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
get_ptr_to_entry(const size_t&& index) {
static_assert (isArray == ArrayStore, "SFINAE");
return &this->values[std::move(index)](0, 0);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool noArray>
T* TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index) {
static_assert (noArray != ArrayStore, "SFINAE");
return &this->values[NbComponents*std::move(index)];
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool isArray>
std::enable_if_t<isArray, const T*>
TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
get_ptr_to_entry(const size_t&& index) const {
static_assert (isArray == ArrayStore, "SFINAE");
return &this->values[std::move(index)](0, 0);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool noArray>
const T* TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index) const {
static_assert (noArray != ArrayStore, "SFINAE");
return &this->values[NbComponents*std::move(index)];
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
void TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
resize(size_t size) {
if (ArrayStore) {
this->values.resize(size);
} else {
this->values.resize(size*NbComponents);
}
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool isArrayStore>
void
TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
push_back(const std::enable_if_t<isArrayStore,StoredType> & value) {
static_assert(isArrayStore == ArrayStore, "SFINAE");
this->values.push_back(value);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool componentStore>
std::enable_if_t<componentStore>
TypedSizedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
push_back(const StoredType & value) {
static_assert(componentStore != ArrayStore, "SFINAE");
for (Dim_t i = 0; i < NbComponents; ++i) {
this->values.push_back(value(i));
}
}
} // internal
/* ---------------------------------------------------------------------- */
//! Factory function, guarantees that only fields get created that are
//! properly registered and linked to a collection.
template<class FieldType, class FieldCollection, typename... Args>
FieldType &
make_field(std::string unique_name,
FieldCollection & collection,
Args&&... args) {
auto ptr = std::unique_ptr<FieldType>{
new FieldType(unique_name, collection, args...)};
auto& retref = *ptr;
collection.register_field(std::move(ptr));
return retref;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
TensorField<FieldCollection, T, order, dim>::
TensorField(std::string unique_name, FieldCollection & collection)
:Parent(unique_name, collection) {}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
Dim_t TensorField<FieldCollection, T, order, dim>::
get_order() const {
return order;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
Dim_t TensorField<FieldCollection, T, order, dim>::
get_dim() const {
return dim;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
MatrixField<FieldCollection, T, NbRow, NbCol>::
MatrixField(std::string unique_name, FieldCollection & collection)
:Parent(unique_name, collection) {}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
Dim_t MatrixField<FieldCollection, T, NbRow, NbCol>::
get_nb_col() const {
return NbCol;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
Dim_t MatrixField<FieldCollection, T, NbRow, NbCol>::
get_nb_row() const {
return NbRow;
}
} // muSpectre
#include "common/field_map.hh"
namespace muSpectre {
namespace internal {
/* ---------------------------------------------------------------------- */
/**
* defines the default mapped type obtained when calling
* `muSpectre::TensorField::get_map()`
*/
template <class FieldCollection, typename T, size_t order, Dim_t dim,
bool ConstMap>
struct tensor_map_type {
};
/// specialisation for vectors
template <class FieldCollection, typename T, Dim_t dim, bool ConstMap>
struct tensor_map_type<FieldCollection, T, firstOrder, dim, ConstMap> {
//! use this type
using type = MatrixFieldMap<FieldCollection, T, dim, 1, ConstMap>;
};
/// specialisation to second-order tensors (matrices)
template <class FieldCollection, typename T, Dim_t dim, bool ConstMap>
struct tensor_map_type<FieldCollection, T, secondOrder, dim, ConstMap> {
//! use this type
using type = MatrixFieldMap<FieldCollection, T, dim, dim, ConstMap>;
};
/// specialisation to fourth-order tensors
template <class FieldCollection, typename T, Dim_t dim, bool ConstMap>
struct tensor_map_type<FieldCollection, T, fourthOrder, dim, ConstMap> {
//! use this type
using type = T4MatrixFieldMap<FieldCollection, T, dim, ConstMap>;
};
/* ---------------------------------------------------------------------- */
/**
* defines the default mapped type obtained when calling
* `muSpectre::MatrixField::get_map()`
*/
template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol,
bool ConstMap>
struct matrix_map_type {
//! mapping type
using type =
MatrixFieldMap<FieldCollection, T, NbRow, NbCol, ConstMap>;
};
- //! spectialisation to scalar fields
+ //! specialisation to scalar fields
template <class FieldCollection, typename T, bool ConstMap>
struct matrix_map_type<FieldCollection, T, oneD, oneD, ConstMap> {
//! mapping type
using type = ScalarFieldMap<FieldCollection, T, ConstMap>;
};
} // internal
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
decltype(auto) TensorField<FieldCollection, T, order, dim>::
get_map() {
constexpr bool map_constness{false};
using RawMap_t =
typename internal::tensor_map_type<FieldCollection, T, order, dim,
map_constness>::type;
return RawMap_t(*this);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
decltype(auto) TensorField<FieldCollection, T, order, dim>::
get_const_map() {
constexpr bool map_constness{true};
using RawMap_t =
typename internal::tensor_map_type<FieldCollection, T, order, dim,
map_constness>::type;
return RawMap_t(*this);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
decltype(auto) TensorField<FieldCollection, T, order, dim>::
get_map() const {
constexpr bool map_constness{true};
using RawMap_t =
typename internal::tensor_map_type<FieldCollection, T, order, dim,
map_constness>::type;
return RawMap_t(*this);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
decltype(auto) MatrixField<FieldCollection, T, NbRow, NbCol>::
get_map() {
constexpr bool map_constness{false};
using RawMap_t =
typename internal::matrix_map_type<FieldCollection, T, NbRow, NbCol,
map_constness>::type;
return RawMap_t(*this);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
decltype(auto) MatrixField<FieldCollection, T, NbRow, NbCol>::
get_const_map() {
constexpr bool map_constness{true};
using RawMap_t =
typename internal::matrix_map_type<FieldCollection, T, NbRow, NbCol,
map_constness>::type;
return RawMap_t(*this);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol>
decltype(auto) MatrixField<FieldCollection, T, NbRow, NbCol>::
get_map() const {
constexpr bool map_constness{true};
using RawMap_t =
typename internal::matrix_map_type<FieldCollection, T, NbRow, NbCol,
map_constness>::type;
return RawMap_t(*this);
}
} // muSpectre
#endif /* FIELD_H */
diff --git a/src/common/field_collection_base.hh b/src/common/field_collection_base.hh
index 7c59c61..d7d13c1 100644
--- a/src/common/field_collection_base.hh
+++ b/src/common/field_collection_base.hh
@@ -1,173 +1,184 @@
/**
* @file field_collection_base.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 05 Nov 2017
*
* @brief Base class for field collections
*
* Copyright © 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_BASE_H
#define FIELD_COLLECTION_BASE_H
#include "common/common.hh"
#include "common/field.hh"
#include <map>
#include <vector>
namespace muSpectre {
/* ---------------------------------------------------------------------- */
- //! DimS spatial dimension (dimension of problem)
- //! DimM material_dimension (dimension of constitutive law)
+ /** `FieldCollectionBase` is the base class for collections of fields. All
+ * fields in a field collection have the same number of pixels. The field
+ * collection is templated with @a DimS is the spatial dimension (i.e.
+ * whether the simulation domain is one, two or three-dimensional) and
+ * @a DimM is the material dimension (i.e., the dimension of constitutive
+ * law; even for e.g. two-dimensional problems the constitutive law could
+ * live in three-dimensional space for e.g. plane strain or stress problems).
+ * All fields within a field collection have a unique string identifier.
+ * A `FieldCollectionBase` is therefore comparable to a dictionary of fields
+ * that live on the same grid.
+ * `FieldCollectionBase` has the specialisations `GlobalFieldCollection` and
+ * `LocalFieldCollection`.
+ */
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
class FieldCollectionBase
{
public:
//! polymorphic base type to store
using Field = internal::FieldBase<FieldCollectionDerived>;
using Field_p = std::unique_ptr<Field>; //!< stored type
using Ccoord = Ccoord_t<DimS>; //!< cell coordinates type
//! Default constructor
FieldCollectionBase();
//! Copy constructor
FieldCollectionBase(const FieldCollectionBase &other) = delete;
//! Move constructor
FieldCollectionBase(FieldCollectionBase &&other) = delete;
//! Destructor
virtual ~FieldCollectionBase() = default;
//! Copy assignment operator
FieldCollectionBase& operator=(const FieldCollectionBase &other) = delete;
//! Move assignment operator
FieldCollectionBase& operator=(FieldCollectionBase &&other) = delete;
//! Register a new field (fields need to be in heap, so I want to keep them
//! as 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_;}
//! check whether a field is present
bool check_field_exists(std::string unique_name);
protected:
std::map<const std::string, Field_p> fields{}; //!< contains the field ptrs
bool is_initialised{false}; //!< to handle double initialisation correctly
const Uint id; //!< unique identifier
static Uint counter; //!< used to assign unique identifiers
size_t size_{0}; //!< holds the number of pixels after initialisation
private:
};
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
Uint FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::counter{0};
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::FieldCollectionBase()
:id(counter++){}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
void FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::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());
}
if (this->is_initialised) {
field->resize(this->size());
}
this->fields[field->get_name()] = std::move(field);
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
constexpr Dim_t FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::
spatial_dim() {
return DimS;
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
Dim_t FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::
get_spatial_dim() const {
return DimS;
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
typename FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::Field&
FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::
operator[](std::string unique_name) {
return *(this->fields[unique_name]);
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
typename FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::Field&
FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::
at(std::string unique_name) {
return *(this->fields.at(unique_name));
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
bool
FieldCollectionBase<DimS, DimM, FieldCollectionDerived>::
check_field_exists(std::string unique_name) {
return this->fields.find(unique_name) != this->fields.end();
}
} // muSpectre
#endif /* FIELD_COLLECTION_BASE_H */
diff --git a/src/common/field_collection_global.hh b/src/common/field_collection_global.hh
index 2af5958..98eb516 100644
--- a/src/common/field_collection_global.hh
+++ b/src/common/field_collection_global.hh
@@ -1,191 +1,193 @@
/**
* @file field_collection_global.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 05 Nov 2017
*
* @brief FieldCollection base-class for global fields
*
* Copyright © 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_GLOBAL_H
#define FIELD_COLLECTION_GLOBAL_H
#include "common/field_collection_base.hh"
#include "common/ccoord_operations.hh"
namespace muSpectre {
/* ---------------------------------------------------------------------- */
- //! DimS spatial dimension (dimension of problem)
- //! DimM material_dimension (dimension of constitutive law)
+ /** `GlobalFieldCollection` derives from `FieldCollectionBase` and stores
+ * global fields that live throughout the whole computational domain, i.e.
+ * are defined for each pixel.
+ */
template<Dim_t DimS, Dim_t DimM>
class GlobalFieldCollection:
public FieldCollectionBase<DimS, DimM, GlobalFieldCollection<DimS, DimM>>
{
public:
using Parent = FieldCollectionBase
<DimS, DimM, GlobalFieldCollection<DimS, DimM>>; //!< base class
using Ccoord = typename Parent::Ccoord; //!< cell coordinates type
using Field_p = typename Parent::Field_p; //!< spatial coordinates type
//! iterator over all pixels contained it the collection
using iterator = typename CcoordOps::Pixels<DimS>::iterator;
//! Default constructor
GlobalFieldCollection();
//! Copy constructor
GlobalFieldCollection(const GlobalFieldCollection &other) = delete;
//! Move constructor
GlobalFieldCollection
(GlobalFieldCollection &&other) = default;
//! Destructor
virtual ~GlobalFieldCollection() = default;
//! Copy assignment operator
GlobalFieldCollection&
operator=(const GlobalFieldCollection &other) = delete;
//! Move assignment operator
GlobalFieldCollection&
operator=(GlobalFieldCollection &&other) = default;
/** allocate memory, etc. At this point, the collection is
informed aboud the size and shape of the domain (through the
sizes parameter). The job of initialise is to make sure that
all fields are either of size 0, in which case they need to be
allocated, or are of the same size as the product of 'sizes'
(if standard strides apply) any field of a different size is
wrong.
TODO: check whether it makes sense to put a runtime check here
**/
inline void initialise(Ccoord sizes);
//! return the pixel sizes
inline const Ccoord & get_sizes() const;
//! returns the linear index corresponding to cell coordinates
template <class CcoordRef>
inline size_t get_index(CcoordRef && ccoord) const;
//! returns the cell coordinates corresponding to a linear index
inline Ccoord get_ccoord(size_t index) const;
inline iterator begin(); //!< returns iterator to first pixel
inline iterator end(); //!< returns iterator past the last pixel
//! return spatial dimension (template parameter)
static constexpr inline Dim_t spatial_dim() {return DimS;}
//! return material dimension (template parameter)
static constexpr inline Dim_t material_dim() {return DimM;}
protected:
//! number of discretisation cells in each of the DimS spatial directions
Ccoord sizes{};
CcoordOps::Pixels<DimS> pixels{}; //!< helper to iterate over the grid
private:
};
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
GlobalFieldCollection<DimS, DimM>::GlobalFieldCollection()
:Parent()
{}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void GlobalFieldCollection<DimS, DimM>::
initialise(Ccoord sizes) {
if (this->is_initialised) {
throw std::runtime_error("double initialisation");
}
this->pixels = CcoordOps::Pixels<DimS>(sizes);
this->size_ = CcoordOps::get_size(sizes);
this->sizes = sizes;
std::for_each(std::begin(this->fields), std::end(this->fields),
[this](auto && item) {
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;
}
//----------------------------------------------------------------------------//
//! return the pixel sizes
template <Dim_t DimS, Dim_t DimM>
const typename GlobalFieldCollection<DimS, DimM>::Ccoord &
GlobalFieldCollection<DimS, DimM>::get_sizes() const {
return this->sizes;
}
//----------------------------------------------------------------------------//
//! returns the cell coordinates corresponding to a linear index
template <Dim_t DimS, Dim_t DimM>
typename GlobalFieldCollection<DimS, DimM>::Ccoord
GlobalFieldCollection<DimS, DimM>::get_ccoord(size_t index) const {
return CcoordOps::get_ccoord(this->get_sizes(), std::move(index));
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
typename GlobalFieldCollection<DimS, DimM>::iterator
GlobalFieldCollection<DimS, DimM>::begin() {
return this->pixels.begin();
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
typename GlobalFieldCollection<DimS, DimM>::iterator
GlobalFieldCollection<DimS, DimM>::end() {
return this->pixels.end();
}
//----------------------------------------------------------------------------//
//! returns the linear index corresponding to cell coordinates
template <Dim_t DimS, Dim_t DimM>
template <class CcoordRef>
size_t
GlobalFieldCollection<DimS, DimM>::get_index(CcoordRef && ccoord) const {
static_assert(std::is_same<
Ccoord,
std::remove_const_t<
std::remove_reference_t<CcoordRef>>>::value,
"can only be called with values or references of Ccoord");
return CcoordOps::get_index(this->get_sizes(),
std::forward<CcoordRef>(ccoord));
}
} // muSpectre
#endif /* FIELD_COLLECTION_GLOBAL_H */
diff --git a/src/common/field_collection_local.hh b/src/common/field_collection_local.hh
index 3eff772..2de7864 100644
--- a/src/common/field_collection_local.hh
+++ b/src/common/field_collection_local.hh
@@ -1,176 +1,181 @@
/**
* @file field_collection_local.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 05 Nov 2017
*
* @brief FieldCollection base-class for local fields
*
* Copyright © 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_LOCAL_H
#define FIELD_COLLECTION_LOCAL_H
#include "common/field_collection_base.hh"
namespace muSpectre {
/* ---------------------------------------------------------------------- */
- //! DimS spatial dimension (dimension of problem)
- //! DimM material_dimension (dimension of constitutive law)
+ /** `LocalFieldCollection` derives from `FieldCollectionBase` and stores
+ * local fields, i.e. fields that are only defined for a subset of all pixels
+ * in the computational domain. The coordinates of these active pixels are
+ * explicitly stored by this field collection.
+ * `LocalFieldCollection::add_pixel` allows to add individual pixels to the
+ * field collection.
+ */
template<Dim_t DimS, Dim_t DimM>
class LocalFieldCollection:
public FieldCollectionBase<DimS, DimM, LocalFieldCollection<DimS, DimM>>
{
public:
//! base class
using Parent = FieldCollectionBase<DimS, DimM,
LocalFieldCollection<DimS, DimM>>;
using Ccoord = typename Parent::Ccoord; //!< cell coordinates type
using Field_p = typename Parent::Field_p; //!< field pointer
using ccoords_container = std::vector<Ccoord>; //!< list of pixels
//! iterator over managed pixels
using iterator = typename ccoords_container::iterator;
//! Default constructor
LocalFieldCollection();
//! Copy constructor
LocalFieldCollection(const LocalFieldCollection &other) = delete;
//! Move constructor
LocalFieldCollection(LocalFieldCollection &&other) = delete;
//! Destructor
virtual ~LocalFieldCollection() = default;
//! Copy assignment operator
LocalFieldCollection& operator=(const LocalFieldCollection &other) = delete;
//! Move assignment operator
LocalFieldCollection& operator=(LocalFieldCollection &&other) = delete;
//! add a pixel/voxel to the field collection
inline void add_pixel(const Ccoord & local_ccoord);
/** allocate memory, etc. at this point, the field collection
knows how many entries it should have from the size of the
coords containes (which grows by one every time add_pixel is
called. The job of initialise is to make sure that all fields
are either of size 0, in which case they need to be allocated,
or are of the same size as the product of 'sizes' any field of
a different size is wrong TODO: check whether it makes sense
to put a runtime check here
**/
inline void initialise();
//! returns the linear index corresponding to cell coordinates
template <class CcoordRef>
inline size_t get_index(CcoordRef && ccoord) const;
//! returns the cell coordinates corresponding to a linear index
inline Ccoord get_ccoord(size_t index) const;
//! iterator to first pixel
inline iterator begin() {return this->ccoords.begin();}
//! iterator past last pixel
inline iterator end() {return this->ccoords.end();}
protected:
//! container of pixel coords for non-global collections
ccoords_container ccoords{};
//! container of indices for non-global collections (slow!)
std::map<Ccoord, size_t> indices{};
private:
};
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
LocalFieldCollection<DimS, DimM>::LocalFieldCollection()
:Parent(){}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void LocalFieldCollection<DimS, DimM>::
add_pixel(const Ccoord & local_ccoord) {
if (this->is_initialised) {
throw FieldCollectionError
("once a field collection has been initialised, you can't add new "
"pixels.");
}
this->indices[local_ccoord] = this->ccoords.size();
this->ccoords.push_back(local_ccoord);
this->size_++;
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void LocalFieldCollection<DimS, DimM>::
initialise() {
if (this->is_initialised) {
throw std::runtime_error("double initialisation");
}
std::for_each(std::begin(this->fields), std::end(this->fields),
[this](auto && item) {
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;
}
//----------------------------------------------------------------------------//
//! returns the linear index corresponding to cell coordinates
template <Dim_t DimS, Dim_t DimM>
template <class CcoordRef>
size_t
LocalFieldCollection<DimS, DimM>::get_index(CcoordRef && ccoord) const {
static_assert(std::is_same<
Ccoord,
std::remove_const_t<
std::remove_reference_t<CcoordRef>>>::value,
"can only be called with values or references of Ccoord");
return this->indices.at(std::forward<CcoordRef>(ccoord));
}
//----------------------------------------------------------------------------//
//! returns the cell coordinates corresponding to a linear index
template <Dim_t DimS, Dim_t DimM>
typename LocalFieldCollection<DimS, DimM>::Ccoord
LocalFieldCollection<DimS, DimM>::get_ccoord(size_t index) const {
return this->ccoords[std::move(index)];
}
} // muSpectre
#endif /* FIELD_COLLECTION_LOCAL_H */
diff --git a/src/common/field_map_base.hh b/src/common/field_map_base.hh
index 891fbfb..827c11c 100644
--- a/src/common/field_map_base.hh
+++ b/src/common/field_map_base.hh
@@ -1,644 +1,654 @@
/**
* @file field_map.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 12 Sep 2017
*
* @brief Defined a strongly defines proxy that iterates efficiently over a field
*
* Copyright © 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 "common/field.hh"
#include "field_collection_base.hh"
#include "common/common.hh"
#include <unsupported/Eigen/CXX11/Tensor>
#include <array>
#include <string>
#include <memory>
#include <type_traits>
namespace muSpectre {
namespace internal {
//----------------------------------------------------------------------------//
//! little helper to automate creation of const maps without duplication
template<class T, bool isConst>
struct const_corrector {
//! non-const type
using type = typename T::reference;
};
//! specialisation for constant case
template<class T>
struct const_corrector<T, true> {
//! const type
using type = typename T::const_reference;
};
//! convenience alias
template<class T, bool isConst>
using const_corrector_t = typename const_corrector<T, isConst>::type;
//----------------------------------------------------------------------------//
+ /**
+ * `FieldMap` provides two central mechanisms:
+ * - Map a field (that knows only about the size of the underlying object,
+ * onto the mathematical object (reprensented by the respective Eigen class)
+ * that provides linear algebra functionality.
+ * - Provide an iterator that allows to iterate over all pixels.
+ * A field is represented by `FieldBase` or a derived class.
+ * `FieldMap` has the specialisations `MatrixLikeFieldMap`,
+ * `ScalarFieldMap` and `TensorFieldMap`.
+ */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
class FieldMap
{
public:
//! number of scalars per entry
constexpr static auto nb_components{NbComponents};
using TypedField_nc = TypedSizedFieldBase
<FieldCollection, T, NbComponents>; //!< non-constant version of field
//! field type as seen from iterator
using TypedField = std::conditional_t<ConstField,
const TypedField_nc,
TypedField_nc>;
using Field = typename TypedField::Base; //!< iterated field type
using size_type = std::size_t; //!< stl conformance
using pointer = std::conditional_t<ConstField,
const T*,
T*>; //!< stl conformance
//! Default constructor
FieldMap() = delete;
//! constructor
template <bool isntConst=!ConstField>
FieldMap(std::enable_if_t<isntConst, Field &> field);
//! constructor
template <bool isConst=ConstField>
FieldMap(std::enable_if_t<isConst, const Field &> field);
//! constructor with run-time cost (for python and debugging)
template<class FC, typename T2, Dim_t NbC>
FieldMap(TypedSizedFieldBase<FC, T2, NbC> & field);
//! Copy constructor
FieldMap(const FieldMap &other) = default;
//! Move constructor
FieldMap(FieldMap &&other) = default;
//! Destructor
virtual ~FieldMap() = default;
//! Copy assignment operator
FieldMap& operator=(const FieldMap &other) = delete;
//! Move assignment operator
FieldMap& operator=(FieldMap &&other) = delete;
//! 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;
//! compile-time compatibility check
template<class TypedField>
struct is_compatible;
/**
* iterates over all pixels in the `muSpectre::FieldCollection`
* and dereferences to an Eigen map to the currently used field.
*/
template <class FullyTypedFieldMap, bool ConstIter=false>
class iterator
{
static_assert(!((ConstIter==false) && (ConstField==true)),
"You can't have a non-const iterator over a const "
"field");
public:
//! stl conformance
using value_type =
const_corrector_t<FullyTypedFieldMap, ConstIter>;
//! stl conformance
using const_value_type =
const_corrector_t<FullyTypedFieldMap, true>;
//! stl conformance
using pointer = typename FullyTypedFieldMap::pointer;
//! stl conformance
using difference_type = std::ptrdiff_t;
//! stl conformance
using iterator_category = std::random_access_iterator_tag;
//! cell coordinates type
using Ccoord = typename FieldCollection::Ccoord;
//! stl conformance
using reference = typename FullyTypedFieldMap::reference;
//! fully typed reference as seen by the iterator
using TypedRef = std::conditional_t<ConstIter,
const FullyTypedFieldMap &,
FullyTypedFieldMap>;
//! Default constructor
iterator() = delete;
//! constructor
inline iterator(TypedRef fieldmap, bool begin=true);
//! constructor for random access
inline iterator(TypedRef fieldmap, size_t index);
//! Copy constructor
iterator(const iterator &other)= default;
//! Move constructor
iterator(iterator &&other) = default;
//! Destructor
virtual ~iterator() = default;
//! Copy assignment operator
iterator& operator=(const iterator &other) = default;
//! Move assignment operator
iterator& operator=(iterator &&other) = default;
//! pre-increment
inline iterator & operator++();
//! post-increment
inline iterator operator++(int);
//! dereference
inline value_type operator*();
//! dereference
inline const_value_type operator*() const;
//! 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);
//! access subscripting
inline const_value_type operator[](const difference_type diff) const;
//! 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;
//! div. comparisons
inline bool operator<=(const iterator & other) const;
//! div. comparisons
inline bool operator>(const iterator & other) const;
//! div. comparisons
inline bool operator>=(const iterator & other) const;
//! additions, subtractions and corresponding assignments
inline iterator operator+(difference_type diff) const;
//! additions, subtractions and corresponding assignments
inline iterator operator-(difference_type diff) const;
//! additions, subtractions and corresponding assignments
inline iterator& operator+=(difference_type diff);
//! additions, subtractions and corresponding assignments
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 (ConstIter) {
os << "const ";
}
os << "iterator on field '"
<< it.fieldmap.get_name()
<< "', entry " << it.index;
return os;
}
protected:
const FieldCollection & collection; //!< collection of the field
TypedRef fieldmap; //!< ref to the field itself
size_t index; //!< index of currently pointed-to pixel
private:
};
protected:
//! raw pointer to entry (for Eigen Map)
inline pointer get_ptr_to_entry(size_t index);
//! raw pointer to entry (for Eigen Map)
inline const T* get_ptr_to_entry(size_t index) const;
const FieldCollection & collection; //!< collection holding Field
TypedField & field; //!< mapped Field
private:
};
} // internal
namespace internal {
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template <bool isntConst>
FieldMap<FieldCollection, T, NbComponents, ConstField>::
FieldMap(std::enable_if_t<isntConst, Field &> field)
:collection(field.get_collection()), field(static_cast<TypedField&>(field)) {
static_assert(NbComponents>0,
"Only fields with more than 0 components allowed");
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template <bool isConst>
FieldMap<FieldCollection, T, NbComponents, ConstField>::
FieldMap(std::enable_if_t<isConst, const Field &> field)
:collection(field.get_collection()), field(static_cast<const TypedField&>(field)) {
static_assert(NbComponents>0,
"Only fields with more than 0 components allowed");
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template <class FC, typename T2, Dim_t NbC>
FieldMap<FieldCollection, T, NbComponents, ConstField>::
FieldMap(TypedSizedFieldBase<FC, T2, NbC> & field)
:collection(field.get_collection()), field(static_cast<TypedField&>(field)) {
static_assert(std::is_same<FC, FieldCollection>::value,
"The field does not have the expected FieldCollection type");
static_assert(std::is_same<T2, T>::value,
"The field does not have the expected Scalar type");
static_assert((NbC == NbComponents),
"The field does not have the expected number of components");
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
void
FieldMap<FieldCollection, T, NbComponents, ConstField>::
check_compatibility() {
if (typeid(T).hash_code() !=
this->field.get_stored_typeid().hash_code()) {
std::string err{"Cannot create a Map of type '" +
this->info_string() +
"' for field '" + this->field.get_name() + "' of type '" +
this->field.get_stored_typeid().name() + "'"};
throw FieldInterpretationError
(err);
}
//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<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
size_t
FieldMap<FieldCollection, T, NbComponents, ConstField>::
size() const {
return this->collection.size();
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template <class myField>
struct FieldMap<FieldCollection, T, NbComponents, ConstField>::is_compatible {
//! creates a more readable compile error
constexpr static bool explain() {
static_assert
(std::is_same<typename myField::collection_t, FieldCollection>::value,
"The field does not have the expected FieldCollection type");
static_assert
(std::is_same<typename myField::Scalar, T>::value,
"The // field does not have the expected Scalar type");
static_assert((TypedField::nb_components == NbComponents),
"The field does not have the expected number of components");
//The static asserts wouldn't pass in the incompatible case, so this is it
return true;
}
//! evaluated compatibility
constexpr static bool value{std::is_base_of<TypedField, myField>::value};
};
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
const std::string &
FieldMap<FieldCollection, T, NbComponents, ConstField>::
get_name() const {
return this->field.get_name();
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
const FieldCollection &
FieldMap<FieldCollection, T, NbComponents, ConstField>::
get_collection() const {
return this->collection;
}
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
// Iterator implementations
//! constructor
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
FieldMap<FieldCollection, T, NbComponents, ConstField>::iterator
<FullyTypedFieldMap, ConstIter>::
iterator(TypedRef fieldmap, bool begin)
:collection(fieldmap.get_collection()), fieldmap(fieldmap),
index(begin ? 0 : fieldmap.field.size()) {}
/* ---------------------------------------------------------------------- */
//! constructor for random access
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
FieldMap<FieldCollection, T, NbComponents, ConstField>::iterator<FullyTypedFieldMap, ConstIter>::
iterator(TypedRef fieldmap, size_t index)
:collection(fieldmap.collection), fieldmap(fieldmap),
index(index) {}
/* ---------------------------------------------------------------------- */
//! pre-increment
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template
iterator<FullyTypedFieldMap, ConstIter> &
FieldMap<FieldCollection, T, NbComponents, ConstField>::iterator<FullyTypedFieldMap, ConstIter>::
operator++() {
this->index++;
return *this;
}
/* ---------------------------------------------------------------------- */
//! post-increment
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template
iterator<FullyTypedFieldMap, ConstIter>
FieldMap<FieldCollection, T, NbComponents, ConstField>::iterator<FullyTypedFieldMap, ConstIter>::
operator++(int) {
iterator current = *this;
this->index++;
return current;
}
/* ---------------------------------------------------------------------- */
//! dereference
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::
template iterator<FullyTypedFieldMap, ConstIter>::value_type
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator*() {
return this->fieldmap.operator[](this->index);
}
/* ---------------------------------------------------------------------- */
//! dereference
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::
template iterator<FullyTypedFieldMap, ConstIter>::const_value_type
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator*() const {
return this->fieldmap.operator[](this->index);
}
/* ---------------------------------------------------------------------- */
//! member of pointer
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FullyTypedFieldMap::pointer
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator->() {
return this->fieldmap.ptr_to_val_t(this->index);
}
/* ---------------------------------------------------------------------- */
//! pre-decrement
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template
iterator<FullyTypedFieldMap, ConstIter> &
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator--() {
this->index--;
return *this;
}
/* ---------------------------------------------------------------------- */
//! post-decrement
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template
iterator<FullyTypedFieldMap, ConstIter>
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator--(int) {
iterator current = *this;
this->index--;
return current;
}
/* ---------------------------------------------------------------------- */
//! Access subscripting
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::
template iterator<FullyTypedFieldMap, ConstIter>::value_type
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator[](difference_type diff) {
return this->fieldmap[this->index+diff];
}
/* ---------------------------------------------------------------------- */
//! Access subscripting
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template iterator<FullyTypedFieldMap, ConstIter>::const_value_type
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator[](const difference_type diff) const {
return this->fieldmap[this->index+diff];
}
/* ---------------------------------------------------------------------- */
//! equality
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
bool
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator==(const iterator & other) const {
return (this->index == other.index);
}
/* ---------------------------------------------------------------------- */
//! inquality
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
bool
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator!=(const iterator & other) const {
return !(*this == other);
}
/* ---------------------------------------------------------------------- */
//! div. comparisons
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
bool
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator<(const iterator & other) const {
return (this->index < other.index);
}
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
bool
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator<=(const iterator & other) const {
return (this->index <= other.index);
}
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
bool
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator>(const iterator & other) const {
return (this->index > other.index);
}
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
bool
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator>=(const iterator & other) const {
return (this->index >= other.index);
}
/* ---------------------------------------------------------------------- */
//! additions, subtractions and corresponding assignments
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template
iterator<FullyTypedFieldMap, ConstIter>
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator+(difference_type diff) const {
return iterator(this->fieldmap, this->index + diff);
}
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template
iterator<FullyTypedFieldMap, ConstIter>
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator-(difference_type diff) const {
return iterator(this->fieldmap, this->index - diff);
}
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template
iterator<FullyTypedFieldMap, ConstIter> &
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator+=(difference_type diff) {
this->index += diff;
return *this;
}
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::template
iterator<FullyTypedFieldMap, ConstIter> &
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
operator-=(difference_type diff) {
this->index -= diff;
return *this;
}
/* ---------------------------------------------------------------------- */
//! get pixel coordinates
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
template<class FullyTypedFieldMap, bool ConstIter>
typename FieldCollection::Ccoord
FieldMap<FieldCollection, T, NbComponents, ConstField>::
iterator<FullyTypedFieldMap, ConstIter>::
get_ccoord() const {
return this->collection.get_ccoord(this->index);
}
////----------------------------------------------------------------------------//
//template<class FieldCollection, typename T, Dim_t NbComponents, class FullyTypedFieldMap>
//std::ostream & operator <<
//(std::ostream &os,
// const typename FieldMap<FieldCollection, T, NbComponents, ConstField>::
// template iterator<FullyTypedFieldMap, ConstIter> & it) {
// os << "iterator on field '"
// << it.field.get_name()
// << "', entry " << it.index;
// return os;
//}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
typename FieldMap<FieldCollection, T, NbComponents, ConstField>::pointer
FieldMap<FieldCollection, T, NbComponents, ConstField>::
get_ptr_to_entry(size_t index) {
return this->field.get_ptr_to_entry(std::move(index));
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t NbComponents, bool ConstField>
const T*
FieldMap<FieldCollection, T, NbComponents, ConstField>::
get_ptr_to_entry(size_t index) const {
return this->field.get_ptr_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 332bc3f..ac4f67e 100644
--- a/src/common/field_map_matrixlike.hh
+++ b/src/common/field_map_matrixlike.hh
@@ -1,402 +1,404 @@
/**
* @file field_map_matrixlike.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 26 Sep 2017
*
* @brief Eigen-Matrix and -Array maps over strongly typed fields
*
* Copyright © 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 "common/field_map_base.hh"
#include "common/T4_map_proxy.hh"
#include <Eigen/Dense>
#include <memory>
namespace muSpectre {
namespace internal {
/* ---------------------------------------------------------------------- */
/**
* lists all matrix-like types consideres by
* `muSpectre::internal::MatrixLikeFieldMap`
*/
enum class Map_t{
Matrix, //!< for wrapping `Eigen::Matrix`
Array, //!< for wrapping `Eigen::Array`
T4Matrix //!< for wrapping `Eigen::T4Matrix`
};
/**
* traits structure to define the name shown when a
* `muSpectre::MatrixLikeFieldMap` output into an ostream
*/
template<Map_t map_type>
struct NameWrapper {
};
/// specialisation for `muSpectre::ArrayFieldMap`
template<>
struct NameWrapper<Map_t::Array> {
//! string to use for printing
static std::string field_info_root() {return "Array";}
};
/// specialisation for `muSpectre::MatrixFieldMap`
template<>
struct NameWrapper<Map_t::Matrix> {
//! string to use for printing
static std::string field_info_root() {return "Matrix";}
};
/// specialisation for `muSpectre::T4MatrixFieldMap`
template<>
struct NameWrapper<Map_t::T4Matrix> {
//! string to use for printing
static std::string field_info_root() {return "T4Matrix";}
};
/* ---------------------------------------------------------------------- */
/*!
- * base class for maps of matrices, arrays and fourth-order
- * tensors mapped onty matrices
+ * A `MatrixLikeFieldMap` is the base class for maps of matrices, arrays and
+ * fourth-order tensors mapped onto matrices.
*
* It should never be necessary to call directly any of the
- * constructors if this class, but rather use the template aliases
- * `muSpectre::ArrayFieldMap`, `muSpectre::MatrixFieldMap`, and
- * `muSpectre::T4MatrixFieldMap`
+ * constructors if this class, but rather use the template aliases:
+ * - `ArrayFieldMap`: iterate in the form of `Eigen::Array<...>`.
+ * - `MatrixFieldMap`: iterate in the form of
+ * `Eigen::Matrix<...>`.
+ * - `T4MatrixFieldMap`: iterate in the form of `T4MatMap`.
*/
template <class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
class MatrixLikeFieldMap: public FieldMap<FieldCollection,
typename EigenArray::Scalar,
EigenArray::SizeAtCompileTime,
ConstField>
{
public:
//! base class
using Parent = FieldMap<FieldCollection,
typename EigenArray::Scalar,
EigenArray::SizeAtCompileTime, ConstField>;
using T_t = EigenPlain; //!< plain Eigen type to map
//! cell coordinates type
using Ccoord = Ccoord_t<FieldCollection::spatial_dim()>;
using value_type = EigenArray; //!< stl conformance
using const_reference = EigenConstArray; //!< stl conformance
//! stl conformance
using reference = std::conditional_t<ConstField,
const_reference,
value_type>; // since it's a resource handle
using size_type = typename Parent::size_type; //!< stl conformance
using pointer = std::unique_ptr<EigenArray>; //!< stl conformance
using TypedField = typename Parent::TypedField; //!< stl conformance
using Field = typename TypedField::Base; //!< stl conformance
//! stl conformance
using const_iterator= typename Parent::template iterator<MatrixLikeFieldMap, true>;
//! stl conformance
using iterator = std::conditional_t<
ConstField,
const_iterator,
typename Parent::template iterator<MatrixLikeFieldMap>>;
using reverse_iterator = std::reverse_iterator<iterator>; //!< stl conformance
//! stl conformance
using const_reverse_iterator = std::reverse_iterator<const_iterator>;
//! stl conformance
friend iterator;
//! Default constructor
MatrixLikeFieldMap() = delete;
/**
* Constructor using a (non-typed) field. Compatibility is enforced at
* runtime. This should not be a performance concern, as this constructor
* will not be called in anny inner loops (if used as intended).
*/
template <bool isntConst=!ConstField>
MatrixLikeFieldMap(std::enable_if_t<isntConst, Field &> field);
/**
* Constructor using a (non-typed) field. Compatibility is enforced at
* runtime. This should not be a performance concern, as this constructor
* will not be called in anny inner loops (if used as intended).
*/
template <bool isConst=ConstField>
MatrixLikeFieldMap(std::enable_if_t<isConst, const Field &> field);
/**
* Constructor using a typed field. Compatibility is enforced
* statically. It is not always possible to call this constructor, as the
* type of the field might not be known at compile time.
*/
template<class FC, typename T2, Dim_t NbC>
MatrixLikeFieldMap(TypedSizedFieldBase<FC, T2, NbC> & field);
//! Copy constructor
MatrixLikeFieldMap(const MatrixLikeFieldMap &other) = default;
//! Move constructor
MatrixLikeFieldMap(MatrixLikeFieldMap &&other) = default;
//! Destructor
virtual ~MatrixLikeFieldMap() = default;
//! Copy assignment operator
MatrixLikeFieldMap& operator=(const MatrixLikeFieldMap &other) = delete;
//! Move assignment operator
MatrixLikeFieldMap& operator=(MatrixLikeFieldMap &&other) = delete;
//! Assign a matrixlike value to every entry
template <class Derived>
inline MatrixLikeFieldMap & operator=(const Eigen::EigenBase<Derived> & val);
//! give human-readable field map type
inline std::string info_string() const override final;
//! member access
inline reference operator[](size_type index);
//! member access
inline reference operator[](const Ccoord& ccoord);
//! member access
inline const_reference operator[](size_type index) const;
//! member access
inline const_reference operator[](const Ccoord& ccoord) const;
//! return an iterator to head of field for ranges
inline iterator begin(){return iterator(*this);}
//! return an iterator to head of field for ranges
inline const_iterator cbegin() const {return const_iterator(*this);}
//! return an iterator to head of field for ranges
inline const_iterator begin() const {return this->cbegin();}
//! return an iterator to tail of field for ranges
inline iterator end(){return iterator(*this, false);};
//! return an iterator to tail of field for ranges
inline const_iterator cend() const {return const_iterator(*this, false);}
//! return an iterator to tail of field for ranges
inline const_iterator end() const {return this->cend();}
//! evaluate the average of the field
inline T_t mean() const;
protected:
//! for sad, legacy iterator use
inline pointer ptr_to_val_t(size_type index);
const static std::string field_info_root; //!< for printing and debug
private:
};
/* ---------------------------------------------------------------------- */
template <class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
const std::string MatrixLikeFieldMap<FieldCollection, EigenArray,
EigenConstArray, EigenPlain, map_type,
ConstField>::
field_info_root{NameWrapper<map_type>::field_info_root()};
/* ---------------------------------------------------------------------- */
template<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
template <bool isntConst>
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
map_type, ConstField>::
MatrixLikeFieldMap(std::enable_if_t<isntConst, Field &> field)
:Parent(field) {
static_assert((isntConst != ConstField),
"let the compiler deduce isntConst, this is a SFINAE "
"parameter");
this->check_compatibility();
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
template <bool isConst>
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
map_type, ConstField>::
MatrixLikeFieldMap(std::enable_if_t<isConst, const Field &> field)
:Parent(field) {
static_assert((isConst == ConstField),
"let the compiler deduce isntConst, this is a SFINAE "
"parameter");
this->check_compatibility();
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
template<class FC, typename T2, Dim_t NbC>
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::
MatrixLikeFieldMap(TypedSizedFieldBase<FC, T2, NbC> & field)
:Parent(field) {
}
/* ---------------------------------------------------------------------- */
//! human-readable field map type
template<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
std::string
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::
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<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::reference
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::
operator[](size_type index) {
return reference(this->get_ptr_to_entry(index));
}
template<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::reference
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
map_type, ConstField>::
operator[](const Ccoord & ccoord) {
size_t index{};
index = this->collection.get_index(ccoord);
return reference(this->get_ptr_to_entry(std::move(index)));
}
/* ---------------------------------------------------------------------- */
//! member access
template<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::const_reference
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
map_type, ConstField>::
operator[](size_type index) const {
return const_reference(this->get_ptr_to_entry(index));
}
template<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::const_reference
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
map_type, ConstField>::
operator[](const Ccoord & ccoord) const{
size_t index{};
index = this->collection.get_index(ccoord);
return const_reference(this->get_ptr_to_entry(std::move(index)));
}
//----------------------------------------------------------------------------//
template <class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
template <class Derived>
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField> &
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::
operator=(const Eigen::EigenBase<Derived> & val) {
for (auto && entry: *this) {
entry = val;
}
return *this;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::T_t
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::
mean() const {
T_t mean{T_t::Zero()};
for (auto && val: *this) {
mean += val;
}
mean *= 1./Real(this->size());
return mean;
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, class EigenArray, class EigenConstArray,
class EigenPlain, Map_t map_type, bool ConstField>
typename MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray,
EigenPlain, map_type, ConstField>::pointer
MatrixLikeFieldMap<FieldCollection, EigenArray, EigenConstArray, EigenPlain,
map_type, ConstField>::
ptr_to_val_t(size_type index) {
return std::make_unique<value_type>
(this->get_ptr_to_entry(std::move(index)));
}
} // internal
/* ---------------------------------------------------------------------- */
//! short-hand for an Eigen matrix map as iterate
template <class FieldCollection, typename T, Dim_t NbRows, Dim_t NbCols,
bool ConstField=false>
using MatrixFieldMap = internal::MatrixLikeFieldMap
<FieldCollection,
Eigen::Map<Eigen::Matrix<T, NbRows, NbCols>>,
Eigen::Map<const Eigen::Matrix<T, NbRows, NbCols>>,
Eigen::Matrix<T, NbRows, NbCols>,
internal::Map_t::Matrix, ConstField>;
/* ---------------------------------------------------------------------- */
//! short-hand for an Eigen matrix map as iterate
template <class FieldCollection, typename T, Dim_t Dim,
bool MapConst=false>
using T4MatrixFieldMap = internal::MatrixLikeFieldMap
<FieldCollection,
T4MatMap<T, Dim, false>,
T4MatMap<T, Dim, true>,
T4Mat<T, Dim>,
internal::Map_t::T4Matrix, MapConst>;
/* ---------------------------------------------------------------------- */
//! short-hand for an Eigen array map as iterate
template <class FieldCollection, typename T, Dim_t NbRows, Dim_t NbCols,
bool ConstField=false>
using ArrayFieldMap = internal::MatrixLikeFieldMap
<FieldCollection,
Eigen::Map<Eigen::Array<T, NbRows, NbCols>>,
Eigen::Map<const Eigen::Array<T, NbRows, NbCols>>,
Eigen::Array<T, NbRows, NbCols>,
internal::Map_t::Array, ConstField>;
} // muSpectre
#endif /* FIELD_MAP_MATRIXLIKE_H */
diff --git a/src/common/field_map_scalar.hh b/src/common/field_map_scalar.hh
index 70e2eaf..3982a39 100644
--- a/src/common/field_map_scalar.hh
+++ b/src/common/field_map_scalar.hh
@@ -1,233 +1,234 @@
/**
* @file field_map_scalar.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 26 Sep 2017
*
* @brief maps over scalar fields
*
* Copyright © 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_SCALAR_H
#define FIELD_MAP_SCALAR_H
#include "common/field_map_base.hh"
namespace muSpectre {
/**
- * implements maps on scalar fields (i.e. material properties,
- * temperatures, etc)
+ * Implements maps on scalar fields (i.e. material properties,
+ * temperatures, etc). Maps onto a `muSpectre::internal::TypedSizedFieldBase`
+ * and lets you iterate over it in the form of the bare type of the field.
*/
template <class FieldCollection, typename T, bool ConstField=false>
class ScalarFieldMap
: public internal::FieldMap<FieldCollection, T, 1, ConstField>
{
public:
//! base class
using parent = internal::FieldMap<FieldCollection, T, 1, ConstField>;
//! cell coordinates type
using Ccoord = Ccoord_t<FieldCollection::spatial_dim()>;
using value_type = T; //!< stl conformance
using const_reference = const value_type &; //!< stl conformance
//! stl conformance
using reference = std::conditional_t<ConstField,
const_reference,
value_type &>;
using size_type = typename parent::size_type; //!< stl conformance
using pointer = T*; //!< stl conformance
using Field = typename parent::Field; //!< stl conformance
//! stl conformance
using const_iterator= typename parent::template iterator<ScalarFieldMap, true>;
//! iterator over a scalar field
using iterator = std::conditional_t
<ConstField,
const_iterator,
typename parent::template iterator<ScalarFieldMap>>;
using reverse_iterator = std::reverse_iterator<iterator>; //!< stl conformance
//! stl conformance
using const_reverse_iterator = std::reverse_iterator<const_iterator>;
//! give access to the protected fields
friend iterator;
//! Default constructor
ScalarFieldMap() = delete;
//! constructor
template <bool isntConst=!ConstField>
ScalarFieldMap(std::enable_if_t<isntConst, Field &> field);
//! constructor
template <bool isConst=ConstField>
ScalarFieldMap(std::enable_if_t<isConst, const Field &> field);
//! Copy constructor
ScalarFieldMap(const ScalarFieldMap &other) = default;
//! Move constructor
ScalarFieldMap(ScalarFieldMap &&other) = default;
//! Destructor
virtual ~ScalarFieldMap() = default;
//! Copy assignment operator
ScalarFieldMap& operator=(const ScalarFieldMap &other) = delete;
//! Move assignment operator
ScalarFieldMap& operator=(ScalarFieldMap &&other) = delete;
//! Assign a value to every entry
ScalarFieldMap& operator=(T val);
//! give human-readable field map type
inline std::string info_string() const override final;
//! member access
inline reference operator[](size_type index);
inline reference operator[](const Ccoord& ccoord);
inline const_reference operator[] (size_type index) const;
inline const_reference operator[] (const Ccoord& ccoord) const;
//! return an iterator to the first pixel of the field
inline iterator begin(){return iterator(*this);}
//! return an iterator to the first pixel of the field
inline const_iterator cbegin() const {return const_iterator(*this);}
//! return an iterator to the first pixel of the field
inline const_iterator begin() const {return this->cbegin();}
//! return an iterator to tail of field for ranges
inline iterator end(){return iterator(*this, false);}
//! return an iterator to tail of field for ranges
inline const_iterator cend() const {return const_iterator(*this, false);}
//! return an iterator to tail of field for ranges
inline const_iterator end() const {return this->cend();}
//! evaluate the average of the field
inline T mean() const;
protected:
//! for sad, legacy iterator use
inline pointer ptr_to_val_t(size_type index);
//! type identifier for printing and debugging
const static std::string field_info_root;
private:
};
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, bool ConstField>
template <bool isntConst>
ScalarFieldMap<FieldCollection, T, ConstField>::
ScalarFieldMap(std::enable_if_t<isntConst, Field &> field)
:parent(field) {
static_assert((isntConst != ConstField),
"let the compiler deduce isntConst, this is a SFINAE "
"parameter");
this->check_compatibility();
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, bool ConstField>
template <bool isConst>
ScalarFieldMap<FieldCollection, T, ConstField>::
ScalarFieldMap(std::enable_if_t<isConst, const Field &> field)
:parent(field) {
static_assert((isConst == ConstField),
"let the compiler deduce isntConst, this is a SFINAE "
"parameter");
this->check_compatibility();
}
/* ---------------------------------------------------------------------- */
//! human-readable field map type
template<class FieldCollection, typename T, bool ConstField>
std::string
ScalarFieldMap<FieldCollection, T, ConstField>::info_string() const {
std::stringstream info;
info << "Scalar(" << typeid(T).name() << ")";
return info.str();
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, bool ConstField>
const std::string ScalarFieldMap<FieldCollection, T, ConstField>::field_info_root{
"Scalar"};
/* ---------------------------------------------------------------------- */
//! member access
template <class FieldCollection, typename T, bool ConstField>
typename ScalarFieldMap<FieldCollection, T, ConstField>::reference
ScalarFieldMap<FieldCollection, T, ConstField>::operator[](size_type index) {
return this->get_ptr_to_entry(std::move(index))[0];
}
/* ---------------------------------------------------------------------- */
//! member access
template <class FieldCollection, typename T, bool ConstField>
typename ScalarFieldMap<FieldCollection, T, ConstField>::reference
ScalarFieldMap<FieldCollection, T, ConstField>::operator[](const Ccoord& ccoord) {
auto && index = this->collection.get_index(std::move(ccoord));
return this->get_ptr_to_entry(std::move(index))[0];
}
/* ---------------------------------------------------------------------- */
//! member access
template <class FieldCollection, typename T, bool ConstField>
typename ScalarFieldMap<FieldCollection, T, ConstField>::const_reference
ScalarFieldMap<FieldCollection, T, ConstField>::
operator[](size_type index) const {
return this->get_ptr_to_entry(std::move(index))[0];
}
/* ---------------------------------------------------------------------- */
//! member access
template <class FieldCollection, typename T, bool ConstField>
typename ScalarFieldMap<FieldCollection, T, ConstField>::const_reference
ScalarFieldMap<FieldCollection, T, ConstField>::
operator[](const Ccoord& ccoord) const {
auto && index = this->collection.get_index(std::move(ccoord));
return this->get_ptr_to_entry(std::move(index))[0];
}
/* ---------------------------------------------------------------------- */
//! Assign a value to every entry
template <class FieldCollection, typename T, bool ConstField>
ScalarFieldMap<FieldCollection, T, ConstField> &
ScalarFieldMap<FieldCollection, T, ConstField>::operator=(T val) {
for (auto & scalar:*this) {
scalar = val;
}
return *this;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, bool ConstField>
T
ScalarFieldMap<FieldCollection, T, ConstField>::
mean() const {
T mean{0};
for (auto && val: *this) {
mean += val;
}
mean /= Real(this->size());
return mean;
}
} // muSpectre
#endif /* FIELD_MAP_SCALAR_H */
diff --git a/src/common/field_map_tensor.hh b/src/common/field_map_tensor.hh
index cab4d16..e112ef7 100644
--- a/src/common/field_map_tensor.hh
+++ b/src/common/field_map_tensor.hh
@@ -1,287 +1,286 @@
/**
* @file field_map_tensor.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 26 Sep 2017
*
* @brief Defines an Eigen-Tensor map over strongly typed fields
*
* Copyright © 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 "common/eigen_tools.hh"
#include "common/field_map_base.hh"
#include <sstream>
#include <memory>
namespace muSpectre {
/* ---------------------------------------------------------------------- */
/**
- * maps onto a `muSpectre::internal::TypedSizedFieldBase` and lets
- * you iterate over it in the form ot
- * `Eigen::TensorMap<TensorFixedSize<...>>`
+ * Maps onto a `muSpectre::internal::TypedSizedFieldBase` and lets
+ * you iterate over it in the form of `Eigen::TensorMap<TensorFixedSize<...>>`
*/
template <class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField=false>
class TensorFieldMap: public
internal::FieldMap<FieldCollection,
T,
SizesByOrder<order, dim>::Sizes::total_size,
ConstField>
{
public:
//! base class
using Parent = internal::FieldMap<FieldCollection, T,
TensorFieldMap::nb_components, ConstField>;
//! cell coordinates type
using Ccoord = Ccoord_t<FieldCollection::spatial_dim()>;
//! static tensor size
using Sizes = typename SizesByOrder<order, dim>::Sizes;
//! plain Eigen type
using T_t = Eigen::TensorFixedSize<T, Sizes>;
using value_type = Eigen::TensorMap<T_t>; //!< stl conformance
using const_reference = Eigen::TensorMap<const T_t>; //!< stl conformance
//! stl conformance
using reference = std::conditional_t<ConstField,
const_reference,
value_type>; // since it's a resource handle
using size_type = typename Parent::size_type; //!< stl conformance
using pointer = std::unique_ptr<value_type>; //!< stl conformance
using TypedField = typename Parent::TypedField; //!< field to map
//! polymorphic base field type (for debug and python)
using Field = typename TypedField::Base;
//! stl conformance
using const_iterator = typename Parent::template iterator<TensorFieldMap, true>;
//! stl conformance
using iterator = std::conditional_t<
ConstField,
const_iterator,
typename Parent::template iterator<TensorFieldMap>>;
//! stl conformance
using reverse_iterator = std::reverse_iterator<iterator>;
//! stl conformance
using const_reverse_iterator = std::reverse_iterator<const_iterator>;
//! stl conformance
friend iterator;
//! Default constructor
TensorFieldMap() = delete;
//! constructor
template <bool isntConst=!ConstField>
TensorFieldMap(std::enable_if_t<isntConst, Field &> field);
//! constructor
template <bool isConst=ConstField>
TensorFieldMap(std::enable_if_t<isConst, const Field &> field);
//! Copy constructor
TensorFieldMap(const TensorFieldMap &other) = default;
//! Move constructor
TensorFieldMap(TensorFieldMap &&other) = default;
//! Destructor
virtual ~TensorFieldMap() = default;
//! Copy assignment operator
TensorFieldMap& operator=(const TensorFieldMap &other) = delete;
//! Assign a matrixlike value to every entry
inline TensorFieldMap & operator=(const T_t & val);
//! Move assignment operator
TensorFieldMap& operator=(TensorFieldMap &&other) = delete;
//! give human-readable field map type
inline std::string info_string() const override final;
//! member access
inline reference operator[](size_type index);
//! member access
inline reference operator[](const Ccoord & ccoord);
//! member access
inline const_reference operator[](size_type index) const;
//! member access
inline const_reference operator[](const Ccoord& ccoord) const;
//! return an iterator to first entry of field
inline iterator begin(){return iterator(*this);}
//! return an iterator to first entry of field
inline const_iterator cbegin() const {return const_iterator(*this);}
//! return an iterator to first entry of field
inline const_iterator begin() const {return this->cbegin();}
//! return an iterator past the last entry of field
inline iterator end(){return iterator(*this, false);}
//! return an iterator past the last entry of field
inline const_iterator cend() const {return const_iterator(*this, false);}
//! return an iterator past the last entry of field
inline const_iterator end() const {return this->cend();}
//! evaluate the average of the field
inline T_t mean() const;
protected:
//! for sad, legacy iterator use
inline pointer ptr_to_val_t(size_type index);
private:
};
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
template <bool isntConst>
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
TensorFieldMap(std::enable_if_t<isntConst, Field &> field)
:Parent(field) {
static_assert((isntConst != ConstField),
"let the compiler deduce isntConst, this is a SFINAE "
"parameter");
this->check_compatibility();
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
template <bool isConst>
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
TensorFieldMap(std::enable_if_t<isConst, const Field &> field)
:Parent(field) {
static_assert((isConst == ConstField),
"let the compiler deduce isntConst, this is a SFINAE "
"parameter");
this->check_compatibility();
}
/* ---------------------------------------------------------------------- */
//! human-readable field map type
template<class FieldCollection, typename T, Dim_t order, Dim_t dim, bool ConstField>
std::string
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::info_string() const {
std::stringstream info;
info << "Tensor(" << typeid(T).name() << ", " << order
<< "_o, " << dim << "_d)";
return info.str();
}
/* ---------------------------------------------------------------------- */
//! member access
template <class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::reference
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
operator[](size_type index) {
auto && lambda = [this, &index](auto&&...tens_sizes) {
return reference(this->get_ptr_to_entry(index), tens_sizes...);
};
return call_sizes<order, dim>(lambda);
}
template<class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::reference
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
operator[](const Ccoord & ccoord) {
auto && index = this->collection.get_index(ccoord);
auto && lambda = [this, &index](auto&&...sizes) {
return reference(this->get_ptr_to_entry(index), sizes...);
};
return call_sizes<order, dim>(lambda);
}
template <class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
const_reference
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
operator[](size_type index) const {
// Warning: due to a inconsistency in Eigen's API, tensor maps
// cannot be constructed from a const ptr, hence this nasty const
// cast :(
auto && lambda = [this, &index](auto&&...tens_sizes) {
return const_reference(const_cast<T*>(this->get_ptr_to_entry(index)),
tens_sizes...);
};
return call_sizes<order, dim>(lambda);
}
template<class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
const_reference
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
operator[](const Ccoord & ccoord) const {
auto && index = this->collection.get_index(ccoord);
auto && lambda = [this, &index](auto&&...sizes) {
return const_reference(const_cast<T*>(this->get_ptr_to_entry(index)),
sizes...);
};
return call_sizes<order, dim>(lambda);
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
TensorFieldMap<FieldCollection, T, order, dim, ConstField> &
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
operator=(const T_t & val) {
for (auto && tens: *this) {
tens = val;
}
return *this;
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::T_t
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
mean() const {
T_t mean{T_t::Zero()};
for (auto && val: *this) {
mean += val;
}
mean *= 1./Real(this->size());
return mean;
}
/* ---------------------------------------------------------------------- */
//! for sad, legacy iterator use. Don't use unless you have to.
template<class FieldCollection, typename T, Dim_t order, Dim_t dim,
bool ConstField>
typename TensorFieldMap<FieldCollection, T, order, dim, ConstField>::pointer
TensorFieldMap<FieldCollection, T, order, dim, ConstField>::
ptr_to_val_t(size_type index) {
auto && lambda = [this, &index](auto&&... tens_sizes) {
return std::make_unique<value_type>(this->get_ptr_to_entry(index),
tens_sizes...);
};
return call_sizes<order, dim>(lambda);
}
} // muSpectre
#endif /* FIELD_MAP_TENSOR_H */