From 7825f61ef992dfebdddecbe836a664c7bfd974a5 Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 6 Dec 1999 09:39:57 +0000 Subject: [PATCH] Add example 2. git-svn-id: https://svn.dealii.org/trunk@1980 0785d39b-7218-0410-832d-ea1e28bc413d --- .../Attic/examples/step-by-step/Makefile | 4 +- .../examples/step-by-step/step-1/step-1.cc | 23 +- .../examples/step-by-step/step-2/Makefile | 120 ++++++ .../examples/step-by-step/step-2/step-2.cc | 360 ++++++++++++++++++ deal.II/examples/Makefile | 4 +- deal.II/examples/step-1/step-1.cc | 23 +- deal.II/examples/step-2/Makefile | 120 ++++++ deal.II/examples/step-2/step-2.cc | 360 ++++++++++++++++++ 8 files changed, 1010 insertions(+), 4 deletions(-) create mode 100644 deal.II/deal.II/Attic/examples/step-by-step/step-2/Makefile create mode 100644 deal.II/deal.II/Attic/examples/step-by-step/step-2/step-2.cc create mode 100644 deal.II/examples/step-2/Makefile create mode 100644 deal.II/examples/step-2/step-2.cc diff --git a/deal.II/deal.II/Attic/examples/step-by-step/Makefile b/deal.II/deal.II/Attic/examples/step-by-step/Makefile index fbf1950096..9c40c4dc1e 100644 --- a/deal.II/deal.II/Attic/examples/step-by-step/Makefile +++ b/deal.II/deal.II/Attic/examples/step-by-step/Makefile @@ -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 diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-1/step-1.cc b/deal.II/deal.II/Attic/examples/step-by-step/step-1/step-1.cc index 8ffd572b54..195b23fee1 100644 --- a/deal.II/deal.II/Attic/examples/step-by-step/step-1/step-1.cc +++ b/deal.II/deal.II/Attic/examples/step-by-step/step-1/step-1.cc @@ -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 index 0000000000..8fb10b1fec --- /dev/null +++ b/deal.II/deal.II/Attic/examples/step-by-step/step-2/Makefile @@ -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 index 0000000000..c952420de8 --- /dev/null +++ b/deal.II/deal.II/Attic/examples/step-by-step/step-2/step-2.cc @@ -0,0 +1,360 @@ +/* $Id$ */ + + // The following includes are just + // like for the previous program, so + // will not be commented further +#include +#include +#include +#include +#include + + // We need this include file for the + // association of degrees of freedom + // ("DoF"s) to vertices, lines, and + // cells. +#include + // 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 + // In the following file, several + // tools for manipulating degrees of + // freedom can be found: +#include + // 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 + // We will want to use a special + // algorithm to renumber degrees of + // freedom. It is declared here: +#include + + // This is needed for C++ output: +#include + + + + // 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); +}; diff --git a/deal.II/examples/Makefile b/deal.II/examples/Makefile index fbf1950096..9c40c4dc1e 100644 --- a/deal.II/examples/Makefile +++ b/deal.II/examples/Makefile @@ -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 diff --git a/deal.II/examples/step-1/step-1.cc b/deal.II/examples/step-1/step-1.cc index 8ffd572b54..195b23fee1 100644 --- a/deal.II/examples/step-1/step-1.cc +++ b/deal.II/examples/step-1/step-1.cc @@ -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 index 0000000000..8fb10b1fec --- /dev/null +++ b/deal.II/examples/step-2/Makefile @@ -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 index 0000000000..c952420de8 --- /dev/null +++ b/deal.II/examples/step-2/step-2.cc @@ -0,0 +1,360 @@ +/* $Id$ */ + + // The following includes are just + // like for the previous program, so + // will not be commented further +#include +#include +#include +#include +#include + + // We need this include file for the + // association of degrees of freedom + // ("DoF"s) to vertices, lines, and + // cells. +#include + // 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 + // In the following file, several + // tools for manipulating degrees of + // freedom can be found: +#include + // 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 + // We will want to use a special + // algorithm to renumber degrees of + // freedom. It is declared here: +#include + + // This is needed for C++ output: +#include + + + + // 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); +}; -- 2.39.5