mathfem.h

Go to the documentation of this file.

/*
 *
 *                 #####    #####   ######  ######  ###   ###
 *               ##   ##  ##   ##  ##      ##      ## ### ##
 *              ##   ##  ##   ##  ####    ####    ##  #  ##
 *             ##   ##  ##   ##  ##      ##      ##     ##
 *            ##   ##  ##   ##  ##      ##      ##     ##
 *            #####    #####   ##      ######  ##     ##
 *
 *
 *             OOFEM : Object Oriented Finite Element Code
 *
 *               Copyright (C) 1993 - 2013   Borek Patzak
 *
 *
 *
 *       Czech Technical University, Faculty of Civil Engineering,
 *   Department of Structural Mechanics, 166 29 Prague, Czech Republic
 *
 *  This library 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 2.1 of the License, or (at your option) any later version.
 *
 *  This program 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 this library; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 */

/*
 * This file contains some functions used in the finite element application.
 * ref : Lippman p 104
 */
#ifndef mathfem_h
#define mathfem_h

#include "error.h"
#include "oofemcfg.h"

#include <cmath>
#include <cfloat> // For _isnan

namespace oofem {
class FloatArray;

#ifndef M_PI
 #define M_PI 3.1415926535897932384626433832795029L
#endif
#ifndef M_LN2
 #define M_LN2 0.6931471805599453094172321214581766L
#endif

inline int min(int i, int j)
{ return ( i <= j ? i : j ); }

inline long min(long i, long j)
{ return ( i <= j ? i : j ); }

inline double min(double i, double j)
{ return ( i <= j ? i : j ); }

inline int max(int i, int j)
{ return ( i >= j ? i : j ); }

inline double clamp(int a, int lower, int upper)
{ return ( a <= lower ? lower : ( a >= upper ? upper : a ) ); }

inline long max(long i, long j)
{ return ( i >= j ? i : j ); }

inline double max(double i, double j)
{ return ( i >= j ? i : j ); }

inline double clamp(double a, double lower, double upper)
{ return ( a <= lower ? lower : ( a >= upper ? upper : a ) ); }

inline double sgn(double i)
{ return ( i < 0. ? -1. : 1. ); }

double signum(double i);

#ifndef HAVE_ISNAN
 #ifdef _MSC_VER
inline bool isnan(double x) { return _isnan(x) != 0; }
 #else
// Last attempt to find a isnan function, rely on C++11
inline bool isnan(double x) { return std :: isnan(x); }
 #endif
#endif

#ifndef HAVE_CBRT
inline double cbrt(double x) { return sgn(x) * pow(fabs(x), 1.0 / 3.0); }
#endif

inline double sqr(double x) { return x * x; }

inline double macbra(double x) { return ( x >= 0 ? x : 0 ); }
inline double negbra(double x) { return ( x <= 0 ? x : 0 ); }

void cubic(double a, double b, double c, double d, double *r1, double *r2, double *r3, int *num);

void cubic3r(double a, double b, double c, double d, double *r1, double *r2, double *r3, int *num);

int iperm(int val, int rank);


#define MATHFEM_C 0.38196601
#define MATHFEM_R ( 1 - MATHFEM_C )
#define MATHFEM_BRENT_MAXITER 100

template< class T > class mem_fun
{
    double ( T :: *pmf )(double);
    T *ptr;
public:
    mem_fun( T * o, double ( T :: *p )(double) ) : pmf(p), ptr(o) { }
    double operator() (double x) const { return ( ptr->*pmf )(x); }
};


class OOFEM_EXPORT c_fun
{
    double ( * func )(double);
public:
    c_fun( double ( * p )(double) ) : func(p) { }
    double operator() (double x) const { return ( * func )( x ); }
};



template< class T > double gss(double ax, double bx, double cx, const T &f,
                               double tol, double &xmin)
{
    int ii = 0;
    double f1, f2, x0, x1, x2, x3;

    x0 = ax;
    x3 = cx;

    // initialization x0-x1 made the smaller segment
    if ( fabs(cx - bx) > fabs(bx - ax) ) {
        x1 = bx;
        x2 = bx + MATHFEM_C * ( cx - bx );
    } else {
        x2 = bx;
        x1 = bx - MATHFEM_C * ( bx - ax );
    }

    f1 = f(x1);
    f2 = f(x2);

    // iteration loop
    while ( fabs(x3 - x0) > tol * ( fabs(x1) + fabs(x2) ) ) {
        if ( f2 < f1 ) {
            // minimum bracketed by (x1,x2,x3) triplet
            x0 = x1;
            x1 = x2;
            x2 = MATHFEM_R * x1 + MATHFEM_C * x3; // x2=x1+C*(x3-x1)

            f1 = f2;
            f2 = f(x2);
        } else {
            // minimum bracketed by (x0,x1,x2) triplet
            x3 = x2;
            x2 = x1;
            x1 = MATHFEM_R * x2 + MATHFEM_C * x0; // x1 = x2+c*(x0-x2)

            f2 = f1;
            f1 = f(x1);
        }

        ii++;
    }

    //printf ("gss: convergence reached in %d iterations\n", ii);
    if ( f1 < f2 ) {
        xmin = x1;
        return f1;
    } else {
        xmin = x2;
        return f2;
    }
}

template< class T > double brent(double ax, double bx, double cx, const T &f,
                                 double tol, double &xmin)
{
    int ii;
    double x_left = ax, x_right = cx;
    double x, x_midpoint, v, w, u, tol1, tol2, p, q, r, e_tmp, d = 0.0, fx, fv, fw, fu;
    double e = 0.0;

    x = v = w = bx;
    fx = fv = fw = f(x);

    for ( ii = 1; ii <= MATHFEM_BRENT_MAXITER; ii++ ) {
        x_midpoint = 0.5 * ( x_left + x_right );

        // check for convergence here
        tol1 = tol * fabs(x) + 1.0e-10;
        tol2 = 2.0 * tol1;
        if ( fabs(x - x_midpoint) <= ( tol2 - 0.5 * ( x_right - x_left ) ) ) {
            //printf ("brent: convergence in %d iterations\n", ii);
            xmin = x;
            return fx;
        }

        if ( fabs(e) > tol1 ) {
            /* fit parabola */
            r = ( x - w ) * ( fx - fv );
            q = ( x - v ) * ( fx - fw );
            p = ( x - v ) * q - ( x - w ) * r;
            q = 2.0 * ( q - r );

            if ( q > 0 ) {
                p = -p;
            } else {
                q = -q;
            }

            e_tmp = e;
            e = d;

            if ( fabs(p) < fabs(0.5 * q * e_tmp) && p < q * ( x - x_left ) && p < q * ( x_right - x ) ) {
                d = p / q;
                u = x + d;
                if ( ( u - x_left ) < tol2 || ( x_right - u ) < tol2 ) {
                    d = ( x < x_midpoint ) ? tol1 : -tol1;
                }
            } else {
                e = ( x < x_midpoint ) ? x_right - x : -( x - x_left );
                d = MATHFEM_C * e;
            }
        } else {
            e = ( x < x_midpoint ) ? x_right - x : -( x - x_left );
            d = MATHFEM_C * e;
        }

        if ( fabs(d) >= tol1 ) {
            u = x + d;
        } else {
            u = x + ( ( d > 0 ) ? tol1 : -tol1 );
        }

        fu = f(u);

        if ( fu <= fx ) {
            if ( u >= x ) {
                x_left = x;
            } else {
                x_right = x;
            }

            v = w;
            w = x;
            x = u;
            fv = fw;
            fw = fx;
            fx = fu;
        } else {
            if ( u < x ) {
                x_left = u;
            } else {
                x_right = u;
            }

            if ( fu <= fw || w == x ) {
                v = w;
                w = u;
                fv = fw;
                fw = fu;
            } else if ( fu <= fv || v == x || v == w ) {
                v = u;
                fv = fu;
            }
        }
    }

    // too many iterations
    OOFEM_LOG_WARNING("brent : too many iterations\n");
    xmin = x;
    return fx;
}

void ls2fit(const FloatArray &x, const FloatArray &y, FloatArray &a);
} // end namespace oofem
#endif // mathfem_h