From d3b85902271dd330373dd690b50265eb76e84070 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 8 Feb 2006 18:12:10 +0000 Subject: [PATCH] Read over first part. git-svn-id: https://svn.dealii.org/trunk@12269 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-8/step-8.cc | 456 ++++++++++++++++++------------ 1 file changed, 270 insertions(+), 186 deletions(-) diff --git a/deal.II/examples/step-8/step-8.cc b/deal.II/examples/step-8/step-8.cc index 8fe13b97b3..1601d595d7 100644 --- a/deal.II/examples/step-8/step-8.cc +++ b/deal.II/examples/step-8/step-8.cc @@ -11,6 +11,8 @@ /* to the file deal.II/doc/license.html for the text and */ /* further information on this license. */ + // @sect3{Include files} + // As usual, the first few include // files are already known, so we // will not comment on them further. @@ -54,12 +56,30 @@ #include + // @sect3{The ``ElasticProblem'' class template} + // The main class is, except for its // name, almost unchanged with - // respect to the step-6 example. The - // only change is the use of a + // respect to the step-6 example. + // + // The only change is the use of a // different class for the ``fe'' - // variable. + // variable: Instead of a concrete + // finite element class such as + // ``FE_Q'', we now use a more + // generic one, ``FESystem''. In + // fact, ``FESystem'' is not really a + // finite element itself in that it + // does not implement shape functions + // of its own. Rather, it is a class + // that can be used to stack several + // other elements together to form + // one vector-valued finite + // element. In our case, we will + // compose the vector-valued element + // of ``FE_Q(1)'' objects, as shown + // below in the constructor of this + // class. template class ElasticProblem { @@ -78,21 +98,6 @@ class ElasticProblem Triangulation triangulation; DoFHandler dof_handler; - // Instead of a concrete finite - // element class such as - // ``FE_Q'', we now use a more - // generic one, ``FESystem''. In - // fact, it is not a finite - // element itself, but rather a - // class that can be used to - // stack several usual elements - // together to form one - // vector-valued finite - // element. In our case, we will - // compose the vector-valued - // element of ``FE_Q(1)'' objects, - // as shown below in the - // constructor of this class. FESystem fe; ConstraintMatrix hanging_node_constraints; @@ -105,6 +110,8 @@ class ElasticProblem }; + // @sect3{Right hand side values} + // Before going over to the // implementation of the main class, // we declare and define the class @@ -112,44 +119,52 @@ class ElasticProblem // side. This time, the right hand // side is vector-valued, as is the // solution, so we will describe the - // new elements in some more detail. + // changes required for this in some + // more detail. + // + // The first thing is that + // vector-valued functions have to + // have a constructor, since they + // need to pass down to the base + // class of how many components the + // function consists. The default + // value in the constructor of the + // base class is one (i.e.: a scalar + // function), which is why we did not + // need not define a constructor for + // the scalar function used in + // previous programs. template class RightHandSide : public Function { public: - // The first thing is that - // vector-valued functions have a - // constructor, since they need - // to pass down to the base class - // of how many components the - // function consists. The default - // value in the constructor of - // the base class is one, so we - // need not define a constructor - // for the usual scalar function. RightHandSide (); - // The next function is a - // replacement for the ``value'' - // function of the previous - // examples. There, a second - // parameter ``component'' was - // given, which denoted which + // The next change is that we + // want a replacement for the + // ``value'' function of the + // previous examples. There, a + // second parameter ``component'' + // was given, which denoted which // component was requested. Here, // we implement a function that // returns the whole vector of // values at the given place at - // once. - virtual void vector_value (const Point &p, - Vector &values) const; - - // Then, in analogy to the + // once, in the second argument + // of the function. The obvious + // name for such a replacement + // function is ``vector_value''. + // + // Secondly, in analogy to the // ``value_list'' function, there // is a function // ``vector_value_list'', which // returns the values of the // vector-valued function at // several points at once: + virtual void vector_value (const Point &p, + Vector &values) const; + virtual void vector_value_list (const std::vector > &points, std::vector > &value_list) const; }; @@ -160,69 +175,110 @@ class RightHandSide : public Function // above, it only passes down to the // base class the number of // components, which is ``dim'' in - // the present case. Note that - // although the implementation is - // very short here, we do not move it - // into the class declaration, since - // our style guides require that - // inside the class declaration only - // declarations have to happen and - // that definitions are always to be - // found outside. + // the present case (one force + // component in each of the ``dim'' + // space directions). + // + // Some people would have moved the + // definition of such a short + // function right into the class + // declaration. We do not do that, as + // a matter of style: the deal.II + // style guides require that class + // declarations contain only + // declarations, and that definitions + // are always to be found + // outside. This is, obviously, as + // much as matter of taste as + // indentation, but we try to be + // consistent in this direction. template -RightHandSide::RightHandSide () : +RightHandSide::RightHandSide () + : Function (dim) {} - // This is the function that returns + // Next the function that returns // the whole vector of values at the - // point ``p'' at once: + // point ``p'' at once. + // + // To prevent cases where the return + // vector has not previously been set + // to the right size we test for this + // case and otherwise throw an + // exception at the beginning of the + // function. Note that enforcing that + // output arguments already have the + // correct size is a convention in + // deal.II, and enforced almost + // everywhere. The reason is that we + // would otherwise have to check at + // the beginning of the function and + // possibly change the size of the + // output vector. This is expensive, + // and would almost always be + // unnecessary (the first call to the + // function would set the vector to + // the right size, and subsequent + // calls would only have to do + // redundant checks). In addition, + // checking and possibly resizing the + // vector is an operation that can + // not be removed if we can't rely on + // the assumption that the vector + // already has the correct size; this + // is in contract to the ``Assert'' + // call that is completely removed if + // the program is compiled in + // optimized mode. + // + // Likewise, if by some accident + // someone tried to compile and run + // the program in only one space + // dimension (in which the elastic + // equations do not make much sense + // since they reduce to the ordinary + // Laplace equation), we terminate + // the program in the second + // assertion. The program will work + // just fine in 3d, however. template inline void RightHandSide::vector_value (const Point &p, Vector &values) const { - // To prevent cases where the - // return value has not previously - // been set to the right size - // (which is kind of a convention - // in the deal.II library), we test - // for this case and otherwise - // throw an exception: Assert (values.size() == dim, ExcDimensionMismatch (values.size(), dim)); - // Likewise, if by some accident - // someone tried to compile and run - // the program in only one space - // dimension (in which the elastic - // equations do not make much sense - // since they reduce to the - // ordinary Laplace equation), we - // terminate the program if the - // dimension is not as expected. - Assert (dim >= 2, ExcInternalError()); + Assert (dim >= 2, ExcNotImplemented()); - // The rest of the function is as - // would probably be expected given - // the form of the right hand side - // function. First we define the - // centers of the two points around - // which are the sources of - // x-displacement, i.e. (0.5,0) and - // (-0.5,0). Note that upon - // construction of the ``Point'' - // objects, all components are set - // to zero. + // The rest of the function + // implements computing force + // values. We will use a constant + // (unit) force in x-direction + // located in two little circles + // (or spheres, in 3d) around + // points (0.5,0) and (-0.5,0), and + // y-force in an area around the + // origin; in 3d, the z-component + // of these centers is zero as + // well. + // + // For this, let us first define + // two objects that denote the + // centers of these areas. Note + // that upon construction of the + // ``Point'' objects, all + // components are set to zero. Point point_1, point_2; point_1(0) = 0.5; point_2(0) = -0.5; - // If now the point ``p'' is in the - // circle of radius 0.2 around one - // of these points, then set the - // force in x-direction to one, - // otherwise to zero: + // If now the point ``p'' is in a + // circle (sphere) of radius 0.2 + // around one of these points, then + // set the force in x-direction to + // one, otherwise to zero: if (((p-point_1).square() < 0.2*0.2) || ((p-point_2).square() < 0.2*0.2)) values(0) = 1; @@ -244,22 +300,24 @@ void RightHandSide::vector_value (const Point &p, // Now, this is the function of the // right hand side class that returns // the values at several points at - // once. + // once. The function starts out with + // checking that the number of input + // and output arguments is equal (the + // sizes of the individual output + // vectors will be checked in the + // function that we call further down + // below). Next, we define an + // abbreviation for the number of + // points which we shall work on, to + // make some things simpler below. template void RightHandSide::vector_value_list (const std::vector > &points, std::vector > &value_list) const { - // First we define an abbreviation - // for the number of points which - // we shall work on: - const unsigned int n_points = points.size(); + Assert (value_list.size() == points.size(), + ExcDimensionMismatch (value_list.size(), points.size())); - // Then we check whether the number - // of output slots has been set - // correctly, i.e. to the number of - // input points: - Assert (value_list.size() == n_points, - ExcDimensionMismatch (value_list.size(), n_points)); + const unsigned int n_points = points.size(); // Finally we treat each of the // points. In one of the previous @@ -274,73 +332,105 @@ void RightHandSide::vector_value_list (const std::vector > &poin // twice, which can lead to // confusion if one function is // changed but the other is - // not. However, we can prevent - // this situation using the - // following construct: + // not. + // + // We can prevent this situation by + // calling + // ``RightHandSide::vector_valued'' + // on each point in the input + // list. Note that by giving the + // full name of the function, + // including the class name, we + // instruct the compiler to + // explicitly call this function, + // and not to use the virtual + // function call mechanism that + // would be used if we had just + // called ``vector_value''. This is + // important, since the compiler + // generally can't make any + // assumptions which function is + // called when using virtual + // functions, and it therefore + // can't inline the called function + // into the site of the call. On + // the contrary, here we give the + // fully qualified name, which + // bypasses the virtual function + // call, and consequently the + // compiler knows exactly which + // function is called and will + // inline above function into the + // present location. (Note that we + // have declared the + // ``vector_value'' function above + // ``inline'', though modern + // compilers are also able to + // inline functions even if they + // have not been declared as + // inline). + // + // It is worth noting why we go to + // such length explaining what we + // do. Using this construct, we + // manage to avoid any + // inconsistency: if we want to + // change the right hand side + // function, it would be difficult + // to always remember that we + // always have to change two + // functions in the same way. Using + // this forwarding mechanism, we + // only have to change a single + // place (the ``vector_value'' + // function), and the second place + // (the ``vector_value_list'' + // function) will always be + // consistent with it. At the same + // time, using virtual function + // call bypassing, the code is no + // less efficient than if we had + // written it twice in the first + // place: for (unsigned int p=0; p::vector_value (points[p], value_list[p]); - // It calls the ``vector_value'' - // function defined above for each - // point, and thus preempts all - // chances for inconsistency. It is - // important to note how the - // function was called: using the - // full class qualification using - // ``RightHandSide::'', since this - // calls the function directly and - // not using the virtual function - // table. The call is thus as fast - // as a call to any non-virtual - // function. In addition, we have - // declared the ``vector_value'' - // function ``inline'', i.e. the - // compiler can remove the function - // call altogether and the - // resulting code can in principle - // be as fast as if we had - // duplicated the code. } + // @sect3{The ``ElasticProblem'' class implementation} + + // @sect4{ElasticProblem::ElasticProblem} + + // Following is the constructor of + // the main class. As said before, we + // would like to construct a + // vector-valued finite element that + // is composed of several scalar + // finite elements (i.e., we want to + // build the vector-valued element so + // that each of its vector components + // consists of the shape functions of + // a scalar element). Of course, the + // number of scalar finite elements we + // would like to stack together + // equals the number of components + // the solution function has, which + // is ``dim'' since we consider + // displacement in each space + // direction. The ``FESystem'' class + // can handle this: we pass it the + // finite element of which we would + // like to compose the system of, and + // how often it shall be repeated: template -ElasticProblem::ElasticProblem () : +ElasticProblem::ElasticProblem () + : dof_handler (triangulation), - // As said before, we - // would like to - // construct one - // vector-valued - // finite element as - // outer product of - // several scalar - // finite - // elements. Of - // course, the number - // of scalar finite - // element we would - // like to stack - // together equals - // the number of - // components the - // solution function - // has, which is - // ``dim'' since we - // consider - // displacement in - // each space - // direction. The - // ``FESystem'' class - // can handle this: - // we pass it the - // finite element of - // which we would - // like to compose - // the system of, and - // how often it shall - // be repeated: fe (FE_Q(1), dim) +{} // In fact, the ``FESystem'' class // has several more constructors // which can perform more complex @@ -349,21 +439,13 @@ ElasticProblem::ElasticProblem () : // elements of the same type into // one; we will get to know these // possibilities in later examples. - // - // It should be noted that the - // ``FESystem'' object so created - // does not actually use the finite - // element which we have passed to it - // as first parameter. We could thus - // use an anonymous object created - // in-place. The ``FESystem'' - // constructor only needs the - // parameter to generate a copy of - // the finite element from this. -{} + // @sect4{ElasticProblem::~ElasticProblem} + + // The destructor, on the other hand, + // is exactly as in step-6: template ElasticProblem::~ElasticProblem () { @@ -371,23 +453,29 @@ ElasticProblem::~ElasticProblem () } + // @sect4{ElasticProblem::setup_system} + // Setting up the system of equations - // is equal to the function used in - // the step-6 example. The + // is identitical to the function + // used in the step-6 example. The // ``DoFHandler'' class and all other - // classes used take care of the - // vector-valuedness of the finite - // element themselves (in fact, the - // do not do so, since they only take - // care how many degrees of freedom - // there are per vertex, line and - // cell, and they do not ask what they - // represent, i.e. whether the finite - // element under consideration is - // vector-valued or whether it is, - // for example, a scalar Hermite - // element with several degrees of - // freedom on each vertex). + // classes used here are fully aware + // that the finite element we want to + // use is vector-valued, and take + // care of the vector-valuedness of + // the finite element themselves. (In + // fact, they do not, but this does + // not need to bother you: since they + // only need to know how many degrees + // of freedom there are per vertex, + // line and cell, and they do not ask + // what they represent, i.e. whether + // the finite element under + // consideration is vector-valued or + // whether it is, for example, a + // scalar Hermite element with + // several degrees of freedom on each + // vertex). template void ElasticProblem::setup_system () { @@ -399,13 +487,6 @@ void ElasticProblem::setup_system () sparsity_pattern.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), dof_handler.max_couplings_between_dofs()); - // When making the sparsity - // pattern, there is some potential - // for optimization if not all - // components couple to all - // others. However, this is not the - // case for the elastic equations, - // so we use the standard call: DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); hanging_node_constraints.condense (sparsity_pattern); @@ -419,6 +500,8 @@ void ElasticProblem::setup_system () } + // @sect4{ElasticProblem::assemble_system} + // The big changes in this program // are in the creation of matrix and // right hand side, since they are @@ -1058,6 +1141,7 @@ void ElasticProblem::run () }; } + // @sect3{The ``main'' function} // The main function is again exactly // like in step-6 (apart from the -- 2.39.5