diff --git a/src/common/field.hh b/src/common/field.hh
index a6802f3..dda09be 100644
--- a/src/common/field.hh
+++ b/src/common/field.hh
@@ -1,535 +1,536 @@
/**
* file field.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 07 Sep 2017
*
* @brief header-only implementation of a field for field collections
*
* @section LICENCE
*
* Copyright (C) 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef FIELD_H
#define FIELD_H
#include <string>
#include <sstream>
#include <utility>
#include <typeinfo>
#include <vector>
#include <algorithm>
#include <Eigen/Dense>
#include <cmath>
#include <memory>
#include <type_traits>
namespace muSpectre {
namespace internal {
/* ---------------------------------------------------------------------- */
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;
//! Copy constructor
FieldBase(const FieldBase &other) = delete;
//! Move constructor
FieldBase(FieldBase &&other) noexcept = delete;
//! Destructor
virtual ~FieldBase() noexcept = default;
//! Copy assignment operator
FieldBase& operator=(const FieldBase &other) = delete;
//! Move assignment operator
FieldBase& operator=(FieldBase &&other) noexcept = delete;
/* ---------------------------------------------------------------------- */
//!Identifying accessors
//! return field name
inline const std::string & get_name() const;
//! return field type
//inline const Field_t & get_type() const;
//! return my collection (for iterating)
inline const FieldCollection & get_collection() const;
//! return my collection (for iterating)
inline const size_t & get_nb_components() const;
//! return type_id of stored type
virtual const std::type_info & get_stored_typeid() const = 0;
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;
+ friend typename FieldCollection::Parent;
protected:
/* ---------------------------------------------------------------------- */
//! allocate memory etc
virtual void resize(size_t size) = 0;
const std::string name;
const size_t nb_components;
const FieldCollection & collection;
private:
};
/* ---------------------------------------------------------------------- */
//! declaraton for friending
template <class FieldCollection, typename T, Dim_t NbComponents, bool isConst>
class FieldMap;
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents,
bool ArrayStore=false>
class TypedFieldBase: public FieldBase<FieldCollection>
{
friend class FieldMap<FieldCollection, T, NbComponents, true>;
friend class FieldMap<FieldCollection, T, NbComponents, false>;
public:
constexpr static auto nb_components{NbComponents};
using Parent = FieldBase<FieldCollection>;
using Parent::collection_t;
using Scalar = T;
using Base = Parent;
//using storage_type = Eigen::Array<T, Eigen::Dynamic, NbComponents>;
using StoredType = Eigen::Array<T, NbComponents, 1>;
using StorageType = std::conditional_t
<ArrayStore,
std::vector<StoredType,
Eigen::aligned_allocator<StoredType>>,
std::vector<T,Eigen::aligned_allocator<T>>>;
TypedFieldBase(std::string unique_name,
FieldCollection& collection);
virtual ~TypedFieldBase() = default;
//! return type_id of stored type
virtual const std::type_info & get_stored_typeid() const;
//! initialise field to zero (do more complicated initialisations through
//! fully typed maps)
inline void set_zero() override final;
template <bool isArrayStore = ArrayStore>
inline void push_back(const
std::enable_if_t<isArrayStore, StoredType> &
value);
template <bool componentStore = !ArrayStore>
inline std::enable_if_t<componentStore> push_back(const StoredType & value);
size_t size() const override final;
static TypedFieldBase & check_ref(Parent & other);
static const TypedFieldBase & check_ref(const Base & parent);
inline T* data() {return this->get_ptr_to_entry(0);}
inline const T* data() const {return this->get_ptr_to_entry(0);}
protected:
template <bool isArray=ArrayStore>
inline std::enable_if_t<isArray, T*> get_ptr_to_entry(const size_t&& index);
template <bool isArray=ArrayStore>
inline std::enable_if_t<isArray, const T*>
get_ptr_to_entry(const size_t&& index) const;
template <bool noArray = !ArrayStore>
inline T* get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index);
template <bool noArray = !ArrayStore>
inline const T*
get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index) const;
inline virtual void resize(size_t size) override final;
StorageType array{};
};
} // internal
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t order, Dim_t dim>
class TensorField: public internal::TypedFieldBase<FieldCollection,
T,
ipow(dim,order)>
{
public:
using Parent = internal::TypedFieldBase<FieldCollection,
T,
ipow(dim,order)>;
using Base = typename Parent::Base;
using Field_p = typename FieldCollection::Field_p;
using component_type = T;
//! Copy constructor
TensorField(const TensorField &other) = delete;
//! Move constructor
TensorField(TensorField &&other) noexcept = delete;
//! Destructor
virtual ~TensorField() noexcept = default;
//! Copy assignment operator
TensorField& operator=(const TensorField &other) = delete;
//! Move assignment operator
TensorField& operator=(TensorField &&other) noexcept = delete;
//! accessors
inline Dim_t get_order() const;
inline Dim_t get_dim() const;
//! factory function
template<class FieldType, class CollectionType, typename... Args>
friend FieldType& make_field(std::string unique_name,
CollectionType & collection,
Args&&... args);
static TensorField & check_ref(Base & other) {
return static_cast<TensorField &>(Parent::check_ref(other));}
static const TensorField & check_ref(const Base & other) {
return static_cast<const TensorField &>(Parent::check_ref(other));}
protected:
//! constructor protected!
TensorField(std::string unique_name, FieldCollection & collection);
private:
};
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbRow, Dim_t NbCol=NbRow>
class MatrixField: public internal::TypedFieldBase<FieldCollection,
T,
NbRow*NbCol>
{
public:
using Parent = internal::TypedFieldBase<FieldCollection,
T,
NbRow*NbCol>;
using Base = typename Parent::Base;
using Field_p = std::unique_ptr<internal::FieldBase<FieldCollection>>;
using component_type = T;
//! Copy constructor
MatrixField(const MatrixField &other) = delete;
//! Move constructor
MatrixField(MatrixField &&other) noexcept = delete;
//! Destructor
virtual ~MatrixField() noexcept = default;
//! Copy assignment operator
MatrixField& operator=(const MatrixField &other) = delete;
//! Move assignment operator
MatrixField& operator=(MatrixField &&other) noexcept = delete;
//! accessors
inline Dim_t get_nb_row() const;
inline Dim_t get_nb_col() const;
//! factory function
template<class FieldType, class CollectionType, typename... Args>
friend FieldType& make_field(std::string unique_name,
CollectionType & collection,
Args&&... args);
static MatrixField & check_ref(Base & other) {
return static_cast<MatrixField &>(Parent::check_ref(other));}
static const MatrixField & check_ref(const Base & other) {
return static_cast<const MatrixField &>(Parent::check_ref(other));}
protected:
//! constructor protected!
MatrixField(std::string unique_name, FieldCollection & collection);
private:
};
/* ---------------------------------------------------------------------- */
//! convenience alias
template <class FieldCollection, typename T>
using ScalarField = TensorField<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, Dim_t NbComponents, bool ArrayStore>
TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
TypedFieldBase(std::string unique_name, FieldCollection & collection)
:FieldBase<FieldCollection>(unique_name, NbComponents, collection){
static_assert
((std::is_arithmetic<T>::value ||
std::is_same<T, Complex>::value),
"Use TypedFieldBase for integer, real or complex scalars for T");
static_assert(NbComponents > 0,
"Only fields with more than 0 components");
}
/* ---------------------------------------------------------------------- */
//! return type_id of stored type
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
const std::type_info & TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
get_stored_typeid() const {
return typeid(T);
}
/* ---------------------------------------------------------------------- */
template<class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
void
TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
set_zero() {
std::fill(this->array.begin(), this->array.end(), T{});
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
size_t TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
size() const {
if (ArrayStore) {
return this->array.size();
} else {
return this->array.size()/NbComponents;
}
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore> &
TypedFieldBase<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<TypedFieldBase&>(other);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
const TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore> &
TypedFieldBase<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 TypedFieldBase&>(other);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool isArray>
std::enable_if_t<isArray, T*>
TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
get_ptr_to_entry(const size_t&& index) {
static_assert (isArray == ArrayStore, "SFINAE");
return &this->array[std::move(index)](0, 0);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool noArray>
T* TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
get_ptr_to_entry(std::enable_if_t<noArray, const size_t&&> index) {
static_assert (noArray != ArrayStore, "SFINAE");
return &this->array[NbComponents*std::move(index)];
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool isArray>
std::enable_if_t<isArray, const T*>
TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
get_ptr_to_entry(const size_t&& index) const {
static_assert (isArray == ArrayStore, "SFINAE");
return &this->array[std::move(index)](0, 0);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool noArray>
const T* TypedFieldBase<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->array[NbComponents*std::move(index)];
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
void TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
resize(size_t size) {
if (ArrayStore) {
this->array.resize(size);
} else {
this->array.resize(size*NbComponents);
}
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool isArrayStore>
void
TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
push_back(const std::enable_if_t<isArrayStore,StoredType> & value) {
static_assert(isArrayStore == ArrayStore, "SFINAE");
this->array.push_back(value);
}
/* ---------------------------------------------------------------------- */
template <class FieldCollection, typename T, Dim_t NbComponents, bool ArrayStore>
template <bool componentStore>
std::enable_if_t<componentStore>
TypedFieldBase<FieldCollection, T, NbComponents, ArrayStore>::
push_back(const StoredType & value) {
static_assert(componentStore != ArrayStore, "SFINAE");
for (Dim_t i = 0; i < NbComponents; ++i) {
this->array.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
#endif /* FIELD_H */
diff --git a/src/common/field_collection_base.hh b/src/common/field_collection_base.hh
index 1c16da9..669c6cc 100644
--- a/src/common/field_collection_base.hh
+++ b/src/common/field_collection_base.hh
@@ -1,185 +1,188 @@
/**
* file field_collection_base.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 05 Nov 2017
*
* @brief Base class for field collections
*
* @section LICENCE
*
* Copyright (C) 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef FIELD_COLLECTION_BASE_H
#define FIELD_COLLECTION_BASE_H
#include <map>
#include <vector>
#include "common/common.hh"
#include "common/field.hh"
namespace muSpectre {
/* ---------------------------------------------------------------------- */
class FieldCollectionError: public std::runtime_error {
public:
explicit FieldCollectionError(const std::string& what)
:std::runtime_error(what){}
explicit FieldCollectionError(const char * what)
:std::runtime_error(what){}
};
class FieldError: public FieldCollectionError {
using Parent = FieldCollectionError;
public:
explicit FieldError(const std::string& what)
:Parent(what){}
explicit FieldError(const char * what)
:Parent(what){}
};
class FieldInterpretationError: public FieldError
{
public:
explicit FieldInterpretationError(const std::string & what)
:FieldError(what){}
explicit FieldInterpretationError(const char * what)
:FieldError(what){}
};
/* ---------------------------------------------------------------------- */
//! DimS spatial dimension (dimension of problem)
//! DimM material_dimension (dimension of constitutive law)
template <Dim_t DimS, Dim_t DimM, class FieldCollectionDerived>
class FieldCollectionBase
{
public:
using Field = internal::FieldBase<FieldCollectionDerived>;
using Field_p = std::unique_ptr<Field>;
using Ccoord = Ccoord_t<DimS>;
//! Default constructor
FieldCollectionBase();
//! Copy constructor
FieldCollectionBase(const FieldCollectionBase &other) = delete;
//! Move constructor
FieldCollectionBase(FieldCollectionBase &&other) noexcept = delete;
//! Destructor
virtual ~FieldCollectionBase() noexcept = default;
//! Copy assignment operator
FieldCollectionBase& operator=(const FieldCollectionBase &other) = delete;
//! Move assignment operator
FieldCollectionBase& operator=(FieldCollectionBase &&other) noexcept = 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_;}
protected:
std::map<const std::string, Field_p> fields{};
bool is_initialised{false};
const Uint id;
static Uint counter;
size_t size_{0};
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));
}
} // muSpectre
#endif /* FIELD_COLLECTION_BASE_H */
diff --git a/src/common/field_collection_global.hh b/src/common/field_collection_global.hh
index d04e0c5..75fd5db 100644
--- a/src/common/field_collection_global.hh
+++ b/src/common/field_collection_global.hh
@@ -1,186 +1,189 @@
/**
* file field_collection_global.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 05 Nov 2017
*
* @brief FieldCollection base-class for global fields
*
* @section LICENCE
*
* Copyright (C) 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef FIELD_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)
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>>;
using Ccoord = typename Parent::Ccoord;
using Field_p = typename Parent::Field_p;
using iterator = typename CcoordOps::Pixels<DimS>::iterator;
//! Default constructor
GlobalFieldCollection();
//! Copy constructor
GlobalFieldCollection(const GlobalFieldCollection &other) = delete;
//! Move constructor
GlobalFieldCollection
(GlobalFieldCollection &&other) noexcept = delete;
//! Destructor
virtual ~GlobalFieldCollection() noexcept = default;
//! Copy assignment operator
GlobalFieldCollection&
operator=(const GlobalFieldCollection &other) = delete;
//! Move assignment operator
GlobalFieldCollection&
operator=(GlobalFieldCollection &&other) noexcept = delete;
/** 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'
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();
inline iterator end();
static constexpr inline Dim_t spatial_dim() {return DimS;}
static constexpr inline Dim_t material_dim() {return DimM;}
protected:
//! number of discretisation cells in each of the DimS spatial directions
Ccoord sizes{0};
CcoordOps::Pixels<DimS> pixels;
private:
};
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
GlobalFieldCollection<DimS, DimM>::GlobalFieldCollection()
:Parent(), pixels{{0}}{}
/* ---------------------------------------------------------------------- */
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_ = std::accumulate(sizes.begin(), sizes.end(), 1,
std::multiplies<Dim_t>());
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 2348238..a8ff2ad 100644
--- a/src/common/field_collection_local.hh
+++ b/src/common/field_collection_local.hh
@@ -1,171 +1,174 @@
/**
* file field_collection_local.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 05 Nov 2017
*
* @brief FieldCollection base-class for local fields
*
* @section LICENCE
*
* Copyright (C) 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef FIELD_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)
template<Dim_t DimS, Dim_t DimM>
class LocalFieldCollection:
public FieldCollectionBase<DimS, DimM, LocalFieldCollection<DimS, DimM>>
{
public:
using Parent = FieldCollectionBase<DimS, DimM,
LocalFieldCollection<DimS, DimM>>;
using Ccoord = typename Parent::Ccoord;
using Field_p = typename Parent::Field_p;
using ccoords_container = std::vector<Ccoord>;
using iterator = typename ccoords_container::iterator;
//! Default constructor
LocalFieldCollection();
//! Copy constructor
LocalFieldCollection(const LocalFieldCollection &other) = delete;
//! Move constructor
LocalFieldCollection(LocalFieldCollection &&other) noexcept = delete;
//! Destructor
virtual ~LocalFieldCollection() noexcept = default;
//! Copy assignment operator
LocalFieldCollection& operator=(const LocalFieldCollection &other) = delete;
//! Move assignment operator
LocalFieldCollection& operator=(LocalFieldCollection &&other) noexcept = 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;
inline iterator begin() {return this->ccoords.begin();}
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);
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void LocalFieldCollection<DimS, DimM>::
initialise() {
+ if (this->is_initialised) {
+ throw std::runtime_error("double initialisation");
+ }
this->size_ = this->ccoords.size();
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/fft/projection_base.hh b/src/fft/projection_base.hh
index 8fc8735..479372b 100644
--- a/src/fft/projection_base.hh
+++ b/src/fft/projection_base.hh
@@ -1,93 +1,100 @@
/**
* file projection_base.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 03 Dec 2017
*
* @brief Base class for Projection operators
*
* @section LICENCE
*
* Copyright (C) 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef PROJECTION_BASE_H
#define PROJECTION_BASE_H
#include "common/common.hh"
#include "common/field_collection.hh"
#include "common/field.hh"
#include "fft/fft_engine_base.hh"
namespace muSpectre {
template<class Projection>
struct Projection_traits {
};
template <Dim_t DimS, Dim_t DimM>
class ProjectionBase
{
public:
using FFT_Engine = FFT_Engine_base<DimS, DimM>;
using Ccoord = typename FFT_Engine::Ccoord;
+ using Rcoord = typename FFT_Engine::Rcoord;
using GFieldCollection_t = typename FFT_Engine::GFieldCollection_t;
using LFieldCollection_t = typename FFT_Engine::LFieldCollection_t;
using Field_t = typename FFT_Engine::Field_t;
using iterator = typename FFT_Engine::iterator;
//! Default constructor
ProjectionBase() = delete;
//! Constructor with system sizes
ProjectionBase(FFT_Engine & engine);
//! Copy constructor
ProjectionBase(const ProjectionBase &other) = delete;
//! Move constructor
ProjectionBase(ProjectionBase &&other) noexcept = default;
//! Destructor
virtual ~ProjectionBase() noexcept = default;
//! Copy assignment operator
ProjectionBase& operator=(const ProjectionBase &other) = delete;
//! Move assignment operator
ProjectionBase& operator=(ProjectionBase &&other) noexcept = default;
//! initialises the fft engine (plan the transform)
void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate);
//! apply the projection operator to a field
void apply_projection(Field_t & field) const;
+ //!
+ const Ccoord & get_resolutions() const {
+ return this->fft_engine.get_resolutions();}
+ const Rcoord & get_lengths() const {
+ return this->fft_engine.get_lengths();}
+
protected:
FFT_Engine & fft_engine;
LFieldCollection_t & projection_container{};
private:
};
} // muSpectre
#endif /* PROJECTION_BASE_H */
diff --git a/src/materials/material_base.hh b/src/materials/material_base.hh
index 1da34a2..7f55cc3 100644
--- a/src/materials/material_base.hh
+++ b/src/materials/material_base.hh
@@ -1,139 +1,143 @@
/**
* file material_base.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 25 Oct 2017
*
* @brief Base class for materials (constitutive models)
*
* @section LICENCE
*
* Copyright (C) 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include <string>
#include "common/common.hh"
#include "common/field_map_matrixlike.hh"
#include "common/field_collection.hh"
#ifndef MATERIAL_BASE_H
#define MATERIAL_BASE_H
namespace muSpectre {
//! DimS spatial dimension (dimension of problem
//! DimM material_dimension (dimension of constitutive law)
template<Dim_t DimS, Dim_t DimM>
class MaterialBase
{
public:
//! typedefs for data handled by this interface
//! global field collection for system-wide fields, like stress, strain, etc
using GFieldCollection_t = FieldCollection<DimS, DimM, true>;
//! field collection for internal variables, such as eigen-strains,
//! plastic strains, damage variables, etc, but also for managing which
//! pixels the material is responsible for
using MFieldCollection_t = FieldCollection<DimS, DimM, false>;
+ using iterator = typename MFieldCollection_t::iterator;
using Field_t = internal::FieldBase<GFieldCollection_t>;
using StressField_t = TensorField<GFieldCollection_t, Real, secondOrder, DimM>;
using StrainField_t = StressField_t;
using TangentField_t = TensorField<GFieldCollection_t, Real, fourthOrder, DimM>;
using Ccoord = Ccoord_t<DimS>;
//! Default constructor
MaterialBase() = delete;
//! Construct by name
MaterialBase(std::string name);
//! Copy constructor
MaterialBase(const MaterialBase &other) = delete;
//! Move constructor
MaterialBase(MaterialBase &&other) noexcept = delete;
//! Destructor
virtual ~MaterialBase() noexcept = default;
//! Copy assignment operator
MaterialBase& operator=(const MaterialBase &other) = delete;
//! Move assignment operator
MaterialBase& operator=(MaterialBase &&other) noexcept = delete;
/**
* take responsibility for a pixel identified by its cell coordinates
* WARNING: this won't work for materials with additional info per pixel
* (as, e.g. for eigenstrain), we need to pass more parameters. Materials
* of this tye need to overload add_pixel
*/
void add_pixel(const Ccoord & ccooord);
//! allocate memory, etc
virtual void initialise(bool stiffness = false) = 0;
//! return the materil's name
const std::string & get_name() const;
//! for static inheritance stuff
constexpr static Dim_t sdim() {return DimS;}
constexpr static Dim_t mdim() {return DimM;}
//! computes stress
virtual void compute_stresses(const StrainField_t & F,
StressField_t & P,
Formulation form) = 0;
/**
* Convenience function to compute stresses, mostly for debugging and
* testing. Has runtime-cost associated with compatibility-checking and
* conversion of the Field_t arguments that can be avoided by using the
* version with strongly typed field references
*/
void compute_stresses(const Field_t & F,
Field_t & P,
Formulation form);
//! computes stress and tangent moduli
virtual void compute_stresses_tangent(const StrainField_t & F,
StressField_t & P,
TangentField_t & K,
Formulation form) = 0;
/**
* Convenience function to compute stresses and tangent moduli, mostly for
* debugging and testing. Has runtime-cost associated with
* compatibility-checking and conversion of the Field_t arguments that can
* be avoided by using the version with strongly typed field references
*/
void compute_stresses_tangent(const Field_t & F,
Field_t & P,
Field_t & K,
Formulation form);
+ inline iterator begin() {return this->internal_fields.begin();}
+ inline iterator end() {return this->internal_fields.end();}
+
protected:
//! members
const std::string name;
MFieldCollection_t internal_fields{};
private:
};
} // muSpectre
#endif /* MATERIAL_BASE_H */
diff --git a/src/system/CMakeLists.txt b/src/system/CMakeLists.txt
index 3ff0f3c..a3c3189 100644
--- a/src/system/CMakeLists.txt
+++ b/src/system/CMakeLists.txt
@@ -1,9 +1,7 @@
set (systems_SRC
system_base.cc
- #fftw_engine_c2c.cc
- #fftw_engine_r2c.cc
)
find_package(FFTW REQUIRED)
add_library(systems ${systems_SRC})
-target_link_libraries(systems Eigen3::Eigen common ${FFTW_LIBRARIES})
+target_link_libraries(systems Eigen3::Eigen common)
diff --git a/src/system/system_base.cc b/src/system/system_base.cc
index 8286acf..4051497 100644
--- a/src/system/system_base.cc
+++ b/src/system/system_base.cc
@@ -1,37 +1,127 @@
/**
* file system_base.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 01 Nov 2017
*
* @brief Implementation for system base class
*
* @section LICENCE
*
* Copyright (C) 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
+#include <sstream>
+#include <algorithm>
+
#include "system/system_base.hh"
+#include "common/ccoord_operations.hh"
namespace muSpectre {
-
+ /* ---------------------------------------------------------------------- */
+ template <Dim_t DimS, Dim_t DimM>
+ SystemBase<DimS, DimM>::SystemBase(Projection_ptr projection)
+ :resolutions{projection->get_resolutions()},
+ lengths{projection->get_lengths()},
+ F{make_field<StrainField_t>("Gradient", this->fields)},
+ P{make_field<StressField_t>("Piola-Kirchhoff-1", this->fields)},
+ K_ptr{nullptr}, projection{std::move(projection)}
+ {}
+
+ /* ---------------------------------------------------------------------- */
+ template <Dim_t DimS, Dim_t DimM>
+ void SystemBase<DimS, DimM>::add_material(Material_ptr mat) {
+ this->materials.push_back(std::move(mat));
+ }
+
+ /* ---------------------------------------------------------------------- */
+ template <Dim_t DimS, Dim_t DimM>
+ typename SystemBase<DimS, DimM>::FullResponse_t
+ SystemBase<DimS, DimM>::evaluate_stress_tangent(StrainField_t & grad) {
+ if (this->is_initialised == false) {
+ this->initialise();
+ }
+ //! High level compatibility checks
+ if (grad.size() != this->F.size()) {
+ throw std::runtime_error("Size mismatch");
+ }
+ if (this->K_ptr == nullptr) {
+ K_ptr = &make_field<TangentField_t>("Tangent Stiffness", this->fields);
+ }
+
+ for (auto & mat: this->materials) {
+ mat->compute_stresses_tangent(grad, this->P, *this->K_ptr, form);
+ }
+ return std::tie(this->P, *this->K_ptr);
+ }
+
+ /* ---------------------------------------------------------------------- */
+ template <Dim_t DimS, Dim_t DimM>
+ void SystemBase<DimS, DimM>::initialise() {
+ this->check_material_coverage();
+ this->fields.initialise(this->resolutions);
+ this->is_initialised = true;
+ }
+
+ /* ---------------------------------------------------------------------- */
+ template <Dim_t DimS, Dim_t DimM>
+ void SystemBase<DimS, DimM>::check_material_coverage() {
+ auto nb_pixels = CcoordOps::get_size(this->resolutions);
+ std::vector<MaterialBase<DimS, DimM>*> assignments(nb_pixels, nullptr);
+ for (auto & mat: this->materials) {
+ for (auto & pixel: *mat) {
+ auto index = CcoordOps::get_index(this->resolutions, pixel);
+ auto& assignment{assignments.at(index)};
+ if (assignment != nullptr) {
+ std::stringstream err{};
+ err << "Pixel " << pixel << "is already assigned to material '"
+ << assignment->get_name()
+ << "' and cannot be reassigned to material '" << mat->get_name();
+ throw std::runtime_error(err.str());
+ } else {
+ assignments[index] = mat.get();
+ }
+ }
+ }
+
+ // find and identify unassigned pixels
+ std::vector<Ccoord> unassigned_pixels;
+ for (size_t i = 0; i < assignments.size(); ++i) {
+ if (assignments[i] == nullptr) {
+ unassigned_pixels.push_back(CcoordOps::get_ccoord(this->resolutions, i));
+ }
+ }
+
+ if (unassigned_pixels.size() != 0) {
+ std::stringstream err {};
+ err << "The following pixels have were not assigned a material: ";
+ for (auto & pixel: unassigned_pixels) {
+ err << pixel << ", ";
+ }
+ err << "and that cannot be handled";
+ throw std::runtime_error(err.str());
+ }
+ }
+
+ template class SystemBase<twoD, twoD>;
+ template class SystemBase<threeD, threeD>;
} // muSpectre
diff --git a/src/system/system_base.hh b/src/system/system_base.hh
index d6afc2c..d35adc9 100644
--- a/src/system/system_base.hh
+++ b/src/system/system_base.hh
@@ -1,46 +1,113 @@
/**
* file system_base.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 01 Nov 2017
*
- * @brief Base class representing a a unit cell system
+ * @brief Base class representing a unit cell system with single
+ * projection operator
*
* @section LICENCE
*
* Copyright (C) 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
-#include <common/common.hh>
+#include <vector>
+#include <memory>
+#include <tuple>
+
+#include "common/common.hh"
+#include "common/field_collection.hh"
+#include "materials/material_base.hh"
+#include "fft/projection_base.hh"
+
#ifndef SYSTEM_BASE_H
#define SYSTEM_BASE_H
namespace muSpectre {
//! DimS spatial dimension (dimension of problem
//! DimM material_dimension (dimension of constitutive law)
template<Dim_t DimS, Dim_t DimM>
class SystemBase
{
+ public:
+ constexpr static Formulation form{Formulation::finite_strain};
+ using Ccoord = Ccoord_t<DimS>;
+ using Rcoord = Rcoord_t<DimS>;
+ using FieldCollection_t = FieldCollection<DimS, DimM>;
+ using Material_ptr = std::unique_ptr<MaterialBase<DimS, DimM>>;
+ using Projection_ptr = std::unique_ptr<ProjectionBase<DimS, DimM>>;
+ using StrainField_t = TensorField<FieldCollection_t, Real, secondOrder, DimM>;
+ using StressField_t = TensorField<FieldCollection_t, Real, secondOrder, DimM>;
+ using TangentField_t = TensorField<FieldCollection_t, Real, secondOrder, DimM>;
+ using FullResponse_t = std::tuple<const StressField_t&, const TangentField_t&>;
+
+ //! Default constructor
+ SystemBase() = delete;
+
+ //! constructor using sizes and resolution
+ SystemBase(Projection_ptr projection);
+
+ //! Copy constructor
+ SystemBase(const SystemBase &other) = delete;
+
+ //! Move constructor
+ SystemBase(SystemBase &&other) noexcept = default;
+
+ //! Destructor
+ virtual ~SystemBase() noexcept = default;
+
+ //! Copy assignment operator
+ SystemBase& operator=(const SystemBase &other) = delete;
+
+ //! Move assignment operator
+ SystemBase& operator=(SystemBase &&other) noexcept = default;
+
+ /**
+ * Materials can only be moved. This is to assure exclusive
+ * ownership of any material by this system
+ */
+ void add_material(Material_ptr mat);
+
+ FullResponse_t evaluate_stress_tangent(StrainField_t & F);
+
+ void initialise();
+ protected:
+ //! make sure that every pixel is assigned to one and only one material
+ void check_material_coverage();
+
+ const Ccoord & resolutions;
+ const Rcoord & lengths;
+ FieldCollection_t fields{};
+ StrainField_t & F;
+ StressField_t & P;
+ //! Tangent field migth not even be required; so this is a
+ //! pointer instead of a ref
+ TangentField_t * K_ptr;
+ std::vector<Material_ptr> materials{};
+ Projection_ptr projection;
+ bool is_initialised{false};
+ private:
};
} // muSpectre
#endif /* SYSTEM_BASE_H */
diff --git a/tests/test_system_base.cc b/tests/test_system_base.cc
new file mode 100644
index 0000000..8870c73
--- /dev/null
+++ b/tests/test_system_base.cc
@@ -0,0 +1,80 @@
+/**
+ * file test_system_base.cc
+ *
+ * @author Till Junge <till.junge@epfl.ch>
+ *
+ * @date 14 Dec 2017
+ *
+ * @brief Tests for the basic system class
+ *
+ * @section LICENCE
+ *
+ * Copyright (C) 2017 Till Junge
+ *
+ * µSpectre is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License as
+ * published by the Free Software Foundation, either version 3, or (at
+ * your option) any later version.
+ *
+ * µSpectre is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with GNU Emacs; see the file COPYING. If not, write to the
+ * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+ * Boston, MA 02111-1307, USA.
+ */
+
+
+#include <boost/mpl/list.hpp>
+
+#include "tests.hh"
+#include "common/common.hh"
+#include "system/system_base.hh"
+#include "fft/projection_finite_strain_fast.hh"
+#include "fft/fftw_engine.hh"
+
+
+namespace muSpectre {
+
+ BOOST_AUTO_TEST_SUITE(system_base);
+ template <Dim_t DimS>
+ struct Sizes {
+ };
+ template<>
+ struct Sizes<twoD> {
+ constexpr static Ccoord_t<twoD> get_resolution() {
+ return Ccoord_t<twoD>{3, 5};}
+ constexpr static Rcoord_t<twoD> get_lengths() {
+ return Rcoord_t<twoD>{3.4, 5.8};}
+ };
+ template<>
+ struct Sizes<threeD> {
+ constexpr static Ccoord_t<threeD> get_resolution() {
+ return Ccoord_t<threeD>{3, 5, 7};}
+ constexpr static Rcoord_t<threeD> get_lengths() {
+ return Rcoord_t<threeD>{3.4, 5.8, 6.7};}
+ };
+
+ template <Dim_t DimS, Dim_t DimM>
+ struct SystemBaseFixture: SystemBase<DimS, DimM> {
+ SystemBaseFixture()
+ :fft_engine{Sizes<DimS>::get_resolution(), Sizes<DimS>::get_lengths()},
+ proj{fft_engine},
+ SystemBase<DimS, DimM>{std::move(proj)}{}
+
+ FFTW_Engine<DimS, DimM> fft_engine;
+ ProjectionFiniteStrainFast<DimS, DimM> proj;
+ };
+
+ using fixlist = boost::mpl::list<SystemBaseFixture<twoD, twoD>,
+ SystemBaseFixture<threeD, threeD>>;
+
+ BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) {
+ }
+
+ BOOST_AUTO_TEST_SUITE_END();
+
+} // muSpectre