]> https://gitweb.dealii.org/ - dealii.git/commitdiff
More support for hanging nodes and more structure for quadrature rules.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 17 Mar 1998 12:44:10 +0000 (12:44 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 17 Mar 1998 12:44:10 +0000 (12:44 +0000)
git-svn-id: https://svn.dealii.org/trunk@72 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/quadrature.h
deal.II/base/include/base/quadrature_lib.h
deal.II/base/source/quadrature_lib.cc
deal.II/deal.II/Attic/examples/Makefile
deal.II/deal.II/Attic/examples/poisson/poisson.cc
deal.II/deal.II/include/dofs/dof_constraints.h
deal.II/deal.II/source/dofs/dof_constraints.cc
tests/big-tests/Makefile
tests/big-tests/poisson/poisson.cc

index 3f432333c8994bad1c05faf2c67eee1179f15708..cf56a1914df1bb4d2f4c977c1a4b4042c900d267 100644 (file)
   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:
index 5d3a147ea7993816c8b808fd6f82a41923fc3eaa..bfd620fc98cc1318dc72c6d7c3b91a9a409d6e25 100644 (file)
@@ -8,6 +8,11 @@
 #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:
@@ -16,14 +21,24 @@ class QGauss2 : public Quadrature<dim> {
 
 
 
+/**
+   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:
@@ -32,6 +47,63 @@ class QGauss4 : public Quadrature<dim> {
 
 
 
+
+/**
+   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:
@@ -40,6 +112,11 @@ class QGauss8 : public Quadrature<dim> {
 
 
 
+
+
+/**
+   First order midpoint quadrature rule.
+*/
 template <int dim>
 class QMidpoint : public Quadrature<dim> {
   public:
@@ -48,6 +125,9 @@ class QMidpoint : public Quadrature<dim> {
 
 
 
+/**
+   Simpson quadrature rule.
+*/
 template <int dim>
 class QSimpson : public Quadrature<dim> {
   public:
@@ -56,6 +136,9 @@ class QSimpson : public Quadrature<dim> {
 
 
 
+/**
+   Trapezoidal quadrature rule.
+*/
 template <int dim>
 class QTrapez : public Quadrature<dim> {
   public:
index 44c574b6b67ae9fa514be27c404589b54426626b..0b38e6e7d24af4bec82310c5381232e8ef30de85 100644 (file)
@@ -1,14 +1,23 @@
 /* $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) 
     {
@@ -19,20 +28,58 @@ QGauss2<1>::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<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) 
     {
@@ -43,18 +90,395 @@ QGauss2x4<1>::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<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) 
     {
@@ -68,21 +492,42 @@ QGauss4<1>::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<n_quadrature_points; ++i) 
     {
@@ -93,6 +538,12 @@ QGauss8<1>::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<n_quadrature_points; ++i) 
+    {
+      quadrature_points[i] = Point<2>(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<n_quadrature_points; ++i) 
     {
       quadrature_points[i] = Point<2>(xpts[i], ypts[i]);
       weights[i] = wts[i];
     };
-  
 };
+
+
+
+
index 5f54731868392616f702fc2020de110f56bb7255..f40c158b70b349ddaec9a4e098df28930c4bdc76 100644 (file)
@@ -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
index 9be2e12b2e7ebbd10b593cdbef1f8cd1e19e35ac..7b6b860d3191643d3c1628fe898dd4ad3cb6f683 100644 (file)
@@ -48,8 +48,9 @@ double PoissonEquation<dim>::right_hand_side (const Point<dim> &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<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)) *
@@ -112,20 +113,26 @@ int main () {
   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;
index a084417d2a9b5e4b81c211d6fc2b975e97e420b0..ec9b0a975d5a5f32170233d9243936b6a08b2ddf 100644 (file)
@@ -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:
 
index 525def890fe81116f7eec56d140f83cb82209ac0..e250b7a28f7973c86d3085c7330bc1700f718e0b 100644 (file)
@@ -4,6 +4,7 @@
 #include <algo.h>
 #include <grid/dof_constraints.h>
 #include <lac/dsmatrix.h>
+#include <lac/dvector.h>
 #include <iostream.h>
 
 
@@ -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<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;
+      };
+};
+  
 
 
 
index 5f54731868392616f702fc2020de110f56bb7255..f40c158b70b349ddaec9a4e098df28930c4bdc76 100644 (file)
@@ -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
index 9be2e12b2e7ebbd10b593cdbef1f8cd1e19e35ac..7b6b860d3191643d3c1628fe898dd4ad3cb6f683 100644 (file)
@@ -48,8 +48,9 @@ double PoissonEquation<dim>::right_hand_side (const Point<dim> &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<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)) *
@@ -112,20 +113,26 @@ int main () {
   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;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.