Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F62398003
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, May 12, 22:18
Size
8 KB
Mime Type
text/x-c++
Expires
Tue, May 14, 22:18 (2 d)
Engine
blob
Format
Raw Data
Handle
17642606
Attached To
rTAMAAS tamaas
influence.hh
View Options
/**
* @file
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Tamaas is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Tamaas 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 Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "loop.hh"
#include "static_types.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
namespace
influence
{
namespace
{
const
Complex
I
(
0
,
1
);
/* < imaginary unit */
constexpr
Real
sign
(
bool
upper
)
{
return
(
upper
)
?
1.
:
-
1.
;
}
}
template
<
bool
conjugate
,
UInt
dim
>
// dim + 1 to help dim template deduction
inline
Vector
<
Complex
,
dim
+
1
>
computeD
(
VectorProxy
<
const
Real
,
dim
>&
q
)
{
constexpr
Real
s
=
-
sign
(
conjugate
);
const
auto
q_norm
=
q
.
l2norm
();
return
{{{
s
*
I
*
q
(
0
)
/
q_norm
,
s
*
I
*
q
(
1
)
/
q_norm
,
1
}}};
}
template
<
UInt
dim
>
// dim + 1 to help dim template deduction
inline
Vector
<
Complex
,
dim
+
1
>
computeD2
(
VectorProxy
<
const
Real
,
dim
>&
q
)
{
const
auto
q_norm
=
q
.
l2norm
();
return
{{{
I
*
q
(
0
)
/
q_norm
,
I
*
q
(
1
)
/
q_norm
,
0
}}};
}
/// Class for the Kelvin tensor
template
<
UInt
dim
,
UInt
derivative_order
>
class
Kelvin
;
/// See details in notebook and Marc Bonnet's notes
template
<>
class
Kelvin
<
3
,
0
>
{
protected
:
constexpr
static
UInt
dim
=
3
;
public
:
Kelvin
(
Real
mu
,
Real
nu
)
:
mu
(
mu
),
b
(
4
*
(
1
-
nu
))
{}
template
<
bool
upper
,
typename
ST
>
inline
Vector
<
Complex
,
dim
>
applyU0
(
VectorProxy
<
const
Real
,
dim
-
1
>&
q
,
StaticVector
<
Complex
,
ST
,
dim
>&
f
)
const
;
template
<
bool
upper
,
typename
ST
>
inline
Vector
<
Complex
,
dim
>
applyU1
(
VectorProxy
<
const
Real
,
dim
-
1
>&
q
,
StaticVector
<
Complex
,
ST
,
dim
>&
f
)
const
;
protected
:
const
Real
mu
,
b
;
};
/* -------------------------------------------------------------------------- */
/* Base Kelvin tensor */
/* -------------------------------------------------------------------------- */
template
<
bool
upper
,
typename
ST
>
//< this is ignored
inline
Vector
<
Complex
,
Kelvin
<
3
,
0
>::
dim
>
Kelvin
<
3
,
0
>::
applyU0
(
VectorProxy
<
const
Real
,
dim
-
1
>&
q
,
StaticVector
<
Complex
,
ST
,
dim
>&
f
)
const
{
auto
d2
{
computeD2
(
q
)};
auto
tmp
{
d2
.
dot
(
f
)
*
d2
};
tmp
+=
b
*
f
;
tmp
(
2
)
-=
f
(
2
);
tmp
*=
1.
/
(
2.
*
mu
*
q
.
l2norm
()
*
b
);
return
tmp
;
}
template
<
bool
upper
,
typename
ST
>
inline
Vector
<
Complex
,
Kelvin
<
3
,
0
>::
dim
>
Kelvin
<
3
,
0
>::
applyU1
(
VectorProxy
<
const
Real
,
dim
-
1
>&
q
,
StaticVector
<
Complex
,
ST
,
dim
>&
f
)
const
{
auto
d
{
computeD
<
upper
>
(
q
)};
auto
tmp
{
d
.
dot
(
f
)
*
d
};
tmp
*=
sign
(
upper
)
/
(
2.
*
mu
*
q
.
l2norm
()
*
b
);
return
tmp
;
}
/* -------------------------------------------------------------------------- */
/* Kelvin first derivative */
/* -------------------------------------------------------------------------- */
template
<>
class
Kelvin
<
3
,
1
>
:
Kelvin
<
3
,
0
>
{
using
parent
=
Kelvin
<
3
,
0
>
;
static
constexpr
UInt
dim
=
parent
::
dim
;
public
:
using
parent
::
parent
;
template
<
bool
upper
,
typename
ST
>
inline
Vector
<
Complex
,
dim
>
applyU0
(
VectorProxy
<
const
Real
,
dim
-
1
>&
q
,
StaticMatrix
<
Complex
,
ST
,
dim
,
dim
>&
f
)
const
;
template
<
bool
upper
,
typename
ST
>
inline
Vector
<
Complex
,
dim
>
applyU1
(
VectorProxy
<
const
Real
,
dim
-
1
>&
q
,
StaticMatrix
<
Complex
,
ST
,
dim
,
dim
>&
f
)
const
;
};
template
<
bool
upper
,
typename
ST
>
inline
Vector
<
Complex
,
Kelvin
<
3
,
1
>::
dim
>
Kelvin
<
3
,
1
>::
applyU0
(
VectorProxy
<
const
Real
,
dim
-
1
>&
q
,
StaticMatrix
<
Complex
,
ST
,
dim
,
dim
>&
f
)
const
{
auto
d_apply
{
f
*
computeD
<
false
>
(
q
)};
auto
tmp
{
f
*
computeD
<
upper
>
(
q
)};
auto
res
{
Kelvin
<
3
,
0
>::
applyU0
<
false
>
(
q
,
tmp
)};
res
*=
-
sign
(
upper
);
res
+=
Kelvin
<
3
,
0
>::
applyU1
<
upper
>
(
q
,
tmp
);
res
*=
q
.
l2norm
();
return
res
;
}
template
<
bool
upper
,
typename
ST
>
inline
Vector
<
Complex
,
Kelvin
<
3
,
1
>::
dim
>
Kelvin
<
3
,
1
>::
applyU1
(
VectorProxy
<
const
Real
,
dim
-
1
>&
q
,
StaticMatrix
<
Complex
,
ST
,
dim
,
dim
>&
f
)
const
{
auto
tmp
{
f
*
computeD
<
upper
>
(
q
)};
auto
res
{
Kelvin
<
3
,
0
>::
applyU1
<
upper
>
(
q
,
tmp
)};
res
*=
-
sign
(
upper
)
*
q
.
l2norm
();
return
res
;
}
/* -------------------------------------------------------------------------- */
/* Implementation */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
class
KelvinIntegrator
{
static
inline
Real
int_exp
(
Real
a
,
Real
z
)
{
return
std
::
exp
(
a
*
z
)
/
a
;
}
static
inline
Real
int_exp_z
(
Real
a
,
Real
z
)
{
const
Real
a_inv
=
1.
/
a
;
return
std
::
exp
(
a
*
z
)
*
(
z
-
a_inv
)
*
a_inv
;
}
static
inline
Real
int_exp_zz
(
Real
a
,
Real
z
)
{
const
Real
a_inv
=
1.
/
a
;
return
std
::
exp
(
a
*
z
)
*
(
z
*
z
-
2
*
a_inv
*
z
+
2
*
a_inv
*
a_inv
)
*
a_inv
;
}
template
<
bool
upper
>
static
inline
Real
g0
(
Real
z
)
{
return
std
::
exp
(
-
sign
(
upper
)
*
z
);
}
template
<
bool
upper
>
static
inline
Real
int_g0
(
const
thrust
::
pair
<
Real
,
Real
>&
bounds
,
Real
alpha
)
{
alpha
*=
-
sign
(
upper
);
return
int_exp
(
alpha
,
bounds
.
second
)
-
int_exp
(
alpha
,
bounds
.
first
);
}
template
<
bool
upper
>
static
inline
Real
int_g0_z
(
const
thrust
::
pair
<
Real
,
Real
>&
bounds
,
Real
alpha
)
{
alpha
*=
-
sign
(
upper
);
return
int_exp_z
(
alpha
,
bounds
.
second
)
-
int_exp_z
(
alpha
,
bounds
.
first
);
}
template
<
bool
upper
>
static
inline
Real
g1
(
Real
z
)
{
return
g0
<
upper
>
(
z
)
*
z
;
}
template
<
bool
upper
>
static
inline
Real
int_g1
(
const
thrust
::
pair
<
Real
,
Real
>&
bounds
,
Real
alpha
)
{
return
alpha
*
int_g0_z
<
upper
>
(
bounds
,
alpha
);
}
template
<
bool
upper
>
static
inline
Real
int_g1_z
(
const
thrust
::
pair
<
Real
,
Real
>&
bounds
,
Real
alpha
)
{
return
alpha
*
(
int_exp_zz
(
-
alpha
*
sign
(
upper
),
bounds
.
second
)
-
int_exp_zz
(
-
alpha
*
sign
(
upper
),
bounds
.
first
));
}
public
:
template
<
Int
yj_xi
,
UInt
dim_q
,
typename
OutType
,
typename
InType
,
typename
KelvinType
>
static
inline
void
integrate
(
OutType
&&
out
,
InType
&&
in
,
KelvinType
&&
kelvin
,
VectorProxy
<
const
Real
,
dim_q
>&
q_
,
Real
dij
,
Real
dl
)
{
const
auto
q
=
q_
.
l2norm
();
const
auto
alpha
=
q
*
dl
,
beta
=
q
*
dij
;
// Integrate on [-1, 0]
const
auto
bounds_low
=
thrust
::
make_pair
(
-
1
,
0
);
constexpr
bool
domain_low
=
yj_xi
>
0
;
auto
g0_mul
{
kelvin
.
template
applyU0
<
domain_low
>
(
q_
,
in
)};
auto
g1_mul
{
kelvin
.
template
applyU1
<
domain_low
>
(
q_
,
in
)};
decltype
(
g0_mul
)
res
;
res
=
g0_mul
*
(
g0
<
domain_low
>
(
beta
)
*
(
int_g0
<
domain_low
>
(
bounds_low
,
alpha
)
+
int_g0_z
<
domain_low
>
(
bounds_low
,
alpha
)));
res
+=
g1_mul
*
(
g1
<
domain_low
>
(
beta
)
*
(
int_g0
<
domain_low
>
(
bounds_low
,
alpha
)
+
int_g0_z
<
domain_low
>
(
bounds_low
,
alpha
))
+
g0
<
domain_low
>
(
beta
)
*
(
int_g1
<
domain_low
>
(
bounds_low
,
alpha
)
+
int_g1_z
<
domain_low
>
(
bounds_low
,
alpha
)));
// Integrate on [0, 1]
const
auto
bounds_up
=
thrust
::
make_pair
(
0
,
1
);
constexpr
bool
domain_up
=
yj_xi
>=
0
;
// Avoid computing
if
(
yj_xi
==
0
)
{
g0_mul
=
kelvin
.
template
applyU0
<
true
>
(
q_
,
in
);
g1_mul
=
kelvin
.
template
applyU1
<
true
>
(
q_
,
in
);
}
res
+=
g0_mul
*
(
g0
<
domain_up
>
(
beta
)
*
(
int_g0
<
domain_up
>
(
bounds_up
,
alpha
)
-
int_g0_z
<
domain_up
>
(
bounds_up
,
alpha
)));
res
+=
g1_mul
*
(
g1
<
domain_up
>
(
beta
)
*
(
int_g0
<
domain_up
>
(
bounds_up
,
alpha
)
-
int_g0_z
<
domain_up
>
(
bounds_up
,
alpha
))
+
g0
<
domain_up
>
(
beta
)
*
(
int_g1
<
domain_up
>
(
bounds_up
,
alpha
)
-
int_g1_z
<
domain_up
>
(
bounds_up
,
alpha
)));
// res *= dl;
out
+=
res
;
}
};
}
// namespace influence
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
Event Timeline
Log In to Comment