]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add Markus Buerg's step-45.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 7 May 2010 12:58:08 +0000 (12:58 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 7 May 2010 12:58:08 +0000 (12:58 +0000)
git-svn-id: https://svn.dealii.org/trunk@21092 0785d39b-7218-0410-832d-ea1e28bc413d

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

diff --git a/deal.II/examples/step-45/Makefile b/deal.II/examples/step-45/Makefile
new file mode 100644 (file)
index 0000000..a2f3f96
--- /dev/null
@@ -0,0 +1,154 @@
+# $Id$
+
+
+# For the small projects Makefile, you basically need to fill in only
+# four fields.
+#
+# The first is the name of the application. It is assumed that the
+# application name is the same as the base file name of the single C++
+# file from which the application is generated.
+target = $(basename $(shell echo step-*.cc))
+
+# The second field determines whether you want to run your program in
+# debug or optimized mode. The latter is significantly faster, but no
+# run-time checking of parameters and internal states is performed, so
+# you should set this value to `on' while you develop your program,
+# and to `off' when running production computations.
+debug-mode = on
+
+
+# As third field, we need to give the path to the top-level deal.II
+# directory. You need to adjust this to your needs. Since this path is
+# probably the most often needed one in the Makefile internals, it is
+# designated by a single-character variable, since that can be
+# reference using $D only, i.e. without the parentheses that are
+# required for most other parameters, as e.g. in $(target).
+D = ../../
+
+
+# The last field specifies the names of data and other files that
+# shall be deleted when calling `make clean'. Object and backup files,
+# executables and the like are removed anyway. Here, we give a list of
+# files in the various output formats that deal.II supports.
+clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk
+
+
+
+
+#
+#
+# Usually, you will not need to change anything beyond this point.
+#
+#
+# The next statement tell the `make' program where to find the
+# deal.II top level directory and to include the file with the global
+# settings
+include $D/common/Make.global_options
+
+
+# Since the whole project consists of only one file, we need not
+# consider difficult dependencies. We only have to declare the
+# libraries which we want to link to the object file, and there need
+# to be two sets of libraries: one for the debug mode version of the
+# application and one for the optimized mode. Here we have selected
+# the versions for 2d. Note that the order in which the libraries are
+# given here is important and that your applications won't link
+# properly if they are given in another order.
+#
+# You may need to augment the lists of libraries when compiling your
+# program for other dimensions, or when using third party libraries
+libs.g   = $(lib-deal2-2d.g) \
+          $(lib-lac.g)      \
+           $(lib-base.g)
+libs.o   = $(lib-deal2-2d.o) \
+          $(lib-lac.o)      \
+           $(lib-base.o)
+
+
+# We now use the variable defined above which switch between debug and
+# optimized mode to select the set of libraries to link with. Included
+# in the list of libraries is the name of the object file which we
+# will produce from the single C++ file. Note that by default we use
+# the extension .g.o for object files compiled in debug mode and .o for
+# object files in optimized mode (or whatever the local default on your
+# system is instead of .o).
+ifeq ($(debug-mode),on)
+  libraries = $(target).g.$(OBJEXT) $(libs.g)
+else
+  libraries = $(target).$(OBJEXT) $(libs.o)
+endif
+
+
+# Now comes the first production rule: how to link the single object
+# file produced from the single C++ file into the executable. Since
+# this is the first rule in the Makefile, it is the one `make' selects
+# if you call it without arguments.
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS)
+
+
+# To make running the application somewhat independent of the actual
+# program name, we usually declare a rule `run' which simply runs the
+# program. You can then run it by typing `make run'. This is also
+# useful if you want to call the executable with arguments which do
+# not change frequently. You may then want to add them to the
+# following rule:
+run: $(target)
+       @echo ============================ Running $<
+       @./$(target)$(EXEEXT)
+
+
+# As a last rule to the `make' program, we define what to do when
+# cleaning up a directory. This usually involves deleting object files
+# and other automatically created files such as the executable itself,
+# backup files, and data files. Since the latter are not usually quite
+# diverse, you needed to declare them at the top of this file.
+clean:
+       -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files)
+
+
+# Since we have not yet stated how to make an object file from a C++
+# file, we should do so now. Since the many flags passed to the
+# compiler are usually not of much interest, we suppress the actual
+# command line using the `at' sign in the first column of the rules
+# and write the string indicating what we do instead.
+./%.g.$(OBJEXT) :
+       @echo ==============debug========= $(<F)
+       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+./%.$(OBJEXT) :
+       @echo ==============optimized===== $(<F)
+       @$(CXX) $(CXXFLAGS.o) -c $< -o $@
+
+
+# The following statement tells make that the rules `run' and `clean'
+# are not expected to produce files of the same name as Makefile rules
+# usually do.
+.PHONY: run clean
+
+
+# Finally there is a rule which you normally need not care much about:
+# since the executable depends on some include files from the library,
+# besides the C++ application file of course, it is necessary to
+# re-generate the executable when one of the files it depends on has
+# changed. The following rule to created a dependency file
+# `Makefile.dep', which `make' uses to determine when to regenerate
+# the executable. This file is automagically remade whenever needed,
+# i.e. whenever one of the cc-/h-files changed. Make detects whether
+# to remake this file upon inclusion at the bottom of this file.
+#
+# If the creation of Makefile.dep fails, blow it away and fail
+Makefile.dep: $(target).cc Makefile \
+              $(shell echo $D/*/include/*/*.h)
+       @echo ============================ Remaking $@
+       @$D/common/scripts/make_dependencies  $(INCLUDE) -B. $(target).cc \
+               > $@ \
+         || (rm -f $@ ; false)
+       @if test -s $@ ; then : else rm $@ ; fi
+
+
+# To make the dependencies known to `make', we finally have to include
+# them:
+include Makefile.dep
+
+
diff --git a/deal.II/examples/step-45/doc/intro.dox b/deal.II/examples/step-45/doc/intro.dox
new file mode 100644 (file)
index 0000000..a799344
--- /dev/null
@@ -0,0 +1,63 @@
+<br>
+
+<i>This program was contributed by Markus B&uuml;rg.</i>
+<br>
+
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
+In this example we consider how to use periodic boundary conditions in
+deal.II. Periodic boundary conditions often occur in computations of photonic
+crystals, because they have a lattice-like structure and, thus, it often
+suffices to do the actual computation on only one cell. To be able to proceed
+this way one has to assume that the computation can be periodically extended
+to the other cells. This requires the solution to be periodic w.r.t. the
+cells. Hence the solution has to obtain the same values on the parts of the
+boundary, which are facing each other. In the figure below we show this
+concept in two space-dimensions. There, all faces with the same color should
+have the same boundary values:
+
+@image html step-45.periodic_cells.png
+
+To keep things simple, in this tutorial we will consider an academic problem,
+which is a little bit easier: We solve the Poisson problem on a domain, where
+the left and right parts of the boundary are identified. Let $\Omega=(0,1)^2$
+and consider the problem 
+@f{align*}
+   -\Delta u &=
+   \pi^2\sin(\pi x)\sin(\pi y)  \qquad &\text{in }\Omega
+  \\ 
+   u(0,y) &= 0 \qquad &\text{for }y\in(0,1)
+  \\
+   u(1,y) &= 0 \qquad &\text{for }y\in(0,1)
+  \\
+   u(x,0) &= u(x,1) \qquad &\text{for }x\in(0,1) 
+@f}
+
+The way one has to see these periodic boundary conditions $u(x,0) = u(x,1)$ is
+as follows: Assume for a moment (as we do in this program) that we have a
+uniformly refined mesh. Then, after discretization there are a number of nodes
+(degrees of freedom) with indices $i \in {\cal I}_b$ on the bottom boundary of
+the domain, and  a second set of nodes at the top boundary $j \in {\cal
+I}_t$. Since we have assumed that the mesh is uniformly refined, there is
+exactly one node $j \in {\cal I}_t$ for each $i \in {\cal I}_b$ so that
+${\mathrm x}_j = {\mathrm x}_i + (0,1)^T$, i.e. the two of them match with
+respect to the periodicity. We will then write that $j=\text{periodic}(i)$
+(and, if you want, $i=\text{periodic}(j)$).
+If now $U_k, k=0,\ldots,N-1$ are the unknowns of our discretized problem, then
+the periodic boundary condition boils down to the following set of
+constraints:
+@f{align*}
+  U_{\text{periodic}(i)} = U_i, \qquad \forall i \in {\cal I}_b. 
+@f}
+Now, this is exactly the sort of constraint that the ConstraintMatrix class
+handles and can enforce in linear system. Consequently, the main point of this
+program is how we fill the ConstraintMatrix object that stores these
+constraints, and how this is applied to the resulting linear system.
+
+The code for solving this problem is simple and based on step-3 since we want
+to focus on the implementation of the periodic boundary conditions. The code
+could be much more sophisticated, of course. For example, we could want to
+enforce periodic boundary conditions for adaptively refined meshes in which
+there is no longer a one-to-one relationship between degrees of freedom. We
+will discuss this at the end of the results section of this program.
diff --git a/deal.II/examples/step-45/doc/results.dox b/deal.II/examples/step-45/doc/results.dox
new file mode 100644 (file)
index 0000000..63961a3
--- /dev/null
@@ -0,0 +1,20 @@
+<h1>Results</h1>
+
+The textual output the program generates is not very surprising. It just
+prints out the usual information on degrees of freedom, active cells and
+iteration steps of the solver: 
+
+@code
+============================ Running step-45
+Number of active cells: 1024
+Degrees of freedom: 1089
+DEAL:cg::Starting value 0.153970
+DEAL:cg::Convergence step 54 value 7.73356e-16
+@endcode
+
+The plot of the solution can be found in the figure below. We can see that the
+solution is constantly zero on the upper and the lower part of the boundary as
+required by the homogeneous Dirichlet boundary conditions. On the left and
+right parts the values coincide with each other, just as we wanted:
+
+@image html step-45.solution.png
diff --git a/deal.II/examples/step-45/doc/step-45.periodic_cells.png b/deal.II/examples/step-45/doc/step-45.periodic_cells.png
new file mode 100644 (file)
index 0000000..3e19509
Binary files /dev/null and b/deal.II/examples/step-45/doc/step-45.periodic_cells.png differ
diff --git a/deal.II/examples/step-45/doc/step-45.solution.png b/deal.II/examples/step-45/doc/step-45.solution.png
new file mode 100644 (file)
index 0000000..20dae47
Binary files /dev/null and b/deal.II/examples/step-45/doc/step-45.solution.png differ
diff --git a/deal.II/examples/step-45/step-45.cc b/deal.II/examples/step-45/step-45.cc
new file mode 100644 (file)
index 0000000..962fb7d
--- /dev/null
@@ -0,0 +1,275 @@
+/* $Id$ */
+/* Author: Markus Buerg, University of Karlsruhe, 2010 */
+
+/*    $Id$       */
+/*                                                                */
+/*    Copyright (C) 2010 by the deal.II authors */
+/*                                                                */
+/*    This file is subject to QPL and may not be  distributed     */
+/*    without copyright and license information. Please refer     */
+/*    to the file deal.II/doc/license.html for the  text  and     */
+/*    further information on this license.                        */
+
+                                // The include files are already known.
+#include <base/function.h>
+#include <base/quadrature_lib.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+#include <grid/grid_generator.h>
+#include <grid/tria.h>
+#include <lac/constraint_matrix.h>
+#include <lac/precondition.h>
+#include <lac/solver_cg.h>
+#include <lac/solver_control.h>
+#include <lac/sparse_matrix.h>
+#include <lac/sparsity_pattern.h>
+#include <numerics/data_out.h>
+#include <numerics/vectors.h>
+
+using namespace dealii;
+                                // The RightHandSide class is a function
+                                // object representing the right-hand side of
+                                // the problem.
+class RightHandSide: public Function<2> {
+  public:
+    RightHandSide ();
+    virtual double value (const Point<2>& p, const unsigned int component = 0) const;
+};
+
+                                // This function returns the value of the
+                                // right-hand side at a given point
+                                // <code>p</code>.
+double RightHandSide::value (const Point<2>&p, const unsigned int) const {
+  return numbers::PI * numbers::PI * std::sin (numbers::PI * p (0)) * std::sin (numbers::PI * p (1));
+}
+
+                                // Here comes the constructor of the class
+                                // <code>RightHandSide</code>.
+RightHandSide::RightHandSide (): Function<2> () {
+}
+
+                                // The class <code>LaplaceProblem</code> is
+                                // the main class, which solves the problem.
+class LaplaceProblem {
+  private:
+    ConstraintMatrix constraints;
+    const RightHandSide right_hand_side;
+    Triangulation<2> triangulation;
+    DoFHandler<2> dof_handler;
+    FE_Q<2> fe;
+    SparseMatrix<double> A;
+    SparsityPattern sparsity_pattern;
+    Vector<double> b;
+    Vector<double> u;
+    void assemble_system ();
+    void output_results ();
+    void setup_system ();
+    void solve ();
+   
+  public:
+    LaplaceProblem ();
+    void run ();
+};
+
+                                // The constructor of the class, where the
+                                // <code>DoFHandler</code> and the finite
+                                // element object are initialized.
+LaplaceProblem::LaplaceProblem (): dof_handler (triangulation), fe (1) {
+}
+
+                                //Assembling the system matrix and the
+                                //right-hand side vector is done as in other
+                                //tutorials before.
+void LaplaceProblem::assemble_system () {
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  QGauss<2> quadrature (2);
+  const unsigned int n_quadrature_points = quadrature.size ();
+  double JxW;
+  FEValues<2> fe_values (fe, quadrature, update_gradients | update_JxW_values | update_quadrature_points | update_values);
+  FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+  std::vector<Point<2> > quadrature_points;
+  std::vector<unsigned int> cell_dof_indices (dofs_per_cell);
+  Vector<double> cell_rhs (dofs_per_cell);
+   
+  for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active (); cell != dof_handler.end (); ++cell) {
+    cell_rhs = 0;
+    fe_values.reinit (cell);
+    quadrature_points = fe_values.get_quadrature_points ();
+      
+    for (unsigned int i = 0; i < dofs_per_cell; ++i)
+      for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) {
+       JxW = fe_values.JxW (q_point);
+       cell_rhs (i) += JxW * fe_values.shape_value (i, q_point) * right_hand_side.value (quadrature_points[q_point]);
+            
+       for (unsigned int j = 0; j < dofs_per_cell; ++j)
+         cell_matrix (i, j) += JxW * fe_values.shape_grad (i, q_point) * fe_values.shape_grad (j, q_point);
+      }
+      
+    cell->get_dof_indices (cell_dof_indices);
+    constraints.distribute_local_to_global (cell_matrix, cell_rhs, cell_dof_indices, A, b);
+  }
+}
+
+void LaplaceProblem::setup_system () {
+  GridGenerator::hyper_cube (triangulation);
+                                  // We change the boundary indicator on the
+                                  // parts of the boundary, where we have
+                                  // Dirichlet boundary conditions, to one
+                                  // such that we can distinguish between the
+                                  // parts of the boundary, where periodic
+                                  // and where Dirichlet boundary conditions
+                                  // hold.
+  Triangulation<2>::active_cell_iterator cell = triangulation.begin_active ();
+   
+  cell->face (2)->set_boundary_indicator (1);
+  cell->face (3)->set_boundary_indicator (1);   
+  triangulation.refine_global (5);
+                                  // Here the degrees of freedom are
+                                  // distributed.
+  dof_handler.distribute_dofs (fe);
+  std::cout << "Number of active cells: " << triangulation.n_active_cells () << std::endl << "Degrees of freedom: " << dof_handler.n_dofs () << std::endl;
+                                  // Now it is the time for the constraint
+                                  // matrix. The first constraints we put in
+                                  // are the periodic boundary
+                                  // conditions. For this let us consider the
+                                  // constraints we have to take care of in
+                                  // more detail first: We want to identify
+                                  // all degrees of freedom located on the
+                                  // right part of the boundary with the ones
+                                  // located on the left part. Thus, first we
+                                  // select a degree of freedom on the right
+                                  // part of the boundary. Then we look for
+                                  // the corresponding one the left-hand side
+                                  // and identify these two. Since we are
+                                  // using finite elements of order 1 here,
+                                  // finding the corresponding degree of
+                                  // freedom on the other side is quite easy:
+                                  // All degrees of freedom are located on
+                                  // vertices and thus we simply can take the
+                                  // degree of freedom on the other side,
+                                  // which is located on the vertex with the
+                                  // same y-component.  Here starts the
+                                  // implementation: First we declare a
+                                  // vector, which stores the global index of
+                                  // the boundary degree of freedom together
+                                  // with the y-component of the vertex on
+                                  // which it is located.
+  constraints.clear ();
+   
+  std::vector<std::pair<unsigned int, double> > dof_locations;
+   
+  dof_locations.reserve (dof_handler.n_boundary_dofs ());
+   
+  unsigned int dofs_per_face = fe.dofs_per_face;
+                                  // Then we loop over all active cells and
+                                  // check, whether the cell is located at
+                                  // the boundary. If this is the case, we
+                                  // check, if it is located on the right
+                                  // part of the boundary, hence, if face 1
+                                  // is at the boundary.
+  for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active (); cell != dof_handler.end (); ++cell)
+    if (cell->at_boundary () && cell->face (1)->at_boundary ()) {
+                                      // Since each degree of freedom of the
+                                      // face is located on one vertex, we
+                                      // can access them directly and store
+                                      // the global index of the degree of
+                                      // freedom togehter with the
+                                      // y-component of the vertex on which
+                                      // it is located.
+      dof_locations.push_back (std::pair<unsigned int, double> (cell->vertex_dof_index (1, 0), cell->vertex (1) (1)));
+      dof_locations.push_back (std::pair<unsigned int, double> (cell->vertex_dof_index (3, 0), cell->vertex (3) (1)));
+    }
+                                  // Now we have to find the corresponding
+                                  // degrees of freedom on the left part of
+                                  // the boundary. Therefore we loop over all
+                                  // cells again and choose the ones, where
+                                  // face 0 is at the boundary.
+  for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active (); cell != dof_handler.end (); ++cell)
+    if (cell->at_boundary () && cell->face (0)->at_boundary ()) {
+                                      // For every degree of freedom on this
+                                      // face we add a new line to the
+                                      // constraint matrix. Then we identify
+                                      // it with the corresponding degree of
+                                      // freedom on the right part of the
+                                      // boundary.
+      constraints.add_line (cell->vertex_dof_index (0, 0));
+         
+      for (unsigned int i = 0; i < dof_locations.size (); ++i)
+       if (dof_locations[i].second == cell->vertex (0) (1)) {
+         constraints.add_entry (cell->vertex_dof_index (0, 0), dof_locations[i].first, 1.0);
+         break;
+       }
+         
+      constraints.add_line (cell->vertex_dof_index (2, 0));
+         
+      for (unsigned int i = 0; i < dof_locations.size (); ++i)
+       if (dof_locations[i].second == cell->vertex (2) (1)) {
+         constraints.add_entry (cell->vertex_dof_index (2, 0), dof_locations[i].first, 1.0);
+         break;
+       }
+    }
+                                  // Finally we have to set the homogeneous
+                                  // Dirichlet boundary conditions on the
+                                  // upper and lower parts of the boundary
+                                  // and close the
+                                  // <code>ConstraintMatrix</code> object.
+  VectorTools::interpolate_boundary_values (dof_handler, 1, ZeroFunction<2> (), constraints);
+  constraints.close ();
+                                  // Then we create the sparsity pattern and
+                                  // the system matrix and initialize the
+                                  // solution and right-hand side vectors.
+  const unsigned int n_dofs = dof_handler.n_dofs ();
+   
+  sparsity_pattern.reinit (n_dofs, n_dofs, dof_handler.max_couplings_between_dofs ());
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern, constraints, false);
+  sparsity_pattern.compress ();
+  A.reinit (sparsity_pattern);
+  b.reinit (n_dofs);
+  u.reinit (n_dofs);
+}
+                                // To solve the linear system of equations
+                                // $Au=b$ we use the CG solver with an
+                                // SSOR-preconditioner.
+void LaplaceProblem::solve () {
+  SolverControl solver_control (dof_handler.n_dofs (), 1e-15);
+  PreconditionSSOR<SparseMatrix<double> > precondition;
+   
+  precondition.initialize (A);
+   
+  SolverCG<> cg (solver_control);
+   
+  cg.solve (A, u, b, precondition);
+  constraints.distribute (u);
+}
+
+void LaplaceProblem::output_results () {
+                                  // As graphical output we create vtk-file
+                                  // of the computed solution.
+  DataOut<2> data_out;
+   
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (u, "u");
+  data_out.build_patches ();
+   
+  std::ofstream output ("solution.vtk");
+   
+  data_out.write_vtk (output);
+}
+                                // This function manages the solving process
+                                // of the problem.
+void LaplaceProblem::run () {
+  setup_system ();
+  assemble_system ();
+  solve ();
+  output_results ();
+}
+                                // And at the end we have the main function
+                                // as usual.
+int main () {
+  LaplaceProblem laplace_problem;
+   
+  laplace_problem.run ();
+}

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.