--- /dev/null
+# $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
+
+
--- /dev/null
+<br>
+
+<i>This program was contributed by Markus Bü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.
--- /dev/null
+/* $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 ();
+}
projects / dealii-svn.git / commitdiff
