Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85281172
field.hh
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Sep 28, 00:00
Size
40 KB
Mime Type
text/x-c
Expires
Mon, Sep 30, 00:00 (2 d)
Engine
blob
Format
Raw Data
Handle
21154287
Attached To
rMUSPECTRE µSpectre
field.hh
View Options
/**
* @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 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
* 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
;
/* ---------------------------------------------------------------------- */
/**
* 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
/* ---------------------------------------------------------------------- */
/**
* 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
));}
/**
* 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
();
/**
* 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
();
/**
* 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
:
};
/* ---------------------------------------------------------------------- */
/**
* 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
));}
/**
* 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
();
/**
* 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
();
/**
* 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
>
;
};
//! 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 */
Event Timeline
Log In to Comment