Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F68758789
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, Jun 28, 19:43
Size
9 KB
Mime Type
text/x-c++
Expires
Sun, Jun 30, 19:43 (2 d)
Engine
blob
Format
Raw Data
Handle
18629389
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
{
/// Mother class for influence functions
template
<
UInt
dim
>
class
Influence
{
public
:
Influence
(
StaticVector
<
const
Real
,
dim
>
domain
)
:
domain
(
domain
)
{}
protected
:
StaticVector
<
const
Real
,
dim
>
domain
;
};
/// Class representing the basic influence functions from Stanley & Kato (1997)
template
<
UInt
dim
>
class
Basic
:
public
Influence
<
dim
>
{
public
:
Basic
(
Real
E_hertz
,
StaticVector
<
const
Real
,
dim
>
domain
)
:
Influence
<
dim
>
(
domain
),
E_hertz
(
E_hertz
)
{}
/// Influence function: computes the influence coefficient F
CUDA_LAMBDA
void
operator
()(
StaticVector
<
Real
,
dim
>&&
q_vec
,
Complex
&
F
)
const
{
q_vec
*=
2
*
M_PI
;
q_vec
/=
this
->
domain
;
F
=
2.
/
(
q_vec
.
l2norm
()
*
E_hertz
);
}
private
:
Real
E_hertz
;
};
template
<
UInt
dim
>
class
Surface
:
public
Influence
<
dim
>
{};
/// Class for the influence with all components of applied tractions (dim = 2)
template
<>
class
Surface
<
2
>
:
public
Influence
<
2
>
{
constexpr
static
UInt
dim
=
2
;
public
:
Surface
(
Real
E
,
Real
nu
,
StaticVector
<
const
Real
,
dim
>
domain
)
:
Influence
<
dim
>
(
domain
),
E
(
E
),
nu
(
nu
)
{}
CUDA_LAMBDA
void
operator
()(
StaticVector
<
Real
,
dim
>&&
q_
,
StaticMatrix
<
Complex
,
dim
+
1
,
dim
+
1
>&&
F
)
const
{
q_
*=
2
*
M_PI
;
q_
/=
this
->
domain
;
const
Real
q
=
q_
.
l2norm
();
// First line of influence matrix
F
(
0
,
0
)
=
1.
/
(
q
*
q
)
*
(
-
nu
*
nu
*
q_
(
0
)
*
q_
(
0
)
-
nu
*
q_
(
1
)
*
q_
(
1
)
+
q
*
q
);
F
(
0
,
1
)
=
-
q_
(
0
)
*
q_
(
1
)
/
(
q
*
q
)
*
(
1
+
nu
)
*
2
*
nu
;
F
(
0
,
2
)
=
I
*
q_
(
0
)
/
q
*
(
1
+
nu
)
*
(
1
-
2
*
nu
);
// Second line of influence matrix
F
(
1
,
0
)
=
-
q_
(
0
)
*
q_
(
1
)
/
(
q
*
q
)
*
(
1
+
nu
)
*
2
*
nu
;
F
(
1
,
1
)
=
1.
/
(
q
*
q
)
*
(
-
nu
*
nu
*
q_
(
1
)
*
q_
(
1
)
-
nu
*
q_
(
0
)
*
q_
(
0
)
+
q
*
q
);
F
(
1
,
2
)
=
I
*
q_
(
1
)
/
q
*
(
1
+
nu
)
*
(
1
-
2
*
nu
);
// Third line of influence matrix
F
(
2
,
0
)
=
-
F
(
0
,
2
);
F
(
2
,
1
)
=
-
F
(
1
,
2
);
F
(
2
,
2
)
=
2
*
(
1
-
nu
*
nu
);
F
*=
1.
/
(
E
*
q
);
}
private
:
Real
E
,
nu
;
const
Complex
I
=
Complex
(
0
,
1
);
};
/// Class for the influence with all components of applied tractions (dim = 2)
template
<>
class
Surface
<
1
>
:
public
Influence
<
1
>
{
constexpr
static
UInt
dim
=
1
;
public
:
Surface
(
Real
E
,
Real
nu
,
StaticVector
<
const
Real
,
dim
>
domain
)
:
Influence
<
dim
>
(
domain
),
E
(
E
),
nu
(
nu
)
{}
CUDA_LAMBDA
void
operator
()(
StaticVector
<
Real
,
dim
>&&
q_
,
StaticMatrix
<
Complex
,
dim
+
1
,
dim
+
1
>&&
F
)
const
{
q_
*=
2
*
M_PI
;
q_
/=
this
->
domain
;
const
Real
q
=
q_
.
l2norm
();
// First line of influence matrix
F
(
0
,
0
)
=
2.
/
q
*
(
1
-
nu
*
nu
);
F
(
0
,
1
)
=
I
/
q_
(
0
)
*
(
1
+
nu
)
*
(
1
-
2
*
nu
);
// Second line of influence matrix
F
(
1
,
0
)
=
-
F
(
0
,
1
);
F
(
1
,
1
)
=
F
(
0
,
0
);
F
*=
1.
/
E
;
}
private
:
Real
E
,
nu
;
const
Complex
I
=
Complex
(
0
,
1
);
};
/// Class for the Kelvin tensor
template
<
UInt
dim
,
UInt
derivative_order
>
class
Kelvin
;
template
<
UInt
derivative_order
>
class
Kelvin
<
3
,
derivative_order
>
:
Influence
<
3
>
{
constexpr
static
UInt
dim
=
3
;
public
:
Kelvin
(
Real
mu
,
Real
nu
,
StaticVector
<
const
Real
,
dim
>
domain
)
:
Influence
<
dim
>
(
domain
),
mu
(
mu
),
nu
(
nu
)
{}
private
:
/// Compute leading linear term A
template
<
bool
upper
>
inline
void
computeA
(
StaticVector
<
const
Real
,
dim
-
1
>&
q_
,
StaticMatrix
<
Complex
,
dim
,
dim
>&
A
)
const
;
/// Compute leading constant term B, which does not depend on y_3
inline
void
computeB
(
StaticVector
<
const
Real
,
dim
-
1
>&
q_
,
StaticMatrix
<
Complex
,
dim
,
dim
>&
B
)
const
;
public
:
/// Compute leading terms
template
<
bool
upper
>
inline
void
computeTerms
(
StaticVector
<
const
Real
,
dim
-
1
>&
q_
,
StaticMatrix
<
Complex
,
dim
,
dim
>&
A
,
StaticMatrix
<
Complex
,
dim
,
dim
>&
B
)
const
{
computeA
<
upper
>
(
q_
,
A
);
computeB
(
q_
,
B
);
}
private
:
Real
mu
,
nu
;
const
Complex
I
=
Complex
(
0
,
1
);
};
/* -------------------------------------------------------------------------- */
/* Base Kelvin tensor */
/* -------------------------------------------------------------------------- */
template
<>
inline
void
Kelvin
<
3
,
0
>::
computeB
(
StaticVector
<
const
Real
,
dim
-
1
>&
q_
,
StaticMatrix
<
Complex
,
dim
,
dim
>&
B
)
const
{
const
Real
q
=
q_
.
l2norm
();
const
Real
bq1
=
q_
(
0
)
/
q
;
const
Real
bq2
=
q_
(
1
)
/
q
;
for
(
UInt
i
=
0
;
i
<
dim
;
++
i
)
B
(
i
,
i
)
=
4
*
(
1
-
nu
);
B
(
2
,
2
)
-=
1
;
B
(
0
,
0
)
-=
bq1
*
bq1
;
B
(
1
,
1
)
-=
bq2
*
bq2
;
B
(
0
,
1
)
-=
bq1
*
bq2
;
B
(
1
,
0
)
-=
bq1
*
bq2
;
B
(
0
,
2
)
=
B
(
1
,
2
)
=
B
(
2
,
0
)
=
B
(
2
,
1
)
=
0
;
B
*=
(
1
/
(
8
*
mu
*
q
*
(
1
-
nu
)));
}
template
<>
template
<
bool
upper
>
inline
void
Kelvin
<
3
,
0
>::
computeA
(
StaticVector
<
const
Real
,
dim
-
1
>&
q_
,
StaticMatrix
<
Complex
,
dim
,
dim
>&
A
)
const
{
const
Real
q
=
q_
.
l2norm
();
const
Real
bq1
=
q_
(
0
)
/
q
;
const
Real
bq2
=
q_
(
1
)
/
q
;
constexpr
Real
sign
=
(
upper
)
?
1.
:
-
1
;
A
(
0
,
0
)
=
sign
*
bq1
*
bq1
;
A
(
1
,
1
)
=
sign
*
bq2
*
bq2
;
A
(
2
,
2
)
=
1
;
A
(
1
,
0
)
=
A
(
0
,
1
)
=
sign
*
bq1
*
bq2
;
A
(
2
,
0
)
=
A
(
0
,
2
)
=
-
I
*
bq1
;
A
(
2
,
1
)
=
A
(
1
,
2
)
=
-
I
*
bq2
;
A
*=
(
1
/
(
8
*
mu
*
(
1
-
nu
)));
}
/* -------------------------------------------------------------------------- */
class
KelvinIntegrator
{
static
inline
Complex
constantPrimitive
(
const
Real
&
x
,
const
Real
&
alpha
,
const
Complex
&
beta
,
const
Complex
&
gamma
)
{
return
std
::
exp
(
alpha
*
x
)
/
alpha
*
(
beta
*
x
+
gamma
-
beta
/
alpha
);
}
static
inline
Complex
linearPrimitive
(
const
Real
&
x
,
const
Real
&
alpha
,
const
Complex
&
beta
,
const
Complex
&
gamma
)
{
auto
coeff
=
gamma
-
2.
*
beta
/
alpha
;
return
std
::
exp
(
alpha
*
x
)
/
alpha
*
(
beta
*
x
*
x
+
coeff
*
x
-
coeff
/
alpha
);
}
static
inline
void
computePolyCoefficients
(
Complex
&
beta
,
Complex
&
gamma
,
const
Real
&
dl
,
const
Real
&
dij
,
const
Complex
&
A
,
const
Complex
&
B
)
{
beta
=
A
*
dl
;
gamma
=
A
*
dij
+
B
;
}
public
:
template
<
UInt
dim
>
static
inline
void
partialIntegration
(
bool
upper
,
StaticMatrix
<
Complex
,
dim
,
dim
>&
F
,
const
Real
&
alpha
,
const
Real
&
dij
,
const
Real
&
dl
,
const
StaticMatrix
<
Complex
,
dim
,
dim
>&
A
,
const
StaticMatrix
<
Complex
,
dim
,
dim
>&
B
)
{
Real
xsup
=
(
upper
)
?
1
:
0
;
Real
xinf
=
(
upper
)
?
0
:-
1
;
Real
multiplier
=
(
upper
)
?
-
1
:
1
;
Complex
beta
=
0
,
gamma
=
0
;
F
=
0
;
for
(
UInt
i
=
0
;
i
<
dim
;
++
i
)
{
for
(
UInt
j
=
0
;
j
<
dim
;
++
j
)
{
computePolyCoefficients
(
beta
,
gamma
,
dl
,
dij
,
A
(
i
,
j
),
B
(
i
,
j
));
F
(
i
,
j
)
+=
constantPrimitive
(
xsup
,
alpha
,
beta
,
gamma
)
-
constantPrimitive
(
xinf
,
alpha
,
beta
,
gamma
);
F
(
i
,
j
)
+=
multiplier
*
(
linearPrimitive
(
xsup
,
alpha
,
beta
,
gamma
)
-
linearPrimitive
(
xinf
,
alpha
,
beta
,
gamma
));
F
(
i
,
j
)
*=
std
::
exp
(
alpha
*
dij
/
dl
);
}
}
}
template
<
UInt
dim
,
UInt
derivative
,
Int
yj_xi
>
static
inline
void
integrate
(
const
Kelvin
<
dim
,
derivative
>&
kelvin
,
Real
yj
,
Real
xi
,
Real
dl
,
StaticVector
<
const
Real
,
dim
-
1
>&
q
,
StaticMatrix
<
Complex
,
dim
,
dim
>&
F
)
{
Complex
amem
[
dim
*
dim
];
Complex
bmem
[
dim
*
dim
];
StaticMatrix
<
Complex
,
dim
,
dim
>
A
{
*
amem
};
StaticMatrix
<
Complex
,
dim
,
dim
>
B
{
*
bmem
};
const
Real
q_norm
=
q
.
l2norm
();
const
Real
dij
=
yj
-
xi
;
constexpr
Real
sign
=
(
yj_xi
>
0
)
?
1.
:
-
1
;
Real
alpha
=
0
;
// Integrate [-1, 0]
kelvin
.
template
computeTerms
<
(
yj_xi
>
0
)
>
(
q
,
A
,
B
);
alpha
=
-
sign
*
q_norm
*
dl
;
partialIntegration
(
false
,
F
,
alpha
,
dij
,
dl
,
A
,
B
);
// Integrate [0, 1]
if
(
yj_xi
==
0
)
{
// avoid computing the same A & B twice
kelvin
.
template
computeTerms
<
true
>
(
q
,
A
,
B
);
alpha
*=
-
1
;
}
partialIntegration
(
true
,
F
,
alpha
,
dij
,
dl
,
A
,
B
);
}
};
}
// namespace influence
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
Event Timeline
Log In to Comment