optimisation.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008,2009 Oliver Gloth                                     +
   7 // +                                                                      +
   8 // + enGrid is free software: you can redistribute it and/or modify       +
   9 // + it under the terms of the GNU General Public License as published by +
  10 // + the Free Software Foundation, either version 3 of the License, or    +
  11 // + (at your option) any later version.                                  +
  12 // +                                                                      +
  13 // + enGrid is distributed in the hope that it will be useful,            +
  14 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  15 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  16 // + GNU General Public License for more details.                         +
  17 // +                                                                      +
  18 // + You should have received a copy of the GNU General Public License    +
  19 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  20 // +                                                                      +
  21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  22 //
  23 #include "optimisation.h"
  24
  25 Optimisation::Optimisation()
  26 {
  27   setDeltas(1e-6);
  28   F = new double**[3];
  29   for (int i = 0; i < 3; ++i) {
  30     F[i] = new double*[3];
  31     for (int j = 0; j < 3; ++j) {
  32       F[i][j] = new double[3];
  33     };
  34   };
  35 };
  36
  37 void Optimisation::computeDerivatives(vec3_t x)
  38 {
  39   F[0][0][0] = func(x[0]-Dx, x[1]-Dy, x[2]-Dz);
  40   F[1][0][0] = func(x[0],    x[1]-Dy, x[2]-Dz);
  41   F[2][0][0] = func(x[0]+Dx, x[1]-Dy, x[2]-Dz);
  42   F[0][1][0] = func(x[0]-Dx, x[1],    x[2]-Dz);
  43   F[1][1][0] = func(x[0],    x[1],    x[2]-Dz);
  44   F[2][1][0] = func(x[0]+Dx, x[1],    x[2]-Dz);
  45   F[0][2][0] = func(x[0]-Dx, x[1]+Dy, x[2]-Dz);
  46   F[1][2][0] = func(x[0],    x[1]+Dy, x[2]-Dz);
  47   F[2][2][0] = func(x[0]+Dx, x[1]+Dy, x[2]-Dz);
  48   F[0][0][1] = func(x[0]-Dx, x[1]-Dy, x[2]);
  49   F[1][0][1] = func(x[0],    x[1]-Dy, x[2]);
  50   F[2][0][1] = func(x[0]+Dx, x[1]-Dy, x[2]);
  51   F[0][1][1] = func(x[0]-Dx, x[1],    x[2]);
  52   F[1][1][1] = func(x[0],    x[1],    x[2]);
  53   F[2][1][1] = func(x[0]+Dx, x[1],    x[2]);
  54   F[0][2][1] = func(x[0]-Dx, x[1]+Dy, x[2]);
  55   F[1][2][1] = func(x[0],    x[1]+Dy, x[2]);
  56   F[2][2][1] = func(x[0]+Dx, x[1]+Dy, x[2]);
  57   F[0][0][2] = func(x[0]-Dx, x[1]-Dy, x[2]+Dz);
  58   F[1][0][2] = func(x[0],    x[1]-Dy, x[2]+Dz);
  59   F[2][0][2] = func(x[0]+Dx, x[1]-Dy, x[2]+Dz);
  60   F[0][1][2] = func(x[0]-Dx, x[1],    x[2]+Dz);
  61   F[1][1][2] = func(x[0],    x[1],    x[2]+Dz);
  62   F[2][1][2] = func(x[0]+Dx, x[1],    x[2]+Dz);
  63   F[0][2][2] = func(x[0]-Dx, x[1]+Dy, x[2]+Dz);
  64   F[1][2][2] = func(x[0],    x[1]+Dy, x[2]+Dz);
  65   F[2][2][2] = func(x[0]+Dx, x[1]+Dy, x[2]+Dz);
  66
  67   grad_f[0] = (F[2][1][1]-F[0][1][1])/(2*Dx);
  68   grad_f[1] = (F[1][2][1]-F[1][0][1])/(2*Dy);
  69   grad_f[2] = (F[1][1][2]-F[1][1][0])/(2*Dz);
  70
  71   J[0][0] = (F[0][1][1]-2*F[1][1][1]+F[2][1][1])/(Dx*Dx);
  72   J[1][1] = (F[1][0][1]-2*F[1][1][1]+F[1][2][1])/(Dy*Dy);
  73   J[2][2] = (F[1][1][0]-2*F[1][1][1]+F[1][1][2])/(Dz*Dz);
  74
  75   J[0][1] = ((F[2][2][1]-F[0][2][1]) - (F[2][0][1]-F[0][0][1]))/(4*Dx*Dy);
  76   J[0][2] = ((F[2][1][2]-F[0][1][2]) - (F[2][1][0]-F[0][1][0]))/(4*Dx*Dz);
  77   J[1][2] = ((F[1][2][2]-F[1][2][2]) - (F[1][0][0]-F[1][0][0]))/(4*Dy*Dz);
  78   J[1][0] = J[0][1];
  79   J[2][0] = J[0][2];
  80   J[2][1] = J[1][2];
  81 };
  82
  83 vec3_t Optimisation::optimise(vec3_t x)
  84 {
  85   computeDerivatives(x);
  86   mat3_t M = J.inverse();
  87   return (-1)*(M*grad_f);
  88 };