Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F90648113
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
Sun, Nov 3, 13:48
Size
10 KB
Mime Type
text/x-c++
Expires
Tue, Nov 5, 13:48 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
22114452
Attached To
rTAMAAS tamaas
influence.hh
View Options
/**
* @file
* LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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
>
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
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
>
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
>
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
>
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
>
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
>
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
>
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
>
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
>
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
>
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
>
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
>
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
>
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
>
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
>
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