From: wolf Date: Wed, 29 Dec 1999 22:05:31 +0000 (+0000) Subject: Add step-5. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=50ee3edd9ed251b37be3c4b74f77d78b71a7b677;p=dealii-svn.git Add step-5. git-svn-id: https://svn.dealii.org/trunk@2126 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-4/step-4.cc b/deal.II/deal.II/Attic/examples/step-by-step/step-4/step-4.cc index cf073c5509..5c95b88071 100644 --- a/deal.II/deal.II/Attic/examples/step-by-step/step-4/step-4.cc +++ b/deal.II/deal.II/Attic/examples/step-by-step/step-4/step-4.cc @@ -252,7 +252,7 @@ void LaplaceProblem::make_grid_and_dofs () cout << " Number of active cells: " << triangulation.n_active_cells() << endl - << " Total number of cells: " + << " Total number of cells: " << triangulation.n_cells() << endl; diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-5/Makefile b/deal.II/deal.II/Attic/examples/step-by-step/step-5/Makefile new file mode 100644 index 0000000000..dc2798f839 --- /dev/null +++ b/deal.II/deal.II/Attic/examples/step-by-step/step-5/Makefile @@ -0,0 +1,94 @@ +# $Id$ +# Copyright W. Bangerth, University of Heidelberg, 1999 + +# Template for makefiles for the examples subdirectory. In principle, +# everything should be done automatically if you set the target file +# here correctly. We get deduce it from the files in the present +# directory: +target = $(basename $(shell echo step-*.cc)) + +# All dependencies between files should be updated by the included +# file Makefile.dep if necessary. Object files are compiled into +# the archives ./Obj.a and ./Obj.g.a. By default, the debug version +# is used to link. It you don't like that, change the following +# variable to "off" +debug-mode = on + + + +############################################################################### +# Internals + +#deal include base path +D = ../../../.. + +include $D/common/Make.global_options + +# get lists of files we need + + +# list of libraries needed to link with +libs = -ldeal_II_2d -llac -lbase +libs.g = -ldeal_II_2d.g -llac.g -lbase.g + + +# check whether we use debug mode or not +ifeq ($(debug-mode),on) + libraries = $(target).go $(libs.g) + flags = $(CXXFLAGS.g) +else + libraries = $(target).go $(libs) + flags = $(CXXFLAGS.o) +endif + + + +# make rule for the target. $^ is the object file $(target).g?o +$(target) : $(libraries) + @echo ============================ Linking $@ + @$(CXX) $(flags) -o $@ $^ + +# rule how to run the program +run: $(target) + @echo ============================ Running $@ + @./$(target) + + +# rule to make object files +%.go : %.cc + @echo ============================ Compiling with debugging information: $< + @$(CXX) $(CXXFLAGS.g) -c $< -o $@ +%.o : %.cc + @echo ============================ Compiling with optimization: $< + @$(CXX) $(CXXFLAGS) -c $< -o $@ + + +clean: + -rm -f *.o *.go *~ Makefile.dep $(target) *gmv *gnuplot *gpl *eps + + + +.PHONY: clean + + +# Rule to generate the dependency file. 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. +# +# Since the script prefixes the output names by lib/g?o, we have to +# strip that again (the script was written for the main libraries and +# large projects where object files are put into subdirs) +Makefile.dep: $(target).cc Makefile \ + $(shell echo $D/base/include/base/*.h \ + $D/lac/include/lac/*.h \ + $D/deal.II/include/*/*.h) + @echo ============================ Remaking Makefile + @perl $D/common/scripts/make_dependencies.pl $(INCLUDE) $(target).cc \ + | perl -pi -e 's!lib/g?o/!!g;' \ + > Makefile.dep + + +include Makefile.dep + diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-5/step-5.cc b/deal.II/deal.II/Attic/examples/step-by-step/step-5/step-5.cc new file mode 100644 index 0000000000..16e2c5c52e --- /dev/null +++ b/deal.II/deal.II/Attic/examples/step-by-step/step-5/step-5.cc @@ -0,0 +1,504 @@ +/* $Id$ */ + + // The first few (many?) include + // files have already been used in + // the previous example, so we will + // not explain their meaning here + // again. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + + +template +class LaplaceProblem +{ + public: + LaplaceProblem (); + void run (); + + private: + void make_grid_and_dofs (const unsigned int refinement); + void assemble_system (); + void solve (); + void output_results (); + void clear (); + + Triangulation triangulation; + FEQ1 fe; + DoFHandler dof_handler; + + SparseMatrixStruct sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + +template +class Coefficient : public Function +{ + public: + virtual double value (const Point &p, + const unsigned int component = 0) const; + virtual void value_list (const vector > &points, + vector &values, + const unsigned int component = 0) const; +}; + + + +template +double Coefficient::value (const Point &p, + const unsigned int) const +{ + if (p.square() < 0.5*0.5) + return 10; + else + return 1; +}; + + + +template +void Coefficient::value_list (const vector > &points, + vector &values, + const unsigned int component) const +{ + const unsigned int n_points = points.size(); + + Assert (values.size() == n_points, + ExcVectorHasWrongSize (values.size(), n_points)); + + Assert (component == 0, + ExcWrongComponent (component, 1)); + + for (unsigned int i=0; i +LaplaceProblem::LaplaceProblem () : + dof_handler (triangulation) +{}; + + + +template +void LaplaceProblem::make_grid_and_dofs (const unsigned int refinement) +{ + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (refinement); + + cout << " Number of active cells: " + << triangulation.n_active_cells() + << endl + << " Total number of cells: " + << triangulation.n_cells() + << endl; + + dof_handler.distribute_dofs (fe); + + cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << endl; + + sparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); + sparsity_pattern.compress(); + + system_matrix.reinit (sparsity_pattern); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); +}; + + + + // As in the previous examples, this + // function is not changed much with + // regard to its functionality, but + // there are still some optimizations + // which we will show. For this, it + // is important to note that if + // efficient solvers are used (such + // as the preconditions CG method), + // assembling the matrix and right + // hand side can take a comparable + // time, and it is worth the effort + // to use one or two optimizations at + // some places. + // + // What we will show here is how we + // can avoid calls to the + // shape_value, shape_grad, and + // quadrature_point functions of the + // FEValues object, and in particular + // optimize away most of the virtual + // function calls of the Function + // object. The way to do so will be + // explained in the following, while + // those parts of this function that + // are not changed with respect to + // the previous example are not + // commented on. +template +void LaplaceProblem::assemble_system () +{ + // This time, we will again use a + // constant right hand side + // function, but a variable + // coefficient. The following + // object will be used for this: + const Coefficient coefficient; + + QGauss3 quadrature_formula; + + FEValues fe_values (fe, quadrature_formula, + UpdateFlags(update_values | + update_gradients | + update_q_points | + update_JxW_values)); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + vector local_dof_indices (dofs_per_cell); + + // ... + vector coefficient_values (n_q_points); + + DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_matrix.clear (); + cell_rhs.clear (); + + // As before, we want the + // FEValues object to compute + // the quantities which we told + // him to compute in the + // constructor using the update + // flags. + fe_values.reinit (cell); + // Now, these quantities are + // stored in arrays in the + // FEValues object. Usually, + // the internals of how and + // where they are stored is not + // something that the outside + // world should know, but since + // this is a time critical + // function we decided to + // publicize these arrays a + // little bit, and provide + // facilities to export the + // address where this data is + // stored. + // + // For example, the values of + // shape function j at + // quadrature point q is stored + // in a matrix, of which we can + // get the address as follows + // (note that this is a + // reference to the matrix, + // symbolized by the ampersand, + // and that it must be a + // constant reference, since + // only read-only access is + // granted): + const FullMatrix + & shape_values = fe_values.get_shape_values(); + // Instead of writing + // fe_values.shape_value(j,q) + // we can now write + // shape_values(j,q), i.e. the + // function call needed + // previously for each access + // has been otimized away. + // + // There are alike functions + // for almost all data elements + // in the FEValues class. The + // gradient are accessed as + // follows: + const vector > > + & shape_grads = fe_values.get_shape_grads(); + // The data type looks a bit + // unwieldy, since each entry + // in the matrix (j,q) now + // needs to be the gradient of + // the shape function, which is + // a vector. + // + // Similarly, access to the + // place where quadrature + // points and the determinants + // of the Jacobian matrices + // times the weights of the + // respective quadrature points + // are stored, can be obtained + // like this: + const vector + & JxW_values = fe_values.get_JxW_values(); + const vector > + & q_points = fe_values.get_quadrature_points(); + // Admittedly, the declarations + // above are not easily + // readable, but they can save + // many function calls in the + // inner loops and can thus + // make assemblage faster. + // + // An additional advantage is + // that the inner loops are + // simpler to read, since the + // fe_values object is no more + // explicitely needed to access + // the different fields (see + // below). Unfortunately, + // things became a bit + // inconsistent, since the + // shape values are accessed + // via the FullMatrix operator + // (), i.e. using parentheses, + // while all the other fields + // are accessed through vector + // operator [], i.e. using + // brackets. This is due to + // historical reasons and + // frequently leads to a bit of + // confusion, but since the + // places where this happens + // are few in well-written + // programs, this is not too + // big a problem. + + // There is one more thing: in + // this example, we want to use + // a non-constant + // coefficient. In the previous + // example, we have called the + // ``value'' function of the + // right hand side object for + // each quadrature + // point. Unfortunately, that + // is a virtual function, so + // calling it is relatively + // expensive. Therefore, we use + // a function of the Function + // class which returns the + // values at all quadrature + // points at once; that + // function is still virtual, + // but it needs to be computed + // once per cell only, not once + // in the inner loop: + coefficient.value_list (q_points, coefficient_values); + // It should be noted that the + // creation of the + // coefficient_values object is + // done outside the loop over + // all cells to avoid memory + // allocation each time we + // visit a new cell. Contrary + // to this, the other variables + // above were created inside + // the loop, but they were only + // references to memory that + // has already been allocated + // (i.e. they are pointers to + // that memory) and therefore, + // no new memory needs to be + // allocated; in particular, by + // declaring the pointers as + // close to their use as + // possible, we give the + // compiler a better choice to + // optimize them away + // altogether, something which + // it definitely can't do with + // the coefficient_values + // object since it is too + // complicated, but mostly + // because it's address is + // passed to a virtual function + // which is not knows at + // compile time. + + // Using the various + // abbreviations, the loops + // then look like this (the + // parentheses around the + // product of the two gradients + // are needed to indicate the + // dot product; we have to + // overrule associativity of + // the operator* here, since + // the compiler would otherwise + // complain about an undefined + // product of double*gradient + // since it parses + // left-to-right): + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + for (unsigned int i=0; i boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); +}; + + + +template +void LaplaceProblem::solve () +{ + SolverControl solver_control (1000, 1e-12); + PrimitiveVectorMemory<> vector_memory; + SolverCG<> cg (solver_control, vector_memory); + + // ... + PreconditionRelaxation<> + preconditioner(system_matrix, + &SparseMatrix::template precondition_SSOR, + 1.2); + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + + cout << " " << solver_control.last_step() + << " CG iterations needed to obtain convergence." + << endl; +}; + + + +template +void LaplaceProblem::output_results () +{ + DataOut data_out; + + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "solution"); + + data_out.build_patches (); + + ofstream output (dim == 2 ? + "solution-2d.gmv" : + "solution-3d.gmv"); + // ... + data_out.write_gnuplot (output); +}; + + + +template +void LaplaceProblem::clear () +{ + system_rhs.reinit (0); + solution.reinit (0); + system_matrix.reinit (); + sparsity_pattern.reinit (0, 0, 0); + dof_handler.clear (); + triangulation.clear (); +}; + + + +template +void LaplaceProblem::run () +{ + cout << "Solving problem in " << dim << " space dimensions." << endl; + + for (unsigned int refinement=0; refinement<7; ++refinement) + { + cout << "Refinement step: " << refinement << endl; + + make_grid_and_dofs(refinement); + assemble_system (); + solve (); + output_results (); + + clear (); + }; +}; + + + +int main () +{ + deallog.depth_console (0); + + LaplaceProblem<2> laplace_problem_2d; + laplace_problem_2d.run (); + + return 0; +}; diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/navbar.html b/deal.II/doc/tutorial/chapter-2.step-by-step/navbar.html index fbb14abf65..0124224506 100644 --- a/deal.II/doc/tutorial/chapter-2.step-by-step/navbar.html +++ b/deal.II/doc/tutorial/chapter-2.step-by-step/navbar.html @@ -48,6 +48,12 @@ Step 4

+ +
  • +

    + Step 5 +

    +
  • Back to the tutorial index diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.results b/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.results index 3201c10650..7514fac593 100644 --- a/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.results +++ b/deal.II/doc/tutorial/chapter-2.step-by-step/step-4.results @@ -1,7 +1,10 @@

    Results

    -The output of the program looks as follows: +The output of the program looks as follows (the number of iterations +may vary by one or two, depending on your computer, since this is +often dependent on the round-off accuracy of floating point +operations, which differs between micro-processors):
     
     Solving problem in 2 space dimensions.
    diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-5.intro b/deal.II/doc/tutorial/chapter-2.step-by-step/step-5.intro
    new file mode 100644
    index 0000000000..1518606e43
    --- /dev/null
    +++ b/deal.II/doc/tutorial/chapter-2.step-by-step/step-5.intro
    @@ -0,0 +1,3 @@
    +
    +

    Introduction

    + diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-5.results b/deal.II/doc/tutorial/chapter-2.step-by-step/step-5.results new file mode 100644 index 0000000000..747bdbbb3d --- /dev/null +++ b/deal.II/doc/tutorial/chapter-2.step-by-step/step-5.results @@ -0,0 +1,3 @@ + +

    Results

    + diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/toc.html b/deal.II/doc/tutorial/chapter-2.step-by-step/toc.html index 32b036a221..0bec85471b 100644 --- a/deal.II/doc/tutorial/chapter-2.step-by-step/toc.html +++ b/deal.II/doc/tutorial/chapter-2.step-by-step/toc.html @@ -65,6 +65,13 @@ At present, the following programs exist: 3D, although the program is exactly the same. Non-constant right hand side function. Non-homogeneous boundary values. + +
    Step 5
    +
    What's new: Computations on successively + refined grids. Some optimizations. Non-constant coefficient in + the elliptic operator. Preconditioning of the CG solver for the + linear system of equations. +
    diff --git a/deal.II/examples/step-4/step-4.cc b/deal.II/examples/step-4/step-4.cc index cf073c5509..5c95b88071 100644 --- a/deal.II/examples/step-4/step-4.cc +++ b/deal.II/examples/step-4/step-4.cc @@ -252,7 +252,7 @@ void LaplaceProblem::make_grid_and_dofs () cout << " Number of active cells: " << triangulation.n_active_cells() << endl - << " Total number of cells: " + << " Total number of cells: " << triangulation.n_cells() << endl; diff --git a/deal.II/examples/step-5/Makefile b/deal.II/examples/step-5/Makefile new file mode 100644 index 0000000000..dc2798f839 --- /dev/null +++ b/deal.II/examples/step-5/Makefile @@ -0,0 +1,94 @@ +# $Id$ +# Copyright W. Bangerth, University of Heidelberg, 1999 + +# Template for makefiles for the examples subdirectory. In principle, +# everything should be done automatically if you set the target file +# here correctly. We get deduce it from the files in the present +# directory: +target = $(basename $(shell echo step-*.cc)) + +# All dependencies between files should be updated by the included +# file Makefile.dep if necessary. Object files are compiled into +# the archives ./Obj.a and ./Obj.g.a. By default, the debug version +# is used to link. It you don't like that, change the following +# variable to "off" +debug-mode = on + + + +############################################################################### +# Internals + +#deal include base path +D = ../../../.. + +include $D/common/Make.global_options + +# get lists of files we need + + +# list of libraries needed to link with +libs = -ldeal_II_2d -llac -lbase +libs.g = -ldeal_II_2d.g -llac.g -lbase.g + + +# check whether we use debug mode or not +ifeq ($(debug-mode),on) + libraries = $(target).go $(libs.g) + flags = $(CXXFLAGS.g) +else + libraries = $(target).go $(libs) + flags = $(CXXFLAGS.o) +endif + + + +# make rule for the target. $^ is the object file $(target).g?o +$(target) : $(libraries) + @echo ============================ Linking $@ + @$(CXX) $(flags) -o $@ $^ + +# rule how to run the program +run: $(target) + @echo ============================ Running $@ + @./$(target) + + +# rule to make object files +%.go : %.cc + @echo ============================ Compiling with debugging information: $< + @$(CXX) $(CXXFLAGS.g) -c $< -o $@ +%.o : %.cc + @echo ============================ Compiling with optimization: $< + @$(CXX) $(CXXFLAGS) -c $< -o $@ + + +clean: + -rm -f *.o *.go *~ Makefile.dep $(target) *gmv *gnuplot *gpl *eps + + + +.PHONY: clean + + +# Rule to generate the dependency file. 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. +# +# Since the script prefixes the output names by lib/g?o, we have to +# strip that again (the script was written for the main libraries and +# large projects where object files are put into subdirs) +Makefile.dep: $(target).cc Makefile \ + $(shell echo $D/base/include/base/*.h \ + $D/lac/include/lac/*.h \ + $D/deal.II/include/*/*.h) + @echo ============================ Remaking Makefile + @perl $D/common/scripts/make_dependencies.pl $(INCLUDE) $(target).cc \ + | perl -pi -e 's!lib/g?o/!!g;' \ + > Makefile.dep + + +include Makefile.dep + diff --git a/deal.II/examples/step-5/step-5.cc b/deal.II/examples/step-5/step-5.cc new file mode 100644 index 0000000000..16e2c5c52e --- /dev/null +++ b/deal.II/examples/step-5/step-5.cc @@ -0,0 +1,504 @@ +/* $Id$ */ + + // The first few (many?) include + // files have already been used in + // the previous example, so we will + // not explain their meaning here + // again. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + + +template +class LaplaceProblem +{ + public: + LaplaceProblem (); + void run (); + + private: + void make_grid_and_dofs (const unsigned int refinement); + void assemble_system (); + void solve (); + void output_results (); + void clear (); + + Triangulation triangulation; + FEQ1 fe; + DoFHandler dof_handler; + + SparseMatrixStruct sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + +template +class Coefficient : public Function +{ + public: + virtual double value (const Point &p, + const unsigned int component = 0) const; + virtual void value_list (const vector > &points, + vector &values, + const unsigned int component = 0) const; +}; + + + +template +double Coefficient::value (const Point &p, + const unsigned int) const +{ + if (p.square() < 0.5*0.5) + return 10; + else + return 1; +}; + + + +template +void Coefficient::value_list (const vector > &points, + vector &values, + const unsigned int component) const +{ + const unsigned int n_points = points.size(); + + Assert (values.size() == n_points, + ExcVectorHasWrongSize (values.size(), n_points)); + + Assert (component == 0, + ExcWrongComponent (component, 1)); + + for (unsigned int i=0; i +LaplaceProblem::LaplaceProblem () : + dof_handler (triangulation) +{}; + + + +template +void LaplaceProblem::make_grid_and_dofs (const unsigned int refinement) +{ + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (refinement); + + cout << " Number of active cells: " + << triangulation.n_active_cells() + << endl + << " Total number of cells: " + << triangulation.n_cells() + << endl; + + dof_handler.distribute_dofs (fe); + + cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << endl; + + sparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); + sparsity_pattern.compress(); + + system_matrix.reinit (sparsity_pattern); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); +}; + + + + // As in the previous examples, this + // function is not changed much with + // regard to its functionality, but + // there are still some optimizations + // which we will show. For this, it + // is important to note that if + // efficient solvers are used (such + // as the preconditions CG method), + // assembling the matrix and right + // hand side can take a comparable + // time, and it is worth the effort + // to use one or two optimizations at + // some places. + // + // What we will show here is how we + // can avoid calls to the + // shape_value, shape_grad, and + // quadrature_point functions of the + // FEValues object, and in particular + // optimize away most of the virtual + // function calls of the Function + // object. The way to do so will be + // explained in the following, while + // those parts of this function that + // are not changed with respect to + // the previous example are not + // commented on. +template +void LaplaceProblem::assemble_system () +{ + // This time, we will again use a + // constant right hand side + // function, but a variable + // coefficient. The following + // object will be used for this: + const Coefficient coefficient; + + QGauss3 quadrature_formula; + + FEValues fe_values (fe, quadrature_formula, + UpdateFlags(update_values | + update_gradients | + update_q_points | + update_JxW_values)); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + vector local_dof_indices (dofs_per_cell); + + // ... + vector coefficient_values (n_q_points); + + DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_matrix.clear (); + cell_rhs.clear (); + + // As before, we want the + // FEValues object to compute + // the quantities which we told + // him to compute in the + // constructor using the update + // flags. + fe_values.reinit (cell); + // Now, these quantities are + // stored in arrays in the + // FEValues object. Usually, + // the internals of how and + // where they are stored is not + // something that the outside + // world should know, but since + // this is a time critical + // function we decided to + // publicize these arrays a + // little bit, and provide + // facilities to export the + // address where this data is + // stored. + // + // For example, the values of + // shape function j at + // quadrature point q is stored + // in a matrix, of which we can + // get the address as follows + // (note that this is a + // reference to the matrix, + // symbolized by the ampersand, + // and that it must be a + // constant reference, since + // only read-only access is + // granted): + const FullMatrix + & shape_values = fe_values.get_shape_values(); + // Instead of writing + // fe_values.shape_value(j,q) + // we can now write + // shape_values(j,q), i.e. the + // function call needed + // previously for each access + // has been otimized away. + // + // There are alike functions + // for almost all data elements + // in the FEValues class. The + // gradient are accessed as + // follows: + const vector > > + & shape_grads = fe_values.get_shape_grads(); + // The data type looks a bit + // unwieldy, since each entry + // in the matrix (j,q) now + // needs to be the gradient of + // the shape function, which is + // a vector. + // + // Similarly, access to the + // place where quadrature + // points and the determinants + // of the Jacobian matrices + // times the weights of the + // respective quadrature points + // are stored, can be obtained + // like this: + const vector + & JxW_values = fe_values.get_JxW_values(); + const vector > + & q_points = fe_values.get_quadrature_points(); + // Admittedly, the declarations + // above are not easily + // readable, but they can save + // many function calls in the + // inner loops and can thus + // make assemblage faster. + // + // An additional advantage is + // that the inner loops are + // simpler to read, since the + // fe_values object is no more + // explicitely needed to access + // the different fields (see + // below). Unfortunately, + // things became a bit + // inconsistent, since the + // shape values are accessed + // via the FullMatrix operator + // (), i.e. using parentheses, + // while all the other fields + // are accessed through vector + // operator [], i.e. using + // brackets. This is due to + // historical reasons and + // frequently leads to a bit of + // confusion, but since the + // places where this happens + // are few in well-written + // programs, this is not too + // big a problem. + + // There is one more thing: in + // this example, we want to use + // a non-constant + // coefficient. In the previous + // example, we have called the + // ``value'' function of the + // right hand side object for + // each quadrature + // point. Unfortunately, that + // is a virtual function, so + // calling it is relatively + // expensive. Therefore, we use + // a function of the Function + // class which returns the + // values at all quadrature + // points at once; that + // function is still virtual, + // but it needs to be computed + // once per cell only, not once + // in the inner loop: + coefficient.value_list (q_points, coefficient_values); + // It should be noted that the + // creation of the + // coefficient_values object is + // done outside the loop over + // all cells to avoid memory + // allocation each time we + // visit a new cell. Contrary + // to this, the other variables + // above were created inside + // the loop, but they were only + // references to memory that + // has already been allocated + // (i.e. they are pointers to + // that memory) and therefore, + // no new memory needs to be + // allocated; in particular, by + // declaring the pointers as + // close to their use as + // possible, we give the + // compiler a better choice to + // optimize them away + // altogether, something which + // it definitely can't do with + // the coefficient_values + // object since it is too + // complicated, but mostly + // because it's address is + // passed to a virtual function + // which is not knows at + // compile time. + + // Using the various + // abbreviations, the loops + // then look like this (the + // parentheses around the + // product of the two gradients + // are needed to indicate the + // dot product; we have to + // overrule associativity of + // the operator* here, since + // the compiler would otherwise + // complain about an undefined + // product of double*gradient + // since it parses + // left-to-right): + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + for (unsigned int i=0; i boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); +}; + + + +template +void LaplaceProblem::solve () +{ + SolverControl solver_control (1000, 1e-12); + PrimitiveVectorMemory<> vector_memory; + SolverCG<> cg (solver_control, vector_memory); + + // ... + PreconditionRelaxation<> + preconditioner(system_matrix, + &SparseMatrix::template precondition_SSOR, + 1.2); + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + + cout << " " << solver_control.last_step() + << " CG iterations needed to obtain convergence." + << endl; +}; + + + +template +void LaplaceProblem::output_results () +{ + DataOut data_out; + + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "solution"); + + data_out.build_patches (); + + ofstream output (dim == 2 ? + "solution-2d.gmv" : + "solution-3d.gmv"); + // ... + data_out.write_gnuplot (output); +}; + + + +template +void LaplaceProblem::clear () +{ + system_rhs.reinit (0); + solution.reinit (0); + system_matrix.reinit (); + sparsity_pattern.reinit (0, 0, 0); + dof_handler.clear (); + triangulation.clear (); +}; + + + +template +void LaplaceProblem::run () +{ + cout << "Solving problem in " << dim << " space dimensions." << endl; + + for (unsigned int refinement=0; refinement<7; ++refinement) + { + cout << "Refinement step: " << refinement << endl; + + make_grid_and_dofs(refinement); + assemble_system (); + solve (); + output_results (); + + clear (); + }; +}; + + + +int main () +{ + deallog.depth_console (0); + + LaplaceProblem<2> laplace_problem_2d; + laplace_problem_2d.run (); + + return 0; +};