Polynomial.hcc

Go to the documentation of this file.

//Copyright (c) 2013-2014, The MercuryDPM Developers Team. All rights reserved.
//For the list of developers, see <http://www.MercuryDPM.org/Team>.
//
//Redistribution and use in source and binary forms, with or without
//modification, are permitted provided that the following conditions are met:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name MercuryDPM nor the
//    names of its contributors may be used to endorse or promote products
//    derived from this software without specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
//ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
//WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
//DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
//LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
//ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
//(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
//SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

template <StatType T>
void NORMALIZED_POLYNOMIAL<T>::set_polynomial( std::vector<double> new_coefficients, unsigned int new_dim) {    
    coefficients = new_coefficients;
    dim = new_dim;
    finish_set_polynomial();
}

template <StatType T>
void NORMALIZED_POLYNOMIAL<T>::set_polynomial(double* new_coefficients, unsigned int num_coeff, unsigned int new_dim) {
    coefficients.resize(num_coeff); 
    for (unsigned int i=0; i<num_coeff; i++)
        coefficients[i] = new_coefficients[i];
    dim = new_dim;
    finish_set_polynomial();
}

template <StatType T>
void NORMALIZED_POLYNOMIAL<T>::finish_set_polynomial() {
    //normalizes the Polynomial.
    double volume = get_volume();
    if (fabs(volume-1.0)>1e-12) 
        std::cout << "Warning: Polynomial has to be normalized by " << volume << std::endl;
    for (unsigned int i=0; i<coefficients.size(); i++) 
        coefficients[i] /= volume;
    
    //sets #averaged_coefficients   
    set_average();
    
    //std::couts the Polynomial
    std::stringstream sstr;
    sstr << *this;
    set_name(sstr.str().c_str());
    std::cout << "p=" << *this << std::endl;
}

template <StatType T>
double NORMALIZED_POLYNOMIAL<T>::get_volume() {
    double volume = 0;
    unsigned int N = coefficients.size();
    if (dim==3) {
        for (unsigned int i=0; i<N; i++)
            volume += coefficients[i]/(2.+N-i);
        volume *= 4.*constants::pi;
    } else if (dim==2) {
        std::cerr << "dim=2 is not working yet" << std::endl;   exit(-1);
        for (unsigned int i=0; i<coefficients.size(); i++)
            volume += coefficients[i]/(1.+N-i);
        volume *= 2.*constants::pi;
    } else if (dim==1) {
        std::cerr << "dim=1 is not working yet" << std::endl;   exit(-1);
        for (unsigned int i=0; i<coefficients.size(); i++)
            volume += coefficients[i]/(0.+N-i);
        volume *= 2.;
    } else {
        std::cerr << "Error in get_volume: dim=" << dim << std::endl; exit(-1);
    }
    return volume;
}

template <StatType T>
void NORMALIZED_POLYNOMIAL<T>::set_average() {std::cout << "set_average" << std::endl;}

template <StatType T>
void NORMALIZED_POLYNOMIAL<T>::set_average_1D() {
    std::cout << "set_average_1D" << std::endl;
    unsigned int N = coefficients.size();
    averaged_coefficients.resize(N+2);
    for (unsigned int i=0; i<N+2; i++) {
        averaged_coefficients[i] = 0;
    }
    for (unsigned int i=0; i<N; i++) {
        averaged_coefficients[N+1] += coefficients[i]*2*constants::pi/(N-1-i+2);
        averaged_coefficients[i] -= coefficients[i]*2*constants::pi/(N-1-i+2);
    }
}

//not needed
template <StatType T>
void NORMALIZED_POLYNOMIAL<T>::set_average_2D() {
}

//Returns the value of the Polynomial \f$p(r)=\sum_{i=0}^{N-1) c_i r^{N-1-i}\f$.
template <StatType T>
double NORMALIZED_POLYNOMIAL<T>::evaluate(double r) {
    unsigned int N = coefficients.size();
    double value = coefficients[0];
    for (unsigned int i=1; i<N; i++)
        value = value*r+coefficients[i];
    return value;
} 

//Returns the value of the Polynomial \f$p(r)=\sum_{i=0}^{N-1) (N-1-i) c_i r^{N-2-i}\f$.
template <StatType T>
double NORMALIZED_POLYNOMIAL<T>::evaluateGradient(double r) {
    unsigned int N = coefficients.size();
    double value = (N-1)*coefficients[0];
    for (unsigned int i=1; i<N-1; i++)
        value = value*r+(N-1-i)*coefficients[i];
    return value;
} 

template <StatType T>
double NORMALIZED_POLYNOMIAL<T>::evaluate_1D(double r) {
    unsigned int N = averaged_coefficients.size();
    double value = averaged_coefficients[0];
    for (unsigned int i=1; i<N; i++)
        value = value*r+averaged_coefficients[i];
    return value;
}

template <StatType T>
double NORMALIZED_POLYNOMIAL<T>::evaluateGradient_1D(double r) {
    unsigned int N = coefficients.size();
    double value = coefficients[0];
    double value2 = coefficients[0];
    for (unsigned int i=1; i<N-1; i++) {
        value = value*r+coefficients[i];
        value2 += coefficients[i];
    }
    return 2.*constants::pi*r*(value2-value*r);
}

template <StatType T>
double NORMALIZED_POLYNOMIAL<T>::evaluateGradient_2D(double) {
    return 0.0;
}


template <StatType T>
double NORMALIZED_POLYNOMIAL<T>::evaluate_2D(double r) {
    //~ std::cout << "eval_2D" << std::endl;
    unsigned int N = coefficients.size();
    double value1 = 2.*coefficients[N-1];
    double value2 = 0;
    double r2 = r*r;
    switch (N) {
        case 7:
            value1 += coefficients[N-7]*2./35.*(5.+r2*(6.+r2*(8*r2+16)));
        case 6:
            value1 += coefficients[N-6]*(8.+(10.+15.*r2)*r2)/24.;
            value2 += coefficients[N-6]*(15.*r2*r2*r2)/24.;
        case 5:
            value1 += coefficients[N-5]*.4*(1.+r2*(4./3.+r2*8./3.));
        case 4:
            value1 += coefficients[N-4]*(2.+r2*3.)/4.;
            value2 += coefficients[N-4]*.75*r2*r2;
        case 3:
            value1 += coefficients[N-3]*2./3.*(1.+r2*2.);
        case 2:
            value1 += coefficients[N-2]*1.;
            value2 += coefficients[N-2]*r2;
        case 1:
            break;
        default:
            std::cerr << "Error: no rules set for high-order polynomials" << std::endl; exit(-1);
    }
    double root = sqrt(1-r2);
    double invz = root/r;
    double arccsch = log(sqrt(1.+invz*invz)+invz);
    return value1*root+value2*arccsch;
}

template <StatType T>
double NORMALIZED_POLYNOMIAL<T>::evaluateIntegral(double a, double b, double t) {
    unsigned int N = coefficients.size();
    double value1 = coefficients[N-1]*(b-a);
    double value2 = 0;
    double t2 = t*t;
    double a2 = a*a;
    double b2 = b*b;
    double roota = sqrt(a*a+t2);
    double rootb = sqrt(b*b+t2);
    switch (N) {
        case 7:
            value1 += coefficients[N-7]*(b*(b2*b2*b2/7.0+0.6*b2*b2*t2+b2*t2*t2+t2*t2*t2)-a*(a2*a2*a2/7.0+0.6*a2*a2*t2+a2*t2*t2+t2*t2*t2));
        case 6:
            value1 += coefficients[N-6]*(b*rootb*(8.*b2*b2+26.*b2*t2+33.*t2*t2)-a*roota*(8.*a2*a2+26.*a2*t2+33.*t2*t2))/48.;
            value2 += coefficients[N-6]*15./48.*t2*t2*t2;
        case 5:
            value1 += coefficients[N-5]*(b*(b2*b2/5.+2./3.*b2*t2+t2*t2)-a*(a2*a2/5.+2./3.*a2*t2+t2*t2));
        case 4:
            value1 += coefficients[N-4]*(b*rootb*(2.*b2+5.*t2)-a*roota*(2.*a2+5.*t2))/8.;
            value2 += coefficients[N-4]*3./8.*t2*t2;
        case 3:
            value1 += coefficients[N-3]*(b*(b2/3.+t2)-a*(a2/3.+t2));
        case 2:
            value1 += coefficients[N-2]*(b*rootb-a*roota)/2.;
            value2 += coefficients[N-2]*t2/2.;
        case 1:
            break;
        default:
            std::cerr << "Error: no rules set for high-order polynomials" << std::endl; exit(-1);
    }
    double A = a/t;
    double B = b/t;
    double arcsinh = t<1e-12 ? 0.0 : log(B+sqrt(1.+B*B))-log(A+sqrt(1.+A*A));
    //~ std::cout 
    //~ << "value1 " << value1 << std::endl
    //~ << "value2 " << value2 << std::endl
    //~ << "arcsinh " << arcsinh << std::endl
    //~ << "a " << a << std::endl
    //~ << "b " << b << std::endl
    //~ << "t " << t << std::endl;
    return value1+value2*arcsinh;
} 
    
template <StatType T>
double NORMALIZED_POLYNOMIAL<T>::evaluateIntegral_2D(double a UNUSED, double b UNUSED, double t UNUSED) {
    std::cerr << "evaluateIntegral_2D is not implemented yet" << std::endl;
    exit(-1);
}

template <StatType T>
double NORMALIZED_POLYNOMIAL<T>::evaluateIntegral_1D(double a, double b, double t) {
    unsigned int N = averaged_coefficients.size();
    double value1 = averaged_coefficients[N-1]*(b-a);
    double value2 = 0;
    double t2 = t*t;
    double a2 = a*a;
    double b2 = b*b;
    double roota = sqrt(a*a+t2);
    double rootb = sqrt(b*b+t2);
    switch (N) {
        case 7:
            value1 += averaged_coefficients[N-7]*(b*(b2*b2*b2/7.0+0.6*b2*b2*t2+b2*t2*t2+t2*t2*t2)-a*(a2*a2*a2/7.0+0.6*a2*a2*t2+a2*t2*t2+t2*t2*t2));
        case 6:
            value1 += averaged_coefficients[N-6]*(b*rootb*(8.*b2*b2+26.*b2*t2+33.*t2*t2)-a*roota*(8.*a2*a2+26.*a2*t2+33.*t2*t2))/48.;
            value2 += averaged_coefficients[N-6]*15./48.*t2*t2*t2;
        case 5:
            value1 += averaged_coefficients[N-5]*(b*(b2*b2/5.+2./3.*b2*t2+t2*t2)-a*(a2*a2/5.+2./3.*a2*t2+t2*t2));
        case 4:
            value1 += averaged_coefficients[N-4]*(b*rootb*(2.*b2+5.*t2)-a*roota*(2.*a2+5.*t2))/8.;
            value2 += averaged_coefficients[N-4]*3./8.*t2*t2;
        case 3:
            value1 += averaged_coefficients[N-3]*(b*(b2/3.+t2)-a*(a2/3.+t2));
        case 2:
            value1 += averaged_coefficients[N-2]*(b*rootb-a*roota)/2.;
            value2 += averaged_coefficients[N-2]*t2/2.;
        case 1:
            break;
        default:
            std::cerr << "Error: no rules set for high-order polynomials" << std::endl; exit(-1);
    }
    double A = a/t;
    double B = b/t;
    double arcsinh = t<1e-12 ? 0.0 : log(B+sqrt(1.+B*B))-log(A+sqrt(1.+A*A));
    return value1+value2*arcsinh;
    
    // 3D Heaviside:  ( b-a );
    // 3D Polynomial: evaluateIntegral(max(a,-wn)/w,min(b,wn)/w,tangential.GetLength()/w)*w;
    // 1D Polynomial: evaluateIntegral(max(a,-wn)/w,min(b,wn)/w,tangential.GetLength()/w)*w;
    // 1D Heaviside:  ((b*pi*(wn2-sqr(b)/3.0)))-((a*pi*(wn2-sqr(a)/3.0)));
        //~ averaged_coefficients[N+1] += coefficients[i]*2*pi/(N-1-i+2);
        //~ averaged_coefficients[i] -= coefficients[i]*2*pi/(N-1-i+2);

} 

template<> void NORMALIZED_POLYNOMIAL<X>::set_average() {set_average_1D();}
template<> void NORMALIZED_POLYNOMIAL<Y>::set_average() {set_average_1D();}
template<> void NORMALIZED_POLYNOMIAL<Z>::set_average() {set_average_1D();}
template<> void NORMALIZED_POLYNOMIAL<XY>::set_average() {set_average_2D();}
template<> void NORMALIZED_POLYNOMIAL<XZ>::set_average() {set_average_2D();}
template<> void NORMALIZED_POLYNOMIAL<YZ>::set_average() {set_average_2D();}

template<> double NORMALIZED_POLYNOMIAL<X>::evaluate(double r) {return evaluate_1D(r);}
template<> double NORMALIZED_POLYNOMIAL<Y>::evaluate(double r) {return evaluate_1D(r);}
template<> double NORMALIZED_POLYNOMIAL<Z>::evaluate(double r) {return evaluate_1D(r);}
template<> double NORMALIZED_POLYNOMIAL<XY>::evaluate(double r) {return evaluate_2D(r);}
template<> double NORMALIZED_POLYNOMIAL<XZ>::evaluate(double r) {return evaluate_2D(r);}
template<> double NORMALIZED_POLYNOMIAL<YZ>::evaluate(double r) {return evaluate_2D(r);}

template<> double NORMALIZED_POLYNOMIAL<X>::evaluateGradient(double r) {return evaluateGradient_1D(r);}
template<> double NORMALIZED_POLYNOMIAL<Y>::evaluateGradient(double r) {return evaluateGradient_1D(r);}
template<> double NORMALIZED_POLYNOMIAL<Z>::evaluateGradient(double r) {return evaluateGradient_1D(r);}
template<> double NORMALIZED_POLYNOMIAL<XY>::evaluateGradient(double r) {return evaluateGradient_2D(r);}
template<> double NORMALIZED_POLYNOMIAL<XZ>::evaluateGradient(double r) {return evaluateGradient_2D(r);}
template<> double NORMALIZED_POLYNOMIAL<YZ>::evaluateGradient(double r) {return evaluateGradient_2D(r);}

template<> double NORMALIZED_POLYNOMIAL<X>::evaluateIntegral(double a, double b, double t) {return evaluateIntegral_1D(a,b,t);}
template<> double NORMALIZED_POLYNOMIAL<Y>::evaluateIntegral(double a, double b, double t) {return evaluateIntegral_1D(a,b,t);}
template<> double NORMALIZED_POLYNOMIAL<Z>::evaluateIntegral(double a, double b, double t) {return evaluateIntegral_1D(a,b,t);}
template<> double NORMALIZED_POLYNOMIAL<XY>::evaluateIntegral(double a, double b, double t) {return evaluateIntegral_2D(a,b,t);}
template<> double NORMALIZED_POLYNOMIAL<XZ>::evaluateIntegral(double a, double b, double t) {return evaluateIntegral_2D(a,b,t);}
template<> double NORMALIZED_POLYNOMIAL<YZ>::evaluateIntegral(double a, double b, double t) {return evaluateIntegral_2D(a,b,t);}