GCC Code Coverage Report

Directory: ./
File: Geometry/Tools/interpolator.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 50 72 69.4%
Branches: 32 132 24.2%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors:
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Geometry/Tools/interpolator.hpp"
21
22 namespace felisce
23 {
24
25 4 Interpolator::Interpolator(GeometricMeshRegion* mesh, const bool allowProjection, const double tolRescaling, const int verbose) :
26 4 Locator(mesh, allowProjection, tolRescaling, verbose)
27 {
28
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 1 times.
 4 if ( m_mesh->domainDim() == GeometricMeshRegion::GeoMesh2D ) {
29
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
 3 if ( m_mesh->bagElementTypeDomain().size() != 1 )
30 ✗ FEL_ERROR("Interpolator 2D works only for domain made of Tria3.");
31
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
 3 if ( m_mesh->bagElementTypeDomain()[0] != GeometricMeshRegion::Tria3 )
32 ✗ FEL_ERROR("Interpolator 2D works only for domain made of Tria3.");
33
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 3 if ( allowProjection ) {
34 ✗ if ( m_mesh->bagElementTypeDomainBoundary().size() != 1 )
35 ✗ FEL_ERROR("Interpolator 2D works only for boundaries made of Seg2.");
36 ✗ if ( m_mesh->bagElementTypeDomainBoundary()[0] != GeometricMeshRegion::Seg2 )
37 ✗ FEL_ERROR("Interpolator 2D works only for boundaries made of Seg2.");
38 }
39 }
40
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 else if ( m_mesh->domainDim() == GeometricMeshRegion::GeoMesh3D ) {
41
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
 1 if ( m_mesh->bagElementTypeDomain().size() != 1 )
42 ✗ FEL_ERROR("Interpolator 3D works only for domain made of Tetra4.");
43
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
 1 if ( m_mesh->bagElementTypeDomain()[0] != GeometricMeshRegion::Tetra4 )
44 ✗ FEL_ERROR("Interpolator 3D works only for domain made of Tetra4.");
45
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if ( allowProjection ) {
46 ✗ if ( m_mesh->bagElementTypeDomainBoundary().size() != 1 )
47 ✗ FEL_ERROR("Interpolator 3D works only for boundaries made of Tria3");
48 ✗ if ( m_mesh->bagElementTypeDomainBoundary()[0] != GeometricMeshRegion::Tria3 )
49 ✗ FEL_ERROR("Interpolator 3D works only for boundaries made of Tria3");
50 }
51 }
52 else {
53 ✗ FEL_ERROR("Interpolator works only in 3D or 2D");
54 }
55 4 }
56
57 /***********************************************************************************/
58 /***********************************************************************************/
59
60 82 void Interpolator::interpolate(std::vector<felInt>& indElt, std::vector<double>& intVal, std::vector<Point>& lstPoints)
61 {
62
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 81 times.
 82 if ( m_mesh->domainDim() == GeometricMeshRegion::GeoMesh3D ){
63
64 1 interpolate3D(indElt, intVal, lstPoints);
65 }
66
1/2
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
 81 else if ( m_mesh->domainDim() == GeometricMeshRegion::GeoMesh2D ) {
67
68 81 interpolate2D(indElt, intVal, lstPoints);
69
70 } else {
71
72 ✗ FEL_ERROR("Interpolator: Not implemented yet for dimension " + std::to_string(m_mesh->domainDim()) + " !" );
73 }
74 82 }
75
76 /***********************************************************************************/
77 /***********************************************************************************/
78
79 81 void Interpolator::interpolate2D(std::vector<felInt>& idxElt, std::vector<double>& intVal, const std::vector<Point>& lstPoints)
80 {
81 Point point;
82 81 seedElement element;
83 std::array<double, 3> bary;
84
85 // Initialize return arguments
86
1/2
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
 81 idxElt.assign(lstPoints.size(), -1);
87
1/2
✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
 81 intVal.assign(3*lstPoints.size(), -1.);
88
89
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
 81 if ( m_verbose > 3 )
90 ✗ std::cout << "Interpolator: process in progress." << std::endl;
91
92 // For each point of the structure mesh
93
2/2
✓ Branch 1 taken 1753 times.
✓ Branch 2 taken 81 times.
 1834 for (std::size_t k = 0; k < lstPoints.size(); ++k) {
94 1753 point = lstPoints[k];
95
96
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1753 times.
 1753 if ( m_verbose > 3 )
97 ✗ std::cout << "Interpolator: Processing point " << k << std::endl;
98
99 // Locate the point and even project it if necessary - when it is inside the box but not inside the mesh
100
2/4
✓ Branch 2 taken 1753 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1753 times.
 1753 if ( !( findPointInMesh(element, lstPoints[k]) ) )
101 // Failure
102 ✗ FEL_ERROR( "Interpolator: point " + std::to_string(k) + " has not been localized even with the projection on the mesh (if it has been allowed)!" );
103
104 // Compute barycentric coordinates
105
1/2
✓ Branch 1 taken 1753 times.
✗ Branch 2 not taken.
 1753 m_mesh->isPointInTri(element.second, point, bary, m_tol);
106
107 // Save results
108 1753 idxElt[k] = element.second;
109
110 1753 intVal[3*k] = bary[1];
111 1753 intVal[3*k+1] = bary[2];
112 1753 intVal[3*k+2] = bary[0];
113 }
114
115
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
 81 if ( m_verbose > 3 )
116 ✗ std::cout << "Interpolator: process completed." << std::endl;
117 81 }
118
119 /***********************************************************************************/
120 /***********************************************************************************/
121
122 1 void Interpolator::interpolate3D(std::vector<felInt>& idxElt, std::vector<double>& intVal, const std::vector<Point>& lstPoints)
123 {
124 Point point;
125 1 seedElement element;
126 std::array<double,4> bary;
127
128 // Initialize return arguments
129
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 idxElt.assign(lstPoints.size(), -1);
130
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 intVal.assign(4*lstPoints.size(), -1.);
131
132
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if ( m_verbose > 3 )
133 ✗ std::cout << "Interpolator: process in progress." << std::endl;
134
135 // For each point of the structure mesh
136
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
 5 for (std::size_t k = 0; k < lstPoints.size(); ++k) {
137 4 point = lstPoints[k];
138
139
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if ( m_verbose > 3 )
140 ✗ std::cout << "Interpolator: Processing point " << k << std::endl;
141
142 // Locate the point and even project it if necessary - when it is inside the box but not inside the mesh
143
2/4
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
 4 if ( !( findPointInMesh(element, lstPoints[k]) ) )
144 // Failure
145 ✗ FEL_ERROR( "Interpolator: point " + std::to_string(k) + " has not been localized even with the projection on the mesh (if it has been allowed)!" );
146
147 // Compute barycentric coordinates
148
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_mesh->isPointInTetra(element.second, point, bary, m_tol);
149
150 // Save results
151 4 idxElt[k] = element.second;
152
153 4 intVal[4*k] = bary[2];
154 4 intVal[4*k+1] = bary[3];
155 4 intVal[4*k+2] = bary[1];
156 4 intVal[4*k+3] = bary[0];
157 }
158
159
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if ( m_verbose > 3 )
160 ✗ std::cout << "Interpolator: process completed." << std::endl;
161 1 }
162
163 }
164