From: bangerth Date: Fri, 7 May 2010 12:58:08 +0000 (+0000) Subject: Add Markus Buerg's step-45. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=15971114992e852037774dee9f02642cd057aef2;p=dealii-svn.git Add Markus Buerg's step-45. git-svn-id: https://svn.dealii.org/trunk@21092 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-45/Makefile b/deal.II/examples/step-45/Makefile new file mode 100644 index 0000000000..a2f3f96ca5 --- /dev/null +++ b/deal.II/examples/step-45/Makefile @@ -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========= $( $@ \ + || (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 index 0000000000..a799344358 --- /dev/null +++ b/deal.II/examples/step-45/doc/intro.dox @@ -0,0 +1,63 @@ +
+ +This program was contributed by Markus Bürg. +
+ + +

Introduction

+ +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 index 0000000000..63961a384d --- /dev/null +++ b/deal.II/examples/step-45/doc/results.dox @@ -0,0 +1,20 @@ +

Results

+ +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 index 0000000000..3e19509c87 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 index 0000000000..20dae47f30 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 index 0000000000..962fb7d18c --- /dev/null +++ b/deal.II/examples/step-45/step-45.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +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 + // p. +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 + // RightHandSide. +RightHandSide::RightHandSide (): Function<2> () { +} + + // The class LaplaceProblem 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 A; + SparsityPattern sparsity_pattern; + Vector b; + Vector u; + void assemble_system (); + void output_results (); + void setup_system (); + void solve (); + + public: + LaplaceProblem (); + void run (); +}; + + // The constructor of the class, where the + // DoFHandler 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 cell_matrix (dofs_per_cell, dofs_per_cell); + std::vector > quadrature_points; + std::vector cell_dof_indices (dofs_per_cell); + Vector 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 > 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 (cell->vertex_dof_index (1, 0), cell->vertex (1) (1))); + dof_locations.push_back (std::pair (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 + // ConstraintMatrix 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 > 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 (); +}