square [0,1]x[0,1], etc. This information is used together with
objects of the \Ref{FiniteElement} class to compute the values stored
in the \Ref{FEValues} objects.
- */
+
+ There are a number of derived classes, denoting concrete integration
+ formulae. These are named by a prefixed #Q#, the name of the formula
+ (e.g. #Gauss#) and finally the order of integration. For example,
+ #QGauss2<dim># denotes a second order Gauss integration formula in
+ any dimension. Second order means that it integrates polynomials of
+ third order exact. In general, a formula of order #n# exactly
+ integrates polynomials of order #2n-1#.
+*/
template <int dim>
class Quadrature {
public:
#include <fe/quadrature.h>
+/**
+ Second order Gauss quadrature formula.
+
+ Reference: Ward Cheney, David Kincaid: Numerical Mathematics and Computing.
+*/
template <int dim>
class QGauss2 : public Quadrature<dim> {
public:
+/**
+ Third order Gauss quadrature formula.
+
+ Reference: Ward Cheney, David Kincaid: Numerical Mathematics and Computing.
+*/
template <int dim>
-class QGauss2x4 : public Quadrature<dim> {
+class QGauss3 : public Quadrature<dim> {
public:
- QGauss2x4 ();
+ QGauss3 ();
};
+/**
+ Fourth order Gauss quadrature formula.
+
+ Reference: Ward Cheney, David Kincaid: Numerical Mathematics and Computing.
+*/
template <int dim>
class QGauss4 : public Quadrature<dim> {
public:
+
+/**
+ Fifth order Gauss quadrature formula.
+
+ Reference: Ward Cheney, David Kincaid: Numerical Mathematics and Computing.
+*/
+template <int dim>
+class QGauss5 : public Quadrature<dim> {
+ public:
+ QGauss5 ();
+};
+
+
+
+/**
+ Sixth order Gauss quadrature formula. I have not found explicite
+ representations of the zeros of the Legendre functions of sixth
+ and higher degree. If anyone finds them, please replace the existing
+ numbers by these expressions.
+
+ Reference: J. E. Akin: Application and Implementation of Finite
+ Element Methods
+*/
+template <int dim>
+class QGauss6 : public Quadrature<dim> {
+ public:
+ QGauss6 ();
+};
+
+
+
+/**
+ Seventh order Gauss quadrature formula. I have not found explicite
+ representations of the zeros of the Legendre functions of sixth
+ and higher degree. If anyone finds them, please replace the existing
+ numbers by these expressions.
+
+ Reference: J. E. Akin: Application and Implementation of Finite
+ Element Methods
+*/
+template <int dim>
+class QGauss7 : public Quadrature<dim> {
+ public:
+ QGauss7 ();
+};
+
+
+
+/**
+ Eighth order Gauss quadrature formula. I have not found explicite
+ representations of the zeros of the Legendre functions of sixth
+ and higher degree. If anyone finds them, please replace the existing
+ numbers by these expressions.
+
+ Reference: J. E. Akin: Application and Implementation of Finite
+ Element Methods
+*/
template <int dim>
class QGauss8 : public Quadrature<dim> {
public:
+
+
+/**
+ First order midpoint quadrature rule.
+*/
template <int dim>
class QMidpoint : public Quadrature<dim> {
public:
+/**
+ Simpson quadrature rule.
+*/
template <int dim>
class QSimpson : public Quadrature<dim> {
public:
+/**
+ Trapezoidal quadrature rule.
+*/
template <int dim>
class QTrapez : public Quadrature<dim> {
public:
/* $Id$ */
#include <fe/quadrature_lib.h>
+#include <cmath>
QGauss2<1>::QGauss2 () :
Quadrature<1> (2)
{
- static const double xpts[] = { 0.288675135, 0.71132486 };
- static const double wts[] = { 0.5, 0.5 };
+ // points on [-1,1]
+ static const double xpts_normal[] = { -sqrt(1./3.), sqrt(1./3.) };
+ // weights on [-1,1]
+ static const double wts_normal[] = { 1., 1. };
+
+ // points and weights on [0,1]
+ static const double xpts[] = { (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2. };
+ static const double wts[] = { wts_normal[0]/2.,
+ wts_normal[1]/2. };
for (unsigned int i=0; i<n_quadrature_points; ++i)
{
-QGauss2x4<1>::QGauss2x4 () :
- Quadrature<1> (8)
+QGauss2<2>::QGauss2 () :
+ Quadrature<2> (4)
{
- static const double G0=0.930568156,
- G1=0.669990522,
- G2=0.330009478,
- G3=0.069431844;
- static const double W0=0.173927423,
- W1=0.326072577;
-
- static const double xpts[] = { 0.5*G0, 0.5*G1, 0.5*G2, 0.5*G3,
- 0.5*G0+0.5, 0.5*G1+0.5, 0.5*G2+0.5, 0.5*G3+0.5 };
- static const double wts[] = { 0.5*W0, 0.5*W1, 0.5*W1, 0.5*W0,
- 0.5*W0, 0.5*W1, 0.5*W1, 0.5*W0 };
+ // points on [-1,1]
+ static const double xpts_normal[] = { -sqrt(1./3.), sqrt(1./3.) };
+ // weights on [-1,1]
+ static const double wts_normal[] = { 1., 1. };
+
+ // points and weights on [0,1]^2
+ static const double xpts[] = { (xpts_normal[0]+1)/2.,
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[1]+1)/2. };
+ static const double ypts[] = { (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2. };
+ static const double wts[] = { wts_normal[0]/4.,
+ wts_normal[1]/4.,
+ wts_normal[0]/4.,
+ wts_normal[1]/4. };
+
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
+ {
+ quadrature_points[i] = Point<2>(xpts[i], ypts[i]);
+ weights[i] = wts[i];
+ };
+};
+
+
+
+
+
+QGauss3<1>::QGauss3 () :
+ Quadrature<1> (3)
+{
+ // points on [-1,1]
+ static const double xpts_normal[] = { -sqrt(3./5.),
+ 0.,
+ sqrt(3./5.) };
+ // weights on [-1,1]
+ static const double wts_normal[] = { 5./9.,
+ 8./9.,
+ 5./9. };
+
+ // points and weights on [0,1]
+ static const double xpts[] = { (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2. };
+ static const double wts[] = { wts_normal[0]/2.,
+ wts_normal[1]/2.,
+ wts_normal[2]/2. };
for (unsigned int i=0; i<n_quadrature_points; ++i)
{
+QGauss3<2>::QGauss3 () :
+ Quadrature<2> (9)
+{
+ // points on [-1,1]
+ static const double xpts_normal[] = { -sqrt(3./5.),
+ 0.,
+ sqrt(3./5.) };
+ // weights on [-1,1]
+ static const double wts_normal[] = { 5./9.,
+ 8./9.,
+ 5./9. };
+
+ // points and weights on [0,1]^2
+ static const double xpts[] = { (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2. };
+ static const double ypts[] = { (xpts_normal[0]+1)/2.,
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[2]+1)/2. };
+ static const double wts[] = { wts_normal[0]/2.*wts_normal[0]/2.,
+ wts_normal[1]/2.*wts_normal[0]/2.,
+ wts_normal[2]/2.*wts_normal[0]/2.,
+ wts_normal[0]/2.*wts_normal[1]/2.,
+ wts_normal[1]/2.*wts_normal[1]/2.,
+ wts_normal[2]/2.*wts_normal[1]/2.,
+ wts_normal[0]/2.*wts_normal[2]/2.,
+ wts_normal[1]/2.*wts_normal[2]/2.,
+ wts_normal[2]/2.*wts_normal[2]/2. };
+
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
+ {
+ quadrature_points[i] = Point<2>(xpts[i], ypts[i]);
+ weights[i] = wts[i];
+ };
+};
+
+
+
+
+
+
QGauss4<1>::QGauss4 () :
Quadrature<1> (4)
{
- static const double G0=0.930568156,
- G1=0.669990522,
- G2=0.330009478,
- G3=0.069431844;
- static const double W0=0.173927423,
- W1=0.326072577;
+ // points on [-1,1]
+ static const double xpts_normal[] = { -sqrt(1./7.*(3-4*sqrt(0.3))),
+ -sqrt(1./7.*(3+4*sqrt(0.3))),
+ +sqrt(1./7.*(3-4*sqrt(0.3))),
+ +sqrt(1./7.*(3+4*sqrt(0.3))) };
+ // weights on [-1,1]
+ static const double wts_normal[] = { 1./2. + 1./12.*sqrt(10./3.),
+ 1./2. - 1./12.*sqrt(10./3.),
+ 1./2. + 1./12.*sqrt(10./3.),
+ 1./2. - 1./12.*sqrt(10./3.) };
+
+ // points and weights on [0,1]
+ static const double xpts[] = { (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[3]+1)/2. };
+ static const double wts[] = { wts_normal[0]/2.,
+ wts_normal[1]/2.,
+ wts_normal[2]/2.,
+ wts_normal[3]/2. };
+
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
+ {
+ quadrature_points[i] = Point<1>(xpts[i]);
+ weights[i] = wts[i];
+ };
+};
+
+
+
+QGauss4<2>::QGauss4 () :
+ Quadrature<2> (16)
+{
+ // points on [-1,1]
+ static const double xpts_normal[] = { -sqrt(1./7.*(3-4*sqrt(0.3))),
+ -sqrt(1./7.*(3+4*sqrt(0.3))),
+ +sqrt(1./7.*(3-4*sqrt(0.3))),
+ +sqrt(1./7.*(3+4*sqrt(0.3))) };
+ // weights on [-1,1]
+ static const double wts_normal[] = { 1./2. + 1./12.*sqrt(10./3.),
+ 1./2. - 1./12.*sqrt(10./3.),
+ 1./2. + 1./12.*sqrt(10./3.),
+ 1./2. - 1./12.*sqrt(10./3.) };
+
+ // points and weights on [0,1]^2
+ static const double xpts[] = { (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[3]+1)/2. };
+ static const double ypts[] = { (xpts_normal[0]+1)/2.,
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[3]+1)/2.};
+ static const double wts[] = { wts_normal[0]/2.*wts_normal[0]/2.,
+ wts_normal[1]/2.*wts_normal[0]/2.,
+ wts_normal[2]/2.*wts_normal[0]/2.,
+ wts_normal[3]/2.*wts_normal[0]/2.,
+ wts_normal[0]/2.*wts_normal[1]/2.,
+ wts_normal[1]/2.*wts_normal[1]/2.,
+ wts_normal[2]/2.*wts_normal[1]/2.,
+ wts_normal[3]/2.*wts_normal[1]/2.,
+ wts_normal[0]/2.*wts_normal[2]/2.,
+ wts_normal[1]/2.*wts_normal[2]/2.,
+ wts_normal[2]/2.*wts_normal[2]/2.,
+ wts_normal[3]/2.*wts_normal[2]/2.,
+ wts_normal[0]/2.*wts_normal[3]/2.,
+ wts_normal[1]/2.*wts_normal[3]/2.,
+ wts_normal[2]/2.*wts_normal[3]/2.,
+ wts_normal[3]/2.*wts_normal[3]/2. };
+
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
+ {
+ quadrature_points[i] = Point<2>(xpts[i], ypts[i]);
+ weights[i] = wts[i];
+ };
+};
+
+
+
+
+QGauss5<1>::QGauss5 () :
+ Quadrature<1> (5)
+{
+ // points on [-1,1]
+ static const double xpts_normal[] = { -sqrt(1./9.*(5.-2*sqrt(10./7.))),
+ -sqrt(1./9.*(5.+2*sqrt(10./7.))),
+ 0,
+ +sqrt(1./9.*(5.-2*sqrt(10./7.))),
+ +sqrt(1./9.*(5.+2*sqrt(10./7.))) };
+ // weights on [-1,1]
+ static const double wts_normal[] = { 0.3*(-0.7+5.*sqrt(0.7))/(-2.+5.*sqrt(0.7)),
+ 0.3*(+0.7+5.*sqrt(0.7))/(+2.+5.*sqrt(0.7)),
+ 128./225.,
+ 0.3*(-0.7+5.*sqrt(0.7))/(-2.+5.*sqrt(0.7)),
+ 0.3*(+0.7+5.*sqrt(0.7))/(+2.+5.*sqrt(0.7)) };
+
+ // points and weights on [0,1]
+ static const double xpts[] = { (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[4]+1)/2. };
+ static const double wts[] = { wts_normal[0]/2.,
+ wts_normal[1]/2.,
+ wts_normal[2]/2.,
+ wts_normal[3]/2.,
+ wts_normal[4]/2. };
+
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
+ {
+ quadrature_points[i] = Point<1>(xpts[i]);
+ weights[i] = wts[i];
+ };
+};
+
+
+
+QGauss5<2>::QGauss5 () :
+ Quadrature<2> (25)
+{
+ // points on [-1,1]
+ static const double xpts_normal[] = { -sqrt(1./9.*(5.-2*sqrt(10./7.))),
+ -sqrt(1./9.*(5.+2*sqrt(10./7.))),
+ 0,
+ +sqrt(1./9.*(5.-2*sqrt(10./7.))),
+ +sqrt(1./9.*(5.+2*sqrt(10./7.))) };
+ // weights on [-1,1]
+ static const double wts_normal[] = { 0.3*(-0.7+5.*sqrt(0.7))/(-2.+5.*sqrt(0.7)),
+ 0.3*(+0.7+5.*sqrt(0.7))/(+2.+5.*sqrt(0.7)),
+ 128./225.,
+ 0.3*(-0.7+5.*sqrt(0.7))/(-2.+5.*sqrt(0.7)),
+ 0.3*(+0.7+5.*sqrt(0.7))/(+2.+5.*sqrt(0.7)) };
+
+ // points and weights on [0,1]^2
+ static const double xpts[] = { (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[4]+1)/2.,
+
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[4]+1)/2.,
+
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[4]+1)/2.,
+
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[4]+1)/2.,
+
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[4]+1)/2. };
+
+ static const double ypts[] = { (xpts_normal[0]+1)/2.,
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[0]+1)/2.,
+ (xpts_normal[0]+1)/2.,
+
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+
+ (xpts_normal[4]+1)/2.,
+ (xpts_normal[4]+1)/2.,
+ (xpts_normal[4]+1)/2.,
+ (xpts_normal[4]+1)/2.,
+ (xpts_normal[4]+1)/2. };
+
+ static const double wts[] = { wts_normal[0]/2.*wts_normal[0]/2.,
+ wts_normal[1]/2.*wts_normal[0]/2.,
+ wts_normal[2]/2.*wts_normal[0]/2.,
+ wts_normal[3]/2.*wts_normal[0]/2.,
+ wts_normal[4]/2.*wts_normal[0]/2.,
+
+ wts_normal[0]/2.*wts_normal[1]/2.,
+ wts_normal[1]/2.*wts_normal[1]/2.,
+ wts_normal[2]/2.*wts_normal[1]/2.,
+ wts_normal[3]/2.*wts_normal[1]/2.,
+ wts_normal[4]/2.*wts_normal[1]/2.,
+
+ wts_normal[0]/2.*wts_normal[2]/2.,
+ wts_normal[1]/2.*wts_normal[2]/2.,
+ wts_normal[2]/2.*wts_normal[2]/2.,
+ wts_normal[3]/2.*wts_normal[2]/2.,
+ wts_normal[4]/2.*wts_normal[2]/2.,
+
+ wts_normal[0]/2.*wts_normal[3]/2.,
+ wts_normal[1]/2.*wts_normal[3]/2.,
+ wts_normal[2]/2.*wts_normal[3]/2.,
+ wts_normal[3]/2.*wts_normal[3]/2.,
+ wts_normal[4]/2.*wts_normal[3]/2.,
+
+ wts_normal[0]/2.*wts_normal[4]/2.,
+ wts_normal[1]/2.*wts_normal[4]/2.,
+ wts_normal[2]/2.*wts_normal[4]/2.,
+ wts_normal[3]/2.*wts_normal[4]/2.,
+ wts_normal[4]/2.*wts_normal[4]/2. };
- static const double xpts[] = { G0, G1, G2, G3 };
- static const double wts[] = { W0, W1, W1, W0 };
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
+ {
+ quadrature_points[i] = Point<2>(xpts[i], ypts[i]);
+ weights[i] = wts[i];
+ };
+};
+
+
+
+QGauss6<1>::QGauss6 () :
+ Quadrature<1> (6)
+{
+ // points on [-1,1]
+ static const double xpts_normal[] = { -0.932469514203152,
+ -0.661209386466265,
+ -0.238619186083197,
+ +0.238619186083197,
+ +0.661209386466265,
+ +0.932469514203152 };
+ // weights on [-1,1]
+ static const double wts_normal[] = { 0.171324492379170,
+ 0.360761573048139,
+ 0.467913934572691,
+ 0.467913934572691,
+ 0.360761573048139,
+ 0.171324492379170 };
+
+ // points and weights on [0,1]
+ static const double xpts[] = { (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[4]+1)/2.,
+ (xpts_normal[5]+1)/2. };
+ static const double wts[] = { wts_normal[0]/2.,
+ wts_normal[1]/2.,
+ wts_normal[2]/2.,
+ wts_normal[3]/2.,
+ wts_normal[4]/2.,
+ wts_normal[5]/2. };
+
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
+ {
+ quadrature_points[i] = Point<1>(xpts[i]);
+ weights[i] = wts[i];
+ };
+};
+
+
+
+QGauss7<1>::QGauss7 () :
+ Quadrature<1> (7)
+{
+ // points on [-1,1]
+ static const double xpts_normal[] = { -0.949107912342759,
+ -0.741531185599394,
+ -0.405845151377397,
+ 0,
+ +0.405845151377397,
+ +0.741531185599394,
+ +0.949107912342759 };
+ // weights on [-1,1]
+ static const double wts_normal[] = { 0.129484966168870,
+ 0.279705391489277,
+ 0.381830050505119,
+ 0.417959183673469,
+ 0.381830050505119,
+ 0.279705391489277,
+ 0.129484966168870 };
+
+ // points and weights on [0,1]
+ static const double xpts[] = { (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[4]+1)/2.,
+ (xpts_normal[5]+1)/2.,
+ (xpts_normal[6]+1)/2. };
+ static const double wts[] = { wts_normal[0]/2.,
+ wts_normal[1]/2.,
+ wts_normal[2]/2.,
+ wts_normal[3]/2.,
+ wts_normal[4]/2.,
+ wts_normal[5]/2.,
+ wts_normal[6]/2. };
for (unsigned int i=0; i<n_quadrature_points; ++i)
{
QGauss8<1>::QGauss8 () :
Quadrature<1> (8)
{
- static const double G0=0.0198550717512321,
- G1=0.1016667612931866,
- G2=0.2372337950418355,
- G3=0.4082826787521749,
- G4=0.5917173212478251,
- G5=0.7627662049581646,
- G6=0.8983332387068134,
- G7=0.9801449282487679;
- static const double W0=0.0506142681451880,
- W1=0.1111905172266870,
- W2=0.1568533229389435,
- W3=0.1813418916891810;
-
- static const double xpts[] = { G0, G1, G2, G3, G4, G5, G6, G7 };
- static const double wts[] = { W0, W1, W2, W3, W3, W2, W1, W0 };
+ // points on [-1,1]
+ static const double xpts_normal[] = { -0.960289856497536,
+ -0.796666477413627,
+ -0.525532409916329,
+ -0.183434642495650,
+ +0.183434642495650,
+ +0.525532409916329,
+ +0.796666477413627,
+ +0.960289856497536 };
+ // weights on [-1,1]
+ static const double wts_normal[] = { 0.101228536200376,
+ 0.222381034453374,
+ 0.313706645877887,
+ 0.362683783378362,
+ 0.362683783378362,
+ 0.313706645877887,
+ 0.222381034453374,
+ 0.101228536200376 };
+
+ // points and weights on [0,1]
+ static const double xpts[] = { (xpts_normal[0]+1)/2.,
+ (xpts_normal[1]+1)/2.,
+ (xpts_normal[2]+1)/2.,
+ (xpts_normal[3]+1)/2.,
+ (xpts_normal[4]+1)/2.,
+ (xpts_normal[5]+1)/2.,
+ (xpts_normal[6]+1)/2.,
+ (xpts_normal[7]+1)/2. };
+ static const double wts[] = { wts_normal[0]/2.,
+ wts_normal[1]/2.,
+ wts_normal[2]/2.,
+ wts_normal[3]/2.,
+ wts_normal[4]/2.,
+ wts_normal[5]/2.,
+ wts_normal[6]/2.,
+ wts_normal[7]/2. };
for (unsigned int i=0; i<n_quadrature_points; ++i)
{
+
+
+
+
+
+
QMidpoint<1>::QMidpoint () :
Quadrature<1>(1)
{
+QMidpoint<2>::QMidpoint () :
+ Quadrature<2>(1)
+{
+ quadrature_points[0] = Point<2> (1./2., 1./2.);
+ weights[0] = 1.0;
+};
+
+
+
+
QSimpson<1>::QSimpson () :
Quadrature<1> (3)
{
+QSimpson<2>::QSimpson () :
+ Quadrature<2> (9)
+{
+ static const double xpts[] = { 0.0, 0.0, 0.0,
+ 0.5, 0.5, 0.5,
+ 1.0, 1.0, 1.0 };
+ static const double ypts[] = { 0.0, 0.5, 1.0,
+ 0.0, 0.5, 1.0,
+ 0.0, 0.5, 1.0 };
+ static const double wts[] = { 1./6.*1./6., 2./3.*1./6., 1./6.*1./6.,
+ 1./6.*2./3., 2./3.*2./3., 1./6.*2./3.,
+ 1./6.*1./6., 2./3.*1./6., 1./6.*1./6. };
+
+ for (unsigned int i=0; i<n_quadrature_points; ++i)
+ {
+ quadrature_points[i] = Point<2>(xpts[i], ypts[i]);
+ weights[i] = wts[i];
+ };
+};
+
+
+
+
+
+
QTrapez<1>::QTrapez () :
Quadrature<1> (2)
{
+QTrapez<2>::QTrapez () :
+ Quadrature<2> (4)
+{
+ static const double xpts[] = { 0.0, 0.0, 1.0, 1.0 };
+ static const double ypts[] = { 0.0, 1.0, 0.0, 1.0 };
+ static const double wts[] = { 0.25, 0.25, 0.25, 0.25 };
-QGauss4<2>::QGauss4 () :
- Quadrature<2> (4) {
- static const double xpts[] = { 0.211324865, 0.788675135,
- 0.211324865, 0.788675135 };
- static const double ypts[] = { 0.211324865, 0.211324865,
- 0.788675135, 0.788675135 };
- static const double wts[] = { 1./4., 1./4., 1./4., 1./4. };
-
for (unsigned int i=0; i<n_quadrature_points; ++i)
{
quadrature_points[i] = Point<2>(xpts[i], ypts[i]);
weights[i] = wts[i];
};
-
};
+
+
+
+
poisson/poisson: poisson/poisson.o
@echo ================= Linking $@
- @$(CXX) $(CXXFLAGS.g) -o $@ $< ../../mia/control.o $(LIBS.g) -lg++
+ @$(CXX) $(CXXFLAGS.g) -g -o $@ $< ../../mia/control.o $(LIBS.g) -lg++
run: run_grid_test run_dof_test
switch (dim)
{
case 1:
- return ((1-4*3.1415926536*3.1415926536) *
- cos(2*3.1415926536*p(0)));
+// return ((1-4*3.1415926536*3.1415926536) *
+// cos(2*3.1415926536*p(0)));
+ return p(0)*p(0)*p(0)-3./2.*p(0)*p(0)-6*p(0)+3;
case 2:
return ((1-3.1415926536*3.1415926536) *
cos(3.1415926536*p(0)) *
{
for (unsigned int j=0; j<fe_values.total_dofs; ++j)
cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point) +
+ fe_values.shape_grad(j,point) /*+
fe_values.shape_value(i,point) *
- fe_values.shape_value(j,point)) *
+ fe_values.shape_value(j,point)*/) *
fe_values.JxW(point);
rhs[0](i) += fe_values.shape_value(i,point) *
right_hand_side(fe_values.quadrature_point(point)) *
PoissonEquation<2> equation;
QGauss4<2> quadrature;
- HyperBallBoundary<2> boundary(Point<2>(2,3), 4);
+
+// HyperBallBoundary<2> boundary(Point<2>(2,3), 4);
+
+ tria.create_hypercube ();
+// tria.create_hyper_ball(Point<2>(2,3),4);
+// tria.set_boundary (&boundary);
- tria.create_hyper_ball(Point<2>(2,3),4);
- tria.set_boundary (&boundary);
+ tria.refine_global (1);
+ tria.begin_active()->set_refine_flag();
+ tria.execute_refinement ();
- tria.refine_global (5);
dof.distribute_dofs (fe);
problem.assemble (equation, quadrature, fe);
- problem.solve ();
+
+// problem.solve ();
DataOut<2> out;
ofstream gnuplot("gnuplot.out.5");
- problem.fill_data (out);
+ problem.fill_data (out);
out.write_gnuplot (gnuplot);
return 0;
class ostream;
class dSMatrix;
class dSMatrixStruct;
+class dVector;
/**
This class provides two sets of #condense# functions: those taking two
arguments refer to the first possibility above, those taking only one do
their job in-place and refer to the second possibility.
+
+ Condensing vectors works exactly as described above for matrices.
+
+ After solving the condensed system of equations, the solution vector has to
+ be redistributed. This is done by the two #distribute# function, one working
+ with two vectors, one working in-place. The operation of distribution undoes
+ the condensation process in some sense, but it should be noted that it is not
+ the inverse operation.
*/
class ConstraintMatrix {
public:
/**
* Condense a given matrix. The associated
* matrix struct should be condensed and
- * compressed.
+ * compressed. It is the user's
+ * responsibility to guarantee that all
+ * entries in the #condensed# matrix be
+ * zero!
*
* The constraint matrix object must be
* closed to call this function.
void condense (const dSMatrix &uncondensed,
dSMatrix &condensed) const;
-
/**
* This function does much the same as
* the above one, except that it condenses
*/
void condense (dSMatrix &matrix) const;
+ /**
+ * Condense the given vector #uncondensed#
+ * into #condensed#. It is the user's
+ * responsibility to guarantee that all
+ * entries of #condensed# be zero!
+ */
+ void condense (const dVector &uncondensed,
+ dVector &condensed) const;
+
+ /**
+ * Condense the given vector in-place.
+ */
+ void condense (dVector &vec) const;
+
+ /**
+ * Re-distribute the elements of the vector
+ * #condensed# to #uncondensed#. It is the user's
+ * responsibility to guarantee that all
+ * entries of #uncondensed# be zero!
+ *
+ * This function undoes the action of
+ * #condense# somehow, but it should be noted
+ * that it is not the inverse of #condense#-
+ */
+ void distribute (const dVector &condensed,
+ dVector &uncondensed) const;
+
+ /**
+ * Re-distribute the elements of the vector
+ * in-place.
+ */
+ void distribute (dVector &vec) const;
+
/**
* Print the constraint lines. Mainly for
* Exception
*/
DeclException0 (ExcMatrixNotSquare);
-
+ /**
+ * Exception
+ */
+ DeclException0 (ExcWrongDimension);
private:
#include <algo.h>
#include <grid/dof_constraints.h>
#include <lac/dsmatrix.h>
+#include <lac/dvector.h>
#include <iostream.h>
{
// this line is constrained
new_line.push_back (-1);
- // not that #lines# is ordered
+ // note that #lines# is ordered
++next_constraint;
++shift;
}
++next_constraint;
};
+
+ condensed.compress();
};
void ConstraintMatrix::condense (const dSMatrix &uncondensed,
- dSMatrix &condensed) const {};
-void ConstraintMatrix::condense (dSMatrix &uncondensed) const {};
+ dSMatrix &condensed) const {
+ const dSMatrixStruct &uncondensed_struct = uncondensed.get_sparsity_pattern ();
+
+ Assert (sorted == true, ExcMatrixNotClosed());
+ Assert (uncondensed_struct.compressed == true, ExcMatrixNotClosed());
+ Assert (condensed.get_sparsity_pattern().compressed == true, ExcMatrixNotClosed());
+ Assert (uncondensed_struct.n_rows() == uncondensed_struct.n_cols(),
+ ExcMatrixNotSquare());
+ Assert (condensed.n() == condensed.m(),
+ ExcMatrixNotSquare());
+ Assert ((unsigned int)condensed.n()+n_constraints() == (unsigned int)uncondensed.n(),
+ ExcWrongDimension());
+
+ // store for each line of the matrix
+ // its new line number
+ // after compression. If the shift is
+ // -1, this line will be condensed away
+ vector<int> new_line;
+
+ new_line.reserve (uncondensed_struct.n_rows());
+
+ vector<ConstraintLine>::const_iterator next_constraint = lines.begin();
+ unsigned int shift = 0;
+ unsigned int n_rows = (unsigned int)uncondensed_struct.n_rows();
+ for (unsigned int row=0; row!=n_rows; ++row)
+ if (row == (*next_constraint).line)
+ {
+ // this line is constrained
+ new_line.push_back (-1);
+ // note that #lines# is ordered
+ ++next_constraint;
+ ++shift;
+ }
+ else
+ new_line.push_back (row-shift);
+
+
+ next_constraint = lines.begin();
+ for (int row=0; row<uncondensed_struct.rows; ++row)
+ if (new_line[row] != -1)
+ // line not constrained
+ // copy entries if column will not
+ // be condensed away, distribute
+ // otherwise
+ for (int j=uncondensed_struct.rowstart[row]; j<uncondensed_struct.rowstart[row+1]; ++j)
+ if (new_line[uncondensed_struct.colnums[j]] != -1)
+ condensed.add (new_line[row], new_line[uncondensed_struct.colnums[j]],
+ uncondensed.val[j]);
+ else
+ {
+ // let c point to the constraint
+ // of this column
+ vector<ConstraintLine>::const_iterator c = lines.begin();
+ while ((*c).line != (unsigned int)uncondensed_struct.colnums[j]) ++c;
+
+ for (unsigned int q=0; q!=(*c).entries.size(); ++q)
+ // distribute to rows with
+ // appropriate weight
+ condensed.add (new_line[row], new_line[(*c).entries[q].first],
+ uncondensed.val[j] * (*c).entries[q].second);
+ }
+ else
+ // line must be distributed
+ {
+ for (int j=uncondensed_struct.rowstart[row]; j<uncondensed_struct.rowstart[row+1]; ++j)
+ // for each column: distribute
+ if (new_line[uncondensed_struct.colnums[j]] != -1)
+ // column is not constrained
+ for (unsigned int q=0; q!=(*next_constraint).entries.size(); ++q)
+ condensed.add (new_line[(*next_constraint).entries[q].first],
+ new_line[uncondensed_struct.colnums[j]],
+ uncondensed.val[j] * (*next_constraint).entries[q].second);
+
+ else
+ // not only this line but
+ // also this col is constrained
+ {
+ // let c point to the constraint
+ // of this column
+ vector<ConstraintLine>::const_iterator c = lines.begin();
+ while ((*c).line != (unsigned int)uncondensed_struct.colnums[j]) ++c;
+
+ for (unsigned int p=0; p!=(*c).entries.size(); ++p)
+ for (unsigned int q=0; q!=(*next_constraint).entries.size(); ++q)
+ condensed.add (new_line[(*next_constraint).entries[q].first],
+ new_line[(*c).entries[p].first],
+ uncondensed.val[j] *
+ (*next_constraint).entries[q].second *
+ (*c).entries[p].second);
+ };
+
+ ++next_constraint;
+ };
+};
+
+
+
+void ConstraintMatrix::condense (dSMatrix &uncondensed) const {
+ const dSMatrixStruct &sparsity = uncondensed.get_sparsity_pattern ();
+
+ Assert (sorted == true, ExcMatrixNotClosed());
+ Assert (sparsity.compressed == true, ExcMatrixNotClosed());
+ Assert (sparsity.n_rows() == sparsity.n_cols(),
+ ExcMatrixNotSquare());
+
+ // store for each index whether it
+ // must be distributed or not. If entry
+ // is -1, no distribution is necessary.
+ // otherwise, the number states which
+ // line in the constraint matrix handles
+ // this index
+ vector<int> distribute;
+ distribute.reserve (sparsity.n_rows());
+ distribute.insert (distribute.end(), sparsity.n_rows(), -1);
+
+ for (int c=0; c<(signed int)lines.size(); ++c)
+ distribute[lines[c].line] = c;
+
+ int n_rows = sparsity.n_rows();
+ for (int row=0; row<n_rows; ++row)
+ if (distribute[row] == -1)
+ // regular line. loop over cols
+ for (int j=sparsity.rowstart[row]; j<sparsity.rowstart[row+1]; ++j)
+ // end of row reached?
+ if (sparsity.colnums[j] == -1)
+ break;
+ else
+ {
+ if (distribute[sparsity.colnums[j]] != -1)
+ // distribute entry at regular
+ // row #row# and irregular column
+ // sparsity.colnums[j]; set old
+ // entry to zero
+ {
+ for (unsigned int q=0;
+ q!=lines[distribute[sparsity.colnums[j]]].entries.size(); ++q)
+ uncondensed.add (row,
+ lines[distribute[sparsity.colnums[j]]].entries[q].first,
+ uncondensed.val[j] *
+ lines[distribute[sparsity.colnums[j]]].entries[q].second);
+ uncondensed.val[j] = 0.;
+ };
+ }
+ else
+ // row must be distributed
+ for (int j=sparsity.rowstart[row]; j<sparsity.rowstart[row+1]; ++j)
+ // end of row reached?
+ if (sparsity.colnums[j] == -1)
+ break;
+ else
+ {
+ if (distribute[sparsity.colnums[j]] == -1)
+ // distribute entry at irregular
+ // row #row# and regular column
+ // sparsity.colnums[j]. set old
+ // entry to zero
+ {
+ for (unsigned int q=0;
+ q!=lines[distribute[row]].entries.size(); ++q)
+ uncondensed.add (lines[distribute[row]].entries[q].first,
+ sparsity.colnums[j],
+ uncondensed.val[j] *
+ lines[distribute[row]].entries[q].second);
+ uncondensed.val[j] = 0.;
+ }
+ else
+ // distribute entry at irregular
+ // row #row# and irregular column
+ // sparsity.colnums[j]
+ // set old entry to one if on main
+ // diagonal, zero otherwise
+ {
+ for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
+ for (unsigned int q=0;
+ q!=lines[distribute[sparsity.colnums[j]]].entries.size(); ++q)
+ uncondensed.add (lines[distribute[row]].entries[p].first,
+ lines[distribute[sparsity.colnums[j]]].entries[q].first,
+ uncondensed.val[j] *
+ lines[distribute[row]].entries[p].second *
+ lines[distribute[sparsity.colnums[j]]].entries[q].second);
+ uncondensed.val[j] = (row == sparsity.colnums[j] ?
+ 1. : 0. );
+ };
+ };
+};
+
+
+
+
+
+
+void ConstraintMatrix::condense (const dVector &uncondensed,
+ dVector &condensed) const {
+ Assert (sorted == true, ExcMatrixNotClosed());
+ Assert ((unsigned int)condensed.n()+n_constraints() == (unsigned int)uncondensed.n(),
+ ExcWrongDimension());
+
+ // store for each line of the vector
+ // its new line number
+ // after compression. If the shift is
+ // -1, this line will be condensed away
+ vector<int> new_line;
+
+ new_line.reserve (uncondensed.n());
+
+ vector<ConstraintLine>::const_iterator next_constraint = lines.begin();
+ unsigned int shift = 0;
+ unsigned int n_rows = (unsigned int)uncondensed.n();
+ for (unsigned int row=0; row!=n_rows; ++row)
+ if (row == (*next_constraint).line)
+ {
+ // this line is constrained
+ new_line.push_back (-1);
+ // note that #lines# is ordered
+ ++next_constraint;
+ ++shift;
+ }
+ else
+ new_line.push_back (row-shift);
+
+
+ next_constraint = lines.begin();
+ for (int row=0; row<uncondensed.n(); ++row)
+ if (new_line[row] != -1)
+ // line not constrained
+ // copy entry
+ condensed(new_line[row]) += uncondensed(row);
+
+ else
+ // line must be distributed
+ {
+ for (unsigned int q=0; q!=(*next_constraint).entries.size(); ++q)
+ condensed(new_line[(*next_constraint).entries[q].first])
+ +=
+ uncondensed(row) * (*next_constraint).entries[q].second;
+
+ ++next_constraint;
+ };
+};
+
+
+
+void ConstraintMatrix::condense (dVector &vec) const {
+ vector<ConstraintLine>::const_iterator next_constraint = lines.begin();
+ for (unsigned int row=0; row<(unsigned int)vec.n(); ++row)
+ if (row == (*next_constraint).line)
+ // line must be distributed
+ {
+ for (unsigned int q=0; q!=(*next_constraint).entries.size(); ++q)
+ vec((*next_constraint).entries[q].first)
+ +=
+ vec(row) * (*next_constraint).entries[q].second;
+ // set entry to zero
+ vec(row) = 0.;
+
+ ++next_constraint;
+ };
+};
+
poisson/poisson: poisson/poisson.o
@echo ================= Linking $@
- @$(CXX) $(CXXFLAGS.g) -o $@ $< ../../mia/control.o $(LIBS.g) -lg++
+ @$(CXX) $(CXXFLAGS.g) -g -o $@ $< ../../mia/control.o $(LIBS.g) -lg++
run: run_grid_test run_dof_test
switch (dim)
{
case 1:
- return ((1-4*3.1415926536*3.1415926536) *
- cos(2*3.1415926536*p(0)));
+// return ((1-4*3.1415926536*3.1415926536) *
+// cos(2*3.1415926536*p(0)));
+ return p(0)*p(0)*p(0)-3./2.*p(0)*p(0)-6*p(0)+3;
case 2:
return ((1-3.1415926536*3.1415926536) *
cos(3.1415926536*p(0)) *
{
for (unsigned int j=0; j<fe_values.total_dofs; ++j)
cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point) +
+ fe_values.shape_grad(j,point) /*+
fe_values.shape_value(i,point) *
- fe_values.shape_value(j,point)) *
+ fe_values.shape_value(j,point)*/) *
fe_values.JxW(point);
rhs[0](i) += fe_values.shape_value(i,point) *
right_hand_side(fe_values.quadrature_point(point)) *
PoissonEquation<2> equation;
QGauss4<2> quadrature;
- HyperBallBoundary<2> boundary(Point<2>(2,3), 4);
+
+// HyperBallBoundary<2> boundary(Point<2>(2,3), 4);
+
+ tria.create_hypercube ();
+// tria.create_hyper_ball(Point<2>(2,3),4);
+// tria.set_boundary (&boundary);
- tria.create_hyper_ball(Point<2>(2,3),4);
- tria.set_boundary (&boundary);
+ tria.refine_global (1);
+ tria.begin_active()->set_refine_flag();
+ tria.execute_refinement ();
- tria.refine_global (5);
dof.distribute_dofs (fe);
problem.assemble (equation, quadrature, fe);
- problem.solve ();
+
+// problem.solve ();
DataOut<2> out;
ofstream gnuplot("gnuplot.out.5");
- problem.fill_data (out);
+ problem.fill_data (out);
out.write_gnuplot (gnuplot);
return 0;