]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add test case for static condensation with FE_TraceQ. 585/head
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 23 Feb 2015 12:04:19 +0000 (13:04 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 23 Feb 2015 12:04:19 +0000 (13:04 +0100)
tests/bits/static_condensation.cc [new file with mode: 0644]
tests/bits/static_condensation.output [new file with mode: 0644]

diff --git a/tests/bits/static_condensation.cc b/tests/bits/static_condensation.cc
new file mode 100644 (file)
index 0000000..14ffae2
--- /dev/null
@@ -0,0 +1,732 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2005 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// a un-hp-ified version of hp/step-7
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <fstream>
+std::ofstream logfile("output");
+
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_trace.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/base/smartpointer.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/base/convergence_table.h>
+#include <deal.II/hp/fe_values.h>
+
+#include <typeinfo>
+#include <fstream>
+#include <iomanip>
+
+
+
+template <int dim>
+class SolutionBase
+{
+protected:
+  static const unsigned int n_source_centers = 3;
+  static const Point<dim>   source_centers[n_source_centers];
+  static const double       width;
+};
+
+
+template <>
+const Point<1>
+SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers]
+  = { Point<1>(-1.0 / 3.0),
+      Point<1>(0.0),
+      Point<1>(+1.0 / 3.0)
+    };
+
+template <>
+const Point<2>
+SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers]
+  = { Point<2>(-0.5, +0.5),
+      Point<2>(-0.5, -0.5),
+      Point<2>(+0.5, -0.5)
+    };
+
+template <>
+const Point<3>
+SolutionBase<3>::source_centers[SolutionBase<3>::n_source_centers]
+  = { Point<3>(-0.5, +0.5, 0.5),
+      Point<3>(-0.5, -0.5, -0.5),
+      Point<3>(+0.5, -0.5, 0)
+    };
+
+template <int dim>
+const double SolutionBase<dim>::width = 1./3.;
+
+
+
+template <int dim>
+class Solution : public Function<dim>,
+  protected SolutionBase<dim>
+{
+public:
+  Solution () : Function<dim>() {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component = 0) const;
+
+  virtual Tensor<1,dim> gradient (const Point<dim>   &p,
+                                  const unsigned int  component = 0) const;
+};
+
+
+template <int dim>
+double Solution<dim>::value (const Point<dim>   &p,
+                             const unsigned int) const
+{
+  double return_value = 0;
+  for (unsigned int i=0; i<this->n_source_centers; ++i)
+    {
+      const Tensor<1,dim> x_minus_xi = p - this->source_centers[i];
+      return_value += std::exp(-x_minus_xi.norm_square() /
+                               (this->width * this->width));
+    }
+
+  return return_value;
+}
+
+
+template <int dim>
+Tensor<1,dim> Solution<dim>::gradient (const Point<dim>   &p,
+                                       const unsigned int) const
+{
+  Tensor<1,dim> return_value;
+
+  for (unsigned int i=0; i<this->n_source_centers; ++i)
+    {
+      const Tensor<1,dim> x_minus_xi = p - this->source_centers[i];
+
+      return_value += (-2 / (this->width * this->width) *
+                       std::exp(-x_minus_xi.norm_square() /
+                                (this->width * this->width)) *
+                       x_minus_xi);
+    }
+
+  return return_value;
+}
+
+
+
+template <int dim>
+class RightHandSide : public Function<dim>,
+  protected SolutionBase<dim>
+{
+public:
+  RightHandSide () : Function<dim>() {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component = 0) const;
+};
+
+
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim>   &p,
+                                  const unsigned int) const
+{
+  double return_value = 0;
+  for (unsigned int i=0; i<this->n_source_centers; ++i)
+    {
+      const Tensor<1,dim> x_minus_xi = p - this->source_centers[i];
+
+      return_value += ((2*dim - 4*x_minus_xi.norm_square()/
+                        (this->width * this->width)) /
+                       (this->width * this->width) *
+                       std::exp(-x_minus_xi.norm_square() /
+                                (this->width * this->width)));
+      return_value += std::exp(-x_minus_xi.norm_square() /
+                               (this->width * this->width));
+    }
+
+  return return_value;
+}
+
+
+
+template <int dim>
+class HelmholtzProblem
+{
+public:
+  enum RefinementMode
+  {
+    global_refinement, adaptive_refinement
+  };
+
+  HelmholtzProblem (const unsigned int   fe_degree,
+                    const RefinementMode refinement_mode);
+
+  void run ();
+
+private:
+  void setup_system ();
+  void assemble_system (const bool do_reconstruct);
+  void solve ();
+  void refine_grid ();
+  void process_solution (const unsigned int cycle);
+
+  Triangulation<dim>                      triangulation;
+  FE_Q<dim>                               fe;
+  DoFHandler<dim>                         dof_handler;
+  ConstraintMatrix                        constraints;
+  SparsityPattern                         sparsity_pattern;
+  SparseMatrix<double>                    system_matrix;
+  Vector<double>                          solution;
+  Vector<double>                          system_rhs;
+
+  FE_TraceQ<dim>                          fe_trace;
+  DoFHandler<dim>                         dof_handler_trace;
+  ConstraintMatrix                        constraints_trace;
+  SparsityPattern                         sparsity_pattern_trace;
+  SparseMatrix<double>                    system_matrix_trace;
+  Vector<double>                          solution_trace_full;
+  Vector<double>                          solution_trace;
+  Vector<double>                          system_rhs_trace;
+
+  const RefinementMode                    refinement_mode;
+
+  ConvergenceTable                        convergence_table;
+};
+
+
+
+
+template <int dim>
+HelmholtzProblem<dim>::HelmholtzProblem (const unsigned int   fe_degree,
+                                         const RefinementMode refinement_mode) :
+  fe (fe_degree),
+  dof_handler (triangulation),
+  fe_trace (fe_degree),
+  dof_handler_trace (triangulation),
+  refinement_mode (refinement_mode)
+{
+  deallog << "Solving with Q" << fe_degree << " elements, "
+          << (refinement_mode == global_refinement ? "global" : "adaptive")
+          << " refinement" << std::endl
+          << "==========================================="
+          << (refinement_mode == global_refinement ? "" : "==") << std::endl
+          << std::endl;
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (fe);
+
+  constraints.clear ();
+  DoFTools::make_hanging_node_constraints (dof_handler,
+                                           constraints);
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                            0,
+                                            Solution<dim>(),
+                                            constraints);
+  constraints.close ();
+
+  {
+    CompressedSimpleSparsityPattern csp (dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false);
+    sparsity_pattern.copy_from(csp);
+  }
+
+  system_matrix.reinit (sparsity_pattern);
+
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+
+
+  dof_handler_trace.distribute_dofs (fe_trace);
+
+  constraints_trace.clear ();
+  DoFTools::make_hanging_node_constraints (dof_handler_trace,
+                                           constraints_trace);
+  VectorTools::interpolate_boundary_values (dof_handler_trace,
+                                            0,
+                                            Solution<dim>(),
+                                            constraints_trace);
+  constraints_trace.close ();
+
+  {
+    CompressedSimpleSparsityPattern csp (dof_handler_trace.n_dofs());
+    DoFTools::make_sparsity_pattern (dof_handler_trace, csp, constraints_trace, false);
+    sparsity_pattern_trace.copy_from(csp);
+  }
+
+  system_matrix_trace.reinit (sparsity_pattern_trace);
+
+  solution_trace.reinit (dof_handler_trace.n_dofs());
+  system_rhs_trace.reinit (dof_handler_trace.n_dofs());
+  solution_trace_full.reinit (dof_handler.n_dofs());
+
+  deallog << "Number of DoFs:       " << dof_handler.n_dofs() << " / "
+          << dof_handler_trace.n_dofs() << std::endl;
+  deallog << "Number of matrix nnz: " << sparsity_pattern.n_nonzero_elements() << " / "
+          << sparsity_pattern_trace.n_nonzero_elements() << std::endl;
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::assemble_system (const bool do_reconstruct)
+{
+  QGauss<dim>   quadrature_formula (fe.degree+1);
+  QGauss<dim-1> face_quadrature_formula (fe.degree+1);
+
+  const unsigned int n_q_points    = quadrature_formula.size();
+  const unsigned int n_face_q_points = face_quadrature_formula.size();
+
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       cell_rhs (dofs_per_cell);
+
+  FullMatrix<double>   trace_matrix (fe_trace.dofs_per_cell, fe_trace.dofs_per_cell);
+  Vector<double>       trace_rhs (fe_trace.dofs_per_cell);
+
+  FullMatrix<double>   eliminate_matrix(cell_matrix.m()-trace_matrix.m(),
+                                        cell_matrix.m()-trace_matrix.m());
+  FullMatrix<double>   temp_matrix(trace_matrix.m(), eliminate_matrix.m());
+
+  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+  std::vector<types::global_dof_index> dof_indices_trace (fe_trace.dofs_per_cell);
+  Vector<double> local_trace(fe_trace.dofs_per_cell);
+  Vector<double> local_interior(dofs_per_cell-fe_trace.dofs_per_cell);
+
+  FEValues<dim>  x_fe_values (fe, quadrature_formula,
+                              update_values   | update_gradients |
+                              update_q_points | update_JxW_values);
+
+  FEFaceValues<dim> x_fe_face_values (fe, face_quadrature_formula,
+                                      update_values         | update_q_points  |
+                                      update_normal_vectors | update_JxW_values);
+
+  const RightHandSide<dim> right_hand_side;
+  std::vector<double>  rhs_values (n_q_points);
+
+  const Solution<dim> exact_solution;
+
+  typename DoFHandler<dim>::active_cell_iterator
+  cell = dof_handler.begin_active(),
+  endc = dof_handler.end(),
+  tracec = dof_handler_trace.begin_active();
+  for (; cell!=endc; ++cell, ++tracec)
+    {
+      cell_rhs = 0;
+
+      x_fe_values.reinit (cell);
+      const FEValues<dim> &fe_values = x_fe_values.get_present_fe_values();
+
+      right_hand_side.value_list (fe_values.get_quadrature_points(),
+                                  rhs_values);
+
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        {
+          for (unsigned int j=0; j<dofs_per_cell; ++j)
+            {
+              double sum = 0;
+              const Tensor<1,dim> *shape_grad_i = &fe_values.shape_grad(i,0);
+              const Tensor<1,dim> *shape_grad_j = &fe_values.shape_grad(j,0);
+              const double *shape_value_i = &fe_values.shape_value(i,0);
+              const double *shape_value_j = &fe_values.shape_value(j,0);
+              for (unsigned int q_index=0; q_index<n_q_points; ++q_index)
+                sum += (shape_grad_i[q_index] * shape_grad_j[q_index] +
+                        shape_value_i[q_index] * shape_value_j[q_index]
+                       ) * fe_values.JxW(q_index);
+              cell_matrix(i,j) = sum;
+            }
+
+          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+            cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+                            rhs_values [q_point] *
+                            fe_values.JxW(q_point));
+        }
+
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+        if (cell->face(face)->at_boundary()
+            &&
+            (cell->face(face)->boundary_indicator() == 1))
+          {
+            x_fe_face_values.reinit (cell, face);
+            const FEFaceValues<dim> &fe_face_values
+              = x_fe_face_values.get_present_fe_values ();
+
+            for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
+              {
+                const double neumann_value
+                  = (exact_solution.gradient (fe_face_values.quadrature_point(q_point)) *
+                     fe_face_values.normal_vector(q_point));
+
+                for (unsigned int i=0; i<dofs_per_cell; ++i)
+                  cell_rhs(i) += (neumann_value *
+                                  fe_face_values.shape_value(i,q_point) *
+                                  fe_face_values.JxW(q_point));
+              }
+          }
+
+      const unsigned int shift = fe_trace.dofs_per_cell;
+      const unsigned int sizes = fe.dofs_per_cell - fe_trace.dofs_per_cell;
+      for (unsigned int i=0; i<sizes; ++i)
+        for (unsigned int j=0; j<sizes; ++j)
+          eliminate_matrix(i,j) = cell_matrix(i+shift, j+shift);
+      if (sizes > 0)
+        eliminate_matrix.gauss_jordan();
+
+      cell->get_dof_indices (local_dof_indices);
+      if (do_reconstruct == false)
+        {
+          cell->get_dof_indices (local_dof_indices);
+          constraints.distribute_local_to_global(cell_matrix, cell_rhs,
+                                                 local_dof_indices,
+                                                 system_matrix, system_rhs);
+
+          for (unsigned int i=0; i<shift; ++i)
+            for (unsigned int j=0; j<sizes; ++j)
+              {
+                double sum = 0;
+                for (unsigned int k=0; k<sizes; ++k)
+                  sum += cell_matrix(i,k+shift) * eliminate_matrix(k,j);
+                temp_matrix(i,j) = sum;
+              }
+          for (unsigned int i=0; i<shift; ++i)
+            for (unsigned int j=0; j<shift; ++j)
+              {
+                double sum = 0;
+                for (unsigned int k=0; k<sizes; ++k)
+                  sum += temp_matrix(i,k) * cell_matrix(shift+k,j);
+                trace_matrix(i,j) = cell_matrix(i,j) - sum;
+              }
+          for (unsigned int i=0; i<shift; ++i)
+            {
+              double sum = 0;
+              for (unsigned int k=0; k<sizes; ++k)
+                sum += temp_matrix(i,k) * cell_rhs(shift+k);
+              trace_rhs(i) = cell_rhs(i) - sum;
+            }
+
+          tracec->get_dof_indices (dof_indices_trace);
+          constraints_trace.distribute_local_to_global(trace_matrix, trace_rhs,
+                                                       dof_indices_trace,
+                                                       system_matrix_trace,
+                                                       system_rhs_trace);
+        }
+      else
+        {
+          tracec->get_interpolated_dof_values(solution_trace, local_trace);
+          for (unsigned int i=0; i<shift; ++i)
+            solution_trace_full(local_dof_indices[i]) = local_trace(i);
+          for (unsigned int i=0; i<sizes; ++i)
+            {
+              double sum = 0;
+              for (unsigned int k=0; k<shift; ++k)
+                sum += cell_matrix(shift+i,k) * local_trace(k);
+              local_interior(i) = cell_rhs(shift+i) - sum;
+            }
+          for (unsigned int i=0; i<sizes; ++i)
+            {
+              double sum = 0;
+              for (unsigned int j=0; j<sizes; ++j)
+                sum += eliminate_matrix(i,j) * local_interior(j);
+              solution_trace_full(local_dof_indices[shift+i]) = sum;
+            }
+        }
+    }
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::solve ()
+{
+  {
+    SolverControl           solver_control (1000, 1e-12);
+    SolverCG<>              cg (solver_control);
+
+    PreconditionSSOR<> preconditioner;
+    preconditioner.initialize(system_matrix, 1.2);
+
+    cg.solve (system_matrix, solution, system_rhs,
+              preconditioner);
+
+    constraints.distribute (solution);
+  }
+  {
+    SolverControl           solver_control (1000, 1e-12);
+    SolverCG<>              cg (solver_control);
+
+    PreconditionSSOR<> preconditioner;
+    preconditioner.initialize(system_matrix_trace, 1.2);
+
+    cg.solve (system_matrix_trace, solution_trace, system_rhs_trace,
+              preconditioner);
+
+    constraints_trace.distribute (solution_trace);
+  }
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::refine_grid ()
+{
+  switch (refinement_mode)
+    {
+    case global_refinement:
+    {
+      triangulation.refine_global (1);
+      break;
+    }
+
+    case adaptive_refinement:
+    {
+      Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+      typename FunctionMap<dim>::type neumann_boundary;
+      KellyErrorEstimator<dim>::estimate (dof_handler,
+                                          QGauss<dim-1>(3),
+                                          neumann_boundary,
+                                          solution,
+                                          estimated_error_per_cell);
+
+      GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                       estimated_error_per_cell,
+                                                       0.3, 0.03);
+
+      triangulation.execute_coarsening_and_refinement ();
+
+      break;
+    }
+
+    default:
+    {
+      Assert (false, ExcNotImplemented());
+    }
+    }
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::process_solution (const unsigned int cycle)
+{
+  {
+    Vector<float> difference_per_cell (triangulation.n_active_cells());
+    VectorTools::integrate_difference (dof_handler,
+                                       solution,
+                                       Solution<dim>(),
+                                       difference_per_cell,
+                                       QGauss<dim>(fe.degree+2),
+                                       VectorTools::L2_norm);
+    const double L2_error = difference_per_cell.l2_norm();
+
+    VectorTools::integrate_difference (dof_handler,
+                                       solution,
+                                       Solution<dim>(),
+                                       difference_per_cell,
+                                       QGauss<dim>(fe.degree+2),
+                                       VectorTools::H1_seminorm);
+    const double H1_error = difference_per_cell.l2_norm();
+
+    const unsigned int n_active_cells=triangulation.n_active_cells();
+    const unsigned int n_dofs=dof_handler.n_dofs();
+
+    convergence_table.add_value("cycle", cycle);
+    convergence_table.add_value("cells", n_active_cells);
+    convergence_table.add_value("dofs", n_dofs);
+    convergence_table.add_value("L2", L2_error);
+    convergence_table.add_value("H1", H1_error);
+  }
+  {
+    Vector<float> difference_per_cell (triangulation.n_active_cells());
+    VectorTools::integrate_difference (dof_handler,
+                                       solution_trace_full,
+                                       Solution<dim>(),
+                                       difference_per_cell,
+                                       QGauss<dim>(fe.degree+2),
+                                       VectorTools::L2_norm);
+    const double L2_error = difference_per_cell.l2_norm();
+
+    VectorTools::integrate_difference (dof_handler,
+                                       solution_trace_full,
+                                       Solution<dim>(),
+                                       difference_per_cell,
+                                       QGauss<dim>(fe.degree+2),
+                                       VectorTools::H1_seminorm);
+    const double H1_error = difference_per_cell.l2_norm();
+
+    const unsigned int n_active_cells=triangulation.n_active_cells();
+    const unsigned int n_dofs=dof_handler.n_dofs();
+
+    convergence_table.add_value("L2 sc", L2_error);
+    convergence_table.add_value("H1 sc", H1_error);
+  }
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::run ()
+{
+  for (unsigned int cycle=0; cycle<6-dim; ++cycle)
+    {
+      if (cycle == 0)
+        {
+          GridGenerator::hyper_cube (triangulation, -1, 1);
+          triangulation.refine_global (1);
+
+          typename Triangulation<dim>::cell_iterator
+          cell = triangulation.begin (),
+          endc = triangulation.end();
+          for (; cell!=endc; ++cell)
+            for (unsigned int face=0;
+                 face<GeometryInfo<dim>::faces_per_cell;
+                 ++face)
+              if ((cell->face(face)->center()(0) == -1)
+                  ||
+                  (cell->face(face)->center()(1) == -1))
+                cell->face(face)->set_boundary_indicator (1);
+        }
+      else
+        refine_grid ();
+
+
+      setup_system ();
+
+      assemble_system (false);
+      solve ();
+      assemble_system (true);
+
+      process_solution (cycle);
+    }
+
+  convergence_table.set_precision("L2", 3);
+  convergence_table.set_precision("H1", 3);
+
+  convergence_table.set_scientific("L2", true);
+  convergence_table.set_scientific("H1", true);
+
+  convergence_table.set_precision("L2 sc", 3);
+  convergence_table.set_precision("H1 sc", 3);
+
+  convergence_table.set_scientific("L2 sc", true);
+  convergence_table.set_scientific("H1 sc", true);
+
+  deallog << std::endl;
+  convergence_table.write_text(deallog.get_file_stream());
+}
+
+
+int main ()
+{
+  deallog << std::setprecision(2);
+  logfile << std::setprecision(2);
+
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      deallog.depth_console (0);
+
+
+      for (unsigned int deg = 1; deg < 5; ++deg)
+        {
+          HelmholtzProblem<2>
+          helmholtz_problem_2d (deg, HelmholtzProblem<2>::adaptive_refinement);
+
+          helmholtz_problem_2d.run ();
+
+          deallog << std::endl;
+        }
+      for (unsigned int deg = 1; deg < 4; ++deg)
+        {
+          HelmholtzProblem<2>
+          helmholtz_problem_2d (deg, HelmholtzProblem<2>::global_refinement);
+
+          helmholtz_problem_2d.run ();
+
+          deallog << std::endl;
+        }
+
+      {
+        HelmholtzProblem<3>
+        helmholtz_problem_3d (2, HelmholtzProblem<3>::adaptive_refinement);
+
+        helmholtz_problem_3d.run ();
+
+        deallog << std::endl;
+      }
+    }
+  catch (std::exception &exc)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Exception on processing: " << std::endl
+              << exc.what() << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      return 1;
+    }
+  catch (...)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Unknown exception!" << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      return 1;
+    }
+
+  return 0;
+}
+
+
+template const double SolutionBase<2>::width;
diff --git a/tests/bits/static_condensation.output b/tests/bits/static_condensation.output
new file mode 100644 (file)
index 0000000..1c8d2ce
--- /dev/null
@@ -0,0 +1,266 @@
+
+DEAL::Solving with Q1 elements, adaptive refinement
+DEAL::=============================================
+DEAL::
+DEAL::Number of DoFs:       9 / 9
+DEAL::Number of matrix nnz: 21 / 21
+DEAL:cg::Starting value 5.5
+DEAL:cg::Convergence step 4 value 0
+DEAL:cg::Starting value 5.5
+DEAL:cg::Convergence step 4 value 0
+DEAL::Number of DoFs:       22 / 22
+DEAL::Number of matrix nnz: 86 / 86
+DEAL:cg::Starting value 6.0
+DEAL:cg::Convergence step 11 value 0
+DEAL:cg::Starting value 6.0
+DEAL:cg::Convergence step 11 value 0
+DEAL::Number of DoFs:       46 / 46
+DEAL::Number of matrix nnz: 222 / 222
+DEAL:cg::Starting value 5.3
+DEAL:cg::Convergence step 14 value 0
+DEAL:cg::Starting value 5.3
+DEAL:cg::Convergence step 14 value 0
+DEAL::Number of DoFs:       84 / 84
+DEAL::Number of matrix nnz: 444 / 444
+DEAL:cg::Starting value 3.7
+DEAL:cg::Convergence step 17 value 0
+DEAL:cg::Starting value 3.7
+DEAL:cg::Convergence step 17 value 0
+DEAL::
+cycle cells dofs    L2        H1       L2 sc     H1 sc   
+0     4     9    4.037e+00 4.590e+00 4.037e+00 4.590e+00 
+1     13    22   5.358e-01 1.430e+00 5.358e-01 1.430e+00 
+2     31    46   3.087e-01 1.258e+00 3.087e-01 1.258e+00 
+3     61    84   7.744e-02 1.096e+00 7.744e-02 1.096e+00 
+DEAL::
+DEAL::Solving with Q2 elements, adaptive refinement
+DEAL::=============================================
+DEAL::
+DEAL::Number of DoFs:       25 / 21
+DEAL::Number of matrix nnz: 153 / 107
+DEAL:cg::Starting value 11.
+DEAL:cg::Convergence step 11 value 0
+DEAL:cg::Starting value 1.5
+DEAL:cg::Convergence step 11 value 0
+DEAL::Number of DoFs:       71 / 58
+DEAL::Number of matrix nnz: 643 / 460
+DEAL:cg::Starting value 3.9
+DEAL:cg::Convergence step 19 value 0
+DEAL:cg::Starting value 3.9
+DEAL:cg::Convergence step 16 value 0
+DEAL::Number of DoFs:       148 / 120
+DEAL::Number of matrix nnz: 1492 / 1062
+DEAL:cg::Starting value 3.4
+DEAL:cg::Convergence step 27 value 0
+DEAL:cg::Starting value 3.6
+DEAL:cg::Convergence step 21 value 0
+DEAL::Number of DoFs:       265 / 210
+DEAL::Number of matrix nnz: 3101 / 2240
+DEAL:cg::Starting value 2.6
+DEAL:cg::Convergence step 39 value 0
+DEAL:cg::Starting value 2.6
+DEAL:cg::Convergence step 28 value 0
+DEAL::
+cycle cells dofs    L2        H1       L2 sc     H1 sc   
+0     4     25   1.471e+00 2.988e+00 1.471e+00 2.988e+00 
+1     13    71   9.942e-02 1.127e+00 9.942e-02 1.127e+00 
+2     28    148  7.311e-02 9.526e-01 7.311e-02 9.526e-01 
+3     55    265  1.091e-02 2.577e-01 1.091e-02 2.577e-01 
+DEAL::
+DEAL::Solving with Q3 elements, adaptive refinement
+DEAL::=============================================
+DEAL::
+DEAL::Number of DoFs:       49 / 33
+DEAL::Number of matrix nnz: 589 / 261
+DEAL:cg::Starting value 3.9
+DEAL:cg::Convergence step 20 value 0
+DEAL:cg::Starting value 0.30
+DEAL:cg::Convergence step 14 value 0
+DEAL::Number of DoFs:       146 / 94
+DEAL::Number of matrix nnz: 2378 / 1122
+DEAL:cg::Starting value 3.7
+DEAL:cg::Convergence step 32 value 0
+DEAL:cg::Starting value 3.3
+DEAL:cg::Convergence step 22 value 0
+DEAL::Number of DoFs:       342 / 218
+DEAL::Number of matrix nnz: 6126 / 2902
+DEAL:cg::Starting value 3.0
+DEAL:cg::Convergence step 47 value 0
+DEAL:cg::Starting value 3.0
+DEAL:cg::Convergence step 30 value 0
+DEAL::Number of DoFs:       630 / 386
+DEAL::Number of matrix nnz: 12314 / 5938
+DEAL:cg::Starting value 1.8
+DEAL:cg::Convergence step 65 value 0
+DEAL:cg::Starting value 2.0
+DEAL:cg::Convergence step 39 value 0
+DEAL::
+cycle cells dofs    L2        H1       L2 sc     H1 sc   
+0     4     49   3.566e-01 1.675e+00 3.566e-01 1.675e+00 
+1     13    146  1.336e-02 2.093e-01 1.336e-02 2.093e-01 
+2     31    342  7.824e-03 1.526e-01 7.824e-03 1.526e-01 
+3     61    630  1.181e-03 3.811e-02 1.181e-03 3.811e-02 
+DEAL::
+DEAL::Solving with Q4 elements, adaptive refinement
+DEAL::=============================================
+DEAL::
+DEAL::Number of DoFs:       81 / 45
+DEAL::Number of matrix nnz: 1617 / 483
+DEAL:cg::Starting value 5.9
+DEAL:cg::Convergence step 30 value 0
+DEAL:cg::Starting value 0.75
+DEAL:cg::Convergence step 17 value 0
+DEAL::Number of DoFs:       247 / 130
+DEAL::Number of matrix nnz: 6311 / 2072
+DEAL:cg::Starting value 2.9
+DEAL:cg::Convergence step 49 value 0
+DEAL:cg::Starting value 2.8
+DEAL:cg::Convergence step 26 value 0
+DEAL::Number of DoFs:       583 / 304
+DEAL::Number of matrix nnz: 16071 / 5334
+DEAL:cg::Starting value 2.5
+DEAL:cg::Convergence step 73 value 0
+DEAL:cg::Starting value 2.5
+DEAL:cg::Convergence step 35 value 0
+DEAL::Number of DoFs:       1113 / 564
+DEAL::Number of matrix nnz: 32325 / 10914
+DEAL:cg::Starting value 1.5
+DEAL:cg::Convergence step 98 value 0
+DEAL:cg::Starting value 1.6
+DEAL:cg::Convergence step 46 value 0
+DEAL::
+cycle cells dofs    L2        H1       L2 sc     H1 sc   
+0     4     81   4.573e-02 5.591e-01 4.573e-02 5.591e-01 
+1     13    247  3.694e-03 8.184e-02 3.694e-03 8.184e-02 
+2     31    583  2.238e-03 5.906e-02 2.238e-03 5.906e-02 
+3     61    1113 1.549e-04 6.438e-03 1.549e-04 6.438e-03 
+DEAL::
+DEAL::Solving with Q1 elements, global refinement
+DEAL::===========================================
+DEAL::
+DEAL::Number of DoFs:       9 / 9
+DEAL::Number of matrix nnz: 21 / 21
+DEAL:cg::Starting value 5.5
+DEAL:cg::Convergence step 4 value 0
+DEAL:cg::Starting value 5.5
+DEAL:cg::Convergence step 4 value 0
+DEAL::Number of DoFs:       25 / 25
+DEAL::Number of matrix nnz: 109 / 109
+DEAL:cg::Starting value 6.0
+DEAL:cg::Convergence step 11 value 0
+DEAL:cg::Starting value 6.0
+DEAL:cg::Convergence step 11 value 0
+DEAL::Number of DoFs:       81 / 81
+DEAL::Number of matrix nnz: 501 / 501
+DEAL:cg::Starting value 3.8
+DEAL:cg::Convergence step 16 value 0
+DEAL:cg::Starting value 3.8
+DEAL:cg::Convergence step 16 value 0
+DEAL::Number of DoFs:       289 / 289
+DEAL::Number of matrix nnz: 2149 / 2149
+DEAL:cg::Starting value 2.3
+DEAL:cg::Convergence step 26 value 0
+DEAL:cg::Starting value 2.3
+DEAL:cg::Convergence step 26 value 0
+DEAL::
+cycle cells dofs    L2        H1       L2 sc     H1 sc   
+0     4     9    4.037e+00 4.590e+00 4.037e+00 4.590e+00 
+1     16    25   4.926e-01 1.409e+00 4.926e-01 1.409e+00 
+2     64    81   7.984e-02 1.130e+00 7.984e-02 1.130e+00 
+3     256   289  2.101e-02 5.828e-01 2.101e-02 5.828e-01 
+DEAL::
+DEAL::Solving with Q2 elements, global refinement
+DEAL::===========================================
+DEAL::
+DEAL::Number of DoFs:       25 / 21
+DEAL::Number of matrix nnz: 153 / 107
+DEAL:cg::Starting value 11.
+DEAL:cg::Convergence step 11 value 0
+DEAL:cg::Starting value 1.5
+DEAL:cg::Convergence step 11 value 0
+DEAL::Number of DoFs:       81 / 65
+DEAL::Number of matrix nnz: 801 / 575
+DEAL:cg::Starting value 3.8
+DEAL:cg::Convergence step 20 value 0
+DEAL:cg::Starting value 3.8
+DEAL:cg::Convergence step 17 value 0
+DEAL::Number of DoFs:       289 / 225
+DEAL::Number of matrix nnz: 3633 / 2639
+DEAL:cg::Starting value 2.6
+DEAL:cg::Convergence step 40 value 0
+DEAL:cg::Starting value 2.6
+DEAL:cg::Convergence step 29 value 0
+DEAL::Number of DoFs:       1089 / 833
+DEAL::Number of matrix nnz: 15441 / 11279
+DEAL:cg::Starting value 1.4
+DEAL:cg::Convergence step 79 value 0
+DEAL:cg::Starting value 1.5
+DEAL:cg::Convergence step 53 value 0
+DEAL::
+cycle cells dofs    L2        H1       L2 sc     H1 sc   
+0     4     25   1.471e+00 2.988e+00 1.471e+00 2.988e+00 
+1     16    81   9.523e-02 1.122e+00 9.523e-02 1.122e+00 
+2     64    289  9.409e-03 2.484e-01 9.409e-03 2.484e-01 
+3     256   1089 1.194e-03 6.228e-02 1.194e-03 6.228e-02 
+DEAL::
+DEAL::Solving with Q3 elements, global refinement
+DEAL::===========================================
+DEAL::
+DEAL::Number of DoFs:       49 / 33
+DEAL::Number of matrix nnz: 589 / 261
+DEAL:cg::Starting value 3.9
+DEAL:cg::Convergence step 20 value 0
+DEAL:cg::Starting value 0.30
+DEAL:cg::Convergence step 14 value 0
+DEAL::Number of DoFs:       169 / 105
+DEAL::Number of matrix nnz: 2941 / 1397
+DEAL:cg::Starting value 3.7
+DEAL:cg::Convergence step 34 value 0
+DEAL:cg::Starting value 3.3
+DEAL:cg::Convergence step 23 value 0
+DEAL::Number of DoFs:       625 / 369
+DEAL::Number of matrix nnz: 13045 / 6381
+DEAL:cg::Starting value 1.8
+DEAL:cg::Convergence step 63 value 0
+DEAL:cg::Starting value 2.0
+DEAL:cg::Convergence step 38 value 0
+DEAL::Number of DoFs:       2401 / 1377
+DEAL::Number of matrix nnz: 54853 / 27197
+DEAL:cg::Starting value 1.1
+DEAL:cg::Convergence step 120 value 0
+DEAL:cg::Starting value 1.2
+DEAL:cg::Convergence step 69 value 0
+DEAL::
+cycle cells dofs    L2        H1       L2 sc     H1 sc   
+0     4     49   3.566e-01 1.675e+00 3.566e-01 1.675e+00 
+1     16    169  9.770e-03 1.813e-01 9.770e-03 1.813e-01 
+2     64    625  9.013e-04 3.472e-02 9.013e-04 3.472e-02 
+3     256   2401 6.648e-05 5.054e-03 6.648e-05 5.054e-03 
+DEAL::
+DEAL::Solving with Q2 elements, adaptive refinement
+DEAL::=============================================
+DEAL::
+DEAL::Number of DoFs:       125 / 117
+DEAL::Number of matrix nnz: 1085 / 893
+DEAL:cg::Starting value 7.8
+DEAL:cg::Convergence step 11 value 0
+DEAL:cg::Starting value 3.6
+DEAL:cg::Convergence step 10 value 0
+DEAL::Number of DoFs:       407 / 378
+DEAL::Number of matrix nnz: 6443 / 5454
+DEAL:cg::Starting value 5.2
+DEAL:cg::Convergence step 18 value 0
+DEAL:cg::Starting value 1.5
+DEAL:cg::Convergence step 15 value 0
+DEAL::Number of DoFs:       1198 / 1099
+DEAL::Number of matrix nnz: 30428 / 26217
+DEAL:cg::Starting value 1.5
+DEAL:cg::Convergence step 26 value 0
+DEAL:cg::Starting value 1.5
+DEAL:cg::Convergence step 21 value 0
+DEAL::
+cycle cells dofs    L2        H1       L2 sc     H1 sc   
+0     8     125  5.433e-01 2.365e+00 5.433e-01 2.365e+00 
+1     29    407  2.488e-01 1.565e+00 2.488e-01 1.565e+00 
+2     99    1198 5.615e-02 7.982e-01 5.615e-02 7.982e-01 
+DEAL::

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.