]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Preparation for systems and mixed discretizations
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 24 Feb 1999 13:19:48 +0000 (13:19 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 24 Feb 1999 13:19:48 +0000 (13:19 +0000)
git-svn-id: https://svn.dealii.org/trunk@893 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/Makefile
deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/fe/fe_system.h
deal.II/deal.II/include/numerics/vectors.h
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_lib.criss_cross.cc
deal.II/deal.II/source/fe/fe_linear_mapping.cc
deal.II/deal.II/source/fe/fe_system.cc
deal.II/deal.II/source/fe/fe_values.cc
deal.II/deal.II/source/numerics/data_io.cc
deal.II/deal.II/source/numerics/vectors.cc

index e277543fe0489bdc7f56a5dd5e5e7b2c9038ac9f..fa5822401695ca1c6e6bf111f069e4825a33b31b 100644 (file)
@@ -23,7 +23,7 @@ examples-clean:
 source-clean:
        cd source ; make clean
 TAGS:
-       etags include/*/* source/*/*
+       etags include/*/*.h source/*/*.cc
 
 
 .PHONY: all deal.II examples TAGS
index 7d29d970579f18ddac4061d1c725c67536bc91b1..b6d47ef7f4d89e7254ae20f4cdbe83e682e6c217 100644 (file)
@@ -92,8 +92,9 @@ class FiniteElementData
                                      */
     const unsigned int n_transform_functions;
 
+    
                                     /**
-                                     * Number of components.
+                                     * Number of components and dimension of the image space.
                                      */
     const unsigned int n_components;
 
@@ -853,7 +854,12 @@ class FiniteElement : public FiniteElementBase<dim> {
                                      *
                                      * The function assumes that the fields
                                      * already have the right number of
-                                     * elements.
+                                     * elements. It has to be
+                                     * guaranteed, that fields that are
+                                     * not requested for update are not changed.
+                                     * This also means, that these
+                                     * fields have to be filled with
+                                     * the correct values beforehand.
                                      *
                                      * This function is more or less an
                                      * interface to the #FEValues# class and
index d7026356834e3b7d0cdb83717d646abb2e8f1188..00c152ec23fb0d1ceb744750ed9d956e3649464f 100644 (file)
@@ -316,6 +316,23 @@ class FESystem : public FiniteElement<dim>
                                     const vector<Point<dim-1> > &unit_points,
                                     vector<Point<dim> >         &normal_vectors) const;
 
+                                    /**
+                                     * Implementation of the corresponding function of #FiniteElement#.
+                                     */
+    virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
+                                const vector<Point<dim> >            &unit_points,
+                                vector<Tensor<2,dim> >               &jacobians,
+                                const bool              compute_jacobians,
+                                vector<Tensor<3,dim> > &jacobians_grad,
+                                const bool              compute_jacobians_grad,
+                                vector<Point<dim> > &support_points,
+                                const bool           compute_support_points,
+                                vector<Point<dim> > &q_points,
+                                const bool           compute_q_points,
+                                const dFMatrix      &shape_values_transform,
+                                const vector<vector<Tensor<1,dim> > > &shape_grad_transform,
+                                const Boundary<dim> &boundary) const;
+    
                                     /** 
                                      * Number of different base
                                      * elements of this object.
@@ -425,6 +442,10 @@ class FESystem : public FiniteElement<dim>
                                      * the #.h# file.
                                      */
     void initialize();
+                                    /**
+                                     *Exception.
+                                     */
+    DeclException0(ExcElementTransformNotEqual);
 };
 
 
@@ -478,6 +499,9 @@ FESystem<dim>::FESystem (const FE1 &fe1, const unsigned int n1,
                base_elements(2),
                component_to_base_table(n_components)
 {
+  Assert(fe1.n_transform_functions == fe2.n_transform_functions,
+        ExcElementTransformNotEqual());
+  
   base_elements[0] = ElementPair(new FE1, n1);
   base_elements[0].first -> subscribe ();
   base_elements[1] = ElementPair(new FE2, n2);
index 8a32f053c3e85533a3f0133d2ae0fe73b5e5f3dc..90720acd888fdc2b20b46d64b32f7e84d9b92e3e 100644 (file)
@@ -19,6 +19,7 @@ template <int dim> class Boundary;
 template <int dim> class StraightBoundary;
 class ConstraintMatrix;
 class dVector;
+class VectorFunction;
 
 
 /**
@@ -376,6 +377,7 @@ class VectorTools {
                                         map<int,double>          &boundary_values);
     
                                     /**
+                                     * Compute the error of the finite element solution.
                                      * Integrate the difference between
                                      * a finite element function and
                                      * the reference function, which
@@ -393,6 +395,20 @@ class VectorTools {
                                      const NormType           &norm,
                                      const Boundary<dim> &boundary=StraightBoundary<dim>());
 
+                                    /**
+                                     * Compute the error for the solution of a system.
+                                     * See the other #integrate_difference#.
+                                     */
+    static void integrate_difference (const DoFHandler<dim>    &dof,
+                                     const dVector            &fe_function,
+                                     const VectorFunction     &exact_solution,
+                                     dVector                  &difference,
+                                     const Quadrature<dim>    &q,
+                                     const FiniteElement<dim> &fe,
+                                     const NormType           &norm,
+                                     const Boundary<dim>      &boundary);
+    
+
                                     /**
                                      * Exception
                                      */
index bf65e97833cff0220c2723abe0daf34e07960c46..32c438a6155cf125f8243eb7297ab7c866a199b1 100644 (file)
@@ -99,7 +99,6 @@ FiniteElementData<dim>::FiniteElementData (const unsigned int dofs_per_vertex,
 };
 
 
-
 template<int dim>
 bool FiniteElementData<dim>::operator== (const FiniteElementData<dim> &f) const
 {
@@ -241,11 +240,13 @@ void FiniteElement<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell,
                                       const dFMatrix       &,
                                       const vector<vector<Tensor<1,1> > > &,
                                       const Boundary<1> &boundary) const {
-  Assert (jacobians.size() == unit_points.size(),
+  Assert ((!compute_jacobians) || (jacobians.size() == unit_points.size()),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
-  Assert (q_points.size() == unit_points.size(),
+  Assert ((!compute_jacobians_grad) || (jacobians_grad.size() == unit_points.size()),
+         ExcWrongFieldDimension(jacobians_grad.size(), unit_points.size()));
+  Assert ((!compute_q_points) || (q_points.size() == unit_points.size()),
          ExcWrongFieldDimension(q_points.size(), unit_points.size()));
-  Assert (support_points.size() == total_dofs,
+  Assert ((!compute_support_points) || (support_points.size() == total_dofs),
          ExcWrongFieldDimension(support_points.size(), total_dofs));
 
 
index 40f21263bdc9a4f4395ac02ec58708157db92077..6000e8c3c3f248dc24ea34513d2ac4278fba5851 100644 (file)
@@ -939,7 +939,7 @@ void FECrissCross<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &ce
                                        vector<Tensor<3,dim> > &jacobians_grad,
                                        const bool              compute_jacobians_grad,
                                        vector<Point<dim> >    &support_points,
-                                       const bool,
+                                       const bool compute_support_points,
                                        vector<Point<dim> >    &q_points,
                                        const bool              compute_q_points,
                                        const dFMatrix         &shape_values_transform,
@@ -957,7 +957,8 @@ void FECrissCross<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &ce
 
                                   // we need the support points in any
                                   // way, wanted or not by the user
-  get_support_points (cell, boundary, support_points);
+  if (compute_support_points)
+    get_support_points (cell, boundary, support_points);
 
   if (compute_q_points) 
     {
index 0d0acd43979a11f24ad08a9d59a74276f532f199..bef427b76a3f667972f7cbb4ec880212ff2f89ca 100644 (file)
@@ -567,12 +567,15 @@ void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator
                                           const bool           compute_q_points,
                                           const dFMatrix      &shape_values_transform,
                                           const vector<vector<Tensor<1,dim> > > &/*shape_grad_transform*/,
-                                          const Boundary<dim> &boundary) const {
-  Assert (jacobians.size() == unit_points.size(),
+                                          const Boundary<dim> &boundary) const
+{
+  Assert ((!compute_jacobians) || (jacobians.size() == unit_points.size()),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
-  Assert (q_points.size() == unit_points.size(),
+  Assert ((!compute_jacobians_grad) || (jacobians_grad.size() == unit_points.size()),
+         ExcWrongFieldDimension(jacobians_grad.size(), unit_points.size()));
+  Assert ((!compute_q_points) || (q_points.size() == unit_points.size()),
          ExcWrongFieldDimension(q_points.size(), unit_points.size()));
-  Assert (support_points.size() == total_dofs,
+  Assert ((!compute_support_points) || (support_points.size() == total_dofs),
          ExcWrongFieldDimension(support_points.size(), total_dofs));
 
   
index 94991e57d4a6e1aae54c73576c7febc24625f68b..6ea9c6f79a452683f433ab61e417dbe07ed7b07f 100644 (file)
@@ -1,5 +1,5 @@
 /* $Id$ */
-/* Copyright W. Bangerth, University of Heidelberg, 1990 */
+/* Copyright W. Bangerth, G. Kanschat University of Heidelberg, 1990 */
 
 
 #include <fe/fe_system.h>
@@ -196,7 +196,7 @@ FESystem<1>::multiply_dof_numbers (const FiniteElementData<1> &fe1,
 {
   return FiniteElementData<1> (fe1.dofs_per_vertex * N1 + fe2.dofs_per_vertex * N2 ,
                               fe1.dofs_per_line * N1 + fe2.dofs_per_line * N2 ,
-                              fe1.n_transform_functions + fe2.n_transform_functions ,
+                              fe1.n_transform_functions,
                               fe1.n_components * N1 + fe2.n_components * N2 );
 };
 
@@ -227,7 +227,7 @@ FESystem<2>::multiply_dof_numbers (const FiniteElementData<2> &fe1,
   return FiniteElementData<2> (fe1.dofs_per_vertex * N1 + fe2.dofs_per_vertex * N2 ,
                               fe1.dofs_per_line * N1 + fe2.dofs_per_line * N2 ,
                               fe1.dofs_per_quad * N1 + fe2.dofs_per_quad * N2 ,
-                              fe1.n_transform_functions + fe2.n_transform_functions ,
+                              fe1.n_transform_functions,
                               fe1.n_components * N1 + fe2.n_components * N2 );
 };
 
@@ -374,82 +374,126 @@ void FESystem<dim>::get_local_mass_matrix (const DoFHandler<dim>::cell_iterator
 
 
 template <int dim>
-double FESystem<dim>::shape_value_transform (const unsigned int /*i*/,
-                                            const Point<dim>  &/*p*/) const
+double FESystem<dim>::shape_value_transform (const unsigned int i,
+                                            const Point<dim>  &p) const
 {
-  Assert(false, ExcNotImplemented());
-  return 0.;
+  return base_elements[0].first->shape_value_transform(i,p);
 };
 
 
 template <int dim>
-Tensor<1,dim> FESystem<dim>::shape_grad_transform (const unsigned int /*i*/,
-                                                  const Point<dim>  &/*p*/) const
+Tensor<1,dim> FESystem<dim>::shape_grad_transform (const unsigned int i,
+                                                  const Point<dim>  &p) const
 {
-  Assert(false, ExcNotImplemented());
-  return Tensor<1,dim>();
-
-//  return base_element->shape_grad_transform (i, p);
+  return base_elements[0].first->shape_grad_transform (i, p);
 };
 
 
 
 template <int dim>
-void FESystem<dim>::get_face_jacobians (const DoFHandler<dim>::face_iterator &/*face*/,
-                                       const Boundary<dim>         &/*boundary*/,
-                                       const vector<Point<dim-1> > &/*unit_points*/,
-                                       vector<double>      &/*face_jacobi_determinants*/) const
+void FESystem<dim>::get_face_jacobians (const DoFHandler<dim>::face_iterator &face,
+                                       const Boundary<dim>         &boundary,
+                                       const vector<Point<dim-1> > &unit_points,
+                                       vector<double>      &face_jacobi_determinants) const
 {
-  Assert(false, ExcNotImplemented());
-
-//  base_element->get_face_jacobians (face, boundary, unit_points, face_jacobi_determinants);
+  base_elements[0].first->get_face_jacobians (face, boundary, unit_points,
+                                       face_jacobi_determinants);
 };
 
 
 
 template <int dim>
-void FESystem<dim>::get_subface_jacobians (const DoFHandler<dim>::face_iterator &/*face*/,
-                                          const unsigned int           /*subface_no*/,
-                                          const vector<Point<dim-1> > &/*unit_points*/,
-                                          vector<double>      &/*face_jacobi_determinants*/) const
+void FESystem<dim>::get_subface_jacobians (const DoFHandler<dim>::face_iterator &face,
+                                          const unsigned int           subface_no,
+                                          const vector<Point<dim-1> > &unit_points,
+                                          vector<double>      &face_jacobi_determinants) const
 {
-  Assert(false, ExcNotImplemented());
-
-//  base_element->get_subface_jacobians (face, subface_no, unit_points, face_jacobi_determinants);
+  base_elements[0].first->get_subface_jacobians (face, subface_no, unit_points,
+                                          face_jacobi_determinants);
 };
 
 
 
 
 template <int dim>
-void FESystem<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &/*cell*/,
-                                       const unsigned int          /*face_no*/,
-                                       const Boundary<dim>         &/*boundary*/,
-                                       const vector<Point<dim-1> > &/*unit_points*/,
-                                       vector<Point<dim> >         &/*normal_vectors*/) const
+void FESystem<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
+                                       const unsigned int          face_no,
+                                       const Boundary<dim>         &boundary,
+                                       const vector<Point<dim-1> > &unit_points,
+                                       vector<Point<dim> >         &normal_vectors) const
 {
-  Assert(false, ExcNotImplemented());
-
-//  base_element->get_normal_vectors (cell, face_no, boundary, unit_points, normal_vectors);
+  base_elements[0].first->get_normal_vectors (cell, face_no, boundary, unit_points,
+                                       normal_vectors);
 };
 
 
 
 template <int dim>
-void FESystem<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &/*cell*/,
-                                       const unsigned int          /*face_no*/,
-                                       const unsigned int          /*subface_no*/,
-                                       const vector<Point<dim-1> > &/*unit_points*/,
-                                       vector<Point<dim> >         &/*normal_vectors*/) const
+void FESystem<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
+                                       const unsigned int          face_no,
+                                       const unsigned int          subface_no,
+                                       const vector<Point<dim-1> > &unit_points,
+                                       vector<Point<dim> >         &normal_vectors) const
 {
-  Assert(false, ExcNotImplemented());
-
-//  base_element->get_normal_vectors (cell, face_no, subface_no, unit_points, normal_vectors);
+  base_elements[0].first->get_normal_vectors (cell, face_no, subface_no, unit_points,
+                                   normal_vectors);
 };
 
 
+template <int dim>
+void
+FESystem<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
+                              const vector<Point<dim> >            &unit_points,
+                              vector<Tensor<2,dim> >               &jacobians,
+                              const bool              compute_jacobians,
+                              vector<Tensor<3,dim> > &jacobians_grad,
+                              const bool              compute_jacobians_grad,
+                              vector<Point<dim> > &support_points,
+                              const bool           compute_support_points,
+                              vector<Point<dim> > &q_points,
+                              const bool           compute_q_points,
+                              const dFMatrix      &shape_values_transform,
+                              const vector<vector<Tensor<1,dim> > > &shape_grad_transform,
+                              const Boundary<dim> &boundary) const
+{
+  vector<Point<dim> > supp(base_elements[0].first->total_dofs);
 
+  base_elements[0].first->fill_fe_values (cell, unit_points, jacobians, compute_jacobians,
+                                         jacobians_grad, compute_jacobians_grad,
+                                         support_points, compute_support_points,
+                                         q_points, compute_q_points,
+                                         shape_values_transform, shape_grad_transform, boundary);
+  
+  if (compute_support_points)
+    {
+      unsigned component = 0;
+      
+      for (unsigned m=0 ; m < element_multiplicity(0) ; ++ m)
+       {
+         for (unsigned i=0 ; i < base_element(0).total_dofs ; ++i)
+           support_points[component_to_system_index(component,i)] = supp[i];
+         ++component;
+       }
+      for (unsigned base=1 ; base < n_base_elements() ; ++base)
+       {
+         supp.resize(base_elements[base].first->total_dofs);
+         base_elements[base].first->fill_fe_values (cell, unit_points, jacobians, false,
+                                                    jacobians_grad, false,
+                                                    supp, true,
+                                                    q_points, false,
+                                                    shape_values_transform, shape_grad_transform, boundary);
+         
+         for (unsigned m=0 ; m < element_multiplicity(base) ; ++ m)
+           {
+             for (unsigned i=0 ; i < base_element(base).total_dofs ; ++i)
+               support_points[component_to_system_index(component,i)] = supp[i];
+             ++component;
+           }
+       }    
+    }
+}
 
+    
 
 // explicit instantiations
 template class FESystem<deal_II_dimension>;
index 8b51c2e287e83d9e7c6a6e2cbfec2515f3b165f3..870594cd7a29178e3a653617ac1956634eb72bc9 100644 (file)
@@ -88,11 +88,42 @@ void FEValuesBase<dim>::get_function_values (const dVector  &fe_function,
 };
 
 
+template <int dim>
+void FEValuesBase<dim>::get_function_values (const dVector  &fe_function,
+                                            vector< vector<double> > &values) const
+{
+  Assert (fe->n_components == values.size(),
+         ExcWrongNoOfComponents());
+  Assert (selected_dataset<shape_values.size(),
+         ExcInvalidIndex (selected_dataset, shape_values.size()));
+  for (unsigned i=0;i<values.size();++i)
+    Assert (values[i].size() == n_quadrature_points,
+           ExcWrongVectorSize(values.size(), n_quadrature_points));
+
+                                  // get function values of dofs
+                                  // on this cell
+  dVector dof_values (total_dofs);
+  present_cell->get_dof_values (fe_function, dof_values);
+
+                                  // initialize with zero
+  for (unsigned i=0;i<values.size();++i)
+    fill_n (values[i].begin(), n_quadrature_points, 0);
+
+                                  // add up contributions of trial
+                                  // 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]
+       += (dof_values(shape_func) * shape_values[selected_dataset](shape_func, point));
+};
+
+
 
 template <int dim>
 const Tensor<1,dim> &
 FEValuesBase<dim>::shape_grad (const unsigned int i,
-                              const unsigned int j) const {
+                              const unsigned int j) const
+{
   Assert (i<shape_gradients.size(),
          ExcInvalidIndex (i, shape_gradients.size()));
   Assert (j<shape_gradients[i].size(),
index 8e3b0de24e53d4ac63c13dfb677c006257674d31..49edbc694cec749a40de54cfeb8a233e51fca346 100644 (file)
@@ -543,7 +543,8 @@ 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<double> *values = new vector<double> [data.size()];
+  vector< vector <vector<double> > >
+    values (data.size(), vector< vector<double> >(dofs->get_fe().n_components, vector<double>(points.n_quadrature_points)));
 
   for (cell=dofs->begin_active(); cell!=endc; ++cell) 
     {
@@ -551,7 +552,7 @@ void DataOut<dim>::write_gnuplot (ostream &out, unsigned int accuracy) const
 
       for (unsigned i=0; i<data.size(); ++i)
       {
-       values[i].resize(points.n_quadrature_points);
+//     values[i].resize(points.n_quadrature_points);
        fe.get_function_values(*data[i].data, values[i]);
       }
       
@@ -568,8 +569,9 @@ void DataOut<dim>::write_gnuplot (ostream &out, unsigned int accuracy) const
                    Point<dim> pt = fe.quadrature_point(supp_pt);
                    out << pt << "  ";
                    for (unsigned int i=0; i!=data.size(); ++i)
-                     out << values[i][supp_pt]
-                         << ' ';
+                     for (unsigned int j=0; j < dofs->get_fe().n_components; ++j)
+                       out << values[i][j][supp_pt]
+                           << ' ';
                    out << endl;
                  };
                
@@ -589,8 +591,9 @@ void DataOut<dim>::write_gnuplot (ostream &out, unsigned int accuracy) const
                    out << pt << "  ";
                    
                    for (unsigned int i=0; i!=data.size(); ++i)
-                     out << values[i][supp_pt]
-                         << ' ';
+                     for (unsigned int j=0; j < dofs->get_fe().n_components; ++j)
+                       out << values[i][j][supp_pt]
+                           << ' ';
                    out << endl;
                  }
                  out << endl;
@@ -621,8 +624,9 @@ void DataOut<dim>::write_gnuplot (ostream &out, unsigned int accuracy) const
                            out << pt << "  ";
                            
                            for (unsigned int i=0; i!=data.size(); ++i)
-                             out << values[i][supp_pt]
-                                 << ' ';
+                             for (unsigned int j=0; j < dofs->get_fe().n_components; ++j)
+                               out << values[i][j][supp_pt]
+                                   << ' ';
                            out << endl;
                          }
                        out << endl;
@@ -638,7 +642,6 @@ void DataOut<dim>::write_gnuplot (ostream &out, unsigned int accuracy) const
        }
     out << endl;
     }
-  delete [] values;
 
   AssertThrow (out, ExcIO());
 }
index 9d8213c380d982632e6a1343c03839d339f3c494..58fe08915ee436aeebe2a80a69e161f683a68f72 100644 (file)
@@ -555,5 +555,193 @@ void VectorTools<dim>::integrate_difference (const DoFHandler<dim>    &dof,
     };
 };
 
+template <int dim>
+void VectorTools<dim>::integrate_difference (const DoFHandler<dim>    &dof,
+                                            const dVector            &fe_function,
+                                            const TensorFunction<1,dim>      &exact_solution,
+                                            dVector                  &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());
+  
+  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;
+         };
+                                            
+         default:
+               Assert (false, ExcNotImplemented());
+       };
+
+      
+                                      // append result of this cell
+                                      // to the end of the vector
+      difference(index) = diff;
+    };
+};
+
 
 template VectorTools<deal_II_dimension>;

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.