Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F84095127
influence.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
Fri, Sep 20, 17:06
Size
11 KB
Mime Type
text/x-c++
Expires
Sun, Sep 22, 17:06 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20939131
Attached To
rTAMAAS tamaas
influence.hh
View Options
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2022 Lucas Frérot
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef INFLUENCE_HH
#define INFLUENCE_HH
/* -------------------------------------------------------------------------- */
#include "loop.hh"
#include "static_types.hh"
#include <type_traits>
/* -------------------------------------------------------------------------- */
namespace
tamaas
{
/* -------------------------------------------------------------------------- */
namespace
influence
{
namespace
{
const
Complex
I
(
0
,
1
);
/// < imaginary unit
constexpr
inline
Real
sign
(
bool
upper
)
{
return
(
upper
)
?
1.
:
-
1.
;
}
}
// namespace
template
<
bool
conjugate
,
UInt
dim_q
>
// dim + 1 to help dim template deduction
inline
Vector
<
Complex
,
dim_q
+
1
>
__device__
__host__
computeD
(
const
VectorProxy
<
const
Real
,
dim_q
>&
q
)
{
constexpr
Real
s
=
-
sign
(
conjugate
);
const
auto
q_norm
=
q
.
l2norm
();
return
{{{
Complex
{
0
,
s
*
q
(
0
)
/
q_norm
},
Complex
{
0
,
s
*
q
(
1
)
/
q_norm
},
Complex
{
1
,
0
}}}};
}
template
<
UInt
dim_q
>
// dim + 1 to help dim template deduction
__device__
__host__
inline
Vector
<
Complex
,
dim_q
+
1
>
computeD2
(
const
VectorProxy
<
const
Real
,
dim_q
>&
q
)
{
const
auto
q_norm
=
(
q
.
l2norm
()
<
1e-16
)
?
1
:
q
.
l2norm
();
return
{
{{
Complex
{
0
,
q
(
0
)
/
q_norm
},
Complex
{
0
,
q
(
1
)
/
q_norm
},
Complex
{
0
,
0
}}}};
}
/* -------------------------------------------------------------------------- */
/// Functor to apply Hooke's tensor
template
<
UInt
dim
>
struct
ElasticHelper
{
ElasticHelper
(
Real
mu
,
Real
nu
)
:
mu
(
mu
),
nu
(
nu
),
lambda
(
2
*
mu
*
nu
/
(
1
-
2
*
nu
))
{}
template
<
typename
DT
,
typename
ST
>
__device__
__host__
Matrix
<
std
::
remove_cv_t
<
DT
>
,
dim
,
dim
>
operator
()(
const
StaticMatrix
<
DT
,
ST
,
dim
,
dim
>&
gradu
)
const
{
auto
trace
=
gradu
.
trace
();
Matrix
<
std
::
remove_cv_t
<
DT
>
,
dim
,
dim
>
sigma
;
for
(
UInt
i
=
0
;
i
<
dim
;
++
i
)
for
(
UInt
j
=
0
;
j
<
dim
;
++
j
)
sigma
(
i
,
j
)
=
(
i
==
j
)
*
lambda
*
trace
+
mu
*
(
gradu
(
i
,
j
)
+
gradu
(
j
,
i
));
return
sigma
;
}
template
<
typename
DT
,
typename
ST
>
__device__
__host__
SymMatrix
<
std
::
remove_cv_t
<
DT
>
,
dim
>
operator
()(
const
StaticSymMatrix
<
DT
,
ST
,
dim
>&
eps
)
const
{
SymMatrix
<
std
::
remove_cv_t
<
DT
>
,
dim
>
sigma
;
auto
trace
=
eps
.
trace
();
for
(
UInt
i
=
0
;
i
<
dim
;
++
i
)
sigma
(
i
)
=
lambda
*
trace
+
2
*
mu
*
eps
(
i
);
for
(
UInt
i
=
dim
;
i
<
voigt_size
<
dim
>::
value
;
++
i
)
sigma
(
i
)
=
2
*
mu
*
eps
(
i
);
return
sigma
;
}
template
<
typename
DT
,
typename
ST
>
__device__
__host__
Matrix
<
std
::
remove_cv_t
<
DT
>
,
dim
,
dim
>
inverse
(
const
StaticMatrix
<
DT
,
ST
,
dim
,
dim
>&
sigma
)
const
{
Matrix
<
std
::
remove_cv_t
<
DT
>
,
dim
,
dim
>
epsilon
;
auto
trace
=
sigma
.
trace
();
auto
c1
=
1.
/
(
2
*
mu
);
auto
c2
=
-
lambda
/
(
3
*
lambda
+
2
*
mu
);
for
(
UInt
i
=
0
;
i
<
dim
;
++
i
)
for
(
UInt
j
=
0
;
j
<
dim
;
++
j
)
epsilon
(
i
,
j
)
=
c1
*
(
sigma
(
i
,
j
)
+
c2
*
trace
*
(
i
==
j
));
}
const
Real
mu
,
nu
,
lambda
;
};
/* -------------------------------------------------------------------------- */
/// Class for the Kelvin tensor
template
<
UInt
dim
,
UInt
derivative_order
>
class
Kelvin
;
/**
* @brief Kelvin base tensor
* See arXiv:1811.11558 for details.
*/
template
<>
class
Kelvin
<
3
,
0
>
{
protected
:
static
constexpr
UInt
dim
=
3
;
static
constexpr
UInt
order
=
0
;
public
:
Kelvin
(
Real
mu
,
Real
nu
)
:
mu
(
mu
),
b
(
4
*
(
1
-
nu
))
{}
template
<
bool
upper
,
bool
apply_q_power
=
false
,
typename
ST
>
__device__
__host__
inline
Vector
<
Complex
,
dim
>
applyU0
(
const
VectorProxy
<
const
Real
,
dim
-
1
>&
q
,
const
StaticVector
<
Complex
,
ST
,
dim
>&
f
)
const
{
// upper is ignored here
auto
d2
=
computeD2
(
q
);
auto
q_power
=
(
apply_q_power
)
?
1.
/
q
.
l2norm
()
:
1.
;
d2
*=
d2
.
dot
(
f
);
d2
+=
b
*
f
;
d2
(
2
)
-=
f
(
2
);
d2
*=
q_power
/
(
2.
*
mu
*
b
);
return
d2
;
}
template
<
bool
upper
,
bool
apply_q_power
=
false
,
typename
ST
>
__device__
__host__
inline
Vector
<
Complex
,
dim
>
applyU1
(
const
VectorProxy
<
const
Real
,
dim
-
1
>&
q
,
const
StaticVector
<
Complex
,
ST
,
dim
>&
f
)
const
{
auto
d
=
computeD
<
upper
>
(
q
);
auto
q_power
=
(
apply_q_power
)
?
1.
/
q
.
l2norm
()
:
1.
;
d
*=
d
.
dot
(
f
);
d
*=
q_power
*
sign
(
upper
)
/
(
2.
*
mu
*
b
);
return
d
;
}
protected
:
const
Real
mu
,
b
;
};
/* -------------------------------------------------------------------------- */
/// Kelvin first derivative
template
<>
class
Kelvin
<
3
,
1
>
:
protected
Kelvin
<
3
,
0
>
{
using
parent
=
Kelvin
<
3
,
0
>
;
protected
:
static
constexpr
UInt
dim
=
parent
::
dim
;
static
constexpr
UInt
order
=
parent
::
order
+
1
;
public
:
using
parent
::
parent
;
template
<
bool
upper
,
bool
apply_q_power
=
false
,
typename
ST
>
__device__
__host__
inline
Vector
<
Complex
,
dim
>
applyU0
(
const
VectorProxy
<
const
Real
,
dim
-
1
>&
q
,
const
StaticMatrix
<
Complex
,
ST
,
dim
,
dim
>&
f
)
const
{
// apply_q_power is ignored here
auto
tmp
=
f
*
computeD
<
upper
>
(
q
);
auto
res
=
Kelvin
<
3
,
0
>::
applyU0
<
false
>
(
q
,
tmp
);
res
*=
-
sign
(
upper
);
Vector
<
Real
,
dim
>
e3
{{{
0
,
0
,
1
}}};
tmp
.
template
mul
<
false
>
(
f
,
e3
);
res
+=
Kelvin
<
3
,
0
>::
applyU1
<
upper
>
(
q
,
tmp
);
return
res
;
}
template
<
bool
upper
,
bool
apply_q_power
=
false
,
typename
ST
>
inline
Vector
<
Complex
,
dim
>
__device__
__host__
applyU1
(
const
VectorProxy
<
const
Real
,
dim
-
1
>&
q
,
const
StaticMatrix
<
Complex
,
ST
,
dim
,
dim
>&
f
)
const
{
// apply_q_power is ignored here
auto
tmp
=
f
*
computeD
<
upper
>
(
q
);
auto
res
=
Kelvin
<
3
,
0
>::
applyU1
<
upper
>
(
q
,
tmp
);
res
*=
-
sign
(
upper
);
return
res
;
}
};
/* -------------------------------------------------------------------------- */
/// Kelvin second derivative
template
<>
class
Kelvin
<
3
,
2
>
:
protected
Kelvin
<
3
,
1
>
{
using
parent
=
Kelvin
<
3
,
1
>
;
static
constexpr
UInt
dim
=
parent
::
dim
;
static
constexpr
UInt
order
=
parent
::
order
+
1
;
public
:
using
parent
::
parent
;
template
<
bool
upper
,
bool
apply_q_power
=
false
,
typename
ST
>
__device__
__host__
inline
Matrix
<
Complex
,
dim
,
dim
>
applyU0
(
const
VectorProxy
<
const
Real
,
dim
-
1
>&
q
,
const
StaticMatrix
<
Complex
,
ST
,
dim
,
dim
>&
f
)
const
{
auto
tmp
=
Kelvin
<
3
,
1
>::
applyU0
<
upper
>
(
q
,
f
);
auto
q_power
=
(
apply_q_power
)
?
q
.
l2norm
()
:
1.
;
tmp
*=
-
sign
(
upper
);
auto
res
=
outer
(
computeD
<
upper
>
(
q
),
tmp
);
Vector
<
Real
,
dim
>
e3
{{{
0
,
0
,
1
}}};
res
+=
outer
(
e3
,
Kelvin
<
3
,
1
>::
applyU1
<
upper
>
(
q
,
f
));
res
*=
-
1.
*
q_power
;
// << -1 for grad_x
return
res
;
}
template
<
bool
upper
,
bool
apply_q_power
=
false
,
typename
ST
>
__device__
__host__
inline
Matrix
<
Complex
,
dim
,
dim
>
applyU1
(
const
VectorProxy
<
const
Real
,
dim
-
1
>&
q
,
const
StaticMatrix
<
Complex
,
ST
,
dim
,
dim
>&
f
)
const
{
auto
tmp
=
Kelvin
<
3
,
1
>::
applyU1
<
upper
>
(
q
,
f
);
auto
q_power
=
(
apply_q_power
)
?
q
.
l2norm
()
:
1.
;
auto
res
=
outer
(
computeD
<
upper
>
(
q
),
tmp
);
res
*=
-
1.
*
-
sign
(
upper
)
*
q_power
;
// << -1 for grad_x
return
res
;
}
template
<
typename
ST
>
__device__
__host__
inline
Matrix
<
Complex
,
dim
,
dim
>
applyDiscontinuityTerm
(
const
StaticMatrix
<
Complex
,
ST
,
dim
,
dim
>&
f
)
const
{
Vector
<
Real
,
dim
>
e3
{{{
0
,
0
,
1
}}};
auto
res
=
(
f
*
e3
);
res
*=
-
this
->
b
;
res
(
2
)
+=
2.
*
f
(
2
,
2
);
e3
(
2
)
*=
-
1
/
(
this
->
mu
*
this
->
b
);
// << -1 for grad_x
return
outer
(
e3
,
res
);
}
};
/* -------------------------------------------------------------------------- */
/// Class for the Boussinesq tensor
template
<
UInt
dim
,
UInt
derivative_order
>
class
Boussinesq
;
template
<>
class
Boussinesq
<
3
,
0
>
{
protected
:
static
constexpr
UInt
dim
=
3
;
static
constexpr
UInt
order
=
0
;
public
:
/// Constructor
Boussinesq
(
Real
mu
,
Real
nu
)
:
mu
(
mu
),
nu
(
nu
),
lambda
(
2
*
mu
*
nu
/
(
1
-
2
*
nu
))
{}
template
<
bool
apply_q_power
=
false
,
typename
ST
>
__device__
__host__
inline
Vector
<
Complex
,
dim
>
applyU0
(
const
StaticVector
<
Complex
,
ST
,
dim
>&
t
,
const
VectorProxy
<
const
Real
,
dim
-
1
>&
q
)
const
{
// This guy has the property Bt(q) = B(-q). Since we know we recieve -q as
// an argument, we can just apply B without transpose
auto
dplus
=
computeD
<
true
>
(
q
),
dminus
=
computeD
<
false
>
(
q
),
d2
=
computeD2
(
q
);
auto
q_power
=
(
apply_q_power
)
?
1.
/
q
.
l2norm
()
:
1.
;
auto
res
=
2.
*
t
;
res
+=
((
1.
-
2.
*
nu
)
*
dminus
.
dot
(
t
))
*
dplus
;
res
+=
d2
.
dot
(
t
)
*
d2
;
res
(
2
)
-=
t
(
2
);
res
*=
q_power
/
(
2.
*
mu
);
return
res
;
}
template
<
bool
apply_q_power
=
false
,
typename
ST
>
__device__
__host__
inline
Vector
<
Complex
,
dim
>
applyU1
(
const
StaticVector
<
Complex
,
ST
,
dim
>&
t
,
const
VectorProxy
<
const
Real
,
dim
-
1
>&
q
)
const
{
// This guy does not have the property that Bt(q) = B(-q), so we have to
// apply Bt(q) for real. This makes us correct (in the line below) for the
// fact that we recieve -q as an argument.
auto
dplus
=
computeD
<
false
>
(
q
);
auto
q_power
=
(
apply_q_power
)
?
1.
/
q
.
l2norm
()
:
1.
;
auto
res
=
q_power
*
(
dplus
.
dot
(
t
)
/
(
2.
*
mu
))
*
dplus
;
return
res
;
}
protected
:
const
Real
mu
,
nu
,
lambda
;
};
/* -------------------------------------------------------------------------- */
/// Boussinesq first gradient
template
<>
class
Boussinesq
<
3
,
1
>
:
protected
Boussinesq
<
3
,
0
>
{
protected
:
using
parent
=
Boussinesq
<
3
,
0
>
;
static
constexpr
UInt
dim
=
parent
::
dim
;
static
constexpr
UInt
order
=
parent
::
order
+
1
;
public
:
using
parent
::
parent
;
template
<
bool
apply_q_power
=
false
,
typename
ST
>
__device__
__host__
inline
Matrix
<
Complex
,
dim
,
dim
>
applyU0
(
const
StaticVector
<
Complex
,
ST
,
dim
>&
t
,
const
VectorProxy
<
const
Real
,
dim
-
1
>&
q
)
const
{
// taking into account K(-q) convolution
auto
dbar
=
computeD
<
false
>
(
q
);
Vector
<
Real
,
dim
>
e3
{{{
0
,
0
,
1
}}};
auto
res
=
outer
(
dbar
,
Boussinesq
<
3
,
0
>::
applyU0
(
t
,
q
));
res
*=
-
1.
;
res
+=
outer
(
e3
,
Boussinesq
<
3
,
0
>::
applyU1
(
t
,
q
));
return
res
;
}
template
<
bool
apply_q_power
=
false
,
typename
ST
>
__device__
__host__
inline
Matrix
<
Complex
,
dim
,
dim
>
applyU1
(
const
StaticVector
<
Complex
,
ST
,
dim
>&
t
,
const
VectorProxy
<
const
Real
,
dim
-
1
>&
q
)
const
{
// taking into account K(-q) convolution
auto
dbar
=
computeD
<
false
>
(
q
);
auto
res
=
outer
(
dbar
,
Boussinesq
<
3
,
0
>::
applyU1
(
t
,
q
));
res
*=
-
1
;
return
res
;
}
};
}
// namespace influence
/* -------------------------------------------------------------------------- */
}
// namespace tamaas
#endif
// INFLUENCE_HH
Event Timeline
Log In to Comment