GCC Code Coverage Report

Directory: ./
File: Geometry/geometricMeshRegion.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 858 1150 74.6%
Branches: 669 1349 49.6%
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: J.Castelneau & J.Foulon & V.Martin
13 //
14
15 // System includes
16 #include <numeric>
17 #include <queue>
18
19 // External includes
20
21 // Project includes
22 #include "Core/array_1d.hpp"
23 #include "Core/felisceParam.hpp"
24 #include "FiniteElement/geoElement.hpp"
25 #include "Geometry/geometricEdges.hpp"
26 #include "Geometry/geometricFaces.hpp"
27 #include "Geometry/geometricMeshRegion.hpp"
28 #include "Geometry/neighFacesOfEdges.hpp"
29 #include "Geometry/neighVolumesOfEdges.hpp"
30 #include "Geometry/neighVolumesOfFaces.hpp"
31 #include "Geometry/geo_utilities.hpp"
32 #include "Tools/math_utilities.hpp"
33
34 namespace felisce
35 {
36
37 const int GeometricMeshRegion::m_numPointsPerElt[] = {
38 1, //0D
39 2,3,3, //1D
40 3,4,6, 4,5,6,8,9, //2D
41 4,5,10, //3D
42 5,13, //3D
43 6,9,15, //3D
44 8,9,20,26,27 //3D
45 };
46
47 const int GeometricMeshRegion::m_numVerticesPerElt[] = {
48 1, //0D
49 2,2,2, //1D
50 3,3,3, 4,4,4,4, //2D
51 4,4,4, //3D
52 5,5, //3D
53 6,6,6, //3D
54 8,8,8,8,8 //3D
55 };
56
57 const int GeometricMeshRegion::m_numEdgesPerElt[] = {
58 0, //0D
59 1,1,1, //1D
60 3,3,3, 4,4,4,4,4, //2D
61 6,6,6, //3D
62 0,0, //3D
63 9,9,9, //3D
64 12,12,12,12,12 //3D
65 };
66
67 const int GeometricMeshRegion::m_numFacesPerElt[] = {
68 0, //0D
69 0,0,0, //1D
70 1,1,1, 1,1,1,1,1, //2D
71 4,4,4, //3D
72 0,0, //3D
73 5,5,5, //3D
74 6,6,6,6,6 //3D
75 };
76
77 GeometricMeshRegion::EnumToFelNameGeoEle_type GeometricMeshRegion::eltEnumToFelNameGeoEle = GeometricMeshRegion::EnumToFelNameGeoEle_type({
78 std::make_pair("Nodes" , &geoElementNode) // 0D
79 , std::make_pair("Segments2", &geoElementSegmentP1) // 1D
80 , std::make_pair("Segments3b", &geoElementSegmentP1b)
81 , std::make_pair("Segments3", &geoElementSegmentP2)
82 , std::make_pair("Triangles3", &geoElementTriangleP1) // 2D
83 , std::make_pair("Triangles4", &geoElementTriangleP1b)
84 , std::make_pair("Triangles6", &geoElementTriangleP2)
85 , std::make_pair("Quadrangles4", &geoElementQuadrangleQ1)
86 , std::make_pair("Quadrangles5", &geoElementQuadrangleQ1b)
87 , std::make_pair("Quadrangles6", &geoElementQuadrangleP1xP2)
88 , std::make_pair("Quadrangles8", &geoElementQuadrangleQ2)
89 , std::make_pair("Quadrangles9", &geoElementQuadrangleQ2c)
90 , std::make_pair("Tetrahedra4", &geoElementTetrahedronP1) // 3D
91 , std::make_pair("Tetrahedra5", &geoElementTetrahedronP1b)
92 , std::make_pair("Tetrahedra10", &geoElementTetrahedronP2)
93 , std::make_pair("Pyramids5", &geoElementNULL)
94 , std::make_pair("Pyramids13", &geoElementNULL)
95 , std::make_pair("Prisms6", &geoElementPrismR1)
96 , std::make_pair("Prisms9", &geoElementPrismP1xP2)
97 , std::make_pair("Prisms15", &geoElementPrismR2)
98 , std::make_pair("Hexahedra8", &geoElementHexahedronQ1)
99 , std::make_pair("Hexahedra9", &geoElementHexahedronQ1b)
100 , std::make_pair("Hexahedra20", &geoElementHexahedronQ2)
101 , std::make_pair("Hexahedra26", &geoElementNULL)
102 , std::make_pair("Hexahedra27", &geoElementHexahedronQ2c)
103 });
104
105 GeometricMeshRegion::StringToElementType_type GeometricMeshRegion::eltFelNameToEnum = {
106 std::make_pair("Nodes" , GeometricMeshRegion::Nod)
107 , std::make_pair("Segments2", GeometricMeshRegion::Seg2)
108 , std::make_pair("Segments3b", GeometricMeshRegion::Seg3b)
109 , std::make_pair("Segments3", GeometricMeshRegion::Seg3)
110 , std::make_pair("Triangles3", GeometricMeshRegion::Tria3)
111 , std::make_pair("Triangles4", GeometricMeshRegion::Tria4)
112 , std::make_pair("Triangles6", GeometricMeshRegion::Tria6)
113 , std::make_pair("Quadrangles4", GeometricMeshRegion::Quad4)
114 , std::make_pair("Quadrangles5", GeometricMeshRegion::Quad5)
115 , std::make_pair("Quadrangles6", GeometricMeshRegion::Quad6)
116 , std::make_pair("Quadrangles8", GeometricMeshRegion::Quad8)
117 , std::make_pair("Quadrangles9", GeometricMeshRegion::Quad9)
118 , std::make_pair("Tetrahedra4", GeometricMeshRegion::Tetra4)
119 , std::make_pair("Tetrahedra5", GeometricMeshRegion::Tetra5)
120 , std::make_pair("Tetrahedra10", GeometricMeshRegion::Tetra10)
121 , std::make_pair("Pyramids5", GeometricMeshRegion::Pyram5)
122 , std::make_pair("Pyramids13", GeometricMeshRegion::Pyram13)
123 , std::make_pair("Prisms6", GeometricMeshRegion::Prism6)
124 , std::make_pair("Prisms9", GeometricMeshRegion::Prism9)
125 , std::make_pair("Prisms15", GeometricMeshRegion::Prism15)
126 , std::make_pair("Hexahedra8", GeometricMeshRegion::Hexa8)
127 , std::make_pair("Hexahedra9", GeometricMeshRegion::Hexa9)
128 , std::make_pair("Hexahedra20", GeometricMeshRegion::Hexa20)
129 , std::make_pair("Hexahedra26", GeometricMeshRegion::Hexa26)
130 , std::make_pair("Hexahedra27", GeometricMeshRegion::Hexa27)
131 };
132
133 GeometricMeshRegion::ElementTypeLinearToDifferentElementType_type GeometricMeshRegion::eltLinearToEltQuad = {
134 std::make_pair(GeometricMeshRegion::Seg2, GeometricMeshRegion::Seg3)
135 , std::make_pair(GeometricMeshRegion::Tria3, GeometricMeshRegion::Tria6)
136 , std::make_pair(GeometricMeshRegion::Quad4, GeometricMeshRegion::Quad8)
137 , std::make_pair(GeometricMeshRegion::Quad6, GeometricMeshRegion::Quad8)
138 , std::make_pair(GeometricMeshRegion::Tetra4, GeometricMeshRegion::Tetra10)
139 , std::make_pair(GeometricMeshRegion::Pyram5, GeometricMeshRegion::Pyram13)
140 , std::make_pair(GeometricMeshRegion::Prism6, GeometricMeshRegion::Prism15)
141 , std::make_pair(GeometricMeshRegion::Prism9, GeometricMeshRegion::Prism15)
142 , std::make_pair(GeometricMeshRegion::Hexa8, GeometricMeshRegion::Hexa20)
143 };
144
145 GeometricMeshRegion::ElementTypeLinearToDifferentElementType_type GeometricMeshRegion::eltLinearToEltBubble = {
146 std::make_pair(GeometricMeshRegion::Seg2, GeometricMeshRegion::Seg3b)
147 , std::make_pair(GeometricMeshRegion::Tria3, GeometricMeshRegion::Tria4)
148 , std::make_pair(GeometricMeshRegion::Quad4, GeometricMeshRegion::Quad5)
149 , std::make_pair(GeometricMeshRegion::Tetra4, GeometricMeshRegion::Tetra5)
150 , std::make_pair(GeometricMeshRegion::Pyram5, GeometricMeshRegion::Pyram5)
151 , std::make_pair(GeometricMeshRegion::Prism6, GeometricMeshRegion::Prism9)
152 , std::make_pair(GeometricMeshRegion::Hexa8, GeometricMeshRegion::Hexa9)
153 };
154
155 std::vector<GeometricMeshRegion::ElementType> GeometricMeshRegion::bagElementType0D = { GeometricMeshRegion::Nod };
156
157 std::vector<GeometricMeshRegion::ElementType> GeometricMeshRegion::bagElementType1D = {
158 GeometricMeshRegion::Seg2
159 , GeometricMeshRegion::Seg3b
160 , GeometricMeshRegion::Seg3
161 };
162
163 std::vector<GeometricMeshRegion::ElementType> GeometricMeshRegion::bagElementType2D = {
164 GeometricMeshRegion::Tria3
165 , GeometricMeshRegion::Tria4
166 , GeometricMeshRegion::Tria6
167 , GeometricMeshRegion::Quad4
168 , GeometricMeshRegion::Quad5
169 , GeometricMeshRegion::Quad6
170 , GeometricMeshRegion::Quad8
171 , GeometricMeshRegion::Quad9
172 };
173
174 std::vector<GeometricMeshRegion::ElementType> GeometricMeshRegion::bagElementType3D = {
175 GeometricMeshRegion::Tetra4
176 , GeometricMeshRegion::Tetra5
177 , GeometricMeshRegion::Tetra10
178 , GeometricMeshRegion::Prism6
179 , GeometricMeshRegion::Prism9
180 , GeometricMeshRegion::Prism15
181 , GeometricMeshRegion::Hexa8
182 , GeometricMeshRegion::Hexa9
183 , GeometricMeshRegion::Hexa20
184 , GeometricMeshRegion::Hexa26
185 , GeometricMeshRegion::Hexa27
186 };
187
188 std::vector<GeometricMeshRegion::ElementType> GeometricMeshRegion::bagElementTypeLinear = {
189 GeometricMeshRegion::Seg2
190 , GeometricMeshRegion::Tria3
191 , GeometricMeshRegion::Quad4
192 , GeometricMeshRegion::Quad6
193 , GeometricMeshRegion::Tetra4
194 , GeometricMeshRegion::Pyram5
195 , GeometricMeshRegion::Prism6
196 , GeometricMeshRegion::Prism9
197 , GeometricMeshRegion::Hexa8
198 };
199
200 GeometricMeshRegion::NameLabelToListLabel_type GeometricMeshRegion::descriptionLineEnsightToListLabel;
201
202 const std::vector< std::vector< std::vector<unsigned short int> > > GeometricMeshRegion::m_eltTypMapVerToFac =
203 {
204 {}, // "Nodes"
205 { {0} }, // "Segments2"
206 { {0} }, // "Segments3b"
207 { {0} }, // "Segments3"
208 { {2,0}, {0,1}, {1,2} }, // "Triangles3"
209 { {2,0}, {0,1}, {1,2} }, // "Triangles4"
210 { {2,0}, {0,1}, {1,2} }, // "Triangles6"
211 { {3,0}, {0,1}, {1,2}, {2,3} }, // "Quadrangles4"
212 { {3,0}, {0,1}, {1,2}, {2,3} }, // "Quadrangles5"
213 { {3,0}, {0,1}, {1,2}, {2,3} }, // "Quadrangles6"
214 { {3,0}, {0,1}, {1,2}, {2,3} }, // "Quadrangles8"
215 { {3,0}, {0,1}, {1,2}, {2,3} }, // "Quadrangles9"
216 { {0,1,3}, {0,1,2}, {0,2,3}, {1,2,3} }, // "Tetrahedra4"
217 { {0,1,3}, {0,1,2}, {0,2,3}, {1,2,3} }, // "Tetrahedra5"
218 { {0,1,3}, {0,1,2}, {0,2,3}, {1,2,3} }, // "Tetrahedra10"
219 {}, // "Pyramids5"
220 {}, // "Pyramids13"
221 {}, // "Prisms6"
222 {}, // "Prisms9"
223 {}, // "Prisms15"
224 {}, // "Hexahedra8"
225 {}, // "Hexahedra9"
226 {}, // "Hexahedra20"
227 {}, // "Hexahedra26"
228 {} // "Hexahedra27"
229 }; // TODO replace this ugly structure defining a class for every element type + polymorphism
230
231 const std::vector< std::vector< std::vector<unsigned short int> > > GeometricMeshRegion::m_eltTypMapFacToVer =
232 {
233 { {0} }, // "Nodes"
234 { {0,1} }, // "Segments2"
235 { {0,1} }, // "Segments3b"
236 { {0,1} }, // "Segments3"
237 { {0,1}, {1,2}, {2,0} }, // "Triangles3"
238 { {0,1}, {1,2}, {2,0} }, // "Triangles4"
239 { {0,1}, {1,2}, {2,0} }, // "Triangles6"
240 { {0,1}, {1,2}, {2,3}, {3,0} }, // "Quadrangles4"
241 { {0,1}, {1,2}, {2,3}, {3,0} }, // "Quadrangles5"
242 { {0,1}, {1,2}, {2,3}, {3,0} }, // "Quadrangles6"
243 { {0,1}, {1,2}, {2,3}, {3,0} }, // "Quadrangles8"
244 { {0,1}, {1,2}, {2,3}, {3,0} }, // "Quadrangles9"
245 { {0,1,2}, {0,3,1}, {1,3,2}, {0,2,3} }, // "Tetrahedra4"
246 { {0,1,2}, {0,3,1}, {1,3,2}, {0,2,3} }, // "Tetrahedra5"
247 { {0,1,2}, {0,3,1}, {1,3,2}, {0,2,3} }, // "Tetrahedra10"
248 {}, // "Pyramids5"
249 {}, // "Pyramids13"
250 {}, // "Prisms6"
251 {}, // "Prisms9"
252 {}, // "Prisms15"
253 {}, // "Hexahedra8"
254 {}, // "Hexahedra9"
255 {}, // "Hexahedra20"
256 {}, // "Hexahedra26"
257 {} // "Hexahedra27"
258 }; // TODO replace this ugly structure defining a class for every element type + polymorphism
259
260 const std::vector< std::vector< std::vector<unsigned short int> > > GeometricMeshRegion::m_eltTypMapEdgToVer =
261 {
262 { }, // "Nodes"
263 { }, // "Segments2"
264 { }, // "Segments3b"
265 { }, // "Segments3"
266 { }, // "Triangles3"
267 { }, // "Triangles4"
268 { }, // "Triangles6"
269 { }, // "Quadrangles4"
270 { }, // "Quadrangles5"
271 { }, // "Quadrangles6"
272 { }, // "Quadrangles8"
273 { }, // "Quadrangles9"
274 { {0,1}, {0,2}, {0,3}, {1,2}, {1,3}, {2,3} }, // "Tetrahedra4"
275 { {0,1}, {0,2}, {0,3}, {1,2}, {1,3}, {2,3} }, // "Tetrahedra5"
276 { {0,1}, {0,2}, {0,3}, {1,2}, {1,3}, {2,3} }, // "Tetrahedra10"
277 {}, // "Pyramids5"
278 {}, // "Pyramids13"
279 {}, // "Prisms6"
280 {}, // "Prisms9"
281 {}, // "Prisms15"
282 {}, // "Hexahedra8"
283 {}, // "Hexahedra9"
284 {}, // "Hexahedra20"
285 {}, // "Hexahedra26"
286 {} // "Hexahedra27"
287 }; // TODO replace this ugly structure defining a class for every element type + polymorphism
288
289
290
291 const std::vector< std::vector< std::vector<unsigned short int> > > GeometricMeshRegion::m_eltTypMapEdgToFac =
292 {
293 { }, // "Nodes"
294 { }, // "Segments2"
295 { }, // "Segments3b"
296 { }, // "Segments3"
297 { {0}, {0}, {0} }, // "Triangles3"
298 { {0}, {0}, {0} }, // "Triangles4"
299 { {0}, {0}, {0} }, // "Triangles6"
300 { {0}, {0}, {0}, {0} }, // "Quadrangles4"
301 { {0}, {0}, {0}, {0} }, // "Quadrangles5"
302 { {0}, {0}, {0}, {0} }, // "Quadrangles6"
303 { {0}, {0}, {0}, {0} }, // "Quadrangles8"
304 { {0}, {0}, {0}, {0} }, // "Quadrangles9"
305 { {0,1}, {0,3}, {1,3}, {0,2}, {1,2}, {2,3} }, // "Tetrahedra4"
306 { {0,1}, {0,3}, {1,3}, {0,2}, {1,2}, {2,3} }, // "Tetrahedra5"
307 { {0,1}, {0,3}, {1,3}, {0,2}, {1,2}, {2,3} }, // "Tetrahedra10"
308 {}, // "Pyramids5"
309 {}, // "Pyramids13"
310 {}, // "Prisms6"
311 {}, // "Prisms9"
312 {}, // "Prisms15"
313 {}, // "Hexahedra8"
314 {}, // "Hexahedra9"
315 {}, // "Hexahedra20"
316 {}, // "Hexahedra26"
317 {} // "Hexahedra27"
318 }; // TODO replace this ugly structure defining a class for every element type + polymorphism
319
320
321
322
323
324
325
326
327
328
329 /***********************************************************************************/
330 /***********************************************************************************/
331
332
2/2
✓ Branch 1 taken 31925 times.
✓ Branch 2 taken 1277 times.
 33202 GeometricMeshRegion::GeometricMeshRegion()
333 {
334 // Initialization of tab m_numElements
335
2/2
✓ Branch 0 taken 31925 times.
✓ Branch 1 taken 1277 times.
 33202 for (int ityp = 0; ityp < m_numTypesOfElement; ++ityp) {
336 31925 m_numElements[ityp] = 0;
337 31925 m_listElements[ityp] = nullptr;
338 }
339 1277 }
340
341 /***********************************************************************************/
342 /***********************************************************************************/
343
344 33202 GeometricMeshRegion::~GeometricMeshRegion()
345 {
346
2/2
✓ Branch 0 taken 31925 times.
✓ Branch 1 taken 1277 times.
 33202 for (int ityp = 0; ityp < m_numTypesOfElement; ++ityp)
347
2/2
✓ Branch 0 taken 2464 times.
✓ Branch 1 taken 29461 times.
 31925 if ( m_numElements[ityp] != 0 )
348
1/2
✓ Branch 0 taken 2464 times.
✗ Branch 1 not taken.
 2464 delete [] m_listElements[ityp];
349
3/4
✓ Branch 11 taken 1277 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 1277 times.
✓ Branch 14 taken 31925 times.
 34479 }
350
351 /***********************************************************************************/
352 /***********************************************************************************/
353
354 515 void GeometricMeshRegion::setDomainDim()
355 {
356
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 515 times.
 515 if ( m_domainDim != GeoMeshUndefined )
357 ✗ FEL_ERROR("Problem: you can std::set the mesh dimension only once!");
358
359
6/6
✓ Branch 4 taken 4377 times.
✓ Branch 5 taken 330 times.
✓ Branch 6 taken 4192 times.
✓ Branch 7 taken 185 times.
✓ Branch 8 taken 4192 times.
✓ Branch 9 taken 515 times.
 4707 for (auto it_eltype = bagElementType3D.begin(); it_eltype != bagElementType3D.end() && m_domainDim == GeoMeshUndefined ; ++it_eltype) {
360
2/2
✓ Branch 1 taken 185 times.
✓ Branch 2 taken 4007 times.
 4192 if ( m_numElements[*it_eltype] > 0 )
361 185 m_domainDim = GeoMesh3D;
362 }
363
364
6/6
✓ Branch 4 taken 1427 times.
✓ Branch 5 taken 81 times.
✓ Branch 6 taken 993 times.
✓ Branch 7 taken 434 times.
✓ Branch 8 taken 993 times.
✓ Branch 9 taken 515 times.
 1508 for (auto it_eltype = bagElementType2D.begin(); it_eltype != bagElementType2D.end() && m_domainDim == GeoMeshUndefined ; ++it_eltype) {
365
2/2
✓ Branch 1 taken 249 times.
✓ Branch 2 taken 744 times.
 993 if ( m_numElements[*it_eltype] > 0 )
366 249 m_domainDim = GeoMesh2D;
367 }
368
369
5/6
✓ Branch 4 taken 596 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 81 times.
✓ Branch 7 taken 515 times.
✓ Branch 8 taken 81 times.
✓ Branch 9 taken 515 times.
 596 for (auto it_eltype = bagElementType1D.begin(); it_eltype != bagElementType1D.end() && m_domainDim == GeoMeshUndefined ; ++it_eltype) {
370
1/2
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
 81 if ( m_numElements[*it_eltype] > 0 )
371 81 m_domainDim = GeoMesh1D;
372 }
373
374
3/6
✓ Branch 4 taken 515 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 515 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 515 times.
 515 for (auto it_eltype = bagElementType0D.begin(); it_eltype != bagElementType0D.end() && m_domainDim == GeoMeshUndefined ; ++it_eltype) {
375 ✗ if ( m_numElements[*it_eltype] > 0 )
376 ✗ m_domainDim = GeoMesh0D;
377 }
378 515 }
379
380 /***********************************************************************************/
381 /***********************************************************************************/
382
383 2818 void GeometricMeshRegion::allocateElements(const ElementType& eltType)
384 {
385
386
1/2
✓ Branch 0 taken 2818 times.
✗ Branch 1 not taken.
 2818 m_listElements[eltType] = new felInt[ m_numElements[eltType] * m_numPointsPerElt[eltType] ];
387 2818 }
388
389 /***********************************************************************************/
390 /***********************************************************************************/
391
392 11 void GeometricMeshRegion::deleteElements()
393 {
394
2/2
✓ Branch 0 taken 275 times.
✓ Branch 1 taken 11 times.
 286 for (int ityp = 0; ityp < m_numTypesOfElement; ++ityp)
395
2/2
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 253 times.
 275 if ( m_numElements[ityp] != 0 )
396
1/2
✓ Branch 0 taken 22 times.
✗ Branch 1 not taken.
 22 delete [] m_listElements[ityp];
397 11 }
398
399 /***********************************************************************************/
400 /***********************************************************************************/
401
402 13010 void GeometricMeshRegion::reorderListElePerRef(const std::map<int, std::vector<felInt>>& RefToElements, const ElementType eltType)
403 {
404 felInt tmpSizeElements, begElement, nbElements;
405 int refNumber;
406
407
2/2
✓ Branch 1 taken 1476 times.
✓ Branch 2 taken 11534 times.
 13010 if ( RefToElements.size() > 0 ) {
408
409 1476 tmpSizeElements = m_numElements[eltType] * m_numPointsPerElt[eltType];
410
411
2/4
✓ Branch 0 taken 1476 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 1476 times.
✗ Branch 5 not taken.
 1476 felInt *tmpElements = new felInt[ tmpSizeElements ];
412 1476 begElement = 0;
413 1476 nbElements = 0;
414
2/2
✓ Branch 3 taken 3118 times.
✓ Branch 4 taken 1476 times.
 4594 for(auto it_ref = RefToElements.begin(); it_ref !=RefToElements.end(); ++it_ref) {
415 3118 refNumber = it_ref->first;
416
417 // Retrieve all elements having the given ref (their IDs are given by it_ref->second)
418
1/2
✓ Branch 5 taken 3118 times.
✗ Branch 6 not taken.
 3118 getElements(eltType, &(it_ref->second[0]), it_ref->second.size(), tmpElements, tmpSizeElements, begElement);
419
420 3118 nbElements = it_ref->second.size();
421
2/4
✓ Branch 1 taken 3118 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3118 times.
✗ Branch 5 not taken.
 3118 intRefToBegEndMaps[eltType][refNumber] = std::make_pair(begElement, nbElements);
422
423 // Shift the starting position
424 3118 begElement += it_ref->second.size() * m_numPointsPerElt[eltType];
425 }
426
1/2
✓ Branch 0 taken 1476 times.
✗ Branch 1 not taken.
 1476 delete[] m_listElements[eltType];
427
428 // swap the lists.
429 1476 m_listElements[eltType] = tmpElements;
430 }
431 13010 }
432
433 /***********************************************************************************/
434 /***********************************************************************************/
435
436 1427 void GeometricMeshRegion::setBagElementTypeDomain()
437 {
438 1427 m_bagElementTypeDomain.clear();
439
440
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 170 times.
✓ Branch 2 taken 538 times.
✓ Branch 3 taken 719 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 1427 switch ( m_domainDim ) {
441
442 ✗ case GeoMesh0D: // 0D
443 ✗ for (auto it_eltype = bagElementType0D.begin(); it_eltype != bagElementType0D.end(); ++it_eltype)
444 ✗ if ( m_numElements[*it_eltype] > 0 )
445 ✗ m_bagElementTypeDomain.push_back(*it_eltype);
446
447 ✗ break;
448
449 170 case GeoMesh1D: // 1D
450
2/2
✓ Branch 4 taken 510 times.
✓ Branch 5 taken 170 times.
 680 for (auto it_eltype = bagElementType1D.begin(); it_eltype != bagElementType1D.end(); ++it_eltype)
451
2/2
✓ Branch 1 taken 170 times.
✓ Branch 2 taken 340 times.
 510 if ( m_numElements[*it_eltype] > 0 )
452
1/2
✓ Branch 2 taken 170 times.
✗ Branch 3 not taken.
 170 m_bagElementTypeDomain.push_back(*it_eltype);
453
454 170 break;
455
456 538 case GeoMesh2D: // 2D
457
2/2
✓ Branch 4 taken 4304 times.
✓ Branch 5 taken 538 times.
 4842 for (auto it_eltype = bagElementType2D.begin(); it_eltype != bagElementType2D.end(); ++it_eltype)
458
2/2
✓ Branch 1 taken 528 times.
✓ Branch 2 taken 3776 times.
 4304 if ( m_numElements[*it_eltype] > 0 )
459
1/2
✓ Branch 2 taken 528 times.
✗ Branch 3 not taken.
 528 m_bagElementTypeDomain.push_back(*it_eltype);
460
461 538 break;
462
463 719 case GeoMesh3D: // 3D
464
2/2
✓ Branch 4 taken 7909 times.
✓ Branch 5 taken 719 times.
 8628 for (auto it_eltype = bagElementType3D.begin(); it_eltype != bagElementType3D.end(); ++it_eltype)
465
2/2
✓ Branch 1 taken 691 times.
✓ Branch 2 taken 7218 times.
 7909 if ( m_numElements[*it_eltype] > 0 )
466
1/2
✓ Branch 2 taken 691 times.
✗ Branch 3 not taken.
 691 m_bagElementTypeDomain.push_back(*it_eltype);
467
468 719 break;
469
470 ✗ case GeoMeshUndefined:
471 ✗ FEL_ERROR("GeometricMeshRegion: you must define the mesh dimension!");
472
473 ✗ break;
474 }
475 1427 }
476
477 /***********************************************************************************/
478 /***********************************************************************************/
479
480 1427 void GeometricMeshRegion::setBagElementTypeDomainBoundary()
481 {
482 1427 const auto& r_instance = FelisceParam::instance();
483 1427 const bool extendToLowerDimensions = r_instance.extendToLowerDimensions;
484 1427 m_bagElementTypeDomainBoundary.clear();
485
486
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 170 times.
✓ Branch 2 taken 538 times.
✓ Branch 3 taken 719 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 1427 switch ( m_domainDim ) {
487
488 ✗ case GeoMesh0D: // 0D
489
490 ✗ break;
491
492 170 case GeoMesh1D: // 1D
493
2/2
✓ Branch 4 taken 170 times.
✓ Branch 5 taken 170 times.
 340 for (auto it_eltype = bagElementType0D.begin(); it_eltype != bagElementType0D.end(); ++it_eltype)
494
2/2
✓ Branch 1 taken 112 times.
✓ Branch 2 taken 58 times.
 170 if ( m_numElements[*it_eltype] > 0 )
495
1/2
✓ Branch 2 taken 112 times.
✗ Branch 3 not taken.
 112 m_bagElementTypeDomainBoundary.push_back(*it_eltype);
496
497 170 break;
498
499 538 case GeoMesh2D: // 2D
500
501 // TODO
502 // if ( extendToLowerDimensions ) {
503 // for (auto it_eltype = bagElementType0D.begin(); it_eltype != bagElementType0D.end(); ++it_eltype)
504 // if ( m_numElements[*it_eltype] > 0 )
505 // m_bagElementTypeDomainBoundary.push_back(*it_eltype);
506 // }
507
508
2/2
✓ Branch 4 taken 1614 times.
✓ Branch 5 taken 538 times.
 2152 for (auto it_eltype = bagElementType1D.begin(); it_eltype != bagElementType1D.end(); ++it_eltype)
509
2/2
✓ Branch 1 taken 447 times.
✓ Branch 2 taken 1167 times.
 1614 if ( m_numElements[*it_eltype] > 0 )
510
1/2
✓ Branch 2 taken 447 times.
✗ Branch 3 not taken.
 447 m_bagElementTypeDomainBoundary.push_back(*it_eltype);
511
512 538 break;
513
514 719 case GeoMesh3D: // 3D
515
516
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 695 times.
 719 if ( extendToLowerDimensions ) {
517
2/2
✓ Branch 4 taken 72 times.
✓ Branch 5 taken 24 times.
 96 for (auto it_eltype = bagElementType1D.begin(); it_eltype != bagElementType1D.end(); ++it_eltype)
518
2/2
✓ Branch 1 taken 15 times.
✓ Branch 2 taken 57 times.
 72 if ( m_numElements[*it_eltype] > 0 )
519
1/2
✓ Branch 2 taken 15 times.
✗ Branch 3 not taken.
 15 m_bagElementTypeDomainBoundary.push_back(*it_eltype);
520 }
521
522
2/2
✓ Branch 4 taken 5752 times.
✓ Branch 5 taken 719 times.
 6471 for (auto it_eltype = bagElementType2D.begin(); it_eltype != bagElementType2D.end(); ++it_eltype)
523
2/2
✓ Branch 1 taken 677 times.
✓ Branch 2 taken 5075 times.
 5752 if ( m_numElements[*it_eltype] > 0 )
524
1/2
✓ Branch 2 taken 677 times.
✗ Branch 3 not taken.
 677 m_bagElementTypeDomainBoundary.push_back(*it_eltype);
525
526 719 break;
527
528 ✗ case GeoMeshUndefined:
529 ✗ FEL_ERROR("GeometricMeshRegion: you must define the mesh dimension!");
530
531 ✗ break;
532 }
533 1427 }
534
535 /***********************************************************************************/
536 /***********************************************************************************/
537
538 548 void GeometricMeshRegion::setLocalMesh(GeometricMeshRegion& meshGlobal, const std::vector<int>& eltPartition,
539 const int rank, std::vector<felInt>& loc2globElem)
540 {
541 548 m_builtEdges = false;
542 548 m_builtFaces = false;
543 548 m_flagFormatMesh = meshGlobal.m_flagFormatMesh;
544 548 m_listPoints.clear();
545 548 m_listNormals.clear();
546 548 m_listTangents.clear();
547 548 m_numCoor = meshGlobal.numCoor();
548 548 m_domainDim = meshGlobal.domainDim();
549 548 m_createNormalTangent = meshGlobal.createNormalTangent();
550
1/2
✓ Branch 1 taken 548 times.
✗ Branch 2 not taken.
 548 m_listTangents.resize(m_domainDim);
551
2/2
✓ Branch 0 taken 13700 times.
✓ Branch 1 taken 548 times.
 14248 for (int ityp = 0; ityp < m_numTypesOfElement; ++ityp) {
552 13700 m_numElements[ityp] = 0;
553 }
554
555 // Copy of all mesh points
556
1/2
✓ Branch 2 taken 548 times.
✗ Branch 3 not taken.
 548 m_listPoints = meshGlobal.listPoints();
557
2/2
✓ Branch 0 taken 94 times.
✓ Branch 1 taken 454 times.
 548 if (m_createNormalTangent) {
558
1/2
✓ Branch 2 taken 94 times.
✗ Branch 3 not taken.
 94 m_listNormals = meshGlobal.listNormals();
559
2/2
✓ Branch 0 taken 168 times.
✓ Branch 1 taken 94 times.
 262 for (int i = 0 ; i<m_domainDim; i++) {
560
1/2
✓ Branch 3 taken 168 times.
✗ Branch 4 not taken.
 168 m_listTangents[i] = meshGlobal.m_listTangents[i];
561 }
562 }
563
564 // Allocation of memory for DOMAIN m_listElements and BOUNDARY m_listElements
565 548 felInt countEltTot = 0;
566 548 auto& r_element_type_domain = meshGlobal.bagElementTypeDomain();
567
1/2
✓ Branch 2 taken 548 times.
✗ Branch 3 not taken.
 548 m_allocateListElementsByRef(meshGlobal.bagElementTypeDomain(), meshGlobal,
568 eltPartition, rank, countEltTot);
569
570 548 auto& r_element_type_domain_boundary = meshGlobal.bagElementTypeDomainBoundary();
571
1/2
✓ Branch 1 taken 548 times.
✗ Branch 2 not taken.
 548 m_allocateListElementsByRef(r_element_type_domain_boundary, meshGlobal,
572 eltPartition, rank, countEltTot);
573
574 // Reset to 0 the element counter
575 548 countEltTot = 0;
576
577 // Reset intRefToBegEndMaps for domain an boundary
578
2/2
✓ Branch 4 taken 548 times.
✓ Branch 5 taken 548 times.
 1096 for (auto it_eltype = r_element_type_domain.begin(); it_eltype != r_element_type_domain.end(); ++it_eltype) {
579 548 ElementType eltType = static_cast<ElementType>(*it_eltype);
580 548 intRefToBegEndMaps[eltType].clear();
581 }
582
2/2
✓ Branch 4 taken 617 times.
✓ Branch 5 taken 548 times.
 1165 for (auto it_eltype = r_element_type_domain_boundary.begin(); it_eltype != r_element_type_domain_boundary.end(); ++it_eltype) {
583 617 ElementType eltType = static_cast<ElementType>(*it_eltype);
584 617 intRefToBegEndMaps[eltType].clear();
585 }
586
587 // Fill of m_listElements and map intRefToBegEndMaps
588
1/2
✓ Branch 1 taken 548 times.
✗ Branch 2 not taken.
 548 m_fillListElementsByRef(r_element_type_domain, meshGlobal, eltPartition, rank, loc2globElem, countEltTot);
589
590 // Fill of m_listElements and map intRefToBegEndMaps
591
1/2
✓ Branch 1 taken 548 times.
✗ Branch 2 not taken.
 548 m_fillListElementsByRef(r_element_type_domain_boundary, meshGlobal, eltPartition, rank, loc2globElem, countEltTot);
592
593 // Set dim and bags
594
1/2
✓ Branch 1 taken 548 times.
✗ Branch 2 not taken.
 548 setBagElementTypeDomain();
595
1/2
✓ Branch 1 taken 548 times.
✗ Branch 2 not taken.
 548 setBagElementTypeDomainBoundary();
596 548 }
597
598 /***********************************************************************************/
599 /***********************************************************************************/
600
601 114 void GeometricMeshRegion::buildEdges()
602 {
603
2/2
✓ Branch 0 taken 93 times.
✓ Branch 1 taken 21 times.
 114 if ( !m_builtEdges ) {
604 93 m_listEdges.clearAndInit( this->numPoints() );
605
606 //! order of the list of edges: first 1d, then 2d then 3d.
607 //! (given edges first, then edges of faces, then edges of volumes.)
608 93 m_buildEdgesPerBag(bagElementType1D, 1);
609 93 m_listEdges.setNumGivenEdges(); // right after having inserted the given edges
610 93 m_buildEdgesPerBag(bagElementType2D, 2);
611 93 m_buildEdgesPerBag(bagElementType3D, 3);
612
613 93 felInt numElem1D = getNumElement1D();
614
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 93 times.
 93 if ( m_listEdges.numGivenEdges() < numElem1D ) {
615 ✗ std::cout << "We found " << m_listEdges.numGivenEdges() << " edges, while " << numElem1D << " edges were provided in the mesh!\n";
616 ✗ FEL_WARNING("Beware: some duplicated edges in the mesh were discarded. (Or there is a problem in the buildEdge!)");
617 }
618
619
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 93 times.
 93 if ( m_listEdges.numGivenEdges() > numElem1D ) {
620 ✗ std::cout << "We found " << m_listEdges.numGivenEdges() << " edges, while " << numElem1D << " edges were provided in the mesh!\n";
621 ✗ FEL_ERROR("There is a problem in the buildEdge: the number of given edges is not correct!");
622 }
623
624 93 m_builtEdges = true;
625 }
626 114 }
627
628 /***********************************************************************************/
629 /***********************************************************************************/
630
631 177 void GeometricMeshRegion::buildFaces()
632 {
633
2/2
✓ Branch 0 taken 133 times.
✓ Branch 1 taken 44 times.
 177 if ( !m_builtFaces ) {
634 133 m_listFaces.clearAndInit(this->numPoints());
635
636 //! order of the list of faces: first 2d, then 3d.
637 //! (given faces first, then faces of volumes.)
638 133 m_buildFacesPerBag(bagElementType2D, 2);
639 133 m_listFaces.setNumGivenFaces();// right after having inserted the given faces
640 133 m_buildFacesPerBag(bagElementType3D, 3);
641
642 133 felInt numElem2D = getNumElement2D();
643
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 133 times.
 133 if ( m_listFaces.numGivenFaces() < numElem2D ) {
644 ✗ std::cout << "We found " << m_listFaces.numGivenFaces() << " faces, while " << numElem2D << " faces were provided in the mesh!\n";
645 ✗ FEL_WARNING("Beware: some duplicated faces in the mesh were discarded. (Or there is a problem in the buildFace!)");
646 }
647
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 133 times.
 133 if ( m_listFaces.numGivenFaces() > numElem2D ) {
648 ✗ std::cout << "We found " << m_listFaces.numGivenFaces() << " faces, while " << numElem2D << " faces were provided in the mesh!\n";
649 ✗ FEL_ERROR("There is a problem in the buildFace: the number of given faces is not correct!");
650 }
651
652 133 m_builtFaces = true;
653 }
654 177 }
655
656 /***********************************************************************************/
657 /***********************************************************************************/
658
659 995 void GeometricMeshRegion::moveMesh(const std::vector<double>& disp, const felReal& coef)
660 {
661
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 995 times.
 995 if ( static_cast<int>(disp.size()) != m_numCoor * this->numPoints() ) {
662 ✗ PetscPrintf(MpiInfo::petscComm(),"Displacement size = %ld vs numCoor * numPoints = %d \n",static_cast<unsigned long>(disp.size()), m_numCoor * numPoints());
663 ✗ FEL_ERROR("GeometricMeshRegion: We can not move the mesh with this displacement\n" );
664 }
665
666
2/2
✓ Branch 0 taken 51 times.
✓ Branch 1 taken 944 times.
 995 if ( !m_moved ) {
667 // We store the reference position of the mesh
668 51 m_moved = true;
669 51 m_listReferencePoints.resize( numPoints() );
670
2/2
✓ Branch 1 taken 233966 times.
✓ Branch 2 taken 51 times.
 234017 for (int i = 0; i < numPoints(); ++i)
671 233966 m_listReferencePoints[i] = m_listPoints[i];
672 }
673
674
2/2
✓ Branch 1 taken 4184373 times.
✓ Branch 2 taken 995 times.
 4185368 for (int i = 0; i < numPoints() ; ++i)
675
2/2
✓ Branch 0 taken 12450186 times.
✓ Branch 1 taken 4184373 times.
 16634559 for (int j = 0; j < m_numCoor; ++j)
676 12450186 m_listPoints[i].coor(j) = m_listReferencePoints[i].coor(j) + coef*disp[j + m_numCoor*i];
677 995 }
678
679 /***********************************************************************************/
680 /***********************************************************************************/
681 // TODO this function should be rewritten introducing a method "computeNormal" for every shape
682 148 void GeometricMeshRegion::computeNormalTangent(const InitialNormalTangentProvided initialNormalTangentProvided)
683 {
684 // Initialization
685 148 felInt numPointPerElement = 0;
686 148 double norm = 0.;
687 148 const std::size_t number_of_points = m_listPoints.size();
688
1/2
✓ Branch 2 taken 148 times.
✗ Branch 3 not taken.
 148 std::vector<double> meas(number_of_points, 0.0);
689 148 double eltMeas = 0.;
690 148 std::vector<felInt> listPt;
691 ElementType eltType;
692 148 ElementType previous_eltType = ELEMENT_TYPE_COUNTER;
693
694 // Allocate values
695
1/2
✓ Branch 1 taken 148 times.
✗ Branch 2 not taken.
 148 allocateListNormalsAndTangents(initialNormalTangentProvided);
696
697
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 146 times.
 148 const auto& r_bag_elements = m_domainDim == 3 ? m_bagElementTypeDomainBoundary : m_bagElementTypeDomain;
698
2/2
✓ Branch 0 taken 19 times.
✓ Branch 1 taken 129 times.
 148 if (m_domainDim == 1) {
699
700
2/2
✓ Branch 0 taken 17 times.
✓ Branch 1 taken 2 times.
 19 if ( m_numCoor == 2 ){
701 // Compute the normals on each elements
702 // Loop over the elements
703
2/2
✓ Branch 1 taken 17 times.
✓ Branch 2 taken 17 times.
 34 for(std::size_t itype=0; itype<r_bag_elements.size(); ++itype) {
704 17 eltType = r_bag_elements[itype];
705 17 numPointPerElement = m_numPointsPerElt[eltType];
706
1/2
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
 17 listPt.resize(numPointPerElement);
707
708
1/2
✓ Branch 0 taken 17 times.
✗ Branch 1 not taken.
 17 if ( eltType == Seg2 ) {
709 17 Point normal_tmp(0.0, 0.0, 0.0);
710 double norm;
711
2/2
✓ Branch 0 taken 590 times.
✓ Branch 1 taken 17 times.
 607 for(felInt iel=0; iel<m_numElements[eltType]; ++iel) {
712
2/2
✓ Branch 0 taken 1180 times.
✓ Branch 1 taken 590 times.
 1770 for(felInt ipt=0; ipt<numPointPerElement; ++ipt) {
713 1180 listPt[ipt] = m_listElements[eltType][numPointPerElement * iel + ipt];
714 }
715 590 const auto& r_point_0 = m_listPoints[listPt[0]];
716 590 const auto& r_point_1 = m_listPoints[listPt[1]];
717 590 normal_tmp.x() = -(r_point_1.y() - r_point_0.y());
718 590 normal_tmp.y() = r_point_1.x() - r_point_0.x();
719
1/2
✓ Branch 1 taken 590 times.
✗ Branch 2 not taken.
 590 norm = normal_tmp.norm();
720
1/2
✓ Branch 1 taken 590 times.
✗ Branch 2 not taken.
 590 if (norm > std::numeric_limits<double>::epsilon()) {
721 590 normal_tmp /= norm;
722 }
723 590 eltMeas = 0.0;
724
2/2
✓ Branch 0 taken 1180 times.
✓ Branch 1 taken 590 times.
 1770 for(std::size_t icoor=0; icoor<2; ++icoor) {
725 1180 eltMeas += std::pow(r_point_0[icoor] - r_point_1[icoor], 2);
726 }
727 590 eltMeas = std::sqrt(eltMeas);
728
2/2
✓ Branch 0 taken 1180 times.
✓ Branch 1 taken 590 times.
 1770 for(felInt ipt=0; ipt<numPointPerElement; ++ipt) {
729 1180 const auto& id = listPt[ipt];
730
2/2
✓ Branch 0 taken 2360 times.
✓ Branch 1 taken 1180 times.
 3540 for(std::size_t icoor=0; icoor<2; ++icoor) {
731 2360 m_listNormals[id][icoor] += normal_tmp[icoor] * eltMeas;
732 }
733 1180 meas[id] += eltMeas;
734 }
735 }
736
737 // Compute the normal at node
738
2/2
✓ Branch 1 taken 606 times.
✓ Branch 2 taken 17 times.
 623 for(std::size_t i=0; i<m_listNormals.size(); ++i) {
739
2/2
✓ Branch 0 taken 1212 times.
✓ Branch 1 taken 606 times.
 1818 for(felInt icoor=0; icoor<2; ++icoor) {
740
1/2
✓ Branch 2 taken 1212 times.
✗ Branch 3 not taken.
 1212 if (meas[i] > std::numeric_limits<double>::epsilon()) {
741 1212 m_listNormals[i][icoor] /= meas[i];
742 }
743 }
744
1/2
✓ Branch 2 taken 606 times.
✗ Branch 3 not taken.
 606 norm = m_listNormals[i].norm();
745
1/2
✓ Branch 1 taken 606 times.
✗ Branch 2 not taken.
 606 if (norm > std::numeric_limits<double>::epsilon()) {
746 606 m_listNormals[i] /= norm;
747 }
748 }
749
750 // Zero z coordinate
751
2/2
✓ Branch 0 taken 606 times.
✓ Branch 1 taken 17 times.
 623 for(std::size_t i=0; i<number_of_points; ++i) {
752 606 m_listNormals[i].z() = 0.0;
753 }
754 // Compute the tangent at node from the normal
755
2/2
✓ Branch 0 taken 606 times.
✓ Branch 1 taken 17 times.
 623 for(std::size_t i=0; i<number_of_points; ++i) {
756 606 m_listTangents[0][i].x() = m_listNormals[i].y();
757 606 m_listTangents[0][i].y() = -m_listNormals[i].x();
758 606 m_listTangents[0][i].z() = 0.0;
759 }
760 }
761 }
762
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 } else if (m_numCoor == 3) {
763 // Compute the normals on each elements
764 // Loop over the elements
765
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
 4 for(std::size_t itype=0; itype<r_bag_elements.size(); ++itype) {
766 2 eltType = r_bag_elements[itype];
767 2 numPointPerElement = m_numPointsPerElt[eltType];
768
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 listPt.resize(numPointPerElement);
769
770
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if ( eltType == Seg2 ) {
771 2 Point tangent_tmp(0.0, 0.0, 0.0);
772 double norm;
773
2/2
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 2 times.
 24 for(felInt iel=0; iel<m_numElements[eltType]; ++iel) {
774
2/2
✓ Branch 0 taken 44 times.
✓ Branch 1 taken 22 times.
 66 for(felInt ipt=0; ipt<numPointPerElement; ++ipt) {
775 44 listPt[ipt] = m_listElements[eltType][numPointPerElement * iel + ipt];
776 }
777 22 const auto& r_point_0 = m_listPoints[listPt[0]];
778 22 const auto& r_point_1 = m_listPoints[listPt[1]];
779 22 tangent_tmp.x() = r_point_1.x() - r_point_0.x();
780 22 tangent_tmp.y() = r_point_1.y() - r_point_0.y();
781 22 tangent_tmp.z() = r_point_1.z() - r_point_0.z();
782
1/2
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
 22 norm = tangent_tmp.norm();
783
1/2
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
 22 if (norm > std::numeric_limits<double>::epsilon()) {
784 22 tangent_tmp /= norm;
785 }
786 22 eltMeas = 0.0;
787
2/2
✓ Branch 0 taken 66 times.
✓ Branch 1 taken 22 times.
 88 for(std::size_t icoor=0; icoor<3; ++icoor) {
788 66 eltMeas += std::pow(r_point_0[icoor] - r_point_1[icoor], 2);
789 }
790 22 eltMeas = std::sqrt(eltMeas);
791
2/2
✓ Branch 0 taken 44 times.
✓ Branch 1 taken 22 times.
 66 for(felInt ipt=0; ipt<numPointPerElement; ++ipt) {
792 44 const auto& id = listPt[ipt];
793
2/2
✓ Branch 0 taken 132 times.
✓ Branch 1 taken 44 times.
 176 for(std::size_t icoor=0; icoor<3; ++icoor) {
794 132 m_listTangents[0][id][icoor] += tangent_tmp[icoor] * eltMeas;
795 }
796 44 meas[id] += eltMeas;
797 }
798 }
799
800 // Compute the tangent at node
801
2/2
✓ Branch 2 taken 24 times.
✓ Branch 3 taken 2 times.
 26 for(std::size_t i=0; i<m_listTangents[0].size(); ++i) {
802
2/2
✓ Branch 0 taken 72 times.
✓ Branch 1 taken 24 times.
 96 for(felInt icoor=0; icoor<3; ++icoor) {
803
1/2
✓ Branch 2 taken 72 times.
✗ Branch 3 not taken.
 72 if (meas[i] > std::numeric_limits<double>::epsilon()) {
804 72 m_listTangents[0][i][icoor] /= meas[i];
805 }
806 }
807
1/2
✓ Branch 3 taken 24 times.
✗ Branch 4 not taken.
 24 norm = m_listTangents[0][i].norm();
808
1/2
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
 24 if (norm > std::numeric_limits<double>::epsilon()) {
809 24 m_listTangents[0][i] /= norm;
810 }
811 }
812 // Compute the orthonormal base
813
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 array_1d<double, 3> normal, tangent_1, tangent_2;
814
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 if (initialNormalTangentProvided == InitialNormalTangentProvided::NONE) { // Compute whole base
815
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 1 times.
 13 for(std::size_t i=0; i<number_of_points; ++i) {
816
1/2
✓ Branch 3 taken 12 times.
✗ Branch 4 not taken.
 12 tangent_1.data() = m_listTangents[0][i].getCoor();
817
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 MathUtilities::OrthonormalBasis(tangent_1, normal, tangent_2);
818
1/2
✓ Branch 3 taken 12 times.
✗ Branch 4 not taken.
 12 m_listNormals[i].getCoor() = normal.data();
819
1/2
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
 12 m_listTangents[1][i].getCoor() = tangent_2.data();
820 }
821
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 } else if (initialNormalTangentProvided == InitialNormalTangentProvided::NORMAL) {
822 ✗ for(std::size_t i=0; i<number_of_points; ++i) {
823 ✗ normal.data() = m_listNormals[i].getCoor();
824 ✗ tangent_1.data() = m_listTangents[0][i].getCoor();
825 ✗ MathUtilities::UnitCrossProduct(tangent_2, tangent_1, normal);
826 ✗ m_listTangents[1][i].getCoor() = tangent_2.data();
827 }
828
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 } else if (initialNormalTangentProvided == InitialNormalTangentProvided::TANGENT_2) {
829
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 1 times.
 13 for(std::size_t i=0; i<number_of_points; ++i) {
830
1/2
✓ Branch 3 taken 12 times.
✗ Branch 4 not taken.
 12 tangent_1.data() = m_listTangents[0][i].getCoor();
831
1/2
✓ Branch 3 taken 12 times.
✗ Branch 4 not taken.
 12 tangent_2.data() = m_listTangents[1][i].getCoor();
832
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 MathUtilities::UnitCrossProduct(normal, tangent_2, tangent_1);
833
1/2
✓ Branch 3 taken 12 times.
✗ Branch 4 not taken.
 12 m_listNormals[i].getCoor() = normal.data();
834 }
835 } else {
836 ✗ FEL_ERROR("Tangent 1 cannot be provided, must be normal or tangent 2")
837 }
838 2 }
839 }
840 }
841 } else {
842
843 // Compute the normals on each elements
844 // Loop over the elements
845
2/2
✓ Branch 1 taken 129 times.
✓ Branch 2 taken 129 times.
 258 for(std::size_t itype=0; itype<r_bag_elements.size(); ++itype) {
846 258 eltType = r_bag_elements[itype];
847 129 numPointPerElement = m_numPointsPerElt[eltType];
848
1/2
✓ Branch 1 taken 129 times.
✗ Branch 2 not taken.
 129 listPt.resize(numPointPerElement);
849
850 // Compute normal and tangents in 3D
851 129 int iq[2][numPointPerElement];
852 129 std::vector<felInt> listPtQ1Elem;
853 std::array<Point, 2> vecTangentialPlan; // < a1,a2 > define the plan gives by the triangle
854
855
2/2
✓ Branch 0 taken 43914 times.
✓ Branch 1 taken 129 times.
 44043 for(felInt iel=0; iel<m_numElements[eltType]; ++iel) {
856
2/2
✓ Branch 0 taken 132090 times.
✓ Branch 1 taken 43914 times.
 176004 for(felInt ipt=0; ipt<numPointPerElement; ++ipt)
857 132090 listPt[ipt] = m_listElements[eltType][numPointPerElement * iel + ipt];
858
859 // ElementTriangleP1 iq = { {1,2,0},{2,0,1} };
860 // ElementQuadrangleQ1 iq = { {1,2,3,0},{3,0,1,2} };
861
3/4
✓ Branch 0 taken 348 times.
✓ Branch 1 taken 43566 times.
✓ Branch 2 taken 348 times.
✗ Branch 3 not taken.
 43914 if (eltType == Tria3 || eltType == Quad4) {
862
2/2
✓ Branch 0 taken 88176 times.
✓ Branch 1 taken 43914 times.
 132090 for (int l=0 ; l<numPointPerElement-1 ; l++) {
863 88176 iq[0][l]=l+1;
864 88176 iq[1][l+1]=l;
865 }
866 43914 iq[0][numPointPerElement-1]=0;
867 43914 iq[1][0]=numPointPerElement-1;
868
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 43914 } else if (eltType == Quad9) {
869 ✗ listPtQ1Elem.resize(16);
870 ✗ for (int i=0 ; i<4 ; i++) {
871 ✗ listPtQ1Elem[4*i] = listPt[i];
872 ✗ listPtQ1Elem[4*i+1] = listPt[4+i];
873 ✗ listPtQ1Elem[4*i+2] = listPt[8];
874 ✗ listPtQ1Elem[4*i+3] = listPt[3+i];
875 }
876 ✗ listPtQ1Elem[3] = listPt[7];
877
878 ✗ for (int l=0 ; l<3 ; l++) {
879 ✗ iq[0][l]=l+1;
880 ✗ iq[1][l+1]=l;
881 }
882 ✗ iq[0][4-1]=0;
883 ✗ iq[1][0]=4-1;
884 } else {
885 ✗ if (eltType != previous_eltType) {
886 ✗ std::cout << "WARNING: Not compatible geometry type: " << eltType << std::endl;
887 }
888 ✗ previous_eltType = eltType;
889 ✗ continue;
890 }
891
892
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 43914 times.
 43914 if (eltType == Quad9) {
893 ✗ for (int k=0 ; k<4 ; k++) {
894 ✗ for (int i=0 ; i<4 ; i++) {
895
896 ✗ for (int j=0 ; j<2 ; j++)
897 ✗ vecTangentialPlan[j] = m_listPoints[ listPtQ1Elem[4*k+iq[j][i]] ] - m_listPoints[ listPtQ1Elem[4*k+i] ];
898
899 ✗ m_listNormals[listPtQ1Elem[4*k+i]].x() += vecTangentialPlan[0].y()*vecTangentialPlan[1].z()-vecTangentialPlan[0].z()*vecTangentialPlan[1].y();
900 ✗ m_listNormals[listPtQ1Elem[4*k+i]].y() += vecTangentialPlan[0].z()*vecTangentialPlan[1].x()-vecTangentialPlan[0].x()*vecTangentialPlan[1].z();
901 ✗ m_listNormals[listPtQ1Elem[4*k+i]].z() += vecTangentialPlan[0].x()*vecTangentialPlan[1].y()-vecTangentialPlan[0].y()*vecTangentialPlan[1].x();
902 }
903 }
904 } else {
905 // Tria3 || Quad4
906
2/2
✓ Branch 0 taken 132090 times.
✓ Branch 1 taken 43914 times.
 176004 for (int i = 0; i < numPointPerElement; ++i) {
907
908
2/2
✓ Branch 0 taken 264180 times.
✓ Branch 1 taken 132090 times.
 396270 for (int j = 0 ; j < 2; ++j)
909
1/2
✓ Branch 5 taken 264180 times.
✗ Branch 6 not taken.
 264180 vecTangentialPlan[j] = m_listPoints[ listPt[ iq[j][i] ] ] - m_listPoints[ listPt[i] ];
910
911 132090 m_listNormals[listPt[i]].x() += vecTangentialPlan[0].y()*vecTangentialPlan[1].z()-vecTangentialPlan[0].z()*vecTangentialPlan[1].y();
912 132090 m_listNormals[listPt[i]].y() += vecTangentialPlan[0].z()*vecTangentialPlan[1].x()-vecTangentialPlan[0].x()*vecTangentialPlan[1].z();
913 132090 m_listNormals[listPt[i]].z() += vecTangentialPlan[0].x()*vecTangentialPlan[1].y()-vecTangentialPlan[0].y()*vecTangentialPlan[1].x();
914 }
915 }
916 }
917 258 }
918
919 // Normalize normals
920
2/2
✓ Branch 1 taken 24687 times.
✓ Branch 2 taken 129 times.
 24816 for(std::size_t i=0; i<m_listNormals.size(); i++) {
921
1/2
✓ Branch 2 taken 24687 times.
✗ Branch 3 not taken.
 24687 const double norm_normal = m_listNormals[i].norm();
922
2/2
✓ Branch 1 taken 23940 times.
✓ Branch 2 taken 747 times.
 24687 if (norm_normal > std::numeric_limits<double>::epsilon()) {
923 23940 m_listNormals[i] /= norm_normal;
924 }
925 }
926 129 std::unordered_map<std::string, int> mapOfTypeDirectionTangent;
927
2/4
✓ Branch 2 taken 129 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 129 times.
✗ Branch 6 not taken.
 129 mapOfTypeDirectionTangent["e1"] = 1;
928
2/4
✓ Branch 2 taken 129 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 129 times.
✗ Branch 6 not taken.
 129 mapOfTypeDirectionTangent["e2"] = 2;
929
2/4
✓ Branch 2 taken 129 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 129 times.
✗ Branch 6 not taken.
 129 mapOfTypeDirectionTangent["e3"] = 3;
930
2/4
✓ Branch 1 taken 129 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 129 times.
✗ Branch 5 not taken.
 129 m_choiceDirectionTangent = mapOfTypeDirectionTangent[FelisceParam::instance().choiceDirectionTangent];
931 std::array<double,3> e;
932
933 // Compute V1
934
2/2
✓ Branch 1 taken 24687 times.
✓ Branch 2 taken 129 times.
 24816 for(std::size_t i=0; i<m_listPoints.size(); ++i) {
935
3/4
✓ Branch 2 taken 24687 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 23940 times.
✓ Branch 6 taken 747 times.
 24687 if (m_listNormals[i].norm() > std::numeric_limits<double>::epsilon()) {
936 23940 bool firstTangentValidated = false;
937 //test the good choice of std::vector e
938
2/2
✓ Branch 0 taken 23956 times.
✓ Branch 1 taken 23940 times.
 47896 while (!firstTangentValidated) {
939
4/4
✓ Branch 0 taken 11697 times.
✓ Branch 1 taken 120 times.
✓ Branch 2 taken 2292 times.
✓ Branch 3 taken 9847 times.
 23956 switch (m_choiceDirectionTangent) {
940 11697 case 1 :
941 11697 e[0] = 1.; e[1] = 0.; e[2] = 0.;
942 11697 break;
943 120 case 2 :
944 120 e[0] = 0.; e[1] = 1.; e[2] = 0.;
945 120 break;
946 2292 case 3 :
947 2292 e[0] = 0.; e[1] = 0.; e[2] = 1.;
948 2292 break;
949 9847 default:
950 9847 e[0] = 1.; e[1] = 0.; e[2] = 0.;
951 9847 break;
952 }
953
954 23956 m_listTangents[0][i].x() = m_listNormals[i].y()*e[2]-m_listNormals[i].z()*e[1];
955 23956 m_listTangents[0][i].y() = m_listNormals[i].z()*e[0]-m_listNormals[i].x()*e[2];
956 23956 m_listTangents[0][i].z() = m_listNormals[i].x()*e[1]-m_listNormals[i].y()*e[0];
957
1/2
✓ Branch 3 taken 23956 times.
✗ Branch 4 not taken.
 23956 norm = m_listTangents[0][i].norm();
958
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 23940 times.
 23956 if (norm<0.001) {
959 16 m_choiceDirectionTangent++;
960 } else {
961 23940 firstTangentValidated = true;
962 }
963 }
964 23940 m_listTangents[0][i] /= norm;
965 //compute V2
966 23940 m_listTangents[1][i].x() = m_listNormals[i].y()*m_listTangents[0][i].z()-m_listNormals[i].z()*m_listTangents[0][i].y();
967 23940 m_listTangents[1][i].y() = m_listNormals[i].z()*m_listTangents[0][i].x()-m_listNormals[i].x()*m_listTangents[0][i].z();
968 23940 m_listTangents[1][i].z() = m_listNormals[i].x()*m_listTangents[0][i].y()-m_listNormals[i].y()*m_listTangents[0][i].x();
969
1/2
✓ Branch 3 taken 23940 times.
✗ Branch 4 not taken.
 23940 m_listTangents[1][i] /= m_listTangents[1][i].norm();
970 }
971 }
972 129 }
973 148 }
974
975 /***********************************************************************************/
976 /***********************************************************************************/
977
978 185 void GeometricMeshRegion::allocateListNormalsAndTangents(const InitialNormalTangentProvided initialNormalTangentProvided)
979 {
980 // TODO: Normals and Tangents should be arrays instead of Points for efficiency
981 185 const std::size_t number_of_points = this->numPoints();
982
2/2
✓ Branch 1 taken 183 times.
✓ Branch 2 taken 2 times.
 185 if (m_listNormals.size() != number_of_points) {
983
1/2
✓ Branch 2 taken 183 times.
✗ Branch 3 not taken.
 183 m_listNormals.resize(number_of_points, Point(0.0));
984
2/2
✓ Branch 2 taken 45 times.
✓ Branch 3 taken 138 times.
 183 const int aux_size = this->domainDim() > this->numCoor() - 1 ? this->domainDim() : this->numCoor() - 1;
985
2/4
✓ Branch 3 taken 183 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 183 times.
✗ Branch 7 not taken.
 183 m_listTangents.resize(aux_size, std::vector<Point>(number_of_points,0.0));
986 } else {
987
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 if (initialNormalTangentProvided == InitialNormalTangentProvided::NONE) {
988
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 std::fill(m_listNormals.begin(), m_listNormals.end(), 0.);
989
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 2 std::for_each(m_listTangents.begin(), m_listTangents.end(), [](auto& listTan){ std::fill(listTan.begin(), listTan.end(), 0.); } );
990
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 } else if (initialNormalTangentProvided == InitialNormalTangentProvided::NORMAL) {
991 ✗ std::for_each(m_listTangents.begin(), m_listTangents.end(), [](auto& listTan){ std::fill(listTan.begin(), listTan.end(), 0.); } );
992
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 } else if (initialNormalTangentProvided == InitialNormalTangentProvided::TANGENT_2) {
993
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 std::fill(m_listNormals.begin(), m_listNormals.end(), 0.);
994
1/2
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 1 std::fill(m_listTangents[0].begin(), m_listTangents[0].end(), 0.);
995 } else{
996 ✗ FEL_ERROR("Tangent 1 cannot be provided, must be normal or tangent 2");
997 }
998 }
999 185 }
1000
1001 /***********************************************************************************/
1002 /***********************************************************************************/
1003
1004 24 void GeometricMeshRegion::computeElementNormal(std::vector<std::vector<Point> >& listElementNormals) const
1005 {
1006 // Initialization
1007 felInt numPointPerElement;
1008 ElementType eltType;
1009 24 std::vector<felInt> veridxOfElt;
1010 24 std::vector<Point*> pointsOfElt;
1011
1012 // Allocate values
1013
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 20 times.
 24 if ( listElementNormals.size() != static_cast<std::size_t>(m_numTypesOfElement) )
1014
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 listElementNormals.resize(m_numTypesOfElement);
1015
1016 // Loop over all the domain elements
1017
2/2
✓ Branch 1 taken 24 times.
✓ Branch 2 taken 24 times.
 48 for(std::size_t itype = 0; itype < m_bagElementTypeDomain.size(); ++itype) {
1018
1019 // get element type
1020 24 eltType = m_bagElementTypeDomain[itype];
1021
1022 // get number point per elementtype
1023 24 numPointPerElement = m_numPointsPerElt[eltType];
1024
1025 // resize vector
1026
1/2
✓ Branch 2 taken 24 times.
✗ Branch 3 not taken.
 24 listElementNormals[eltType].resize(m_numElements[eltType]);
1027
1/2
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
 24 veridxOfElt.resize(numPointPerElement);
1028
1/2
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
 24 pointsOfElt.resize(numPointPerElement);
1029
1030 // Compute normal for every "eltType" element
1031
2/2
✓ Branch 0 taken 432 times.
✓ Branch 1 taken 24 times.
 456 for(felInt iel = 0; iel < m_numElements[eltType]; ++iel) {
1032
1033 // Get element vertices
1034
1/2
✓ Branch 1 taken 432 times.
✗ Branch 2 not taken.
 432 getOneElement(eltType, iel, veridxOfElt, pointsOfElt, 0);
1035
1036 // Compute normal
1037
1/2
✓ Branch 3 taken 432 times.
✗ Branch 4 not taken.
 432 GeoUtilities::computeElementNormal(eltType, pointsOfElt, listElementNormals[eltType][iel]);
1038 }
1039 }
1040 24 }
1041
1042 /***********************************************************************************/
1043 /***********************************************************************************/
1044
1045 3118 void GeometricMeshRegion::getElements(const ElementType& eltType, const felInt* indexList, felInt numIndexList,
1046 felInt* array, int sizeArray, felInt startPosArray) const
1047 {
1048
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3118 times.
 3118 if ( m_numElements[eltType] == 0 ) {
1049 ✗ std::cout << "Warning: trying to access a type of element with no elements: " << eltEnumToFelNameGeoEle[eltType].first << "!" << std::endl;
1050 ✗ return;
1051 }
1052
1053 3118 const int nVpEl = m_numPointsPerElt[eltType];
1054
1055 // Check if all the elements can fit in array
1056
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3118 times.
 3118 FEL_ASSERT( sizeArray >= numIndexList*nVpEl + startPosArray );
1057
1058
2/2
✓ Branch 0 taken 2391211 times.
✓ Branch 1 taken 3118 times.
 2394329 for (felInt ii = 0; ii < numIndexList; ++ii) {
1059 2391211 felInt kind = indexList[ii];
1060
1061
2/4
✓ Branch 0 taken 2391211 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2391211 times.
 2391211 FEL_ASSERT(kind >= 0 && kind < m_numElements[eltType]);
1062
1063
2/2
✓ Branch 0 taken 8826786 times.
✓ Branch 1 taken 2391211 times.
 11217997 for (int jj = 0; jj < nVpEl; ++jj)
1064 8826786 array[ startPosArray + ii*nVpEl + jj ] = m_listElements[eltType][ kind*nVpEl + jj];
1065 }
1066 }
1067
1068 /***********************************************************************************/
1069 /***********************************************************************************/
1070
1071 ✗ void GeometricMeshRegion::setElements(const ElementType& eltType, const felInt* indexList,
1072 felInt numIndexList, const felInt* array, int sizeArray)
1073 {
1074 ✗ const int nVpEl = m_numPointsPerElt[eltType];
1075
1076 ✗ FEL_ASSERT(sizeArray == numIndexList*nVpEl);
1077
1078 ✗ for (felInt ii = 0; ii < numIndexList; ++ii) {
1079 ✗ felInt kind = indexList[ii];
1080
1081 ✗ FEL_ASSERT(kind >= 0 && kind < m_numElements[eltType]);
1082
1083 ✗ for (felInt jj = 0; jj< m_numPointsPerElt[eltType]; ++jj)
1084 ✗ m_listElements[eltType][ kind*nVpEl + jj ] = array[ ii*nVpEl + jj ];
1085 }
1086 }
1087
1088 /***********************************************************************************/
1089 /***********************************************************************************/
1090
1091 47979429 void GeometricMeshRegion::getOneElement(const ElementType& eltType, const felInt index, std::vector<felInt>& array,
1092 const felInt startindex, const bool flagIncrement) const
1093 {
1094 47979429 const int nVpEl = m_numPointsPerElt[eltType];
1095
1096
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 47979429 times.
 47979429 FEL_ASSERT( array.size() == static_cast<std::size_t>(nVpEl) );
1097
2/4
✓ Branch 0 taken 47979429 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 47979429 times.
 47979429 FEL_ASSERT( index >= 0 && index < m_numElements[eltType]);
1098
1099
2/2
✓ Branch 0 taken 181284013 times.
✓ Branch 1 taken 47979429 times.
 229263442 for (int jj = 0; jj < nVpEl; ++jj)
1100 181284013 array[jj] = m_listElements[eltType][ startindex + index*nVpEl + jj ] + flagIncrement;
1101 47979429 }
1102
1103 /***********************************************************************************/
1104 /***********************************************************************************/
1105
1106 203725 void GeometricMeshRegion::getOneElement(const ElementType& eltType, const felInt index, std::vector<felInt>& array,
1107 std::vector<Point*>& rPoints, const felInt startindex, const bool flagIncrement) const
1108 {
1109 203725 getOneElement(eltType, index, array, startindex, flagIncrement);
1110
1111 203725 const int numPointsPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
1112
2/2
✓ Branch 0 taken 804813 times.
✓ Branch 1 taken 203725 times.
 1008538 for (int i_point = 0; i_point < numPointsPerElt; ++i_point)
1113 804813 rPoints[i_point] = const_cast<Point*>(&(m_listPoints[array[i_point]]));
1114 203725 }
1115
1116 /***********************************************************************************/
1117 /***********************************************************************************/
1118
1119 418784 void GeometricMeshRegion::getOneElement(const felInt idElem, std::vector<felInt>& array, const bool flagIncrement) const
1120 {
1121 ElementType eltType;
1122 418784 felInt idByType = -1;
1123
1124
1/2
✓ Branch 1 taken 418784 times.
✗ Branch 2 not taken.
 418784 getTypeElemFromIdElem(idElem, eltType, idByType);
1125
1/2
✓ Branch 1 taken 418784 times.
✗ Branch 2 not taken.
 418784 getOneElement(eltType, idByType, array, 0, flagIncrement);
1126 418784 }
1127
1128 /***********************************************************************************/
1129 /***********************************************************************************/
1130
1131 3459873 void GeometricMeshRegion::setOneElement(const ElementType& eltType, const felInt index,
1132 const std::vector<felInt>& array, const bool flagDecrement)
1133 {
1134 3459873 const int nVpEl = m_numPointsPerElt[eltType];
1135
1136
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 3459873 times.
 3459873 FEL_ASSERT( array.size() == static_cast<std::size_t>(nVpEl) );
1137
2/4
✓ Branch 0 taken 3459873 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3459873 times.
 3459873 FEL_ASSERT( index >= 0 && index < m_numElements[eltType] );
1138
1139
2/2
✓ Branch 0 taken 12964973 times.
✓ Branch 1 taken 3459873 times.
 16424846 for (int jj = 0; jj < nVpEl; jj++)
1140 12964973 m_listElements[eltType][ index*nVpEl + jj ] = array[jj] - flagDecrement;
1141 3459873 }
1142
1143 /***********************************************************************************/
1144 /***********************************************************************************/
1145
1146 ✗ int GeometricMeshRegion::getElementsPerRef(const ElementType& eltType, const int ref, std::vector<felInt>& array) const
1147 {
1148 ✗ if ( m_numElements[eltType] == 0 ) {
1149 ✗ std::cout << "Warning: trying to access a type of element with no elements: " << eltEnumToFelNameGeoEle[eltType].first << "!" << std::endl;
1150 ✗ return 0;
1151 }
1152
1153 ✗ auto it_ref = intRefToBegEndMaps[eltType].find(ref);
1154 ✗ if ( it_ref == intRefToBegEndMaps[eltType].end() ) {
1155 ✗ std::cout << "Warning: trying to access an unexisting reference: " << ref << " for element type "<< eltType << " !" << std::endl;
1156 ✗ return 0;
1157 }
1158
1159 ✗ const int nVpEl = m_numPointsPerElt[eltType];
1160 ✗ const felInt startindex = it_ref->second.first;
1161 ✗ const felInt numRefEle = it_ref->second.second;
1162
1163 ✗ array.resize(nVpEl*numRefEle);
1164
1165 ✗ for (felInt iel = 0 ; iel < numRefEle; ++iel)
1166 ✗ for (int jvert = 0; jvert < nVpEl; ++jvert)
1167 ✗ array[ iel*nVpEl + jvert] = m_listElements[eltType][ startindex + iel*nVpEl + jvert];
1168
1169 ✗ return numRefEle;
1170 }
1171
1172 /***********************************************************************************/
1173 /***********************************************************************************/
1174
1175 // TODO felInt& out_startindex, felInt& out_numRefEle
1176 8 void GeometricMeshRegion::getElementsIDPerRef(const ElementType& eltType, const int ref, felInt* out_startindex, felInt* out_numRefEle) const
1177 {
1178
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
 8 if ( m_numElements[eltType] == 0 )
1179 ✗ std::cout << "Warning: trying to access a type of element with no elements: " << eltEnumToFelNameGeoEle[eltType].first << "!" << std::endl;
1180
1181
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 auto it_ref = intRefToBegEndMaps[eltType].find(ref);
1182
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
 8 if ( it_ref == intRefToBegEndMaps[eltType].end() )
1183 ✗ std::cout << "Warning: trying to access an unexisting reference: " << ref << " for element type "<< eltType << " !" << std::endl;
1184
1185 8 const int nVpEl = m_numPointsPerElt[eltType];
1186 8 const felInt startindex = it_ref->second.first / nVpEl;
1187 8 const felInt numRefEle = it_ref->second.second;
1188
1189 8 *out_startindex = startindex;
1190 8 *out_numRefEle = numRefEle;
1191 8 }
1192
1193 /***********************************************************************************/
1194 /***********************************************************************************/
1195
1196 1398482 void GeometricMeshRegion::getTypeElemFromIdElem(const felInt idElem, ElementType& eltType, felInt& idByType) const
1197 {
1198 felInt numEltPerLabel;
1199 1398482 felInt countElt = 0;
1200
1201 1398482 std::array<const std::vector<ElementType>*,2> bags{ &m_bagElementTypeDomain, &m_bagElementTypeDomainBoundary };
1202
1203
3/6
✓ Branch 2 taken 1398482 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 1398482 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1398482 times.
 1398482 if( idElem >= getNumElementByBag(*bags[0]) + getNumElementByBag(*bags[1]) )
1204 ✗ FEL_ERROR("GeometricMeshRegion: element not found!");
1205
1206
1/2
✓ Branch 1 taken 1398534 times.
✗ Branch 2 not taken.
 1398534 for(std::size_t ibag = 0; ibag < bags.size(); ++ibag) {
1207
3/4
✓ Branch 2 taken 1398534 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1398482 times.
✓ Branch 5 taken 52 times.
 1398534 if( idElem < countElt + getNumElementByBag(*bags[ibag]) ) {
1208
1209
1/2
✓ Branch 2 taken 1398482 times.
✗ Branch 3 not taken.
 1398482 for (std::size_t itype = 0; itype < bags[ibag]->size(); ++itype) {
1210 1398482 eltType = (*bags[ibag])[itype];
1211 1398482 idByType = 0;
1212
1213
1/2
✓ Branch 0 taken 1398482 times.
✗ Branch 1 not taken.
 1398482 if( idElem < countElt + m_numElements[eltType] ) {
1214
1/2
✓ Branch 3 taken 1398534 times.
✗ Branch 4 not taken.
 1398534 for (auto itRef = intRefToBegEndMaps[eltType].begin(); itRef != intRefToBegEndMaps[eltType].end(); ++itRef) {
1215 1398534 numEltPerLabel = itRef->second.second;
1216
1217
2/2
✓ Branch 0 taken 1398482 times.
✓ Branch 1 taken 52 times.
 1398534 if( idElem < countElt + numEltPerLabel ) {
1218 1398482 idByType += idElem - countElt;
1219 1398482 countElt = idElem;
1220
1221 1398482 return;
1222 } else {
1223 52 idByType += numEltPerLabel;
1224 52 countElt += numEltPerLabel;
1225 }
1226 }
1227 } else {
1228 ✗ idByType += m_numElements[eltType];
1229 ✗ countElt += m_numElements[eltType];
1230 }
1231 }
1232 } else {
1233
1/2
✓ Branch 2 taken 52 times.
✗ Branch 3 not taken.
 52 countElt += getNumElementByBag(*bags[ibag]);
1234 }
1235 }
1236 }
1237
1238 /***********************************************************************************/
1239 /***********************************************************************************/
1240
1241 33937583 void GeometricMeshRegion::getIdElemFromTypeElemAndIdByType(const ElementType eltType, const felInt iElByType, felInt& idElt) const
1242 {
1243 ElementType eltTypeTest;
1244
1245 33937583 std::array<const std::vector<ElementType>*,2> bags{ &m_bagElementTypeDomain, &m_bagElementTypeDomainBoundary };
1246
1247 // first, look in the bagDomain for the element type, then look in the bagDomainBoundary for the element type
1248 33937583 idElt = 0;
1249
1/2
✓ Branch 1 taken 35287581 times.
✗ Branch 2 not taken.
 35287581 for(std::size_t ibag = 0; ibag < bags.size(); ++ibag) {
1250
2/2
✓ Branch 2 taken 35289501 times.
✓ Branch 3 taken 1349998 times.
 36639499 for (std::size_t ieltype = 0; ieltype < bags[ibag]->size(); ++ieltype) {
1251 35289501 eltTypeTest = (*bags[ibag])[ieltype];
1252
1253
2/2
✓ Branch 0 taken 33937583 times.
✓ Branch 1 taken 1351918 times.
 35289501 if ( eltTypeTest == eltType ) {
1254 33937583 idElt += iElByType;
1255 33937583 return;
1256 } else {
1257 1351918 idElt += numElements(eltTypeTest);
1258 }
1259 }
1260 }
1261
1262 // if we are here, it means that this element type is not known
1263 ✗ FEL_ERROR("GeometricMeshRegion: We did not find the type " + eltEnumToFelNameGeoEle[eltType].first + " in the mesh!");
1264 }
1265
1266 /***********************************************************************************/
1267 /***********************************************************************************/
1268
1269 193677 void GeometricMeshRegion::getElementFromBdElem(const std::vector<felInt>& facePts, std::vector<felInt>& elemIdPointVol,
1270 std::vector< Point*>& elemPointVol, int& idLocFace, felInt& idElemVol)
1271 {
1272 ElementType eltType;
1273 int nVpEl;
1274
1275
2/3
✓ Branch 1 taken 193437 times.
✓ Branch 2 taken 240 times.
✗ Branch 3 not taken.
 193677 switch ( facePts.size() ) {
1276
1277 193437 case 3:
1278 case 4: {
1279
1280 193437 Face foundFace;
1281
1/2
✓ Branch 1 taken 193437 times.
✗ Branch 2 not taken.
 193437 getFace(facePts, foundFace);
1282
1283
1/4
✗ Branch 2 not taken.
✓ Branch 3 taken 193437 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 193437 FEL_ASSERT( ( foundFace.listNeighVolumes().size() ) == 1 ); // face on boundary!
1284
1285 193437 idElemVol = ( foundFace.listNeighVolumes()[0] )->idVolume();
1286 193437 idLocFace = ( foundFace.listNeighVolumes()[0] )->idLocalFace();
1287
1288 193437 eltType = static_cast<ElementType>( ( foundFace.listNeighVolumes()[0] )->typeVolume() );
1289
1290 193437 nVpEl = m_numPointsPerElt[eltType];
1291
1/2
✓ Branch 1 taken 193437 times.
✗ Branch 2 not taken.
 193437 elemIdPointVol.resize(nVpEl, 0);
1292
1/2
✓ Branch 1 taken 193437 times.
✗ Branch 2 not taken.
 193437 elemPointVol.resize(nVpEl, nullptr);
1293
1/2
✓ Branch 1 taken 193437 times.
✗ Branch 2 not taken.
 193437 getOneElement( eltType, idElemVol, elemIdPointVol, elemPointVol, 0);
1294
1295 193437 break;
1296 193437 }
1297
1298 240 case 2: {
1299
1300 240 Edge foundEdge;
1301
1/2
✓ Branch 1 taken 240 times.
✗ Branch 2 not taken.
 240 getEdge(facePts, foundEdge);
1302
1303
1/4
✗ Branch 2 not taken.
✓ Branch 3 taken 240 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 240 FEL_ASSERT( ( foundEdge.listNeighFacesOfEdges().size() ) == 1 ); // Edge on boundary!
1304
1305 240 idElemVol = ( foundEdge.listNeighFacesOfEdges()[0] )->idFace();
1306 240 idLocFace = ( foundEdge.listNeighFacesOfEdges()[0] )->idLocalEdge();
1307
1308 240 eltType = static_cast<ElementType>( ( foundEdge.listNeighFacesOfEdges()[0] )->typeFace() );
1309
1310 240 nVpEl = m_numPointsPerElt[eltType];
1311
1/2
✓ Branch 1 taken 240 times.
✗ Branch 2 not taken.
 240 elemIdPointVol.resize(nVpEl, 0);
1312
1/2
✓ Branch 1 taken 240 times.
✗ Branch 2 not taken.
 240 elemPointVol.resize(nVpEl, nullptr);
1313
1/2
✓ Branch 1 taken 240 times.
✗ Branch 2 not taken.
 240 getOneElement( eltType, idElemVol, elemIdPointVol, elemPointVol, 0);
1314
1315 240 break;
1316 240 }
1317
1318 ✗ default:
1319 ✗ FEL_ERROR("GeometricMeshRegion: wrong number of poins!");
1320 }
1321 193677 }
1322
1323 /***********************************************************************************/
1324 /***********************************************************************************/
1325
1326 1698793 bool GeometricMeshRegion::getAdjElement(ElementType& eltType, felInt& idElem, int localIndexOfEdgeOrFace) const
1327 {
1328 ElementType eltTypeTmp;
1329 felInt idElemTmp;
1330
1331 1698793 const GeoElement* geoEle = eltEnumToFelNameGeoEle[eltType].second;
1332 1698793 const auto domain_dim = domainDim();
1333
1334 // TODO is really stupid to use two differnt names (i.e. edge and face) for the boundaries of an element...
1335 // for every dimension we must write exactly the same thing twice...
1336
1337
2/2
✓ Branch 0 taken 1698703 times.
✓ Branch 1 taken 90 times.
 1698793 if ( domain_dim == 2 ) {
1338
1339
1/2
✓ Branch 4 taken 1698703 times.
✗ Branch 5 not taken.
 1698703 std::vector<felInt> edges(geoEle->shape().numEdge(),0);
1340
1341
1/4
✗ Branch 2 not taken.
✓ Branch 3 taken 1698703 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 1698703 FEL_ASSERT( localIndexOfEdgeOrFace <= geoEle->shape().numEdge() );
1342
1343
1/2
✓ Branch 1 taken 1698703 times.
✗ Branch 2 not taken.
 1698703 getAllEdgeOfElement(eltType, idElem, edges);
1344
1345
1/2
✓ Branch 5 taken 1698703 times.
✗ Branch 6 not taken.
 1698703 auto neighFaces = m_listEdges.list()[edges[localIndexOfEdgeOrFace]]->listNeighFacesOfEdges();
1346 1698703 std::size_t numNeigh = neighFaces.size();
1347
1348 // case when going out of the domain
1349
2/2
✓ Branch 0 taken 33050 times.
✓ Branch 1 taken 1665653 times.
 1698703 if ( numNeigh < 2 )
1350 33050 return false;
1351
1352
1/2
✓ Branch 0 taken 2498778 times.
✗ Branch 1 not taken.
 2498778 for(std::size_t iNeigh = 0; iNeigh < numNeigh; ++iNeigh) {
1353 2498778 eltTypeTmp = static_cast<ElementType>( neighFaces[iNeigh]->typeFace() );
1354 2498778 idElemTmp = neighFaces[iNeigh]->idFace();
1355
1356 // found new element (not the element itself)
1357
3/4
✓ Branch 0 taken 2498778 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1665653 times.
✓ Branch 3 taken 833125 times.
 2498778 if ( eltTypeTmp != eltType || idElemTmp != idElem ) {
1358 1665653 eltType = eltTypeTmp;
1359 1665653 idElem = idElemTmp;
1360
1361 1665653 return true;
1362 }
1363 }
1364
1365 ✗ return false;
1366
1/2
✓ Branch 2 taken 90 times.
✗ Branch 3 not taken.
 1698793 } else if ( domain_dim == 3 ) {
1367
1368
1/2
✓ Branch 4 taken 90 times.
✗ Branch 5 not taken.
 90 std::vector<felInt> faces(geoEle->shape().numFace(),0);
1369
1370
1/4
✗ Branch 2 not taken.
✓ Branch 3 taken 90 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 90 FEL_ASSERT( localIndexOfEdgeOrFace <= geoEle->shape().numFace() );
1371
1372
1/2
✓ Branch 1 taken 90 times.
✗ Branch 2 not taken.
 90 getAllFaceOfElement(eltType, idElem, faces);
1373
1374 // small std::vector: copy it for legibility...
1375
1/2
✓ Branch 5 taken 90 times.
✗ Branch 6 not taken.
 90 auto neighVols = m_listFaces.list()[faces[localIndexOfEdgeOrFace]]->listNeighVolumes();
1376 90 std::size_t numNeigh = neighVols.size();
1377
1378 // case when going out of the domain
1379
2/2
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 40 times.
 90 if ( numNeigh < 2 )
1380 50 return false;
1381
1382
1/2
✓ Branch 0 taken 60 times.
✗ Branch 1 not taken.
 60 for(std::size_t iNeigh = 0; iNeigh < numNeigh; ++iNeigh) {
1383 60 eltTypeTmp = static_cast<ElementType>( neighVols[iNeigh]->typeVolume() );
1384 60 idElemTmp = neighVols[iNeigh]->idVolume();
1385
1386 // found new element (not the element itself)
1387
3/4
✓ Branch 0 taken 60 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 40 times.
✓ Branch 3 taken 20 times.
 60 if ( eltTypeTmp != eltType || idElemTmp != idElem ) {
1388 40 eltType = eltTypeTmp;
1389 40 idElem = idElemTmp;
1390
1391 40 return true;
1392 }
1393 }
1394
1395 ✗ return false;
1396 90 }
1397
1398 ✗ FEL_ERROR("GeometricMeshRegion: the method should have been already exited at this point");
1399
1400 ✗ return false;
1401 }
1402
1403 /***********************************************************************************/
1404 /***********************************************************************************/
1405
1406 429960 std::pair<bool, GeometricMeshRegion::ElementType> GeometricMeshRegion::getAdjElement(felInt& idElem, int localIndexOfEdgeOrFace) const
1407 {
1408 ElementType eltType;
1409 429960 felInt idByType = -1;
1410
1411
1/2
✓ Branch 1 taken 429960 times.
✗ Branch 2 not taken.
 429960 getTypeElemFromIdElem(idElem, eltType, idByType);
1412
1413
1/2
✓ Branch 1 taken 429960 times.
✗ Branch 2 not taken.
 429960 bool isfound = getAdjElement(eltType, idByType, localIndexOfEdgeOrFace);
1414 429960 idElem = idByType;
1415
1/2
✓ Branch 1 taken 429960 times.
✗ Branch 2 not taken.
 859920 return std::make_pair(isfound, eltType);
1416 }
1417
1418 /***********************************************************************************/
1419 /***********************************************************************************/
1420
1421 ✗ void GeometricMeshRegion::setListNormal(int iNode, Point vecNormal)
1422 {
1423 ✗ m_listNormals[iNode] = vecNormal;
1424 }
1425
1426 /***********************************************************************************/
1427 /***********************************************************************************/
1428
1429 ✗ void GeometricMeshRegion::setListTangent(int idim, int iNode, Point vecTangent)
1430 {
1431 ✗ m_listTangents[idim][iNode] = vecTangent;
1432 }
1433
1434 /***********************************************************************************/
1435 /***********************************************************************************/
1436
1437 4200426 felInt GeometricMeshRegion::getNumElementByBag(const std::vector<ElementType>& bagElem) const
1438 {
1439 4200426 felInt numElem = 0;
1440
2/2
✓ Branch 3 taken 4204298 times.
✓ Branch 4 taken 4200426 times.
 8404724 for (auto it_eltype = bagElem.begin(); it_eltype != bagElem.end(); ++it_eltype)
1441 4204298 numElem += numElements( *it_eltype );
1442
1443 4200426 return numElem;
1444 }
1445
1446 /***********************************************************************************/
1447 /***********************************************************************************/
1448
1449 2489766 void GeometricMeshRegion::getAllEdgeOfElement(const ElementType eltType, const felInt iel, std::vector<felInt>& array) const
1450 {
1451 2489766 const auto& r_edges = eltEnumToFelNameGeoEle[eltType].second;
1452 2489766 const felInt numEdge = r_edges->numEdge();
1453
1454
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 2489766 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 2489766 FEL_ASSERT( array.size() == static_cast<std::size_t>(numEdge) );
1455
1456
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2489766 times.
 2489766 if ( !m_builtEdges )
1457 ✗ FEL_ERROR("GeometricMeshRegion: Edges not build!");
1458
1459 2489766 felInt id_pt0 = 0, id_pt1 = 0;
1460 2489766 felInt pt_edge_beg = 0, pt_edge_end = 0;
1461
1/2
✓ Branch 2 taken 2489766 times.
✗ Branch 3 not taken.
 2489766 std::vector<felInt> the_elem( m_numPointsPerElt[eltType], 0 );
1462
1463 // Get the ids of the points in the iel-th elem
1464
1/2
✓ Branch 1 taken 2489766 times.
✗ Branch 2 not taken.
 2489766 getOneElement(eltType, iel, the_elem, 0);
1465
1466
2/2
✓ Branch 0 taken 7732223 times.
✓ Branch 1 taken 2489766 times.
 10221989 for (felInt iedge = 0; iedge < numEdge; ++iedge) {
1467
1468 // the first point in an edge has the smallest id
1469
1/2
✓ Branch 1 taken 7732223 times.
✗ Branch 2 not taken.
 7732223 id_pt0 = the_elem[ r_edges->pointOfEdge(iedge,0) ];
1470
1/2
✓ Branch 1 taken 7732223 times.
✗ Branch 2 not taken.
 7732223 id_pt1 = the_elem[ r_edges->pointOfEdge(iedge,1) ];
1471
2/2
✓ Branch 0 taken 3913625 times.
✓ Branch 1 taken 3818598 times.
 7732223 if ( id_pt0 < id_pt1 ) {
1472 3913625 pt_edge_beg = id_pt0;
1473 3913625 pt_edge_end = id_pt1;
1474 } else {
1475 3818598 pt_edge_beg = id_pt1;
1476 3818598 pt_edge_end = id_pt0;
1477 }
1478
1479 // find the edge
1480
1/2
✓ Branch 1 taken 7732223 times.
✗ Branch 2 not taken.
 7732223 array[iedge] = m_listEdges.findEdge(pt_edge_beg, pt_edge_end);
1481
1482
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 7732223 times.
 7732223 if( array[iedge] == -1 )
1483 ✗ FEL_ERROR("GeometricMeshRegion: This edge doesn't exist");
1484 }
1485 2489766 }
1486
1487 /***********************************************************************************/
1488 /***********************************************************************************/
1489
1490 548768 void GeometricMeshRegion::getAllEdgeOfElement(const felInt iel, std::vector<felInt>& array) const
1491 {
1492 ElementType eltType;
1493 548768 felInt idByType = -1;
1494
1495
1/2
✓ Branch 1 taken 548768 times.
✗ Branch 2 not taken.
 548768 getTypeElemFromIdElem(iel, eltType, idByType);
1496
1/2
✓ Branch 1 taken 548768 times.
✗ Branch 2 not taken.
 548768 getAllEdgeOfElement(eltType, idByType, array);
1497 548768 }
1498
1499 /***********************************************************************************/
1500 /***********************************************************************************/
1501
1502 110 void GeometricMeshRegion::getAllFaceOfElement(const ElementType eltType, const felInt iel, std::vector<felInt>& array) const
1503 {
1504 110 const auto& r_faces = eltEnumToFelNameGeoEle[eltType].second;
1505 110 const felInt numFace = r_faces->numFace();
1506
1507
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 110 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 110 FEL_ASSERT( array.size() == static_cast<std::size_t>(numFace) );
1508
1509
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 110 times.
 110 if ( !m_builtFaces )
1510 ✗ FEL_ERROR("GeometricMeshRegion: Faces not build!");
1511
1512
1/2
✓ Branch 2 taken 110 times.
✗ Branch 3 not taken.
 110 std::vector<felInt> the_elem( m_numPointsPerElt[eltType], 0);
1513 110 std::vector<felInt> ptOfFace;
1514
1515 // Get the ids of the points in the iel-th elem
1516
1/2
✓ Branch 1 taken 110 times.
✗ Branch 2 not taken.
 110 getOneElement(eltType, iel, the_elem, 0);
1517
2/2
✓ Branch 0 taken 448 times.
✓ Branch 1 taken 110 times.
 558 for (felInt iface = 0; iface < numFace; ++iface) {
1518
2/4
✓ Branch 2 taken 448 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 448 times.
✗ Branch 6 not taken.
 448 ptOfFace.resize( eltEnumToFelNameGeoEle[eltType].second->numPointPerFace( iface ), 0);
1519
1520
2/2
✓ Branch 1 taken 1368 times.
✓ Branch 2 taken 448 times.
 1816 for (std::size_t iptface = 0; iptface < ptOfFace.size(); iptface++ )
1521
1/2
✓ Branch 1 taken 1368 times.
✗ Branch 2 not taken.
 1368 ptOfFace[iptface] = the_elem[ r_faces->pointOfFace(iface, iptface) ];
1522
1523 // look for current face
1524
1/2
✓ Branch 1 taken 448 times.
✗ Branch 2 not taken.
 448 array[iface] = m_listFaces.findFace(ptOfFace);
1525
1526
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 448 times.
 448 if( array[iface] == -1 )
1527 ✗ FEL_ERROR("GeometricMeshRegion: This face doesn't exist");
1528 }
1529 110 }
1530
1531 /***********************************************************************************/
1532 /***********************************************************************************/
1533
1534 ✗ void GeometricMeshRegion::getAllFaceOfElement(const felInt iel, std::vector<felInt>& array) const
1535 {
1536 ElementType eltType;
1537 ✗ felInt idByType = -1;
1538
1539 ✗ getTypeElemFromIdElem(iel, eltType, idByType);
1540 ✗ getAllFaceOfElement(eltType, idByType, array);
1541 }
1542
1543 /***********************************************************************************/
1544 /***********************************************************************************/
1545
1546 2446 void GeometricMeshRegion::getEdge(const std::vector<felInt>& edgePts, Edge& resultEdge) const
1547 {
1548
1/2
✓ Branch 0 taken 2446 times.
✗ Branch 1 not taken.
 2446 if ( m_builtEdges ) {
1549 2446 m_listEdges.searchEdge(edgePts, resultEdge);
1550 } else {
1551 ✗ FEL_ERROR("GeometricMeshRegion: Edges not build!");
1552 }
1553 2446 }
1554
1555 /***********************************************************************************/
1556 /***********************************************************************************/
1557
1558 193437 void GeometricMeshRegion::getFace(const std::vector<felInt>& facePts, Face& resultFace) const
1559 {
1560
1/2
✓ Branch 0 taken 193437 times.
✗ Branch 1 not taken.
 193437 if ( m_builtFaces ) {
1561 193437 m_listFaces.searchFace(facePts, resultFace);
1562 } else {
1563 ✗ FEL_ERROR("GeometricMeshRegion: Faces not build!");
1564 }
1565 193437 }
1566
1567 /***********************************************************************************/
1568 /***********************************************************************************/
1569
1570 391 void GeometricMeshRegion::print(int verbose, std::ostream& outstr) const
1571 {
1572 int rank;
1573
1/2
✓ Branch 1 taken 391 times.
✗ Branch 2 not taken.
 391 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
1574
1575
2/2
✓ Branch 0 taken 333 times.
✓ Branch 1 taken 58 times.
 391 if ( verbose ) {
1576
5/10
✓ Branch 1 taken 333 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 333 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 333 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 333 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 333 times.
✗ Branch 14 not taken.
 333 outstr << "[" << rank << "] Dimension: " << m_domainDim << std::endl;
1577
5/10
✓ Branch 1 taken 333 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 333 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 333 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 333 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 333 times.
✗ Branch 15 not taken.
 333 outstr << "[" << rank << "] Number of points: " << m_listPoints.size() << std::endl;
1578
1579 // print only the num of elements:
1580
6/12
✓ Branch 1 taken 333 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 333 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 333 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 333 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 333 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 333 times.
✗ Branch 17 not taken.
 333 outstr << "\n[" << rank << "] Domain elements: " << getNumDomainElement() << std::endl;
1581
1/2
✓ Branch 1 taken 333 times.
✗ Branch 2 not taken.
 333 printElements(m_bagElementTypeDomain, /*verbose=*/ 1, outstr);
1582
6/12
✓ Branch 1 taken 333 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 333 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 333 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 333 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 333 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 333 times.
✗ Branch 17 not taken.
 333 outstr << "[" << rank << "] Boundary elements: " << getNumBoundaryElement() << std::endl;
1583
1/2
✓ Branch 1 taken 333 times.
✗ Branch 2 not taken.
 333 printElements(m_bagElementTypeDomainBoundary, 1, outstr);
1584
1585
3/6
✓ Branch 1 taken 333 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 333 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 333 times.
✗ Branch 8 not taken.
 333 outstr << "[" << rank << "] Domain Map: ref -> <startPos, numElements>\n";
1586
1/2
✓ Branch 1 taken 333 times.
✗ Branch 2 not taken.
 333 printRefToBegEnd(m_bagElementTypeDomain, verbose, outstr);
1587
3/6
✓ Branch 1 taken 333 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 333 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 333 times.
✗ Branch 8 not taken.
 333 outstr << "[" << rank << "] Boundary Map: ref -> <startPos, numElements>\n";
1588
1/2
✓ Branch 1 taken 333 times.
✗ Branch 2 not taken.
 333 printRefToBegEnd(m_bagElementTypeDomainBoundary, verbose, outstr);
1589 }
1590
1591 // print all vertex/elem info.
1592
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 391 times.
 391 if ( verbose > 28 ) {
1593 ✗ printVertices(outstr);
1594
1595 ✗ outstr << "[" << rank << "] List of domain elements:\n";
1596 ✗ printElements(m_bagElementTypeDomain, verbose, outstr, false);
1597 ✗ outstr << "[" << rank << "] List of boundary elements:\n";
1598 ✗ printElements(m_bagElementTypeDomainBoundary, verbose, outstr, false);
1599 }
1600 391 }
1601
1602 /***********************************************************************************/
1603 /***********************************************************************************/
1604
1605 ✗ void GeometricMeshRegion::printVertices(std::ostream& outstr) const
1606 {
1607 int rank;
1608 ✗ MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
1609
1610 ✗ outstr << "[" << rank << "] List of Points: " << m_listPoints.size() << "\n";
1611 ✗ outstr << "[" << rank << "] ";
1612
1613 ✗ for (auto& r_point : m_listPoints)
1614 ✗ r_point.print(1, outstr);
1615
1616 ✗ outstr << std::endl;
1617 }
1618
1619 /***********************************************************************************/
1620 /***********************************************************************************/
1621
1622 666 void GeometricMeshRegion::printElements(const std::vector<ElementType>& bagElem, const int verbose,
1623 std::ostream& outstr, const bool incrementPointID) const
1624 {
1625 666 int rank = 0;
1626 // MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
1627
1628 666 std::vector<felInt> elem;
1629
2/2
✓ Branch 4 taken 723 times.
✓ Branch 5 taken 666 times.
 1389 for (auto it_eltype = bagElem.begin(); it_eltype != bagElem.end(); ++it_eltype) {
1630 723 ElementType eltType = static_cast<ElementType>(*it_eltype);
1631
1632
7/14
✓ Branch 1 taken 723 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 723 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 723 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 723 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 723 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 723 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 723 times.
✗ Branch 21 not taken.
 723 outstr << "[" << rank << "] " << eltEnumToFelNameGeoEle[eltType].first << ": " << m_numElements[eltType] << std::endl;
1633
1634 // print all elem info.
1635
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 723 times.
 723 if (verbose > 28) {
1636 ✗ int nVpEl = m_numPointsPerElt[eltType];
1637 ✗ elem.resize(nVpEl, 0);
1638
1639 ✗ for ( felInt iel = 0; iel < m_numElements[eltType]; ++iel) {
1640 ✗ getOneElement(eltType, iel, elem, 0, incrementPointID);
1641
1642 ✗ outstr << "[" << rank << "] ";
1643 ✗ for (int iv = 0; iv < nVpEl; ++iv)
1644 ✗ outstr << elem[iv] << " ";
1645 ✗ outstr << "\n";
1646 }
1647 }
1648
1649
1/2
✓ Branch 1 taken 723 times.
✗ Branch 2 not taken.
 723 outstr << std::endl;
1650 }
1651 666 }
1652
1653 /***********************************************************************************/
1654 /***********************************************************************************/
1655
1656 666 void GeometricMeshRegion::printRefToBegEnd(const std::vector<ElementType>& bagElem, const int verbose, std::ostream& outstr) const
1657 {
1658 int rank;
1659
1/2
✓ Branch 1 taken 666 times.
✗ Branch 2 not taken.
 666 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
1660
1661 666 std::vector<felInt> elem;
1662
2/2
✓ Branch 4 taken 723 times.
✓ Branch 5 taken 666 times.
 1389 for (auto it_eltype = bagElem.begin(); it_eltype != bagElem.end(); ++it_eltype) {
1663 723 ElementType eltType = static_cast<ElementType>(*it_eltype);
1664
1665
7/14
✓ Branch 1 taken 723 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 723 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 723 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 723 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 723 times.
✗ Branch 15 not taken.
✓ Branch 18 taken 723 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 723 times.
✗ Branch 22 not taken.
 723 outstr << "[" << rank << "] " << eltEnumToFelNameGeoEle[eltType].first << ": " << intRefToBegEndMaps[eltType].size() << std::endl;
1666
1667 // print all info.
1668
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 713 times.
 723 if (verbose > 10) {
1669 10 int nVpEl = m_numPointsPerElt[eltType];
1670
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 elem.resize(nVpEl, 0);
1671
1672
2/2
✓ Branch 4 taken 16 times.
✓ Branch 5 taken 10 times.
 26 for (auto it_ref = intRefToBegEndMaps[eltType].begin(); it_ref != intRefToBegEndMaps[eltType].end(); ++it_ref)
1673
5/10
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 16 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 16 times.
✗ Branch 15 not taken.
 16 outstr << "[" << rank << "] ref = " << it_ref->first << " \t-->\t < startPos = "
1674
4/8
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 16 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 16 times.
✗ Branch 13 not taken.
 16 << it_ref->second.first << " ; numElements = " << it_ref->second.second << " >\n";
1675 }
1676
1677
1/2
✓ Branch 1 taken 723 times.
✗ Branch 2 not taken.
 723 outstr << std::endl;
1678 }
1679 666 }
1680
1681 /***********************************************************************************/
1682 /***********************************************************************************/
1683
1684 1096 void GeometricMeshRegion::m_allocateListElementsByRef(const std::vector<ElementType>& theElementBag, GeometricMeshRegion& meshGlobal,
1685 const std::vector<int>& eltPartition, const int rank, felInt& countEltTot)
1686 {
1687
2/2
✓ Branch 4 taken 1165 times.
✓ Branch 5 taken 1096 times.
 2261 for (auto it_eltype = theElementBag.begin(); it_eltype != theElementBag.end(); ++it_eltype) {
1688 1165 ElementType eltType = static_cast<ElementType>(*it_eltype);
1689
1690 // Order the elements by reference
1691
2/2
✓ Branch 4 taken 2675 times.
✓ Branch 5 taken 1165 times.
 3840 for (auto itRef = meshGlobal.intRefToBegEndMaps[eltType].begin(); itRef != meshGlobal.intRefToBegEndMaps[eltType].end(); ++itRef) {
1692
1693 2675 felInt numElemPerRefGlob = itRef->second.second;
1694
1695
2/2
✓ Branch 0 taken 3704477 times.
✓ Branch 1 taken 2675 times.
 3707152 for (felInt iel = 0; iel < numElemPerRefGlob; iel++) {
1696
2/2
✓ Branch 1 taken 1039942 times.
✓ Branch 2 taken 2664535 times.
 3704477 if (eltPartition[ countEltTot ] == rank )
1697 1039942 m_numElements[eltType]++;
1698
1699 3704477 countEltTot++;
1700 }
1701 }
1702
1703
1/2
✓ Branch 1 taken 1165 times.
✗ Branch 2 not taken.
 1165 allocateElements(eltType);
1704 }
1705 1096 }
1706
1707 /***********************************************************************************/
1708 /***********************************************************************************/
1709
1710 1096 void GeometricMeshRegion::m_fillListElementsByRef(const std::vector<ElementType>& theElementBag, GeometricMeshRegion& meshGlobal,
1711 const std::vector<int>& eltPartition, const int rank, std::vector<felInt>& loc2globElem, felInt& countEltTot)
1712 {
1713 1096 std::vector<felInt> elem;
1714
1715 // for each element type in the bag
1716
2/2
✓ Branch 4 taken 1165 times.
✓ Branch 5 taken 1096 times.
 2261 for (auto it_eltype = theElementBag.begin(); it_eltype != theElementBag.end(); ++it_eltype) {
1717 1165 ElementType eltType = static_cast<ElementType>(*it_eltype);
1718 1165 felInt neTotLoc = 0;
1719 1165 int nVpEl = m_numPointsPerElt[eltType];
1720
1721
1/2
✓ Branch 1 taken 1165 times.
✗ Branch 2 not taken.
 1165 elem.resize(nVpEl, 0);
1722
1723 // select the elements in the partition (ordered by reference)
1724 // for each label in the global mesh for the current element type
1725
2/2
✓ Branch 4 taken 2675 times.
✓ Branch 5 taken 1165 times.
 3840 for (auto itRef = meshGlobal.intRefToBegEndMaps[eltType].begin(); itRef != meshGlobal.intRefToBegEndMaps[eltType].end(); ++itRef) {
1726 2675 const int theRef = itRef->first;
1727 2675 const felInt numElemPerRefGlobal = itRef->second.second;
1728
1729 2675 felInt numElemPerRefLocal = 0;
1730 2675 const felInt posRefElemLocal = neTotLoc*nVpEl;
1731 2675 const felInt startindex = itRef->second.first;
1732
1733 // for each element with the current label
1734
2/2
✓ Branch 0 taken 3704477 times.
✓ Branch 1 taken 2675 times.
 3707152 for (felInt iel = 0; iel < numElemPerRefGlobal; iel++) {
1735
2/2
✓ Branch 1 taken 1039942 times.
✓ Branch 2 taken 2664535 times.
 3704477 if ( eltPartition[ countEltTot ] == rank ) {
1736 // fill the local to global mapping
1737
1/2
✓ Branch 1 taken 1039942 times.
✗ Branch 2 not taken.
 1039942 loc2globElem.push_back(countEltTot);
1738
1739 // get the coordinate of the cooresponding global element
1740
1/2
✓ Branch 1 taken 1039942 times.
✗ Branch 2 not taken.
 1039942 meshGlobal.getOneElement(eltType, iel, elem, startindex);
1741
1742 // std::set the local element wi
1743
1/2
✓ Branch 1 taken 1039942 times.
✗ Branch 2 not taken.
 1039942 setOneElement(eltType, neTotLoc, elem, false);
1744
1745 // increment the number of elements
1746 1039942 neTotLoc++;
1747 1039942 numElemPerRefLocal++;
1748 }
1749 3704477 countEltTot++;
1750 }
1751
1752 // Update intRefToBegEndMapsb
1753
2/2
✓ Branch 0 taken 1716 times.
✓ Branch 1 taken 959 times.
 2675 if ( numElemPerRefLocal != 0 ) {
1754
2/4
✓ Branch 1 taken 1716 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1716 times.
✗ Branch 5 not taken.
 1716 intRefToBegEndMaps[eltType][theRef] = std::make_pair(posRefElemLocal, numElemPerRefLocal);
1755
2/4
✓ Branch 1 taken 1716 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1716 times.
✗ Branch 5 not taken.
 1716 m_intRefToEnum[ theRef ].insert( eltType );
1756 }
1757 }
1758 }
1759 1096 }
1760
1761 /***********************************************************************************/
1762 /***********************************************************************************/
1763
1764 279 void GeometricMeshRegion::m_buildEdgesPerBag(const std::vector<ElementType>& theBagElt, int theDim)
1765 {
1766 279 felInt pt_edge_beg = 0, pt_edge_end = 0;
1767 279 felInt id_pt0 = 0, id_pt1 = 0;
1768
1769
2/2
✓ Branch 4 taken 2046 times.
✓ Branch 5 taken 279 times.
 2325 for (auto it_elttype = theBagElt.begin(); it_elttype != theBagElt.end(); ++it_elttype) {
1770 2046 ElementType eltType = (ElementType)*it_elttype;
1771
1/2
✓ Branch 2 taken 2046 times.
✗ Branch 3 not taken.
 2046 std::vector<felInt> the_elem( m_numPointsPerElt[eltType] , m_listEdges.m_null_int );
1772
1773
2/2
✓ Branch 0 taken 145079 times.
✓ Branch 1 taken 2046 times.
 147125 for ( felInt iel = 0; iel < m_numElements[eltType]; iel++ ) {
1774
1/2
✓ Branch 1 taken 145079 times.
✗ Branch 2 not taken.
 145079 getOneElement(eltType, iel, the_elem, 0);
1775
1776
2/2
✓ Branch 2 taken 561748 times.
✓ Branch 3 taken 145079 times.
 706827 for (felInt iedge = 0; iedge < eltEnumToFelNameGeoEle[eltType].second->numEdge(); iedge ++ ) {
1777 561748 NeighFacesOfEdges::Pointer the_FaceNeighbour = nullptr;
1778 561748 NeighVolumesOfEdges::Pointer the_VolumeNeighbour = nullptr;
1779
1780
3/4
✓ Branch 0 taken 6428 times.
✓ Branch 1 taken 281720 times.
✓ Branch 2 taken 273600 times.
✗ Branch 3 not taken.
 561748 switch (theDim) {
1781
1782 6428 case 1:
1783 6428 break;
1784
1785 281720 case 2:
1786
1/2
✓ Branch 1 taken 281720 times.
✗ Branch 2 not taken.
 281720 the_FaceNeighbour = felisce::make_shared<NeighFacesOfEdges>();
1787 281720 the_FaceNeighbour->idFace() = iel;
1788 281720 the_FaceNeighbour->typeFace() = eltType;
1789 281720 the_FaceNeighbour->idLocalEdge() = iedge;
1790 281720 break;
1791
1792 273600 case 3:
1793
1/2
✓ Branch 1 taken 273600 times.
✗ Branch 2 not taken.
 273600 the_VolumeNeighbour = felisce::make_shared<NeighVolumesOfEdges>();
1794 273600 the_VolumeNeighbour->idVolume() = iel;
1795 273600 the_VolumeNeighbour->typeVolume() = eltType;
1796 273600 the_VolumeNeighbour->idLocalEdge() = iedge;
1797 273600 break;
1798
1799 ✗ default:
1800 ✗ FEL_ERROR( "ERROR: dimension for edge building should be 1, 2 or 3. Exiting." );
1801 }
1802
1803
1/2
✓ Branch 2 taken 561748 times.
✗ Branch 3 not taken.
 561748 id_pt0 = the_elem[ eltEnumToFelNameGeoEle[eltType].second->pointOfEdge(iedge, 0) ];
1804
1/2
✓ Branch 2 taken 561748 times.
✗ Branch 3 not taken.
 561748 id_pt1 = the_elem[ eltEnumToFelNameGeoEle[eltType].second->pointOfEdge(iedge, 1) ];
1805
1806
2/2
✓ Branch 0 taken 278110 times.
✓ Branch 1 taken 283638 times.
 561748 if ( id_pt0 < id_pt1 ) {
1807 278110 pt_edge_beg = id_pt0;
1808 278110 pt_edge_end = id_pt1;
1809 } else {
1810 283638 pt_edge_beg = id_pt1;
1811 283638 pt_edge_end = id_pt0;
1812 }
1813
1814
1/2
✓ Branch 3 taken 561748 times.
✗ Branch 4 not taken.
 561748 m_listEdges.searchAddEdge(pt_edge_beg, pt_edge_end,
1815 the_FaceNeighbour, the_VolumeNeighbour);
1816 561748 }
1817 }
1818 2046 }
1819 279 }
1820
1821 /***********************************************************************************/
1822 /***********************************************************************************/
1823
1824 266 void GeometricMeshRegion::m_buildFacesPerBag(const std::vector<ElementType>& theBagElt, int theDim)
1825 {
1826 266 bool check_surface_mesh_flag = false;
1827
1828
2/2
✓ Branch 4 taken 2527 times.
✓ Branch 5 taken 266 times.
 2793 for (auto it_elttype = theBagElt.begin(); it_elttype != theBagElt.end(); ++it_elttype) {
1829 2527 ElementType eltType = static_cast<ElementType>(*it_elttype);
1830
1/2
✓ Branch 2 taken 2527 times.
✗ Branch 3 not taken.
 2527 std::vector <felInt> the_elem( m_numPointsPerElt[eltType] , m_listFaces.m_null_int );
1831 2527 std::vector <felInt> ptOfFace;
1832
1833
2/2
✓ Branch 0 taken 881010 times.
✓ Branch 1 taken 2527 times.
 883537 for ( felInt iel = 0; iel < m_numElements[eltType]; iel++ ) {
1834
1/2
✓ Branch 1 taken 881010 times.
✗ Branch 2 not taken.
 881010 getOneElement(eltType, iel, the_elem, 0);
1835
1836
2/2
✓ Branch 2 taken 3025590 times.
✓ Branch 3 taken 881010 times.
 3906600 for (felInt iface = 0; iface < eltEnumToFelNameGeoEle[eltType].second->numFace(); iface ++) {
1837 3025590 NeighVolumesOfFaces* the_VolumeNeighbour = nullptr;
1838
1839
2/3
✓ Branch 0 taken 166646 times.
✓ Branch 1 taken 2858944 times.
✗ Branch 2 not taken.
 3025590 switch (theDim) {
1840
1841 166646 case 2:
1842 166646 break;
1843
1844 2858944 case 3:
1845
1/2
✓ Branch 1 taken 2858944 times.
✗ Branch 2 not taken.
 2858944 the_VolumeNeighbour = new NeighVolumesOfFaces;
1846 2858944 the_VolumeNeighbour->idVolume() = iel;
1847 2858944 the_VolumeNeighbour->typeVolume() = eltType;
1848 2858944 the_VolumeNeighbour->idLocalFace() = iface;
1849 2858944 break;
1850
1851 ✗ default:
1852 ✗ FEL_ERROR( "ERROR: dimension for face building should be 2 or 3. Exiting." );
1853 }
1854
1855 //! a priori: for each element, a face can switch to quadrangle or triangle (think of prisms).
1856
1/2
✓ Branch 2 taken 3025590 times.
✗ Branch 3 not taken.
 3025590 int numPointPerFace = eltEnumToFelNameGeoEle[eltType].second->numPointPerFace( iface ); //std::cout <<" test face " << std::endl;
1857
1/2
✓ Branch 1 taken 3025590 times.
✗ Branch 2 not taken.
 3025590 ptOfFace.resize( numPointPerFace, 0);
1858
1859 // global point numbering of the local face points:
1860
2/2
✓ Branch 0 taken 9083592 times.
✓ Branch 1 taken 3025590 times.
 12109182 for (int iptface = 0; iptface < numPointPerFace; iptface ++) {
1861
1/2
✓ Branch 2 taken 9083592 times.
✗ Branch 3 not taken.
 9083592 ptOfFace[iptface] = the_elem[ eltEnumToFelNameGeoEle[eltType].second->pointOfFace(iface, iptface) ];
1862 }
1863
2/4
✓ Branch 1 taken 3025590 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3025590 times.
✗ Branch 5 not taken.
 3025590 m_listFaces.searchAddFace(ptOfFace, the_VolumeNeighbour, check_surface_mesh_flag);
1864 }
1865 }
1866 2527 }
1867 266 }
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900 /***********************************************************************************/
1901 /***********************************************************************************/
1902
1903 ✗ void GeometricMeshRegion::extractVerticesFromBoundary(const std::vector<int>& label, std::set<felInt>& list)
1904 {
1905 int currentLabel;
1906 ✗ std::vector<felInt> elmList;
1907
1908 ✗ for (std::size_t iEltType = 0; iEltType < bagElementTypeDomainBoundary().size(); ++iEltType) {
1909 ✗ ElementType eltType = bagElementTypeDomainBoundary()[iEltType];
1910
1911 ✗ auto& listRefPerType = intRefToBegEndMaps[eltType];
1912 ✗ for (auto itRef = listRefPerType.begin(); itRef != listRefPerType.end(); ++itRef) {
1913 ✗ currentLabel = itRef->first;
1914
1915 // Skip boundaries not included in label
1916 ✗ if ( std::find(label.begin(), label.end(), currentLabel) == label.end() )
1917 ✗ continue;
1918
1919 // Get elements and fill the set with points
1920 ✗ getElementsPerRef(eltType, currentLabel, elmList);
1921
1922 // Store ids
1923 ✗ list.insert(elmList.begin(), elmList.end());
1924 }
1925 }
1926 }
1927
1928 /***********************************************************************************/
1929 /***********************************************************************************/
1930 // TODO this function should be split into two different functions
1931 // 1) compute local coordinates
1932 // 2) isPointInElement
1933 3179 bool GeometricMeshRegion::findLocalCoord(ElementType& ityp, felInt& iel, const Point& point, Point& pointLoc)
1934 {
1935 3179 const GeoElement* geoEl = eltEnumToFelNameGeoEle[ityp].second;
1936 // const int nbrVer = m_numPointsPerElt[ityp];
1937 3179 const int nbrVer = geoEl->numPoint();
1938
1939
1/2
✓ Branch 2 taken 3179 times.
✗ Branch 3 not taken.
 3179 std::vector<felInt> veridxOfElt( nbrVer, 0);
1940
1/2
✓ Branch 2 taken 3179 times.
✗ Branch 3 not taken.
 3179 std::vector<Point*> pointsOfElt( nbrVer, nullptr);
1941
1942
1/2
✓ Branch 1 taken 3179 times.
✗ Branch 2 not taken.
 3179 getOneElement(ityp, iel, veridxOfElt, pointsOfElt, 0);
1943
1944 // Point localPt(0.0);
1945 3179 pointLoc.clear();
1946
1947
1/2
✓ Branch 1 taken 3179 times.
✗ Branch 2 not taken.
 3179 UBlasVector RHS(m_numCoor);
1948
1/2
✓ Branch 1 taken 3179 times.
✗ Branch 2 not taken.
 3179 UBlasVector X(m_numCoor);
1949
1/2
✓ Branch 1 taken 3179 times.
✗ Branch 2 not taken.
 3179 UBlasMatrix Jac(m_numCoor, m_numCoor);
1950
1/2
✓ Branch 1 taken 3179 times.
✗ Branch 2 not taken.
 3179 UBlasMatrix invJac(m_numCoor, m_numCoor);
1951 double det;
1952
1953 // Determines the local coordinate
1954 3179 std::size_t iter = 0;
1955 3179 double resid = NWTOLL + 1.;
1956
3/4
✓ Branch 0 taken 6347 times.
✓ Branch 1 taken 3179 times.
✓ Branch 2 taken 6347 times.
✗ Branch 3 not taken.
 9526 while ( resid > NWTOLL && iter < NMAX ) {
1957
1958 // Assemble matrix and RHS
1959
1/2
✓ Branch 1 taken 6347 times.
✗ Branch 2 not taken.
 6347 Jac.clear();
1960
2/2
✓ Branch 0 taken 12729 times.
✓ Branch 1 taken 6347 times.
 19076 for (int iDim = 0; iDim < m_numCoor; ++iDim) {
1961
1/2
✓ Branch 2 taken 12729 times.
✗ Branch 3 not taken.
 12729 RHS(iDim) = - point.coor(iDim);
1962
1963
2/2
✓ Branch 0 taken 38394 times.
✓ Branch 1 taken 12729 times.
 51123 for(int iPoint = 0; iPoint < nbrVer; ++iPoint) {
1964
2/4
✓ Branch 2 taken 38394 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 38394 times.
✗ Branch 8 not taken.
 38394 RHS(iDim) += geoEl->basisFunction().phi( iPoint, pointLoc) * pointsOfElt[iPoint]->coor(iDim);
1965
1966
2/2
✓ Branch 0 taken 77292 times.
✓ Branch 1 taken 38394 times.
 115686 for (int jDim = 0; jDim < m_numCoor; ++jDim)
1967
2/4
✓ Branch 2 taken 77292 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 77292 times.
✗ Branch 8 not taken.
 77292 Jac(iDim,jDim) -= geoEl->basisFunction().dPhi( iPoint, jDim, pointLoc) * pointsOfElt[iPoint]->coor(iDim);
1968 }
1969 }
1970
1971
1/2
✓ Branch 1 taken 6347 times.
✗ Branch 2 not taken.
 6347 MathUtilities::InvertMatrix(Jac, invJac, det);
1972
1973
3/6
✓ Branch 1 taken 6347 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6347 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6347 times.
✗ Branch 8 not taken.
 6347 noalias(X) = prod(invJac, RHS);
1974
1975 6347 resid = 0.0;
1976
2/2
✓ Branch 0 taken 12729 times.
✓ Branch 1 taken 6347 times.
 19076 for(int iCoor = 0; iCoor < m_numCoor; ++iCoor) {
1977
1/2
✓ Branch 1 taken 12729 times.
✗ Branch 2 not taken.
 12729 resid += std::pow( X[iCoor], 2);
1978
1/2
✓ Branch 1 taken 12729 times.
✗ Branch 2 not taken.
 12729 pointLoc.coor(iCoor) += X[iCoor];
1979 }
1980
1981 6347 iter++;
1982 }
1983
1984
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3179 times.
 3179 if ( iter >= NMAX ) {
1985 ✗ fprintf(stderr,"WARNING: maximum number of iterations %d reached for element %d of type...",static_cast<int>(NMAX),iel);
1986 }
1987
1988 /*
1989 To access at the edge defined by the 2 vertices (i,j):
1990 - go to the i-th index in the list of pointer
1991 - browse the chained list of edges where 1st vertex is i
1992 */
1993
1994 // Check if point is interior or not
1995 felInt tmpNeigh;
1996 ElementType tmptyp;
1997
1998 3179 TypeShape elShape = static_cast<TypeShape>( geoEl->shape().typeShape() );
1999
1/2
✓ Branch 2 taken 3179 times.
✗ Branch 3 not taken.
 3179 std::vector<felInt> faces(0,0);
2000
1/2
✓ Branch 2 taken 3179 times.
✗ Branch 3 not taken.
 3179 std::vector<felInt> edges(0,0);
2001
2002
2/2
✓ Branch 0 taken 3159 times.
✓ Branch 1 taken 20 times.
 3179 if ( m_numCoor == 2 ) {
2003
1/2
✓ Branch 3 taken 3159 times.
✗ Branch 4 not taken.
 3159 edges.resize( geoEl->shape().numEdge(), 0);
2004
1/2
✓ Branch 1 taken 3159 times.
✗ Branch 2 not taken.
 3159 getAllEdgeOfElement(ityp, iel, edges);
2005 }
2006
1/2
✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
 20 else if ( m_numCoor == 3 ) {
2007
1/2
✓ Branch 3 taken 20 times.
✗ Branch 4 not taken.
 20 faces.resize( geoEl->shape().numFace(), 0);
2008
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 getAllFaceOfElement(ityp, iel, faces);
2009 }
2010
2011
4/7
✗ Branch 0 not taken.
✓ Branch 1 taken 3154 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 16 times.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 3179 switch (elShape) {
2012
2013 ✗ case Segment:
2014
2015 ✗ if ( pointLoc[0] < - 1 - GeoUtilities::GEO_TOL )
2016 ✗ return false;
2017
2018 ✗ if ( pointLoc[0] > + 1 + GeoUtilities::GEO_TOL )
2019 ✗ return false;
2020
2021 ✗ return true;
2022
2023 break;
2024
2025 3154 case Triangle:
2026
2027 // edge 0 ->adjacency
2028
2/2
✓ Branch 1 taken 1010 times.
✓ Branch 2 taken 2144 times.
 3154 if ( pointLoc[1] < - GeoUtilities::GEO_TOL ) {
2029
2030 1010 std::size_t numNeigh = m_listEdges.list()[edges[0]]->listNeighFacesOfEdges().size();
2031
1/2
✓ Branch 0 taken 1010 times.
✗ Branch 1 not taken.
 1010 for (std::size_t i = 0; i < numNeigh; ++i) {
2032 1010 tmptyp = static_cast<ElementType>(m_listEdges.list()[edges[0]]->listNeighFacesOfEdges()[i]->typeFace());
2033 1010 tmpNeigh = m_listEdges.list()[edges[0]]->listNeighFacesOfEdges()[i]->idFace();
2034
2035
2/4
✓ Branch 0 taken 1010 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1010 times.
✗ Branch 3 not taken.
 1010 if ( tmptyp != ityp || tmpNeigh != iel ) {
2036 1010 ityp = tmptyp;
2037 1010 iel = tmpNeigh;
2038 1010 break;
2039 }
2040 }
2041
2042 1010 return false;
2043 }
2044
2045 // edge 2 ->adjacency
2046
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2144 times.
 2144 if ( pointLoc[0] < - GeoUtilities::GEO_TOL ) {
2047
2048 ✗ std::size_t numNeigh = m_listEdges.list()[edges[2]]->listNeighFacesOfEdges().size();
2049 ✗ for (std::size_t i = 0; i < numNeigh; ++i) {
2050 ✗ tmptyp = static_cast<ElementType>(m_listEdges.list()[edges[2]]->listNeighFacesOfEdges()[i]->typeFace());
2051 ✗ tmpNeigh = m_listEdges.list()[edges[2]]->listNeighFacesOfEdges()[i]->idFace();
2052
2053 ✗ if ( tmptyp != ityp || tmpNeigh != iel ) {
2054 ✗ ityp = tmptyp;
2055 ✗ iel = tmpNeigh;
2056 ✗ break;
2057 }
2058 }
2059
2060 ✗ return false;
2061 }
2062
2063 // edge 1 ->adjacency
2064
2/2
✓ Branch 2 taken 380 times.
✓ Branch 3 taken 1764 times.
 2144 if ( 1 - pointLoc[0] - pointLoc[1] < - GeoUtilities::GEO_TOL ) {
2065
2066 380 std::size_t numNeigh = m_listEdges.list()[edges[1]]->listNeighFacesOfEdges().size();
2067
1/2
✓ Branch 0 taken 520 times.
✗ Branch 1 not taken.
 520 for (std::size_t i = 0; i < numNeigh; ++i) {
2068 520 tmptyp = static_cast<ElementType>(m_listEdges.list()[edges[1]]->listNeighFacesOfEdges()[i]->typeFace());
2069 520 tmpNeigh = m_listEdges.list()[edges[1]]->listNeighFacesOfEdges()[i]->idFace();
2070
2071
3/4
✓ Branch 0 taken 520 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 380 times.
✓ Branch 3 taken 140 times.
 520 if ( tmptyp != ityp || tmpNeigh != iel ) {
2072 380 ityp = tmptyp;
2073 380 iel = tmpNeigh;
2074 380 break;
2075 }
2076 }
2077
2078 380 return false;
2079 }
2080
2081 1764 return true;
2082
2083 break;
2084
2085 5 case Quadrilateral:
2086
2087 // edge 3 ->adjacency
2088
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
 5 if ( pointLoc[0] < - 1 - GeoUtilities::GEO_TOL ) {
2089
2090 ✗ std::size_t numNeigh = m_listEdges.list()[edges[3]]->listNeighFacesOfEdges().size();
2091 ✗ for (std::size_t i = 0; i < numNeigh; ++i) {
2092 ✗ tmptyp = static_cast<ElementType>(m_listEdges.list()[edges[3]]->listNeighFacesOfEdges()[i]->typeFace());
2093 ✗ tmpNeigh = m_listEdges.list()[edges[3]]->listNeighFacesOfEdges()[i]->idFace();
2094
2095 ✗ if ( tmptyp != ityp || tmpNeigh != iel ) {
2096 ✗ ityp = tmptyp;
2097 ✗ iel = tmpNeigh;
2098 ✗ break;
2099 }
2100 }
2101
2102 ✗ return false;
2103 }
2104
2105 // edge 1 ->adjacency
2106
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
 5 if ( pointLoc[0] > + 1 + GeoUtilities::GEO_TOL ) {
2107
2108 ✗ std::size_t numNeigh = m_listEdges.list()[edges[1]]->listNeighFacesOfEdges().size();
2109 ✗ for (std::size_t i=0; i<numNeigh; ++i) {
2110 ✗ tmptyp = static_cast<ElementType>(m_listEdges.list()[edges[1]]->listNeighFacesOfEdges()[i]->typeFace());
2111 ✗ tmpNeigh = m_listEdges.list()[edges[1]]->listNeighFacesOfEdges()[i]->idFace();
2112
2113 ✗ if ( tmptyp != ityp || tmpNeigh != iel ) {
2114 ✗ ityp = tmptyp;
2115 ✗ iel = tmpNeigh;
2116 ✗ break;
2117 }
2118 }
2119
2120 ✗ return false;
2121 }
2122
2123 // edge 0 ->adjacency
2124
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
 5 if ( pointLoc[1] < - 1 - GeoUtilities::GEO_TOL ) {
2125
2126 ✗ std::size_t numNeigh = m_listEdges.list()[edges[0]]->listNeighFacesOfEdges().size();
2127 ✗ for (std::size_t i = 0; i < numNeigh; ++i) {
2128 ✗ tmptyp = static_cast<ElementType>(m_listEdges.list()[edges[0]]->listNeighFacesOfEdges()[i]->typeFace());
2129 ✗ tmpNeigh = m_listEdges.list()[edges[0]]->listNeighFacesOfEdges()[i]->idFace();
2130
2131 ✗ if ( tmptyp != ityp || tmpNeigh != iel ) {
2132 ✗ ityp = tmptyp;
2133 ✗ iel = tmpNeigh;
2134 ✗ break;
2135 }
2136 }
2137
2138 ✗ return false;
2139 }
2140
2141 // edge 2 ->adjacency
2142
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
 5 if ( pointLoc[1] > + 1 + GeoUtilities::GEO_TOL ) {
2143
2144 ✗ std::size_t numNeigh = m_listEdges.list()[edges[2]]->listNeighFacesOfEdges().size();
2145 ✗ for (std::size_t i = 0; i < numNeigh; ++i) {
2146 ✗ tmptyp = static_cast<ElementType>(m_listEdges.list()[edges[2]]->listNeighFacesOfEdges()[i]->typeFace());
2147 ✗ tmpNeigh = m_listEdges.list()[edges[2]]->listNeighFacesOfEdges()[i]->idFace();
2148
2149 ✗ if ( tmptyp != ityp || tmpNeigh != iel ) {
2150 ✗ ityp = tmptyp;
2151 ✗ iel = tmpNeigh;
2152 ✗ break;
2153 }
2154 }
2155
2156 ✗ return false;
2157 }
2158
2159 5 return true;
2160
2161 break;
2162
2163 16 case Tetrahedron:
2164
2165 // face 0 ->adjacency
2166
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 16 times.
 16 if ( pointLoc[2] < - GeoUtilities::GEO_TOL ) {
2167
2168 ✗ std::size_t numNeigh = m_listFaces.list()[faces[0]]->listNeighVolumes().size();
2169 ✗ for (std::size_t i = 0; i < numNeigh; ++i) {
2170 ✗ tmptyp = static_cast<ElementType>(m_listFaces.list()[faces[0]]->listNeighVolumes()[i]->typeVolume());
2171 ✗ tmpNeigh = m_listFaces.list()[faces[0]]->listNeighVolumes()[i]->idVolume();
2172
2173 ✗ if ( tmptyp != ityp || tmpNeigh != iel ) {
2174 ✗ ityp = tmptyp;
2175 ✗ iel = tmpNeigh;
2176 ✗ break;
2177 }
2178 }
2179
2180 ✗ return false;
2181 }
2182
2183 // face 1 ->adjacency
2184
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 16 times.
 16 if ( pointLoc[1] < - GeoUtilities::GEO_TOL ) {
2185
2186 ✗ std::size_t numNeigh = m_listFaces.list()[faces[1]]->listNeighVolumes().size();
2187 ✗ for (std::size_t i = 0; i < numNeigh; ++i) {
2188 ✗ tmptyp = static_cast<ElementType>(m_listFaces.list()[faces[1]]->listNeighVolumes()[i]->typeVolume());
2189 ✗ tmpNeigh = m_listFaces.list()[faces[1]]->listNeighVolumes()[i]->idVolume();
2190
2191 ✗ if ( tmptyp != ityp || tmpNeigh != iel ) {
2192 ✗ ityp = tmptyp;
2193 ✗ iel = tmpNeigh;
2194 ✗ break;
2195 }
2196 }
2197
2198 ✗ return false;
2199 }
2200
2201 // face 3 ->adjacency
2202
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 16 times.
 16 if ( pointLoc[0] < - GeoUtilities::GEO_TOL ) {
2203
2204 ✗ std::size_t numNeigh = m_listFaces.list()[faces[3]]->listNeighVolumes().size();
2205 ✗ for (std::size_t i = 0; i < numNeigh; ++i) {
2206 ✗ tmptyp = static_cast<ElementType>(m_listFaces.list()[faces[3]]->listNeighVolumes()[i]->typeVolume());
2207 ✗ tmpNeigh = m_listFaces.list()[faces[3]]->listNeighVolumes()[i]->idVolume();
2208
2209 ✗ if ( tmptyp != ityp || tmpNeigh != iel ) {
2210 ✗ ityp = tmptyp;
2211 ✗ iel = tmpNeigh;
2212 ✗ break;
2213 }
2214 }
2215
2216 ✗ return false;
2217 }
2218
2219 // face 2 ->adjacency
2220
2/2
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 15 times.
 16 if ( 1 - pointLoc[0] - pointLoc[1] - pointLoc[2] < - GeoUtilities::GEO_TOL ) {
2221
2222 1 std::size_t numNeigh = m_listFaces.list()[faces[2]]->listNeighVolumes().size();
2223
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 for (std::size_t i = 0; i < numNeigh; ++i) {
2224 2 tmptyp = static_cast<ElementType>(m_listFaces.list()[faces[2]]->listNeighVolumes()[i]->typeVolume());
2225 2 tmpNeigh = m_listFaces.list()[faces[2]]->listNeighVolumes()[i]->idVolume();
2226
2227
3/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 1 times.
 2 if ( tmptyp != ityp || tmpNeigh != iel ) {
2228 1 ityp = tmptyp;
2229 1 iel = tmpNeigh;
2230 1 break;
2231 }
2232 }
2233
2234 1 return false;
2235 }
2236
2237 15 return true;
2238
2239 break;
2240
2241 4 case Hexahedron:
2242
2243 // face 0 ->adjacency
2244
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
 4 if ( pointLoc[2] < - 1 - GeoUtilities::GEO_TOL ) {
2245
2246 ✗ std::size_t numNeigh = m_listFaces.list()[faces[0]]->listNeighVolumes().size();
2247 ✗ for (std::size_t i = 0; i < numNeigh; ++i) {
2248 ✗ tmptyp = static_cast<ElementType>(m_listFaces.list()[faces[0]]->listNeighVolumes()[i]->typeVolume());
2249 ✗ tmpNeigh = m_listFaces.list()[faces[0]]->listNeighVolumes()[i]->idVolume();
2250
2251 ✗ if ( tmptyp != ityp || tmpNeigh != iel ) {
2252 ✗ ityp = tmptyp;
2253 ✗ iel = tmpNeigh;
2254 ✗ break;
2255 }
2256 }
2257
2258 ✗ return false;
2259 }
2260
2261 // face 3->adjacency
2262
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
 4 if ( pointLoc[2] > + 1 + GeoUtilities::GEO_TOL ) {
2263
2264 ✗ std::size_t numNeigh = m_listFaces.list()[faces[3]]->listNeighVolumes().size();
2265 ✗ for (std::size_t i = 0; i < numNeigh; ++i) {
2266 ✗ tmptyp = static_cast<ElementType>(m_listFaces.list()[faces[3]]->listNeighVolumes()[i]->typeVolume());
2267 ✗ tmpNeigh = m_listFaces.list()[faces[3]]->listNeighVolumes()[i]->idVolume();
2268
2269 ✗ if ( tmptyp != ityp || tmpNeigh != iel ) {
2270 ✗ ityp = tmptyp;
2271 ✗ iel = tmpNeigh;
2272 ✗ break;
2273 }
2274 }
2275
2276 ✗ return false;
2277 }
2278
2279 // face 2 ->adjacency
2280
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
 4 if ( pointLoc[1] > + 1 + GeoUtilities::GEO_TOL ) {
2281
2282 ✗ std::size_t numNeigh = m_listFaces.list()[faces[2]]->listNeighVolumes().size();
2283 ✗ for (std::size_t i = 0; i < numNeigh; ++i) {
2284 ✗ tmptyp = static_cast<ElementType>(m_listFaces.list()[faces[2]]->listNeighVolumes()[i]->typeVolume());
2285 ✗ tmpNeigh = m_listFaces.list()[faces[2]]->listNeighVolumes()[i]->idVolume();
2286
2287 ✗ if ( tmptyp != ityp || tmpNeigh != iel ) {
2288 ✗ ityp = tmptyp;
2289 ✗ iel = tmpNeigh;
2290 ✗ break;
2291 }
2292 }
2293
2294 ✗ return false;
2295 }
2296
2297 // face 5 ->adjacency
2298
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
 4 if ( pointLoc[1] < - 1 - GeoUtilities::GEO_TOL ) {
2299
2300 ✗ std::size_t numNeigh = m_listFaces.list()[faces[5]]->listNeighVolumes().size();
2301 ✗ for (std::size_t i = 0; i < numNeigh; ++i) {
2302 ✗ tmptyp = static_cast<ElementType>(m_listFaces.list()[faces[5]]->listNeighVolumes()[i]->typeVolume());
2303 ✗ tmpNeigh = m_listFaces.list()[faces[5]]->listNeighVolumes()[i]->idVolume();
2304
2305 ✗ if ( tmptyp != ityp || tmpNeigh != iel ) {
2306 ✗ ityp = tmptyp;
2307 ✗ iel = tmpNeigh;
2308 ✗ break;
2309 }
2310 }
2311
2312 ✗ return false;
2313 }
2314
2315 // face 1 ->adjacency
2316
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
 4 if ( pointLoc[0] < - 1 - GeoUtilities::GEO_TOL ) {
2317
2318 ✗ std::size_t numNeigh = m_listFaces.list()[faces[1]]->listNeighVolumes().size();
2319 ✗ for (std::size_t i = 0; i < numNeigh; ++i) {
2320 ✗ tmptyp = static_cast<ElementType>(m_listFaces.list()[faces[1]]->listNeighVolumes()[i]->typeVolume());
2321 ✗ tmpNeigh = m_listFaces.list()[faces[1]]->listNeighVolumes()[i]->idVolume();
2322
2323 ✗ if ( tmptyp != ityp || tmpNeigh != iel ) {
2324 ✗ ityp = tmptyp;
2325 ✗ iel = tmpNeigh;
2326 ✗ break;
2327 }
2328 }
2329
2330 ✗ return false;
2331 }
2332
2333 // face 4 ->adjacency
2334
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
 4 if ( pointLoc[0] > + 1 + GeoUtilities::GEO_TOL ) {
2335
2336 ✗ std::size_t numNeigh = m_listFaces.list()[faces[4]]->listNeighVolumes().size();
2337 ✗ for (std::size_t i = 0; i < numNeigh; ++i) {
2338 ✗ tmptyp = static_cast<ElementType>(m_listFaces.list()[faces[4]]->listNeighVolumes()[i]->typeVolume());
2339 ✗ tmpNeigh = m_listFaces.list()[faces[4]]->listNeighVolumes()[i]->idVolume();
2340
2341 ✗ if ( tmptyp != ityp || tmpNeigh != iel ) {
2342 ✗ ityp = tmptyp;
2343 ✗ iel = tmpNeigh;
2344 ✗ break;
2345 }
2346 }
2347
2348 ✗ return false;
2349 }
2350
2351 4 return true;
2352
2353 break;
2354
2355 ✗ case NullShape:
2356 case Node:
2357 case Prism:
2358 case Pyramid: {
2359 ✗ FEL_ERROR("unknown (or not implemented) element Shape")
2360 ✗ break;
2361 }
2362 }
2363
2364 ✗ ityp = ELEMENT_TYPE_COUNTER;
2365 ✗ iel = -1;
2366
2367 ✗ return false;
2368 3179 }
2369
2370 /***********************************************************************************/
2371 /***********************************************************************************/
2372
2373 4 bool GeometricMeshRegion::isPointInTetra(const felInt iel, const Point& p, std::array<double, 4>& bary, const double tol)
2374 {
2375
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 std::vector<felInt> verOfTetra(4, -1);
2376
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 std::vector<Point*> pntOfTetra(4, nullptr);
2377
2378 // Get points of the tetrahedra
2379
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 getOneElement(felisce::GeometricMeshRegion::Tetra4, iel, verOfTetra, pntOfTetra, 0);
2380
2381 // Check if the point is inside
2382
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 return GeoUtilities::isPointInTetra(pntOfTetra, p, bary, tol);
2383 4 }
2384
2385 /***********************************************************************************/
2386 /***********************************************************************************/
2387
2388 1753 bool GeometricMeshRegion::isPointInTri(const felInt iel, const Point& p, std::array<double, 3>& bary, const double tol)
2389 {
2390
1/2
✓ Branch 2 taken 1753 times.
✗ Branch 3 not taken.
 1753 std::vector<felInt> verOfTri(3, -1);
2391
1/2
✓ Branch 2 taken 1753 times.
✗ Branch 3 not taken.
 1753 std::vector<Point*> pntOfTri(3, nullptr);
2392
2393 // Get nodes of the triangle
2394
1/2
✓ Branch 1 taken 1753 times.
✗ Branch 2 not taken.
 1753 getOneElement(felisce::GeometricMeshRegion::Tria3, iel, verOfTri, pntOfTri, 0);
2395
2396 // Check if the point is inside
2397
1/2
✓ Branch 1 taken 1753 times.
✗ Branch 2 not taken.
 3506 return GeoUtilities::isPointInTri(pntOfTri, p, bary, tol);
2398 1753 }
2399
2400 /***********************************************************************************/
2401 /***********************************************************************************/
2402
2403 ✗ bool GeometricMeshRegion::isPointInSeg(const felInt iel, const Point& p, std::array<double, 2>& bary, const double tol)
2404 {
2405 ✗ std::vector<felInt> verOfSeg(2, -1);
2406 ✗ std::vector<Point*> pntOfSeg(2, nullptr);
2407
2408 // Get nodes of the segment
2409 ✗ getOneElement(felisce::GeometricMeshRegion::Seg2, iel, verOfSeg, pntOfSeg, 0);
2410
2411 // Check if the point is inside
2412 ✗ return GeoUtilities::isPointInSeg(pntOfSeg, p, bary, tol);
2413 }
2414
2415 /***********************************************************************************/
2416 /***********************************************************************************/
2417
2418 18912 bool GeometricMeshRegion::isVertexOfElt(const felInt idxVer, const ElementType eltTyp, const felInt eltIdx, felInt &verPos) const
2419 {
2420
1/2
✓ Branch 2 taken 18912 times.
✗ Branch 3 not taken.
 18912 std::vector<felInt> eltIdxVer(m_numPointsPerElt[eltTyp]);
2421
2422
1/2
✓ Branch 1 taken 18912 times.
✗ Branch 2 not taken.
 18912 getOneElement(eltTyp, eltIdx, eltIdxVer, 0, 0);
2423
2424
1/2
✓ Branch 1 taken 37825 times.
✗ Branch 2 not taken.
 37825 for(felInt i = 0; i < static_cast<felInt>(eltIdxVer.size()); ++i) {
2425
2/2
✓ Branch 1 taken 18912 times.
✓ Branch 2 taken 18913 times.
 37825 if ( eltIdxVer[i] == idxVer ) {
2426 18912 verPos = i;
2427 18912 return true;
2428 }
2429 }
2430
2431 ✗ return false;
2432 18912 }
2433
2434 /***********************************************************************************/
2435 /***********************************************************************************/
2436
2437 30 bool GeometricMeshRegion::isEdgeOfElt(const std::vector<felInt>& idxVer, const ElementType eltTyp, const felInt eltIdx, felInt &verPos) const
2438 {
2439 std::array<felInt,2> edgVer;
2440
1/2
✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
 30 std::vector<felInt> eltIdxVer(m_numPointsPerElt[eltTyp]);
2441 30 auto& mapEdgToVer = m_eltTypMapEdgToVer[eltTyp];
2442
2443
1/2
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
 30 getOneElement(eltTyp, eltIdx, eltIdxVer, 0, 0);
2444
2445
1/2
✓ Branch 1 taken 114 times.
✗ Branch 2 not taken.
 114 for(felInt i = 0; i < static_cast<felInt>(mapEdgToVer.size()); ++i) {
2446 114 edgVer[0] = eltIdxVer[mapEdgToVer[i][0]];
2447 114 edgVer[1] = eltIdxVer[mapEdgToVer[i][1]];
2448
2449
3/4
✓ Branch 4 taken 114 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 30 times.
✓ Branch 7 taken 84 times.
 114 if ( std::is_permutation(edgVer.begin(), edgVer.end(), idxVer.begin()) ) {
2450 30 verPos = i;
2451 30 return true;
2452 }
2453 }
2454
2455 ✗ return false;
2456 30 }
2457
2458 /***********************************************************************************/
2459 /***********************************************************************************/
2460
2461 3312 void GeometricMeshRegion::getVertexBall(const felInt idxVer, ElementType eltTyp, const felInt eltGrm, BallMap& ball) const
2462 {
2463 ElementType ngbTyp;
2464 felInt eltIdx, ngbIdx;
2465 felInt verPos, ngbPos;
2466
2467 3312 std::set<std::pair<ElementType, felInt>> mark;
2468
2469 3312 ball.clear();
2470
2471 // Get position of vertex in element and insert in ball
2472
2/4
✓ Branch 1 taken 3312 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3312 times.
✗ Branch 4 not taken.
 3312 if ( isVertexOfElt(idxVer, eltTyp, eltGrm, verPos) ) { // TODO warning here this methos should work with shapes and not with geoelements
2473
1/2
✓ Branch 2 taken 3312 times.
✗ Branch 3 not taken.
 3312 ball.push_back( {eltTyp, eltGrm, verPos} );
2474
1/2
✓ Branch 2 taken 3312 times.
✗ Branch 3 not taken.
 3312 mark.insert( {eltTyp, eltGrm} );
2475 }
2476 else {
2477 ✗ FEL_ERROR("Intersection: Provided vertex is not in element!");
2478 }
2479
2480 // Loop over the neighbors to set the ball
2481
2/2
✓ Branch 1 taken 18912 times.
✓ Branch 2 taken 3312 times.
 22224 for (std::size_t iel = 0; iel < ball.size(); ++iel) {
2482 18912 std::tie(eltTyp, eltIdx, verPos) = ball[iel];
2483
2484 // Loop over edges/faces/volumes sharing that vertex
2485 18912 auto& verToFac = m_eltTypMapVerToFac[eltTyp][verPos];
2486
2/2
✓ Branch 1 taken 37834 times.
✓ Branch 2 taken 18912 times.
 56746 for (std::size_t ifa = 0; ifa < verToFac.size(); ++ifa) {
2487
2488 // For every edge/face/volume surrounding the vertex, add neighbor element in ball
2489 37834 ngbTyp = eltTyp;
2490 37834 ngbIdx = eltIdx;
2491
3/4
✓ Branch 2 taken 37834 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 634 times.
✓ Branch 5 taken 37200 times.
 37834 if ( !getAdjElement(ngbTyp, ngbIdx, verToFac[ifa]) )
2492 634 continue;
2493
2494
3/4
✓ Branch 2 taken 37200 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 21600 times.
✓ Branch 5 taken 15600 times.
 37200 if ( !(mark.insert({ngbTyp, ngbIdx})).second )
2495 21600 continue;
2496
2497 // Get local position of the edge/face/volume
2498
1/2
✓ Branch 1 taken 15600 times.
✗ Branch 2 not taken.
 15600 isVertexOfElt(idxVer, ngbTyp, ngbIdx, ngbPos);
2499
2500 // Insert new element in ball
2501
1/2
✓ Branch 2 taken 15600 times.
✗ Branch 3 not taken.
 15600 ball.push_back( {ngbTyp, ngbIdx, ngbPos} );
2502 }
2503 }
2504 3312 }
2505
2506 /***********************************************************************************/
2507 /***********************************************************************************/
2508
2509 10 void GeometricMeshRegion::getEdgeShell(const std::vector<felInt>& idxVer, ElementType eltTyp, const felInt eltGrm, BallMap& shell) const
2510 {
2511 ElementType ngbTyp;
2512 felInt eltIdx, ngbIdx;
2513 felInt edgPos, ngbPos;
2514
2515 10 std::set<std::pair<ElementType, felInt>> mark;
2516
2517 10 shell.clear();
2518
2519 // Get position of edge in element and insert in shell
2520
2/4
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
 10 if ( isEdgeOfElt(idxVer, eltTyp, eltGrm, edgPos) ) { // TODO warning here this methos should work with shapes and not with geoelements
2521
1/2
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
 10 shell.push_back( {eltTyp, eltGrm, edgPos} );
2522
1/2
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
 10 mark.insert( {eltTyp, eltGrm} );
2523 }
2524 else {
2525 ✗ FEL_ERROR("Intersection: Provided vertex is not in element!");
2526 }
2527
2528 // Loop over the neighbors to set the shell
2529
2/2
✓ Branch 1 taken 30 times.
✓ Branch 2 taken 10 times.
 40 for (std::size_t iel = 0; iel < shell.size(); ++iel) {
2530 30 std::tie(eltTyp, eltIdx, edgPos) = shell[iel];
2531
2532 // Loop over edges/faces/volumes sharing that vertex
2533 30 auto& edgToFac = m_eltTypMapEdgToFac[eltTyp][edgPos];
2534
2/2
✓ Branch 1 taken 60 times.
✓ Branch 2 taken 30 times.
 90 for (std::size_t ifa = 0; ifa < edgToFac.size(); ++ifa) {
2535
2536 // For every edge/face/volume surrounding the vertex, add neighbor element in shell
2537 60 ngbTyp = eltTyp;
2538 60 ngbIdx = eltIdx;
2539
3/4
✓ Branch 2 taken 60 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 20 times.
✓ Branch 5 taken 40 times.
 60 if ( !getAdjElement(ngbTyp, ngbIdx, edgToFac[ifa]) )
2540 20 continue;
2541
2542
3/4
✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 20 times.
✓ Branch 5 taken 20 times.
 40 if ( !(mark.insert({ngbTyp, ngbIdx})).second )
2543 20 continue;
2544
2545 // Get local position of the edge/face/volume
2546
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 isEdgeOfElt(idxVer, ngbTyp, ngbIdx, ngbPos);
2547
2548 // Insert new element in shell
2549
1/2
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
 20 shell.push_back( {ngbTyp, ngbIdx, ngbPos} );
2550 }
2551 }
2552 10 }
2553
2554 /***********************************************************************************/
2555 /***********************************************************************************/
2556
2557 1100 void GeometricMeshRegion::getElementPatch(const ElementType typGrm, const felInt eltGrm, std::vector<seedElement>& elementsInPatch) const
2558 {
2559 1100 BallMap ball;
2560 ElementType eltTyp;
2561 felInt eltIdx, verPos;
2562
1/2
✓ Branch 2 taken 1100 times.
✗ Branch 3 not taken.
 1100 std::vector<felInt> veridxOfElt(GeometricMeshRegion::m_numPointsPerElt[typGrm], 0);
2563
2564 // Clean vector
2565 1100 elementsInPatch.clear();
2566
2567 // Get element vertices
2568
1/2
✓ Branch 1 taken 1100 times.
✗ Branch 2 not taken.
 1100 getOneElement(typGrm, eltGrm, veridxOfElt, 0);
2569
2570 // Loop over the shape vertices
2571 1100 const int numVer = eltEnumToFelNameGeoEle[typGrm].second->shape().numVertex();
2572
2/2
✓ Branch 0 taken 3300 times.
✓ Branch 1 taken 1100 times.
 4400 for (int ive = 0; ive < numVer; ++ive) {
2573
2574 // Get ball
2575
1/2
✓ Branch 2 taken 3300 times.
✗ Branch 3 not taken.
 3300 getVertexBall(veridxOfElt[ive], typGrm, eltGrm, ball);
2576
2577 // Save elements
2578
2/2
✓ Branch 1 taken 18900 times.
✓ Branch 2 taken 3300 times.
 22200 for (std::size_t iel = 0; iel < ball.size(); ++iel) {
2579 18900 std::tie(eltTyp, eltIdx, verPos) = ball[iel];
2580
2581 // Insert elements in vector
2582
1/2
✓ Branch 2 taken 18900 times.
✗ Branch 3 not taken.
 18900 elementsInPatch.push_back( {eltTyp, eltIdx} );
2583 }
2584 }
2585
2586
1/2
✓ Branch 3 taken 1100 times.
✗ Branch 4 not taken.
 1100 std::sort( elementsInPatch.begin(), elementsInPatch.end() );
2587
2/4
✓ Branch 5 taken 1100 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 1100 times.
✗ Branch 10 not taken.
 1100 elementsInPatch.erase( std::unique( elementsInPatch.begin(), elementsInPatch.end() ), elementsInPatch.end() );
2588 1100 }
2589
2590 /***********************************************************************************/
2591 /***********************************************************************************/
2592
2593 // std::size_t GeometricMeshRegion::computeConnectedComponents(std::map<ElementType,std::vector<std::size_t>>& mapCco) const
2594 // {
2595
2596 // // Get the array to use to know the number of neighbor TODO D.C. this should be done using polymorphism
2597 // const int* numNeigh;
2598 // if ( m_domainDim == GeoMesh2D ) {
2599 // numNeigh = m_numEdgesPerElt;
2600 // }
2601 // else if ( m_domainDim == GeoMesh3D ) {
2602 // numNeigh = m_numFacesPerElt;
2603 // }
2604 // else {
2605 // FEL_ERROR("Not implemented for this case");
2606 // }
2607
2608 // std::queue<std::pair<ElementType, felInt>> queue;
2609 // ElementType eltTyp, neiTyp;
2610 // felInt prvGlbIdx;
2611 // felInt eltIdx, neiIdx;
2612 // felInt nbrElt, glbNbrElt = 0;
2613
2614 // // Global number of elements and initialize mapCco
2615 // for (std::size_t i = 0; i < m_bagElementTypeDomain.size(); ++i) {
2616 // eltTyp = m_bagElementTypeDomain[i];
2617 // nbrElt = m_numElements[eltTyp];
2618
2619 // if ( nbrElt > 0 )
2620 // mapCco[eltTyp].assign(nbrElt,0);
2621
2622 // glbNbrElt += nbrElt;
2623 // }
2624
2625 // // Set counter to zero
2626 // felInt cnt = 0;
2627 // // Set connected component progressive idx to 0
2628 // std::size_t ccoIdx = 0;
2629
2630 // // Find all connected components
2631 // NextCco:
2632 // ccoIdx++;
2633
2634 // // Find connected component germ
2635 // for(felInt eltGlbIdx = prvGlbIdx; eltGlbIdx < glbNbrElt; ++eltGlbIdx) {
2636
2637 // // Get element type and id by element
2638 // getTypeElemFromIdElem(eltGlbIdx, eltTyp, eltIdx);
2639
2640 // // If element not marked yet
2641 // if ( mapCco[neiTyp][neiIdx] == 0 ) {
2642
2643 // // Update
2644 // prvGlbIdx = eltGlbIdx+1;
2645
2646 // // Insert in the queue
2647 // queue.emplace(std::make_pair(eltTyp, eltIdx));
2648
2649 // // Set connected component
2650 // mapCco[neiTyp][neiIdx] = ccoIdx;
2651
2652 // // Increase counter
2653 // cnt++;
2654
2655 // break;
2656 // }
2657 // }
2658
2659 // // Color the current connected component
2660 // while ( cnt < glbNbrElt ) {
2661
2662 // // Get front element
2663 // eltTyp = queue.front().first;
2664 // eltIdx = queue.front().second;
2665
2666 // //- For every neighbor of iElt:
2667 // for(int ine = 0; ine < numNeigh[eltTyp]; ++ine) {
2668
2669 // // Get the neighbor element if it exists
2670 // neiTyp = eltTyp;
2671 // neiIdx = eltIdx;
2672 // if ( !getAdjElement(neiTyp, neiIdx, ine) )
2673 // continue;
2674
2675 // //- If the neighbor exists
2676 // if ( mapCco[neiTyp][neiIdx] == 0 ) {
2677 // // Insert in the queue
2678 // queue.emplace(std::make_pair(neiTyp, neiIdx));
2679 // // Set connected component
2680 // mapCco[neiTyp][neiIdx] = ccoIdx;
2681 // // Increase counter
2682 // cnt++;
2683 // }
2684 // }
2685
2686 // // Remove first element from list
2687 // queue.pop();
2688
2689 // // Check if we found all the connected component elements
2690 // if ( queue.empty() )
2691 // goto NextCco;
2692 // }
2693
2694 // //- Check if all elements have been marked
2695 // FEL_CHECK( cnt == glbNbrElt, "Problem detecting connected components! ");
2696
2697 // // Return number of connected components
2698 // return ccoIdx;
2699 // }
2700
2701 }
2702