GCC Code Coverage Report

Directory: ./
File: FiniteElement/quadratureRule.hpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 54 59 91.5%
Branches: 6 14 42.9%
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 #ifndef QUADRATURE_RULE_HPP
16 #define QUADRATURE_RULE_HPP
17
18 // System includes
19 #include <cassert>
20 #include <iostream>
21
22 // External includes
23
24 // Project includes
25 #include "Geometry/point.hpp"
26 #include "Core/felisce.hpp"
27 #include "FiniteElement/refShape.hpp"
28
29 ///@name felisce Globals
30 ///@{
31
32 ///@}
33 ///@name Type Definitions
34 ///@{
35
36 ///@}
37 ///@name Enum's
38 ///@{
39
40 ///@}
41 ///@name Functions
42 ///@{
43
44 ///@}
45 ///@name felisce Classes
46 ///@{
47 namespace felisce
48 {
49 /**
50 * @class QuadraturePoint
51 * @brief The class for the quadrature points i.e. (x,y,z,weight)
52 * @authors J-F. Gerbeau
53 * @date 12/2009
54 */
55 class QuadraturePoint
56 : public Point
57 {
58 public:
59 ///@name Type Definitions
60 ///@{
61
62 ///@}
63 ///@name Life Cycle
64 ///@{
65
66 /// Default constructor
67 128592 QuadraturePoint(){};
68
69 /**
70 * @brief Constructor with three coordinates and one weight
71 */
72 244776 QuadraturePoint(
73 const double x,
74 const double y,
75 const double z,
76 const double w
77 244776 ) : Point(x,y,z),
78 244776 m_weight( w )
79 {
80 244776 }
81
82 /**
83 * @brief Constructor with two coordinates and one weight
84 */
85 21432 QuadraturePoint(
86 const double x,
87 const double y,
88 const double w
89 21432 ) : Point(x,y,0.),
90 21432 m_weight( w )
91 {
92 21432 }
93
94 /**
95 * @brief Constructor with one coordinate and one weight
96 */
97 3948 QuadraturePoint(
98 const double x,
99 const double w
100 3948 ) : Point(x,0.,0.),
101 3948 m_weight( w )
102 {
103 3948 }
104
105 /**
106 * @brief Constructor with point and one weight
107 */
108 QuadraturePoint(
109 const Point& rPoint,
110 const double w
111 ) : Point(rPoint),
112 m_weight( w )
113 {
114 }
115
116 /*! Default destructor */
117 ~QuadraturePoint() = default;
118
119 //! Copy constructor, to avoid warning due to user-declared destructor.
120 QuadraturePoint(const QuadraturePoint& P)
121 : Point(P.m_coor),
122 m_weight(P.m_weight)
123 {}
124
125 ///@}
126 ///@name Operators
127 ///@{
128
129 // Assignment.
130 128592 QuadraturePoint& operator=(const QuadraturePoint& P) {
131 128592 m_coor = P.m_coor;
132 128592 m_weight = P.m_weight;
133
134 128592 return *this;
135 }
136
137 ///@}
138 ///@name Operations
139 ///@{
140
141 ///@}
142 ///@name Access
143 ///@{
144
145 /// Retrieves the weight of the QuadraturePoint
146 230792497 inline double weight() const {
147 230792497 return m_weight;
148 }
149
150 ///@}
151 ///@name Inquiry
152 ///@{
153
154 ///@}
155 ///@name Input and output
156 ///@{
157
158 ✗ void print(int verbose,std::ostream& c=std::cout) const {
159 ✗ if(verbose) {
160 ✗ c << "[" << m_coor[0] << ", " << m_coor[1] << ", " << m_coor[2] << "]" << ", weight = " << m_weight;
161 }
162 }
163
164 ///@}
165 ///@name Friends
166 ///@{
167
168 ///@}
169 protected:
170 ///@name Protected static Member Variables
171 ///@{
172
173 ///@}
174 ///@name Protected member Variables
175 ///@{
176 ///@}
177 ///@name Protected Operators
178 ///@{
179
180 ///@}
181 ///@name Protected Operations
182 ///@{
183
184 ///@}
185 ///@name Protected Access
186 ///@{
187
188 ///@}
189 ///@name Protected Inquiry
190 ///@{
191
192 ///@}
193 ///@name Protected LifeCycle
194 ///@{
195
196 ///@}
197 private:
198 ///@name Private static Member Variables
199 ///@{
200
201 ///@}
202 ///@name Private member Variables
203 ///@{
204
205 double m_weight = 0.0; /// Weight associated to the quadrature point
206
207 ///@}
208 ///@name Private Operators
209 ///@{
210
211 ///@}
212 ///@name Private Operations
213 ///@{
214
215 ///@}
216 ///@name Private Access
217 ///@{
218
219 ///@}
220 ///@name Private Inquiry
221 ///@{
222
223 ///@}
224 ///@name Private LifeCycle
225 ///@{
226
227 ///@}
228 };
229
230 /**
231 * @class QuadratureRule
232 * @brief The class of quadrature rules
233 * @authors J-F. Gerbeau
234 * @date 12/2009
235 * This class contains the quadrature points and the associated weights.
236 * @par How to add a new quadrature rule ?
237 * Suppose you want to add the quadrature rule qr_Pipo on a tetrahedron:
238 * in the file quadratureRule.hpp, add
239 * @code
240 * extern const QuadratureRule quadratureRulePipo;
241 * @endcode
242 * In the file definitionGlobalVariables.cpp, add the quadrature rule (follow the example of another one) and add it to the array of all quadrature rule on a tetrahedron, namely :
243 * @code
244 * static const QuadratureRule quad_rule_on_tetra[] =
245 * {quadratureRuleTetra1pt,quadratureRuleTetra4pt,quadratureRuleTetra5pt,quadRulePipo,
246 * quadratureRuleTetra15pt,quadratureRuleTetra64pt}
247 * @endcode
248 * Be careful : it is mandatory that the quadrature are ordered by increasing degree of exactness, and two different
249 * quadrature rules of this list cannot share the same degree of exactness. The reason of this limitation is due to the way
250 * ListOfQuadratureRule::quadratureRuleByExactness is built.
251 * Finally, do not forget to increase the number of quadrature in the file definitionGlobalVariables.cpp
252 * @code
253 * const ListOfQuadratureRule listQuadratureRuleTetrahedron("listQuadratureRuleTetrahedron",
254 * 6,quad_rule_on_tetra);
255 * @endcode
256 */
257 class QuadratureRule
258 {
259 public:
260 ///@name Type Definitions
261 ///@{
262
263 ///@}
264 ///@name Life Cycle
265 ///@{
266
267 /**
268 * @brief Default constructor
269 */
270
3/6
✓ Branch 2 taken 15228 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 15228 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 15228 times.
✗ Branch 10 not taken.
 15228 QuadratureRule(){};
271
272 /**
273 * @brief Constructor with arguments
274 */
275 QuadratureRule(
276 const QuadraturePoint* pt,
277 const std::string name,
278 const RefShape& shape,
279 const int numQuadraturePoint,
280 const int defOfExact
281 );
282
283 /**
284 * @brief Constructor with arguments (combination)
285 */
286 QuadratureRule(
287 const QuadratureRule& rQuadratureRule1,
288 const QuadratureRule& rQuadratureRule2,
289 const std::string name,
290 const RefShape& shape,
291 const double Coefficient = 1.0
292 );
293
294 /*! Default destructor */
295 46248 ~QuadratureRule() = default;
296
297 //! Copy constructor, to avoid warning due to user-declared destructor.
298 15228 QuadratureRule(const QuadratureRule& QR)
299 15228 : m_point(QR.m_point),
300 15228 m_shape(QR.m_shape),
301 15228 m_name(QR.m_name),
302 15228 m_numQuadraturePoint(QR.m_numQuadraturePoint),
303 15228 m_degreeOfExactness(QR.m_degreeOfExactness)
304 {
305 15228 }
306
307 ///@}
308 ///@name Operators
309 ///@{
310
311 // Assignment.
312 QuadratureRule& operator=(const QuadratureRule& QR);
313
314 ///@}
315 ///@name Operations
316 ///@{
317
318 ///@}
319 ///@name Access
320 ///@{
321
322 /*! Return the name of the quadrature rule */
323 ✗ std::string name() const {
324 ✗ return m_name;
325 }
326
327 /*! Return the shape of the quadrature rule */
328 20304 const RefShape& refShape() const {
329 20304 return *m_shape;
330 }
331
332 /*! Return the number of quadrature points */
333 1341704 inline int numQuadraturePoint(const std::size_t Index = 0) const {
334 1341704 return m_numQuadraturePoint[Index];
335 }
336
337 /*! quadraturePoint(ig) is the ig-th quadrature point */
338 30456 inline const QuadraturePoint* quadraturePoints() const {
339 30456 return m_point;
340 }
341
342 /*! quadraturePoint(ig) is the ig-th quadrature point */
343 5996791 inline const QuadraturePoint& quadraturePoint(const int ig) const {
344
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5996791 times.
 5996791 FEL_ASSERT( ig < m_numQuadraturePoint[0] );
345 5996791 return m_point[ ig ];
346 }
347
348 /*! weight(ig) is the ig-th quadrature weight */
349 230535313 inline double weight(const int ig) const {
350
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 230535313 times.
 230535313 FEL_ASSERT( ig < m_numQuadraturePoint[0] );
351 230535313 return m_point[ ig ].weight();
352 }
353
354 /*! weight(ig, jg) is the ig-th/jg-th quadrature weight */
355 inline double weight(const int ig, const int jg) const {
356 FEL_ASSERT( ig < m_numQuadraturePoint[1] );
357 FEL_ASSERT( jg < m_numQuadraturePoint[2] );
358 return m_point[ig * m_numQuadraturePoint[2] + jg].weight();
359 }
360
361 /*! \c quadraturePointCoor(ig) is the array of coordinates of the quadrature point ig */
362 inline const std::array<double, 3>& quadraturePointgetCoor(const int ig) const {
363 FEL_ASSERT( ig < m_numQuadraturePoint[0] );
364 return m_point[ ig ].getCoor();
365 }
366
367 /*! \c quadraturePointCoor(ig) is the array of coordinates of the quadrature point ig/jg */
368 inline const std::array<double, 3>& quadraturePointgetCoor(const int ig, const int jg) const {
369 FEL_ASSERT( ig < m_numQuadraturePoint[1] );
370 FEL_ASSERT( jg < m_numQuadraturePoint[2] );
371 return m_point[ig * m_numQuadraturePoint[2] + jg].getCoor();
372 }
373
374 /*! \c quadraturePointCoor(ig,icoor) is the coordinate icoor of the quadrature point ig */
375 2170355 inline double quadraturePointCoor(const int ig, const int icoor) const {
376
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2170355 times.
 2170355 FEL_ASSERT( ig < m_numQuadraturePoint[0] );
377 2170355 return m_point[ ig ].coor( icoor );
378 }
379
380 /*! \c quadraturePointCoor(ig,icoor) is the coordinate icoor of the quadrature point ig/jg */
381 inline double quadraturePointCoor(const int ig, const int jg, const int icoor) const {
382 FEL_ASSERT( ig < m_numQuadraturePoint[1] );
383 FEL_ASSERT( jg < m_numQuadraturePoint[2] );
384 return m_point[ig * m_numQuadraturePoint[2] + jg].coor(icoor);
385 }
386
387 /*! Return the degree of exactness of the quadrature rule */
388 179628 inline int degreeOfExactness(const std::size_t Index = 0) const {
389 179628 return m_degreeOfExactness[Index];
390 }
391
392 ///@}
393 ///@name Inquiry
394 ///@{
395
396 ///@}
397 ///@name Input and output
398 ///@{
399
400 void print(int verbose,std::ostream& c=std::cout) const;
401
402 ///@}
403 ///@name Friends
404 ///@{
405
406 ///@}
407 protected:
408 ///@name Protected static Member Variables
409 ///@{
410
411 ///@}
412 ///@name Protected member Variables
413 ///@{
414
415 const QuadraturePoint* m_point; /// List of quadrature points
416 const RefShape* m_shape; /// Geometric shape on which the quadrature rule can be used
417 std::string m_name = ""; /// Name of the quadrature rule
418 int* m_numQuadraturePoint= new int[1]; /// Number of quadrature points
419 int* m_degreeOfExactness= new int[1]; /// Degree of exactness of the quadrature rule
420
421 ///@}
422 ///@name Protected Operators
423 ///@{
424
425 ///@}
426 ///@name Protected Operations
427 ///@{
428
429 ///@}
430 ///@name Protected Access
431 ///@{
432
433 ///@}
434 ///@name Protected Inquiry
435 ///@{
436
437 ///@}
438 ///@name Protected LifeCycle
439 ///@{
440
441 ///@}
442 private:
443 ///@name Private static Member Variables
444 ///@{
445
446 ///@}
447 ///@name Private member Variables
448 ///@{
449
450 ///@}
451 ///@name Private Operations
452 ///@{
453
454 ///@}
455 ///@name Private Access
456 ///@{
457
458 ///@}
459 ///@name Private Inquiry
460 ///@{
461
462 ///@}
463 ///@name Private LifeCycle
464 ///@{
465
466 ///@}
467 };
468
469 /**
470 * @brief This is an enum defined in order to define the degree of exactness
471 */
472 enum DegreeOfExactness {
473 DegreeOfExactness_0 = 0,
474 DegreeOfExactness_1 = 1,
475 DegreeOfExactness_2 = 2,
476 DegreeOfExactness_3 = 3,
477 DegreeOfExactness_4 = 4,
478 DegreeOfExactness_5 = 5,
479 DegreeOfExactness_6 = 6,
480 DegreeOfExactness_7 = 7
481 };
482
483 /**
484 * @class ListOfQuadratureRule
485 * @brief List of all quadrature rules on a given geometrical shape
486 * @authors J-F. Gerbeau
487 * @date 02/2010
488 * This class contains the list of quadrature rules on a given geometrical shape.
489 * To add a new quadrature rule see the documentation of \c QuadratureRule.
490 * @see QuadratureRule
491 */
492 class ListOfQuadratureRule
493 {
494 public:
495
496 /*! Default constructor */
497 ListOfQuadratureRule(
498 const std::string name,
499 const int numQuadratureRule,
500 const QuadratureRule* quadratureRule
501 );
502
503 /*! Default constructor (combination) */
504 ListOfQuadratureRule(
505 const std::string name,
506 const ListOfQuadratureRule& listQuadratureRule1,
507 const ListOfQuadratureRule& listQuadratureRule2,
508 const double Coefficient = 1.0
509 );
510
511 /*! Default destructor */
512 5640 ~ListOfQuadratureRule() = default;
513
514 ///@}
515 ///@name Operators
516 ///@{
517
518 ///@}
519 ///@name Operations
520 ///@{
521
522 /**
523 * @brief Return the cheapest quadrature rule of the list ensuring that polynoms of degree \c deg are computed exactly.
524 * @details For the sake of clarity, it is recommended to use the enum \c DegreeOfExactness instead of \c int. Thus prefer for example deg=DegreeOfExactness_2 instead of deg=2.
525 */
526 const QuadratureRule& quadratureRuleByExactness(const DegreeOfExactness deg) const;
527
528 /**
529 * @brief Return the cheapest combined quadrature rule of the list ensuring that polynoms of degree \c deg are computed exactly.
530 * @details For the sake of clarity, it is recommended to use the enum \c DegreeOfExactness instead of \c int. Thus prefer for example deg=DegreeOfExactness_2 instead of deg=2.
531 */
532 const QuadratureRule& quadratureRuleByExactness(const DegreeOfExactness degi, const DegreeOfExactness degj) const;
533
534 ///@}
535 ///@name Access
536 ///@{
537
538 /**
539 * @brief Retrieves the maximum degree of exactness
540 * @return The maximum degree of exactness
541 */
542 2256 int getMaxDegreeOfExactness(const std::size_t Index = 0) const
543 {
544 2256 return m_maxDegreeOfExactness[Index];
545 }
546
547 /**
548 * @brief Retrieves the number of quadrature rules
549 * @return The number of quadrature rules
550 */
551 3384 int getNumQuadratureRule(const std::size_t Index = 0) const
552 {
553 3384 return m_numQuadratureRule[Index];
554 }
555
556 /**
557 * @brief Retrieves the quadrature rule
558 * @return The quadrature rule
559 */
560 3384 const QuadratureRule* getQuadratureRule() const
561 {
562 3384 return m_quadratureRule;
563 }
564
565 /**
566 * @brief Retrieves the quadrature rule by degree of exactness
567 * @return The quadrature rule by degree of exactness
568 */
569 2256 const int* getquadratureByDegreeOfExactness() const
570 {
571 2256 return m_quadratureByDegreeOfExactness;
572 }
573
574 ///@}
575 ///@name Inquiry
576 ///@{
577
578 ///@}
579 ///@name Input and output
580 ///@{
581
582 void print(int verbose,std::ostream& c=std::cout) const;
583
584 ///@}
585 ///@name Friends
586 ///@{
587
588 ///@}
589 protected:
590 ///@name Protected static Member Variables
591 ///@{
592
593 ///@}
594 ///@name Protected member Variables
595 ///@{
596
597 protected:
598 /// TODO: Most of these variables can be computed from m_quadratureRule, therefore reduce the quantity of memory consumed
599 std::string m_name; /// The name of the quadrature rule
600 int* m_numQuadratureRule = new int[1]; /// The number of quadrature rules
601 const QuadratureRule* m_quadratureRule; /// List of quadrature rule in strict increasing order of exactness (compulsory)
602 int* m_maxDegreeOfExactness = new int[1]; /// The maximum level of exactness
603 int* m_quadratureByDegreeOfExactness; /// m_quadratureByDegreeOfExactness[i] is a quadrature rule integrating exactly polynoms of degree i
604
605 ///@}
606 ///@name Protected Operators
607 ///@{
608
609 ///@}
610 ///@name Protected Operations
611 ///@{
612
613 ///@}
614 ///@name Protected Access
615 ///@{
616
617 ///@}
618 ///@name Protected Inquiry
619 ///@{
620
621 ///@}
622 ///@name Protected LifeCycle
623 ///@{
624
625 ///@}
626 private:
627 ///@name Private static Member Variables
628 ///@{
629
630 ///@}
631 ///@name Private member Variables
632 ///@{
633
634 ///@}
635 ///@name Private Operations
636 ///@{
637
638 ///@}
639 ///@name Private Access
640 ///@{
641
642 ///@}
643 ///@name Private Inquiry
644 ///@{
645
646 ///@}
647 ///@name Private LifeCycle
648 ///@{
649
650 ///@}
651 };
652
653 ///@}
654 ///@name Type Definitions
655 ///@{
656
657 extern const QuadratureRule quadratureRuleNULL;
658 extern const ListOfQuadratureRule listQuadratureRuleNULL;
659
660 extern const QuadratureRule quadratureRuleNode;
661 extern const ListOfQuadratureRule listQuadratureRuleNode;
662
663 extern const QuadratureRule quadratureRuleSeg1pt;
664 extern const QuadratureRule quadratureRuleSeg2pt;
665 extern const QuadratureRule quadratureRuleSeg3pt;
666 extern const ListOfQuadratureRule listQuadratureRuleSegment;
667
668 extern const QuadratureRule quadratureRuleTria1pt;
669 extern const QuadratureRule quadratureRuleTria3pt;
670 extern const QuadratureRule quadratureRuleMitc3;
671 extern const QuadratureRule quadratureRuleTria4pt;
672 extern const QuadratureRule quadratureRuleTria6pt;
673 extern const QuadratureRule quadratureRuleTria7pt;
674 extern const ListOfQuadratureRule listQuadratureRuleTriangle;
675
676 extern const QuadratureRule quadratureRuleQuad1pt;
677 extern const QuadratureRule quadratureRuleQuad4pt;
678 extern const QuadratureRule quadratureRuleQuad9pt;
679 extern const ListOfQuadratureRule listQuadratureRuleQuadrilateral;
680
681 extern const QuadratureRule quadratureRuleTetra1pt;
682 extern const QuadratureRule quadratureRuleTetra4pt;
683 extern const QuadratureRule quadratureRuleTetra5pt;
684 extern const QuadratureRule quadratureRuleTetra15pt;
685 extern const QuadratureRule quadratureRuleTetra64pt;
686 extern const ListOfQuadratureRule listQuadratureRuleTetrahedron;
687
688 extern const QuadratureRule quadratureRuleHexa1pt;
689 extern const QuadratureRule quadratureRuleHexa8pt;
690 extern const QuadratureRule quadratureRuleHexa27pt;
691 extern const ListOfQuadratureRule listQuadratureRuleHexahedron;
692 extern const ListOfQuadratureRule listQuadratureRuleHexahedronCombined;
693
694 extern const QuadratureRule quadratureRulePrism6pt;
695 extern const QuadratureRule quadratureRulePrism12ptComposite;
696 extern const QuadratureRule quadratureRulePrism21pt;
697 extern const QuadratureRule quadratureRulePrism42ptComposite;
698 extern const ListOfQuadratureRule listQuadratureRulePrism;
699 extern const ListOfQuadratureRule listQuadratureRulePrismCombined;
700
701 ///@}
702 } /* namespace felisce.*/
703 #endif /* QUADRATURE_RULE_HPP */
704