2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4 // + This file is part of enGrid. +
6 // + Copyright 2008,2009 Oliver Gloth +
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. +
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. +
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/>. +
21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
23 #include "optimisation.h"
25 Optimisation::Optimisation()
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];
37 void Optimisation::computeDerivatives(vec3_t x
)
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
);
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
);
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
);
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
);
83 vec3_t
Optimisation::optimise(vec3_t x
)
85 computeDerivatives(x
);
86 mat3_t M
= J
.inverse();
87 return (-1)*(M
*grad_f
);