Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F87705932
tensor_algebra.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
Mon, Oct 14, 09:35
Size
9 KB
Mime Type
text/x-c
Expires
Wed, Oct 16, 09:35 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
21495102
Attached To
rMUSPECTRE µSpectre
tensor_algebra.hh
View Options
/**
* @file tensor_algebra.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 05 Nov 2017
*
* @brief collection of compile-time quantities and algrebraic functions for
* tensor operations
*
* 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 TENSOR_ALGEBRA_H
#define TENSOR_ALGEBRA_H
#include "common/T4_map_proxy.hh"
#include "common/common.hh"
#include "common/eigen_tools.hh"
#include <Eigen/Dense>
#include <unsupported/Eigen/CXX11/Tensor>
#include <type_traits>
namespace
muSpectre
{
namespace
Tensors
{
//! second-order tensor representation
template
<
Dim_t
dim
>
using
Tens2_t
=
Eigen
::
TensorFixedSize
<
Real
,
Eigen
::
Sizes
<
dim
,
dim
>>
;
//! fourth-order tensor representation
template
<
Dim_t
dim
>
using
Tens4_t
=
Eigen
::
TensorFixedSize
<
Real
,
Eigen
::
Sizes
<
dim
,
dim
,
dim
,
dim
>>
;
//----------------------------------------------------------------------------//
//! compile-time second-order identity
template
<
Dim_t
dim
>
constexpr
inline
Tens2_t
<
dim
>
I2
()
{
Tens2_t
<
dim
>
T
;
using
Mat_t
=
Eigen
::
Matrix
<
Real
,
dim
,
dim
>
;
Eigen
::
Map
<
Mat_t
>
(
&
T
(
0
,
0
))
=
Mat_t
::
Identity
();
return
T
;
}
/* ---------------------------------------------------------------------- */
//! Check whether a given expression represents a Tensor specified order
template
<
class
T
,
Dim_t
order
>
struct
is_tensor
{
//! evaluated test
constexpr
static
bool
value
=
(
std
::
is_convertible
<
T
,
Eigen
::
Tensor
<
Real
,
order
>>::
value
||
std
::
is_convertible
<
T
,
Eigen
::
Tensor
<
Int
,
order
>>::
value
||
std
::
is_convertible
<
T
,
Eigen
::
Tensor
<
Complex
,
order
>>::
value
);
};
/* ---------------------------------------------------------------------- */
/** compile-time outer tensor product as defined by Curnier
* R_ijkl = A_ij.B_klxx
* 0123 01 23
*/
template
<
Dim_t
dim
,
typename
T1
,
typename
T2
>
constexpr
inline
decltype
(
auto
)
outer
(
T1
&&
A
,
T2
&&
B
)
{
// Just make sure that the right type of parameters have been given
constexpr
Dim_t
order
{
2
};
static_assert
(
is_tensor
<
T1
,
order
>::
value
,
"T1 needs to be convertible to a second order Tensor"
);
static_assert
(
is_tensor
<
T2
,
order
>::
value
,
"T2 needs to be convertible to a second order Tensor"
);
// actual function
std
::
array
<
Eigen
::
IndexPair
<
Dim_t
>
,
0
>
dims
{};
return
A
.
contract
(
B
,
dims
);
}
/* ---------------------------------------------------------------------- */
/** compile-time underlined outer tensor product as defined by Curnier
* R_ijkl = A_ik.B_jlxx
* 0123 02 13
* 0213 01 23 <- this defines the shuffle order
*/
template
<
Dim_t
dim
,
typename
T1
,
typename
T2
>
constexpr
inline
decltype
(
auto
)
outer_under
(
T1
&&
A
,
T2
&&
B
)
{
constexpr
size_t
order
{
4
};
return
outer
<
dim
>
(
A
,
B
).
shuffle
(
std
::
array
<
Dim_t
,
order
>
{{
0
,
2
,
1
,
3
}});
}
/* ---------------------------------------------------------------------- */
/** compile-time overlined outer tensor product as defined by Curnier
* R_ijkl = A_il.B_jkxx
* 0123 03 12
* 0231 01 23 <- this defines the shuffle order
*/
template
<
Dim_t
dim
,
typename
T1
,
typename
T2
>
constexpr
inline
decltype
(
auto
)
outer_over
(
T1
&&
A
,
T2
&&
B
)
{
constexpr
size_t
order
{
4
};
return
outer
<
dim
>
(
A
,
B
).
shuffle
(
std
::
array
<
Dim_t
,
order
>
{{
0
,
2
,
3
,
1
}});
}
//! compile-time fourth-order symmetrising identity
template
<
Dim_t
dim
>
constexpr
inline
Tens4_t
<
dim
>
I4S
()
{
auto
I
=
I2
<
dim
>
();
return
0.5
*
(
outer_under
<
dim
>
(
I
,
I
)
+
outer_over
<
dim
>
(
I
,
I
));
}
}
// Tensors
namespace
Matrices
{
//! second-order tensor representation
template
<
Dim_t
dim
>
using
Tens2_t
=
Eigen
::
Matrix
<
Real
,
dim
,
dim
>
;
//! fourth-order tensor representation
template
<
Dim_t
dim
>
using
Tens4_t
=
T4Mat
<
Real
,
dim
>
;
//----------------------------------------------------------------------------//
//! compile-time second-order identity
template
<
Dim_t
dim
>
constexpr
inline
Tens2_t
<
dim
>
I2
()
{
return
Tens2_t
<
dim
>::
Identity
();
}
/* ---------------------------------------------------------------------- */
/** compile-time outer tensor product as defined by Curnier
* R_ijkl = A_ij.B_klxx
* 0123 01 23
*/
template
<
typename
T1
,
typename
T2
>
constexpr
inline
decltype
(
auto
)
outer
(
T1
&&
A
,
T2
&&
B
)
{
// Just make sure that the right type of parameters have been given
constexpr
Dim_t
dim
{
EigenCheck
::
tensor_dim
<
T1
>::
value
};
static_assert
((
dim
==
EigenCheck
::
tensor_dim
<
T2
>::
value
),
"A and B do not have the same dimension"
);
Tens4_t
<
dim
>
product
;
for
(
Dim_t
i
=
0
;
i
<
dim
;
++
i
)
{
for
(
Dim_t
j
=
0
;
j
<
dim
;
++
j
)
{
for
(
Dim_t
k
=
0
;
k
<
dim
;
++
k
)
{
for
(
Dim_t
l
=
0
;
l
<
dim
;
++
l
)
{
get
(
product
,
i
,
j
,
k
,
l
)
=
A
(
i
,
j
)
*
B
(
k
,
l
);
}
}
}
}
return
product
;
}
/* ---------------------------------------------------------------------- */
/** compile-time underlined outer tensor product as defined by Curnier
* R_ijkl = A_ik.B_jlxx
* 0123 02 13
* 0213 01 23 <- this defines the shuffle order
*/
template
<
typename
T1
,
typename
T2
>
constexpr
inline
decltype
(
auto
)
outer_under
(
T1
&&
A
,
T2
&&
B
)
{
// Just make sure that the right type of parameters have been given
constexpr
Dim_t
dim
{
EigenCheck
::
tensor_dim
<
T1
>::
value
};
static_assert
((
dim
==
EigenCheck
::
tensor_dim
<
T2
>::
value
),
"A and B do not have the same dimension"
);
Tens4_t
<
dim
>
product
;
for
(
Dim_t
i
=
0
;
i
<
dim
;
++
i
)
{
for
(
Dim_t
j
=
0
;
j
<
dim
;
++
j
)
{
for
(
Dim_t
k
=
0
;
k
<
dim
;
++
k
)
{
for
(
Dim_t
l
=
0
;
l
<
dim
;
++
l
)
{
get
(
product
,
i
,
j
,
k
,
l
)
=
A
(
i
,
k
)
*
B
(
j
,
l
);
}
}
}
}
return
product
;
}
/* ---------------------------------------------------------------------- */
/** compile-time overlined outer tensor product as defined by Curnier
* R_ijkl = A_il.B_jkxx
* 0123 03 12
* 0231 01 23 <- this defines the shuffle order
*/
template
<
typename
T1
,
typename
T2
>
constexpr
inline
decltype
(
auto
)
outer_over
(
T1
&&
A
,
T2
&&
B
)
{
// Just make sure that the right type of parameters have been given
constexpr
Dim_t
dim
{
EigenCheck
::
tensor_dim
<
T1
>::
value
};
static_assert
((
dim
==
EigenCheck
::
tensor_dim
<
T2
>::
value
),
"A and B do not have the same dimension"
);
Tens4_t
<
dim
>
product
;
for
(
Dim_t
i
=
0
;
i
<
dim
;
++
i
)
{
for
(
Dim_t
j
=
0
;
j
<
dim
;
++
j
)
{
for
(
Dim_t
k
=
0
;
k
<
dim
;
++
k
)
{
for
(
Dim_t
l
=
0
;
l
<
dim
;
++
l
)
{
get
(
product
,
i
,
j
,
k
,
l
)
=
A
(
i
,
l
)
*
B
(
j
,
k
);
}
}
}
}
return
product
;
}
/**
* Standart tensor multiplication
*/
template
<
typename
T4
,
typename
T2
>
constexpr
inline
decltype
(
auto
)
tensmult
(
T4
&&
A
,
T2
&&
B
)
{
constexpr
Dim_t
dim
{
EigenCheck
::
tensor_4_dim
<
T4
>::
value
};
static_assert
((
dim
==
EigenCheck
::
tensor_dim
<
T2
>::
value
),
"Dimensionality check failed. Expects A to be a fourth-"
"order tensor of dimension dim and B to be a second-order "
"Tensor of same dimension"
);
Tens2_t
<
dim
>
result
;
result
.
setZero
();
for
(
Dim_t
i
=
0
;
i
<
dim
;
++
i
)
{
for
(
Dim_t
j
=
0
;
j
<
dim
;
++
j
)
{
for
(
Dim_t
k
=
0
;
k
<
dim
;
++
k
)
{
for
(
Dim_t
l
=
0
;
l
<
dim
;
++
l
)
{
result
(
i
,
j
)
+=
get
(
A
,
i
,
j
,
k
,
l
)
*
B
(
k
,
l
);
}
}
}
}
return
result
;
}
//! compile-time fourth-order tracer
template
<
Dim_t
dim
>
constexpr
inline
Tens4_t
<
dim
>
Itrac
()
{
auto
I
=
I2
<
dim
>
();
return
outer
(
I
,
I
);
}
//! compile-time fourth-order identity
template
<
Dim_t
dim
>
constexpr
inline
Tens4_t
<
dim
>
Iiden
()
{
auto
I
=
I2
<
dim
>
();
return
outer_under
(
I
,
I
);
}
//! compile-time fourth-order transposer
template
<
Dim_t
dim
>
constexpr
inline
Tens4_t
<
dim
>
Itrns
()
{
auto
I
=
I2
<
dim
>
();
return
outer_over
(
I
,
I
);
}
//! compile-time fourth-order symmetriser
template
<
Dim_t
dim
>
constexpr
inline
Tens4_t
<
dim
>
Isymm
()
{
auto
I
=
I2
<
dim
>
();
return
0.5
*
(
outer_under
(
I
,
I
)
+
outer_over
(
I
,
I
));
}
}
// Matrices
}
// muSpectre
#endif
/* TENSOR_ALGEBRA_H */
Event Timeline
Log In to Comment