?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
updated to modern VTK
[engrid-github.git] / src / libengrid / updatedesiredmeshdensity.cpp
blobee6a7f0e509d78fdba502ded8ab69c649f8ec210
1 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 // + +
3 // + This file is part of enGrid. +
4 // + +
5 // + Copyright 2008-2014 enGits GmbH +
6 // + +
7 // + enGrid is free software: you can redistribute it and/or modify +
8 // + it under the terms of the GNU General Public License as published by +
9 // + the Free Software Foundation, either version 3 of the License, or +
10 // + (at your option) any later version. +
11 // + +
12 // + enGrid is distributed in the hope that it will be useful, +
13 // + but WITHOUT ANY WARRANTY; without even the implied warranty of +
14 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +
15 // + GNU General Public License for more details. +
16 // + +
17 // + You should have received a copy of the GNU General Public License +
18 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>. +
19 // + +
20 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
21 #include "updatedesiredmeshdensity.h"
22 #include "engrid.h"
23 #include "guimainwindow.h"
24 #include "pointfinder.h"
25 #include "cgaltricadinterface.h"
26
27 #include <vtkSignedCharArray.h>
28
29 UpdateDesiredMeshDensity::UpdateDesiredMeshDensity() : SurfaceOperation()
30 {
31 EG_TYPENAME;
32 m_MaxEdgeLength = 1e99;
33 m_NodesPerQuarterCircle = 0;
34 m_OnlySurfaceCells = true;
35
36 m_GrowthFactor = 0.0;
37 m_MinEdgeLength = 0.0;
38 m_FeatureResolution2D = 0;
39 m_FeatureResolution3D = 0;
40
41 getSet("surface meshing", "threshold angle for feature resolution", 20, m_FeatureThresholdAngle);
42 m_FeatureThresholdAngle = deg2rad(m_FeatureThresholdAngle);
43 m_Relaxation = true;
44
45 }
46
47 double UpdateDesiredMeshDensity::computeSearchDistance(vtkIdType id_face)
48 {
49 EG_GET_CELL(id_face, m_Grid);
50 QVector<vec3_t> x(num_pts + 1);
51 for (int i = 0; i < num_pts; ++i) {
52 m_Grid->GetPoint(pts[i], x[i].data());
53 }
54 x[num_pts] = x[0];
55 double L = 0;
56 for (int i = 0; i < num_pts; ++i) {
57 L = max(L, (x[i] - x[i+1]).abs());
58 }
59 return L;
60 }
61
62 void UpdateDesiredMeshDensity::computeExistingLengths()
63 {
64 QSet<int> all_bcs = GuiMainWindow::pointer()->getAllBoundaryCodes();
65 QSet<int> fixed_bcs = all_bcs - m_BoundaryCodes;
66 QVector<double> edge_length(m_Grid->GetNumberOfPoints(), 1e99);
67 QVector<int> edge_count(m_Grid->GetNumberOfPoints(), 0);
68 m_Fixed.fill(false, m_Grid->GetNumberOfPoints());
69 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
70 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
71 if (isSurface(id_cell, m_Grid)) {
72 if (fixed_bcs.contains(cell_code->GetValue(id_cell))) {
73 EG_GET_CELL(id_cell, m_Grid);
74 QVector<vec3_t> x(num_pts);
75 for (int i = 0; i < num_pts; ++i) {
76 m_Grid->GetPoint(pts[i], x[i].data());
77 m_Fixed[pts[i]] = true;
78 }
79 for (int i = 0; i < num_pts; ++i) {
80 int j = i + 1;
81 if (j >= num_pts) {
82 j = 0;
83 }
84 double L = (x[i] - x[j]).abs();
85 edge_length[pts[i]] = min(edge_length[pts[i]], L);
86 edge_length[pts[j]] = min(edge_length[pts[j]], L);
87 ++edge_count[pts[i]];
88 ++edge_count[pts[j]];
89 }
90 }
91 }
92 }
93 EG_VTKDCN(vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired");
94 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
95 if (edge_count[id_node] > 0) {
96 if (edge_length[id_node] > 1e98) {
97 EG_BUG;
98 }
99 characteristic_length_desired->SetValue(id_node, edge_length[id_node]);
100 }
101 }
102 }
103
104 void UpdateDesiredMeshDensity::computeFeature(const QList<point_t> points, QVector<double> &cl_pre, double res, int restriction_type)
105 {
106 int N = 0;
107 QVector<vec3_t> pts(points.size());
108 for (int i = 0; i < points.size(); ++i) {
109 pts[i] = points[i].x;
110 }
111 PointFinder pfind;
112 pfind.setMaxNumPoints(5000);
113 pfind.setPoints(pts);
114 for (int i = 0; i < points.size(); ++i) {
115 double h = -1;
116 foreach (int i_points, points[i].idx) {
117 h = max(h, cl_pre[i_points]);
118 }
119 if (h < 0) {
120 h = 1e99;
121 }
122 QVector<int> close_points;
123 pfind.getClosePoints(points[i].x, close_points, res*points[i].L);
124 foreach (int j, close_points) {
125 if (i != j) {
126 ++N;
127 vec3_t x1 = points[i].x;
128 vec3_t x2 = points[j].x;
129 vec3_t n1 = points[i].n;
130 vec3_t n2 = points[j].n;
131 vec3_t v = x2 - x1;
132 if (n1*n2 < 0) {
133 if (n1*v > 0) {
134 if (fabs(GeometryTools::angle(n1, (-1)*n2)) <= m_FeatureThresholdAngle) {
135 double l = v.abs()/fabs(n1*n2);
136 if (l/res < h) {
137 // check topological distance
138 int search_dist = int(ceil(res));
139 bool topo_dist_ok = true;
140 foreach (int k, points[i].idx) {
141 vtkIdType id_node1 = m_Part.globalNode(k);
142 foreach (int l, points[j].idx) {
143 vtkIdType id_node2 = m_Part.globalNode(l);
144 if (m_Part.computeTopoDistance(id_node1, id_node2, search_dist, restriction_type) < search_dist) {
145 topo_dist_ok = false;
146 break;
147 }
148 if (!topo_dist_ok) {
149 break;
150 }
151 }
152 }
153 if (topo_dist_ok) {
154 h = l/res;
155 }
156 }
157 }
158 }
159 }
160 }
161 }
162 foreach (int i_points, points[i].idx) {
163 cl_pre[i_points] = min(h, cl_pre[i_points]);
164 }
165 }
166 }
167
168 void UpdateDesiredMeshDensity::computeFeature2D(QVector<double> &cl_pre)
169 {
170 if (m_FeatureResolution2D < 1e-3) {
171 return;
172 }
173 QSet<int> bcs = getAllBoundaryCodes(m_Grid);
174 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
175 foreach (int bc, bcs) {
176 QList<point_t> points;
177 for (vtkIdType id_face = 0; id_face < m_Grid->GetNumberOfCells(); ++id_face) {
178 if (isSurface(id_face, m_Grid)) {
179 if (cell_code->GetValue(id_face) == bc) {
180 vec3_t xc = cellCentre(m_Grid, id_face);
181 EG_GET_CELL(id_face, m_Grid);
182 QVector<vec3_t> xn(num_pts + 1);
183 QVector<vtkIdType> idx(num_pts + 1);
184 for (int i_pts = 0; i_pts < num_pts; ++i_pts) {
185 m_Grid->GetPoint(pts[i_pts], xn[i_pts].data());
186 idx[i_pts] = m_Part.localNode(pts[i_pts]);
187 }
188 xn[num_pts] = xn[0];
189 idx[num_pts] = idx[0];
190 for (int i_neigh = 0; i_neigh < m_Part.c2cGSize(id_face); ++i_neigh) {
191 vtkIdType id_neigh = m_Part.c2cGG(id_face, i_neigh);
192 if (id_neigh != -1) {
193 if (cell_code->GetValue(id_neigh) != bc) {
194 point_t P;
195 P.x = 0.5*(xn[i_neigh] + xn[i_neigh + 1]);
196 P.n = xc - xn[i_neigh];
197 vec3_t v = xn[i_neigh + 1] - xn[i_neigh];
198 v.normalise();
199 P.n -= (P.n*v)*v;
200 P.n.normalise();
201 P.idx.append(idx[i_neigh]);
202 P.idx.append(idx[i_neigh+1]);
203 P.L = computeSearchDistance(id_face);
204 points.append(P);
205 }
206 }
207 }
208 }
209 }
210 }
211 computeFeature(points, cl_pre, m_FeatureResolution2D, 2);
212 }
213 }
214
215 void UpdateDesiredMeshDensity::computeFeature3D(QVector<double> &cl_pre)
216 {
217 // add mesh density radiation to 3D feature resolution
218 //EG_STOPDATE("2015-06-01");
219
220 if (m_FeatureResolution3D < 1e-3) {
221 return;
222 }
223 CgalTriCadInterface cad(m_Grid);
224
225 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
226 vec3_t n0 = m_Part.globalNormal(id_node);
227 vec3_t x0;
228 m_Grid->GetPoint(id_node, x0.data());
229 QVector<QPair<vec3_t, vtkIdType> > inters1;
230 cad.computeIntersections(x0, n0, inters1);
231 QVector<QPair<vec3_t, vtkIdType> > inters2;
232 cad.computeIntersections(x0, (-1)*n0, inters2);
233 QVector<QPair<vec3_t, vtkIdType> > intersections = inters1 + inters2;
234 double d_min = EG_LARGE_REAL;
235 vtkIdType id_tri_min = -1;
236 for (int i = 0; i < intersections.size(); ++i) {
237 vec3_t x = intersections[i].first;
238 vtkIdType id_tri = intersections[i].second;
239 vec3_t n = GeometryTools::cellNormal(m_Grid, id_tri);
240 n.normalise();
241 vec3_t dx = x - x0;
242 if (dx*n0 < 0 && n*n0 < 0) {
243 double d = dx.abs();
244 if (d < d_min) {
245 d_min = d;
246 id_tri_min = id_tri;
247 }
248 }
249 }
250 if (id_tri_min > -1) {
251 double h = d_min/m_FeatureResolution3D;
252 int max_topo_dist = int(ceil(m_FeatureResolution3D));
253 bool topo_dist_ok = true;
254 EG_GET_CELL(id_tri_min, m_Grid);
255 for (int i = 0; i < num_pts; ++i) {
256 int topo_dist = m_Part.computeTopoDistance(pts[i], id_node, max_topo_dist, 0);
257 if (topo_dist < max_topo_dist) {
258 topo_dist_ok = false;
259 break;
260 }
261 }
262 if (topo_dist_ok) {
263 cl_pre[id_node] = h;
264 for (int i = 0; i < num_pts; ++i) {
265 cl_pre[pts[i]] = h;
266 }
267 }
268 }
269 }
270 }
271
272 void UpdateDesiredMeshDensity::readSettings()
273 {
274 readVMD();
275 QString buffer = GuiMainWindow::pointer()->getXmlSection("engrid/surface/settings").replace("\n", " ");
276 QTextStream in(&buffer, QIODevice::ReadOnly);
277 in >> m_MaxEdgeLength;
278 in >> m_MinEdgeLength;
279 in >> m_GrowthFactor;
280 in >> m_NodesPerQuarterCircle;
281 int num_bcs;
282 in >> num_bcs;
283 QVector<int> tmp_bcs;
284 GuiMainWindow::pointer()->getAllBoundaryCodes(tmp_bcs);
285 m_BoundaryCodes.clear();
286 if (num_bcs == tmp_bcs.size()) {
287 foreach (int bc, tmp_bcs) {
288 int check_state;
289 in >> check_state;
290 if (check_state == 1) {
291 m_BoundaryCodes.insert(bc);
292 }
293 }
294 }
295 if (!in.atEnd()) {
296 in >> m_FeatureResolution2D;
297 in >> m_FeatureResolution3D;
298 }
299 }
300
301 void UpdateDesiredMeshDensity::operate()
302 {
303 m_ELSManager.read();
304 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
305
306 if (m_OnlySurfaceCells) {
307 setAllSurfaceCells();
308 } else {
309 setAllCells();
310 }
311 l2g_t nodes = getPartNodes();
312 l2g_t cells = getPartCells();
313 l2l_t n2n = getPartN2N();
314
315 EG_VTKDCN(vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired");
316 EG_VTKDCN(vtkIntArray, characteristic_length_specified, m_Grid, "node_specified_density");
317
318 QVector<vec3_t> normals(cells.size(), vec3_t(0,0,0));
319 QVector<vec3_t> centres(cells.size(), vec3_t(0,0,0));
320 QVector<double> cl_pre(nodes.size(), 1e99);
321
322 computeExistingLengths();
323 if (m_BoundaryCodes.size() == 0) {
324 return;
325 }
326
327 if (m_NodesPerQuarterCircle > 1e-3) {
328 QVector<double> R(nodes.size(), 1e99);
329
330 //#pragma omp parallel for
331 for (int i_cell = 0; i_cell < cells.size(); ++i_cell) {
332 vtkIdType id_cell = cells[i_cell];
333 EG_GET_CELL(id_cell, m_Grid);
334 int bc = cell_code->GetValue(id_cell);
335 for (int i = 0; i < num_pts; ++i) {
336 int i_nodes = m_Part.localNode(pts[i]);
337 R[i_nodes] = min(R[i_nodes], fabs(GuiMainWindow::pointer()->getCadInterface(bc)->getRadius(pts[i])));
338 }
339 }
340
341 for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
342 if (cl_pre[i_nodes] == 0) {
343 EG_BUG;
344 }
345 cl_pre[i_nodes] = max(m_MinEdgeLength, min(cl_pre[i_nodes], 0.5*R[i_nodes]*M_PI/m_NodesPerQuarterCircle));
346 }
347 }
348
349 // cells across branches
350 computeFeature2D(cl_pre);
351 computeFeature3D(cl_pre);
352
353 // set everything to desired mesh density and find maximal mesh-density
354 double cl_min = 1e99;
355 int i_nodes_min = -1;
356 for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
357 vtkIdType id_node = nodes[i_nodes];
358 double cl = m_MaxEdgeLength;
359 if (m_BoundaryCodes.size() > 0) {
360 int idx = characteristic_length_specified->GetValue(id_node);
361 if ((idx >= 0) && (idx < m_VMDvector.size())) {
362 cl = m_VMDvector[idx].density;
363 }
364 }
365 if (m_Fixed[id_node]) {
366 cl = characteristic_length_desired->GetValue(id_node);
367 }
368 cl = min(cl_pre[i_nodes], cl);
369 vec3_t x;
370 m_Grid->GetPoint(id_node, x.data());
371 double cl_src = m_ELSManager.minEdgeLength(x);
372 if (cl_src > 0) {
373 cl = min(cl, cl_src);
374 }
375
376 cl = max(m_MinEdgeLength, cl);
377
378 if(cl == 0) {
379 EG_BUG;
380 }
381 characteristic_length_desired->SetValue(id_node, cl);
382
383 if (cl < cl_min) {
384 for (int i = 0; i < m_Part.n2bcGSize(id_node); ++i) {
385 if (m_BoundaryCodes.contains(m_Part.n2bcG(id_node, i))) {
386 cl_min = cl;
387 i_nodes_min = i_nodes;
388 break;
389 }
390 }
391 }
392 }
393 if (i_nodes_min == -1) {
394 EG_ERR_RETURN("There are no edges that need improving.")
395 }
396
397
398 if (m_Relaxation) {
399
400 // start from smallest characteristic length and loop as long as nodes are updated
401 int num_updated = 0;
402
403 do {
404 num_updated = 0;
405 for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
406 double cli = characteristic_length_desired->GetValue(nodes[i_nodes]);
407 if (cli <= cl_min) {
408 vec3_t xi;
409 m_Grid->GetPoint(nodes[i_nodes], xi.data());
410 for (int j = 0; j < n2n[i_nodes].size(); ++j) {
411 int j_nodes = n2n[i_nodes][j];
412 double clj = characteristic_length_desired->GetValue(nodes[j_nodes]);
413 if (clj > cli && clj > cl_min) {
414 vec3_t xj;
415 m_Grid->GetPoint(nodes[j_nodes], xj.data());
416 ++num_updated;
417 double L_new = min(m_MaxEdgeLength, cli * m_GrowthFactor);
418 if (!m_Fixed[nodes[j_nodes]]) {
419
420 double cl_min = min(characteristic_length_desired->GetValue(nodes[j_nodes]), L_new);
421 if(cl_min==0) {
422 qWarning()<<"m_MaxEdgeLength="<<m_MaxEdgeLength;
423 qWarning()<<"cli="<<cli;
424 qWarning()<<"m_GrowthFactor="<<m_GrowthFactor;
425 qWarning()<<"characteristic_length_desired->GetValue(nodes[j_nodes])="<<characteristic_length_desired->GetValue(nodes[j_nodes]);
426 qWarning()<<"L_new="<<L_new;
427 EG_BUG;
428 }
429 characteristic_length_desired->SetValue(nodes[j_nodes], cl_min);
430
431 }
432 }
433 }
434 }
435 }
436 cl_min *= m_GrowthFactor;
437 } while (num_updated > 0);
438 }
439 }
This is a mirror repository of the new location for engrid@github