?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
attempt to keep existing node types
[engrid.git] / src / libengrid / surfaceoperation.cpp
blob70be49c1a2a7d41cb0fcda406f0ba6dbc020467e
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, EG_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) == EG_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) == EG_BOUNDARY_EDGE_VERTEX) ||
206 ( edge && node_type->GetValue(id_node1) == EG_FEATURE_EDGE_VERTEX) ||
207 (!edge && node_type->GetValue(id_node1) == EG_SIMPLE_VERTEX))
208 {
209 m_PotentialSnapPoints[id_node1].push_back( id_node2 );
210 if ( node_type->GetValue( id_node1 ) && edge == EG_BOUNDARY_EDGE_VERTEX ) {
211 if ( update_node_types ) node_type->SetValue( id_node1, EG_BOUNDARY_EDGE_VERTEX );//EG_BOUNDARY_EDGE_VERTEX has priority over EG_FEATURE_EDGE_VERTEX
212 }
213 }
214
215 if ( edge && node_type->GetValue( id_node2 ) == EG_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 ) == EG_BOUNDARY_EDGE_VERTEX ) ||
221 ( edge && node_type->GetValue( id_node2 ) == EG_FEATURE_EDGE_VERTEX ) ||
222 ( !edge && node_type->GetValue( id_node2 ) == EG_SIMPLE_VERTEX ) ) {
223 m_PotentialSnapPoints[id_node2].push_back( id_node1 );
224 if ( node_type->GetValue( id_node2 ) && edge == EG_BOUNDARY_EDGE_VERTEX ) {
225 if ( update_node_types ) node_type->SetValue( id_node2, EG_BOUNDARY_EDGE_VERTEX );//EG_BOUNDARY_EDGE_VERTEX has priority over EG_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 ) == EG_FEATURE_EDGE_VERTEX || node_type->GetValue( id_node ) == EG_BOUNDARY_EDGE_VERTEX ) { //see how many edges; if two, what the angle is
240
241 if ( !this->m_BoundarySmoothing && node_type->GetValue( id_node ) == EG_BOUNDARY_EDGE_VERTEX ) {
242 if ( update_node_types ) node_type->SetValue( id_node, EG_FIXED_VERTEX );
243 } else if ( m_PotentialSnapPoints[id_node].size() != 2 ) {
244 if ( update_node_types ) node_type->SetValue( id_node, EG_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, EG_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 EG_FIXED_VERTEX;
278 }
279 }
280
281 // experimental ...
282 EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
283 char type = node_type->GetValue(id_node);
284 if (type != EG_SIMPLE_VERTEX) {
285 return type;
286 }
287
288 // initialize default value
289 //char type = EG_SIMPLE_VERTEX;
290
291 // loop through edges around id_node
292
293 QVector <vtkIdType> edges;
294 double cos_edge_angle = cos(this->m_EdgeAngle);
295
296 // safety switch to make sure edge angle fixing is disabled
297 if (cos_edge_angle < -0.99) {
298 cos_edge_angle = -2;
299 }
300
301 foreach (int i_node2, n2n[_nodes[id_node]]) {
302
303 vtkIdType id_node2 = nodes[i_node2];
304
305 // determine edge type
306 char edge = getEdgeType(id_node2, id_node, fix_unselected);
307
308 // determine node type pre-processing
309 // count number of complex edges if the node is complex,
310 // otherwise, just count the number of edges
311 //
312 if (edge != EG_SIMPLE_VERTEX && type == EG_SIMPLE_VERTEX) {
313
314 edges.clear();
315 edges.push_back(id_node2);
316 type = edge;
317
318 } else if ((edge != EG_SIMPLE_VERTEX && type == EG_BOUNDARY_EDGE_VERTEX) ||
319 (edge != EG_SIMPLE_VERTEX && type == EG_FEATURE_EDGE_VERTEX) ||
320 (edge == EG_SIMPLE_VERTEX && type == EG_SIMPLE_VERTEX))
321 {
322 edges.push_back(id_node2);
323 if (type != EG_SIMPLE_VERTEX && edge == EG_BOUNDARY_EDGE_VERTEX) {
324 // EG_BOUNDARY_EDGE_VERTEX has priority over EG_FEATURE_EDGE_VERTEX
325 type = EG_BOUNDARY_EDGE_VERTEX;
326 }
327 }
328 }
329
330 // determine node type post-processing
331 if (type == EG_FEATURE_EDGE_VERTEX || type == EG_BOUNDARY_EDGE_VERTEX) {
332
333 // check how how many edges; if two, what is the angle?
334
335 if (!this->m_BoundarySmoothing && type == EG_BOUNDARY_EDGE_VERTEX) {
336 type = EG_FIXED_VERTEX;
337 } else if (edges.size() != 2) {
338 //type = EG_FIXED_VERTEX;
339 } else {
340
341 // check angle between edges
342 vec3_t x1, x2, x3, l1, l2;
343 m_Grid->GetPoint(edges[0], x1.data());
344 m_Grid->GetPoint(id_node, x2.data());
345 m_Grid->GetPoint(edges[1], x3.data());
346 l1 = x2 - x1;
347 l2 = x3 - x2;
348 if (l1.abs() >= 0.0 && l2.abs() >= 0.0 && l1*l2 < cos_edge_angle) {
349 //type = EG_FIXED_VERTEX;
350 }
351 } // if along edge
352 } // if edge vertex
353
354 return type;
355 }
356
357 int SurfaceOperation::getEdgeCells(vtkIdType id_node1, vtkIdType id_node2, QVector <vtkIdType> &EdgeCells)
358 {
359 g2l_t _nodes = getPartLocalNodes();
360 l2g_t cells = getPartCells();
361 l2l_t n2c = getPartN2C();
362
363 QSet<vtkIdType> S1;
364 foreach (int i, n2c[_nodes[id_node1]]) {
365 S1.insert(cells[i]);
366 }
367
368 QSet<vtkIdType> S2;
369 foreach( int i, n2c[_nodes[id_node2]] ) {
370 S2.insert(cells[i]);
371 }
372
373 S2.intersect(S1);
374 EdgeCells = Set2Vector(S2, false);
375 return EdgeCells.size();
376 }
377
378 int SurfaceOperation::getEdgeCells( vtkIdType id_node1, vtkIdType id_node2, QSet <vtkIdType> &EdgeCells )
379 {
380 g2l_t _nodes = getPartLocalNodes();
381 l2g_t cells = getPartCells();
382 l2l_t n2c = getPartN2C();
383
384 QSet<vtkIdType> S1;
385 foreach( int i, n2c[_nodes[id_node1]] ) {
386 S1.insert( cells[i] );
387 }
388
389 QSet<vtkIdType> S2;
390 foreach( int i, n2c[_nodes[id_node2]] ) {
391 S2.insert( cells[i] );
392 }
393
394 EdgeCells = S2.intersect( S1 );
395 return EdgeCells.size();
396 }
397
398 char SurfaceOperation::getEdgeType(vtkIdType a_node1, vtkIdType a_node2, bool fix_unselected)
399 {
400 double cos_feature_angle = cos(this->m_FeatureAngle);
401 bool feature_edges_disabled = m_FeatureAngle >= M_PI;
402
403 // compute number of cells around edge [a_node,p2] and put them into neighbour_cells
404 QVector <vtkIdType> neighbour_cells;
405 int numNei = getEdgeCells(a_node1, a_node2, neighbour_cells) - 1;
406
407 // set default value
408 char edge = EG_SIMPLE_VERTEX;
409
410 if (numNei == 0) {
411 edge = EG_BOUNDARY_EDGE_VERTEX;
412 } else if (numNei >= 2) {
413 edge = EG_BOUNDARY_EDGE_VERTEX;
414 } else if (numNei == 1) {
415
416 // check angle between cell1 and cell2 against FeatureAngle
417 if (cosAngle(m_Grid, neighbour_cells[0], neighbour_cells[1] ) <= cos_feature_angle && !feature_edges_disabled) {
418 edge = EG_FEATURE_EDGE_VERTEX;
419 }
420
421 // check the boundary codes
422 EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
423 int cell_code_0 = cell_code->GetValue( neighbour_cells[0] );
424 int cell_code_1 = cell_code->GetValue( neighbour_cells[1] );
425 if (cell_code_0 != cell_code_1) {
426 edge = EG_BOUNDARY_EDGE_VERTEX;
427 }
428
429 if (fix_unselected) {
430 if (!m_BoundaryCodes.contains(cell_code_0) || !m_BoundaryCodes.contains(cell_code_1)) {
431 // does not make sense, but should make the points of the edge fixed
432 edge = EG_FIXED_VERTEX;
433 }
434 }
435 }
436
437 return edge;
438 }
439
440 VertexMeshDensity SurfaceOperation::getVMD( vtkIdType id_node )
441 {
442 g2l_t _nodes = getPartLocalNodes();
443 l2g_t cells = getPartCells();
444 l2l_t n2c = getPartN2C();
445
446 EG_VTKDCN( vtkCharArray, node_type, m_Grid, "node_type" );
447 EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
448
449 VertexMeshDensity VMD;
450 VMD.type = node_type->GetValue( id_node );
451 VMD.density = 0;
452 VMD.CurrentNode = id_node;
453
454 foreach( int i_cell, n2c[_nodes[id_node]] ) {
455 vtkIdType id_cell = cells[i_cell];
456 VMD.BCmap[cell_code->GetValue( id_cell )] = 2;
457 }
458 return( VMD );
459 }
460
461 //////////////////////////////////////////////
462 double SurfaceOperation::currentVertexAvgDist( vtkIdType id_node )
463 {
464 l2g_t nodes = getPartNodes();
465 g2l_t _nodes = getPartLocalNodes();
466 l2l_t n2n = getPartN2N();
467
468 double total_dist = 0;
469 double avg_dist = 0;
470 int N = n2n[_nodes[id_node]].size();
471 vec3_t C;
472 m_Grid->GetPoint( id_node, C.data() );
473 foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
474 vtkIdType id_node_neighbour = nodes[i_node_neighbour];
475 vec3_t M;
476 m_Grid->GetPoint( id_node_neighbour, M.data() );
477 total_dist += ( M - C ).abs();
478 }
479 avg_dist = total_dist / ( double )N;
480 return( avg_dist );
481 }
482
483 double SurfaceOperation::CurrentMeshDensity( vtkIdType id_node )
484 {
485 return 1.0 / currentVertexAvgDist( id_node );
486 }
487
488 ///\todo change meshdensity fields to edgelength fields since this is what is mostly used?
489
490 /// desired edge length for id_node
491 double SurfaceOperation::desiredEdgeLength( vtkIdType id_node )
492 {
493 EG_VTKDCN( vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired" );
494 return( 1.0 / characteristic_length_desired->GetValue( id_node ) );
495 }
496
497 /// mean desired edge length for id_cell
498 double SurfaceOperation::meanDesiredEdgeLength( vtkIdType id_cell )
499 {
500 vtkIdType num_pts, *pts;
501 m_Grid->GetCellPoints( id_cell, num_pts, pts );
502 int total = 0;
503 for ( int i = 0; i < num_pts; i++ ) {
504 total += desiredEdgeLength( pts[i] );
505 }
506 return total / ( double )num_pts;
507 }
508
509 QVector <vtkIdType> SurfaceOperation::getPotentialSnapPoints( vtkIdType id_node )
510 {
511 if ((id_node < 0) || (id_node >= m_PotentialSnapPoints.size())) {
512 // UpdatePotentialSnapPoints should probably be called before using this function.
513 EG_BUG;
514 }
515 return m_PotentialSnapPoints[id_node];
516 }
517
518 bool SurfaceOperation::isCell(vtkIdType id_node1, vtkIdType id_node2, vtkIdType id_node3)
519 {
520 QVector <vtkIdType> EdgeCells_12;
521 QVector <vtkIdType> EdgeCells_13;
522 QVector <vtkIdType> inter;
523
524 getEdgeCells( id_node1, id_node2, EdgeCells_12 );
525 getEdgeCells( id_node1, id_node3, EdgeCells_13 );
526 qcontIntersection( EdgeCells_12, EdgeCells_13, inter );
527 if(inter.size()>1) {
528 qWarning()<<"(id_node1, id_node2, id_node3)="<<"("<<id_node1<<", "<<id_node2<<", "<<id_node3<<")";
529 qWarning()<<"EdgeCells_12="<<EdgeCells_12;
530 qWarning()<<"EdgeCells_13="<<EdgeCells_13;
531 qWarning()<<"inter="<<inter;
532 writeGrid(m_Grid, "abort");
533 EG_BUG;// multiple cells in the same place
534 }
535 if(DebugLevel>100) {
536 qDebug()<<"(id_node1, id_node2, id_node3)="<<"("<<id_node1<<", "<<id_node2<<", "<<id_node3<<")";
537 qDebug()<<"EdgeCells_12="<<EdgeCells_12;
538 qDebug()<<"EdgeCells_13="<<EdgeCells_13;
539 qDebug()<<"inter="<<inter;
540 }
541 return(inter.size()>0);
542 }
543
544 void SurfaceOperation::computeNormals()
545 {
546 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
547 m_NodeNormal.fill(vec3_t(0,0,0), m_Grid->GetNumberOfPoints());
548 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
549 if (m_Part.localNode(id_node) != -1) {
550 QSet<int> bcs;
551 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
552 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
553 if (isSurface(id_cell, m_Grid)) {
554 int bc = cell_code->GetValue(id_cell);
555 if (m_BoundaryCodes.contains(bc)) {
556 bcs.insert(bc);
557 }
558 }
559 }
560 int num_bcs = bcs.size();
561 QVector<vec3_t> normal(num_bcs, vec3_t(0,0,0));
562 QMap<int,int> bcmap;
563 int i_bc = 0;
564 foreach (int bc, bcs) {
565 bcmap[bc] = i_bc;
566 ++i_bc;
567 }
568 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
569 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
570 if (isSurface(id_cell, m_Grid)) {
571 int bc = cell_code->GetValue(id_cell);
572 if (m_BoundaryCodes.contains(bc)) {
573 vtkIdType N_pts, *pts;
574 m_Grid->GetCellPoints(id_cell, N_pts, pts);
575 vec3_t a, b, c;
576 for (int j = 0; j < N_pts; ++j) {
577 if (pts[j] == id_node) {
578 m_Grid->GetPoint(pts[j], a.data());
579 if (j > 0) {
580 m_Grid->GetPoint(pts[j-1], b.data());
581 } else {
582 m_Grid->GetPoint(pts[N_pts-1], b.data());
583 }
584 if (j < N_pts - 1) {
585 m_Grid->GetPoint(pts[j+1], c.data());
586 } else {
587 m_Grid->GetPoint(pts[0], c.data());
588 }
589 }
590 }
591 vec3_t u = b - a;
592 vec3_t v = c - a;
593 double alpha = GeometryTools::angle(u, v);
594 vec3_t n = u.cross(v);
595 n.normalise();
596 if (checkVector(n)) {
597 normal[bcmap[bc]] += alpha*n;
598 }
599 }
600 }
601 }
602 for (int i = 0; i < num_bcs; ++i) {
603 normal[i].normalise();
604 }
605 if (num_bcs > 0) {
606 if (num_bcs > 1) {
607 if (num_bcs == 3) {
608 for (int i = 0; i < num_bcs; ++i) {
609 for (int j = i + 1; j < num_bcs; ++j) {
610 vec3_t n = normal[i] + normal[j];
611 n.normalise();
612 m_NodeNormal[id_node] += n;
613 }
614 }
615 } else {
616 for (int i = 0; i < num_bcs; ++i) {
617 m_NodeNormal[id_node] += normal[i];
618 }
619 }
620 } else {
621 m_NodeNormal[id_node] = normal[0];
622 }
623 m_NodeNormal[id_node].normalise();
624 }
625 }
626 }
627 }
628
629 double SurfaceOperation::normalIrregularity(vtkIdType id_node)
630 {
631 double nirr = 0;
632 QVector<vec3_t> nc(m_Part.n2cGSize(id_node));
633 for (int i = 0; i < nc.size(); ++i) {
634 nc[i] = GeometryTools::cellNormal(m_Grid, m_Part.n2cGG(id_node, i));
635 nc[i].normalise();
636 }
637 for (int i = 0; i < nc.size(); ++i) {
638 for (int j = i + 1; j < nc.size(); ++j) {
639 nirr += 1.0 - nc[i]*nc[j];
640 }
641 }
642 return nirr;
643 }
Clone of the official engrid repository