?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
feature magic defaults to 1 now
[engrid.git] / src / libengrid / surfaceoperation.cpp
blob1b8fbc6932af768f3511b2fbf6d02315d9099aaf
1 //
2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 // + +
4 // + This file is part of enGrid. +
5 // + +
6 // + Copyright 2008-2013 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 #include "meshqualityfaceorientation.h"
34
35 using namespace GeometryTools;
36
37 SurfaceOperation::SurfaceOperation() : Operation()
38 {
39 //default values for determining node types and for smoothing operations
40 getSet("surface meshing", "edge angle to determine fixed vertices", 180, m_EdgeAngle);
41 getSet("surface meshing", "feature angle", 180, m_FeatureAngle);
42 getSet("surface meshing", "boundary codes define features", true, m_BCodeFeatureDefinition);
43 getSet("surface meshing", "number of steps to protect node types", 2, m_TypeProtectionCount);
44 getSet("surface meshing", "threshold for face orientation quality", 0.1, m_FaceOrientationThreshold);
45 m_FeatureAngle = GeometryTools::deg2rad(m_FeatureAngle);
46 m_EdgeAngle = GeometryTools::deg2rad(m_EdgeAngle);
47 setEdgeAngle(m_EdgeAngle);
48 m_StretchingFactor = 0;
49 m_UniformSnapPoints = false;
50 m_StrictFeatureSnap = true;
51 }
52
53 void SurfaceOperation::operate()
54 {
55
56 }
57
58 ostream& operator<<(ostream &out, stencil_t S)
59 {
60 out << "S.id_cell = " << S.id_cell << " ";
61 out << "S.id_node = " << S.id_node << " ";
62 out << "S.sameBC = " << S.sameBC << " ";
63 out << "S.type = " << S.type_cell << " ";
64 out << "S.p1 = " << S.p1 << " ";
65 out << "S.p2 = " << S.p2 << " ";
66 return(out);
67 }
68
69 stencil_t SurfaceOperation::getStencil(vtkIdType id_cell1, int j1)
70 {
71 stencil_t S;
72 {
73 vtkIdType N_pts, *pts;
74 m_Grid->GetCellPoints(id_cell1, N_pts, pts);
75 S.p1 = pts[j1];
76 S.p2 = pts[0];
77 if (j1 < N_pts - 1) {
78 S.p2 = pts[j1 + 1];
79 }
80 }
81 QSet<vtkIdType> cells_p1;
82 for (int i = 0; i < m_Part.n2cGSize(S.p1); ++i) {
83 vtkIdType id_cell = m_Part.n2cGG(S.p1, i);
84 if (id_cell != id_cell1) {
85 cells_p1.insert(id_cell);
86 }
87 }
88 QSet<vtkIdType> cells_p2;
89 for (int i = 0; i < m_Part.n2cGSize(S.p2); ++i) {
90 vtkIdType id_cell = m_Part.n2cGG(S.p2, i);
91 if (id_cell != id_cell1) {
92 cells_p2.insert(id_cell);
93 }
94 }
95 QSet<vtkIdType> cells = cells_p1.intersect(cells_p2);
96 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
97 S.sameBC = true;
98 S.id_cell.resize(1);
99 S.id_cell[0] = id_cell1;
100 foreach (vtkIdType id_cell, cells) {
101 if (isSurface(id_cell, m_Grid)) {
102 S.id_cell.push_back(id_cell);
103 if (cell_code->GetValue(id_cell) != cell_code->GetValue(id_cell1)) {
104 S.sameBC = false;
105 }
106 }
107 }
108 S.id_node.resize(S.id_cell.size());
109 S.type_cell.resize(S.id_cell.size());
110 for (int i = 0; i < S.id_cell.size(); ++i) {
111 vtkIdType N_pts, *pts;
112 m_Grid->GetCellPoints(S.id_cell[i], N_pts, pts);
113 S.type_cell[i] = m_Grid->GetCellType(S.id_cell[i]);
114 for (int j = 0; j < N_pts; ++j) {
115 if (pts[j] != S.p1 && pts[j] != S.p2) {
116 S.id_node[i] = pts[j];
117 break;
118 }
119 }
120 }
121 return S;
122 }
123
124 void SurfaceOperation::setBCodesFeatureDefinition(bool flag)
125 {
126 m_BCodeFeatureDefinition = flag;
127 if (m_BCodeFeatureDefinition) {
128 m_FeatureAngle = deg2rad(180);
129 }
130 }
131
132 int SurfaceOperation::UpdateCurrentMeshDensity()
133 {
134 if ( DebugLevel > 0 ) {
135 cout << "===UpdateMeshDensity START===" << endl;
136 }
137 QVector<vtkIdType> cells;
138 getAllSurfaceCells( cells, m_Grid );
139 EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
140 EG_VTKDCN( vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired" );
141 setGrid( m_Grid );
142 setCells( cells );
143 if ( DebugLevel > 5 ) {
144 cout << "cells.size()=" << cells.size() << endl;
145 }
146 EG_VTKDCN( vtkDoubleArray, node_meshdensity_current, m_Grid, "node_meshdensity_current" );
147 l2g_t nodes = getPartNodes();
148 foreach( vtkIdType node, nodes ) {
149 node_meshdensity_current->SetValue( node, CurrentMeshDensity( node ) );
150 }
151 if ( DebugLevel > 0 ) {
152 cout << "===UpdateMeshDensity END===" << endl;
153 }
154 return( 0 ); ///\todo what for???
155 }
156
157 void SurfaceOperation::readVMD()
158 {
159 QString buffer = GuiMainWindow::pointer()->getXmlSection("engrid/surface/table").replace("\n", " ");
160 int row_count = 0;
161 int column_count = 0;
162 m_VMDvector.clear();
163
164 if(!buffer.isEmpty()) {
165 QTextStream in(&buffer, QIODevice::ReadOnly);
166 in >> row_count >> column_count;
167 QVector<int> tmp_bcs;
168 GuiMainWindow::pointer()->getAllBoundaryCodes(tmp_bcs);
169 if (column_count == tmp_bcs.size() + 3) {
170 m_VMDvector.fill(VertexMeshDensity(), row_count);
171 for (int i = 0; i < row_count; ++i) {
172 int row, column;
173 QString formula;
174 foreach (int bc, tmp_bcs) {
175 in >> row >> column >> formula;
176 m_VMDvector[row].BCmap[bc] = formula.toInt();
177 }
178 in >> row >> column >> formula;
179 m_VMDvector[row].type = Str2VertexType(formula);
180 in >> row >> column >> formula;
181 if (formula == "{{{empty}}}") {
182 formula = "";
183 }
184 m_VMDvector[i].setNodes(formula);
185 in >> row >> column >> formula;
186 m_VMDvector[i].density = formula.toDouble();
187 }
188 } else {
189 EG_ERR_RETURN(QObject::tr("Mismatch of number of boundary codes!"));
190 }
191 }
192 }
193
194 void SurfaceOperation::updateNodeInfo()
195 {
196 setAllCells();
197 readVMD();
198 l2g_t nodes = getPartNodes();
199 computeNormals();
200 EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
201 EG_VTKDCN(vtkIntArray, node_type_counter, m_Grid, "node_type_counter");
202 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
203 foreach (vtkIdType id_node, nodes) {
204
205 char old_type = node_type->GetValue(id_node);
206 char new_type = getNodeType(id_node, true);
207
208 if (old_type == EG_FIXED_VERTEX && new_type != EG_FIXED_VERTEX) {
209 //EG_BUG;
210 }
211
212 if ((old_type != EG_FEATURE_CORNER_VERTEX && old_type != EG_FEATURE_EDGE_VERTEX) || m_BCodeFeatureDefinition) {
213 node_type->SetValue(id_node, new_type);
214 } else {
215
216 bool done = false;
217
218 // EG_FEATURE_CORNER_VERTEX -> any other type
219 if (old_type == EG_FEATURE_CORNER_VERTEX && new_type != EG_FEATURE_CORNER_VERTEX) {
220 if (node_type_counter->GetValue(id_node) > m_TypeProtectionCount) {
221 node_type->SetValue(id_node, new_type);
222 } else {
223 node_type_counter->SetValue(id_node, node_type_counter->GetValue(id_node) + 1);
224 }
225 done = true;
226 }
227
228 // EG_FEATURE_EDGE_VERTEX -> EG_SIMPLE_VERTEX
229 if (old_type == EG_FEATURE_EDGE_VERTEX && new_type == EG_SIMPLE_VERTEX) {
230 if (node_type_counter->GetValue(id_node) > m_TypeProtectionCount) {
231 node_type->SetValue(id_node, new_type);
232 } else {
233 node_type_counter->SetValue(id_node, node_type_counter->GetValue(id_node) + 1);
234 }
235 done = true;
236 }
237
238 if (!done) {
239 node_type->SetValue(id_node, new_type);
240 node_type_counter->SetValue(id_node, 0);
241 }
242
243 }
244
245 //density index from table
246 EG_VTKDCN(vtkIntArray, node_specified_density, m_Grid, "node_specified_density");
247
248 VertexMeshDensity nodeVMD = getVMD(id_node);
249 int idx = nodeVMD.findSmallestVMD(m_VMDvector);
250 node_specified_density->SetValue(id_node, idx);
251 }
252
253 // mesh quality
254 if (!m_BCodeFeatureDefinition) {
255 MeshQualityFaceOrientation mesh_quality;
256 mesh_quality();
257 EG_VTKDCN(vtkDoubleArray, node_mesh_quality, m_Grid, "node_mesh_quality");
258 EG_FORALL_NODES(id_node, m_Grid) {
259 if (node_mesh_quality->GetValue(id_node) < m_FaceOrientationThreshold) {
260 node_type->SetValue(id_node, EG_SIMPLE_VERTEX);
261 }
262 }
263 }
264
265 /*
266 // look for feature edge triple stars
267 QList<vtkIdType> new_corners;
268 EG_FORALL_NODES(id_node1, m_Grid) {
269 if (node_type->GetValue(id_node1) == EG_FEATURE_EDGE_VERTEX) {
270 int N = 0;
271 for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
272 vtkIdType id_node2 = m_Part.n2cGG(id_node1, i);
273 if (node_type->GetValue(id_node2) == EG_FEATURE_EDGE_VERTEX) {
274 ++N;
275 }
276 }
277 if (N == 3) {
278 new_corners.append(id_node1);
279 }
280 }
281 }
282 foreach (vtkIdType id_node, new_corners) {
283 node_type->SetValue(id_node, EG_FEATURE_CORNER_VERTEX);
284 }
285 */
286
287 updatePotentialSnapPoints();
288 }
289
290 bool SurfaceOperation::checkSnapPointPairForBcMatch(vtkIdType id_node1, vtkIdType id_node2)
291 {
292 QSet<int> bcs1, bcs2;
293 for (int i = 0; i < m_Part.n2bcGSize(id_node1); ++i) {
294 bcs1.insert(m_Part.n2bcG(id_node1, i));
295 }
296 for (int i = 0; i < m_Part.n2bcGSize(id_node2); ++i) {
297 bcs2.insert(m_Part.n2bcG(id_node2, i));
298 }
299 if (bcs1 == bcs2) {
300 return true;
301 }
302 return false;
303 }
304
305 void SurfaceOperation::updatePotentialSnapPoints()
306 {
307 setAllSurfaceCells();
308 l2g_t nodes = getPartNodes();
309
310 m_PotentialSnapPoints.resize(m_Grid->GetNumberOfPoints());
311
312 if (m_UniformSnapPoints) {
313 m_PotentialSnapPoints.resize(m_Grid->GetNumberOfPoints());
314 foreach( vtkIdType id_node, nodes ) {
315 for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
316 m_PotentialSnapPoints[id_node].append(m_Part.n2nGG(id_node, i));
317 }
318 }
319 return;
320 }
321
322 EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
323
324 foreach( vtkIdType id_node1, nodes ) {
325 m_PotentialSnapPoints[id_node1].clear();
326 char type1 = node_type->GetValue(id_node1);
327 if (type1 != EG_FIXED_VERTEX) { // fixed vertices do not have any snap-points
328 QSet<vtkIdType> exclude_nodes;
329 if (type1 == EG_FEATURE_EDGE_VERTEX || type1 == EG_BOUNDARY_EDGE_VERTEX) {
330 for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
331 vtkIdType id_node2 = m_Part.n2nGG(id_node1, i);
332 char type2 = node_type->GetValue(id_node2);
333 if (type2 == EG_FEATURE_CORNER_VERTEX || type2 == EG_FIXED_VERTEX) {
334 for (int j = 0; j < m_Part.n2nGSize(id_node2); ++j) {
335 vtkIdType id_node3 = m_Part.n2nGG(id_node2, j);
336 exclude_nodes.insert(id_node3);
337 }
338 }
339 }
340 }
341 for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
342 vtkIdType id_node2 = m_Part.n2nGG(id_node1, i);
343 char type2 = node_type->GetValue(id_node2);
344 if (m_StrictFeatureSnap) {
345 if ( (type1 == EG_SIMPLE_VERTEX)
346 || (type1 == EG_FEATURE_EDGE_VERTEX && (type2 == EG_FEATURE_EDGE_VERTEX || type2 == EG_FEATURE_CORNER_VERTEX))
347 || (type1 == EG_FEATURE_CORNER_VERTEX && type2 == EG_FEATURE_CORNER_VERTEX)
348 || (type1 == EG_BOUNDARY_EDGE_VERTEX && (type2 == EG_BOUNDARY_EDGE_VERTEX || type2 == EG_FIXED_VERTEX)))
349 {
350 if (!exclude_nodes.contains(id_node2)) {
351 if (checkSnapPointPairForBcMatch(id_node1, id_node2)) {
352 m_PotentialSnapPoints[id_node1].append(id_node2);
353 }
354 }
355 }
356 } else {
357 if ( (type1 == EG_SIMPLE_VERTEX)
358 || (type1 == EG_FEATURE_EDGE_VERTEX)
359 || (type1 == EG_BOUNDARY_EDGE_VERTEX && (type2 == EG_BOUNDARY_EDGE_VERTEX || type2 == EG_FIXED_VERTEX)))
360 {
361 if (checkSnapPointPairForBcMatch(id_node1, id_node2)) {
362 m_PotentialSnapPoints[id_node1].append(id_node2);
363 }
364 }
365 }
366 }
367
368 // make sure feature edge vertices have at least two snap points ...
369 /*
370 if (type1 == EG_FEATURE_EDGE_VERTEX) {
371 if (m_PotentialSnapPoints[id_node1].size() < 2) {
372 m_PotentialSnapPoints[id_node1].clear();
373 }
374 }
375 */
376 }
377 }
378 }
379
380 char SurfaceOperation::getNodeType(vtkIdType id_node, bool fix_unselected)
381 {
382 if (m_Part.n2nGSize(id_node) == 0) {
383 return EG_FIXED_VERTEX;
384 }
385
386 //char type = EG_SIMPLE_VERTEX;
387 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
388
389 QSet<int> bcs;
390 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
391 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
392
393 // fix all vertices that are part of a volume element or a quad
394 if (m_Grid->GetCellType(id_cell) != VTK_TRIANGLE) {
395 return EG_FIXED_VERTEX;
396 }
397
398 int bc = cell_code->GetValue(id_cell);
399
400 // fix nodes which belong to faces with unselected boundary codes
401 if (!m_BoundaryCodes.contains(bc) && fix_unselected) {
402 return EG_FIXED_VERTEX;
403 }
404
405 bcs.insert(bc);
406 }
407
408 if (bcs.size() >= 3 || bcs.size() == 0) {
409 return EG_FIXED_VERTEX;
410 }
411 if (bcs.size() == 2) {
412 return EG_BOUNDARY_EDGE_VERTEX;
413 }
414
415 CadInterface* cad_interface = GuiMainWindow::pointer()->getCadInterface(*bcs.begin(), true);
416 if (cad_interface && !m_BCodeFeatureDefinition) {
417
418 vec3_t x, x0;
419 m_Grid->GetPoint(id_node, x0.data());
420
421 // compute average edge length
422 double L = 0;
423 for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
424 vec3_t xn;
425 m_Grid->GetPoint(m_Part.n2nGG(id_node, i), xn.data());
426 L += (x0 - xn).abs();
427 }
428 L /= m_Part.n2nGSize(id_node);
429 L *= 0.1;
430 bool convex = isConvexNode(id_node);
431
432 x0 = cad_interface->projectNode(id_node, x0, m_NodeNormal[id_node]);
433 if (convex) {
434 x = x0 - L*m_NodeNormal[id_node];
435 } else {
436 x = x0 + L*m_NodeNormal[id_node];
437 }
438 vec3_t n = cad_interface->getLastNormal();
439 if (GeometryTools::angle(n, m_NodeNormal[id_node]) > 0.5*m_FeatureAngle && !cad_interface->failed()) {
440 double d = L/tan(0.5*m_FeatureAngle);
441 static const int num_steps = 36;
442 double D_alpha = 2*M_PI/num_steps;
443 vec3_t v;
444
445 v = GeometryTools::orthogonalVector(m_NodeNormal[id_node]);
446 int num_miss = 0;
447 for (int i = 0; i < num_steps; ++i) {
448 v = GeometryTools::rotate(v, m_NodeNormal[id_node], D_alpha);
449 vec3_t xp = cad_interface->projectNode(id_node, x, v, true);
450 if (cad_interface->failed()) {
451 ++num_miss;
452 } else {
453 double l = (x - xp).abs();
454 if (l >= d) {
455 ++num_miss;
456 }
457 }
458 }
459 if (num_miss == 0) {
460 return EG_FEATURE_CORNER_VERTEX;
461 }
462
463 v = GeometryTools::orthogonalVector(n);
464 int num_hit = 0;
465 for (int i = 0; i < num_steps; ++i) {
466 v = GeometryTools::rotate(v, n, D_alpha);
467 vec3_t xp = cad_interface->projectNode(id_node, x, v, true);
468 if (!cad_interface->failed()) {
469 double l = (x - xp).abs();
470 if (l < d) {
471 ++num_hit;
472 }
473 }
474 }
475 if (num_hit > 0) {
476 return EG_FEATURE_EDGE_VERTEX;
477 }
478 }
479 }
480
481 return EG_SIMPLE_VERTEX;
482 }
483
484 int SurfaceOperation::getEdgeCells(vtkIdType id_node1, vtkIdType id_node2, QVector <vtkIdType> &EdgeCells)
485 {
486 g2l_t _nodes = getPartLocalNodes();
487 l2g_t cells = getPartCells();
488 l2l_t n2c = getPartN2C();
489
490 QSet<vtkIdType> S1;
491 foreach (int i, n2c[_nodes[id_node1]]) {
492 S1.insert(cells[i]);
493 }
494
495 QSet<vtkIdType> S2;
496 foreach( int i, n2c[_nodes[id_node2]] ) {
497 S2.insert(cells[i]);
498 }
499
500 S2.intersect(S1);
501 EdgeCells = Set2Vector(S2, false);
502 return EdgeCells.size();
503 }
504
505 int SurfaceOperation::getEdgeCells( vtkIdType id_node1, vtkIdType id_node2, QSet <vtkIdType> &EdgeCells )
506 {
507 g2l_t _nodes = getPartLocalNodes();
508 l2g_t cells = getPartCells();
509 l2l_t n2c = getPartN2C();
510
511 QSet<vtkIdType> S1;
512 foreach( int i, n2c[_nodes[id_node1]] ) {
513 S1.insert( cells[i] );
514 }
515
516 QSet<vtkIdType> S2;
517 foreach( int i, n2c[_nodes[id_node2]] ) {
518 S2.insert( cells[i] );
519 }
520
521 EdgeCells = S2.intersect( S1 );
522 return EdgeCells.size();
523 }
524
525 char SurfaceOperation::getEdgeType(vtkIdType a_node1, vtkIdType a_node2, bool fix_unselected)
526 {
527 double cos_feature_angle = cos(this->m_FeatureAngle);
528 bool feature_edges_disabled = m_FeatureAngle >= M_PI;
529
530 // compute number of cells around edge [a_node,p2] and put them into neighbour_cells
531 QVector <vtkIdType> neighbour_cells;
532 int numNei = getEdgeCells(a_node1, a_node2, neighbour_cells) - 1;
533
534 // set default value
535 char edge = EG_SIMPLE_VERTEX;
536
537 if (numNei == 0) {
538 edge = EG_BOUNDARY_EDGE_VERTEX;
539 } else if (numNei >= 2) {
540 edge = EG_BOUNDARY_EDGE_VERTEX;
541 } else if (numNei == 1) {
542
543 // check angle between cell1 and cell2 against FeatureAngle
544 double cosa = cosAngle(m_Grid, neighbour_cells[0], neighbour_cells[1]);
545 if (cosa <= cos_feature_angle && !feature_edges_disabled) {
546 edge = EG_FEATURE_EDGE_VERTEX;
547 }
548
549 // check the boundary codes
550 EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
551 int cell_code_0 = cell_code->GetValue( neighbour_cells[0] );
552 int cell_code_1 = cell_code->GetValue( neighbour_cells[1] );
553 if (cell_code_0 != cell_code_1) {
554 edge = EG_BOUNDARY_EDGE_VERTEX;
555 }
556
557 if (fix_unselected) {
558 if (!m_BoundaryCodes.contains(cell_code_0) || !m_BoundaryCodes.contains(cell_code_1)) {
559 // does not make sense, but should make the points of the edge fixed
560 edge = EG_FIXED_VERTEX;
561 }
562 }
563 }
564
565 return edge;
566 }
567
568 VertexMeshDensity SurfaceOperation::getVMD( vtkIdType id_node )
569 {
570 g2l_t _nodes = getPartLocalNodes();
571 l2g_t cells = getPartCells();
572 l2l_t n2c = getPartN2C();
573
574 EG_VTKDCN( vtkCharArray, node_type, m_Grid, "node_type" );
575 EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
576
577 VertexMeshDensity VMD;
578 VMD.type = node_type->GetValue( id_node );
579 VMD.density = 0;
580 VMD.CurrentNode = id_node;
581
582 foreach( int i_cell, n2c[_nodes[id_node]] ) {
583 vtkIdType id_cell = cells[i_cell];
584 VMD.BCmap[cell_code->GetValue( id_cell )] = 2;
585 }
586 return( VMD );
587 }
588
589 //////////////////////////////////////////////
590 double SurfaceOperation::currentVertexAvgDist( vtkIdType id_node )
591 {
592 l2g_t nodes = getPartNodes();
593 g2l_t _nodes = getPartLocalNodes();
594 l2l_t n2n = getPartN2N();
595
596 double total_dist = 0;
597 double avg_dist = 0;
598 int N = n2n[_nodes[id_node]].size();
599 vec3_t C;
600 m_Grid->GetPoint( id_node, C.data() );
601 foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
602 vtkIdType id_node_neighbour = nodes[i_node_neighbour];
603 vec3_t M;
604 m_Grid->GetPoint( id_node_neighbour, M.data() );
605 total_dist += ( M - C ).abs();
606 }
607 avg_dist = total_dist / ( double )N;
608 return( avg_dist );
609 }
610
611 double SurfaceOperation::CurrentMeshDensity( vtkIdType id_node )
612 {
613 return 1.0 / currentVertexAvgDist( id_node );
614 }
615
616 ///\todo change meshdensity fields to edgelength fields since this is what is mostly used?
617
618 /// desired edge length for id_node
619 double SurfaceOperation::desiredEdgeLength( vtkIdType id_node )
620 {
621 EG_VTKDCN( vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired" );
622 return( 1.0 / characteristic_length_desired->GetValue( id_node ) );
623 }
624
625 /// mean desired edge length for id_cell
626 double SurfaceOperation::meanDesiredEdgeLength( vtkIdType id_cell )
627 {
628 vtkIdType num_pts, *pts;
629 m_Grid->GetCellPoints( id_cell, num_pts, pts );
630 int total = 0;
631 for ( int i = 0; i < num_pts; i++ ) {
632 total += desiredEdgeLength( pts[i] );
633 }
634 return total / ( double )num_pts;
635 }
636
637 QVector <vtkIdType> SurfaceOperation::getPotentialSnapPoints( vtkIdType id_node )
638 {
639 if ((id_node < 0) || (id_node >= m_PotentialSnapPoints.size())) {
640 // UpdatePotentialSnapPoints should probably be called before using this function.
641 EG_BUG;
642 }
643 return m_PotentialSnapPoints[id_node];
644 }
645
646 bool SurfaceOperation::isCell(vtkIdType id_node1, vtkIdType id_node2, vtkIdType id_node3)
647 {
648 QVector <vtkIdType> EdgeCells_12;
649 QVector <vtkIdType> EdgeCells_13;
650 QVector <vtkIdType> inter;
651
652 getEdgeCells( id_node1, id_node2, EdgeCells_12 );
653 getEdgeCells( id_node1, id_node3, EdgeCells_13 );
654 qcontIntersection( EdgeCells_12, EdgeCells_13, inter );
655 if(inter.size()>1) {
656 qWarning()<<"(id_node1, id_node2, id_node3)="<<"("<<id_node1<<", "<<id_node2<<", "<<id_node3<<")";
657 qWarning()<<"EdgeCells_12="<<EdgeCells_12;
658 qWarning()<<"EdgeCells_13="<<EdgeCells_13;
659 qWarning()<<"inter="<<inter;
660 writeGrid(m_Grid, "abort");
661 EG_BUG;// multiple cells in the same place
662 }
663 if(DebugLevel>100) {
664 qDebug()<<"(id_node1, id_node2, id_node3)="<<"("<<id_node1<<", "<<id_node2<<", "<<id_node3<<")";
665 qDebug()<<"EdgeCells_12="<<EdgeCells_12;
666 qDebug()<<"EdgeCells_13="<<EdgeCells_13;
667 qDebug()<<"inter="<<inter;
668 }
669 return(inter.size()>0);
670 }
671
672 void SurfaceOperation::computeNormals()
673 {
674 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
675 m_NodeNormal.fill(vec3_t(0,0,0), m_Grid->GetNumberOfPoints());
676 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
677 if (m_Part.localNode(id_node) != -1) {
678 QSet<int> bcs;
679 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
680 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
681 if (isSurface(id_cell, m_Grid)) {
682 int bc = cell_code->GetValue(id_cell);
683 if (m_BoundaryCodes.contains(bc)) {
684 bcs.insert(bc);
685 }
686 }
687 }
688 int num_bcs = bcs.size();
689 QVector<vec3_t> normal(num_bcs, vec3_t(0,0,0));
690 QMap<int,int> bcmap;
691 int i_bc = 0;
692 foreach (int bc, bcs) {
693 bcmap[bc] = i_bc;
694 ++i_bc;
695 }
696 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
697 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
698 if (isSurface(id_cell, m_Grid)) {
699 int bc = cell_code->GetValue(id_cell);
700 if (m_BoundaryCodes.contains(bc)) {
701 vtkIdType N_pts, *pts;
702 m_Grid->GetCellPoints(id_cell, N_pts, pts);
703 vec3_t a, b, c;
704 for (int j = 0; j < N_pts; ++j) {
705 if (pts[j] == id_node) {
706 m_Grid->GetPoint(pts[j], a.data());
707 if (j > 0) {
708 m_Grid->GetPoint(pts[j-1], b.data());
709 } else {
710 m_Grid->GetPoint(pts[N_pts-1], b.data());
711 }
712 if (j < N_pts - 1) {
713 m_Grid->GetPoint(pts[j+1], c.data());
714 } else {
715 m_Grid->GetPoint(pts[0], c.data());
716 }
717 }
718 }
719 vec3_t u = b - a;
720 vec3_t v = c - a;
721 double alpha = GeometryTools::angle(u, v);
722 vec3_t n = u.cross(v);
723 n.normalise();
724 if (checkVector(n)) {
725 normal[bcmap[bc]] -= alpha*n;
726 }
727 }
728 }
729 }
730 for (int i = 0; i < num_bcs; ++i) {
731 normal[i].normalise();
732 }
733 if (num_bcs > 0) {
734 if (num_bcs > 1) {
735 if (num_bcs == 3) {
736 for (int i = 0; i < num_bcs; ++i) {
737 for (int j = i + 1; j < num_bcs; ++j) {
738 vec3_t n = normal[i] + normal[j];
739 n.normalise();
740 m_NodeNormal[id_node] += n;
741 }
742 }
743 } else {
744 for (int i = 0; i < num_bcs; ++i) {
745 m_NodeNormal[id_node] += normal[i];
746 }
747 }
748 } else {
749 m_NodeNormal[id_node] = normal[0];
750 }
751 m_NodeNormal[id_node].normalise();
752 }
753 }
754 }
755 }
756
757 double SurfaceOperation::normalIrregularity(vtkIdType id_node)
758 {
759 double nirr = 0;
760 QVector<vec3_t> nc(m_Part.n2cGSize(id_node));
761 for (int i = 0; i < nc.size(); ++i) {
762 nc[i] = GeometryTools::cellNormal(m_Grid, m_Part.n2cGG(id_node, i));
763 nc[i].normalise();
764 }
765 for (int i = 0; i < nc.size(); ++i) {
766 for (int j = i + 1; j < nc.size(); ++j) {
767 nirr += 1.0 - nc[i]*nc[j];
768 }
769 }
770 return nirr;
771 }
772
773 bool SurfaceOperation::isConvexNode(vtkIdType id_node)
774 {
775 int N = m_Part.n2nGSize(id_node);
776 if (N == 0) {
777 return false;
778 }
779 vec3_t x1, x2(0,0,0);
780 m_Grid->GetPoint(id_node, x1.data());
781 for (int i = 0; i < N; ++i) {
782 vec3_t x;
783 m_Grid->GetPoint(m_Part.n2nGG(id_node, i), x.data());
784 x2 += x;
785 }
786 x2 *= 1.0/N;
787 if ((x1 - x2)*m_NodeNormal[id_node] > 0) {
788 return true;
789 }
790 return false;
791 }
792
793 double SurfaceOperation::getSurfaceDeviation(vtkIdType id_node)
794 {
795 vec3_t x;
796 m_Grid->GetPoint(id_node, x.data());
797 double d = 0;
798 for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
799 vec3_t xe;
800 m_Grid->GetPoint(m_Part.n2nGG(id_node, i), xe.data());
801 d += (x - xe)*m_NodeNormal[id_node];
802 }
803 if (m_Part.n2nGSize(id_node) > 0) {
804 d *= 1.0/m_Part.n2nGSize(id_node);
805 }
806 return d;
807 }
808
809 bool SurfaceOperation::isFeatureNode(vtkIdType id_node)
810 {
811 EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
812 char type = node_type->GetValue(id_node);
813 if (type == EG_FEATURE_CORNER_VERTEX || type == EG_FEATURE_EDGE_VERTEX) {
814 return true;
815 }
816 return false;
817 }
Clone of the official engrid repository