?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
all points are snap points for smoothing with the 1 BC method
[engrid.git] / src / libengrid / surfaceoperation.cpp
blob463873450fd9b221b8e0e87591f278791461b7e4
1 //
2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 // + +
4 // + This file is part of enGrid. +
5 // + +
6 // + Copyright 2008-2012 enGits GmbH +
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 "surfaceoperation.h"
24
25 #include "guimainwindow.h"
26
27 #include <vtkCharArray.h>
28 #include <vtkMath.h>
29 #include <vtkCellArray.h>
30 #include <vtkPolygon.h>
31
32 #include "geometrytools.h"
33 using namespace GeometryTools;
34
35 SurfaceOperation::SurfaceOperation() : Operation()
36 {
37 //default values for determining node types and for smoothing operations
38 getSet("surface meshing", "edge angle to determine fixed vertices", 180, m_EdgeAngle);
39 getSet("surface meshing", "feature angle", 180, m_FeatureAngle);
40 m_FeatureAngle = GeometryTools::deg2rad(m_FeatureAngle);
41 m_EdgeAngle = GeometryTools::deg2rad(m_EdgeAngle);
42 setEdgeAngle(m_EdgeAngle);
43 m_BoundarySmoothing = 1;
44 m_StretchingFactor = 0;
45 m_UniformSnapPoints = false;
46 }
47
48 void SurfaceOperation::operate()
49 {
50
51 }
52
53 ostream& operator<<(ostream &out, stencil_t S)
54 {
55 out << "S.id_cell = " << S.id_cell << " ";
56 out << "S.id_node = " << S.id_node << " ";
57 out << "S.sameBC = " << S.sameBC << " ";
58 out << "S.type = " << S.type_cell << " ";
59 out << "S.p1 = " << S.p1 << " ";
60 out << "S.p2 = " << S.p2 << " ";
61 return(out);
62 }
63
64 stencil_t SurfaceOperation::getStencil(vtkIdType id_cell1, int j1)
65 {
66 stencil_t S;
67 {
68 vtkIdType N_pts, *pts;
69 m_Grid->GetCellPoints(id_cell1, N_pts, pts);
70 S.p1 = pts[j1];
71 S.p2 = pts[0];
72 if (j1 < N_pts - 1) {
73 S.p2 = pts[j1 + 1];
74 }
75 }
76 QSet<vtkIdType> cells_p1;
77 for (int i = 0; i < m_Part.n2cGSize(S.p1); ++i) {
78 vtkIdType id_cell = m_Part.n2cGG(S.p1, i);
79 if (id_cell != id_cell1) {
80 cells_p1.insert(id_cell);
81 }
82 }
83 QSet<vtkIdType> cells_p2;
84 for (int i = 0; i < m_Part.n2cGSize(S.p2); ++i) {
85 vtkIdType id_cell = m_Part.n2cGG(S.p2, i);
86 if (id_cell != id_cell1) {
87 cells_p2.insert(id_cell);
88 }
89 }
90 QSet<vtkIdType> cells = cells_p1.intersect(cells_p2);
91 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
92 S.sameBC = true;
93 S.id_cell.resize(1);
94 S.id_cell[0] = id_cell1;
95 foreach (vtkIdType id_cell, cells) {
96 if (isSurface(id_cell, m_Grid)) {
97 S.id_cell.push_back(id_cell);
98 if (cell_code->GetValue(id_cell) != cell_code->GetValue(id_cell1)) {
99 S.sameBC = false;
100 }
101 }
102 }
103 S.id_node.resize(S.id_cell.size());
104 S.type_cell.resize(S.id_cell.size());
105 for (int i = 0; i < S.id_cell.size(); ++i) {
106 vtkIdType N_pts, *pts;
107 m_Grid->GetCellPoints(S.id_cell[i], N_pts, pts);
108 S.type_cell[i] = m_Grid->GetCellType(S.id_cell[i]);
109 for (int j = 0; j < N_pts; ++j) {
110 if (pts[j] != S.p1 && pts[j] != S.p2) {
111 S.id_node[i] = pts[j];
112 break;
113 }
114 }
115 }
116 return S;
117 }
118
119 int SurfaceOperation::UpdateCurrentMeshDensity()
120 {
121 if ( DebugLevel > 0 ) {
122 cout << "===UpdateMeshDensity START===" << endl;
123 }
124 QVector<vtkIdType> cells;
125 getAllSurfaceCells( cells, m_Grid );
126 EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
127 EG_VTKDCN( vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired" );
128 setGrid( m_Grid );
129 setCells( cells );
130 if ( DebugLevel > 5 ) {
131 cout << "cells.size()=" << cells.size() << endl;
132 }
133 EG_VTKDCN( vtkDoubleArray, node_meshdensity_current, m_Grid, "node_meshdensity_current" );
134 l2g_t nodes = getPartNodes();
135 foreach( vtkIdType node, nodes ) {
136 node_meshdensity_current->SetValue( node, CurrentMeshDensity( node ) );
137 }
138 if ( DebugLevel > 0 ) {
139 cout << "===UpdateMeshDensity END===" << endl;
140 }
141 return( 0 ); ///\todo what for???
142 }
143
144 int SurfaceOperation::UpdatePotentialSnapPoints(bool update_node_types, bool fix_unselected)
145 {
146 setAllSurfaceCells();
147
148 l2g_t nodes = getPartNodes();
149 l2g_t cells = getPartCells();
150
151 m_PotentialSnapPoints.resize(m_Grid->GetNumberOfPoints());
152
153 if (m_UniformSnapPoints) {
154 m_PotentialSnapPoints.resize(m_Grid->GetNumberOfPoints());
155 EG_FORALL_NODES(id_node, m_Grid) {
156 for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
157 m_PotentialSnapPoints[id_node].append(m_Part.n2nGG(id_node, i));
158 }
159 }
160 return 0;
161 }
162
163 //initialize default values
164 EG_VTKDCN( vtkCharArray, node_type, m_Grid, "node_type" );
165 foreach( vtkIdType id_node, nodes ) {
166 if ( update_node_types ) node_type->SetValue( id_node, VTK_SIMPLE_VERTEX );
167 m_PotentialSnapPoints[id_node].clear();
168 }
169
170 int num_edges = 0;
171
172 // loop through edges
173 foreach (vtkIdType id_cell, cells) {
174 vtkIdType *pts, Npts;
175 m_Grid->GetCellPoints( id_cell, Npts, pts );
176 for ( int i = 0; i < Npts; i++ ) {
177
178 QVector <vtkIdType> EdgeCells;
179 int Ncells = getEdgeCells( pts[i], pts[(i+1)%Npts], EdgeCells );
180
181 bool already_visited_edge = false;
182 for(int i_cell=0;i_cell<Ncells;i_cell++) {
183 if(EdgeCells[i_cell]<id_cell) {
184 already_visited_edge = true;
185 break;
186 }
187 }
188
189 if ( already_visited_edge ) continue;//already visited edge
190 num_edges++;
191
192 vtkIdType id_node1 = pts[i];
193 vtkIdType id_node2 = pts[( i+1 )%Npts];
194
195 // determine edge type
196 char edge = getEdgeType( id_node2, id_node1, fix_unselected );
197
198 // determine node type pre-processing (count nb of complex edges if the node is complex, otherwise, just count the nb of edges)
199 if (edge && node_type->GetValue(id_node1) == VTK_SIMPLE_VERTEX)
200 {
201 m_PotentialSnapPoints[id_node1].clear();
202 m_PotentialSnapPoints[id_node1].push_back( id_node2 );
203 if ( update_node_types ) node_type->SetValue( id_node1, edge );
204
205 } else if (( edge && node_type->GetValue(id_node1) == VTK_BOUNDARY_EDGE_VERTEX) ||
206 ( edge && node_type->GetValue(id_node1) == VTK_FEATURE_EDGE_VERTEX) ||
207 (!edge && node_type->GetValue(id_node1) == VTK_SIMPLE_VERTEX))
208 {
209 m_PotentialSnapPoints[id_node1].push_back( id_node2 );
210 if ( node_type->GetValue( id_node1 ) && edge == VTK_BOUNDARY_EDGE_VERTEX ) {
211 if ( update_node_types ) node_type->SetValue( id_node1, VTK_BOUNDARY_EDGE_VERTEX );//VTK_BOUNDARY_EDGE_VERTEX has priority over VTK_FEATURE_EDGE_VERTEX
212 }
213 }
214
215 if ( edge && node_type->GetValue( id_node2 ) == VTK_SIMPLE_VERTEX ) {
216 m_PotentialSnapPoints[id_node2].clear();
217 m_PotentialSnapPoints[id_node2].push_back( id_node1 );
218 if ( update_node_types ) node_type->SetValue( id_node2, edge );
219 }
220 else if (( edge && node_type->GetValue( id_node2 ) == VTK_BOUNDARY_EDGE_VERTEX ) ||
221 ( edge && node_type->GetValue( id_node2 ) == VTK_FEATURE_EDGE_VERTEX ) ||
222 ( !edge && node_type->GetValue( id_node2 ) == VTK_SIMPLE_VERTEX ) ) {
223 m_PotentialSnapPoints[id_node2].push_back( id_node1 );
224 if ( node_type->GetValue( id_node2 ) && edge == VTK_BOUNDARY_EDGE_VERTEX ) {
225 if ( update_node_types ) node_type->SetValue( id_node2, VTK_BOUNDARY_EDGE_VERTEX );//VTK_BOUNDARY_EDGE_VERTEX has priority over VTK_FEATURE_EDGE_VERTEX
226 }
227 }
228 }
229 }
230
231 //cout<<"num_edges="<<num_edges<<endl;
232
233 //-----------------------
234 //determine node type post-processing
235 double CosEdgeAngle = cos(this->m_EdgeAngle);
236 //cout<<"===post-processing==="<<endl;
237 //This time, we loop through nodes
238 foreach( vtkIdType id_node, nodes ) {
239 if ( node_type->GetValue( id_node ) == VTK_FEATURE_EDGE_VERTEX || node_type->GetValue( id_node ) == VTK_BOUNDARY_EDGE_VERTEX ) { //see how many edges; if two, what the angle is
240
241 if ( !this->m_BoundarySmoothing && node_type->GetValue( id_node ) == VTK_BOUNDARY_EDGE_VERTEX ) {
242 if ( update_node_types ) node_type->SetValue( id_node, VTK_FIXED_VERTEX );
243 } else if ( m_PotentialSnapPoints[id_node].size() != 2 ) {
244 if ( update_node_types ) node_type->SetValue( id_node, VTK_FIXED_VERTEX );
245 }
246 else { //check angle between edges
247 double x1[3], x2[3], x3[3], l1[3], l2[3];
248 m_Grid->GetPoint( m_PotentialSnapPoints[id_node][0], x1 );
249 m_Grid->GetPoint( id_node, x2 );
250 m_Grid->GetPoint( m_PotentialSnapPoints[id_node][1], x3 );
251 for ( int k = 0; k < 3; k++ ) {
252 l1[k] = x2[k] - x1[k];
253 l2[k] = x3[k] - x2[k];
254 }
255 if ( vtkMath::Normalize( l1 ) >= 0.0 &&
256 vtkMath::Normalize( l2 ) >= 0.0 &&
257 vtkMath::Dot( l1, l2 ) < CosEdgeAngle ) {
258 if ( update_node_types ) node_type->SetValue( id_node, VTK_FIXED_VERTEX );
259 }
260 }//if along edge
261 }//if edge vertex
262 }
263 //cout<<"m_PotentialSnapPoints.size()="<<m_PotentialSnapPoints.size()<<endl;
264 //cout<<"=== UpdatePotentialSnapPoints END ==="<<endl;
265 return( 0 );
266 }
267
268 char SurfaceOperation::getNodeType(vtkIdType id_node, bool fix_unselected)
269 {
270 l2g_t nodes = getPartNodes();
271 g2l_t _nodes = getPartLocalNodes();
272 l2l_t n2n = getPartN2N();
273
274 // fix all vertices that are part of a volume element or a quad
275 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
276 if (m_Grid->GetCellType(m_Part.n2cGG(id_node, i)) != VTK_TRIANGLE) {
277 return VTK_FIXED_VERTEX;
278 }
279 }
280
281 // initialize default value
282 char type = VTK_SIMPLE_VERTEX;
283
284 // loop through edges around id_node
285
286 QVector <vtkIdType> edges;
287 double cos_edge_angle = cos(this->m_EdgeAngle);
288
289 // safety switch to make sure edge angle fixing is disabled
290 if (cos_edge_angle < -0.99) {
291 cos_edge_angle = -2;
292 }
293
294 foreach (int i_node2, n2n[_nodes[id_node]]) {
295
296 vtkIdType id_node2 = nodes[i_node2];
297
298 // determine edge type
299 char edge = getEdgeType(id_node2, id_node, fix_unselected);
300
301 // determine node type pre-processing
302 // count number of complex edges if the node is complex,
303 // otherwise, just count the number of edges
304 //
305 if (edge != VTK_SIMPLE_VERTEX && type == VTK_SIMPLE_VERTEX) {
306
307 edges.clear();
308 edges.push_back(id_node2);
309 type = edge;
310
311 } else if ((edge != VTK_SIMPLE_VERTEX && type == VTK_BOUNDARY_EDGE_VERTEX) ||
312 (edge != VTK_SIMPLE_VERTEX && type == VTK_FEATURE_EDGE_VERTEX) ||
313 (edge == VTK_SIMPLE_VERTEX && type == VTK_SIMPLE_VERTEX))
314 {
315 edges.push_back(id_node2);
316 if (type != VTK_SIMPLE_VERTEX && edge == VTK_BOUNDARY_EDGE_VERTEX) {
317 // VTK_BOUNDARY_EDGE_VERTEX has priority over VTK_FEATURE_EDGE_VERTEX
318 type = VTK_BOUNDARY_EDGE_VERTEX;
319 }
320 }
321 }
322
323 // determine node type post-processing
324 if (type == VTK_FEATURE_EDGE_VERTEX || type == VTK_BOUNDARY_EDGE_VERTEX) {
325
326 // check how how many edges; if two, what is the angle?
327
328 if (!this->m_BoundarySmoothing && type == VTK_BOUNDARY_EDGE_VERTEX) {
329 type = VTK_FIXED_VERTEX;
330 } else if (edges.size() != 2) {
331 //type = VTK_FIXED_VERTEX;
332 } else {
333
334 // check angle between edges
335 vec3_t x1, x2, x3, l1, l2;
336 m_Grid->GetPoint(edges[0], x1.data());
337 m_Grid->GetPoint(id_node, x2.data());
338 m_Grid->GetPoint(edges[1], x3.data());
339 l1 = x2 - x1;
340 l2 = x3 - x2;
341 if (l1.abs() >= 0.0 && l2.abs() >= 0.0 && l1*l2 < cos_edge_angle) {
342 //type = VTK_FIXED_VERTEX;
343 }
344 } // if along edge
345 } // if edge vertex
346
347 return type;
348 }
349
350 int SurfaceOperation::getEdgeCells(vtkIdType id_node1, vtkIdType id_node2, QVector <vtkIdType> &EdgeCells)
351 {
352 g2l_t _nodes = getPartLocalNodes();
353 l2g_t cells = getPartCells();
354 l2l_t n2c = getPartN2C();
355
356 QSet<vtkIdType> S1;
357 foreach (int i, n2c[_nodes[id_node1]]) {
358 S1.insert(cells[i]);
359 }
360
361 QSet<vtkIdType> S2;
362 foreach( int i, n2c[_nodes[id_node2]] ) {
363 S2.insert(cells[i]);
364 }
365
366 S2.intersect(S1);
367 EdgeCells = Set2Vector(S2, false);
368 return EdgeCells.size();
369 }
370
371 int SurfaceOperation::getEdgeCells( vtkIdType id_node1, vtkIdType id_node2, QSet <vtkIdType> &EdgeCells )
372 {
373 g2l_t _nodes = getPartLocalNodes();
374 l2g_t cells = getPartCells();
375 l2l_t n2c = getPartN2C();
376
377 QSet<vtkIdType> S1;
378 foreach( int i, n2c[_nodes[id_node1]] ) {
379 S1.insert( cells[i] );
380 }
381
382 QSet<vtkIdType> S2;
383 foreach( int i, n2c[_nodes[id_node2]] ) {
384 S2.insert( cells[i] );
385 }
386
387 EdgeCells = S2.intersect( S1 );
388 return EdgeCells.size();
389 }
390
391 char SurfaceOperation::getEdgeType(vtkIdType a_node1, vtkIdType a_node2, bool fix_unselected)
392 {
393 double cos_feature_angle = cos(this->m_FeatureAngle);
394 bool feature_edges_disabled = m_FeatureAngle >= M_PI;
395
396 // compute number of cells around edge [a_node,p2] and put them into neighbour_cells
397 QVector <vtkIdType> neighbour_cells;
398 int numNei = getEdgeCells(a_node1, a_node2, neighbour_cells) - 1;
399
400 // set default value
401 char edge = VTK_SIMPLE_VERTEX;
402
403 if (numNei == 0) {
404 edge = VTK_BOUNDARY_EDGE_VERTEX;
405 } else if (numNei >= 2) {
406 edge = VTK_BOUNDARY_EDGE_VERTEX;
407 } else if (numNei == 1) {
408
409 // check angle between cell1 and cell2 against FeatureAngle
410 if (cosAngle(m_Grid, neighbour_cells[0], neighbour_cells[1] ) <= cos_feature_angle && !feature_edges_disabled) {
411 edge = VTK_FEATURE_EDGE_VERTEX;
412 }
413
414 // check the boundary codes
415 EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
416 int cell_code_0 = cell_code->GetValue( neighbour_cells[0] );
417 int cell_code_1 = cell_code->GetValue( neighbour_cells[1] );
418 if (cell_code_0 != cell_code_1) {
419 edge = VTK_BOUNDARY_EDGE_VERTEX;
420 }
421
422 if (fix_unselected) {
423 if (!m_BoundaryCodes.contains(cell_code_0) || !m_BoundaryCodes.contains(cell_code_1)) {
424 // does not make sense, but should make the points of the edge fixed
425 edge = VTK_FIXED_VERTEX;
426 }
427 }
428 }
429
430 return edge;
431 }
432
433 VertexMeshDensity SurfaceOperation::getVMD( vtkIdType id_node )
434 {
435 g2l_t _nodes = getPartLocalNodes();
436 l2g_t cells = getPartCells();
437 l2l_t n2c = getPartN2C();
438
439 EG_VTKDCN( vtkCharArray, node_type, m_Grid, "node_type" );
440 EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
441
442 VertexMeshDensity VMD;
443 VMD.type = node_type->GetValue( id_node );
444 VMD.density = 0;
445 VMD.CurrentNode = id_node;
446
447 foreach( int i_cell, n2c[_nodes[id_node]] ) {
448 vtkIdType id_cell = cells[i_cell];
449 VMD.BCmap[cell_code->GetValue( id_cell )] = 2;
450 }
451 return( VMD );
452 }
453
454 //////////////////////////////////////////////
455 double SurfaceOperation::currentVertexAvgDist( vtkIdType id_node )
456 {
457 l2g_t nodes = getPartNodes();
458 g2l_t _nodes = getPartLocalNodes();
459 l2l_t n2n = getPartN2N();
460
461 double total_dist = 0;
462 double avg_dist = 0;
463 int N = n2n[_nodes[id_node]].size();
464 vec3_t C;
465 m_Grid->GetPoint( id_node, C.data() );
466 foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
467 vtkIdType id_node_neighbour = nodes[i_node_neighbour];
468 vec3_t M;
469 m_Grid->GetPoint( id_node_neighbour, M.data() );
470 total_dist += ( M - C ).abs();
471 }
472 avg_dist = total_dist / ( double )N;
473 return( avg_dist );
474 }
475
476 double SurfaceOperation::CurrentMeshDensity( vtkIdType id_node )
477 {
478 return 1.0 / currentVertexAvgDist( id_node );
479 }
480
481 ///\todo change meshdensity fields to edgelength fields since this is what is mostly used?
482
483 /// desired edge length for id_node
484 double SurfaceOperation::desiredEdgeLength( vtkIdType id_node )
485 {
486 EG_VTKDCN( vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired" );
487 return( 1.0 / characteristic_length_desired->GetValue( id_node ) );
488 }
489
490 /// mean desired edge length for id_cell
491 double SurfaceOperation::meanDesiredEdgeLength( vtkIdType id_cell )
492 {
493 vtkIdType num_pts, *pts;
494 m_Grid->GetCellPoints( id_cell, num_pts, pts );
495 int total = 0;
496 for ( int i = 0; i < num_pts; i++ ) {
497 total += desiredEdgeLength( pts[i] );
498 }
499 return total / ( double )num_pts;
500 }
501
502 QVector <vtkIdType> SurfaceOperation::getPotentialSnapPoints( vtkIdType id_node )
503 {
504 if ((id_node < 0) || (id_node >= m_PotentialSnapPoints.size())) {
505 // UpdatePotentialSnapPoints should probably be called before using this function.
506 EG_BUG;
507 }
508 return m_PotentialSnapPoints[id_node];
509 }
510
511 bool SurfaceOperation::isCell(vtkIdType id_node1, vtkIdType id_node2, vtkIdType id_node3)
512 {
513 QVector <vtkIdType> EdgeCells_12;
514 QVector <vtkIdType> EdgeCells_13;
515 QVector <vtkIdType> inter;
516
517 getEdgeCells( id_node1, id_node2, EdgeCells_12 );
518 getEdgeCells( id_node1, id_node3, EdgeCells_13 );
519 qcontIntersection( EdgeCells_12, EdgeCells_13, inter );
520 if(inter.size()>1) {
521 qWarning()<<"(id_node1, id_node2, id_node3)="<<"("<<id_node1<<", "<<id_node2<<", "<<id_node3<<")";
522 qWarning()<<"EdgeCells_12="<<EdgeCells_12;
523 qWarning()<<"EdgeCells_13="<<EdgeCells_13;
524 qWarning()<<"inter="<<inter;
525 writeGrid(m_Grid, "abort");
526 EG_BUG;// multiple cells in the same place
527 }
528 if(DebugLevel>100) {
529 qDebug()<<"(id_node1, id_node2, id_node3)="<<"("<<id_node1<<", "<<id_node2<<", "<<id_node3<<")";
530 qDebug()<<"EdgeCells_12="<<EdgeCells_12;
531 qDebug()<<"EdgeCells_13="<<EdgeCells_13;
532 qDebug()<<"inter="<<inter;
533 }
534 return(inter.size()>0);
535 }
536
537 void SurfaceOperation::computeNormals()
538 {
539 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
540 m_NodeNormal.fill(vec3_t(0,0,0), m_Grid->GetNumberOfPoints());
541 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
542 if (m_Part.localNode(id_node) != -1) {
543 QSet<int> bcs;
544 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
545 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
546 if (isSurface(id_cell, m_Grid)) {
547 int bc = cell_code->GetValue(id_cell);
548 if (m_BoundaryCodes.contains(bc)) {
549 bcs.insert(bc);
550 }
551 }
552 }
553 int num_bcs = bcs.size();
554 QVector<vec3_t> normal(num_bcs, vec3_t(0,0,0));
555 QMap<int,int> bcmap;
556 int i_bc = 0;
557 foreach (int bc, bcs) {
558 bcmap[bc] = i_bc;
559 ++i_bc;
560 }
561 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
562 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
563 if (isSurface(id_cell, m_Grid)) {
564 int bc = cell_code->GetValue(id_cell);
565 if (m_BoundaryCodes.contains(bc)) {
566 vtkIdType N_pts, *pts;
567 m_Grid->GetCellPoints(id_cell, N_pts, pts);
568 vec3_t a, b, c;
569 for (int j = 0; j < N_pts; ++j) {
570 if (pts[j] == id_node) {
571 m_Grid->GetPoint(pts[j], a.data());
572 if (j > 0) {
573 m_Grid->GetPoint(pts[j-1], b.data());
574 } else {
575 m_Grid->GetPoint(pts[N_pts-1], b.data());
576 }
577 if (j < N_pts - 1) {
578 m_Grid->GetPoint(pts[j+1], c.data());
579 } else {
580 m_Grid->GetPoint(pts[0], c.data());
581 }
582 }
583 }
584 vec3_t u = b - a;
585 vec3_t v = c - a;
586 double alpha = GeometryTools::angle(u, v);
587 vec3_t n = u.cross(v);
588 n.normalise();
589 if (checkVector(n)) {
590 normal[bcmap[bc]] += alpha*n;
591 }
592 }
593 }
594 }
595 for (int i = 0; i < num_bcs; ++i) {
596 normal[i].normalise();
597 }
598 if (num_bcs > 0) {
599 if (num_bcs > 1) {
600 if (num_bcs == 3) {
601 for (int i = 0; i < num_bcs; ++i) {
602 for (int j = i + 1; j < num_bcs; ++j) {
603 vec3_t n = normal[i] + normal[j];
604 n.normalise();
605 m_NodeNormal[id_node] += n;
606 }
607 }
608 } else {
609 for (int i = 0; i < num_bcs; ++i) {
610 m_NodeNormal[id_node] += normal[i];
611 }
612 }
613 } else {
614 m_NodeNormal[id_node] = normal[0];
615 }
616 m_NodeNormal[id_node].normalise();
617 }
618 }
619 }
620 }
621
622 double SurfaceOperation::normalIrregularity(vtkIdType id_node)
623 {
624 double nirr = 0;
625 QVector<vec3_t> nc(m_Part.n2cGSize(id_node));
626 for (int i = 0; i < nc.size(); ++i) {
627 nc[i] = GeometryTools::cellNormal(m_Grid, m_Part.n2cGG(id_node, i));
628 nc[i].normalise();
629 }
630 for (int i = 0; i < nc.size(); ++i) {
631 for (int j = i + 1; j < nc.size(); ++j) {
632 nirr += 1.0 - nc[i]*nc[j];
633 }
634 }
635 return nirr;
636 }
Clone of the official engrid repository