From 05621823f1b863c0f8993a7337ff027d4f0c1e3c Mon Sep 17 00:00:00 2001 From: wolf Date: Fri, 14 Sep 2001 09:28:50 +0000 Subject: [PATCH] Add yet another example program. git-svn-id: https://svn.dealii.org/trunk@5001 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-11/.cvsignore | 3 + deal.II/examples/step-11/Makefile | 159 +++++++++++ deal.II/examples/step-11/step-11.cc | 400 ++++++++++++++++++++++++++++ 3 files changed, 562 insertions(+) create mode 100644 deal.II/examples/step-11/.cvsignore create mode 100644 deal.II/examples/step-11/Makefile create mode 100644 deal.II/examples/step-11/step-11.cc diff --git a/deal.II/examples/step-11/.cvsignore b/deal.II/examples/step-11/.cvsignore new file mode 100644 index 0000000000..1f752bee9c --- /dev/null +++ b/deal.II/examples/step-11/.cvsignore @@ -0,0 +1,3 @@ +*.o *.go Makefile.dep *.gnuplot *.gmv *.eps +step-11 +*.dat diff --git a/deal.II/examples/step-11/Makefile b/deal.II/examples/step-11/Makefile new file mode 100644 index 0000000000..c02ef0f148 --- /dev/null +++ b/deal.II/examples/step-11/Makefile @@ -0,0 +1,159 @@ +# $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 something 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 correct compiler flags and 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 .go for object files +# compiled in debug mode and .o for object files in optimized mode. +ifeq ($(debug-mode),on) + libraries = $(target).go $(libs.g) + flags = $(CXXFLAGS.g) +else + libraries = $(target).o $(libs.o) + flags = $(CXXFLAGS.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 $@ $^ $(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) + + +# 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 *.o *.go *~ Makefile.dep $(target) $(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. +%.go : %.cc + @echo ==============debug========= $( Makefile.dep + +# To make the dependencies known to `make', we finally have to include +# them: +include Makefile.dep + + diff --git a/deal.II/examples/step-11/step-11.cc b/deal.II/examples/step-11/step-11.cc new file mode 100644 index 0000000000..a4f38dd082 --- /dev/null +++ b/deal.II/examples/step-11/step-11.cc @@ -0,0 +1,400 @@ +/* $Id$ */ +/* Author: Wolfgang Bangerth, University of Heidelberg, 2001 */ + + // As usual, the program starts with + // a rather long list of include + // files which you are probably + // already used to by now: +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + // Just this one is new: it declares + // a class + // ``CompressedSparsityPattern'', + // which we will use and explain + // further down below. +#include + + // We will make use of the std::find + // algorithm of the C++ standard + // library, so we have to include the + // following file for its + // declaration: +#include + + +template +double measure (const DoFHandler &dof_handler, + const Mapping &mapping) +{ + QGauss4 quadrature; + FEValues fe_values (mapping, dof_handler.get_fe(), quadrature, + update_JxW_values); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + double measure = 0; + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + for (unsigned int i=0; i +double measure (const Triangulation &triangulation, + const Mapping &mapping) +{ + FE_Q dummy_fe(1); + DoFHandler dof_handler (const_cast&>(triangulation)); + dof_handler.distribute_dofs(dummy_fe); + return measure (dof_handler, mapping); +}; + + +template +double surface (const DoFHandler &dof_handler, + const Mapping &mapping) +{ + QGauss4 quadrature; + FEFaceValues fe_values (mapping, dof_handler.get_fe(), quadrature, + update_JxW_values); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + double surface = 0; + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->at_boundary()) + { + fe_values.reinit (cell, face); + for (unsigned int i=0; i +double surface (const Triangulation &triangulation, + const Mapping &mapping) +{ + FE_Q dummy_fe(1); + DoFHandler dof_handler (const_cast&>(triangulation)); + dof_handler.distribute_dofs(dummy_fe); + return surface (dof_handler, mapping); +}; + + +template double surface (const Triangulation<2> &, const Mapping<2> &); +template double measure (const Triangulation<2> &, const Mapping<2> &); + + + + + // Then we declare a class which + // represents the solution of a + // Laplace problem. As this example + // program is based on step-5, the + // class looks rather the same, with + // the sole structural difference + // that we have merged the functions + // ``assemble_system'' and ``solve'', + // and the output function was + // dropped since the solution + // function is so boring that it is + // not worth being viewed. + // + // The only other noteworthy change + // is that the constructor takes a + // value representing the polynomial + // degree of the mapping to be used + // later on, and that it has another + // member variable representing + // exactly this mapping. In general, + // this variable will occur in real + // applications at the same places + // where the finite element is + // declared or used. +template +class LaplaceProblem +{ + public: + LaplaceProblem (const unsigned int mapping_degree); + void run (); + + private: + void setup_system (); + void assemble_and_solve (); + void solve (); + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + MappingQ mapping; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + ConstraintMatrix mean_value_constraints; + + Vector solution; + Vector system_rhs; +}; + + + + // Construct such an object, by + // initializing the variables. Here, + // we use linear finite elements (the + // argument to the ``fe'' variable + // denotes the polynomial degree), + // and mappings of given order. Print + // to screen what we are about to do. +template +LaplaceProblem::LaplaceProblem (const unsigned int mapping_degree) : + fe (1), + dof_handler (triangulation), + mapping (mapping_degree) +{ + std::cout << "Using mapping with degree " << mapping_degree << ":" + << std::endl + << "============================" + << std::endl; +}; + + + + // The first task is to set up the + // variables for this problem. This + // includes generating a valid + // ``DoFHandler'' object, as well as + // the sparsity patterns for the + // matrix, and the object + // representing the constraints that + // the mean value of the degrees of + // freedom on the boundary be zero. +template +void LaplaceProblem::setup_system () +{ + // The first task is trivial: + // generate an enumeration of the + // degrees of freedom: + dof_handler.distribute_dofs (fe); + + // Next task is to construct the + // object representing the + // constraint that the mean value + // of the degrees of freedom on the + // boundary shall be zero. For + // this, we first want a list of + // those nodes which are actually + // at the boundary. The + // ``DoFTools'' class has a + // function that returns an array + // of boolean values where ``true'' + // indicates that the node is at + // the boundary. The second + // argument denotes a mask + // selecting which components of + // vector valued finite elements we + // want to be considered. Since we + // have a scalar finite element + // anyway, this mask consists of + // only one entry, and its value + // must be ``true''. + std::vector boundary_dofs (dof_handler.n_dofs(), false); + DoFTools::extract_boundary_dofs (dof_handler, std::vector(1,true), + boundary_dofs); + + // Let us first pick out the first + // boundary node from this list. We + // do that by searching for the + // first ``true'' value in the + // array (note that ``std::find'' + // returns an iterator to this + // element), and computing its + // distance to the overall first + // element in the array to get its + // index: + const unsigned int first_boundary_dof + = std::distance (std::find (boundary_dofs.begin(), + boundary_dofs.end(), + true), + boundary_dofs.begin()); + + mean_value_constraints.clear (); + mean_value_constraints.add_line (first_boundary_dof); + for (unsigned int i=first_boundary_dof+1; i +void LaplaceProblem::assemble_and_solve () +{ + QGauss2 cell_quadrature; + QGauss2 face_quadrature; + MatrixTools::create_laplace_matrix (mapping, dof_handler, + cell_quadrature, + system_matrix); + VectorTools::create_right_hand_side (mapping, dof_handler, + cell_quadrature, + ConstantFunction(-2), + system_rhs); + + Vector tmp (system_rhs.size()); + VectorTools::create_boundary_right_hand_side (mapping, dof_handler, + face_quadrature, + ConstantFunction(1), + tmp); + system_rhs += tmp; + + mean_value_constraints.condense (system_matrix); + mean_value_constraints.condense (system_rhs); + + solve (); + mean_value_constraints.distribute (solution); + + Vector difference_per_cell (triangulation.n_active_cells()); + VectorTools::integrate_difference (mapping, dof_handler, + solution, + ZeroFunction(), + difference_per_cell, + QGauss3(), + H1_seminorm); + std::cout << " " << triangulation.n_active_cells() << " cells: " + << " |u|_1=" + << difference_per_cell.l2_norm() + << ", error=" + << fabs(difference_per_cell.l2_norm()-sqrt(3.14159265358/2)) + << std::endl; +}; + + + +template +void LaplaceProblem::solve () +{ + SolverControl solver_control (1000, 1e-12); + PrimitiveVectorMemory<> vector_memory; + SolverCG<> cg (solver_control, vector_memory); + + PreconditionSSOR<> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); +}; + + + +template +void LaplaceProblem::run () +{ + GridGenerator::hyper_ball (triangulation); + static const HyperBallBoundary boundary; + triangulation.set_boundary (0, boundary); + + for (unsigned int cycle=0; cycle<6; ++cycle, triangulation.refine_global(1)) + { + setup_system (); + assemble_and_solve (); + }; +}; + + + + // Finally the main function. It's + // structure is the same as that used + // in several of the previous + // examples, so probably needs no + // more explanation. +int main () +{ + try + { + deallog.depth_console (0); + std::cout.precision(5); + + // This is the main loop, doing + // the computations with + // mappings of linear through + // cubic mappings. Note that + // since we need the object of + // type ``LaplaceProblem<2>'' + // only once, we do not even + // name it, but create an + // unnamed such object and call + // the ``run'' function of it, + // subsequent to which it is + // immediately destroyed again. + for (unsigned int mapping_degree=1; mapping_degree<=3; ++mapping_degree) + LaplaceProblem<2>(mapping_degree).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; +}; -- 2.39.5