?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
moved gridsmoother.* from master branch
[engrid.git] / src / gridsmoother.cpp
blob4482696ee582862be68bf5018d02f2cfd692b385
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 if (GeometryTools::tetraVol(xtet[0], xtet[1], xtet[2], xtet[3]) < 0) {
173 ok = false;
174 }
175 }
176 if (ok) {
177 break;
178 }
179 }
180 if (!ok) {
181 move = false;
182 }
183 }
184 }
185
186
187 if (!move) {
188 m_Grid->GetPoints()->SetPoint(id_node, x_old.data());
189 }
190 return move;
191 }
192
193 void GridSmoother::correctDx(int i_nodes, vec3_t &Dx)
194 {
195 l2g_t nodes = m_Part.getNodes();
196 l2l_t n2c = m_Part.getN2C();
197 for (int i_boundary_correction = 0; i_boundary_correction < m_NumBoundaryCorrections; ++i_boundary_correction) {
198 foreach (vtkIdType id_cell, n2c[i_nodes]) {
199 if (isSurface(id_cell, m_Grid)) {
200 double A = GeometryTools::cellVA(m_Grid, id_cell);
201 if (A > 1e-20) {
202 vec3_t n = GeometryTools::cellNormal(m_Grid, id_cell);
203 n.normalise();
204 Dx -= (n*Dx)*n;
205 }
206 } else {
207 if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
208 vtkIdType N_pts, *pts;
209 m_Grid->GetCellPoints(id_cell, N_pts, pts);
210 vtkIdType id_surf_node = -1;
211 if (pts[3] == nodes[i_nodes]) id_surf_node = pts[0];
212 if (pts[4] == nodes[i_nodes]) id_surf_node = pts[1];
213 if (pts[5] == nodes[i_nodes]) id_surf_node = pts[2];
214 if (id_surf_node != -1) {
215 vec3_t x0,x1,x2;
216 m_Grid->GetPoint(pts[0],x0.data());
217 m_Grid->GetPoint(pts[1],x1.data());
218 m_Grid->GetPoint(pts[2],x2.data());
219 vec3_t a = x1-x0;
220 vec3_t b = x2-x0;
221 vec3_t c = b-a;
222 double L = (a.abs()+b.abs()+c.abs())/3.0;
223 vec3_t n = b.cross(a);
224 n.normalise();
225 vec3_t x_old;
226 m_Grid->GetPoint(nodes[i_nodes],x_old.data());
227 vec3_t x_new = x_old + Dx - x0;
228 if ( (n*x_new) <= 0 ) {
229 x_new -= (x_new*n)*n;
230 x_new += 1e-4*L*n;
231 x_new += x0;
232 Dx = x_new - x_old;
233 }
234 }
235 }
236 }
237 }
238 }
239 }
240
241 bool GridSmoother::moveNode(int i_nodes, vec3_t &Dx)
242 {
243 l2g_t nodes = m_Part.getNodes();
244 vtkIdType id_node = nodes[i_nodes];
245 vec3_t x_old;
246 m_Grid->GetPoint(id_node, x_old.data());
247 bool moved = false;
248 for (int i_relaxation = 0; i_relaxation < m_NumRelaxations; ++i_relaxation) {
249 if (setNewPosition(id_node, x_old + Dx)) {
250 moved = true;
251 break;
252 }
253 Dx *= 0.5;
254 }
255 return moved;
256 }
257
258 double GridSmoother::errThickness(double x)
259 {
260 if (x > 1) x = 2 - x;
261 double delta = 0.01;
262 double a = 5.0;
263 double err = 0;
264 if (x < 0) err = 1.0/delta - 1.0/(a+delta);
265 else if (x < 1) err = 1/(a*x+delta) - 1.0/(a+delta);
266 return err;
267 }
268
269 double GridSmoother::func(vec3_t x)
270 {
271 l2g_t nodes = m_Part.getNodes();
272 l2g_t cells = m_Part.getCells();
273 l2l_t n2c = m_Part.getN2C();
274 l2l_t c2c = m_Part.getC2C();
275
276 vec3_t x_old;
277 m_Grid->GetPoint(nodes[m_INodesOpt], x_old.data());
278 m_Grid->GetPoints()->SetPoint(nodes[m_INodesOpt], x.data());
279 double f = 0;
280 double f13 = 1.0/3.0;
281 double f14 = 0.25;
282
283 vec3_t n_node(1,0,0);
284 QList<vec3_t> n_pri;
285
286 foreach (int i_cells, n2c[m_INodesOpt]) {
287 vtkIdType id_cell = cells[i_cells];
288 if (isVolume(id_cell, m_Grid)) {
289 vtkIdType type_cell = m_Grid->GetCellType(id_cell);
290 vtkIdType N_pts, *pts;
291 m_Grid->GetCellPoints(id_cell, N_pts, pts);
292 QVector<vec3_t> xn(N_pts);
293 for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
294 m_Grid->GetPoint(pts[i_pts],xn[i_pts].data());
295 }
296 if (type_cell == VTK_TETRA) {
297 double L_max = 1e-20;
298 L_max = max((xn[0]-xn[1]).abs(), L_max);
299 L_max = max((xn[0]-xn[2]).abs(), L_max);
300 L_max = max((xn[0]-xn[3]).abs(), L_max);
301 L_max = max((xn[1]-xn[2]).abs(), L_max);
302 L_max = max((xn[1]-xn[3]).abs(), L_max);
303 L_max = max((xn[2]-xn[3]).abs(), L_max);
304 double V1 = GeometryTools::cellVA(m_Grid, id_cell, true);
305 double V2 = 0.1178511301977579207*L_max*L_max*L_max;
306 double e = fabs(V1-V2)/V2;
307 f += m_WTet*pow(e, m_ETet);
308 }
309 if (type_cell == VTK_WEDGE) {
310 double L = 0;
311 L += (xn[0]-xn[1]).abs();
312 L += (xn[0]-xn[2]).abs();
313 L += (xn[1]-xn[2]).abs();
314 L *= m_H/3.0;
315 vec3_t a = xn[2]-xn[0];
316 vec3_t b = xn[1]-xn[0];
317 vec3_t c = xn[5]-xn[3];
318 vec3_t d = xn[4]-xn[3];
319 vec3_t n_face[5];
320 vec3_t x_face[5];
321 n_face[0] = GeometryTools::triNormal(xn[0],xn[1],xn[2]);
322 n_face[1] = GeometryTools::triNormal(xn[3],xn[5],xn[4]);
323 n_face[2] = GeometryTools::quadNormal(xn[0],xn[3],xn[4],xn[1]);
324 n_face[3] = GeometryTools::quadNormal(xn[1],xn[4],xn[5],xn[2]);
325 n_face[4] = GeometryTools::quadNormal(xn[0],xn[2],xn[5],xn[3]);
326 x_face[0] = f13*(xn[0]+xn[1]+xn[2]);
327 x_face[1] = f13*(xn[3]+xn[4]+xn[5]);
328 x_face[2] = f14*(xn[0]+xn[3]+xn[4]+xn[1]);
329 x_face[3] = f14*(xn[1]+xn[4]+xn[5]+xn[2]);
330 x_face[4] = f14*(xn[0]+xn[2]+xn[5]+xn[3]);
331
332 double A1 = 0.5*n_face[0].abs();
333 double A2 = 0.5*n_face[1].abs();
334 for (int i_face = 0; i_face < 5; ++i_face) {
335 n_face[i_face].normalise();
336 }
337 m_IdFoot[pts[3]] = pts[0];
338 m_IdFoot[pts[4]] = pts[1];
339 m_IdFoot[pts[5]] = pts[2];
340 if (nodes[m_INodesOpt] == pts[3]) {
341 n_node = xn[3]-xn[0];
342 n_pri.append(n_face[1]);
343 m_L[nodes[m_INodesOpt]] = L;
344 }
345 if (nodes[m_INodesOpt] == pts[4]) {
346 n_node = xn[4]-xn[1];
347 n_pri.append(n_face[1]);
348 m_L[nodes[m_INodesOpt]] = L;
349 }
350 if (nodes[m_INodesOpt] == pts[5]) {
351 n_node = xn[5]-xn[2];
352 n_pri.append(n_face[1]);
353 m_L[nodes[m_INodesOpt]] = L;
354 }
355 vec3_t v0 = xn[0]-xn[3];
356 vec3_t v1 = xn[1]-xn[4];
357 vec3_t v2 = xn[2]-xn[5];
358
359 double h0 = v0*n_face[0];
360 double h1 = v1*n_face[0];
361 double h2 = v2*n_face[0];
362 if (h0 > 0.5*L) h0 = max(v0.abs(),h0);
363 if (h1 > 0.5*L) h1 = max(v1.abs(),h1);
364 if (h2 > 0.5*L) h2 = max(v2.abs(),h2);
365
366 {
367
368 double e1 = errThickness(h0/L);
369 double e2 = errThickness(h1/L);
370 double e3 = errThickness(h2/L);
371 double e = max(e1,max(e2,e3));
372 //f += w_h*f13*e1;
373 //f += w_h*f13*e2;
374 //f += w_h*f13*e3;
375 //err_pria->SetValue(id_cell,f13*(e1+e2+e3));
376
377 f += m_WH*e;
378 }
379 if ((h0 > 0.01*L) && (h1 > 0.01*L) && (h2 > 0.01*L)) {
380 v0.normalise();
381 v1.normalise();
382 v2.normalise();
383 double e1 = f13*(1-v0*v1);
384 double e2 = f13*(1-v0*v2);
385 double e3 = f13*(1-v1*v2);
386 f += m_WPar*e1;
387 f += m_WPar*e2;
388 f += m_WPar*e3;
389 }
390 if ((h0 > 0.01*L) && (h1 > 0.01*L) && (h2 > 0.01*L)) {
391 double e = (1+n_face[0]*n_face[1]);
392 f += m_WN*e;
393 }
394 if ((h0 > 0.01*L) && (h1 > 0.01*L) && (h2 > 0.01*L)) {
395 double e = sqr(0.5*(A1-A2)/A1);
396 f += m_WA*e;
397 }
398 if ((h0 > 0.01*L) && (h1 > 0.01*L) && (h2 > 0.01*L)) {
399 double e_skew = 0;
400 double e_orth = 0;
401 int N = 0;
402 vec3_t xc = cellCentre(m_Grid, id_cell);
403 for (int i_face = 0; i_face < 5; ++i_face) {
404 int i_cells_neigh = c2c[i_cells][i_face];
405 if (i_cells_neigh != -1) {
406 vtkIdType id_neigh_cell = cells[i_cells_neigh];
407 if (isVolume(id_neigh_cell, m_Grid)) {
408 vec3_t vc = cellCentre(m_Grid, id_neigh_cell) - xc;
409 vec3_t vf = x_face[i_face] - xc;
410 vc.normalise();
411 vf.normalise();
412 e_skew += (1-vc*vf);
413 e_orth += (1-vc*n_face[i_face]);
414 ++N;
415 }
416 }
417 }
418 e_skew /= N;
419 e_orth /= N;
420 f += m_WSkew*e_skew + m_WOrth*e_orth;
421 }
422
423 double f_sharp2 = 0;
424 for (int j = 2; j <= 4; ++j) {
425 vtkIdType id_ncell = c2c[id_cell][j];
426 if (id_ncell != -1) {
427 if (m_Grid->GetCellType(id_ncell) == VTK_WEDGE) {
428 vtkIdType N_pts, *pts;
429 m_Grid->GetCellPoints(id_ncell, N_pts, pts);
430 QVector<vec3_t> x(3);
431 for (int i_pts = 3; i_pts <= 5; ++i_pts) {
432 m_Grid->GetPoint(pts[i_pts],x[i_pts-3].data());
433 }
434 vec3_t n = GeometryTools::triNormal(x[0],x[2],x[1]);
435 n.normalise();
436 f_sharp2 += pow(fabs(1-n_face[1]*n), m_ESharp2);
437 }
438 }
439 }
440 f += m_WSharp2*f_sharp2;
441 }
442 }
443 }
444 m_Grid->GetPoints()->SetPoint(nodes[m_INodesOpt], x_old.data());
445 n_node.normalise();
446 {
447 double f_sharp1 = 0;
448 foreach (vec3_t n, n_pri) {
449 f_sharp1 += pow(fabs(1-n_node*n), m_ESharp1);
450 }
451 f += m_WSharp1*f_sharp1;
452 }
453 return f;
454 }
455
456 void GridSmoother::resetStencil()
457 {
458 m_Stencil.clear();
459 }
460
461 void GridSmoother::addToStencil(double C, vec3_t x)
462 {
463 stencil_node_t sn;
464 sn.x = x;
465 sn.C = C;
466 m_Stencil.append(sn);
467 }
468
469 void GridSmoother::operateOptimisation()
470 {
471 markNodes();
472 findCriticalTetras();
473
474 l2g_t nodes = m_Part.getNodes();
475 l2g_t cells = m_Part.getCells();
476 g2l_t _nodes = m_Part.getLocalNodes();
477 l2l_t n2c = m_Part.getN2C();
478 l2l_t n2n = m_Part.getN2N();
479
480 EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
481 EG_VTKDCN(vtkIntArray, node_status, m_Grid, "node_status");
482 EG_VTKDCN(vtkIntArray, node_layer, m_Grid, "node_layer");
483
484 m_IdFoot.fill(-1, m_Grid->GetNumberOfPoints());
485 m_L.fill(0, m_Grid->GetNumberOfPoints());
486
487 QVector<QSet<int> > n2bc(nodes.size());
488 QVector<bool> prism_node(nodes.size(),false);
489 foreach (vtkIdType id_cell, cells) {
490 if (isSurface(id_cell, m_Grid)) {
491 vtkIdType N_pts, *pts;
492 m_Grid->GetCellPoints(id_cell, N_pts, pts);
493 for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
494 n2bc[_nodes[pts[i_pts]]].insert(bc->GetValue(id_cell));
495 }
496 }
497 if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
498 vtkIdType N_pts, *pts;
499 m_Grid->GetCellPoints(id_cell, N_pts, pts);
500 for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
501 if (_nodes[pts[i_pts]] != -1) {
502 prism_node[_nodes[pts[i_pts]]] = true;
503 }
504 }
505 }
506 }
507
508 m_FOld = 0;
509 m_FMaxOld = 0;
510 setPrismWeighting();
511 for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
512 if (prism_node[i_nodes]) {
513 vec3_t x;
514 m_INodesOpt = i_nodes;
515 m_Grid->GetPoint(nodes[i_nodes], x.data());
516 double f = func(x);
517 m_FOld += f;
518 m_FMaxOld = max(m_FMaxOld, f);
519 }
520 }
521 setAllWeighting();
522
523 cout << "\nsmoothing volume mesh (" << m_NumMarkedNodes << " nodes)" << endl;
524 for (int i_iterations = 0; i_iterations < m_NumIterations; ++i_iterations) {
525 cout << "iteration " << i_iterations+1 << "/" << m_NumIterations << endl;
526 int N_blocked = 0;
527 int m_NumSearched = 0;
528 int N_illegal = 0;
529 int N1 = 0;
530 int N2 = 0;
531 int N3 = 0;
532 int N4 = 0;
533 QTime start = QTime::currentTime();
534 for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
535 bool smooth = m_NodeMarked[nodes[i_nodes]];
536 if (smooth) {
537 foreach (int bc, n2bc[i_nodes]) {
538 if (m_BoundaryCodes.contains(bc) || (m_BoundaryCodes.size() == 0)) {
539 smooth = false;
540 }
541 }
542 foreach (int i_cells, n2c[i_nodes]) {
543 vtkIdType type_cell = m_Grid->GetCellType(cells[i_cells]);
544 if ((type_cell == VTK_WEDGE) && !m_SmoothPrisms) {
545 smooth = false;
546 }
547 }
548 }
549 if (smooth) {
550 vtkIdType id_node = nodes[i_nodes];
551 vtkIdType id_foot = m_IdFoot[id_node];
552 bool use_simple = false;
553 if (id_foot != -1) {
554 if (m_IsSharpNode[id_foot]) {
555 use_simple = true;
556 }
557 }
558 if (use_simple) {
559 simpleNodeMovement(i_nodes);
560 ++N4;
561 } else {
562 vec3_t xn;
563 vec3_t x_old;
564 m_Grid->GetPoint(id_node, x_old.data());
565 resetStencil();
566 bool is_surf = n2bc[i_nodes].size() > 0;
567 int N = 0;
568 foreach (int j_nodes, n2n[i_nodes]) {
569 if (!is_surf || (n2bc[j_nodes].size() > 0)) {
570 vtkIdType id_neigh_node = nodes[j_nodes];
571 m_Grid->GetPoint(id_neigh_node, xn.data());
572 addToStencil(1.0, xn);
573 ++N;
574 }
575 }
576 if (N == 0) {
577 EG_BUG;
578 }
579 vec3_t x_new1 = vec3_t(0,0,0);
580 m_SumC = 0;
581 m_L0 = 0;
582 double L_min = 1e99;
583 foreach (stencil_node_t sn, m_Stencil) {
584 m_SumC += sn.C;
585 x_new1 += sn.C*sn.x;
586 double L = (sn.x - x_old).abs();
587 m_L0 += sn.C*L;
588 L_min = min(L_min, L);
589 }
590 m_L0 /= m_SumC;
591 x_new1 *= 1.0/m_SumC;
592 vec3_t Dx1 = x_new1 - x_old;
593 setDeltas(1e-3*m_L0);
594 m_INodesOpt = i_nodes;
595 vec3_t Dx2(0,0,0);
596 Dx2 = optimise(x_old);
597 vec3_t Dx3 = (-10e-4/func(x_old))*grad_f;
598 correctDx(i_nodes, Dx1);
599 correctDx(i_nodes, Dx2);
600 correctDx(i_nodes, Dx3);
601 vec3_t Dx = Dx1;
602 double f = func(x_old + Dx1);
603 int _N1 = 1;
604 int _N2 = 0;
605 int _N3 = 0;
606 if (f > func(x_old + Dx2)) {
607 Dx = Dx2;
608 f = func(x_old + Dx2);
609 _N1 = 0;
610 _N2 = 1;
611 _N3 = 0;
612 }
613 if (f > func(x_old + Dx3)) {
614 Dx = Dx3;
615 _N1 = 0;
616 _N2 = 0;
617 _N3 = 1;
618 }
619 if (!moveNode(i_nodes, Dx)) {
620 // search for a better place
621 vec3_t x_save = x_old;
622 vec3_t ex(m_LSearch*m_L0/m_NumSearch, 0, 0);
623 vec3_t ey(0,m_LSearch*m_L0/m_NumSearch, 0);
624 vec3_t ez(0,0,m_LSearch*m_L0/m_NumSearch);
625 vec3_t x_best = x_old;
626 double f_min = func(x_old);
627 bool found = false;
628 bool illegal = false;
629 if (!setNewPosition(id_node,x_old)) {
630 illegal = true;
631 }
632 for (int i = -m_NumSearch; i <= m_NumSearch; ++i) {
633 for (int j = -m_NumSearch; j <= m_NumSearch; ++j) {
634 for (int k = -m_NumSearch; k <= m_NumSearch; ++k) {
635 if ((i != 0) || (j != 0) || (k != 0)) {
636 vec3_t Dx = double(i)*ex + double(j)*ey + double(k)*ez;
637 correctDx(i_nodes, Dx);
638 vec3_t x = x_old + Dx;
639 if (setNewPosition(id_node, x)) {
640 m_Grid->GetPoints()->SetPoint(id_node, x_old.data());
641 double f = func(x);
642 if (f < f_min) {
643 f_min = f;
644 x_best = x;
645 found = true;
646 illegal = false;
647 }
648 }
649 }
650 }
651 }
652 }
653 if (found) {
654 m_Grid->GetPoints()->SetPoint(id_node, x_best.data());
655 ++m_NumSearched;
656 } else {
657 ++N_blocked;
658 if (illegal) {
659 ++N_illegal;
660 }
661 }
662 } else {
663 N1 += _N1;
664 N2 += _N2;
665 N3 += _N3;
666 }
667 }
668 }
669 }
670
671 cout << N1 << " type 1 movements (simple)" << endl;
672 cout << N2 << " type 2 movements (Newton)" << endl;
673 cout << N3 << " type 3 movements (gradient)" << endl;
674 cout << N4 << " type 4 movements (sharp edges)" << endl;
675 //cout << N_blocked << " movements blocked" << endl;
676 cout << m_NumSearched << " type 5 movements (search)" << endl;
677 //cout << N_illegal << " nodes in illegal positions" << endl;
678
679 cout << start.secsTo(QTime::currentTime()) << " seconds elapsed" << endl;
680 m_FNew = 0;
681 m_FMaxNew = 0;
682 setPrismWeighting();
683 for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
684 if (prism_node[i_nodes]) {
685 vec3_t x;
686 m_INodesOpt = i_nodes;
687 m_Grid->GetPoint(nodes[i_nodes], x.data());
688 double f = func(x);
689 m_FNew += f;
690 m_FMaxNew = max(m_FMaxNew, f);
691 }
692 }
693 setAllWeighting();
694 cout << "total prism error (old) = " << m_FOld << endl;
695 cout << "total prism error (new) = " << m_FNew << endl;
696 double f_old = max(1e-10, m_FOld);
697 cout << "total prism improvement = " << 100*(1 - m_FNew/f_old) << "%" << endl;
698 }
699 cout << "done" << endl;
700 }
701
702 double GridSmoother::improvement()
703 {
704 double f_old = max(1e-10, m_FOld);
705 return 1-m_FNew/f_old;
706 }
707
708 void GridSmoother::computeNormals()
709 {
710 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
711 m_NodeNormal.fill(vec3_t(0,0,0), m_Grid->GetNumberOfPoints());
712 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
713 QSet<int> bcs;
714 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
715 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
716 if (isSurface(id_cell, m_Grid)) {
717 int bc = cell_code->GetValue(id_cell);
718 if (m_BoundaryCodes.contains(bc)) {
719 bcs.insert(bc);
720 }
721 }
722 }
723 int num_bcs = bcs.size();
724 QVector<vec3_t> normal(num_bcs, vec3_t(0,0,0));
725 QMap<int,int> bcmap;
726 int i_bc = 0;
727 foreach (int bc, bcs) {
728 bcmap[bc] = i_bc;
729 ++i_bc;
730 }
731 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
732 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
733 if (isSurface(id_cell, m_Grid)) {
734 int bc = cell_code->GetValue(id_cell);
735 if (m_BoundaryCodes.contains(bc)) {
736 vtkIdType N_pts, *pts;
737 m_Grid->GetCellPoints(id_cell, N_pts, pts);
738 vec3_t a, b, c;
739 for (int j = 0; j < N_pts; ++j) {
740 if (pts[j] == id_node) {
741 m_Grid->GetPoint(pts[j], a.data());
742 if (j > 0) {
743 m_Grid->GetPoint(pts[j-1], b.data());
744 } else {
745 m_Grid->GetPoint(pts[N_pts-1], b.data());
746 }
747 if (j < N_pts - 1) {
748 m_Grid->GetPoint(pts[j+1], c.data());
749 } else {
750 m_Grid->GetPoint(pts[0], c.data());
751 }
752 }
753 }
754 vec3_t u = b - a;
755 vec3_t v = c - a;
756 double alpha = GeometryTools::angle(u, v);
757 vec3_t n = u.cross(v);
758 n.normalise();
759 normal[bcmap[bc]] += alpha*n;
760 }
761 }
762 }
763 for (int i = 0; i < num_bcs; ++i) {
764 normal[i].normalise();
765 }
766 if (num_bcs > 0) {
767 if (num_bcs > 1) {
768 for (int i = 0; i < num_bcs; ++i) {
769 for (int j = i + 1; j < num_bcs; ++j) {
770 vec3_t n = normal[i] + normal[j];
771 n.normalise();
772 m_NodeNormal[id_node] += n;
773 }
774 }
775 } else {
776 m_NodeNormal[id_node] = normal[0];
777 }
778 m_NodeNormal[id_node].normalise();
779 }
780 }
781 m_IsSharpNode.fill(false, m_Grid->GetNumberOfPoints());
782 for (int i = 0; i < m_Part.getNumberOfCells(); ++i) {
783 vtkIdType id_cell1 = m_Part.globalCell(i);
784 if (isSurface(id_cell1, m_Grid)) {
785 if (m_BoundaryCodes.contains(cell_code->GetValue(id_cell1))) {
786 vec3_t n1 = cellNormal(m_Grid, id_cell1);
787 n1.normalise();
788 vtkIdType N_pts, *pts;
789 m_Grid->GetCellPoints(id_cell1, N_pts, pts);
790 for (int j = 0; j < m_Part.c2cLSize(i); ++j) {
791 vtkIdType id_cell2 = m_Part.c2cLG(i, j);
792 if (m_BoundaryCodes.contains(cell_code->GetValue(id_cell2))) {
793 vec3_t n2 = cellNormal(m_Grid, id_cell2);
794 n2.normalise();
795 if (GeometryTools::angle(n1, n2) > m_CritAngle) {
796 vtkIdType id_node1 = pts[j];
797 vtkIdType id_node2 = pts[0];
798 if (j < m_Part.c2cLSize(i) - 1) {
799 id_node2 = pts[j+1];
800 }
801 m_IsSharpNode[id_node1] = true;
802 m_IsSharpNode[id_node2] = true;
803 }
804 }
805 }
806 }
807 }
808 }
809 m_IsTripleNode.fill(false, m_Grid->GetNumberOfPoints());
810 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
811 if (m_IsSharpNode[id_node]) {
812 QVector<vec3_t> normals;
813 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
814 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
815 if (isSurface(id_cell, m_Grid)) {
816 if (m_BoundaryCodes.contains(cell_code->GetValue(id_cell))) {
817 vec3_t n = cellNormal(m_Grid, id_cell);
818 n.normalise();
819 normals.push_back(n);
820 }
821 }
822 }
823 int N = 0;
824 foreach (vec3_t n1, normals) {
825 foreach (vec3_t n2, normals) {
826 if (GeometryTools::angle(n1, n2) > m_CritAngle) {
827 ++N;
828 }
829 }
830 }
831 }
832 }
833 }
834
835 void GridSmoother::computeFeet()
836 {
837 m_IdFoot.fill(-1, m_Grid->GetNumberOfPoints());
838 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
839 if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
840 vtkIdType N_pts, *pts;
841 m_Grid->GetCellPoints(id_cell, N_pts, pts);
842 m_IdFoot[pts[3]] = pts[0];
843 m_IdFoot[pts[4]] = pts[1];
844 m_IdFoot[pts[5]] = pts[2];
845 }
846 }
847 }
848
849 void GridSmoother::simpleNodeMovement(int i_nodes)
850 {
851 EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired" );
852 vtkIdType id_node = m_Part.globalNode(i_nodes);
853 vtkIdType id_foot = m_IdFoot[id_node];
854 if (id_foot != -1) {
855 vec3_t x_foot, x_node;
856 m_Grid->GetPoint(id_foot, x_foot.data());
857 m_Grid->GetPoint(id_node, x_node.data());
858 double L = cl->GetValue(id_foot);
859 vec3_t x_new = x_foot + m_RelativeHeight*L*m_NodeNormal[id_foot];
860 vec3_t Dx = x_new - x_node;
861 correctDx(i_nodes, Dx);
862 m_L[id_node] = L;
863 moveNode(i_nodes, Dx);
864 }
865 }
866
867 void GridSmoother::operateSimple()
868 {
869 cout << "performing simple boundary layer adjustment" << endl;
870 computeFeet();
871 m_L.fill(0, m_Grid->GetNumberOfPoints());
872 l2g_t nodes = m_Part.getNodes();
873 for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
874 simpleNodeMovement(i_nodes);
875 }
876 }
877
878 void GridSmoother::operatePostSmoothing()
879 {
880 QVector<bool> top_node(m_Grid->GetNumberOfPoints(), false);
881 QVector<vec3_t> x_new(m_Grid->GetNumberOfPoints());
882 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
883 if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
884 vtkIdType N_pts, *pts;
885 m_Grid->GetCellPoints(id_cell, N_pts, pts);
886 top_node[pts[3]] = true;
887 top_node[pts[4]] = true;
888 top_node[pts[5]] = true;
889 }
890 }
891 for (vtkIdType id_node1 = 0; id_node1 < m_Grid->GetNumberOfPoints(); ++id_node1) {
892 m_Grid->GetPoint(id_node1, x_new[id_node1].data());
893 if (top_node[id_node1]) {
894 double w = m_PostSmoothingStrength;
895 vec3_t x_av(0,0,0);
896 vec3_t x_old;
897 m_Grid->GetPoint(id_node1, x_old.data());
898 QSet<vtkIdType> ids;
899 for (int i = 0; i < m_Part.n2cGSize(id_node1); ++i) {
900 vtkIdType id_cell = m_Part.n2cGG(id_node1, i);
901 if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
902 vtkIdType N_pts, *pts;
903 m_Grid->GetCellPoints(id_cell, N_pts, pts);
904 for (int j = 0; j < 6; ++j) {
905 if (top_node[pts[j]] && (pts[j] != id_node1)) {
906 ids.insert(pts[j]);
907 }
908 }
909 }
910 }
911 foreach (vtkIdType id_node2, ids) {
912 vec3_t x;
913 m_Grid->GetPoint(id_node2, x.data());
914 x_av += x;
915 }
916 if (ids.size() > 0) {
917 x_av *= 1.0/ids.size();
918 x_new[id_node1] = w*x_av + (1-w)*x_old;
919 }
920 }
921 }
922 l2g_t nodes = m_Part.getNodes();
923 for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
924 vtkIdType id_node = m_Part.globalNode(i_nodes);
925 vec3_t x_old;
926 m_Grid->GetPoint(id_node, x_old.data());
927 vec3_t Dx = x_new[id_node] - x_old;
928 correctDx(i_nodes, Dx);
929 moveNode(i_nodes, Dx);
930 }
931 }
932
933 void GridSmoother::operate()
934 {
935 markNodes();
936 computeNormals();
937 if (m_SimpleOperation) {
938 operateSimple();
939 } else if (m_PostOperation) {
940 operatePostSmoothing();
941 } else {
942 operateOptimisation();
943 }
944 }
945
Clone of the official engrid repository