From 08e10f34d45b62f854c8e63cd48003ab96aa2105 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sun, 17 Jul 2016 15:28:30 -0400 Subject: [PATCH] Use a plain function, not a Function, in step-8. --- examples/step-8/step-8.cc | 171 +++++++------------------------------- 1 file changed, 28 insertions(+), 143 deletions(-) diff --git a/examples/step-8/step-8.cc b/examples/step-8/step-8.cc index 1dc6ed1dfa..9df1bba196 100644 --- a/examples/step-8/step-8.cc +++ b/examples/step-8/step-8.cc @@ -24,7 +24,7 @@ // comment on them further. #include #include -#include +#include #include #include #include @@ -111,62 +111,10 @@ namespace Step8 // @sect3{Right hand side values} // Before going over to the implementation of the main class, we declare and - // define the class which describes the right hand side. This time, the + // define the function which describes the right hand side. This time, the // right hand side is vector-valued, as is the solution, so we will describe // the 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: - RightHandSide (); - - // 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, 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; - }; - - - // This is the constructor of the right hand side class. As said above, it - // only passes down to the base class the number of components, which is - // dim in 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 () - : - Function (dim) - {} - - - // Next the function that returns the whole vector of values at the 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 @@ -187,12 +135,11 @@ namespace Step8 // 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 + void right_hand_side (const std::vector > &points, + std::vector > &values) { - Assert (values.size() == dim, - ExcDimensionMismatch (values.size(), dim)); + Assert (values.size() == points.size(), + ExcDimensionMismatch (values.size(), points.size())); Assert (dim >= 2, ExcNotImplemented()); // The rest of the function implements computing force values. We will use @@ -208,77 +155,24 @@ namespace Step8 point_1(0) = 0.5; point_2(0) = -0.5; - // 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).norm_square() < 0.2*0.2) || - ((p-point_2).norm_square() < 0.2*0.2)) - values(0) = 1; - else - values(0) = 0; - - // Likewise, if p is in the vicinity of the origin, then set - // the y-force to 1, otherwise to zero: - if (p.norm_square() < 0.2*0.2) - values(1) = 1; - else - values(1) = 0; - } - - + for (unsigned int point_n = 0; point_n < points.size(); ++point_n) + { + // If points[point_n] 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 (((points[point_n]-point_1).norm_square() < 0.2*0.2) || + ((points[point_n]-point_2).norm_square() < 0.2*0.2)) + values[point_n][0] = 1.0; + else + values[point_n][0] = 0.0; - // Now, this is the function of the right hand side class that returns the - // values at several points at 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 - { - Assert (value_list.size() == points.size(), - ExcDimensionMismatch (value_list.size(), points.size())); - - const unsigned int n_points = points.size(); - - // Finally we treat each of the points. In one of the previous examples, - // we have explained why the - // value_list/vector_value_list function had - // been introduced: to prevent us from calling virtual functions too - // frequently. On the other hand, we now need to implement the same - // function twice, which can lead to confusion if one function is changed - // but the other is not. - // - // We can prevent this situation by calling - // RightHandSide::vector_value 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]); + // Likewise, if points[point_n] is in the vicinity of the + // origin, then set the y-force to one, otherwise to zero: + if (points[point_n].norm_square() < 0.2*0.2) + values[point_n][1] = 1.0; + else + values[point_n][1] = 0.0; + } } @@ -405,16 +299,9 @@ namespace Step8 // demonstration. ConstantFunction lambda(1.), mu(1.); - // Then again, we need to have the same for the right hand side. This is - // exactly as before in previous examples. However, we now have a - // vector-valued right hand side, which is why the data type of the - // rhs_values array is changed. We initialize it by - // n_q_points elements, each of which is a - // Vector@ with dim elements. - RightHandSide right_hand_side; - std::vector > rhs_values (n_q_points, - Vector(dim)); - + // Like the two constant functions above, we will call the function + // right_hand_side just once per cell to make things simpler. + std::vector > rhs_values (n_q_points); // Now we can begin with the loop over all cells: typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), @@ -430,9 +317,7 @@ namespace Step8 // points. Likewise for the right hand side: lambda.value_list (fe_values.get_quadrature_points(), lambda_values); mu.value_list (fe_values.get_quadrature_points(), mu_values); - - right_hand_side.vector_value_list (fe_values.get_quadrature_points(), - rhs_values); + right_hand_side (fe_values.get_quadrature_points(), rhs_values); // Then assemble the entries of the local stiffness matrix and right // hand side vector. This follows almost one-to-one the pattern @@ -522,7 +407,7 @@ namespace Step8 for (unsigned int q_point=0; q_point