From 4a9d2cadd10fbb757829ff4c9b3cb30625856f8b Mon Sep 17 00:00:00 2001 From: wolf Date: Sun, 27 Feb 2000 22:12:08 +0000 Subject: [PATCH] Almost finish program, some text. git-svn-id: https://svn.dealii.org/trunk@2497 0785d39b-7218-0410-832d-ea1e28bc413d --- .../examples/step-by-step/step-7/step-7.cc | 706 +++++++++++++++--- deal.II/examples/step-7/step-7.cc | 706 +++++++++++++++--- 2 files changed, 1242 insertions(+), 170 deletions(-) diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-7/step-7.cc b/deal.II/deal.II/Attic/examples/step-by-step/step-7/step-7.cc index 0caa12938f..46b0754f77 100644 --- a/deal.II/deal.II/Attic/examples/step-by-step/step-7/step-7.cc +++ b/deal.II/deal.II/Attic/examples/step-by-step/step-7/step-7.cc @@ -1,6 +1,10 @@ /* $Id$ */ /* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */ + // These first include files have all + // been treated in previous examples, + // so we won't explain what is in + // them again. #include #include #include @@ -11,63 +15,81 @@ #include #include #include -#include #include #include #include #include +#include +#include #include #include -#include -#include -#include #include -#include +#include #include +#include + // In this example, we will not use + // the numeration scheme which is + // used per default by the + // ``DoFHandler'' class, but will + // renumber them using the + // Cuthill-McKee algorithm. The + // necessary functions are declared + // in the following file: #include + // Then we will show a little trick + // how we can make sure that objects + // are not deleted while they are + // still in use. For this purpose, + // there is the ``SmartPointer'' + // helper class, which is declared in + // this file: #include + // Then we will want to use the + // ``integrate_difference'' function + // mentioned in the introduction. It + // comes from this file: +#include + // And finally, we need to use the + // ``FEFaceValues'' class, which is + // declare in the same file as the + // ``FEValues'' class: +#include #include -template -class LaplaceProblem -{ - public: - enum RefinementMode { - global_refinement, adaptive_refinement - }; - - LaplaceProblem (const FiniteElement &fe, - const RefinementMode refinement_mode); - ~LaplaceProblem (); - void run (); - - private: - void setup_system (); - void assemble_system (); - void solve (); - void refine_grid (); - void process_solution (const unsigned int cycle) const; - - Triangulation triangulation; - DoFHandler dof_handler; - //... - SmartPointer > fe; - ConstraintMatrix hanging_node_constraints; - - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; - - Vector solution; - Vector system_rhs; - - RefinementMode refinement_mode; -}; - - + // Since we want to compare the + // exactly known continuous solution + // to the computed one, we need a + // function object which represents + // the continuous solution. On the + // other hand, we need the right hand + // side function, and that one of + // course shares some characteristics + // with the solution. In order to + // reduce dependencies which arise if + // we have to change something in + // both classes at the same time, we + // exclude the common characteristics + // of both functions into a base + // class. + // + // The common characteristics for the + // given solution, which as explained + // in the introduction is a sum of + // three exponentials, are here: the + // number of exponentials, their + // centers, and their half width. We + // declare them in the following + // class. Since the number of + // exponentials is a constant scalar + // integral quantity, C++ allows its + // definition (i.e. assigning a + // value) right at the place of + // declaration (i.e. where we declare + // that such a variable exists). template class SolutionBase { @@ -78,30 +100,24 @@ class SolutionBase }; -template -class Solution : public Function, - protected SolutionBase -{ - public: - virtual double value (const Point &p, - const unsigned int component = 0) const; - virtual Tensor<1,dim> gradient (const Point &p, - const unsigned int component = 0) const; -}; - - - -template -class RightHandSide : public Function, - protected SolutionBase -{ - public: - virtual double value (const Point &p, - const unsigned int component = 0) const; -}; - - - + // The variables which denote the + // centers and the width of the + // exponentials have just been + // declared, now we still need to + // assign values to them. Here, we + // can show another small piece of + // template sourcery, namely how we + // can assign different values to + // these variables depending on the + // dimension. We will only use the 2d + // case in the program, but we show + // the 1d case for exposition of a + // useful technique. + // + // First we assign values to the + // centers for the 1d case, where we + // place the centers equidistanly at + // -1/3, 0, and 1/3: template <> const Point<1> SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers] @@ -109,6 +125,8 @@ SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers] Point<1>(0.0), Point<1>(+1.0 / 3.0) }; + // Then we place the centers for the + // 2d case as follows: template <> const Point<2> SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers] @@ -116,11 +134,56 @@ SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers] Point<2>(-0.5, -0.5), Point<2>(+0.5, -0.5) }; + // There remains to assign a value to + // the half-width of the + // exponentials. We would like to use + // the same value for all dimensions, + // so here is how that works: template -const double SolutionBase::width = 0.15; - +const double SolutionBase::width = 1./3.; + + + + // After declaring and defining the + // characteristics of solution and + // right hand side, we can declare + // the classes representing these + // two. They both represent + // continuous functions, so they are + // derived from the ``Function'' + // base class, and they also inherit + // the characteristics defined in the + // ``SolutionBase'' class. + // + // The actual classes are declared in + // the following. Note that in order + // to compute the error of the + // numerical solution against the + // continuous one in the L2 and H1 + // norms, we have to export value and + // gradient of the exact solution, + // which is done by overloading the + // respective virtual member + // functions in the ``Function'' base + // class. +template +class Solution : public Function, + protected SolutionBase +{ + public: + virtual double value (const Point &p, + const unsigned int component = 0) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; +}; + // The actual definition of the + // values and gradients of the exact + // solution class is according to + // their mathematical definition and + // probably needs not much + // explanation. template double Solution::value (const Point &p, const unsigned int) const @@ -128,8 +191,18 @@ double Solution::value (const Point &p, double return_value = 0; for (unsigned int i=0; i shifted_point = p-source_centers[i]; + // The ``Point'' class + // offers a member function + // ``square'' that does what + // it's name suggests. return_value += exp(-shifted_point.square() / (width*width)); }; @@ -142,11 +215,40 @@ template Tensor<1,dim> Solution::gradient (const Point &p, const unsigned int) const { + // In order to accumulate the + // gradient from the contributions + // of the exponentials, we allocate + // an object which denotes the + // mathematical quantity of a + // tensor of rank ``1'' and + // dimension ``dim''. Its default + // constructor sets it to the + // vector containing only zeroes, + // so we need not explicitely care + // for its initialization. Tensor<1,dim> return_value; + // Note that we could as well have + // taken the type of the object to + // be ``Point''. Tensors of + // rank 1 and points are almost + // exchangeable, and have only very + // slightly different mathematical + // meanings. In fact, the + // ``Point'' class is derived + // from the ``Tensor<1,dim>'' + // class, which makes up for their + // mutual exchangeability. + for (unsigned int i=0; i shifted_point = p-source_centers[i]; + // For the gradient, note that + // it's direction is along + // (x-x_i), so we add up + // multiples of this distance + // vector, where the factor is + // given by the exponentials. return_value += (-2 / (width*width) * exp(-shifted_point.square() / (width*width)) * shifted_point); @@ -157,6 +259,33 @@ Tensor<1,dim> Solution::gradient (const Point &p, + // Besides the function that + // represents the exact solution, we + // also need a function which we can + // use as right hand side when + // assembling the linear system of + // discretized equations. This is + // accomplished using the following + // class and the following definition + // of its function. Note that here we + // only need the value of the + // function, not its gradients or + // higher derivatives. +template +class RightHandSide : public Function, + protected SolutionBase +{ + public: + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + // The value of the right hand side + // is given by the negative Laplacian + // of the solution plus the solution + // itself, since we wanted to solve + // Helmholtz's equation: template double RightHandSide::value (const Point &p, const unsigned int) const @@ -166,8 +295,14 @@ double RightHandSide::value (const Point &p, { const Point shifted_point = p-source_centers[i]; - return_value += ((2*dim - 4*shifted_point.square()/(width*width)) / (width*width) * + // The first contribution is + // the Laplacian: + return_value += ((2*dim - 4*shifted_point.square()/(width*width)) / + (width*width) * exp(-shifted_point.square() / (width*width))); + // And the second is the + // solution itself: + return_value += exp(-shifted_point.square() / (width*width)); }; return return_value; @@ -175,6 +310,47 @@ double RightHandSide::value (const Point &p, + // Then we need the class that does + // all the work. +//....................... +template +class LaplaceProblem +{ + public: + enum RefinementMode { + global_refinement, adaptive_refinement + }; + + LaplaceProblem (const FiniteElement &fe, + const RefinementMode refinement_mode); + ~LaplaceProblem (); + void run (); + + private: + void setup_system (); + void assemble_system (); + void solve (); + void refine_grid (); + void process_solution (const unsigned int cycle) const; + + Triangulation triangulation; + DoFHandler dof_handler; + //... + SmartPointer > fe; + ConstraintMatrix hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + + RefinementMode refinement_mode; +}; + + + + template LaplaceProblem::LaplaceProblem (const FiniteElement &fe, const RefinementMode refinement_mode) : @@ -192,13 +368,42 @@ LaplaceProblem::~LaplaceProblem () }; - + // The following function sets up the + // degrees of freedom, sizes of + // matrices and vectors, etc. Most of + // its functionality has been showed + // in previous examples, the only + // difference being the renumbering + // step. template void LaplaceProblem::setup_system () { dof_handler.distribute_dofs (*fe); - // Renumber the degrees of freedom... + // Renumbering the degrees of + // freedom is not overly difficult, + // as long as you use one of the + // algorithms included in the + // library. It requires just one + // line of code, namely the + // following: DoFRenumbering::Cuthill_McKee (dof_handler); + // Note, however, that when you + // renumber the degrees of freedom, + // you must do so immediately after + // distributing them, since such + // things as hanging nodes, the + // sparsity pattern etc. depend on + // the absolute numbers which are + // altered by renumbering. + // + // Renumbering does not serve any + // specific purpose in this + // example, it is done only for + // exposition of the technique. To + // see the effect of renumbering on + // the sparsity pattern of the + // matrix, refer to the second + // example program. hanging_node_constraints.clear (); DoFTools::make_hanging_node_constraints (dof_handler, @@ -220,19 +425,101 @@ void LaplaceProblem::setup_system () + // Assembling the system of equations + // for the problem at hand is mostly + // as for the example programs + // before. However, some things have + // changed anyway, so we comment on + // this function fairly extensively. template void LaplaceProblem::assemble_system () { - QGauss3 quadrature_formula; - FEValues fe_values (*fe, quadrature_formula, - UpdateFlags(update_values | - update_gradients | - update_q_points | - update_JxW_values)); - + // First we need to define objects + // which will be used as quadrature + // formula for domain and face + // integrals. + // + // Note the way in which we define + // a quadrature rule for the faces: + // it is simply a quadrature rule + // for one dimension less! + QGauss3 quadrature_formula; + QGauss3 face_quadrature_formula; + // For simpler use later on, we + // alias the number of quadrature + // points to local variables: + const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + const unsigned int n_face_q_points = face_quadrature_formula.n_quadrature_points; + + // Then we need objects which can + // evaluate the values, gradients, + // etc of the shape functions at + // the quadrature points. While it + // seems that it should be feasible + // to do it with one object for + // both domain and face integrals, + // there is a subtle difference + // since the weights in the domain + // integrals include the measure of + // the cell in the domain, while + // the face integral quadrature + // requires the measure of the face + // in a lower-dimensional + // mannifold. Internally these two + // classes are rooted on a common + // base class which does most of + // the work; that, however, is + // something that you need not + // worry about. + // + // For the domain integrals in the + // bilinear form for Helmholtz's + // equation, we need to compute the + // values and gradients, as well as + // the weights at the quadrature + // points. Furthermore, we need the + // quadrature points on the real + // cell (rather than on the unit + // cell) to evaluate the right hand + // side function. + FEValues fe_values (*fe, quadrature_formula, + UpdateFlags(update_values | + update_gradients | + update_q_points | + update_JxW_values)); + + // For the face integrals, we only + // need the values of the shape + // functions, as well as the + // weights. We also need the normal + // vectors and quadrature points on + // the real cell since we want to + // determine the Neumann values + // from the exact solution object + // (see below). + FEFaceValues fe_face_values (*fe, face_quadrature_formula, + UpdateFlags(update_values | + update_q_points | + update_normal_vectors | + update_JxW_values)); + + // In order to make programming + // more readable below, we alias + // the number of degrees of freedom + // per cell to a local variable, as + // already done for the number of + // quadrature points above: const unsigned int dofs_per_cell = fe->dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + // Then we need some objects + // already known from previous + // examples: An object denoting the + // right hand side function, its + // values at the quadrature points + // on a cell, the cell matrix and + // right hand side, and the indices + // of the degrees of freedom on a + // cell. RightHandSide right_hand_side; vector rhs_values (n_q_points); @@ -241,6 +528,28 @@ void LaplaceProblem::assemble_system () vector local_dof_indices (dofs_per_cell); + // Then we define an object + // denoting the exact solution + // function. We will use it to + // compute the Neumann values at + // the boundary from it. Usually, + // one would of course do so using + // a separate object, in particular + // since the exact solution is not + // known while the Neumann values + // are prescribed. We will, + // however, be a little bit lazy + // and use what we already have in + // information. Real-life programs + // would to go other ways here, of + // course. + Solution exact_solution; + + // Now for the main loop over all + // cells. This is mostly unchanged + // from previous examples, so we + // only comment on the things that + // have changed. DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) @@ -264,16 +573,134 @@ void LaplaceProblem::assemble_system () for (unsigned int i=0; i::faces_per_cell; ++face) + if (cell->face(face)->boundary_indicator() == 1) + { + // If we came into here, + // then we have found an + // external face + // belonging to + // Gamma2. Next, we have + // to compute the values + // of the shape functions + // and the other + // quantities which we + // will need for the + // computation of the + // contour integral. This + // is done using the + // ``reinit'' function + // which we already know + // from the ``FEValue'' + // class: + fe_face_values.reinit (cell, face); + + // Then, for simpler + // access, we alias the + // various quantities to + // local variables: + const FullMatrix + & face_shape_values = fe_face_values.get_shape_values(); + const vector + & face_JxW_values = fe_face_values.get_JxW_values(); + const vector > + & face_q_points = fe_face_values.get_quadrature_points(); + const vector > + & face_normal_vectors = fe_face_values.get_normal_vectors (); + + // And we can then + // perform the + // integration by using a + // loop over all + // quadrature points. + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); for (unsigned int i=0; i::assemble_system () }; }; + // The rest of the function has + // also been shown previously: hanging_node_constraints.condense (system_matrix); hanging_node_constraints.condense (system_rhs); + // Only with the interpolation of + // boundary values, there is one + // notable thing, namely that now + // the boundary indicator for which + // we interpolate boundary values + // (denoted by the second parameter + // to + // ``interpolate_boundary_values'') + // does not represent the whole + // boundary an more. Rather, it is + // that portion of the boundary + // which we have not assigned + // another indicator (see + // below). The degrees of freedom + // at the boundary that do not + // belong to Gamma1 are therefore + // excluded from the interpolation + // of boundary values. map boundary_values; VectorTools::interpolate_boundary_values (dof_handler, 0, @@ -301,7 +748,8 @@ void LaplaceProblem::assemble_system () }; - + // Solving the system of equations is + // done in the same way as before. template void LaplaceProblem::solve () { @@ -355,8 +803,6 @@ void LaplaceProblem::refine_grid () }; -#include - template void LaplaceProblem::process_solution (const unsigned int cycle) const { @@ -402,15 +848,105 @@ void LaplaceProblem::process_solution (const unsigned int cycle) const + // The following function is the main + // one which controls the flow of + // execution. The basic layout is as + // in previous examples: an outer + // loop over successively refined + // grids, and in this loop first + // problem setup, assemblage of the + // linear system, solution, and + // postprocessing. template void LaplaceProblem::run () { - for (unsigned int cycle=0; cycle<12; ++cycle) + for (unsigned int cycle=0; cycle<9; ++cycle) { + // The first action in each + // iteration of the outer loop + // is setting up the grid on + // which we will solve in this + // iteration. In the first + // iteration, the coarsest grid + // is generated, in later + // iterations it is refined, + // for which we call the + // ``refine_grid'' function. if (cycle == 0) { + // Setting up the coarse + // grid is done as in + // previous examples: we + // first create an initial + // grid, which is the unit + // square [-1,1]x[-1,1] in + // the present case. Then + // we refine it globally a + // specific number of + // times. GridGenerator::hyper_cube (triangulation, -1, 1); triangulation.refine_global (1); + + // However, here we have to + // do something else in + // addition: mark those + // faces that belong to the + // different components of + // the boundary, Gamma1 and + // Gamma2. We will use the + // following convention: + // Faces belonging to + // Gamma1 will have the + // boundary indicator ``0'' + // (which is the default, + // so we don't have to set + // it explicitely), and + // faces belonging to + // Gamma2 will use ``1'' as + // boundary indicator. + // + // To set these values, we + // loop over all cells, + // then over all faces of a + // given cell, check + // whether it belongs to + // the boundary Gamma2, and + // if so set its boundary + // indicator to ``1''. + // + // It is worth noting that + // we have to loop over all + // cells here, not only the + // active ones. The reason + // is that upon refinement, + // newly created faces + // inherit the boundary + // indicator of their + // parent face. If we now + // only set the boundary + // indicator for active + // faces, coarsen some + // cells and refine them + // later on, they will + // again have the boundary + // indicator of the parent + // cell which we have not + // modified, instead of the + // one we + // intended. Therefore, we + // have to change the + // boundary indicators of + // all faces on Gamma2, + // irrespective whether + // they are active or not. + Triangulation::cell_iterator cell = triangulation.begin (), + endc = triangulation.end(); + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if ((cell->face(face)->center()(0) == -1) + || + (cell->face(face)->center()(1) == -1)) + cell->face(face)->set_boundary_indicator (1); } else refine_grid (); diff --git a/deal.II/examples/step-7/step-7.cc b/deal.II/examples/step-7/step-7.cc index 0caa12938f..46b0754f77 100644 --- a/deal.II/examples/step-7/step-7.cc +++ b/deal.II/examples/step-7/step-7.cc @@ -1,6 +1,10 @@ /* $Id$ */ /* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */ + // These first include files have all + // been treated in previous examples, + // so we won't explain what is in + // them again. #include #include #include @@ -11,63 +15,81 @@ #include #include #include -#include #include #include #include #include +#include +#include #include #include -#include -#include -#include #include -#include +#include #include +#include + // In this example, we will not use + // the numeration scheme which is + // used per default by the + // ``DoFHandler'' class, but will + // renumber them using the + // Cuthill-McKee algorithm. The + // necessary functions are declared + // in the following file: #include + // Then we will show a little trick + // how we can make sure that objects + // are not deleted while they are + // still in use. For this purpose, + // there is the ``SmartPointer'' + // helper class, which is declared in + // this file: #include + // Then we will want to use the + // ``integrate_difference'' function + // mentioned in the introduction. It + // comes from this file: +#include + // And finally, we need to use the + // ``FEFaceValues'' class, which is + // declare in the same file as the + // ``FEValues'' class: +#include #include -template -class LaplaceProblem -{ - public: - enum RefinementMode { - global_refinement, adaptive_refinement - }; - - LaplaceProblem (const FiniteElement &fe, - const RefinementMode refinement_mode); - ~LaplaceProblem (); - void run (); - - private: - void setup_system (); - void assemble_system (); - void solve (); - void refine_grid (); - void process_solution (const unsigned int cycle) const; - - Triangulation triangulation; - DoFHandler dof_handler; - //... - SmartPointer > fe; - ConstraintMatrix hanging_node_constraints; - - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; - - Vector solution; - Vector system_rhs; - - RefinementMode refinement_mode; -}; - - + // Since we want to compare the + // exactly known continuous solution + // to the computed one, we need a + // function object which represents + // the continuous solution. On the + // other hand, we need the right hand + // side function, and that one of + // course shares some characteristics + // with the solution. In order to + // reduce dependencies which arise if + // we have to change something in + // both classes at the same time, we + // exclude the common characteristics + // of both functions into a base + // class. + // + // The common characteristics for the + // given solution, which as explained + // in the introduction is a sum of + // three exponentials, are here: the + // number of exponentials, their + // centers, and their half width. We + // declare them in the following + // class. Since the number of + // exponentials is a constant scalar + // integral quantity, C++ allows its + // definition (i.e. assigning a + // value) right at the place of + // declaration (i.e. where we declare + // that such a variable exists). template class SolutionBase { @@ -78,30 +100,24 @@ class SolutionBase }; -template -class Solution : public Function, - protected SolutionBase -{ - public: - virtual double value (const Point &p, - const unsigned int component = 0) const; - virtual Tensor<1,dim> gradient (const Point &p, - const unsigned int component = 0) const; -}; - - - -template -class RightHandSide : public Function, - protected SolutionBase -{ - public: - virtual double value (const Point &p, - const unsigned int component = 0) const; -}; - - - + // The variables which denote the + // centers and the width of the + // exponentials have just been + // declared, now we still need to + // assign values to them. Here, we + // can show another small piece of + // template sourcery, namely how we + // can assign different values to + // these variables depending on the + // dimension. We will only use the 2d + // case in the program, but we show + // the 1d case for exposition of a + // useful technique. + // + // First we assign values to the + // centers for the 1d case, where we + // place the centers equidistanly at + // -1/3, 0, and 1/3: template <> const Point<1> SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers] @@ -109,6 +125,8 @@ SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers] Point<1>(0.0), Point<1>(+1.0 / 3.0) }; + // Then we place the centers for the + // 2d case as follows: template <> const Point<2> SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers] @@ -116,11 +134,56 @@ SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers] Point<2>(-0.5, -0.5), Point<2>(+0.5, -0.5) }; + // There remains to assign a value to + // the half-width of the + // exponentials. We would like to use + // the same value for all dimensions, + // so here is how that works: template -const double SolutionBase::width = 0.15; - +const double SolutionBase::width = 1./3.; + + + + // After declaring and defining the + // characteristics of solution and + // right hand side, we can declare + // the classes representing these + // two. They both represent + // continuous functions, so they are + // derived from the ``Function'' + // base class, and they also inherit + // the characteristics defined in the + // ``SolutionBase'' class. + // + // The actual classes are declared in + // the following. Note that in order + // to compute the error of the + // numerical solution against the + // continuous one in the L2 and H1 + // norms, we have to export value and + // gradient of the exact solution, + // which is done by overloading the + // respective virtual member + // functions in the ``Function'' base + // class. +template +class Solution : public Function, + protected SolutionBase +{ + public: + virtual double value (const Point &p, + const unsigned int component = 0) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; +}; + // The actual definition of the + // values and gradients of the exact + // solution class is according to + // their mathematical definition and + // probably needs not much + // explanation. template double Solution::value (const Point &p, const unsigned int) const @@ -128,8 +191,18 @@ double Solution::value (const Point &p, double return_value = 0; for (unsigned int i=0; i shifted_point = p-source_centers[i]; + // The ``Point'' class + // offers a member function + // ``square'' that does what + // it's name suggests. return_value += exp(-shifted_point.square() / (width*width)); }; @@ -142,11 +215,40 @@ template Tensor<1,dim> Solution::gradient (const Point &p, const unsigned int) const { + // In order to accumulate the + // gradient from the contributions + // of the exponentials, we allocate + // an object which denotes the + // mathematical quantity of a + // tensor of rank ``1'' and + // dimension ``dim''. Its default + // constructor sets it to the + // vector containing only zeroes, + // so we need not explicitely care + // for its initialization. Tensor<1,dim> return_value; + // Note that we could as well have + // taken the type of the object to + // be ``Point''. Tensors of + // rank 1 and points are almost + // exchangeable, and have only very + // slightly different mathematical + // meanings. In fact, the + // ``Point'' class is derived + // from the ``Tensor<1,dim>'' + // class, which makes up for their + // mutual exchangeability. + for (unsigned int i=0; i shifted_point = p-source_centers[i]; + // For the gradient, note that + // it's direction is along + // (x-x_i), so we add up + // multiples of this distance + // vector, where the factor is + // given by the exponentials. return_value += (-2 / (width*width) * exp(-shifted_point.square() / (width*width)) * shifted_point); @@ -157,6 +259,33 @@ Tensor<1,dim> Solution::gradient (const Point &p, + // Besides the function that + // represents the exact solution, we + // also need a function which we can + // use as right hand side when + // assembling the linear system of + // discretized equations. This is + // accomplished using the following + // class and the following definition + // of its function. Note that here we + // only need the value of the + // function, not its gradients or + // higher derivatives. +template +class RightHandSide : public Function, + protected SolutionBase +{ + public: + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + // The value of the right hand side + // is given by the negative Laplacian + // of the solution plus the solution + // itself, since we wanted to solve + // Helmholtz's equation: template double RightHandSide::value (const Point &p, const unsigned int) const @@ -166,8 +295,14 @@ double RightHandSide::value (const Point &p, { const Point shifted_point = p-source_centers[i]; - return_value += ((2*dim - 4*shifted_point.square()/(width*width)) / (width*width) * + // The first contribution is + // the Laplacian: + return_value += ((2*dim - 4*shifted_point.square()/(width*width)) / + (width*width) * exp(-shifted_point.square() / (width*width))); + // And the second is the + // solution itself: + return_value += exp(-shifted_point.square() / (width*width)); }; return return_value; @@ -175,6 +310,47 @@ double RightHandSide::value (const Point &p, + // Then we need the class that does + // all the work. +//....................... +template +class LaplaceProblem +{ + public: + enum RefinementMode { + global_refinement, adaptive_refinement + }; + + LaplaceProblem (const FiniteElement &fe, + const RefinementMode refinement_mode); + ~LaplaceProblem (); + void run (); + + private: + void setup_system (); + void assemble_system (); + void solve (); + void refine_grid (); + void process_solution (const unsigned int cycle) const; + + Triangulation triangulation; + DoFHandler dof_handler; + //... + SmartPointer > fe; + ConstraintMatrix hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + + RefinementMode refinement_mode; +}; + + + + template LaplaceProblem::LaplaceProblem (const FiniteElement &fe, const RefinementMode refinement_mode) : @@ -192,13 +368,42 @@ LaplaceProblem::~LaplaceProblem () }; - + // The following function sets up the + // degrees of freedom, sizes of + // matrices and vectors, etc. Most of + // its functionality has been showed + // in previous examples, the only + // difference being the renumbering + // step. template void LaplaceProblem::setup_system () { dof_handler.distribute_dofs (*fe); - // Renumber the degrees of freedom... + // Renumbering the degrees of + // freedom is not overly difficult, + // as long as you use one of the + // algorithms included in the + // library. It requires just one + // line of code, namely the + // following: DoFRenumbering::Cuthill_McKee (dof_handler); + // Note, however, that when you + // renumber the degrees of freedom, + // you must do so immediately after + // distributing them, since such + // things as hanging nodes, the + // sparsity pattern etc. depend on + // the absolute numbers which are + // altered by renumbering. + // + // Renumbering does not serve any + // specific purpose in this + // example, it is done only for + // exposition of the technique. To + // see the effect of renumbering on + // the sparsity pattern of the + // matrix, refer to the second + // example program. hanging_node_constraints.clear (); DoFTools::make_hanging_node_constraints (dof_handler, @@ -220,19 +425,101 @@ void LaplaceProblem::setup_system () + // Assembling the system of equations + // for the problem at hand is mostly + // as for the example programs + // before. However, some things have + // changed anyway, so we comment on + // this function fairly extensively. template void LaplaceProblem::assemble_system () { - QGauss3 quadrature_formula; - FEValues fe_values (*fe, quadrature_formula, - UpdateFlags(update_values | - update_gradients | - update_q_points | - update_JxW_values)); - + // First we need to define objects + // which will be used as quadrature + // formula for domain and face + // integrals. + // + // Note the way in which we define + // a quadrature rule for the faces: + // it is simply a quadrature rule + // for one dimension less! + QGauss3 quadrature_formula; + QGauss3 face_quadrature_formula; + // For simpler use later on, we + // alias the number of quadrature + // points to local variables: + const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + const unsigned int n_face_q_points = face_quadrature_formula.n_quadrature_points; + + // Then we need objects which can + // evaluate the values, gradients, + // etc of the shape functions at + // the quadrature points. While it + // seems that it should be feasible + // to do it with one object for + // both domain and face integrals, + // there is a subtle difference + // since the weights in the domain + // integrals include the measure of + // the cell in the domain, while + // the face integral quadrature + // requires the measure of the face + // in a lower-dimensional + // mannifold. Internally these two + // classes are rooted on a common + // base class which does most of + // the work; that, however, is + // something that you need not + // worry about. + // + // For the domain integrals in the + // bilinear form for Helmholtz's + // equation, we need to compute the + // values and gradients, as well as + // the weights at the quadrature + // points. Furthermore, we need the + // quadrature points on the real + // cell (rather than on the unit + // cell) to evaluate the right hand + // side function. + FEValues fe_values (*fe, quadrature_formula, + UpdateFlags(update_values | + update_gradients | + update_q_points | + update_JxW_values)); + + // For the face integrals, we only + // need the values of the shape + // functions, as well as the + // weights. We also need the normal + // vectors and quadrature points on + // the real cell since we want to + // determine the Neumann values + // from the exact solution object + // (see below). + FEFaceValues fe_face_values (*fe, face_quadrature_formula, + UpdateFlags(update_values | + update_q_points | + update_normal_vectors | + update_JxW_values)); + + // In order to make programming + // more readable below, we alias + // the number of degrees of freedom + // per cell to a local variable, as + // already done for the number of + // quadrature points above: const unsigned int dofs_per_cell = fe->dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + // Then we need some objects + // already known from previous + // examples: An object denoting the + // right hand side function, its + // values at the quadrature points + // on a cell, the cell matrix and + // right hand side, and the indices + // of the degrees of freedom on a + // cell. RightHandSide right_hand_side; vector rhs_values (n_q_points); @@ -241,6 +528,28 @@ void LaplaceProblem::assemble_system () vector local_dof_indices (dofs_per_cell); + // Then we define an object + // denoting the exact solution + // function. We will use it to + // compute the Neumann values at + // the boundary from it. Usually, + // one would of course do so using + // a separate object, in particular + // since the exact solution is not + // known while the Neumann values + // are prescribed. We will, + // however, be a little bit lazy + // and use what we already have in + // information. Real-life programs + // would to go other ways here, of + // course. + Solution exact_solution; + + // Now for the main loop over all + // cells. This is mostly unchanged + // from previous examples, so we + // only comment on the things that + // have changed. DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) @@ -264,16 +573,134 @@ void LaplaceProblem::assemble_system () for (unsigned int i=0; i::faces_per_cell; ++face) + if (cell->face(face)->boundary_indicator() == 1) + { + // If we came into here, + // then we have found an + // external face + // belonging to + // Gamma2. Next, we have + // to compute the values + // of the shape functions + // and the other + // quantities which we + // will need for the + // computation of the + // contour integral. This + // is done using the + // ``reinit'' function + // which we already know + // from the ``FEValue'' + // class: + fe_face_values.reinit (cell, face); + + // Then, for simpler + // access, we alias the + // various quantities to + // local variables: + const FullMatrix + & face_shape_values = fe_face_values.get_shape_values(); + const vector + & face_JxW_values = fe_face_values.get_JxW_values(); + const vector > + & face_q_points = fe_face_values.get_quadrature_points(); + const vector > + & face_normal_vectors = fe_face_values.get_normal_vectors (); + + // And we can then + // perform the + // integration by using a + // loop over all + // quadrature points. + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); for (unsigned int i=0; i::assemble_system () }; }; + // The rest of the function has + // also been shown previously: hanging_node_constraints.condense (system_matrix); hanging_node_constraints.condense (system_rhs); + // Only with the interpolation of + // boundary values, there is one + // notable thing, namely that now + // the boundary indicator for which + // we interpolate boundary values + // (denoted by the second parameter + // to + // ``interpolate_boundary_values'') + // does not represent the whole + // boundary an more. Rather, it is + // that portion of the boundary + // which we have not assigned + // another indicator (see + // below). The degrees of freedom + // at the boundary that do not + // belong to Gamma1 are therefore + // excluded from the interpolation + // of boundary values. map boundary_values; VectorTools::interpolate_boundary_values (dof_handler, 0, @@ -301,7 +748,8 @@ void LaplaceProblem::assemble_system () }; - + // Solving the system of equations is + // done in the same way as before. template void LaplaceProblem::solve () { @@ -355,8 +803,6 @@ void LaplaceProblem::refine_grid () }; -#include - template void LaplaceProblem::process_solution (const unsigned int cycle) const { @@ -402,15 +848,105 @@ void LaplaceProblem::process_solution (const unsigned int cycle) const + // The following function is the main + // one which controls the flow of + // execution. The basic layout is as + // in previous examples: an outer + // loop over successively refined + // grids, and in this loop first + // problem setup, assemblage of the + // linear system, solution, and + // postprocessing. template void LaplaceProblem::run () { - for (unsigned int cycle=0; cycle<12; ++cycle) + for (unsigned int cycle=0; cycle<9; ++cycle) { + // The first action in each + // iteration of the outer loop + // is setting up the grid on + // which we will solve in this + // iteration. In the first + // iteration, the coarsest grid + // is generated, in later + // iterations it is refined, + // for which we call the + // ``refine_grid'' function. if (cycle == 0) { + // Setting up the coarse + // grid is done as in + // previous examples: we + // first create an initial + // grid, which is the unit + // square [-1,1]x[-1,1] in + // the present case. Then + // we refine it globally a + // specific number of + // times. GridGenerator::hyper_cube (triangulation, -1, 1); triangulation.refine_global (1); + + // However, here we have to + // do something else in + // addition: mark those + // faces that belong to the + // different components of + // the boundary, Gamma1 and + // Gamma2. We will use the + // following convention: + // Faces belonging to + // Gamma1 will have the + // boundary indicator ``0'' + // (which is the default, + // so we don't have to set + // it explicitely), and + // faces belonging to + // Gamma2 will use ``1'' as + // boundary indicator. + // + // To set these values, we + // loop over all cells, + // then over all faces of a + // given cell, check + // whether it belongs to + // the boundary Gamma2, and + // if so set its boundary + // indicator to ``1''. + // + // It is worth noting that + // we have to loop over all + // cells here, not only the + // active ones. The reason + // is that upon refinement, + // newly created faces + // inherit the boundary + // indicator of their + // parent face. If we now + // only set the boundary + // indicator for active + // faces, coarsen some + // cells and refine them + // later on, they will + // again have the boundary + // indicator of the parent + // cell which we have not + // modified, instead of the + // one we + // intended. Therefore, we + // have to change the + // boundary indicators of + // all faces on Gamma2, + // irrespective whether + // they are active or not. + Triangulation::cell_iterator cell = triangulation.begin (), + endc = triangulation.end(); + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if ((cell->face(face)->center()(0) == -1) + || + (cell->face(face)->center()(1) == -1)) + cell->face(face)->set_boundary_indicator (1); } else refine_grid (); -- 2.39.5