]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Factor out common code.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 25 Aug 2006 17:23:48 +0000 (17:23 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 25 Aug 2006 17:23:48 +0000 (17:23 +0000)
git-svn-id: https://svn.dealii.org/trunk@13748 0785d39b-7218-0410-832d-ea1e28bc413d

tests/fe/abf_projection_01.cc
tests/fe/abf_projection_01/cmp/generic
tests/fe/deformed_projection.h [new file with mode: 0644]
tests/fe/fe_sys_projection_01.cc
tests/fe/rt_projection_01.cc
tests/fe/rt_projection_01/cmp/generic

index d48ab215978b19f293c4a3b091e4e04cbf8a9dc7..3b4decec2fb7453b9c80204cbc393dbadd6e5de3 100644 (file)
@@ -1,5 +1,5 @@
-//----------------------------  abf_projection_01.cc  ---------------------------
-//    abf_projection_01.cc,v 1.3 2003/06/09 16:00:38 wolf Exp
+//----------------------------  abfprojection_01.cc  ---------------------------
+//    abfprojection_01.cc,v 1.3 2003/06/09 16:00:38 wolf Exp
 //    Version: 
 //
 //    Copyright (C) 2003, 2005, 2006 by the deal.II authors
@@ -9,7 +9,7 @@
 //    to the file deal.II/doc/license.html for the  text  and
 //    further information on this license.
 //
-//----------------------------  abf_projection_01.cc  ---------------------------
+//----------------------------  abfprojection_01.cc  ---------------------------
 
 /*
  * Project the function [1,1] onto a deformed grid and see whether the ABF
 
 
 
-#include "../tests.h"
-#include <base/logstream.h>
+char logname[] = "abf_projection_01/output";
+#include "deformed_projection.h"
 
-#define PRECISION 2
 
-#include <fstream>
-
-std::ofstream logfile ("abf_projection_01/output");
-
-char buf[1000];
-#include <fstream>
-
-#include <grid/grid_generator.h>
-#include <grid/grid_in.h>
-#include <grid/grid_out.h>
-#include <grid/tria_boundary_lib.h>
-
-#include <fe/mapping_q.h>
-#include <fe/mapping_q1_eulerian.h>
-
-#include <numerics/vectors.h>
-#include <numerics/matrices.h>
-#include <numerics/data_out.h>
-#include <lac/vector.h>
-#include <lac/solver_cg.h>
-#include <lac/precondition.h>
-#include <lac/vector_memory.h>
-#include <lac/sparsity_pattern.h>
-#include <lac/sparse_matrix.h>
-
-#include <base/function.h>
-#include <base/quadrature_lib.h>
-#include <dofs/dof_constraints.h>
-#include <dofs/dof_tools.h>
-
-#include <fe/fe_system.h>
-#include <fe/fe.h>
-#include <fe/fe_q.h>
-#include <fe/fe_values.h>
-#include <fe/fe_raviart_thomas.h>
-#include <fe/fe_abf.h>
-#include <fe/fe_dgq.h>
-
-
-template <int dim>
-class TestMap1 : public Function<dim>
-{
-  public:
-    TestMap1 (const unsigned int n_components) :
-                   Function<dim> (n_components)
-      {};
-  
-    virtual ~TestMap1 () {};
-  
-    virtual double value (const Point<dim>   &p,
-                         const unsigned int  component = 0) const;
-  
-    void vector_value (const Point<dim> &p,
-                      Vector<double>   &return_value) const;
-};
-
-
-
-template <int dim>
-double
-TestMap1<dim>::value (const Point<dim>   &,
-                     const unsigned int  ) const 
-{
-  return (1.0);
-}
-
-
-template <int dim>
-void TestMap1<dim>::vector_value (const Point<dim> &p,
-                                 Vector<double>   &return_value) const
-{
-  Assert (return_value.size() == this->n_components,
-         ExcDimensionMismatch (return_value.size(), this->n_components));
-  
-                                  // Parabolic inflow profile
-  for (unsigned int iCount = 0; iCount < this->n_components; iCount++)
-    return_value (iCount) = value (p, iCount);
-}
-
-
-/*
- * Check the value of the derivative field.
- */
-
-void EvaluateDerivative (DoFHandler<2> *dof_handler,
-                        Vector<double> &solution)
-{
-                                  // This quadrature rule determines the points, where the
-                                  // derivative will be evaluated.
-  QGauss<2> quad (3);
-  FEValues<2> fe_values (dof_handler->get_fe (), quad, 
-                        UpdateFlags(update_values    |
-                                    update_q_points  |
-                                    update_gradients |
-                                    update_JxW_values));
-
-  const unsigned int   n_q_points    = quad.n_quadrature_points;
-  const unsigned int   n_components   = dof_handler->get_fe().n_components();
-  const unsigned int dofs_per_cell = dof_handler->get_fe().dofs_per_cell;
-
-                                  // Cell iterators
-  DoFHandler<2>::active_cell_iterator cell = dof_handler->begin_active(),
-                                     endc = dof_handler->end();
-
-  double err_l2 = 0,
-       err_hdiv = 0;
-
-  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-
-  for (; cell!=endc; ++cell)
-    {
-      cell->get_dof_indices (local_dof_indices);
-    
-      fe_values.reinit (cell);
-
-                                      // Get function values
-      std::vector<Vector<double> > this_value(n_q_points,
-                                             Vector<double>(n_components));
-      fe_values.get_function_values (solution, this_value);
-
-                                      // Get values from solution vector (For Trap.Rule)
-      std::vector<std::vector<Tensor<1,2> > >
-       grads_here (n_q_points, std::vector<Tensor<1,2> > (n_components));
-      fe_values.get_function_grads (solution, grads_here);
-
-      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-       {
-                                          //     double u0 = this_value[q_point](0);
-                                          //double v0 = this_value[q_point](1);
-
-         double u0 = 0;
-         double v0 = 0;
-         for (unsigned int i = 0; i < dofs_per_cell; ++i)
-           {
-             u0 += (solution (local_dof_indices[i]) *
-                    fe_values.shape_value_component(i, q_point, 0));
-             v0 += (solution (local_dof_indices[i]) *
-                    fe_values.shape_value_component(i, q_point, 1));
-           }
-         u0 -= 1.0;
-         v0 -= 1.0;
-
-         double dudx = grads_here[q_point][0][0];
-         double dvdy = grads_here[q_point][1][1];
-
-                                          /*
-                                            sprintf (buf, "QP %i %e %e %e %e %e %e Div %e\n", q_point, u0, v0,
-                                            dudx, dudy, dvdx, dvdy, dudx + dvdy);
-                                          */
-
-         err_l2 += (u0 * u0 + v0 * v0) * fe_values.JxW (q_point);
-         err_hdiv += (dudx + dvdy) * (dudx + dvdy) * fe_values.JxW (q_point);
-       }
-    }
-
-  sprintf (buf, "L2-Err %e  Hdiv-Err %e\n", pow (err_l2, 0.5),
-          pow (err_hdiv, 0.5));
-  deallog << buf;
-}
-
-
-template <int dim>
-void create_mass_matrix (const Mapping<dim>       &mapping,
-                        const DoFHandler<dim>    &dof,
-                        const Quadrature<dim>    &q,
-                        SparseMatrix<double>     &matrix,
-                        const Function<dim>   &rhs_function,
-                        Vector<double>        &rhs_vector,
-                        const Function<dim> * const coefficient = 0)
+void test ()
 {
-  UpdateFlags update_flags = UpdateFlags(update_values | update_JxW_values | update_q_points);
-  if (coefficient != 0)
-    update_flags = UpdateFlags (update_flags | update_q_points);
-  
-  FEValues<dim> fe_values (mapping, dof.get_fe(), q, update_flags);
-    
-  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
-                    n_q_points    = fe_values.n_quadrature_points;
-  const FiniteElement<dim>    &fe  = fe_values.get_fe();
-  const unsigned int n_components  = fe.n_components();
-
-  Assert(coefficient == 0 ||
-        coefficient->n_components==1 ||
-        coefficient->n_components==n_components, ExcInternalError());
-  
-  FullMatrix<double>  cell_matrix (dofs_per_cell, dofs_per_cell);
-  Vector<double> cell_vector (dofs_per_cell);
-  std::vector<double> coefficient_values (n_q_points);
-  std::vector<Vector<double> > coefficient_vector_values (n_q_points,
-                                                         Vector<double> (n_components));
-  
-  std::vector<unsigned int> dof_indices (dofs_per_cell);
-
-  std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
-  
-  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active (),
-                                                endc = dof.end ();
-  for (; cell!=endc; ++cell)
-    {
-      fe_values.reinit (cell);
-      
-      cell_matrix = 0;
-      cell->get_dof_indices (dof_indices);
-
-      const std::vector<double> &weights   = fe_values.get_JxW_values ();
-      rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
-      cell_vector = 0;
-      
-      if (coefficient != 0)
-       {
-         if (coefficient->n_components==1)
-           {
-             coefficient->value_list (fe_values.get_quadrature_points(),
-                                      coefficient_values);
-             for (unsigned int point=0; point<n_q_points; ++point)
-               {
-                 const double weight = fe_values.JxW(point);
-                 for (unsigned int i=0; i<dofs_per_cell; ++i)
-                   {
-                     const double v = fe_values.shape_value(i,point);
-                     for (unsigned int j=0; j<dofs_per_cell; ++j)
-                       {
-                         const double u = fe_values.shape_value(j,point);
-                         
-                         if ((n_components==1) ||
-                             (fe.system_to_component_index(i).first ==
-                              fe.system_to_component_index(j).first))
-                           cell_matrix(i,j) +=
-                             (u * v * weight * coefficient_values[point]);
-                       }
-                   }
-               }
-           }
-         else
-           {
-             coefficient->vector_value_list (fe_values.get_quadrature_points(),
-                                             coefficient_vector_values);
-             for (unsigned int point=0; point<n_q_points; ++point)
-               {
-                 const double weight = fe_values.JxW(point);
-                 for (unsigned int i=0; i<dofs_per_cell; ++i)
-                   {
-                     const double v = fe_values.shape_value(i,point);
-                     const unsigned int component_i=
-                       fe.system_to_component_index(i).first;
-                     for (unsigned int j=0; j<dofs_per_cell; ++j)
-                       {
-                         const double u = fe_values.shape_value(j,point);
-                         if ((n_components==1) ||
-                             (fe.system_to_component_index(j).first == component_i))
-                           cell_matrix(i,j) +=
-                             (u * v * weight *
-                              coefficient_vector_values[point](component_i));
-                       }
-                   }
-               }
-           }
-       }
-      else
-       {
-                                          // Compute eventual sign changes depending on the neighborhood
-                                          // between two faces.
-         std::vector<double> sign_change (dofs_per_cell, 1.0);
-         const unsigned int dofs_per_face = fe.dofs_per_face;
-         std::vector<unsigned int> face_dof_indices (dofs_per_face);
-         
-         for (unsigned int f = 0; f < 2; ++f)
-           {
-             typename DoFHandler<dim>::active_face_iterator face = cell->face (f);
-             if (!face->at_boundary ())
-               {
-                 unsigned int nn = cell->neighbor_of_neighbor (f);
-                 sprintf (buf, "Face %i  NeigNeig %i\n", f, nn);
-                 deallog << buf;
-                 if (nn > 1)
-                   {
-                     face->get_dof_indices (face_dof_indices);
-                     for (unsigned int j = 0; j < dofs_per_face; ++j)
-                       {
-                         sign_change[face_dof_indices[j]] = -1.0;
-                         sprintf (buf, "DoF %i\n", face_dof_indices[j]);
-                         deallog << buf;
-                       }
-                   }
-               }
-           }
-         
-         for (unsigned int point=0; point<n_q_points; ++point)
-           {
-             const double weight = fe_values.JxW(point);
-                                              //           const double weight = q.weight(point);
-             
-             std::vector<Vector<double> > val_vector (dofs_per_cell,
-                                                      Vector<double> (n_components));
-             
-                                              // Precompute the component values
-             for (unsigned int i=0; i < dofs_per_cell; ++i)
-               for (unsigned int comp_i = 0; comp_i < fe.n_components (); 
-                    ++comp_i)
-                 {
-                   val_vector[i](comp_i) = sign_change[i] * 
-                                           fe_values.shape_value_component(i,point,comp_i);
-                 }
-                                              // Now eventually switch the sign of some of the ansatzfunctions.
-                                              // TODO
-             
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
-               for (unsigned int comp_i = 0; comp_i < fe.n_components (); 
-                    ++comp_i)
-                 if (fe.get_nonzero_components(i)[comp_i] == true)
-                   {
-                     const double v = val_vector[i](comp_i);
-                     for (unsigned int j=0; j<dofs_per_cell; ++j)
-                       for (unsigned int comp_j = 0;
-                            comp_j < fe.n_components (); ++comp_j)
-                         if (fe.get_nonzero_components(j)[comp_j] == true)
-                           {
-                             const double u = val_vector[j](comp_j);
-                             if ((n_components==1) ||
-                                 (comp_i == comp_j))
-                               cell_matrix(i,j) += (u * v * weight);
-                           }
-                   }
-             
-             
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
-               for (unsigned int comp_i = 0; comp_i < fe.n_components (); 
-                    ++comp_i)
-                 if (fe.get_nonzero_components(i)[comp_i] == true)
-                   {
-                     cell_vector(i) += rhs_values[point](comp_i) *
-                                       val_vector[i](comp_i) * weights[point];
-                   }
-           }
-       }
-                                      // transfer everything into the
-                                      // global object. lock the
-                                      // matrix meanwhile
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       for (unsigned int j=0; j<dofs_per_cell; ++j)
-         if ((n_components==1) ||
-             (cell_matrix (i,j) != 0.0))
-/*
-  (fe.system_to_component_index(i).first ==
-  fe.system_to_component_index(j).first))
-*/
-           matrix.add (dof_indices[i], dof_indices[j],
-                       cell_matrix(i,j));
-
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       rhs_vector(dof_indices[i]) += cell_vector(i);
-    }
-}
-
-
-template <int dim>
-void create_right_hand_side (const Mapping<dim>    &mapping,
-                            const DoFHandler<dim> &dof_handler,
-                            const Quadrature<dim> &quadrature,
-                            const Function<dim>   &rhs_function,
-                            Vector<double>        &rhs_vector)
-{
-  const FiniteElement<dim> &fe  = dof_handler.get_fe();
-  Assert (fe.n_components() == rhs_function.n_components,
-         ExcInternalError());
-  Assert (rhs_vector.size() == dof_handler.n_dofs(),
-         ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
-  rhs_vector = 0;
-  
-  UpdateFlags update_flags = UpdateFlags(update_values   |
-                                        update_q_points |
-                                        update_JxW_values);
-  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,
-                    n_components  = fe.n_components();
-  
-  std::vector<unsigned int> dofs (dofs_per_cell);
-  Vector<double> cell_vector (dofs_per_cell);
-
-  typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
-                                                endc = dof_handler.end();
-
-  if (n_components==1)
-    {
-      std::vector<double> rhs_values(n_q_points);
-      
-      for (; cell!=endc; ++cell) 
-       {
-         fe_values.reinit(cell);
-         
-         const std::vector<double> &weights   = fe_values.get_JxW_values ();
-         rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values);
-         
-         cell_vector = 0;
-         for (unsigned int point=0; point<n_q_points; ++point)
-           for (unsigned int i=0; i<dofs_per_cell; ++i) 
-             cell_vector(i) += rhs_values[point] *
-                               fe_values.shape_value(i,point) *
-                               weights[point];
-       
-         cell->get_dof_indices (dofs);
-         
-         for (unsigned int i=0; i<dofs_per_cell; ++i)
-           rhs_vector(dofs[i]) += cell_vector(i);
-       }
-    }
-  else
-    {
-      std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
-      
-      for (; cell!=endc; ++cell) 
-       {
-         fe_values.reinit(cell);
-         
-         const std::vector<double> &weights   = fe_values.get_JxW_values ();
-         rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
-         
-         cell_vector = 0;
-         for (unsigned int point=0; point<n_q_points; ++point)
-           for (unsigned int i=0; i<dofs_per_cell; ++i)
-             for (unsigned int comp_i = 0; comp_i < fe.n_components (); 
-                  ++comp_i)
-                                                //                 if (fe.get_nonzero_components(i)[comp_i] == true)
-               {
-                 cell_vector(i) += rhs_values[point](comp_i) *
-                                   fe_values.shape_value_component(i,point,comp_i) *
-                                   weights[point];
-               }
-         
-         cell->get_dof_indices (dofs);
-         
-         for (unsigned int i=0; i<dofs_per_cell; ++i)
-           rhs_vector(dofs[i]) += cell_vector(i);
-       }
-    }
-}
-
-
-
-//
-// This function replaces the deal.II implementation of the projection.
-// The purpose is to have more freedom in assembling the matrix.
-//
-
-template <int dim>
-void project (const Mapping<dim>       &mapping,
-             const DoFHandler<dim>    &dof,
-             const ConstraintMatrix   &constraints,
-             const Quadrature<dim>    &quadrature,
-             const Function<dim>      &function,
-             Vector<double>           &vec,
-             const bool                enforce_zero_boundary = false,
-             const Quadrature<dim-1>  & = QGauss2<dim-1>(),
-             const bool                project_to_boundary_first = false)
-{
-  Assert (dof.get_fe().n_components() == function.n_components,
-         ExcInternalError());
-  
-  const FiniteElement<dim> &fe = dof.get_fe();
-
-                                  // make up boundary values
-  std::map<unsigned int,double> boundary_values;
-
-  if (enforce_zero_boundary == true) 
-                                    // no need to project boundary
-                                    // values, but enforce
-                                    // homogeneous boundary values
-                                    // anyway
-    {
-                                      // loop over all boundary faces
-                                      // to get all dof indices of
-                                      // dofs on the boundary. note
-                                      // that in 3d there are cases
-                                      // where a face is not at the
-                                      // boundary, yet one of its
-                                      // lines is, and we should
-                                      // consider the degrees of
-                                      // freedom on it as boundary
-                                      // nodes. likewise, in 2d and
-                                      // 3d there are cases where a
-                                      // cell is only at the boundary
-                                      // by one vertex. nevertheless,
-                                      // since we do not support
-                                      // boundaries with dimension
-                                      // less or equal to dim-2, each
-                                      // such boundary dof is also
-                                      // found from some other face
-                                      // that is actually wholly on
-                                      // the boundary, not only by
-                                      // one line or one vertex
-      typename DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
-                                                    endf = dof.end_face();
-      std::vector<unsigned int> face_dof_indices (fe.dofs_per_face);
-      for (; face!=endf; ++face)
-       if (face->at_boundary())
-         {
-           face->get_dof_indices (face_dof_indices);
-           for (unsigned int i=0; i<fe.dofs_per_face; ++i)
-                                              // enter zero boundary values
-                                              // for all boundary nodes
-                                              //
-                                              // we need not care about
-                                              // vector valued elements here,
-                                              // since we set all components
-             boundary_values[face_dof_indices[i]] = 0.;
-         };
-    }
-  else
-                                    // no homogeneous boundary values
-    if (project_to_boundary_first == true)
-                                      // boundary projection required
-      {
-/*
-                                // set up a list of boundary functions for
-                                                                 // the different boundary parts. We want the
-                                                                                                  // @p{function} to hold on all parts of the
-                                                                                                                                   // boundary
-                                                                                                                                   typename FunctionMap<dim>::type boundary_functions;
-                                                                                                                                   for (unsigned char c=0; c<255; ++c)
-                                                                                                                                   boundary_functions[c] = &function;
-                                                                                                                                   project_boundary_values (dof, boundary_functions, q_boundary,
-                                                                                                                                   boundary_values);
-*/
-      }
-
-
-                                  // set up mass matrix and right hand side
-  vec.reinit (dof.n_dofs());
-  SparsityPattern sparsity(dof.n_dofs(),
-                          dof.n_dofs(),
-                          dof.max_couplings_between_dofs());
-  DoFTools::make_sparsity_pattern (dof, sparsity);
-  constraints.condense (sparsity);
-  
-  SparseMatrix<double> mass_matrix (sparsity);
-  Vector<double> tmp (mass_matrix.n());
-
-  create_mass_matrix (mapping, dof, quadrature, mass_matrix, function, tmp);
-  sprintf (buf, "MM created\n");
-  deallog << buf;
-  create_right_hand_side (mapping, dof, quadrature, function, tmp);
-  sprintf (buf, "RHS created\n");
-  deallog << buf;
-
-  constraints.condense (mass_matrix);
-  constraints.condense (tmp);
-  if (boundary_values.size() != 0)
-    MatrixTools::apply_boundary_values (boundary_values,
-                                       mass_matrix, vec, tmp,
-                                       true);
-
-  SolverControl           control(1000,1e-16);
-  PrimitiveVectorMemory<> memory;
-  SolverCG<>              cg(control,memory);
-
-  PreconditionSSOR<> prec;
-  prec.initialize(mass_matrix, 1.2);
-                                  // solve
-  cg.solve (mass_matrix, vec, tmp, prec);
-  
-                                  // distribute solution
-  constraints.distribute (vec);
-}
-
-
-int create_tria (unsigned int elm, Triangulation<2> &tria)
-{
-  std::vector<Point<2> > points_glob;
-  std::vector<Point<2> > points;
-
-  points_glob.push_back (Point<2> (0.0, 0.0));
-  points_glob.push_back (Point<2> (1.0, 0.0));
-  points_glob.push_back (Point<2> (1.0, 0.5));
-  points_glob.push_back (Point<2> (1.0, 1.0));
-  points_glob.push_back (Point<2> (0.6, 0.5));
-  points_glob.push_back (Point<2> (0.5, 1.0));
-  points_glob.push_back (Point<2> (0.0, 1.0));
-
-                                  // Prepare cell data
-  std::vector<CellData<2> > cells (1);
-
-  switch (elm)
-    {
-      case 0:
-           cells[0].vertices[0] = 0;
-           cells[0].vertices[1] = 1;
-           cells[0].vertices[2] = 4;
-           cells[0].vertices[3] = 2;
-           cells[0].material_id = 0;
-           break;
-      case 1:
-           cells[0].vertices[0] = 4;
-           cells[0].vertices[1] = 2;
-           cells[0].vertices[2] = 5;
-           cells[0].vertices[3] = 3;
-           cells[0].material_id = 0;
-           break;
-      case 2:
-           cells[0].vertices[0] = 0;
-           cells[0].vertices[1] = 4;
-           cells[0].vertices[2] = 6;
-           cells[0].vertices[3] = 5;
-           cells[0].material_id = 0;
-           break;
-    }
-
-  for (unsigned int i = 0; i < 4; ++i)
-    {
-      points.push_back (points_glob[cells[0].vertices[i]]);
-      cells[0].vertices[i] = i;
-    }
-
-  tria.create_triangulation (points, cells, SubCellData());
-
-  return (0);
-}
-
-
-int plot_shapes (DoFHandler<2> &dof_handler)
-{
-  Vector<double> solution (dof_handler.n_dofs ());
-  std::set<unsigned int> face_dofs;
-  
-                                  // Create set of all DoFs, which are on the boundary.
-  DoFHandler<2>::active_face_iterator face = dof_handler.begin_active_face(),
-                                     endf = dof_handler.end_face();
-  std::vector<unsigned int> face_dof_indices (dof_handler.get_fe().dofs_per_face);
-  for (; face!=endf; ++face)
-    {
-      face->get_dof_indices (face_dof_indices);
-      for (unsigned int i=0; i<dof_handler.get_fe().dofs_per_face; ++i)
-       face_dofs.insert (face_dof_indices[i]);
-    }
-  
-                                  // Now set a 1 at the place of the different DoFs and
-                                  // output the solution.
-  std::set<unsigned int>::iterator face_dof_iter;
-  for (face_dof_iter = face_dofs.begin ();
-       face_dof_iter != face_dofs.end (); ++face_dof_iter)
-    {
-      unsigned int dof = *face_dof_iter;
-
-      sprintf (buf, "%i\n", dof);
-      deallog << buf;
-
-      solution (dof) = 1.0;
-
-                                      // Test the core functionality2
-      DataOut<2> *data_out = new DataOut<2>;
-      data_out->attach_dof_handler (dof_handler);
-      data_out->add_data_vector (solution, "solution");
-      data_out->build_patches (4, 0);
-
-      data_out->write_gnuplot (deallog.get_file_stream());
-
-      delete data_out;
-
-      solution (dof) = 0.0;
-    }
-
-  return (0);
-}
-
-
-int main (int /*argc*/, char **/*argv*/)
-{
-  logfile.precision (PRECISION);
-  logfile.setf(std::ios::fixed);  
-  deallog.attach(logfile);
-  deallog.depth_console(0);
-  deallog.threshold_double(1.e-10);
-
-  Triangulation<2> tria_test;
-  DoFHandler<2> *dof_handler;
-
   FE_ABF<2> fe (0);
-  deallog << "Dofs/cell " << fe.dofs_per_cell
-         << "Dofs/face " << fe.dofs_per_face << std::endl;
-
-  for (unsigned int elm = 0; elm < 3; ++elm)
-    {
-      create_tria (elm, tria_test);
-
-                                      // Create a DoFHandler
-      dof_handler = new DoFHandler<2> (tria_test);
-      dof_handler->distribute_dofs (fe);
-
-      deallog << "Dofs total " << dof_handler->n_dofs () << std::endl;
-
-                                      // Alloc some DoFs
-      Vector<double> solution;
-      solution.reinit (dof_handler->n_dofs ());
-
-                                      // Project solution onto FE field
-      ConstraintMatrix     hn_constraints;
-      hn_constraints.clear ();
-      DoFTools::make_hanging_node_constraints (*dof_handler, 
-                                              hn_constraints);
-      hn_constraints.close ();
-
-      MappingQ1<2> map_default;
-      
-      project (map_default, *dof_handler, hn_constraints,
-              QGauss6<2> (), TestMap1<2>(2),
-              solution);
-
-                                      // Test the core functionality
-      DataOut<2> *data_out = new DataOut<2>;
-      data_out->attach_dof_handler (*dof_handler);
-      data_out->add_data_vector (solution, "solution");
-      data_out->build_patches (4);
-
-      data_out->write_gnuplot (deallog.get_file_stream());
-
-                                      // Now write face DoFs ..
-      for (unsigned int i = 0; i < 4; ++i)
-       {
-         sprintf (buf, "%i %e\n", i, solution (i));
-         deallog << buf;
-       }
-      
-                                      // Clean up ...
-      delete data_out;
-      delete (dof_handler);
-      tria_test.clear ();
-    }
-
-  return (0);
+  check (fe);
 }
index 84f06527b87d972656f1e5638cec915249b878db..0f500d410f4504f1b3b488c243e6553fde79715c 100644 (file)
@@ -2,8 +2,8 @@
 DEAL::Dofs/cell 6Dofs/face 1
 DEAL::Dofs total 6
 DEAL::MM created
-RHS created
-Starting value 1.99
+DEAL::RHS created
+DEAL:cg::Starting value 1.99
 DEAL:cg::Convergence step 9 value 0
 # This file was generated by the deal.II library.
 
@@ -43,14 +43,14 @@ DEAL:cg::Convergence step 9 value 0
 1.00 0.50 1.00 1.00 
 
 
-DEAL::0 -1.000000e-01
-1 5.000000e-01
-2 1.000000e+00
-3 4.000000e-01
-Dofs total 6
+DEAL::-0.10
+DEAL::0.50
+DEAL::1.00
+DEAL::0.40
+DEAL::Dofs total 6
 DEAL::MM created
-RHS created
-Starting value 1.19
+DEAL::RHS created
+DEAL:cg::Starting value 1.19
 DEAL:cg::Convergence step 8 value 0
 # This file was generated by the deal.II library.
 
@@ -90,14 +90,14 @@ DEAL:cg::Convergence step 8 value 0
 1.00 1.00 1.00 1.00 
 
 
-DEAL::0 6.000000e-01
-1 5.000000e-01
-2 4.000000e-01
-3 5.000000e-01
-Dofs total 6
+DEAL::0.60
+DEAL::0.50
+DEAL::0.40
+DEAL::0.50
+DEAL::Dofs total 6
 DEAL::MM created
-RHS created
-Starting value 1.99
+DEAL::RHS created
+DEAL:cg::Starting value 1.99
 DEAL:cg::Convergence step 9 value 0
 # This file was generated by the deal.II library.
 
@@ -137,7 +137,7 @@ DEAL:cg::Convergence step 9 value 0
 0.50 1.00 1.00 1.00 
 
 
-DEAL::0 1.000000e+00
-1 6.000000e-01
-2 1.000000e-01
-3 5.000000e-01
+DEAL::1.00
+DEAL::0.60
+DEAL::0.10
+DEAL::0.50
diff --git a/tests/fe/deformed_projection.h b/tests/fe/deformed_projection.h
new file mode 100644 (file)
index 0000000..03e6274
--- /dev/null
@@ -0,0 +1,753 @@
+//----------------------------  deformed_projection.h  ---------------------------
+//    deformed_projection.h,v 1.3 2003/06/09 16:00:38 wolf Exp
+//    Version: 
+//
+//    Copyright (C) 2003, 2005, 2006 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
+//    further information on this license.
+//
+//----------------------------  deformed_projection.h  ---------------------------
+
+/*
+ * Project the function [1,1] onto a deformed grid and see whether the element
+ * elements can represent it exactly. This is the common code for the various
+ * *_projection_* tests.
+ */
+
+
+
+#include "../tests.h"
+#include <base/logstream.h>
+
+#define PRECISION 2
+
+#include <fstream>
+
+#include <grid/grid_generator.h>
+#include <grid/grid_in.h>
+#include <grid/grid_out.h>
+#include <grid/tria_boundary_lib.h>
+
+#include <fe/mapping_q.h>
+#include <fe/mapping_q1_eulerian.h>
+
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <lac/vector.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <lac/vector_memory.h>
+#include <lac/sparsity_pattern.h>
+#include <lac/sparse_matrix.h>
+
+#include <base/function.h>
+#include <base/quadrature_lib.h>
+#include <dofs/dof_constraints.h>
+#include <dofs/dof_tools.h>
+
+#include <fe/fe_system.h>
+#include <fe/fe.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+#include <fe/fe_raviart_thomas.h>
+#include <fe/fe_abf.h>
+#include <fe/fe_dgq.h>
+
+void test();
+
+
+template <int dim>
+class TestMap1 : public Function<dim>
+{
+  public:
+    TestMap1 (const unsigned int n_components) :
+                   Function<dim> (n_components)
+      {};
+  
+    virtual ~TestMap1 () {};
+  
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+  
+    void vector_value (const Point<dim> &p,
+                      Vector<double>   &return_value) const;
+};
+
+
+
+template <int dim>
+double
+TestMap1<dim>::value (const Point<dim>   &,
+                     const unsigned int  ) const 
+{
+  return (1.0);
+}
+
+
+template <int dim>
+void TestMap1<dim>::vector_value (const Point<dim> &p,
+                                 Vector<double>   &return_value) const
+{
+  Assert (return_value.size() == this->n_components,
+         ExcDimensionMismatch (return_value.size(), this->n_components));
+  
+                                  // Parabolic inflow profile
+  for (unsigned int iCount = 0; iCount < this->n_components; iCount++)
+    return_value (iCount) = value (p, iCount);
+}
+
+
+/*
+ * Check the value of the derivative field.
+ */
+
+void EvaluateDerivative (DoFHandler<2> *dof_handler,
+                        Vector<double> &solution)
+{
+                                  // This quadrature rule determines the points, where the
+                                  // derivative will be evaluated.
+  QGauss<2> quad (3);
+  FEValues<2> fe_values (dof_handler->get_fe (), quad, 
+                        UpdateFlags(update_values    |
+                                    update_q_points  |
+                                    update_gradients |
+                                    update_JxW_values));
+
+  const unsigned int   n_q_points    = quad.n_quadrature_points;
+  const unsigned int   n_components   = dof_handler->get_fe().n_components();
+  const unsigned int dofs_per_cell = dof_handler->get_fe().dofs_per_cell;
+
+                                  // Cell iterators
+  DoFHandler<2>::active_cell_iterator cell = dof_handler->begin_active(),
+                                     endc = dof_handler->end();
+
+  double err_l2 = 0,
+       err_hdiv = 0;
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+  for (; cell!=endc; ++cell)
+    {
+      cell->get_dof_indices (local_dof_indices);
+    
+      fe_values.reinit (cell);
+
+                                      // Get function values
+      std::vector<Vector<double> > this_value(n_q_points,
+                                             Vector<double>(n_components));
+      fe_values.get_function_values (solution, this_value);
+
+                                      // Get values from solution vector (For Trap.Rule)
+      std::vector<std::vector<Tensor<1,2> > >
+       grads_here (n_q_points, std::vector<Tensor<1,2> > (n_components));
+      fe_values.get_function_grads (solution, grads_here);
+
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+       {
+                                          //     double u0 = this_value[q_point](0);
+                                          //double v0 = this_value[q_point](1);
+
+         double u0 = 0;
+         double v0 = 0;
+         for (unsigned int i = 0; i < dofs_per_cell; ++i)
+           {
+             u0 += (solution (local_dof_indices[i]) *
+                    fe_values.shape_value_component(i, q_point, 0));
+             v0 += (solution (local_dof_indices[i]) *
+                    fe_values.shape_value_component(i, q_point, 1));
+           }
+         u0 -= 1.0;
+         v0 -= 1.0;
+
+         double dudx = grads_here[q_point][0][0];
+         double dvdy = grads_here[q_point][1][1];
+
+                                          /*
+                                            sprintf (buf, "QP %i %e %e %e %e %e %e Div %e\n", q_point, u0, v0,
+                                            dudx, dudy, dvdx, dvdy, dudx + dvdy);
+                                          */
+
+         err_l2 += (u0 * u0 + v0 * v0) * fe_values.JxW (q_point);
+         err_hdiv += (dudx + dvdy) * (dudx + dvdy) * fe_values.JxW (q_point);
+       }
+    }
+
+  deallog << "L2-Err " << pow (err_l2, 0.5)
+         << "  Hdiv-Err " << pow (err_hdiv, 0.5)
+         << std::endl;
+}
+
+
+template <int dim>
+void create_mass_matrix (const Mapping<dim>       &mapping,
+                        const DoFHandler<dim>    &dof,
+                        const Quadrature<dim>    &q,
+                        SparseMatrix<double>     &matrix,
+                        const Function<dim>   &rhs_function,
+                        Vector<double>        &rhs_vector,
+                        const Function<dim> * const coefficient = 0)
+{
+  UpdateFlags update_flags = UpdateFlags(update_values | update_JxW_values | update_q_points);
+  if (coefficient != 0)
+    update_flags = UpdateFlags (update_flags | update_q_points);
+  
+  FEValues<dim> fe_values (mapping, dof.get_fe(), q, update_flags);
+    
+  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+                    n_q_points    = fe_values.n_quadrature_points;
+  const FiniteElement<dim>    &fe  = fe_values.get_fe();
+  const unsigned int n_components  = fe.n_components();
+
+  Assert(coefficient == 0 ||
+        coefficient->n_components==1 ||
+        coefficient->n_components==n_components, ExcInternalError());
+  
+  FullMatrix<double>  cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double> cell_vector (dofs_per_cell);
+  std::vector<double> coefficient_values (n_q_points);
+  std::vector<Vector<double> > coefficient_vector_values (n_q_points,
+                                                         Vector<double> (n_components));
+  
+  std::vector<unsigned int> dof_indices (dofs_per_cell);
+
+  std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
+  
+  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active (),
+                                                endc = dof.end ();
+  for (; cell!=endc; ++cell)
+    {
+      fe_values.reinit (cell);
+      
+      cell_matrix = 0;
+      cell->get_dof_indices (dof_indices);
+
+      const std::vector<double> &weights   = fe_values.get_JxW_values ();
+      rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
+      cell_vector = 0;
+      
+      if (coefficient != 0)
+       {
+         if (coefficient->n_components==1)
+           {
+             coefficient->value_list (fe_values.get_quadrature_points(),
+                                      coefficient_values);
+             for (unsigned int point=0; point<n_q_points; ++point)
+               {
+                 const double weight = fe_values.JxW(point);
+                 for (unsigned int i=0; i<dofs_per_cell; ++i)
+                   {
+                     const double v = fe_values.shape_value(i,point);
+                     for (unsigned int j=0; j<dofs_per_cell; ++j)
+                       {
+                         const double u = fe_values.shape_value(j,point);
+                         
+                         if ((n_components==1) ||
+                             (fe.system_to_component_index(i).first ==
+                              fe.system_to_component_index(j).first))
+                           cell_matrix(i,j) +=
+                             (u * v * weight * coefficient_values[point]);
+                       }
+                   }
+               }
+           }
+         else
+           {
+             coefficient->vector_value_list (fe_values.get_quadrature_points(),
+                                             coefficient_vector_values);
+             for (unsigned int point=0; point<n_q_points; ++point)
+               {
+                 const double weight = fe_values.JxW(point);
+                 for (unsigned int i=0; i<dofs_per_cell; ++i)
+                   {
+                     const double v = fe_values.shape_value(i,point);
+                     const unsigned int component_i=
+                       fe.system_to_component_index(i).first;
+                     for (unsigned int j=0; j<dofs_per_cell; ++j)
+                       {
+                         const double u = fe_values.shape_value(j,point);
+                         if ((n_components==1) ||
+                             (fe.system_to_component_index(j).first == component_i))
+                           cell_matrix(i,j) +=
+                             (u * v * weight *
+                              coefficient_vector_values[point](component_i));
+                       }
+                   }
+               }
+           }
+       }
+      else
+       {
+                                          // Compute eventual sign changes depending on the neighborhood
+                                          // between two faces.
+         std::vector<double> sign_change (dofs_per_cell, 1.0);
+         const unsigned int dofs_per_face = fe.dofs_per_face;
+         std::vector<unsigned int> face_dof_indices (dofs_per_face);
+         
+         for (unsigned int f = 0; f < 2; ++f)
+           {
+             typename DoFHandler<dim>::active_face_iterator face = cell->face (f);
+             if (!face->at_boundary ())
+               {
+                 unsigned int nn = cell->neighbor_of_neighbor (f);
+                 deallog << "Face " << f << " NeigNeig " << nn << std::endl;
+                 if (nn > 1)
+                   {
+                     face->get_dof_indices (face_dof_indices);
+                     for (unsigned int j = 0; j < dofs_per_face; ++j)
+                       {
+                         sign_change[face_dof_indices[j]] = -1.0;
+                         deallog << "DoF " << face_dof_indices[j] << std::endl;
+                       }
+                   }
+               }
+           }
+         
+         for (unsigned int point=0; point<n_q_points; ++point)
+           {
+             const double weight = fe_values.JxW(point);
+                                              //           const double weight = q.weight(point);
+             
+             std::vector<Vector<double> > val_vector (dofs_per_cell,
+                                                      Vector<double> (n_components));
+             
+                                              // Precompute the component values
+             for (unsigned int i=0; i < dofs_per_cell; ++i)
+               for (unsigned int comp_i = 0; comp_i < fe.n_components (); 
+                    ++comp_i)
+                 {
+                   val_vector[i](comp_i) = sign_change[i] * 
+                                           fe_values.shape_value_component(i,point,comp_i);
+                 }
+                                              // Now eventually switch the sign of some of the ansatzfunctions.
+                                              // TODO
+             
+             for (unsigned int i=0; i<dofs_per_cell; ++i)
+               for (unsigned int comp_i = 0; comp_i < fe.n_components (); 
+                    ++comp_i)
+                 if (fe.get_nonzero_components(i)[comp_i] == true)
+                   {
+                     const double v = val_vector[i](comp_i);
+                     for (unsigned int j=0; j<dofs_per_cell; ++j)
+                       for (unsigned int comp_j = 0;
+                            comp_j < fe.n_components (); ++comp_j)
+                         if (fe.get_nonzero_components(j)[comp_j] == true)
+                           {
+                             const double u = val_vector[j](comp_j);
+                             if ((n_components==1) ||
+                                 (comp_i == comp_j))
+                               cell_matrix(i,j) += (u * v * weight);
+                           }
+                   }
+             
+             
+             for (unsigned int i=0; i<dofs_per_cell; ++i)
+               for (unsigned int comp_i = 0; comp_i < fe.n_components (); 
+                    ++comp_i)
+                 if (fe.get_nonzero_components(i)[comp_i] == true)
+                   {
+                     cell_vector(i) += rhs_values[point](comp_i) *
+                                       val_vector[i](comp_i) * weights[point];
+                   }
+           }
+       }
+                                      // transfer everything into the
+                                      // global object. lock the
+                                      // matrix meanwhile
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       for (unsigned int j=0; j<dofs_per_cell; ++j)
+         if ((n_components==1) ||
+             (cell_matrix (i,j) != 0.0))
+/*
+  (fe.system_to_component_index(i).first ==
+  fe.system_to_component_index(j).first))
+*/
+           matrix.add (dof_indices[i], dof_indices[j],
+                       cell_matrix(i,j));
+
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       rhs_vector(dof_indices[i]) += cell_vector(i);
+    }
+}
+
+
+template <int dim>
+void create_right_hand_side (const Mapping<dim>    &mapping,
+                            const DoFHandler<dim> &dof_handler,
+                            const Quadrature<dim> &quadrature,
+                            const Function<dim>   &rhs_function,
+                            Vector<double>        &rhs_vector)
+{
+  const FiniteElement<dim> &fe  = dof_handler.get_fe();
+  Assert (fe.n_components() == rhs_function.n_components,
+         ExcInternalError());
+  Assert (rhs_vector.size() == dof_handler.n_dofs(),
+         ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
+  rhs_vector = 0;
+  
+  UpdateFlags update_flags = UpdateFlags(update_values   |
+                                        update_q_points |
+                                        update_JxW_values);
+  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,
+                    n_components  = fe.n_components();
+  
+  std::vector<unsigned int> dofs (dofs_per_cell);
+  Vector<double> cell_vector (dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+                                                endc = dof_handler.end();
+
+  if (n_components==1)
+    {
+      std::vector<double> rhs_values(n_q_points);
+      
+      for (; cell!=endc; ++cell) 
+       {
+         fe_values.reinit(cell);
+         
+         const std::vector<double> &weights   = fe_values.get_JxW_values ();
+         rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values);
+         
+         cell_vector = 0;
+         for (unsigned int point=0; point<n_q_points; ++point)
+           for (unsigned int i=0; i<dofs_per_cell; ++i) 
+             cell_vector(i) += rhs_values[point] *
+                               fe_values.shape_value(i,point) *
+                               weights[point];
+       
+         cell->get_dof_indices (dofs);
+         
+         for (unsigned int i=0; i<dofs_per_cell; ++i)
+           rhs_vector(dofs[i]) += cell_vector(i);
+       }
+    }
+  else
+    {
+      std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
+      
+      for (; cell!=endc; ++cell) 
+       {
+         fe_values.reinit(cell);
+         
+         const std::vector<double> &weights   = fe_values.get_JxW_values ();
+         rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
+         
+         cell_vector = 0;
+         for (unsigned int point=0; point<n_q_points; ++point)
+           for (unsigned int i=0; i<dofs_per_cell; ++i)
+             for (unsigned int comp_i = 0; comp_i < fe.n_components (); 
+                  ++comp_i)
+                                                //                 if (fe.get_nonzero_components(i)[comp_i] == true)
+               {
+                 cell_vector(i) += rhs_values[point](comp_i) *
+                                   fe_values.shape_value_component(i,point,comp_i) *
+                                   weights[point];
+               }
+         
+         cell->get_dof_indices (dofs);
+         
+         for (unsigned int i=0; i<dofs_per_cell; ++i)
+           rhs_vector(dofs[i]) += cell_vector(i);
+       }
+    }
+}
+
+
+
+//
+// This function replaces the deal.II implementation of the projection.
+// The purpose is to have more freedom in assembling the matrix.
+//
+
+template <int dim>
+void project (const Mapping<dim>       &mapping,
+             const DoFHandler<dim>    &dof,
+             const ConstraintMatrix   &constraints,
+             const Quadrature<dim>    &quadrature,
+             const Function<dim>      &function,
+             Vector<double>           &vec,
+             const bool                enforce_zero_boundary = false,
+             const Quadrature<dim-1>  & = QGauss2<dim-1>(),
+             const bool                project_to_boundary_first = false)
+{
+  Assert (dof.get_fe().n_components() == function.n_components,
+         ExcInternalError());
+  
+  const FiniteElement<dim> &fe = dof.get_fe();
+
+                                  // make up boundary values
+  std::map<unsigned int,double> boundary_values;
+
+  if (enforce_zero_boundary == true) 
+                                    // no need to project boundary
+                                    // values, but enforce
+                                    // homogeneous boundary values
+                                    // anyway
+    {
+                                      // loop over all boundary faces
+                                      // to get all dof indices of
+                                      // dofs on the boundary. note
+                                      // that in 3d there are cases
+                                      // where a face is not at the
+                                      // boundary, yet one of its
+                                      // lines is, and we should
+                                      // consider the degrees of
+                                      // freedom on it as boundary
+                                      // nodes. likewise, in 2d and
+                                      // 3d there are cases where a
+                                      // cell is only at the boundary
+                                      // by one vertex. nevertheless,
+                                      // since we do not support
+                                      // boundaries with dimension
+                                      // less or equal to dim-2, each
+                                      // such boundary dof is also
+                                      // found from some other face
+                                      // that is actually wholly on
+                                      // the boundary, not only by
+                                      // one line or one vertex
+      typename DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
+                                                    endf = dof.end_face();
+      std::vector<unsigned int> face_dof_indices (fe.dofs_per_face);
+      for (; face!=endf; ++face)
+       if (face->at_boundary())
+         {
+           face->get_dof_indices (face_dof_indices);
+           for (unsigned int i=0; i<fe.dofs_per_face; ++i)
+                                              // enter zero boundary values
+                                              // for all boundary nodes
+                                              //
+                                              // we need not care about
+                                              // vector valued elements here,
+                                              // since we set all components
+             boundary_values[face_dof_indices[i]] = 0.;
+         };
+    }
+  else
+                                    // no homogeneous boundary values
+    if (project_to_boundary_first == true)
+                                      // boundary projection required
+      {
+/*
+                                // set up a list of boundary functions for
+                                                                 // the different boundary parts. We want the
+                                                                                                  // @p{function} to hold on all parts of the
+                                                                                                                                   // boundary
+                                                                                                                                   typename FunctionMap<dim>::type boundary_functions;
+                                                                                                                                   for (unsigned char c=0; c<255; ++c)
+                                                                                                                                   boundary_functions[c] = &function;
+                                                                                                                                   project_boundary_values (dof, boundary_functions, q_boundary,
+                                                                                                                                   boundary_values);
+*/
+      }
+
+
+                                  // set up mass matrix and right hand side
+  vec.reinit (dof.n_dofs());
+  SparsityPattern sparsity(dof.n_dofs(),
+                          dof.n_dofs(),
+                          dof.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern (dof, sparsity);
+  constraints.condense (sparsity);
+  
+  SparseMatrix<double> mass_matrix (sparsity);
+  Vector<double> tmp (mass_matrix.n());
+
+  create_mass_matrix (mapping, dof, quadrature, mass_matrix, function, tmp);
+  deallog << "MM created" << std::endl;
+  create_right_hand_side (mapping, dof, quadrature, function, tmp);
+  deallog << "RHS created" << std::endl;
+
+  constraints.condense (mass_matrix);
+  constraints.condense (tmp);
+  if (boundary_values.size() != 0)
+    MatrixTools::apply_boundary_values (boundary_values,
+                                       mass_matrix, vec, tmp,
+                                       true);
+
+  SolverControl           control(1000,1e-16);
+  PrimitiveVectorMemory<> memory;
+  SolverCG<>              cg(control,memory);
+
+  PreconditionSSOR<> prec;
+  prec.initialize(mass_matrix, 1.2);
+                                  // solve
+  cg.solve (mass_matrix, vec, tmp, prec);
+  
+                                  // distribute solution
+  constraints.distribute (vec);
+}
+
+
+int create_tria (unsigned int elm, Triangulation<2> &tria)
+{
+  std::vector<Point<2> > points_glob;
+  std::vector<Point<2> > points;
+
+  points_glob.push_back (Point<2> (0.0, 0.0));
+  points_glob.push_back (Point<2> (1.0, 0.0));
+  points_glob.push_back (Point<2> (1.0, 0.5));
+  points_glob.push_back (Point<2> (1.0, 1.0));
+  points_glob.push_back (Point<2> (0.6, 0.5));
+  points_glob.push_back (Point<2> (0.5, 1.0));
+  points_glob.push_back (Point<2> (0.0, 1.0));
+
+                                  // Prepare cell data
+  std::vector<CellData<2> > cells (1);
+
+  switch (elm)
+    {
+      case 0:
+           cells[0].vertices[0] = 0;
+           cells[0].vertices[1] = 1;
+           cells[0].vertices[2] = 4;
+           cells[0].vertices[3] = 2;
+           cells[0].material_id = 0;
+           break;
+      case 1:
+           cells[0].vertices[0] = 4;
+           cells[0].vertices[1] = 2;
+           cells[0].vertices[2] = 5;
+           cells[0].vertices[3] = 3;
+           cells[0].material_id = 0;
+           break;
+      case 2:
+           cells[0].vertices[0] = 0;
+           cells[0].vertices[1] = 4;
+           cells[0].vertices[2] = 6;
+           cells[0].vertices[3] = 5;
+           cells[0].material_id = 0;
+           break;
+    }
+
+  for (unsigned int i = 0; i < 4; ++i)
+    {
+      points.push_back (points_glob[cells[0].vertices[i]]);
+      cells[0].vertices[i] = i;
+    }
+
+  tria.create_triangulation (points, cells, SubCellData());
+
+  return (0);
+}
+
+
+int plot_shapes (DoFHandler<2> &dof_handler)
+{
+  Vector<double> solution (dof_handler.n_dofs ());
+  std::set<unsigned int> face_dofs;
+  
+                                  // Create set of all DoFs, which are on the boundary.
+  DoFHandler<2>::active_face_iterator face = dof_handler.begin_active_face(),
+                                     endf = dof_handler.end_face();
+  std::vector<unsigned int> face_dof_indices (dof_handler.get_fe().dofs_per_face);
+  for (; face!=endf; ++face)
+    {
+      face->get_dof_indices (face_dof_indices);
+      for (unsigned int i=0; i<dof_handler.get_fe().dofs_per_face; ++i)
+       face_dofs.insert (face_dof_indices[i]);
+    }
+  
+                                  // Now set a 1 at the place of the different DoFs and
+                                  // output the solution.
+  std::set<unsigned int>::iterator face_dof_iter;
+  for (face_dof_iter = face_dofs.begin ();
+       face_dof_iter != face_dofs.end (); ++face_dof_iter)
+    {
+      unsigned int dof = *face_dof_iter;
+
+      deallog << dof << std::endl;
+
+      solution (dof) = 1.0;
+
+                                      // Test the core functionality2
+      DataOut<2> *data_out = new DataOut<2>;
+      data_out->attach_dof_handler (dof_handler);
+      data_out->add_data_vector (solution, "solution");
+      data_out->build_patches (4, 0);
+
+      data_out->write_gnuplot (deallog.get_file_stream());
+
+      delete data_out;
+
+      solution (dof) = 0.0;
+    }
+
+  return (0);
+}
+
+
+int main (int /*argc*/, char **/*argv*/)
+{
+  std::ofstream logfile (logname);
+  logfile.precision (PRECISION);
+  logfile.setf(std::ios::fixed);  
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test ();
+}
+
+
+
+void check (const FiniteElement<2> &fe)
+{
+  Triangulation<2> tria_test;
+  DoFHandler<2> *dof_handler;
+
+  deallog << "Dofs/cell " << fe.dofs_per_cell
+         << "Dofs/face " << fe.dofs_per_face << std::endl;
+
+  for (unsigned int elm = 0; elm < 3; ++elm)
+    {
+      create_tria (elm, tria_test);
+
+                                      // Create a DoFHandler
+      dof_handler = new DoFHandler<2> (tria_test);
+      dof_handler->distribute_dofs (fe);
+
+      deallog << "Dofs total " << dof_handler->n_dofs () << std::endl;
+
+                                      // Alloc some DoFs
+      Vector<double> solution;
+      solution.reinit (dof_handler->n_dofs ());
+
+                                      // Project solution onto FE field
+      ConstraintMatrix     hn_constraints;
+      hn_constraints.clear ();
+      DoFTools::make_hanging_node_constraints (*dof_handler, 
+                                              hn_constraints);
+      hn_constraints.close ();
+
+      MappingQ1<2> map_default;
+      
+      project (map_default, *dof_handler, hn_constraints,
+              QGauss6<2> (), TestMap1<2>(2),
+              solution);
+
+                                      // Test the core functionality
+      DataOut<2> *data_out = new DataOut<2>;
+      data_out->attach_dof_handler (*dof_handler);
+      data_out->add_data_vector (solution, "solution");
+      data_out->build_patches (4);
+
+      data_out->write_gnuplot (deallog.get_file_stream());
+
+                                      // Now write face DoFs ..
+      for (unsigned int i = 0; i < 4; ++i)
+       {
+         deallog << solution(i) << std::endl;
+       }
+      
+                                      // Clean up ...
+      delete data_out;
+      delete (dof_handler);
+      tria_test.clear ();
+    }
+}
index 229dcb44513523723fff799cb4dcbcbd2dcb14b4..107787871f4f9122936ddf11925ec04015673f93 100644 (file)
@@ -1,5 +1,5 @@
-//----------------------------  fe_sys_projection_01.cc  ---------------------------
-//    fe_sys_projection_01.cc,v 1.3 2003/06/09 16:00:38 wolf Exp
+//----------------------------  abfprojection_01.cc  ---------------------------
+//    abfprojection_01.cc,v 1.3 2003/06/09 16:00:38 wolf Exp
 //    Version: 
 //
 //    Copyright (C) 2003, 2005, 2006 by the deal.II authors
@@ -9,7 +9,7 @@
 //    to the file deal.II/doc/license.html for the  text  and
 //    further information on this license.
 //
-//----------------------------  fe_sys_projection_01.cc  ---------------------------
+//----------------------------  abfprojection_01.cc  ---------------------------
 
 /*
  * Project the function [1,1] onto a deformed grid and see whether the
 
 
 
-#include "../tests.h"
-#include <base/logstream.h>
+char logname[] = "abf_projection_01/output";
+#include "deformed_projection.h"
 
-#define PRECISION 2
 
-#include <fstream>
-
-std::ofstream logfile ("fe_sys_projection_01/output");
-
-char buf[1000];
-#include <fstream>
-
-#include <grid/grid_generator.h>
-#include <grid/grid_in.h>
-#include <grid/grid_out.h>
-#include <grid/tria_boundary_lib.h>
-
-#include <fe/mapping_q.h>
-#include <fe/mapping_q1_eulerian.h>
-
-#include <numerics/vectors.h>
-#include <numerics/matrices.h>
-#include <numerics/data_out.h>
-#include <lac/vector.h>
-#include <lac/solver_cg.h>
-#include <lac/precondition.h>
-#include <lac/vector_memory.h>
-#include <lac/sparsity_pattern.h>
-#include <lac/sparse_matrix.h>
-
-#include <base/function.h>
-#include <base/quadrature_lib.h>
-#include <dofs/dof_constraints.h>
-#include <dofs/dof_tools.h>
-
-#include <fe/fe_system.h>
-#include <fe/fe.h>
-#include <fe/fe_q.h>
-#include <fe/fe_values.h>
-#include <fe/fe_raviart_thomas.h>
-#include <fe/fe_abf.h>
-#include <fe/fe_dgq.h>
-
-
-template <int dim>
-class TestMap1 : public Function<dim>
-{
-  public:
-    TestMap1 (const unsigned int n_components) :
-                   Function<dim> (n_components)
-      {};
-  
-    virtual ~TestMap1 () {};
-  
-    virtual double value (const Point<dim>   &p,
-                         const unsigned int  component = 0) const;
-  
-    void vector_value (const Point<dim> &p,
-                      Vector<double>   &return_value) const;
-};
-
-
-
-template <int dim>
-double
-TestMap1<dim>::value (const Point<dim>   &,
-                     const unsigned int  ) const 
-{
-  return (1.0);
-}
-
-
-template <int dim>
-void TestMap1<dim>::vector_value (const Point<dim> &p,
-                                 Vector<double>   &return_value) const
-{
-  Assert (return_value.size() == this->n_components,
-         ExcDimensionMismatch (return_value.size(), this->n_components));
-  
-                                  // Parabolic inflow profile
-  for (unsigned int iCount = 0; iCount < this->n_components; iCount++)
-    return_value (iCount) = value (p, iCount);
-}
-
-
-/*
- * Check the value of the derivative field.
- */
-
-void EvaluateDerivative (DoFHandler<2> *dof_handler,
-                        Vector<double> &solution)
-{
-                                  // This quadrature rule determines the points, where the
-                                  // derivative will be evaluated.
-  QGauss<2> quad (3);
-  FEValues<2> fe_values (dof_handler->get_fe (), quad, 
-                        UpdateFlags(update_values    |
-                                    update_q_points  |
-                                    update_gradients |
-                                    update_JxW_values));
-
-  const unsigned int   n_q_points    = quad.n_quadrature_points;
-  const unsigned int   n_components   = dof_handler->get_fe().n_components();
-  const unsigned int dofs_per_cell = dof_handler->get_fe().dofs_per_cell;
-
-                                  // Cell iterators
-  DoFHandler<2>::active_cell_iterator cell = dof_handler->begin_active(),
-                                     endc = dof_handler->end();
-
-  double err_l2 = 0,
-       err_hdiv = 0;
-
-  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-
-  for (; cell!=endc; ++cell)
-    {
-      cell->get_dof_indices (local_dof_indices);
-    
-      fe_values.reinit (cell);
-
-                                      // Get function values
-      std::vector<Vector<double> > this_value(n_q_points,
-                                             Vector<double>(n_components));
-      fe_values.get_function_values (solution, this_value);
-
-                                      // Get values from solution vector (For Trap.Rule)
-      std::vector<std::vector<Tensor<1,2> > >
-       grads_here (n_q_points, std::vector<Tensor<1,2> > (n_components));
-      fe_values.get_function_grads (solution, grads_here);
-
-      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-       {
-                                          //     double u0 = this_value[q_point](0);
-                                          //double v0 = this_value[q_point](1);
-
-         double u0 = 0;
-         double v0 = 0;
-         for (unsigned int i = 0; i < dofs_per_cell; ++i)
-           {
-             u0 += (solution (local_dof_indices[i]) *
-                    fe_values.shape_value_component(i, q_point, 0));
-             v0 += (solution (local_dof_indices[i]) *
-                    fe_values.shape_value_component(i, q_point, 1));
-           }
-         u0 -= 1.0;
-         v0 -= 1.0;
-
-         double dudx = grads_here[q_point][0][0];
-         double dvdy = grads_here[q_point][1][1];
-
-                                          /*
-                                            sprintf (buf, "QP %i %e %e %e %e %e %e Div %e\n", q_point, u0, v0,
-                                            dudx, dudy, dvdx, dvdy, dudx + dvdy);
-                                          */
-
-         err_l2 += (u0 * u0 + v0 * v0) * fe_values.JxW (q_point);
-         err_hdiv += (dudx + dvdy) * (dudx + dvdy) * fe_values.JxW (q_point);
-       }
-    }
-
-  sprintf (buf, "L2-Err %e  Hdiv-Err %e\n", pow (err_l2, 0.5),
-          pow (err_hdiv, 0.5));
-  deallog << buf;
-}
-
-
-template <int dim>
-void create_mass_matrix (const Mapping<dim>       &mapping,
-                        const DoFHandler<dim>    &dof,
-                        const Quadrature<dim>    &q,
-                        SparseMatrix<double>     &matrix,
-                        const Function<dim>   &rhs_function,
-                        Vector<double>        &rhs_vector,
-                        const Function<dim> * const coefficient = 0)
+void test ()
 {
-  UpdateFlags update_flags = UpdateFlags(update_values | update_JxW_values | update_q_points);
-  if (coefficient != 0)
-    update_flags = UpdateFlags (update_flags | update_q_points);
-  
-  FEValues<dim> fe_values (mapping, dof.get_fe(), q, update_flags);
-    
-  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
-                    n_q_points    = fe_values.n_quadrature_points;
-  const FiniteElement<dim>    &fe  = fe_values.get_fe();
-  const unsigned int n_components  = fe.n_components();
-
-  Assert(coefficient == 0 ||
-        coefficient->n_components==1 ||
-        coefficient->n_components==n_components, ExcInternalError());
-  
-  FullMatrix<double>  cell_matrix (dofs_per_cell, dofs_per_cell);
-  Vector<double> cell_vector (dofs_per_cell);
-  std::vector<double> coefficient_values (n_q_points);
-  std::vector<Vector<double> > coefficient_vector_values (n_q_points,
-                                                         Vector<double> (n_components));
-  
-  std::vector<unsigned int> dof_indices (dofs_per_cell);
-
-  std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
-  
-  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active (),
-                                                endc = dof.end ();
-  for (; cell!=endc; ++cell)
-    {
-      fe_values.reinit (cell);
-      
-      cell_matrix = 0;
-      cell->get_dof_indices (dof_indices);
-
-      const std::vector<double> &weights   = fe_values.get_JxW_values ();
-      rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
-      cell_vector = 0;
-      
-      if (coefficient != 0)
-       {
-         if (coefficient->n_components==1)
-           {
-             coefficient->value_list (fe_values.get_quadrature_points(),
-                                      coefficient_values);
-             for (unsigned int point=0; point<n_q_points; ++point)
-               {
-                 const double weight = fe_values.JxW(point);
-                 for (unsigned int i=0; i<dofs_per_cell; ++i)
-                   {
-                     const double v = fe_values.shape_value(i,point);
-                     for (unsigned int j=0; j<dofs_per_cell; ++j)
-                       {
-                         const double u = fe_values.shape_value(j,point);
-                         
-                         if ((n_components==1) ||
-                             (fe.system_to_component_index(i).first ==
-                              fe.system_to_component_index(j).first))
-                           cell_matrix(i,j) +=
-                             (u * v * weight * coefficient_values[point]);
-                       }
-                   }
-               }
-           }
-         else
-           {
-             coefficient->vector_value_list (fe_values.get_quadrature_points(),
-                                             coefficient_vector_values);
-             for (unsigned int point=0; point<n_q_points; ++point)
-               {
-                 const double weight = fe_values.JxW(point);
-                 for (unsigned int i=0; i<dofs_per_cell; ++i)
-                   {
-                     const double v = fe_values.shape_value(i,point);
-                     const unsigned int component_i=
-                       fe.system_to_component_index(i).first;
-                     for (unsigned int j=0; j<dofs_per_cell; ++j)
-                       {
-                         const double u = fe_values.shape_value(j,point);
-                         if ((n_components==1) ||
-                             (fe.system_to_component_index(j).first == component_i))
-                           cell_matrix(i,j) +=
-                             (u * v * weight *
-                              coefficient_vector_values[point](component_i));
-                       }
-                   }
-               }
-           }
-       }
-      else
-       {
-                                          // Compute eventual sign changes depending on the neighborhood
-                                          // between two faces.
-         std::vector<double> sign_change (dofs_per_cell, 1.0);
-         const unsigned int dofs_per_face = fe.dofs_per_face;
-         std::vector<unsigned int> face_dof_indices (dofs_per_face);
-         
-         for (unsigned int f = 0; f < 2; ++f)
-           {
-             typename DoFHandler<dim>::active_face_iterator face = cell->face (f);
-             if (!face->at_boundary ())
-               {
-                 unsigned int nn = cell->neighbor_of_neighbor (f);
-                 sprintf (buf, "Face %i  NeigNeig %i\n", f, nn);
-                 deallog << buf;
-                 if (nn > 1)
-                   {
-                     face->get_dof_indices (face_dof_indices);
-                     for (unsigned int j = 0; j < dofs_per_face; ++j)
-                       {
-                         sign_change[face_dof_indices[j]] = -1.0;
-                         sprintf (buf, "DoF %i\n", face_dof_indices[j]);
-                         deallog << buf;
-                       }
-                   }
-               }
-           }
-         
-         for (unsigned int point=0; point<n_q_points; ++point)
-           {
-             const double weight = fe_values.JxW(point);
-                                              //           const double weight = q.weight(point);
-             
-             std::vector<Vector<double> > val_vector (dofs_per_cell,
-                                                      Vector<double> (n_components));
-             
-                                              // Precompute the component values
-             for (unsigned int i=0; i < dofs_per_cell; ++i)
-               for (unsigned int comp_i = 0; comp_i < fe.n_components (); 
-                    ++comp_i)
-                 {
-                   val_vector[i](comp_i) = sign_change[i] * 
-                                           fe_values.shape_value_component(i,point,comp_i);
-                 }
-                                              // Now eventually switch the sign of some of the ansatzfunctions.
-                                              // TODO
-             
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
-               for (unsigned int comp_i = 0; comp_i < fe.n_components (); 
-                    ++comp_i)
-                 if (fe.get_nonzero_components(i)[comp_i] == true)
-                   {
-                     const double v = val_vector[i](comp_i);
-                     for (unsigned int j=0; j<dofs_per_cell; ++j)
-                       for (unsigned int comp_j = 0;
-                            comp_j < fe.n_components (); ++comp_j)
-                         if (fe.get_nonzero_components(j)[comp_j] == true)
-                           {
-                             const double u = val_vector[j](comp_j);
-                             if ((n_components==1) ||
-                                 (comp_i == comp_j))
-                               cell_matrix(i,j) += (u * v * weight);
-                           }
-                   }
-             
-             
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
-               for (unsigned int comp_i = 0; comp_i < fe.n_components (); 
-                    ++comp_i)
-                 if (fe.get_nonzero_components(i)[comp_i] == true)
-                   {
-                     cell_vector(i) += rhs_values[point](comp_i) *
-                                       val_vector[i](comp_i) * weights[point];
-                   }
-           }
-       }
-                                      // transfer everything into the
-                                      // global object. lock the
-                                      // matrix meanwhile
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       for (unsigned int j=0; j<dofs_per_cell; ++j)
-         if ((n_components==1) ||
-             (cell_matrix (i,j) != 0.0))
-/*
-  (fe.system_to_component_index(i).first ==
-  fe.system_to_component_index(j).first))
-*/
-           matrix.add (dof_indices[i], dof_indices[j],
-                       cell_matrix(i,j));
-
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       rhs_vector(dof_indices[i]) += cell_vector(i);
-    }
-}
-
-
-template <int dim>
-void create_right_hand_side (const Mapping<dim>    &mapping,
-                            const DoFHandler<dim> &dof_handler,
-                            const Quadrature<dim> &quadrature,
-                            const Function<dim>   &rhs_function,
-                            Vector<double>        &rhs_vector)
-{
-  const FiniteElement<dim> &fe  = dof_handler.get_fe();
-  Assert (fe.n_components() == rhs_function.n_components,
-         ExcInternalError());
-  Assert (rhs_vector.size() == dof_handler.n_dofs(),
-         ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
-  rhs_vector = 0;
-  
-  UpdateFlags update_flags = UpdateFlags(update_values   |
-                                        update_q_points |
-                                        update_JxW_values);
-  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,
-                    n_components  = fe.n_components();
-  
-  std::vector<unsigned int> dofs (dofs_per_cell);
-  Vector<double> cell_vector (dofs_per_cell);
-
-  typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
-                                                endc = dof_handler.end();
-
-  if (n_components==1)
-    {
-      std::vector<double> rhs_values(n_q_points);
-      
-      for (; cell!=endc; ++cell) 
-       {
-         fe_values.reinit(cell);
-         
-         const std::vector<double> &weights   = fe_values.get_JxW_values ();
-         rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values);
-         
-         cell_vector = 0;
-         for (unsigned int point=0; point<n_q_points; ++point)
-           for (unsigned int i=0; i<dofs_per_cell; ++i) 
-             cell_vector(i) += rhs_values[point] *
-                               fe_values.shape_value(i,point) *
-                               weights[point];
-       
-         cell->get_dof_indices (dofs);
-         
-         for (unsigned int i=0; i<dofs_per_cell; ++i)
-           rhs_vector(dofs[i]) += cell_vector(i);
-       }
-    }
-  else
-    {
-      std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
-      
-      for (; cell!=endc; ++cell) 
-       {
-         fe_values.reinit(cell);
-         
-         const std::vector<double> &weights   = fe_values.get_JxW_values ();
-         rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
-         
-         cell_vector = 0;
-         for (unsigned int point=0; point<n_q_points; ++point)
-           for (unsigned int i=0; i<dofs_per_cell; ++i)
-             for (unsigned int comp_i = 0; comp_i < fe.n_components (); 
-                  ++comp_i)
-                                                //                 if (fe.get_nonzero_components(i)[comp_i] == true)
-               {
-                 cell_vector(i) += rhs_values[point](comp_i) *
-                                   fe_values.shape_value_component(i,point,comp_i) *
-                                   weights[point];
-               }
-         
-         cell->get_dof_indices (dofs);
-         
-         for (unsigned int i=0; i<dofs_per_cell; ++i)
-           rhs_vector(dofs[i]) += cell_vector(i);
-       }
-    }
-}
-
-
-
-//
-// This function replaces the deal.II implementation of the projection.
-// The purpose is to have more freedom in assembling the matrix.
-//
-
-template <int dim>
-void project (const Mapping<dim>       &mapping,
-             const DoFHandler<dim>    &dof,
-             const ConstraintMatrix   &constraints,
-             const Quadrature<dim>    &quadrature,
-             const Function<dim>      &function,
-             Vector<double>           &vec,
-             const bool                enforce_zero_boundary = false,
-             const Quadrature<dim-1>  & = QGauss2<dim-1>(),
-             const bool                project_to_boundary_first = false)
-{
-  Assert (dof.get_fe().n_components() == function.n_components,
-         ExcInternalError());
-  
-  const FiniteElement<dim> &fe = dof.get_fe();
-
-                                  // make up boundary values
-  std::map<unsigned int,double> boundary_values;
-
-  if (enforce_zero_boundary == true) 
-                                    // no need to project boundary
-                                    // values, but enforce
-                                    // homogeneous boundary values
-                                    // anyway
-    {
-                                      // loop over all boundary faces
-                                      // to get all dof indices of
-                                      // dofs on the boundary. note
-                                      // that in 3d there are cases
-                                      // where a face is not at the
-                                      // boundary, yet one of its
-                                      // lines is, and we should
-                                      // consider the degrees of
-                                      // freedom on it as boundary
-                                      // nodes. likewise, in 2d and
-                                      // 3d there are cases where a
-                                      // cell is only at the boundary
-                                      // by one vertex. nevertheless,
-                                      // since we do not support
-                                      // boundaries with dimension
-                                      // less or equal to dim-2, each
-                                      // such boundary dof is also
-                                      // found from some other face
-                                      // that is actually wholly on
-                                      // the boundary, not only by
-                                      // one line or one vertex
-      typename DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
-                                                    endf = dof.end_face();
-      std::vector<unsigned int> face_dof_indices (fe.dofs_per_face);
-      for (; face!=endf; ++face)
-       if (face->at_boundary())
-         {
-           face->get_dof_indices (face_dof_indices);
-           for (unsigned int i=0; i<fe.dofs_per_face; ++i)
-                                              // enter zero boundary values
-                                              // for all boundary nodes
-                                              //
-                                              // we need not care about
-                                              // vector valued elements here,
-                                              // since we set all components
-             boundary_values[face_dof_indices[i]] = 0.;
-         };
-    }
-  else
-                                    // no homogeneous boundary values
-    if (project_to_boundary_first == true)
-                                      // boundary projection required
-      {
-/*
-                                // set up a list of boundary functions for
-                                                                 // the different boundary parts. We want the
-                                                                                                  // @p{function} to hold on all parts of the
-                                                                                                                                   // boundary
-                                                                                                                                   typename FunctionMap<dim>::type boundary_functions;
-                                                                                                                                   for (unsigned char c=0; c<255; ++c)
-                                                                                                                                   boundary_functions[c] = &function;
-                                                                                                                                   project_boundary_values (dof, boundary_functions, q_boundary,
-                                                                                                                                   boundary_values);
-*/
-      }
-
-
-                                  // set up mass matrix and right hand side
-  vec.reinit (dof.n_dofs());
-  SparsityPattern sparsity(dof.n_dofs(),
-                          dof.n_dofs(),
-                          dof.max_couplings_between_dofs());
-  DoFTools::make_sparsity_pattern (dof, sparsity);
-  constraints.condense (sparsity);
-  
-  SparseMatrix<double> mass_matrix (sparsity);
-  Vector<double> tmp (mass_matrix.n());
-
-  create_mass_matrix (mapping, dof, quadrature, mass_matrix, function, tmp);
-  sprintf (buf, "MM created\n");
-  deallog << buf;
-  create_right_hand_side (mapping, dof, quadrature, function, tmp);
-  sprintf (buf, "RHS created\n");
-  deallog << buf;
-
-  constraints.condense (mass_matrix);
-  constraints.condense (tmp);
-  if (boundary_values.size() != 0)
-    MatrixTools::apply_boundary_values (boundary_values,
-                                       mass_matrix, vec, tmp,
-                                       true);
-
-  SolverControl           control(1000,1e-16);
-  PrimitiveVectorMemory<> memory;
-  SolverCG<>              cg(control,memory);
-
-  PreconditionSSOR<> prec;
-  prec.initialize(mass_matrix, 1.2);
-                                  // solve
-  cg.solve (mass_matrix, vec, tmp, prec);
-  
-                                  // distribute solution
-  constraints.distribute (vec);
-}
-
-
-int create_tria (unsigned int elm, Triangulation<2> &tria)
-{
-  std::vector<Point<2> > points_glob;
-  std::vector<Point<2> > points;
-
-  points_glob.push_back (Point<2> (0.0, 0.0));
-  points_glob.push_back (Point<2> (1.0, 0.0));
-  points_glob.push_back (Point<2> (1.0, 0.5));
-  points_glob.push_back (Point<2> (1.0, 1.0));
-  points_glob.push_back (Point<2> (0.6, 0.5));
-  points_glob.push_back (Point<2> (0.5, 1.0));
-  points_glob.push_back (Point<2> (0.0, 1.0));
-
-                                  // Prepare cell data
-  std::vector<CellData<2> > cells (1);
-
-  switch (elm)
-    {
-      case 0:
-           cells[0].vertices[0] = 0;
-           cells[0].vertices[1] = 1;
-           cells[0].vertices[2] = 4;
-           cells[0].vertices[3] = 2;
-           cells[0].material_id = 0;
-           break;
-      case 1:
-           cells[0].vertices[0] = 4;
-           cells[0].vertices[1] = 2;
-           cells[0].vertices[2] = 5;
-           cells[0].vertices[3] = 3;
-           cells[0].material_id = 0;
-           break;
-      case 2:
-           cells[0].vertices[0] = 0;
-           cells[0].vertices[1] = 4;
-           cells[0].vertices[2] = 6;
-           cells[0].vertices[3] = 5;
-           cells[0].material_id = 0;
-           break;
-    }
-
-  for (unsigned int i = 0; i < 4; ++i)
-    {
-      points.push_back (points_glob[cells[0].vertices[i]]);
-      cells[0].vertices[i] = i;
-    }
-
-  tria.create_triangulation (points, cells, SubCellData());
-
-  return (0);
-}
-
-
-int plot_shapes (DoFHandler<2> &dof_handler)
-{
-  Vector<double> solution (dof_handler.n_dofs ());
-  std::set<unsigned int> face_dofs;
-  
-                                  // Create set of all DoFs, which are on the boundary.
-  DoFHandler<2>::active_face_iterator face = dof_handler.begin_active_face(),
-                                     endf = dof_handler.end_face();
-  std::vector<unsigned int> face_dof_indices (dof_handler.get_fe().dofs_per_face);
-  for (; face!=endf; ++face)
-    {
-      face->get_dof_indices (face_dof_indices);
-      for (unsigned int i=0; i<dof_handler.get_fe().dofs_per_face; ++i)
-       face_dofs.insert (face_dof_indices[i]);
-    }
-  
-                                  // Now set a 1 at the place of the different DoFs and
-                                  // output the solution.
-  std::set<unsigned int>::iterator face_dof_iter;
-  for (face_dof_iter = face_dofs.begin ();
-       face_dof_iter != face_dofs.end (); ++face_dof_iter)
-    {
-      unsigned int dof = *face_dof_iter;
-
-      sprintf (buf, "%i\n", dof);
-      deallog << buf;
-
-      solution (dof) = 1.0;
-
-                                      // Test the core functionality2
-      DataOut<2> *data_out = new DataOut<2>;
-      data_out->attach_dof_handler (dof_handler);
-      data_out->add_data_vector (solution, "solution");
-      data_out->build_patches (4, 0);
-
-      data_out->write_gnuplot (deallog.get_file_stream());
-
-      delete data_out;
-
-      solution (dof) = 0.0;
-    }
-
-  return (0);
-}
-
-
-int main (int /*argc*/, char **/*argv*/)
-{
-  logfile.precision (PRECISION);
-  logfile.setf(std::ios::fixed);  
-  deallog.attach(logfile);
-  deallog.depth_console(0);
-  deallog.threshold_double(1.e-10);
-
-  Triangulation<2> tria_test;
-  DoFHandler<2> *dof_handler;
-
   FESystem<2> fe (FE_Q<2>(3), 2);
-  deallog << "Dofs/cell " << fe.dofs_per_cell
-         << "Dofs/face " << fe.dofs_per_face << std::endl;
-
-  for (unsigned int elm = 0; elm < 3; ++elm)
-    {
-      create_tria (elm, tria_test);
-
-                                      // Create a DoFHandler
-      dof_handler = new DoFHandler<2> (tria_test);
-      dof_handler->distribute_dofs (fe);
-
-      deallog << "Dofs total " << dof_handler->n_dofs () << std::endl;
-
-                                      // Alloc some DoFs
-      Vector<double> solution;
-      solution.reinit (dof_handler->n_dofs ());
-
-                                      // Project solution onto FE field
-      ConstraintMatrix     hn_constraints;
-      hn_constraints.clear ();
-      DoFTools::make_hanging_node_constraints (*dof_handler, 
-                                              hn_constraints);
-      hn_constraints.close ();
-
-      MappingQ1<2> map_default;
-      
-      project (map_default, *dof_handler, hn_constraints,
-              QGauss6<2> (), TestMap1<2>(2),
-              solution);
-
-                                      // Test the core functionality
-      DataOut<2> *data_out = new DataOut<2>;
-      data_out->attach_dof_handler (*dof_handler);
-      data_out->add_data_vector (solution, "solution");
-      data_out->build_patches (4);
-
-      data_out->write_gnuplot (deallog.get_file_stream());
-
-                                      // Now write face DoFs ..
-      for (unsigned int i = 0; i < 4; ++i)
-       {
-         sprintf (buf, "%i %e\n", i, solution (i));
-         deallog << buf;
-       }
-      
-                                      // Clean up ...
-      delete data_out;
-      delete (dof_handler);
-      tria_test.clear ();
-    }
-
-  return (0);
+  check (fe);
 }
index a28673a88d2d37421cfc9dc6e4eed1936ff9f236..22d0e320e7fb1c97d1477183b0f286331c383f16 100644 (file)
 
 
 
-#include "../tests.h"
-#include <base/logstream.h>
+char logname[] = "rt_projection_01/output";
+#include "deformed_projection.h"
 
-#define PRECISION 2
 
-#include <fstream>
-
-std::ofstream logfile ("rt_projection_01/output");
-
-char buf[1000];
-#include <fstream>
-
-#include <grid/grid_generator.h>
-#include <grid/grid_in.h>
-#include <grid/grid_out.h>
-#include <grid/tria_boundary_lib.h>
-
-#include <fe/mapping_q.h>
-#include <fe/mapping_q1_eulerian.h>
-
-#include <numerics/vectors.h>
-#include <numerics/matrices.h>
-#include <numerics/data_out.h>
-#include <lac/vector.h>
-#include <lac/solver_cg.h>
-#include <lac/precondition.h>
-#include <lac/vector_memory.h>
-#include <lac/sparsity_pattern.h>
-#include <lac/sparse_matrix.h>
-
-#include <base/function.h>
-#include <base/quadrature_lib.h>
-#include <dofs/dof_constraints.h>
-#include <dofs/dof_tools.h>
-
-#include <fe/fe_system.h>
-#include <fe/fe.h>
-#include <fe/fe_q.h>
-#include <fe/fe_values.h>
-#include <fe/fe_raviart_thomas.h>
-#include <fe/fe_abf.h>
-#include <fe/fe_dgq.h>
-
-
-template <int dim>
-class TestMap1 : public Function<dim>
-{
-  public:
-    TestMap1 (const unsigned int n_components) :
-                   Function<dim> (n_components)
-      {};
-  
-    virtual ~TestMap1 () {};
-  
-    virtual double value (const Point<dim>   &p,
-                         const unsigned int  component = 0) const;
-  
-    void vector_value (const Point<dim> &p,
-                      Vector<double>   &return_value) const;
-};
-
-
-
-template <int dim>
-double
-TestMap1<dim>::value (const Point<dim>   &,
-                     const unsigned int  ) const 
-{
-  return (1.0);
-}
-
-
-template <int dim>
-void TestMap1<dim>::vector_value (const Point<dim> &p,
-                                 Vector<double>   &return_value) const
-{
-  Assert (return_value.size() == this->n_components,
-         ExcDimensionMismatch (return_value.size(), this->n_components));
-  
-                                  // Parabolic inflow profile
-  for (unsigned int iCount = 0; iCount < this->n_components; iCount++)
-    return_value (iCount) = value (p, iCount);
-}
-
-
-/*
- * Check the value of the derivative field.
- */
-
-void EvaluateDerivative (DoFHandler<2> *dof_handler,
-                        Vector<double> &solution)
-{
-                                  // This quadrature rule determines the points, where the
-                                  // derivative will be evaluated.
-  QGauss<2> quad (3);
-  FEValues<2> fe_values (dof_handler->get_fe (), quad, 
-                        UpdateFlags(update_values    |
-                                    update_q_points  |
-                                    update_gradients |
-                                    update_JxW_values));
-
-  const unsigned int   n_q_points    = quad.n_quadrature_points;
-  const unsigned int   n_components   = dof_handler->get_fe().n_components();
-  const unsigned int dofs_per_cell = dof_handler->get_fe().dofs_per_cell;
-
-                                  // Cell iterators
-  DoFHandler<2>::active_cell_iterator cell = dof_handler->begin_active(),
-                                     endc = dof_handler->end();
-
-  double err_l2 = 0,
-       err_hdiv = 0;
-
-  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-
-  for (; cell!=endc; ++cell)
-    {
-      cell->get_dof_indices (local_dof_indices);
-    
-      fe_values.reinit (cell);
-
-                                      // Get function values
-      std::vector<Vector<double> > this_value(n_q_points,
-                                             Vector<double>(n_components));
-      fe_values.get_function_values (solution, this_value);
-
-                                      // Get values from solution vector (For Trap.Rule)
-      std::vector<std::vector<Tensor<1,2> > >
-       grads_here (n_q_points, std::vector<Tensor<1,2> > (n_components));
-      fe_values.get_function_grads (solution, grads_here);
-
-      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-       {
-                                          //     double u0 = this_value[q_point](0);
-                                          //double v0 = this_value[q_point](1);
-
-         double u0 = 0;
-         double v0 = 0;
-         for (unsigned int i = 0; i < dofs_per_cell; ++i)
-           {
-             u0 += (solution (local_dof_indices[i]) *
-                    fe_values.shape_value_component(i, q_point, 0));
-             v0 += (solution (local_dof_indices[i]) *
-                    fe_values.shape_value_component(i, q_point, 1));
-           }
-         u0 -= 1.0;
-         v0 -= 1.0;
-
-         double dudx = grads_here[q_point][0][0];
-         double dvdy = grads_here[q_point][1][1];
-
-                                          /*
-                                            sprintf (buf, "QP %i %e %e %e %e %e %e Div %e\n", q_point, u0, v0,
-                                            dudx, dudy, dvdx, dvdy, dudx + dvdy);
-                                          */
-
-         err_l2 += (u0 * u0 + v0 * v0) * fe_values.JxW (q_point);
-         err_hdiv += (dudx + dvdy) * (dudx + dvdy) * fe_values.JxW (q_point);
-       }
-    }
-
-  sprintf (buf, "L2-Err %e  Hdiv-Err %e\n", pow (err_l2, 0.5),
-          pow (err_hdiv, 0.5));
-  deallog << buf;
-}
-
-
-template <int dim>
-void create_mass_matrix (const Mapping<dim>       &mapping,
-                        const DoFHandler<dim>    &dof,
-                        const Quadrature<dim>    &q,
-                        SparseMatrix<double>     &matrix,
-                        const Function<dim>   &rhs_function,
-                        Vector<double>        &rhs_vector,
-                        const Function<dim> * const coefficient = 0)
+void test ()
 {
-  UpdateFlags update_flags = UpdateFlags(update_values | update_JxW_values | update_q_points);
-  if (coefficient != 0)
-    update_flags = UpdateFlags (update_flags | update_q_points);
-  
-  FEValues<dim> fe_values (mapping, dof.get_fe(), q, update_flags);
-    
-  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
-                    n_q_points    = fe_values.n_quadrature_points;
-  const FiniteElement<dim>    &fe  = fe_values.get_fe();
-  const unsigned int n_components  = fe.n_components();
-
-  Assert(coefficient == 0 ||
-        coefficient->n_components==1 ||
-        coefficient->n_components==n_components, ExcInternalError());
-  
-  FullMatrix<double>  cell_matrix (dofs_per_cell, dofs_per_cell);
-  Vector<double> cell_vector (dofs_per_cell);
-  std::vector<double> coefficient_values (n_q_points);
-  std::vector<Vector<double> > coefficient_vector_values (n_q_points,
-                                                         Vector<double> (n_components));
-  
-  std::vector<unsigned int> dof_indices (dofs_per_cell);
-
-  std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
-  
-  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active (),
-                                                endc = dof.end ();
-  for (; cell!=endc; ++cell)
-    {
-      fe_values.reinit (cell);
-      
-      cell_matrix = 0;
-      cell->get_dof_indices (dof_indices);
-
-      const std::vector<double> &weights   = fe_values.get_JxW_values ();
-      rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
-      cell_vector = 0;
-      
-      if (coefficient != 0)
-       {
-         if (coefficient->n_components==1)
-           {
-             coefficient->value_list (fe_values.get_quadrature_points(),
-                                      coefficient_values);
-             for (unsigned int point=0; point<n_q_points; ++point)
-               {
-                 const double weight = fe_values.JxW(point);
-                 for (unsigned int i=0; i<dofs_per_cell; ++i)
-                   {
-                     const double v = fe_values.shape_value(i,point);
-                     for (unsigned int j=0; j<dofs_per_cell; ++j)
-                       {
-                         const double u = fe_values.shape_value(j,point);
-                         
-                         if ((n_components==1) ||
-                             (fe.system_to_component_index(i).first ==
-                              fe.system_to_component_index(j).first))
-                           cell_matrix(i,j) +=
-                             (u * v * weight * coefficient_values[point]);
-                       }
-                   }
-               }
-           }
-         else
-           {
-             coefficient->vector_value_list (fe_values.get_quadrature_points(),
-                                             coefficient_vector_values);
-             for (unsigned int point=0; point<n_q_points; ++point)
-               {
-                 const double weight = fe_values.JxW(point);
-                 for (unsigned int i=0; i<dofs_per_cell; ++i)
-                   {
-                     const double v = fe_values.shape_value(i,point);
-                     const unsigned int component_i=
-                       fe.system_to_component_index(i).first;
-                     for (unsigned int j=0; j<dofs_per_cell; ++j)
-                       {
-                         const double u = fe_values.shape_value(j,point);
-                         if ((n_components==1) ||
-                             (fe.system_to_component_index(j).first == component_i))
-                           cell_matrix(i,j) +=
-                             (u * v * weight *
-                              coefficient_vector_values[point](component_i));
-                       }
-                   }
-               }
-           }
-       }
-      else
-       {
-                                          // Compute eventual sign changes depending on the neighborhood
-                                          // between two faces.
-         std::vector<double> sign_change (dofs_per_cell, 1.0);
-         const unsigned int dofs_per_face = fe.dofs_per_face;
-         std::vector<unsigned int> face_dof_indices (dofs_per_face);
-         
-         for (unsigned int f = 0; f < 2; ++f)
-           {
-             typename DoFHandler<dim>::active_face_iterator face = cell->face (f);
-             if (!face->at_boundary ())
-               {
-                 unsigned int nn = cell->neighbor_of_neighbor (f);
-                 sprintf (buf, "Face %i  NeigNeig %i\n", f, nn);
-                 deallog << buf;
-                 if (nn > 1)
-                   {
-                     face->get_dof_indices (face_dof_indices);
-                     for (unsigned int j = 0; j < dofs_per_face; ++j)
-                       {
-                         sign_change[face_dof_indices[j]] = -1.0;
-                         sprintf (buf, "DoF %i\n", face_dof_indices[j]);
-                         deallog << buf;
-                       }
-                   }
-               }
-           }
-         
-         for (unsigned int point=0; point<n_q_points; ++point)
-           {
-             const double weight = fe_values.JxW(point);
-                                              //           const double weight = q.weight(point);
-             
-             std::vector<Vector<double> > val_vector (dofs_per_cell,
-                                                      Vector<double> (n_components));
-             
-                                              // Precompute the component values
-             for (unsigned int i=0; i < dofs_per_cell; ++i)
-               for (unsigned int comp_i = 0; comp_i < fe.n_components (); 
-                    ++comp_i)
-                 {
-                   val_vector[i](comp_i) = sign_change[i] * 
-                                           fe_values.shape_value_component(i,point,comp_i);
-                 }
-                                              // Now eventually switch the sign of some of the ansatzfunctions.
-                                              // TODO
-             
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
-               for (unsigned int comp_i = 0; comp_i < fe.n_components (); 
-                    ++comp_i)
-                 if (fe.get_nonzero_components(i)[comp_i] == true)
-                   {
-                     const double v = val_vector[i](comp_i);
-                     for (unsigned int j=0; j<dofs_per_cell; ++j)
-                       for (unsigned int comp_j = 0;
-                            comp_j < fe.n_components (); ++comp_j)
-                         if (fe.get_nonzero_components(j)[comp_j] == true)
-                           {
-                             const double u = val_vector[j](comp_j);
-                             if ((n_components==1) ||
-                                 (comp_i == comp_j))
-                               cell_matrix(i,j) += (u * v * weight);
-                           }
-                   }
-             
-             
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
-               for (unsigned int comp_i = 0; comp_i < fe.n_components (); 
-                    ++comp_i)
-                 if (fe.get_nonzero_components(i)[comp_i] == true)
-                   {
-                     cell_vector(i) += rhs_values[point](comp_i) *
-                                       val_vector[i](comp_i) * weights[point];
-                   }
-           }
-       }
-                                      // transfer everything into the
-                                      // global object. lock the
-                                      // matrix meanwhile
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       for (unsigned int j=0; j<dofs_per_cell; ++j)
-         if ((n_components==1) ||
-             (cell_matrix (i,j) != 0.0))
-/*
-  (fe.system_to_component_index(i).first ==
-  fe.system_to_component_index(j).first))
-*/
-           matrix.add (dof_indices[i], dof_indices[j],
-                       cell_matrix(i,j));
-
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       rhs_vector(dof_indices[i]) += cell_vector(i);
-    }
-}
-
-
-template <int dim>
-void create_right_hand_side (const Mapping<dim>    &mapping,
-                            const DoFHandler<dim> &dof_handler,
-                            const Quadrature<dim> &quadrature,
-                            const Function<dim>   &rhs_function,
-                            Vector<double>        &rhs_vector)
-{
-  const FiniteElement<dim> &fe  = dof_handler.get_fe();
-  Assert (fe.n_components() == rhs_function.n_components,
-         ExcInternalError());
-  Assert (rhs_vector.size() == dof_handler.n_dofs(),
-         ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
-  rhs_vector = 0;
-  
-  UpdateFlags update_flags = UpdateFlags(update_values   |
-                                        update_q_points |
-                                        update_JxW_values);
-  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,
-                    n_components  = fe.n_components();
-  
-  std::vector<unsigned int> dofs (dofs_per_cell);
-  Vector<double> cell_vector (dofs_per_cell);
-
-  typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
-                                                endc = dof_handler.end();
-
-  if (n_components==1)
-    {
-      std::vector<double> rhs_values(n_q_points);
-      
-      for (; cell!=endc; ++cell) 
-       {
-         fe_values.reinit(cell);
-         
-         const std::vector<double> &weights   = fe_values.get_JxW_values ();
-         rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values);
-         
-         cell_vector = 0;
-         for (unsigned int point=0; point<n_q_points; ++point)
-           for (unsigned int i=0; i<dofs_per_cell; ++i) 
-             cell_vector(i) += rhs_values[point] *
-                               fe_values.shape_value(i,point) *
-                               weights[point];
-       
-         cell->get_dof_indices (dofs);
-         
-         for (unsigned int i=0; i<dofs_per_cell; ++i)
-           rhs_vector(dofs[i]) += cell_vector(i);
-       }
-    }
-  else
-    {
-      std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
-      
-      for (; cell!=endc; ++cell) 
-       {
-         fe_values.reinit(cell);
-         
-         const std::vector<double> &weights   = fe_values.get_JxW_values ();
-         rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
-         
-         cell_vector = 0;
-         for (unsigned int point=0; point<n_q_points; ++point)
-           for (unsigned int i=0; i<dofs_per_cell; ++i)
-             for (unsigned int comp_i = 0; comp_i < fe.n_components (); 
-                  ++comp_i)
-                                                //                 if (fe.get_nonzero_components(i)[comp_i] == true)
-               {
-                 cell_vector(i) += rhs_values[point](comp_i) *
-                                   fe_values.shape_value_component(i,point,comp_i) *
-                                   weights[point];
-               }
-         
-         cell->get_dof_indices (dofs);
-         
-         for (unsigned int i=0; i<dofs_per_cell; ++i)
-           rhs_vector(dofs[i]) += cell_vector(i);
-       }
-    }
-}
-
-
-
-//
-// This function replaces the deal.II implementation of the projection.
-// The purpose is to have more freedom in assembling the matrix.
-//
-
-template <int dim>
-void project (const Mapping<dim>       &mapping,
-             const DoFHandler<dim>    &dof,
-             const ConstraintMatrix   &constraints,
-             const Quadrature<dim>    &quadrature,
-             const Function<dim>      &function,
-             Vector<double>           &vec,
-             const bool                enforce_zero_boundary = false,
-             const Quadrature<dim-1>  & = QGauss2<dim-1>(),
-             const bool                project_to_boundary_first = false)
-{
-  Assert (dof.get_fe().n_components() == function.n_components,
-         ExcInternalError());
-  
-  const FiniteElement<dim> &fe = dof.get_fe();
-
-                                  // make up boundary values
-  std::map<unsigned int,double> boundary_values;
-
-  if (enforce_zero_boundary == true) 
-                                    // no need to project boundary
-                                    // values, but enforce
-                                    // homogeneous boundary values
-                                    // anyway
-    {
-                                      // loop over all boundary faces
-                                      // to get all dof indices of
-                                      // dofs on the boundary. note
-                                      // that in 3d there are cases
-                                      // where a face is not at the
-                                      // boundary, yet one of its
-                                      // lines is, and we should
-                                      // consider the degrees of
-                                      // freedom on it as boundary
-                                      // nodes. likewise, in 2d and
-                                      // 3d there are cases where a
-                                      // cell is only at the boundary
-                                      // by one vertex. nevertheless,
-                                      // since we do not support
-                                      // boundaries with dimension
-                                      // less or equal to dim-2, each
-                                      // such boundary dof is also
-                                      // found from some other face
-                                      // that is actually wholly on
-                                      // the boundary, not only by
-                                      // one line or one vertex
-      typename DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
-                                                    endf = dof.end_face();
-      std::vector<unsigned int> face_dof_indices (fe.dofs_per_face);
-      for (; face!=endf; ++face)
-       if (face->at_boundary())
-         {
-           face->get_dof_indices (face_dof_indices);
-           for (unsigned int i=0; i<fe.dofs_per_face; ++i)
-                                              // enter zero boundary values
-                                              // for all boundary nodes
-                                              //
-                                              // we need not care about
-                                              // vector valued elements here,
-                                              // since we set all components
-             boundary_values[face_dof_indices[i]] = 0.;
-         };
-    }
-  else
-                                    // no homogeneous boundary values
-    if (project_to_boundary_first == true)
-                                      // boundary projection required
-      {
-/*
-                                // set up a list of boundary functions for
-                                                                 // the different boundary parts. We want the
-                                                                                                  // @p{function} to hold on all parts of the
-                                                                                                                                   // boundary
-                                                                                                                                   typename FunctionMap<dim>::type boundary_functions;
-                                                                                                                                   for (unsigned char c=0; c<255; ++c)
-                                                                                                                                   boundary_functions[c] = &function;
-                                                                                                                                   project_boundary_values (dof, boundary_functions, q_boundary,
-                                                                                                                                   boundary_values);
-*/
-      }
-
-
-                                  // set up mass matrix and right hand side
-  vec.reinit (dof.n_dofs());
-  SparsityPattern sparsity(dof.n_dofs(),
-                          dof.n_dofs(),
-                          dof.max_couplings_between_dofs());
-  DoFTools::make_sparsity_pattern (dof, sparsity);
-  constraints.condense (sparsity);
-  
-  SparseMatrix<double> mass_matrix (sparsity);
-  Vector<double> tmp (mass_matrix.n());
-
-  create_mass_matrix (mapping, dof, quadrature, mass_matrix, function, tmp);
-  sprintf (buf, "MM created\n");
-  deallog << buf;
-  create_right_hand_side (mapping, dof, quadrature, function, tmp);
-  sprintf (buf, "RHS created\n");
-  deallog << buf;
-
-  constraints.condense (mass_matrix);
-  constraints.condense (tmp);
-  if (boundary_values.size() != 0)
-    MatrixTools::apply_boundary_values (boundary_values,
-                                       mass_matrix, vec, tmp,
-                                       true);
-
-  SolverControl           control(1000,1e-16);
-  PrimitiveVectorMemory<> memory;
-  SolverCG<>              cg(control,memory);
-
-  PreconditionSSOR<> prec;
-  prec.initialize(mass_matrix, 1.2);
-                                  // solve
-  cg.solve (mass_matrix, vec, tmp, prec);
-  
-                                  // distribute solution
-  constraints.distribute (vec);
-}
-
-
-int create_tria (unsigned int elm, Triangulation<2> &tria)
-{
-  std::vector<Point<2> > points_glob;
-  std::vector<Point<2> > points;
-
-  points_glob.push_back (Point<2> (0.0, 0.0));
-  points_glob.push_back (Point<2> (1.0, 0.0));
-  points_glob.push_back (Point<2> (1.0, 0.5));
-  points_glob.push_back (Point<2> (1.0, 1.0));
-  points_glob.push_back (Point<2> (0.6, 0.5));
-  points_glob.push_back (Point<2> (0.5, 1.0));
-  points_glob.push_back (Point<2> (0.0, 1.0));
-
-                                  // Prepare cell data
-  std::vector<CellData<2> > cells (1);
-
-  switch (elm)
-    {
-      case 0:
-           cells[0].vertices[0] = 0;
-           cells[0].vertices[1] = 1;
-           cells[0].vertices[2] = 4;
-           cells[0].vertices[3] = 2;
-           cells[0].material_id = 0;
-           break;
-      case 1:
-           cells[0].vertices[0] = 4;
-           cells[0].vertices[1] = 2;
-           cells[0].vertices[2] = 5;
-           cells[0].vertices[3] = 3;
-           cells[0].material_id = 0;
-           break;
-      case 2:
-           cells[0].vertices[0] = 0;
-           cells[0].vertices[1] = 4;
-           cells[0].vertices[2] = 6;
-           cells[0].vertices[3] = 5;
-           cells[0].material_id = 0;
-           break;
-    }
-
-  for (unsigned int i = 0; i < 4; ++i)
-    {
-      points.push_back (points_glob[cells[0].vertices[i]]);
-      cells[0].vertices[i] = i;
-    }
-
-  tria.create_triangulation (points, cells, SubCellData());
-
-  return (0);
-}
-
-
-int plot_shapes (DoFHandler<2> &dof_handler)
-{
-  Vector<double> solution (dof_handler.n_dofs ());
-  std::set<unsigned int> face_dofs;
-  
-                                  // Create set of all DoFs, which are on the boundary.
-  DoFHandler<2>::active_face_iterator face = dof_handler.begin_active_face(),
-                                     endf = dof_handler.end_face();
-  std::vector<unsigned int> face_dof_indices (dof_handler.get_fe().dofs_per_face);
-  for (; face!=endf; ++face)
-    {
-      face->get_dof_indices (face_dof_indices);
-      for (unsigned int i=0; i<dof_handler.get_fe().dofs_per_face; ++i)
-       face_dofs.insert (face_dof_indices[i]);
-    }
-  
-                                  // Now set a 1 at the place of the different DoFs and
-                                  // output the solution.
-  std::set<unsigned int>::iterator face_dof_iter;
-  for (face_dof_iter = face_dofs.begin ();
-       face_dof_iter != face_dofs.end (); ++face_dof_iter)
-    {
-      unsigned int dof = *face_dof_iter;
-
-      sprintf (buf, "%i\n", dof);
-      deallog << buf;
-
-      solution (dof) = 1.0;
-
-                                      // Test the core functionality2
-      DataOut<2> *data_out = new DataOut<2>;
-      data_out->attach_dof_handler (dof_handler);
-      data_out->add_data_vector (solution, "solution");
-      data_out->build_patches (4, 0);
-
-      data_out->write_gnuplot (deallog.get_file_stream());
-
-      delete data_out;
-
-      solution (dof) = 0.0;
-    }
-
-  return (0);
-}
-
-
-int main (int /*argc*/, char **/*argv*/)
-{
-  logfile.precision (PRECISION);
-  logfile.setf(std::ios::fixed);  
-  deallog.attach(logfile);
-  deallog.depth_console(0);
-  deallog.threshold_double(1.e-10);
-
-  Triangulation<2> tria_test;
-  DoFHandler<2> *dof_handler;
-
   FE_RaviartThomas<2> fe (0);
-                                  //FESystem<2> fe (FE_Q<2>(3), 2);
-                                  //FE_ABF<2> fe (0);
-  deallog << "Dofs/cell " << fe.dofs_per_cell
-         << "Dofs/face " << fe.dofs_per_face << std::endl;
-
-  for (unsigned int elm = 0; elm < 3; ++elm)
-    {
-      create_tria (elm, tria_test);
-
-                                      // Create a DoFHandler
-      dof_handler = new DoFHandler<2> (tria_test);
-      dof_handler->distribute_dofs (fe);
-
-      deallog << "Dofs total " << dof_handler->n_dofs () << std::endl;
-
-                                      // Alloc some DoFs
-      Vector<double> solution;
-      solution.reinit (dof_handler->n_dofs ());
-
-                                      // Project solution onto FE field
-      ConstraintMatrix     hn_constraints;
-      hn_constraints.clear ();
-      DoFTools::make_hanging_node_constraints (*dof_handler, 
-                                              hn_constraints);
-      hn_constraints.close ();
-
-      MappingQ1<2> map_default;
-      
-      project (map_default, *dof_handler, hn_constraints,
-              QGauss6<2> (), TestMap1<2>(2),
-              solution);
-
-                                      // Test the core functionality
-      DataOut<2> *data_out = new DataOut<2>;
-      data_out->attach_dof_handler (*dof_handler);
-      data_out->add_data_vector (solution, "solution");
-      data_out->build_patches (4);
-
-      data_out->write_gnuplot (deallog.get_file_stream());
-
-                                      // Now write face DoFs ..
-      for (unsigned int i = 0; i < 4; ++i)
-       {
-         sprintf (buf, "%i %e\n", i, solution (i));
-         deallog << buf;
-       }
-      
-                                      // Clean up ...
-      delete data_out;
-      delete (dof_handler);
-      tria_test.clear ();
-    }
-
-  return (0);
+  check (fe);
 }
index a6f3e7bae912e45f0701e0bc7fa7f30f8a5d4d37..f6800890d1b4ad52e1e899558600a8106a5c3fe7 100644 (file)
@@ -2,8 +2,8 @@
 DEAL::Dofs/cell 4Dofs/face 1
 DEAL::Dofs total 4
 DEAL::MM created
-RHS created
-Starting value 0.75
+DEAL::RHS created
+DEAL:cg::Starting value 0.75
 DEAL:cg::Convergence step 4 value 0
 # This file was generated by the deal.II library.
 
@@ -43,14 +43,14 @@ DEAL:cg::Convergence step 4 value 0
 1.00 0.50 1.00 1.00 
 
 
-DEAL::0 -1.000000e-01
-1 5.000000e-01
-2 1.000000e+00
-3 4.000000e-01
-Dofs total 4
+DEAL::-0.10
+DEAL::0.50
+DEAL::1.00
+DEAL::0.40
+DEAL::Dofs total 4
 DEAL::MM created
-RHS created
-Starting value 0.45
+DEAL::RHS created
+DEAL:cg::Starting value 0.45
 DEAL:cg::Convergence step 4 value 0
 # This file was generated by the deal.II library.
 
@@ -90,14 +90,14 @@ DEAL:cg::Convergence step 4 value 0
 1.00 1.00 1.00 1.00 
 
 
-DEAL::0 6.000000e-01
-1 5.000000e-01
-2 4.000000e-01
-3 5.000000e-01
-Dofs total 4
+DEAL::0.60
+DEAL::0.50
+DEAL::0.40
+DEAL::0.50
+DEAL::Dofs total 4
 DEAL::MM created
-RHS created
-Starting value 0.75
+DEAL::RHS created
+DEAL:cg::Starting value 0.75
 DEAL:cg::Convergence step 4 value 0
 # This file was generated by the deal.II library.
 
@@ -137,7 +137,7 @@ DEAL:cg::Convergence step 4 value 0
 0.50 1.00 1.00 1.00 
 
 
-DEAL::0 1.000000e+00
-1 6.000000e-01
-2 1.000000e-01
-3 5.000000e-01
+DEAL::1.00
+DEAL::0.60
+DEAL::0.10
+DEAL::0.50

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.