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 "gridsmoother.h"
24 #include "guimainwindow.h"
30 GridSmoother::GridSmoother()
34 N_boundary_corrections
= 20;
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
);
58 void GridSmoother::markNodes()
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
) {
72 for (int i_pts
= 0; i_pts
< N_pts
; ++i_pts
) {
73 if (node_marked
[pts
[i_pts
]]) {
79 for (int i_pts
= 0; i_pts
< N_pts
; ++i_pts
) {
80 new_mark
[pts
[i_pts
]] = true;
84 qCopy(new_mark
.begin(),new_mark
.end(),node_marked
.begin());
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
]) {
96 bool GridSmoother::setNewPosition(vtkIdType id_node
, vec3_t x_new
)
98 using namespace GeometryTools
;
101 grid
->GetPoint(id_node
, x_old
.data());
102 grid
->GetPoints()->SetPoint(id_node
, x_new
.data());
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) {
112 //if (dbg) cout << id_node << " : tetra negative" << endl;
115 if (type_cell
== VTK_WEDGE
) {
116 vtkIdType N_pts
, *pts
;
118 grid
->GetCellPoints(id_cell
, N_pts
, pts
);
120 for (int i
= 0; i
< 4; ++i
) { // variation
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());
126 if (GeometryTools::tetraVol(xtet
[0], xtet
[1], xtet
[2], xtet
[3]) < 0) {
128 //if (dbg) cout << id_node << " : prism negative" << endl;
141 grid
->GetPoints()->SetPoint(id_node
, x_old
.data());
146 void GridSmoother::correctDx(int i_nodes
, vec3_t
&Dx
)
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
);
153 vec3_t n
= GeometryTools::cellNormal(grid
, id_cell
);
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) {
167 grid
->GetPoint(pts
[0],x0
.data());
168 grid
->GetPoint(pts
[1],x1
.data());
169 grid
->GetPoint(pts
[2],x2
.data());
173 double L
= (a
.abs()+b
.abs()+c
.abs())/3.0;
174 vec3_t n
= b
.cross(a
);
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
;
192 bool GridSmoother::moveNode(int i_nodes
, vec3_t
&Dx
)
194 vtkIdType id_node
= nodes
[i_nodes
];
196 grid
->GetPoint(id_node
, x_old
.data());
198 for (int i_relaxation
= 0; i_relaxation
< N_relaxations
; ++i_relaxation
) {
199 if (setNewPosition(id_node
, x_old
+ Dx
)) {
208 double GridSmoother::errThickness(double x
)
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;
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]);
221 cout << x << ',' << err <<endl;
223 cout << (x-X[i]) << endl;
224 cout << (X[i+1]-X[i]) << endl;
225 cout << (Y[i+1]-Y[i]) << endl;
230 if (x
> 1) x
= 2 - x
;
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
);
240 double GridSmoother::func(vec3_t x
)
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");
251 grid
->GetPoint(nodes
[i_nodes_opt
], x_old
.data());
252 grid
->GetPoints()->SetPoint(nodes
[i_nodes_opt
], x
.data());
254 double f13
= 1.0/3.0;
257 vec3_t
n_node(1,0,0);
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());
275 if (type_cell
== VTK_TETRA
) {
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();
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
);
288 err_tet
->SetValue(id_cell
,e
);
290 if (type_cell
== VTK_WEDGE
) {
292 L
+= (xn
[0]-xn
[1]).abs();
293 L
+= (xn
[0]-xn
[2]).abs();
294 L
+= (xn
[1]-xn
[2]).abs();
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];
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]);
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();
318 if (nodes
[i_nodes_opt
] == pts
[3]) {
319 n_node
= xn
[3]-xn
[0];
320 n_pri
.append(n_face
[1]);
322 if (nodes
[i_nodes_opt
] == pts
[4]) {
323 n_node
= xn
[4]-xn
[1];
324 n_pri
.append(n_face
[1]);
326 if (nodes
[i_nodes_opt
] == pts
[5]) {
327 n_node
= xn
[5]-xn
[2];
328 n_pri
.append(n_face
[1]);
330 vec3_t v0
= xn
[0]-xn
[3];
331 vec3_t v1
= xn
[1]-xn
[4];
332 vec3_t v2
= xn
[2]-xn
[5];
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
);
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
));
350 //err_pria->SetValue(id_cell,f13*(e1+e2+e3));
353 err_pria
->SetValue(id_cell
,0);
355 if ((h0
> 0.01*L
) && (h1
> 0.01*L
) && (h2
> 0.01*L
)) {
359 double e1
= f13
*(1-v0
*v1
);
360 double e2
= f13
*(1-v0
*v2
);
361 double e3
= f13
*(1-v1
*v2
);
365 err_pric
->SetValue(id_cell
,f13
*(e1
+e2
+e3
));
367 if ((h0
> 0.01*L
) && (h1
> 0.01*L
) && (h2
> 0.01*L
)) {
368 double e
= (1+n_face
[0]*n_face
[1]);
370 err_prib
->SetValue(id_cell
,e
);
372 if ((h0
> 0.01*L
) && (h1
> 0.01*L
) && (h2
> 0.01*L
)) {
373 double e
= sqr((A1
-A2
)/(A1
+A2
));
375 err_prid
->SetValue(id_cell
,e
);
377 if ((h0
> 0.01*L
) && (h1
> 0.01*L
) && (h2
> 0.01*L
)) {
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
;
392 e_orth
+= (1-vc
*n_face
[i_face
]);
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
);
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());
415 vec3_t n
= GeometryTools::triNormal(x
[0],x
[2],x
[1]);
417 f_sharp2
+= pow(fabs(1-n_face
[1]*n
), e_sharp2
);
421 f
+= w_sharp2
*f_sharp2
;
425 grid
->GetPoints()->SetPoint(nodes
[i_nodes_opt
], x_old
.data());
429 foreach (vec3_t n
, n_pri
) {
430 f_sharp1
+= pow(fabs(1-n_node
*n
), e_sharp1
);
432 f
+= w_sharp1
*f_sharp1
;
437 void GridSmoother::resetStencil()
442 void GridSmoother::addToStencil(double C
, vec3_t x
)
450 void GridSmoother::operate()
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");
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
));
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;
482 for (int i_nodes
= 0; i_nodes
< nodes
.size(); ++i_nodes
) {
483 if (prism_node
[i_nodes
]) {
485 i_nodes_opt
= i_nodes
;
486 grid
->GetPoint(nodes
[i_nodes
], x
.data());
489 F_max_old
= max(F_max_old
,f
);
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
;
503 QTime start
= QTime::currentTime();
504 for (int i_nodes
= 0; i_nodes
< nodes
.size(); ++i_nodes
) {
506 bool smooth
= node_marked
[nodes
[i_nodes
]];
508 foreach (int bc
, n2bc
[i_nodes
]) {
509 if (boundary_codes
.contains(bc
) || (boundary_codes
.size() ==0)) {
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
) {
521 vtkIdType id_node
= nodes
[i_nodes
];
524 grid
->GetPoint(id_node
, x_old
.data());
526 bool is_surf
= n2bc
[i_nodes
].size() > 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
);
539 vec3_t x_new1
= vec3_t(0,0,0);
543 foreach (stencil_node_t sn
, stencil
) {
546 double L
= (sn
.x
- x_old
).abs();
548 L_min
= min(L_min
, L
);
552 vec3_t Dx1
= x_new1
- x_old
;
554 i_nodes_opt
= i_nodes
;
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
);
562 double f
= func(x_old
+ Dx1
);
566 if (f
> func(x_old
+ Dx2
)) {
568 f
= func(x_old
+ Dx2
);
573 if (f
> func(x_old
+ Dx3
)) {
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
);
588 bool illegal
= false;
589 if (!setNewPosition(id_node
,x_old
)) {
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());
614 grid
->GetPoints()->SetPoint(id_node
, x_best
.data());
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;
637 cout
<< start
.secsTo(QTime::currentTime()) << " seconds elapsed" << endl
;
641 for (int i_nodes
= 0; i_nodes
< nodes
.size(); ++i_nodes
) {
642 if (prism_node
[i_nodes
]) {
644 i_nodes_opt
= i_nodes
;
645 grid
->GetPoint(nodes
[i_nodes
], x
.data());
648 F_max_new
= max(F_max_new
,f
);
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
;
659 cout
<< "done" << endl
;
662 double GridSmoother::improvement()
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
;