src/ODE/findRoot/BisectionRoot/BisectionRoot.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Class
  26     BisectionRoot
  27
  28 Description
  29     Bisection root Based on Numerical Recipes in C++, Section 9.1, page 358.
  30     Function is provided as a template parameter function object, evaluated
  31     using operator()(const scalar x)
  32
  33 Author
  34     Hrvoje Jasak, Wikki Ltd.  All rights reserved.
  35
  36 \*----------------------------------------------------------------------------*/
  37
  38 #include "BisectionRoot.H"
  39
  40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  41
  42 template<class Func>
  43 const Foam::label Foam::BisectionRoot<Func>::maxIter = 60;
  44
  45
  46 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  47
  48 template<class Func>
  49 Foam::BisectionRoot<Func>::BisectionRoot(const Func& f, const scalar eps)
  50 :
  51     f_(f),
  52     eps_(eps)
  53 {}
  54
  55 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  56
  57 template<class Func>
  58 Foam::scalar Foam::BisectionRoot<Func>::root
  59 (
  60     const scalar x0,
  61     const scalar x1
  62 ) const
  63 {
  64     scalar f, fMid, dx, rtb, xMid;
  65
  66     f = f_(x0);
  67     fMid = f_(x1);
  68
  69     if (f*fMid >= 0)
  70     {
  71         FatalErrorIn
  72         (
  73             "Foam::scalar Foam::BisectionRoot<Func>::root\n"
  74             "(\n"
  75             "    const scalar x0,\n"
  76             "    const scalar x1\n"
  77             ") const"
  78         )   << "Root is not bracketed.  f(x0) = " << f << " f(x1) = " << fMid
  79             << abort(FatalError);
  80     }
  81
  82     // Orient the search such that f > 0 lies at x + dx
  83     if (f < 0)
  84     {
  85         dx = x1 - x0;
  86         rtb = x0;
  87     }
  88     else
  89     {
  90         dx = x0 - x1;
  91         rtb = x1;
  92     }
  93
  94     for (label nIter = 0; nIter < maxIter; nIter++)
  95     {
  96         dx *= 0.5;
  97         xMid = rtb + dx;
  98
  99         fMid = f_(xMid);
 100
 101         if (fMid <= 0)
 102         {
 103             rtb = xMid;
 104         }
 105
 106         if (mag(dx) < eps_ || mag(fMid) < SMALL)
 107         {
 108             return rtb;
 109         }
 110     }
 111
 112     FatalErrorIn
 113     (
 114         "Foam::scalar Foam::BisectionRoot<Func>::root\n"
 115         "(\n"
 116         "    const scalar x0,\n"
 117         "    const scalar x1\n"
 118         ") const"
 119     )   << "Maximum number of iterations exceeded" << abort(FatalError);
 120
 121     // Dummy return to keep compiler happy
 122     return x0;
 123 }
 124
 125
 126 // ************************************************************************* //