From: Wolfgang Bangerth Date: Mon, 22 Jun 2009 14:07:52 +0000 (+0000) Subject: Add the skeleton of an eigenvalue/function solver using SLEPc. X-Git-Tag: v8.0.0~7603 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9d23073c009102c08357e9aa678d41b40cf8e78d;p=dealii.git Add the skeleton of an eigenvalue/function solver using SLEPc. git-svn-id: https://svn.dealii.org/trunk@18947 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-36/Makefile b/deal.II/examples/step-36/Makefile new file mode 100644 index 0000000000..b753dea123 --- /dev/null +++ b/deal.II/examples/step-36/Makefile @@ -0,0 +1,174 @@ +# $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 + + + + +# +# +# 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 + +################################################################ +# This example program will only work if PETSc is installed. If this +# is not the case, then simply redefine the main targets to to nothing +ifneq ($(USE_CONTRIB_PETSC),yes) +default run clean: + @echo + @echo "===========================================================" + @echo "= This program cannot be compiled without PETSc. Make =" + @echo "= sure you have PETSc installed and detected during =" + @echo "= configuration of deal.II =" + @echo "===========================================================" + @echo + +else +# +################################################################ + + + + +# 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 + + +endif # CONTRIB_USE_PETSC diff --git a/deal.II/examples/step-36/doc/intro.dox b/deal.II/examples/step-36/doc/intro.dox new file mode 100644 index 0000000000..928f558739 --- /dev/null +++ b/deal.II/examples/step-36/doc/intro.dox @@ -0,0 +1,10 @@ +
+ +This program was contributed by Toby Young and Wolfgang +Bangerth. + + + + +

Introduction

+ diff --git a/deal.II/examples/step-36/doc/results.dox b/deal.II/examples/step-36/doc/results.dox new file mode 100644 index 0000000000..b5eaba9377 --- /dev/null +++ b/deal.II/examples/step-36/doc/results.dox @@ -0,0 +1,2 @@ +

Results

+ diff --git a/deal.II/examples/step-36/step-36.cc b/deal.II/examples/step-36/step-36.cc new file mode 100644 index 0000000000..f9f6d047bb --- /dev/null +++ b/deal.II/examples/step-36/step-36.cc @@ -0,0 +1,254 @@ +/* $Id$ */ +/* Author: Toby Young, Polish Academy of Sciences, */ +/* Wolfgang Bangerth, Texas A&M University */ +/* $Id$ */ +/* */ +/* Copyright (C) 2009 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. */ + + // @sect3{Include files} + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + + // The final step, as in previous + // programs, is to import all the + // deal.II class and function names + // into the global namespace: +using namespace dealii; + + // @sect3{The EigenvalueProblem class template} + +template +class EigenvalueProblem +{ + public: + EigenvalueProblem (); + void run (); + + private: + void make_grid_and_dofs (); + void assemble_system (); + void solve (); + void output_results () const; + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + + PETScWrappers::SparseMatrix stiffness_matrix, mass_matrix; + std::vector eigenfunctions; + + ParameterHandler parameters; +}; + + + + // @sect3{Implementation of the EigenvalueProblem class} + + // @sect4{EigenvalueProblem::EigenvalueProblem} + +template +EigenvalueProblem::EigenvalueProblem () + : + fe (1), + dof_handler (triangulation) +{ + parameters.declare_entry ("Number of eigenvalues/eigenfunctions", "10", + Patterns::Integer (0, 100), + "The number of eigenvalues/eigenfunctions " + "to be computed."); + parameters.declare_entry ("Potential", "0", + Patterns::Anything(), + "A functional description of the potential."); +} + + + // @sect4{EigenvalueProblem::make_grid_and_dofs} +template +void EigenvalueProblem::make_grid_and_dofs () +{ + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (4); + dof_handler.distribute_dofs (fe); + + CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(), + dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, csp); + stiffness_matrix.reinit (csp); + mass_matrix.reinit (csp); + + eigenfunctions + .resize (parameters.get_integer ("Number of eigenvalues/eigenfunctions")); + for (unsigned int i=0; i +void EigenvalueProblem::assemble_system () +{ + QGauss quadrature_formula(2); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_stiffness_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix cell_mass_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + FunctionParser potential; + potential.initialize ((dim == 2 ? + "x,y" : + "x,y,z"), + parameters.get ("Potential"), + typename FunctionParser::ConstMap()); + + ConstraintMatrix constraints; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(), + constraints); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + cell_stiffness_matrix = 0; + cell_mass_matrix = 0; + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + + constraints.distribute_local_to_global (cell_stiffness_matrix, + local_dof_indices, + stiffness_matrix); + constraints.distribute_local_to_global (cell_mass_matrix, + local_dof_indices, + mass_matrix); + } +} + + + // @sect4{EigenvalueProblem::solve} +template +void EigenvalueProblem::solve () +{ +// do whatever is necessary here +} + + + // @sect4{EigenvalueProblem::output_results} + +template +void EigenvalueProblem::output_results () const +{ +// do whatever is necessary here +} + + + + // @sect4{EigenvalueProblem::run} + + // This is the function which has the + // top-level control over everything. It is + // the same as for the previous example. +template +void EigenvalueProblem::run () +{ + std::cout << "Solving problem in " << dim << " space dimensions." << std::endl; + + make_grid_and_dofs(); + assemble_system (); + solve (); + output_results (); +} + + + // @sect3{The main function} + +int main () +{ + try + { + deallog.depth_console (0); + EigenvalueProblem<2> laplace_problem_2d; + laplace_problem_2d.run (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +}