From 1ae52b8bf25b3de494b5061f58ba5b33162994b4 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 19 Aug 2011 04:39:39 +0000 Subject: [PATCH] Do like in step-7: put everything into a namespace. git-svn-id: https://svn.dealii.org/trunk@24113 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-8/step-8.cc | 2157 +++++++++++++++-------------- 1 file changed, 1083 insertions(+), 1074 deletions(-) diff --git a/deal.II/examples/step-8/step-8.cc b/deal.II/examples/step-8/step-8.cc index aecef3b529..d617f803df 100644 --- a/deal.II/examples/step-8/step-8.cc +++ b/deal.II/examples/step-8/step-8.cc @@ -3,7 +3,7 @@ /* $Id$ */ /* */ -/* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2010 by the deal.II authors */ +/* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2010, 2011 by the deal.II authors */ /* */ /* This file is subject to QPL and may not be distributed */ /* without copyright and license information. Please refer */ @@ -55,1097 +55,1106 @@ #include // The last step is as in previous - // programs: -using namespace dealii; - - // @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 - // different class for the fe - // 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 + // programs. In particular, just like in + // step-7, we pack everything that's specific + // to this program into a namespace of its + // own. +namespace Step8 { - public: - ElasticProblem (); - ~ElasticProblem (); - void run (); - - private: - void setup_system (); - void assemble_system (); - void solve (); - void refine_grid (); - void output_results (const unsigned int cycle) const; - - Triangulation triangulation; - DoFHandler dof_handler; - - FESystem fe; - - ConstraintMatrix hanging_node_constraints; - - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; - - Vector solution; - Vector system_rhs; -}; - - - // @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 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 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 -{ - Assert (values.size() == dim, - ExcDimensionMismatch (values.size(), dim)); - Assert (dim >= 2, ExcNotImplemented()); - - // 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 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; - 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.square() < 0.2*0.2) - values(1) = 1; - else - values(1) = 0; -} + using namespace dealii; + // @sect3{The ElasticProblem class template} - - // 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. + // The main class is, except for its + // name, almost unchanged with + // respect to the step-6 example. // - // 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). + // The only change is the use of a + // different class for the fe + // 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 + { + public: + ElasticProblem (); + ~ElasticProblem (); + void run (); + + private: + void setup_system (); + void assemble_system (); + void solve (); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + DoFHandler dof_handler; + + FESystem fe; + + ConstraintMatrix hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + }; + + + // @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 right hand + // side is vector-valued, as is the + // solution, so we will describe the + // changes required for this in some + // more detail. // - // 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]); -} - - - - // @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 () - : - dof_handler (triangulation), - fe (FE_Q(1), dim) -{} - // In fact, the FESystem class - // has several more constructors - // which can perform more complex - // operations than just stacking - // together several scalar finite - // elements of the same type into - // one; we will get to know these - // possibilities in later examples. - - - - // @sect4{ElasticProblem::~ElasticProblem} - - // The destructor, on the other hand, - // is exactly as in step-6: -template -ElasticProblem::~ElasticProblem () -{ - dof_handler.clear (); -} - - - // @sect4{ElasticProblem::setup_system} - - // Setting up the system of equations - // is identitical to the function - // used in the step-6 example. The - // DoFHandler class and all other - // 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 () -{ - dof_handler.distribute_dofs (fe); - hanging_node_constraints.clear (); - DoFTools::make_hanging_node_constraints (dof_handler, - hanging_node_constraints); - hanging_node_constraints.close (); - 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); - - hanging_node_constraints.condense (sparsity_pattern); - - sparsity_pattern.compress(); - - system_matrix.reinit (sparsity_pattern); - - solution.reinit (dof_handler.n_dofs()); - system_rhs.reinit (dof_handler.n_dofs()); -} - - - // @sect4{ElasticProblem::assemble_system} - - // The big changes in this program - // are in the creation of matrix and - // right hand side, since they are - // problem-dependent. We will go - // through that process step-by-step, - // since it is a bit more complicated - // than in previous examples. - // - // The first parts of this function - // are the same as before, however: - // setting up a suitable quadrature - // formula, initializing an - // FEValues object for the - // (vector-valued) finite element we - // use as well as the quadrature - // object, and declaring a number of - // auxiliary arrays. In addition, we - // declare the ever same two - // abbreviations: n_q_points and - // dofs_per_cell. The number of - // degrees of freedom per cell we now - // obviously ask from the composed - // finite element rather than from - // the underlying scalar Q1 - // element. Here, it is dim times - // the number of degrees of freedom - // per cell of the Q1 element, though - // this is not explicit knowledge we - // need to care about: -template -void ElasticProblem::assemble_system () -{ - QGauss quadrature_formula(2); - - FEValues fe_values (fe, quadrature_formula, - update_values | update_gradients | - update_quadrature_points | update_JxW_values); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - Vector cell_rhs (dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - // As was shown in previous - // examples as well, we need a - // place where to store the values - // of the coefficients at all the - // quadrature points on a cell. In - // the present situation, we have - // two coefficients, lambda and mu. - std::vector lambda_values (n_q_points); - std::vector mu_values (n_q_points); - - // Well, we could as well have - // omitted the above two arrays - // since we will use constant - // coefficients for both lambda and - // mu, which can be declared like - // this. They both represent - // functions always returning the - // constant value 1.0. Although we - // could omit the respective - // factors in the assemblage of the - // matrix, we use them here for - // purpose of 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)); - - - // Now we can begin with the loop - // over all cells: - typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - cell_matrix = 0; - cell_rhs = 0; - - fe_values.reinit (cell); - - // Next we get the values of - // the coefficients at the - // quadrature 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); - - // Then assemble the entries of - // the local stiffness matrix - // and right hand side - // vector. This follows almost - // one-to-one the pattern - // described in the - // introduction of this - // example. One of the few - // comments in place is that we - // can compute the number - // comp(i), i.e. the index - // of the only nonzero vector - // component of shape function - // i using the - // fe.system_to_component_index(i).first - // function call below. + // 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. // - // (By accessing the - // first variable of - // the return value of the - // system_to_component_index - // function, you might - // already have guessed - // that there is more in - // it. In fact, the - // function returns a - // std::pair@, of - // which the first element - // is comp(i) and the - // second is the value - // base(i) also noted - // in the introduction, i.e. - // the index - // of this shape function - // within all the shape - // functions that are nonzero - // in this component, - // i.e. base(i) in the - // diction of the - // introduction. This is not a - // number that we are usually - // interested in, however.) - // - // With this knowledge, we can - // assemble the local matrix - // contributions: - for (unsigned int i=0; ishape_grad(i,q_point) - // returns the - // gradient of - // the only - // nonzero - // component of - // the i-th shape - // function at - // quadrature - // point - // q_point. The - // component - // comp(i) of - // the gradient, - // which is the - // derivative of - // this only - // nonzero vector - // component of - // the i-th shape - // function with - // respect to the - // comp(i)th - // coordinate is - // accessed by - // the appended - // brackets. - ( - (fe_values.shape_grad(i,q_point)[component_i] * - fe_values.shape_grad(j,q_point)[component_j] * - lambda_values[q_point]) - + - (fe_values.shape_grad(i,q_point)[component_j] * - fe_values.shape_grad(j,q_point)[component_i] * - mu_values[q_point]) - + - // The second term is - // (mu nabla u_i, nabla v_j). - // We need not - // access a - // specific + // 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 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 + { + Assert (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + Assert (dim >= 2, ExcNotImplemented()); + + // 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 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; + 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.square() < 0.2*0.2) + values(1) = 1; + else + values(1) = 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]); + } + + + + // @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 () + : + dof_handler (triangulation), + fe (FE_Q(1), dim) + {} + // In fact, the FESystem class + // has several more constructors + // which can perform more complex + // operations than just stacking + // together several scalar finite + // elements of the same type into + // one; we will get to know these + // possibilities in later examples. + + + + // @sect4{ElasticProblem::~ElasticProblem} + + // The destructor, on the other hand, + // is exactly as in step-6: + template + ElasticProblem::~ElasticProblem () + { + dof_handler.clear (); + } + + + // @sect4{ElasticProblem::setup_system} + + // Setting up the system of equations + // is identitical to the function + // used in the step-6 example. The + // DoFHandler class and all other + // 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 () + { + dof_handler.distribute_dofs (fe); + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, + hanging_node_constraints); + hanging_node_constraints.close (); + 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); + + hanging_node_constraints.condense (sparsity_pattern); + + sparsity_pattern.compress(); + + system_matrix.reinit (sparsity_pattern); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + } + + + // @sect4{ElasticProblem::assemble_system} + + // The big changes in this program + // are in the creation of matrix and + // right hand side, since they are + // problem-dependent. We will go + // through that process step-by-step, + // since it is a bit more complicated + // than in previous examples. + // + // The first parts of this function + // are the same as before, however: + // setting up a suitable quadrature + // formula, initializing an + // FEValues object for the + // (vector-valued) finite element we + // use as well as the quadrature + // object, and declaring a number of + // auxiliary arrays. In addition, we + // declare the ever same two + // abbreviations: n_q_points and + // dofs_per_cell. The number of + // degrees of freedom per cell we now + // obviously ask from the composed + // finite element rather than from + // the underlying scalar Q1 + // element. Here, it is dim times + // the number of degrees of freedom + // per cell of the Q1 element, though + // this is not explicit knowledge we + // need to care about: + template + void ElasticProblem::assemble_system () + { + QGauss quadrature_formula(2); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + // As was shown in previous + // examples as well, we need a + // place where to store the values + // of the coefficients at all the + // quadrature points on a cell. In + // the present situation, we have + // two coefficients, lambda and mu. + std::vector lambda_values (n_q_points); + std::vector mu_values (n_q_points); + + // Well, we could as well have + // omitted the above two arrays + // since we will use constant + // coefficients for both lambda and + // mu, which can be declared like + // this. They both represent + // functions always returning the + // constant value 1.0. Although we + // could omit the respective + // factors in the assemblage of the + // matrix, we use them here for + // purpose of 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)); + + + // Now we can begin with the loop + // over all cells: + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_matrix = 0; + cell_rhs = 0; + + fe_values.reinit (cell); + + // Next we get the values of + // the coefficients at the + // quadrature 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); + + // Then assemble the entries of + // the local stiffness matrix + // and right hand side + // vector. This follows almost + // one-to-one the pattern + // described in the + // introduction of this + // example. One of the few + // comments in place is that we + // can compute the number + // comp(i), i.e. the index + // of the only nonzero vector + // component of shape function + // i using the + // fe.system_to_component_index(i).first + // function call below. + // + // (By accessing the + // first variable of + // the return value of the + // system_to_component_index + // function, you might + // already have guessed + // that there is more in + // it. In fact, the + // function returns a + // std::pair@, of + // which the first element + // is comp(i) and the + // second is the value + // base(i) also noted + // in the introduction, i.e. + // the index + // of this shape function + // within all the shape + // functions that are nonzero + // in this component, + // i.e. base(i) in the + // diction of the + // introduction. This is not a + // number that we are usually + // interested in, however.) + // + // With this knowledge, we can + // assemble the local matrix + // contributions: + for (unsigned int i=0; ishape_grad(i,q_point) + // returns the + // gradient of + // the only + // nonzero // component of - // the - // gradient, - // since we - // only have to - // compute the - // scalar - // product of - // the two - // gradients, - // of which an - // overloaded - // version of - // the - // operator* - // takes care, - // as in - // previous - // examples. - // - // Note that by - // using the ?: - // operator, we - // only do this - // if comp(i) - // equals - // comp(j), - // otherwise a - // zero is - // added (which - // will be - // optimized - // away by the - // compiler). - ((component_i == component_j) ? - (fe_values.shape_grad(i,q_point) * - fe_values.shape_grad(j,q_point) * - mu_values[q_point]) : - 0) - ) - * - fe_values.JxW(q_point); - } - } - } - - // Assembling the right hand - // side is also just as - // discussed in the - // introduction: - for (unsigned int i=0; iget_dof_indices (local_dof_indices); - for (unsigned int i=0; iZeroFunction - // constructor accepts a parameter - // that tells it that it shall - // represent a vector valued, - // constant zero function with that - // many components. By default, - // this parameter is equal to one, - // in which case the - // ZeroFunction object would - // represent a scalar - // function. Since the solution - // vector has dim components, - // we need to pass dim as - // number of components to the zero - // function as well. - std::map boundary_values; - VectorTools::interpolate_boundary_values (dof_handler, - 0, - ZeroFunction(dim), - boundary_values); - MatrixTools::apply_boundary_values (boundary_values, - system_matrix, - solution, - system_rhs); -} - - - - // @sect4{ElasticProblem::solve} - - // The solver does not care about - // where the system of equations - // comes, as long as it stays - // positive definite and symmetric - // (which are the requirements for - // the use of the CG solver), which - // the system indeed is. Therefore, - // we need not change anything. -template -void ElasticProblem::solve () -{ - SolverControl solver_control (1000, 1e-12); - SolverCG<> cg (solver_control); - - PreconditionSSOR<> preconditioner; - preconditioner.initialize(system_matrix, 1.2); - - cg.solve (system_matrix, solution, system_rhs, - preconditioner); - - hanging_node_constraints.distribute (solution); -} - - - // @sect4{ElasticProblem::refine_grid} - - // The function that does the - // refinement of the grid is the same - // as in the step-6 example. The - // quadrature formula is adapted to - // the linear elements again. Note - // that the error estimator by - // default adds up the estimated - // obtained from all components of - // the finite element solution, i.e., - // it uses the displacement in all - // directions with the same - // weight. If we would like the grid - // to be adapted to the - // x-displacement only, we could pass - // the function an additional - // parameter which tells it to do so - // and do not consider the - // displacements in all other - // directions for the error - // indicators. However, for the - // current problem, it seems - // appropriate to consider all - // displacement components with equal - // weight. -template -void ElasticProblem::refine_grid () -{ - Vector estimated_error_per_cell (triangulation.n_active_cells()); - - typename FunctionMap::type neumann_boundary; - KellyErrorEstimator::estimate (dof_handler, - QGauss(2), - neumann_boundary, - solution, - estimated_error_per_cell); - - GridRefinement::refine_and_coarsen_fixed_number (triangulation, - estimated_error_per_cell, - 0.3, 0.03); - - triangulation.execute_coarsening_and_refinement (); -} - - - // @sect4{ElasticProblem::output_results} - - // The output happens mostly as has - // been shown in previous examples - // already. The only difference is - // that the solution function is - // vector valued. The DataOut - // class takes care of this - // automatically, but we have to give - // each component of the solution - // vector a different name. -template -void ElasticProblem::output_results (const unsigned int cycle) const -{ - std::string filename = "solution-"; - filename += ('0' + cycle); - Assert (cycle < 10, ExcInternalError()); - - filename += ".gmv"; - std::ofstream output (filename.c_str()); - - DataOut data_out; - data_out.attach_dof_handler (dof_handler); - - - - // As said above, we need a - // different name for each - // component of the solution - // function. To pass one name for - // each component, a vector of - // strings is used. Since the - // number of components is the same - // as the number of dimensions we - // are working in, the following - // switch statement is used. + // the i-th shape + // function at + // quadrature + // point + // q_point. The + // component + // comp(i) of + // the gradient, + // which is the + // derivative of + // this only + // nonzero vector + // component of + // the i-th shape + // function with + // respect to the + // comp(i)th + // coordinate is + // accessed by + // the appended + // brackets. + ( + (fe_values.shape_grad(i,q_point)[component_i] * + fe_values.shape_grad(j,q_point)[component_j] * + lambda_values[q_point]) + + + (fe_values.shape_grad(i,q_point)[component_j] * + fe_values.shape_grad(j,q_point)[component_i] * + mu_values[q_point]) + + + // The second term is + // (mu nabla u_i, nabla v_j). + // We need not + // access a + // specific + // component of + // the + // gradient, + // since we + // only have to + // compute the + // scalar + // product of + // the two + // gradients, + // of which an + // overloaded + // version of + // the + // operator* + // takes care, + // as in + // previous + // examples. + // + // Note that by + // using the ?: + // operator, we + // only do this + // if comp(i) + // equals + // comp(j), + // otherwise a + // zero is + // added (which + // will be + // optimized + // away by the + // compiler). + ((component_i == component_j) ? + (fe_values.shape_grad(i,q_point) * + fe_values.shape_grad(j,q_point) * + mu_values[q_point]) : + 0) + ) + * + fe_values.JxW(q_point); + } + } + } + + // Assembling the right hand + // side is also just as + // discussed in the + // introduction: + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + for (unsigned int i=0; iZeroFunction + // constructor accepts a parameter + // that tells it that it shall + // represent a vector valued, + // constant zero function with that + // many components. By default, + // this parameter is equal to one, + // in which case the + // ZeroFunction object would + // represent a scalar + // function. Since the solution + // vector has dim components, + // we need to pass dim as + // number of components to the zero + // function as well. + std::map boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(dim), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); + } + + + + // @sect4{ElasticProblem::solve} + + // The solver does not care about + // where the system of equations + // comes, as long as it stays + // positive definite and symmetric + // (which are the requirements for + // the use of the CG solver), which + // the system indeed is. Therefore, + // we need not change anything. + template + void ElasticProblem::solve () + { + SolverControl solver_control (1000, 1e-12); + SolverCG<> cg (solver_control); + + PreconditionSSOR<> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + + hanging_node_constraints.distribute (solution); + } + + + // @sect4{ElasticProblem::refine_grid} + + // The function that does the + // refinement of the grid is the same + // as in the step-6 example. The + // quadrature formula is adapted to + // the linear elements again. Note + // that the error estimator by + // default adds up the estimated + // obtained from all components of + // the finite element solution, i.e., + // it uses the displacement in all + // directions with the same + // weight. If we would like the grid + // to be adapted to the + // x-displacement only, we could pass + // the function an additional + // parameter which tells it to do so + // and do not consider the + // displacements in all other + // directions for the error + // indicators. However, for the + // current problem, it seems + // appropriate to consider all + // displacement components with equal + // weight. + template + void ElasticProblem::refine_grid () + { + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + typename FunctionMap::type neumann_boundary; + KellyErrorEstimator::estimate (dof_handler, + QGauss(2), + neumann_boundary, + solution, + estimated_error_per_cell); + + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.3, 0.03); + + triangulation.execute_coarsening_and_refinement (); + } + + + // @sect4{ElasticProblem::output_results} + + // The output happens mostly as has + // been shown in previous examples + // already. The only difference is + // that the solution function is + // vector valued. The DataOut + // class takes care of this + // automatically, but we have to give + // each component of the solution + // vector a different name. + template + void ElasticProblem::output_results (const unsigned int cycle) const + { + std::string filename = "solution-"; + filename += ('0' + cycle); + Assert (cycle < 10, ExcInternalError()); + + filename += ".gmv"; + std::ofstream output (filename.c_str()); + + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + + + + // As said above, we need a + // different name for each + // component of the solution + // function. To pass one name for + // each component, a vector of + // strings is used. Since the + // number of components is the same + // as the number of dimensions we + // are working in, the following + // switch statement is used. + // + // We note that some graphics + // programs have restriction as to + // what characters are allowed in + // the names of variables. The + // library therefore supports only + // the minimal subset of these + // characters that is supported by + // all programs. Basically, these + // are letters, numbers, + // underscores, and some other + // characters, but in particular no + // whitespace and minus/hyphen. The + // library will throw an exception + // otherwise, at least if in debug + // mode. + // + // After listing the 1d, 2d, and 3d + // case, it is good style to let + // the program die if we run upon a + // case which we did not + // consider. Remember that the + // Assert macro generates an + // exception if the condition in + // the first parameter is not + // satisfied. Of course, the + // condition false can never be + // satisfied, so the program will + // always abort whenever it gets to + // the default statement: + std::vector solution_names; + switch (dim) + { + case 1: + solution_names.push_back ("displacement"); + break; + case 2: + solution_names.push_back ("x_displacement"); + solution_names.push_back ("y_displacement"); + break; + case 3: + solution_names.push_back ("x_displacement"); + solution_names.push_back ("y_displacement"); + solution_names.push_back ("z_displacement"); + break; + default: + Assert (false, ExcNotImplemented()); + } + + // After setting up the names for + // the different components of the + // solution vector, we can add the + // solution vector to the list of + // data vectors scheduled for + // output. Note that the following + // function takes a vector of + // strings as second argument, + // whereas the one which we have + // used in all previous examples + // accepted a string there. In + // fact, the latter function is + // only a shortcut for the function + // which we call here: it puts the + // single string that is passed to + // it into a vector of strings with + // only one element and forwards + // that to the other function. + data_out.add_data_vector (solution, solution_names); + data_out.build_patches (); + data_out.write_gmv (output); + } + + + + // @sect4{ElasticProblem::run} + + // The run function does the same + // things as in step-6, for + // example. This time, we use the + // square [-1,1]^d as domain, and we + // refine it twice globally before + // starting the first iteration. // - // We note that some graphics - // programs have restriction as to - // what characters are allowed in - // the names of variables. The - // library therefore supports only - // the minimal subset of these - // characters that is supported by - // all programs. Basically, these - // are letters, numbers, - // underscores, and some other - // characters, but in particular no - // whitespace and minus/hyphen. The - // library will throw an exception - // otherwise, at least if in debug - // mode. + // The reason is the following: we + // use the Gauss quadrature + // formula with two points in each + // direction for integration of the + // right hand side; that means that + // there are four quadrature points + // on each cell (in 2D). If we only + // refine the initial grid once + // globally, then there will be only + // four quadrature points in each + // direction on the domain. However, + // the right hand side function was + // chosen to be rather localized and + // in that case all quadrature points + // lie outside the support of the + // right hand side function. The + // right hand side vector will then + // contain only zeroes and the + // solution of the system of + // equations is the zero vector, + // i.e. a finite element function + // that it zero everywhere. We should + // not be surprised about such things + // happening, since we have chosen an + // initial grid that is totally + // unsuitable for the problem at + // hand. // - // After listing the 1d, 2d, and 3d - // case, it is good style to let - // the program die if we run upon a - // case which we did not - // consider. Remember that the - // Assert macro generates an - // exception if the condition in - // the first parameter is not - // satisfied. Of course, the - // condition false can never be - // satisfied, so the program will - // always abort whenever it gets to - // the default statement: - std::vector solution_names; - switch (dim) - { - case 1: - solution_names.push_back ("displacement"); - break; - case 2: - solution_names.push_back ("x_displacement"); - solution_names.push_back ("y_displacement"); - break; - case 3: - solution_names.push_back ("x_displacement"); - solution_names.push_back ("y_displacement"); - solution_names.push_back ("z_displacement"); - break; - default: - Assert (false, ExcNotImplemented()); - } - - // After setting up the names for - // the different components of the - // solution vector, we can add the - // solution vector to the list of - // data vectors scheduled for - // output. Note that the following - // function takes a vector of - // strings as second argument, - // whereas the one which we have - // used in all previous examples - // accepted a string there. In - // fact, the latter function is - // only a shortcut for the function - // which we call here: it puts the - // single string that is passed to - // it into a vector of strings with - // only one element and forwards - // that to the other function. - data_out.add_data_vector (solution, solution_names); - data_out.build_patches (); - data_out.write_gmv (output); -} - - - - // @sect4{ElasticProblem::run} - - // The run function does the same - // things as in step-6, for - // example. This time, we use the - // square [-1,1]^d as domain, and we - // refine it twice globally before - // starting the first iteration. - // - // The reason is the following: we - // use the Gauss quadrature - // formula with two points in each - // direction for integration of the - // right hand side; that means that - // there are four quadrature points - // on each cell (in 2D). If we only - // refine the initial grid once - // globally, then there will be only - // four quadrature points in each - // direction on the domain. However, - // the right hand side function was - // chosen to be rather localized and - // in that case all quadrature points - // lie outside the support of the - // right hand side function. The - // right hand side vector will then - // contain only zeroes and the - // solution of the system of - // equations is the zero vector, - // i.e. a finite element function - // that it zero everywhere. We should - // not be surprised about such things - // happening, since we have chosen an - // initial grid that is totally - // unsuitable for the problem at - // hand. - // - // The unfortunate thing is that if - // the discrete solution is constant, - // then the error indicators computed - // by the KellyErrorEstimator - // class are zero for each cell as - // well, and the call to - // refine_and_coarsen_fixed_number - // on the triangulation object - // will not flag any cells for - // refinement (why should it if the - // indicated error is zero for each - // cell?). The grid in the next - // iteration will therefore consist - // of four cells only as well, and - // the same problem occurs again. - // - // The conclusion needs to be: while - // of course we will not choose the - // initial grid to be well-suited for - // the accurate solution of the - // problem, we must at least choose - // it such that it has the chance to - // capture the most striking features - // of the solution. In this case, it - // needs to be able to see the right - // hand side. Thus, we refine twice - // globally. (Note that the - // refine_global function is not - // part of the GridRefinement - // class in which - // refine_and_coarsen_fixed_number - // is declared, for example. The - // reason is first that it is not an - // algorithm that computed refinement - // flags from indicators, but more - // importantly that it actually - // performs the refinement, in - // contrast to the functions in - // GridRefinement that only flag - // cells without actually refining - // the grid.) -template -void ElasticProblem::run () -{ - for (unsigned int cycle=0; cycle<8; ++cycle) - { - std::cout << "Cycle " << cycle << ':' << std::endl; - - if (cycle == 0) - { - GridGenerator::hyper_cube (triangulation, -1, 1); - triangulation.refine_global (2); - } - else - refine_grid (); - - std::cout << " Number of active cells: " - << triangulation.n_active_cells() - << std::endl; - - setup_system (); - - std::cout << " Number of degrees of freedom: " - << dof_handler.n_dofs() - << std::endl; - - assemble_system (); - solve (); - output_results (cycle); - } + // The unfortunate thing is that if + // the discrete solution is constant, + // then the error indicators computed + // by the KellyErrorEstimator + // class are zero for each cell as + // well, and the call to + // refine_and_coarsen_fixed_number + // on the triangulation object + // will not flag any cells for + // refinement (why should it if the + // indicated error is zero for each + // cell?). The grid in the next + // iteration will therefore consist + // of four cells only as well, and + // the same problem occurs again. + // + // The conclusion needs to be: while + // of course we will not choose the + // initial grid to be well-suited for + // the accurate solution of the + // problem, we must at least choose + // it such that it has the chance to + // capture the most striking features + // of the solution. In this case, it + // needs to be able to see the right + // hand side. Thus, we refine twice + // globally. (Note that the + // refine_global function is not + // part of the GridRefinement + // class in which + // refine_and_coarsen_fixed_number + // is declared, for example. The + // reason is first that it is not an + // algorithm that computed refinement + // flags from indicators, but more + // importantly that it actually + // performs the refinement, in + // contrast to the functions in + // GridRefinement that only flag + // cells without actually refining + // the grid.) + template + void ElasticProblem::run () + { + for (unsigned int cycle=0; cycle<8; ++cycle) + { + std::cout << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (2); + } + else + refine_grid (); + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + setup_system (); + + std::cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + assemble_system (); + solve (); + output_results (cycle); + } + } } // @sect3{The main function} - // The main function is again exactly - // like in step-6 (apart from the - // changed class names, of course). + // After closing the Step8 + // namespace in the last line above, the + // following is the main function of the + // program and is again exactly like in + // step-6 (apart from the changed class + // names, of course). int main () { try { - deallog.depth_console (0); + dealii::deallog.depth_console (0); - ElasticProblem<2> elastic_problem_2d; + Step8::ElasticProblem<2> elastic_problem_2d; elastic_problem_2d.run (); } catch (std::exception &exc) -- 2.39.5