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