]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Merge with mapping-branch.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 6 Mar 2001 10:30:36 +0000 (10:30 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 6 Mar 2001 10:30:36 +0000 (10:30 +0000)
git-svn-id: https://svn.dealii.org/trunk@4115 0785d39b-7218-0410-832d-ea1e28bc413d

20 files changed:
deal.II/deal.II/include/grid/grid_out.h
deal.II/deal.II/include/numerics/error_estimator.h
deal.II/deal.II/include/numerics/matrices.h
deal.II/deal.II/source/grid/grid_out.all_dimensions.cc
deal.II/deal.II/source/grid/grid_out.cc
deal.II/deal.II/source/numerics/assembler.cc
deal.II/deal.II/source/numerics/base.cc
deal.II/deal.II/source/numerics/data_out.cc
deal.II/deal.II/source/numerics/data_out_faces.cc
deal.II/deal.II/source/numerics/data_out_rotation.cc
deal.II/deal.II/source/numerics/data_out_stack.cc
deal.II/deal.II/source/numerics/derivative_approximation.cc
deal.II/deal.II/source/numerics/error_estimator.cc
deal.II/deal.II/source/numerics/matrices.cc
deal.II/deal.II/source/numerics/vectors.cc
deal.II/examples/step-2/step-2.cc
deal.II/examples/step-3/step-3.cc
deal.II/examples/step-4/step-4.cc
deal.II/examples/step-5/step-5.cc
deal.II/examples/step-6/step-6.cc

index b9eb33ecb9336f42e079bde74b22850cd1f3f6c0..f42be24657dc6de96bc5b7b0f7275cd4c9a9ea43 100644 (file)
@@ -18,6 +18,7 @@
 #include <string>
 
 template <int dim> class Triangulation;
+template <int dim> class Mapping;
 
 /**
  * This class provides a means to output a triangulation to a file in different
@@ -172,10 +173,21 @@ class GridOut
                                          * Default: @p{false}.
                                          */
        bool write_cell_numbers;
+
+                                        /**
+                                         * This is the number of
+                                         * points on a boundary face,
+                                         * that are ploted
+                                         * additionally to the
+                                         * vertices of the face.
+                                         */
+       unsigned int n_boundary_face_points;
+       
                                         /**
                                          * Constructor.
                                          */
-       GnuplotFlags (const bool write_cell_number = false);
+       GnuplotFlags (const bool         write_cell_number = false,
+                     const unsigned int n_boundary_face_points = 2);
     };
 
                                     /**
@@ -232,6 +244,20 @@ class GridOut
                                          * be drawn in a different color?
                                          */
        bool color_lines_on_user_flag;
+
+                                        /**
+                                         * This is the number of
+                                         * points on a boundary face,
+                                         * that are ploted
+                                         * additionally to the
+                                         * vertices of the face.
+                                         *
+                                         * This is used if the
+                                         * mapping used is not the
+                                         * standard @p{MappingQ1}
+                                         * mapping.
+                                         */
+       unsigned int n_boundary_face_points;
        
                                         /**
                                          * Constructor.
@@ -239,7 +265,8 @@ class GridOut
        EpsFlagsBase (const SizeType     size_type  = width,
                      const unsigned int size       = 300,
                      const double       line_width = 0.5,
-                     bool color_lines_on_user_flag = false);
+                     const bool color_lines_on_user_flag = false,
+                     const unsigned int n_boundary_face_points = 2);
     };
 
 
@@ -252,13 +279,49 @@ class GridOut
                                      * listed in the base class
                                      */
     template <int dim>
-    struct EpsFlags : public EpsFlagsBase 
+    struct EpsFlags;
+//        : public EpsFlagsBase
+//      {
+       
+//                                      /**
+//                                       * Constructor.
+//                                       */
+//     EpsFlags (const SizeType     size_type  = width,
+//               const unsigned int size       = 300,
+//               const double       line_width = 0.5,
+//               bool color_lines_on_user_flag = false,
+//               const unsigned int n_boundary_face_points = 2);
+//      };
+
+//  //#if (__GNUC__==2) && (__GNUC_MINOR__ < 95)
+//  //    template <>
+    
+    struct EpsFlags<1> : public EpsFlagsBase 
     {
+       
+                                        /**
+                                         * Constructor.
+                                         */
+       EpsFlags (const SizeType     size_type  = width,
+                 const unsigned int size       = 300,
+                 const double       line_width = 0.5,
+                 bool color_lines_on_user_flag = false,
+                 const unsigned int n_boundary_face_points = 2);
     };
 
-//#if (__GNUC__==2) && (__GNUC_MINOR__ < 95)
-//    template <>
-
+    
+    struct EpsFlags<2> : public EpsFlagsBase 
+    {
+       
+                                        /**
+                                         * Constructor.
+                                         */
+       EpsFlags (const SizeType     size_type  = width,
+                 const unsigned int size       = 300,
+                 const double       line_width = 0.5,
+                 bool color_lines_on_user_flag = false,
+                 const unsigned int n_boundary_face_points = 2);
+    };
                                     /**
                                      * Flags specific to the output of
                                      * grids in three space dimensions.
@@ -292,11 +355,16 @@ class GridOut
                                         /**
                                          * Constructor.
                                          */
-       EpsFlags (const double        azimut_angle    = 60,
-                 const double        turn_angle      = 30) :
-                       azimut_angle (azimut_angle),
-                       turn_angle (turn_angle)      {};
+       EpsFlags (const SizeType     size_type  = width,
+                 const unsigned int size       = 300,
+                 const double       line_width = 0.5,
+                 const bool color_lines_on_user_flag = false,
+                 const unsigned int n_boundary_face_points = 2,
+                 const double        azimut_angle    = 60,
+                 const double        turn_angle      = 30);
     };
+
+
 //#else
 //#   warning Not implemented for gcc2.95
 //#endif
@@ -355,6 +423,12 @@ class GridOut
                                      * plots of the different levels
                                      * of grid.
                                      *
+                                     * @p{mapping} is a pointer to a
+                                     * mapping used for the
+                                     * transformation of cells at the
+                                     * boundary. If zero, then use
+                                     * standard Q1 mapping.
+                                     *
                                      * Names and values of additional
                                      * flags controlling the output
                                      * can be found in the
@@ -363,7 +437,8 @@ class GridOut
                                      */
     template <int dim>
     void write_gnuplot (const Triangulation<dim> &tria,
-                       std::ostream             &out);
+                       std::ostream             &out,
+                       const Mapping<dim> *mapping=0);
 
                                     /**
                                      * Write the triangulation in the
@@ -444,6 +519,12 @@ class GridOut
                                      * file and can be changed there
                                      * according to need.
                                      *
+                                     * @p{mapping} is a pointer to a
+                                     * mapping used for the
+                                     * transformation of cells at the
+                                     * boundary. If zero, then use
+                                     * standard Q1 mapping.
+                                     *
                                      * Names and values of additional
                                      * flags controlling the output
                                      * can be found in the
@@ -456,7 +537,8 @@ class GridOut
                                      */
     template <int dim>
     void write_eps (const Triangulation<dim> &tria,
-                   std::ostream             &out);
+                   std::ostream             &out,
+                   const Mapping<dim>       *mapping=0);
     
                                     /**
                                      * Write data and grid to @p{out} according
@@ -467,7 +549,8 @@ class GridOut
     template <int dim>
     void write (const Triangulation<dim> &tria,
                std::ostream             &out,
-               const OutputFormat        output_format);
+               const OutputFormat        output_format,
+               const Mapping<dim>       *mapping=0);
 
                                     /**
                                      * Set the flags to be used for output
index 4f28bcdb4d93139b66442a5b4be03d5bd4eeb00c..96f53add2fc76574fc02b5575a229dbe1430622b 100644 (file)
@@ -1,5 +1,9 @@
 //----------------------------  error_estimator.h  ---------------------------
+//    $Id$
 //    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+//
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
 //    to the file deal.II/doc/license.html for the  text  and
index c81d70f4bdc8c68ec495e444599ed2dde77ea683..32e186486daee71caeabe925a4f774f15ef8ce9c 100644 (file)
@@ -197,15 +197,7 @@ class MatrixCreator
                                     /**
                                      * Assemble the mass matrix. If no 
                                      * coefficient is given, it is assumed
-                                     * to be constant one.
-                                     * 
-                                     * If the coefficient is constant, it
-                                     * may be more adequate to use the
-                                     * functions assembling the mass matrix
-                                     * without quadrature. However, the
-                                     * two functions have different effects
-                                     * for finite elements composed of
-                                     * several subobjects.
+                                     * to be unity.
                                      * 
                                      * See the general doc of this class
                                      * for more information.
@@ -219,16 +211,8 @@ class MatrixCreator
                                      * Assemble the mass matrix and a right
                                      * hand side vector. If no 
                                      * coefficient is given, it is assumed
-                                     * to be constant one.
+                                     * to be unity.
                                      *
-                                     * If the coefficient is constant, it
-                                     * may be more adequate to use the
-                                     * functions assembling the mass matrix
-                                     * without quadrature. However, the
-                                     * two functions have different effects
-                                     * for finite elements composed of
-                                     * several subobjects.
-                                     * 
                                      * See the general doc of this class
                                      * for more information.
                                      */
@@ -239,28 +223,6 @@ class MatrixCreator
                                    Vector<double>           &rhs_vector,
                                    const Function<dim>      *a = 0);
 
-                                    /**
-                                     * Create the mass matrix by exact
-                                     * evaluation without using a quadrature
-                                     * formula.
-                                     *
-                                     * No right hand side may be created using
-                                     * this function. See the general doc of
-                                     * this class for more information.
-                                     *
-                                     * It is assumed that the matrix already
-                                     * has the right size. The mass matrix
-                                     * elements are summed up to the values
-                                     * previously in the matrix, so if you want
-                                     * the pure mass matrix, you have to clear
-                                     * the matrix beforehand.
-                                     * 
-                                     * See the general doc of this class
-                                     * for more information.
-                                     */
-    static void create_mass_matrix (const DoFHandler<dim>    &dof,
-                                   SparseMatrix<double>     &matrix);
-    
                                     /**
                                      * Assemble the mass matrix and a right
                                      * hand side vector along the boundary.
index fb3a56742279f53a122b238cb71c57a5b6945df7..812385cf807411699e66f2b8e882c7b11ab9932a 100644 (file)
@@ -22,31 +22,62 @@ GridOut::UcdFlags::UcdFlags (const bool write_preamble,
 {}
 
 
-GridOut::GnuplotFlags::GnuplotFlags (const bool write_cell_numbers) :
-               write_cell_numbers (write_cell_numbers)
+GridOut::GnuplotFlags::GnuplotFlags (const bool write_cell_numbers,
+                                    const unsigned int n_boundary_face_points) :
+               write_cell_numbers (write_cell_numbers),
+               n_boundary_face_points(n_boundary_face_points)
 {}
 
 
 GridOut::EpsFlagsBase::EpsFlagsBase (const SizeType     size_type,
                                     const unsigned int size,
                                     const double       line_width,
-                                    bool color_lines_on_user_flag) :
+                                    const bool color_lines_on_user_flag,
+                                    const unsigned int n_boundary_face_points) :
                size_type (size_type),
                size (size),
                line_width (line_width),
-               color_lines_on_user_flag(color_lines_on_user_flag)
+               color_lines_on_user_flag(color_lines_on_user_flag),
+               n_boundary_face_points(n_boundary_face_points)
 {}
 
 
-// egcs 1.1.2 does not understand this, so I made it inlined in the
-// class declaration
-//
-// template <>
-// GridOut::EpsFlags<3>::EpsFlags (const double        azimut_angle,
-//                             const double        turn_angle) :
-//             azimut_angle (azimut_angle),
-//             turn_angle (turn_angle)
-// {};
+
+GridOut::EpsFlags<1>::EpsFlags (const SizeType     size_type,
+                               const unsigned int size,
+                               const double       line_width,
+                               const bool color_lines_on_user_flag,
+                               const unsigned int n_boundary_face_points):
+                       EpsFlagsBase(size_type, size, line_width,
+                                    color_lines_on_user_flag,
+                                    n_boundary_face_points)
+{}
+
+
+GridOut::EpsFlags<2>::EpsFlags (const SizeType     size_type,
+                               const unsigned int size,
+                               const double       line_width,
+                               const bool color_lines_on_user_flag,
+                               const unsigned int n_boundary_face_points):
+                       EpsFlagsBase(size_type, size, line_width,
+                                    color_lines_on_user_flag,
+                                    n_boundary_face_points)
+{}
+
+
+GridOut::EpsFlags<3>::EpsFlags (const SizeType     size_type,
+                               const unsigned int size,
+                               const double       line_width,
+                               const bool color_lines_on_user_flag,
+                               const unsigned int n_boundary_face_points,
+                               const double        azimut_angle,
+                               const double        turn_angle):
+                       EpsFlagsBase(size_type, size, line_width,
+                                    color_lines_on_user_flag,
+                                    n_boundary_face_points),
+                       azimut_angle (azimut_angle),
+                       turn_angle (turn_angle)
+{}
 
 
 void GridOut::set_flags (const UcdFlags &flags) 
index 1eda7a75a6717085ecf083e364cda3090daee15f..0720790031426e4c24775d4b025eba980f3c8f0d 100644 (file)
 
 
 #include <base/point.h>
+#include <base/quadrature.h>
 #include <grid/grid_out.h>
 #include <grid/tria.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
+#include <fe/mapping.h>
 
 #include <iomanip>
 #include <algorithm>
@@ -222,16 +224,68 @@ void GridOut::write_ucd_faces (const Triangulation<dim> &tria,
       };         
 };
 
+#if deal_II_dimension==1
+
+template <>
+void GridOut::write_gnuplot (const Triangulation<1> &tria,
+                            std::ostream           &out,
+                            const Mapping<1>       *) 
+{
+  AssertThrow (out, ExcIO());
+
+  Triangulation<1>::active_cell_iterator        cell=tria.begin_active();
+  const Triangulation<1>::active_cell_iterator  endc=tria.end();
+  for (; cell!=endc; ++cell)
+    {
+      if (gnuplot_flags.write_cell_numbers)
+       out << "# cell " << cell << std::endl;
+//TODO: plot level and material according to switches      
+
+      out << cell->vertex(0)
+         << ' ' << cell->level()
+         << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+         << cell->vertex(1)
+                   << ' ' << cell->level()
+         << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+         << std::endl;
+      break;
+    }
+  AssertThrow (out, ExcIO());
+}
+
+
+#else
 
 
 template <int dim>
 void GridOut::write_gnuplot (const Triangulation<dim> &tria,
-                            std::ostream             &out) 
+                            std::ostream             &out,
+                            const Mapping<dim>       *mapping) 
 {
   AssertThrow (out, ExcIO());
 
+  const unsigned int n_points=
+    gnuplot_flags.n_boundary_face_points;
+
   typename Triangulation<dim>::active_cell_iterator        cell=tria.begin_active();
   const typename Triangulation<dim>::active_cell_iterator  endc=tria.end();
+
+                                  // if we are to treat curved
+                                  // boundaries, then generate a
+                                  // quadrature formula which will be
+                                  // used to probe boundary points at
+                                  // curved faces
+  QProjector<dim> *q_projector=0;
+  if (mapping!=0)
+    {
+      typename std::vector<Point<dim-1> > boundary_points(n_points);
+      for (unsigned int i=0; i<n_points; ++i)
+       boundary_points[i](0)= 1.*(i+1)/(n_points+1);
+
+      Quadrature<dim-1> quadrature(boundary_points, std::vector<double> (n_points, 1));
+      q_projector = new QProjector<dim> (quadrature, false);
+    }
+  
   for (; cell!=endc; ++cell)
     {
       if (gnuplot_flags.write_cell_numbers)
@@ -240,16 +294,24 @@ void GridOut::write_gnuplot (const Triangulation<dim> &tria,
       switch (dim)
        {
          case 1:
-               out << cell->vertex(0)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << cell->vertex(1)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << std::endl;
-               break;
-
+         {
+           Assert(false, ExcInternalError());
+           break;
+         };
+          
          case 2:
+         {
+           if (mapping==0 || !cell->at_boundary())
+             {
+                                                // write out the four
+                                                // sides of this cell
+                                                // by putting the
+                                                // four points (+ the
+                                                // initial point
+                                                // again) in a row
+                                                // and lifting the
+                                                // drawing pencil at
+                                                // the end
                out << cell->vertex(0)
                    << ' ' << cell->level()
                    << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
@@ -267,81 +329,153 @@ void GridOut::write_gnuplot (const Triangulation<dim> &tria,
                    << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
                    << std::endl  // double new line for gnuplot 3d plots
                    << std::endl;
-               break;
-
+             }
+           else
+                                              // cell is at boundary
+                                              // and we are to treat
+                                              // curved
+                                              // boundaries. so loop
+                                              // over all faces and
+                                              // draw them as small
+                                              // pieces of lines
+             {
+               for (unsigned int face_no=0;
+                    face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+                 {
+                   const typename Triangulation<dim>::face_iterator
+                     face = cell->face(face_no);
+                   if (face->at_boundary())
+                     {
+                       out << face->vertex(0)
+                           << ' ' << cell->level()
+                           << ' ' << static_cast<unsigned int>(cell->material_id())
+                           << std::endl;
+                       
+                       const unsigned int offset=face_no*n_points;
+                       for (unsigned int i=0; i<n_points; ++i)
+                         out << (mapping->transform_unit_to_real_cell
+                                 (cell, q_projector->point(offset+i)))
+                             << ' ' << cell->level()
+                             << ' ' << static_cast<unsigned int>(cell->material_id())
+                             << std::endl;
+                       
+                       out << face->vertex(1)
+                           << ' ' << cell->level()
+                           << ' ' << static_cast<unsigned int>(cell->material_id())
+                           << std::endl
+                           << std::endl
+                           << std::endl;
+                     }
+                   else
+                     {
+                                                        // if,
+                                                        // however,
+                                                        // the face
+                                                        // is not at
+                                                        // the
+                                                        // boundary,
+                                                        // then draw
+                                                        // it as
+                                                        // usual
+                       out << face->vertex(0)
+                           << ' ' << cell->level()
+                           << ' ' << static_cast<unsigned int>(cell->material_id())
+                           << std::endl
+                           << face->vertex(1)
+                           << ' ' << cell->level()
+                           << ' ' << static_cast<unsigned int>(cell->material_id())
+                           << std::endl
+                           << std::endl
+                           << std::endl;
+                     };
+                 };
+             };
+           
+           break;
+         };
+          
          case 3:
-                                                // front face
-               out << cell->vertex(0)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << cell->vertex(1)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << cell->vertex(2)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << cell->vertex(3)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << cell->vertex(0)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << std::endl;
-                                                // back face
-               out << cell->vertex(4)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << cell->vertex(5)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << cell->vertex(6)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << cell->vertex(7)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << cell->vertex(4)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << std::endl;
-
-                                                // now for the four connecting lines
-               out << cell->vertex(0)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << cell->vertex(4)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << std::endl;
-               out << cell->vertex(1)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << cell->vertex(5)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << std::endl;
-               out << cell->vertex(2)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << cell->vertex(6)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << std::endl;
-               out << cell->vertex(3)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << cell->vertex(7)
-                   << ' ' << cell->level()
-                   << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
-                   << std::endl;
-               break;
+         {
+//TODO: curved boundaries in 3d gnuplot not supported      
+           Assert (mapping == 0, ExcNotImplemented());
+           
+                                            // front face
+           out << cell->vertex(0)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << cell->vertex(1)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << cell->vertex(2)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << cell->vertex(3)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << cell->vertex(0)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << std::endl;
+                                            // back face
+           out << cell->vertex(4)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << cell->vertex(5)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << cell->vertex(6)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << cell->vertex(7)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << cell->vertex(4)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << std::endl;
+           
+                                            // now for the four connecting lines
+           out << cell->vertex(0)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << cell->vertex(4)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << std::endl;
+           out << cell->vertex(1)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << cell->vertex(5)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << std::endl;
+           out << cell->vertex(2)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << cell->vertex(6)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << std::endl;
+           out << cell->vertex(3)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << cell->vertex(7)
+               << ' ' << cell->level()
+               << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+               << std::endl;
+           break;
+         };
        };
     };
+
+  if (q_projector != 0)
+    delete q_projector;
+  
   
   AssertThrow (out, ExcIO());
 };
 
-
+#endif
 
 struct LineEntry
 {
@@ -355,10 +489,23 @@ struct LineEntry
 };
 
 
+#if deal_II_dimension==1
+
+template <>
+void GridOut::write_eps (const Triangulation<1> &,
+                        std::ostream &,
+                        const Mapping<1> *) 
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+#else
 
 template <int dim>
 void GridOut::write_eps (const Triangulation<dim> &tria,
-                        std::ostream             &out) 
+                        std::ostream             &out,
+                        const Mapping<dim>       *mapping) 
 {
   
   typedef std::list<LineEntry> LineList;
@@ -368,15 +515,14 @@ void GridOut::write_eps (const Triangulation<dim> &tria,
                                   // in order to avoid the recurring
                                   // distinctions between
                                   // eps_flags_1, eps_flags_2, ...
-  const EpsFlagsBase &eps_flags_base = (dim==1 ?
-                                       (EpsFlagsBase&)eps_flags_1 :
-                                       (dim==2 ?
-                                        (EpsFlagsBase&)eps_flags_2 :
-                                        (dim==3 ?
-                                         (EpsFlagsBase&)eps_flags_3 :
-                                         *(EpsFlagsBase*)0)));
+  const EpsFlagsBase &eps_flags_base = (dim==2 ?
+                                       static_cast<EpsFlagsBase&>(eps_flags_2) :
+                                       (dim==3 ?
+                                        static_cast<EpsFlagsBase&>(eps_flags_3) :
+                                        *static_cast<EpsFlagsBase*>(0)));
   
   AssertThrow (out, ExcIO());
+  const unsigned int n_points = eps_flags_base.n_boundary_face_points;
 
                                   // make up a list of lines by which
                                   // we will construct the triangulation
@@ -393,36 +539,120 @@ void GridOut::write_eps (const Triangulation<dim> &tria,
 
   switch (dim)
     {
+      case 1:
+      {
+       Assert(false, ExcInternalError());
+       break;
+      };
+       
       case 2:
       {
        Triangulation<dim>::active_line_iterator line   =tria.begin_active_line ();
        Triangulation<dim>::active_line_iterator endline=tria.end_line ();
-       
+
+                                        // first treat all interior
+                                        // lines and make up a list
+                                        // of them. if curved lines
+                                        // shall not be supported
+                                        // (i.e. no mapping is
+                                        // provided), then also treat
+                                        // all other lines
        for (; line!=endline; ++line)
-                                          // one would expect
-                                          // make_pair(line->vertex(0),
-                                          //           line->vertex(1))
-                                          // here, but that is not
-                                          // dimension independent, since
-                                          // vertex(i) is Point<dim>,
-                                          // but we want a Point<2>.
-                                          // in fact, whenever we're here,
-                                          // the vertex is a Point<dim>,
-                                          // but the compiler does not
-                                          // know this. hopefully, the
-                                          // compiler will optimize away
-                                          // this little kludge
-         line_list.push_back (LineEntry(Point<2>(line->vertex(0)(0),
-                                                 line->vertex(0)(1)),
-                                        Point<2>(line->vertex(1)(0),
-                                                 line->vertex(1)(1)),
-                                        line->user_flag_set()));
+         if (mapping==0 || !line->at_boundary())
+                                            // one would expect
+                                            // make_pair(line->vertex(0),
+                                            //           line->vertex(1))
+                                            // here, but that is not
+                                            // dimension independent, since
+                                            // vertex(i) is Point<dim>,
+                                            // but we want a Point<2>.
+                                            // in fact, whenever we're here,
+                                            // the vertex is a Point<dim>,
+                                            // but the compiler does not
+                                            // know this. hopefully, the
+                                            // compiler will optimize away
+                                            // this little kludge
+           line_list.push_back (LineEntry(Point<2>(line->vertex(0)(0),
+                                                   line->vertex(0)(1)),
+                                          Point<2>(line->vertex(1)(0),
+                                                   line->vertex(1)(1)),
+                                          line->user_flag_set()));
+       
+                                        // next if we are to treat
+                                        // curved boundaries
+                                        // specially, then add lines
+                                        // to the list consisting of
+                                        // pieces of the boundary
+                                        // lines
+       if (mapping!=0)
+         {
+                                            // to do so, first
+                                            // generate a sequence of
+                                            // points on a face and
+                                            // project them onto the
+                                            // faces of a unit cell
+           typename std::vector<Point<dim-1> > boundary_points (n_points);
+           
+           for (unsigned int i=0; i<n_points; ++i)
+             boundary_points[i](0) = 1.*(i+1)/(n_points+1);
+           
+           Quadrature<dim-1> quadrature (boundary_points,
+                                         std::vector<double> (n_points, 1));
+           QProjector<dim> q_projector (quadrature, false);
+
+                                            // next loop over all
+                                            // boundary faces and
+                                            // generate the info from
+                                            // them
+           typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active ();
+           const typename Triangulation<dim>::active_cell_iterator end=tria.end ();
+           for (; cell!=end; ++cell)
+             for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+               {
+                 const typename Triangulation<dim>::face_iterator
+                   face = cell->face(face_no);
+                 
+                 if (face->at_boundary())
+                   {
+                     Point<dim> p0_dim(face->vertex(0));
+                     Point<2>   p0    (p0_dim(0), p0_dim(1));
+
+                                                      // loop over
+                                                      // all pieces
+                                                      // of the line
+                                                      // and generate
+                                                      // line-lets
+                     const unsigned int offset=face_no*n_points;
+                     for (unsigned int i=0; i<n_points; ++i)
+                       {
+                         const Point<dim> p1_dim (mapping->transform_unit_to_real_cell
+                                                  (cell, q_projector.point(offset+i)));
+                         const Point<2>   p1     (p1_dim(0), p1_dim(1));
+                         
+                         line_list.push_back (LineEntry(p0, p1,
+                                                        face->user_flag_set()));
+                         p0=p1;
+                       }
+
+                                                      // generate last piece
+                     const Point<dim> p1_dim (face->vertex(1));
+                     const Point<2>   p1     (p1_dim(0), p1_dim(1));
+                     line_list.push_back (LineEntry(p0, p1,
+                                                    face->user_flag_set()));
+                   };
+               };
+         };
        
        break;
       };
        
       case 3:
       {
+                                        // curved boundary output
+                                        // presently not supported
+//TODO: curved boundaries in eps for 3d        
+       Assert (mapping == 0, ExcNotImplemented());
+       
        Triangulation<dim>::active_line_iterator line   =tria.begin_active_line ();
        Triangulation<dim>::active_line_iterator endline=tria.end_line ();
        
@@ -442,14 +672,8 @@ void GridOut::write_eps (const Triangulation<dim> &tria,
                                         // is in direction of the viewer, but
                                         // I am too tired at present to fix
                                         // this
-//#if !((__GNUC__==2) && (__GNUC_MINOR__==95))
        const double z_angle    = eps_flags_3.azimut_angle;
        const double turn_angle = eps_flags_3.turn_angle;
-//#else
-//#   warning Not implemented for gcc2.95
-//     const double z_angle    = 60;
-//     const double turn_angle = 30;
-//#endif
        const double pi = 3.1415926536;
        const Point<dim> view_direction(-sin(z_angle * 2.*pi / 360.) * sin(turn_angle * 2.*pi / 360.),
                                        +sin(z_angle * 2.*pi / 360.) * cos(turn_angle * 2.*pi / 360.),
@@ -488,6 +712,7 @@ void GridOut::write_eps (const Triangulation<dim> &tria,
     };
 
 
+
                                   // find out minimum and maximum x and
                                   // y coordinates to compute offsets
                                   // and scaling factors
@@ -545,9 +770,9 @@ void GridOut::write_eps (const Triangulation<dim> &tria,
                                         // lower left corner
          << "0 0 "
                                         // upper right corner
-         << static_cast<unsigned int>( (x_max-x_min) * scale )
+         << static_cast<unsigned int>( (x_max-x_min) * scale )+1
          << ' '
-         << static_cast<unsigned int>( (y_max-y_min) * scale )
+         << static_cast<unsigned int>( (y_max-y_min) * scale )+1
          << std::endl;
 
                                       // define some abbreviations to keep
@@ -582,12 +807,14 @@ void GridOut::write_eps (const Triangulation<dim> &tria,
   AssertThrow (out, ExcIO());
 };
 
+#endif
 
 
 template <int dim>
 void GridOut::write (const Triangulation<dim> &tria,
                     std::ostream             &out,
-                    OutputFormat              output_format)
+                    OutputFormat              output_format,
+                    const Mapping<dim>       *mapping)
 {
   switch (output_format)
     {
@@ -596,11 +823,11 @@ void GridOut::write (const Triangulation<dim> &tria,
            return;
 
       case gnuplot:
-           write_gnuplot (tria, out);
+           write_gnuplot (tria, out, mapping);
            return;
 
       case eps:
-           write_eps (tria, out);
+           write_eps (tria, out, mapping);
            return;
     };
 
@@ -612,10 +839,14 @@ void GridOut::write (const Triangulation<dim> &tria,
 // explicit instantiations
 template void GridOut::write_ucd (const Triangulation<deal_II_dimension> &,
                                  std::ostream &);
+#if deal_II_dimension != 1
 template void GridOut::write_gnuplot (const Triangulation<deal_II_dimension> &,
-                                     std::ostream &);
+                                     std::ostream &,
+                                     const Mapping<deal_II_dimension> *);
 template void GridOut::write_eps (const Triangulation<deal_II_dimension> &,
-                                 std::ostream &);
+                                 std::ostream &,
+                                 const Mapping<deal_II_dimension> *);
+#endif
 template void GridOut::write (const Triangulation<deal_II_dimension> &,
                              std::ostream &,
                              OutputFormat);
index a1001adde3a828a96a74e2a0ce7b48fdec4e9460..794690cbc5eacb04fc62b1be43451a0fae304cb0 100644 (file)
@@ -19,6 +19,9 @@
 #include <lac/full_matrix.h>
 #include <lac/vector.h>
 #include <lac/sparse_matrix.h>
+#include <base/quadrature.h>
+#include <fe/mapping_q1.h>
+
 
 // if necessary try to work around a bug in the IBM xlC compiler
 #ifdef XLC_WORK_AROUND_STD_BUG
@@ -26,6 +29,9 @@ using namespace std;
 #endif
 
 
+
+static MappingQ1<deal_II_dimension> mapping;
+
 template <int dim>
 Assembler<dim>::AssemblerData::AssemblerData (const DoFHandler<dim>    &dof,
                                              const bool                assemble_matrix,
@@ -56,7 +62,7 @@ Assembler<dim>::Assembler (Triangulation<dim>  *tria,
                assemble_rhs (local_data->assemble_rhs),
                matrix(local_data->matrix),
                rhs_vector(local_data->rhs_vector),
-               fe_values (dof_handler->get_fe(),
+               fe_values (mapping, dof_handler->get_fe(),
                           local_data->quadrature,
                           local_data->update_flags)
 {
index 0c6e160f7da285b93f69579f09e8a1beae572fd0..a94cdf9920c0fd93888fb83b7f4f70e6cac7e180 100644 (file)
 #include <algorithm>
 #include <cmath>
 
-// if necessary try to work around a bug in the IBM xlC compiler
-#ifdef XLC_WORK_AROUND_STD_BUG
-using namespace std;
-#endif
-
 
 template <int dim>
 ProblemBase<dim>::ProblemBase () :
@@ -170,7 +165,8 @@ void ProblemBase<dim>::solve ()
 
 
 template <int dim>
-void ProblemBase<dim>::fill_data (DataOut<dim> &out) const {
+void ProblemBase<dim>::fill_data (DataOut<dim> &out) const
+{
   Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
   
   out.clear ();
@@ -181,7 +177,8 @@ void ProblemBase<dim>::fill_data (DataOut<dim> &out) const {
 
 
 template <int dim>
-std::string ProblemBase<dim>::get_solution_name () const {
+std::string ProblemBase<dim>::get_solution_name () const
+{
   return "solution";
 };
 
index aff5d97d6fcc72363ab7231d4aec7c39e8da5873..7902c30df498e2451b898c78c77de3bf7a2217fc 100644 (file)
 #include <grid/tria_iterator.h>
 #include <fe/fe.h>
 #include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+
+//TODO: Do this more clever
+static MappingQ1<deal_II_dimension> mapping;
 
 #ifdef DEAL_II_USE_MT
 #include <base/thread_management.h>
@@ -275,7 +279,7 @@ void DataOut<dim>::build_some_patches (Data data)
   QTrapez<1>     q_trapez;
   QIterated<dim> patch_points (q_trapez, data.n_subdivisions);
   
-  FEValues<dim> fe_patch_values(dofs->get_fe(),
+  FEValues<dim> fe_patch_values(mapping, dofs->get_fe(),
                                patch_points,
                                update_values);
 
index bea1beb2a6cdc0d7b0c3321138109e2b9668b38d..f29a1f763d04fcd4b957dd491fac70984937c3af 100644 (file)
 #include <fe/fe.h>
 #include <fe/fe_values.h>
 
+#include <fe/mapping_q1.h>
+
+//TODO: Do this more clever
+static MappingQ1<deal_II_dimension> mapping;
+
 #ifdef DEAL_II_USE_MT
 #include <base/thread_management.h>
 #endif
@@ -36,7 +41,8 @@ void DataOutFaces<dim>::build_some_patches (Data data)
   QTrapez<1>        q_trapez;
   QIterated<dim-1>  patch_points (q_trapez, data.n_subdivisions);
   
-  FEFaceValues<dim> fe_patch_values(dofs->get_fe(),
+  FEFaceValues<dim> fe_patch_values(mapping,
+                                   dofs->get_fe(),
                                    patch_points,
                                    update_values);
 
index 2709ad609b8b9432df9206dfca4a69100dc865a1..9142d73b4ae542c72aa7e7038e4adc0f05559d86 100644 (file)
 #include <dofs/dof_handler.h>
 #include <dofs/dof_accessor.h>
 #include <grid/tria_iterator.h>
+#include <fe/mapping_q1.h>
 #include <fe/fe.h>
 #include <fe/fe_values.h>
 
+//TODO: Do this more clever
+static MappingQ1<deal_II_dimension> mapping;
+
 #ifdef DEAL_II_USE_MT
 #include <base/thread_management.h>
 #endif
@@ -35,15 +39,14 @@ using namespace std;
 #endif
 
 
-
-
 template <int dim>
 void DataOutRotation<dim>::build_some_patches (Data data)
 {
   QTrapez<1>     q_trapez;
   QIterated<dim> patch_points (q_trapez, data.n_subdivisions);
   
-  FEValues<dim> fe_patch_values(dofs->get_fe(),
+  FEValues<dim> fe_patch_values(mapping,
+                               dofs->get_fe(),
                                patch_points,
                                update_values);
 
index 0bfbd648755ab97dab8a0b610bd3ab3901902b02..6aeee12a0ee264558d118c09a4d9cb7b38294a2e 100644 (file)
@@ -21,6 +21,8 @@
 #include <grid/tria_iterator.h>
 #include <fe/fe.h>
 #include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+
 
 
 // if necessary try to work around a bug in the IBM xlC compiler
@@ -29,6 +31,9 @@ using namespace std;
 #endif
 
 
+//TODO: Comment?? Use proper mapping!
+static MappingQ1<deal_II_dimension> mapping;
+
 
 template <int dim>
 unsigned int
@@ -77,8 +82,8 @@ void DataOutStack<dim>::attach_dof_handler (const DoFHandler<dim> &dof)
 
 
 template <int dim>
-void DataOutStack<dim>::declare_data_vector (const std::string     &name,
-                                            const VectorType  vector_type)
+void DataOutStack<dim>::declare_data_vector (const std::string &name,
+                                            const VectorType   vector_type)
 {
   std::vector<std::string> names;
   names.push_back (name);
@@ -88,7 +93,7 @@ void DataOutStack<dim>::declare_data_vector (const std::string     &name,
 
 template <int dim>
 void DataOutStack<dim>::declare_data_vector (const std::vector<std::string> &names,
-                                            const VectorType      vector_type)
+                                            const VectorType    vector_type)
 {
                                   // make sure this function is
                                   // not called after some parameter
@@ -129,7 +134,7 @@ void DataOutStack<dim>::declare_data_vector (const std::vector<std::string> &nam
 template <int dim>
 template <typename number>
 void DataOutStack<dim>::add_data_vector (const Vector<number> &vec,
-                                        const std::string         &name)
+                                        const std::string    &name)
 {
   std::vector<std::string> names;
   names.push_back (name);
@@ -165,7 +170,7 @@ void DataOutStack<dim>::add_data_vector (const Vector<number> &vec,
          {
            data_vector->data.reinit (vec.size());
            std::copy (vec.begin(), vec.end(),
-                 data_vector->data.begin());
+                      data_vector->data.begin());
            break;
          };
       Assert (data_vector != dof_data.end(),
@@ -218,13 +223,13 @@ void DataOutStack<dim>::build_patches (const unsigned int n_subdivisions)
                                   // cell to these points
   QTrapez<1>     q_trapez;
   QIterated<dim> patch_points (q_trapez, n_subdivisions);
-  FEValues<dim>  fe_patch_values (dof_handler->get_fe(),
+  FEValues<dim>  fe_patch_values (mapping, dof_handler->get_fe(),
                                  patch_points,
                                  update_values);
   const unsigned int n_q_points = patch_points.n_quadrature_points;
   std::vector<double>          patch_values (n_q_points);
   std::vector<Vector<double> > patch_values_system (n_q_points,
-                                              Vector<double>(n_components));
+                                                   Vector<double>(n_components));
 
                                   // add the required number of patches
   DataOutBase::Patch<dim+1>  default_patch;
@@ -368,7 +373,7 @@ std::vector<std::string> DataOutStack<dim>::get_dataset_names () const
 // explicit instantiations
 template class DataOutStack<deal_II_dimension>;
 template void DataOutStack<deal_II_dimension>::add_data_vector (const Vector<double> &vec,
-                                                               const std::string         &name);
+                                                               const std::string    &name);
 template void DataOutStack<deal_II_dimension>::add_data_vector (const Vector<float>  &vec,
-                                                               const std::string         &name);
+                                                               const std::string    &name);
 
index d8f400fe11d1c316833b5916a27023bcdb93b3aa..7f534b38e4c296615b5d3887816d6945dbc4713d 100644 (file)
@@ -1,4 +1,4 @@
-//----------------------------  derivative_approximation.cc  ---------------------------
+//----------------------------  gradient_estimator.cc  ---------------------------
 //    $Id$
 //    Version: $Name$
 //
 #include <dofs/dof_handler.h>
 #include <fe/fe.h>
 #include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
 #include <numerics/derivative_approximation.h>
 
 #include <cmath>
 
 
+//TODO: Comment? Use proper mapping!
+static MappingQ1<deal_II_dimension> mapping;
+
 
 template <typename T>
 static T sqr (const T t)
@@ -257,7 +261,7 @@ DerivativeApproximation::approximate (const DoFHandler<dim> &dof_handler,
                                      Vector<float>         &derivative_norm)
 {
   QMidpoint<dim> midpoint_rule;
-  FEValues<dim>  fe_midpoint_value (dof_handler.get_fe(),
+  FEValues<dim>  fe_midpoint_value (mapping, dof_handler.get_fe(),
                                    midpoint_rule,
                                    UpdateFlags(DerivativeDescription::update_flags |
                                                update_q_points));
index a026e12f070b8416ec86a78238220e606e339614..d66742d6c12581abb68a8283726e5561a66446de 100644 (file)
 //----------------------------  error_estimator.cc  ---------------------------
 
 
-#include <fe/fe.h>
-#include <fe/fe_values.h>
-#include <fe/fe_update_flags.h>
+#include <base/timer.h>
 #include <base/thread_management.h>
 #include <base/quadrature.h>
 #include <base/quadrature_lib.h>
-#include <numerics/error_estimator.h>
-#include <dofs/dof_handler.h>
+#include <lac/vector.h>
 #include <grid/tria_iterator.h>
-#include <dofs/dof_accessor.h>
 #include <grid/geometry_info.h>
-#include <lac/vector.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe.h>
+#include <fe/fe_values.h>
+#include <fe/fe_update_flags.h>
+#include <fe/mapping_q1.h>
+#include <numerics/error_estimator.h>
 
 #include <numeric>
 #include <algorithm>
 #include <cmath>
 #include <vector>
-#include <base/timer.h>
-
 
 // if necessary try to work around a bug in the IBM xlC compiler
 #ifdef XLC_WORK_AROUND_STD_BUG
@@ -38,6 +38,12 @@ using namespace std;
 #endif
 
 
+
+//TODO: Comment?? Proper Mapping
+static MappingQ1<deal_II_dimension> mapping;
+
+
+
 static
 inline
 double sqr (const double x)
@@ -46,8 +52,6 @@ double sqr (const double x)
 };
 
 
-
-
 #if deal_II_dimension == 1
 
 template <>
@@ -104,8 +108,7 @@ KellyErrorEstimator<dim>::Data::Data(const DoFHandler<dim>               &dof_ha
          ExcInvalidBoundaryIndicator());
   
   for (typename FunctionMap::const_iterator i=neumann_bc.begin(); i!=neumann_bc.end(); ++i)
-    Assert (i->second->n_components == n_components, ExcInvalidBoundaryFunction());
-  
+    Assert (i->second->n_components == n_components, ExcInvalidBoundaryFunction());  
   
                                   // Init the size of a lot of vectors
                                   // needed in the calculations once
@@ -151,7 +154,7 @@ void KellyErrorEstimator<dim>::estimate (const DoFHandler<dim>   &dof_handler,
                                         const FunctionMap       &neumann_bc,
                                         const Vector<double>    &solution,
                                         Vector<float>           &error,
-                                        const std::vector<bool>      &component_mask,
+                                        const std::vector<bool> &component_mask,
                                         const Function<dim>     *coefficients,
                                         unsigned int             n_threads)
 {
@@ -201,8 +204,8 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1>                 &dof_
   
                                   // if no mask given: treat all components
   std::vector<bool> component_mask ((component_mask_.size() == 0)    ?
-                              std::vector<bool>(n_components, true) :
-                              component_mask_);
+                                   std::vector<bool>(n_components, true) :
+                                   component_mask_);
   Assert (component_mask.size() == n_components, ExcInvalidComponentMask());
   Assert (count(component_mask.begin(), component_mask.end(), true) > 0,
          ExcInvalidComponentMask());
@@ -254,7 +257,7 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1>                 &dof_
                                   // two contributions from the two
                                   // vertices of each cell.
   QTrapez<1> quadrature;
-  FEValues<1> fe_values (dof_handler.get_fe(), quadrature, update_gradients);
+  FEValues<1> fe_values (mapping, dof_handler.get_fe(), quadrature, update_gradients);
   active_cell_iterator cell = dof_handler.begin_active();
   for (unsigned int cell_index=0; cell != dof_handler.end(); ++cell, ++cell_index)
     {
@@ -373,7 +376,8 @@ void KellyErrorEstimator<dim>::estimate_some (Data               &data,
                                   // need not compute all values on the
                                   // neighbor cells, so using two objects
                                   // gives us a performance gain).
-  FEFaceValues<dim> fe_face_values_cell (data.dof_handler.get_fe(),
+  FEFaceValues<dim> fe_face_values_cell (mapping,
+                                        data.dof_handler.get_fe(),
                                         data.quadrature,
                                         UpdateFlags(update_gradients      |
                                                     update_JxW_values     |
@@ -381,10 +385,12 @@ void KellyErrorEstimator<dim>::estimate_some (Data               &data,
                                                       (data.coefficients != 0))  ?
                                                      update_q_points : 0) |
                                                     update_normal_vectors)); 
-  FEFaceValues<dim> fe_face_values_neighbor (data.dof_handler.get_fe(),
+  FEFaceValues<dim> fe_face_values_neighbor (mapping,
+                                            data.dof_handler.get_fe(),
                                             data.quadrature,
                                             update_gradients);
-  FESubfaceValues<dim> fe_subface_values (data.dof_handler.get_fe(),
+  FESubfaceValues<dim> fe_subface_values (mapping,
+                                         data.dof_handler.get_fe(),
                                          data.quadrature,
                                          update_gradients);
 
@@ -564,7 +570,7 @@ void KellyErrorEstimator<dim>::estimate (const DoFHandler<dim>               &do
                                   solutions,
                                   ((component_mask.size() == 0)    ?
                                    std::vector<bool>(dof_handler.get_fe().n_components(),
-                                                true) :
+                                                     true) :
                                    component_mask),
                                   coefficients,
                                   n_threads,
@@ -989,8 +995,8 @@ integrate_over_irregular_face (Data                       &data,
              data.face_integrals.end(),
              ExcInternalError());
       Assert (data.face_integrals[face->child(subface_no)][0] >= 0,
-             ExcInternalError());
-
+             ExcInternalError());
+      
       for (unsigned int n=0; n<n_solution_vectors; ++n)
        sum[n] += data.face_integrals[face->child(subface_no)][n];
     };
index ec2156ee8f629df93de1a0db85ebaf3bb5a4ce2a..b6a0cee4b93e7037bb829281b5022e1ea76c6017 100644 (file)
 #include <lac/block_vector.h>
 #include <lac/sparse_matrix.h>
 #include <lac/block_sparse_matrix.h>
+#include <fe/mapping_q1.h>
 
 #include <algorithm>
 #include <set>
 #include <cmath>
 
+
 // if necessary try to work around a bug in the IBM xlC compiler
 #ifdef XLC_WORK_AROUND_STD_BUG
 using namespace std;
 #endif
 
 
+
+
+//TODO: Comment?? Use proper mapping!
+static MappingQ1<deal_II_dimension> mapping;
+
+
+
+//TODO: re-create the create_mass_matrix function with 2 args
+
 template <int dim>
 void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
                                             const Quadrature<dim>    &q,
@@ -96,42 +107,16 @@ void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
 
 
 
-template <int dim>
-void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
-                                            SparseMatrix<double>     &matrix)
-{
-  const FiniteElement<dim> &fe = dof.get_fe();
-
-  const unsigned int dofs_per_cell = fe.dofs_per_cell;
-  
-  FullMatrix<double>   local_mass_matrix (dofs_per_cell, dofs_per_cell);
-  std::vector<unsigned int> dofs_on_this_cell (dofs_per_cell);
-  
-  DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
-                                       endc = dof.end();
-  for (; cell!=endc; ++cell) 
-    {
-      cell->get_dof_indices (dofs_on_this_cell);
-      fe.get_local_mass_matrix (cell, local_mass_matrix);
-      
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       for (unsigned int j=0; j<dofs_per_cell; ++j)
-         matrix.add (dofs_on_this_cell[i], dofs_on_this_cell[j],
-                     local_mass_matrix(i,j));
-    };
-};
-
-
 #if deal_II_dimension == 1
 
 template <>
-void MatrixCreator<1>::create_boundary_mass_matrix (const DoFHandler<1>    &,
-                                                   const Quadrature<0>    &,
-                                                   SparseMatrix<double>   &,
-                                                   const FunctionMap      &,
-                                                   Vector<double>         &,
-                                                   std::vector<unsigned int>   &,
-                                                   const Function<1>      *)
+void MatrixCreator<1>::create_boundary_mass_matrix (const DoFHandler<1>       &,
+                                                   const Quadrature<0>       &,
+                                                   SparseMatrix<double>      &,
+                                                   const FunctionMap         &,
+                                                   Vector<double>            &,
+                                                   std::vector<unsigned int> &,
+                                                   const Function<1>         *)
 {
   Assert (false, ExcNotImplemented());
 };
@@ -140,13 +125,13 @@ void MatrixCreator<1>::create_boundary_mass_matrix (const DoFHandler<1>    &,
 
 
 template <int dim>
-void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>    &dof,
-                                                     const Quadrature<dim-1>  &q,
-                                                     SparseMatrix<double>     &matrix,
-                                                     const FunctionMap        &rhs,
-                                                     Vector<double>           &rhs_vector,
-                                                     std::vector<unsigned int>     &dof_to_boundary_mapping,
-                                                     const Function<dim>      *a)
+void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>     &dof,
+                                                     const Quadrature<dim-1>   &q,
+                                                     SparseMatrix<double>      &matrix,
+                                                     const FunctionMap         &rhs,
+                                                     Vector<double>            &rhs_vector,
+                                                     std::vector<unsigned int> &dof_to_boundary_mapping,
+                                                     const Function<dim>       *a)
 {
   const FiniteElement<dim> &fe = dof.get_fe();
   const unsigned int n_components  = fe.n_components();
@@ -184,7 +169,7 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>    &
   UpdateFlags update_flags = UpdateFlags (update_values     |
                                          update_JxW_values |
                                          update_q_points);
-  FEFaceValues<dim> fe_values (fe, q, update_flags);
+  FEFaceValues<dim> fe_values (mapping, fe, q, update_flags);
 
                                   // two variables for the coefficient,
                                   // one for the two cases indicated in
@@ -452,6 +437,8 @@ void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
 
 
 
+
+//TODO: recreate this function
 /*
 
 template <int dim>
@@ -488,6 +475,8 @@ void MatrixCreator<dim>::create_level_laplace_matrix (unsigned int level,
 
 
 
+
+
 template <int dim>
 void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
                                                const Quadrature<dim>    &q,
@@ -525,10 +514,10 @@ template <int dim>
 template <typename number>
 void
 MatrixTools<dim>::apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
-                                        SparseMatrix<number> &matrix,
-                                        Vector<number>       &solution,
-                                        Vector<number>       &right_hand_side,
-                                        const bool            preserve_symmetry)
+                                        SparseMatrix<number>  &matrix,
+                                        Vector<number>   &solution,
+                                        Vector<number>   &right_hand_side,
+                                        const bool        preserve_symmetry)
 {
   Assert (matrix.n() == matrix.m(),
          ExcDimensionsDontMatch(matrix.n(), matrix.m()));
@@ -542,8 +531,8 @@ MatrixTools<dim>::apply_boundary_values (const std::map<unsigned int,double> &bo
     return;
 
 
-  std::map<unsigned int,double>::const_iterator  dof  = boundary_values.begin(),
-                                           endd = boundary_values.end();
+  std::map<unsigned int,double>::const_iterator dof  = boundary_values.begin(),
+                                               endd = boundary_values.end();
   const unsigned int n_dofs             = matrix.m();
   const SparsityPattern    &sparsity    = matrix.get_sparsity_pattern();
   const unsigned int *sparsity_rowstart = sparsity.get_rowstart_indices();
@@ -695,7 +684,7 @@ MatrixTools<dim>::apply_boundary_values (const std::map<unsigned int,double> &bo
                                         BlockSparseMatrix<double>  &matrix,
                                         BlockVector<double>   &solution,
                                         BlockVector<double>   &right_hand_side,
-                                        const bool                    preserve_symmetry)
+                                        const bool             preserve_symmetry)
 {
   const unsigned int blocks = matrix.n_block_rows();
   
@@ -724,8 +713,8 @@ MatrixTools<dim>::apply_boundary_values (const std::map<unsigned int,double> &bo
     return;
 
 
-  std::map<unsigned int,double>::const_iterator  dof  = boundary_values.begin(),
-                                           endd = boundary_values.end();
+  std::map<unsigned int,double>::const_iterator dof  = boundary_values.begin(),
+                                               endd = boundary_values.end();
   const unsigned int n_dofs = matrix.m();
   const BlockSparsityPattern &
     sparsity_pattern = matrix.get_sparsity_pattern();
@@ -949,17 +938,17 @@ MatrixTools<dim>::apply_boundary_values (const std::map<unsigned int,double> &bo
                            [this_sparsity.get_rowstart_indices()[row]];
                      else
                        p = std::lower_bound(&this_sparsity.get_column_numbers()
-                                       [this_sparsity.get_rowstart_indices()[row]+1],
-                                       &this_sparsity.get_column_numbers()
-                                       [this_sparsity.get_rowstart_indices()[row+1]],
-                                       block_index.second);
+                                            [this_sparsity.get_rowstart_indices()[row]+1],
+                                            &this_sparsity.get_column_numbers()
+                                            [this_sparsity.get_rowstart_indices()[row+1]],
+                                            block_index.second);
                    }
                  else
                    p = std::lower_bound(&this_sparsity.get_column_numbers()
-                                   [this_sparsity.get_rowstart_indices()[row]],
-                                   &this_sparsity.get_column_numbers()
-                                   [this_sparsity.get_rowstart_indices()[row+1]],
-                                   block_index.second);
+                                        [this_sparsity.get_rowstart_indices()[row]],
+                                        &this_sparsity.get_column_numbers()
+                                        [this_sparsity.get_rowstart_indices()[row+1]],
+                                        block_index.second);
 
                                                   // check whether this line has
                                                   // an entry in the regarding column
@@ -1205,6 +1194,8 @@ void MassMatrix<dim>::assemble (Vector<double>      &rhs,
 
 
 
+
+
 template <int dim>
 LaplaceMatrix<dim>::LaplaceMatrix (const Function<dim> * const rhs,
                                   const Function<dim> * const a) :
index b99173ebc91da0f9eb1229e50f3ebb5f11009441..69f15b5e67e2811a34a4c1ea7dd06974130f3955 100644 (file)
@@ -29,6 +29,7 @@
 #include <lac/precondition.h>
 #include <lac/solver_cg.h>
 #include <lac/vector_memory.h>
+#include <fe/mapping_q1.h>
 
 #include <numeric>
 #include <algorithm>
@@ -42,6 +43,12 @@ using namespace std;
 #endif
 
 
+
+//TODO: Comment?? Use proper mapping!
+static MappingQ1<deal_II_dimension> mapping;
+
+
+
 static inline double sqr (const double x)
 {
   return x*x;
@@ -57,6 +64,7 @@ inline double sqr_point (const Tensor<1,dim> &p)
 
 
 
+
 template <int dim>
 void VectorTools::interpolate (const DoFHandler<dim> &dof,
                               const Function<dim>   &function,
@@ -83,89 +91,10 @@ void VectorTools::interpolate (const DoFHandler<dim> &dof,
                                   // avoided to evaluate
                                   // the vectorfunction multiply at
                                   // the same point on a cell.
-  std::vector<Point<dim> > unit_support_points (fe.dofs_per_cell);
+  std::vector<Point<dim> > unit_support_points;
   fe.get_unit_support_points(unit_support_points);
-
-                                  // The following works well
-                                  // if @p{dofs_per_x<=1 (x=vertex,line,cell)}
-                                  // as then
-                                  // the multiple support_points
-                                  // are placed one after another.
-
-                                  // find the support points 
-                                  // on a cell that
-                                  // are multiply mentioned in 
-                                  // @p{unit_support_points}.
-                                  // Mark the first representative
-                                  // of each multiply mentioned
-                                  // support point by setting
-                                  // @p{true} in the boolean vector 
-                                  // @p{is_representative_point}.
-//   vector<bool>  is_representative_point(fe.dofs_per_cell, false);
-//   is_representative_point[0]=true;
-//   unsigned int n_rep_points=1;
-//   for (unsigned int last_rep_point=0, i=1; i<fe.dofs_per_cell; ++i)
-//     {
-//       if (unit_support_points[i] != unit_support_points[last_rep_point])
-//     {
-//       is_representative_point[i] = true;
-//       last_rep_point=i;
-//       ++n_rep_points;
-//     }
-//    };
-
-//   vector<unsigned int> dofs_on_cell (fe.dofs_per_cell);
-//   vector<Point<dim> >  support_points (fe.dofs_per_cell);
-
-//   vector<Point<dim> > rep_points (n_rep_points);
-//   vector<Vector<double> > function_values_at_rep_points (
-//     n_rep_points, Vector<double>(fe.n_components()));
-
-//   for (; cell!=endc; ++cell)
-//     {
-//                                    // for each cell:
-//                                    // get location of finite element
-//                                    // off-points (support_points)
-//       fe.get_support_points (cell, support_points);
-
-//                                    // pick out the representative
-//                                    // support points
-//       unsigned int j=0;
-//       for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-//     if (is_representative_point[i])
-//       rep_points[j++]=support_points[i];
-//       Assert(j == n_rep_points, ExcInternalError());
-
-//                                    // get function values at these points
-//       vectorfunction.value_list (rep_points, function_values_at_rep_points);
-  
-//                                          // get indices of the dofs on this cell
-//       cell->get_dof_indices (dofs_on_cell);
-
-//                                    // distribute function values to the
-//                                    // whole vector
-//       int last_rep_point = -1;
-//                                    // it holds `is_representative_point[0]=true'
-//                                    // therefore the first @p{last_rep_point} is 0
-//                                    // and we need to start with
-//                                    // `last_rep_point = -1'
-//       for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-//     {
-//       if (is_representative_point[i])
-//         ++last_rep_point;
-
-//       const unsigned int component
-//         = fe.system_to_component_index(i).first;
-//       vec(dofs_on_cell[i])
-//         = function_values_at_rep_points[last_rep_point](component);
-//     } 
-//     }
-
-                                  // The following is more general.
-                                  // It also works if @p{dofs_per_x>1},
-                                  // i.e. it is usable also for systems
-                                  // including
-                                  // FEQ3, FEQ4, FEDG_Qx.
+  Assert (unit_support_points.size() != 0,
+         ExcNonInterpolatingFE());
 
                                   // Find the support points 
                                   // on a cell that
@@ -222,7 +151,7 @@ void VectorTools::interpolate (const DoFHandler<dim> &dof,
   Assert(dof_to_rep_index_table.size()==fe.dofs_per_cell, ExcInternalError());
 
   std::vector<unsigned int> dofs_on_cell (fe.dofs_per_cell);
-  std::vector<Point<dim> >  support_points (fe.dofs_per_cell);
+  std::vector<double>       dummy_weights (fe.dofs_per_cell);
 
   std::vector<Point<dim> >  rep_points (n_rep_points);
 
@@ -236,12 +165,23 @@ void VectorTools::interpolate (const DoFHandler<dim> &dof,
   std::vector<Vector<double> > function_values_system (n_rep_points,
                                                  Vector<double>(fe.n_components()));
 
+                                  // Make a quadrature rule from support points
+                                  // to feed it into FEValues
+  Quadrature<dim> support_quadrature(unit_support_points, dummy_weights);
+
+//TODO: Higher order mapping?  
+                                  // Transformed support points are computed by
+                                  // FEValues
+  FEValues<dim> fe_values (mapping, fe, support_quadrature, update_q_points);
+  
   for (; cell!=endc; ++cell)
     {
                                       // for each cell:
                                       // get location of finite element
-                                      // off-points (support_points)
-      fe.get_support_points (cell, support_points);
+                                      // support_points
+      fe_values.reinit(cell);
+      const std::vector<Point<dim> >& support_points =
+       fe_values.get_quadrature_points();
       
                                       // pick out the representative
                                       // support points
@@ -289,6 +229,7 @@ void VectorTools::interpolate (const DoFHandler<dim> &dof,
 
 
 
+
 template <int dim> void
 VectorTools::interpolate (const DoFHandler<dim>           &dof_1,
                          const DoFHandler<dim>           &dof_2,
@@ -452,17 +393,7 @@ void VectorTools::project (const DoFHandler<dim>    &dof,
   SparseMatrix<double> mass_matrix (sparsity);
   Vector<double> tmp (mass_matrix.n());
 
-                                  // try to assemble the mass matrix by exact
-                                  // integration. if this is not supported,
-                                  // then use quadrature
-  try 
-    {
-      MatrixCreator<dim>::create_mass_matrix (dof, mass_matrix);
-    }
-  catch (FiniteElement<dim>::ExcComputationNotUseful)
-    {
-      MatrixCreator<dim>::create_mass_matrix (dof, quadrature, mass_matrix);
-    };
+  MatrixCreator<dim>::create_mass_matrix (dof, quadrature, mass_matrix);
   
   VectorTools::create_right_hand_side (dof, quadrature, function, tmp);
 
@@ -504,7 +435,7 @@ void VectorTools::create_right_hand_side (const DoFHandler<dim>    &dof_handler,
   UpdateFlags update_flags = UpdateFlags(update_values   |
                                         update_q_points |
                                         update_JxW_values);
-  FEValues<dim> fe_values (fe, quadrature, update_flags);
+  FEValues<dim> fe_values (mapping, fe, quadrature, update_flags);
 
   const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
                     n_q_points    = fe_values.n_quadrature_points,
@@ -513,8 +444,8 @@ void VectorTools::create_right_hand_side (const DoFHandler<dim>    &dof_handler,
   std::vector<unsigned int> dofs (dofs_per_cell);
   Vector<double> cell_vector (dofs_per_cell);
 
-  DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
-                                       endc = dof_handler.end();
+  typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+                                                endc = dof_handler.end();
 
   if (n_components==1)
     {
@@ -524,8 +455,8 @@ void VectorTools::create_right_hand_side (const DoFHandler<dim>    &dof_handler,
        {
          fe_values.reinit(cell);
          
-         const FullMatrix<double> &values    = fe_values.get_shape_values ();
-         const std::vector<double>     &weights   = fe_values.get_JxW_values ();
+         const FullMatrix<double>  &values    = fe_values.get_shape_values ();
+         const std::vector<double> &weights   = fe_values.get_JxW_values ();
          rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values);
          
          cell_vector.clear();
@@ -550,8 +481,8 @@ void VectorTools::create_right_hand_side (const DoFHandler<dim>    &dof_handler,
        {
          fe_values.reinit(cell);
          
-         const FullMatrix<double> &values    = fe_values.get_shape_values ();
-         const std::vector<double>     &weights   = fe_values.get_JxW_values ();
+         const FullMatrix<double>  &values    = fe_values.get_shape_values ();
+         const std::vector<double> &weights   = fe_values.get_JxW_values ();
          rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
          
          cell_vector.clear();
@@ -597,8 +528,8 @@ VectorTools::interpolate_boundary_values (const DoFHandler<1>      &dof,
                                   // the original value or a vector
                                   // of @p{true}s
   const std::vector<bool> component_mask ((component_mask_.size() == 0) ?
-                                    std::vector<bool> (fe.n_components(), true) :
-                                    component_mask_);
+                                         std::vector<bool> (fe.n_components(), true) :
+                                         component_mask_);
   Assert (count(component_mask.begin(), component_mask.end(), true) > 0,
          ExcComponentMismatch());
   
@@ -675,14 +606,15 @@ VectorTools::interpolate_boundary_values (const DoFHandler<dim>    &dof,
                                   // the original value or a vector
                                   // of @p{true}s
   const std::vector<bool> component_mask ((component_mask_.size() == 0) ?
-                                    std::vector<bool> (fe.n_components(), true) :
-                                    component_mask_);
+                                         std::vector<bool> (fe.n_components(), true) :
+                                         component_mask_);
   Assert (count(component_mask.begin(), component_mask.end(), true) > 0,
          ExcComponentMismatch());
 
-                                  // field to store the indices of dofs
+                                  // field to store the indices
   std::vector<unsigned int> face_dofs (fe.dofs_per_face, -1);
   std::vector<Point<dim> >  dof_locations (face_dofs.size(), Point<dim>());
+  
                                   // array to store the values of
                                   // the boundary function at the
                                   // boundary points. have to arrays
@@ -691,55 +623,78 @@ VectorTools::interpolate_boundary_values (const DoFHandler<dim>    &dof,
                                   // respectively
   std::vector<double>          dof_values_scalar (fe.dofs_per_face);
   std::vector<Vector<double> > dof_values_system (fe.dofs_per_face,
-                                            Vector<double>(fe.n_components()));
-       
-  DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
-                                       endf = dof.end_face();
-  for (; face!=endf; ++face)
-    if (boundary_component == face->boundary_indicator()) 
-                                      // face is of the right component
-      {
-                                        // get indices, physical location and
-                                        // boundary values of dofs on this
-                                        // face
-       face->get_dof_indices (face_dofs);
-       fe.get_face_support_points (face, dof_locations);
+                                                 Vector<double>(fe.n_components()));
 
-       if (fe_is_system)
+                                  // next generate a quadrature rule
+                                  // on the face from the unit
+                                  // support points. this wil be used
+                                  // to obtain the quadrature points
+                                  // on the real cell's face
+  typename std::vector<Point<dim-1> > unit_support_points;
+  fe.get_unit_face_support_points(unit_support_points);
+                                  // check whether there are support
+                                  // points on the face, if not, then
+                                  // this FE does not allow to be
+                                  // used in this function
+  Assert (unit_support_points.size() != 0, ExcNonInterpolatingFE());
+
+  std::vector<double> dummy_weights (unit_support_points.size());
+  Quadrature<dim-1> aux_quad (unit_support_points, dummy_weights);
+  FEFaceValues<dim> fe_values (mapping, fe, aux_quad, update_q_points);
+
+  DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
+                                       endc = dof.end();
+  for (; cell!=endc; ++cell)
+    for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
+        ++face_no)
+      {
+       DoFHandler<dim>::face_iterator face = cell->face(face_no);
+       if (boundary_component == face->boundary_indicator()) 
+                                          // face is of the right component
          {
-           boundary_function.vector_value_list (dof_locations, dof_values_system);
-           
-                                            // enter into list
+                                            // get indices, physical location and
+                                            // boundary values of dofs on this
+                                            // face
+           face->get_dof_indices (face_dofs);
+           fe_values.reinit(cell, face_no);
+           const std::vector<Point<dim> > &dof_locations = fe_values.get_quadrature_points ();
            
-           for (unsigned int i=0; i<face_dofs.size(); ++i)
-             if (component_mask[fe.face_system_to_component_index(i).first])
-               boundary_values[face_dofs[i]]
-                 = dof_values_system[i](fe.face_system_to_component_index(i).first);
+           if (fe_is_system)
+             {
+               boundary_function.vector_value_list (dof_locations, dof_values_system);
+               
+                                                // enter into list
+               
+               for (unsigned int i=0; i<face_dofs.size(); ++i)
+                 if (component_mask[fe.face_system_to_component_index(i).first])
+                   boundary_values[face_dofs[i]]
+                     = dof_values_system[i](fe.face_system_to_component_index(i).first);
+             }
+           else
+                                              // fe has only one component,
+                                              // so save some computations
+             {
+                                                // get only the one component that
+                                                // this function has
+               boundary_function.value_list (dof_locations,
+                                             dof_values_scalar,
+                                             0);
+               
+                                                // enter into list
+               
+               for (unsigned int i=0; i<face_dofs.size(); ++i)
+                 boundary_values[face_dofs[i]] = dof_values_scalar[i];
+             }
          }
-       else
-                                          // fe has only one component,
-                                          // so save some computations
-         {
-                                            // get only the one component that
-                                            // this function has
-           boundary_function.value_list (dof_locations,
-                                         dof_values_scalar,
-                                         0);
-           
-                                            // enter into list
-           
-           for (unsigned int i=0; i<face_dofs.size(); ++i)
-             boundary_values[face_dofs[i]] = dof_values_scalar[i];
-         };
-      };
+      }
 }
 
-
-
+  
+  
 template <int dim>
 void
 VectorTools::project_boundary_values (const DoFHandler<dim>    &dof,
-                                     const typename std::map<unsigned char,const Function<dim>*> &boundary_functions,
+                                     const std::map<unsigned char,const Function<dim>*> &boundary_functions,
                                      const Quadrature<dim-1>  &q,
                                      std::map<unsigned int,double> &boundary_values)
 {
@@ -858,19 +813,19 @@ VectorTools::integrate_difference (const DoFHandler<dim>    &dof,
   if ((norm==H1_seminorm) || (norm==H1_norm))
     update_flags = UpdateFlags (update_flags | update_gradients);
   
-  FEValues<dim> fe_values(fe, q, update_flags);
+  FEValues<dim> fe_values(mapping, fe, q, update_flags);
 
   std::vector< Vector<double> >        function_values (n_q_points,
-                                                  Vector<double>(n_components));
+                                                       Vector<double>(n_components));
   std::vector<std::vector<Tensor<1,dim> > > function_grads (n_q_points,
-                                                 std::vector<Tensor<1,dim> >(n_components));
+                                                           std::vector<Tensor<1,dim> >(n_components));
   std::vector<double> weight_values (n_q_points);
   std::vector<Vector<double> > weight_vectors (n_q_points, n_components);
   
   std::vector<Vector<double> >         psi_values (n_q_points,
-                                             Vector<double>(n_components));
+                                                  Vector<double>(n_components));
   std::vector<std::vector<Tensor<1,dim> > > psi_grads (n_q_points,
-                                            std::vector<Tensor<1,dim> >(n_components));
+                                                      std::vector<Tensor<1,dim> >(n_components));
   std::vector<double> psi_scalar (n_q_points);
                                   // tmp vector when we use the
                                   // Function<dim> functions for
@@ -913,15 +868,15 @@ VectorTools::integrate_difference (const DoFHandler<dim>    &dof,
                                             // the function really has only
                                             // one component
            if (fe_is_system)
-               exact_solution.vector_value_list (fe_values.get_quadrature_points(),
-                                                 psi_values);
+             exact_solution.vector_value_list (fe_values.get_quadrature_points(),
+                                               psi_values);
            else
              {
                exact_solution.value_list (fe_values.get_quadrature_points(),
                                           tmp_values);
                for (unsigned int i=0; i<n_q_points; ++i)
                  psi_values[i](0) = tmp_values[i];
-             }
+             };
            
                                             // then subtract finite element
                                             // fe_function
@@ -1109,7 +1064,7 @@ VectorTools::integrate_difference (const DoFHandler<dim>    &dof,
                psi_grads[q][k] -= function_grads[q][k];
 
                                             // take square of integrand
-           fill_n (psi_scalar.begin(), n_q_points, 0.0);
+           std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
 
            for (unsigned int k=0; k<n_components; ++k)
              if (weight != 0)
@@ -1134,9 +1089,9 @@ VectorTools::integrate_difference (const DoFHandler<dim>    &dof,
 
                                             // add seminorm to L_2 norm or
                                             // to zero
-           diff += inner_product (psi_scalar.begin(), psi_scalar.end(),
-                                  fe_values.get_JxW_values().begin(),
-                                  0.0);
+           diff += std::inner_product (psi_scalar.begin(), psi_scalar.end(),
+                                       fe_values.get_JxW_values().begin(),
+                                       0.0);
            diff = sqrt(diff);
 
            break;
@@ -1165,20 +1120,18 @@ VectorTools::compute_mean_value (const DoFHandler<dim> &dof,
   Assert (component < dof.get_fe().n_components(),
          ExcIndexRange(component, 0, dof.get_fe().n_components()));
   
-  FEValues<dim> fe(dof.get_fe(), quadrature,
+  FEValues<dim> fe(mapping, dof.get_fe(), quadrature,
                   UpdateFlags(update_JxW_values
                               | update_values));
 
-  DoFHandler<dim>::active_cell_iterator c;
+  typename DoFHandler<dim>::active_cell_iterator c;
   std::vector<Vector<double> > values(quadrature.n_quadrature_points,
-                                Vector<double> (dof.get_fe().n_components()));
+                                     Vector<double> (dof.get_fe().n_components()));
   
   double mean = 0.;
   double area = 0.;
                                   // Compute mean value
-  for (c = dof.begin_active();
-       c != dof.end();
-       ++c)
+  for (c = dof.begin_active(); c != dof.end(); ++c)
     {
       fe.reinit (c);
       fe.get_function_values(v, values);
@@ -1186,10 +1139,13 @@ VectorTools::compute_mean_value (const DoFHandler<dim> &dof,
        {
          mean += fe.JxW(k) * values[k](component);
          area += fe.JxW(k);
-       }
-    }
+       };
+    };
+  
   return (mean/area);
-}
+};
+
+
 
 
 
@@ -1231,8 +1187,6 @@ double VectorTools::compute_mean_value (const DoFHandler<deal_II_dimension> &dof
                                        const unsigned int component);
 
 
-
-
 // the following two functions are not derived from a template in 1d
 // and thus need no explicit instantiation
 #if deal_II_dimension > 1
index 9cec48bc60a42be93b921c4643951780400ddf81..02f9183595de140dc4d67c6322beeeb01021fa00 100644 (file)
@@ -15,6 +15,7 @@
                                 // ("DoF"s) to vertices, lines, and
                                 // cells.
 #include <dofs/dof_handler.h>
+
                                 // The following include contains the
                                 // description of the bilinear finite
                                 // element, including the facts that
@@ -30,7 +31,7 @@
                                 // elements, but not only for two
                                 // space dimensions, but also for one
                                 // and three dimensions.
-#include <fe/fe_lib.lagrange.h>
+#include <fe/fe_q.h>
                                 // In the following file, several
                                 // tools for manipulating degrees of
                                 // freedom can be found:
@@ -144,21 +145,33 @@ void distribute_dofs (DoFHandler<2> &dof_handler)
                                   // object which describes how many
                                   // degrees of freedom are to be
                                   // associated to each of these
-                                  // objects. For (bi-, tri-)linear
-                                  // finite elements, this is done
-                                  // using the FEQ1 class, which
+                                  // objects. This is done using the
+                                  // ``FE_Q'' class. Giving a
+                                   // constructor argument one
                                   // states that one degree of
                                   // freedom is to be assigned to
                                   // each vertex, while there are
                                   // none on lines and inside the
-                                  // quadrilateral. We first need to
-                                  // create an object of this class
-                                  // and use it to distribute the
-                                  // degrees of freedom. Note that
-                                  // the DoFHandler object will store
-                                  // a reference to this object, so
-                                  // we need to make it static as
-                                  // well, in order to prevent its
+                                  // quadrilateral. In fact, the
+                                  // argument denotes the polynomial
+                                  // degree, in this case, we get
+                                  // bilinear finite elements in two
+                                  // space dimensions; a value of,
+                                  // say, three would give us
+                                  // bi-cubic ones. In general, ``FE_Q''
+                                  // denotes the family of continuous
+                                  // elements with complete polynomials
+                                  // (i.e. tensor-product polynomials)
+                                  // up to the specified order 
+                                  //
+                                  // We first need to create an
+                                  // object of this class and use it
+                                  // to distribute the degrees of
+                                  // freedom. Note that the
+                                  // DoFHandler object will store a
+                                  // reference to this object, so we
+                                  // need to make it static as well,
+                                  // in order to prevent its
                                   // preemptive
                                   // destruction. (However, the
                                   // library would warn us about this
@@ -166,7 +179,7 @@ void distribute_dofs (DoFHandler<2> &dof_handler)
                                   // occured. You can check this, if
                                   // you want, by removing the
                                   // 'static' declaration.)
-  static const FEQ1<2> finite_element;
+  static const FE_Q<2> finite_element(1);
   dof_handler.distribute_dofs (finite_element);
 
                                   // Now we have associated a number
index 9c6fdc2d4b4093a297fb816a012c15cc2f609f53..f1f7379396be8a3eac1479d0f0de3ea0df1b99ae 100644 (file)
@@ -23,7 +23,7 @@
 
                                 // In this file are the finite
                                 // element descriptions.
-#include <fe/fe_lib.lagrange.h>
+#include <fe/fe_q.h>
 
                                 // And this file is needed for the
                                 // creation of sparsity patterns of
@@ -118,11 +118,13 @@ class LaplaceProblem
                                     // variables. There are variables
                                     // describing the triangulation
                                     // and the numbering of the
-                                    // degrees of freedom...
-                                    // (FEQ1<2> is the space with
-                                    // shape functions {1,x,y,xy})
+                                    // degrees of freedom (we will
+                                    // specify the exact polynomial
+                                    // degree of the finite element
+                                    // in the constructor of this
+                                    // class)...
     Triangulation<2>     triangulation;
-    FEQ1<2>              fe;
+    FE_Q<2>              fe;
     DoFHandler<2>        dof_handler;
 
                                     // ...variables for the sparsity
@@ -142,14 +144,20 @@ class LaplaceProblem
 
 
                                 // Here comes the constructor. It
-                                // does not much more than associate
-                                // the dof_handler variable to the
+                                // does not much more than first to
+                                // specify that we want bi-linear
+                                // elements (denoted by the parameter
+                                // to the finite element object,
+                                // which specifies the polynomial
+                                // degree), and to associate the
+                                // dof_handler variable to the
                                 // triangulation we use. All the
                                 // other member variables of the
                                 // LaplaceProblem class have a
                                 // default constructor which does all
                                 // we want.
 LaplaceProblem::LaplaceProblem () :
+                fe (1),
                dof_handler (triangulation)
 {};
 
index 8566ce580dd60ecb450aee961e4e02bf91db758f..754ddc63019670d9f84477b73e09cedf3e49831d 100644 (file)
@@ -12,7 +12,7 @@
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
 #include <dofs/dof_accessor.h>
-#include <fe/fe_lib.lagrange.h>
+#include <fe/fe_q.h>
 #include <dofs/dof_tools.h>
 #include <fe/fe_values.h>
 #include <base/quadrature_lib.h>
@@ -72,7 +72,7 @@ class LaplaceProblem
     void output_results () const;
 
     Triangulation<dim>   triangulation;
-    FEQ1<dim>            fe;
+    FE_Q<dim>            fe;
     DoFHandler<dim>      dof_handler;
 
     SparsityPattern      sparsity_pattern;
@@ -199,12 +199,15 @@ double BoundaryValues<dim>::value (const Point<dim> &p,
 
 
                                 // This is the constructor of the
-                                // LaplaceProblem class. It
-                                // associates the DoFHandler to the
+                                // LaplaceProblem class. It specifies
+                                // the desired polynomial degree of
+                                // the finite elements and associates
+                                // the DoFHandler to the
                                 // triangulation just as in the
                                 // previous example.
 template <int dim>
 LaplaceProblem<dim>::LaplaceProblem () :
+                fe (1),
                dof_handler (triangulation)
 {};
 
@@ -414,8 +417,7 @@ void LaplaceProblem<dim>::assemble_system ()
                            right_hand_side.value (fe_values.quadrature_point (q_point)) *
                            fe_values.JxW (q_point));
          };
-
-
+      
                                       // The transfer into the global
                                       // matrix and right hand side
                                       // is done exactly as before,
index 5477e710f346375df58955815891ca47461c9f1a..8a6fcb63f7707ceffc2904ee0ca45ffde2b7a1c7 100644 (file)
@@ -19,7 +19,7 @@
 #include <grid/tria_iterator.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
-#include <fe/fe_lib.lagrange.h>
+#include <fe/fe_q.h>
 #include <fe/fe_values.h>
 #include <numerics/vectors.h>
 #include <numerics/matrices.h>
@@ -69,7 +69,7 @@ class LaplaceProblem
     void output_results (const unsigned int cycle) const;
 
     Triangulation<dim>   triangulation;
-    FEQ1<dim>            fe;
+    FE_Q<dim>            fe;
     DoFHandler<dim>      dof_handler;
 
     SparsityPattern      sparsity_pattern;
@@ -257,6 +257,7 @@ void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
                                 // This function is as before.
 template <int dim>
 LaplaceProblem<dim>::LaplaceProblem () :
+                fe (1),
                dof_handler (triangulation)
 {};
 
index f617d9843e784976ad37203bba979969f45f9650..3682e235205e74983db3cbd9d99d07b7bdf70e97 100644 (file)
 #include <numerics/matrices.h>
 #include <numerics/data_out.h>
                                 // From the following include file we
-                                // will import the declaration of the
-                                // quadratic finite element class,
-                                // which in analogy to ``FEQ1'' for
-                                // the linear element is called
-                                // ``FEQ2''. The Lagrange elements of
-                                // polynomial degrees one through four
-                                // are all declared in this file.
-#include <fe/fe_lib.lagrange.h>
+                                // will import the declaration of
+                                // H1-conforming finite element shape
+                                // functions. This family of 
+                                // finite elements is called ``FE_Q''.
+#include <fe/fe_q.h>
                                 // We will not read the grid from a
                                 // file as in the previous example,
                                 // but generate it using a function
@@ -114,13 +111,14 @@ class LaplaceProblem
     void output_results (const unsigned int cycle) const;
 
     Triangulation<dim>   triangulation;
-    DoFHandler<dim>      dof_handler;
 
-                                    // In order to use the quadratic
-                                    // element, we only have to
-                                    // replace the declaration of the
-                                    // ``fe'' variable like this:
-    FEQ2<dim>            fe;
+                                    // We need a finite element
+                                    // again. This time, we will want
+                                    // to use quadratic polynomials
+                                    // (but this is only specified in
+                                    // the constructor):
+    FE_Q<dim>            fe;
+    DoFHandler<dim>      dof_handler;
 
                                     // This is the new variable in
                                     // the main class. We need an
@@ -187,8 +185,17 @@ void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
 };
 
 
+                                // This is mostly the same as before,
+                                // but this time we want to use the
+                                // quadratic element. To do so, we
+                                // only have to replace the
+                                // constructor argument (which was
+                                // ``1'' in all previous examples) by
+                                // the desired polynomial degree
+                                // (here ``2''):
 template <int dim>
 LaplaceProblem<dim>::LaplaceProblem () :
+                fe (2),
                dof_handler (triangulation)
 {};
 

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.