From 8d5b7b9f0d757776bf8a239de95fb6e0a9ce9de9 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 8 Feb 2006 20:44:53 +0000 Subject: [PATCH] Read over the rest. git-svn-id: https://svn.dealii.org/trunk@12270 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-8/step-8.cc | 405 +++++++++++++++--------------- 1 file changed, 198 insertions(+), 207 deletions(-) diff --git a/deal.II/examples/step-8/step-8.cc b/deal.II/examples/step-8/step-8.cc index 1601d595d7..3b5fdf0a56 100644 --- a/deal.II/examples/step-8/step-8.cc +++ b/deal.II/examples/step-8/step-8.cc @@ -509,39 +509,37 @@ void ElasticProblem::setup_system () // 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 () { - // First thing: the quadrature - // formula does not need - // modification since we still deal - // with bilinear functions. QGauss quadrature_formula(2); - // Also, the ``FEValues'' objects - // takes care of everything for us - // (or better: it does not really - // so; as in the comment in the - // function setting up the system, - // here as well the ``FEValues'' - // object computes the same data on - // each cell, but it has some - // functionality to access data - // stored inside the finite element - // where they are precomputed upon - // construction). + FEValues fe_values (fe, quadrature_formula, update_values | update_gradients | update_q_points | update_JxW_values); - // 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, but this is not - // something we need to care about. const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.n_quadrature_points; @@ -604,50 +602,65 @@ void ElasticProblem::assemble_system () // Next we get the values of // the coefficients at the - // quadrature points: + // 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 - // and will not comment much on - // this. + // 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; i'', of - // which the first element - // is ``comp(i)'' and the - // second is the value - // ``base(i)'' also noted - // in the text. You will - // rather seldom need to - // access this second - // value, but the first is - // important when using - // vector valued elements. for (unsigned int j=0; j::assemble_system () for (unsigned int q_point=0; q_point::assemble_system () ) * fe_values.JxW(q_point); - }; - }; - }; + } + } + } // Assembling the right hand // side is also just as // discussed in the - // introduction. We will - // therefore not discuss it - // further. - right_hand_side.vector_value_list (fe_values.get_quadrature_points(), - rhs_values); + // introduction: for (unsigned int i=0; i::assemble_system () cell_rhs(i) += fe_values.shape_value(i,q_point) * rhs_values[q_point](component_i) * fe_values.JxW(q_point); - }; + } // The transfer from local // degrees of freedom into the @@ -780,7 +787,12 @@ void ElasticProblem::assemble_system () // on the equation under // consideration, and is thus // the same as in all previous - // examples. + // examples. The same holds for + // the elimination of hanging + // nodes from the matrix and + // right hand side, once we are + // done with assembling the + // entire linear system: cell->get_dof_indices (local_dof_indices); for (unsigned int i=0; i::assemble_system () cell_matrix(i,j)); system_rhs(local_dof_indices[i]) += cell_rhs(i); - }; - }; + } + } hanging_node_constraints.condense (system_matrix); hanging_node_constraints.condense (system_rhs); @@ -800,7 +812,7 @@ void ElasticProblem::assemble_system () // boundary values needs a small // modification: since the solution // function is vector-valued, so - // needs to be the boundary + // need to be the boundary // values. The ``ZeroFunction'' // constructor accepts a parameter // that tells it that it shall @@ -829,14 +841,16 @@ void ElasticProblem::assemble_system () + // @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 is. Therefore, we need - // not change anything. + // the system indeed is. Therefore, + // we need not change anything. template void ElasticProblem::solve () { @@ -853,6 +867,7 @@ void ElasticProblem::solve () } + // @sect4{ElasticProblem::refine_grid} // The function that does the // refinement of the grid is the same @@ -862,8 +877,8 @@ void ElasticProblem::solve () // that the error estimator by // default adds up the estimated // obtained from all components of - // the finite element solution, that - // is it uses the displacement in all + // 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 @@ -873,7 +888,11 @@ void ElasticProblem::solve () // and do not consider the // displacements in all other // directions for the error - // indicators. + // indicators. However, for the + // current problem, it seems + // appropriate to consider all + // displacement components with equal + // weight. template void ElasticProblem::refine_grid () { @@ -894,6 +913,8 @@ void ElasticProblem::refine_grid () } + // @sect4{ElasticProblem::output_results} + // The output happens mostly as has // been shown in previous examples // already. The only difference is @@ -944,6 +965,20 @@ void ElasticProblem::output_results (const unsigned int cycle) const // 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) { @@ -959,26 +994,9 @@ void ElasticProblem::output_results (const unsigned int cycle) const solution_names.push_back ("y_displacement"); solution_names.push_back ("z_displacement"); break; - // 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 throws 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 this statement: default: - Assert (false, ExcInternalError()); - }; + Assert (false, ExcNotImplemented()); + } // After setting up the names for // the different components of the @@ -1005,6 +1023,84 @@ void ElasticProblem::output_results (const unsigned int cycle) const + // @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 () { @@ -1014,112 +1110,7 @@ void ElasticProblem::run () if (cycle == 0) { - // As in previous examples, - // we use the unit square - // (or cube) as domain. GridGenerator::hyper_cube (triangulation, -1, 1); - // This time, we have to - // refine the coarse grid - // twice before we first - // solve on it. 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.) triangulation.refine_global (2); } else @@ -1138,7 +1129,7 @@ void ElasticProblem::run () assemble_system (); solve (); output_results (cycle); - }; + } } // @sect3{The ``main'' function} @@ -1178,7 +1169,7 @@ int main () << "----------------------------------------------------" << std::endl; return 1; - }; + } return 0; } -- 2.39.5