Branch data Line data Source code
1 : : #include "meshkit/AF2Front.hpp"
2 : :
3 : : // C++
4 : : #include <algorithm>
5 : : #include <cmath>
6 : :
7 : : // MeshKit
8 : : #include "meshkit/Error.hpp"
9 : : #include "meshkit/AF2LocalTransform.hpp"
10 : :
11 : 339619 : bool EndPointLess::operator()(const AF2Edge3D* const & oneEdge,
12 : : const AF2Edge3D* const & otherEdge) const
13 : : {
14 [ + + ]: 339619 : if (oneEdge->getStart()->getLocalId() < otherEdge->getStart()->getLocalId())
15 : : {
16 : 155597 : return true;
17 : : }
18 [ + + ]: 184022 : if (oneEdge->getStart()->getLocalId() > otherEdge->getStart()->getLocalId())
19 : : {
20 : 121174 : return false;
21 : : }
22 : : // otherwise the start endpoints are equal
23 [ + + ]: 62848 : if (oneEdge->getEnd()->getLocalId() < otherEdge->getEnd()->getLocalId())
24 : : {
25 : 18987 : return true;
26 : : }
27 : 43861 : return false;
28 : : }
29 : :
30 [ + - ][ + - ]: 92 : AF2Front::AF2Front() : points(), edges(), qualityCount(4, 0u)
31 : : {
32 : 46 : }
33 : :
34 : 103 : AF2Front::~AF2Front()
35 : : {
36 : : typedef std::set<AF2Edge3D*>::iterator EdgeSetItr;
37 [ + + ]: 142 : for (EdgeSetItr edgeItr = edges.begin(); edgeItr != edges.end(); ++edgeItr)
38 : : {
39 : 96 : delete *edgeItr;
40 : : }
41 [ - + ]: 57 : }
42 : :
43 [ # # ][ # # ]: 0 : AF2Front::AF2Front(const AF2Front & toCopy)
44 : : {
45 [ # # ]: 0 : MeshKit::Error notImpl(MeshKit::MK_NOT_IMPLEMENTED);
46 [ # # ]: 0 : notImpl.set_string("AF2Front copy construction is not supported.");
47 [ # # ]: 0 : throw notImpl;
48 : : }
49 : :
50 : 0 : AF2Front& AF2Front::operator=(const AF2Front & rhs)
51 : : {
52 [ # # ]: 0 : MeshKit::Error notImpl(MeshKit::MK_NOT_IMPLEMENTED);
53 [ # # ]: 0 : notImpl.set_string("AF2Front assignment operator is not supported.");
54 [ # # ]: 0 : throw notImpl;
55 : : }
56 : :
57 : 3544 : void AF2Front::qualityDecreased(const AF2Edge3D* const & anEdge)
58 : : {
59 : 3544 : unsigned int curQual = anEdge->getQualityLevel();
60 [ + + ]: 3544 : if (curQual > qualityCount.size())
61 : : {
62 [ + - ]: 54 : qualityCount.resize(curQual, 0u);
63 : : }
64 : 3544 : ++qualityCount[curQual - 1];
65 : 3544 : --qualityCount[curQual - 2];
66 : 3544 : }
67 : :
68 : 5366 : void AF2Front::addPoint(AF2Point3D* pointPtr)
69 : : {
70 [ + - ][ + - ]: 5366 : if (points.find(pointPtr) != points.end())
[ - + ]
71 : : {
72 [ # # ]: 0 : MeshKit::Error duplicate(MeshKit::MK_BAD_INPUT);
73 [ # # ]: 0 : duplicate.set_string("The point is already on the advancing front.");
74 [ # # ]: 0 : throw duplicate;
75 : : }
76 : 5366 : points[pointPtr] = 0;
77 : 5366 : }
78 : :
79 : 9317 : void AF2Front::advanceFront(std::list<AF2Edge3D*> edgeList)
80 : : {
81 : : typedef std::map<AF2Point3D*, int>::iterator PointCountItr;
82 : : typedef std::set<AF2Edge3D*>::iterator EdgeSetItr;
83 : : typedef std::list<AF2Edge3D*>::iterator EdgeListItr;
84 : :
85 : : // verify validity of endpoints, which could throw exception,
86 : : // before doing any other processing
87 [ + - ][ + - ]: 38445 : for (EdgeListItr itr = edgeList.begin(); itr != edgeList.end(); ++itr)
[ + + ]
88 : : {
89 [ + - ][ + - ]: 29128 : PointCountItr startPntItr = points.find((*itr)->getStart());
[ + - ]
90 [ + - ][ + - ]: 29128 : PointCountItr endPntItr = points.find((*itr)->getEnd());
[ + - ]
91 [ + - ][ + + ]: 29128 : if (startPntItr == points.end() || endPntItr == points.end())
[ + - ][ - + ]
[ + + ][ + - ]
[ + + # #
# # ]
92 : : {
93 [ + - ]: 1 : MeshKit::Error badArg(MeshKit::MK_BAD_INPUT);
94 : : badArg.set_string(
95 [ + - ]: 1 : "One or both endpoints of an edge are not on the advancing front.");
96 [ + - ]: 1 : throw badArg;
97 : : }
98 [ + - ][ + - ]: 29127 : if ((*itr)->getQualityLevel() != 1u)
[ + + ]
99 : : {
100 [ + - ]: 1 : MeshKit::Error badArg(MeshKit::MK_BAD_INPUT);
101 : : badArg.set_string(
102 [ + - ]: 1 : "Edges added to the advancing front must have quality level 1.");
103 [ + - ]: 2 : throw badArg;
104 : : }
105 : : }
106 : :
107 : : // decide which half edges need to be added
108 : : // and which (reverse) half edges need to be removed
109 [ + - ]: 9315 : std::list<AF2Edge3D*> edgesToAdd;
110 [ + - ]: 18630 : std::list<EdgeSetItr> edgesToRemove;
111 : :
112 [ + - ][ + - ]: 38437 : for (EdgeListItr itr = edgeList.begin(); itr != edgeList.end(); ++itr)
[ + + ]
113 : : {
114 [ + - ][ + - ]: 29122 : AF2Edge3D reverseEdge((*itr)->getEnd(), (*itr)->getStart());
[ + - ][ + - ]
[ + - ]
115 [ + - ]: 29122 : EdgeSetItr revEdgeItr = edges.find(&reverseEdge);
116 [ + - ][ + + ]: 29122 : if (revEdgeItr != edges.end())
117 : : {
118 : : // prepare for later processing of the reverse half edge that was
119 : : // allocated on the heatp to remove it from the edge set and delete it
120 : : // don't delete it now, since that would make it more difficult to
121 : : // determine whether its endpoints should be deleted
122 [ + - ]: 14513 : edgesToRemove.push_back(revEdgeItr);
123 : : // right now delete the forward half edge that was allocated on the heap
124 [ + - ]: 14513 : delete *itr;
125 : : }
126 : : else
127 : : {
128 : : // prepare for later processing of the forward half edge
129 : : // don't add it now, since that would prevent adding both half edges
130 : : // of an "isolated" edge
131 [ + - ][ + - ]: 14609 : edgesToAdd.push_back(*itr);
132 : : }
133 : : }
134 : :
135 : : // add the half edges that need to be added
136 [ + - ][ + - ]: 23924 : for (EdgeListItr itr = edgesToAdd.begin(); itr != edgesToAdd.end(); ++itr)
[ + + ]
137 : : {
138 : : // add it to the set of edges
139 [ + - ][ + - ]: 14609 : edges.insert(*itr);
140 : : // start tracking the edge's quality level
141 [ + - ]: 14609 : ++qualityCount[0];
142 [ + - ][ + - ]: 14609 : (*itr)->setObserver(this);
143 : : // increase the count of the number of edges incident to each endpoint
144 [ + - ][ + - ]: 14609 : ++points[(*itr)->getStart()];
[ + - ]
145 [ + - ][ + - ]: 14609 : ++points[(*itr)->getEnd()];
[ + - ]
146 : : // update the distance to the boundary (locally)
147 : : // * Note * The distance to the boundary may be globally inaccurate
148 : : // after this update, but this agrees with how NetGen did the update.
149 : : // It would be possible to cascade the distance to boundary updates
150 : : // and ensure globally accuracy, but that might affect a lot of points.
151 : 14609 : unsigned int maxDist = 1 + std::min(
152 [ + - ][ + - ]: 14609 : (*itr)->getStart()->getDistanceToBoundary(),
[ + - ]
153 [ + - ][ + - ]: 29218 : (*itr)->getEnd()->getDistanceToBoundary());
[ + - ][ + - ]
154 [ + - ][ + - ]: 14609 : (*itr)->getStart()->limitDistanceToBoundary(maxDist);
[ + - ]
155 [ + - ][ + - ]: 14609 : (*itr)->getEnd()->limitDistanceToBoundary(maxDist);
[ + - ]
156 : : }
157 : :
158 : : // remove the half edges that need to be removed and any endpoints of
159 : : // those edges that have no incident edges after the edge is removed
160 [ + - + - ]: 47656 : for (std::list<EdgeSetItr>::iterator itr = edgesToRemove.begin();
[ + + ]
161 : 23828 : itr != edgesToRemove.end(); ++itr)
162 : : {
163 [ + - ]: 14513 : EdgeSetItr edgeItr = *itr;
164 : : // decrease the count of the number of edges incident to each endpoint
165 [ + - ][ + - ]: 14513 : --points[(*edgeItr)->getStart()];
[ + - ]
166 [ + - ][ + - ]: 14513 : --points[(*edgeItr)->getEnd()];
[ + - ]
167 : : // remove endpoints from points if they now have zero incident edges
168 [ + - ][ + - ]: 14513 : if (points[(*edgeItr)->getStart()] == 0u)
[ + - ][ + + ]
169 : : {
170 [ + - ][ + - ]: 1341 : points.erase((*edgeItr)->getStart());
[ + - ]
171 : : }
172 [ + - ][ + - ]: 14513 : if (points[(*edgeItr)->getEnd()] == 0u)
[ + - ][ + + ]
173 : : {
174 [ + - ][ + - ]: 3921 : points.erase((*edgeItr)->getEnd());
[ + - ]
175 : : }
176 : : // stop tracking the quality level of the edge
177 [ + - ][ + - ]: 14513 : --qualityCount[(*edgeItr)->getQualityLevel() - 1u];
[ + - ]
178 : : // delete the edge object that was allocated on the heap
179 [ + - ]: 14513 : delete *edgeItr;
180 : : // remove the (now dangling) pointer to the edge from the edge set
181 [ + - ]: 14513 : edges.erase(edgeItr);
182 : 9315 : }
183 : 9315 : }
184 : :
185 : 10055 : bool AF2Front::isEmpty() const
186 : : {
187 [ + + ][ + + ]: 10055 : return edges.empty() && points.empty();
188 : : }
189 : :
190 : 9990 : unsigned int AF2Front::getMaximumQuality() const
191 : : {
192 [ + - ]: 30134 : for (unsigned int qualityLevel = 1;
193 : 15067 : qualityLevel <= qualityCount.size(); ++qualityLevel)
194 : : {
195 [ + + ]: 15067 : if (qualityCount[qualityLevel - 1] > 0)
196 : : {
197 : 9990 : return qualityLevel;
198 : : }
199 : : }
200 : :
201 : 0 : return 0u;
202 : : }
203 : :
204 : 9990 : AF2Neighborhood* AF2Front::selectNeighborhood(
205 : : const AF2LocalTransformMaker* const & transformMaker) const
206 : : {
207 : : typedef std::set<AF2Edge3D*>::const_iterator EdgeSetItr;
208 : : typedef std::map<AF2Point3D*, int>::const_iterator PointCountItr;
209 : :
210 : : // Select a baseline edge that minimizes a metric value for the metric
211 : : // edge quality level + sum of endpoint distances to boundary.
212 : : // The current implementation loops through all of the edges
213 : : // to find an edge with minimal metric value. This could be made faster
214 : : // using a C++ priority queue, updating each time the metric value changes.
215 : 9990 : AF2Edge3D* baselineEdge = NULL;
216 : : unsigned int minMetricVal;
217 [ + - ][ + - ]: 544749 : for (EdgeSetItr edgeItr = edges.begin(); edgeItr != edges.end(); ++edgeItr)
[ + + ]
218 : : {
219 [ + + ]: 534759 : if (baselineEdge == NULL)
220 : : {
221 [ + - ][ + - ]: 9988 : minMetricVal = (*edgeItr)->getQualityLevel() +
222 [ + - ][ + - ]: 9988 : (*edgeItr)->getStart()->getDistanceToBoundary() +
[ + - ]
223 [ + - ][ + - ]: 9988 : (*edgeItr)->getEnd()->getDistanceToBoundary();
[ + - ]
224 [ + - ]: 9988 : baselineEdge = *edgeItr;
225 : : }
226 : : else
227 : : {
228 [ + - ][ + - ]: 524771 : unsigned int metricVal = (*edgeItr)->getQualityLevel() +
229 [ + - ][ + - ]: 524771 : (*edgeItr)->getStart()->getDistanceToBoundary() +
[ + - ]
230 [ + - ][ + - ]: 524771 : (*edgeItr)->getEnd()->getDistanceToBoundary();
[ + - ]
231 [ + + ]: 524771 : if (metricVal < minMetricVal)
232 : : {
233 : 192436 : metricVal = minMetricVal;
234 [ + - ]: 192436 : baselineEdge = *edgeItr;
235 : : }
236 : : }
237 : : }
238 : :
239 : : // throw an exception if no baseline edge was found
240 [ + + ]: 9990 : if (baselineEdge == NULL)
241 : : {
242 [ + - ]: 2 : MeshKit::Error fail(MeshKit::MK_FAILURE);
243 [ + - ]: 2 : fail.set_string("The advancing front has no edges.");
244 [ + - ]: 2 : throw fail;
245 : : }
246 : :
247 : : // determine the neighborhood size from the length and quality
248 : : // level of the baseline edge
249 [ + - ]: 9988 : AF2Point3D* baseStart = baselineEdge->getStart();
250 [ + - ]: 9988 : AF2Point3D* baseEnd = baselineEdge->getEnd();
251 [ + - ][ + - ]: 9988 : double xDiff = baseEnd->getX() - baseStart->getX();
252 [ + - ][ + - ]: 9988 : double yDiff = baseEnd->getY() - baseStart->getY();
253 [ + - ][ + - ]: 9988 : double zDiff = baseEnd->getZ() - baseStart->getZ();
254 : 9988 : double sqLen = xDiff * xDiff + yDiff * yDiff + zDiff * zDiff;
255 [ + - ][ + - ]: 9988 : unsigned int qLevelFactor = std::min(baselineEdge->getQualityLevel(), 7u);
256 : 9988 : double ngbhdSize = 4.0 * (3u + qLevelFactor) * (3u + qLevelFactor) * sqLen;
257 : :
258 : : // create lists to hold the points and (other) edges in the neighborhood
259 [ + - ]: 9988 : std::list<AF2Point3D*> ngbhdPoints;
260 [ + - ]: 19976 : std::list<const AF2Edge3D*> ngbhdEdges;
261 : : // create a set to track points included in the neighorhood and avoid
262 : : // duplicating them if they are endpoints of multiple neighborhood edges
263 [ + - ]: 19976 : std::set<AF2Point3D*> ngbhdPointSet;
264 : :
265 : : // find edges that are close enough to the baseline edge
266 : : // to be included in the neighborhood
267 : 9988 : AF2Point3D* edgeStart = NULL;
268 : 9988 : AF2Point3D* edgeEnd = NULL;
269 : 9988 : double edgeXDiff = 0.0;
270 : 9988 : double edgeYDiff = 0.0;
271 : 9988 : double edgeZDiff = 0.0;
272 : 9988 : double edgeSqLen = 0.0;
273 : 9988 : double dotProd = 0.0;
274 [ + - ][ + - ]: 544747 : for (EdgeSetItr edgeItr = edges.begin(); edgeItr != edges.end(); ++edgeItr)
[ + + ]
275 : : {
276 : : // skip the baseline edge to avoid duplicating it in the neighborhood
277 [ + - ][ + + ]: 534759 : if ((*edgeItr) == baselineEdge)
278 : : {
279 : 9988 : continue;
280 : : }
281 : :
282 : : // calculate components of vector along the direction of the edge
283 : : // and the squared length of the edge
284 [ + - ][ + - ]: 524771 : edgeStart = (*edgeItr)->getStart();
285 [ + - ][ + - ]: 524771 : edgeEnd = (*edgeItr)->getEnd();
286 [ + - ][ + - ]: 524771 : edgeXDiff = edgeEnd->getX() - edgeStart->getX();
287 [ + - ][ + - ]: 524771 : edgeYDiff = edgeEnd->getY() - edgeStart->getY();
288 [ + - ][ + - ]: 524771 : edgeZDiff = edgeEnd->getZ() - edgeStart->getZ();
289 : 1049542 : edgeSqLen = edgeXDiff * edgeXDiff +
290 : 1049542 : edgeYDiff * edgeYDiff + edgeZDiff * edgeZDiff;
291 : :
292 : : // calculate components of vector from edge start point to base start point
293 [ + - ][ + - ]: 524771 : xDiff = baseStart->getX() - edgeStart->getX();
294 [ + - ][ + - ]: 524771 : yDiff = baseStart->getY() - edgeStart->getY();
295 [ + - ][ + - ]: 524771 : zDiff = baseStart->getZ() - edgeStart->getZ();
296 : :
297 : : // calculate the dot product of the two vectors
298 : 524771 : dotProd = edgeXDiff * xDiff + edgeYDiff * yDiff + edgeZDiff * zDiff;
299 : :
300 : : // calculate the components of the vector from the baseline edge
301 : : // start point to the closest point on the edge, placing the results
302 : : // in (and reusing the variables) xDiff, yDiff, and zDiff
303 [ + + ]: 524771 : if (dotProd < 0)
304 : : {
305 : : // the orthogonal projection of the baseline edge start point onto
306 : : // the line that contains the current edge lies outside the edge
307 : : // and the edge start point is closest to the baseline edge start point
308 : :
309 : : // the vector from the edge start point to the baseline edge start
310 : : // point is correct, and xDiff, yDiff, and zDiff already has that vector
311 : : // --> do nothing
312 : : }
313 [ + + ]: 283082 : else if (dotProd > edgeSqLen)
314 : : {
315 : : // the orthogonal projection of the baseline edge start point onto
316 : : // the line that contains the current edge lies outside the edge
317 : : // and the edge end point is closest to the baseline edge start point
318 : :
319 : : // the vector from the edge end point to the baseline edge start
320 : : // point is the one that is needed
321 [ + - ][ + - ]: 241592 : xDiff = baseStart->getX() - edgeEnd->getX();
322 [ + - ][ + - ]: 241592 : yDiff = baseStart->getY() - edgeEnd->getY();
323 [ + - ][ + - ]: 241592 : zDiff = baseStart->getZ() - edgeEnd->getZ();
324 : : }
325 : : else
326 : : {
327 : : // the orthogonal projection of the baseline edge start point onto
328 : : // the line that contains the current edge lies on the edge
329 : : // and is the closest point on the edge to the baseline edge start point
330 : :
331 : : // calculate the components of the projection of the vector from the
332 : : // edge start point to the baseline edge start point onto the vector
333 : : // in the direction of the edge, reusing edge difference variables
334 : 41490 : edgeXDiff = dotProd * edgeXDiff / edgeSqLen;
335 : 41490 : edgeYDiff = dotProd * edgeYDiff / edgeSqLen;
336 : 41490 : edgeZDiff = dotProd * edgeZDiff / edgeSqLen;
337 : :
338 : : // calculate the components of the projection of the vector from the
339 : : // edge start point to the baseline edge start point onto the plane
340 : : // orthogonal to the direction of the edge
341 : 41490 : xDiff = xDiff - edgeXDiff;
342 : 41490 : yDiff = yDiff - edgeYDiff;
343 : 41490 : zDiff = zDiff - edgeZDiff;
344 : : }
345 : :
346 : : // calculate the square of the length of said vector, which is also the
347 : : // square of the distance from the edge to the baseline edge start point
348 : 524771 : sqLen = xDiff * xDiff + yDiff * yDiff + zDiff * zDiff;
349 : :
350 : : // if the edge is close enough, add it and its endpoints
351 : : // to the neighborhood
352 [ + + ]: 524771 : if (sqLen < ngbhdSize)
353 : : {
354 [ + - ][ + - ]: 270293 : ngbhdEdges.push_back(*edgeItr);
355 [ + - ][ + - ]: 270293 : if (ngbhdPointSet.find(edgeStart) == ngbhdPointSet.end())
[ + + ]
356 : : {
357 [ + - ]: 137463 : ngbhdPointSet.insert(edgeStart);
358 [ + - ]: 137463 : ngbhdPoints.push_back(edgeStart);
359 : : }
360 [ + - ][ + - ]: 270293 : if (ngbhdPointSet.find(edgeEnd) == ngbhdPointSet.end())
[ + + ]
361 : : {
362 [ + - ]: 150857 : ngbhdPointSet.insert(edgeEnd);
363 [ + - ]: 150857 : ngbhdPoints.push_back(edgeEnd);
364 : : }
365 : : }
366 : : }
367 : :
368 : : // find points that are close enough to the baseline edge
369 : : // to be included in the neighborhood
370 [ + - ][ + - ]: 541562 : for (PointCountItr itr = points.begin(); itr != points.end(); ++itr)
[ + + ]
371 : : {
372 [ + - ]: 531574 : AF2Point3D* frontPoint = itr->first;
373 : :
374 : : // calculate the square of the distance from the baseline edge start point
375 [ + - ][ + - ]: 531574 : xDiff = frontPoint->getX() - baseStart->getX();
376 [ + - ][ + - ]: 531574 : yDiff = frontPoint->getY() - baseStart->getY();
377 [ + - ][ + - ]: 531574 : zDiff = frontPoint->getZ() - baseStart->getZ();
378 : 531574 : sqLen = xDiff * xDiff + yDiff * yDiff + zDiff * zDiff;
379 : :
380 : : // if the distance is small enough and the point has not been added
381 : : // to the neighborhood yet, add it to the neighborhood
382 [ + + ][ + + ]: 1596190 : if (sqLen < ngbhdSize &&
[ + + ]
383 [ + - ][ + - ]: 1064616 : ngbhdPointSet.find(frontPoint) == ngbhdPointSet.end())
[ + + ][ + + ]
[ # # # # ]
384 : : {
385 [ + - ]: 2 : ngbhdPointSet.insert(frontPoint);
386 [ + - ]: 2 : ngbhdPoints.push_back(frontPoint);
387 : : }
388 : : }
389 : :
390 : : // use the transform maker to make a local transform appropriate
391 : : // to the neighborhood
392 : : AF2LocalTransform* localTransform =
393 [ + - ]: 9988 : transformMaker->makeLocalTransform(ngbhdPoints, baselineEdge, ngbhdEdges);
394 : :
395 : : // build and return the neighborhood object
396 : : return new AF2Neighborhood(
397 [ + - ][ + - ]: 19976 : ngbhdPoints, baselineEdge, ngbhdEdges, localTransform);
398 : : }
|
