GCC Code Coverage Report

Directory: ./
File: Tools/fe_utilities.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 67 67 100.0%
Branches: 35 50 70.0%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors: Vicente Mataix Ferrandiz
13 //
14
15 // System includes
16 #include <iostream>
17 #include <unordered_set>
18 #include <unordered_map>
19
20 // External includes
21
22 // Project includes
23 #include "Tools/fe_utilities.hpp"
24 #include "FiniteElement/refElement.hpp"
25 #include "FiniteElement/currentFiniteElement.hpp"
26
27 namespace felisce
28 {
29
30 namespace FEUtilities
31 {
32
33 // List of the FE triple product compatible elements
34 static std::unordered_set<GeometricMeshRegion::ElementType> compatible_triple_product_fe =
35 {
36 GeometricMeshRegion::Tetra4,
37 GeometricMeshRegion::Hexa8
38 };
39 static std::unordered_map<GeometricMeshRegion::ElementType, std::vector<std::array<int, 4>>> vertex_tetrahedra =
40 {
41 {GeometricMeshRegion::Tetra4, {{0,1,2,3}}},
42 {GeometricMeshRegion::Hexa8, { // http://www.math.udel.edu/~szhang/research/p/subtettest.pdf
43 {0,1,2,6}, {1,2,3,7}, {2,3,0,4}, {3,0,1,5}, {4,7,6,2}, {7,6,5,1},
44 {6,5,4,0}, {5,4,7,3}, {1,0,4,7}, {3,2,6,5}, {0,3,7,6}, {2,1,5,4},
45 {0,1,2,4}, {1,2,3,5}, {2,3,0,6}, {3,0,1,7}, {4,7,6,0}, {7,6,5,3},
46 {6,5,4,2}, {5,4,7,1}, {4,0,3,5}, {2,6,7,1}, {3,7,4,2}, {1,5,6,0},
47 {0,1,3,4}, {1,2,0,5}, {2,3,1,6}, {3,0,2,7}, {4,7,5,0}, {5,4,6,1}, {6,5,7,2}, {7,6,4,3}
48 }}
49 };
50
51 /***********************************************************************************/
52 /***********************************************************************************/
53
54 4 double volumeTetrahedra(const std::vector<Point*>& elemPoint, const std::vector<std::array<int, 4>>& vertexTetrahedra)
55 {
56 // Auxiliary double
57 4 double auxiliary_double = 0.0;
58
59 // Loop over the tetrahedra
60
2/2
✓ Branch 1 taken 35 times.
✓ Branch 2 taken 2 times.
 37 for(std::size_t i = 0; i < vertexTetrahedra.size(); ++i) {
61 // Get the vertices
62 35 const auto& r_array = vertexTetrahedra[i];
63 35 const Point* p0 = elemPoint[r_array[0]];
64 35 const Point* p1 = elemPoint[r_array[1]];
65 35 const Point* p2 = elemPoint[r_array[2]];
66 35 const Point* p3 = elemPoint[r_array[3]];
67
68 35 const double x10 = p1->x() - p0->x();
69 35 const double y10 = p1->y() - p0->y();
70 35 const double z10 = p1->z() - p0->z();
71
72 35 const double x20 = p2->x() - p0->x();
73 35 const double y20 = p2->y() - p0->y();
74 35 const double z20 = p2->z() - p0->z();
75
76 35 const double x30 = p3->x() - p0->x();
77 35 const double y30 = p3->y() - p0->y();
78 35 const double z30 = p3->z() - p0->z();
79
80 // Compute the volume
81 35 const double volume = (x10 * y20 * z30 - x10 * y30 * z20 + y10 * z20 * x30 - y10 * x20 * z30 + z10 * x20 * y30 - z10 * y20 * x30)/6.0;
82
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 33 times.
 35 if (volume < 0.0) return -1.0;
83 33 auxiliary_double += volume;
84 }
85
86 // Return the volume
87 2 return auxiliary_double;
88 }
89
90 /***********************************************************************************/
91 /***********************************************************************************/
92
93 8 double volumeFE(GeometricMeshRegion& rMesh, GeometricMeshRegion::ElementType& eltType, felInt iel, std::vector<Point*>& elemPoint, std::vector<felInt>& elemIdPoint, felInt& ielSupportDof)
94 {
95 IGNORE_UNUSED_ARGUMENT(ielSupportDof);
96 // Generate element from previous mesh
97 8 const GeoElement* geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
98
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 const RefElement* refEle = geoEle->defineFiniteEle(eltType, 0, rMesh);
99
1/2
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
 8 std::vector<DegreeOfExactness> degree_of_exactness(1, DegreeOfExactness::DegreeOfExactness_1);
100
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 auto element = CurrentFiniteElement(*refEle,*geoEle, degree_of_exactness);
101
102
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 rMesh.getOneElement(eltType, iel, elemIdPoint, elemPoint);
103
104 // Compute integration points
105
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 element.updateMeas(0, elemPoint);
106
107
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 16 return element.measure();
108 8 }
109
110 /***********************************************************************************/
111 /***********************************************************************************/
112
113 12 bool CheckVolumeIsInverted(
114 GeometricMeshRegion& rMesh,
115 const bool jacobianCheck
116 )
117 {
118 /* Loop function */
119 // Execute the loop
120 12 std::unordered_set<GeometricMeshRegion::ElementType> aux_bag_set;
121 12 const auto& elements_bag = rMesh.bagElementTypeDomain();
122
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 12 times.
✗ Branch 8 not taken.
 12 std::copy(elements_bag.begin(), elements_bag.end(), std::inserter(aux_bag_set, aux_bag_set.begin()));
123
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 6 times.
 12 if (!jacobianCheck) {
124
2/2
✓ Branch 5 taken 11 times.
✓ Branch 6 taken 4 times.
 15 for (auto& compatible_eltType : compatible_triple_product_fe) {
125 11 const int num_elem = rMesh.numElements(compatible_eltType);
126
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 7 times.
 11 if (num_elem > 0) {
127
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 aux_bag_set.erase(compatible_eltType);
128
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 const std::vector<GeometricMeshRegion::ElementType> aux_bag(1, compatible_eltType);
129
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 const auto& connectivity = vertex_tetrahedra[compatible_eltType];
130 8 auto function_triple_product = [&rMesh,&connectivity](GeometricMeshRegion::ElementType& eltType, felInt iel, std::vector<Point*>& elemPoint, std::vector<felInt>& elemIdPoint, felInt& ielSupportDof) {
131 IGNORE_UNUSED_ARGUMENT(ielSupportDof);
132 4 rMesh.getOneElement(eltType, iel, elemIdPoint, elemPoint);
133 4 const double volume = volumeTetrahedra(elemPoint, connectivity);
134
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 4 if (volume < 0.0) {
135 2 std::cout << "WARNING:: The element " << iel << ", which type is " << eltType << " is inverted. Check connectivity" << std::endl;
136 2 return 1.0;
137 } else {
138 2 return 0.0;
139 }
140 4 };
141
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 const double aux = FEUtilities::LoopOverElements(rMesh, aux_bag, &function_triple_product);
142
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 4 if (aux > 0.0) {
143 2 return true;
144 }
145
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
 4 }
146 }
147 }
148 10 std::vector<GeometricMeshRegion::ElementType> effective_elements_bag;
149
2/4
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 10 times.
✗ Branch 7 not taken.
 10 std::copy(aux_bag_set.begin(), aux_bag_set.end(), std::back_inserter(effective_elements_bag));
150
151 // Generic jacobian check
152 /* Loop function */
153 8 auto function_jacobian = [&rMesh](GeometricMeshRegion::ElementType& eltType, felInt iel, std::vector<Point*>& elemPoint, std::vector<felInt>& elemIdPoint, felInt& ielSupportDof) {
154 8 const double volume = volumeFE(rMesh, eltType, iel, elemPoint, elemIdPoint, ielSupportDof);
155
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 8 if (volume < 0.0) {
156 4 std::cout << "WARNING:: The element " << iel << ", which type is " << eltType << " is inverted. Check connectivity" << std::endl;
157 4 return 1.0;
158 } else {
159 4 return 0.0;
160 }
161 10 };
162
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 const double aux = FEUtilities::LoopOverElements(rMesh, effective_elements_bag, &function_jacobian);
163
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 6 times.
 10 if (aux > 0.0) {
164 4 return true;
165 } else {
166 6 return false;
167 }
168 12 }
169
170 } // namespace FEUtilities
171
172 } // namespace felisce
173