]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Mostly bugfixes in the assembler and fe part.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 11 Mar 1998 16:14:20 +0000 (16:14 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 11 Mar 1998 16:14:20 +0000 (16:14 +0000)
git-svn-id: https://svn.dealii.org/trunk@56 0785d39b-7218-0410-832d-ea1e28bc413d

20 files changed:
deal.II/base/include/base/quadrature.h
deal.II/base/source/quadrature.cc [new file with mode: 0644]
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/Todo
deal.II/deal.II/include/dofs/dof_constraints.h
deal.II/deal.II/include/dofs/dof_handler.h
deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/numerics/assembler.h
deal.II/deal.II/include/numerics/base.h
deal.II/deal.II/source/dofs/dof_constraints.cc
deal.II/deal.II/source/dofs/dof_handler.cc
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_lib.linear.cc
deal.II/deal.II/source/numerics/assembler.cc
deal.II/deal.II/source/numerics/base.cc
tests/big-tests/Makefile
tests/big-tests/Makefile.dep
tests/big-tests/poisson/poisson.cc

index 7f06f7c545674f6d8cc0848e77ce0daa0b1c610f..3f432333c8994bad1c05faf2c67eee1179f15708 100644 (file)
@@ -41,6 +41,14 @@ class Quadrature {
                                      */
     double weight (const unsigned int i) const;
 
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcInvalidIndex,
+                   int, int,
+                   << "The index " << arg1
+                   << " is out of range, it should be less than " << arg2);
+
   protected:
                                     /**
                                      * List of quadrature points. To be filled
diff --git a/deal.II/base/source/quadrature.cc b/deal.II/base/source/quadrature.cc
new file mode 100644 (file)
index 0000000..187cd73
--- /dev/null
@@ -0,0 +1,34 @@
+/* $Id$ */
+
+#include <fe/quadrature.h>
+
+
+template <int dim>
+Quadrature<dim>::Quadrature (const unsigned int n_q) :
+               n_quadrature_points(n_q),
+               quadrature_points (n_q, Point<dim>()),
+               weights (n_q, 0) {};
+
+
+
+template <int dim>
+const Point<dim> & Quadrature<dim>::quad_point (const unsigned int i) const {
+  Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
+  return quadrature_points[i];
+};
+
+
+
+template <int dim>
+double Quadrature<dim>::weight (const unsigned int i) const {
+  Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
+  return weights[i];
+};
+
+
+
+
+
+// explicite instantiations
+template class Quadrature<1>;
+template class Quadrature<2>;
index 0f21f502fbab79b0bff2348637ec525c5799e331..44c574b6b67ae9fa514be27c404589b54426626b 100644 (file)
@@ -12,8 +12,8 @@ QGauss2<1>::QGauss2 () :
 
   for (unsigned int i=0; i<n_quadrature_points; ++i) 
     {
-      quadrature_points.push_back (Point<1> (xpts[i]));
-      weights.push_back (wts[i]);
+      quadrature_points[i] = Point<1>(xpts[i]);
+      weights[i] = wts[i];
     };
 };
 
@@ -36,8 +36,8 @@ QGauss2x4<1>::QGauss2x4 () :
 
   for (unsigned int i=0; i<n_quadrature_points; ++i) 
     {
-      quadrature_points.push_back (Point<1> (xpts[i]));
-      weights.push_back (wts[i]);
+      quadrature_points[i] = Point<1>(xpts[i]);
+      weights[i] = wts[i];
     };
 };
 
@@ -58,8 +58,8 @@ QGauss4<1>::QGauss4 () :
 
   for (unsigned int i=0; i<n_quadrature_points; ++i) 
     {
-      quadrature_points.push_back (Point<1> (xpts[i]));
-      weights.push_back (wts[i]);
+      quadrature_points[i] = Point<1>(xpts[i]);
+      weights[i] = wts[i];
     };
 };
 
@@ -86,8 +86,8 @@ QGauss8<1>::QGauss8 () :
 
   for (unsigned int i=0; i<n_quadrature_points; ++i) 
     {
-      quadrature_points.push_back (Point<1> (xpts[i]));
-      weights.push_back (wts[i]);
+      quadrature_points[i] = Point<1>(xpts[i]);
+      weights[i] = wts[i];
     };
 };
 
@@ -96,8 +96,8 @@ QGauss8<1>::QGauss8 () :
 QMidpoint<1>::QMidpoint () :
                Quadrature<1>(1)
 {
-  quadrature_points.push_back (0.5);
-  weights.push_back (1.0);
+  quadrature_points[0] = 0.5;
+  weights[0] = 1.0;
 };
 
 
@@ -110,8 +110,8 @@ QSimpson<1>::QSimpson () :
 
   for (unsigned int i=0; i<n_quadrature_points; ++i) 
     {
-      quadrature_points.push_back (Point<1> (xpts[i]));
-      weights.push_back (wts[i]);
+      quadrature_points[i] = Point<1>(xpts[i]);
+      weights[i] = wts[i];
     };
 };
 
@@ -125,7 +125,26 @@ QTrapez<1>::QTrapez () :
 
   for (unsigned int i=0; i<n_quadrature_points; ++i) 
     {
-      quadrature_points.push_back (Point<1> (xpts[i]));
-      weights.push_back (wts[i]);
+      quadrature_points[i] = Point<1>(xpts[i]);
+      weights[i] = wts[i];
     };
 };
+
+
+
+
+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 719e01c5d08b40bfd3613bf5c251b51a9fc4bf3c..e874cc83e65b4515b61b83ef3cb43343a45789cb 100644 (file)
@@ -5,8 +5,9 @@ INCPATH   = -I../include -I../../lac/include -I../../base/include
 CXXFLAGS.g= -DDEBUG -g -Wall $(INCPATH)
 CXXFLAGS  = -O3 $(INCPATH)
 
-LIBS.g    = -L../lib -lbasic.g -lgrid.g -lfe.g -lnumerics.g -llac.g
-LIBS      = -L../lib -lbasic -lgrid -lfe -lnumerics 
+LIBPATH   = -L../lib -L../../lac/lib
+LIBS.g    = $(LIBPATH) -lbasic.g -lgrid.g -lfe.g -lnumerics.g -llac.g
+LIBS      = $(LIBPATH) -lbasic -lgrid -lfe -lnumerics 
 
 ifeq ($(shell uname),SunOS)
 CXXFLAGS  := -V2.7.2.3 $(CXXFLAGS)
@@ -14,17 +15,22 @@ CXXFLAGS.g:= -V2.7.2.3 $(CXXFLAGS.g)
 endif
 
 
-cc-files = grid/grid_test.cc
+cc-files = grid/grid_test.cc dof/dof_test.cc poisson/poisson.cc
 o-files  = $(cc-files:.cc=.o)
 h-files  = $(wildcard ../include/*.h)
 
-examples: grid/grid_test dof/dof_test
+examples: $(cc-files:.cc=)
 
+%.o : %.cc
+       $(CXX) $(CXXFLAGS) -o $@ -c $<
 
-grid/grid_test: grid/grid_test.o ../lib/libgrid.a
+grid/grid_test: grid/grid_test.o
        $(CXX) $(CXXFLAGS) -o $@ $< $(LIBS)
 
-dof/dof_test: dof/dof_test.o ../lib/libgrid.g.a ../lib/libdeal.g.a
+dof/dof_test: dof/dof_test.o
+       $(CXX) $(CXXFLAGS.g) -o $@ $< $(LIBS.g) -lg++
+
+poisson/poisson: poisson/poisson.o
        $(CXX) $(CXXFLAGS.g) -o $@ $< $(LIBS.g) -lg++
 
 
@@ -52,6 +58,6 @@ clean:
 #to remake this file upon inclusion at the bottom
 #of this file.
 Makefile.dep: $(cc-files) $(h-files)
-       $(CXX) -I../include -M $(cc-files) > Makefile.dep
+       $(CXX) $(INCPATH) -M $(cc-files) > Makefile.dep
 
 include Makefile.dep
index f4f0530b0581d7492e0837f8fd2955f516c7a3ad..d6d6321fd973bf188e3adb634a5acf3ea71ea1af 100644 (file)
@@ -4,13 +4,15 @@
 #include <grid/dof.h>
 #include <grid/tria.h>
 #include <fe/fe_lib.h>
-#include <grid/tria_iterator.h>
 #include <grid/tria_accessor.h>
+#include <grid/dof_accessor.h>
+#include <grid/tria_iterator.h>
 #include <lac/dsmatrix.h>
 #include <grid/dof_constraints.h>
 #include <basic/data_io.h>
 #include <numerics/base.h>
 #include <numerics/assembler.h>
+#include <fe/quadrature_lib.h>
 
 #include <fstream.h>
 #include <cmath>
@@ -20,6 +22,7 @@ extern "C" {
 
 extern TriaActiveIterator<1,CellAccessor<1> > x;
 extern TriaActiveIterator<2,CellAccessor<2> > y;
+extern TriaRawIterator<1,DoFLineAccessor<1,LineAccessor<1> > > z;
 
 
 
@@ -33,6 +36,16 @@ class PoissonEquation :  public Equation<dim> {
                           vector<dVector>     &rhs,
                           const FEValues<dim> &fe_values,
                           const Triangulation<dim>::cell_iterator &cell) const;
+    double right_hand_side (const Point<dim> &) const;
+};
+
+
+
+template <int dim>
+inline
+double PoissonEquation<dim>::right_hand_side (const Point<dim> &) const {
+                                  // this yields as solution u=x^2+y^2+...=\vec x^2
+  return -2*dim;
 };
 
 
@@ -45,26 +58,49 @@ void PoissonEquation<1>::assemble (dFMatrix            &cell_matrix,
     for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
       {
        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.JxW(point);
-       rhs[0] = fe_values.shape_value(i,point) *
-                right_hand_side(fe_values.quadrature_point(point)) *
-                fe_values.JxW(point);
+         cell_matrix(i,j) += fe_values.shape_grad(i,point) *
+                             fe_values.shape_grad(j,point) *
+                             fe_values.JxW(point);
+       rhs[0](i) += fe_values.shape_value(i,point) *
+                    right_hand_side(fe_values.quadrature_point(point)) *
+                    fe_values.JxW(point);
+      };
+};
+
+
+
+void PoissonEquation<2>::assemble (dFMatrix            &cell_matrix,
+                                  vector<dVector>     &rhs,
+                                  const FEValues<2>   &fe_values,
+                                  const Triangulation<2>::cell_iterator &) const {
+  for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+    for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+      {
+       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.JxW(point);
+       rhs[0](i) += fe_values.shape_value(i,point) *
+                    right_hand_side(fe_values.quadrature_point(point)) *
+                    fe_values.JxW(point);
       };
 };
 
   
 
 
-void main () {
-  Triangulation<1> tria;
-  DoFHandler<1>    dof(tria);
-  FELinear<1>      fe;
-  ProblemBase<1>   problem(&tria, &dof);
-  Quadrature<1>    q;
+int main () {
+  Triangulation<2>   tria;
+  DoFHandler<2>      dof(&tria);
+  FELinear<2>        fe;
+  ProblemBase<2>     problem(&tria, &dof);
+  PoissonEquation<2> equation;
+  QGauss4<2>         quadrature;
   
-  tria.create_hypercube();
+  tria.create_unit_hypercube();
+  tria.refine_global (1);
   dof.distribute_dofs (fe);
-  problem.assemble (q);
+  problem.assemble (equation, quadrature, fe);
+  
+  return 0;
 };
index eedc18636659eea8481cdb354159fb5185c36549..42d502cd440c5e5b52b7f12e31045121173ef1f8 100644 (file)
@@ -34,12 +34,19 @@ Unify CellAccessor<1> and <2> by renaming
 
 Add support for fast assemblage of mass matrices.
 
+Make AssemblerData a local class to Assembler again if gcc2.8 supports it.
 
 
 
 DEAL:
 make dSMatrixStruct warn if max # of entries per row is overrun!
 
+let dSMatrixStruct have an empty constructor to allow for member
+  objects being initialized later
+
 why do loops in dFMatrix run backwards?
 
 Let an empty dfmatrix assume it has dimensions zero rather than one!
+
+Is dvector::operator= (double) really useful or does it make more
+  confusion than it helps?
index 305b46f1c4c3a153b0ddb1e1381205bd4ca7866f..a084417d2a9b5e4b81c211d6fc2b975e97e420b0 100644 (file)
@@ -12,6 +12,7 @@
 
 
 class ostream;
+class dSMatrix;
 class dSMatrixStruct;
 
 
@@ -173,6 +174,28 @@ class ConstraintMatrix {
                                      * end of the function.
                                      */
     void condense (dSMatrixStruct &sparsity) const;
+
+                                    /**
+                                     * Condense a given matrix. The associated
+                                     * matrix struct should be condensed and
+                                     * compressed.
+                                     *
+                                     * 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
+                                     * the matrix 'in-place'. See the general
+                                     * documentation of this class for more
+                                     * detailed information.
+                                     */
+    void condense (dSMatrix &matrix) const;
+
     
                                     /**
                                      * Print the constraint lines. Mainly for
index 8a72bd24151e63ea26ca1656ae4e58d15f167d31..78a868da31c12ae8f78393dd57df784093caa972 100644 (file)
@@ -469,9 +469,10 @@ class DoFHandler : public DoFDimensionInfo<dim> {
                                      * one dimension, the functions returns
                                      * immediately after clearing the
                                      * constraint matrix.
-                                     *
                                      * For more than one dimension, the matrix
-                                     * is cleared before usage.
+                                     * is cleared before usage. The constraint
+                                     * matrix is closed anyway, no matter of the
+                                     * dimension.
                                      *
                                      * To condense a given sparsity pattern,
                                      * use #ConstraintMatrix::condense#.
index ca758eb7e6207a0c7ba05a71c49f69b5cbbcb976..d23143ef7e9de8a6681c6d65f761ca8ce7777288 100644 (file)
@@ -354,6 +354,10 @@ class FiniteElementBase {
                                      * Exception
                                      */
     DeclException0 (ExcPureFunctionCalled);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcNotImplemented);
     
   protected:
                                     /**
index 4f6c6ac9876d8dd15252eec867a2e20f66cf93a4..a7d2e036c1f83997024bf97c1484323258142419 100644 (file)
@@ -89,6 +89,62 @@ class Equation {
 
 
 
+/**
+ * Structure to be passed upon
+ * construction of an assembler object
+ * through the iterator object. See
+ * \Ref{TriaRawIterator} for a discussion
+ * of this mechanism.
+ */
+template <int dim>
+struct AssemblerData {
+                                    /**
+                                     * Constructor.
+                                     */
+    AssemblerData (const DoFHandler<dim>    &dof,
+                  const bool                assemble_matrix,
+                  const bool                assemble_rhs,
+                  const unsigned int        n_rhs,
+                  ProblemBase<dim>         &problem,
+                  const Quadrature<dim>    &quadrature,
+                  const FiniteElement<dim> &fe);
+    
+                                    /**
+                                     * Pointer to the dof handler object
+                                     * to be used to iterate on.
+                                     */
+    const DoFHandler<dim>  &dof;
+                                    /**
+                                     * Flags whether to assemble the matrix
+                                     * and right hand sides.
+                                     */
+    const bool              assemble_matrix, assemble_rhs;
+                                    /**
+                                     * How many right hand sides are there.
+                                     */
+    const unsigned int      n_rhs;
+                                    /**
+                                     * Pointer to the #ProblemBase# object
+                                     * of which the global matrices and
+                                     * vectors are to be assembled.
+                                     */
+    ProblemBase<dim>       &problem;
+                                    /**
+                                     * Pointer to a quadrature object to be
+                                     * used for this assemblage process.
+                                     */
+    const Quadrature<dim>  &quadrature;
+                                    /**
+                                     * Use this FE type for the assemblage
+                                     * process. The FE object must match that
+                                     * used to construct the #DoFHandler#
+                                     * object.
+                                     */
+    const FiniteElement<dim> &fe;
+};
+
+
+
 
 /**
   An #Assembler# is a specialized version of a #DoFCellAccessor# which adds
@@ -97,38 +153,6 @@ class Equation {
 template <int dim>
 class Assembler : public DoFCellAccessor<dim> {
   public:
-                                    /**
-                                     * Structure to be passed upon
-                                     * construction of an assembler object.
-                                     */
-    struct AssemblerData {
-                                        /**
-                                         * Pointer to the dof handler object
-                                         * to be used to iterate on.
-                                         */
-       DoFHandler<dim>  *dof;
-                                        /**
-                                         * Flags whether to assemble the matrix
-                                         * and right hand sides.
-                                         */
-       bool              assemble_matrix, assemble_rhs;
-                                        /**
-                                         * How many right hand sides are there.
-                                         */
-       unsigned int      n_rhs;
-                                        /**
-                                         * Pointer to the #ProblemBase# object
-                                         * of which the global matrices and
-                                         * vectors are to be assembled.
-                                         */
-       ProblemBase<dim> *problem;
-                                        /**
-                                         * Pointer to a quadrature object to be
-                                         * used for this assemblage process.
-                                         */
-       Quadrature<dim>  *quadrature;
-    };
-
                                     /**
                                      * Default constructor, unused thus not
                                      * implemented.
@@ -196,7 +220,23 @@ class Assembler : public DoFCellAccessor<dim> {
                                      * of which we are to use the matrices
                                      * and vectors.
                                      */
-    ProblemBase<dim> *problem;
+    ProblemBase<dim> &problem;
+
+                                    /**
+                                     * Pointer to the finite element used for
+                                     * the assemblage process. We store a
+                                     * pointer to the finite element, not a
+                                     * copy, but this is not a problem, since
+                                     * #Assembler# objects are usually created
+                                     * and deleted by the
+                                     * #ProblemBase::assemble# routine, which
+                                     * itself gets passed the finite element;
+                                     * the lifetime of the #Assembler# object
+                                     * is thus less than the lifetime of the
+                                     * finite element object and storing
+                                     * a pointer is not risky.
+                                     */
+    const FiniteElement<dim> &fe;
 
                                     /**
                                      * The finite element evaluated at the
index a55cb9098944763c705e6554702d41fb0c32ce37..ec9f7e6548cda6d213ee20ecc34ef60aa5450eb9 100644 (file)
@@ -6,7 +6,7 @@
 
 
 #include <lac/dsmatrix.h>
-
+#include <base/exceptions.h>
 
 // forward declaration
 template <int dim> class Triangulation;
@@ -17,24 +17,103 @@ template <int dim> class Quadrature;
 
 
 
+/**
+  Base class for user problems. This class stores the system matrix and right
+  hand side vectors as well as a solution vector. It initiates the assemblage
+  process of matrix and vectors and so on.
+
+  This class is not extremely versatile as could certainly be. For example
+  it presently only supports sparse matrices and has no multigrid features.
+  However, all these things depend strongly on the problem and it seems
+  best to implement many of these things yourself. Thus, this class is more
+  a display of concept haw to work with deal.II.
+
 
+  {\bf Assemblage}
+
+  The #assemble# member function does the assemblage of the system matrix and
+  the given number of right hand sides. It does the following steps:
+  \begin{itemize}
+    \item Create sparsity pattern of the system matrix and condense it with
+      the constraints induced by hanging nodes.
+    \item Initialize an assembler object.
+    \item Loop over all cells and assemble matrix and vectors using the given
+      quadrature formula and the equation object which contains the weak
+      formulation of the equation.
+    \item Condense the system matrix with the constraints induced by hanging
+      nodes.
+  \end{itemize}
+  */
 template <int dim>
 class ProblemBase {
   public:
+                                    /**
+                                     * Constructor. Use this triangulation and
+                                     * degree of freedom object during the
+                                     * lifetime of this object. The dof
+                                     * object must refer to the given
+                                     * triangulation. The number of right hand
+                                     * sides being assembled defaults to one.
+                                     */
     ProblemBase (Triangulation<dim> *tria,
-                DoFHandler<dim>    *dof_handler);
-
-    virtual void assemble (const Equation<dim> &equation,
-                          const Quadrature<dim> &q);
-
+                DoFHandler<dim>    *dof_handler,
+                const unsigned int  n_rhs=1);
+
+                                    /**
+                                     * Initiate the process of assemblage of
+                                     * vectors and system matrix. Use the
+                                     * given equation object and the given
+                                     * quadrature formula.
+                                     *
+                                     * For what exactly happens here, refer to
+                                     * the general doc of this class.
+                                     */
+    virtual void assemble (const Equation<dim>      &equation,
+                          const Quadrature<dim>    &q,
+                          const FiniteElement<dim> &fe);
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcDofAndTriaDontMatch);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcNoMemory);
+    
   protected:
+                                    /**
+                                     * Pointer to the triangulation to work on.
+                                     */
     Triangulation<dim> *tria;
+
+                                    /**
+                                     * Pointer to the degree of freedom handler
+                                     * to work with.
+                                     */
     DoFHandler<dim>    *dof_handler;
 
-    dSMatrix            system_matrix;
+                                    /**
+                                     * Sparsity pattern of the system matrix.
+                                     */
     dSMatrixStruct      system_sparsity;
 
+                                    /**
+                                     * System matrix.
+                                     */
+    dSMatrix            system_matrix;
+
+                                    /**
+                                     * Pointers to the right hand sides to the
+                                     * problem. Usually, one only needs one
+                                     * right hand side, but simetimes one
+                                     * wants more.
+                                     */
     vector<dVector*>    right_hand_sides;
+
+                                    /**
+                                     * Solution vector.
+                                     */
     dVector             solution;
 
   friend class Assembler<dim>;
index 115730e9d4cbb0a694c791159074213853902f40..525def890fe81116f7eec56d140f83cb82209ac0 100644 (file)
@@ -235,6 +235,11 @@ void ConstraintMatrix::condense (dSMatrixStruct &sparsity) const {
 
 
 
+void ConstraintMatrix::condense (const dSMatrix &uncondensed,
+                                dSMatrix       &condensed) const {};
+void ConstraintMatrix::condense (dSMatrix &uncondensed) const {};
+
+
 
 unsigned int ConstraintMatrix::n_constraints () const {
   return lines.size();
index ef6bf77ba5e3fd1cfd1e140bd94402e2c149a879..e88aa3bdd7d4772efd81f4b67dccf12a5eb72d69 100644 (file)
@@ -834,6 +834,7 @@ void DoFHandler<2>::do_renumbering (const vector<int> &new_numbers) {
 
 void DoFHandler<1>::make_constraint_matrix (ConstraintMatrix &cm) const {
   cm.clear ();
+  cm.close ();
 };
 
 
index 6d70b92e3b3b3f5c7f8ffd3171f88aa6a7ced401..c2e14343593a5f1721bd39b74eb2eea36c1e1219 100644 (file)
@@ -26,7 +26,8 @@ FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
                JxW_values(quadrature.n_quadrature_points, 0),
                quadrature_points(quadrature.n_quadrature_points, Point<dim>()),
                unit_quadrature_points(quadrature.n_quadrature_points, Point<dim>()),
-               jacobi_matrices (quadrature.n_quadrature_points)
+               jacobi_matrices (quadrature.n_quadrature_points,
+                                dFMatrix(dim,dim))
 {
   for (unsigned int i=0; i<fe.total_dofs; ++i)
     for (unsigned int j=0; j<n_quadrature_points; ++j) 
@@ -95,6 +96,7 @@ void FEValues<1>::reinit (const Triangulation<1>::cell_iterator &cell,
                     unit_quadrature_points,
                     jacobi_matrices,
                     quadrature_points);
+
                                   // compute gradients on real element
   for (unsigned int i=0; i<fe.total_dofs; ++i)
     for (unsigned int j=0; j<n_quadrature_points; ++j)
index f1b7be3b181eeb1a2df4243a049c2be44027afd8..16c685166e419f884a6e9ea2dd4ac0ffca2d8a07 100644 (file)
@@ -51,7 +51,7 @@ FELinear<1>::FELinear () :
 
 double
 FELinear<1>::shape_value(const unsigned int i,
-                              const Point<1>& p) const
+                        const Point<1>& p) const
 {
   Assert((i<total_dofs), ExcInvalidIndex(i));
   switch (i)
@@ -226,8 +226,7 @@ FELinear<2>::shape_grad (const unsigned int i,
 
 
 // this function may be generalised to three or more dimensions with gcc2.8
-// you will have to change the number of vertices and the setting up of the
-// jacobi matrices
+// you will have to change the number of vertices
 void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell,
                                  const vector<Point<2> >               &unit_points,
                                  vector<dFMatrix>  &jacobians,
@@ -332,6 +331,40 @@ FEQuadratic<2>::FEQuadratic () :
 
 
 
+template <int dim>
+double
+FEQuadratic<dim>::shape_value (const unsigned int i,
+                              const Point<dim> &) const
+{
+  Assert (i<total_dofs, ExcInvalidIndex(i));
+  Assert (false, ExcNotImplemented());
+  return 0.;
+};
+
+
+
+template <int dim>
+Point<dim>
+FEQuadratic<dim>::shape_grad (const unsigned int i,
+                             const Point<dim> &) const
+{
+  Assert (i<total_dofs, ExcInvalidIndex(i));
+  Assert (false, ExcNotImplemented());
+  return Point<dim> ();
+};
+
+
+
+void FEQuadratic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
+                                    const vector<Point<2> >               &,
+                                    vector<dFMatrix>  &,
+                                    vector<Point<2> > &) const {
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+
 
 
 FECubic<1>::FECubic () :
@@ -353,3 +386,51 @@ void FECubic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
 FECubic<2>::FECubic () :
                FiniteElement<2> (1, 2, 4) {};
 
+
+
+
+template <int dim>
+double
+FECubic<dim>::shape_value (const unsigned int i,
+                          const Point<dim> &) const
+{
+  Assert (i<total_dofs, ExcInvalidIndex(i));
+  Assert (false, ExcNotImplemented());
+  return 0.;
+};
+
+
+
+template <int dim>
+Point<dim>
+FECubic<dim>::shape_grad (const unsigned int i,
+                         const Point<dim> &) const
+{
+  Assert (i<total_dofs, ExcInvalidIndex(i));
+  Assert (false, ExcNotImplemented());
+  return Point<dim> ();
+};
+
+
+
+void FECubic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
+                                const vector<Point<2> >               &,
+                                vector<dFMatrix>  &,
+                                vector<Point<2> > &) const {
+  Assert (false, ExcNotImplemented());
+};
+
+
+
+
+
+
+// explicite instantiations
+
+template class FELinear<1>;
+template class FELinear<2>;
+template class FEQuadratic<1>;
+template class FEQuadratic<2>;
+template class FECubic<1>;
+template class FECubic<2>;
+
index 1a1da666a930c4edf2ba35674294ce2c5b230a69..3b99a37c698128a2f2391d3a69fa74c8d85226da 100644 (file)
@@ -3,6 +3,7 @@
 #include <numerics/assembler.h>
 #include <numerics/base.h>
 #include <grid/tria_iterator.h>
+#include <grid/tria_iterator.templates.h>
 #include <lac/dfmatrix.h>
 #include <lac/dvector.h>
 
@@ -67,18 +68,20 @@ void Equation<2>::assemble (vector<dVector>   &,
 
 
 
-/*
+
 template <int dim>
-Assembler<dim>::AssemblerData::AssemblerData (DoFHandler<dim> *dof,
-                                             const bool assemble_matrix,
-                                             const bool assemble_rhs,
-                                             const unsigned int n_rhs,
-                                             ProblemBase<dim> *problem,
-                                             const Quadrature<dim> *quadrature) :
+AssemblerData<dim>::AssemblerData (const DoFHandler<dim>   &dof,
+                                  const bool               assemble_matrix,
+                                  const bool               assemble_rhs,
+                                  const unsigned int       n_rhs,
+                                  ProblemBase<dim>         &problem,
+                                  const Quadrature<dim>    &quadrature,
+                                  const FiniteElement<dim> &fe) :
                dof(dof), assemble_matrix(assemble_matrix),
                assemble_rhs(assemble_rhs), n_rhs(n_rhs),
-               problem(problem), quadrature(quadrature) {};
-               */
+               problem(problem), quadrature(quadrature),
+               fe(fe) {};
+
 
 
 
@@ -88,31 +91,31 @@ Assembler<dim>::Assembler (Triangulation<dim> *tria,
                           const int           index,
                           const void         *local_data) :
                DoFCellAccessor<dim> (tria,level,index,
-                                     ((AssemblerData*)local_data)->dof),
+                                     &((AssemblerData<dim>*)local_data)->dof),
                cell_matrix (dof_handler->get_selected_fe().total_dofs),
-               cell_vectors (((AssemblerData*)local_data)->n_rhs,
+               cell_vectors (((AssemblerData<dim>*)local_data)->n_rhs,
                              dVector(dof_handler->get_selected_fe().total_dofs)),
-               assemble_matrix (((AssemblerData*)local_data)->assemble_matrix),
-               assemble_rhs (((AssemblerData*)local_data)->assemble_rhs),
-               problem (((AssemblerData*)local_data)->problem),
-               fe_values (dof_handler->get_selected_fe(),
-                          *((AssemblerData*)local_data)->quadrature)
+               assemble_matrix (((AssemblerData<dim>*)local_data)->assemble_matrix),
+               assemble_rhs (((AssemblerData<dim>*)local_data)->assemble_rhs),
+               problem (((AssemblerData<dim>*)local_data)->problem),
+               fe(((AssemblerData<dim>*)local_data)->fe),
+               fe_values (((AssemblerData<dim>*)local_data)->fe,
+                          ((AssemblerData<dim>*)local_data)->quadrature)
 {
-  Assert (((AssemblerData*)local_data)->dof != 0,
+  Assert (((AssemblerData<dim>*)local_data)->n_rhs != 0,
          ExcInvalidData());
-  Assert (((AssemblerData*)local_data)->n_rhs != 0,
+  Assert ((unsigned int)problem.system_matrix.m() == dof_handler->n_dofs(),
          ExcInvalidData());
-  Assert (((AssemblerData*)local_data)->problem != 0,
+  Assert ((unsigned int)problem.system_matrix.n() == dof_handler->n_dofs(),
          ExcInvalidData());
-
-  Assert ((unsigned int)problem->system_matrix.m() == dof_handler->n_dofs(),
+  Assert (problem.right_hand_sides.size() == cell_vectors.size(),
          ExcInvalidData());
-  Assert ((unsigned int)problem->system_matrix.n() == dof_handler->n_dofs(),
-         ExcInvalidData());
-  Assert (problem->right_hand_sides.size() == cell_vectors.size(),
+  Assert (((AssemblerData<dim>*)local_data)->fe == dof_handler->get_selected_fe(),
          ExcInvalidData());
+
+
   for (unsigned int i=0; i<cell_vectors.size(); ++i)
-    Assert ((unsigned int)problem->right_hand_sides[i]->n() == dof_handler->n_dofs(),
+    Assert ((unsigned int)problem.right_hand_sides[i]->n() == dof_handler->n_dofs(),
            ExcInvalidData());
 };
 
@@ -124,7 +127,7 @@ void Assembler<dim>::assemble (const Equation<dim> &equation) {
   fe_values.reinit (Triangulation<dim>::cell_iterator (tria,
                                                       present_level,
                                                       present_index),
-                   dof_handler->get_selected_fe());
+                   fe);
   
   if (assemble_matrix)
                                     // clear cell matrix
@@ -170,13 +173,13 @@ void Assembler<dim>::assemble (const Equation<dim> &equation) {
                                     // distribute cell matrix
     for (unsigned int i=0; i<dof_handler->get_selected_fe().total_dofs; ++i)
       for (unsigned int j=0; j<dof_handler->get_selected_fe().total_dofs; ++j)
-       problem->system_matrix.set(dofs[i], dofs[j], cell_matrix(i,j));
+       problem.system_matrix.add(dofs[i], dofs[j], cell_matrix(i,j));
 
   if (assemble_rhs)
                                     // distribute cell vector
     for (unsigned int vec=0; vec<cell_vectors.size(); ++vec)
       for (unsigned int j=0; j<dof_handler->get_selected_fe().total_dofs; ++j)
-       (*problem->right_hand_sides[vec])(dofs[j]) = cell_vectors[vec](j);
+       (*problem.right_hand_sides[vec])(dofs[j]) += cell_vectors[vec](j);
 };
 
                
@@ -189,3 +192,15 @@ template class Equation<2>;
 
 template class Assembler<1>;
 template class Assembler<2>;
+
+
+template class AssemblerData<1>;
+template class AssemblerData<2>;
+
+
+template class TriaRawIterator<1,Assembler<1> >;
+template class TriaIterator<1,Assembler<1> >;
+template class TriaActiveIterator<1,Assembler<1> >;
+template class TriaRawIterator<2,Assembler<2> >;
+template class TriaIterator<2,Assembler<2> >;
+template class TriaActiveIterator<2,Assembler<2> >;
index a250604a3aeb2d0a2bb4d61dcfcfb883658f8dc6..1bd73903bfe7579a255d6c7819996622971acc8e 100644 (file)
@@ -18,7 +18,7 @@ ProblemBase<dim>::ProblemBase (Triangulation<dim> *tria,
                tria(tria),
                dof_handler(dof),
                system_sparsity(1,1,1),        // dummy initialisation, is later reinit'd
-               system_matrix(system_sparsity),// dummy initialisation, is later reinit'd
+               system_matrix(),               // dummy initialisation, is later reinit'd
                right_hand_sides (n_rhs, (dVector*)0)   
 {
   Assert (tria == &dof->get_tria(), ExcDofAndTriaDontMatch());
@@ -35,7 +35,8 @@ ProblemBase<dim>::ProblemBase (Triangulation<dim> *tria,
 
 template <int dim>
 void ProblemBase<dim>::assemble (const Equation<dim>   &equation,
-                                const Quadrature<dim> &quadrature) {
+                                const Quadrature<dim> &quadrature,
+                                const FiniteElement<dim> &fe) {
   system_sparsity.reinit (dof_handler->n_dofs(),
                          dof_handler->n_dofs(),
                          dof_handler->max_couplings_between_dofs());
@@ -55,11 +56,12 @@ void ProblemBase<dim>::assemble (const Equation<dim>   &equation,
   system_matrix.reinit (system_sparsity);
 
                                   // create assembler
-  Assembler<dim>::AssemblerData data (dof_handler,
-                                     true, true, //assemble matrix and rhs
-                                     right_hand_sides.size(),   //one rhs
-                                     this,
-                                     &quadrature);
+  AssemblerData<dim> data (*dof_handler,
+                          true, true, //assemble matrix and rhs
+                          right_hand_sides.size(),   //one rhs
+                          *this,
+                          quadrature,
+                          fe);
   TriaActiveIterator<dim, Assembler<dim> > assembler (tria,
                                                      tria->begin_active()->level(),
                                                      tria->begin_active()->index(),
index 719e01c5d08b40bfd3613bf5c251b51a9fc4bf3c..e874cc83e65b4515b61b83ef3cb43343a45789cb 100644 (file)
@@ -5,8 +5,9 @@ INCPATH   = -I../include -I../../lac/include -I../../base/include
 CXXFLAGS.g= -DDEBUG -g -Wall $(INCPATH)
 CXXFLAGS  = -O3 $(INCPATH)
 
-LIBS.g    = -L../lib -lbasic.g -lgrid.g -lfe.g -lnumerics.g -llac.g
-LIBS      = -L../lib -lbasic -lgrid -lfe -lnumerics 
+LIBPATH   = -L../lib -L../../lac/lib
+LIBS.g    = $(LIBPATH) -lbasic.g -lgrid.g -lfe.g -lnumerics.g -llac.g
+LIBS      = $(LIBPATH) -lbasic -lgrid -lfe -lnumerics 
 
 ifeq ($(shell uname),SunOS)
 CXXFLAGS  := -V2.7.2.3 $(CXXFLAGS)
@@ -14,17 +15,22 @@ CXXFLAGS.g:= -V2.7.2.3 $(CXXFLAGS.g)
 endif
 
 
-cc-files = grid/grid_test.cc
+cc-files = grid/grid_test.cc dof/dof_test.cc poisson/poisson.cc
 o-files  = $(cc-files:.cc=.o)
 h-files  = $(wildcard ../include/*.h)
 
-examples: grid/grid_test dof/dof_test
+examples: $(cc-files:.cc=)
 
+%.o : %.cc
+       $(CXX) $(CXXFLAGS) -o $@ -c $<
 
-grid/grid_test: grid/grid_test.o ../lib/libgrid.a
+grid/grid_test: grid/grid_test.o
        $(CXX) $(CXXFLAGS) -o $@ $< $(LIBS)
 
-dof/dof_test: dof/dof_test.o ../lib/libgrid.g.a ../lib/libdeal.g.a
+dof/dof_test: dof/dof_test.o
+       $(CXX) $(CXXFLAGS.g) -o $@ $< $(LIBS.g) -lg++
+
+poisson/poisson: poisson/poisson.o
        $(CXX) $(CXXFLAGS.g) -o $@ $< $(LIBS.g) -lg++
 
 
@@ -52,6 +58,6 @@ clean:
 #to remake this file upon inclusion at the bottom
 #of this file.
 Makefile.dep: $(cc-files) $(h-files)
-       $(CXX) -I../include -M $(cc-files) > Makefile.dep
+       $(CXX) $(INCPATH) -M $(cc-files) > Makefile.dep
 
 include Makefile.dep
index 1ce93aa62f560cd08c03d2be94b7af9f96413ea7..7dafc011520047d240b41794b68e68e6801d6489 100644 (file)
@@ -12,7 +12,7 @@ grid_test.o: grid/grid_test.cc ../include/grid/tria_iterator.h \
  /usr/local/lib/g++-include/new.h /usr/local/lib/g++-include/std/new.h \
  /usr/local/lib/g++-include/std/cstddef.h \
  /usr/local/include/g++/stl_raw_storage_iter.h \
- ../include/basic/exceptions.h ../include/grid/point.h \
+ ../../base/include/base/exceptions.h ../include/grid/point.h \
  ../include/grid/tria_accessor.h ../include/grid/tria.h \
  /usr/local/include/g++/vector.h /usr/local/include/g++/algobase.h \
  /usr/local/include/g++/pair.h /usr/local/include/g++/stl_pair.h \
@@ -31,3 +31,85 @@ grid_test.o: grid/grid_test.cc ../include/grid/tria_iterator.h \
  /usr/include/math.h /usr/include/floatingpoint.h /usr/include/stdio.h \
  /usr/include/sys/va_list.h /usr/include/sys/ieeefp.h \
  /usr/local/include/g++/fstream.h
+dof_test.o: dof/dof_test.cc ../include/grid/dof.h \
+ /usr/local/include/g++/vector.h /usr/local/include/g++/algobase.h \
+ /usr/local/include/g++/pair.h /usr/local/include/g++/stl_config.h \
+ /usr/local/include/g++/stl_relops.h /usr/local/include/g++/stl_pair.h \
+ /usr/local/include/g++/iterator.h /usr/local/include/g++/function.h \
+ /usr/local/lib/gcc-lib/sparc-sun-solaris2.5.1/2.8.0/include/stddef.h \
+ /usr/local/include/g++/stl_function.h \
+ /usr/local/include/g++/iostream.h /usr/local/include/g++/streambuf.h \
+ /usr/local/include/g++/libio.h /usr/local/lib/g++-include/_G_config.h \
+ /usr/local/include/g++/stl_iterator.h \
+ /usr/local/include/g++/type_traits.h \
+ /usr/local/include/g++/stl_construct.h \
+ /usr/local/lib/g++-include/new.h /usr/local/lib/g++-include/std/new.h \
+ /usr/local/lib/g++-include/std/cstddef.h \
+ /usr/local/include/g++/stl_raw_storage_iter.h \
+ /usr/local/include/g++/stl_algobase.h /usr/include/string.h \
+ /usr/include/sys/feature_tests.h \
+ /usr/local/lib/gcc-lib/sparc-sun-solaris2.5.1/2.8.0/include/limits.h \
+ /usr/local/lib/gcc-lib/sparc-sun-solaris2.5.1/2.8.0/include/syslimits.h \
+ /usr/include/limits.h /usr/include/sys/isa_defs.h \
+ /usr/include/stdlib.h /usr/local/include/g++/stl_uninitialized.h \
+ /usr/local/include/g++/alloc.h /usr/local/include/g++/stl_alloc.h \
+ /usr/local/sparc-sun-solaris2.5.1/include/assert.h \
+ /usr/local/include/g++/stl_vector.h ../include/fe/fe.h \
+ ../../lac/include/lac/dfmatrix.h /usr/include/stdio.h \
+ /usr/include/sys/va_list.h ../../lac/include/lac/dvector.h \
+ ../../base/include/base/types.h ../../lac/include/lac/vectorbase.h \
+ ../../base/include/base/exceptions.h ../../lac/include/lac/ivector.h \
+ ../include/grid/point.h ../include/grid/tria.h \
+ ../include/grid/tria_line.h ../include/grid/tria_quad.h \
+ ../include/grid/tria_boundary.h ../include/fe/fe_lib.h \
+ ../include/grid/tria_iterator.h ../include/grid/tria_accessor.h \
+ ../../lac/include/lac/dsmatrix.h \
+ /usr/local/lib/gcc-lib/sparc-sun-solaris2.5.1/2.8.0/include/math.h \
+ /usr/include/math.h /usr/include/floatingpoint.h \
+ /usr/include/sys/ieeefp.h ../include/basic/parameter_handler.h \
+ /usr/local/include/g++/map.h /usr/local/include/g++/tree.h \
+ /usr/local/include/g++/stl_tree.h /usr/local/include/g++/stl_map.h \
+ /usr/local/lib/g++-include/String.h \
+ /usr/local/lib/g++-include/Regex.h ../include/grid/dof_constraints.h \
+ /usr/local/include/g++/fstream.h /usr/local/include/g++/cmath
+poisson.o: poisson/poisson.cc ../include/grid/dof.h \
+ /usr/local/include/g++/vector.h /usr/local/include/g++/algobase.h \
+ /usr/local/include/g++/pair.h /usr/local/include/g++/stl_config.h \
+ /usr/local/include/g++/stl_relops.h /usr/local/include/g++/stl_pair.h \
+ /usr/local/include/g++/iterator.h /usr/local/include/g++/function.h \
+ /usr/local/lib/gcc-lib/sparc-sun-solaris2.5.1/2.8.0/include/stddef.h \
+ /usr/local/include/g++/stl_function.h \
+ /usr/local/include/g++/iostream.h /usr/local/include/g++/streambuf.h \
+ /usr/local/include/g++/libio.h /usr/local/lib/g++-include/_G_config.h \
+ /usr/local/include/g++/stl_iterator.h \
+ /usr/local/include/g++/type_traits.h \
+ /usr/local/include/g++/stl_construct.h \
+ /usr/local/lib/g++-include/new.h /usr/local/lib/g++-include/std/new.h \
+ /usr/local/lib/g++-include/std/cstddef.h \
+ /usr/local/include/g++/stl_raw_storage_iter.h \
+ /usr/local/include/g++/stl_algobase.h /usr/include/string.h \
+ /usr/include/sys/feature_tests.h \
+ /usr/local/lib/gcc-lib/sparc-sun-solaris2.5.1/2.8.0/include/limits.h \
+ /usr/local/lib/gcc-lib/sparc-sun-solaris2.5.1/2.8.0/include/syslimits.h \
+ /usr/include/limits.h /usr/include/sys/isa_defs.h \
+ /usr/include/stdlib.h /usr/local/include/g++/stl_uninitialized.h \
+ /usr/local/include/g++/alloc.h /usr/local/include/g++/stl_alloc.h \
+ /usr/local/sparc-sun-solaris2.5.1/include/assert.h \
+ /usr/local/include/g++/stl_vector.h ../include/fe/fe.h \
+ ../../lac/include/lac/dfmatrix.h /usr/include/stdio.h \
+ /usr/include/sys/va_list.h ../../lac/include/lac/dvector.h \
+ ../../base/include/base/types.h ../../lac/include/lac/vectorbase.h \
+ ../../base/include/base/exceptions.h ../../lac/include/lac/ivector.h \
+ ../include/grid/point.h ../include/grid/tria.h \
+ ../include/grid/tria_line.h ../include/grid/tria_quad.h \
+ ../include/grid/tria_boundary.h ../include/fe/fe_lib.h \
+ ../include/grid/tria_accessor.h ../include/grid/dof_accessor.h \
+ ../include/grid/tria_iterator.h ../../lac/include/lac/dsmatrix.h \
+ /usr/local/lib/gcc-lib/sparc-sun-solaris2.5.1/2.8.0/include/math.h \
+ /usr/include/math.h /usr/include/floatingpoint.h \
+ /usr/include/sys/ieeefp.h ../include/grid/dof_constraints.h \
+ ../include/basic/data_io.h /usr/local/lib/g++-include/String.h \
+ /usr/local/lib/g++-include/Regex.h ../include/numerics/base.h \
+ ../include/numerics/assembler.h ../include/fe/quadrature_lib.h \
+ ../include/fe/quadrature.h /usr/local/include/g++/fstream.h \
+ /usr/local/include/g++/cmath
index f4f0530b0581d7492e0837f8fd2955f516c7a3ad..d6d6321fd973bf188e3adb634a5acf3ea71ea1af 100644 (file)
@@ -4,13 +4,15 @@
 #include <grid/dof.h>
 #include <grid/tria.h>
 #include <fe/fe_lib.h>
-#include <grid/tria_iterator.h>
 #include <grid/tria_accessor.h>
+#include <grid/dof_accessor.h>
+#include <grid/tria_iterator.h>
 #include <lac/dsmatrix.h>
 #include <grid/dof_constraints.h>
 #include <basic/data_io.h>
 #include <numerics/base.h>
 #include <numerics/assembler.h>
+#include <fe/quadrature_lib.h>
 
 #include <fstream.h>
 #include <cmath>
@@ -20,6 +22,7 @@ extern "C" {
 
 extern TriaActiveIterator<1,CellAccessor<1> > x;
 extern TriaActiveIterator<2,CellAccessor<2> > y;
+extern TriaRawIterator<1,DoFLineAccessor<1,LineAccessor<1> > > z;
 
 
 
@@ -33,6 +36,16 @@ class PoissonEquation :  public Equation<dim> {
                           vector<dVector>     &rhs,
                           const FEValues<dim> &fe_values,
                           const Triangulation<dim>::cell_iterator &cell) const;
+    double right_hand_side (const Point<dim> &) const;
+};
+
+
+
+template <int dim>
+inline
+double PoissonEquation<dim>::right_hand_side (const Point<dim> &) const {
+                                  // this yields as solution u=x^2+y^2+...=\vec x^2
+  return -2*dim;
 };
 
 
@@ -45,26 +58,49 @@ void PoissonEquation<1>::assemble (dFMatrix            &cell_matrix,
     for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
       {
        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.JxW(point);
-       rhs[0] = fe_values.shape_value(i,point) *
-                right_hand_side(fe_values.quadrature_point(point)) *
-                fe_values.JxW(point);
+         cell_matrix(i,j) += fe_values.shape_grad(i,point) *
+                             fe_values.shape_grad(j,point) *
+                             fe_values.JxW(point);
+       rhs[0](i) += fe_values.shape_value(i,point) *
+                    right_hand_side(fe_values.quadrature_point(point)) *
+                    fe_values.JxW(point);
+      };
+};
+
+
+
+void PoissonEquation<2>::assemble (dFMatrix            &cell_matrix,
+                                  vector<dVector>     &rhs,
+                                  const FEValues<2>   &fe_values,
+                                  const Triangulation<2>::cell_iterator &) const {
+  for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+    for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+      {
+       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.JxW(point);
+       rhs[0](i) += fe_values.shape_value(i,point) *
+                    right_hand_side(fe_values.quadrature_point(point)) *
+                    fe_values.JxW(point);
       };
 };
 
   
 
 
-void main () {
-  Triangulation<1> tria;
-  DoFHandler<1>    dof(tria);
-  FELinear<1>      fe;
-  ProblemBase<1>   problem(&tria, &dof);
-  Quadrature<1>    q;
+int main () {
+  Triangulation<2>   tria;
+  DoFHandler<2>      dof(&tria);
+  FELinear<2>        fe;
+  ProblemBase<2>     problem(&tria, &dof);
+  PoissonEquation<2> equation;
+  QGauss4<2>         quadrature;
   
-  tria.create_hypercube();
+  tria.create_unit_hypercube();
+  tria.refine_global (1);
   dof.distribute_dofs (fe);
-  problem.assemble (q);
+  problem.assemble (equation, quadrature, fe);
+  
+  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.