]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add example 2.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 6 Dec 1999 09:39:57 +0000 (09:39 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 6 Dec 1999 09:39:57 +0000 (09:39 +0000)
git-svn-id: https://svn.dealii.org/trunk@1980 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/Attic/examples/step-by-step/Makefile
deal.II/deal.II/Attic/examples/step-by-step/step-1/step-1.cc
deal.II/deal.II/Attic/examples/step-by-step/step-2/Makefile [new file with mode: 0644]
deal.II/deal.II/Attic/examples/step-by-step/step-2/step-2.cc [new file with mode: 0644]
deal.II/examples/Makefile
deal.II/examples/step-1/step-1.cc
deal.II/examples/step-2/Makefile [new file with mode: 0644]
deal.II/examples/step-2/step-2.cc [new file with mode: 0644]

index fbf1950096430e2fc769ba0bf63e60b71b2d8fba..9c40c4dc1e605b579f625dc799a7e6e1a7a68cb3 100644 (file)
@@ -4,11 +4,13 @@
 
 default:
        cd step-1 ; $(MAKE)
+       cd step-2 ; $(MAKE)
 
 
 run:
        cd step-1 ; $(MAKE) run
-
+       cd step-2 ; $(MAKE) run
 
 clean:
        cd step-1 ; $(MAKE) clean
+       cd step-2 ; $(MAKE) clean
index 8ffd572b546afb39afe658a7ca3f0f14fa082770..195b23fee1dced710ad17c4e9e86d9433b22f7be 100644 (file)
@@ -32,7 +32,28 @@ void first_grid ()
 {
                                   // Define an object for a
                                   // triangulation of a
-                                  // two-dimensional domain
+                                  // two-dimensional domain. Here and
+                                  // in many following cases, the
+                                  // string "<2>" after a class name
+                                  // indicates that this is an object
+                                  // that shall work in two space
+                                  // dimensions. Likewise, there are
+                                  // version working in one ("<1>")
+                                  // and three ("<3>") space
+                                  // dimension, or for all
+                                  // dimensions. We will see such
+                                  // constructs in later examples,
+                                  // where we show how to program
+                                  // dimension independently.
+                                  // (At present, only one through
+                                  // three space dimensions are
+                                  // supported, but that is not a
+                                  // restriction. In case someone
+                                  // would like to implement four
+                                  // dimensional finite elements, for
+                                  // example for general relativity,
+                                  // this would be a straightforward
+                                  // thing.)
   Triangulation<2> triangulation;
   
                                   // Fill it with a square
diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-2/Makefile b/deal.II/deal.II/Attic/examples/step-by-step/step-2/Makefile
new file mode 100644 (file)
index 0000000..8fb10b1
--- /dev/null
@@ -0,0 +1,120 @@
+# $Id$
+# Copyright W. Bangerth, University of Heidelberg, 1998
+
+# Template for makefiles for the examples subdirectory. In principle,
+# everything should be done automatically if you set the target file
+# here correctly:
+target   = step-2
+
+# All dependencies between files should be updated by the included
+# file Makefile.dep if necessary. Object files are compiled into
+# the archives ./Obj.a and ./Obj.g.a. By default, the debug version
+# is used to link. It you don't like that, change the following
+# variable to "off"
+debug-mode = on
+
+# If you want your program to be linked with extra object or library
+# files, specify them here:
+user-libs =
+
+# To run the program, use "make run"; to give parameters to the program,
+# give the parameters to the following variable:
+run-parameters  = 
+
+# To execute additional action apart from running the program, fill
+# in this list:
+additional-run-action =
+
+# To specify which files are to be deleted by "make clean" (apart from
+# the usual ones: object files, executables, backups, etc), fill in the
+# following list
+delete-files = *gnuplot *inp *history
+
+
+
+
+###############################################################################
+# Internals
+
+#deal include base path
+D = ../../../..
+
+include $D/common/Make.global_options
+
+# get lists of files we need
+cc-files    = $(filter-out *%, $(shell echo *.cc))
+o-files     = $(cc-files:.cc=.o)
+go-files    = $(cc-files:.cc=.go)
+h-files     = $(filter-out *%, $(shell echo *.h))
+lib-h-files = $(filter-out *%, $(shell echo ../../include/*/*.h))
+
+# list of libraries needed to link with
+libs     = ./Obj.a   -ldeal_II_2d  -llac -lbase
+libs.g   = ./Obj.g.a -ldeal_II_2d.g -llac.g -lbase.g
+
+
+# check whether we use debug mode or not
+ifeq ($(debug-mode),on)
+libraries = $(libs.g)
+flags     = $(CXXFLAGS.g)
+endif
+
+ifeq ($(debug-mode),off)
+libraries = $(libs)
+flags     = $(CXXFLAGS)
+endif
+
+
+
+# make rule for the target
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) $(flags) -o $@ $^ $(user-libs)
+
+# rule how to run the program
+run: $(target)
+       $(target) $(run-parameters)
+       $(additional-run-action)
+
+
+# rule to make object files
+%.go : %.cc
+       @echo ============================ Compiling with debugging information:   $<
+       @echo $(CXX) ... -c $< -o $@
+       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc
+       @echo ============================ Compiling with optimization:   $<
+       @echo $(CXX) ... -c $< -o $@
+       @$(CXX) $(CXXFLAGS) -c $< -o $@
+
+
+# rules which files the libraries depend upon
+Obj.a: ./Obj.a($(o-files))
+Obj.g.a: ./Obj.g.a($(go-files))
+
+
+clean:
+       -rm -f *.o *.go *~ Makefile.dep Obj.a Obj.g.a $(target) $(delete-files)
+
+
+
+.PHONY: clean
+
+
+#Rule to generate the dependency file. 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.
+#
+#use perl to generate rules for the .go files as well
+#as to make rules not for tria.o and the like, but
+#rather for libnumerics.a(tria.o)
+Makefile.dep: $(cc-files) $(h-files) $(lib-h-files)
+       @echo ============================ Remaking Makefile
+       @perl $D/common/scripts/Make_dep.pl ./Obj $(INCLUDE) $(cc-files) \
+               > Makefile.dep
+
+
+include Makefile.dep
+
diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-2/step-2.cc b/deal.II/deal.II/Attic/examples/step-by-step/step-2/step-2.cc
new file mode 100644 (file)
index 0000000..c952420
--- /dev/null
@@ -0,0 +1,360 @@
+/* $Id$ */
+
+                                // The following includes are just
+                                // like for the previous program, so
+                                // will not be commented further
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+
+                                // We need this include file for the
+                                // association of degrees of freedom
+                                // ("DoF"s) to vertices, lines, and
+                                // cells.
+#include <grid/dof.h>
+                                // The following include contains the
+                                // description of the bilinear finite
+                                // element, including the facts that
+                                // it has one degree of freedom on
+                                // each vertex of the triangulation,
+                                // but none on faces and none in the
+                                // interior of the cells.
+                                //
+                                // In fact, the file contains the
+                                // description of several more finite
+                                // elements as well, such as
+                                // biquadratic, bicubic and biquartic
+                                // elements, but not only for two
+                                // space dimensions, but also for one
+                                // and three dimensions.
+#include <fe/fe_lib.lagrange.h>
+                                // In the following file, several
+                                // tools for manipulating degrees of
+                                // freedom can be found:
+#include <basic/dof_tools.h>
+                                // We will use a sparse matrix to
+                                // visualize the pattern of nonzero
+                                // entries resulting from the
+                                // distribution of degrees of freedom
+                                // on the grid. That class can be
+                                // found here:
+#include <lac/sparsematrix.h>
+                                // We will want to use a special
+                                // algorithm to renumber degrees of
+                                // freedom. It is declared here:
+#include <numerics/dof_renumbering.h>
+
+                                // This is needed for C++ output:
+#include <fstream>
+
+
+
+                                // This is the function that produced
+                                // the circular grid in the previous
+                                // example. The sole difference is
+                                // that it returns the grid it
+                                // produces via its argument.
+                                //
+                                // We won't comment on the internals
+                                // of this function, since this has
+                                // been done in the previous
+                                // example. If you don't understand
+                                // what is happening here, look
+                                // there.
+void make_grid (Triangulation<2> &triangulation)
+{
+  const Point<2> center (1,0);
+  const double inner_radius = 0.5,
+              outer_radius = 1.0;
+  GridGenerator::hyper_shell (triangulation,
+                             center, inner_radius, outer_radius);
+
+                                  // This is the single difference to
+                                  // the respetive function in the
+                                  // previous program: since we want
+                                  // to export the triangulation
+                                  // through this function's
+                                  // parameter, we need to make sure
+                                  // that the boundary object lives
+                                  // at least as long as the
+                                  // triangulation does. However,
+                                  // since the boundary object is a
+                                  // local variable, it would be
+                                  // deleted at the end of this
+                                  // function, which is too early; by
+                                  // declaring it 'static', we can
+                                  // assure that it lives until the
+                                  // end of the program.
+  static const HyperShellBoundary<2> boundary_description(center);
+  triangulation.set_boundary (0, boundary_description);
+  
+  for (unsigned int step=0; step<5; ++step)
+    {
+      Triangulation<2>::active_cell_iterator cell, endc;
+      cell = triangulation.begin_active();
+      endc = triangulation.end();
+
+      for (; cell!=endc; ++cell)
+       for (unsigned int vertex=0;
+            vertex < GeometryInfo<2>::vertices_per_cell;
+            ++vertex)
+         {
+           const Point<2> vector_to_center
+             = (cell->vertex(vertex) - center);
+           const double distance_from_center
+             = sqrt(vector_to_center.square());
+           
+           if (fabs(distance_from_center - inner_radius) < 1e-10)
+             {
+               cell->set_refine_flag ();
+               break;
+             };
+         };
+
+      triangulation.execute_coarsening_and_refinement ();
+    };
+};
+
+
+                                // up to now, we only have a grid,
+                                // i.e. some geometrical (the
+                                // position of the vertices and which
+                                // vertices make up which cell) and
+                                // some topological information
+                                // (neighborhoods of cells). To use
+                                // numerical algorithms, one needs
+                                // some logic information in addition
+                                // to that: we would like to
+                                // associate degree of freedom
+                                // numbers to each vertex (or line,
+                                // or cell, in case we were using
+                                // higher order elements) to later
+                                // generate matrices and vectors
+                                // which describe a finite element
+                                // field on the triangulation.
+void distribute_dofs (DoFHandler<2> &dof_handler) 
+{
+                                  // In order to associate degrees of
+                                  // freedom with features of a
+                                  // triangulation (vertices, lines,
+                                  // quadrilaterals), we need an
+                                  // object which describes how many
+                                  // degrees of freedom are to be
+                                  // associated to each of these
+                                  // objects. For (bi-, tri-)linear
+                                  // finite elements, this is done
+                                  // using the FEQ1 class, which
+                                  // states that one degree of
+                                  // freedom is to be assigned to
+                                  // each vertex, while there are
+                                  // none on lines and inside the
+                                  // quadrilateral. We first need to
+                                  // create an object of this class
+                                  // and use it to distribute the
+                                  // degrees of freedom. Note that
+                                  // the DoFHandler object will store
+                                  // a reference to this object, so
+                                  // we need to make it static as
+                                  // well, in order to prevent its
+                                  // preemptive
+                                  // destruction. (However, the
+                                  // library would warn us about this
+                                  // and exit the program if that
+                                  // occured. You can check this, if
+                                  // you want, by removing the
+                                  // 'static' declaration.)
+  static const FEQ1<2> finite_element;
+  dof_handler.distribute_dofs (finite_element);
+
+                                  // Now we have associated a number
+                                  // to each vertex, but how can we
+                                  // visualize this? Unfortunately,
+                                  // presently there is no way
+                                  // implemented to directly show the
+                                  // DoF number associated with each
+                                  // vertex. However, such
+                                  // information would hardly ever be
+                                  // truly important, since the
+                                  // numbering itself is more or less
+                                  // arbitrary. There are more
+                                  // important factors, of which we
+                                  // will visualize one in the
+                                  // following.
+                                  //
+                                  // Associated with each vertex of
+                                  // the triangulation is a shape
+                                  // function. Assume we want to
+                                  // solve something like Laplace's
+                                  // equation, then the different
+                                  // matrix entries will be the
+                                  // integrals over the gradient of
+                                  // each two such shape
+                                  // functions. Obviously, since the
+                                  // shape functions are not equal to
+                                  // zero only on the cells adjacent
+                                  // to the vertex they are
+                                  // associated to, matrix entries
+                                  // will be nonzero only of the
+                                  // supports of the shape functions
+                                  // associated to the column and row
+                                  // numbers intersect. This is only
+                                  // the case for adjacent shape
+                                  // functions, and therefore only
+                                  // for adjacent vertices. Now,
+                                  // since the vertices are numbered
+                                  // more or less randomly be the
+                                  // above function
+                                  // (distribute_dofs), the pattern
+                                  // of nonzero entries in the matrix
+                                  // will be somewhat ragged, and we
+                                  // will take a look at it now.
+                                  //
+                                  // First we have to create a
+                                  // structure which we use to store
+                                  // the places of nonzero
+                                  // elements. We have to give it the
+                                  // size of the matrix, which in our
+                                  // case will be square with that
+                                  // many rows and columns as there
+                                  // are degrees of freedom on the
+                                  // grid:
+  SparseMatrixStruct sparsity_pattern (dof_handler.n_dofs(),
+                                      dof_handler.n_dofs());
+                                  // We fill it with the places where
+                                  // nonzero elements will be located
+                                  // given the present numbering of
+                                  // degrees of freedom:
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+                                  // Before further work can be done
+                                  // on the object, we have to allow
+                                  // for some internal
+                                  // reorganization:
+  sparsity_pattern.compress ();
+
+                                  // Now write the results to a file
+  ofstream out ("sparsity_pattern.1");
+  sparsity_pattern.print_gnuplot (out);
+                                  // The result is in GNUPLOT format,
+                                  // where in each line of the output
+                                  // file, the coordinates of one
+                                  // nonzero entry are listed. The
+                                  // output will be shown below.
+                                  //
+                                  // If you look at it, you will note
+                                  // that the sparsity pattern is
+                                  // symmetric, which is quite often
+                                  // so, unless you have a rather
+                                  // special equation you want to
+                                  // solve. You will also note that
+                                  // it has several distinct region,
+                                  // which stem from the fact that
+                                  // the numbering starts from the
+                                  // coarsest cells and moves on to
+                                  // the finer ones; since they are
+                                  // all distributed symmetrically
+                                  // around the origin, this shows up
+                                  // again in the sparsity pattern.
+};
+
+
+
+                                // In the sparsity pattern produced
+                                // above, the nonzero entries
+                                // extended quite far off from the
+                                // diagonal. For some algorithms,
+                                // this is unfavorable, and we will
+                                // show a simple way how to improve
+                                // this situation.
+                                //
+                                // Remember that for an entry (i,j)
+                                // in the matrix to be nonzero, the
+                                // supports of the shape functions i
+                                // and j needed to intersect
+                                // (otherwise in the integral, the
+                                // integrand would be zero everywhere
+                                // since either the one or the other
+                                // shape function is zero at some
+                                // point). However, the supports of
+                                // shape functions intersected only
+                                // of they were adjacent to each
+                                // other, so in order to have the
+                                // nonzero entries clustered around
+                                // the diagonal (where i equals j),
+                                // we would like to have adjacent
+                                // shape functions to be numbered
+                                // with indices (DoF numbers) that
+                                // differ not too much.
+                                //
+                                // This can be accomplished by a
+                                // simple front marching algorithm,
+                                // where one starts at a given vertex
+                                // and gives it the index zero. Then,
+                                // its neighbors are numbered
+                                // successively, making their indices
+                                // close to the original one. Then,
+                                // their neighbors, if not yet
+                                // numbered, are numbered, and so
+                                // on. One such algorithm is the one
+                                // by Cuthill and McKee, which is a
+                                // little more complicated, but works
+                                // along the same lines. We will use
+                                // it to renumber the degrees of
+                                // freedom such that the resulting
+                                // sparsity pattern is more localized
+                                // around the diagonal.
+void renumber_dofs (DoFHandler<2> &dof_handler) 
+{
+                                  // Renumber the degrees of freedom...
+  DoFRenumbering::Cuthill_McKee (dof_handler);
+                                  // ...regenerate the sparsity pattern...
+  SparseMatrixStruct sparsity_pattern (dof_handler.n_dofs(),
+                                      dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+  sparsity_pattern.compress ();
+                                  // ...and output the result:
+  ofstream out ("sparsity_pattern.2");
+  sparsity_pattern.print_gnuplot (out);
+                                  // Again, the output is shown
+                                  // below. Note that the nonzero
+                                  // entries are clustered far better
+                                  // around the diagonal than
+                                  // before. This effect is even more
+                                  // distinguished for larger
+                                  // matrices (the present one has
+                                  // 1260 rows and columns, but large
+                                  // matrices often have several
+                                  // 100,000s).
+};
+
+
+
+
+                                // This is the main program, which
+                                // only calls the other functions in
+                                // their respective order.
+int main () 
+{
+                                  // Allocate space for a triangulation...
+  Triangulation<2> triangulation;
+                                  // ...and create it
+  make_grid (triangulation);
+
+                                  // A variable that will hold the
+                                  // information which vertex has
+                                  // which number. The geometric
+                                  // information is passed as
+                                  // parameter and a pointer to the
+                                  // triangulation will be stored
+                                  // inside the DoFHandler object.
+  DoFHandler<2> dof_handler (triangulation);
+                                  // Associate vertices and degrees
+                                  // of freedom.
+  distribute_dofs (dof_handler);
+
+                                  // Show the effect of renumbering
+                                  // of degrees of freedom to the
+                                  // sparsity pattern of the matrix.
+  renumber_dofs (dof_handler);
+};
index fbf1950096430e2fc769ba0bf63e60b71b2d8fba..9c40c4dc1e605b579f625dc799a7e6e1a7a68cb3 100644 (file)
@@ -4,11 +4,13 @@
 
 default:
        cd step-1 ; $(MAKE)
+       cd step-2 ; $(MAKE)
 
 
 run:
        cd step-1 ; $(MAKE) run
-
+       cd step-2 ; $(MAKE) run
 
 clean:
        cd step-1 ; $(MAKE) clean
+       cd step-2 ; $(MAKE) clean
index 8ffd572b546afb39afe658a7ca3f0f14fa082770..195b23fee1dced710ad17c4e9e86d9433b22f7be 100644 (file)
@@ -32,7 +32,28 @@ void first_grid ()
 {
                                   // Define an object for a
                                   // triangulation of a
-                                  // two-dimensional domain
+                                  // two-dimensional domain. Here and
+                                  // in many following cases, the
+                                  // string "<2>" after a class name
+                                  // indicates that this is an object
+                                  // that shall work in two space
+                                  // dimensions. Likewise, there are
+                                  // version working in one ("<1>")
+                                  // and three ("<3>") space
+                                  // dimension, or for all
+                                  // dimensions. We will see such
+                                  // constructs in later examples,
+                                  // where we show how to program
+                                  // dimension independently.
+                                  // (At present, only one through
+                                  // three space dimensions are
+                                  // supported, but that is not a
+                                  // restriction. In case someone
+                                  // would like to implement four
+                                  // dimensional finite elements, for
+                                  // example for general relativity,
+                                  // this would be a straightforward
+                                  // thing.)
   Triangulation<2> triangulation;
   
                                   // Fill it with a square
diff --git a/deal.II/examples/step-2/Makefile b/deal.II/examples/step-2/Makefile
new file mode 100644 (file)
index 0000000..8fb10b1
--- /dev/null
@@ -0,0 +1,120 @@
+# $Id$
+# Copyright W. Bangerth, University of Heidelberg, 1998
+
+# Template for makefiles for the examples subdirectory. In principle,
+# everything should be done automatically if you set the target file
+# here correctly:
+target   = step-2
+
+# All dependencies between files should be updated by the included
+# file Makefile.dep if necessary. Object files are compiled into
+# the archives ./Obj.a and ./Obj.g.a. By default, the debug version
+# is used to link. It you don't like that, change the following
+# variable to "off"
+debug-mode = on
+
+# If you want your program to be linked with extra object or library
+# files, specify them here:
+user-libs =
+
+# To run the program, use "make run"; to give parameters to the program,
+# give the parameters to the following variable:
+run-parameters  = 
+
+# To execute additional action apart from running the program, fill
+# in this list:
+additional-run-action =
+
+# To specify which files are to be deleted by "make clean" (apart from
+# the usual ones: object files, executables, backups, etc), fill in the
+# following list
+delete-files = *gnuplot *inp *history
+
+
+
+
+###############################################################################
+# Internals
+
+#deal include base path
+D = ../../../..
+
+include $D/common/Make.global_options
+
+# get lists of files we need
+cc-files    = $(filter-out *%, $(shell echo *.cc))
+o-files     = $(cc-files:.cc=.o)
+go-files    = $(cc-files:.cc=.go)
+h-files     = $(filter-out *%, $(shell echo *.h))
+lib-h-files = $(filter-out *%, $(shell echo ../../include/*/*.h))
+
+# list of libraries needed to link with
+libs     = ./Obj.a   -ldeal_II_2d  -llac -lbase
+libs.g   = ./Obj.g.a -ldeal_II_2d.g -llac.g -lbase.g
+
+
+# check whether we use debug mode or not
+ifeq ($(debug-mode),on)
+libraries = $(libs.g)
+flags     = $(CXXFLAGS.g)
+endif
+
+ifeq ($(debug-mode),off)
+libraries = $(libs)
+flags     = $(CXXFLAGS)
+endif
+
+
+
+# make rule for the target
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) $(flags) -o $@ $^ $(user-libs)
+
+# rule how to run the program
+run: $(target)
+       $(target) $(run-parameters)
+       $(additional-run-action)
+
+
+# rule to make object files
+%.go : %.cc
+       @echo ============================ Compiling with debugging information:   $<
+       @echo $(CXX) ... -c $< -o $@
+       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc
+       @echo ============================ Compiling with optimization:   $<
+       @echo $(CXX) ... -c $< -o $@
+       @$(CXX) $(CXXFLAGS) -c $< -o $@
+
+
+# rules which files the libraries depend upon
+Obj.a: ./Obj.a($(o-files))
+Obj.g.a: ./Obj.g.a($(go-files))
+
+
+clean:
+       -rm -f *.o *.go *~ Makefile.dep Obj.a Obj.g.a $(target) $(delete-files)
+
+
+
+.PHONY: clean
+
+
+#Rule to generate the dependency file. 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.
+#
+#use perl to generate rules for the .go files as well
+#as to make rules not for tria.o and the like, but
+#rather for libnumerics.a(tria.o)
+Makefile.dep: $(cc-files) $(h-files) $(lib-h-files)
+       @echo ============================ Remaking Makefile
+       @perl $D/common/scripts/Make_dep.pl ./Obj $(INCLUDE) $(cc-files) \
+               > Makefile.dep
+
+
+include Makefile.dep
+
diff --git a/deal.II/examples/step-2/step-2.cc b/deal.II/examples/step-2/step-2.cc
new file mode 100644 (file)
index 0000000..c952420
--- /dev/null
@@ -0,0 +1,360 @@
+/* $Id$ */
+
+                                // The following includes are just
+                                // like for the previous program, so
+                                // will not be commented further
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+
+                                // We need this include file for the
+                                // association of degrees of freedom
+                                // ("DoF"s) to vertices, lines, and
+                                // cells.
+#include <grid/dof.h>
+                                // The following include contains the
+                                // description of the bilinear finite
+                                // element, including the facts that
+                                // it has one degree of freedom on
+                                // each vertex of the triangulation,
+                                // but none on faces and none in the
+                                // interior of the cells.
+                                //
+                                // In fact, the file contains the
+                                // description of several more finite
+                                // elements as well, such as
+                                // biquadratic, bicubic and biquartic
+                                // elements, but not only for two
+                                // space dimensions, but also for one
+                                // and three dimensions.
+#include <fe/fe_lib.lagrange.h>
+                                // In the following file, several
+                                // tools for manipulating degrees of
+                                // freedom can be found:
+#include <basic/dof_tools.h>
+                                // We will use a sparse matrix to
+                                // visualize the pattern of nonzero
+                                // entries resulting from the
+                                // distribution of degrees of freedom
+                                // on the grid. That class can be
+                                // found here:
+#include <lac/sparsematrix.h>
+                                // We will want to use a special
+                                // algorithm to renumber degrees of
+                                // freedom. It is declared here:
+#include <numerics/dof_renumbering.h>
+
+                                // This is needed for C++ output:
+#include <fstream>
+
+
+
+                                // This is the function that produced
+                                // the circular grid in the previous
+                                // example. The sole difference is
+                                // that it returns the grid it
+                                // produces via its argument.
+                                //
+                                // We won't comment on the internals
+                                // of this function, since this has
+                                // been done in the previous
+                                // example. If you don't understand
+                                // what is happening here, look
+                                // there.
+void make_grid (Triangulation<2> &triangulation)
+{
+  const Point<2> center (1,0);
+  const double inner_radius = 0.5,
+              outer_radius = 1.0;
+  GridGenerator::hyper_shell (triangulation,
+                             center, inner_radius, outer_radius);
+
+                                  // This is the single difference to
+                                  // the respetive function in the
+                                  // previous program: since we want
+                                  // to export the triangulation
+                                  // through this function's
+                                  // parameter, we need to make sure
+                                  // that the boundary object lives
+                                  // at least as long as the
+                                  // triangulation does. However,
+                                  // since the boundary object is a
+                                  // local variable, it would be
+                                  // deleted at the end of this
+                                  // function, which is too early; by
+                                  // declaring it 'static', we can
+                                  // assure that it lives until the
+                                  // end of the program.
+  static const HyperShellBoundary<2> boundary_description(center);
+  triangulation.set_boundary (0, boundary_description);
+  
+  for (unsigned int step=0; step<5; ++step)
+    {
+      Triangulation<2>::active_cell_iterator cell, endc;
+      cell = triangulation.begin_active();
+      endc = triangulation.end();
+
+      for (; cell!=endc; ++cell)
+       for (unsigned int vertex=0;
+            vertex < GeometryInfo<2>::vertices_per_cell;
+            ++vertex)
+         {
+           const Point<2> vector_to_center
+             = (cell->vertex(vertex) - center);
+           const double distance_from_center
+             = sqrt(vector_to_center.square());
+           
+           if (fabs(distance_from_center - inner_radius) < 1e-10)
+             {
+               cell->set_refine_flag ();
+               break;
+             };
+         };
+
+      triangulation.execute_coarsening_and_refinement ();
+    };
+};
+
+
+                                // up to now, we only have a grid,
+                                // i.e. some geometrical (the
+                                // position of the vertices and which
+                                // vertices make up which cell) and
+                                // some topological information
+                                // (neighborhoods of cells). To use
+                                // numerical algorithms, one needs
+                                // some logic information in addition
+                                // to that: we would like to
+                                // associate degree of freedom
+                                // numbers to each vertex (or line,
+                                // or cell, in case we were using
+                                // higher order elements) to later
+                                // generate matrices and vectors
+                                // which describe a finite element
+                                // field on the triangulation.
+void distribute_dofs (DoFHandler<2> &dof_handler) 
+{
+                                  // In order to associate degrees of
+                                  // freedom with features of a
+                                  // triangulation (vertices, lines,
+                                  // quadrilaterals), we need an
+                                  // object which describes how many
+                                  // degrees of freedom are to be
+                                  // associated to each of these
+                                  // objects. For (bi-, tri-)linear
+                                  // finite elements, this is done
+                                  // using the FEQ1 class, which
+                                  // states that one degree of
+                                  // freedom is to be assigned to
+                                  // each vertex, while there are
+                                  // none on lines and inside the
+                                  // quadrilateral. We first need to
+                                  // create an object of this class
+                                  // and use it to distribute the
+                                  // degrees of freedom. Note that
+                                  // the DoFHandler object will store
+                                  // a reference to this object, so
+                                  // we need to make it static as
+                                  // well, in order to prevent its
+                                  // preemptive
+                                  // destruction. (However, the
+                                  // library would warn us about this
+                                  // and exit the program if that
+                                  // occured. You can check this, if
+                                  // you want, by removing the
+                                  // 'static' declaration.)
+  static const FEQ1<2> finite_element;
+  dof_handler.distribute_dofs (finite_element);
+
+                                  // Now we have associated a number
+                                  // to each vertex, but how can we
+                                  // visualize this? Unfortunately,
+                                  // presently there is no way
+                                  // implemented to directly show the
+                                  // DoF number associated with each
+                                  // vertex. However, such
+                                  // information would hardly ever be
+                                  // truly important, since the
+                                  // numbering itself is more or less
+                                  // arbitrary. There are more
+                                  // important factors, of which we
+                                  // will visualize one in the
+                                  // following.
+                                  //
+                                  // Associated with each vertex of
+                                  // the triangulation is a shape
+                                  // function. Assume we want to
+                                  // solve something like Laplace's
+                                  // equation, then the different
+                                  // matrix entries will be the
+                                  // integrals over the gradient of
+                                  // each two such shape
+                                  // functions. Obviously, since the
+                                  // shape functions are not equal to
+                                  // zero only on the cells adjacent
+                                  // to the vertex they are
+                                  // associated to, matrix entries
+                                  // will be nonzero only of the
+                                  // supports of the shape functions
+                                  // associated to the column and row
+                                  // numbers intersect. This is only
+                                  // the case for adjacent shape
+                                  // functions, and therefore only
+                                  // for adjacent vertices. Now,
+                                  // since the vertices are numbered
+                                  // more or less randomly be the
+                                  // above function
+                                  // (distribute_dofs), the pattern
+                                  // of nonzero entries in the matrix
+                                  // will be somewhat ragged, and we
+                                  // will take a look at it now.
+                                  //
+                                  // First we have to create a
+                                  // structure which we use to store
+                                  // the places of nonzero
+                                  // elements. We have to give it the
+                                  // size of the matrix, which in our
+                                  // case will be square with that
+                                  // many rows and columns as there
+                                  // are degrees of freedom on the
+                                  // grid:
+  SparseMatrixStruct sparsity_pattern (dof_handler.n_dofs(),
+                                      dof_handler.n_dofs());
+                                  // We fill it with the places where
+                                  // nonzero elements will be located
+                                  // given the present numbering of
+                                  // degrees of freedom:
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+                                  // Before further work can be done
+                                  // on the object, we have to allow
+                                  // for some internal
+                                  // reorganization:
+  sparsity_pattern.compress ();
+
+                                  // Now write the results to a file
+  ofstream out ("sparsity_pattern.1");
+  sparsity_pattern.print_gnuplot (out);
+                                  // The result is in GNUPLOT format,
+                                  // where in each line of the output
+                                  // file, the coordinates of one
+                                  // nonzero entry are listed. The
+                                  // output will be shown below.
+                                  //
+                                  // If you look at it, you will note
+                                  // that the sparsity pattern is
+                                  // symmetric, which is quite often
+                                  // so, unless you have a rather
+                                  // special equation you want to
+                                  // solve. You will also note that
+                                  // it has several distinct region,
+                                  // which stem from the fact that
+                                  // the numbering starts from the
+                                  // coarsest cells and moves on to
+                                  // the finer ones; since they are
+                                  // all distributed symmetrically
+                                  // around the origin, this shows up
+                                  // again in the sparsity pattern.
+};
+
+
+
+                                // In the sparsity pattern produced
+                                // above, the nonzero entries
+                                // extended quite far off from the
+                                // diagonal. For some algorithms,
+                                // this is unfavorable, and we will
+                                // show a simple way how to improve
+                                // this situation.
+                                //
+                                // Remember that for an entry (i,j)
+                                // in the matrix to be nonzero, the
+                                // supports of the shape functions i
+                                // and j needed to intersect
+                                // (otherwise in the integral, the
+                                // integrand would be zero everywhere
+                                // since either the one or the other
+                                // shape function is zero at some
+                                // point). However, the supports of
+                                // shape functions intersected only
+                                // of they were adjacent to each
+                                // other, so in order to have the
+                                // nonzero entries clustered around
+                                // the diagonal (where i equals j),
+                                // we would like to have adjacent
+                                // shape functions to be numbered
+                                // with indices (DoF numbers) that
+                                // differ not too much.
+                                //
+                                // This can be accomplished by a
+                                // simple front marching algorithm,
+                                // where one starts at a given vertex
+                                // and gives it the index zero. Then,
+                                // its neighbors are numbered
+                                // successively, making their indices
+                                // close to the original one. Then,
+                                // their neighbors, if not yet
+                                // numbered, are numbered, and so
+                                // on. One such algorithm is the one
+                                // by Cuthill and McKee, which is a
+                                // little more complicated, but works
+                                // along the same lines. We will use
+                                // it to renumber the degrees of
+                                // freedom such that the resulting
+                                // sparsity pattern is more localized
+                                // around the diagonal.
+void renumber_dofs (DoFHandler<2> &dof_handler) 
+{
+                                  // Renumber the degrees of freedom...
+  DoFRenumbering::Cuthill_McKee (dof_handler);
+                                  // ...regenerate the sparsity pattern...
+  SparseMatrixStruct sparsity_pattern (dof_handler.n_dofs(),
+                                      dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+  sparsity_pattern.compress ();
+                                  // ...and output the result:
+  ofstream out ("sparsity_pattern.2");
+  sparsity_pattern.print_gnuplot (out);
+                                  // Again, the output is shown
+                                  // below. Note that the nonzero
+                                  // entries are clustered far better
+                                  // around the diagonal than
+                                  // before. This effect is even more
+                                  // distinguished for larger
+                                  // matrices (the present one has
+                                  // 1260 rows and columns, but large
+                                  // matrices often have several
+                                  // 100,000s).
+};
+
+
+
+
+                                // This is the main program, which
+                                // only calls the other functions in
+                                // their respective order.
+int main () 
+{
+                                  // Allocate space for a triangulation...
+  Triangulation<2> triangulation;
+                                  // ...and create it
+  make_grid (triangulation);
+
+                                  // A variable that will hold the
+                                  // information which vertex has
+                                  // which number. The geometric
+                                  // information is passed as
+                                  // parameter and a pointer to the
+                                  // triangulation will be stored
+                                  // inside the DoFHandler object.
+  DoFHandler<2> dof_handler (triangulation);
+                                  // Associate vertices and degrees
+                                  // of freedom.
+  distribute_dofs (dof_handler);
+
+                                  // Show the effect of renumbering
+                                  // of degrees of freedom to the
+                                  // sparsity pattern of the matrix.
+  renumber_dofs (dof_handler);
+};

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.