From 7499f781ac2fd91b0e5c48538bb122420ab6921d Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 17 Mar 1998 12:44:10 +0000 Subject: [PATCH] More support for hanging nodes and more structure for quadrature rules. git-svn-id: https://svn.dealii.org/trunk@72 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/quadrature.h | 10 +- deal.II/base/include/base/quadrature_lib.h | 87 ++- deal.II/base/source/quadrature_lib.cc | 581 ++++++++++++++++-- deal.II/deal.II/Attic/examples/Makefile | 2 +- .../deal.II/Attic/examples/poisson/poisson.cc | 27 +- .../deal.II/include/dofs/dof_constraints.h | 53 +- .../deal.II/source/dofs/dof_constraints.cc | 266 +++++++- tests/big-tests/Makefile | 2 +- tests/big-tests/poisson/poisson.cc | 27 +- 9 files changed, 977 insertions(+), 78 deletions(-) diff --git a/deal.II/base/include/base/quadrature.h b/deal.II/base/include/base/quadrature.h index 3f432333c8..cf56a1914d 100644 --- a/deal.II/base/include/base/quadrature.h +++ b/deal.II/base/include/base/quadrature.h @@ -16,7 +16,15 @@ 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# 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 class Quadrature { public: diff --git a/deal.II/base/include/base/quadrature_lib.h b/deal.II/base/include/base/quadrature_lib.h index 5d3a147ea7..bfd620fc98 100644 --- a/deal.II/base/include/base/quadrature_lib.h +++ b/deal.II/base/include/base/quadrature_lib.h @@ -8,6 +8,11 @@ #include +/** + Second order Gauss quadrature formula. + + Reference: Ward Cheney, David Kincaid: Numerical Mathematics and Computing. +*/ template class QGauss2 : public Quadrature { public: @@ -16,14 +21,24 @@ class QGauss2 : public Quadrature { +/** + Third order Gauss quadrature formula. + + Reference: Ward Cheney, David Kincaid: Numerical Mathematics and Computing. +*/ template -class QGauss2x4 : public Quadrature { +class QGauss3 : public Quadrature { public: - QGauss2x4 (); + QGauss3 (); }; +/** + Fourth order Gauss quadrature formula. + + Reference: Ward Cheney, David Kincaid: Numerical Mathematics and Computing. +*/ template class QGauss4 : public Quadrature { public: @@ -32,6 +47,63 @@ class QGauss4 : public Quadrature { + +/** + Fifth order Gauss quadrature formula. + + Reference: Ward Cheney, David Kincaid: Numerical Mathematics and Computing. +*/ +template +class QGauss5 : public Quadrature { + 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 +class QGauss6 : public Quadrature { + 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 +class QGauss7 : public Quadrature { + 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 class QGauss8 : public Quadrature { public: @@ -40,6 +112,11 @@ class QGauss8 : public Quadrature { + + +/** + First order midpoint quadrature rule. +*/ template class QMidpoint : public Quadrature { public: @@ -48,6 +125,9 @@ class QMidpoint : public Quadrature { +/** + Simpson quadrature rule. +*/ template class QSimpson : public Quadrature { public: @@ -56,6 +136,9 @@ class QSimpson : public Quadrature { +/** + Trapezoidal quadrature rule. +*/ template class QTrapez : public Quadrature { public: diff --git a/deal.II/base/source/quadrature_lib.cc b/deal.II/base/source/quadrature_lib.cc index 44c574b6b6..0b38e6e7d2 100644 --- a/deal.II/base/source/quadrature_lib.cc +++ b/deal.II/base/source/quadrature_lib.cc @@ -1,14 +1,23 @@ /* $Id$ */ #include +#include 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::QGauss2 () : -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(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::QGauss2x4 () : +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(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(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(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(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(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(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::QGauss4 () : 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::QGauss8 () : + + + + + + QMidpoint<1>::QMidpoint () : Quadrature<1>(1) { @@ -102,6 +553,16 @@ QMidpoint<1>::QMidpoint () : +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) { @@ -117,6 +578,31 @@ QSimpson<1>::QSimpson () : +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(xpts[i], ypts[i]); + weights[i] = wts[i]; + }; +}; + + + + + + QTrapez<1>::QTrapez () : Quadrature<1> (2) { @@ -132,19 +618,20 @@ QTrapez<1>::QTrapez () : +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(xpts[i], ypts[i]); weights[i] = wts[i]; }; - }; + + + + diff --git a/deal.II/deal.II/Attic/examples/Makefile b/deal.II/deal.II/Attic/examples/Makefile index 5f54731868..f40c158b70 100644 --- a/deal.II/deal.II/Attic/examples/Makefile +++ b/deal.II/deal.II/Attic/examples/Makefile @@ -35,7 +35,7 @@ dof/dof_test: dof/dof_test.o 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 diff --git a/deal.II/deal.II/Attic/examples/poisson/poisson.cc b/deal.II/deal.II/Attic/examples/poisson/poisson.cc index 9be2e12b2e..7b6b860d31 100644 --- a/deal.II/deal.II/Attic/examples/poisson/poisson.cc +++ b/deal.II/deal.II/Attic/examples/poisson/poisson.cc @@ -48,8 +48,9 @@ double PoissonEquation::right_hand_side (const Point &p) const { 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)) * @@ -91,9 +92,9 @@ void PoissonEquation<2>::assemble (dFMatrix &cell_matrix, { for (unsigned int j=0; j 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; diff --git a/deal.II/deal.II/include/dofs/dof_constraints.h b/deal.II/deal.II/include/dofs/dof_constraints.h index a084417d2a..ec9b0a975d 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.h @@ -14,6 +14,7 @@ class ostream; class dSMatrix; class dSMatrixStruct; +class dVector; /** @@ -83,6 +84,14 @@ class dSMatrixStruct; 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: @@ -178,7 +187,10 @@ class ConstraintMatrix { /** * 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. @@ -186,7 +198,6 @@ class ConstraintMatrix { void condense (const dSMatrix &uncondensed, dSMatrix &condensed) const; - /** * This function does much the same as * the above one, except that it condenses @@ -196,6 +207,39 @@ class ConstraintMatrix { */ 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 @@ -238,7 +282,10 @@ class ConstraintMatrix { * Exception */ DeclException0 (ExcMatrixNotSquare); - + /** + * Exception + */ + DeclException0 (ExcWrongDimension); private: diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc index 525def890f..e250b7a28f 100644 --- a/deal.II/deal.II/source/dofs/dof_constraints.cc +++ b/deal.II/deal.II/source/dofs/dof_constraints.cc @@ -4,6 +4,7 @@ #include #include #include +#include #include @@ -106,7 +107,7 @@ void ConstraintMatrix::condense (const dSMatrixStruct &uncondensed, { // this line is constrained new_line.push_back (-1); - // not that #lines# is ordered + // note that #lines# is ordered ++next_constraint; ++shift; } @@ -162,6 +163,8 @@ void ConstraintMatrix::condense (const dSMatrixStruct &uncondensed, ++next_constraint; }; + + condensed.compress(); }; @@ -236,8 +239,265 @@ void ConstraintMatrix::condense (dSMatrixStruct &sparsity) const { 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 new_line; + + new_line.reserve (uncondensed_struct.n_rows()); + + vector::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::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::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 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 new_line; + + new_line.reserve (uncondensed.n()); + + vector::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::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; + }; +}; + diff --git a/tests/big-tests/Makefile b/tests/big-tests/Makefile index 5f54731868..f40c158b70 100644 --- a/tests/big-tests/Makefile +++ b/tests/big-tests/Makefile @@ -35,7 +35,7 @@ dof/dof_test: dof/dof_test.o 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 diff --git a/tests/big-tests/poisson/poisson.cc b/tests/big-tests/poisson/poisson.cc index 9be2e12b2e..7b6b860d31 100644 --- a/tests/big-tests/poisson/poisson.cc +++ b/tests/big-tests/poisson/poisson.cc @@ -48,8 +48,9 @@ double PoissonEquation::right_hand_side (const Point &p) const { 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)) * @@ -91,9 +92,9 @@ void PoissonEquation<2>::assemble (dFMatrix &cell_matrix, { for (unsigned int j=0; j 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; -- 2.39.5