Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65043296
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, May 31, 08:11
Size
15 KB
Mime Type
text/x-c++
Expires
Sun, Jun 2, 08:11 (2 d)
Engine
blob
Format
Raw Data
Handle
17994243
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
(
VectorProxy
<
const
Real
,
dim
>
domain
)
:
domain
(
domain
)
{}
protected
:
VectorProxy
<
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
,
VectorProxy
<
const
Real
,
dim
>
domain
)
:
Influence
<
dim
>
(
domain
),
E_hertz
(
E_hertz
)
{}
/// Influence function: computes the influence coefficient F
CUDA_LAMBDA
void
operator
()(
VectorProxy
<
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
,
VectorProxy
<
const
Real
,
dim
>
domain
)
:
Influence
<
dim
>
(
domain
),
E
(
E
),
nu
(
nu
)
{}
CUDA_LAMBDA
void
operator
()(
VectorProxy
<
Real
,
dim
>&&
q_
,
MatrixProxy
<
Complex
,
dim
+
1
,
dim
+
1
>&&
F
)
const
{
q_
*=
2
*
M_PI
;
q_
/=
this
->
domain
;
const
Real
q
=
q_
.
l2norm
();
q_
/=
q
;
// Diagonal
F
(
0
,
0
)
=
2
*
(
1
+
nu
)
*
(
1
-
nu
*
q_
(
0
)
*
q_
(
0
));
F
(
1
,
1
)
=
2
*
(
1
+
nu
)
*
(
1
-
nu
*
q_
(
1
)
*
q_
(
1
));
F
(
2
,
2
)
=
2
*
(
1
-
nu
*
nu
);
F
(
0
,
1
)
=
F
(
1
,
0
)
=
-
q_
(
0
)
*
q_
(
1
)
*
2
*
nu
*
(
1
+
nu
);
F
(
0
,
2
)
=
I
*
q_
(
0
)
*
(
1
+
nu
)
*
(
1
-
2
*
nu
);
F
(
1
,
2
)
=
I
*
q_
(
1
)
*
(
1
+
nu
)
*
(
1
-
2
*
nu
);
F
(
2
,
0
)
=
-
F
(
0
,
2
);
F
(
2
,
1
)
=
-
F
(
1
,
2
);
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
,
VectorProxy
<
const
Real
,
dim
>
domain
)
:
Influence
<
dim
>
(
domain
),
E
(
E
),
nu
(
nu
)
{}
CUDA_LAMBDA
void
operator
()(
VectorProxy
<
Real
,
dim
>&&
q_
,
MatrixProxy
<
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
,
VectorProxy
<
const
Real
,
dim
>
domain
)
:
Influence
<
dim
>
(
domain
),
mu
(
mu
),
nu
(
nu
)
{}
/// Compute leading linear term A and leading constant B
template
<
bool
upper
>
CUDA_LAMBDA
inline
void
computeTerms
(
VectorProxy
<
Real
,
dim
-
1
>&&
q_
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
A
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
B
)
const
;
/// Compute linear integral version of operator
template
<
Int
yj_xi
>
CUDA_LAMBDA
inline
void
linearIntegrate
(
Real
yj
,
Real
xi
,
Real
dl
,
VectorProxy
<
Real
,
dim
-
1
>&&
q
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
F
)
const
;
private
:
Real
mu
,
nu
;
const
Complex
I
=
Complex
(
0
,
1
);
};
/* -------------------------------------------------------------------------- */
/// Kelvin tensor for upper domain
template
<>
template
<>
inline
void
Kelvin
<
3
,
0
>::
computeTerms
<
true
>
(
VectorProxy
<
Real
,
dim
-
1
>&&
q_
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
A
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
B
)
const
{
q_
*=
2
*
M_PI
;
q_
/=
VectorProxy
<
const
Real
,
dim
-
1
>
(
this
->
domain
(
0
));
const
Real
q
=
q_
.
l2norm
();
const
Real
bq1
=
q_
(
0
)
/
q
;
const
Real
bq2
=
q_
(
1
)
/
q
;
A
(
0
,
0
)
=
-
bq1
*
bq1
;
B
(
0
,
0
)
=
3
*
bq1
*
bq1
+
4
*
bq2
*
bq2
-
4
*
nu
;
A
(
0
,
1
)
=
A
(
1
,
0
)
=
-
bq1
*
bq2
;
B
(
0
,
1
)
=
B
(
1
,
0
)
=
-
bq1
*
bq2
;
A
(
0
,
2
)
=
A
(
2
,
0
)
=
-
I
*
bq1
;
B
(
0
,
2
)
=
B
(
2
,
0
)
=
0
;
A
(
1
,
1
)
=
-
bq2
*
bq2
;
B
(
1
,
1
)
=
4
*
bq1
*
bq1
+
3
*
bq2
*
bq2
-
4
*
nu
;
A
(
1
,
2
)
=
A
(
2
,
1
)
=
-
I
*
bq2
;
B
(
1
,
2
)
=
B
(
2
,
1
)
=
0
;
A
(
2
,
2
)
=
1
;
B
(
2
,
2
)
=
3
-
4
*
nu
;
A
*=
(
1
/
(
8
*
mu
*
(
1
-
nu
)));
B
*=
(
1
/
(
8
*
mu
*
q
*
(
1
-
nu
)));
}
/// Kelvin tensor for lower domain
template
<>
template
<>
inline
void
Kelvin
<
3
,
0
>::
computeTerms
<
false
>
(
VectorProxy
<
Real
,
dim
-
1
>&&
q_
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
A
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
B
)
const
{
q_
*=
2
*
M_PI
;
q_
/=
VectorProxy
<
const
Real
,
dim
-
1
>
(
this
->
domain
(
0
));
const
Real
q
=
q_
.
l2norm
();
const
Real
bq1
=
q_
(
0
)
/
q
;
const
Real
bq2
=
q_
(
1
)
/
q
;
A
(
0
,
0
)
=
bq1
*
bq1
;
B
(
0
,
0
)
=
3
*
bq1
*
bq1
+
4
*
bq2
*
bq2
-
4
*
nu
;
A
(
0
,
1
)
=
A
(
1
,
0
)
=
bq1
*
bq2
;
B
(
0
,
1
)
=
B
(
1
,
0
)
=
-
bq1
*
bq2
;
A
(
0
,
2
)
=
A
(
2
,
0
)
=
-
I
*
bq1
;
B
(
0
,
2
)
=
B
(
2
,
0
)
=
0
;
A
(
1
,
1
)
=
bq2
*
bq2
;
B
(
1
,
1
)
=
4
*
bq1
*
bq1
+
3
*
bq2
*
bq2
-
4
*
nu
;
A
(
1
,
2
)
=
A
(
2
,
1
)
=
-
I
*
bq2
;
B
(
1
,
2
)
=
B
(
2
,
1
)
=
0
;
A
(
2
,
2
)
=
-
1
;
B
(
2
,
2
)
=
3
-
4
*
nu
;
A
*=
(
1
/
(
8
*
mu
*
(
1
-
nu
)));
B
*=
(
1
/
(
8
*
mu
*
q
*
(
1
-
nu
)));
}
/* -------------------------------------------------------------------------- */
/// First derivative of Kelvin tensor (upper domain)
template
<>
template
<>
inline
void
Kelvin
<
3
,
1
>::
computeTerms
<
true
>
(
VectorProxy
<
Real
,
dim
-
1
>&&
q_
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
A
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
B
)
const
{
q_
*=
2
*
M_PI
;
q_
/=
VectorProxy
<
const
Real
,
dim
-
1
>
(
this
->
domain
(
0
));
const
Real
q
=
q_
.
l2norm
();
const
Real
bq1
=
q_
(
0
)
/
q
;
const
Real
bq2
=
q_
(
1
)
/
q
;
A
(
0
,
0
)
=
bq1
*
bq1
;
B
(
0
,
0
)
=
-
4
*
(
1
-
nu
);
A
(
0
,
1
)
=
A
(
1
,
0
)
=
bq1
*
bq2
;
B
(
0
,
1
)
=
B
(
1
,
0
)
=
0
;
A
(
0
,
2
)
=
A
(
2
,
0
)
=
I
*
bq1
;
B
(
0
,
2
)
=
B
(
2
,
0
)
=
-
I
*
bq1
;
A
(
1
,
1
)
=
bq2
*
bq2
;
B
(
1
,
1
)
=
-
4
*
(
1
-
nu
);
A
(
1
,
2
)
=
A
(
2
,
1
)
=
I
*
bq2
;
B
(
1
,
2
)
=
B
(
2
,
1
)
=
-
I
*
bq2
;
A
(
2
,
2
)
=
-
1
;
B
(
2
,
2
)
=
4
*
nu
-
2
;
A
*=
(
q
/
(
8
*
mu
*
(
1
-
nu
)));
B
*=
(
1
/
(
8
*
mu
*
(
1
-
nu
)));
}
/// First derivative of Kelvin tensor (lower domain)
template
<>
template
<>
inline
void
Kelvin
<
3
,
1
>::
computeTerms
<
false
>
(
VectorProxy
<
Real
,
dim
-
1
>&&
q_
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
A
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
B
)
const
{
q_
*=
2
*
M_PI
;
q_
/=
VectorProxy
<
const
Real
,
dim
-
1
>
(
this
->
domain
(
0
));
const
Real
q
=
q_
.
l2norm
();
const
Real
bq1
=
q_
(
0
)
/
q
;
const
Real
bq2
=
q_
(
1
)
/
q
;
A
(
0
,
0
)
=
bq1
*
bq1
;
B
(
0
,
0
)
=
4
*
(
1
-
nu
);
A
(
0
,
1
)
=
A
(
1
,
0
)
=
bq1
*
bq2
;
B
(
0
,
1
)
=
B
(
1
,
0
)
=
0
;
A
(
0
,
2
)
=
A
(
2
,
0
)
=
-
I
*
bq1
;
B
(
0
,
2
)
=
B
(
2
,
0
)
=
-
I
*
bq1
;
A
(
1
,
1
)
=
bq2
*
bq2
;
B
(
1
,
1
)
=
4
*
(
1
-
nu
);
A
(
1
,
2
)
=
A
(
2
,
1
)
=
-
I
*
bq2
;
B
(
1
,
2
)
=
B
(
2
,
1
)
=
-
I
*
bq2
;
A
(
2
,
2
)
=
-
1
;
B
(
2
,
2
)
=
-
(
4
*
nu
-
2
);
A
*=
(
q
/
(
8
*
mu
*
(
1
-
nu
)));
B
*=
(
1
/
(
8
*
mu
*
(
1
-
nu
)));
}
/* -------------------------------------------------------------------------- */
/// Regular second derivative of Kelvin tensor (upper domain)
template
<>
template
<>
inline
void
Kelvin
<
3
,
2
>::
computeTerms
<
true
>
(
VectorProxy
<
Real
,
dim
-
1
>&&
q_
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
A
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
B
)
const
{
q_
*=
2
*
M_PI
;
q_
/=
VectorProxy
<
const
Real
,
dim
-
1
>
(
this
->
domain
(
0
));
const
Real
q
=
q_
.
l2norm
();
const
Real
bq1
=
q_
(
0
)
/
q
;
const
Real
bq2
=
q_
(
1
)
/
q
;
A
(
0
,
0
)
=
-
bq1
*
bq1
;
B
(
0
,
0
)
=
bq1
*
bq1
-
4.
*
nu
+
4.
;
A
(
0
,
1
)
=
A
(
1
,
0
)
=
-
bq1
*
bq2
;
B
(
0
,
1
)
=
B
(
1
,
0
)
=
bq1
*
bq2
;
A
(
0
,
2
)
=
A
(
2
,
0
)
=
-
I
*
bq1
;
B
(
0
,
2
)
=
B
(
2
,
0
)
=
2.
*
I
*
bq1
;
A
(
1
,
1
)
=
-
bq2
*
bq2
;
B
(
1
,
1
)
=
bq2
*
bq2
-
4.
*
nu
+
4.
;
A
(
1
,
2
)
=
A
(
2
,
1
)
=
-
I
*
bq2
;
B
(
1
,
2
)
=
B
(
2
,
1
)
=
2.
*
I
*
bq2
;
A
(
2
,
2
)
=
1.
;
B
(
2
,
2
)
=
1.
-
4.
*
nu
;
A
*=
(
q
*
q
/
(
8.
*
mu
*
(
1.
-
nu
)));
B
*=
(
q
/
(
8.
*
mu
*
(
1.
-
nu
)));
}
/// Regular second derivative of Kelvin tensor (lower domain)
template
<>
template
<>
inline
void
Kelvin
<
3
,
2
>::
computeTerms
<
false
>
(
VectorProxy
<
Real
,
dim
-
1
>&&
q_
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
A
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
B
)
const
{
q_
*=
2
*
M_PI
;
q_
/=
VectorProxy
<
const
Real
,
dim
-
1
>
(
this
->
domain
(
0
));
const
Real
q
=
q_
.
l2norm
();
const
Real
bq1
=
q_
(
0
)
/
q
;
const
Real
bq2
=
q_
(
1
)
/
q
;
A
(
0
,
0
)
=
bq1
*
bq1
;
B
(
0
,
0
)
=
bq1
*
bq1
-
4.
*
nu
+
4.
;
A
(
0
,
1
)
=
A
(
1
,
0
)
=
bq1
*
bq2
;
B
(
0
,
1
)
=
B
(
1
,
0
)
=
bq1
*
bq2
;
A
(
0
,
2
)
=
A
(
2
,
0
)
=
-
I
*
bq1
;
B
(
0
,
2
)
=
B
(
2
,
0
)
=
-
2.
*
I
*
bq1
;
A
(
1
,
1
)
=
bq2
*
bq2
;
B
(
1
,
1
)
=
bq2
*
bq2
-
4.
*
nu
+
4.
;
A
(
1
,
2
)
=
A
(
2
,
1
)
=
-
I
*
bq2
;
B
(
1
,
2
)
=
B
(
2
,
1
)
=
-
2.
*
I
*
bq2
;
A
(
2
,
2
)
=
-
1.
;
B
(
2
,
2
)
=
1.
-
4.
*
nu
;
A
*=
(
q
*
q
/
(
8.
*
mu
*
(
1.
-
nu
)));
B
*=
(
q
/
(
8.
*
mu
*
(
1.
-
nu
)));
}
/* -------------------------------------------------------------------------- */
template
<
UInt
dim
,
UInt
derivative
,
Int
yi_xi
>
struct
KelvinIntegrator
;
inline
Real
kelvin_primitive
(
Real
x
,
Real
alpha
,
Real
beta
,
Real
gamma
,
Real
lambda
)
{
return
std
::
exp
(
alpha
*
x
)
/
alpha
*
(
beta
*
x
*
x
+
(
gamma
-
2
*
beta
/
alpha
)
*
x
+
lambda
+
2
*
beta
/
(
alpha
*
alpha
));
}
/// Integrate for yj > xi
template
<
UInt
dim
,
UInt
derivative
>
struct
KelvinIntegrator
<
dim
,
derivative
,
1
>
{
CUDA_LAMBDA
inline
void
integrate
(
const
Kelvin
<
dim
,
derivative
>&
kelvin
,
Real
yj
,
Real
xi
,
Real
dl
,
VectorProxy
<
Real
,
dim
-
1
>&&
q
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
F
)
{
Complex
amem
[
dim
*
dim
],
bmem
[
dim
*
dim
];
MatrixProxy
<
Complex
,
dim
,
dim
>
A
(
*
amem
),
B
(
*
bmem
);
const
Real
dij
=
yj
-
xi
;
kelvin
.
computeTerms
<
true
>
(
q
,
A
,
B
);
// Integrate [-1, 0] (upper domain)
Real
alpha
=
-
q
*
dl
,
beta
=
0
,
gamma
=
0
,
lambda
=
0
;
for
(
UInt
i
=
0
;
i
<
dim
;
++
i
)
{
for
(
UInt
j
=
0
;
j
<
dim
;
++
j
)
{
beta
=
A
(
i
,
j
)
*
dl
;
gamma
=
A
(
i
,
j
)
*
(
dl
+
dij
)
+
B
(
i
,
j
);
lambda
=
A
(
i
,
j
)
*
dij
+
B
(
i
,
j
);
F
(
i
,
j
)
=
kelvin_primitive
(
0
,
alpha
,
beta
,
gamma
,
lambda
)
-
kelvin_primitive
(
-
1
,
alpha
,
beta
,
gamma
,
lambda
);
}
}
// Integrate [0, 1] (also upper domain)
for
(
UInt
i
=
0
;
i
<
dim
;
++
i
)
{
for
(
UInt
j
=
0
;
j
<
dim
;
++
j
)
{
beta
=
-
A
(
i
,
j
)
*
dl
;
gamma
=
A
(
i
,
j
)
*
(
dl
-
dij
)
-
B
(
i
,
j
);
lambda
=
A
(
i
,
j
)
*
dij
+
B
(
i
,
j
);
F
(
i
,
j
)
=
kelvin_primitive
(
0
,
alpha
,
beta
,
gamma
,
lambda
)
-
kelvin_primitive
(
-
1
,
alpha
,
beta
,
gamma
,
lambda
);
}
}
F
*=
std
::
exp
(
dij
);
}
};
/* -------------------------------------------------------------------------- */
/// Integrate for yj = xi
template
<
UInt
dim
,
UInt
derivative
>
struct
KelvinIntegrator
<
dim
,
derivative
,
0
>
{
CUDA_LAMBDA
inline
void
integrate
(
const
Kelvin
<
dim
,
derivative
>&
kelvin
,
Real
yj
,
Real
xi
,
Real
dl
,
VectorProxy
<
Real
,
dim
-
1
>&&
q
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
F
)
{
Complex
amem
[
dim
*
dim
],
bmem
[
dim
*
dim
];
MatrixProxy
<
Complex
,
dim
,
dim
>
A
(
*
amem
),
B
(
*
bmem
);
const
Real
dij
=
yj
-
xi
;
kelvin
.
computeTerms
<
false
>
(
q
,
A
,
B
);
// Integrate [-1, 0] (lower domain)
Real
alpha
=
q
*
dl
,
beta
=
0
,
gamma
=
0
,
lambda
=
0
;
for
(
UInt
i
=
0
;
i
<
dim
;
++
i
)
{
for
(
UInt
j
=
0
;
j
<
dim
;
++
j
)
{
beta
=
A
(
i
,
j
)
*
dl
;
gamma
=
A
(
i
,
j
)
*
(
dl
+
dij
)
+
B
(
i
,
j
);
lambda
=
A
(
i
,
j
)
*
dij
+
B
(
i
,
j
);
F
(
i
,
j
)
=
kelvin_primitive
(
0
,
alpha
,
beta
,
gamma
,
lambda
)
-
kelvin_primitive
(
-
1
,
alpha
,
beta
,
gamma
,
lambda
);
}
}
// Integrate [0, 1] (upper domain)
kelvin
.
computeTerms
<
true
>
(
q
,
A
,
B
);
alpha
=
-
q
*
dl
;
for
(
UInt
i
=
0
;
i
<
dim
;
++
i
)
{
for
(
UInt
j
=
0
;
j
<
dim
;
++
j
)
{
beta
=
-
A
(
i
,
j
)
*
dl
;
gamma
=
A
(
i
,
j
)
*
(
dl
-
dij
)
-
B
(
i
,
j
);
lambda
=
A
(
i
,
j
)
*
dij
+
B
(
i
,
j
);
F
(
i
,
j
)
=
kelvin_primitive
(
0
,
alpha
,
beta
,
gamma
,
lambda
)
-
kelvin_primitive
(
-
1
,
alpha
,
beta
,
gamma
,
lambda
);
}
}
F
*=
std
::
exp
(
dij
);
}
};
/* -------------------------------------------------------------------------- */
/// Integrate for yj > xi
template
<
UInt
dim
,
UInt
derivative
>
struct
KelvinIntegrator
<
dim
,
derivative
,
-
1
>
{
CUDA_LAMBDA
inline
void
integrate
(
const
Kelvin
<
dim
,
derivative
>&
kelvin
,
Real
yj
,
Real
xi
,
Real
dl
,
VectorProxy
<
Real
,
dim
-
1
>&&
q
,
MatrixProxy
<
Complex
,
dim
,
dim
>&&
F
)
{
Complex
amem
[
dim
*
dim
],
bmem
[
dim
*
dim
];
MatrixProxy
<
Complex
,
dim
,
dim
>
A
(
*
amem
),
B
(
*
bmem
);
const
Real
dij
=
yj
-
xi
;
kelvin
.
computeTerms
<
false
>
(
q
,
A
,
B
);
// Integrate [-1, 0] (lower domain)
Real
alpha
=
q
*
dl
,
beta
=
0
,
gamma
=
0
,
lambda
=
0
;
for
(
UInt
i
=
0
;
i
<
dim
;
++
i
)
{
for
(
UInt
j
=
0
;
j
<
dim
;
++
j
)
{
beta
=
A
(
i
,
j
)
*
dl
;
gamma
=
A
(
i
,
j
)
*
(
dl
+
dij
)
+
B
(
i
,
j
);
lambda
=
A
(
i
,
j
)
*
dij
+
B
(
i
,
j
);
F
(
i
,
j
)
=
kelvin_primitive
(
0
,
alpha
,
beta
,
gamma
,
lambda
)
-
kelvin_primitive
(
-
1
,
alpha
,
beta
,
gamma
,
lambda
);
}
}
// Integrate [0, 1] (also lower domain)
for
(
UInt
i
=
0
;
i
<
dim
;
++
i
)
{
for
(
UInt
j
=
0
;
j
<
dim
;
++
j
)
{
beta
=
-
A
(
i
,
j
)
*
dl
;
gamma
=
A
(
i
,
j
)
*
(
dl
-
dij
)
-
B
(
i
,
j
);
lambda
=
A
(
i
,
j
)
*
dij
+
B
(
i
,
j
);
F
(
i
,
j
)
=
kelvin_primitive
(
0
,
alpha
,
beta
,
gamma
,
lambda
)
-
kelvin_primitive
(
-
1
,
alpha
,
beta
,
gamma
,
lambda
);
}
}
F
*=
std
::
exp
(
dij
);
}
};
}
// namespace influence
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
Event Timeline
Log In to Comment