GCC Code Coverage Report

Directory: ./
File: Geometry/Tools/locator.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 99 148 66.9%
Branches: 76 196 38.8%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors:
13 //
14
15 // System includes
16 #include <numeric>
17
18 // External includes
19
20 // Project includes
21 #include "Geometry/Tools/locator.hpp"
22 #include "Geometry/geo_utilities.hpp"
23
24 namespace felisce
25 {
26
27 17 Locator::Locator(GeometricMeshRegion* mesh, const bool allowProjection, const double tolRescaling, const int verbose):
28 17 m_allowProjection(allowProjection), m_verbose(verbose)
29 {
30 // Set tolerance
31 17 m_tol = GeoUtilities::GEO_TOL * tolRescaling;
32
33 // Set mesh pointer
34 17 m_mesh = mesh;
35
36 // Allocate and Initialize array m_markElement
37
2/2
✓ Branch 1 taken 425 times.
✓ Branch 2 taken 17 times.
 442 for (std::size_t i_eltyp = 0; i_eltyp < m_markElement.size(); ++i_eltyp)
38 425 m_markElement[i_eltyp] = nullptr;
39
40 ElementType eltType;
41
2/2
✓ Branch 2 taken 17 times.
✓ Branch 3 taken 17 times.
 34 for (std::size_t i = 0; i < m_mesh->bagElementTypeDomain().size(); ++i) {
42 17 eltType = m_mesh->bagElementTypeDomain()[i];
43
44
2/4
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 17 times.
✗ Branch 6 not taken.
 17 m_markElement[eltType] = new std::size_t[m_mesh->numElements(eltType)];
45
2/2
✓ Branch 1 taken 7884 times.
✓ Branch 2 taken 17 times.
 7901 for (felInt k = 0; k < m_mesh->numElements(eltType); ++k)
46 7884 m_markElement[eltType][k] = 0;
47 }
48 17 m_mark = 0;
49
50 // Build edges or faces
51
2/2
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 7 times.
 17 if ( m_mesh->domainDim() == 2 ) {
52
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 if ( !m_mesh->statusEdges() )
53
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 m_mesh->buildEdges();
54
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 } else if ( m_mesh->domainDim() == 3 ) {
55
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 if ( !m_mesh->statusFaces() )
56
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 m_mesh->buildFaces();
57 }
58
59 // Build bounding box
60
1/2
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
 17 m_bbox = felisce::make_shared<BoundingBox>(*m_mesh);
61
1/2
✓ Branch 3 taken 17 times.
✗ Branch 4 not taken.
 17 m_bbox->computeDelta(mesh->numCoor());
62
63 // If default hashtab size < numPoints reduce it
64
2/2
✓ Branch 1 taken 15 times.
✓ Branch 2 taken 2 times.
 17 if ( static_cast<std::size_t>(m_mesh->numPoints()) < HashTable::HASHSIZE )
65
1/2
✓ Branch 3 taken 15 times.
✗ Branch 4 not taken.
 15 m_hash = felisce::make_shared<HashTable>(m_mesh, static_cast<std::size_t>(m_mesh->numPoints()), m_bbox.get());
66 else
67
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 m_hash = felisce::make_shared<HashTable>(m_mesh, HashTable::HASHSIZE, m_bbox.get());
68 17 }
69
70 /***********************************************************************************/
71 /***********************************************************************************/
72
73 17 Locator::~Locator()
74 {
75
2/2
✓ Branch 1 taken 425 times.
✓ Branch 2 taken 17 times.
 442 for (std::size_t i_eltyp = 0; i_eltyp < m_markElement.size(); ++i_eltyp) {
76
2/2
✓ Branch 1 taken 17 times.
✓ Branch 2 taken 408 times.
 425 if( m_markElement[i_eltyp] != nullptr ) {
77
78
1/2
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
 17 delete [] m_markElement[i_eltyp];
79 17 m_markElement[i_eltyp] = nullptr;
80 }
81 }
82 17 }
83
84 /***********************************************************************************/
85 /***********************************************************************************/
86
87 13 void Locator::findListPointsInMesh(std::vector<felInt>& indElt, std::vector<ElementType>& typElt, const std::vector<Point>& lstPoints)
88 {
89 13 std::vector<seedElement> elt;
90
91
1/2
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
 13 findListPointsInMesh(elt, lstPoints);
92
93
2/2
✓ Branch 1 taken 41 times.
✓ Branch 2 taken 13 times.
 54 for (std::size_t i = 0; i < lstPoints.size(); ++i) {
94 41 typElt[i] = elt[i].first;
95 41 indElt[i] = elt[i].second;
96 }
97 13 }
98
99 /***********************************************************************************/
100 /***********************************************************************************/
101
102 13 void Locator::findListPointsInMesh(std::vector<seedElement>& elt, const std::vector<Point>& lstPoints)
103 {
104
1/2
✓ Branch 3 taken 13 times.
✗ Branch 4 not taken.
 13 elt.assign(lstPoints.size(), seedElement(GeometricMeshRegion::ELEMENT_TYPE_COUNTER, -1) );
105
106 // This method requires the hashtable
107
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 13 times.
 13 if ( !m_hash ) FEL_ERROR("Localization: hashtable not allocated");
108
109 // For each point of the structure mesh
110
2/2
✓ Branch 1 taken 41 times.
✓ Branch 2 taken 13 times.
 54 for (std::size_t k = 0; k < lstPoints.size(); ++k) {
111
112
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 41 times.
 41 if ( m_verbose > 3 )
113 ✗ std::cout << "Localization: Processing point " << k << " = " << lstPoints[k] <<std::endl;
114
115 // Point has not been localised yet
116
1/2
✓ Branch 1 taken 41 times.
✗ Branch 2 not taken.
 41 if( elt[k].second == -1 ) {
117
118 // Locate the point and even project it if necessary - when it is inside the box but not inside the mesh
119
2/2
✓ Branch 3 taken 10 times.
✓ Branch 4 taken 31 times.
 41 if ( !(findPointInMesh(elt[k], lstPoints[k])) ) {
120 // Failure
121
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 if ( m_allowProjection )
122 ✗ FEL_ERROR( "Localization: point " + std::to_string(k) + " has not been localized even with the projection on the mesh!" );
123
124
4/8
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 10 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 10 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 10 times.
✗ Branch 12 not taken.
 10 std::cout << "Localization: point " + std::to_string(k) + " has not been localized!" << std::endl;
125 }
126 }
127 }
128
129
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13 times.
 13 if ( m_verbose > 3 )
130 ✗ std::cout << "Localization: all points in the list have been analysed!" << std::endl;
131 13 }
132
133 /***********************************************************************************/
134 /***********************************************************************************/
135
136 1805 bool Locator::findPointInMesh(seedElement& elt, const Point& point)
137 {
138 ElementType eltType;
139 1805 seedElement tmpSeed;
140 Point pointInRefElm;
141 1805 std::size_t mark = ++getMark();
142
143 // If the point is in the bounding box of the mesh
144
2/2
✓ Branch 3 taken 1788 times.
✓ Branch 4 taken 17 times.
 1805 if( m_bbox->isInBBox(m_mesh->numCoor(), point, m_tol) ) {
145
146 // Try to locate the point in an optimized way with the hashtable
147
1/2
✓ Branch 3 taken 1788 times.
✗ Branch 4 not taken.
 1788 tmpSeed = m_hash->findSeed(point);
148
1/2
✓ Branch 1 taken 1788 times.
✗ Branch 2 not taken.
 1788 tmpSeed = localizePoint(tmpSeed, point, pointInRefElm, mark);
149
150 // If succeeded, return True
151
2/4
✓ Branch 0 taken 1788 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1788 times.
✗ Branch 3 not taken.
 1788 if ( tmpSeed.first != GeometricMeshRegion::ELEMENT_TYPE_COUNTER && tmpSeed.second > -1 ) {
152
153
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1788 times.
 1788 if ( m_verbose > 3 )
154 ✗ std::cout << "Localization: success, point located using the hash table!" << std::endl;
155
156 // Copy information from the hashtable
157 1788 elt = tmpSeed;
158
159 1788 return true;
160 }
161
162 // If failure,
163 ✗ if ( m_verbose > 3 )
164 ✗ std::cout << "Localization: failure, point cannot be located using the hash table!" << std::endl;
165
166 // Go through all type of elements of mesh
167 ✗ mark = ++getMark();
168 ✗ for (std::size_t l = 0; l < m_mesh->bagElementTypeDomain().size(); ++l) {
169 ✗ eltType = m_mesh->bagElementTypeDomain()[l];
170
171 ✗ const int nbrPnt = GeometricMeshRegion::m_numPointsPerElt[eltType];
172 ✗ std::vector<felInt> veridxOfElt( nbrPnt, 0);
173 ✗ std::vector<Point*> pointsOfElt( nbrPnt, nullptr);
174
175 // Go through all elements of this type
176 ✗ for (felInt jel = 0; jel < m_mesh->numElements(eltType); ++jel) {
177 ✗ tmpSeed.first = eltType;
178 ✗ tmpSeed.second = jel;
179
180 // Skip already visited elements
181 ✗ if ( markElement(eltType, jel) == mark )
182 ✗ continue;
183
184 // Check bounding box
185 ✗ m_mesh->getOneElement(eltType, jel, veridxOfElt, pointsOfElt, 0);
186 ✗ BoundingBox box(m_mesh->numCoor(), pointsOfElt);
187
188 ✗ if ( !box.isInBBox(m_mesh->numCoor(), point, m_tol) )
189 ✗ continue;
190
191 // Check if point inside current element
192 ✗ if ( m_mesh->findLocalCoord(tmpSeed.first, tmpSeed.second, point, pointInRefElm) ) {
193
194 ✗ if ( m_verbose > 3 )
195 ✗ std::cout << "Localization: success, the point has been located in one of the element of the mesh!" << std::endl;
196
197 // Copy information from the hashtable
198 ✗ elt = tmpSeed;
199
200 ✗ return true;
201 }
202 }
203
204 // If failure
205 ✗ if ( m_verbose > 3 )
206 ✗ std::cout << "Localization: failure, point cannot be found into the mesh, projection needed!" << std::endl;
207 }
208 }
209
210 // If the point is outside the bounding box or it needs to be projected
211
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 10 times.
 17 if ( m_allowProjection ) {
212
213 // If the projection is successful
214
2/4
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 7 times.
✗ Branch 4 not taken.
 7 if( projectPointOnMesh(elt, point) ) {
215
216
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
 7 if ( m_verbose > 3 )
217 ✗ std::cout << "Localization: success, the projection has been located in one of the element of the mesh!" << std::endl;
218
219 7 return true;
220 }
221 }
222
223 // If failure
224
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 if ( m_verbose > 3 )
225 ✗ std::cout << "Localization: failure, the point (as well as its projection if enable) has not been localised into the mesh!" << std::endl;
226
227 10 return false;
228 }
229
230 /***********************************************************************************/
231 /***********************************************************************************/
232
233 3179 Locator::seedElement Locator::localizePoint(seedElement seed, const Point& point, Point& pointInRefElm, std::size_t mark)
234 {
235
236 // Localize point in mesh
237 while (1) {
238
239 // Boundary reached
240
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3179 times.
 3179 if ( seed.second < 0 ) {
241 ✗ seed.first = GeometricMeshRegion::ELEMENT_TYPE_COUNTER;
242 ✗ seed.second = -1;
243
244 ✗ return seed;
245 }
246
247 // Element already marked -> looping
248
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 3179 times.
 3179 if ( markElement(seed.first, seed.second) == mark ) {
249 ✗ seed.first = GeometricMeshRegion::ELEMENT_TYPE_COUNTER;
250 ✗ seed.second = -1;
251
252 ✗ return seed;
253 }
254
255 // Mark current element
256 3179 markElement(seed.first, seed.second) = mark;
257
258 // If point inside element, seed fuond. Otherwise go to next element
259
2/2
✓ Branch 1 taken 1788 times.
✓ Branch 2 taken 1391 times.
 3179 if ( m_mesh->findLocalCoord(seed.first, seed.second, point, pointInRefElm) )
260 1788 return seed;
261 }
262
263 FEL_ERROR("Localization: this point should never be reached seed found, or stuck on boundary or cycle");
264
265 return seed;
266 }
267
268 /***********************************************************************************/
269 /***********************************************************************************/
270
271 7 bool Locator::projectPointOnMesh(seedElement& seed, const Point& pntIn)
272 {
273 bool isProjectionInMesh, isProjectionFound;
274 int numPointPerElt;
275 double distProjPlane;
276 double distToElt, minDistToElt;
277 Point projection_tmp, projection;
278
279 7 std::vector<felInt> veridxOfElt;
280 7 std::vector<Point*> pointsOfElt;
281
282 ElementType eltTypeBD;
283
284 7 minDistToElt = std::numeric_limits<double>::max();
285
286 // Go through all the boundary elements of the mesh
287 7 isProjectionFound = false;
288 7 const std::vector<ElementType>& bagElementTypeBoundary = m_mesh->bagElementTypeDomainBoundary();
289
2/2
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 7 times.
 14 for (std::size_t i = 0; i < bagElementTypeBoundary.size(); ++i) {
290 7 eltTypeBD = bagElementTypeBoundary[i];
291
292 // Resize vectors
293 7 numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltTypeBD];
294
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 veridxOfElt.assign(numPointPerElt, 0);
295
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 pointsOfElt.assign(numPointPerElt, nullptr);
296
297
2/2
✓ Branch 1 taken 52 times.
✓ Branch 2 taken 7 times.
 59 for (felInt ii = 0; ii < m_mesh->numElements(eltTypeBD); ++ii) {
298
299 // Get vertices of the current element
300
1/2
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
 52 m_mesh->getOneElement(eltTypeBD, ii, veridxOfElt, pointsOfElt, 0);
301
302 // We compute the element normal
303
1/2
✓ Branch 2 taken 52 times.
✗ Branch 3 not taken.
 52 GeoUtilities::projectOnGeomEntity( static_cast<felInt>(m_mesh->domainDim()), pointsOfElt, pntIn, projection_tmp);
304
305 // Compute the distance from the point to its projection on the plane
306 52 distProjPlane = projection_tmp.dist(pntIn);
307
308 // If the new projected point is closer than the previous ones
309
2/2
✓ Branch 0 taken 31 times.
✓ Branch 1 taken 21 times.
 52 if ( distProjPlane < minDistToElt ) {
310
311 // Get closest projection
312
1/2
✓ Branch 1 taken 31 times.
✗ Branch 2 not taken.
 31 GeoUtilities::getClosestProjectionOnElement(eltTypeBD, pointsOfElt, projection_tmp, projection_tmp, m_tol);
313
314 // Compute total distance from the point and its projection in the current considered element
315 31 distToElt = projection_tmp.dist(pntIn);
316
317 // Looking for the minimal total distance
318
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 21 times.
 31 if ( distToElt < minDistToElt ) {
319
320 // Update minimal distance
321 10 minDistToElt = distToElt;
322
323 // Save coordinates of the point which minimizes the total distance
324 10 projection = projection_tmp;
325
326 // Projection found
327 10 isProjectionFound = true;
328 }
329 }
330 }
331 }
332
333
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
 7 if ( !isProjectionFound ) {
334
335 ✗ if ( m_verbose > 3 )
336 ✗ std::cout << "Localization: failure, no valid projection could be found!" << std::endl;
337
338 ✗ seed.first = GeometricMeshRegion::ELEMENT_TYPE_COUNTER;
339 ✗ seed.second = -1;
340
341 ✗ return false;
342 }
343
344
345
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
 7 if ( m_verbose > 3 )
346 ✗ std::cout << "Localization: Projected point is " << projection << std::endl;
347
348 // Turn off projection to avoid infinity loop
349 7 m_allowProjection = false;
350
351 // Try to localize the projected point
352
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 isProjectionInMesh = findPointInMesh(seed, projection);
353
354 // Turn on projection
355 7 m_allowProjection = true;
356
357 7 return isProjectionInMesh;
358 7 }
359
360 /***********************************************************************************/
361 /***********************************************************************************/
362
363 ✗ void Locator::resetMarkElement()
364 {
365 // Reset m_markElement
366 ElementType eltType;
367 ✗ for (std::size_t i = 0; i < m_mesh->bagElementTypeDomain().size(); ++i) {
368 ✗ eltType = m_mesh->bagElementTypeDomain()[i];
369
370 ✗ if ( m_markElement[eltType] ) {
371
372 // Reset arrays
373 ✗ for (felInt k = 0; k < m_mesh->numElements(eltType); ++k)
374 ✗ m_markElement[eltType][k] = 0;
375 }
376 }
377
378 ✗ m_mark = 0;
379 }
380
381 }
382