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