]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
add
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 27 Jan 2008 04:39:35 +0000 (04:39 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 27 Jan 2008 04:39:35 +0000 (04:39 +0000)
git-svn-id: https://svn.dealii.org/trunk@15691 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-31/Makefile [new file with mode: 0644]
deal.II/examples/step-31/doc/intro.dox [new file with mode: 0644]
deal.II/examples/step-31/doc/results.dox [new file with mode: 0644]
deal.II/examples/step-31/step-31.cc [new file with mode: 0644]

diff --git a/deal.II/examples/step-31/Makefile b/deal.II/examples/step-31/Makefile
new file mode 100644 (file)
index 0000000..26266ad
--- /dev/null
@@ -0,0 +1,156 @@
+# $Id: Makefile,v 1.4 2006/02/10 17:53:05 wolf Exp $
+
+
+# For the small projects Makefile, you basically need to fill in only
+# four fields.
+#
+# The first is the name of the application. It is assumed that the
+# application name is the same as the base file name of the single C++
+# file from which the application is generated.
+target = $(basename $(shell echo step-*.cc))
+
+# The second field determines whether you want to run your program in
+# debug or optimized mode. The latter is significantly faster, but no
+# run-time checking of parameters and internal states is performed, so
+# you should set this value to `on' while you develop your program,
+# and to `off' when running production computations.
+debug-mode = off
+
+
+# As third field, we need to give the path to the top-level deal.II
+# directory. You need to adjust this to your needs. Since this path is
+# probably the most often needed one in the Makefile internals, it is
+# designated by a single-character variable, since that can be
+# reference using $D only, i.e. without the parentheses that are
+# required for most other parameters, as e.g. in $(target).
+D = ../../
+
+
+# The last field specifies the names of data and other files that
+# shall be deleted when calling `make clean'. Object and backup files,
+# executables and the like are removed anyway. Here, we give a list of
+# files in the various output formats that deal.II supports.
+clean-up-files = *gmv *gnuplot *gpl *eps *pov
+
+
+
+
+#
+#
+# Usually, you will not need to change something beyond this point.
+#
+#
+# The next statement tell the `make' program where to find the
+# deal.II top level directory and to include the file with the global
+# settings
+include $D/common/Make.global_options
+
+
+# Since the whole project consists of only one file, we need not
+# consider difficult dependencies. We only have to declare the
+# libraries which we want to link to the object file, and there need
+# to be two sets of libraries: one for the debug mode version of the
+# application and one for the optimized mode. Here we have selected
+# the versions for 2d. Note that the order in which the libraries are
+# given here is important and that your applications won't link
+# properly if they are given in another order.
+#
+# You may need to augment the lists of libraries when compiling your
+# program for other dimensions, or when using third party libraries
+libs.g   = $(lib-deal2-2d.g) \
+          $(lib-deal2-3d.g) \
+          $(lib-lac.g)      \
+           $(lib-base.g)
+libs.o   = $(lib-deal2-2d.o) \
+          $(lib-deal2-3d.o) \
+          $(lib-lac.o)      \
+           $(lib-base.o)
+
+
+# We now use the variable defined above which switch between debug and
+# optimized mode to select the set of libraries to link with. Included
+# in the list of libraries is the name of the object file which we
+# will produce from the single C++ file. Note that by default we use
+# the extension .g.o for object files compiled in debug mode and .o for
+# object files in optimized mode (or whatever the local default on your
+# system is instead of .o).
+ifeq ($(debug-mode),on)
+  libraries = $(target).g.$(OBJEXT) $(libs.g)
+else
+  libraries = $(target).$(OBJEXT) $(libs.o)
+endif
+
+
+# Now comes the first production rule: how to link the single object
+# file produced from the single C++ file into the executable. Since
+# this is the first rule in the Makefile, it is the one `make' selects
+# if you call it without arguments.
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS)
+
+
+# To make running the application somewhat independent of the actual
+# program name, we usually declare a rule `run' which simply runs the
+# program. You can then run it by typing `make run'. This is also
+# useful if you want to call the executable with arguments which do
+# not change frequently. You may then want to add them to the
+# following rule:
+run: $(target)
+       @echo ============================ Running $<
+       @./$(target)$(EXEEXT)
+
+
+# As a last rule to the `make' program, we define what to do when
+# cleaning up a directory. This usually involves deleting object files
+# and other automatically created files such as the executable itself,
+# backup files, and data files. Since the latter are not usually quite
+# diverse, you needed to declare them at the top of this file.
+clean:
+       -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files)
+
+
+# Since we have not yet stated how to make an object file from a C++
+# file, we should do so now. Since the many flags passed to the
+# compiler are usually not of much interest, we suppress the actual
+# command line using the `at' sign in the first column of the rules
+# and write the string indicating what we do instead.
+./%.g.$(OBJEXT) :
+       @echo ==============debug========= $(<F)
+       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+./%.$(OBJEXT) :
+       @echo ==============optimized===== $(<F)
+       @$(CXX) $(CXXFLAGS.o) -c $< -o $@
+
+
+# The following statement tells make that the rules `run' and `clean'
+# are not expected to produce files of the same name as Makefile rules
+# usually do.
+.PHONY: run clean
+
+
+# Finally there is a rule which you normally need not care much about:
+# since the executable depends on some include files from the library,
+# besides the C++ application file of course, it is necessary to
+# re-generate the executable when one of the files it depends on has
+# changed. The following rule to created a dependency file
+# `Makefile.dep', which `make' uses to determine when to regenerate
+# the executable. This file is automagically remade whenever needed,
+# i.e. whenever one of the cc-/h-files changed. Make detects whether
+# to remake this file upon inclusion at the bottom of this file.
+#
+# If the creation of Makefile.dep fails, blow it away and fail
+Makefile.dep: $(target).cc Makefile \
+              $(shell echo $D/*/include/*/*.h)
+       @echo ============================ Remaking $@
+       @$D/common/scripts/make_dependencies  $(INCLUDE) -B. $(target).cc \
+               > $@ \
+         || (rm -f $@ ; false)
+       @if test -s $@ ; then : else rm $@ ; fi
+
+
+# To make the dependencies known to `make', we finally have to include
+# them:
+include Makefile.dep
+
+
diff --git a/deal.II/examples/step-31/doc/intro.dox b/deal.II/examples/step-31/doc/intro.dox
new file mode 100644 (file)
index 0000000..c59cfdf
--- /dev/null
@@ -0,0 +1,2 @@
+<a name="Intro"></a>
+<h1>Introduction</h1>
diff --git a/deal.II/examples/step-31/doc/results.dox b/deal.II/examples/step-31/doc/results.dox
new file mode 100644 (file)
index 0000000..f4c6fee
--- /dev/null
@@ -0,0 +1 @@
+<h1>Results</h1>
diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc
new file mode 100644 (file)
index 0000000..1fea326
--- /dev/null
@@ -0,0 +1,792 @@
+/* $Id: step-22.cc 15679 2008-01-24 23:28:37Z bangerth $ */
+/* Author: Wolfgang Bangerth, Texas A&M University, 2008 */
+
+/*    $Id: step-22.cc 15679 2008-01-24 23:28:37Z bangerth $       */
+/*    Version: $Name$                                          */
+/*                                                                */
+/*    Copyright (C) 2008 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.                        */
+
+
+
+#include <base/quadrature_lib.h>
+#include <base/logstream.h>
+#include <base/function.h>
+#include <base/utilities.h>
+
+#include <lac/block_vector.h>
+#include <lac/full_matrix.h>
+#include <lac/block_sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <lac/sparse_direct.h>
+
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/grid_tools.h>
+#include <grid/grid_refinement.h>
+
+#include <dofs/dof_handler.h>
+#include <dofs/dof_renumbering.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <dofs/dof_constraints.h>
+
+#include <fe/fe_q.h>
+#include <fe/fe_dgq.h>
+#include <fe/fe_system.h>
+#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+#include <fe/mapping_c1.h>
+
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <numerics/error_estimator.h>
+#include <numerics/solution_transfer.h>
+
+#include <fstream>
+#include <sstream>
+
+using namespace dealii;
+
+
+                                 
+template <int dim>
+class StokesProblem 
+{
+  public:
+    StokesProblem (const unsigned int degree);
+    void run ();
+    
+  private:
+    void setup_dofs ();
+    void assemble_system ();
+    void solve ();
+    void output_results (const unsigned int refinement_cycle) const;
+    void refine_mesh ();
+    
+    const unsigned int   degree;
+    
+    Triangulation<dim>   triangulation;
+    FESystem<dim>        fe;
+    DoFHandler<dim>      dof_handler;
+
+    ConstraintMatrix     hanging_node_constraints;
+    
+    BlockSparsityPattern      sparsity_pattern;
+    BlockSparseMatrix<double> system_matrix;
+
+    BlockVector<double> solution;
+    BlockVector<double> system_rhs;
+
+    boost::shared_ptr<SparseDirectUMFPACK> A_preconditioner;
+};
+
+
+
+
+
+template <int dim>
+class PressureBoundaryValues : public Function<dim> 
+{
+  public:
+    PressureBoundaryValues () : Function<dim>(1) {}
+    
+    virtual double value (const Point<dim>   &p,
+                          const unsigned int  component = 0) const;
+};
+
+
+template <int dim>
+double
+PressureBoundaryValues<dim>::value (const Point<dim>  &/*p*/,
+                                    const unsigned int /*component*/) const 
+{
+  return 0;
+}
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim> 
+{
+  public:
+    BoundaryValues () : Function<dim>(dim+1) {}
+    
+    virtual double value (const Point<dim>   &p,
+                          const unsigned int  component = 0) const;
+
+    virtual void vector_value (const Point<dim> &p, 
+                               Vector<double>   &value) const;
+};
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value (const Point<dim>  &p,
+                           const unsigned int component) const 
+{
+  if (component == 0)
+    return (p[0] < 0 ? -1 : (p[0] > 0 ? 1 : 0));
+  return 0;
+}
+
+
+template <int dim>
+void
+BoundaryValues<dim>::vector_value (const Point<dim> &p,
+                                  Vector<double>   &values) const 
+{
+  for (unsigned int c=0; c<this->n_components; ++c)
+    values(c) = BoundaryValues<dim>::value (p, c);
+}
+
+
+
+template <int dim>
+class RightHandSide : public Function<dim> 
+{
+  public:
+    RightHandSide () : Function<dim>(dim+1) {}
+    
+    virtual double value (const Point<dim>   &p,
+                          const unsigned int  component = 0) const;
+
+    virtual void vector_value (const Point<dim> &p, 
+                               Vector<double>   &value) const;
+};
+
+
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim>  &p,
+                           const unsigned int component) const 
+{
+  return 0;
+}
+
+
+template <int dim>
+void
+RightHandSide<dim>::vector_value (const Point<dim> &p,
+                                  Vector<double>   &values) const 
+{
+  for (unsigned int c=0; c<this->n_components; ++c)
+    values(c) = RightHandSide<dim>::value (p, c);
+}
+
+
+
+
+
+
+template <int dim>
+Tensor<1,dim>
+extract_u (const FEValuesBase<dim> &fe_values,
+           const unsigned int i,
+           const unsigned int q)
+{
+  Tensor<1,dim> tmp;
+
+  const unsigned int component
+    = fe_values.get_fe().system_to_component_index(i).first;
+
+  if (component < dim)
+    tmp[component] = fe_values.shape_value (i,q);
+
+  return tmp;
+}
+
+
+
+template <int dim>
+double
+extract_div_u (const FEValuesBase<dim> &fe_values,
+               const unsigned int i,
+               const unsigned int q)
+{
+  const unsigned int component
+    = fe_values.get_fe().system_to_component_index(i).first;
+
+  if (component < dim)
+    return fe_values.shape_grad (i,q)[component];
+  else
+    return 0;
+}
+
+
+
+template <int dim>
+Tensor<2,dim>
+extract_grad_s_u (const FEValuesBase<dim> &fe_values,
+                 const unsigned int i,
+                 const unsigned int q)
+{
+  Tensor<2,dim> tmp;
+
+  const unsigned int component
+    = fe_values.get_fe().system_to_component_index(i).first;
+  
+  if (component < dim)
+    {
+      const Tensor<1,dim> grad_phi_over_2 = fe_values.shape_grad (i,q) / 2;
+      
+      for (unsigned int e=0; e<dim; ++e)
+       tmp[component][e] += grad_phi_over_2[e];
+      for (unsigned int d=0; d<dim; ++d)
+       tmp[d][component] += grad_phi_over_2[d];
+    }
+  
+  return tmp;
+}
+
+
+  
+template <int dim>
+double extract_p (const FEValuesBase<dim> &fe_values,
+                  const unsigned int i,
+                  const unsigned int q)
+{
+  const unsigned int component
+    = fe_values.get_fe().system_to_component_index(i).first;
+
+  if (component == dim)
+    return fe_values.shape_value (i,q);
+  else
+    return 0;
+}
+
+
+
+
+template <class Matrix, class Preconditioner>
+class InverseMatrix : public Subscriptor
+{
+  public:
+    InverseMatrix (const Matrix         &m,
+                  const Preconditioner &preconditioner);
+
+    void vmult (Vector<double>       &dst,
+                const Vector<double> &src) const;
+
+  private:
+    const SmartPointer<const Matrix> matrix;
+    const Preconditioner &preconditioner;
+
+    mutable GrowingVectorMemory<> vector_memory;    
+};
+
+
+template <class Matrix, class Preconditioner>
+InverseMatrix<Matrix,Preconditioner>::InverseMatrix (const Matrix &m,
+                                                    const Preconditioner &preconditioner)
+                :
+                matrix (&m),
+               preconditioner (preconditioner)
+{}
+
+
+                                 
+template <class Matrix, class Preconditioner>
+void InverseMatrix<Matrix,Preconditioner>::vmult (Vector<double>       &dst,
+                                                 const Vector<double> &src) const
+{
+  SolverControl solver_control (src.size(), 1e-6*src.l2_norm());
+  SolverCG<> cg (solver_control, vector_memory);
+
+  dst = 0;
+
+  try
+    {
+      cg.solve (*matrix, dst, src, preconditioner);
+    }
+  catch (std::exception &e)
+    {
+      Assert (false, ExcMessage(e.what()));
+    }
+}
+
+
+
+template <class Preconditioner>
+class SchurComplement : public Subscriptor
+{
+  public:
+    SchurComplement (const BlockSparseMatrix<double> &A,
+                     const InverseMatrix<SparseMatrix<double>,Preconditioner> &Minv);
+
+    void vmult (Vector<double>       &dst,
+                const Vector<double> &src) const;
+
+  private:
+    const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
+    const SmartPointer<const InverseMatrix<SparseMatrix<double>,Preconditioner> > m_inverse;
+    
+    mutable Vector<double> tmp1, tmp2;
+};
+
+
+
+template <class Preconditioner>
+SchurComplement<Preconditioner>::
+SchurComplement (const BlockSparseMatrix<double> &A,
+                 const InverseMatrix<SparseMatrix<double>,Preconditioner> &Minv)
+                :
+                system_matrix (&A),
+                m_inverse (&Minv),
+                tmp1 (A.block(0,0).m()),
+                tmp2 (A.block(0,0).m())
+{}
+
+
+template <class Preconditioner>
+void SchurComplement<Preconditioner>::vmult (Vector<double>       &dst,
+                                            const Vector<double> &src) const
+{
+  system_matrix->block(0,1).vmult (tmp1, src);
+  m_inverse->vmult (tmp2, tmp1);
+  system_matrix->block(1,0).vmult (dst, tmp2);
+}
+
+
+
+template <int dim>
+StokesProblem<dim>::StokesProblem (const unsigned int degree)
+                :
+                degree (degree),
+               triangulation (Triangulation<dim>::maximum_smoothing),
+                fe (FE_Q<dim>(degree+1), dim,
+                    FE_Q<dim>(degree), 1),
+                dof_handler (triangulation)
+{}
+
+
+
+
+template <int dim>
+void StokesProblem<dim>::setup_dofs ()
+{  
+  dof_handler.distribute_dofs (fe); 
+  DoFRenumbering::component_wise (dof_handler);
+
+  hanging_node_constraints.clear ();
+  DoFTools::make_hanging_node_constraints (dof_handler,
+                                          hanging_node_constraints);
+  hanging_node_constraints.close ();
+
+  std::vector<unsigned int> dofs_per_component (dim+1);
+  DoFTools::count_dofs_per_component (dof_handler, dofs_per_component);  
+  const unsigned int n_u = dofs_per_component[0] * dim,
+                     n_p = dofs_per_component[dim];
+
+  std::cout << "   Number of active cells: "
+            << triangulation.n_active_cells()
+            << std::endl
+            << "   Number of degrees of freedom: "
+            << dof_handler.n_dofs()
+            << " (" << n_u << '+' << n_p <<')'
+            << std::endl;
+  
+  const unsigned int
+    n_couplings = dof_handler.max_couplings_between_dofs();
+
+  system_matrix.clear ();
+      
+  sparsity_pattern.reinit (2,2);
+  sparsity_pattern.block(0,0).reinit (n_u, n_u, n_couplings);
+  sparsity_pattern.block(1,0).reinit (n_p, n_u, n_couplings);
+  sparsity_pattern.block(0,1).reinit (n_u, n_p, n_couplings);
+  sparsity_pattern.block(1,1).reinit (n_p, n_p, n_couplings);
+  
+  sparsity_pattern.collect_sizes();    
+
+  
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+  hanging_node_constraints.condense (sparsity_pattern);
+  sparsity_pattern.compress();
+  
+  system_matrix.reinit (sparsity_pattern);
+                                   
+  solution.reinit (2);
+  solution.block(0).reinit (n_u);
+  solution.block(1).reinit (n_p);
+  solution.collect_sizes ();
+  
+  system_rhs.reinit (2);
+  system_rhs.block(0).reinit (n_u);
+  system_rhs.block(1).reinit (n_p);
+  system_rhs.collect_sizes ();
+}
+
+
+template <int dim>
+double
+scalar_product (const Tensor<2,dim> &a,
+               const Tensor<2,dim> &b)
+{
+  double tmp = 0;
+  for (unsigned int i=0; i<dim; ++i)
+    for (unsigned int j=0; j<dim; ++j)
+      tmp += a[i][j] * b[i][j];
+  return tmp;
+}
+
+
+
+template <int dim>
+void StokesProblem<dim>::assemble_system () 
+{
+  system_matrix=0;
+  system_rhs=0;
+  
+  QGauss<dim>   quadrature_formula(degree+2); 
+  QGauss<dim-1> face_quadrature_formula(degree+2);
+
+  FEValues<dim> fe_values (fe, quadrature_formula,
+                          update_values    |
+                          update_quadrature_points  |
+                          update_JxW_values |
+                          update_gradients);
+  FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula, 
+                                    update_values    | update_normal_vectors |
+                                    update_quadrature_points  | update_JxW_values);
+
+  const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+  
+  const unsigned int   n_q_points      = quadrature_formula.size();
+  const unsigned int   n_face_q_points = face_quadrature_formula.size();
+
+  FullMatrix<double>   local_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       local_rhs (dofs_per_cell);
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+  
+  const PressureBoundaryValues<dim> pressure_boundary_values;
+  
+  std::vector<double>               boundary_values (n_face_q_points);
+  
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    { 
+      fe_values.reinit (cell);
+      local_matrix = 0;
+      local_rhs = 0;
+
+      for (unsigned int q=0; q<n_q_points; ++q)
+       {
+         for (unsigned int i=0; i<dofs_per_cell; ++i)
+           {
+
+             const Tensor<1,dim> phi_i_u      = extract_u (fe_values, i, q);
+
+             const Tensor<2,dim> phi_i_grads_u= extract_grad_s_u (fe_values, i, q);
+             const double        div_phi_i_u  = extract_div_u (fe_values, i, q);
+             const double        phi_i_p      = extract_p (fe_values, i, q);
+            
+             for (unsigned int j=0; j<dofs_per_cell; ++j)
+               {
+                 const Tensor<2,dim> phi_j_grads_u     = extract_grad_s_u (fe_values, j, q);
+                 const double        div_phi_j_u = extract_div_u (fe_values, j, q);
+                 const double        phi_j_p     = extract_p (fe_values, j, q);
+                
+                 local_matrix(i,j) += (scalar_product(phi_i_grads_u, phi_j_grads_u)
+                                       - div_phi_i_u * phi_j_p
+                                       - phi_i_p * div_phi_j_u
+                                       + phi_i_p * phi_j_p)
+                                      * fe_values.JxW(q);     
+               }
+          }
+       }
+      
+
+      for (unsigned int face_no=0;
+           face_no<GeometryInfo<dim>::faces_per_cell;
+           ++face_no)
+        if (cell->at_boundary(face_no))
+          {
+            fe_face_values.reinit (cell, face_no);
+            
+            pressure_boundary_values
+              .value_list (fe_face_values.get_quadrature_points(),
+                           boundary_values);
+
+            for (unsigned int q=0; q<n_face_q_points; ++q) 
+              for (unsigned int i=0; i<dofs_per_cell; ++i)
+                {
+                  const Tensor<1,dim>
+                    phi_i_u = extract_u (fe_face_values, i, q);
+
+                  local_rhs(i) += -(phi_i_u *
+                                    fe_face_values.normal_vector(q) *
+                                    boundary_values[q] *
+                                    fe_face_values.JxW(q));
+                }
+          }
+
+      cell->get_dof_indices (local_dof_indices);
+
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       for (unsigned int j=0; j<dofs_per_cell; ++j)
+         system_matrix.add (local_dof_indices[i],
+                            local_dof_indices[j],
+                            local_matrix(i,j));
+      
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        system_rhs(local_dof_indices[i]) += local_rhs(i);
+    }
+
+  hanging_node_constraints.condense (system_matrix);
+  hanging_node_constraints.condense (system_rhs);  
+
+  {
+    std::map<unsigned int,double> boundary_values;
+    std::vector<bool> component_mask (dim+1, true);
+    component_mask[dim] = false;
+    VectorTools::interpolate_boundary_values (dof_handler,
+                                             3,
+                                             BoundaryValues<dim>(),
+                                             boundary_values,
+                                             component_mask);
+
+    MatrixTools::apply_boundary_values (boundary_values,
+                                       system_matrix,
+                                       solution,
+                                       system_rhs);
+  }
+  
+  std::cout << "   Computing preconditioner..." << std::flush;
+      
+  A_preconditioner
+    = boost::shared_ptr<SparseDirectUMFPACK>(new SparseDirectUMFPACK());
+  A_preconditioner->initialize (system_matrix.block(0,0));
+
+  std::cout << std::endl;
+}
+
+
+
+template <int dim>
+void StokesProblem<dim>::solve () 
+{
+  const InverseMatrix<SparseMatrix<double>,SparseDirectUMFPACK>
+    A_inverse (system_matrix.block(0,0), *A_preconditioner);
+  Vector<double> tmp (solution.block(0).size());
+  Vector<double> schur_rhs (solution.block(1).size());
+  
+  {
+    A_inverse.vmult (tmp, system_rhs.block(0));
+    system_matrix.block(1,0).vmult (schur_rhs, tmp);
+    schur_rhs -= system_rhs.block(1);
+
+    SchurComplement<SparseDirectUMFPACK>
+      schur_complement (system_matrix, A_inverse);
+    
+    SolverControl solver_control (system_matrix.block(0,0).m(),
+                                  1e-6*schur_rhs.l2_norm());
+    SolverCG<>    cg (solver_control);
+    
+    PreconditionSSOR<> preconditioner;
+    preconditioner.initialize (system_matrix.block(1,1), 1.2);
+
+    InverseMatrix<SparseMatrix<double>,PreconditionSSOR<> >
+      m_inverse (system_matrix.block(1,1), preconditioner);
+    
+    try
+      {
+       cg.solve (schur_complement, solution.block(1), schur_rhs,
+                 m_inverse);
+      }
+    catch (...)
+      {
+       abort ();
+      }
+
+                                    // produce a consistent flow field
+    hanging_node_constraints.distribute (solution);
+  
+    std::cout << "   "
+              << solver_control.last_step()
+              << " CG Schur complement iterations for pressure."
+              << std::endl;    
+  }
+
+  {
+    system_matrix.block(0,1).vmult (tmp, solution.block(1));
+    tmp *= -1;
+    tmp += system_rhs.block(0);
+
+    A_inverse.vmult (solution.block(0), tmp);
+
+                                    // produce a consistent pressure field
+    hanging_node_constraints.distribute (solution);
+  }
+}
+                                 
+
+
+template <int dim>
+void
+StokesProblem<dim>::output_results (const unsigned int refinement_cycle)  const
+{
+  std::vector<std::string> solution_names (dim, "velocity");
+  solution_names.push_back ("p");
+  
+  DataOut<dim> data_out;
+
+  data_out.attach_dof_handler (dof_handler);
+
+  std::vector<DataComponentInterpretation::DataComponentInterpretation>
+    data_component_interpretation
+    (dim+1, DataComponentInterpretation::component_is_scalar);
+  for (unsigned int i=0; i<dim; ++i)
+    data_component_interpretation[i]
+      = DataComponentInterpretation::component_is_part_of_vector;
+  
+  data_out.add_data_vector (solution, solution_names,
+                           DataOut<dim>::type_dof_data,
+                           data_component_interpretation);
+  
+  data_out.build_patches ();
+  
+  std::ostringstream filename;
+  filename << "solution-"
+          << Utilities::int_to_string (refinement_cycle, 2)
+          << ".vtk";
+
+  std::ofstream output (filename.str().c_str());
+  data_out.write_vtk (output);
+}
+
+
+
+template <int dim>
+void
+StokesProblem<dim>::refine_mesh () 
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+  std::vector<bool> component_mask (dim+1, false);
+  component_mask[dim] = true;
+  KellyErrorEstimator<dim>::estimate (dof_handler,
+                                     QGauss<dim-1>(degree+1),
+                                     typename FunctionMap<dim>::type(),
+                                     solution,
+                                     estimated_error_per_cell,
+                                     component_mask);
+
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                  estimated_error_per_cell,
+                                                  0.3, 0.0);
+  triangulation.execute_coarsening_and_refinement ();
+}
+
+
+template <int dim>
+void StokesProblem<dim>::run () 
+{
+  switch (dim)
+    {
+      case 2:
+      {
+       std::vector<unsigned int> subdivisions (dim, 1);
+       subdivisions[0] = 4;
+       
+       GridGenerator::subdivided_hyper_rectangle (triangulation,
+                                                  subdivisions,
+                                                  Point<dim>(-2,-1),
+                                                  Point<dim>(2,0),
+                                                  true);
+       
+       triangulation.refine_global (2);
+
+       break;
+      }
+
+      case 3:
+      {
+       GridGenerator::hyper_shell (triangulation,
+                                   Point<dim>(), 0.5, 1.0);
+       
+       static HyperShellBoundary<dim> boundary;
+       triangulation.set_boundary (0, boundary);
+       
+       triangulation.refine_global (2);
+
+       break;
+      }
+
+      default:
+           Assert (false, ExcNotImplemented());
+    }
+  
+
+  for (unsigned int refinement_cycle = 0; refinement_cycle<9;
+       ++refinement_cycle)
+    {
+      std::cout << "Refinement cycle " << refinement_cycle << std::endl;
+      
+      if (refinement_cycle > 0)
+       refine_mesh ();
+      
+      setup_dofs ();
+
+      std::cout << "   Assembling..." << std::endl;
+      assemble_system ();      
+
+      std::cout << "   Solving..." << std::endl;
+      solve ();
+      
+      output_results (refinement_cycle);
+
+      std::cout << std::endl;
+    }
+}
+
+    
+
+int main () 
+{
+  try
+    {
+      deallog.depth_console (0);
+
+      StokesProblem<2> flow_problem(1);
+      flow_problem.run ();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      
+      return 1;
+    }
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+
+  return 0;
+}

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.