?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
fixed inside sharp feature problem
[engrid.git] / src / gridsmoother.cpp
blobbcf302384853fc8699658eb9f890b10153fd3c1e
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 m_NumIterations = 5;
33 m_NumRelaxations = 5;
34 m_ReductionFactor = 0.2;
35 m_NumBoundaryCorrections = 20;
36 m_NumSearch = 10;
37 m_LSearch = 0.5;
38 m_SmoothPrisms = true;
39 m_FOld = 0;
40 m_FNew = 0;
41 m_H = 1.5;
42
43 getSet("boundary layer", "tetra weighting", 0.1, m_WTet);
44 getSet("boundary layer", "tetra exponent", 1.2, m_ETet);
45 getSet("boundary layer", "layer height weighting", 1.0, m_WH);
46 getSet("boundary layer", "parallel edges weighting", 3.0, m_WPar);
47 getSet("boundary layer", "parallel faces weighting", 5.0, m_WN);
48 getSet("boundary layer", "similar face area weighting", 5.0, m_WA);
49 getSet("boundary layer", "skewness weighting", 0.0, m_WSkew);
50 getSet("boundary layer", "orthogonality weighting", 0.0, m_WOrth);
51 getSet("boundary layer", "sharp features on nodes weighting", 8.0, m_WSharp1);
52 getSet("boundary layer", "sharp features on nodes exponent", 1.4, m_ESharp1);
53 getSet("boundary layer", "sharp features on edges weighting", 3.0, m_WSharp2);
54 getSet("boundary layer", "sharp features on edges exponent", 1.3, m_ESharp2);
55 getSet("boundary layer", "number of smoothing sub-iterations", 5, m_NumIterations);
56 getSet("boundary layer", "angle for sharp features", 45.00, m_CritAngle);
57 getSet("boundary layer", "post smoothing strength", 0.1, m_PostSmoothingStrength);
58
59 getSet("boundary layer", "use strict prism checking", false, m_StrictPrismChecking);
60
61 m_CritAngle = GeometryTools::deg2rad(m_CritAngle);
62 m_SimpleOperation = false;
63 m_PostOperation = false;
64
65 }
66
67 void GridSmoother::markNodes()
68 {
69 m_NodeMarked.fill(false,m_Grid->GetNumberOfPoints());
70 QVector<bool> new_mark(m_Grid->GetNumberOfPoints());
71 for (int i_iterations = 0; i_iterations < 4; ++i_iterations) {
72 qCopy(m_NodeMarked.begin(),m_NodeMarked.end(),new_mark.begin());
73 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
74 bool mark_cell = false;
75 vtkIdType type_cell, N_pts, *pts;
76 type_cell = m_Grid->GetCellType(id_cell);
77 m_Grid->GetCellPoints(id_cell, N_pts, pts);
78 if (type_cell == VTK_WEDGE) {
79 mark_cell = true;
80 } else {
81 for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
82 if (m_NodeMarked[pts[i_pts]]) {
83 mark_cell = true;
84 }
85 }
86 }
87 if (mark_cell) {
88 for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
89 new_mark[pts[i_pts]] = true;
90 }
91 }
92 }
93 qCopy(new_mark.begin(),new_mark.end(),m_NodeMarked.begin());
94 }
95 m_NumMarkedNodes = 0;
96 QVector<vtkIdType> nodes = m_Part.getNodes();
97 foreach (vtkIdType id_node, nodes) {
98 if (id_node < 0) EG_BUG;
99 if (id_node > m_Grid->GetNumberOfPoints()) EG_BUG;
100 if (m_NodeMarked[id_node]) {
101 ++m_NumMarkedNodes;
102 }
103 }
104 }
105
106 void GridSmoother::findCriticalTetras()
107 {
108 m_CriticalTetra.fill(false, m_Grid->GetNumberOfCells());
109 QVector<bool> pure_tetra_node(m_Grid->GetNumberOfPoints(), true);
110 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
111 if (isVolume(id_cell, m_Grid)) {
112 vtkIdType type_cell = m_Grid->GetCellType(id_cell);
113 if (type_cell != VTK_TETRA) {
114 vtkIdType N_pts, *pts;
115 m_Grid->GetCellPoints(id_cell, N_pts, pts);
116 for (int i = 0; i < N_pts; ++i) {
117 pure_tetra_node[pts[i]] = false;
118 }
119 }
120 }
121 }
122 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
123 vtkIdType type_cell = m_Grid->GetCellType(id_cell);
124 if (type_cell == VTK_TETRA) {
125 vtkIdType N_pts, *pts;
126 bool crit = true;
127 m_Grid->GetCellPoints(id_cell, N_pts, pts);
128 for (int i = 0; i < N_pts; ++i) {
129 if (pure_tetra_node[pts[i]]) {
130 crit = false;
131 break;
132 }
133 }
134 m_CriticalTetra[id_cell] = crit;
135 }
136 }
137 }
138
139 bool GridSmoother::setNewPosition(vtkIdType id_node, vec3_t x_new)
140 {
141 using namespace GeometryTools;
142
143 vec3_t x_old;
144 m_Grid->GetPoint(id_node, x_old.data());
145 m_Grid->GetPoints()->SetPoint(id_node, x_new.data());
146 bool move = true;
147 Elements E;
148
149 l2g_t cells = m_Part.getCells();
150 l2l_t n2c = m_Part.getN2C();
151
152 foreach (int i_cells, n2c[id_node]) {
153 vtkIdType id_cell = cells[i_cells];
154 vtkIdType type_cell = m_Grid->GetCellType(id_cell);
155 if (type_cell == VTK_TETRA) {
156 if (GeometryTools::cellVA(m_Grid, id_cell) < 0) {
157 move = false;
158 }
159 }
160
161 if (type_cell == VTK_WEDGE && m_StrictPrismChecking) {
162 vtkIdType N_pts, *pts;
163 vec3_t xtet[4];
164 m_Grid->GetCellPoints(id_cell, N_pts, pts);
165 bool ok = true;
166 for (int i = 0; i < 4; ++i) { // variation
167 ok = true;
168 for (int j = 0; j < 3; ++j) { // tetrahedron
169 for (int k = 0; k < 4; ++k) { // node
170 m_Grid->GetPoint(pts[E.priTet(i,j,k)], xtet[k].data());
171 }
172 double V = GeometryTools::tetraVol(xtet[0], xtet[1], xtet[2], xtet[3]);
173 if (V <= 0) {
174 ok = false;
175 }
176 }
177 if (ok) {
178 break;
179 }
180 }
181 if (!ok) {
182 move = false;
183 }
184 }
185 }
186
187
188 if (!move) {
189 m_Grid->GetPoints()->SetPoint(id_node, x_old.data());
190 }
191 return move;
192 }
193
194 void GridSmoother::correctDx(int i_nodes, vec3_t &Dx)
195 {
196 l2g_t nodes = m_Part.getNodes();
197 l2l_t n2c = m_Part.getN2C();
198 for (int i_boundary_correction = 0; i_boundary_correction < m_NumBoundaryCorrections; ++i_boundary_correction) {
199 foreach (vtkIdType id_cell, n2c[i_nodes]) {
200 if (isSurface(id_cell, m_Grid)) {
201 double A = GeometryTools::cellVA(m_Grid, id_cell);
202 if (A > 1e-20) {
203 vec3_t n = GeometryTools::cellNormal(m_Grid, id_cell);
204 n.normalise();
205 Dx -= (n*Dx)*n;
206 }
207 } else {
208 if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
209 vtkIdType N_pts, *pts;
210 m_Grid->GetCellPoints(id_cell, N_pts, pts);
211 vtkIdType id_surf_node = -1;
212 if (pts[3] == nodes[i_nodes]) id_surf_node = pts[0];
213 if (pts[4] == nodes[i_nodes]) id_surf_node = pts[1];
214 if (pts[5] == nodes[i_nodes]) id_surf_node = pts[2];
215 if (id_surf_node != -1) {
216 vec3_t x0,x1,x2;
217 m_Grid->GetPoint(pts[0],x0.data());
218 m_Grid->GetPoint(pts[1],x1.data());
219 m_Grid->GetPoint(pts[2],x2.data());
220 vec3_t a = x1-x0;
221 vec3_t b = x2-x0;
222 vec3_t c = b-a;
223 double L = (a.abs()+b.abs()+c.abs())/3.0;
224 vec3_t n = b.cross(a);
225 n.normalise();
226 vec3_t x_old;
227 m_Grid->GetPoint(nodes[i_nodes],x_old.data());
228 vec3_t x_new = x_old + Dx - x0;
229 if ( (n*x_new) <= 0 ) {
230 x_new -= (x_new*n)*n;
231 x_new += 1e-4*L*n;
232 x_new += x0;
233 Dx = x_new - x_old;
234 }
235 }
236 }
237 }
238 }
239 }
240 }
241
242 bool GridSmoother::moveNode(int i_nodes, vec3_t &Dx)
243 {
244 l2g_t nodes = m_Part.getNodes();
245 vtkIdType id_node = nodes[i_nodes];
246 vec3_t x_old;
247 m_Grid->GetPoint(id_node, x_old.data());
248 bool moved = false;
249 for (int i_relaxation = 0; i_relaxation < m_NumRelaxations; ++i_relaxation) {
250 if (setNewPosition(id_node, x_old + Dx)) {
251 moved = true;
252 break;
253 }
254 Dx *= 0.5;
255 }
256 return moved;
257 }
258
259 double GridSmoother::errThickness(double x)
260 {
261 if (x > 1) x = 2 - x;
262 double delta = 0.01;
263 double a = 5.0;
264 double err = 0;
265 if (x < 0) err = 1.0/delta - 1.0/(a+delta);
266 else if (x < 1) err = 1/(a*x+delta) - 1.0/(a+delta);
267 return err;
268 }
269
270 double GridSmoother::func(vec3_t x)
271 {
272 l2g_t nodes = m_Part.getNodes();
273 l2g_t cells = m_Part.getCells();
274 l2l_t n2c = m_Part.getN2C();
275 l2l_t c2c = m_Part.getC2C();
276
277 vec3_t x_old;
278 m_Grid->GetPoint(nodes[m_INodesOpt], x_old.data());
279 m_Grid->GetPoints()->SetPoint(nodes[m_INodesOpt], x.data());
280 double f = 0;
281 double f13 = 1.0/3.0;
282 double f14 = 0.25;
283
284 vec3_t n_node(1,0,0);
285 QList<vec3_t> n_pri;
286
287 foreach (int i_cells, n2c[m_INodesOpt]) {
288 vtkIdType id_cell = cells[i_cells];
289 if (isVolume(id_cell, m_Grid)) {
290 vtkIdType type_cell = m_Grid->GetCellType(id_cell);
291 vtkIdType N_pts, *pts;
292 m_Grid->GetCellPoints(id_cell, N_pts, pts);
293 QVector<vec3_t> xn(N_pts);
294 for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
295 m_Grid->GetPoint(pts[i_pts],xn[i_pts].data());
296 }
297 if (type_cell == VTK_TETRA) {
298 double L_max = 1e-20;
299 L_max = max((xn[0]-xn[1]).abs(), L_max);
300 L_max = max((xn[0]-xn[2]).abs(), L_max);
301 L_max = max((xn[0]-xn[3]).abs(), L_max);
302 L_max = max((xn[1]-xn[2]).abs(), L_max);
303 L_max = max((xn[1]-xn[3]).abs(), L_max);
304 L_max = max((xn[2]-xn[3]).abs(), L_max);
305 double V1 = GeometryTools::cellVA(m_Grid, id_cell, true);
306 double V2 = 0.1178511301977579207*L_max*L_max*L_max;
307 double e = fabs(V1-V2)/V2;
308 f += m_WTet*pow(e, m_ETet);
309 }
310 if (type_cell == VTK_WEDGE) {
311 double L = 0;
312 L += (xn[0]-xn[1]).abs();
313 L += (xn[0]-xn[2]).abs();
314 L += (xn[1]-xn[2]).abs();
315 L *= m_H/3.0;
316 vec3_t a = xn[2]-xn[0];
317 vec3_t b = xn[1]-xn[0];
318 vec3_t c = xn[5]-xn[3];
319 vec3_t d = xn[4]-xn[3];
320 vec3_t n_face[5];
321 vec3_t x_face[5];
322 n_face[0] = GeometryTools::triNormal(xn[0],xn[1],xn[2]);
323 n_face[1] = GeometryTools::triNormal(xn[3],xn[5],xn[4]);
324 n_face[2] = GeometryTools::quadNormal(xn[0],xn[3],xn[4],xn[1]);
325 n_face[3] = GeometryTools::quadNormal(xn[1],xn[4],xn[5],xn[2]);
326 n_face[4] = GeometryTools::quadNormal(xn[0],xn[2],xn[5],xn[3]);
327 x_face[0] = f13*(xn[0]+xn[1]+xn[2]);
328 x_face[1] = f13*(xn[3]+xn[4]+xn[5]);
329 x_face[2] = f14*(xn[0]+xn[3]+xn[4]+xn[1]);
330 x_face[3] = f14*(xn[1]+xn[4]+xn[5]+xn[2]);
331 x_face[4] = f14*(xn[0]+xn[2]+xn[5]+xn[3]);
332
333 double A1 = 0.5*n_face[0].abs();
334 double A2 = 0.5*n_face[1].abs();
335 for (int i_face = 0; i_face < 5; ++i_face) {
336 n_face[i_face].normalise();
337 }
338 m_IdFoot[pts[3]] = pts[0];
339 m_IdFoot[pts[4]] = pts[1];
340 m_IdFoot[pts[5]] = pts[2];
341 if (nodes[m_INodesOpt] == pts[3]) {
342 n_node = xn[3]-xn[0];
343 n_pri.append(n_face[1]);
344 m_L[nodes[m_INodesOpt]] = L;
345 }
346 if (nodes[m_INodesOpt] == pts[4]) {
347 n_node = xn[4]-xn[1];
348 n_pri.append(n_face[1]);
349 m_L[nodes[m_INodesOpt]] = L;
350 }
351 if (nodes[m_INodesOpt] == pts[5]) {
352 n_node = xn[5]-xn[2];
353 n_pri.append(n_face[1]);
354 m_L[nodes[m_INodesOpt]] = L;
355 }
356 vec3_t v0 = xn[0]-xn[3];
357 vec3_t v1 = xn[1]-xn[4];
358 vec3_t v2 = xn[2]-xn[5];
359
360 double h0 = v0*n_face[0];
361 double h1 = v1*n_face[0];
362 double h2 = v2*n_face[0];
363 if (h0 > 0.5*L) h0 = max(v0.abs(),h0);
364 if (h1 > 0.5*L) h1 = max(v1.abs(),h1);
365 if (h2 > 0.5*L) h2 = max(v2.abs(),h2);
366
367 {
368
369 double e1 = errThickness(h0/L);
370 double e2 = errThickness(h1/L);
371 double e3 = errThickness(h2/L);
372 double e = max(e1,max(e2,e3));
373 //f += w_h*f13*e1;
374 //f += w_h*f13*e2;
375 //f += w_h*f13*e3;
376 //err_pria->SetValue(id_cell,f13*(e1+e2+e3));
377
378 f += m_WH*e;
379 }
380 if ((h0 > 0.01*L) && (h1 > 0.01*L) && (h2 > 0.01*L)) {
381 v0.normalise();
382 v1.normalise();
383 v2.normalise();
384 double e1 = f13*(1-v0*v1);
385 double e2 = f13*(1-v0*v2);
386 double e3 = f13*(1-v1*v2);
387 f += m_WPar*e1;
388 f += m_WPar*e2;
389 f += m_WPar*e3;
390 }
391 if ((h0 > 0.01*L) && (h1 > 0.01*L) && (h2 > 0.01*L)) {
392 double e = (1+n_face[0]*n_face[1]);
393 f += m_WN*e;
394 }
395 if ((h0 > 0.01*L) && (h1 > 0.01*L) && (h2 > 0.01*L)) {
396 double e = sqr(0.5*(A1-A2)/A1);
397 f += m_WA*e;
398 }
399 if ((h0 > 0.01*L) && (h1 > 0.01*L) && (h2 > 0.01*L)) {
400 double e_skew = 0;
401 double e_orth = 0;
402 int N = 0;
403 vec3_t xc = cellCentre(m_Grid, id_cell);
404 for (int i_face = 0; i_face < 5; ++i_face) {
405 int i_cells_neigh = c2c[i_cells][i_face];
406 if (i_cells_neigh != -1) {
407 vtkIdType id_neigh_cell = cells[i_cells_neigh];
408 if (isVolume(id_neigh_cell, m_Grid)) {
409 vec3_t vc = cellCentre(m_Grid, id_neigh_cell) - xc;
410 vec3_t vf = x_face[i_face] - xc;
411 vc.normalise();
412 vf.normalise();
413 e_skew += (1-vc*vf);
414 e_orth += (1-vc*n_face[i_face]);
415 ++N;
416 }
417 }
418 }
419 e_skew /= N;
420 e_orth /= N;
421 f += m_WSkew*e_skew + m_WOrth*e_orth;
422 }
423
424 double f_sharp2 = 0;
425 for (int j = 2; j <= 4; ++j) {
426 vtkIdType id_ncell = c2c[id_cell][j];
427 if (id_ncell != -1) {
428 if (m_Grid->GetCellType(id_ncell) == VTK_WEDGE) {
429 vtkIdType N_pts, *pts;
430 m_Grid->GetCellPoints(id_ncell, N_pts, pts);
431 QVector<vec3_t> x(3);
432 for (int i_pts = 3; i_pts <= 5; ++i_pts) {
433 m_Grid->GetPoint(pts[i_pts],x[i_pts-3].data());
434 }
435 vec3_t n = GeometryTools::triNormal(x[0],x[2],x[1]);
436 n.normalise();
437 f_sharp2 += pow(fabs(1-n_face[1]*n), m_ESharp2);
438 }
439 }
440 }
441 f += m_WSharp2*f_sharp2;
442 }
443 }
444 }
445 m_Grid->GetPoints()->SetPoint(nodes[m_INodesOpt], x_old.data());
446 n_node.normalise();
447 {
448 double f_sharp1 = 0;
449 foreach (vec3_t n, n_pri) {
450 f_sharp1 += pow(fabs(1-n_node*n), m_ESharp1);
451 }
452 f += m_WSharp1*f_sharp1;
453 }
454 return f;
455 }
456
457 void GridSmoother::resetStencil()
458 {
459 m_Stencil.clear();
460 }
461
462 void GridSmoother::addToStencil(double C, vec3_t x)
463 {
464 stencil_node_t sn;
465 sn.x = x;
466 sn.C = C;
467 m_Stencil.append(sn);
468 }
469
470 void GridSmoother::operateOptimisation()
471 {
472 markNodes();
473 findCriticalTetras();
474
475 l2g_t nodes = m_Part.getNodes();
476 l2g_t cells = m_Part.getCells();
477 g2l_t _nodes = m_Part.getLocalNodes();
478 l2l_t n2c = m_Part.getN2C();
479 l2l_t n2n = m_Part.getN2N();
480
481 EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
482 EG_VTKDCN(vtkIntArray, node_status, m_Grid, "node_status");
483 EG_VTKDCN(vtkIntArray, node_layer, m_Grid, "node_layer");
484
485 m_IdFoot.fill(-1, m_Grid->GetNumberOfPoints());
486 m_L.fill(0, m_Grid->GetNumberOfPoints());
487
488 QVector<QSet<int> > n2bc(nodes.size());
489 QVector<bool> prism_node(nodes.size(),false);
490 foreach (vtkIdType id_cell, cells) {
491 if (isSurface(id_cell, m_Grid)) {
492 vtkIdType N_pts, *pts;
493 m_Grid->GetCellPoints(id_cell, N_pts, pts);
494 for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
495 n2bc[_nodes[pts[i_pts]]].insert(bc->GetValue(id_cell));
496 }
497 }
498 if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
499 vtkIdType N_pts, *pts;
500 m_Grid->GetCellPoints(id_cell, N_pts, pts);
501 for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
502 if (_nodes[pts[i_pts]] != -1) {
503 prism_node[_nodes[pts[i_pts]]] = true;
504 }
505 }
506 }
507 }
508
509 m_FOld = 0;
510 m_FMaxOld = 0;
511 setPrismWeighting();
512 for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
513 if (prism_node[i_nodes]) {
514 vec3_t x;
515 m_INodesOpt = i_nodes;
516 m_Grid->GetPoint(nodes[i_nodes], x.data());
517 double f = func(x);
518 m_FOld += f;
519 m_FMaxOld = max(m_FMaxOld, f);
520 }
521 }
522 setAllWeighting();
523
524 cout << "\nsmoothing volume mesh (" << m_NumMarkedNodes << " nodes)" << endl;
525 for (int i_iterations = 0; i_iterations < m_NumIterations; ++i_iterations) {
526 cout << "iteration " << i_iterations+1 << "/" << m_NumIterations << endl;
527 int N_blocked = 0;
528 int m_NumSearched = 0;
529 int N_illegal = 0;
530 int N1 = 0;
531 int N2 = 0;
532 int N3 = 0;
533 int N4 = 0;
534 QTime start = QTime::currentTime();
535 for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
536 bool smooth = m_NodeMarked[nodes[i_nodes]];
537 if (smooth) {
538 foreach (int bc, n2bc[i_nodes]) {
539 if (m_BoundaryCodes.contains(bc) || (m_BoundaryCodes.size() == 0)) {
540 smooth = false;
541 }
542 }
543 foreach (int i_cells, n2c[i_nodes]) {
544 vtkIdType type_cell = m_Grid->GetCellType(cells[i_cells]);
545 if ((type_cell == VTK_WEDGE) && !m_SmoothPrisms) {
546 smooth = false;
547 }
548 }
549 }
550 if (smooth) {
551 vtkIdType id_node = nodes[i_nodes];
552 vtkIdType id_foot = m_IdFoot[id_node];
553 bool use_simple = false;
554 if (id_foot != -1) {
555 if (m_IsSharpNode[id_foot]) {
556 use_simple = true;
557 }
558 }
559 if (use_simple) {
560 simpleNodeMovement(i_nodes);
561 ++N4;
562 } else {
563 vec3_t xn;
564 vec3_t x_old;
565 m_Grid->GetPoint(id_node, x_old.data());
566 resetStencil();
567 bool is_surf = n2bc[i_nodes].size() > 0;
568 int N = 0;
569 foreach (int j_nodes, n2n[i_nodes]) {
570 if (!is_surf || (n2bc[j_nodes].size() > 0)) {
571 vtkIdType id_neigh_node = nodes[j_nodes];
572 m_Grid->GetPoint(id_neigh_node, xn.data());
573 addToStencil(1.0, xn);
574 ++N;
575 }
576 }
577 if (N == 0) {
578 EG_BUG;
579 }
580 vec3_t x_new1 = vec3_t(0,0,0);
581 m_SumC = 0;
582 m_L0 = 0;
583 double L_min = 1e99;
584 foreach (stencil_node_t sn, m_Stencil) {
585 m_SumC += sn.C;
586 x_new1 += sn.C*sn.x;
587 double L = (sn.x - x_old).abs();
588 m_L0 += sn.C*L;
589 L_min = min(L_min, L);
590 }
591 m_L0 /= m_SumC;
592 x_new1 *= 1.0/m_SumC;
593 vec3_t Dx1 = x_new1 - x_old;
594 setDeltas(1e-3*m_L0);
595 m_INodesOpt = i_nodes;
596 vec3_t Dx2(0,0,0);
597 Dx2 = optimise(x_old);
598 vec3_t Dx3 = (-10e-4/func(x_old))*grad_f;
599 correctDx(i_nodes, Dx1);
600 correctDx(i_nodes, Dx2);
601 correctDx(i_nodes, Dx3);
602 vec3_t Dx = Dx1;
603 double f = func(x_old + Dx1);
604 int _N1 = 1;
605 int _N2 = 0;
606 int _N3 = 0;
607 if (f > func(x_old + Dx2)) {
608 Dx = Dx2;
609 f = func(x_old + Dx2);
610 _N1 = 0;
611 _N2 = 1;
612 _N3 = 0;
613 }
614 if (f > func(x_old + Dx3)) {
615 Dx = Dx3;
616 _N1 = 0;
617 _N2 = 0;
618 _N3 = 1;
619 }
620 if (!moveNode(i_nodes, Dx)) {
621 // search for a better place
622 vec3_t x_save = x_old;
623 vec3_t ex(m_LSearch*m_L0/m_NumSearch, 0, 0);
624 vec3_t ey(0,m_LSearch*m_L0/m_NumSearch, 0);
625 vec3_t ez(0,0,m_LSearch*m_L0/m_NumSearch);
626 vec3_t x_best = x_old;
627 double f_min = func(x_old);
628 bool found = false;
629 bool illegal = false;
630 if (!setNewPosition(id_node,x_old)) {
631 illegal = true;
632 }
633 for (int i = -m_NumSearch; i <= m_NumSearch; ++i) {
634 for (int j = -m_NumSearch; j <= m_NumSearch; ++j) {
635 for (int k = -m_NumSearch; k <= m_NumSearch; ++k) {
636 if ((i != 0) || (j != 0) || (k != 0)) {
637 vec3_t Dx = double(i)*ex + double(j)*ey + double(k)*ez;
638 correctDx(i_nodes, Dx);
639 vec3_t x = x_old + Dx;
640 if (setNewPosition(id_node, x)) {
641 m_Grid->GetPoints()->SetPoint(id_node, x_old.data());
642 double f = func(x);
643 if (f < f_min) {
644 f_min = f;
645 x_best = x;
646 found = true;
647 illegal = false;
648 }
649 }
650 }
651 }
652 }
653 }
654 if (found) {
655 m_Grid->GetPoints()->SetPoint(id_node, x_best.data());
656 ++m_NumSearched;
657 } else {
658 ++N_blocked;
659 if (illegal) {
660 ++N_illegal;
661 }
662 }
663 } else {
664 N1 += _N1;
665 N2 += _N2;
666 N3 += _N3;
667 }
668 }
669 }
670 }
671
672 cout << N1 << " type 1 movements (simple)" << endl;
673 cout << N2 << " type 2 movements (Newton)" << endl;
674 cout << N3 << " type 3 movements (gradient)" << endl;
675 cout << N4 << " type 4 movements (sharp edges)" << endl;
676 //cout << N_blocked << " movements blocked" << endl;
677 cout << m_NumSearched << " type 5 movements (search)" << endl;
678 //cout << N_illegal << " nodes in illegal positions" << endl;
679
680 cout << start.secsTo(QTime::currentTime()) << " seconds elapsed" << endl;
681 m_FNew = 0;
682 m_FMaxNew = 0;
683 setPrismWeighting();
684 for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
685 if (prism_node[i_nodes]) {
686 vec3_t x;
687 m_INodesOpt = i_nodes;
688 m_Grid->GetPoint(nodes[i_nodes], x.data());
689 double f = func(x);
690 m_FNew += f;
691 m_FMaxNew = max(m_FMaxNew, f);
692 }
693 }
694 setAllWeighting();
695 cout << "total prism error (old) = " << m_FOld << endl;
696 cout << "total prism error (new) = " << m_FNew << endl;
697 double f_old = max(1e-10, m_FOld);
698 cout << "total prism improvement = " << 100*(1 - m_FNew/f_old) << "%" << endl;
699
700 //writeDebugFile("gridsmoother");
701 //EG_BUG;
702
703 }
704 cout << "done" << endl;
705 }
706
707 double GridSmoother::improvement()
708 {
709 double f_old = max(1e-10, m_FOld);
710 return 1-m_FNew/f_old;
711 }
712
713 void GridSmoother::computeNormals()
714 {
715 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
716 m_NodeNormal.fill(vec3_t(0,0,0), m_Grid->GetNumberOfPoints());
717 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
718 QSet<int> bcs;
719 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
720 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
721 if (isSurface(id_cell, m_Grid)) {
722 int bc = cell_code->GetValue(id_cell);
723 if (m_BoundaryCodes.contains(bc)) {
724 bcs.insert(bc);
725 }
726 }
727 }
728 int num_bcs = bcs.size();
729 QVector<vec3_t> normal(num_bcs, vec3_t(0,0,0));
730 QMap<int,int> bcmap;
731 int i_bc = 0;
732 foreach (int bc, bcs) {
733 bcmap[bc] = i_bc;
734 ++i_bc;
735 }
736 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
737 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
738 if (isSurface(id_cell, m_Grid)) {
739 int bc = cell_code->GetValue(id_cell);
740 if (m_BoundaryCodes.contains(bc)) {
741 vtkIdType N_pts, *pts;
742 m_Grid->GetCellPoints(id_cell, N_pts, pts);
743 vec3_t a, b, c;
744 for (int j = 0; j < N_pts; ++j) {
745 if (pts[j] == id_node) {
746 m_Grid->GetPoint(pts[j], a.data());
747 if (j > 0) {
748 m_Grid->GetPoint(pts[j-1], b.data());
749 } else {
750 m_Grid->GetPoint(pts[N_pts-1], b.data());
751 }
752 if (j < N_pts - 1) {
753 m_Grid->GetPoint(pts[j+1], c.data());
754 } else {
755 m_Grid->GetPoint(pts[0], c.data());
756 }
757 }
758 }
759 vec3_t u = b - a;
760 vec3_t v = c - a;
761 double alpha = GeometryTools::angle(u, v);
762 vec3_t n = u.cross(v);
763 n.normalise();
764 normal[bcmap[bc]] += alpha*n;
765 }
766 }
767 }
768 for (int i = 0; i < num_bcs; ++i) {
769 normal[i].normalise();
770 }
771 if (num_bcs > 0) {
772 if (num_bcs > 1) {
773 if (num_bcs == 3) {
774 for (int i = 0; i < num_bcs; ++i) {
775 for (int j = i + 1; j < num_bcs; ++j) {
776 vec3_t n = normal[i] + normal[j];
777 n.normalise();
778 m_NodeNormal[id_node] += n;
779 }
780 }
781 } else {
782 for (int i = 0; i < num_bcs; ++i) {
783 m_NodeNormal[id_node] += normal[i];
784 }
785 }
786 } else {
787 m_NodeNormal[id_node] = normal[0];
788 }
789 m_NodeNormal[id_node].normalise();
790 }
791 }
792 m_IsSharpNode.fill(false, m_Grid->GetNumberOfPoints());
793 for (int i = 0; i < m_Part.getNumberOfCells(); ++i) {
794 vtkIdType id_cell1 = m_Part.globalCell(i);
795 if (isSurface(id_cell1, m_Grid)) {
796 if (m_BoundaryCodes.contains(cell_code->GetValue(id_cell1))) {
797 vec3_t n1 = cellNormal(m_Grid, id_cell1);
798 n1.normalise();
799 vtkIdType N_pts, *pts;
800 m_Grid->GetCellPoints(id_cell1, N_pts, pts);
801 for (int j = 0; j < m_Part.c2cLSize(i); ++j) {
802 vtkIdType id_cell2 = m_Part.c2cLG(i, j);
803 if (m_BoundaryCodes.contains(cell_code->GetValue(id_cell2))) {
804 vec3_t n2 = cellNormal(m_Grid, id_cell2);
805 n2.normalise();
806 if (GeometryTools::angle(n1, n2) > m_CritAngle) {
807 vtkIdType id_node1 = pts[j];
808 vtkIdType id_node2 = pts[0];
809 if (j < m_Part.c2cLSize(i) - 1) {
810 id_node2 = pts[j+1];
811 }
812 if (id_node1 == 7 || id_node2 == 7) {
813 cout << "Hello" << endl;
814 }
815 vec3_t x_node1, x_node2;
816 m_Grid->GetPoint(id_node1, x_node1.data());
817 m_Grid->GetPoint(id_node2, x_node2.data());
818 vec3_t x_corner = 0.5*(x_node1 + x_node2);
819 vec3_t x_cell1 = cellCentre(m_Grid, id_cell1);
820 vec3_t x_cell2 = cellCentre(m_Grid, id_cell2);
821 vec3_t v_cell1 = x_cell2 - x_cell1;
822 vec3_t v_cell2 = x_cell1 - x_cell2;
823 if (n1*v_cell1 > 0 || n2*v_cell2 > 0) {
824 m_IsSharpNode[id_node1] = true;
825 m_IsSharpNode[id_node2] = true;
826 }
827 }
828 }
829 }
830 }
831 }
832 }
833 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
834 if (m_IsSharpNode[id_node]) {
835 QVector<vec3_t> normals;
836 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
837 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
838 if (isSurface(id_cell, m_Grid)) {
839 if (m_BoundaryCodes.contains(cell_code->GetValue(id_cell))) {
840 vec3_t n = cellNormal(m_Grid, id_cell);
841 n.normalise();
842 normals.push_back(n);
843 }
844 }
845 }
846 int N = 0;
847 foreach (vec3_t n1, normals) {
848 foreach (vec3_t n2, normals) {
849 if (GeometryTools::angle(n1, n2) > m_CritAngle) {
850 ++N;
851 }
852 }
853 }
854 }
855 }
856 }
857
858 void GridSmoother::computeFeet()
859 {
860 m_IdFoot.fill(-1, m_Grid->GetNumberOfPoints());
861 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
862 if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
863 vtkIdType N_pts, *pts;
864 m_Grid->GetCellPoints(id_cell, N_pts, pts);
865 m_IdFoot[pts[3]] = pts[0];
866 m_IdFoot[pts[4]] = pts[1];
867 m_IdFoot[pts[5]] = pts[2];
868 }
869 }
870 }
871
872 void GridSmoother::simpleNodeMovement(int i_nodes)
873 {
874 EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired" );
875 vtkIdType id_node = m_Part.globalNode(i_nodes);
876 vtkIdType id_foot = m_IdFoot[id_node];
877 if (id_foot != -1) {
878 vec3_t x_foot, x_node;
879 m_Grid->GetPoint(id_foot, x_foot.data());
880 m_Grid->GetPoint(id_node, x_node.data());
881 double L = cl->GetValue(id_foot);
882 vec3_t x_new = x_foot + m_RelativeHeight*L*m_NodeNormal[id_foot];
883 vec3_t Dx = x_new - x_node;
884 correctDx(i_nodes, Dx);
885 m_L[id_node] = L;
886 moveNode(i_nodes, Dx);
887 }
888 }
889
890 void GridSmoother::operateSimple()
891 {
892 cout << "performing simple boundary layer adjustment" << endl;
893 computeFeet();
894 m_L.fill(0, m_Grid->GetNumberOfPoints());
895 l2g_t nodes = m_Part.getNodes();
896 for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
897 simpleNodeMovement(i_nodes);
898 }
899 }
900
901 void GridSmoother::operatePostSmoothing()
902 {
903 QVector<bool> top_node(m_Grid->GetNumberOfPoints(), false);
904 QVector<vec3_t> x_new(m_Grid->GetNumberOfPoints());
905 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
906 if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
907 vtkIdType N_pts, *pts;
908 m_Grid->GetCellPoints(id_cell, N_pts, pts);
909 top_node[pts[3]] = true;
910 top_node[pts[4]] = true;
911 top_node[pts[5]] = true;
912 }
913 }
914 for (vtkIdType id_node1 = 0; id_node1 < m_Grid->GetNumberOfPoints(); ++id_node1) {
915 m_Grid->GetPoint(id_node1, x_new[id_node1].data());
916 if (top_node[id_node1]) {
917 double w = m_PostSmoothingStrength;
918 vec3_t x_av(0,0,0);
919 vec3_t x_old;
920 m_Grid->GetPoint(id_node1, x_old.data());
921 QSet<vtkIdType> ids;
922 for (int i = 0; i < m_Part.n2cGSize(id_node1); ++i) {
923 vtkIdType id_cell = m_Part.n2cGG(id_node1, i);
924 if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
925 vtkIdType N_pts, *pts;
926 m_Grid->GetCellPoints(id_cell, N_pts, pts);
927 for (int j = 0; j < 6; ++j) {
928 if (top_node[pts[j]] && (pts[j] != id_node1)) {
929 ids.insert(pts[j]);
930 }
931 }
932 }
933 }
934 foreach (vtkIdType id_node2, ids) {
935 vec3_t x;
936 m_Grid->GetPoint(id_node2, x.data());
937 x_av += x;
938 }
939 if (ids.size() > 0) {
940 x_av *= 1.0/ids.size();
941 x_new[id_node1] = w*x_av + (1-w)*x_old;
942 }
943 }
944 }
945 l2g_t nodes = m_Part.getNodes();
946 for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
947 vtkIdType id_node = m_Part.globalNode(i_nodes);
948 vec3_t x_old;
949 m_Grid->GetPoint(id_node, x_old.data());
950 vec3_t Dx = x_new[id_node] - x_old;
951 correctDx(i_nodes, Dx);
952 moveNode(i_nodes, Dx);
953 }
954 }
955
956 void GridSmoother::operate()
957 {
958 markNodes();
959 computeNormals();
960 if (m_SimpleOperation) {
961 operateSimple();
962 } else if (m_PostOperation) {
963 operatePostSmoothing();
964 } else {
965 operateOptimisation();
966 }
967 }
968 void GridSmoother::writeDebugFile(QString file_name)
969 {
970 QVector<vtkIdType> bcells;
971 getSurfaceCells(m_BoundaryCodes, bcells, m_Grid);
972 MeshPartition bpart(m_Grid);
973 bpart.setCells(bcells);
974 /*
975 QVector<vtkIdType> foot2field(bpart.getNumberOfNodes(), -1);
976 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
977 vtkIdType id_foot = m_IdFoot[id_node];
978 if (id_foot != -1) {
979 if (bpart.localNode(id_foot) == -1) {
980 EG_BUG;
981 }
982 foot2field[bpart.localNode(id_foot)] = id_node;
983 }
984 }
985 */
986 file_name = GuiMainWindow::pointer()->getCwd() + "/" + file_name + ".vtk";
987 QFile file(file_name);
988 file.open(QIODevice::WriteOnly);
989 QTextStream f(&file);
990 f << "# vtk DataFile Version 2.0\n";
991 f << "m_NodeNormal\n";
992 f << "ASCII\n";
993 f << "DATASET UNSTRUCTURED_GRID\n";
994 f << "POINTS " << bpart.getNumberOfNodes() << " float\n";
995 for (int i = 0; i < bpart.getNumberOfNodes(); ++i) {
996 vec3_t x;
997 vtkIdType id_node = bpart.globalNode(i);
998 m_Grid->GetPoint(id_node, x.data());
999 f << x[0] << " " << x[1] << " " << x[2] << "\n";
1000 }
1001 f << "CELLS " << bpart.getNumberOfCells() << " " << 4*bpart.getNumberOfCells() << "\n";
1002 for (int i = 0; i < bpart.getNumberOfCells(); ++ i) {
1003 vtkIdType id_cell = bpart.globalCell(i);
1004 vtkIdType N_pts, *pts;
1005 if (m_Grid->GetCellType(id_cell) != VTK_TRIANGLE) {
1006 EG_BUG;
1007 }
1008 m_Grid->GetCellPoints(id_cell, N_pts, pts);
1009 QVector<int> nds(3);
1010 for (int j = 0; j < 3; ++j) {
1011 nds[j] = bpart.localNode(pts[j]);
1012 }
1013 f << "3 " << nds[0] << " " << nds[1] << " " << nds[2] << "\n";
1014 }
1015
1016 f << "CELL_TYPES " << bpart.getNumberOfCells() << "\n";
1017 for (int i = 0; i < bpart.getNumberOfCells(); ++ i) {
1018 f << VTK_TRIANGLE << "\n";
1019 }
1020 f << "POINT_DATA " << bpart.getNumberOfNodes() << "\n";
1021 f << "VECTORS N float\n";
1022
1023 for (int i = 0; i < bpart.getNumberOfNodes(); ++i) {
1024 vtkIdType id_node = bpart.globalNode(i);
1025 f << m_NodeNormal[id_node][0] << " " << m_NodeNormal[id_node][1] << " " << m_NodeNormal[id_node][2] << "\n";
1026 }
1027 }
1028
1029
Clone of the official engrid repository