GCC Code Coverage Report

Directory: ./
File: Geometry/Tools/intersector.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 251 429 58.5%
Branches: 177 803 22.0%
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
17 // External includes
18
19 // Project includes
20 #include "Geometry/geometricEdges.hpp"
21 #include "Geometry/geometricFaces.hpp"
22 #include "Geometry/neighFacesOfEdges.hpp"
23 #include "Geometry/neighVolumesOfFaces.hpp"
24 #include "Geometry/Tools/intersector.hpp"
25
26 namespace felisce
27 {
28
29 4 Intersector::Intersector(GeometricMeshRegion* volMesh, GeometricMeshRegion* strMesh, const double tolRescaling, const int verbose):
30 4 m_volMesh(volMesh), m_strMesh(strMesh), m_verbose(verbose)
31 {
32 // Set tolerance
33 4 m_tol = GeoUtilities::GEO_TOL * tolRescaling;
34
35
4/14
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
 4 FEL_CHECK( m_volMesh->domainDim() == 2 || m_volMesh->domainDim() == 3,"Intersector: cannot be called for a mesh with domainDim = " + std::to_string(m_volMesh->domainDim()) + " different from 2 or 3!\n" );
36
5/8
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 3 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 4 times.
 4 FEL_CHECK( m_volMesh->bagElementTypeDomain().size() == 1 &&
37 (m_volMesh->bagElementTypeDomain()[0] == GeometricMeshRegion::Tria3 ||
38 m_volMesh->bagElementTypeDomain()[0] == GeometricMeshRegion::Tetra4 ), "Intersector: volume mesh must contain only Tria3 or Tetra4!\n");
39
5/8
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 3 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 4 times.
 4 FEL_CHECK( m_strMesh->bagElementTypeDomain().size() == 1 &&
40 (m_strMesh->bagElementTypeDomain()[0] == GeometricMeshRegion::Seg2 ||
41 m_strMesh->bagElementTypeDomain()[0] == GeometricMeshRegion::Tria3 ), "Intersector: interface mesh must contain only Seg2 or Tria3!\n");
42 4 }
43
44 /***********************************************************************************/
45 /***********************************************************************************/
46
47 52 void Intersector::intersectMeshes(const std::vector<felInt>& listVerToElt, std::vector<felInt>& listOfIntElt) const
48 {
49
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 52 times.
 52 FEL_ASSERT ( m_strMesh->numPoints() == static_cast<int>(listVerToElt.size()) );
50
51
2/2
✓ Branch 1 taken 51 times.
✓ Branch 2 taken 1 times.
 52 if ( m_volMesh->domainDim() == 2 )
52 51 findIntersectedElements2D(listVerToElt, listOfIntElt);
53 else /*if ( m_volMesh->domainDim() == 3 )*/
54 1 findIntersectedElements3D(listVerToElt, listOfIntElt);
55 52 }
56
57 /***********************************************************************************/
58 /***********************************************************************************/
59
60 51 void Intersector::findIntersectedElements2D(const std::vector<felInt>& listVerToElt, std::vector<felInt>& listOfIntElt) const
61 {
62
1/2
✓ Branch 2 taken 51 times.
✗ Branch 3 not taken.
 51 std::vector<Point*> pntOfTria(3, nullptr);
63
1/2
✓ Branch 2 taken 51 times.
✗ Branch 3 not taken.
 51 std::vector<felInt> verOfTria(3, 0);
64
1/2
✓ Branch 2 taken 51 times.
✗ Branch 3 not taken.
 51 std::vector<Point> pntOfEdge(2);
65
1/2
✓ Branch 2 taken 51 times.
✗ Branch 3 not taken.
 51 std::vector<felInt> verOfEdge(2, -1);
66
1/2
✓ Branch 2 taken 51 times.
✗ Branch 3 not taken.
 51 std::vector<Point> intPoint(2);
67 51 std::array<Location,2> intLocation;
68 std::array<double, 3> bary;
69 bool isInside0, isInside1;
70 51 felInt idxOfPreElt, idxOfNxtElt, idxOfElt = -1;
71
72 // Set of intersected elements
73 51 std::set<felInt> setOfIntElt;
74
75 // Build the edges of the volume mesh
76
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 51 times.
 51 if ( !m_volMesh->statusEdges() )
77 ✗ m_volMesh->buildEdges();
78
79 // For each edge of the structure mesh
80
3/4
✓ Branch 1 taken 1304 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1253 times.
✓ Branch 4 taken 51 times.
 1304 for (felInt i = 0; i < m_strMesh->getNumElement1D(); ++i) {
81
82 // Get nodes of the current structure edge
83
1/2
✓ Branch 1 taken 1253 times.
✗ Branch 2 not taken.
 1253 m_strMesh->getOneElement(GeometricMeshRegion::Seg2, i, verOfEdge, 0);
84 1253 pntOfEdge[0] = m_strMesh->listPoints()[verOfEdge[0]];
85 1253 pntOfEdge[1] = m_strMesh->listPoints()[verOfEdge[1]];
86
87
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1253 times.
 1253 if ( m_verbose > 3 )
88 ✗ std::cout << "Intersector: Processing edge " << i << " = " << verOfEdge[0] << " " << verOfEdge[1] << std::endl;
89
90 // Find volume element to start the intersection
91
1/2
✓ Branch 2 taken 1253 times.
✗ Branch 3 not taken.
 1253 if ( listVerToElt[ verOfEdge[0] ] != -1 ) {
92 1253 idxOfElt = listVerToElt[ verOfEdge[0] ];
93 }
94 ✗ else if ( listVerToElt[ verOfEdge[1] ] != -1 ) {
95 ✗ idxOfElt = listVerToElt[ verOfEdge[1] ];
96
97 ✗ std::swap(verOfEdge[0], verOfEdge[1]);
98 ✗ std::swap(pntOfEdge[0], pntOfEdge[1]);
99 }
100 else {
101 ✗ idxOfElt = getGermElement2D(verOfEdge);
102 }
103
104
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1253 times.
 1253 if ( m_verbose > 3 )
105 ✗ std::cout << "Interface edge " << i << " = " << verOfEdge[0] << " " << verOfEdge[1] << ", vertex " << verOfEdge[0] << " is in volume element " << idxOfElt << std::endl;
106
107 // Reset previous intersected element to -1
108 1253 idxOfPreElt = -1;
109
110 // Look for intersection until pntOfEdge[1] is inside idxOfElt
111 do {
112
113 // Get nodes and coordinates of the current triangle
114
1/2
✓ Branch 1 taken 2304 times.
✗ Branch 2 not taken.
 2304 m_volMesh->getOneElement(GeometricMeshRegion::Tria3, idxOfElt, verOfTria, pntOfTria, 0);
115
116
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2304 times.
 2304 if ( m_verbose > 9 )
117 ✗ std::cout << "Volume element " << idxOfElt << " = " << verOfTria[0] << " " << verOfTria[1] << " " << verOfTria[2] << "." << std::endl;
118
119 // Check the position of pntOfEdge[0] wrt to the current element
120
1/2
✓ Branch 3 taken 2304 times.
✗ Branch 4 not taken.
 2304 isInside0 = GeoUtilities::isPointInTri(pntOfTria, pntOfEdge[0], bary, intLocation[0], m_tol);
121
122
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2304 times.
 2304 if ( !isInside0 ){
123 ✗ FEL_ERROR("Intersector::findIntersectedElements2D : pntOfEdge[0] must be inside idxOfElt !")
124 }
125
126
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2304 times.
 2304 if ( m_verbose > 9 )
127 ✗ std::cout << "Point0 " << pntOfEdge[0] << " is in element " << idxOfElt << " with case = " << static_cast<std::underlying_type_t<Position>>(intLocation[0].first) << " and position = " << static_cast<std::underlying_type_t<Position>>(intLocation[0].second) << "." << std::endl;
128
129 // Add the intersected elements in the list
130
1/2
✓ Branch 2 taken 2304 times.
✗ Branch 3 not taken.
 2304 addIntersetedEltToList2D(idxOfElt, intLocation[0], setOfIntElt);
131
132 // Check the position of pntOfEdge[0] wrt to the current element
133
1/2
✓ Branch 3 taken 2304 times.
✗ Branch 4 not taken.
 2304 isInside1 = GeoUtilities::isPointInTri( pntOfTria, pntOfEdge[1], bary, intLocation[1], m_tol);
134
135 // If pntOfEdge[1] is inside the current element, add the intersected elements to the list
136
2/2
✓ Branch 0 taken 202 times.
✓ Branch 1 taken 2102 times.
 2304 if ( isInside1 ) {
137
138
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 202 times.
 202 if ( m_verbose > 9 )
139 ✗ std::cout << "Point1 " << pntOfEdge[1] << " is in element " << idxOfElt << " with case = " << static_cast<std::underlying_type_t<Position>>(intLocation[1].first) << " and position = " << static_cast<std::underlying_type_t<Position>>(intLocation[1].second) << "." << std::endl;
140
141 // Add the intersected elements in the list
142
1/2
✓ Branch 2 taken 202 times.
✗ Branch 3 not taken.
 202 addIntersetedEltToList2D(idxOfElt, intLocation[1], setOfIntElt);
143 }
144 // If pntOfEdge[1] is outside the current element, we must detect the intersections
145 else {
146
147
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2102 times.
 2102 if ( m_verbose > 9 )
148 ✗ std::cout << "Point1 " << pntOfEdge[1] << " is not in the volume element " << idxOfElt << " so there is an intersection to detect." << std::endl;
149
150 2102 bool intFnd = false;
151 // If verOfEdge[0] on a vertex
152
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2102 times.
 2102 if ( intLocation[0].first == Position::OnVertex )
153 ✗ intFnd = computeIntersectionVertex2D(idxOfElt, verOfTria[intLocation[0].second], idxOfPreElt, idxOfNxtElt, pntOfEdge, intPoint, isInside1);
154 // If verOfEdge[0] on an edge
155
2/2
✓ Branch 1 taken 1052 times.
✓ Branch 2 taken 1050 times.
 2102 else if ( intLocation[0].first == Position::OnEdge ) {
156
1/2
✓ Branch 2 taken 1052 times.
✗ Branch 3 not taken.
 1052 intFnd = computeIntersectionEdge2D(idxOfElt, intLocation[0].second, idxOfPreElt, idxOfNxtElt, pntOfEdge, intPoint, isInside1);
157 }
158 else { // if ( intLocation[0].first == Position::Inside ) {
159
1/2
✓ Branch 1 taken 1050 times.
✗ Branch 2 not taken.
 1050 intFnd = computeIntersectionInside2D(idxOfElt, idxOfNxtElt, pntOfEdge, intPoint);
160 }
161
162
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2102 times.
 2102 if ( !intFnd )
163 ✗ FEL_ERROR("Intersector::findIntersectedElements2D : there must be at least one intersection point !");
164
165
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2102 times.
 2102 if ( m_verbose > 9 )
166 ✗ std::cout << "Intersection detected " << intPoint[0] << "." << std::endl;
167
168 // Update pntOfEdge
169 2102 pntOfEdge[0] = intPoint[0]; // TODO is it possible to have two intPoints???
170
171 // Update previous intersected element
172 2102 idxOfPreElt = idxOfElt;
173
174 // Update triangle for new intersections
175 2102 idxOfElt = idxOfNxtElt;
176 }
177
2/2
✓ Branch 0 taken 1051 times.
✓ Branch 1 taken 1253 times.
 2304 } while ( !isInside1 );
178 }
179
180 // Fill the vector of intersected elements
181
1/2
✓ Branch 3 taken 51 times.
✗ Branch 4 not taken.
 51 listOfIntElt.assign(setOfIntElt.begin(), setOfIntElt.end());
182
183
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 51 times.
 51 if ( m_verbose > 9 )
184 ✗ std::cout << "Intersector: all the edges of the structure have been considered!" << std::endl;
185 51 }
186
187 /***********************************************************************************/
188 /***********************************************************************************/
189
190 2506 void Intersector::addIntersetedEltToList2D(const felInt idxOfElt, const Location& location, std::set<felInt>& setOfIntElt) const
191 {
192 felInt eltIdx;
193
1/2
✓ Branch 2 taken 2506 times.
✗ Branch 3 not taken.
 2506 std::vector<felInt> verOfTria(3, 0);
194
195 // Get nodes and coordinates of the current triangle
196
1/2
✓ Branch 1 taken 2506 times.
✗ Branch 2 not taken.
 2506 m_volMesh->getOneElement(GeometricMeshRegion::Tria3, idxOfElt, verOfTria, 0);
197
198
3/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1154 times.
✓ Branch 2 taken 1350 times.
✗ Branch 3 not taken.
 2506 switch ( location.first ) {
199
200 // Vertex is on a vertex of the triangle
201 2 case Position::OnVertex : {
202
203 // Variables
204 GeometricMeshRegion::ElementType eltTyp;
205 2 GeometricMeshRegion::BallMap ball;
206 felInt verPosInBall;
207
208 // Get ball
209
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 m_volMesh->getVertexBall(verOfTria[location.second], GeometricMeshRegion::Tria3, idxOfElt, ball);
210
211 // Loop over ball
212
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
 4 for (std::size_t iel = 0; iel < ball.size(); ++iel) {
213 2 std::tie(eltTyp, eltIdx, verPosInBall) = ball[iel];
214
215 // Add element to the set of intersected elements
216
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 setOfIntElt.insert(eltIdx);
217 }
218 2 break;
219 2 }
220 // Vertex is on an edge of the triangle
221 1154 case Position::OnEdge : {
222
223 // Variables
224
1/2
✓ Branch 2 taken 1154 times.
✗ Branch 3 not taken.
 1154 std::vector<felInt> verEdgOfTria(2,-1);
225 1154 Edge edgOfElt;
226
227 1154 auto& mapEltEdgToVer = GeometricMeshRegion::m_eltTypMapFacToVer[GeometricMeshRegion::Tria3];
228
229 // Get the ied-th edge of the triangle
230 1154 verEdgOfTria[0] = verOfTria[ mapEltEdgToVer[location.second][0] ];
231 1154 verEdgOfTria[1] = verOfTria[ mapEltEdgToVer[location.second][1] ];
232
1/2
✓ Branch 1 taken 1154 times.
✗ Branch 2 not taken.
 1154 m_volMesh->getEdge(verEdgOfTria, edgOfElt);
233
234 // Loop over elements sharing the edge
235 1154 auto& listNeigh = edgOfElt.listNeighFacesOfEdges();
236
2/2
✓ Branch 1 taken 2205 times.
✓ Branch 2 taken 1154 times.
 3359 for (std::size_t iel = 0; iel < listNeigh.size(); ++iel) {
237 2205 eltIdx = listNeigh[iel]->idFace();
238
239 // Add element to the set of intersected elements
240
1/2
✓ Branch 1 taken 2205 times.
✗ Branch 2 not taken.
 2205 setOfIntElt.insert(eltIdx);
241 }
242 1154 break;
243 1154 }
244 // Vertex is inside the triangle
245 1350 case Position::Inside : {
246
247 // Add element to the set of intersected elements
248
1/2
✓ Branch 1 taken 1350 times.
✗ Branch 2 not taken.
 1350 setOfIntElt.insert(idxOfElt);
249 1350 break;
250 }
251
252 ✗ default : {
253 ✗ FEL_ERROR("Intersector::addIntersetedEltToList2D : the vertex must be on a vertex, edge or inside the triangle !");
254 }
255 }
256 2506 }
257
258 /***********************************************************************************/
259 /***********************************************************************************/
260
261 ✗ bool Intersector::computeIntersectionVertex2D(const felInt idxOfElt, const felInt idxOfVer, const felInt idxOfPreElt, felInt& idxNxtElt, const std::vector<Point>& pntOfEdge, std::vector<Point>& intPoint, bool& isInside) const
262 {
263 GeometricMeshRegion::ElementType eltTyp;
264 felInt eltIdx;
265 felInt verPos;
266 unsigned short int ied;
267 int nbrIntPnt;
268 ✗ std::vector<Point*> pntOfElt(3, nullptr);
269 ✗ std::vector<felInt> verOfElt(3, 0);
270 std::array<double,3> bary;
271 std::array<Point*,2> pntEdgOfElt;
272
273 ✗ GeometricMeshRegion::BallMap ball;
274
275 ✗ auto& mapEltEdgToVer = GeometricMeshRegion::m_eltTypMapFacToVer[GeometricMeshRegion::Tria3];
276
277 // Get ball
278 ✗ m_volMesh->getVertexBall(idxOfVer, GeometricMeshRegion::Tria3, idxOfElt, ball);
279
280 // Loop over ball
281 ✗ for (std::size_t iel = 0; iel < ball.size(); ++iel) {
282 ✗ std::tie(eltTyp, eltIdx, verPos) = ball[iel];
283
284 ✗ if ( eltIdx == idxOfPreElt )
285 ✗ continue;
286
287 // Get nodes and coordinates of the current triangle
288 ✗ m_volMesh->getOneElement(GeometricMeshRegion::Tria3, eltIdx, verOfElt, pntOfElt, 0);
289
290 // Check if pntOfEdge[1] is inside the current element
291 // If it is not inside, look for an intersection with the opposite edge
292 ✗ if ( !GeoUtilities::isPointInTri(pntOfElt, pntOfEdge[1], bary, m_tol) ) {
293
294 // Get the points of the edge opposite to the vertex
295 ✗ ied = mapEltEdgToVer[verPos][1];
296 ✗ pntEdgOfElt[0] = pntOfElt[ mapEltEdgToVer[ied][0] ];
297 ✗ pntEdgOfElt[1] = pntOfElt[ mapEltEdgToVer[ied][1] ];
298
299 // Compute intersection with opposite edge
300 ✗ nbrIntPnt = GeoUtilities::computeSegmentsIntersection(pntOfEdge[0], pntOfEdge[1], *pntEdgOfElt[0], *pntEdgOfElt[1], intPoint[0], intPoint[1], m_tol);
301
302 // If intersection
303 ✗ if ( nbrIntPnt > 0 ) {
304
305 ✗ if ( m_verbose > 9 )
306 ✗ std::cout << "computeIntersectionVertex2D: Intersection detected in " << intPoint[0] << "." << std::endl;
307
308 // Return next element to intersect
309 ✗ idxNxtElt = eltIdx;
310
311 ✗ return true;
312 }
313 }
314 // If it is inside, nothing to do
315 else {
316
317 ✗ if ( m_verbose > 9 )
318 ✗ std::cout << "computeIntersectionVertex2D: Point " << pntOfEdge[1] << " has been found in element " << eltIdx << "." << std::endl;
319
320 // Set intersection
321 ✗ intPoint[0] = pntOfEdge[1];
322
323 // Return next element to intersect
324 ✗ idxNxtElt = eltIdx;
325
326 // Point inside element
327 ✗ isInside = true;
328
329 ✗ return true;
330 }
331 }
332
333 ✗ return false;
334 }
335
336 /***********************************************************************************/
337 /***********************************************************************************/
338
339 1052 bool Intersector::computeIntersectionEdge2D(const felInt idxOfElt, const int idxEdgInt, const felInt idxOfPreElt, felInt& idxNxtElt, const std::vector<Point>& pntOfEdge, std::vector<Point>& intPoint, bool& isInside) const
340 {
341 felInt eltIdx;
342 int nbrIntPnt;
343
1/2
✓ Branch 2 taken 1052 times.
✗ Branch 3 not taken.
 1052 std::vector<Point*> pntOfNxtElt(3, nullptr);
344
1/2
✓ Branch 2 taken 1052 times.
✗ Branch 3 not taken.
 1052 std::vector<felInt> verOfNxtElt(3, 0);
345
1/2
✓ Branch 2 taken 1052 times.
✗ Branch 3 not taken.
 1052 std::vector<felInt> verOfElt(3, 0);
346
1/2
✓ Branch 2 taken 1052 times.
✗ Branch 3 not taken.
 1052 std::vector<felInt> verEdgOfElt(2,-1);
347 std::array<double,3> bary;
348 std::array<felInt,2> verEdgOfNxtElt;
349 std::array<Point*,2> pntEdgOfNxtElt;
350 1052 Edge edgOfElt;
351
352 1052 auto& mapEltEdgToVer = GeometricMeshRegion::m_eltTypMapFacToVer[GeometricMeshRegion::Tria3];
353
354 // Get nodes and coordinates of idxOfElt
355
1/2
✓ Branch 1 taken 1052 times.
✗ Branch 2 not taken.
 1052 m_volMesh->getOneElement(GeometricMeshRegion::Tria3, idxOfElt, verOfElt, 0);
356
357 // Get the ied-th edge of the triangle
358 1052 verEdgOfElt[0] = verOfElt[ mapEltEdgToVer[idxEdgInt][0] ];
359 1052 verEdgOfElt[1] = verOfElt[ mapEltEdgToVer[idxEdgInt][1] ];
360
1/2
✓ Branch 1 taken 1052 times.
✗ Branch 2 not taken.
 1052 m_volMesh->getEdge(verEdgOfElt, edgOfElt);
361
362 // Loop over elements sharing the edge
363 1052 auto& listNeigh = edgOfElt.listNeighFacesOfEdges();
364
1/2
✓ Branch 1 taken 2103 times.
✗ Branch 2 not taken.
 2103 for (std::size_t iel = 0; iel < listNeigh.size(); ++iel) {
365 2103 eltIdx = listNeigh[iel]->idFace();
366
367
2/2
✓ Branch 0 taken 1051 times.
✓ Branch 1 taken 1052 times.
 2103 if ( eltIdx == idxOfPreElt )
368 1051 continue;
369
370 // Get nodes and coordinates of the current triangle
371
1/2
✓ Branch 1 taken 1052 times.
✗ Branch 2 not taken.
 1052 m_volMesh->getOneElement(GeometricMeshRegion::Tria3, eltIdx, verOfNxtElt, pntOfNxtElt, 0);
372
373 // Check if pntOfEdge[1] is inside the current element
374 // If it is not inside, look for an intersection with the opposite edges
375
3/4
✓ Branch 2 taken 1052 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 1051 times.
 1052 if ( !GeoUtilities::isPointInTri(pntOfNxtElt, pntOfEdge[1], bary, m_tol) ) {
376
377 // Loop over the two opposite edges
378
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 for (std::size_t ied = 0; ied < 3; ++ied) {
379
380 // Get vertices of the ied-th edge
381 2 verEdgOfNxtElt[0] = verOfNxtElt[ mapEltEdgToVer[ied][0] ];
382 2 verEdgOfNxtElt[1] = verOfNxtElt[ mapEltEdgToVer[ied][1] ];
383
384 // Skip edge
385
3/4
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 1 times.
 2 if ( std::is_permutation(verEdgOfElt.begin(), verEdgOfElt.end(), verEdgOfNxtElt.begin()) )
386 1 continue;
387
388 // Get the edge
389 1 pntEdgOfNxtElt[0] = pntOfNxtElt[ mapEltEdgToVer[ied][0] ];
390 1 pntEdgOfNxtElt[1] = pntOfNxtElt[ mapEltEdgToVer[ied][1] ];
391
392 // Compute intersection with the edge
393
1/2
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
 1 nbrIntPnt = GeoUtilities::computeSegmentsIntersection(pntOfEdge[0], pntOfEdge[1], *pntEdgOfNxtElt[0], *pntEdgOfNxtElt[1], intPoint[0], intPoint[1], m_tol);
394
395 // If intersection
396
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if ( nbrIntPnt > 0 ) {
397
398
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if ( m_verbose > 9 )
399 ✗ std::cout << "computeIntersectionEdge2D: Intersection detected in " << intPoint[0] << "." << std::endl;
400
401 // Return next element to intersect
402 1 idxNxtElt = eltIdx;
403
404 1 return true;
405 }
406 }
407 }
408 // If it is inside, nothing to do
409 else {
410
411
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1051 times.
 1051 if ( m_verbose > 9 )
412 ✗ std::cout << "computeIntersectionEdge2D: Point " << pntOfEdge[1] << " has been found in element " << eltIdx << "." << std::endl;
413
414 // Set intersection
415 1051 intPoint[0] = pntOfEdge[1];
416
417 // Return next element to intersect
418 1051 idxNxtElt = eltIdx;
419
420 // Point inside element
421 1051 isInside = true;
422
423 1051 return true;
424 }
425 }
426
427 ✗ return false;
428 1052 }
429
430 /***********************************************************************************/
431 /***********************************************************************************/
432
433 1050 bool Intersector::computeIntersectionInside2D(const felInt idxOfElt, felInt& idxNxtElt, const std::vector<Point>& pntOfEdge, std::vector<Point>& intPoint) const
434 {
435 int nbrIntPnt;
436
1/2
✓ Branch 2 taken 1050 times.
✗ Branch 3 not taken.
 1050 std::vector<Point*> pntOfElt(3, nullptr);
437
1/2
✓ Branch 2 taken 1050 times.
✗ Branch 3 not taken.
 1050 std::vector<felInt> verOfElt(3, 0);
438 std::array<Point*,2> pntEdgOfElt;
439
440 1050 auto& mapEltEdgToVer = GeometricMeshRegion::m_eltTypMapFacToVer[GeometricMeshRegion::Tria3];
441
442 // Get nodes and coordinates of the current triangle
443
1/2
✓ Branch 1 taken 1050 times.
✗ Branch 2 not taken.
 1050 m_volMesh->getOneElement(GeometricMeshRegion::Tria3, idxOfElt, verOfElt, pntOfElt, 0);
444
445 // Loop over the edges
446
1/2
✓ Branch 0 taken 2100 times.
✗ Branch 1 not taken.
 2100 for (std::size_t ied = 0; ied < 3; ++ied) {
447
448 // Get the edge
449 2100 pntEdgOfElt[0] = pntOfElt[ mapEltEdgToVer[ied][0] ];
450 2100 pntEdgOfElt[1] = pntOfElt[ mapEltEdgToVer[ied][1] ];
451
452 // Compute intersection with the edge
453
1/2
✓ Branch 7 taken 2100 times.
✗ Branch 8 not taken.
 2100 nbrIntPnt = GeoUtilities::computeSegmentsIntersection(pntOfEdge[0], pntOfEdge[1], *pntEdgOfElt[0], *pntEdgOfElt[1], intPoint[0], intPoint[1], m_tol);
454
455 // If intersection
456
2/2
✓ Branch 0 taken 1050 times.
✓ Branch 1 taken 1050 times.
 2100 if ( nbrIntPnt > 0 ) {
457
458 // Return next element to intersect
459 1050 idxNxtElt = idxOfElt;
460
461 1050 return true;
462 }
463 }
464
465 ✗ return false;
466 1050 }
467
468 /***********************************************************************************/
469 /***********************************************************************************/
470
471 ✗ felInt Intersector::getGermElement2D(const std::vector<felInt>& verOfEdge) const
472 {
473 ✗ std::vector<Point*> pntOfTria(3, nullptr);
474 ✗ std::vector<felInt> verOfTria(3, 0);
475 std::array<double, 3> bary;
476
477 // Go through all elements of the mesh
478 ✗ for (felInt k = 0; k < m_volMesh->numElements(GeometricMeshRegion::Tria3); ++k) {
479
480 // Get nodes of the current triangle
481 ✗ m_volMesh->getOneElement(GeometricMeshRegion::Tria3, k, verOfTria, pntOfTria, 0);
482
483 ✗ if( GeoUtilities::isPointInTri( pntOfTria, m_strMesh->listPoints()[verOfEdge[0]], bary, m_tol) )
484 ✗ return k;
485 }
486
487 ✗ FEL_ERROR("Intersector::getGermElement2D cannot localize the point into the volume mesh.");
488 ✗ return -1;
489 }
490
491 /***********************************************************************************/
492 /***********************************************************************************/
493
494 1 void Intersector::findIntersectedElements3D(const std::vector<felInt>& listVerToElt, std::vector<felInt>& listOfIntElt) const
495 {
496
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 std::vector<Point*> pntOfTetra(4, nullptr);
497
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 std::vector<felInt> verOfTetra(4, 0);
498
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 std::vector<Point> pntOfEdge(2);
499
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 std::vector<felInt> verOfEdge(2, -1);
500
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 std::vector<Point> intPoint(2);
501 1 std::array<Location,2> intLocation;
502 std::array<double, 4> bary;
503 bool isInside0, isInside1;
504 1 felInt idxOfPreElt, idxOfNxtElt, idxOfElt = -1;
505
506 // Set of intersected elements
507 1 std::set<felInt> setOfIntElt;
508
509 // Build the edges of the interface mesh
510
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 if ( !m_strMesh->statusEdges() )
511
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_strMesh->buildEdges();
512
513 // Build the faces of the volume mesh
514
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 1 if ( !m_volMesh->statusFaces() )
515 ✗ m_volMesh->buildFaces();
516
517 // For each edge of the structure mesh
518
2/2
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 1 times.
 6 for (felInt i = 0; i < m_strMesh->numEdges(); ++i) {
519
520 // Get nodes of the current structure edge
521 5 verOfEdge[0] = m_strMesh->listEdges()[i]->idBeg();
522 5 verOfEdge[1] = m_strMesh->listEdges()[i]->idEnd();
523 5 pntOfEdge[0] = m_strMesh->listPoints()[verOfEdge[0]];
524 5 pntOfEdge[1] = m_strMesh->listPoints()[verOfEdge[1]];
525
526
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if ( m_verbose > 3 )
527 ✗ std::cout << "Intersector: Processing edge " << i << " = " << verOfEdge[0] << " " << verOfEdge[1] << std::endl;
528
529 // Find volume element to start the intersection
530
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 if ( listVerToElt[ verOfEdge[0] ] != -1 ) {
531 5 idxOfElt = listVerToElt[ verOfEdge[0] ];
532 }
533 ✗ else if ( listVerToElt[ verOfEdge[1] ] != -1 ) {
534 ✗ idxOfElt = listVerToElt[ verOfEdge[1] ];
535
536 ✗ std::swap(verOfEdge[0], verOfEdge[1]);
537 ✗ std::swap(pntOfEdge[0], pntOfEdge[1]);
538 }
539 else {
540 ✗ idxOfElt = getGermElement3D(verOfEdge);
541 }
542
543
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if ( m_verbose > 3 )
544 ✗ std::cout << "Interface edge " << i << " = " << verOfEdge[0] << " " << verOfEdge[1] << ", vertex " << verOfEdge[0] << " is in volume element " << idxOfElt << std::endl;
545
546 // Reset previous intersected element to -1
547 5 idxOfPreElt = -1;
548
549 // Look for intersection until pntOfEdge[1] is inside idxOfElt
550 do {
551
552 // Get nodes and coordinates of the current tetrahedron
553
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 m_volMesh->getOneElement(GeometricMeshRegion::Tetra4, idxOfElt, verOfTetra, pntOfTetra, 0);
554
555
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 if ( m_verbose > 9 )
556 ✗ std::cout << "Volume element " << idxOfElt << " = " << verOfTetra[0] << " " << verOfTetra[1] << " " << verOfTetra[2] << " " << verOfTetra[3] << "." << std::endl;
557
558 // Check the position of pntOfEdge[0] wrt to the current element
559
1/2
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
 10 isInside0 = GeoUtilities::isPointInTetra(pntOfTetra, pntOfEdge[0], bary, intLocation[0], m_tol);
560
561
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 if ( !isInside0 ){
562 ✗ FEL_ERROR("Intersector::findIntersectedElements3D : pntOfEdge[0] must be inside idxOfElt !")
563 }
564
565
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 if ( m_verbose > 9 )
566 ✗ std::cout << "Point0 " << pntOfEdge[0] << " is in element " << idxOfElt << " with case = " << static_cast<std::underlying_type_t<Position>>(intLocation[0].first) << " and position = " << static_cast<std::underlying_type_t<Position>>(intLocation[0].second) << "." << std::endl;
567
568 // Add the intersected elements in the list
569
1/2
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
 10 addIntersetedEltToList3D(idxOfElt, intLocation[0], setOfIntElt);
570
571 // Check the position of pntOfEdge[0] wrt to the current element
572
1/2
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
 10 isInside1 = GeoUtilities::isPointInTetra(pntOfTetra, pntOfEdge[1], bary, intLocation[1], m_tol);
573
574 // If pntOfEdge[1] is inside the current element, add the intersected elements to the list
575
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 if ( isInside1 ) {
576
577 ✗ if ( m_verbose > 9 )
578 ✗ std::cout << "Point1 " << pntOfEdge[1] << " is in element " << idxOfElt << " with case = " << static_cast<std::underlying_type_t<Position>>(intLocation[1].first) << " and position = " << static_cast<std::underlying_type_t<Position>>(intLocation[1].second) << "." << std::endl;
579
580 // Add the intersected elements in the list
581 ✗ addIntersetedEltToList3D(idxOfElt, intLocation[1], setOfIntElt);
582 }
583 // If pntOfEdge[1] is outside the current element, we must detect the intersections
584 else {
585
586
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 if ( m_verbose > 9 )
587 ✗ std::cout << "Point1 " << pntOfEdge[1] << " is not in the volume element " << idxOfElt << " so there is an intersection to detect." << std::endl;
588
589 10 bool intFnd = false;
590 // If verOfEdge[0] on a vertex
591
2/2
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
 10 if ( intLocation[0].first == Position::OnVertex ){
592
1/2
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
 5 intFnd = computeIntersectionVertex3D(idxOfElt, verOfTetra[intLocation[0].second], idxOfPreElt, idxOfNxtElt, pntOfEdge, intPoint, isInside1);
593 }
594 // If verOfEdge[0] on an edge
595
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 else if ( intLocation[0].first == Position::OnEdge ) {
596 5 auto& edge = GeometricMeshRegion::m_eltTypMapEdgToVer[GeometricMeshRegion::Tetra4][intLocation[0].second];
597
2/4
✓ Branch 6 taken 5 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 5 times.
✗ Branch 10 not taken.
 5 intFnd = computeIntersectionEdge3D(idxOfElt, {verOfTetra[edge[0]], verOfTetra[edge[1]]}, idxOfPreElt, idxOfNxtElt, pntOfEdge, intPoint, isInside1);
598 }
599 ✗ else if ( intLocation[0].first == Position::OnFace ) {
600 ✗ intFnd = computeIntersectionFace3D(idxOfElt, intLocation[0].second, idxOfPreElt, idxOfNxtElt, pntOfEdge, intPoint, isInside1);
601 }
602 else { // if ( intLocation[0].first == Position::Inside ) {
603 ✗ intFnd = computeIntersectionInside3D(idxOfElt, idxOfNxtElt, pntOfEdge, intPoint);
604 }
605
606
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 if ( !intFnd )
607 ✗ FEL_ERROR("Intersector::findIntersectedElements3D : there must be at least one intersection point !");
608
609
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 if ( m_verbose > 9 )
610 ✗ std::cout << "Intersection detected " << intPoint[0] << "." << std::endl;
611
612 // Update pntOfEdge
613 10 pntOfEdge[0] = intPoint[0]; // TODO is it possible to have two intPoints???
614
615 // Update previous intersected element
616 10 idxOfPreElt = idxOfElt;
617
618 // Update element for new intersections
619 10 idxOfElt = idxOfNxtElt;
620 }
621
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 10 } while ( !isInside1 );
622 }
623
624 // Fill the vector of intersected elements
625
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 listOfIntElt.assign(setOfIntElt.begin(), setOfIntElt.end());
626
627
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if ( m_verbose > 9 )
628 ✗ std::cout << "Intersector: all the edges of the structure have been considered!" << std::endl;
629 1 }
630
631 /***********************************************************************************/
632 /***********************************************************************************/
633
634 10 void Intersector::addIntersetedEltToList3D(const felInt idxOfElt, const Location& location, std::set<felInt>& setOfIntElt) const
635 {
636 felInt eltIdx;
637
1/2
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
 10 std::vector<felInt> verOfTetra(4, 0);
638
639 // Get nodes and coordinates of the current element
640
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 m_volMesh->getOneElement(GeometricMeshRegion::Tetra4, idxOfElt, verOfTetra, 0);
641
642
2/5
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 10 switch ( location.first ) {
643
644 // Vertex is on a vertex of the element
645 5 case Position::OnVertex : {
646
647 // Variables
648 GeometricMeshRegion::ElementType eltTyp;
649 5 GeometricMeshRegion::BallMap ball;
650 felInt verPosInBall;
651
652 // Get ball
653
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 m_volMesh->getVertexBall(verOfTetra[location.second], GeometricMeshRegion::Tetra4, idxOfElt, ball);
654
655 // Loop over ball
656
2/2
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
 10 for (std::size_t iel = 0; iel < ball.size(); ++iel) {
657 5 std::tie(eltTyp, eltIdx, verPosInBall) = ball[iel];
658
659 // Add element to the set of intersected elements
660
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 setOfIntElt.insert(eltIdx);
661 }
662 5 break;
663 5 }
664 // Vertex is on an edge of the element
665 5 case Position::OnEdge : {
666
667 // Variables
668 GeometricMeshRegion::ElementType eltTyp;
669 5 GeometricMeshRegion::BallMap shell;
670 felInt edgPosInShell;
671
672 // Get map edge to ver
673 5 auto& edge = GeometricMeshRegion::m_eltTypMapEdgToVer[GeometricMeshRegion::Tetra4][location.second];
674
675 // Get shell
676
2/4
✓ Branch 6 taken 5 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 5 times.
✗ Branch 10 not taken.
 5 m_volMesh->getEdgeShell({verOfTetra[edge[0]], verOfTetra[edge[1]]}, GeometricMeshRegion::Tetra4, idxOfElt, shell);
677
678 // Loop over the shell
679
2/2
✓ Branch 1 taken 15 times.
✓ Branch 2 taken 5 times.
 20 for (std::size_t iel = 0; iel < shell.size(); ++iel) {
680 15 std::tie(eltTyp, eltIdx, edgPosInShell) = shell[iel];
681
682 // Add element to the set of intersected elements
683
1/2
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
 15 setOfIntElt.insert(eltIdx);
684 }
685 5 break;
686 5 }
687 // Vertex is on a face of the element
688 ✗ case Position::OnFace : {
689
690 // Variables
691 ✗ std::vector<felInt> verFacOfElt(3,-1);
692 ✗ Face facOfElt;
693
694 ✗ auto& mapEltFacToVer = GeometricMeshRegion::m_eltTypMapFacToVer[GeometricMeshRegion::Tetra4];
695
696 // Get the ied-th edge of the element
697 ✗ verFacOfElt[0] = verOfTetra[ mapEltFacToVer[location.second][0] ];
698 ✗ verFacOfElt[1] = verOfTetra[ mapEltFacToVer[location.second][1] ];
699 ✗ verFacOfElt[2] = verOfTetra[ mapEltFacToVer[location.second][2] ];
700 ✗ m_volMesh->getFace(verFacOfElt, facOfElt);
701
702 // Loop over elements sharing the face
703 ✗ auto& listNeigh = facOfElt.listNeighVolumes();
704 ✗ for (std::size_t iel = 0; iel < listNeigh.size(); ++iel) {
705 ✗ eltIdx = listNeigh[iel]->idVolume();
706
707 // Add element to the set of intersected elements
708 ✗ setOfIntElt.insert(eltIdx);
709 }
710 ✗ break;
711 }
712 // Vertex is inside the element
713 ✗ case Position::Inside : {
714
715 // Add element to the set of intersected elements
716 ✗ setOfIntElt.insert(idxOfElt);
717 ✗ break;
718 }
719 ✗ default : {
720 ✗ FEL_ERROR("Intersector::addIntersetedEltToList3D : the vertex must be on a vertex, edge or inside the triangle !");
721 }
722 }
723 10 }
724
725 /***********************************************************************************/
726 /***********************************************************************************/
727
728 5 bool Intersector::computeIntersectionVertex3D(const felInt idxOfElt, const felInt idxOfVer, const felInt idxOfPreElt, felInt& idxOfNxtElt, const std::vector<Point>& pntOfEdge, std::vector<Point>& intPoint, bool& isInside) const
729 {
730 GeometricMeshRegion::ElementType eltTyp;
731 felInt eltIdx;
732 felInt verPos;
733 unsigned short int ifa;
734 int nbrIntPnt;
735
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 std::vector<Point*> pntOfNxtElt(4, nullptr);
736
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 std::vector<felInt> verOfNxtElt(4, 0);
737 std::array<double,4> bary;
738 std::array<Point*,3> pntFacOfElt;
739
740 5 GeometricMeshRegion::BallMap ball;
741
742 5 const std::array<const unsigned short int,4> mapVerToOppFac{2, 3, 1, 0};
743
744 5 auto& mapEltFacToVer = GeometricMeshRegion::m_eltTypMapFacToVer[GeometricMeshRegion::Tetra4];
745
746 // Get ball
747
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 m_volMesh->getVertexBall(idxOfVer, GeometricMeshRegion::Tetra4, idxOfElt, ball);
748
749 // Loop over ball
750
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 for (std::size_t iel = 0; iel < ball.size(); ++iel) {
751 5 std::tie(eltTyp, eltIdx, verPos) = ball[iel];
752
753
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if ( eltIdx == idxOfPreElt )
754 ✗ continue;
755
756 // Get nodes and coordinates of the current element
757
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 m_volMesh->getOneElement(GeometricMeshRegion::Tetra4, eltIdx, verOfNxtElt, pntOfNxtElt, 0);
758
759 // Check if pntOfEdge[1] is inside the current element
760 // If it is not inside, look for an intersection with the opposite face
761
2/4
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
 5 if ( !GeoUtilities::isPointInTetra(pntOfNxtElt, pntOfEdge[1], bary, m_tol) ) {
762
763 // Get the points of the face opposite to the vertex
764 5 ifa = mapVerToOppFac[verPos];
765
2/2
✓ Branch 1 taken 15 times.
✓ Branch 2 taken 5 times.
 20 for (std::size_t ive = 0; ive < pntFacOfElt.size(); ++ive)
766 15 pntFacOfElt[ive] = pntOfNxtElt[ mapEltFacToVer[ifa][ive] ];
767
768 // Compute intersection with opposite face
769
1/2
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
 5 nbrIntPnt = GeoUtilities::computeSegmentTriangleIntersection(pntOfEdge[0], pntOfEdge[1], *pntFacOfElt[0], *pntFacOfElt[1], *pntFacOfElt[2], intPoint[0], m_tol);
770
771 // If intersection
772
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 5 if ( nbrIntPnt > 0 ) {
773
774
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if ( m_verbose > 9 )
775 ✗ std::cout << "computeIntersectionVertex3D: Intersection detected in " << intPoint[0] << "." << std::endl;
776
777 // Return next element to intersect
778 5 idxOfNxtElt = eltIdx;
779
780 5 return true;
781 }
782 }
783 // If it is inside, nothing to do
784 else {
785
786 ✗ if ( m_verbose > 9 )
787 ✗ std::cout << "computeIntersectionVertex3D: Point " << pntOfEdge[1] << " has been found in element " << eltIdx << "." << std::endl;
788
789 // Set intersection
790 ✗ intPoint[0] = pntOfEdge[1];
791
792 // Return next element to intersect
793 ✗ idxOfNxtElt = eltIdx;
794
795 // Point inside element
796 ✗ isInside = true;
797
798 ✗ return true;
799 }
800 }
801
802 ✗ return false;
803 5 }
804
805 /***********************************************************************************/
806 /***********************************************************************************/
807
808 5 bool Intersector::computeIntersectionEdge3D(const felInt idxOfElt, const std::vector<felInt>& idxOfVer, const felInt idxOfPreElt, felInt& idxOfNxtElt, const std::vector<Point>& pntOfEdge, std::vector<Point>& intPoint, bool& isInside) const
809 {
810 GeometricMeshRegion::ElementType eltTyp;
811 felInt eltIdx;
812 felInt edgPos;
813 int nbrIntPnt;
814
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 std::vector<Point*> pntOfNxtElt(4, nullptr);
815
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 std::vector<felInt> verOfNxtElt(4, 0);
816 std::array<double,4> bary;
817 std::array<Point*,3> pntFacOfNxtElt;
818
819 5 GeometricMeshRegion::BallMap shell;
820
821 5 auto& mapEltFacToVer = GeometricMeshRegion::m_eltTypMapFacToVer[GeometricMeshRegion::Tetra4];
822 5 auto& mapEltEdgToFac = GeometricMeshRegion::m_eltTypMapEdgToFac[GeometricMeshRegion::Tetra4];
823
824 // Get shell
825
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 m_volMesh->getEdgeShell(idxOfVer, GeometricMeshRegion::Tetra4, idxOfElt, shell);
826
827 // Loop over the shell
828
1/2
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
 15 for (std::size_t iel = 0; iel < shell.size(); ++iel) {
829 15 std::tie(eltTyp, eltIdx, edgPos) = shell[iel];
830
831
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 10 times.
 15 if ( eltIdx == idxOfPreElt )
832 5 continue;
833
834 // Get nodes and coordinates of the current element
835
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 m_volMesh->getOneElement(GeometricMeshRegion::Tetra4, eltIdx, verOfNxtElt, pntOfNxtElt, 0);
836
837 // Check if pntOfEdge[1] is inside the current element
838 // If it is not inside, look for an intersection with the opposite edges
839
3/4
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 5 times.
✓ Branch 5 taken 5 times.
 10 if ( !GeoUtilities::isPointInTetra(pntOfNxtElt, pntOfEdge[1], bary, m_tol) ) {
840
841 // Loop over the two opposite faces
842
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 5 times.
 25 for (std::size_t ifa = 0; ifa < 4; ++ifa) {
843
844 // Skip face
845
6/6
✓ Branch 2 taken 15 times.
✓ Branch 3 taken 5 times.
✓ Branch 6 taken 5 times.
✓ Branch 7 taken 10 times.
✓ Branch 8 taken 10 times.
✓ Branch 9 taken 10 times.
 20 if ( ifa == mapEltEdgToFac[edgPos][0] || ifa == mapEltEdgToFac[edgPos][1] )
846 10 continue;
847
848 // Get the edge
849 10 pntFacOfNxtElt[0] = pntOfNxtElt[ mapEltFacToVer[ifa][0] ];
850 10 pntFacOfNxtElt[1] = pntOfNxtElt[ mapEltFacToVer[ifa][1] ];
851 10 pntFacOfNxtElt[2] = pntOfNxtElt[ mapEltFacToVer[ifa][2] ];
852
853 // Compute intersection with opposite edge
854
1/2
✓ Branch 7 taken 10 times.
✗ Branch 8 not taken.
 10 nbrIntPnt = GeoUtilities::computeSegmentTriangleIntersection(pntOfEdge[0], pntOfEdge[1], *pntFacOfNxtElt[0], *pntFacOfNxtElt[1], *pntFacOfNxtElt[2], intPoint[0], m_tol);
855
856 // If intersection
857
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 if ( nbrIntPnt > 0 ) {
858
859 ✗ if ( m_verbose > 9 )
860 ✗ std::cout << "computeIntersectionEdge3D: Intersection detected in " << intPoint[0] << "." << std::endl;
861
862 // Return next element to intersect
863 ✗ idxOfNxtElt = eltIdx;
864
865 ✗ return true;
866 }
867 }
868 }
869 // If it is inside, nothing to do
870 else {
871
872
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if ( m_verbose > 9 )
873 ✗ std::cout << "computeIntersectionEdge3D: Point " << pntOfEdge[1] << " has been found in element " << eltIdx << "." << std::endl;
874
875 // Set intersection
876 5 intPoint[0] = pntOfEdge[1];
877
878 // Return next element to intersect
879 5 idxOfNxtElt = eltIdx;
880
881 // Point inside element
882 5 isInside = true;
883
884 5 return true;
885 }
886 }
887
888 ✗ return false;
889 5 }
890
891 /***********************************************************************************/
892 /***********************************************************************************/
893
894 ✗ bool Intersector::computeIntersectionFace3D(const felInt idxOfElt, const int idxFacInt, const felInt idxOfPreElt, felInt& idxNxtElt, const std::vector<Point>& pntOfEdge, std::vector<Point>& intPoint, bool& isInside) const
895 {
896 felInt eltIdx;
897 int nbrIntPnt;
898 ✗ std::vector<Point*> pntOfNxtElt(4, nullptr);
899 ✗ std::vector<felInt> verOfNxtElt(4, 0);
900 ✗ std::vector<felInt> verOfElt(4, 0);
901 ✗ std::vector<felInt> verFacOfElt(3,-1);
902 std::array<double,4> bary;
903 std::array<felInt,3> verFacOfNxtElt;
904 std::array<Point*,3> pntFacOfNxtElt;
905 ✗ Face facOfElt;
906
907 ✗ auto& mapEltFacToVer = GeometricMeshRegion::m_eltTypMapFacToVer[GeometricMeshRegion::Tetra4];
908
909 // Get nodes and coordinates of idxOfElt
910 ✗ m_volMesh->getOneElement(GeometricMeshRegion::Tetra4, idxOfElt, verOfElt, 0);
911
912 // Get the ied-th edge of the triangle
913 ✗ verFacOfElt[0] = verOfElt[ mapEltFacToVer[idxFacInt][0] ];
914 ✗ verFacOfElt[1] = verOfElt[ mapEltFacToVer[idxFacInt][1] ];
915 ✗ verFacOfElt[2] = verOfElt[ mapEltFacToVer[idxFacInt][2] ];
916 ✗ m_volMesh->getFace(verFacOfElt, facOfElt);
917
918 // Loop over elements sharing the face
919 ✗ auto& listNeigh = facOfElt.listNeighVolumes();
920 ✗ for (std::size_t iel = 0; iel < listNeigh.size(); ++iel) {
921 ✗ eltIdx = listNeigh[iel]->idVolume();
922
923 ✗ if ( eltIdx == idxOfPreElt )
924 ✗ continue;
925
926 // Get nodes and coordinates of the current element
927 ✗ m_volMesh->getOneElement(GeometricMeshRegion::Tetra4, eltIdx, verOfNxtElt, pntOfNxtElt, 0);
928
929 // Check if pntOfEdge[1] is inside the current element
930 // If it is not inside, look for an intersection with the opposite edges
931 ✗ if ( !GeoUtilities::isPointInTetra(pntOfNxtElt, pntOfEdge[1], bary, m_tol) ) {
932
933 // Loop over the three opposite faces
934 ✗ for (std::size_t ifa = 0; ifa < 4; ++ifa) {
935
936 // Get vertices of the ied-th face
937 ✗ verFacOfNxtElt[0] = verOfNxtElt[ mapEltFacToVer[ifa][0] ];
938 ✗ verFacOfNxtElt[1] = verOfNxtElt[ mapEltFacToVer[ifa][1] ];
939 ✗ verFacOfNxtElt[2] = verOfNxtElt[ mapEltFacToVer[ifa][2] ];
940
941 // Skip edge
942 ✗ if ( std::is_permutation(verFacOfElt.begin(), verFacOfElt.end(), verFacOfNxtElt.begin()) )
943 ✗ continue;
944
945 // Get the face
946 ✗ pntFacOfNxtElt[0] = pntOfNxtElt[ mapEltFacToVer[ifa][0] ];
947 ✗ pntFacOfNxtElt[1] = pntOfNxtElt[ mapEltFacToVer[ifa][1] ];
948 ✗ pntFacOfNxtElt[2] = pntOfNxtElt[ mapEltFacToVer[ifa][2] ];
949
950 // Compute intersection with the edge
951 ✗ nbrIntPnt = GeoUtilities::computeSegmentTriangleIntersection(pntOfEdge[0], pntOfEdge[1], *pntFacOfNxtElt[0], *pntFacOfNxtElt[1], *pntFacOfNxtElt[2], intPoint[0], m_tol);
952
953 // If intersection
954 ✗ if ( nbrIntPnt > 0 ) {
955
956 ✗ if ( m_verbose > 9 )
957 ✗ std::cout << "computeIntersectionFace3D: Intersection detected in " << intPoint[0] << "." << std::endl;
958
959 // Return next element to intersect
960 ✗ idxNxtElt = eltIdx;
961
962 ✗ return true;
963 }
964 }
965 }
966 // If it is inside, nothing to do
967 else {
968
969 ✗ if ( m_verbose > 9 )
970 ✗ std::cout << "computeIntersectionFace3D: Point " << pntOfEdge[1] << " has been found in element " << eltIdx << "." << std::endl;
971
972 // Set intersection
973 ✗ intPoint[0] = pntOfEdge[1];
974
975 // Return next element to intersect
976 ✗ idxNxtElt = eltIdx;
977
978 // Point inside element
979 ✗ isInside = true;
980
981 ✗ return true;
982 }
983 }
984
985 ✗ return false;
986 }
987
988 /***********************************************************************************/
989 /***********************************************************************************/
990
991 ✗ bool Intersector::computeIntersectionInside3D(const felInt idxOfElt, felInt& idxNxtElt, const std::vector<Point>& pntOfEdge, std::vector<Point>& intPoint) const
992 {
993 int nbrIntPnt;
994 ✗ std::vector<Point*> pntOfNxtElt(4, nullptr);
995 ✗ std::vector<felInt> verOfNxtElt(4, 0);
996 std::array<Point*,3> pntFacOfElt;
997
998 ✗ auto& mapEltFacToVer = GeometricMeshRegion::m_eltTypMapFacToVer[GeometricMeshRegion::Tetra4];
999
1000 // Get nodes and coordinates of the current triangle
1001 ✗ m_volMesh->getOneElement(GeometricMeshRegion::Tetra4, idxOfElt, verOfNxtElt, pntOfNxtElt, 0);
1002
1003 // Loop over the faces
1004 ✗ for (std::size_t ifa = 0; ifa < 4; ++ifa) {
1005
1006 // Get the face
1007 ✗ pntFacOfElt[0] = pntOfNxtElt[ mapEltFacToVer[ifa][0] ];
1008 ✗ pntFacOfElt[1] = pntOfNxtElt[ mapEltFacToVer[ifa][1] ];
1009 ✗ pntFacOfElt[2] = pntOfNxtElt[ mapEltFacToVer[ifa][2] ];
1010
1011 // Compute intersection with the face
1012 ✗ nbrIntPnt = GeoUtilities::computeSegmentTriangleIntersection(pntOfEdge[0], pntOfEdge[1], *pntFacOfElt[0], *pntFacOfElt[1], *pntFacOfElt[2], intPoint[0], m_tol);
1013
1014 // If intersection
1015 ✗ if ( nbrIntPnt > 0 ) {
1016
1017 // Return next element to intersect
1018 ✗ idxNxtElt = idxOfElt;
1019
1020 ✗ return true;
1021 }
1022 }
1023
1024 ✗ return false;
1025 }
1026
1027 /***********************************************************************************/
1028 /***********************************************************************************/
1029
1030 ✗ felInt Intersector::getGermElement3D(const std::vector<felInt>& verOfEdge) const
1031 {
1032 ✗ std::vector<Point*> pntOfTetra(4, nullptr);
1033 ✗ std::vector<felInt> verOfTetra(4, 0);
1034 std::array<double, 4> bary;
1035
1036 // Go through all elements of the mesh
1037 ✗ for (felInt k = 0; k < m_volMesh->numElements(GeometricMeshRegion::Tetra4); ++k) {
1038
1039 // Get nodes of the current triangle
1040 ✗ m_volMesh->getOneElement(GeometricMeshRegion::Tetra4, k, verOfTetra, pntOfTetra, 0);
1041
1042 ✗ if( GeoUtilities::isPointInTetra( pntOfTetra, m_strMesh->listPoints()[verOfEdge[0]], bary, m_tol) )
1043 ✗ return k;
1044 }
1045
1046 ✗ FEL_ERROR("Intersector::getGermElement3D cannot localize the point into the volume mesh.");
1047 ✗ return -1;
1048 }
1049
1050 }
1051