GCC Code Coverage Report

Directory: ./
File: FiniteElement/geoElement.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 83 128 64.8%
Branches: 50 245 20.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 <iostream>
17
18 // External includes
19
20 // Project includes
21 #include "FiniteElement/geoElement.hpp"
22 #include "FiniteElement/basisFunction.hpp"
23
24 namespace felisce
25 {
26 12972 GeoElement::GeoElement(const std::string& name,const RefShape& shape,const BasisFunction& basisFunction,const double* pointCoor,const int numEdge,const int numFace,
27 const int* pointOfEdge,const int* pointOfFace,const int* edgeOfFace,const bool* orientationOfEdge,
28 12972 const GeoElement** boundaryGeoElement,const GeoElement** boundaryBoundaryGeoElement,const RefElement& defaultFiniteEle):
29 12972 m_name(name),m_shape(shape),m_basisFunction(basisFunction),m_numCoor(m_basisFunction.numCoor()),
30 12972 m_numPoint(m_basisFunction.size()),m_numEdge(numEdge),m_numFace(numFace),
31
6/6
✓ Branch 0 taken 7332 times.
✓ Branch 1 taken 5640 times.
✓ Branch 2 taken 4512 times.
✓ Branch 3 taken 2820 times.
✓ Branch 4 taken 1692 times.
✓ Branch 5 taken 1128 times.
 12972 m_numBdEle( m_numCoor == 3 ? numFace : ( m_numCoor == 2 ? m_numEdge : (m_numCoor == 1 ? 2 : 0) ) ),
32 12972 m_pointCoor(pointCoor),
33 12972 m_ptOfEd(pointOfEdge),m_ptOfFa(pointOfFace),m_edOfFa(edgeOfFace),m_orientEd(orientationOfEdge),
34 12972 m_boundaryGeoElement(boundaryGeoElement),m_boundaryBoundaryGeoElement(boundaryBoundaryGeoElement),
35 12972 m_defaultFiniteEle(defaultFiniteEle)
36 {
37 FEL_DEBUG("Constructing " << m_name << std::endl);
38
39
2/4
✓ Branch 0 taken 12972 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 12972 times.
✗ Branch 5 not taken.
 12972 m_numPointPerEdge = new int[m_numEdge];
40
2/4
✓ Branch 0 taken 12972 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 12972 times.
✗ Branch 5 not taken.
 12972 m_numPointPerFace = new int[m_numFace];
41
2/4
✓ Branch 0 taken 12972 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 12972 times.
✗ Branch 5 not taken.
 12972 m_numEdgePerFace = new int[m_numFace];
42
2/4
✓ Branch 0 taken 12972 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 12972 times.
✗ Branch 5 not taken.
 12972 m_indexPointPerEdge = new int[m_numEdge];
43
2/4
✓ Branch 0 taken 12972 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 12972 times.
✗ Branch 5 not taken.
 12972 m_indexPointPerFace = new int[m_numFace];
44
2/4
✓ Branch 0 taken 12972 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 12972 times.
✗ Branch 5 not taken.
 12972 m_indexEdgePerFace = new int[m_numFace];
45
46
4/5
✓ Branch 0 taken 1128 times.
✓ Branch 1 taken 1692 times.
✓ Branch 2 taken 4512 times.
✓ Branch 3 taken 5640 times.
✗ Branch 4 not taken.
 12972 switch (m_numCoor) {
47 1128 case 0: // test for geoElement{NULL,Node} vm 2012/09
48 // m_numEdge = m_numFace = 0
49 1128 break;
50 1692 case 1:
51 // 1D edge = the element
52 1692 m_numPointPerEdge[0] = m_numPoint;
53 1692 m_indexPointPerEdge[0] = 0;
54 //! create/update boundary element arrays: (bdEle=2 points)
55
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1692 times.
 1692 FEL_ASSERT( m_numBdEle == 2);
56 1692 m_ptOfBdEle = m_ptOfEd; //={0,1,...}: but only the extremal points {0,1} are used.
57
1/2
✓ Branch 1 taken 1692 times.
✗ Branch 2 not taken.
 1692 m_numPointPerBdEle = new int[2]; // 2 vertices per segment!...
58 1692 m_numPointPerBdEle[0] = 1;
59 1692 m_numPointPerBdEle[1] = 1;
60
1/2
✓ Branch 1 taken 1692 times.
✗ Branch 2 not taken.
 1692 m_indexPointPerBdEle = new int[2];
61 1692 m_indexPointPerBdEle[0] = 0;
62 1692 m_indexPointPerBdEle[1] = 1;
63 1692 break;
64 4512 case 2:
65 // 2D Boundary = edges
66
2/2
✓ Branch 0 taken 16356 times.
✓ Branch 1 taken 4512 times.
 20868 for(int i=0; i<m_numEdge; i++) {
67 16356 m_numPointPerEdge[i] = m_boundaryGeoElement[i]->numPoint();
68 }
69 4512 m_indexPointPerEdge[0] = 0;
70
2/2
✓ Branch 0 taken 11844 times.
✓ Branch 1 taken 4512 times.
 16356 for(int i=1; i<m_numEdge; i++) {
71 11844 m_indexPointPerEdge[i] = m_indexPointPerEdge[i-1] + m_boundaryGeoElement[i-1]->numPoint();
72 }
73 // 2D face = the element
74 4512 m_numPointPerFace[0] = m_numPoint;
75 4512 m_indexPointPerFace[0] = 0;
76 4512 m_numEdgePerFace[0] = m_numEdge;
77 4512 m_indexEdgePerFace[0] = 0;
78 //! update boundary element pointers: (bdEle=edges)
79
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4512 times.
 4512 FEL_ASSERT( m_numBdEle == m_numEdge );
80 4512 m_ptOfBdEle = m_ptOfEd;
81 4512 m_numPointPerBdEle = m_numPointPerEdge;
82 4512 m_indexPointPerBdEle = m_indexPointPerEdge;
83 4512 break;
84 5640 case 3:
85 // 3D Boundary of boundary = edges
86
2/2
✓ Branch 0 taken 52452 times.
✓ Branch 1 taken 5640 times.
 58092 for(int i=0; i<m_numEdge; i++) {
87 52452 m_numPointPerEdge[i] = m_boundaryBoundaryGeoElement[i]->numPoint();
88 }
89 5640 m_indexPointPerEdge[0] = 0;
90
2/2
✓ Branch 0 taken 46812 times.
✓ Branch 1 taken 5640 times.
 52452 for(int i=1; i<m_numEdge; i++) {
91 46812 m_indexPointPerEdge[i] = m_indexPointPerEdge[i-1] + m_boundaryBoundaryGeoElement[i-1]->numPoint();
92 }
93 // 3D Boundary = faces
94
2/2
✓ Branch 0 taken 28764 times.
✓ Branch 1 taken 5640 times.
 34404 for(int i=0; i<m_numFace; i++) {
95 28764 m_numPointPerFace[i] = m_boundaryGeoElement[i]->numPoint();
96 28764 m_numEdgePerFace[i] = m_boundaryGeoElement[i]->numEdge();
97 }
98 5640 m_indexPointPerFace[0] = 0;
99 5640 m_indexEdgePerFace[0] = 0;
100
2/2
✓ Branch 0 taken 23124 times.
✓ Branch 1 taken 5640 times.
 28764 for(int i=1; i<m_numFace; i++) {
101 23124 m_indexPointPerFace[i] = m_indexPointPerFace[i-1] + m_boundaryGeoElement[i-1]->numPoint();
102 23124 m_indexEdgePerFace[i] = m_indexEdgePerFace[i-1] + m_boundaryGeoElement[i-1]->numEdge();
103 }
104 //! update boundary element pointers: (bdEle=edges)
105
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5640 times.
 5640 FEL_ASSERT( m_numBdEle == m_numFace );
106 5640 m_ptOfBdEle = m_ptOfFa;
107 5640 m_numPointPerBdEle = m_numPointPerFace;
108 5640 m_indexPointPerBdEle = m_indexPointPerFace;
109 5640 break;
110 ✗ default:
111 ✗ std::cout << "m_numCoor = " << m_numCoor << std::endl;
112 ✗ felisce_error("Incorrect number of coordinates", __FILE__, __LINE__);
113 break;
114 }
115 12972 }
116
117 /***********************************************************************************/
118 /***********************************************************************************/
119
120 25944 GeoElement::~GeoElement()
121 {
122
1/2
✓ Branch 0 taken 12972 times.
✗ Branch 1 not taken.
 12972 delete [] m_numPointPerEdge;
123
1/2
✓ Branch 0 taken 12972 times.
✗ Branch 1 not taken.
 12972 delete [] m_numPointPerFace;
124
1/2
✓ Branch 0 taken 12972 times.
✗ Branch 1 not taken.
 12972 delete [] m_numEdgePerFace;
125
1/2
✓ Branch 0 taken 12972 times.
✗ Branch 1 not taken.
 12972 delete [] m_indexPointPerEdge;
126
1/2
✓ Branch 0 taken 12972 times.
✗ Branch 1 not taken.
 12972 delete [] m_indexPointPerFace;
127
1/2
✓ Branch 0 taken 12972 times.
✗ Branch 1 not taken.
 12972 delete [] m_indexEdgePerFace;
128 12972 }
129
130 /***********************************************************************************/
131 /***********************************************************************************/
132
133 80655 const RefElement* GeoElement::defineFiniteEle(const GeometricMeshRegion::ElementType eltType, const int typeOfFiniteElement, GeometricMeshRegion& mesh) const
134 {
135
3/4
✓ Branch 0 taken 78005 times.
✓ Branch 1 taken 1682 times.
✓ Branch 2 taken 968 times.
✗ Branch 3 not taken.
 80655 switch (typeOfFiniteElement) {
136 78005 case 0:
137 78005 return &GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second->defaultFiniteEle();
138 1682 case 1:
139 1682 return &GeometricMeshRegion::eltEnumToFelNameGeoEle[GeometricMeshRegion::eltLinearToEltQuad[eltType]].second->defaultFiniteEle();
140 968 case 2:
141 // for bubble finite element, the boundaries are P1, we must not map them.
142
2/2
✓ Branch 1 taken 488 times.
✓ Branch 2 taken 480 times.
 968 if(mesh.domainDim() == m_numCoor)
143 488 return &GeometricMeshRegion::eltEnumToFelNameGeoEle[GeometricMeshRegion::eltLinearToEltBubble[eltType]].second->defaultFiniteEle();
144 else
145 480 return &GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second->defaultFiniteEle();
146 // case 3:
147 // return &refElementSegmentP3H;
148 ✗ default:
149 ✗ FEL_ERROR("Not yet implemented for this type of finite element option");
150 break;
151 }
152 return nullptr;
153 }
154
155 /***********************************************************************************/
156 /***********************************************************************************/
157
158 ✗ void GeoElement::print(int verbose,std::ostream& c) const
159 {
160 ✗ if(verbose) {
161 ✗ c << "GeoElement: " << std::endl;
162 ✗ c << " name: " << m_name << std::endl;
163 ✗ c << " shape: " << m_shape.name() << std::endl;
164 ✗ c << " numCoor: " << m_numCoor << std::endl;
165 ✗ c << " numPoint: " << m_numPoint << std::endl;
166 ✗ c << " numEdge: " << m_numEdge << std::endl;
167 ✗ c << " numFace: " << m_numFace << std::endl;
168 ✗ c << " numBdEle: " << m_numBdEle << std::endl;
169 ✗ c << " numPointPerEdge: ";
170 ✗ for(int ie=0; ie<m_numEdge; ie++) c << m_numPointPerEdge[ie] << " ";
171 ✗ c << std::endl;
172 ✗ c << " numPointPerFace: ";
173 ✗ for(int jf=0; jf<m_numFace; jf++) c << m_numPointPerFace[jf] << " ";
174 ✗ c << std::endl;
175 ✗ c << " numEdgePerFace: ";
176 ✗ for(int jf=0; jf<m_numFace; jf++) c << m_numEdgePerFace[jf] << " ";
177 ✗ c << std::endl;
178 ✗ c << " Point coordinates: " << std::endl;
179 ✗ for(int ip=0; ip<m_numPoint; ip++) {
180 ✗ c << " ";
181 ✗ std::cout << pointCoor(ip, 0) << ", " << pointCoor(ip, 1) << ", " << pointCoor(ip, 2) << std::endl;
182 }
183 ✗ c << " Points of Edges: " << std::endl;
184 ✗ for(int ie=0; ie<m_numEdge; ie++) {
185 ✗ c << " ";
186 ✗ for(int ip=0; ip<m_numPointPerEdge[ie]; ip++) c << pointOfEdge(ie,ip) << " ";
187 ✗ c << std::endl;
188 }
189 ✗ c << " Points of Faces: " << std::endl;
190 ✗ for(int jf=0; jf<m_numFace; jf++) {
191 ✗ c << " ";
192 ✗ for(int ip=0; ip<m_numPointPerFace[jf]; ip++) c << pointOfFace(jf,ip) << " ";
193 ✗ c << std::endl;
194 }
195 ✗ c << " Edges of Faces (orientation in parenthesis): " << std::endl;
196 bool orient;
197 ✗ for(int jf=0; jf<m_numFace; jf++) {
198 ✗ c << " ";
199 ✗ for(int ie=0; ie<m_numEdgePerFace[jf]; ie++) {
200 ✗ int edge = edgeOfFace(jf,ie,orient);
201 ✗ c << edge << " (" << orient << ") ";
202 }
203 ✗ c << std::endl;
204 }
205 }
206 }
207
208 /***********************************************************************************/
209 /***********************************************************************************/
210
211 #define OUTJFG
212 #ifndef OUTJFG
213
214 /************************************************************************
215 * QuadraticSegment
216 *
217 * 0-----2-----1
218 *
219 *************************************************************************/
220 static const int m_ptOfEdQuadSeg[3] = {
221 0,1,2
222 }
223
224 static const double refcoor_P2_1D[ 9 ] = {
225 0. , 0. , 0.,
226 1. , 0. , 0.,
227 0.5 , 0. , 0.
228 }
229
230 const GeoElement quadSeg("Quadratic Segment",SEGMENT,2,3,1,0,3,0,0,refcoor_P2_1D,m_ptOfEdQuadSeg,NULL,NULL,NULL,basisFunctionP2Seg,NULL);
231
232
233 /************************************************************************
234 * QuadraticTriangle
235 *
236 * 2
237 * / \
238 * / \
239 * 5 4
240 * / \
241 * / \
242 * 0-----3-----1
243 *
244 *************************************************************************/
245 static const int m_ptOfEdQuadTria[9] = {
246 0,1,3, 1,2,4, 2,0,5
247 }
248 static const int m_ptOfFaQuadTria[6] = {
249 0,1,2,3,4,5
250 }
251 static const double refcoor_P2_2D[ 18 ] = {
252 0. , 0. , 0.,
253 1. , 0. , 0.,
254 0. , 1. , 0.,
255 0.5 , 0. , 0.,
256 0.5 , 0.5 , 0.,
257 0. , 0.5 , 0.
258 }
259 const GeoElement quadraticTriangle("Quadratic Triangle",TRIANGLE,2,6,3,1,3,6,3,refcoor_P2_2D,m_ptOfEdQuadTria,m_ptOfFaQuadTria,NULL,NULL,
260 basisFunctionP2Tria,&quadSeg);
261
262 /* GeoElement(ReferenceShape shape,int numCoor,int numPoint,int numEdge,int numFace,
263 int numPointPerEdge,int numPointPerFace,int numEdgePerFace,
264 const int* eToP,const int* fToP,const int* fToE,const BasisFunction& basisFunction):*/
265
266 static const double refcoor_P2_3D[ 30 ] = {
267 0. , 0. , 0. ,
268 1. , 0. , 0. ,
269 0. , 1. , 0. ,
270 0. , 0. , 1. ,
271 0.5 , 0. , 0. ,
272 0.5, 0.5 , 0. ,
273 0. , 0.5 , 0. ,
274 0. , 0. , 0.5,
275 0.5, 0. , 0.5,
276 0. , 0.5 , 0.5
277 }
278
279 static const double refcoor_Q1_3D[ 24 ] = {
280 0. , 0. , 0. ,
281 1. , 0. , 0. ,
282 1. , 1. , 0. ,
283 0. , 1. , 0. ,
284 0. , 0. , 1. ,
285 1. , 0. , 1. ,
286 1. , 1. , 1. ,
287 0. , 1. , 1.
288 }
289
290 //===================================================================================================
291 /*
292 -- LinearQuad
293
294 3-------2
295 ! !
296 ! !
297 0-------1
298
299
300 */
301 static const int m_ptOfEdBilinearQuad[ 8 ] = {
302 0, 1, 1, 2, 2, 3, 3, 0
303 }
304 static const int m_ptOfFaBilinearQuad[4] = {
305 0,1,2,3
306 }
307
308 static const double refcoor_Q1_2D[ 12 ] = {
309 0. , 0. , 0.,
310 1. , 0. , 0.,
311 1. , 1. , 0.,
312 0. , 1. , 0.
313 }
314
315 const GeoElement bilinearQuadrangle("Bilinear Quadrangle",QUAD,2,4,4,1,2,4,4,refcoor_Q1_2D,m_ptOfEdBilinearQuad,m_ptOfFaBilinearQuad,NULL,NULL,
316 basisFunctionQ1Quad,&geoElementSegmentP1);
317
318 /*
319 -- QuadraticQuad
320
321 3---6---2
322 ! !
323 7 8 5
324 ! !
325 0---4---1
326
327 */
328 static const int m_ptOfEdBiquadQuad[ 12 ] = {
329 0, 1, 4, 1, 2, 5, 2, 3, 6, 3, 0, 7
330 }
331 static const int m_ptOfFaBiquadQuad[ 9 ] = {
332 0, 1, 2, 3, 4, 5, 6, 7, 8
333 }
334
335 static const double refcoor_Q2_2D[ 27 ] = {
336 0. , 0. , 0.,
337 1. , 0. , 0.,
338 1. , 1. , 0.,
339 0. , 1. , 0.,
340 0.5 , 0. , 0.,
341 1. , 0.5 , 0.,
342 0.5 , 1. , 0.,
343 0. , 0.5 , 0.,
344 0.5 , 0.5 , 0.
345 }
346
347 const GeoElement biquadQuadrangle("Biquadratic Quadrangle",QUAD,2,9,4,1,3,9,4,refcoor_Q2_2D,m_ptOfEdBiquadQuad,m_ptOfFaBiquadQuad,NULL,NULL,
348 basisFunctionQ2Quad,&quadSeg);
349
350
351 //===================================================================================================
352 /*
353 -- Prism
354
355 5
356 / \
357 / \
358 3 ---- 4
359
360 2
361 / \
362 / \
363 0 ---- 1
364
365 */
366
367 static const int m_ptOfEdPrismP1P1[18] = {
368 0,1, 1,2, 2,0, 3,4, 4,5, 5,3, 0,3, 1,4, 2,5
369 }
370
371 static const int m_ptOfFaPrismP1P1[18] = {
372 0,1,2, 3,4,5, 0,1,4,3, 1,4,5,2, 0,3,5,2
373 }
374
375 static const double refcoor_PrismP1P1[ 18 ] = {
376 0. , 0. , 0.,
377 1. , 0. , 0.,
378 0. , 1. , 0.,
379 0. , 0. , 1.,
380 1. , 0. , 1.,
381 0. , 1. , 1.,
382 }
383 /*
384 static const GeoElement* bdEle_Prism[ 5 ] =
385 {
386 &linearTriangle, &linearTriangle, &bilinearQuadrangle,&bilinearQuadrangle,&bilinearQuadrangle
387 }
388
389 */
390 //const GeoElement prismP1P1("Prism P1 P1",PRISM,2,6,3,1,3,6,3,refcoor_PrismP1P1,m_ptOfEdPrismP1P1,m_ptOfFaPrismP1P1,NULL,NULL,
391 // basisFunctionP2Tria,bdEle_Prism);
392 #endif
393 }
394