Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F68717287
kelvin_helper.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, Jun 28, 15:36
Size
9 KB
Mime Type
text/x-c++
Expires
Sun, Jun 30, 15:36 (2 d)
Engine
blob
Format
Raw Data
Handle
18586947
Attached To
rTAMAAS tamaas
kelvin_helper.hh
View Options
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-19 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 KELVIN_HELPER_HH
#define KELVIN_HELPER_HH
/* -------------------------------------------------------------------------- */
#include "grid.hh"
#include "grid_hermitian.hh"
#include "influence.hh"
#include "integration/accumulator.hh"
#include "model.hh"
#include "model_type.hh"
#include <tuple>
/* -------------------------------------------------------------------------- */
namespace
tamaas
{
namespace
detail
{
/* -------------------------------------------------------------------------- */
/// Trait for kelvin local types
template
<
typename
T
>
struct
KelvinTrait
;
template
<
UInt
dim
>
struct
KelvinTrait
<
influence
::
Kelvin
<
dim
,
0
>>
{
using
source_t
=
VectorProxy
<
Complex
,
dim
>
;
using
out_t
=
VectorProxy
<
Complex
,
dim
>
;
};
template
<
UInt
dim
>
struct
KelvinTrait
<
influence
::
Kelvin
<
dim
,
1
>>
{
using
source_t
=
SymMatrixProxy
<
Complex
,
dim
>
;
using
out_t
=
VectorProxy
<
Complex
,
dim
>
;
};
template
<
UInt
dim
>
struct
KelvinTrait
<
influence
::
Kelvin
<
dim
,
2
>>
{
using
source_t
=
SymMatrixProxy
<
Complex
,
dim
>
;
using
out_t
=
SymMatrixProxy
<
Complex
,
dim
>
;
};
template
<
UInt
dim
>
struct
KelvinTrait
<
influence
::
Boussinesq
<
dim
,
0
>>
{
using
source_t
=
VectorProxy
<
Complex
,
dim
>
;
using
out_t
=
VectorProxy
<
Complex
,
dim
>
;
};
template
<
UInt
dim
>
struct
KelvinTrait
<
influence
::
Boussinesq
<
dim
,
1
>>
{
using
source_t
=
VectorProxy
<
Complex
,
dim
>
;
using
out_t
=
SymMatrixProxy
<
Complex
,
dim
>
;
};
/* -------------------------------------------------------------------------- */
/// Helper to apply integral representation on output
template
<
model_type
type
,
typename
kelvin_t
,
typename
=
typename
KelvinTrait
<
kelvin_t
>::
source_t
>
struct
KelvinHelper
{
using
trait
=
model_type_traits
<
type
>
;
static
constexpr
UInt
dim
=
trait
::
dimension
;
static
constexpr
UInt
bdim
=
trait
::
boundary_dimension
;
using
BufferType
=
GridHermitian
<
Real
,
bdim
>
;
using
source_t
=
typename
KelvinTrait
<
kelvin_t
>::
source_t
;
using
out_t
=
typename
KelvinTrait
<
kelvin_t
>::
out_t
;
virtual
~
KelvinHelper
()
=
default
;
/**
@brief Apply the regular part of Kelvin to source and sum into output
This function performs the linear integration algorithm using the
accumulator.
*/
void
applyIntegral
(
std
::
vector
<
BufferType
>&
source
,
std
::
vector
<
BufferType
>&
out
,
const
Grid
<
Real
,
bdim
>&
wavevectors
,
Real
domain_size
,
const
kelvin_t
&
kelvin
)
{
// Sanity check
if
(
source
.
size
()
!=
out
.
size
())
TAMAAS_EXCEPTION
(
"Linear integration requires source and out of same sizes"
);
accumulator
.
makeUniformMesh
(
source
.
size
(),
domain_size
);
auto
waverange
=
range
<
VectorProxy
<
const
Real
,
bdim
>>
(
wavevectors
);
for
(
auto
&&
tuple
:
accumulator
.
forward
(
source
,
wavevectors
))
{
auto
&&
l
=
std
::
get
<
0
>
(
tuple
);
auto
&&
xl_
=
std
::
get
<
1
>
(
tuple
);
auto
&&
acc_g0
=
std
::
get
<
2
>
(
tuple
);
auto
&&
acc_g1
=
std
::
get
<
3
>
(
tuple
);
Loop
::
loop
(
[
xl
=
xl_
,
&
kelvin
](
auto
qv
,
auto
u
,
auto
acc_g0
,
auto
acc_g1
)
{
Real
q
=
qv
.
l2norm
();
auto
dense_G0
=
dense
(
acc_g0
),
dense_G1
=
dense
(
acc_g1
);
auto
tmp
=
kelvin
.
template
applyU1
<
false
,
true
>
(
qv
,
dense_G0
);
tmp
*=
-
q
*
xl
;
tmp
+=
kelvin
.
template
applyU0
<
false
,
true
>
(
qv
,
dense_G0
);
tmp
+=
kelvin
.
template
applyU1
<
false
,
true
>
(
qv
,
dense_G1
);
tmp
*=
std
::
exp
(
-
q
*
xl
);
u
+=
tmp
;
},
waverange
.
headless
(),
range
<
out_t
>
(
out
[
l
]).
headless
(),
range
<
source_t
>
(
acc_g0
).
headless
(),
range
<
source_t
>
(
acc_g1
).
headless
());
}
for
(
auto
&&
tuple
:
accumulator
.
backward
(
source
,
wavevectors
))
{
auto
&&
l
=
std
::
get
<
0
>
(
tuple
);
auto
&&
xl_
=
std
::
get
<
1
>
(
tuple
);
auto
&&
acc_g0
=
std
::
get
<
2
>
(
tuple
);
auto
&&
acc_g1
=
std
::
get
<
3
>
(
tuple
);
Loop
::
loop
(
[
xl
=
xl_
,
&
kelvin
](
auto
qv
,
auto
u
,
auto
acc_g0
,
auto
acc_g1
)
{
Real
q
=
qv
.
l2norm
();
auto
dense_G0
=
dense
(
acc_g0
),
dense_G1
=
dense
(
acc_g1
);
auto
tmp
=
kelvin
.
template
applyU1
<
true
,
true
>
(
qv
,
dense_G0
);
tmp
*=
-
q
*
xl
;
tmp
+=
kelvin
.
template
applyU0
<
true
,
true
>
(
qv
,
dense_G0
);
tmp
+=
kelvin
.
template
applyU1
<
true
,
true
>
(
qv
,
dense_G1
);
tmp
*=
std
::
exp
(
q
*
xl
);
u
+=
tmp
;
// should automatically symmetrize if needed
},
waverange
.
headless
(),
range
<
out_t
>
(
out
[
l
]).
headless
(),
range
<
source_t
>
(
acc_g0
).
headless
(),
range
<
source_t
>
(
acc_g1
).
headless
());
}
}
/**
@brief Apply the regular part of Kelvin to source and sum into output
This function performs the cutoff integration algorithm. Not to be confused
with its overload above.
*/
void
applyIntegral
(
std
::
vector
<
BufferType
>&
/*source*/
,
BufferType
&
/*out*/
,
UInt
/*layer*/
,
const
Grid
<
Real
,
bdim
>&
/*wavevectors*/
,
Real
/*domain_size*/
,
const
kelvin_t
&
/*kelvin*/
)
{
TAMAAS_EXCEPTION
(
"Not implemented"
);
}
/// Apply free term of distribution derivative
void
applyFreeTerm
(
std
::
vector
<
BufferType
>&
source
,
std
::
vector
<
BufferType
>&
out
,
const
kelvin_t
&
kelvin
)
{
for
(
UInt
l
:
Loop
::
range
(
source
.
size
()))
{
applyFreeTerm
(
source
[
l
],
out
[
l
],
kelvin
);
}
}
/// Apply free term of distribution derivative (layer)
void
applyFreeTerm
(
BufferType
&
/*source*/
,
BufferType
&
/*out*/
,
const
kelvin_t
&
/*kelvin*/
)
{}
/// Making the output free of communist NaN
void
makeFundamentalGreatAgain
(
std
::
vector
<
BufferType
>&
out
)
{
for
(
auto
&&
layer
:
out
)
makeFundamentalGreatAgain
(
layer
);
}
/// Making the output free of communist NaN (layer)
void
makeFundamentalGreatAgain
(
BufferType
&
/*out_layer*/
)
{}
protected
:
Accumulator
<
type
,
source_t
>
accumulator
;
};
/// Applying free term for double gradient of Kelvin
template
<>
inline
void
KelvinHelper
<
model_type
::
volume_2d
,
influence
::
Kelvin
<
3
,
2
>>::
applyFreeTerm
(
BufferType
&
source
,
BufferType
&
out
,
const
influence
::
Kelvin
<
3
,
2
>&
kelvin
)
{
// No need for headless loops as discontinuity term does not depend on q
Loop
::
loop
(
[
&
kelvin
](
auto
u
,
auto
w
)
{
u
+=
kelvin
.
applyDiscontinuityTerm
(
dense
(
w
));
},
range
<
out_t
>
(
out
),
range
<
source_t
>
(
source
));
}
/// Setting fundamental mode to zero: we've got no other choice here
template
<>
inline
void
KelvinHelper
<
model_type
::
volume_2d
,
influence
::
Kelvin
<
3
,
0
>>::
makeFundamentalGreatAgain
(
BufferType
&
layer
)
{
out_t
u
(
&
layer
(
0
));
u
=
0
;
}
/* -------------------------------------------------------------------------- */
template
<
model_type
type
>
struct
SurfaceTractionHelper
{
using
trait
=
model_type_traits
<
type
>
;
static
constexpr
UInt
dim
=
trait
::
dimension
;
static
constexpr
UInt
bdim
=
trait
::
boundary_dimension
;
using
BufferType
=
GridHermitian
<
Real
,
bdim
>
;
using
kelvin_t
=
influence
::
Kelvin
<
3
,
2
>
;
using
source_t
=
typename
KelvinTrait
<
kelvin_t
>::
source_t
;
using
out_t
=
typename
KelvinTrait
<
kelvin_t
>::
out_t
;
/// Compute surface tractions due to eigenstress distribution
template
<
bool
apply_q_power
>
void
computeSurfaceTractions
(
std
::
vector
<
BufferType
>&
source
,
BufferType
&
tractions
,
const
Grid
<
Real
,
bdim
>&
wavevectors
,
Real
domain_size
,
const
kelvin_t
&
kelvin
,
const
influence
::
ElasticHelper
<
dim
>&
el
)
{
accumulator
.
makeUniformMesh
(
source
.
size
(),
domain_size
);
// only doing backward loop because forward loop gives 0 for l = 0
for
(
auto
&&
tuple
:
accumulator
.
backward
(
source
,
wavevectors
))
{
auto
&&
l
=
std
::
get
<
0
>
(
tuple
);
auto
&&
acc_g0
=
std
::
get
<
2
>
(
tuple
);
auto
&&
acc_g1
=
std
::
get
<
3
>
(
tuple
);
// Just skip all layers and let the accumulator accumulate
if
(
l
!=
0
)
continue
;
Loop
::
loop
(
[
&
kelvin
,
&
el
](
auto
qv
,
auto
t
,
auto
acc_g0
,
auto
acc_g1
)
{
// computing strains
auto
dense_G0
=
dense
(
acc_g0
),
dense_G1
=
dense
(
acc_g1
);
auto
strains
=
kelvin
.
template
applyU0
<
true
,
apply_q_power
>
(
qv
,
dense_G0
);
strains
+=
kelvin
.
template
applyU1
<
true
,
apply_q_power
>
(
qv
,
dense_G1
);
// computing traction vector
Vector
<
Real
,
dim
>
normal
{{{
0
,
0
,
-
1
}}};
t
=
el
(
strains
)
*
normal
;
},
range
<
VectorProxy
<
const
Real
,
bdim
>>
(
wavevectors
).
headless
(),
range
<
VectorProxy
<
Complex
,
dim
>>
(
tractions
).
headless
(),
range
<
source_t
>
(
acc_g0
).
headless
(),
range
<
source_t
>
(
acc_g1
).
headless
());
}
}
protected
:
Accumulator
<
type
,
source_t
>
accumulator
;
};
}
// namespace detail
}
// namespace tamaas
#endif
// KELVIN_HELPER_HH
Event Timeline
Log In to Comment