Homec4science
Diffusion LibMultiScale (master) Browse

libmultiscale/src/fmm/radial_function.hh
b05d5f2aba78master

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