Homec4science
Diffusion LibMultiScale (features/cadd) Browse

libmultiscale/src/fmm/radial_function.hh
88abcabc76bafeatures/cadd

Last Change
Raw File

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__ */
c4science · Help