From: wolf Date: Sun, 28 Nov 1999 19:47:59 +0000 (+0000) Subject: Various updates. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6620f83212e8b620dcd09f92d6726e920de8a627;p=dealii-svn.git Various updates. git-svn-id: https://svn.dealii.org/trunk@1954 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/tutorial/chapter-1.elements/dofs.html b/deal.II/doc/tutorial/chapter-1.elements/dofs.html index b86fc91d26..e5fb3575af 100644 --- a/deal.II/doc/tutorial/chapter-1.elements/dofs.html +++ b/deal.II/doc/tutorial/chapter-1.elements/dofs.html @@ -20,9 +20,12 @@

Degrees of Freedom

-Degrees of freedom - short DoFs -are virtually the most important entities you will -encounter in deal.II. Most of the task of setting up -your problem is accomplished by manipulating them. +Degrees of freedom - short DoFs - are needed to associate a numering +to the shape functions living on our triangulation. Using this +numbering, we can then interpret the discretized operator as a matrix, +where the matrix row and column numbers are identified with the DoF +numbers, and likewise the index of a vector denotes the number of the +respective degree of freedom.

For now we will only discuss how and when to number and renumber the @@ -49,7 +52,7 @@ afterwards the structure will be invalid.

How to Distribute and Renumber Degrees of Freedom

The degrees of freedom are distributed on a triangulation using the function
-void DoFHandler::distribute_dofs(const FiniteElement<dim> &)
+void DoFHandler<dim>::distribute_dofs(const FiniteElement<dim> &)
The argument to distribute_dofs is a finite element that describes how many degrees of freedom are located on vertices, lines etc. (FiniteElement objects have much more functionality than only storing @@ -74,42 +77,53 @@ The finite elements used are bilinear. #include <grid/grid_generator.h> #include <fe/fe_lib.lagrange.h> -Triangulation<dim> tr; -DoFHandler<dim> dof; -FEQ1<dim> fe; +// first build up the triangulation +Triangulation<dim> triangulation; +GridGenerator::hyper_cube (triangulation,-1,1); + +// you may want to refine the triangulation a bit here +... -GridGenerator<dim>::hyper_cube(tr,-1,1); -dof.distribute_dofs(fe); +// then associate a DoFHandler object with this triangulation... +DoFHandler<dim> dof_handler (&tria); + +// ...and distribute degrees of freedom for linear elements +FEQ1<dim> fe; +dof_handler.distribute_dofs(fe);

-For a number of solution algorithms this way of numbering is way below optimal. +For a number of solution algorithms this way of numbering is way below +optimal. In particular, as mentioned above, the degrees of freedom are +numbered starting on the active cells of coarser levels and going to +finer levels, but depending on how your triangulation was refined, the +order in which the cells are traversed will have no geometrical +pattern. Likewise, the distribution of DoF numbers will look mostly +random. +

+

+For some problems, it is better if the numbering of DoFs is roughly +into a certain direction, for example from inflow to outflow in +hyperbolic problems. To achieve a kind of numbering better suited to your problem the deal.II -offers the function
-void DoFHandler::renumber_dofs(const RenumberingMethod method,const bool use_constraints=false, const vector<int> &starting_points=vector<int>())
-or alternatively, the easier to use
-void DoFHandler::renumber_dofs(const RenumberingMethod method)
-The methods available are Cuthill_McKee and reverse_Cuthill_McKee. Both algorithms are usually better than the one used -by distribute_dofs and neither is optimal (there is no algorithm +offers the class +DoFRenumbering with several functions implementing +different ways of renumbering. +The methods presently available are Cuthill-McKee and reverse +Cuthill-McKee, and some other methods you will only want to ask for in +more advanced applications. The first two algorithms are usually much +better than the one used +by distribute_dofs, but neither is optimal (there is no algorithm that accomplishes optimal results within a reasonable amount of time). -Both algorithms require a good index to start with to achieve good results. -Whenever you know a good set of starting indices you might wish to use the -first way of calling renumber_dofs, -in every other case the second way is more advisable. -

-Example: We use the definitions from the -first example given above. We renumber our -degrees of freedom with the -Cuthill_McKee method. +Both algorithms require a good point to start numbering from (for +example one or more degrees of freedom on the inflow boundary) to +achieve good results.

-
-
-dof.renumber_dofs(Cuthill_McKee);
-
-
- -

Access to triangulation cells, faces and vertices

+

+The documentation of the DoFRenumbering gives a good +overview of how the methods mentioned above work, and how the +functions performing the re-enumerations have to be called. diff --git a/deal.II/doc/tutorial/chapter-1.elements/grid_creation.html b/deal.II/doc/tutorial/chapter-1.elements/grid_creation.html index d62a59f8f9..2723502a42 100644 --- a/deal.II/doc/tutorial/chapter-1.elements/grid_creation.html +++ b/deal.II/doc/tutorial/chapter-1.elements/grid_creation.html @@ -17,30 +17,40 @@

Different kinds of grids

-All numerics are done on a grid. +The numerics for which this library was written is always done on a grid. Choosing the right kind of grid can be essential for solving your problem. A grid should be chosen to fit your problem space in the best way possible. This chapter describes coarse grids for the triangulation of your domain.

-

Please remember throughout this chapter: -A triangulation must be void before it is initialized ! -

-deal.II offers three fundamental grid types: +deal.II offers various ways to produce grids for +standard domains, such as a hypercube, a hyperball -and a hyper_L. Furthermore it is possible to +and a hyper_L. There may be more in the meantime (for +example I can now see a ring domain), and you may want to check out +the documentation of the GridGenerator class for this. +

+Furthermore it is possible to read a triangulation from a ucd-file or generate a custom grid from within a program.

+

We note before actually going to some code snippets that +a triangulation object must be empty before it is initialized with a +grid. However, you need not fear now, since the library will tell you +if you try to initialize a triangulation which is already initialized; +in general, the library will warn you about almost all errors you do, +since in the debug version there are very many checks. +

Hypercube

A hypercube can be created using the function
-void GridGenerator<dim>::hyper_cube(Triangulation<dim> &tria,const double left=0.,const double right=1.)
+void GridGenerator::hyper_cube(Triangulation<dim> &tria,const double left=0.,const double right=1.)
The cube created is the tensor product of the left and right edges which default to 0 and 1, creating the unit hypercube. The hypercube will consist of -exactly one cell. In two dimensions, this amounts to a unit square, in three, to a unit cube. +exactly one cell. In two dimensions, this amounts to a unit square, in +three, to a unit cube. In one dimension, you simply get a line.

@@ -60,10 +70,10 @@ function calls needed. Be sure to use them in their appropriate places. #include <grid/tria.h> #include <grid/grid_generator.h> -const unsigned int dim=2; // Two dimensions; to create a cube set to three -Triangulation<dim> tr; +const unsigned int dim=2; // Two dimensions; to create a cube set to dim=3 +Triangulation<dim> triangulation; -GridGenerator<dim>::hyper_cube(tr,-1,1); +GridGenerator::hyper_cube (triangulation, -1, 1); @@ -71,7 +81,7 @@ GridGenerator<dim>::hyper_cube(tr,-1,1);

Hyperball

A hyperball can be created using the function
-void GridGenerator<dim>::hyper_ball(Triangulation<dim> &tria,const Point<dim> center=0.,const double radius=1.)
+void GridGenerator::hyper_ball(Triangulation<dim> &tria,const Point<dim> center=0.,const double radius=1.)
This will create a hyperball of given centre and radius where the location of the centre defaults to the origin and the radius to unity.

@@ -95,14 +105,21 @@ This example will create a hyperball with unit radius centred on (1,0). #include <base/point.h> const unsigned int dim=2; // For example -Triangulation<dim> tr; +Triangulation<dim> triangulation; Point<dim> centre(1,0); // Taking (1,0) as the centre of the ball -GridGenerator<dim>::hyper_ball(tr,centre,1); +GridGenerator::hyper_ball (triangulation, centre, 1); +As can be seen above, the unrefined ball (or circle) does not look +very much like a ball (or circle). After several refinement +steps, it looks much better, though. However, along the rays pointing +outward in the initial grid, the cells ar distorted even in the fine +grid and if you want to do high-accuracy computations on circles, you +may want to use a better coarse grid than the one on the left. +

Hyper-L

@@ -123,7 +140,11 @@ can be used to test the efficiency of error estimators.

- + @@ -141,9 +162,9 @@ This example will create the default hyper-L. #include <grid/grid_generator.h> const unsigned int dim=2; // For example -Triangulation<dim> tr; +Triangulation<dim> triangulation; -GridGenerator<dim>::hyper_L(tr,-1,1); +GridGenerator::hyper_L (triangulation, -1, 1); @@ -151,14 +172,21 @@ GridGenerator<dim>::hyper_L(tr,-1,1);

Reading a grid from a file

In many problems taken from real life you will find those grid types -insufficient and you will need to create your own grid. Fortunately, -this is possible. deal.II offers the possibility of +insufficient and you will need to create your own grid. +This is possible: deal.II offers the possibility of reading a complete triangulation from a file in the ucd-format used by avs. A ucd-file can be read with the function
void DataIn::read_ucd(istream&)

-

Vertex numbering in input files: +

+Since the creation of input grids is probably not something you need to do on a +daily basis, we will only very briefly describe the format of the +input file. More information can be found in the documentation of the +GridIn class. +

+

+Vertex numbering in input files: The vertex numbering must start at the vertex with the lowest number for lines and be counterclockwise for quads. Also, faces between adjacent cells must face in the same direction. This must be ensured by you, or by the @@ -181,9 +209,10 @@ for creating grid cells and their properties. This is done in the order

  • boundaries
  • create a triangulation from this information
  • -You first create the vertices, then create cells by assigning vertices -to them, then you set the boundary conditions. Last you use this information -to create a triangulation. This is probably best illustrated by an example. +In fact, this is the way the functions in the +GridGenerator class use to produce the grids we have +shown above. +

    Example: Below we show the includes, definitions and @@ -229,20 +258,23 @@ for (unsigned int i=0; i<3; ++i) // to distinguish them from cell data like vertices and material. SubCellData boundary_info; -// We are using a boundary of cell number 1 +// We want to describe a line in 2d, which would be a cell in 1d, +// which is why we can use the object for a 1d cell +// (CellData<1>) for its description: boundary_info.boundary_lines.push_back (CellData<1>()); -// The boundary gets a material id of -boundary_info.boundary_lines.back().material_id = 1; +// The boundary gets a material id of 1. The boundary id of lines in 2d +// will be interpreted as the boundary id +boundary_info.boundary_lines[0].material_id = 1; -// The boundary is between vertices number 1 and 3 +// The boundary is between vertices number 0 and 3 boundary_info.boundary_lines[0].vertices[0] = 0; boundary_info.boundary_lines[0].vertices[1] = 3; // From this information the triangulation is created coarse_grid->create_triangulation (vector<Point<2> >(&vertices[0], &vertices[8]), - cells, boundary_info); + cells, boundary_info); diff --git a/deal.II/doc/tutorial/chapter-1.elements/hanging_nodes.html b/deal.II/doc/tutorial/chapter-1.elements/hanging_nodes.html index 16edd8c19d..e324d1905b 100644 --- a/deal.II/doc/tutorial/chapter-1.elements/hanging_nodes.html +++ b/deal.II/doc/tutorial/chapter-1.elements/hanging_nodes.html @@ -19,16 +19,17 @@

    Hanging Nodes

    -This chapter discusses locally refined grids, which are the only ones with -hanging nodes. A node is called hanging if its value is subject to a -constraint. -In particular, nodes on the boundary of your domain -are not hanging. +This chapter discusses one aspect of locally refined grids. Using the +approach chosen in deal.II, grids are produced with +so-called hanging nodes. +A node is called hanging if its value is subject to a +constraint induced by the requirement of globally continuous +functions. We will make this clearer in the following. -

    When do you get hanging nodes ?

    +

    When do you get hanging nodes?

    -You get hanging nodes after each refinement. Hanging nodes are those that +You get hanging nodes after some refinements. Hanging nodes are those that are situated on the boundary of a larger active cell. The global function is assumed to be linear on the boundary of each active cell, therefore the value @@ -63,6 +64,9 @@ with function values x1 and x2. If you get a hanging node on this boundary during refinement it will be situated in the middle and therefore its value - the value of the function on the boundary at this point - will be x3 = 1/2 (x1+x2). +The value x3 is therefore not a free degree of freedom, +but it is bound by this constraining equation and we will have to take +this into account.

    Mathematical implications, or
    @@ -107,7 +111,7 @@ or, taking

    It is not possible to generate ~A directly, but -w can generate A and condense it to ~A +we can generate A and condense it to ~A using the constraints, solve the system and obtain u by u=Cy.

    diff --git a/deal.II/doc/tutorial/chapter-1.elements/matrix_structure.html b/deal.II/doc/tutorial/chapter-1.elements/matrix_structure.html index 61838520ff..d2f2d86817 100644 --- a/deal.II/doc/tutorial/chapter-1.elements/matrix_structure.html +++ b/deal.II/doc/tutorial/chapter-1.elements/matrix_structure.html @@ -28,6 +28,28 @@ very large, and therefore storage is expensive in terms of memory. providing means of generating memory efficient structures for sparse matrices where only non-zero elements are stored.

    +

    +To show only two example of matrices which are quite typical of finite +element discretizations of partial differential equations, we present +the following two graphics, in which each dot represents a non-zero +entry of the matrix: +

    +

    + Sparsity 1 + Sparsity 2 +

    +

    +The first is a 300x300 matrix, the second has approximately 4000 rows +and columns. As can be seen, the vast majority of entries is zero, and +we can also guess that the proportion of zeroes is even growing as the +matrix size is increased (this is in fact so, but can not so clearly +be seen from the pictures, admittedly). The equation discretized to +yield these matrices is rather complex, so the number of non-zero +entries per row is quite large for finite element matrices. +Both matrices were derived after renumbering of the degrees of freedom +using the algorithm of Cuthill-McKee, which concentrates the entries +around the diagonal of the matrix. +

    Ways of dealing with matrix structures provided by deal.II

    @@ -44,16 +66,49 @@ tell the matrix for which elements to allocate memory. The deal.II function initializing a standard sparse matrix structure is void SparseMatrixStruct::reinit(const unsigned int m, const unsigned int n, const unsigned int max_per_row). It takes as its arguments the size of the matrix in question and the maximum number of non-zero elements -in one row. -This number can be calculated with -int DoFHandler::max_couplings_between_dofs(). -Beware: This function was written to cope with refined grids only ! -If your grid consists exclusively of cells next to boundaries its -return value will be wrong. -This matrix structure can then be used to generate a matrix using -void SparseMatrix::reinit(const SparseMatrixStruct &sparsity). -In cases of locally refined grids you will also need to take care of the constraints to your -matrix given by the hanging nodes. +in one row. This first sets up a blank sparsity pattern which, +however, has as yet no information which entries are non-zero. +

    +

    +Before going on with showing how to find out which elements are +non-zero, we have to talk about how to find out about the maximum +number of non-zero entries per row which we will need for our matrix. +Since this number depends of the actual finite element used (higher +order elements need more non-zero entries than linear elements, for +example), of the equation under consideration, and of some properties +of the triangulation on which we want to do our computations. Since +all these factors can't be considered by the user all the time, there +is a library function which does the necessary computations for you: +unsigned int DoFHandler::max_couplings_between_dofs(). +We note that this function was written to cope with refined grids, and +if you only have structured grids without local refinement and hanging +nodes, the result of this function will overestimate the maximal +number of non-zero entries per row. +

    +

    +After having fixed the basic properties of the sparsity structure, we +can tell it which entries actually are non-zero. This is again done +using a library function, +DoFHandler::make_sparsity_pattern. Then we have to take +into account that constraints due to hanging nodes (see the previous +chapter) add some more non-zero entries, which is again done using +some library functions (see the example below). After this, we can +instruct the sparsity pattern to throw out all unused elements +(in those rows where the actual number of non-zero elements is less +than the maximal number given to reinit, we can discard +the unused one to save some more memory). This is discussed in the +next subsection. +

    +

    +This matrix structure can then be used to serve as sparsity pattern to +one or several matrices using +void SparseMatrix::reinit(const SparseMatrixStruct&sparsity). +We note the distinction that the SparseMatrixStruct only +stores the places of the non-zero elements, while the +SparseMatrix actually stores the values of these +elements. The separation was made since often several matrices are +needed with the same sparsity pattern, but different entry values; +this way, we can save some memory space.

    @@ -68,25 +123,29 @@ appropriate places. This example initializes a sparse square matrix structure. const unsigned int dim=2; // For example -SparseMatrixStruct<double> sparsity_pattern; -DoFHandler<dim> dof; -ConstraintMatrix hanging_nodes; // Only necessary for locally refined grids +SparseMatrixStruct sparsity_pattern; +DoFHandler<dim> dof_handler; -// Your degrees of freedom must already be distributed +// set up a triangulation, refine it, associate the DoFHandler object +// shown above to it, and distribute the degrees of freedom +... -sparsity_pattern.reinit(dof.n_dofs(),dof.n_dofs(),dof.max_couplings_between_dofs()); +// now fix the size of the matrix and the maximal number +// of non-zero entries: +sparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + +// tell the sparsity pattern which entries are non-zero: +dof_handler.make_sparsity_pattern (sparsity_pattern); // Condense the hanging nodes into -// the matrix, taking into account the constraints to the hanging nodes. -// This must be done before the matrix structure is compressed, but only -// on locally refined grids. +// the matrix structure, taking into account the constraints due to +// the hanging nodes. // You can skip this part if your grid is not locally refined. - -dof.make_sparsity_pattern(sparsity_pattern); -hanging_nodes.clear(); -dof.make_constraint_matrix(hanging_nodes); -hanging_nodes.condense(sparsity_pattern); - +ConstraintMatrix hanging_nodes; +dof_handler.make_constraint_matrix(hanging_nodes); +hanging_nodes.condense (sparsity_pattern); @@ -127,6 +186,19 @@ sparsity_pattern.compress(); +

    More complex matrices

    + +Finally, we note that in some cases not all variables in a system of +partial differential equations couple to each other, or that the +operators have a special structure. For this reason, there are more +specialized version of the +DoFHandler::make_sparsity_pattern and +SparseMatrixStruct::reinit functions, which allow you to +more flexibly determine the actual number of non-zero elements of all +or individual rows. See the documentation of the respective classes +for more information. + +
    Hyper-L in three dimensions +Hyper-L in three dimensions (not a good picture, since the proportions +of short and long lines is not exactly one half, but you might get a +clue about it) +
    Hyper-L
    diff --git a/deal.II/doc/tutorial/chapter-1.elements/title.html b/deal.II/doc/tutorial/chapter-1.elements/title.html index 52150dac85..08e3cc5c6a 100644 --- a/deal.II/doc/tutorial/chapter-1.elements/title.html +++ b/deal.II/doc/tutorial/chapter-1.elements/title.html @@ -19,7 +19,8 @@ -

    The deal.II tutorial - Chapter 0 - Structure of a Finite Element Program

    +

    The deal.II tutorial - Chapter 0 - Basic +Elements of a Finite Element Program

    diff --git a/deal.II/doc/tutorial/chapter-1.elements/toc.html b/deal.II/doc/tutorial/chapter-1.elements/toc.html index 94b5781858..48b793463f 100644 --- a/deal.II/doc/tutorial/chapter-1.elements/toc.html +++ b/deal.II/doc/tutorial/chapter-1.elements/toc.html @@ -11,17 +11,18 @@ -

    Structure of a Finite Element Program

    +

    Basic Elements of a Finite Element Program

    About this tutorial

    This part of the deal.II documentation will discuss -how to build a finite element program. We will take a look at the general +some building blocks of a finite element program. We will take a look at the general structure of programs of this kind and at the various possibilities in that respect offered by deal.II. The order of the table of contents below reflects the order in which you need to perform the tasks. Some things may be left out according to -your needs. +your needs. The application in actual programs can be found in other +chapters of this tutorial.

    Table of contents

    diff --git a/deal.II/doc/tutorial/images/sparsity_1.jpg b/deal.II/doc/tutorial/images/sparsity_1.jpg new file mode 100644 index 0000000000..f0f1c070f1 Binary files /dev/null and b/deal.II/doc/tutorial/images/sparsity_1.jpg differ diff --git a/deal.II/doc/tutorial/images/sparsity_2.jpg b/deal.II/doc/tutorial/images/sparsity_2.jpg new file mode 100644 index 0000000000..aef82c191b Binary files /dev/null and b/deal.II/doc/tutorial/images/sparsity_2.jpg differ diff --git a/deal.II/doc/tutorial/index.html b/deal.II/doc/tutorial/index.html index 96be3330a3..8f66c1ea3e 100644 --- a/deal.II/doc/tutorial/index.html +++ b/deal.II/doc/tutorial/index.html @@ -9,24 +9,11 @@ - Deal Logo + Deal Logo

    The deal.II Tutorial

    or deal.II in several easy and some more hard lessons

    -

    Prerequisites

    -

    - This tutorial is written in HTML - 4.0 which your browser may or may not support. If you run into problems - you are advised to try a recent browser or one that is able to parse Document - Type Definitions. The reason for this is the math support of HTML 4.0. Most - people should, however, not run into any - problems.

    -

    - You yourself should have some background knowledge on numerical methods, finite - elements and C++. -

    -

    Usage of this document

    This tutorial is split into several parts, starting with an outline of @@ -43,28 +30,37 @@

    On deal.II

    deal.II is short for Differential Equations Analysis Library, version II, - which says it all, really. Its purpose it the solution of partial differential - equations on unstructured grids with error control. You can find more + which says it all, really. Its purpose is the solution of partial differential + equations on unstructured, locally refined grids with error control. You can find more information on deal.II and some examples at - - http://gaia.iwr.uni-heidelberg.de/~deal/. deal.II was written - by the numerics group at the IWR, University of Heidelberg. + the + deal.II homepage . + deal.II was written + by the numerics group at the IWR, University of Heidelberg, + mainly by Wolfgang Bangerth and in parts by Guido Kanschat, but + several others have contributed as well.

    -

    Chapter 0: Structure of a - Finite Element Program

    + +

    + The structure of this tutorial is as follows: +

    + +

    Chapter 0: Basic Elements of a + Finite Element Program

    (also available in a version without frames)

    We will discuss the various elements needed for a working - finite element program and some of the possibilitier + finite element program and some of the possibilities deal.II offers in this respect.

    -

    Chapter 1: The Laplace - Problem

    +

    Chapter 1: The very simple Laplace + Problem

    (also available in a version without frames)

    + *********** We will solve the stationary Laplace problem on a square grid with boundary conditions of the form cos(2*PI*x)cos(2*PI*y). You can also obtain the source code. @@ -87,7 +83,7 @@


    - Jan Schrage
    + The deal.II group

    Last modified: $Date$