Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #include <fstream>
34 : #include <queue>
35 : #include "lib_grid/lg_base.h"
36 : #include "common/profiler/profiler.h"
37 : #include "simple_grid.h"
38 : #include "edge_length_adjustment.h"
39 : #include "lib_grid/refinement/regular_refinement.h"
40 : #include "common/node_tree/node_tree.h"
41 : #include "lib_grid/algorithms/trees/octree.h"
42 : #include "lib_grid/callbacks/callbacks.h"
43 : #include "../trees/octree.h"
44 :
45 : using namespace std;
46 :
47 : namespace ug
48 : {
49 :
50 : /// only for debugging purposes!!!
51 : /** Output value pairs to gnuplot...
52 : * \{ */
53 : //#define EDGE_LENGTH_ADJUSTMENT__GPLOT_ENABLED
54 : #ifdef EDGE_LENGTH_ADJUSTMENT__GPLOT_ENABLED
55 : typedef vector<pair<number, number> > GnuplotData;
56 : static GnuplotData gplotLengthFac;
57 : static GnuplotData gplotMinCurvature;
58 : static GnuplotData gplotAverageCurvature;
59 :
60 : void WriteGnuplotData(const char* filename, const GnuplotData& data)
61 : {
62 : ofstream out(filename);
63 : if(!out)
64 : return;
65 :
66 : for(size_t i = 0; i < data.size(); ++i)
67 : out << data[i].first << " " << data[i].second << endl;
68 :
69 : out.close();
70 : }
71 :
72 : #define GPLOTPOINT(dataName, x, y) dataName.push_back(make_pair<number, number>((x), (y)));
73 : #define GPLOTSAVE() {WriteGnuplotData("length_fac.gplot", gplotLengthFac);\
74 : WriteGnuplotData("min_curvature.gplot", gplotMinCurvature);\
75 : WriteGnuplotData("average_curvature.gplot", gplotAverageCurvature);}
76 : #else
77 : // do nothing if EDGE_LENGTH_ADJUSTMENT__GPLOT_ENABLED is false
78 : #define GPLOTPOINT(dataName, x, y)
79 : #define GPLOTSAVE()
80 : #endif
81 : /** \} */
82 :
83 :
84 : /*
85 : vector3 PNTrianglePos(const vector3& p0, const vector3& p1, const vector3& p2,
86 : const vector3& n0, const vector3& n1, const vector3& n2);
87 :
88 : vector3 PNTriangleNorm(const vector3& p0, const vector3& p1, const vector3& p2,
89 : const vector3& n0, const vector3& n1, const vector3& n2);
90 :
91 : vector3 PNCTrianglePos(const vector3& p0, const vector3& p1, const vector3& p2,
92 : const vector3& n0, const vector3& n1, const vector3& n2,
93 : const vector3& cn0, const vector3& cn1, const vector3& cn2);
94 :
95 : vector3 PNCTriangleNorm(const vector3& p0, const vector3& p1, const vector3& p2,
96 : const vector3& n0, const vector3& n1, const vector3& n2,
97 : const vector3& cn0, const vector3& cn1, const vector3& cn2);
98 :
99 :
100 : vector3 pos(number bc0, number bc1);
101 : vector3 norm(number bc0, number bc1);
102 : };
103 :
104 : class ProjectedVertex{
105 : Vertex* vertex();
106 : vector3 vertex_position();
107 : vector3 vertex_normal();
108 : vector3 surface_normal();
109 : vector3 surface_position();
110 : size_t num_barycentric_coords();
111 : number barycentric_coord(size_t index);
112 : };
113 :
114 : class SurfaceRepresentation{
115 : public:
116 :
117 : };
118 : */
119 :
120 : ////////////////////////////////////////////////////////////////////////
121 0 : static void AssignFixedVertices(Grid& grid, SubsetHandler& shMarks)
122 : {
123 0 : grid.begin_marking();
124 :
125 : // mark all vertices that are not regular crease-vertices as fixed
126 : for(EdgeIterator iter = shMarks.begin<Edge>(REM_CREASE);
127 0 : iter != shMarks.end<Edge>(REM_CREASE); ++iter)
128 : {
129 : Edge* e = *iter;
130 0 : for(size_t i = 0; i < 2; ++i){
131 0 : Vertex* vrt = e->vertex(i);
132 0 : if(!grid.is_marked(vrt)){
133 : grid.mark(vrt);
134 : int counter = 0;
135 0 : for(Grid::AssociatedEdgeIterator nbIter = grid.associated_edges_begin(vrt);
136 0 : nbIter != grid.associated_edges_end(vrt); ++nbIter)
137 : {
138 0 : if(shMarks.get_subset_index(*nbIter) != REM_NONE)
139 0 : ++counter;
140 : }
141 :
142 0 : if(counter != 2)
143 0 : shMarks.assign_subset(vrt, REM_FIXED);
144 : }
145 : }
146 : }
147 :
148 0 : grid.end_marking();
149 :
150 : // mark all vertices that lie on a fixed edge as fixed vertex
151 : for(EdgeIterator iter = shMarks.begin<Edge>(REM_FIXED);
152 0 : iter != shMarks.end<Edge>(REM_FIXED); ++iter)
153 : {
154 : Edge* e = *iter;
155 0 : shMarks.assign_subset(e->vertex(0), REM_FIXED);
156 0 : shMarks.assign_subset(e->vertex(1), REM_FIXED);
157 : }
158 0 : }
159 :
160 0 : static void AssignCreaseVertices(Grid& grid, SubsetHandler& shMarks)
161 : {
162 : // mark all vertices that lie on a crease and which are not fixed
163 : // as crease vertices.
164 0 : if(shMarks.num_subsets() <= REM_CREASE)
165 : return;
166 :
167 : for(EdgeIterator iter = shMarks.begin<Edge>(REM_CREASE);
168 0 : iter != shMarks.end<Edge>(REM_CREASE); ++iter)
169 : {
170 : Edge* e = *iter;
171 0 : for(uint i = 0; i < 2; ++i)
172 0 : if(shMarks.get_subset_index(e->vertex(i)) != REM_FIXED)
173 0 : shMarks.assign_subset(e->vertex(i), REM_CREASE);
174 : }
175 : }
176 :
177 :
178 : ////////////////////////////////////////////////////////////////////////
179 : template <class TVertexPositionAccessor>
180 : number CalculateNormalDot(TriangleDescriptor& td1, TriangleDescriptor& td2,
181 : TVertexPositionAccessor& aaPos)
182 : {
183 : vector3 n1;
184 : CalculateTriangleNormal(n1, aaPos[td1.vertex(0)],
185 : aaPos[td1.vertex(1)], aaPos[td1.vertex(2)]);
186 : vector3 n2;
187 : CalculateTriangleNormal(n2, aaPos[td2.vertex(0)],
188 : aaPos[td2.vertex(1)], aaPos[td2.vertex(2)]);
189 : return VecDot(n1, n2);
190 : }
191 :
192 :
193 : ////////////////////////////////////////////////////////////////////////
194 : // CalculateCurvature
195 : template <class TAAPosVRT>
196 0 : number CalculateMinCurvature(Grid& grid, SubsetHandler& shMarks,
197 : Vertex* vrt, TAAPosVRT& aaPos)
198 : {
199 : //TODO: check whether static vNormals brings any benefits.
200 : //TODO: special cases for crease vertices
201 :
202 : // face normals
203 0 : static vector<vector3> vNormals;
204 : // vertex normal (mean face normal)
205 : vector3 n(0, 0, 0);
206 :
207 : // calculate the normals of associated faces
208 : vNormals.clear();
209 0 : Grid::AssociatedFaceIterator iterEnd = grid.associated_faces_end(vrt);
210 0 : for(Grid::AssociatedFaceIterator iter = grid.associated_faces_begin(vrt);
211 0 : iter != iterEnd; ++iter)
212 : {
213 : vector3 nTmp;
214 0 : CalculateNormal(nTmp, *iter, aaPos);
215 0 : vNormals.push_back(nTmp);
216 : VecAdd(n, n, nTmp);
217 : }
218 :
219 : // the vertex normal
220 0 : VecNormalize(n, n);
221 :
222 : // get the min dot-product of the vertex normal with associated faces-normals
223 0 : number minDot = 1;
224 0 : for(size_t i = 0; i < vNormals.size(); ++i)
225 0 : minDot = std::min(minDot, VecDot(n, vNormals[i]));
226 :
227 : //todo: Think about converting to radiants.
228 : //minDot = acos(minDot) / (0.5 * PI);
229 : //...
230 : // done
231 : GPLOTPOINT(gplotMinCurvature, 0, minDot);
232 0 : return minDot;
233 : }
234 :
235 : ////////////////////////////////////////////////////////////////////////
236 : template <class TAAPosVRT>
237 0 : number CalculateAverageCurvature(Grid& grid, SubsetHandler& shMarks,
238 : Edge* e, TAAPosVRT& aaPos)
239 : {
240 0 : number avCurv = 0.5 * (CalculateMinCurvature(grid, shMarks, e->vertex(0), aaPos)
241 0 : + CalculateMinCurvature(grid, shMarks, e->vertex(1), aaPos));
242 : GPLOTPOINT(gplotAverageCurvature, 0.5, avCurv);
243 0 : return avCurv;
244 : }
245 :
246 : ////////////////////////////////////////////////////////////////////////
247 : template <class TAAPosVRT>
248 : number CalculateLengthFac(Grid& grid, SubsetHandler& shMarks,
249 : Edge* e, TAAPosVRT& aaPos)
250 : {
251 0 : number lenFac = CalculateAverageCurvature(grid, shMarks, e, aaPos);
252 0 : lenFac = (lenFac - 0.95) / 0.05;
253 0 : lenFac = max(number(0.25), lenFac);
254 : GPLOTPOINT(gplotLengthFac, 0.5, lenFac);
255 : return lenFac;
256 : }
257 :
258 : ////////////////////////////////////////////////////////////////////////
259 : ////////////////////////////////////////////////////////////////////////
260 : template <class TAAPosVRT, class TAANormVRT, class TAAIntVRT>
261 0 : bool TrySwap(Grid& grid, Edge* e, TAAPosVRT& aaPos, TAANormVRT& aaNorm,
262 : TAAIntVRT& aaInt, SubsetHandler* pshMarks = NULL,
263 : EdgeSelector* pCandidates = NULL)
264 : {
265 : // swaps are neither allowed for crease edges nor for fixed edges
266 0 : if(pshMarks){
267 0 : if(pshMarks->get_subset_index(e) == REM_FIXED ||
268 : pshMarks->get_subset_index(e) == REM_CREASE)
269 0 : return false;
270 : }
271 :
272 : // get the associated faces. we need two of them
273 : Face* f[2];
274 0 : if(GetAssociatedFaces(f, grid, e, 2) != 2)
275 : return false;
276 :
277 : // make sure both faces are triangles
278 0 : if((f[0]->num_vertices() != 3) || (f[1]->num_vertices() != 3))
279 0 : return false;
280 :
281 : // create a simple grid
282 : SimpleGrid sg;
283 0 : if(!ObtainSimpleGrid(sg, grid, e->vertex(0), e->vertex(1), 0,
284 : aaPos, aaNorm, aaInt))
285 : {
286 : LOG("ObtainSimpleGrid failed. ignoring edge...\n");
287 : return false;
288 : }
289 :
290 : // calculate geometric-approximation-degree and triangle quality
291 : //number approxDeg = GeometricApproximationDegree(sg);
292 0 : number shapeDeg = ShapeQualityDegree(sg);
293 :
294 : number smoothDeg = VecDot(sg.triangleNormals[0], sg.triangleNormals[1]);
295 :
296 : // this is a new test. the idea is that each edge should be orthogonal to
297 : // the normals of its endpoints - at least in a perfectly smooth surface.
298 :
299 : number approxDeg;
300 : number newApproxDeg;
301 : {
302 : vector3 v;
303 : VecSubtract(v, sg.vertices[1], sg.vertices[0]);
304 0 : VecNormalize(v, v);
305 0 : approxDeg = 1. - 0.5*(fabs(VecDot(v, sg.vertexNormals[0])) +
306 0 : fabs(VecDot(v, sg.vertexNormals[1])));
307 :
308 : VecSubtract(v, sg.vertices[3], sg.vertices[2]);
309 0 : VecNormalize(v, v);
310 0 : newApproxDeg = 1. - 0.5*(fabs(VecDot(v, sg.vertexNormals[2])) +
311 0 : fabs(VecDot(v, sg.vertexNormals[3])));
312 : }
313 :
314 : // perform a swap on the simple grid
315 0 : if(!SwapEdge(sg))
316 : {
317 : LOG("swap edge failed...\n");
318 : return false;
319 : }
320 :
321 : // calculate new geometric-approximation-degree and triangle quality
322 : //number newApproxDeg = GeometricApproximationDegree(sg);
323 0 : number newShapeDeg = ShapeQualityDegree(sg);
324 : number newSmoothDeg = VecDot(sg.triangleNormals[0], sg.triangleNormals[1]);
325 :
326 : // neither the shapeDeg nor the approxDeg may get too bad.
327 0 : if((newApproxDeg < 0.5 * approxDeg) || (newShapeDeg < 0.5 * shapeDeg))
328 : return false;
329 :
330 : // make sure that the swap does not destroy the smoothness
331 0 : if(newSmoothDeg < 0.1 * smoothDeg)
332 : return false;
333 :
334 0 : if(0.2 * (newApproxDeg - approxDeg) + 0.8 * (newShapeDeg - shapeDeg) > 0)
335 : //if(newShapeDeg > shapeDeg)
336 : //if(newApproxDeg > approxDeg)
337 : {
338 : // swap the edge
339 0 : e = SwapEdge(grid, e);
340 :
341 0 : if(e){
342 : // swap was successful
343 : // if pCandidates was specified then add new candidates
344 0 : if(pCandidates){
345 0 : for(int i = 0; i < 2; ++i)
346 0 : pCandidates->select(grid.associated_edges_begin(e->vertex(i)),
347 0 : grid.associated_edges_end(e->vertex(i)));
348 : // e was selected but is not really a candidate
349 0 : pCandidates->deselect(e);
350 : }
351 :
352 0 : return true;
353 : }
354 : }
355 :
356 : return false;
357 0 : }
358 :
359 : ////////////////////////////////////////////////////////////////////////
360 : template <class TAAPosVRT, class TAANormVRT, class TAAIntVRT>
361 0 : bool PerformSwaps(Grid& grid, SubsetHandler& shMarks, EdgeSelector& esel,
362 : TAAPosVRT& aaPos, TAANormVRT& aaNorm, TAAIntVRT& aaInt)
363 : {
364 : PROFILE_FUNC();
365 : LOG(" performing swaps\n");
366 : int numSwaps = 0;
367 0 : int maxNumSwaps = esel.num<Edge>() * 2;
368 :
369 0 : while(!esel.empty())
370 : {
371 : Edge* e = *esel.begin<Edge>();
372 0 : esel.deselect(e);
373 :
374 0 : if(TrySwap(grid, e, aaPos, aaNorm, aaInt, &shMarks, &esel)){
375 0 : ++numSwaps;
376 0 : if(numSwaps > maxNumSwaps){
377 : UG_LOG(" aborting since maxNumSwaps was reached...");
378 0 : esel.clear();
379 : }
380 : }
381 : }
382 0 : LOG(" swaps performed: " << numSwaps << endl);
383 :
384 0 : return true;
385 : }
386 :
387 : /** returns the resulting vertex or NULL, if no collapse was performed.*/
388 : template <class TAAPosVRT, class TAANormVRT, class TAAIntVRT>
389 0 : Vertex* TryCollapse(Grid& grid, Edge* e,
390 : TAAPosVRT& aaPos, TAANormVRT& aaNorm,
391 : TAAIntVRT& aaInt, SubsetHandler* pshMarks = NULL,
392 : EdgeSelector* pCandidates = NULL)
393 : {
394 0 : if(pshMarks)
395 : {
396 : SubsetHandler& shMarks = *pshMarks;
397 : // collapses are not allowed for fixed edges
398 0 : if(shMarks.get_subset_index(e) == REM_FIXED)
399 : return NULL;
400 :
401 : // if both endpoints of are fixed vertices then
402 : // we may not collapse
403 : int vrtSI[2];
404 0 : vrtSI[0] = shMarks.get_subset_index(e->vertex(0));
405 0 : vrtSI[1] = shMarks.get_subset_index(e->vertex(1));
406 :
407 0 : if((vrtSI[0] == REM_FIXED) && (vrtSI[1] == REM_FIXED))
408 : return NULL;
409 :
410 : // if both endpoints are somehow marked, e has to be a
411 : // crease edge
412 0 : if((vrtSI[0] != REM_NONE) && (vrtSI[1] != REM_NONE)
413 0 : && (shMarks.get_subset_index(e) != REM_CREASE))
414 : return NULL;
415 : }
416 :
417 : // check whether the edge can be collapsed
418 0 : if(EdgeCollapseIsValid(grid, e))
419 : {
420 : // test the collapse using a simple-grid
421 : SimpleGrid sgSrc;
422 0 : if(!ObtainSimpleGrid(sgSrc, grid, e->vertex(0), e->vertex(1), 1,
423 : aaPos, aaNorm, aaInt))
424 : {
425 : LOG("ObtainSimpleGrid failed. ignoring edge...\n");
426 : return NULL;
427 : }
428 :
429 : // calculate geometric-approximation-degree and triangle quality
430 0 : number approxDeg = GeometricApproximationDegree(sgSrc);
431 0 : number shapeDeg = ShapeQualityDegree(sgSrc);
432 :
433 : // perform a collapse on the simple grid
434 : SimpleGrid sgDest;
435 0 : if(!ObtainSimpleGrid_CollapseEdge(sgDest, grid, e,
436 : 1, aaPos, aaNorm, aaInt))
437 : {
438 : LOG("collapse edge failed...\n");
439 : return NULL;
440 : }
441 :
442 : // get the positions of the old endpoints
443 : static const int numTestPositions = 3;
444 : int newInd = 0;
445 0 : vector3 v[numTestPositions];
446 0 : v[0] = aaPos[e->vertex(0)];
447 0 : v[1] = aaPos[e->vertex(1)];
448 : v[2] = sgDest.vertices[newInd];
449 :
450 : // we'll compare 3 approximation degrees and three shape degrees
451 : number newApproxDeg[numTestPositions];
452 : number newShapeDeg[numTestPositions];
453 :
454 : // check which position is the best
455 : int bestIndex = -1;
456 : // the vertex subset index is used to support marks (crease and fixed vertices)
457 : int vrtSI[2];
458 : vrtSI[0] = vrtSI[1] = REM_NONE;
459 :
460 0 : if(pshMarks)
461 : {
462 0 : vrtSI[0] = pshMarks->get_subset_index(e->vertex(0));
463 0 : vrtSI[1] = pshMarks->get_subset_index(e->vertex(1));
464 :
465 0 : if((vrtSI[0] == REM_FIXED) || ((vrtSI[0] != REM_NONE) && (vrtSI[1] == REM_NONE))){
466 : bestIndex = 0;
467 : sgDest.vertices[newInd] = v[0];
468 0 : newApproxDeg[0] = GeometricApproximationDegree(sgDest);
469 0 : newShapeDeg[0] = ShapeQualityDegree(sgDest);
470 : }
471 0 : else if((vrtSI[1] == REM_FIXED) || ((vrtSI[1] != REM_NONE) && (vrtSI[0] == REM_NONE))){
472 : bestIndex = 1;
473 : sgDest.vertices[newInd] = v[1];
474 0 : newApproxDeg[1] = GeometricApproximationDegree(sgDest);
475 0 : newShapeDeg[1] = ShapeQualityDegree(sgDest);
476 : }
477 : }
478 :
479 : if(bestIndex == -1){
480 : // check all three approximation degrees
481 0 : for(int i = 0; i < numTestPositions; ++i){
482 : // we'll compute all qualities with the averaged normal
483 : sgDest.vertices[newInd] = v[i];
484 0 : CalculateTriangleNormals(sgDest);
485 0 : newApproxDeg[i] = GeometricApproximationDegree(sgDest);
486 0 : newShapeDeg[i] = ShapeQualityDegree(sgDest);
487 : }
488 : // get the best one
489 : bestIndex = 0;
490 : /*
491 : for(int i = 1; i < numTestPositions; ++i){
492 : if(newApproxDeg[i] > newApproxDeg[bestIndex])
493 : bestIndex = i;
494 : }
495 : */
496 0 : for(int i = 1; i < numTestPositions; ++i){
497 0 : if(newShapeDeg[i] > newShapeDeg[bestIndex])
498 : bestIndex = i;
499 : }
500 : }
501 :
502 : // if the shape-degree of the collapsed region is too bad, we'll skip the collapse
503 0 : if(newShapeDeg[bestIndex] < 0.5 * shapeDeg)
504 : return NULL;
505 : /*
506 : // the approximation degree is only interesting if both endpoints of the
507 : // edge are regular surface vertices
508 : bool regularNeighbourhood = IsRegularSurfaceVertex(grid, e->vertex(0)) &&
509 : IsRegularSurfaceVertex(grid, e->vertex(1));
510 : */
511 :
512 : // if the best approximation degree is not too bad, we'll perform the collapse
513 0 : if(/*!regularNeighbourhood || */(newApproxDeg[bestIndex] > 0.95 * approxDeg))
514 : {
515 : // pick one of the endpoints to be the one that resides
516 : // This has to be done with care, since the residing vertex
517 : // determines which edges will be removed and which will remain
518 : // (this is important, since we want crease edges to remain!).
519 : // If only one of the endpoints lies on a crease or is fixed, the
520 : // decision is easy (the marked one has to remain).
521 : // If both are marked, we have to investigate the triangles which
522 : // are adjacent to the collapsed edge (quadrilaterals don't make
523 : // problems here).
524 : // - If all three corner vertices are marked we have to distinguish:
525 : // (a) all three edges are marked -> abort
526 : // (b) two edges are marked -> the vertex connecting both remains.
527 : // (c) one edge is marked -> the chosen vertex doesn't affect edge-marks.
528 : // - If only the two are marked they don't affect the edge-marks.
529 : //
530 : // Even more care has to be taken if a fixed vertex is involved.
531 : // If case (b) applies and the creas-vertex is has to remain, it
532 : // has to be moved to the position of the fixed-vertex and has to
533 : // be marked fixed itself.
534 :
535 : // choose the vertex that shall remain.
536 0 : Vertex* vrt = e->vertex(0);
537 :
538 0 : if(vrtSI[0] != REM_FIXED && vrtSI[1] != REM_NONE)
539 : {
540 0 : vrt = e->vertex(1);
541 : }
542 :
543 0 : if(vrtSI[0] != REM_NONE && vrtSI[1] != REM_NONE){
544 : // both are marked. Things are getting complicated now!
545 : // get adjacent faces of e
546 : vector<Face*> vFaces;
547 0 : vFaces.reserve(2);
548 0 : CollectFaces(vFaces, grid, e);
549 :
550 : vector<Edge*> vEdges;
551 0 : vEdges.reserve(4);
552 0 : for(size_t i = 0; i < vFaces.size(); ++i){
553 0 : Face* f = vFaces[i];
554 : // only triangles are interesting
555 0 : if(f->num_edges() == 3){
556 : // get the number of marked edges
557 0 : CollectEdges(vEdges, grid, f);
558 : // count the marked edges
559 : int numMarked = 0;
560 0 : for(size_t j = 0; j < vEdges.size(); ++j){
561 : // note that pshMarks exists since vrtSI != REM_NONE
562 0 : if(pshMarks->get_subset_index(vEdges[j]) != REM_NONE)
563 0 : ++numMarked;
564 : }
565 :
566 0 : if(numMarked == 3){
567 : return NULL; // case (a) applies
568 : }
569 0 : else if(numMarked == 2){
570 : // check which of the vrts is connected to two
571 : // marked edges of the triangle
572 0 : for(size_t j = 0; j < 2; ++j){
573 : int numMarked = 0;
574 0 : for(size_t k = 0; k < vEdges.size(); ++k){
575 0 : if(pshMarks->get_subset_index(vEdges[k]) != REM_NONE){
576 0 : if(EdgeContains(vEdges[k], e->vertex(j)))
577 0 : ++numMarked;
578 : }
579 : }
580 : // if numMarked == 2 we found it
581 0 : if(numMarked == 2){
582 : // the connected edge has to be marked as a crease
583 0 : Edge* ce = GetConnectedEdge(grid, e->vertex(j), f);
584 0 : if(ce)
585 0 : pshMarks->assign_subset(ce, REM_CREASE);
586 : // we're done. break
587 : break;
588 : }
589 : }
590 : }
591 : else{
592 : }
593 : }
594 : }
595 0 : }
596 :
597 : // collapse the edge
598 0 : CollapseEdge(grid, e, vrt);
599 :
600 : // assign best position
601 : aaPos[vrt] = v[bestIndex];
602 : // assign the normal
603 : aaNorm[vrt] = sgDest.vertexNormals[newInd];
604 : /*
605 : if(pCandidates){
606 : //TODO: all edges that belong to associated faces are new candidates.
607 : // associated edges of vrt are possible new candidates
608 : pCandidates->select(grid.associated_edges_begin(vrt),
609 : grid.associated_edges_end(vrt));
610 : }*/
611 :
612 0 : return vrt;
613 : }
614 0 : }
615 :
616 : return NULL;
617 : }
618 :
619 : ////////////////////////////////////////////////////////////////////////
620 : template <class TAAPosVRT, class TAANormVRT, class TAAIntVRT>
621 0 : bool PerformCollapses(Grid& grid, SubsetHandler& shMarks, EdgeSelector& esel,
622 : number minEdgeLen, TAAPosVRT& aaPos, TAANormVRT& aaNorm,
623 : TAAIntVRT& aaInt, bool adaptive = true)
624 : {
625 : PROFILE_FUNC();
626 : LOG(" performing collapses\n");
627 : vector<Edge*> edges;
628 : int numCollapses = 0;
629 : // compare squares
630 0 : minEdgeLen *= minEdgeLen;
631 :
632 0 : while(!esel.empty())
633 : {
634 : Edge* e = *esel.begin<Edge>();
635 0 : esel.deselect(e);
636 :
637 : // the higher the curvature the smaller the maxEdgeLen.
638 : // minimal lenFac is 0.1
639 : number lenFac = 1.0;
640 0 : if(adaptive)
641 : lenFac = CalculateLengthFac(grid, shMarks, e, aaPos);
642 :
643 : // check whether the edge is short enough
644 0 : if(VecDistanceSq(aaPos[e->vertex(0)], aaPos[e->vertex(1)]) < lenFac * minEdgeLen)
645 : {
646 0 : Vertex* vrt = TryCollapse(grid, e, aaPos, aaNorm, aaInt, &shMarks, &esel);
647 0 : if(vrt){
648 0 : ++numCollapses;
649 :
650 : // we'll deselect associated edges of vrt, to avoid cascade-collapses
651 : CollectAssociated(edges, grid, vrt);
652 : esel.deselect(edges.begin(), edges.end());
653 : }
654 : }
655 : }
656 0 : LOG(" collapses performed: " << numCollapses << endl);
657 :
658 0 : return true;
659 0 : }
660 :
661 : ////////////////////////////////////////////////////////////////////////
662 : template <class TAAPosVRT, class TAANormVRT>
663 0 : bool TrySplit(Grid& grid, Edge* e, TAAPosVRT& aaPos, TAANormVRT& aaNorm,
664 : EdgeSelector* pCandidates = NULL, SubsetHandler* pshMarks = NULL)
665 : {
666 : bool bCreaseEdge = false;
667 : // splits are not allowed for fixed edges
668 0 : if(pshMarks){
669 0 : if(pshMarks->get_subset_index(e) == REM_FIXED)
670 : return false;
671 0 : else if(pshMarks->get_subset_index(e) == REM_CREASE)
672 : bCreaseEdge = true;
673 : }
674 :
675 : // get the center of the edges
676 0 : vector3 vCenter = CalculateCenter(e, aaPos);
677 :
678 : // the new normal
679 : vector3 n;
680 0 : if(bCreaseEdge){
681 : // interpolating the normal can cause severe problems at highly
682 : // irregular vertices or if one vertecs lies on a very
683 : // sharp edge (the normals of the endpoints thus point
684 : // in different directions.)
685 : // VecAdd(n, aaNorm[e->vertex(0)], aaNorm[e->vertex(1)]);
686 : // VecNormalize(n, n);
687 0 : CalculateNormal(n, grid, e, aaPos);
688 : }
689 :
690 : // split the edge
691 : RegularVertex* vrt = SplitEdge<RegularVertex>(grid, e, false);
692 :
693 : // assign the new position
694 : aaPos[vrt] = vCenter;
695 :
696 : // assign the new normal. calculate it if required
697 0 : if(!bCreaseEdge)
698 0 : CalculateVertexNormal(n, grid, vrt, aaPos);
699 :
700 : aaNorm[vrt] = n;
701 :
702 : // associated edges of vrt are candidates again
703 0 : if(pCandidates)
704 0 : pCandidates->select(grid.associated_edges_begin(vrt),
705 : grid.associated_edges_end(vrt));
706 :
707 : return true;
708 : }
709 :
710 : ////////////////////////////////////////////////////////////////////////
711 : template <class TAAPosVRT, class TAANormVRT>
712 0 : bool PerformSplits(Grid& grid, SubsetHandler& shMarks, EdgeSelector& esel,
713 : number maxEdgeLen, TAAPosVRT& aaPos, TAANormVRT& aaNorm,
714 : bool adaptive = true)
715 : {
716 : PROFILE_FUNC();
717 : // compare squares
718 0 : maxEdgeLen *= maxEdgeLen;
719 :
720 : LOG(" performing splits\n");
721 : int numSplits = 0;
722 :
723 0 : while(!esel.empty())
724 : {
725 : // get an edge
726 : Edge* e = *esel.begin<Edge>();
727 0 : esel.deselect(e);
728 :
729 : // the higher the curvature the smaller the maxEdgeLen.
730 : // minimal lenFac is 0.1
731 : number lenFac = 1.0;
732 0 : if(adaptive)
733 : lenFac = CalculateLengthFac(grid, shMarks, e, aaPos);
734 :
735 : // check whether the edge should be splitted
736 0 : if(VecDistanceSq(aaPos[e->vertex(0)], aaPos[e->vertex(1)]) > lenFac * maxEdgeLen)
737 : {
738 0 : if(TrySplit(grid, e, aaPos, aaNorm, &esel, &shMarks))
739 0 : ++numSplits;
740 : }
741 : }
742 :
743 0 : LOG(" splits performed: " << numSplits << endl);
744 0 : return true;
745 : }
746 :
747 : ////////////////////////////////////////////////////////////////////////
748 : // relocate point by smoothing
749 : // static void RelocatePointBySmoothing(vector3& vOut, const vector3&v,
750 : // std::vector<vector3>& vNodes,
751 : // size_t numIterations, number stepSize,
752 : // std::vector<number>* pvWeights = NULL)
753 : // {
754 : // if(vNodes.empty())
755 : // return;
756 :
757 : // //TODO: incorporate weights
758 : // // iterate through all nodes, calculate the direction
759 : // // from v to each node, and add the scaled sum to v.
760 : // vector3 t, vOld;
761 : // vOut = v;
762 : // stepSize /= (number)vNodes.size();
763 :
764 : // for(size_t j = 0; j < numIterations; ++j){
765 : // vOld = vOut;
766 : // for(size_t i = 0; i < vNodes.size(); ++i){
767 : // VecSubtract(t, vNodes[i], vOld);
768 : // VecScale(t, t, stepSize);
769 : // VecAdd(vOut, vOut, t);
770 : // }
771 : // }
772 : // }
773 :
774 : ////////////////////////////////////////////////////////////////////////
775 : // FixBadTriangles
776 : // template <class TAAPosVRT, class TAANormVRT>
777 : // static bool FixBadTriangles(Grid& grid, SubsetHandler& shMarks, EdgeSelector& esel,
778 : // TAAPosVRT& aaPos, TAANormVRT& aaNorm,
779 : // number qualityThreshold)
780 : // {
781 : // LOG(" fixing bad triangles... ");
782 : // // iterate through marked edges and check whether adjacent triangles are
783 : // // badly shaped. If this is the case, then split the non-marked
784 : // // edges of the triangle.
785 : // // store newly inserted vertices and smooth them at the end of the algo
786 : // vector<Vertex*> vNewVrts;
787 : // vector<Face*> vFaces;
788 : // vector<Edge*> vEdges;
789 :
790 : // // we wont assign this selector to the grid until it is clear that
791 : // // we'll need it.
792 : // Selector sel;
793 :
794 :
795 : // EdgeIterator iter = esel.begin<Edge>();
796 : // while(iter != esel.end<Edge>())
797 : // {
798 : // // store edge and increase iterator immediatly
799 : // Edge* e = *iter;
800 : // ++iter;
801 :
802 : // // get the adjacent faces
803 : // CollectFaces(vFaces, grid, e);
804 :
805 : // // check whether one of them is degenerated.
806 : // for(size_t i = 0; i < vFaces.size(); ++i){
807 : // Face* f = vFaces[i];
808 : // number q = FaceQuality(f, aaPos);
809 : // // if the quality is smaller than threshold, we have to
810 : // // do something
811 : // if(q < qualityThreshold){
812 : // // make sure that the selector is connected to the grid
813 : // if(sel.grid() == NULL)
814 : // sel.assign_grid(grid);
815 :
816 : // // get associated edges and mark them for refinement
817 : // CollectEdges(vEdges, grid, f);
818 :
819 : // for(size_t j = 0; j < vEdges.size(); ++j){
820 : // if(shMarks.get_subset_index(vEdges[j]) == REM_NONE)
821 : // sel.select(vEdges[j]);
822 : // }
823 : // }
824 : // }
825 : // }
826 :
827 : // // if edges have been selected, we'll now call refine.
828 : // if(sel.grid() != NULL){
829 : // if(sel.num<Edge>() > 0){
830 : // if(Refine(grid, sel)){
831 : // LOG(sel.num<Vertex>() << " new vertices... ");
832 : // // retriangulate surface
833 : // if(grid.num<Quadrilateral>() > 0)
834 : // Triangulate(grid, grid.begin<Quadrilateral>(), grid.end<Quadrilateral>());
835 :
836 : // // calculate normals, then
837 : // // smooth new vertices (all vertices selected in sel).
838 : // vector<Vertex*> vNeighbours;
839 : // vector<vector3> vNodes;
840 :
841 : // // calculate normals
842 : // for(VertexIterator iter = sel.begin<Vertex>();
843 : // iter != sel.end<Vertex>(); ++iter)
844 : // {
845 : // // collect neighbour nodes
846 : // Vertex* vrt = *iter;
847 : // CollectNeighbors(vNeighbours, grid, vrt);
848 :
849 : // // sum their normals and interpolate it
850 : // // make sure to only add valid normals
851 : // vector3 n(0, 0, 0);
852 : // for(size_t i = 0; i < vNeighbours.size(); ++i){
853 : // if(!sel.is_selected(vNeighbours[i]))
854 : // VecAdd(n, n, aaNorm[vNeighbours[i]]);
855 : // }
856 : // VecNormalize(aaNorm[vrt], n);
857 : // }
858 :
859 : // // repeat smoothing.
860 : // for(size_t i = 0; i < 5; ++i){
861 : // for(VertexIterator iter = sel.begin<Vertex>();
862 : // iter != sel.end<Vertex>(); ++iter)
863 : // {
864 : // Vertex* vrt = *iter;
865 :
866 : // // collect the neighbours and project them to the plane
867 : // // that is defined by vrt and its normal
868 : // vector3 v = aaPos[vrt];
869 : // vector3 n = aaNorm[vrt];
870 :
871 : // CollectNeighbors(vNeighbours, grid, vrt);
872 : // vNodes.resize(vNeighbours.size());
873 :
874 : // for(size_t j = 0; j < vNodes.size(); ++j)
875 : // ProjectPointToPlane(vNodes[j], aaPos[vNeighbours[j]], v, n);
876 :
877 : // // perform point relocation
878 : // RelocatePointBySmoothing(aaPos[vrt], v, vNodes, 5, 0.05);
879 : // }
880 : // }
881 : // }
882 : // else{
883 : // LOG("refine failed!\n");
884 : // return false;
885 : // }
886 : // }
887 : // }
888 :
889 : // LOG("done\n");
890 : // return true;
891 : // }
892 :
893 : // ////////////////////////////////////////////////////////////////////////
894 : // // PerformSmoothing
895 : // template <class TAAPosVRT, class TAANormVRT>
896 : // static void PerformSmoothing(Grid& grid, SubsetHandler& shMarks,
897 : // TAAPosVRT& aaPos, TAANormVRT& aaNorm,
898 : // size_t numIterations, number stepSize)
899 : // {
900 : // vector<vector3> vNodes;
901 : // vector<Vertex*> vNeighbours;
902 : // for(size_t i = 0; i < numIterations; ++i){
903 : // CalculateVertexNormals(grid, aaPos, aaNorm);
904 : // for(VertexIterator iter = grid.begin<Vertex>();
905 : // iter != grid.end<Vertex>(); ++iter)
906 : // {
907 : // Vertex* vrt = *iter;
908 : // // if the vertex is fixed then leave it where it is.
909 : // if(shMarks.get_subset_index(vrt) == REM_FIXED)
910 : // continue;
911 : // /*
912 : // if(shMarks.get_subset_index(vrt) == REM_CREASE)
913 : // continue;
914 : // */
915 : // // collect the neighbours and project them to the plane
916 : // // that is defined by vrt and its normal
917 : // vector3 v = aaPos[vrt];
918 : // vector3 n = aaNorm[vrt];
919 :
920 : // bool bProjectPointsToPlane = true;
921 :
922 : // if(shMarks.get_subset_index(vrt) == REM_CREASE){
923 : // CollectNeighbors(vNeighbours, grid, vrt, NHT_EDGE_NEIGHBORS,
924 : // IsInSubset(shMarks, REM_CREASE));
925 :
926 : // // we have to choose a special normal
927 : // if(vNeighbours.size() != 2){
928 : // UG_LOG("n"<<vNeighbours.size());
929 : // continue;
930 : // }
931 : // else{
932 : // vector3 v1, v2;
933 : // VecSubtract(v1, v, aaPos[vNeighbours[0]]);
934 : // VecSubtract(v2, v, aaPos[vNeighbours[1]]);
935 : // VecNormalize(v1, v1);
936 : // VecNormalize(v2, v2);
937 : // VecAdd(n, v1, v2);
938 : // if(VecLengthSq(n) > SMALL)
939 : // VecNormalize(n, n);
940 : // else {
941 : // // both edges have the same direction.
942 : // // don't project normals
943 : // bProjectPointsToPlane = false;
944 : // }
945 : // }
946 : // }
947 : // else{
948 : // CollectNeighbors(vNeighbours, grid, vrt);
949 : // }
950 :
951 : // vNodes.resize(vNeighbours.size());
952 :
953 : // if(bProjectPointsToPlane){
954 : // for(size_t j = 0; j < vNodes.size(); ++j)
955 : // ProjectPointToPlane(vNodes[j], aaPos[vNeighbours[j]], v, n);
956 : // }
957 : // else{
958 : // for(size_t j = 0; j < vNodes.size(); ++j)
959 : // vNodes[j] = aaPos[vNeighbours[j]];
960 : // }
961 : // // perform point relocation
962 : // RelocatePointBySmoothing(aaPos[vrt], v, vNodes, 5, stepSize);
963 : // }
964 : // }
965 : // }
966 :
967 : /** Make sure that elements in gridOut directly correspond to
968 : * elements in gridIn*/
969 : template <class TGeomObj>
970 : void CopySelectionStatus(Selector& selOut, Grid& gridOut,
971 : Selector& selIn, Grid& gridIn)
972 : {
973 : typedef typename geometry_traits<TGeomObj>::iterator GeomObjIter;
974 : GeomObjIter endOut = gridOut.end<TGeomObj>();
975 : GeomObjIter endIn = gridIn.end<TGeomObj>();
976 : GeomObjIter iterOut = gridOut.begin<TGeomObj>();
977 : GeomObjIter iterIn = gridIn.begin<TGeomObj>();
978 :
979 : for(; (iterOut != endOut) && (iterIn != endIn); ++iterOut, ++iterIn)
980 : {
981 : if(selIn.is_selected(*iterIn))
982 : selOut.select(*iterOut);
983 : }
984 : }
985 :
986 : ////////////////////////////////////////////////////////////////////////
987 0 : bool AdjustEdgeLength(Grid& grid, SubsetHandler& shMarks,
988 : number minEdgeLen, number maxEdgeLen, int numIterations,
989 : bool projectPoints, bool adaptive)
990 : {
991 : PROFILE_FUNC();
992 :
993 : //TODO: replace this by a template parameter
994 : APosition aPos = aPosition;
995 :
996 : // we have to make sure that the mesh consist of triangles only,
997 : // since the algorithm would produce bad results if not.
998 0 : if(grid.num<Quadrilateral>() > 0)
999 : {
1000 : UG_LOG(" INFO: grid contains quadrilaterals. Converting to triangles...\n");
1001 :
1002 : //TODO: not gridIn but a copy of gridIn (pRefGrid) should be triangulated.
1003 0 : Triangulate(grid, grid.begin<Quadrilateral>(),
1004 0 : grid.end<Quadrilateral>());
1005 : }
1006 :
1007 : // make sure that grid and pRefGrid have position-attachments
1008 0 : if(!grid.has_vertex_attachment(aPos)){
1009 : UG_LOG(" vertex-position-attachment missing in AdjustEdgeLength. Aborting.\n");
1010 : return false;
1011 : }
1012 :
1013 : // make sure that faces create associated edges
1014 0 : if(!grid.option_is_enabled(FACEOPT_AUTOGENERATE_EDGES))
1015 : {
1016 : LOG(" INFO: auto-enabling FACEOPT_AUTOGENERATE_EDGES in AdjustEdgeLength.\n");
1017 0 : grid.enable_options(FACEOPT_AUTOGENERATE_EDGES);
1018 : }
1019 :
1020 : // position attachment
1021 : Grid::VertexAttachmentAccessor<APosition> aaPos(grid, aPos);
1022 :
1023 : // we need an integer attachment (a helper for ObtainSimpleGrid)
1024 : AInt aInt;
1025 0 : grid.attach_to<Vertex>(aInt);
1026 : Grid::AttachmentAccessor<Vertex, AInt> aaInt(grid, aInt);
1027 :
1028 : // TODO: normals shouldn't be created here but instead be passed to the method.
1029 : // attach the vertex normals.
1030 : ANormal aNorm;
1031 0 : grid.attach_to<Vertex>(aNorm);
1032 : Grid::AttachmentAccessor<Vertex, ANormal> aaNorm(grid, aNorm);
1033 0 : CalculateVertexNormals(grid, aPos, aNorm);
1034 :
1035 : // assign vertex marks
1036 0 : AssignFixedVertices(grid, shMarks);
1037 0 : AssignCreaseVertices(grid, shMarks);
1038 :
1039 : // we need an selector that holds all edges that are candidates for a operation
1040 0 : EdgeSelector esel(grid);
1041 0 : esel.enable_selection_inheritance(false);
1042 :
1043 : // sort the triangles of pRefGrid into a octree to speed-up projection performance
1044 : SPOctree octree;
1045 0 : node_tree::Traverser_ProjectPoint pojectionTraverser;
1046 0 : if(projectPoints){
1047 : //PROFILE_BEGIN(octree_construction);
1048 0 : octree = CreateOctree(grid, grid.begin<Triangle>(),
1049 0 : grid.end<Triangle>(),
1050 0 : 10, 30, false, aPos);
1051 :
1052 : //PROFILE_END();
1053 0 : if(!octree.valid()){
1054 : UG_LOG(" Octree creation failed in AdjustEdgeLength. Aborting.\n");
1055 : return false;
1056 : }
1057 : }
1058 :
1059 : // start the main iteration
1060 0 : for(int iteration = 0; iteration < numIterations; ++iteration)
1061 : {
1062 : // perform splits
1063 0 : esel.select(grid.begin<Edge>(), grid.end<Edge>());
1064 0 : if(!PerformSplits(grid, shMarks, esel, maxEdgeLen, aaPos, aaNorm, adaptive))
1065 : return false;
1066 :
1067 : // perform collapses
1068 0 : esel.select(grid.begin<Edge>(), grid.end<Edge>());
1069 0 : if(!PerformCollapses(grid, shMarks, esel, minEdgeLen, aaPos, aaNorm, aaInt, adaptive))
1070 : return false;
1071 :
1072 : // perform swaps
1073 0 : esel.select(grid.begin<Edge>(), grid.end<Edge>());
1074 0 : if(!PerformSwaps(grid, shMarks, esel, aaPos, aaNorm, aaInt))
1075 : return false;
1076 :
1077 : /*
1078 : // This is commented out, since it didn't help with the problems encountered
1079 : // in the geometries at that time.
1080 : // The algorithm however works and may prove useful in the future.
1081 : // fix bad triangles
1082 : // adjacent to crease-edges badly shaped triangles may occur,
1083 : // which have to be treated somehow.
1084 : esel.clear();
1085 : esel.select(shMarks.begin<Edge>(REM_CREASE), shMarks.end<Edge>(REM_CREASE));
1086 : esel.select(shMarks.begin<Edge>(REM_FIXED), shMarks.end<Edge>(REM_FIXED));
1087 : FixBadTriangles(grid, shMarks, esel, aaPos, aaNorm, 0.1);
1088 : */
1089 : // relocate points
1090 : /*LOG(" smoothing points...");
1091 : PerformSmoothing(grid, shMarks, aaPos, aaNorm, 10, 0.1);
1092 : LOG(" done\n");*/
1093 :
1094 : LOG(" updating normals...\n");
1095 0 : CalculateVertexNormals(grid, aPos, aNorm);
1096 :
1097 : // project points back on the surface
1098 0 : if(projectPoints)
1099 : {
1100 : LOG(" projecting points...");
1101 : //PROFILE_BEGIN(projecting_points);
1102 : for(VertexIterator iter = grid.vertices_begin();
1103 0 : iter != grid.vertices_end(); ++iter)
1104 : {
1105 : //TODO: project crease vertices onto creases only! Don't project fixed vertices
1106 0 : if(shMarks.get_subset_index(*iter) != REM_FIXED){
1107 : vector3 vNew;
1108 0 : if(pojectionTraverser.project(aaPos[*iter], octree/*, &aaNorm[*iter]*/)){
1109 : aaPos[*iter] = pojectionTraverser.get_closest_point();
1110 : }
1111 : else{
1112 : LOG("f");
1113 : }
1114 : }
1115 : }
1116 : //PROFILE_END();
1117 : LOG(" done\n");
1118 : }
1119 :
1120 0 : if(iteration < numIterations - 1){
1121 : LOG(" updating normals...");
1122 0 : CalculateVertexNormals(grid, aPos, aNorm);
1123 : }
1124 : }
1125 :
1126 :
1127 : // detach
1128 : grid.detach_from<Vertex>(aInt);
1129 : grid.detach_from<Vertex>(aNorm);
1130 :
1131 : GPLOTSAVE();
1132 : return true;
1133 0 : }
1134 :
1135 :
1136 :
1137 :
1138 : /*
1139 : ////////////////////////////////////////////////////////////////////////
1140 : // This is an alternative version for PerformSplits.
1141 : // It uses the Refine method.
1142 : // It is not yet optimized for maximum speed.
1143 : // While it performs a little less splits, overall runtime of
1144 : // AdjustEdgeLength is not better than with the original
1145 : // PerformSplits method.
1146 : template <class TAAPosVRT, class TAANormVRT>
1147 : bool PerformSplits(Grid& grid, SubsetHandler& shMarks, EdgeSelector& esel,
1148 : number maxEdgeLen, TAAPosVRT& aaPos, TAANormVRT& aaNorm)
1149 : {
1150 : AInt aInt;
1151 : grid.attach_to_edges(aInt);
1152 :
1153 : // compare squares
1154 : maxEdgeLen *= maxEdgeLen;
1155 :
1156 : Selector sel(grid);
1157 : sel.enable_autoselection(true);
1158 : sel.select(esel.begin<Edge>(), esel.end<Edge>());
1159 :
1160 : LOG(" performing splits\n");
1161 : int numSplits = 0;
1162 :
1163 : while(!sel.empty()){
1164 : // deselect all vertices and faces
1165 : sel.clear_selection<Vertex>();
1166 : sel.clear_selection<Face>();
1167 : sel.clear_selection<Volume>();
1168 :
1169 : // deselect all edges that shall not be splitted
1170 : EdgeIterator iter = sel.begin<Edge>();
1171 : while(iter != sel.end<Edge>()){
1172 : Edge* e = *iter;
1173 : ++iter;
1174 :
1175 : // the higher the curvature the smaller the maxEdgeLen.
1176 : // minimal lenFac is 0.1
1177 : number lenFac = CalculateLengthFac(grid, shMarks, e, aaPos);
1178 :
1179 : // fixed edges will not be refined
1180 : if(shMarks.get_subset_index(e) == REM_FIXED)
1181 : sel.deselect(e);
1182 : else if(VecDistanceSq(aaPos[e->vertex(0)], aaPos[e->vertex(1)]) < lenFac * maxEdgeLen)
1183 : sel.deselect(e);
1184 : }
1185 :
1186 : // refine the grid
1187 : Refine(grid, sel, aInt);
1188 :
1189 : // new vertices are selected
1190 : numSplits += sel.num<Vertex>();
1191 :
1192 : // re-triangulate
1193 : Triangulate(grid, &aaPos);
1194 :
1195 : // calculate normal for new vertices
1196 : //TODO: be careful with crease edges
1197 : for(VertexIterator iter = sel.begin<Vertex>();
1198 : iter != sel.end<Vertex>(); ++iter)
1199 : {
1200 : CalculateVertexNormal(aaNorm[*iter], grid, *iter, aaPos);
1201 : }
1202 :
1203 : sel.clear_selection<Vertex>();
1204 : sel.clear_selection<Face>();
1205 : sel.clear_selection<Volume>();
1206 : }
1207 :
1208 : grid.detach_from_edges(aInt);
1209 : LOG(" splits performed: " << numSplits << endl);
1210 : return true;
1211 : }
1212 : */
1213 :
1214 : }// end of namespace
|
