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-2024 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 
42 namespace 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  //=========================================================
48  class Integral
49  {
50  public:
51  /// Default constructor (empty)
52  Integral(){};
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)
650  TGauss(){};
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)
696  TGauss(){};
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)
743  TGauss(){};
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)
782  TGauss(){};
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)
830  TGauss(){};
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)
876  TGauss(){};
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)
923  TGauss(){};
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)
967  TGauss(){};
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)
1012  TGauss(){};
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 #ifdef PARANOID
1059  /// A flag to track whether we have warned the user about the exterior knots
1060  /// and negative weights in this scheme.
1062 #endif
1063 
1064  public:
1065  /// Default constructor (empty)
1067  {
1068 #ifdef PARANOID
1069  // If the user has not been warned about the exterior knots and negative
1070  // weights that this integration scheme uses then do so and set the static
1071  // flag so that it does not happen again.
1072  if (!User_has_been_warned)
1073  {
1074  std::string warning_string =
1075  "The TGauss<2,16> integration scheme uses a high order Dunavant\n"
1076  "scheme which results in a couple of knots (slightly) outside of\n"
1077  "the element as well as some negative weights. These may be\n"
1078  "undesirable features depending on the use case and may also break\n"
1079  "some oomph-lib routines which expect points to be within the\n"
1080  "bounds of an element (locate_zeta). Please ensure that the\n"
1081  "integrand can be evaluated just outside the boundary of this\n"
1082  "element.";
1083  User_has_been_warned = true;
1085  warning_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1086  }
1087 #endif
1088  };
1089 
1090  /// Broken copy constructor
1091  TGauss(const TGauss& dummy) = delete;
1092 
1093  /// Broken assignment operator
1094  void operator=(const TGauss&) = delete;
1095 
1096  /// Number of integration points of the scheme
1097  unsigned nweight() const
1098  {
1099  return Npts;
1100  }
1101 
1102  /// Return coordinate x[j] of integration point i
1103  double knot(const unsigned& i, const unsigned& j) const
1104  {
1105  return Knot[i][j];
1106  }
1107 
1108  /// Return weight of integration point i
1109  double weight(const unsigned& i) const
1110  {
1111  return Weight[i];
1112  }
1113  };
1114 
1115  //------------------------------------------------------------
1116  //"Full integration" weights for 2D triangles
1117  // accurate up to order 15
1118  // http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_tri/quadrature_rules_tri.html
1119  //------------------------------------------------------------
1120 
1121  template<>
1122  class TGauss<2, 5> : public Integral
1123  {
1124  private:
1125  /// Number of integration points in the scheme
1126  static const unsigned Npts = 64;
1127  /// Array to hold the weights and knots (defined in cc file)
1128  static const double Knot[64][2], Weight[64];
1129 
1130  public:
1131  /// Default constructor (empty)
1132  TGauss(){};
1133 
1134  /// Broken copy constructor
1135  TGauss(const TGauss& dummy) = delete;
1136 
1137  /// Broken assignment operator
1138  void operator=(const TGauss&) = delete;
1139 
1140  /// Number of integration points of the scheme
1141  unsigned nweight() const
1142  {
1143  return Npts;
1144  }
1145 
1146  /// Return coordinate x[j] of integration point i
1147  double knot(const unsigned& i, const unsigned& j) const
1148  {
1149  return Knot[i][j];
1150  }
1151 
1152  /// Return weight of integration point i
1153  double weight(const unsigned& i) const
1154  {
1155  return Weight[i];
1156  }
1157  };
1158 
1159  //=========================================================
1160  /// 3D Gaussian integration class for tets.
1161  /// Four integration points. This integration scheme can
1162  /// integrate up to second-order polynomials exactly and
1163  /// is therefore a suitable "full" integration scheme
1164  /// for linear (four-node) elements in which the
1165  /// highest-order polynomial is quadratic.
1166  //=========================================================
1167  template<>
1168  class TGauss<3, 2> : public Integral
1169  {
1170  private:
1171  /// Number of integration points in the scheme
1172  static const unsigned Npts = 4;
1173  /// Array to hold the weights and knots (defined in cc file)
1174  static const double Knot[4][3], Weight[4];
1175 
1176  public:
1177  /// Default constructor (empty)
1178  TGauss(){};
1179 
1180  /// Broken copy constructor
1181  TGauss(const TGauss& dummy) = delete;
1182 
1183  /// Broken assignment operator
1184  void operator=(const TGauss&) = delete;
1185 
1186  /// Number of integration points of the scheme
1187  unsigned nweight() const
1188  {
1189  return Npts;
1190  }
1191 
1192  /// Return coordinate x[j] of integration point i
1193  double knot(const unsigned& i, const unsigned& j) const
1194  {
1195  return Knot[i][j];
1196  }
1197 
1198  /// Return weight of integration point i
1199  double weight(const unsigned& i) const
1200  {
1201  return Weight[i];
1202  }
1203  };
1204 
1205 
1206  //=========================================================
1207  /// 3D Gaussian integration class for tets.
1208  /// Eleven integration points. This integration scheme can
1209  /// integrate up to fourth-order polynomials exactly and
1210  /// is therefore a suitable "full" integration scheme
1211  /// for quadratic (ten-node) elements in which the
1212  /// highest-order polynomial is fourth order.
1213  /// The numbers are from Keast CMAME 55 pp339-348 (1986)
1214  //=========================================================
1215  template<>
1216  class TGauss<3, 3> : public Integral
1217  {
1218  private:
1219  /// Number of integration points in the scheme
1220  static const unsigned Npts = 11;
1221  /// Array to hold the weights and knots (defined in cc file)
1222  static const double Knot[11][3], Weight[11];
1223 
1224  public:
1225  /// Default constructor (empty)
1226  TGauss(){};
1227 
1228  /// Broken copy constructor
1229  TGauss(const TGauss& dummy) = delete;
1230 
1231  /// Broken assignment operator
1232  void operator=(const TGauss&) = delete;
1233 
1234  /// Number of integration points of the scheme
1235  unsigned nweight() const
1236  {
1237  return Npts;
1238  }
1239 
1240  /// Return coordinate x[j] of integration point i
1241  double knot(const unsigned& i, const unsigned& j) const
1242  {
1243  return Knot[i][j];
1244  }
1245 
1246  /// Return weight of integration point i
1247  double weight(const unsigned& i) const
1248  {
1249  return Weight[i];
1250  }
1251  };
1252 
1253 
1254  //=========================================================
1255  /// 3D Gaussian integration class for tets.
1256  /// 45 integration points. This integration scheme can
1257  /// integrate up to eighth-order polynomials exactly and
1258  /// is therefore a suitable "full" integration scheme
1259  /// for quartic elements in which the
1260  /// highest-order polynomial is fourth order.
1261  /// The numbers are from Keast CMAME 55 pp339-348 (1986)
1262  //=========================================================
1263  template<>
1264  class TGauss<3, 5> : public Integral
1265  {
1266  private:
1267  /// Number of integration points in the scheme
1268  static const unsigned Npts = 45;
1269  /// Array to hold the weights and knots (defined in cc file)
1270  static const double Knot[45][3], Weight[45];
1271 
1272  public:
1273  /// Default constructor (empty)
1274  TGauss(){};
1275 
1276  /// Broken copy constructor
1277  TGauss(const TGauss& dummy) = delete;
1278 
1279  /// Broken assignment operator
1280  void operator=(const TGauss&) = delete;
1281 
1282  /// Number of integration points of the scheme
1283  unsigned nweight() const
1284  {
1285  return Npts;
1286  }
1287 
1288  /// Return coordinate x[j] of integration point i
1289  double knot(const unsigned& i, const unsigned& j) const
1290  {
1291  return Knot[i][j];
1292  }
1293 
1294  /// Return weight of integration point i
1295  double weight(const unsigned& i) const
1296  {
1297  return Weight[i];
1298  }
1299  };
1300 
1301 
1302  //===================================================================
1303  /// Class for multidimensional Gauss Lobatto Legendre integration
1304  /// rules
1305  /// empty - just establishes template parameters
1306  //===================================================================
1307  template<unsigned DIM, unsigned NPTS_1D>
1309  {
1310  };
1311 
1312  //===================================================================
1313  /// 1D Gauss Lobatto Legendre integration class
1314  //===================================================================
1315  template<unsigned NPTS_1D>
1316  class GaussLobattoLegendre<1, NPTS_1D> : public Integral
1317  {
1318  private:
1319  /// Number of integration points in scheme
1320  static const unsigned Npts = NPTS_1D;
1321  /// Array to hold weight and knot points
1322  // These are not constant, because they are calculated on the fly
1323  double Knot[NPTS_1D][1], Weight[NPTS_1D];
1324 
1325  public:
1326  /// Deafault constructor. Calculates and stores GLL nodes
1328 
1329  /// Number of integration points of the scheme
1330  unsigned nweight() const
1331  {
1332  return Npts;
1333  }
1334 
1335  /// Return coordinate s[j] (j=0) of integration point i
1336  double knot(const unsigned& i, const unsigned& j) const
1337  {
1338  return Knot[i][j];
1339  }
1340 
1341  /// Return weight of integration point i
1342  double weight(const unsigned& i) const
1343  {
1344  return Weight[i];
1345  }
1346  };
1347 
1348 
1349  //=============================================================
1350  /// Calculate positions and weights for the 1D Gauss Lobatto
1351  /// Legendre integration class
1352  //=============================================================
1353  template<unsigned NPTS_1D>
1355  {
1356  // Temporary storage for the integration points
1357  Vector<double> s(NPTS_1D), w(NPTS_1D);
1358  // call the function that calculates the points
1359  Orthpoly::gll_nodes(NPTS_1D, s, w);
1360  // Populate the arrays
1361  for (unsigned i = 0; i < NPTS_1D; i++)
1362  {
1363  Knot[i][0] = s[i];
1364  Weight[i] = w[i];
1365  }
1366  }
1367 
1368 
1369  //===================================================================
1370  /// 2D Gauss Lobatto Legendre integration class
1371  //===================================================================
1372  template<unsigned NPTS_1D>
1373  class GaussLobattoLegendre<2, NPTS_1D> : public Integral
1374  {
1375  private:
1376  /// Number of integration points in scheme
1377  static const unsigned long int Npts = NPTS_1D * NPTS_1D;
1378 
1379  /// Array to hold weight and knot points
1380  double Knot[NPTS_1D * NPTS_1D][2],
1381  Weight[NPTS_1D * NPTS_1D]; // COULDN'T MAKE THESE
1382  // const BECAUSE THEY ARE CALCULATED (at least currently)
1383 
1384  public:
1385  /// Deafault constructor. Calculates and stores GLL nodes
1387 
1388  /// Number of integration points of the scheme
1389  unsigned nweight() const
1390  {
1391  return Npts;
1392  }
1393 
1394  /// Return coordinate s[j] (j=0) of integration point i
1395  double knot(const unsigned& i, const unsigned& j) const
1396  {
1397  return Knot[i][j];
1398  }
1399 
1400  /// Return weight of integration point i
1401  double weight(const unsigned& i) const
1402  {
1403  return Weight[i];
1404  }
1405  };
1406 
1407  //=============================================================
1408  /// Calculate positions and weights for the 2D Gauss Lobatto
1409  /// Legendre integration class
1410  //=============================================================
1411 
1412  template<unsigned NPTS_1D>
1414  {
1415  // Tempoarary storage for the 1D knots and weights
1416  Vector<double> s(NPTS_1D), w(NPTS_1D);
1417  // Call the function to populate the array
1418  Orthpoly::gll_nodes(NPTS_1D, s, w);
1419  int i_fast = 0, i_slow = 0;
1420  for (unsigned i = 0; i < NPTS_1D * NPTS_1D; i++)
1421  {
1422  if (i_fast == NPTS_1D)
1423  {
1424  i_fast = 0;
1425  i_slow++;
1426  }
1427  Knot[i][0] = s[i_fast];
1428  Knot[i][1] = s[i_slow];
1429  Weight[i] = w[i_fast] * w[i_slow];
1430  i_fast++;
1431  }
1432  }
1433 
1434 
1435  //===================================================================
1436  /// 3D Gauss Lobatto Legendre integration class
1437  //===================================================================
1438  template<unsigned NPTS_1D>
1439  class GaussLobattoLegendre<3, NPTS_1D> : public Integral
1440  {
1441  private:
1442  /// Number of integration points in scheme
1443  static const unsigned long int Npts = NPTS_1D * NPTS_1D * NPTS_1D;
1444 
1445  /// Array to hold weight and knot points
1446  double Knot[NPTS_1D * NPTS_1D * NPTS_1D][3],
1447  Weight[NPTS_1D * NPTS_1D * NPTS_1D]; // COULDN'T MAKE THESE
1448  // const BECAUSE THEY ARE CALCULATED (at least currently)
1449 
1450  public:
1451  /// Deafault constructor. Calculates and stores GLL nodes
1453 
1454  /// Number of integration points of the scheme
1455  unsigned nweight() const
1456  {
1457  return Npts;
1458  }
1459 
1460  /// Return coordinate s[j] (j=0) of integration point i
1461  double knot(const unsigned& i, const unsigned& j) const
1462  {
1463  return Knot[i][j];
1464  }
1465 
1466  /// Return weight of integration point i
1467  double weight(const unsigned& i) const
1468  {
1469  return Weight[i];
1470  }
1471  };
1472 
1473  //=============================================================
1474  /// Calculate positions and weights for the 3D Gauss Lobatto
1475  /// Legendre integration class
1476  //=============================================================
1477 
1478  template<unsigned NPTS_1D>
1480  {
1481  // Tempoarary storage for the 1D knots and weights
1482  Vector<double> s(NPTS_1D), w(NPTS_1D);
1483  // Call the function to populate the array
1484  Orthpoly::gll_nodes(NPTS_1D, s, w);
1485  for (unsigned k = 0; k < NPTS_1D; k++)
1486  {
1487  for (unsigned j = 0; j < NPTS_1D; j++)
1488  {
1489  for (unsigned i = 0; i < NPTS_1D; i++)
1490  {
1491  unsigned index = NPTS_1D * NPTS_1D * k + NPTS_1D * j + i;
1492  Knot[index][0] = s[i];
1493  Knot[index][1] = s[j];
1494  Knot[index][2] = s[k];
1495  Weight[index] = w[i] * w[j] * w[k];
1496  }
1497  }
1498  }
1499  }
1500 
1501  /// /////////////////////////////////////////////////////////////////
1502  /// /////////////////////////////////////////////////////////////////
1503  /// /////////////////////////////////////////////////////////////////
1504 
1505 
1506  //===================================================================
1507  /// Class for multidimensional Gauss Legendre integration
1508  /// rules
1509  /// empty - just establishes template parameters
1510  //===================================================================
1511  template<unsigned DIM, unsigned NPTS_1D>
1513  {
1514  };
1515 
1516  //===================================================================
1517  /// 1D Gauss Legendre integration class
1518  //===================================================================
1519  template<unsigned NPTS_1D>
1520  class GaussLegendre<1, NPTS_1D> : public Integral
1521  {
1522  private:
1523  /// Number of integration points in scheme
1524  static const unsigned Npts = NPTS_1D;
1525  /// Array to hold weight and knot points
1526  // These are not constant, because they are calculated on the fly
1527  double Knot[NPTS_1D][1], Weight[NPTS_1D];
1528 
1529  public:
1530  /// Deafault constructor. Calculates and stores GLL nodes
1531  GaussLegendre();
1532 
1533  /// Number of integration points of the scheme
1534  unsigned nweight() const
1535  {
1536  return Npts;
1537  }
1538 
1539  /// Return coordinate s[j] (j=0) of integration point i
1540  double knot(const unsigned& i, const unsigned& j) const
1541  {
1542  return Knot[i][j];
1543  }
1544 
1545  /// Return weight of integration point i
1546  double weight(const unsigned& i) const
1547  {
1548  return Weight[i];
1549  }
1550  };
1551 
1552 
1553  //=============================================================
1554  /// Calculate positions and weights for the 1D Gauss
1555  /// Legendre integration class
1556  //=============================================================
1557  template<unsigned NPTS_1D>
1559  {
1560  // Temporary storage for the integration points
1561  Vector<double> s(NPTS_1D), w(NPTS_1D);
1562  // call the function that calculates the points
1563  Orthpoly::gl_nodes(NPTS_1D, s, w);
1564  // Populate the arrays
1565  for (unsigned i = 0; i < NPTS_1D; i++)
1566  {
1567  Knot[i][0] = s[i];
1568  Weight[i] = w[i];
1569  }
1570  }
1571 
1572 
1573  //===================================================================
1574  /// 2D Gauss Legendre integration class
1575  //===================================================================
1576  template<unsigned NPTS_1D>
1577  class GaussLegendre<2, NPTS_1D> : public Integral
1578  {
1579  private:
1580  /// Number of integration points in scheme
1581  static const unsigned long int Npts = NPTS_1D * NPTS_1D;
1582 
1583  /// Array to hold weight and knot points
1584  double Knot[NPTS_1D * NPTS_1D][2],
1585  Weight[NPTS_1D * NPTS_1D]; // COULDN'T MAKE THESE
1586  // const BECAUSE THEY ARE CALCULATED (at least currently)
1587 
1588  public:
1589  /// Deafault constructor. Calculates and stores GLL nodes
1590  GaussLegendre();
1591 
1592  /// Number of integration points of the scheme
1593  unsigned nweight() const
1594  {
1595  return Npts;
1596  }
1597 
1598  /// Return coordinate s[j] (j=0) of integration point i
1599  double knot(const unsigned& i, const unsigned& j) const
1600  {
1601  return Knot[i][j];
1602  }
1603 
1604  /// Return weight of integration point i
1605  double weight(const unsigned& i) const
1606  {
1607  return Weight[i];
1608  }
1609  };
1610 
1611  //=============================================================
1612  /// Calculate positions and weights for the 2D Gauss
1613  /// Legendre integration class
1614  //=============================================================
1615 
1616  template<unsigned NPTS_1D>
1618  {
1619  // Tempoarary storage for the 1D knots and weights
1620  Vector<double> s(NPTS_1D), w(NPTS_1D);
1621  // Call the function to populate the array
1622  Orthpoly::gl_nodes(NPTS_1D, s, w);
1623  int i_fast = 0, i_slow = 0;
1624  for (unsigned i = 0; i < NPTS_1D * NPTS_1D; i++)
1625  {
1626  if (i_fast == NPTS_1D)
1627  {
1628  i_fast = 0;
1629  i_slow++;
1630  }
1631  Knot[i][0] = s[i_fast];
1632  Knot[i][1] = s[i_slow];
1633  Weight[i] = w[i_fast] * w[i_slow];
1634  i_fast++;
1635  }
1636  }
1637 
1638 
1639  //===================================================================
1640  /// 3D Gauss Legendre integration class
1641  //===================================================================
1642  template<unsigned NPTS_1D>
1643  class GaussLegendre<3, NPTS_1D> : public Integral
1644  {
1645  private:
1646  /// Number of integration points in scheme
1647  static const unsigned long int Npts = NPTS_1D * NPTS_1D * NPTS_1D;
1648 
1649  /// Array to hold weight and knot points
1650  double Knot[NPTS_1D * NPTS_1D * NPTS_1D][3],
1651  Weight[NPTS_1D * NPTS_1D * NPTS_1D]; // COULDN'T MAKE THESE
1652  // const BECAUSE THEY ARE CALCULATED (at least currently)
1653 
1654  public:
1655  /// Deafault constructor. Calculates and stores GLL nodes
1656  GaussLegendre();
1657 
1658  /// Number of integration points of the scheme
1659  unsigned nweight() const
1660  {
1661  return Npts;
1662  }
1663 
1664  /// Return coordinate s[j] (j=0) of integration point i
1665  double knot(const unsigned& i, const unsigned& j) const
1666  {
1667  return Knot[i][j];
1668  }
1669 
1670  /// Return weight of integration point i
1671  double weight(const unsigned& i) const
1672  {
1673  return Weight[i];
1674  }
1675  };
1676 
1677  //=============================================================
1678  /// Calculate positions and weights for the 3D Gauss
1679  /// Legendre integration class
1680  //=============================================================
1681 
1682  template<unsigned NPTS_1D>
1684  {
1685  // Tempoarary storage for the 1D knots and weights
1686  Vector<double> s(NPTS_1D), w(NPTS_1D);
1687  // Call the function to populate the array
1688  Orthpoly::gl_nodes(NPTS_1D, s, w);
1689  for (unsigned k = 0; k < NPTS_1D; k++)
1690  {
1691  for (unsigned j = 0; j < NPTS_1D; j++)
1692  {
1693  for (unsigned i = 0; i < NPTS_1D; i++)
1694  {
1695  unsigned index = NPTS_1D * NPTS_1D * k + NPTS_1D * j + i;
1696  Knot[index][0] = s[i];
1697  Knot[index][1] = s[j];
1698  Knot[index][2] = s[k];
1699  Weight[index] = w[i] * w[j] * w[k];
1700  }
1701  }
1702  }
1703  }
1704 
1705 
1706 } // namespace oomph
1707 
1708 #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:1546
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] (j=0) of integration point i.
Definition: integral.h:1540
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1534
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1593
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1605
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] (j=0) of integration point i.
Definition: integral.h:1599
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1659
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1671
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] (j=0) of integration point i.
Definition: integral.h:1665
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
Definition: integral.h:1513
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1342
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1330
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] (j=0) of integration point i.
Definition: integral.h:1336
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1389
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1401
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] (j=0) of integration point i.
Definition: integral.h:1395
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1467
double knot(const unsigned &i, const unsigned &j) const
Return coordinate s[j] (j=0) of integration point i.
Definition: integral.h:1461
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1455
Class for multidimensional Gauss Lobatto Legendre integration rules empty - just establishes template...
Definition: integral.h:1309
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 Vector< double > knot(const unsigned &i) const
Return local coordinates of i-th intergration point. Broken virtual.
Definition: integral.h:70
virtual ~Integral()
Virtual destructor (empty)
Definition: integral.h:61
Integral()
Default constructor (empty)
Definition: integral.h:52
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....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
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
static bool User_has_been_warned
A flag to track whether we have warned the user about the exterior knots and negative weights in this...
Definition: integral.h:1061
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:1103
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1097
void operator=(const TGauss &)=delete
Broken assignment operator.
TGauss()
Default constructor (empty)
Definition: integral.h:1066
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1109
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:1147
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1141
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1153
void operator=(const TGauss &)=delete
Broken assignment operator.
TGauss()
Default constructor (empty)
Definition: integral.h:1132
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:1178
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1199
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:1193
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1187
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1247
void operator=(const TGauss &)=delete
Broken assignment operator.
TGauss()
Default constructor (empty)
Definition: integral.h:1226
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1235
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:1241
void operator=(const TGauss &)=delete
Broken assignment operator.
unsigned nweight() const
Number of integration points of the scheme.
Definition: integral.h:1283
double knot(const unsigned &i, const unsigned &j) const
Return coordinate x[j] of integration point i.
Definition: integral.h:1289
TGauss(const TGauss &dummy)=delete
Broken copy constructor.
double weight(const unsigned &i) const
Return weight of integration point i.
Definition: integral.h:1295
TGauss()
Default constructor (empty)
Definition: integral.h:1274
Class for Gaussian integration rules for triangles/tets.
Definition: integral.h:626
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...