Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F90032617
radial_function.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
Mon, Oct 28, 16:20
Size
10 KB
Mime Type
text/x-c++
Expires
Wed, Oct 30, 16:20 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
21995792
Attached To
rLIBMULTISCALE LibMultiScale
radial_function.hh
View Options
/**
* @file radial_function.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Thu Jan 10 12:24:40 2013
*
* @brief Generate a radial function from a 1d function
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale 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.
*
* LibMultiScale 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 LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef __LIBMULTISCALE_RADIAL_FUNCTION_HH__
#define __LIBMULTISCALE_RADIAL_FUNCTION_HH__
/* -------------------------------------------------------------------------- */
#include "../math/function_interface.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
template
<
UInt
Dim
,
typename
myFunc
>
class
RadialFunction
:
public
FunctionInterface
{
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public
:
RadialFunction
(
myFunc
&
f
)
:
myF
(
f
){};
~
RadialFunction
(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
//! compute the value of the function derived partially d_x, d_y and d_z times
//! in x,y,z
inline
Real
compute
(
Real
*
x
,
const
UInt
*
d_X
);
//! compute the value of the function itself
inline
Real
compute
(
Real
*
x
);
private
:
//! compute first order derivative
inline
Real
derivationOrder1
(
Real
*
x
,
UInt
i
);
//! compute second order derivative
inline
Real
derivationOrder2
(
Real
*
x
,
UInt
i
);
//! compute second order derivative
inline
Real
derivationOrder11
(
Real
*
x
,
UInt
i
,
UInt
j
);
//! compute third order derivative
inline
Real
derivationOrder3
(
Real
*
x
,
UInt
i
);
//! compute third order derivative
inline
Real
derivationOrder21
(
Real
*
x
,
UInt
i
,
UInt
j
);
//! compute third order derivative
inline
Real
derivationOrder111
(
Real
*
x
);
//! compute the norm of vector (x,y,z)
inline
Real
computeNorm
(
Real
*
x
);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private
:
myFunc
myF
;
};
/* -------------------------------------------------------------------------- */
template
<
UInt
Dim
,
typename
myFunc
>
inline
Real
RadialFunction
<
Dim
,
myFunc
>::
computeNorm
(
Real
*
x
)
{
UInt
res
=
0.0
;
for
(
UInt
i
=
0
;
i
<
Dim
;
++
i
)
{
res
+=
x
[
i
]
*
x
[
i
];
}
return
sqrt
(
res
);
}
/* -------------------------------------------------------------------------- */
template
<
UInt
Dim
,
typename
myFunc
>
inline
Real
RadialFunction
<
Dim
,
myFunc
>::
compute
(
Real
*
x
)
{
Real
r
=
computeNorm
(
x
);
return
myF
.
template
compute
<
0
>
(
r
);
}
/* -------------------------------------------------------------------------- */
template
<
UInt
Dim
,
typename
myFunc
>
inline
Real
RadialFunction
<
Dim
,
myFunc
>::
compute
(
Real
*
x
,
const
UInt
*
d_X
)
{
// order 0
if
(
d_X
[
0
]
==
0
&&
d_X
[
1
]
==
0
&&
d_X
[
2
]
==
0
)
return
compute
(
x
);
// order 1
else
if
(
d_X
[
0
]
==
1
&&
d_X
[
1
]
==
0
&&
d_X
[
2
]
==
0
)
return
derivationOrder1
(
x
,
0
);
else
if
(
d_X
[
0
]
==
0
&&
d_X
[
1
]
==
1
&&
d_X
[
2
]
==
0
)
return
derivationOrder1
(
x
,
1
);
else
if
(
d_X
[
0
]
==
0
&&
d_X
[
1
]
==
0
&&
d_X
[
2
]
==
1
)
return
derivationOrder1
(
x
,
2
);
// order 2
else
if
(
d_X
[
0
]
==
2
&&
d_X
[
1
]
==
0
&&
d_X
[
2
]
==
0
)
return
derivationOrder2
(
x
,
0
);
else
if
(
d_X
[
0
]
==
1
&&
d_X
[
1
]
==
1
&&
d_X
[
2
]
==
0
)
return
derivationOrder11
(
x
,
0
,
1
);
else
if
(
d_X
[
0
]
==
1
&&
d_X
[
1
]
==
0
&&
d_X
[
2
]
==
1
)
return
derivationOrder11
(
x
,
0
,
2
);
else
if
(
d_X
[
0
]
==
0
&&
d_X
[
1
]
==
2
&&
d_X
[
2
]
==
0
)
return
derivationOrder2
(
x
,
1
);
else
if
(
d_X
[
0
]
==
0
&&
d_X
[
1
]
==
1
&&
d_X
[
2
]
==
1
)
return
derivationOrder11
(
x
,
1
,
1
);
else
if
(
d_X
[
0
]
==
0
&&
d_X
[
1
]
==
0
&&
d_X
[
2
]
==
2
)
return
derivationOrder2
(
x
,
2
);
// order 3
else
if
(
d_X
[
0
]
==
3
&&
d_X
[
1
]
==
0
&&
d_X
[
2
]
==
0
)
return
derivationOrder3
(
x
,
0
);
else
if
(
d_X
[
0
]
==
2
&&
d_X
[
1
]
==
1
&&
d_X
[
2
]
==
0
)
return
derivationOrder21
(
x
,
0
,
1
);
else
if
(
d_X
[
0
]
==
2
&&
d_X
[
1
]
==
0
&&
d_X
[
2
]
==
1
)
return
derivationOrder21
(
x
,
0
,
2
);
else
if
(
d_X
[
0
]
==
1
&&
d_X
[
1
]
==
2
&&
d_X
[
2
]
==
0
)
return
derivationOrder21
(
x
,
1
,
0
);
else
if
(
d_X
[
0
]
==
1
&&
d_X
[
1
]
==
1
&&
d_X
[
2
]
==
1
)
return
derivationOrder111
(
x
);
else
if
(
d_X
[
0
]
==
1
&&
d_X
[
1
]
==
0
&&
d_X
[
2
]
==
2
)
return
derivationOrder21
(
x
,
2
,
0
);
else
if
(
d_X
[
0
]
==
0
&&
d_X
[
1
]
==
3
&&
d_X
[
2
]
==
0
)
return
derivationOrder3
(
x
,
1
);
else
if
(
d_X
[
0
]
==
0
&&
d_X
[
1
]
==
2
&&
d_X
[
2
]
==
1
)
return
derivationOrder21
(
x
,
1
,
2
);
else
if
(
d_X
[
0
]
==
0
&&
d_X
[
1
]
==
1
&&
d_X
[
2
]
==
2
)
return
derivationOrder21
(
x
,
2
,
1
);
else
if
(
d_X
[
0
]
==
0
&&
d_X
[
1
]
==
0
&&
d_X
[
2
]
==
3
)
return
derivationOrder3
(
x
,
2
);
else
LM_FATAL
(
"this combination of partial derivatives still needs to be "
"implemented"
);
}
/* -------------------------------------------------------------------------- */
/** @f[
* \frac{\partial \phi(r)}{\partial x_i} = \phi'(r) \frac{x_i}{r}
* @f]
*/
template
<
UInt
Dim
,
typename
myFunc
>
inline
Real
RadialFunction
<
Dim
,
myFunc
>::
derivationOrder1
(
Real
*
x
,
UInt
i
)
{
if
(
Dim
<
3
&&
i
==
2
)
return
0.0
;
if
(
Dim
==
1
&&
i
==
1
)
return
0.0
;
Real
r
=
computeNorm
(
x
);
return
x
[
i
]
/
r
*
myF
.
template
compute
<
1
>
(
r
);
}
/* -------------------------------------------------------------------------- */
/** @f[
* \frac{\partial^2 \phi(r)}{\partial x_i^2} =
* \phi''(r) \frac{x_i^2}{r^2} + \frac{\phi'(r)}{r} (1 - \frac{x_i^2}{r^2})
* @f]
*/
template
<
UInt
Dim
,
typename
myFunc
>
inline
Real
RadialFunction
<
Dim
,
myFunc
>::
derivationOrder2
(
Real
*
x
,
UInt
i
)
{
if
(
Dim
<
3
&&
i
==
2
)
return
0.0
;
if
(
Dim
==
1
&&
i
==
1
)
return
0.0
;
Real
r
=
computeNorm
(
x
);
Real
_1_r
=
1.
/
r
;
Real
_1_r2
=
_1_r
*
_1_r
;
Real
_x2
=
x
[
i
]
*
x
[
i
];
Real
_x2_r2
=
_x2
*
_1_r2
;
return
(
myF
.
template
compute
<
2
>
(
r
)
*
_x2_r2
+
myF
.
template
compute
<
1
>
(
r
)
*
_1_r
*
(
1
-
_x2_r2
));
}
/* -------------------------------------------------------------------------- */
/** @f[
* \frac{\partial^2 \phi(r)}{\partial x_i \partial x_j} =
* \phi''(r) \frac{x_i x_j}{r^2} - \phi'(r) \frac{x_i x_j}{r^3}
* @f]
*/
template
<
UInt
Dim
,
typename
myFunc
>
inline
Real
RadialFunction
<
Dim
,
myFunc
>::
derivationOrder11
(
Real
*
x
,
UInt
i
,
UInt
j
)
{
if
(
Dim
<
3
&&
(
i
==
2
||
j
==
2
))
return
0.0
;
if
(
Dim
==
1
&&
(
i
==
1
||
j
==
1
))
return
0.0
;
Real
r
=
computeNorm
(
x
);
Real
_1_r
=
1.
/
r
;
Real
_1_r2
=
_1_r
*
_1_r
;
Real
_xy
=
x
[
i
]
*
x
[
j
];
Real
_xy_r2
=
_xy
*
_1_r2
;
Real
_xy_r3
=
_xy_r2
*
_1_r
;
return
(
myF
.
template
compute
<
2
>
(
r
)
*
_xy_r2
-
myF
.
template
compute
<
1
>
(
r
)
*
_xy_r3
);
}
/* -------------------------------------------------------------------------- */
/** @f[
* \frac{\partial^3 \phi(r)}{\partial x_i^3} =
* \phi'''(r) \frac{x_i^3}{r^3} + 3\phi''(r) (\frac{x_i}{r^2} -
* \frac{x_i^3}{r^4})
* + 3 \phi'(r) ( \frac{x_i^3}{r^5} - \frac{x_i}{r^3} )
* @f]
* which can be factorized as
* @f[
* \frac{\partial^3 \phi(r)}{\partial x_i^3} =
* \phi'''(r) \frac{x_i^3}{r^3} + 3(1-\frac{x_i^2}{r^2}) ( \phi''(r)
* \frac{x_i}{r^2} - \phi'(r) \frac{x_i}{r^3} ) =
* \phi'''(r) \frac{x_i^3}{r^3} + 3(1-\frac{x_i^2}{r^2}) \frac{x_i}{r^2}(
* \phi''(r) - \phi'(r) \frac{1}{r} )
* @f]
*/
template
<
UInt
Dim
,
typename
myFunc
>
inline
Real
RadialFunction
<
Dim
,
myFunc
>::
derivationOrder3
(
Real
*
x
,
UInt
i
)
{
if
(
Dim
<
3
&&
i
==
2
)
return
0.0
;
if
(
Dim
==
1
&&
i
==
1
)
return
0.0
;
Real
r
=
computeNorm
(
x
);
Real
_1_r
=
1.
/
r
;
Real
_1_r2
=
_1_r
*
_1_r
;
Real
_1_r3
=
_1_r
*
_1_r2
;
Real
_x
=
x
[
i
];
Real
_x2
=
_x
*
_x
;
Real
_x3
=
_x
*
_x2
;
return
myF
.
template
compute
<
3
>
(
r
)
*
_x3
*
_1_r3
+
3.
*
(
1.
-
_x2
*
_1_r2
)
*
_x
*
_1_r2
*
(
myF
.
template
compute
<
2
>
(
r
)
-
myF
.
template
compute
<
1
>
(
r
)
*
_1_r
);
}
/* -------------------------------------------------------------------------- */
/** @f[
* \frac{\partial^3 \phi(r)}{\partial x_i^2 \partial x_j} =
* \phi'''(r) \frac{x_i^2 x_j}{r^3} + \phi''(r) (-3 \frac{x_i^2 x_j}{r^4} +
* \frac{x_j}{r^2})
* + \phi'(r) (3 \frac{x_i^2 x_j}{r^5} - \frac{x_j}{r^3})
* @f]
* which can be factorized as
* @f[
* \frac{\partial^3 \phi(r)}{\partial x_i^2 \partial x_j} =
* \phi'''(r) \frac{x_i^2 x_j}{r^3} + (1-3\frac{x_i^2}{r^2}) (\phi''(r)
* \frac{x_j}{r^2} - \phi'(r) \frac{x_j}{r^3} )
* = \phi'''(r) \frac{x_i^2 x_j}{r^3} + (1-3\frac{x_i^2}{r^2}) \frac{x_j}{r^2}
* (\phi''(r) - \phi'(r) \frac{1}{r} )
* @f]
*/
template
<
UInt
Dim
,
typename
myFunc
>
inline
Real
RadialFunction
<
Dim
,
myFunc
>::
derivationOrder21
(
Real
*
x
,
UInt
i
,
UInt
j
)
{
if
(
Dim
<
3
&&
(
i
==
2
||
j
==
2
))
return
0.0
;
if
(
Dim
==
1
&&
(
i
==
1
||
j
==
1
))
return
0.0
;
Real
r
=
computeNorm
(
x
);
Real
_1_r
=
1.
/
r
;
Real
_1_r2
=
_1_r
*
_1_r
;
Real
_1_r3
=
_1_r
*
_1_r2
;
Real
_x
=
x
[
i
];
Real
_y
=
x
[
j
];
Real
_x2
=
_x
*
_x
;
return
myF
.
template
compute
<
3
>
(
r
)
*
_x2
*
_y
*
_1_r3
+
(
1.
-
3.
*
_x2
*
_1_r2
)
*
_y
*
_1_r2
*
(
myF
.
template
compute
<
2
>
(
r
)
-
myF
.
template
compute
<
1
>
(
r
)
*
_1_r
);
}
/* -------------------------------------------------------------------------- */
/** @f[
* \frac{\partial^3 \phi(r)}{\partial x \partial y \partial z} =
* \frac{xyz}{r^3} ( \phi'''(r) -3 \phi''(r) \frac{1}{r} + 3 \phi'(r)
* \frac{1}{r^2})
* @f]
*/
template
<
UInt
Dim
,
typename
myFunc
>
inline
Real
RadialFunction
<
Dim
,
myFunc
>::
derivationOrder111
(
Real
*
x
)
{
if
(
Dim
<
3
)
return
0.0
;
Real
r
=
computeNorm
(
x
);
Real
_1_r
=
1.
/
r
;
Real
_1_r2
=
_1_r
*
_1_r
;
Real
_1_r3
=
_1_r
*
_1_r2
;
Real
_xyz
=
x
[
0
]
*
x
[
1
]
*
x
[
2
];
return
_xyz
*
_1_r3
*
(
myF
.
template
compute
<
3
>
(
r
)
-
3.
*
myF
.
template
compute
<
2
>
(
r
)
*
_1_r
+
3.
*
myF
.
template
compute
<
1
>
(
r
)
*
_1_r2
);
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif
/* __LIBMULTISCALE_RADIAL_FUNCTION_HH__ */
Event Timeline
Log In to Comment