integral.h
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Header file for numerical integration routines based on quadrature
27
28// Include guards to prevent multiple inclusions of the header
29#ifndef OOMPH_INTEGRAL_HEADER
30#define OOMPH_INTEGRAL_HEADER
31
32
33// Config header generated by autoconfig
34#ifdef HAVE_CONFIG_H
35#include <oomph-lib-config.h>
36#endif
37
38// oomph-lib headers
39#include "oomph_utilities.h"
40#include "orthpoly.h"
41
42namespace oomph
43{
44 //=========================================================
45 /// Generic class for numerical integration schemes:
46 /// \f[ \int f(x_0,x_1...)\ dx_0 \ dx_1... = \sum_{i_{int}=0}^{\mbox{\tt nweight()} } f( x_j=\mbox{\tt knot($i_{int}$,j)} ) \ \ \ \mbox{\tt weight($i_{int}$)} \f]
47 //=========================================================
49 {
50 public:
51 /// Default constructor (empty)
53
54 /// Broken copy constructor
55 Integral(const Integral& dummy) = delete;
56
57 /// Broken assignment operator
58 void operator=(const Integral&) = delete;
59
60 /// Virtual destructor (empty)
61 virtual ~Integral(){};
62
63 /// Return the number of integration points of the scheme.
64 virtual unsigned nweight() const = 0;
65
66 /// Return local coordinate s[j] of i-th integration point.
67 virtual double knot(const unsigned& i, const unsigned& j) const = 0;
68
69 /// Return local coordinates of i-th intergration point. Broken virtual.
70 virtual Vector<double> knot(const unsigned& i) const
71 {
72 throw OomphLibError("Not implemented for this integration scheme (yet?).",
73 OOMPH_CURRENT_FUNCTION,
74 OOMPH_EXCEPTION_LOCATION);
75 }
76
77 /// Return weight of i-th integration point.
78 virtual double weight(const unsigned& i) const = 0;
79 };
80
81 //=============================================================================
82 /// Broken pseudo-integration scheme for points elements: Iit's not clear
83 /// in general what this integration scheme is supposed to. It probably
84 /// ought to evaluate integrals to zero but we're not sure in what
85 /// context this may be used. Replace by your own integration scheme
86 /// that does what you want!
87 //=============================================================================
88 class PointIntegral : public Integral
89 {
90 public:
91 /// Default constructor (empty)
93
94 /// Broken copy constructor
95 PointIntegral(const PointIntegral& dummy) = delete;
96
97 /// Broken assignment operator
98 void operator=(const PointIntegral&) = delete;
99
100 /// Number of integration points of the scheme
101 unsigned nweight() const
102 {
103 return 1;
104 }
105
106 /// Return coordinate s[j] (j=0) of integration point i --
107 /// deliberately broken!
108 double knot(const unsigned& i, const unsigned& j) const
109 {
110 throw OomphLibError("Local coordinate vector is of size zero, so this "
111 "should never be called.",
112 OOMPH_CURRENT_FUNCTION,
113 OOMPH_EXCEPTION_LOCATION);
114
115 // Dummy return
116 return 0.0;
117 }
118
119 /// Return weight of integration point i
120 double weight(const unsigned& i) const
121 {
122 return 1.0;
123 }
124 };
125
126 //=========================================================
127 /// Class for multidimensional Gaussian integration rules
128 ///
129 /// Empty -- just establishes the template parameters.
130 ///
131 /// General logic: The template parameters correspond to those
132 /// of the QElement family so that the integration scheme
133 /// Gauss<DIM,NNODE_1D> provides the default ("full") integration
134 /// scheme for QElement<DIM,NNODE_1D>. "Full" integration
135 /// means that for linear PDEs that are discretised on a uniform
136 /// mesh, all integrals arising in the Galerkin weak form of the PDE
137 /// are evaluated exactly. In such problems the highest-order
138 /// polynomials arise from the products of the undifferentiated
139 /// shape and test functions so a 4 node quad needs
140 /// an integration scheme that can integrate fourth-order
141 /// polynomials exactly etc.
142 //=========================================================
143 template<unsigned DIM, unsigned NPTS_1D>
144 class Gauss
145 {
146 };
147
148
149 //=========================================================
150 /// 1D Gaussian integration class.
151 /// Two integration points. This integration scheme can
152 /// integrate up to third-order polynomials exactly and
153 /// is therefore a suitable "full" integration scheme
154 /// for linear (two-node) elements in which the
155 /// highest-order polynomial is quadratic.
156 //=========================================================
157 template<>
158 class Gauss<1, 2> : public Integral
159 {
160 private:
161 /// Number of integration points in the scheme
162 static const unsigned Npts = 2;
163 /// Array to hold weights and knot points (defined in cc file)
164 static const double Knot[2][1], Weight[2];
165
166 public:
167 /// Default constructor (empty)
168 Gauss(){};
169
170 /// Broken copy constructor
171 Gauss(const Gauss& dummy) = delete;
172
173 /// Broken assignment operator
174 void operator=(const Gauss&) = delete;
175
176 /// Number of integration points of the scheme
177 unsigned nweight() const
178 {
179 return Npts;
180 }
181
182 /// Return coordinate s[j] (j=0) of integration point i
183 double knot(const unsigned& i, const unsigned& j) const
184 {
185 return Knot[i][j];
186 }
187
188 /// Return weight of integration point i
189 double weight(const unsigned& i) const
190 {
191 return Weight[i];
192 }
193 };
194
195 //=========================================================
196 /// 1D Gaussian integration class.
197 /// Three integration points. This integration scheme can
198 /// integrate up to fifth-order polynomials exactly and
199 /// is therefore a suitable "full" integration scheme
200 /// for quadratic (three-node) elements in which the
201 /// highest-order polynomial is fourth order.
202 //=========================================================
203 template<>
204 class Gauss<1, 3> : public Integral
205 {
206 private:
207 /// Number of integration points in the scheme
208 static const unsigned Npts = 3;
209 /// Array to hold weights and knot points (defined in cc file)
210 static const double Knot[3][1], Weight[3];
211
212 public:
213 /// Default constructor (empty)
214 Gauss(){};
215
216 /// Broken copy constructor
217 Gauss(const Gauss& dummy) = delete;
218
219 /// Broken assignment operator
220 void operator=(const Gauss&) = delete;
221
222 /// Number of integration points of the scheme
223 unsigned nweight() const
224 {
225 return Npts;
226 }
227
228 /// Return coordinate s[j] (j=0) of integration point i
229 double knot(const unsigned& i, const unsigned& j) const
230 {
231 return Knot[i][j];
232 }
233
234 /// Return weight of integration point i
235 double weight(const unsigned& i) const
236 {
237 return Weight[i];
238 }
239 };
240
241 //=========================================================
242 /// 1D Gaussian integration class
243 /// Four integration points. This integration scheme can
244 /// integrate up to seventh-order polynomials exactly and
245 /// is therefore a suitable "full" integration scheme
246 /// for cubic (four-node) elements in which the
247 /// highest-order polynomial is sixth order.
248 //=========================================================
249 template<>
250 class Gauss<1, 4> : public Integral
251 {
252 private:
253 /// Number of integration points in the scheme
254 static const unsigned Npts = 4;
255 /// Array to hold weight and knot points (defined in cc file)
256 static const double Knot[4][1], Weight[4];
257
258 public:
259 /// Default constructor (empty)
260 Gauss(){};
261
262 /// Broken copy constructor
263 Gauss(const Gauss& dummy) = delete;
264
265 /// Broken assignment operator
266 void operator=(const Gauss&) = delete;
267
268 /// Number of integration points of the scheme
269 unsigned nweight() const
270 {
271 return Npts;
272 }
273
274 /// Return coordinate x[j] (j=0) of integration point i
275 double knot(const unsigned& i, const unsigned& j) const
276 {
277 return Knot[i][j];
278 }
279
280 /// Return weight of integration point i
281 double weight(const unsigned& i) const
282 {
283 return Weight[i];
284 }
285 };
286
287 //=========================================================
288 /// 2D Gaussian integration class.
289 /// 2x2 integration points. This integration scheme can
290 /// integrate up to third-order polynomials exactly and
291 /// is therefore a suitable "full" integration scheme
292 /// for linear (four-node) elements in which the
293 /// highest-order polynomial is quadratic.
294 //=========================================================
295 template<>
296 class Gauss<2, 2> : public Integral
297 {
298 private:
299 /// Number of integration points in the scheme
300 static const unsigned Npts = 4;
301 /// Array to hold the weight and know points (defined in cc file)
302 static const double Knot[4][2], Weight[4];
303
304 public:
305 /// Default constructor (empty)
306 Gauss(){};
307
308 /// Broken copy constructor
309 Gauss(const Gauss& dummy) = delete;
310
311 /// Broken assignment operator
312 void operator=(const Gauss&) = delete;
313
314 /// Number of integration points of the scheme
315 unsigned nweight() const
316 {
317 return Npts;
318 }
319
320 /// Return coordinate x[j] of integration point i
321 double knot(const unsigned& i, const unsigned& j) const
322 {
323 return Knot[i][j];
324 }
325
326 /// Return weight of integration point i
327 double weight(const unsigned& i) const
328 {
329 return Weight[i];
330 }
331 };
332
333 //=========================================================
334 /// 2D Gaussian integration class.
335 /// 3x3 integration points. This integration scheme can
336 /// integrate up to fifth-order polynomials exactly and
337 /// is therefore a suitable "full" integration scheme
338 /// for quadratic (nine-node) elements in which the
339 /// highest-order polynomial is fourth order.
340 //=========================================================
341 template<>
342 class Gauss<2, 3> : public Integral
343 {
344 private:
345 /// Number of integration points in the scheme
346 static const unsigned Npts = 9;
347 /// Array to hold the weights and knots (defined in cc file)
348 static const double Knot[9][2], Weight[9];
349
350 public:
351 /// Default constructor (empty)
352 Gauss(){};
353
354 /// Broken copy constructor
355 Gauss(const Gauss& dummy) = delete;
356
357 /// Broken assignment operator
358 void operator=(const Gauss&) = delete;
359
360 /// Number of integration points of the scheme
361 unsigned nweight() const
362 {
363 return Npts;
364 }
365
366 /// Return coordinate s[j] of integration point i
367 double knot(const unsigned& i, const unsigned& j) const
368 {
369 return Knot[i][j];
370 }
371
372 /// Return weight of integration point i
373 double weight(const unsigned& i) const
374 {
375 return Weight[i];
376 }
377 };
378
379 //=========================================================
380 /// 2D Gaussian integration class.
381 /// 4x4 integration points. This integration scheme can
382 /// integrate up to seventh-order polynomials exactly and
383 /// is therefore a suitable "full" integration scheme
384 /// for cubic (sixteen-node) elements in which the
385 /// highest-order polynomial is sixth order.
386 //=========================================================
387 template<>
388 class Gauss<2, 4> : public Integral
389 {
390 private:
391 /// Number of integration points in the scheme
392 static const unsigned Npts = 16;
393 /// Array to hold the weights and knots (defined in cc file)
394 static const double Knot[16][2], Weight[16];
395
396
397 public:
398 /// Default constructor (empty)
399 Gauss(){};
400
401 /// Broken copy constructor
402 Gauss(const Gauss& dummy) = delete;
403
404 /// Broken assignment operator
405 void operator=(const Gauss&) = delete;
406
407 /// Number of integration points of the scheme
408 unsigned nweight() const
409 {
410 return Npts;
411 }
412
413 /// Return coordinate s[j] of integration point i
414 double knot(const unsigned& i, const unsigned& j) const
415 {
416 return Knot[i][j];
417 }
418
419 /// Return weight of integration point i
420 double weight(const unsigned& i) const
421 {
422 return Weight[i];
423 }
424 };
425
426 //=========================================================
427 /// 3D Gaussian integration class
428 /// 2x2x2 integration points. This integration scheme can
429 /// integrate up to third-order polynomials exactly and
430 /// is therefore a suitable "full" integration scheme
431 /// for linear (eight-node) elements in which the
432 /// highest-order polynomial is quadratic.
433 //=========================================================
434 template<>
435 class Gauss<3, 2> : public Integral
436 {
437 private:
438 /// Number of integration points in the scheme
439 static const unsigned Npts = 8;
440 /// Array to hold the weights and knots (defined in cc file)
441 static const double Knot[8][3], Weight[8];
442
443 public:
444 /// Default constructor (empty)
445 Gauss(){};
446
447 /// Broken copy constructor
448 Gauss(const Gauss& dummy) = delete;
449
450 /// Broken assignment operator
451 void operator=(const Gauss&) = delete;
452
453 /// Number of integration points of the scheme
454 unsigned nweight() const
455 {
456 return Npts;
457 }
458
459 /// Return coordinate s[j] of integration point i
460 double knot(const unsigned& i, const unsigned& j) const
461 {
462 return Knot[i][j];
463 }
464
465 /// Return weight of integration point i
466 double weight(const unsigned& i) const
467 {
468 return Weight[i];
469 }
470 };
471
472 //=========================================================
473 /// 3D Gaussian integration class
474 /// 3x3x3 integration points. This integration scheme can
475 /// integrate up to fifth-order polynomials exactly and
476 /// is therefore a suitable "full" integration scheme
477 /// for quadratic (27-node) elements in which the
478 /// highest-order polynomial is fourth order.
479 //=========================================================
480 template<>
481 class Gauss<3, 3> : public Integral
482 {
483 private:
484 /// Number of integration points in the scheme
485 static const unsigned Npts = 27;
486 /// Array to hold the weights and knots (defined in cc file)
487 static const double Knot[27][3], Weight[27];
488
489 public:
490 /// Default constructor (empty)
491 Gauss(){};
492
493 /// Broken copy constructor
494 Gauss(const Gauss& dummy) = delete;
495
496 /// Broken assignment operator
497 void operator=(const Gauss&) = delete;
498
499 /// Number of integration points of the scheme
500 unsigned nweight() const
501 {
502 return Npts;
503 }
504
505 /// Return coordinate x[j] of integration point i
506 double knot(const unsigned& i, const unsigned& j) const
507 {
508 return Knot[i][j];
509 }
510
511 /// Return weight of integration point i
512 double weight(const unsigned& i) const
513 {
514 return Weight[i];
515 }
516 };
517
518 //=========================================================
519 /// 3D Gaussian integration class.
520 /// 4x4x4 integration points. This integration scheme can
521 /// integrate up to seventh-order polynomials exactly and
522 /// is therefore a suitable "full" integration scheme
523 /// for cubic (64-node) elements in which the
524 /// highest-order polynomial is sixth order.
525 //=========================================================
526 template<>
527 class Gauss<3, 4> : public Integral
528 {
529 private:
530 /// Number of integration points in the scheme
531 static const unsigned Npts = 64;
532 /// Array to hold the weights and knots (defined in cc file)
533 static const double Knot[64][3], Weight[64];
534
535 public:
536 /// Default constructor (empty)
537 Gauss(){};
538
539 /// Broken copy constructor
540 Gauss(const Gauss& dummy) = delete;
541
542 /// Broken assignment operator
543 void operator=(const Gauss&) = delete;
544
545 /// Number of integration points of the scheme
546 unsigned nweight() const
547 {
548 return Npts;
549 }
550
551 /// Return coordinate x[j] of integration point i
552 double knot(const unsigned& i, const unsigned& j) const
553 {
554 return Knot[i][j];
555 }
556
557 /// Return weight of integration point i
558 double weight(const unsigned& i) const
559 {
560 return Weight[i];
561 }
562 };
563
564 //=========================================================
565 /// Class for multidimensional Gaussian integration rules,
566 /// over intervals other than -1 to 1, all intervals are
567 /// rescaled in this case
568 //=========================================================
569 template<unsigned DIM, unsigned NPTS_1D>
570 class Gauss_Rescaled : public Gauss<DIM, NPTS_1D>
571 {
572 private:
573 /// Store for the lower and upper limits of integration and the range
574 double Lower, Upper, Range;
575
576 public:
577 /// Default constructor (empty)
579
580 /// Broken copy constructor
581 Gauss_Rescaled(const Gauss_Rescaled& dummy) = delete;
582
583 /// Broken assignment operator
584 void operator=(const Gauss_Rescaled&) = delete;
585
586 /// The constructor in this case takes the lower and upper arguments
587 Gauss_Rescaled(double lower, double upper) : Lower(lower), Upper(upper)
588 {
589 // Set the range of integration
590 Range = upper - lower;
591 }
592
593 /// Return the rescaled knot values s[j] at integration point i
594 double knot(const unsigned& i, const unsigned& j) const
595 {
596 return (0.5 * (Gauss<DIM, NPTS_1D>::knot(i, j) * Range + Lower + Upper));
597 }
598
599 /// Return the rescaled weight at integration point i
600 double weight(const unsigned& i) const
601 {
603 pow(0.5 * Range, static_cast<int>(DIM));
604 }
605 };
606
607 //=========================================================
608 /// Class for Gaussian integration rules for triangles/tets.
609 ///
610 /// Empty -- just establishes the template parameters
611 ///
612 /// General logic: The template parameters correspond to those
613 /// of the TElement family so that the integration scheme
614 /// TGauss<DIM,NNODE_1D> provides the default ("full") integration
615 /// scheme for TElement<DIM,NNODE_1D>. "Full" integration
616 /// means that for linear PDEs that are discretised on a uniform
617 /// mesh, all integrals arising in the Galerkin weak form of the PDE
618 /// are evaluated exactly. In such problems the highest-order
619 /// polynomials arise from the products of the undifferentiated
620 /// shape and test functions so a three node triangle needs
621 /// an integration scheme that can integrate quadratic
622 /// polynomials exactly etc.
623 //=========================================================
624 template<unsigned DIM, unsigned NPTS_1D>
625 class TGauss
626 {
627 };
628
629
630 //=========================================================
631 /// 1D Gaussian integration class for linear "triangular" elements.
632 /// Two integration points. This integration scheme can
633 /// integrate up to second-order polynomials exactly and
634 /// is therefore a suitable "full" integration scheme
635 /// for linear (two-node) elements in which the
636 /// highest-order polynomial is quadratic.
637 //=========================================================
638 template<>
639 class TGauss<1, 2> : public Integral
640 {
641 private:
642 /// Number of integration points in the scheme
643 static const unsigned Npts = 2;
644
645 /// Array to hold the weights and knots (defined in cc file)
646 static const double Knot[2][1], Weight[2];
647
648 public:
649 /// Default constructor (empty)
651
652 /// Broken copy constructor
653 TGauss(const TGauss& dummy) = delete;
654
655 /// Broken assignment operator
656 void operator=(const TGauss&) = delete;
657
658 /// Number of integration points of the scheme
659 unsigned nweight() const
660 {
661 return Npts;
662 }
663
664 /// Return coordinate x[j] of integration point i
665 double knot(const unsigned& i, const unsigned& j) const
666 {
667 return Knot[i][j];
668 }
669
670 /// Return weight of integration point i
671 double weight(const unsigned& i) const
672 {
673 return Weight[i];
674 }
675 };
676
677 //=========================================================
678 /// 1D Gaussian integration class for quadratic "triangular" elements.
679 /// Three integration points. This integration scheme can
680 /// integrate up to fifth-order polynomials exactly and
681 /// is therefore a suitable "full" integration scheme
682 /// for quadratic (three-node) elements in which the
683 /// highest-order polynomial is fourth order.
684 //=========================================================
685 template<>
686 class TGauss<1, 3> : public Integral
687 {
688 private:
689 /// Number of integration points in the scheme
690 static const unsigned Npts = 3;
691 /// Array to hold the weights and knots (defined in cc file)
692 static const double Knot[3][1], Weight[3];
693
694 public:
695 /// Default constructor (empty)
697
698 /// Broken copy constructor
699 TGauss(const TGauss& dummy) = delete;
700
701 /// Broken assignment operator
702 void operator=(const TGauss&) = delete;
703
704 /// Number of integration points of the scheme
705 unsigned nweight() const
706 {
707 return Npts;
708 }
709
710 /// Return coordinate x[j] of integration point i
711 double knot(const unsigned& i, const unsigned& j) const
712 {
713 return Knot[i][j];
714 }
715
716 /// Return weight of integration point i
717 double weight(const unsigned& i) const
718 {
719 return Weight[i];
720 }
721 };
722
723
724 //=========================================================
725 /// 1D Gaussian integration class for cubic "triangular" elements.
726 /// Four integration points. This integration scheme can
727 /// integrate up to seventh-order polynomials exactly and
728 /// is therefore a suitable "full" integration scheme
729 /// for cubic (ten-node) elements in which the
730 /// highest-order polynomial is sixth order.
731 //=========================================================
732 template<>
733 class TGauss<1, 4> : public Integral
734 {
735 private:
736 /// Number of integration points in the scheme
737 static const unsigned Npts = 4;
738 /// Array to hold the weights and knots (defined in cc file)
739 static const double Knot[4][1], Weight[4];
740
741 public:
742 /// Default constructor (empty)
744
745 /// Broken copy constructor
746 TGauss(const TGauss& dummy) = delete;
747
748 /// Broken assignment operator
749 void operator=(const TGauss&) = delete;
750
751 /// Number of integration points of the scheme
752 unsigned nweight() const
753 {
754 return Npts;
755 }
756
757 /// Return coordinate x[j] of integration point i
758 double knot(const unsigned& i, const unsigned& j) const
759 {
760 return Knot[i][j];
761 }
762
763 /// Return weight of integration point i
764 double weight(const unsigned& i) const
765 {
766 return Weight[i];
767 }
768 };
769
770 //=========================================================
771 template<>
772 class TGauss<1, 5> : public Integral
773 {
774 private:
775 /// Number of integration points in the scheme
776 static const unsigned Npts = 5;
777 /// Array to hold the weights and knots (defined in cc file)
778 static const double Knot[5][1], Weight[5];
779
780 public:
781 /// Default constructor (empty)
783
784 /// Broken copy constructor
785 TGauss(const TGauss& dummy) = delete;
786
787 /// Broken assignment operator
788 void operator=(const TGauss&) = delete;
789
790 /// Number of integration points of the scheme
791 unsigned nweight() const
792 {
793 return Npts;
794 }
795
796 /// Return coordinate x[j] of integration point i
797 double knot(const unsigned& i, const unsigned& j) const
798 {
799 return Knot[i][j];
800 }
801
802 /// Return weight of integration point i
803 double weight(const unsigned& i) const
804 {
805 return Weight[i];
806 }
807 };
808
809
810 //=========================================================
811 /// 2D Gaussian integration class for linear triangles.
812 /// Three integration points. This integration scheme can
813 /// integrate up to second-order polynomials exactly and
814 /// is therefore a suitable "full" integration scheme
815 /// for linear (three-node) elements in which the
816 /// highest-order polynomial is quadratic.
817 //=========================================================
818 template<>
819 class TGauss<2, 2> : public Integral
820 {
821 private:
822 /// Number of integration points in the scheme
823 static const unsigned Npts = 3;
824
825 /// Array to hold the weights and knots (defined in cc file)
826 static const double Knot[3][2], Weight[3];
827
828 public:
829 /// Default constructor (empty)
831
832 /// Broken copy constructor
833 TGauss(const TGauss& dummy) = delete;
834
835 /// Broken assignment operator
836 void operator=(const TGauss&) = delete;
837
838 /// Number of integration points of the scheme
839 unsigned nweight() const
840 {
841 return Npts;
842 }
843
844 /// Return coordinate x[j] of integration point i
845 double knot(const unsigned& i, const unsigned& j) const
846 {
847 return Knot[i][j];
848 }
849
850 /// Return weight of integration point i
851 double weight(const unsigned& i) const
852 {
853 return Weight[i];
854 }
855 };
856
857 //=========================================================
858 /// 2D Gaussian integration class for quadratic triangles.
859 /// Seven integration points. This integration scheme can
860 /// integrate up to fifth-order polynomials exactly and
861 /// is therefore a suitable "full" integration scheme
862 /// for quadratic (six-node) elements in which the
863 /// highest-order polynomial is fourth order.
864 //=========================================================
865 template<>
866 class TGauss<2, 3> : public Integral
867 {
868 private:
869 /// Number of integration points in the scheme
870 static const unsigned Npts = 7;
871 /// Array to hold the weights and knots (defined in cc file)
872 static const double Knot[7][2], Weight[7];
873
874 public:
875 /// Default constructor (empty)
877
878 /// Broken copy constructor
879 TGauss(const TGauss& dummy) = delete;
880
881 /// Broken assignment operator
882 void operator=(const TGauss&) = delete;
883
884 /// Number of integration points of the scheme
885 unsigned nweight() const
886 {
887 return Npts;
888 }
889
890 /// Return coordinate x[j] of integration point i
891 double knot(const unsigned& i, const unsigned& j) const
892 {
893 return Knot[i][j];
894 }
895
896 /// Return weight of integration point i
897 double weight(const unsigned& i) const
898 {
899 return Weight[i];
900 }
901 };
902
903
904 //=========================================================
905 /// 2D Gaussian integration class for cubic triangles.
906 /// Thirteen integration points. This integration scheme can
907 /// integrate up to seventh-order polynomials exactly and
908 /// is therefore a suitable "full" integration scheme
909 /// for cubic (ten-node) elements in which the
910 /// highest-order polynomial is sixth order.
911 //=========================================================
912 template<>
913 class TGauss<2, 4> : public Integral
914 {
915 private:
916 /// Number of integration points in the scheme
917 static const unsigned Npts = 13;
918 /// Array to hold the weights and knots (defined in cc file)
919 static const double Knot[13][2], Weight[13];
920
921 public:
922 /// Default constructor (empty)
924
925 /// Broken copy constructor
926 TGauss(const TGauss& dummy) = delete;
927
928 /// Broken assignment operator
929 void operator=(const TGauss&) = delete;
930
931 /// Number of integration points of the scheme
932 unsigned nweight() const
933 {
934 return Npts;
935 }
936
937 /// Return coordinate x[j] of integration point i
938 double knot(const unsigned& i, const unsigned& j) const
939 {
940 return Knot[i][j];
941 }
942
943 /// Return weight of integration point i
944 double weight(const unsigned& i) const
945 {
946 return Weight[i];
947 }
948 };
949
950 //------------------------------------------------------------
951 //"Full integration" weights for 2D triangles
952 // accurate up to order 11
953 // http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_tri/
954 // quadrature_rules_tri.html
955 //------------------------------------------------------------
956 template<>
957 class TGauss<2, 13> : public Integral
958 {
959 private:
960 /// Number of integration points in the scheme
961 static const unsigned Npts = 37;
962 /// Array to hold the weights and knots (defined in cc file)
963 static const double Knot[37][2], Weight[37];
964
965 public:
966 /// Default constructor (empty)
968
969 /// Broken copy constructor
970 TGauss(const TGauss& dummy) = delete;
971
972 /// Broken assignment operator
973 void operator=(const TGauss&) = delete;
974
975 /// Number of integration points of the scheme
976 unsigned nweight() const
977 {
978 return Npts;
979 }
980
981 /// Return coordinate x[j] of integration point i
982 double knot(const unsigned& i, const unsigned& j) const
983 {
984 return Knot[i][j];
985 }
986
987 /// Return weight of integration point i
988 double weight(const unsigned& i) const
989 {
990 return Weight[i];
991 }
992 };
993
994 //------------------------------------------------------------
995 //"Full integration" weights for 2D triangles
996 // accurate up to points 19, degree of precision 8, a rule from ACM TOMS
997 // algorithm #584.
998 // http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_tri/quadrature_rules_tri.html
999 // TOMS589_19 a 19 point Gauss rule accurate to order 8
1000 //------------------------------------------------------------
1001 template<>
1002 class TGauss<2, 9> : public Integral
1003 {
1004 private:
1005 /// Number of integration points in the scheme
1006 static const unsigned Npts = 19;
1007 /// Array to hold the weights and knots (defined in cc file)
1008 static const double Knot[19][2], Weight[19];
1009
1010 public:
1011 /// Default constructor (empty)
1013
1014 /// Broken copy constructor
1015 TGauss(const TGauss& dummy) = delete;
1016
1017 /// Broken assignment operator
1018 void operator=(const TGauss&) = delete;
1019
1020 /// Number of integration points of the scheme
1021 unsigned nweight() const
1022 {
1023 return Npts;
1024 }
1025
1026 /// Return coordinate x[j] of integration point i
1027 double knot(const unsigned& i, const unsigned& j) const
1028 {
1029 return Knot[i][j];
1030 }
1031
1032 /// Return weight of integration point i
1033 double weight(const unsigned& i) const
1034 {
1035 return Weight[i];
1036 }
1037 };
1038
1039 //------------------------------------------------------------
1040 //"Full integration" weights for 2D triangles
1041 // accurate up to order 16
1042 // http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_tri/
1043 // quadrature_rules_tri.html
1044 // This uses the order 16 Dunavant rule, computed using software available
1045 // from
1046 // https://people.sc.fsu.edu/~jburkardt/cpp_src/triangle_dunavant_rule/triangle_dunavant_rule.html
1047 // This integegration scheme is accurate up to order 16 and uses 52 points
1048 //------------------------------------------------------------
1049 template<>
1050 class TGauss<2, 16> : public Integral
1051 {
1052 private:
1053 /// Number of integration points in the scheme
1054 static const unsigned Npts = 52;
1055 /// Array to hold the weights and knots (defined in cc file)
1056 static const double Knot[52][2], Weight[52];
1057
1058 public:
1059 /// Default constructor (empty)
1061
1062 /// Broken copy constructor
1063 TGauss(const TGauss& dummy) = delete;
1064
1065 /// Broken assignment operator
1066 void operator=(const TGauss&) = delete;
1067
1068 /// Number of integration points of the scheme
1069 unsigned nweight() const
1070 {
1071 return Npts;
1072 }
1073
1074 /// Return coordinate x[j] of integration point i
1075 double knot(const unsigned& i, const unsigned& j) const
1076 {
1077 return Knot[i][j];
1078 }
1079
1080 /// Return weight of integration point i
1081 double weight(const unsigned& i) const
1082 {
1083 return Weight[i];
1084 }
1085 };
1086
1087 //------------------------------------------------------------
1088 //"Full integration" weights for 2D triangles
1089 // accurate up to order 15
1090 // http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_tri/quadrature_rules_tri.html
1091 //------------------------------------------------------------
1092
1093 template<>
1094 class TGauss<2, 5> : public Integral
1095 {
1096 private:
1097 /// Number of integration points in the scheme
1098 static const unsigned Npts = 64;
1099 /// Array to hold the weights and knots (defined in cc file)
1100 static const double Knot[64][2], Weight[64];
1101
1102 public:
1103 /// Default constructor (empty)
1105
1106 /// Broken copy constructor
1107 TGauss(const TGauss& dummy) = delete;
1108
1109 /// Broken assignment operator
1110 void operator=(const TGauss&) = delete;
1111
1112 /// Number of integration points of the scheme
1113 unsigned nweight() const
1114 {
1115 return Npts;
1116 }
1117
1118 /// Return coordinate x[j] of integration point i
1119 double knot(const unsigned& i, const unsigned& j) const
1120 {
1121 return Knot[i][j];
1122 }
1123
1124 /// Return weight of integration point i
1125 double weight(const unsigned& i) const
1126 {
1127 return Weight[i];
1128 }
1129 };
1130
1131 //=========================================================
1132 /// 3D Gaussian integration class for tets.
1133 /// Four integration points. This integration scheme can
1134 /// integrate up to second-order polynomials exactly and
1135 /// is therefore a suitable "full" integration scheme
1136 /// for linear (four-node) elements in which the
1137 /// highest-order polynomial is quadratic.
1138 //=========================================================
1139 template<>
1140 class TGauss<3, 2> : public Integral
1141 {
1142 private:
1143 /// Number of integration points in the scheme
1144 static const unsigned Npts = 4;
1145 /// Array to hold the weights and knots (defined in cc file)
1146 static const double Knot[4][3], Weight[4];
1147
1148 public:
1149 /// Default constructor (empty)
1151
1152 /// Broken copy constructor
1153 TGauss(const TGauss& dummy) = delete;
1154
1155 /// Broken assignment operator
1156 void operator=(const TGauss&) = delete;
1157
1158 /// Number of integration points of the scheme
1159 unsigned nweight() const
1160 {
1161 return Npts;
1162 }
1163
1164 /// Return coordinate x[j] of integration point i
1165 double knot(const unsigned& i, const unsigned& j) const
1166 {
1167 return Knot[i][j];
1168 }
1169
1170 /// Return weight of integration point i
1171 double weight(const unsigned& i) const
1172 {
1173 return Weight[i];
1174 }
1175 };
1176
1177
1178 //=========================================================
1179 /// 3D Gaussian integration class for tets.
1180 /// Eleven integration points. This integration scheme can
1181 /// integrate up to fourth-order polynomials exactly and
1182 /// is therefore a suitable "full" integration scheme
1183 /// for quadratic (ten-node) elements in which the
1184 /// highest-order polynomial is fourth order.
1185 /// The numbers are from Keast CMAME 55 pp339-348 (1986)
1186 //=========================================================
1187 template<>
1188 class TGauss<3, 3> : public Integral
1189 {
1190 private:
1191 /// Number of integration points in the scheme
1192 static const unsigned Npts = 11;
1193 /// Array to hold the weights and knots (defined in cc file)
1194 static const double Knot[11][3], Weight[11];
1195
1196 public:
1197 /// Default constructor (empty)
1199
1200 /// Broken copy constructor
1201 TGauss(const TGauss& dummy) = delete;
1202
1203 /// Broken assignment operator
1204 void operator=(const TGauss&) = delete;
1205
1206 /// Number of integration points of the scheme
1207 unsigned nweight() const
1208 {
1209 return Npts;
1210 }
1211
1212 /// Return coordinate x[j] of integration point i
1213 double knot(const unsigned& i, const unsigned& j) const
1214 {
1215 return Knot[i][j];
1216 }
1217
1218 /// Return weight of integration point i
1219 double weight(const unsigned& i) const
1220 {
1221 return Weight[i];
1222 }
1223 };
1224
1225
1226 //=========================================================
1227 /// 3D Gaussian integration class for tets.
1228 /// 45 integration points. This integration scheme can
1229 /// integrate up to eighth-order polynomials exactly and
1230 /// is therefore a suitable "full" integration scheme
1231 /// for quartic elements in which the
1232 /// highest-order polynomial is fourth order.
1233 /// The numbers are from Keast CMAME 55 pp339-348 (1986)
1234 //=========================================================
1235 template<>
1236 class TGauss<3, 5> : public Integral
1237 {
1238 private:
1239 /// Number of integration points in the scheme
1240 static const unsigned Npts = 45;
1241 /// Array to hold the weights and knots (defined in cc file)
1242 static const double Knot[45][3], Weight[45];
1243
1244 public:
1245 /// Default constructor (empty)
1247
1248 /// Broken copy constructor
1249 TGauss(const TGauss& dummy) = delete;
1250
1251 /// Broken assignment operator
1252 void operator=(const TGauss&) = delete;
1253
1254 /// Number of integration points of the scheme
1255 unsigned nweight() const
1256 {
1257 return Npts;
1258 }
1259
1260 /// Return coordinate x[j] of integration point i
1261 double knot(const unsigned& i, const unsigned& j) const
1262 {
1263 return Knot[i][j];
1264 }
1265
1266 /// Return weight of integration point i
1267 double weight(const unsigned& i) const
1268 {
1269 return Weight[i];
1270 }
1271 };
1272
1273
1274 //===================================================================
1275 /// Class for multidimensional Gauss Lobatto Legendre integration
1276 /// rules
1277 /// empty - just establishes template parameters
1278 //===================================================================
1279 template<unsigned DIM, unsigned NPTS_1D>
1281 {
1282 };
1283
1284 //===================================================================
1285 /// 1D Gauss Lobatto Legendre integration class
1286 //===================================================================
1287 template<unsigned NPTS_1D>
1288 class GaussLobattoLegendre<1, NPTS_1D> : public Integral
1289 {
1290 private:
1291 /// Number of integration points in scheme
1292 static const unsigned Npts = NPTS_1D;
1293 /// Array to hold weight and knot points
1294 // These are not constant, because they are calculated on the fly
1295 double Knot[NPTS_1D][1], Weight[NPTS_1D];
1296
1297 public:
1298 /// Deafault constructor. Calculates and stores GLL nodes
1300
1301 /// Number of integration points of the scheme
1302 unsigned nweight() const
1303 {
1304 return Npts;
1305 }
1306
1307 /// Return coordinate s[j] (j=0) of integration point i
1308 double knot(const unsigned& i, const unsigned& j) const
1309 {
1310 return Knot[i][j];
1311 }
1312
1313 /// Return weight of integration point i
1314 double weight(const unsigned& i) const
1315 {
1316 return Weight[i];
1317 }
1318 };
1319
1320
1321 //=============================================================
1322 /// Calculate positions and weights for the 1D Gauss Lobatto
1323 /// Legendre integration class
1324 //=============================================================
1325 template<unsigned NPTS_1D>
1327 {
1328 // Temporary storage for the integration points
1329 Vector<double> s(NPTS_1D), w(NPTS_1D);
1330 // call the function that calculates the points
1331 Orthpoly::gll_nodes(NPTS_1D, s, w);
1332 // Populate the arrays
1333 for (unsigned i = 0; i < NPTS_1D; i++)
1334 {
1335 Knot[i][0] = s[i];
1336 Weight[i] = w[i];
1337 }
1338 }
1339
1340
1341 //===================================================================
1342 /// 2D Gauss Lobatto Legendre integration class
1343 //===================================================================
1344 template<unsigned NPTS_1D>
1345 class GaussLobattoLegendre<2, NPTS_1D> : public Integral
1346 {
1347 private:
1348 /// Number of integration points in scheme
1349 static const unsigned long int Npts = NPTS_1D * NPTS_1D;
1350
1351 /// Array to hold weight and knot points
1352 double Knot[NPTS_1D * NPTS_1D][2],
1353 Weight[NPTS_1D * NPTS_1D]; // COULDN'T MAKE THESE
1354 // const BECAUSE THEY ARE CALCULATED (at least currently)
1355
1356 public:
1357 /// Deafault constructor. Calculates and stores GLL nodes
1359
1360 /// Number of integration points of the scheme
1361 unsigned nweight() const
1362 {
1363 return Npts;
1364 }
1365
1366 /// Return coordinate s[j] (j=0) of integration point i
1367 double knot(const unsigned& i, const unsigned& j) const
1368 {
1369 return Knot[i][j];
1370 }
1371
1372 /// Return weight of integration point i
1373 double weight(const unsigned& i) const
1374 {
1375 return Weight[i];
1376 }
1377 };
1378
1379 //=============================================================
1380 /// Calculate positions and weights for the 2D Gauss Lobatto
1381 /// Legendre integration class
1382 //=============================================================
1383
1384 template<unsigned NPTS_1D>
1386 {
1387 // Tempoarary storage for the 1D knots and weights
1388 Vector<double> s(NPTS_1D), w(NPTS_1D);
1389 // Call the function to populate the array
1390 Orthpoly::gll_nodes(NPTS_1D, s, w);
1391 int i_fast = 0, i_slow = 0;
1392 for (unsigned i = 0; i < NPTS_1D * NPTS_1D; i++)
1393 {
1394 if (i_fast == NPTS_1D)
1395 {
1396 i_fast = 0;
1397 i_slow++;
1398 }
1399 Knot[i][0] = s[i_fast];
1400 Knot[i][1] = s[i_slow];
1401 Weight[i] = w[i_fast] * w[i_slow];
1402 i_fast++;
1403 }
1404 }
1405
1406
1407 //===================================================================
1408 /// 3D Gauss Lobatto Legendre integration class
1409 //===================================================================
1410 template<unsigned NPTS_1D>
1411 class GaussLobattoLegendre<3, NPTS_1D> : public Integral
1412 {
1413 private:
1414 /// Number of integration points in scheme
1415 static const unsigned long int Npts = NPTS_1D * NPTS_1D * NPTS_1D;
1416
1417 /// Array to hold weight and knot points
1418 double Knot[NPTS_1D * NPTS_1D * NPTS_1D][3],
1419 Weight[NPTS_1D * NPTS_1D * NPTS_1D]; // COULDN'T MAKE THESE
1420 // const BECAUSE THEY ARE CALCULATED (at least currently)
1421
1422 public:
1423 /// Deafault constructor. Calculates and stores GLL nodes
1425
1426 /// Number of integration points of the scheme
1427 unsigned nweight() const
1428 {
1429 return Npts;
1430 }
1431
1432 /// Return coordinate s[j] (j=0) of integration point i
1433 double knot(const unsigned& i, const unsigned& j) const
1434 {
1435 return Knot[i][j];
1436 }
1437
1438 /// Return weight of integration point i
1439 double weight(const unsigned& i) const
1440 {
1441 return Weight[i];
1442 }
1443 };
1444
1445 //=============================================================
1446 /// Calculate positions and weights for the 3D Gauss Lobatto
1447 /// Legendre integration class
1448 //=============================================================
1449
1450 template<unsigned NPTS_1D>
1452 {
1453 // Tempoarary storage for the 1D knots and weights
1454 Vector<double> s(NPTS_1D), w(NPTS_1D);
1455 // Call the function to populate the array
1456 Orthpoly::gll_nodes(NPTS_1D, s, w);
1457 for (unsigned k = 0; k < NPTS_1D; k++)
1458 {
1459 for (unsigned j = 0; j < NPTS_1D; j++)
1460 {
1461 for (unsigned i = 0; i < NPTS_1D; i++)
1462 {
1463 unsigned index = NPTS_1D * NPTS_1D * k + NPTS_1D * j + i;
1464 Knot[index][0] = s[i];
1465 Knot[index][1] = s[j];
1466 Knot[index][2] = s[k];
1467 Weight[index] = w[i] * w[j] * w[k];
1468 }
1469 }
1470 }
1471 }
1472
1473 /// /////////////////////////////////////////////////////////////////
1474 /// /////////////////////////////////////////////////////////////////
1475 /// /////////////////////////////////////////////////////////////////
1476
1477
1478 //===================================================================
1479 /// Class for multidimensional Gauss Legendre integration
1480 /// rules
1481 /// empty - just establishes template parameters
1482 //===================================================================
1483 template<unsigned DIM, unsigned NPTS_1D>
1485 {
1486 };
1487
1488 //===================================================================
1489 /// 1D Gauss Legendre integration class
1490 //===================================================================
1491 template<unsigned NPTS_1D>
1492 class GaussLegendre<1, NPTS_1D> : public Integral
1493 {
1494 private:
1495 /// Number of integration points in scheme
1496 static const unsigned Npts = NPTS_1D;
1497 /// Array to hold weight and knot points
1498 // These are not constant, because they are calculated on the fly
1499 double Knot[NPTS_1D][1], Weight[NPTS_1D];
1500
1501 public:
1502 /// Deafault constructor. Calculates and stores GLL nodes
1503 GaussLegendre();
1504
1505 /// Number of integration points of the scheme
1506 unsigned nweight() const
1507 {
1508 return Npts;
1509 }
1510
1511 /// Return coordinate s[j] (j=0) of integration point i
1512 double knot(const unsigned& i, const unsigned& j) const
1513 {
1514 return Knot[i][j];
1515 }
1516
1517 /// Return weight of integration point i
1518 double weight(const unsigned& i) const
1519 {
1520 return Weight[i];
1521 }
1522 };
1523
1524
1525 //=============================================================
1526 /// Calculate positions and weights for the 1D Gauss
1527 /// Legendre integration class
1528 //=============================================================
1529 template<unsigned NPTS_1D>
1531 {
1532 // Temporary storage for the integration points
1533 Vector<double> s(NPTS_1D), w(NPTS_1D);
1534 // call the function that calculates the points
1535 Orthpoly::gl_nodes(NPTS_1D, s, w);
1536 // Populate the arrays
1537 for (unsigned i = 0; i < NPTS_1D; i++)
1538 {
1539 Knot[i][0] = s[i];
1540 Weight[i] = w[i];
1541 }
1542 }
1543
1544
1545 //===================================================================
1546 /// 2D Gauss Legendre integration class
1547 //===================================================================
1548 template<unsigned NPTS_1D>
1549 class GaussLegendre<2, NPTS_1D> : public Integral
1550 {
1551 private:
1552 /// Number of integration points in scheme
1553 static const unsigned long int Npts = NPTS_1D * NPTS_1D;
1554
1555 /// Array to hold weight and knot points
1556 double Knot[NPTS_1D * NPTS_1D][2],
1557 Weight[NPTS_1D * NPTS_1D]; // COULDN'T MAKE THESE
1558 // const BECAUSE THEY ARE CALCULATED (at least currently)
1559
1560 public:
1561 /// Deafault constructor. Calculates and stores GLL nodes
1562 GaussLegendre();
1563
1564 /// Number of integration points of the scheme
1565 unsigned nweight() const
1566 {
1567 return Npts;
1568 }
1569
1570 /// Return coordinate s[j] (j=0) of integration point i
1571 double knot(const unsigned& i, const unsigned& j) const
1572 {
1573 return Knot[i][j];
1574 }
1575
1576 /// Return weight of integration point i
1577 double weight(const unsigned& i) const
1578 {
1579 return Weight[i];
1580 }
1581 };
1582
1583 //=============================================================
1584 /// Calculate positions and weights for the 2D Gauss
1585 /// Legendre integration class
1586 //=============================================================
1587
1588 template<unsigned NPTS_1D>
1590 {
1591 // Tempoarary storage for the 1D knots and weights
1592 Vector<double> s(NPTS_1D), w(NPTS_1D);
1593 // Call the function to populate the array
1594 Orthpoly::gl_nodes(NPTS_1D, s, w);
1595 int i_fast = 0, i_slow = 0;
1596 for (unsigned i = 0; i < NPTS_1D * NPTS_1D; i++)
1597 {
1598 if (i_fast == NPTS_1D)
1599 {
1600 i_fast = 0;
1601 i_slow++;
1602 }
1603 Knot[i][0] = s[i_fast];
1604 Knot[i][1] = s[i_slow];
1605 Weight[i] = w[i_fast] * w[i_slow];
1606 i_fast++;
1607 }
1608 }
1609
1610
1611 //===================================================================
1612 /// 3D Gauss Legendre integration class
1613 //===================================================================
1614 template<unsigned NPTS_1D>
1615 class GaussLegendre<3, NPTS_1D> : public Integral
1616 {
1617 private:
1618 /// Number of integration points in scheme
1619 static const unsigned long int Npts = NPTS_1D * NPTS_1D * NPTS_1D;
1620
1621 /// Array to hold weight and knot points
1622 double Knot[NPTS_1D * NPTS_1D * NPTS_1D][3],
1623 Weight[NPTS_1D * NPTS_1D * NPTS_1D]; // COULDN'T MAKE THESE
1624 // const BECAUSE THEY ARE CALCULATED (at least currently)
1625
1626 public:
1627 /// Deafault constructor. Calculates and stores GLL nodes
1628 GaussLegendre();
1629
1630 /// Number of integration points of the scheme
1631 unsigned nweight() const
1632 {
1633 return Npts;
1634 }
1635
1636 /// Return coordinate s[j] (j=0) of integration point i
1637 double knot(const unsigned& i, const unsigned& j) const
1638 {
1639 return Knot[i][j];
1640 }
1641
1642 /// Return weight of integration point i
1643 double weight(const unsigned& i) const
1644 {
1645 return Weight[i];
1646 }
1647 };
1648
1649 //=============================================================
1650 /// Calculate positions and weights for the 3D Gauss
1651 /// Legendre integration class
1652 //=============================================================
1653
1654 template<unsigned NPTS_1D>
1656 {
1657 // Tempoarary storage for the 1D knots and weights
1658 Vector<double> s(NPTS_1D), w(NPTS_1D);
1659 // Call the function to populate the array
1660 Orthpoly::gl_nodes(NPTS_1D, s, w);
1661 for (unsigned k = 0; k < NPTS_1D; k++)
1662 {
1663 for (unsigned j = 0; j < NPTS_1D; j++)
1664 {
1665 for (unsigned i = 0; i < NPTS_1D; i++)
1666 {
1667 unsigned index = NPTS_1D * NPTS_1D * k + NPTS_1D * j + i;
1668 Knot[index][0] = s[i];
1669 Knot[index][1] = s[j];
1670 Knot[index][2] = s[k];
1671 Weight[index] = w[i] * w[j] * w[k];
1672 }
1673 }
1674 }
1675 }
1676
1677
1678} // namespace oomph
1679
1680#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1518
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] (j=0) of integration point i.
Definition: integral.h:1512
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1506
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1565
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1577
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] (j=0) of integration point i.
Definition: integral.h:1571
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1631
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1643
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] (j=0) of integration point i.
Definition: integral.h:1637
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
Definition: integral.h:1485
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1314
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1302
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] (j=0) of integration point i.
Definition: integral.h:1308
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1361
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1373
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] (j=0) of integration point i.
Definition: integral.h:1367
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1439
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] (j=0) of integration point i.
Definition: integral.h:1433
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1427
Class for multidimensional Gauss Lobatto Legendre integration rules empty - just establishes template...
Definition: integral.h:1281
Gauss()
Default constructor (empty)
Definition: integral.h:168
void operator=(const Gauss &)=delete
Broken assignment operator.
Gauss(const Gauss &dummy)=delete
Broken copy constructor.
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:177
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] (j=0) of integration point i.
Definition: integral.h:183
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:189
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:235
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] (j=0) of integration point i.
Definition: integral.h:229
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:223
Gauss()
Default constructor (empty)
Definition: integral.h:214
void operator=(const Gauss &)=delete
Broken assignment operator.
Gauss(const Gauss &dummy)=delete
Broken copy constructor.
void operator=(const Gauss &)=delete
Broken assignment operator.
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:269
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:281
Gauss()
Default constructor (empty)
Definition: integral.h:260
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] (j=0) of integration point i.
Definition: integral.h:275
Gauss(const Gauss &dummy)=delete
Broken copy constructor.
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:321
Gauss(const Gauss &dummy)=delete
Broken copy constructor.
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:327
void operator=(const Gauss &)=delete
Broken assignment operator.
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:315
Gauss()
Default constructor (empty)
Definition: integral.h:306
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] of integration point i.
Definition: integral.h:367
void operator=(const Gauss &)=delete
Broken assignment operator.
Gauss(const Gauss &dummy)=delete
Broken copy constructor.
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:373
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:361
Gauss()
Default constructor (empty)
Definition: integral.h:352
Gauss(const Gauss &dummy)=delete
Broken copy constructor.
Gauss()
Default constructor (empty)
Definition: integral.h:399
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] of integration point i.
Definition: integral.h:414
void operator=(const Gauss &)=delete
Broken assignment operator.
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:408
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:420
Gauss(const Gauss &dummy)=delete
Broken copy constructor.
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:454
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] of integration point i.
Definition: integral.h:460
void operator=(const Gauss &)=delete
Broken assignment operator.
Gauss()
Default constructor (empty)
Definition: integral.h:445
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:466
void operator=(const Gauss &)=delete
Broken assignment operator.
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:512
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:506
Gauss()
Default constructor (empty)
Definition: integral.h:491
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:500
Gauss(const Gauss &dummy)=delete
Broken copy constructor.
Gauss()
Default constructor (empty)
Definition: integral.h:537
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:546
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:552
void operator=(const Gauss &)=delete
Broken assignment operator.
Gauss(const Gauss &dummy)=delete
Broken copy constructor.
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:558
Class for multidimensional Gaussian integration rules, over intervals other than -1 to 1,...
Definition: integral.h:571
void operator=(const Gauss_Rescaled &)=delete
Broken assignment operator.
double Lower
Store for the lower and upper limits of integration and the range.
Definition: integral.h:574
Gauss_Rescaled()
Default constructor (empty)
Definition: integral.h:578
double knot(const unsigned &i, const unsigned &j) const
Return the rescaled knot values s[j] at integration point i.
Definition: integral.h:594
Gauss_Rescaled(double lower, double upper)
The constructor in this case takes the lower and upper arguments.
Definition: integral.h:587
double weight(const unsigned &i) const
Return the rescaled weight at integration point i.
Definition: integral.h:600
Gauss_Rescaled(const Gauss_Rescaled &dummy)=delete
Broken copy constructor.
Class for multidimensional Gaussian integration rules.
Definition: integral.h:145
Generic class for numerical integration schemes:
Definition: integral.h:49
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
Integral(const Integral &dummy)=delete
Broken copy constructor.
virtual ~Integral()
Virtual destructor (empty)
Definition: integral.h:61
Integral()
Default constructor (empty)
Definition: integral.h:52
virtual Vector< double > knot(const unsigned &i) const
Return local coordinates of i-th intergration point. Broken virtual.
Definition: integral.h:70
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
void operator=(const Integral &)=delete
Broken assignment operator.
An OomphLibError object which should be thrown when an run-time error is encountered....
Broken pseudo-integration scheme for points elements: Iit's not clear in general what this integratio...
Definition: integral.h:89
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] (j=0) of integration point i – deliberately broken!
Definition: integral.h:108
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:101
PointIntegral(const PointIntegral &dummy)=delete
Broken copy constructor.
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:120
PointIntegral()
Default constructor (empty)
Definition: integral.h:92
void operator=(const PointIntegral &)=delete
Broken assignment operator.
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:671
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:659
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:665
TGauss()
Default constructor (empty)
Definition: integral.h:650
void operator=(const TGauss &)=delete
Broken assignment operator.
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:711
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:717
void operator=(const TGauss &)=delete
Broken assignment operator.
TGauss()
Default constructor (empty)
Definition: integral.h:696
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:705
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:758
TGauss()
Default constructor (empty)
Definition: integral.h:743
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
void operator=(const TGauss &)=delete
Broken assignment operator.
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:752
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:764
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:803
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:797
void operator=(const TGauss &)=delete
Broken assignment operator.
TGauss()
Default constructor (empty)
Definition: integral.h:782
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:791
TGauss()
Default constructor (empty)
Definition: integral.h:967
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:988
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:976
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
void operator=(const TGauss &)=delete
Broken assignment operator.
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:982
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:1075
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1069
void operator=(const TGauss &)=delete
Broken assignment operator.
TGauss()
Default constructor (empty)
Definition: integral.h:1060
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1081
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:845
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:851
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:839
void operator=(const TGauss &)=delete
Broken assignment operator.
TGauss()
Default constructor (empty)
Definition: integral.h:830
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:885
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:891
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:897
void operator=(const TGauss &)=delete
Broken assignment operator.
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
TGauss()
Default constructor (empty)
Definition: integral.h:876
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:944
TGauss()
Default constructor (empty)
Definition: integral.h:923
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:932
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:938
void operator=(const TGauss &)=delete
Broken assignment operator.
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:1119
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1113
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1125
void operator=(const TGauss &)=delete
Broken assignment operator.
TGauss()
Default constructor (empty)
Definition: integral.h:1104
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:1027
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1021
TGauss()
Default constructor (empty)
Definition: integral.h:1012
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1033
void operator=(const TGauss &)=delete
Broken assignment operator.
void operator=(const TGauss &)=delete
Broken assignment operator.
TGauss()
Default constructor (empty)
Definition: integral.h:1150
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1171
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:1165
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1159
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1219
void operator=(const TGauss &)=delete
Broken assignment operator.
TGauss()
Default constructor (empty)
Definition: integral.h:1198
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1207
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:1213
void operator=(const TGauss &)=delete
Broken assignment operator.
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1255
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:1261
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1267
TGauss()
Default constructor (empty)
Definition: integral.h:1246
Class for Gaussian integration rules for triangles/tets.
Definition: integral.h:626
void gl_nodes(const unsigned &Nnode, Vector< double > &x)
Definition: orthpoly.cc:105
void gll_nodes(const unsigned &Nnode, Vector< double > &x)
Calculates the Gauss Lobatto Legendre abscissas for degree p = NNode-1.
Definition: orthpoly.cc:33
//////////////////////////////////////////////////////////////////// ////////////////////////////////...