?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
preparing for 1.0 release
[engrid.git] / gridsmoother.cpp
blob72152df95d60b760cd42296bf11be23b2fcbb256
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 "gridsmoother.h"
24 #include "guimainwindow.h"
25
26 #include "elements.h"
27
28 #include <QTime>
29
30 GridSmoother::GridSmoother()
31 {
32 N_iterations = 5;
33 N_relaxations = 5;
34 N_boundary_corrections = 20;
35 N_search = 10;
36 L_search = 0.5;
37 smooth_prisms = true;
38 dbg = false;
39 F_old = 0;
40 F_new = 0;
41
42 getSet("boundary layer", "tetra weighting", 1.0, w_tet);
43 getSet("boundary layer", "layer height weighting", 1.0, w_h);
44 getSet("boundary layer", "parallel edges weighting", 3.0, w_par);
45 getSet("boundary layer", "parallel faces weighting", 5.0, w_n);
46 getSet("boundary layer", "similar face area weighting", 5.0, w_A);
47 getSet("boundary layer", "skewness weighting", 0.0, w_skew);
48 getSet("boundary layer", "orthogonality weighting", 0.0, w_orth);
49 getSet("boundary layer", "sharp features on nodes weighting", 8.0, w_sharp1);
50 getSet("boundary layer", "sharp features on nodes exponent", 1.3, e_sharp1);
51 getSet("boundary layer", "sharp features on edges weighting", 3.0, w_sharp2);
52 getSet("boundary layer", "sharp features on edges exponent", 1.3, e_sharp2);
53 getSet("boundary layer", "relative height of boundary layer", 1.5, H);
54 };
55
56 void GridSmoother::markNodes()
57 {
58 node_marked.fill(false,grid->GetNumberOfPoints());
59 QVector<bool> new_mark(grid->GetNumberOfPoints());
60 for (int i_iterations = 0; i_iterations < 2; ++i_iterations) {
61 qCopy(node_marked.begin(),node_marked.end(),new_mark.begin());
62 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
63 bool mark_cell = false;
64 vtkIdType type_cell, N_pts, *pts;
65 type_cell = grid->GetCellType(id_cell);
66 grid->GetCellPoints(id_cell, N_pts, pts);
67 if (type_cell == VTK_WEDGE) {
68 mark_cell = true;
69 } else {
70 for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
71 if (node_marked[pts[i_pts]]) {
72 mark_cell = true;
73 };
74 };
75 };
76 if (mark_cell) {
77 for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
78 new_mark[pts[i_pts]] = true;
79 };
80 };
81 };
82 qCopy(new_mark.begin(),new_mark.end(),node_marked.begin());
83 };
84 N_marked_nodes = 0;
85 foreach (vtkIdType id_node, nodes) {
86 if (id_node < 0) EG_BUG;
87 if (id_node > grid->GetNumberOfPoints()) EG_BUG;
88 if (node_marked[id_node]) {
89 ++N_marked_nodes;
90 };
91 };
92 };
93
94 bool GridSmoother::setNewPosition(vtkIdType id_node, vec3_t x_new)
95 {
96 using namespace GeometryTools;
97
98 vec3_t x_old;
99 grid->GetPoint(id_node, x_old.data());
100 grid->GetPoints()->SetPoint(id_node, x_new.data());
101 bool move = true;
102 Elements E;
103
104 foreach (int i_cells, n2c[id_node]) {
105 vtkIdType id_cell = cells[i_cells];
106 vtkIdType type_cell = grid->GetCellType(id_cell);
107 if (type_cell == VTK_TETRA) {
108 if (GeometryTools::cellVA(grid, id_cell) < 0) {
109 move = false;
110 //if (dbg) cout << id_node << " : tetra negative" << endl;
111 };
112 };
113 if (type_cell == VTK_WEDGE) {
114 vtkIdType N_pts, *pts;
115 vec3_t xtet[4];
116 grid->GetCellPoints(id_cell, N_pts, pts);
117 bool ok = true;
118 for (int i = 0; i < 4; ++i) { // variation
119 ok = true;
120 for (int j = 0; j < 3; ++j) { // tetrahedron
121 for (int k = 0; k < 4; ++k) { // node
122 grid->GetPoint(pts[E.priTet(i,j,k)], xtet[k].data());
123 };
124 if (GeometryTools::tetraVol(xtet[0], xtet[1], xtet[2], xtet[3]) < 0) {
125 ok = false;
126 //if (dbg) cout << id_node << " : prism negative" << endl;
127 };
128 };
129 if (ok) {
130 break;
131 };
132 };
133 if (!ok) {
134 move = false;
135 };
136 };
137 };
138 if (!move) {
139 grid->GetPoints()->SetPoint(id_node, x_old.data());
140 };
141 return move;
142 };
143
144 void GridSmoother::correctDx(int i_nodes, vec3_t &Dx)
145 {
146 for (int i_boundary_correction = 0; i_boundary_correction < N_boundary_corrections; ++i_boundary_correction) {
147 foreach (vtkIdType id_cell, n2c[i_nodes]) {
148 if (isSurface(id_cell, grid)) {
149 double A = GeometryTools::cellVA(grid, id_cell);
150 if (A > 1e-20) {
151 vec3_t n = GeometryTools::cellNormal(grid, id_cell);
152 n.normalise();
153 Dx -= (n*Dx)*n;
154 };
155 } else {
156 if (grid->GetCellType(id_cell) == VTK_WEDGE) {
157 vtkIdType N_pts, *pts;
158 grid->GetCellPoints(id_cell, N_pts, pts);
159 vtkIdType id_surf_node = -1;
160 if (pts[3] == nodes[i_nodes]) id_surf_node = pts[0];
161 if (pts[4] == nodes[i_nodes]) id_surf_node = pts[1];
162 if (pts[5] == nodes[i_nodes]) id_surf_node = pts[2];
163 if (id_surf_node != -1) {
164 vec3_t x0,x1,x2;
165 grid->GetPoint(pts[0],x0.data());
166 grid->GetPoint(pts[1],x1.data());
167 grid->GetPoint(pts[2],x2.data());
168 vec3_t a = x1-x0;
169 vec3_t b = x2-x0;
170 vec3_t c = b-a;
171 double L = (a.abs()+b.abs()+c.abs())/3.0;
172 vec3_t n = b.cross(a);
173 n.normalise();
174 vec3_t x_old;
175 grid->GetPoint(nodes[i_nodes],x_old.data());
176 vec3_t x_new = x_old + Dx - x0;
177 if (n*x_new <= 0) {
178 x_new -= (x_new*n)*n;
179 x_new += 1e-4*L*n;
180 x_new += x0;
181 Dx = x_new - x_old;
182 };
183 };
184 };
185 };
186 };
187 };
188 };
189
190 bool GridSmoother::moveNode(int i_nodes, vec3_t &Dx)
191 {
192 vtkIdType id_node = nodes[i_nodes];
193 vec3_t x_old;
194 grid->GetPoint(id_node, x_old.data());
195 bool moved = false;
196 for (int i_relaxation = 0; i_relaxation < N_relaxations; ++i_relaxation) {
197 if (setNewPosition(id_node, x_old + Dx)) {
198 moved = true;
199 break;
200 };
201 Dx *= 0.5;
202 };
203 return moved;
204 };
205
206 double GridSmoother::errThickness(double x)
207 {
208 /*
209 double X[5], Y[5];
210 X[0] = 0.0; Y[0] = 10000;
211 X[1] = 0.1; Y[1] = 1.0;
212 X[2] = 1.0; Y[2] = 0.0;
213 X[3] = 2.0; Y[3] = 1.0;
214 X[4] = 3.0; Y[4] = 2.0;
215 int i = 0;
216 while ((i < 3) && (x > X[i+1])) ++i;
217 double err = Y[i] + (x-X[i])*(Y[i+1]-Y[i])/(X[i+1]-X[i]);
218 if (err < 0) {
219 cout << x << ',' << err <<endl;
220 cout << i << endl;
221 cout << (x-X[i]) << endl;
222 cout << (X[i+1]-X[i]) << endl;
223 cout << (Y[i+1]-Y[i]) << endl;
224 EG_BUG;
225 };
226 */
227
228 if (x > 1) x = 2 - x;
229 double delta = 0.01;
230 double a = 5.0;
231 double err = 0;
232 if (x < 0) err = 1.0/delta - 1.0/(a+delta);
233 else if (x < 1) err = 1/(a*x+delta) - 1.0/(a+delta);
234
235 return err;
236 };
237
238 double GridSmoother::func(vec3_t x)
239 {
240 EG_VTKDCC(vtkDoubleArray, err_tet, grid, "cell_err_tet");
241 EG_VTKDCC(vtkDoubleArray, err_pria, grid, "cell_err_pria");
242 EG_VTKDCC(vtkDoubleArray, err_prib, grid, "cell_err_prib");
243 EG_VTKDCC(vtkDoubleArray, err_pric, grid, "cell_err_pric");
244 EG_VTKDCC(vtkDoubleArray, err_prid, grid, "cell_err_prid");
245 EG_VTKDCC(vtkDoubleArray, err_prie, grid, "cell_err_prie");
246 EG_VTKDCC(vtkDoubleArray, err_prif, grid, "cell_err_prif");
247
248 vec3_t x_old;
249 grid->GetPoint(nodes[i_nodes_opt], x_old.data());
250 grid->GetPoints()->SetPoint(nodes[i_nodes_opt], x.data());
251 double f = 0;
252 double f13 = 1.0/3.0;
253 double f14 = 0.25;
254
255 vec3_t n_node(1,0,0);
256 QList<vec3_t> n_pri;
257
258 foreach (int i_cells, n2c[i_nodes_opt]) {
259 vtkIdType id_cell = cells[i_cells];
260 err_tet->SetValue(id_cell,0);
261 err_pria->SetValue(id_cell,0);
262 err_prib->SetValue(id_cell,0);
263 err_pric->SetValue(id_cell,0);
264 err_prid->SetValue(id_cell,0);
265 if (isVolume(grid, id_cell)) {
266 vtkIdType type_cell = grid->GetCellType(id_cell);
267 vtkIdType N_pts, *pts;
268 grid->GetCellPoints(id_cell, N_pts, pts);
269 QVector<vec3_t> xn(N_pts);
270 for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
271 grid->GetPoint(pts[i_pts],xn[i_pts].data());
272 };
273 if (type_cell == VTK_TETRA) {
274 double L = 0;
275 L += (xn[0]-xn[1]).abs();
276 L += (xn[0]-xn[2]).abs();
277 L += (xn[0]-xn[3]).abs();
278 L += (xn[1]-xn[2]).abs();
279 L += (xn[1]-xn[3]).abs();
280 L += (xn[2]-xn[3]).abs();
281 L /= 6;
282 double V1 = GeometryTools::cellVA(grid, id_cell, true);
283 double V2 = sqrt(1.0/72.0)*L*L*L;
284 double e = sqr((V1-V2)/V2);
285 f += w_tet*e;
286 err_tet->SetValue(id_cell,e);
287 };
288 if (type_cell == VTK_WEDGE) {
289 double L = 0;
290 L += (xn[0]-xn[1]).abs();
291 L += (xn[0]-xn[2]).abs();
292 L += (xn[1]-xn[2]).abs();
293 L *= H/3.0;
294 vec3_t a = xn[2]-xn[0];
295 vec3_t b = xn[1]-xn[0];
296 vec3_t c = xn[5]-xn[3];
297 vec3_t d = xn[4]-xn[3];
298 vec3_t n_face[5];
299 vec3_t x_face[5];
300 n_face[0] = GeometryTools::triNormal(xn[0],xn[1],xn[2]);
301 n_face[1] = GeometryTools::triNormal(xn[3],xn[5],xn[4]);
302 n_face[2] = GeometryTools::quadNormal(xn[0],xn[3],xn[4],xn[1]);
303 n_face[3] = GeometryTools::quadNormal(xn[1],xn[4],xn[5],xn[2]);
304 n_face[4] = GeometryTools::quadNormal(xn[0],xn[2],xn[5],xn[3]);
305 x_face[0] = f13*(xn[0]+xn[1]+xn[2]);
306 x_face[1] = f13*(xn[3]+xn[4]+xn[5]);
307 x_face[2] = f14*(xn[0]+xn[3]+xn[4]+xn[1]);
308 x_face[3] = f14*(xn[1]+xn[4]+xn[5]+xn[2]);
309 x_face[4] = f14*(xn[0]+xn[2]+xn[5]+xn[3]);
310
311 double A1 = 0.5*n_face[0].abs();
312 double A2 = 0.5*n_face[1].abs();
313 for (int i_face = 0; i_face < 5; ++i_face) {
314 n_face[i_face].normalise();
315 };
316 if (nodes[i_nodes_opt] == pts[3]) {
317 n_node = xn[3]-xn[0];
318 n_pri.append(n_face[1]);
319 };
320 if (nodes[i_nodes_opt] == pts[4]) {
321 n_node = xn[4]-xn[1];
322 n_pri.append(n_face[1]);
323 };
324 if (nodes[i_nodes_opt] == pts[5]) {
325 n_node = xn[5]-xn[2];
326 n_pri.append(n_face[1]);
327 };
328 vec3_t v0 = xn[0]-xn[3];
329 vec3_t v1 = xn[1]-xn[4];
330 vec3_t v2 = xn[2]-xn[5];
331
332 double h0 = v0*n_face[0];
333 double h1 = v1*n_face[0];
334 double h2 = v2*n_face[0];
335 if (h0 > 0.5*L) h0 = max(v0.abs(),h0);
336 if (h1 > 0.5*L) h1 = max(v1.abs(),h1);
337 if (h2 > 0.5*L) h2 = max(v2.abs(),h2);
338
339 {
340
341 double e1 = errThickness(h0/L);
342 double e2 = errThickness(h1/L);
343 double e3 = errThickness(h2/L);
344 double e = max(e1,max(e2,e3));
345 //f += w_h*f13*e1;
346 //f += w_h*f13*e2;
347 //f += w_h*f13*e3;
348 //err_pria->SetValue(id_cell,f13*(e1+e2+e3));
349
350 f += w_h*e;
351 err_pria->SetValue(id_cell,0);
352 };
353 if ((h0 > 0.01*L) && (h1 > 0.01*L) && (h2 > 0.01*L)) {
354 v0.normalise();
355 v1.normalise();
356 v2.normalise();
357 double e1 = f13*(1-v0*v1);
358 double e2 = f13*(1-v0*v2);
359 double e3 = f13*(1-v1*v2);
360 f += w_par*e1;
361 f += w_par*e2;
362 f += w_par*e3;
363 err_pric->SetValue(id_cell,f13*(e1+e2+e3));
364 };
365 if ((h0 > 0.01*L) && (h1 > 0.01*L) && (h2 > 0.01*L)) {
366 double e = (1+n_face[0]*n_face[1]);
367 f += w_n*e;
368 err_prib->SetValue(id_cell,e);
369 };
370 if ((h0 > 0.01*L) && (h1 > 0.01*L) && (h2 > 0.01*L)) {
371 double e = sqr((A1-A2)/(A1+A2));
372 f += w_A*e;
373 err_prid->SetValue(id_cell,e);
374 };
375 if ((h0 > 0.01*L) && (h1 > 0.01*L) && (h2 > 0.01*L)) {
376 double e_skew = 0;
377 double e_orth = 0;
378 int N = 0;
379 vec3_t xc = cellCentre(grid, id_cell);
380 for (int i_face = 0; i_face < 5; ++i_face) {
381 int i_cells_neigh = c2c[i_cells][i_face];
382 if (i_cells_neigh != -1) {
383 vtkIdType id_neigh_cell = cells[i_cells_neigh];
384 if (isVolume(grid, id_neigh_cell)) {
385 vec3_t vc = cellCentre(grid, id_neigh_cell) - xc;
386 vec3_t vf = x_face[i_face] - xc;
387 vc.normalise();
388 vf.normalise();
389 e_skew += (1-vc*vf);
390 e_orth += (1-vc*n_face[i_face]);
391 ++N;
392 };
393 };
394 };
395 e_skew /= N;
396 e_orth /= N;
397 f += w_skew*e_skew + w_orth*e_orth;
398 err_prie->SetValue(id_cell,e_skew);
399 err_prif->SetValue(id_cell,e_orth);
400 };
401
402 double f_sharp2 = 0;
403 for (int j = 2; j <= 4; ++j) {
404 vtkIdType id_ncell = c2c[id_cell][j];
405 if (id_ncell != -1) {
406 if (grid->GetCellType(id_ncell) == VTK_WEDGE) {
407 vtkIdType N_pts, *pts;
408 grid->GetCellPoints(id_ncell, N_pts, pts);
409 QVector<vec3_t> x(3);
410 for (int i_pts = 3; i_pts <= 5; ++i_pts) {
411 grid->GetPoint(pts[i_pts],x[i_pts-3].data());
412 };
413 vec3_t n = GeometryTools::triNormal(x[0],x[2],x[1]);
414 n.normalise();
415 f_sharp2 += pow(fabs(1-n_face[1]*n), e_sharp2);
416 };
417 };
418 };
419 f += w_sharp2*f_sharp2;
420 };
421 };
422 };
423 grid->GetPoints()->SetPoint(nodes[i_nodes_opt], x_old.data());
424 n_node.normalise();
425 {
426 double f_sharp1 = 0;
427 foreach (vec3_t n, n_pri) {
428 f_sharp1 += pow(fabs(1-n_node*n), e_sharp1);
429 };
430 f += w_sharp1*f_sharp1;
431 };
432 return f;
433 };
434
435 void GridSmoother::resetStencil()
436 {
437 stencil.clear();
438 };
439
440 void GridSmoother::addToStencil(double C, vec3_t x)
441 {
442 stencil_node_t sn;
443 sn.x = x;
444 sn.C = C;
445 stencil.append(sn);
446 };
447
448 void GridSmoother::operate()
449 {
450 markNodes();
451
452 EG_VTKDCC(vtkIntArray, bc, grid, "cell_code");
453 EG_VTKDCN(vtkIntArray, node_status, grid, "node_status");
454 EG_VTKDCN(vtkIntArray, node_layer, grid, "node_layer");
455
456 QVector<QSet<int> > n2bc(nodes.size());
457 QVector<bool> prism_node(nodes.size(),false);
458 foreach (vtkIdType id_cell, cells) {
459 if (isSurface(grid, id_cell)) {
460 vtkIdType N_pts, *pts;
461 grid->GetCellPoints(id_cell, N_pts, pts);
462 for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
463 n2bc[_nodes[pts[i_pts]]].insert(bc->GetValue(id_cell));
464 };
465 };
466 if (grid->GetCellType(id_cell) == VTK_WEDGE) {
467 vtkIdType N_pts, *pts;
468 grid->GetCellPoints(id_cell, N_pts, pts);
469 for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
470 if (_nodes[pts[i_pts]] != -1) {
471 prism_node[_nodes[pts[i_pts]]] = true;
472 };
473 };
474 };
475 };
476
477 F_old = 0;
478 setPrismWeighting();
479 for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
480 if (prism_node[i_nodes]) {
481 vec3_t x;
482 i_nodes_opt = i_nodes;
483 grid->GetPoint(nodes[i_nodes], x.data());
484 F_old += func(x);
485 };
486 };
487 setAllWeighting();
488
489 cout << "\nsmoothing volume mesh (" << N_marked_nodes << " nodes)" << endl;
490 for (int i_iterations = 0; i_iterations < N_iterations; ++i_iterations) {
491 cout << "iteration " << i_iterations+1 << "/" << N_iterations << endl;
492 int N_blocked = 0;
493 int N_searched = 0;
494 int N_illegal = 0;
495 int N1 = 0;
496 int N2 = 0;
497 int N3 = 0;
498 QTime start = QTime::currentTime();
499 for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
500 dbg = false;
501 bool smooth = node_marked[nodes[i_nodes]];
502 if (smooth) {
503 foreach (int bc, n2bc[i_nodes]) {
504 if (boundary_codes.contains(bc) || (boundary_codes.size() ==0)) {
505 smooth = false;
506 };
507 };
508 foreach (int i_cells, n2c[i_nodes]) {
509 vtkIdType type_cell = grid->GetCellType(cells[i_cells]);
510 if ((type_cell == VTK_WEDGE) && !smooth_prisms) {
511 smooth = false;
512 };
513 };
514 };
515 if (smooth) {
516 vtkIdType id_node = nodes[i_nodes];
517 vec3_t xn;
518 vec3_t x_old;
519 grid->GetPoint(id_node, x_old.data());
520 resetStencil();
521 bool is_surf = n2bc[i_nodes].size() > 0;
522 int N = 0;
523 foreach (int j_nodes, n2n[i_nodes]) {
524 if (!is_surf || (n2bc[j_nodes].size() > 0)) {
525 vtkIdType id_neigh_node = nodes[j_nodes];
526 grid->GetPoint(id_neigh_node, xn.data());
527 addToStencil(1.0, xn);
528 ++N;
529 };
530 };
531 if (N == 0) {
532 EG_BUG;
533 };
534 vec3_t x_new1 = vec3_t(0,0,0);
535 sum_C = 0;
536 L0 = 0;
537 double L_min = 1e99;
538 foreach (stencil_node_t sn, stencil) {
539 sum_C += sn.C;
540 x_new1 += sn.C*sn.x;
541 double L = (sn.x - x_old).abs();
542 L0 += sn.C*L;
543 L_min = min(L_min, L);
544 };
545 L0 /= sum_C;
546 x_new1 *= 1.0/sum_C;
547 vec3_t Dx1 = x_new1 - x_old;
548 setDeltas(1e-3*L0);
549 i_nodes_opt = i_nodes;
550 vec3_t Dx2(0,0,0);
551 Dx2 = optimise(x_old);
552 vec3_t Dx3 = (-10e-4/func(x_old))*grad_f;
553 correctDx(i_nodes, Dx1);
554 correctDx(i_nodes, Dx2);
555 correctDx(i_nodes, Dx3);
556 vec3_t Dx = Dx1;
557 double f = func(x_old + Dx1);
558 int _N1 = 1;
559 int _N2 = 0;
560 int _N3 = 0;
561 if (f > func(x_old + Dx2)) {
562 Dx = Dx2;
563 f = func(x_old + Dx2);
564 _N1 = 0;
565 _N2 = 1;
566 _N3 = 0;
567 };
568 if (f > func(x_old + Dx3)) {
569 Dx = Dx3;
570 _N1 = 0;
571 _N2 = 0;
572 _N3 = 1;
573 };
574 if (!moveNode(i_nodes, Dx)) {
575 // search for a better place
576 vec3_t x_save = x_old;
577 vec3_t ex(L_search*L0/N_search,0,0);
578 vec3_t ey(0,L_search*L0/N_search,0);
579 vec3_t ez(0,0,L_search*L0/N_search);
580 vec3_t x_best = x_old;
581 double f_min = func(x_old);
582 bool found = false;
583 bool illegal = false;
584 if (!setNewPosition(id_node,x_old)) {
585 illegal = true;
586 };
587 for (int i = -N_search; i <= N_search; ++i) {
588 for (int j = -N_search; j <= N_search; ++j) {
589 for (int k = -N_search; k <= N_search; ++k) {
590 if ((i != 0) || (j != 0) || (k != 0)) {
591 vec3_t Dx = double(i)*ex + double(j)*ey + double(k)*ez;
592 correctDx(i_nodes, Dx);
593 vec3_t x = x_old + x;
594 if (setNewPosition(id_node, x)) {
595 grid->GetPoints()->SetPoint(id_node, x_old.data());
596 double f = func(x);
597 if (f < f_min) {
598 f_min = f;
599 x_best = x;
600 found = true;
601 illegal = false;
602 };
603 };
604 };
605 };
606 };
607 };
608 if (found) {
609 grid->GetPoints()->SetPoint(id_node, x_best.data());
610 ++N_searched;
611 } else {
612 ++N_blocked;
613 if (illegal) {
614 ++N_illegal;
615 };
616 };
617 } else {
618 N1 += _N1;
619 N2 += _N2;
620 N3 += _N3;
621 };
622 };
623 };
624
625 cout << N1 << " type 1 movements (simple)" << endl;
626 cout << N2 << " type 2 movements (Newton)" << endl;
627 cout << N3 << " type 3 movements (gradient)" << endl;
628 //cout << N_blocked << " movements blocked" << endl;
629 //cout << N_searched << " movements by search" << endl;
630 //cout << N_illegal << " nodes in illegal positions" << endl;
631
632 cout << start.secsTo(QTime::currentTime()) << " seconds elapsed" << endl;
633 F_new = 0;
634 setPrismWeighting();
635 for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
636 if (prism_node[i_nodes]) {
637 vec3_t x;
638 i_nodes_opt = i_nodes;
639 grid->GetPoint(nodes[i_nodes], x.data());
640 F_new += func(x);
641 };
642 };
643 setAllWeighting();
644 cout << "total prism error (old) = " << F_old << endl;
645 cout << "total prism error (new) = " << F_new << endl;
646 double f_old = max(1e-10,F_old);
647 cout << "total prism improvement = " << 100*(1-F_new/f_old) << "%" << endl;
648 };
649 cout << "done" << endl;
650 };
651
652 double GridSmoother::improvement()
653 {
654 double f_old = max(1e-10,F_old);
655 return 1-F_new/f_old;
656 };
657
Clone of the official engrid repository