2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4 // + This file is part of enGrid. +
6 // + Copyright 2008-2010 enGits GmbH +
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 m_NumBoundaryCorrections
= 50;
36 getSet("boundary layer", "number of smoothing sub-iterations", 5, m_NumIterations
);
37 getSet("boundary layer", "prism layer clearance", 0.2, m_LayerClearance
);
38 getSet("boundary layer", "use strict prism checking", false, m_StrictPrismChecking
);
39 getSet("boundary layer", "number of normal vector relax iterations", 10, m_NumNormalRelaxations
);
40 getSet("boundary layer", "number of layer height relax iterations", 3, m_NumHeightRelaxations
);
41 getSet("boundary layer", "radar angle", 30, m_RadarAngle
);
42 getSet("boundary layer", "maximal layer height in gaps", 0.2, m_MaxHeightInGaps
);
44 //m_CritAngle = GeometryTools::deg2rad(m_CritAngle);
47 void GridSmoother::markNodes()
49 m_NodeMarked
.fill(false,m_Grid
->GetNumberOfPoints());
50 QVector
<bool> new_mark(m_Grid
->GetNumberOfPoints());
51 for (int i_iterations
= 0; i_iterations
< 4; ++i_iterations
) {
52 qCopy(m_NodeMarked
.begin(),m_NodeMarked
.end(),new_mark
.begin());
53 for (vtkIdType id_cell
= 0; id_cell
< m_Grid
->GetNumberOfCells(); ++id_cell
) {
54 bool mark_cell
= false;
55 vtkIdType type_cell
, N_pts
, *pts
;
56 type_cell
= m_Grid
->GetCellType(id_cell
);
57 m_Grid
->GetCellPoints(id_cell
, N_pts
, pts
);
58 if (type_cell
== VTK_WEDGE
) {
61 for (int i_pts
= 0; i_pts
< N_pts
; ++i_pts
) {
62 if (m_NodeMarked
[pts
[i_pts
]]) {
68 for (int i_pts
= 0; i_pts
< N_pts
; ++i_pts
) {
69 new_mark
[pts
[i_pts
]] = true;
73 qCopy(new_mark
.begin(),new_mark
.end(),m_NodeMarked
.begin());
76 QVector
<vtkIdType
> nodes
= m_Part
.getNodes();
77 foreach (vtkIdType id_node
, nodes
) {
78 if (id_node
< 0) EG_BUG
;
79 if (id_node
> m_Grid
->GetNumberOfPoints()) EG_BUG
;
80 if (m_NodeMarked
[id_node
]) {
86 bool GridSmoother::noCollision(vtkIdType id_node
)
89 vtkIdType id_foot
= m_IdFoot
[id_node
];
91 QSet
<vtkIdType
> front_neigh
;
92 QSet
<vtkIdType
> own_prisms
;
93 for (int i
= 0; i
< m_Part
.n2nGSize(id_node
); ++i
) {
94 vtkIdType id_neigh
= m_Part
.n2nGG(id_node
, i
);
95 if (m_IdFoot
[id_neigh
] != -1) {
96 front_neigh
.insert(id_neigh
);
98 QSet
<vtkIdType
> neigh_prisms
;
99 for (int j
= 0; j
< m_Part
.n2cGSize(id_neigh
); ++j
) {
100 vtkIdType id_cell
= m_Part
.n2cGG(id_neigh
, j
);
101 if (m_Grid
->GetCellType(id_cell
) == VTK_WEDGE
) {
102 neigh_prisms
.insert(id_cell
);
105 if (neigh_prisms
.size() == 0) {
106 for (int j
= 0; j
< m_Part
.n2nGSize(id_neigh
); ++j
) {
107 vtkIdType id_next_neigh
= m_Part
.n2nGG(id_neigh
, j
);
108 if (m_IdFoot
[id_next_neigh
] != -1) {
109 front_neigh
.insert(id_next_neigh
);
115 for (int i
= 0; i
< m_Part
.n2cGSize(id_node
); ++i
) {
116 vtkIdType id_cell
= m_Part
.n2cGG(id_node
, i
);
117 if (m_Grid
->GetCellType(id_cell
) == VTK_WEDGE
) {
118 own_prisms
.insert(id_cell
);
121 foreach (vtkIdType id_neigh
, front_neigh
) {
122 QSet
<vtkIdType
> neigh_prisms
;
123 for (int i
= 0; i
< m_Part
.n2cGSize(id_neigh
); ++i
) {
124 vtkIdType id_cell
= m_Part
.n2cGG(id_neigh
, i
);
125 if (m_Grid
->GetCellType(id_cell
) == VTK_WEDGE
) {
126 neigh_prisms
.insert(id_cell
);
129 QSet
<vtkIdType
> all_prisms
= neigh_prisms
+ own_prisms
;
130 if (all_prisms
.size() == neigh_prisms
.size() + own_prisms
.size()) { // here is a collision candidate!
132 foreach (vtkIdType id_cell
, own_prisms
) {
133 vtkIdType N_pts
, *pts
;
134 m_Grid
->GetCellPoints(id_cell
, N_pts
, pts
);
136 m_Grid
->GetPoint(pts
[0], a
.data());
137 m_Grid
->GetPoint(pts
[2], b
.data());
138 m_Grid
->GetPoint(pts
[1], c
.data());
141 vec3_t n
= u
.cross(v
);
143 double alpha
= GeometryTools::angle(u
, v
);
148 foreach (vtkIdType id_cell
, neigh_prisms
) {
149 vtkIdType N_pts
, *pts
;
150 m_Grid
->GetCellPoints(id_cell
, N_pts
, pts
);
152 m_Grid
->GetPoint(pts
[0], a
.data());
153 m_Grid
->GetPoint(pts
[2], b
.data());
154 m_Grid
->GetPoint(pts
[1], c
.data());
157 vec3_t n
= u
.cross(v
);
159 double alpha
= GeometryTools::angle(u
, v
);
164 vec3_t x1
, x2
, xf1
, xf2
;
165 m_Grid
->GetPoint(id_node
, x1
.data());
166 m_Grid
->GetPoint(id_neigh
, x2
.data());
167 m_Grid
->GetPoint(m_IdFoot
[id_node
], xf1
.data());
168 m_Grid
->GetPoint(m_IdFoot
[id_neigh
], xf2
.data());
169 double l1
= (x1
-xf1
)*n1
;
170 double l2
= (x2
-xf2
)*n2
;
171 if (n1
*(x2
-x1
) < m_LayerClearance
*l1
) {
174 if (n2
*(x1
-x2
) < m_LayerClearance
*l2
) {
182 m_CollisionDetected
= true;
187 bool GridSmoother::setNewPosition(vtkIdType id_node
, vec3_t x_new
)
189 using namespace GeometryTools
;
192 m_Grid
->GetPoint(id_node
, x_old
.data());
193 m_Grid
->GetPoints()->SetPoint(id_node
, x_new
.data());
194 bool move
= noCollision(id_node
);
199 l2g_t cells
= m_Part
.getCells();
200 l2l_t n2c
= m_Part
.getN2C();
202 foreach (int i_cells
, n2c
[id_node
]) {
203 vtkIdType id_cell
= cells
[i_cells
];
204 vtkIdType type_cell
= m_Grid
->GetCellType(id_cell
);
205 if (type_cell
== VTK_TETRA
) {
206 if (GeometryTools::cellVA(m_Grid
, id_cell
) < 0) {
211 if (type_cell
== VTK_WEDGE
&& m_StrictPrismChecking
) {
212 vtkIdType N_pts
, *pts
;
214 m_Grid
->GetCellPoints(id_cell
, N_pts
, pts
);
216 for (int i
= 0; i
< 4; ++i
) { // variation
218 for (int j
= 0; j
< 3; ++j
) { // tetrahedron
219 for (int k
= 0; k
< 4; ++k
) { // node
220 m_Grid
->GetPoint(pts
[E
.priTet(i
,j
,k
)], xtet
[k
].data());
222 double V
= GeometryTools::tetraVol(xtet
[0], xtet
[1], xtet
[2], xtet
[3]);
238 m_Grid
->GetPoints()->SetPoint(id_node
, x_old
.data());
243 void GridSmoother::correctDx(int i_nodes
, vec3_t
&Dx
)
245 l2g_t nodes
= m_Part
.getNodes();
246 l2l_t n2c
= m_Part
.getN2C();
247 for (int i_boundary_correction
= 0; i_boundary_correction
< m_NumBoundaryCorrections
; ++i_boundary_correction
) {
248 foreach (vtkIdType id_cell
, n2c
[i_nodes
]) {
249 if (isSurface(id_cell
, m_Grid
)) {
250 double A
= GeometryTools::cellVA(m_Grid
, id_cell
);
252 vec3_t n
= GeometryTools::cellNormal(m_Grid
, id_cell
);
257 if (m_Grid
->GetCellType(id_cell
) == VTK_WEDGE
) {
258 vtkIdType N_pts
, *pts
;
259 m_Grid
->GetCellPoints(id_cell
, N_pts
, pts
);
260 vtkIdType id_surf_node
= -1;
261 if (pts
[3] == nodes
[i_nodes
]) id_surf_node
= pts
[0];
262 if (pts
[4] == nodes
[i_nodes
]) id_surf_node
= pts
[1];
263 if (pts
[5] == nodes
[i_nodes
]) id_surf_node
= pts
[2];
264 if (id_surf_node
!= -1) {
266 m_Grid
->GetPoint(pts
[0],x0
.data());
267 m_Grid
->GetPoint(pts
[1],x1
.data());
268 m_Grid
->GetPoint(pts
[2],x2
.data());
272 double L
= (a
.abs()+b
.abs()+c
.abs())/3.0;
273 vec3_t n
= b
.cross(a
);
276 m_Grid
->GetPoint(nodes
[i_nodes
],x_old
.data());
277 vec3_t x_new
= x_old
+ Dx
- x0
;
278 if ( (n
*x_new
) <= 0 ) {
279 x_new
-= (x_new
*n
)*n
;
291 bool GridSmoother::moveNode(int i_nodes
, vec3_t
&Dx
)
293 m_CollisionDetected
= false;
294 l2g_t nodes
= m_Part
.getNodes();
295 vtkIdType id_node
= nodes
[i_nodes
];
297 m_Grid
->GetPoint(id_node
, x_old
.data());
299 for (int i_relaxation
= 0; i_relaxation
< m_NumRelaxations
; ++i_relaxation
) {
300 if (setNewPosition(id_node
, x_old
+ Dx
)) {
309 void GridSmoother::computeNormals()
311 EG_VTKDCC(vtkIntArray
, cell_code
, m_Grid
, "cell_code");
312 m_NodeNormal
.fill(vec3_t(0,0,0), m_Grid
->GetNumberOfPoints());
313 for (vtkIdType id_node
= 0; id_node
< m_Grid
->GetNumberOfPoints(); ++id_node
) {
315 for (int i
= 0; i
< m_Part
.n2cGSize(id_node
); ++i
) {
316 vtkIdType id_cell
= m_Part
.n2cGG(id_node
, i
);
317 if (isSurface(id_cell
, m_Grid
)) {
318 int bc
= cell_code
->GetValue(id_cell
);
319 if (m_BoundaryCodes
.contains(bc
)) {
324 int num_bcs
= bcs
.size();
325 QVector
<vec3_t
> normal(num_bcs
, vec3_t(0,0,0));
328 foreach (int bc
, bcs
) {
332 for (int i
= 0; i
< m_Part
.n2cGSize(id_node
); ++i
) {
333 vtkIdType id_cell
= m_Part
.n2cGG(id_node
, i
);
334 if (isSurface(id_cell
, m_Grid
)) {
335 int bc
= cell_code
->GetValue(id_cell
);
336 if (m_BoundaryCodes
.contains(bc
)) {
337 vtkIdType N_pts
, *pts
;
338 m_Grid
->GetCellPoints(id_cell
, N_pts
, pts
);
340 for (int j
= 0; j
< N_pts
; ++j
) {
341 if (pts
[j
] == id_node
) {
342 m_Grid
->GetPoint(pts
[j
], a
.data());
344 m_Grid
->GetPoint(pts
[j
-1], b
.data());
346 m_Grid
->GetPoint(pts
[N_pts
-1], b
.data());
349 m_Grid
->GetPoint(pts
[j
+1], c
.data());
351 m_Grid
->GetPoint(pts
[0], c
.data());
357 double alpha
= GeometryTools::angle(u
, v
);
358 vec3_t n
= u
.cross(v
);
360 normal
[bcmap
[bc
]] += alpha
*n
;
364 for (int i
= 0; i
< num_bcs
; ++i
) {
365 normal
[i
].normalise();
370 for (int i
= 0; i
< num_bcs
; ++i
) {
371 for (int j
= i
+ 1; j
< num_bcs
; ++j
) {
372 vec3_t n
= normal
[i
] + normal
[j
];
374 m_NodeNormal
[id_node
] += n
;
378 for (int i
= 0; i
< num_bcs
; ++i
) {
379 m_NodeNormal
[id_node
] += normal
[i
];
383 m_NodeNormal
[id_node
] = normal
[0];
385 m_NodeNormal
[id_node
].normalise();
389 m_IsSharpNode.fill(false, m_Grid->GetNumberOfPoints());
390 for (int i = 0; i < m_Part.getNumberOfCells(); ++i) {
391 vtkIdType id_cell1 = m_Part.globalCell(i);
392 if (isSurface(id_cell1, m_Grid)) {
393 if (m_BoundaryCodes.contains(cell_code->GetValue(id_cell1))) {
394 vec3_t n1 = cellNormal(m_Grid, id_cell1);
396 vtkIdType N_pts, *pts;
397 m_Grid->GetCellPoints(id_cell1, N_pts, pts);
398 for (int j = 0; j < m_Part.c2cLSize(i); ++j) {
399 vtkIdType id_cell2 = m_Part.c2cLG(i, j);
400 if (m_BoundaryCodes.contains(cell_code->GetValue(id_cell2))) {
401 vec3_t n2 = cellNormal(m_Grid, id_cell2);
403 if (GeometryTools::angle(n1, n2) > m_CritAngle) {
404 vtkIdType id_node1 = pts[j];
405 vtkIdType id_node2 = pts[0];
406 if (j < m_Part.c2cLSize(i) - 1) {
409 vec3_t x_node1, x_node2;
410 m_Grid->GetPoint(id_node1, x_node1.data());
411 m_Grid->GetPoint(id_node2, x_node2.data());
412 vec3_t x_corner = 0.5*(x_node1 + x_node2);
413 vec3_t x_cell1 = cellCentre(m_Grid, id_cell1);
414 vec3_t x_cell2 = cellCentre(m_Grid, id_cell2);
415 vec3_t v_cell1 = x_cell2 - x_cell1;
416 vec3_t v_cell2 = x_cell1 - x_cell2;
417 if (n1*v_cell1 > 0 || n2*v_cell2 > 0) {
418 m_IsSharpNode[id_node1] = true;
419 m_IsSharpNode[id_node2] = true;
429 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
430 if (m_IsSharpNode[id_node]) {
431 QVector<vec3_t> normals;
432 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
433 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
434 if (isSurface(id_cell, m_Grid)) {
435 if (m_BoundaryCodes.contains(cell_code->GetValue(id_cell))) {
436 vec3_t n = cellNormal(m_Grid, id_cell);
438 normals.push_back(n);
443 foreach (vec3_t n1, normals) {
444 foreach (vec3_t n2, normals) {
445 if (GeometryTools::angle(n1, n2) > m_CritAngle) {
454 relaxNormalVectors();
458 void GridSmoother::relaxNormalVectors()
460 m_SurfNode
.fill(false, m_Grid
->GetNumberOfPoints());
461 for (vtkIdType id_node
= 0; id_node
< m_Grid
->GetNumberOfPoints(); ++id_node
) {
462 for (int i
= 0; i
< m_Part
.n2cGSize(id_node
); ++i
) {
463 if (isSurface(m_Part
.n2cGG(id_node
,i
), m_Grid
)) {
464 m_SurfNode
[id_node
] = true;
469 EG_VTKDCC(vtkIntArray
, bc
, m_Grid
, "cell_code");
470 QVector
<QSet
<int> > n2bc(m_Grid
->GetNumberOfPoints());
471 for (vtkIdType id_cell
= 0; id_cell
< m_Grid
->GetNumberOfCells(); ++id_cell
) {
472 if (isSurface(id_cell
, m_Grid
)) {
473 vtkIdType N_pts
, *pts
;
474 m_Grid
->GetCellPoints(id_cell
, N_pts
, pts
);
475 for (int i
= 0; i
< N_pts
; ++i
) {
476 if (m_SurfNode
[pts
[i
]]) {
477 n2bc
[pts
[i
]].insert(bc
->GetValue(id_cell
));
482 QVector
<int> num_bcs(m_Grid
->GetNumberOfPoints(),0);
483 for (vtkIdType id_node
= 0; id_node
< m_Grid
->GetNumberOfPoints(); ++id_node
) {
484 if (m_SurfNode
[id_node
]) {
485 QList
<vtkIdType
> snap_points
;
486 foreach (int bc
, n2bc
[id_node
]) {
487 if (m_BoundaryCodes
.contains(bc
)) {
493 for (vtkIdType id_node
= 0; id_node
< m_Grid
->GetNumberOfPoints(); ++id_node
) {
494 if (num_bcs
[id_node
] == 0) {
495 m_SurfNode
[id_node
] = false;
496 n2bc
[id_node
].clear();
499 for (int iter
= 0; iter
< m_NumNormalRelaxations
; ++iter
) {
500 QVector
<vec3_t
> n_new(m_NodeNormal
.size(), vec3_t(0,0,0));
501 for (vtkIdType id_node
= 0; id_node
< m_Grid
->GetNumberOfPoints(); ++id_node
) {
502 if (m_SurfNode
[id_node
] ) {
503 QList
<vtkIdType
> snap_points
;
504 for (int i
= 0; i
< m_Part
.n2nGSize(id_node
); ++i
) {
505 vtkIdType id_neigh
= m_Part
.n2nGG(id_node
,i
);
506 if (num_bcs
[id_node
] <= num_bcs
[id_neigh
]) {
507 if (!m_SurfNode
[id_neigh
]) {
508 cout
<< id_node
<< ',' << m_Part
.n2nGG(id_node
,i
) << ',' << num_bcs
[id_node
] << ',' << num_bcs
[id_neigh
] << endl
;
511 snap_points
.append(id_neigh
);
514 if (snap_points
.size() > 0) {
515 n_new
[id_node
] = vec3_t(0,0,0);
516 foreach (vtkIdType id_snap
, snap_points
) {
517 n_new
[id_node
] += m_NodeNormal
[id_snap
];
519 n_new
[id_node
].normalise();
521 n_new
[id_node
] = m_NodeNormal
[id_node
];
523 if (n_new
[id_node
].abs() < 0.1) {
524 cout
<< id_node
<< ',' << n_new
[id_node
] << ',' << num_bcs
[id_node
] << ',' << n2bc
[id_node
] << endl
;
529 m_NodeNormal
= n_new
;
530 correctNormalVectors();
533 num
= "normals_" + num
;
537 //writeDebugFile("normals");
541 void GridSmoother::computeHeights()
543 // first pass (intial heigh)
544 EG_VTKDCN(vtkDoubleArray
, cl
, m_Grid
, "node_meshdensity_desired" );
545 m_Height
.fill(0, m_Grid
->GetNumberOfPoints());
546 for (vtkIdType id_node
= 0; id_node
< m_Grid
->GetNumberOfPoints(); ++id_node
) {
547 if (m_SurfNode
[id_node
]) {
548 m_Height
[id_node
] = cl
->GetValue(id_node
);
549 // if undefined: compute height from surrounding edges
550 if (m_Height
[id_node
] < 1e-99) {
551 m_Height
[id_node
] = 0;
554 m_Grid
->GetPoint(id_node
, x
.data());
555 for (int i
= 0; i
< m_Part
.n2nGSize(id_node
); ++i
) {
556 vtkIdType id_neigh_node
= m_Part
.n2nGG(id_node
, i
);
557 if (m_SurfNode
[id_neigh_node
]) {
560 m_Grid
->GetPoint(id_neigh_node
, xn
.data());
561 m_Height
[id_node
] += (x
- xn
).abs();
567 m_Height
[id_node
] /= N
;
572 // second pass (correct with absolute height if required)
573 for (vtkIdType id_node
= 0; id_node
< m_Grid
->GetNumberOfPoints(); ++id_node
) {
574 if (m_SurfNode
[id_node
]) {
575 m_Height
[id_node
] = m_Blending
*m_AbsoluteHeight
+ (1.0-m_Blending
)*m_RelativeHeight
*m_Height
[id_node
];
580 QList
<vtkIdType
> surf_nodes
;
581 for (vtkIdType id_node
= 0; id_node
< m_Grid
->GetNumberOfPoints(); ++id_node
) {
582 if (m_SurfNode
[id_node
]) {
583 surf_nodes
.append(id_node
);
586 foreach (vtkIdType id_node1
, surf_nodes
) {
587 foreach (vtkIdType id_node2
, surf_nodes
) {
588 if (id_node1
!= id_node2
) {
589 const vec3_t
& n1
= m_NodeNormal
[id_node1
];
590 const vec3_t
& n2
= m_NodeNormal
[id_node2
];
592 m_Grid
->GetPoint(id_node1
, x1
.data());
593 m_Grid
->GetPoint(id_node2
, x2
.data());
598 double alpha
= 180.0/M_PI
*acos(a
/b
);
599 if (alpha
< m_RadarAngle
) {
600 m_Height
[id_node1
] = min(m_Height
[id_node1
], m_MaxHeightInGaps
*a
);
607 // fourth pass (smoothing)
608 for (int iter
= 0; iter
< m_NumHeightRelaxations
; ++iter
) {
609 QVector
<double> h_new
= m_Height
;
610 for (vtkIdType id_node
= 0; id_node
< m_Grid
->GetNumberOfPoints(); ++id_node
) {
611 if (m_SurfNode
[id_node
]) {
614 for (int i
= 0; i
< m_Part
.n2nGSize(id_node
); ++i
) {
615 vtkIdType id_neigh
= m_Part
.n2nGG(id_node
,i
);
616 if (m_SurfNode
[id_neigh
]) {
617 h_new
[id_node
] += m_Height
[id_neigh
];
631 void GridSmoother::correctNormalVectors()
633 EG_VTKDCC(vtkIntArray
, bc
, m_Grid
, "cell_code");
634 for (vtkIdType id_node
= 0; id_node
< m_Grid
->GetNumberOfPoints(); ++id_node
) {
635 if (m_NodeNormal
[id_node
].abs() > 0.1) {
636 for (int iter
= 0; iter
< m_NumBoundaryCorrections
; ++iter
) {
637 for (int i
= 0; i
< m_Part
.n2cGSize(id_node
); ++i
) {
638 vtkIdType id_cell
= m_Part
.n2cGG(id_node
,i
);
639 if (isSurface(id_cell
, m_Grid
)) {
640 if (!m_BoundaryCodes
.contains(bc
->GetValue(id_cell
))) {
641 vec3_t v
= m_NodeNormal
[id_node
];
642 vec3_t n
= GeometryTools::cellNormal(m_Grid
, id_cell
);
644 v
-= (n
*m_NodeNormal
[id_node
])*n
;
646 m_NodeNormal
[id_node
] = v
;
655 void GridSmoother::computeFeet()
657 m_IdFoot
.fill(-1, m_Grid
->GetNumberOfPoints());
658 for (vtkIdType id_cell
= 0; id_cell
< m_Grid
->GetNumberOfCells(); ++id_cell
) {
659 if (m_Grid
->GetCellType(id_cell
) == VTK_WEDGE
) {
660 vtkIdType N_pts
, *pts
;
661 m_Grid
->GetCellPoints(id_cell
, N_pts
, pts
);
662 m_IdFoot
[pts
[3]] = pts
[0];
663 m_IdFoot
[pts
[4]] = pts
[1];
664 m_IdFoot
[pts
[5]] = pts
[2];
665 m_NodeMarked
[pts
[0]] = false;
666 m_NodeMarked
[pts
[1]] = false;
667 m_NodeMarked
[pts
[2]] = false;
672 void GridSmoother::simpleNodeMovement(int i_nodes
)
674 vtkIdType id_node
= m_Part
.globalNode(i_nodes
);
675 vtkIdType id_foot
= m_IdFoot
[id_node
];
676 vec3_t
x_new(0,0,0), x_node
;
677 m_Grid
->GetPoint(id_node
, x_node
.data());
679 vec3_t
x_surf(0,0,0);
680 if (m_SurfNode
[id_node
]) {
682 for (int i
= 0; i
< m_Part
.n2nLSize(i_nodes
); ++i
) {
683 vtkIdType id_neigh
= m_Part
.n2nLG(i_nodes
,i
);
684 if (m_SurfNode
[id_neigh
]) {
686 m_Grid
->GetPoint(id_neigh
, x_neigh
.data());
699 vec3_t x_foot
, x_node
;
700 m_Grid
->GetPoint(id_foot
, x_foot
.data());
701 //double H = m_Blending*m_AbsoluteHeight + (1.0-m_Blending)*m_RelativeHeight*m_Height[id_foot];
702 double H
= m_Height
[id_foot
];
703 x_new
= x_foot
+ H
*m_NodeNormal
[id_foot
];
705 if (m_SurfNode
[id_node
]) {
708 if (m_Part
.n2nLSize(i_nodes
) == 0) {
711 for (int i
= 0; i
< m_Part
.n2nLSize(i_nodes
); ++i
) {
712 vtkIdType id_neigh
= m_Part
.n2nLG(i_nodes
,i
);
714 m_Grid
->GetPoint(id_neigh
, x_neigh
.data());
717 x_new
*= 1.0/m_Part
.n2nLSize(i_nodes
);
720 vec3_t Dx
= x_new
- x_node
;
721 correctDx(i_nodes
, Dx
);
722 moveNode(i_nodes
, Dx
);
725 void GridSmoother::operate()
731 l2g_t nodes
= m_Part
.getNodes();
732 for (int i_nodes
= 0; i_nodes
< nodes
.size(); ++i_nodes
) {
733 if (m_NodeMarked
[nodes
[i_nodes
]]) {
734 simpleNodeMovement(i_nodes
);
739 void GridSmoother::writeDebugFile(QString file_name
)
741 QVector
<vtkIdType
> bcells
;
742 getSurfaceCells(m_BoundaryCodes
, bcells
, m_Grid
);
743 MeshPartition
bpart(m_Grid
);
744 bpart
.setCells(bcells
);
746 QVector<vtkIdType> foot2field(bpart.getNumberOfNodes(), -1);
747 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
748 vtkIdType id_foot = m_IdFoot[id_node];
750 if (bpart.localNode(id_foot) == -1) {
753 foot2field[bpart.localNode(id_foot)] = id_node;
757 file_name
= GuiMainWindow::pointer()->getCwd() + "/" + file_name
+ ".vtk";
758 QFile
file(file_name
);
759 file
.open(QIODevice::WriteOnly
);
760 QTextStream
f(&file
);
761 f
<< "# vtk DataFile Version 2.0\n";
762 f
<< "m_NodeNormal\n";
764 f
<< "DATASET UNSTRUCTURED_GRID\n";
765 f
<< "POINTS " << bpart
.getNumberOfNodes() << " float\n";
766 for (int i
= 0; i
< bpart
.getNumberOfNodes(); ++i
) {
768 vtkIdType id_node
= bpart
.globalNode(i
);
769 m_Grid
->GetPoint(id_node
, x
.data());
770 f
<< x
[0] << " " << x
[1] << " " << x
[2] << "\n";
772 f
<< "CELLS " << bpart
.getNumberOfCells() << " " << 4*bpart
.getNumberOfCells() << "\n";
773 for (int i
= 0; i
< bpart
.getNumberOfCells(); ++ i
) {
774 vtkIdType id_cell
= bpart
.globalCell(i
);
775 vtkIdType N_pts
, *pts
;
776 if (m_Grid
->GetCellType(id_cell
) != VTK_TRIANGLE
) {
779 m_Grid
->GetCellPoints(id_cell
, N_pts
, pts
);
781 for (int j
= 0; j
< 3; ++j
) {
782 nds
[j
] = bpart
.localNode(pts
[j
]);
784 f
<< "3 " << nds
[0] << " " << nds
[1] << " " << nds
[2] << "\n";
787 f
<< "CELL_TYPES " << bpart
.getNumberOfCells() << "\n";
788 for (int i
= 0; i
< bpart
.getNumberOfCells(); ++ i
) {
789 f
<< VTK_TRIANGLE
<< "\n";
791 f
<< "POINT_DATA " << bpart
.getNumberOfNodes() << "\n";
792 f
<< "VECTORS N float\n";
794 for (int i
= 0; i
< bpart
.getNumberOfNodes(); ++i
) {
795 vtkIdType id_node
= bpart
.globalNode(i
);
796 f
<< m_NodeNormal
[id_node
][0] << " " << m_NodeNormal
[id_node
][1] << " " << m_NodeNormal
[id_node
][2] << "\n";