GCC Code Coverage Report

Directory: ./
File: FiniteElement/refElement.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 124 218 56.9%
Branches: 78 193 40.4%
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-F. Gerbeau
13 //
14
15 // System includes
16 #include <cmath>
17
18 // External includes
19
20 // Project includes
21 #include "refElement.hpp"
22
23 namespace felisce
24 {
25 13536 RefElement::RefElement(const std::string& name,
26 const RefShape& refShape,
27 const BasisFunction& basisFunction,
28 const ListOfQuadratureRule& listOfQuadratureRule,
29 const DegreeOfFreedomType* typeDOF,
30 const DegreeOfFreedomSupport* supportDOF,
31 const int* idDOF,
32 const Point* node,int numNodeVertex,int numNodeEdge,int numNodeFace,int numNodeVolume,
33 const RefElement** boundaryRefElement,
34 const RefElement** boundaryBoundaryRefElement,
35 13536 bool edgesCanHaveDifferentNumberOfDof):
36
1/2
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
 13536 m_name(name),
37 13536 m_refShape(refShape),
38 13536 m_basisFunction(basisFunction),
39 // m_listOfQuadratureRule(listOfQuadratureRule),
40 13536 m_typeDOF(typeDOF),
41 13536 m_supportDOF(supportDOF),
42 13536 m_idSupportOfDOF(idDOF),
43 13536 m_node(node),
44 13536 m_numNodeVertex(numNodeVertex),
45 13536 m_numNodeEdge(numNodeEdge),
46 13536 m_numNodeFace(numNodeFace),
47 13536 m_numNodeVolume(numNodeVolume),
48 13536 m_numNode(m_numNodeVertex + m_numNodeEdge + m_numNodeFace + m_numNodeVolume),
49 13536 m_numDOF(basisFunction.size()),
50 13536 m_numVertex(m_refShape.numVertex()),
51 13536 m_numEdge(m_refShape.numEdge()),
52 13536 m_numFace(m_refShape.numFace()),
53 13536 m_numCoor(m_refShape.numCoor()),
54 13536 m_boundaryRefElement(boundaryRefElement),
55 27072 m_edgesCanHaveDifferentNumberOfDof(edgesCanHaveDifferentNumberOfDof)//false by default
56 // m_boundaryBoundaryRefElement(boundaryBoundaryRefElement)
57 {
58 IGNORE_UNUSED_ARGUMENT(listOfQuadratureRule);
59 IGNORE_UNUSED_ARGUMENT(boundaryBoundaryRefElement);
60
61 FEL_DEBUG("Constructing " << m_name << std::endl);
62
1/2
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
 13536 initStaticMember();
63 13536 m_numDOFNodeVertex = 0;
64 13536 m_numDOFNodeEdge = 0;
65 13536 m_numDOFNodeFace = 0;
66 13536 m_numDOFNodeVolume = 0;
67 13536 m_numDOFEdge = 0;
68 13536 m_numDOFFace = 0;
69 13536 m_numDOFVolume = 0;
70
2/2
✓ Branch 0 taken 98700 times.
✓ Branch 1 taken 13536 times.
 112236 for(int i=0; i<m_numDOF; i++) {
71
5/8
✓ Branch 0 taken 57528 times.
✓ Branch 1 taken 32148 times.
✓ Branch 2 taken 5076 times.
✓ Branch 3 taken 1692 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2256 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 98700 switch (m_supportDOF[i]) {
72 57528 case DOF_NODE_VERTEX:
73 57528 m_numDOFNodeVertex++;
74 57528 break;
75 32148 case DOF_NODE_EDGE:
76 32148 m_numDOFNodeEdge++;
77 32148 m_numDOFEdge++;
78 32148 break;
79 5076 case DOF_NODE_FACE:
80 5076 m_numDOFNodeFace++;
81 5076 m_numDOFFace++;
82 5076 break;
83 1692 case DOF_NODE_VOLUME:
84 1692 m_numDOFNodeVolume++;
85 1692 m_numDOFVolume++;
86 1692 break;
87 ✗ case DOF_EDGE:
88 ✗ m_numDOFEdge++;
89 ✗ break;
90 2256 case DOF_FACE:
91 2256 m_numDOFFace++;
92 2256 break;
93 ✗ case DOF_VOLUME:
94 ✗ m_numDOFVolume++;
95 ✗ break;
96 // Default case should appear with a warning at compile time instead of an error in runtime
97 // (that's truly the point of using enums as switch cases)
98 // default:
99 //felisce_error("Something is wrong is typeDOF", __FILE__, __LINE__);
100 //break;
101 }
102 }
103 13536 m_numDOFNode = m_numDOFNodeVertex + m_numDOFNodeEdge + m_numDOFNodeFace + m_numDOFNodeVolume;
104
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 13536 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 13536 FEL_ASSERT( m_numDOF == (m_numDOFNodeVertex + m_numDOFEdge + m_numDOFFace + m_numDOFVolume) );
105 13536 m_numDOFPerEdge = 0;
106 13536 m_numDOFPerFace = 0;
107
4/4
✓ Branch 0 taken 12972 times.
✓ Branch 1 taken 564 times.
✓ Branch 2 taken 11844 times.
✓ Branch 3 taken 1128 times.
 13536 if ( m_numEdge > 0 && !m_edgesCanHaveDifferentNumberOfDof ) {
108
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 11844 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 11844 FEL_ASSERT( m_numDOFEdge % m_numEdge == 0 );
109 11844 m_numDOFPerEdge = (int) (m_numDOFEdge / m_numEdge);
110
2/2
✓ Branch 0 taken 1128 times.
✓ Branch 1 taken 564 times.
 1692 } else if ( m_numEdge > 0 ){
111
1/2
✓ Branch 1 taken 1128 times.
✗ Branch 2 not taken.
 1128 m_numDOFPerSingleEdge.resize(m_numEdge,0);
112
2/2
✓ Branch 0 taken 2820 times.
✓ Branch 1 taken 1128 times.
 3948 for (int i(0); i<m_numDOFEdge; ++i) {
113 2820 m_numDOFPerSingleEdge[m_idSupportOfDOF[m_numDOFNodeVertex + i]]++;
114 }
115 }
116
2/2
✓ Branch 0 taken 10716 times.
✓ Branch 1 taken 2820 times.
 13536 if ( m_numFace > 0 ) {
117
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10716 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 10716 FEL_ASSERT( m_numDOFFace % m_numFace == 0 );
118 10716 m_numDOFPerFace = (int) (m_numDOFFace / m_numFace);
119 }
120
1/2
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
 13536 m_idDOFVertex.resize(m_numVertex);
121
1/2
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
 13536 m_idDOFEdge.resize(m_numEdge);
122
1/2
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
 13536 m_idDOFFace.resize(m_numFace);
123
2/2
✓ Branch 0 taken 98700 times.
✓ Branch 1 taken 13536 times.
 112236 for(int i=0; i<m_numDOF; i++) {
124
4/5
✓ Branch 0 taken 57528 times.
✓ Branch 1 taken 32148 times.
✓ Branch 2 taken 7332 times.
✓ Branch 3 taken 1692 times.
✗ Branch 4 not taken.
 98700 switch (m_supportDOF[i]) {
125 57528 case DOF_NODE_VERTEX:
126
1/2
✓ Branch 2 taken 57528 times.
✗ Branch 3 not taken.
 57528 m_idDOFVertex[ m_idSupportOfDOF[i] ].push_back(i);
127 57528 break;
128 32148 case DOF_NODE_EDGE:
129 case DOF_EDGE:
130 // in the case of m_edgesCanHaveDifferentNumberOfDof=true some of those vectors will be empty
131
1/2
✓ Branch 2 taken 32148 times.
✗ Branch 3 not taken.
 32148 m_idDOFEdge[ m_idSupportOfDOF[i] ].push_back(i);
132 32148 break;
133 7332 case DOF_NODE_FACE:
134 case DOF_FACE:
135
1/2
✓ Branch 2 taken 7332 times.
✗ Branch 3 not taken.
 7332 m_idDOFFace[ m_idSupportOfDOF[i] ].push_back(i);
136 7332 break;
137 1692 case DOF_NODE_VOLUME:
138 case DOF_VOLUME:
139
1/2
✓ Branch 1 taken 1692 times.
✗ Branch 2 not taken.
 1692 m_idDOFVolume.push_back(i);
140 1692 break;
141 // Default case should appear with a warning at compile time instead of an error in runtime
142 // (that's truly the point of using enums as switch cases)
143 // default:
144
145 // felisce_error("Something is wrong is typeDOF", __FILE__, __LINE__);
146 // break;
147 }
148 }
149
150 // This loop is not very general, in case m_numEdge = 2*m_numDOFEdge it does not work anymore maybe replace m_numEdge with m_numDOFEdge
151
1/2
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
 13536 m_nodeDOF.resize(m_numDOF);
152
2/2
✓ Branch 0 taken 12408 times.
✓ Branch 1 taken 1128 times.
 13536 if ( ! m_edgesCanHaveDifferentNumberOfDof ) {
153
2/2
✓ Branch 0 taken 90240 times.
✓ Branch 1 taken 12408 times.
 102648 for(int i=0; i<m_numDOF; i++) {
154
5/6
✓ Branch 0 taken 51888 times.
✓ Branch 1 taken 29328 times.
✓ Branch 2 taken 5076 times.
✓ Branch 3 taken 1692 times.
✓ Branch 4 taken 2256 times.
✗ Branch 5 not taken.
 90240 switch (m_supportDOF[i]) {
155 51888 case DOF_NODE_VERTEX:
156 51888 m_nodeDOF[i] = &(m_node[m_idSupportOfDOF[i]]); // vertices come first in the nodes list
157 51888 break;
158 29328 case DOF_NODE_EDGE:
159 29328 m_nodeDOF[i] = &(m_node[m_numVertex+m_idSupportOfDOF[i]]); // nodes-edges comes just after vertices in the nodes list
160 29328 break;
161 5076 case DOF_NODE_FACE:
162 5076 m_nodeDOF[i] = &(m_node[m_numVertex+m_numEdge+m_idSupportOfDOF[i]]);
163 5076 break;
164 1692 case DOF_NODE_VOLUME:
165 1692 m_nodeDOF[i] = &(m_node[m_numVertex+m_numEdge+m_numFace+m_idSupportOfDOF[i]]);
166 1692 break;
167 2256 case DOF_EDGE:
168 case DOF_FACE:
169 case DOF_VOLUME:
170 2256 m_nodeDOF[i] = nullptr;
171 2256 break;
172 }
173 }
174 } else {
175
2/2
✓ Branch 0 taken 8460 times.
✓ Branch 1 taken 1128 times.
 9588 for(int i=0; i<m_numDOF; i++) {
176
2/6
✓ Branch 0 taken 5640 times.
✓ Branch 1 taken 2820 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 8460 switch (m_supportDOF[i]) {
177 5640 case DOF_NODE_VERTEX:
178 5640 m_nodeDOF[i] = &(m_node[m_idSupportOfDOF[i]]); // vertices come first in the nodes list
179 5640 break;
180 2820 case DOF_NODE_EDGE:
181 2820 m_nodeDOF[i] = &(m_node[i]);// We assume that the points are listed in the same order as the dofs. This is used for instance in PrismsP1xP2
182 2820 break;
183 ✗ case DOF_NODE_FACE:
184 ✗ m_nodeDOF[i] = &(m_node[m_numVertex+m_numDOFEdge+m_idSupportOfDOF[i]]);
185 ✗ break;
186 ✗ case DOF_NODE_VOLUME:
187 ✗ m_nodeDOF[i] = &(m_node[m_numVertex+m_numDOFEdge+m_numDOFFace+m_idSupportOfDOF[i]]);
188 ✗ break;
189 ✗ case DOF_EDGE:
190 case DOF_FACE:
191 case DOF_VOLUME:
192 ✗ m_nodeDOF[i] = nullptr;
193 ✗ break;
194 }
195 }
196 }
197 13536 }
198
199 /***********************************************************************************/
200 /***********************************************************************************/
201
202 13536 void RefElement::initStaticMember()
203 {
204 // static bool first_time = true;
205 // if(first_time){
206 FEL_DEBUG("Initializing static members of RefElement" << std::endl);
207
208
2/4
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13536 times.
✗ Branch 5 not taken.
 13536 m_nameSupportDOF.insert(std::pair<DegreeOfFreedomSupport,std::string>(DOF_NODE_VERTEX,"DOF_NODE_VERTEX"));
209
2/4
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13536 times.
✗ Branch 5 not taken.
 13536 m_nameSupportDOF.insert(std::pair<DegreeOfFreedomSupport,std::string>(DOF_NODE_EDGE,"DOF_NODE_EDGE"));
210
2/4
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13536 times.
✗ Branch 5 not taken.
 13536 m_nameSupportDOF.insert(std::pair<DegreeOfFreedomSupport,std::string>(DOF_NODE_FACE,"DOF_NODE_FACE"));
211
2/4
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13536 times.
✗ Branch 5 not taken.
 13536 m_nameSupportDOF.insert(std::pair<DegreeOfFreedomSupport,std::string>(DOF_NODE_VOLUME,"DOF_NODE_VOLUME"));
212
2/4
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13536 times.
✗ Branch 5 not taken.
 13536 m_nameSupportDOF.insert(std::pair<DegreeOfFreedomSupport,std::string>(DOF_EDGE,"DOF_EDGE"));
213
2/4
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13536 times.
✗ Branch 5 not taken.
 13536 m_nameSupportDOF.insert(std::pair<DegreeOfFreedomSupport,std::string>(DOF_FACE,"DOF_FACE"));
214
2/4
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13536 times.
✗ Branch 5 not taken.
 13536 m_nameSupportDOF.insert(std::pair<DegreeOfFreedomSupport,std::string>(DOF_VOLUME,"DOF_VOLUME"));
215
216
2/4
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13536 times.
✗ Branch 5 not taken.
 13536 m_nameTypeDOF.insert(std::pair<DegreeOfFreedomType,std::string>(DOF_VALUE,"DOF_VALUE"));
217
2/4
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13536 times.
✗ Branch 5 not taken.
 13536 m_nameTypeDOF.insert(std::pair<DegreeOfFreedomType,std::string>(DOF_X_DERIVATIVE,"DOF_X_DERIVATIVE"));
218
2/4
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13536 times.
✗ Branch 5 not taken.
 13536 m_nameTypeDOF.insert(std::pair<DegreeOfFreedomType,std::string>(DOF_Y_DERIVATIVE,"DOF_Y_DERIVATIVE"));
219
2/4
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13536 times.
✗ Branch 5 not taken.
 13536 m_nameTypeDOF.insert(std::pair<DegreeOfFreedomType,std::string>(DOF_Z_DERIVATIVE,"DOF_Z_DERIVATIVE"));
220
2/4
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13536 times.
✗ Branch 5 not taken.
 13536 m_nameTypeDOF.insert(std::pair<DegreeOfFreedomType,std::string>(DOF_MEAN,"DOF_MEAN"));
221
2/4
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13536 times.
✗ Branch 5 not taken.
 13536 m_nameTypeDOF.insert(std::pair<DegreeOfFreedomType,std::string>(DOF_FLUX,"DOF_FLUX"));
222
2/4
✓ Branch 1 taken 13536 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13536 times.
✗ Branch 5 not taken.
 13536 m_nameTypeDOF.insert(std::pair<DegreeOfFreedomType,std::string>(DOF_CIRCULATION,"DOF_CIRCULATION"));
223
224 // }
225 //first_time = false;
226 13536 }
227
228 /***********************************************************************************/
229 /***********************************************************************************/
230
231 ✗ void RefElement::check() const
232 {
233 ✗ std::cout << "Check Reference Element " << m_name << std::endl;
234 ✗ if(m_numDOFNode) std::cout << " Nodal DOF (the following matrix must be identity)" << std::endl;
235 ✗ double delta_ij=0.;
236
237 ✗ for(int i=0; i<m_numDOF; i++) {
238 ✗ if(dofIsANode(i)) {
239 ✗ std::cout << " ";
240 ✗ for(int j=0; j<m_numDOF; j++) {
241 ✗ if(dofIsANode(j)) {
242 ✗ FEL_ASSERT(m_nodeDOF[j]); // check if dof j is really a node
243 ✗ switch(m_typeDOF[i]) {
244 ✗ case DOF_VALUE:
245 ✗ delta_ij = m_basisFunction.phi(i,*(m_nodeDOF[j]));
246 ✗ break;
247 ✗ case DOF_X_DERIVATIVE:
248 ✗ delta_ij = m_basisFunction.dPhi(i,0,*(m_nodeDOF[j]));
249 ✗ break;
250 ✗ case DOF_Y_DERIVATIVE:
251 ✗ delta_ij = m_basisFunction.dPhi(i,1,*(m_nodeDOF[j]));
252 ✗ break;
253 ✗ case DOF_Z_DERIVATIVE:
254 ✗ delta_ij = m_basisFunction.dPhi(i,2,*(m_nodeDOF[j]));
255 ✗ break;
256 ✗ case DOF_NORMAL_DERIVATIVE:
257 case DOF_MEAN:
258 case DOF_FLUX:
259 case DOF_CIRCULATION:
260 ✗ break;
261 }
262 ✗ std::cout << std::fabs(delta_ij) << " "; // fabs is just to avoid the minus sign for -0 (due rounding off errors)
263 }
264 }
265 ✗ std::cout << std::endl;
266 }
267 }
268 }
269
270 /***********************************************************************************/
271 /***********************************************************************************/
272
273 ✗ void RefElement::print(int verbose,std::ostream& c) const
274 {
275 ✗ if(verbose) {
276 ✗ c << "RefElement: " << std::endl;
277 ✗ c << " name: " << m_name << std::endl;
278 ✗ c << " shape: " << m_refShape.name() << std::endl;
279 ✗ c << " basisFunction: " << m_basisFunction.name() << std::endl;
280 ✗ c << " numDOFNode: " << m_numDOFNode ;
281 ✗ c << " ( " << m_numDOFNodeVertex << " on vertices, ";
282 ✗ c << m_numDOFNodeEdge << " on edges, ";
283 ✗ c << m_numDOFNodeFace << " on faces, ";
284 ✗ c << m_numDOFNodeVolume << " on volume) " << std::endl;
285 ✗ if ( ! m_edgesCanHaveDifferentNumberOfDof ) {
286 ✗ c << " numDOFEdge: " << m_numDOFEdge << " (numDOFPerEdge=" << m_numDOFPerEdge << ")" << std::endl;
287 } else {
288 ✗ c << " numDOFEdge: " << m_numDOFEdge << std::endl;
289 ✗ for (int i(0); i<m_numEdge; ++i)
290 ✗ c << "Edge id: "<<i<<" number of dof supported by this edge " << m_idDOFEdge[i].size() << std::endl;
291 }
292 ✗ c << " numDOFFace: " << m_numDOFFace << " (numDOFPerFace=" << m_numDOFPerFace << ")" << std::endl;
293 ✗ c << " numDOFVolume: " << m_numDOFVolume << std::endl;
294 ✗ if(m_numNode) c << " Nodes coordinates : " << std::endl;
295
296 ✗ for(int ip=0; ip<m_numNode; ip++) {
297 ✗ c << " Node " << ip << ": ";
298 ✗ std::cout << m_node[ip].x() << ", " << m_node[ip].y() << ", " << m_node[ip].z() << std::endl;
299 }
300
301 ✗ c << " Degree Of Freedom (DOF Type / Support type / Support id) : " << std::endl;
302 ✗ for(int i=0; i<m_basisFunction.size(); i++) {
303 ✗ c << " DOF " << i << " : " << m_nameTypeDOF[m_typeDOF[i]] << " / "
304 ✗ << m_nameSupportDOF[m_supportDOF[i]] << " / " << m_idSupportOfDOF[i] << std::endl;
305 }
306
307 ✗ if(m_numDOFNodeVertex) {
308 ✗ c << " List of DOF by vertex " << std::endl;
309 ✗ for(int i=0; i<m_numVertex; i++) {
310 ✗ c << " Vertex " << i << " : ";
311 ✗ for(std::size_t j=0; j<m_idDOFVertex[i].size(); j++) {
312 ✗ c << m_idDOFVertex[i][j] << " ";
313 }
314 ✗ c << std::endl;
315 }
316 }
317 ✗ if(m_numDOFEdge) {
318 ✗ c << " List of DOF by edge " << std::endl;
319 ✗ for(int i=0; i<m_numEdge; i++) {
320 ✗ if ( m_idDOFEdge[i].size() > 0 ){
321 ✗ c << " Edge " << i << " : ";
322 ✗ for(std::size_t j=0; j<m_idDOFEdge[i].size(); j++) {
323 ✗ c << m_idDOFEdge[i][j] << " ";
324 }
325 ✗ c << std::endl;
326 }
327 }
328 }
329 ✗ if(m_numDOFFace) {
330 ✗ c << " List of DOF by face " << std::endl;
331 ✗ for(int i=0; i<m_numFace; i++) {
332 ✗ c << " Face " << i << " : ";
333 ✗ for(std::size_t j=0; j<m_idDOFFace[i].size(); j++) {
334 ✗ c << m_idDOFFace[i][j] << " ";
335 }
336 ✗ c << std::endl;
337 }
338 }
339 ✗ if(m_numDOFVolume) {
340 ✗ c << " List of DOF on the volume " << std::endl;
341 ✗ for(int i=0; i<m_numDOFVolume; i++) {
342 ✗ c << m_idDOFVolume[i] << " ";
343 }
344 }
345 }
346 }
347 }
348