]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Errors according to Tensors and new LAC removed
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 26 Feb 1999 13:56:52 +0000 (13:56 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 26 Feb 1999 13:56:52 +0000 (13:56 +0000)
git-svn-id: https://svn.dealii.org/trunk@916 0785d39b-7218-0410-832d-ea1e28bc413d

14 files changed:
deal.II/base/include/base/tensor.h
deal.II/base/include/base/tensor_base.h
deal.II/base/include/base/tensor_function.h
deal.II/base/source/tensor.cc
deal.II/base/source/tensor_function.cc
deal.II/deal.II/include/fe/fe_values.h
deal.II/deal.II/include/grid/tria.h
deal.II/deal.II/include/numerics/vectors.h
deal.II/deal.II/source/fe/fe_system.cc
deal.II/deal.II/source/fe/fe_values.cc
deal.II/deal.II/source/grid/tria.cc
deal.II/deal.II/source/numerics/data_io.cc
deal.II/deal.II/source/numerics/vectors.cc
deal.II/doc/auto/scripts/index.pl

index a9ca5b8d2158dedb5e34cd351127d2ea4ac0ab38..93ae01a7579f7c4736a15948ae504b874768067e 100644 (file)
@@ -138,7 +138,7 @@ class Tensor {
                                      * usual in C++, the rightmost
                                      * index marches fastest.
                                      */
-    void unroll(vector<double> & result) const;
+    void unroll(Vector<double> & result) const;
     
 
                                     /**
@@ -146,18 +146,18 @@ class Tensor {
                                      */
     void clear ();
 
-                                    /**
-                                     * Exception
-                                     */
-    DeclException2(ExcWrongVectorSize, unsigned, int, << "Tensor has " << arg1
-                  << " entries, but vector has size " << arg2);
+//                                  /**
+//                                   * Exception
+//                                   */
+//     DeclException2(ExcWrongVectorSize, unsigned, int, << "Tensor has " << arg1
+//                << " entries, but vector has size " << arg2);
     
-                                    /**
-                                     *  Exception
-                                     */
-    DeclException1 (ExcInvalidIndex,
-                   int,
-                   << "Invalid index " << arg1);
+//                                  /**
+//                                   *  Exception
+//                                   */
+//     DeclException1 (ExcInvalidIndex,
+//                 int,
+//                 << "Invalid index " << arg1);
   private:
                                     /**
                                      * Array of tensors holding the
@@ -168,10 +168,10 @@ class Tensor {
                                     /**
                                      * Help function for unroll.
                                      */
-    void unroll_recursion(vector<double> & result, unsigned& start_index) const;
+    void unroll_recursion(Vector<double> & result, unsigned& start_index) const;
 
     template<>
-    friend void Tensor<rank_+1,dim>::unroll_recursion(vector<double> &, unsigned&) const;
+    friend void Tensor<rank_+1,dim>::unroll_recursion(Vector<double> &, unsigned&) const;
 };
 
 
index af991fbed6ad33dcd6332c8b2cdb17877ba4544f..49b6fec21e8c4ab3fbbafb4a2a59d58f3ca98916 100644 (file)
@@ -13,8 +13,7 @@
 #include <base/exceptions.h>
 #include <iostream>
 #include <vector>
-
-
+#include <lac/forward-declarations.h>
 
 
 // general template; specialized for rank==1; the general template is in
@@ -171,7 +170,7 @@ class Tensor<1,dim> {
                                      * usual in C++, the rightmost
                                      * index marches fastest.
                                      */
-    void unroll(vector<double> & result) const;
+    void unroll(Vector<double> & result) const;
      
   protected:
                                     /**
@@ -182,15 +181,17 @@ class Tensor<1,dim> {
                                     /**
                                      * Help function for unroll.
                                      */
-    void unroll_recursion(vector<double> & result, unsigned& start_index) const;
+    void unroll_recursion(Vector<double> & result, unsigned& start_index) const;
 
     template<>
-    friend void Tensor<2,dim>::unroll_recursion(vector<double> &, unsigned&) const;
+    friend void Tensor<2,dim>::unroll_recursion(Vector<double> &, unsigned&) const;
 };
 
-
-
-
+                                /**
+                                 * Exception
+                                 */
+DeclException2(ExcWrongVectorSize, unsigned, int, << "Tensor has " << arg1
+              << " entries, but vector has size " << arg2);
 DeclException1 (ExcDimTooSmall,
                int,
                << "Given dimensions must be >= 1, but was " << arg1);
index f79825e2f2f812549a7b6fa654650e7d901cafff..810a2d630009936ed684ee4a5f1f262d24cee259 100644 (file)
@@ -11,7 +11,7 @@
 #include <base/point.h>
 #include <base/functiontime.h>
 #include <base/tensorindex.h>
-
+#include <lac/forward-declarations.h>
 
 template<int rank_, int dim>
 class Tensor;
@@ -53,7 +53,7 @@ class VectorFunction :
                                      * the same size as the #n_components#
                                      * array.
                                      */
-    virtual void value (const Point<dim>  &p, vector<double> &values) const;
+    virtual void value (const Point<dim>  &p, Vector<double> &values) const;
 
                                     /**
                                      * Set #values# to the point values
@@ -64,7 +64,7 @@ class VectorFunction :
                                      * array.
                                      */
     virtual void value_list (const vector<Point<dim> > &points,
-                            vector<vector<double> > &values) const;
+                            vector<Vector<double> > &values) const;
 
                                     /**
                                      * Set #gradients# to the gradients of
@@ -172,11 +172,17 @@ class TensorFunction :
     virtual void gradient_list (const vector<Point<dim> > &points,
                                vector<Tensor<rank_+1,dim> > &gradients) const;
 
+                                    /**
+                                     * See #VectorFunction#.
+                                     */  
+    virtual void value (const Point<dim> &points,
+                            Vector<double> &values) const;
+
                                     /**
                                      * See #VectorFunction#.
                                      */  
     virtual void value_list (const vector<Point<dim> > &points,
-                            vector<vector<double> > &values) const;
+                            vector<Vector<double> > &values) const;
 
                                     /**
                                      * See #VectorFunction#.
index 3a4bac21a1db27650e01f93d817f1f3bbe44d3c7..d1035b500899ad3ee99a903a060e5952eb796725 100644 (file)
@@ -2,17 +2,24 @@
 
 #include <base/tensor.h>
 #include <cmath>
+#include <lac/vector.h>
 
 
-template <int dim> void
-Tensor<1,dim>::unroll( vector<double>& result) const
+template <int dim>
+void
+Tensor<1,dim>::unroll( Vector<double>& result) const
 {
-  Assert(false,         ExcNotImplemented());
+  Assert(result.size()==dim,
+        ExcWrongVectorSize(dim, result.size()));
+
+  unsigned index = 0;
+  unroll_recursion(result,index);
 }
 
 
-template <int rank_, int dim> void
-Tensor<rank_, dim>::unroll( vector<double>& result) const
+template <int rank_, int dim>
+void
+Tensor<rank_, dim>::unroll( Vector<double>& result) const
 {
   Assert(result.size()==pow(dim,rank_),
         ExcWrongVectorSize(pow(dim,rank_), result.size()));
@@ -22,8 +29,9 @@ Tensor<rank_, dim>::unroll( vector<double>& result) const
 }
 
 
-template <int rank_, int dim> void
-Tensor<rank_, dim>::unroll_recursion( vector<double>& result, unsigned& index) const
+template <int rank_, int dim>
+void
+Tensor<rank_, dim>::unroll_recursion( Vector<double>& result, unsigned& index) const
 {
   for (unsigned i=0; i<dim; ++i)
     {
@@ -33,13 +41,14 @@ Tensor<rank_, dim>::unroll_recursion( vector<double>& result, unsigned& index) c
 }
 
 
-template<int dim> void
-Tensor<1,dim>::unroll_recursion( vector<double>& result, unsigned& index) const
+template<int dim>
+void
+Tensor<1,dim>::unroll_recursion( Vector<double>& result, unsigned& index) const
 {
   for (unsigned i=0; i<dim; ++i)
     {
-      cerr << "[" << index << ',' << operator[](i) << ']' << endl;
-      result[index++] = operator[](i);
+//      cerr << "[" << index << ',' << operator[](i) << ']' << endl;
+      result(index++) = operator[](i);
     }
   
 }
index 4372c6ceb8a2a9fbc15602f144a18dd2bb845ccf..b2dc8127bf74b5c47d08413d81f576024c40a48b 100644 (file)
@@ -5,6 +5,7 @@
 #include <vector>
 #include <base/tensor.h>
 #include <cmath>
+#include <lac/vector.h>
 
 template <int dim>
 VectorFunction<dim>::VectorFunction(unsigned n_components, const double initial_time)
@@ -29,7 +30,7 @@ VectorFunction<dim>::operator () (const Point<dim> &, unsigned) const
 */
 
 template <int dim> void
-VectorFunction<dim>::value (const Point<dim>  &, vector<double> &) const
+VectorFunction<dim>::value (const Point<dim>  &, Vector<double> &) const
 {
   Assert (false, ExcPureFunctionCalled());
 }
@@ -37,7 +38,7 @@ VectorFunction<dim>::value (const Point<dim>  &, vector<double> &) const
 
 template <int dim> void
 VectorFunction<dim>::value_list (const vector<Point<dim> > &,
-                                vector<vector<double> > &) const
+                                vector<Vector<double> > &) const
 {
   Assert (false, ExcPureFunctionCalled());
 }
@@ -139,7 +140,7 @@ TensorFunction<rank_, dim>::value (const Point<dim>  &p, vector<double> &erg) co
 
 template <int rank_, int dim> void
 TensorFunction<rank_, dim>::value_list (const vector<Point<dim> > & points,
-                                vector<vector<double> > & values) const
+                                vector<Vector<double> > & values) const
 {
   Assert (values.size() == points.size(),
          ExcVectorHasWrongSize(values.size(), points.size()));
index b57150d59f7c58efaf2af882f5323b1cb4652cd4..d55dc026c285a07747aff69198fae6babc751a38 100644 (file)
@@ -304,7 +304,7 @@ class FEValuesBase {
                                    * elements.
                                    */
     void get_function_values (const Vector<double>    &fe_function,
-                             vector<vector<double> > &values) const;
+                             vector<Vector<double> > &values) const;
 
                                     /**
                                      * Return the gradient of the #i#th shape
index e3d6de9137d68df74fe7ae2de5613fd0dd531f90..cddb2210c29891b54e396f6d6e037fd8df79a449 100644 (file)
@@ -1676,7 +1676,8 @@ class Triangulation : public TriaDimensionInfo<dim>, public Subscriptor {
                                      * much needed here and saving memory may
                                      * be a good goal when using many cells.
                                      */
-    void refine (const Vector<float> &criteria,
+    template <typename number>
+    void refine (const Vector<number> &criteria,
                 const double         threshold);
 
                                     /**
@@ -1692,7 +1693,8 @@ class Triangulation : public TriaDimensionInfo<dim>, public Subscriptor {
                                      * much needed here and saving memory may
                                      * be a good goal when using many cells.
                                      */
-    void coarsen (const Vector<float> &criteria,
+    template <typename number>
+    void coarsen (const Vector<number> &criteria,
                  const double         threshold);
     
                                     /**
@@ -1716,7 +1718,8 @@ class Triangulation : public TriaDimensionInfo<dim>, public Subscriptor {
                                      * much needed here and saving memory may
                                      * be a good goal when using many cells.
                                      */
-    void refine_and_coarsen_fixed_number (const Vector<float> &criteria,
+    template <typename number>
+    void refine_and_coarsen_fixed_number (const Vector<number> &criteria,
                                          const double         top_fraction_of_cells,
                                          const double         bottom_fraction_of_cells);
 
@@ -1741,7 +1744,8 @@ class Triangulation : public TriaDimensionInfo<dim>, public Subscriptor {
                                      * much needed here and saving memory may
                                      * be a good goal when using many cells.
                                      */
-    void refine_and_coarsen_fixed_fraction (const Vector<float> &criteria,
+    template<typename number>
+    void refine_and_coarsen_fixed_fraction (const Vector<number> &criteria,
                                            const double         top_fraction,
                                            const double         bottom_fraction);
     
index a4567f8a75f0b021e30f3797f72801c45371bfa9..f3b9c6e57c37d8f30ff7aadedd90b8c5162cefa2 100644 (file)
@@ -398,14 +398,14 @@ class VectorTools {
                                      * Compute the error for the solution of a system.
                                      * See the other #integrate_difference#.
                                      */
-//     static void integrate_difference (const DoFHandler<dim>    &dof,
-//                                   const Vector<double>     &fe_function,
-//                                   const TensorFunction<1,dim>&exact_solution,
-//                                   Vector<float>            &difference,
-//                                   const Quadrature<dim>    &q,
-//                                   const FiniteElement<dim> &fe,
-//                                   const NormType           &norm,
-//                                   const Boundary<dim>      &boundary);
+     static void integrate_difference (const DoFHandler<dim>   &dof,
+                                     const Vector<double>     &fe_function,
+                                     const VectorFunction<dim>&exact_solution,
+                                     Vector<float>            &difference,
+                                     const Quadrature<dim>    &q,
+                                     const FiniteElement<dim> &fe,
+                                     const NormType           &norm,
+                                     const Boundary<dim>      &boundary);
     
 
                                     /**
index b9ae5a99b249a79c31f34e0d98bf49f570c1eb8e..4415ea60f76c05d60e046cb255d09dcd3ad6a0a4 100644 (file)
@@ -36,7 +36,7 @@ void FESystem<dim>::initialize ()
   
                                   // 1. Vertices
       total_index = 0;
-  for (unsigned vertex_number= 0 ; vertex_number < GeometryInfo<2>::vertices_per_cell ;
+  for (unsigned vertex_number= 0 ; vertex_number < GeometryInfo<dim>::vertices_per_cell ;
        ++vertex_number)
     {
       unsigned comp_start = 0;
@@ -60,7 +60,7 @@ void FESystem<dim>::initialize ()
     }
   
                                   // 2. Lines
-  for (unsigned line_number= 0 ; line_number < GeometryInfo<2>::lines_per_cell ;
+  for (unsigned line_number= 0 ; line_number < GeometryInfo<dim>::lines_per_cell ;
        ++line_number)
     {
       unsigned comp_start = 0;
@@ -85,7 +85,7 @@ void FESystem<dim>::initialize ()
     }
   
                                   // 3. Quads
-  for (unsigned quad_number= 0 ; quad_number < GeometryInfo<2>::quads_per_cell ;
+  for (unsigned quad_number= 0 ; quad_number < GeometryInfo<dim>::quads_per_cell ;
        ++quad_number)
     {
       unsigned comp_start = 0;
@@ -110,7 +110,7 @@ void FESystem<dim>::initialize ()
     }
   
                                   // 4. Hex
-  for (unsigned hex_number= 0 ; hex_number < GeometryInfo<2>::hexes_per_cell ;
+  for (unsigned hex_number= 0 ; hex_number < GeometryInfo<dim>::hexes_per_cell ;
        ++hex_number)
     {
       unsigned comp_start = 0;
index 67c57e16fa9c4fa13dff6fe11a83ad68f8728b73..31568f053c17fa671ce07cb1abab71bc98fd26b8 100644 (file)
@@ -90,7 +90,7 @@ void FEValuesBase<dim>::get_function_values (const Vector<double> &fe_function,
 
 template <int dim>
 void FEValuesBase<dim>::get_function_values (const Vector<double>     &fe_function,
-                                            vector< vector<double> > &values) const
+                                            vector< Vector<double> > &values) const
 {
   Assert (fe->n_components == values.size(),
          ExcWrongNoOfComponents());
@@ -113,7 +113,7 @@ void FEValuesBase<dim>::get_function_values (const Vector<double>     &fe_functi
                                   // functions
   for (unsigned int point=0; point<n_quadrature_points; ++point)
     for (unsigned int shape_func=0; shape_func<total_dofs; ++shape_func)
-      values[fe->system_to_component_index(shape_func).first][point]
+      values[point](fe->system_to_component_index(shape_func).first)
        += (dof_values(shape_func) * shape_values[selected_dataset](shape_func, point));
 };
 
index 605ac44c11390070c116f7fa042d3d94b066bd15..85108aead8148de549fbf4c63334438583541261 100644 (file)
@@ -3761,9 +3761,10 @@ void Triangulation<3>::print_gnuplot (ostream &out,
 
 
 
-template <int dim>
-void Triangulation<dim>::refine (const Vector<float> &criteria,
-                                const double         threshold) {
+template <int dim> template <typename number>
+void Triangulation<dim>::refine (const Vector<number> &criteria,
+                                const double         threshold)
+{
   Assert (criteria.size() == n_active_cells(),
          ExcInvalidVectorSize(criteria.size(), n_active_cells()));
 
@@ -3785,9 +3786,10 @@ void Triangulation<dim>::refine (const Vector<float> &criteria,
 
 
 
-template <int dim>
-void Triangulation<dim>::coarsen (const Vector<float> &criteria,
-                                 const double         threshold) {
+template <int dim> template <typename number>
+void Triangulation<dim>::coarsen (const Vector<number> &criteria,
+                                 const double         threshold)
+{
   Assert (criteria.size() == n_active_cells(),
          ExcInvalidVectorSize(criteria.size(), n_active_cells()));
 
@@ -3801,10 +3803,12 @@ void Triangulation<dim>::coarsen (const Vector<float> &criteria,
 
 
 
-template <int dim>
-void Triangulation<dim>::refine_and_coarsen_fixed_number (const Vector<float> &criteria,
-                                                         const double         top_fraction,
-                                                         const double         bottom_fraction) {
+template <int dim> template <typename number>
+void
+Triangulation<dim>::refine_and_coarsen_fixed_number (const Vector<number> &criteria,
+                                                    const double         top_fraction,
+                                                    const double         bottom_fraction)
+{
                                   // correct number of cells is
                                   // checked in #refine#
   Assert ((top_fraction>=0) && (top_fraction<=1), ExcInvalidParameterValue());
@@ -3817,7 +3821,7 @@ void Triangulation<dim>::refine_and_coarsen_fixed_number (const Vector<float> &c
   const int coarsen_cells = max(static_cast<int>(bottom_fraction*criteria.size()),
                                1);
   
-  Vector<float> tmp(criteria);  
+  Vector<number> tmp(criteria);  
   nth_element (tmp.begin(), tmp.begin()+refine_cells,
               tmp.end(),
               greater<double>());
@@ -3839,9 +3843,9 @@ double sqr(double a) {
 
 
 
-template <int dim>
+template <int dim> template <typename number>
 void
-Triangulation<dim>::refine_and_coarsen_fixed_fraction (const Vector<float> &criteria,
+Triangulation<dim>::refine_and_coarsen_fixed_fraction (const Vector<number> &criteria,
                                                       const double         top_fraction,
                                                       const double         bottom_fraction) {
                                   // correct number of cells is
@@ -3854,7 +3858,7 @@ Triangulation<dim>::refine_and_coarsen_fixed_fraction (const Vector<float> &crit
                                   // error, which is what we have to sum
                                   // up and compare with
                                   // #fraction_of_error*total_error#.
-  Vector<float> tmp(criteria);
+  Vector<number> tmp(criteria);
   const double total_error = tmp.l1_norm();
 
   Vector<float> partial_sums(criteria.size());
@@ -7048,3 +7052,36 @@ void Triangulation<dim>::read_bool_vector (const unsigned int  magic_number1,
 // explicit instantiations
 template class Triangulation<deal_II_dimension>;
 
+template void Triangulation<deal_II_dimension>
+::refine (const Vector<float> &, const double);
+
+template void Triangulation<deal_II_dimension>
+::refine (const Vector<double> &, const double);
+
+template void Triangulation<deal_II_dimension>
+::coarsen (const Vector<float> &, const double);
+
+template void Triangulation<deal_II_dimension>
+::coarsen (const Vector<double> &, const double);
+
+
+template void Triangulation<deal_II_dimension>
+::refine_and_coarsen_fixed_number (const Vector<double> &,
+                                  const double         top_fraction,
+                                  const double         bottom_fraction);
+
+template void Triangulation<deal_II_dimension>
+::refine_and_coarsen_fixed_number (const Vector<float> &criteria,
+                                  const double         top_fraction,
+                                  const double         bottom_fraction);
+
+template void Triangulation<deal_II_dimension>
+::refine_and_coarsen_fixed_fraction (const Vector<double> &criteria,
+                                    const double         top_fraction,
+                                    const double         bottom_fraction);
+
+template void Triangulation<deal_II_dimension>
+::refine_and_coarsen_fixed_fraction (const Vector<float> &criteria,
+                                    const double         top_fraction,
+                                    const double         bottom_fraction);
+
index da9a132577044c1aa7ecf02a93f3adb4e4cf1898..b34ce01796ba7707a782afccaa40cb869651a6b5 100644 (file)
@@ -543,8 +543,11 @@ void DataOut<dim>::write_gnuplot (ostream &out, unsigned int accuracy) const
   
   FEValues<dim> fe(dofs->get_fe(), points, UpdateFlags(update_q_points));
   const StraightBoundary<dim> boundary;
-  vector< vector <vector<double> > >
-    values (data.size(), vector< vector<double> >(dofs->get_fe().n_components, vector<double>(points.n_quadrature_points)));
+  vector< vector <Vector<double> > >
+    values (data.size(),
+           vector< Vector<double> >(points.n_quadrature_points,
+                                    Vector<double>(dofs->get_fe().n_components
+                                    )));
 
   for (cell=dofs->begin_active(); cell!=endc; ++cell) 
     {
@@ -570,7 +573,7 @@ void DataOut<dim>::write_gnuplot (ostream &out, unsigned int accuracy) const
                    out << pt << "  ";
                    for (unsigned int i=0; i!=data.size(); ++i)
                      for (unsigned int j=0; j < dofs->get_fe().n_components; ++j)
-                       out << values[i][j][supp_pt]
+                       out << values[i][supp_pt](j)
                            << ' ';
                    out << endl;
                  };
@@ -592,7 +595,7 @@ void DataOut<dim>::write_gnuplot (ostream &out, unsigned int accuracy) const
                    
                    for (unsigned int i=0; i!=data.size(); ++i)
                      for (unsigned int j=0; j < dofs->get_fe().n_components; ++j)
-                       out << values[i][j][supp_pt]
+                       out << values[i][supp_pt](j)
                            << ' ';
                    out << endl;
                  }
@@ -625,7 +628,7 @@ void DataOut<dim>::write_gnuplot (ostream &out, unsigned int accuracy) const
                            
                            for (unsigned int i=0; i!=data.size(); ++i)
                              for (unsigned int j=0; j < dofs->get_fe().n_components; ++j)
-                               out << values[i][j][supp_pt]
+                               out << values[i][supp_pt](j)
                                    << ' ';
                            out << endl;
                          }
index 31330e7b628d7d2ebba12eda1c1e1b5f7a94bc38..68785d54dd4060b2c1c1d3f4c23c789d696bfa28 100644 (file)
@@ -4,6 +4,7 @@
 
 
 #include <base/function.h>
+#include <base/tensorfunction.h>
 #include <grid/dof.h>
 #include <grid/dof_accessor.h>
 #include <grid/tria_iterator.h>
@@ -557,193 +558,196 @@ void VectorTools<dim>::integrate_difference (const DoFHandler<dim>    &dof,
 
 
 
-// template <int dim>
-// void VectorTools<dim>::integrate_difference (const DoFHandler<dim>    &dof,
-//                                          const Vector<double>     &fe_function,
-//                                          const TensorFunction<1,dim>&exact_solution,
-//                                          Vector<float>            &difference,
-//                                          const Quadrature<dim>    &q,
-//                                          const FiniteElement<dim> &fe,
-//                                          const NormType           &norm,
-//                                          const Boundary<dim>      &boundary) {
-//   Assert (fe == dof.get_fe(), ExcInvalidFE());
+template <int dim>
+void
+VectorTools<dim>::integrate_difference (const DoFHandler<dim>   &dof,
+                                       const Vector<double>     &fe_function,
+                                       const VectorFunction<dim>&exact_solution,
+                                       Vector<float>            &difference,
+                                       const Quadrature<dim>    &q,
+                                       const FiniteElement<dim> &fe,
+                                       const NormType           &norm,
+                                       const Boundary<dim>      &boundary)
+{
+   Assert (fe == dof.get_fe(), ExcInvalidFE());
   
-//   difference.reinit (dof.get_tria().n_active_cells());
+   difference.reinit (dof.get_tria().n_active_cells());
   
-//   UpdateFlags update_flags = UpdateFlags (update_q_points  |
-//                                       update_JxW_values);
-//   if ((norm==H1_seminorm) || (norm==H1_norm))
-//     update_flags = UpdateFlags (update_flags | update_gradients);
-//   FEValues<dim> fe_values(fe, q, update_flags);
+   UpdateFlags update_flags = UpdateFlags (update_q_points  |
+                                         update_JxW_values);
+   if ((norm==H1_seminorm) || (norm==H1_norm))
+     update_flags = UpdateFlags (update_flags | update_gradients);
+   FEValues<dim> fe_values(fe, q, update_flags);
   
-//                                // loop over all cells
-//   DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
-//                                     endc = dof.end();
-//   for (unsigned int index=0; cell != endc; ++cell, ++index)
-//     {
-//       double diff=0;
-//                                    // initialize for this cell
-//       fe_values.reinit (cell, boundary);
-
-//       switch (norm) 
-//     {
-//       case mean:
-//       case L1_norm:
-//       case L2_norm:
-//       case Linfty_norm:
-//       case H1_norm:
-//       {
-//                                          // we need the finite element
-//                                          // function \psi at the different
-//                                          // integration points. Compute
-//                                          // it like this:
-//                                          // \psi(x_j)=\sum_i v_i \phi_i(x_j)
-//                                          // with v_i the nodal values of the
-//                                          // fe_function and \phi_i(x_j) the
-//                                          // matrix of the trial function
-//                                          // values at the integration point
-//                                          // x_j. Then the vector
-//                                          // of the \psi(x_j) is v*Phi with
-//                                          // v being the vector of nodal
-//                                          // values on this cell and Phi
-//                                          // the matrix.
-//                                          //
-//                                          // we then need the difference:
-//                                          // reference_function(x_j)-\psi_j
-//                                          // and assign that to the vector
-//                                          // \psi.
-//         const unsigned int n_q_points = q.n_quadrature_points;
-//         vector<double>   psi (n_q_points);
-
-//                                          // in praxi: first compute
-//                                          // exact fe_function vector
-//         exact_solution.value_list (fe_values.get_quadrature_points(),
-//                                    psi);
-//                                          // then subtract finite element
-//                                          // fe_function
-//         if (true) 
-//           {
-//             vector< vector<double> > function_values (fe.n_components, vector<double>(n_q_points, 0));
-//             fe_values.get_function_values (fe_function, function_values);
-
-//             transform (psi.begin(), psi.end(),
-//                        function_values.begin(),
-//                        psi.begin(),
-//                        minus<double>());
-//           };            
-
-//                                          // for L1_norm and Linfty_norm:
-//                                          // take absolute
-//                                          // value, for the L2_norm take
-//                                          // square of psi
-//         switch (norm) 
-//           {
-//             case mean:
-//                   break;
-//             case L1_norm:
-//             case Linfty_norm:
-//                   transform (psi.begin(), psi.end(),
-//                              psi.begin(), ptr_fun(fabs));
-//                   break;
-//             case L2_norm:
-//             case H1_norm:
-//                   transform (psi.begin(), psi.end(),
-//                              psi.begin(), ptr_fun(sqr));
-//                   break;
-//             default:
-//                   Assert (false, ExcNotImplemented());
-//           };
-
-//                                          // ok, now we have the integrand,
-//                                          // let's compute the integral,
-//                                          // which is
-//                                          // sum_j psi_j JxW_j
-//                                          // (or |psi_j| or |psi_j|^2
-//         switch (norm) 
-//           {
-//             case mean:
-//             case L1_norm:
-//                   diff = inner_product (psi.begin(), psi.end(),
-//                                         fe_values.get_JxW_values().begin(),
-//                                         0.0);
-//                   break;
-//             case L2_norm:
-//             case H1_norm:
-//                   diff = sqrt(inner_product (psi.begin(), psi.end(),
-//                                              fe_values.get_JxW_values().begin(),
-//                                              0.0));
-//                   break;
-//             case Linfty_norm:
-//                   diff = *max_element (psi.begin(), psi.end());
-//                   break;
-//             default:
-//                   Assert (false, ExcNotImplemented());
-//           };
-
-//                                          // note: the H1_norm uses the result
-//                                          // of the L2_norm and control goes
-//                                          // over to the next case statement!
-//         if (norm != H1_norm)
-//           break;
-//       };
-
-//       case H1_seminorm:
-//       {
-//                                          // note: the computation of the
-//                                          // H1_norm starts at the previous
-//                                          // case statement, but continues
-//                                          // here!
-
-//                                          // for H1_norm: re-square L2_norm.
-//         diff = sqr(diff);
-
-//                                          // same procedure as above, but now
-//                                          // psi is a vector of gradients
-//         const unsigned int n_q_points = q.n_quadrature_points;
-//         vector<Tensor<1,dim> >   psi (n_q_points);
-
-//                                          // in praxi: first compute
-//                                          // exact fe_function vector
-//         exact_solution.gradient_list (fe_values.get_quadrature_points(),
-//                                       psi);
-           
-//                                          // then subtract finite element
-//                                          // fe_function
-//         if (true) 
-//           {
-//             vector<Tensor<1,dim> > function_grads (n_q_points, Tensor<1,dim>());
-//             fe_values.get_function_grads (fe_function, function_grads);
-
-//             transform (psi.begin(), psi.end(),
-//                        function_grads.begin(),
-//                        psi.begin(),
-//                        minus<Tensor<1,dim> >());
-//           };
-//                                          // take square of integrand
-//         vector<double> psi_square (psi.size(), 0.0);
-//         for (unsigned int i=0; i<n_q_points; ++i)
-//           psi_square[i] = sqr_point(psi[i]);
-
-//                                          // add seminorm to L_2 norm or
-//                                          // to zero
-//         diff += inner_product (psi_square.begin(), psi_square.end(),
-//                                fe_values.get_JxW_values().begin(),
-//                                0.0);
-//         diff = sqrt(diff);
-
-//         break;
-//       };
+                                  // loop over all cells
+   DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
+                                       endc = dof.end();
+   for (unsigned int index=0; cell != endc; ++cell, ++index)
+     {
+       double diff=0;
+                                      // initialize for this cell
+       fe_values.reinit (cell, boundary);
+
+       switch (norm) 
+       {
+         case mean:
+         case L1_norm:
+         case L2_norm:
+         case Linfty_norm:
+         case H1_norm:
+         {
+                                            // we need the finite element
+                                            // function \psi at the different
+                                            // integration points. Compute
+                                            // it like this:
+                                            // \psi(x_j)=\sum_i v_i \phi_i(x_j)
+                                            // with v_i the nodal values of the
+                                            // fe_function and \phi_i(x_j) the
+                                            // matrix of the trial function
+                                            // values at the integration point
+                                            // x_j. Then the vector
+                                            // of the \psi(x_j) is v*Phi with
+                                            // v being the vector of nodal
+                                            // values on this cell and Phi
+                                            // the matrix.
+                                            //
+                                            // we then need the difference:
+                                            // reference_function(x_j)-\psi_j
+                                            // and assign that to the vector
+                                            // \psi.
+           const unsigned int n_q_points = q.n_quadrature_points;
+           vector<Vector<double> >  psi (n_q_points);
+
+                                            // in praxi: first compute
+                                            // exact fe_function vector
+           exact_solution.value_list (fe_values.get_quadrature_points(),
+                                      psi);
+                                            // then subtract finite element
+                                            // fe_function
+           if (true) 
+             {
+               vector< Vector<double> > function_values (n_q_points,
+                                                         Vector<double>(fe.n_components));
+               fe_values.get_function_values (fe_function, function_values);
+
+/*             transform (psi.begin(), psi.end(),
+                          function_values.begin(),
+                          psi.begin(),
+                          minus<double>());
+*/           };            
+
+                                            // for L1_norm and Linfty_norm:
+                                            // take absolute
+                                            // value, for the L2_norm take
+                                            // square of psi
+/*         switch (norm) 
+             {
+               case mean:
+                     break;
+               case L1_norm:
+               case Linfty_norm:
+                     transform (psi.begin(), psi.end(),
+                                psi.begin(), ptr_fun(fabs));
+                     break;
+               case L2_norm:
+               case H1_norm:
+                     transform (psi.begin(), psi.end(),
+                                psi.begin(), ptr_fun(sqr));
+                     break;
+               default:
+                     Assert (false, ExcNotImplemented());
+             };
+*/
+                                            // ok, now we have the integrand,
+                                            // let's compute the integral,
+                                            // which is
+                                            // sum_j psi_j JxW_j
+                                            // (or |psi_j| or |psi_j|^2
+/*         switch (norm) 
+             {
+               case mean:
+               case L1_norm:
+                     diff = inner_product (psi.begin(), psi.end(),
+                                           fe_values.get_JxW_values().begin(),
+                                           0.0);
+                     break;
+               case L2_norm:
+               case H1_norm:
+                     diff = sqrt(inner_product (psi.begin(), psi.end(),
+                                                fe_values.get_JxW_values().begin(),
+                                                0.0));
+                     break;
+               case Linfty_norm:
+                     diff = *max_element (psi.begin(), psi.end());
+                     break;
+               default:
+                     Assert (false, ExcNotImplemented());
+             };
+*/
+                                            // note: the H1_norm uses the result
+                                            // of the L2_norm and control goes
+                                            // over to the next case statement!
+           if (norm != H1_norm)
+             break;
+         };
+
+         case H1_seminorm:
+         {
+                                            // note: the computation of the
+                                            // H1_norm starts at the previous
+                                            // case statement, but continues
+                                            // here!
+
+                                            // for H1_norm: re-square L2_norm.
+           diff = sqr(diff);
+
+                                            // same procedure as above, but now
+                                            // psi is a vector of gradients
+           const unsigned int n_q_points = q.n_quadrature_points;
+           vector<Tensor<1,dim> >   psi (n_q_points);
+
+                                            // in praxi: first compute
+                                            // exact fe_function vector
+/*         exact_solution.gradient_list (fe_values.get_quadrature_points(),
+                                         psi);
+*/         
+                                            // then subtract finite element
+                                            // fe_function
+           if (true) 
+             {
+               vector<Tensor<1,dim> > function_grads (n_q_points, Tensor<1,dim>());
+               fe_values.get_function_grads (fe_function, function_grads);
+
+/*             transform (psi.begin(), psi.end(),
+                          function_grads.begin(),
+                          psi.begin(),
+                          minus<Tensor<1,dim> >());
+*/           };
+                                            // take square of integrand
+           vector<double> psi_square (psi.size(), 0.0);
+           for (unsigned int i=0; i<n_q_points; ++i)
+             psi_square[i] = sqr_point(psi[i]);
+
+                                            // add seminorm to L_2 norm or
+                                            // to zero
+/*         diff += inner_product (psi_square.begin(), psi_square.end(),
+                                  fe_values.get_JxW_values().begin(),
+                                  0.0);
+           diff = sqrt(diff);
+*/
+           break;
+         };
                                             
-//       default:
-//             Assert (false, ExcNotImplemented());
-//     };
+         default:
+               Assert (false, ExcNotImplemented());
+       };
 
       
-//                                    // append result of this cell
-//                                    // to the end of the vector
-//       difference(index) = diff;
-//     };
-// };
+                                      // append result of this cell
+                                      // to the end of the vector
+       difference(index) = diff;
+     };
+ };
 
 
 template VectorTools<deal_II_dimension>;
index 737c09121792c885a861c65bf1b3abe6890151dc..85a4086f55b82454965406d82ebfc164ee3e30e0 100644 (file)
@@ -4,7 +4,9 @@ while(<>)
 {
     next if (m/DeclException/);
     chop;
-    if (m/(.+)=(.+)/)
+    s/\</\&lt\;/g;
+    s/\>/\&gt\;/g;
+    if (m/(.+)=([^=]+)=?$/)
     {
        @n = split "::", $1;
 
@@ -13,7 +15,7 @@ while(<>)
        {
            $entry .= "+" . pop @n;
        }
-       $entry .= '=' . $2 . '=' . $library;
+       $entry .= ',' . $2 . ',' . $library;
        push @entries, $entry;
     }
     else
@@ -34,7 +36,7 @@ EOT
     ;
 foreach $entry (sort @entries)
 {
-    @l = split "=", $entry;
+    @l = split ",", $entry;
     @c = split '\+', @l[0];
     @f = split "#", @l[1];
 

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.