From 7be3be25153a0cb12ba6997f71ac3580029be8f3 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 8 Aug 2007 04:41:10 +0000 Subject: [PATCH] Generate more documentation. git-svn-id: https://svn.dealii.org/trunk@14920 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-27/doc/intro.dox | 22 ++ deal.II/examples/step-27/doc/results.dox | 31 +++ deal.II/examples/step-27/step-27.cc | 244 +++++++++++++++++++---- 3 files changed, 260 insertions(+), 37 deletions(-) diff --git a/deal.II/examples/step-27/doc/intro.dox b/deal.II/examples/step-27/doc/intro.dox index 7ea6fbd59c..2934903d98 100644 --- a/deal.II/examples/step-27/doc/intro.dox +++ b/deal.II/examples/step-27/doc/intro.dox @@ -680,3 +680,25 @@ Using this technique of eliminating constrained nodes already when transferring local contributions into the global objects, we avoid the problem of having to go back later and change these objects. Timing these operations shows that this makes the overall algorithms faster. + + +

The test case

+ +The test case we will solve with this program is a re-take of the one we +already look at in @ref step_14 "step-14": we solve the Laplace equation +@[ + -\Delta u = f +@] +in 2d, with $f=(x+1)(y+1)$, and with zero Dirichlet boundary values for +$u$. We do so on the domain $[-1,1]^2\backslash[-\frac 12,\frac 12]^2$, i.e. a +square with a square hole in the middle. + +The difference to @ref step_14 "step-14" is of course that we use $hp$ finite +elements for the solution. The testcase is of interest because it has +re-entrant corners in the corners of the hole, at which the solution has +singularities. We therefore expect that the solution will be smooth in the +interior of the domain, and rough in the vicinity of the singularities. The +hope is that our refinement and smoothness indicators will be able to see this +behavior and refine the mesh close to the singularities, while the polynomial +degree is increased away from it. As we will see in the results section, this +is indeed the case. diff --git a/deal.II/examples/step-27/doc/results.dox b/deal.II/examples/step-27/doc/results.dox index e67fccc1c1..957c48b172 100644 --- a/deal.II/examples/step-27/doc/results.dox +++ b/deal.II/examples/step-27/doc/results.dox @@ -1,2 +1,33 @@

Results

+ + +
+examples/step-27> make run
+============================ Running step-27
+Cycle 0:
+   Number of active cells:       768
+   Number of degrees of freedom: 3264
+   Number of constraints       : 0
+Cycle 1:
+   Number of active cells:       996
+   Number of degrees of freedom: 5327
+   Number of constraints       : 494
+Cycle 2:
+   Number of active cells:       1335
+   Number of degrees of freedom: 8947
+   Number of constraints       : 1496
+Cycle 3:
+   Number of active cells:       1626
+   Number of degrees of freedom: 12994
+   Number of constraints       : 2707
+Cycle 4:
+   Number of active cells:       1911
+   Number of degrees of freedom: 17988
+   Number of constraints       : 4112
+Cycle 5:
+   Number of active cells:       2577
+   Number of degrees of freedom: 26936
+   Number of constraints       : 6196
+
+ diff --git a/deal.II/examples/step-27/step-27.cc b/deal.II/examples/step-27/step-27.cc index bb6052088d..e99f14542a 100644 --- a/deal.II/examples/step-27/step-27.cc +++ b/deal.II/examples/step-27/step-27.cc @@ -274,6 +274,33 @@ void LaplaceProblem::setup_system () // global objects. Both of these are // explained in detail in the introduction of // this program. + // + // One other slight complication is the fact + // that because we use different polynomial + // degrees on different cells, the matrices + // and vectors holding local contributions do + // not have the same size on all cells. At + // the beginning of the loop over all cells, + // we therefore each time have to resize them + // to the correct size (given by + // dofs_per_cell). Because these + // classes are implement in such a way that + // reducing the size of a matrix or vector + // does not release the currently allocated + // memory (unless the new size is zero), the + // process of resizing at the beginning of + // the loop will only require re-allocation + // of memory during the first few + // iterations. Once we have found in a cell + // with the maximal finite element degree, no + // more re-allocations will happen because + // all subsequent reinit calls + // will only set the size to something that + // fits the currently allocated memory. This + // is important since allocating memory is + // expensive, and doing so every time we + // visit a new cell would take significant + // compute time. template void LaplaceProblem::assemble_system () { @@ -284,18 +311,22 @@ void LaplaceProblem::assemble_system () const RightHandSide rhs_function; + FullMatrix cell_matrix; + Vector cell_rhs; + + std::vector local_dof_indices; + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) { const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - Vector cell_rhs (dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); + cell_matrix.reinit (dofs_per_cell, dofs_per_cell); cell_matrix = 0; + + cell_rhs.reinit (dofs_per_cell); cell_rhs = 0; hp_fe_values.reinit (cell); @@ -321,6 +352,7 @@ void LaplaceProblem::assemble_system () fe_values.JxW(q_point)); } + local_dof_indices.resize (dofs_per_cell); cell->get_dof_indices (local_dof_indices); hanging_node_constraints @@ -380,11 +412,32 @@ void LaplaceProblem::solve () // @sect4{LaplaceProblem::postprocess} + // After solving the linear system, we will + // want to postprocess the solution. Here, + // all we do is to estimate the error, + // estimate the local smoothness of the + // solution as described in the introduction, + // then write graphical output, and finally + // refine the mesh in both $h$ and $p$ + // according to the indicators computed + // before. We do all this in the same + // function because we want the estimated + // error and smoothness indicators not only + // for refinement, but also include them in + // the graphical output. template void LaplaceProblem::postprocess (const unsigned int cycle) { - Assert (cycle < 10, ExcNotImplemented()); - + // Let us start with computing estimated + // error and smoothness indicators, which + // each are one number for each active cell + // of our triangulation. For the error + // indicator, we use the + // KellyErrorEstimator class as + // always. Estimating the smoothness is + // done in the respective function of this + // class; that function is discussed + // further down below: Vector estimated_error_per_cell (triangulation.n_active_cells()); KellyErrorEstimator::estimate (dof_handler, face_quadrature_collection, @@ -394,39 +447,115 @@ void LaplaceProblem::postprocess (const unsigned int cycle) Vector smoothness_indicators (triangulation.n_active_cells()); estimate_smoothness (smoothness_indicators); - + + // Next we want to generate graphical + // output. In addition to the two estimated + // quantities derived above, we would also + // like to output the polynomial degree of + // the finite elements used on each of the + // elements on the mesh. + // + // The way to do that requires that we loop + // over all cells and poll the active + // finite element index of them using + // cell-@>active_fe_index(). We + // then use the result of this operation + // and query the finite element collection + // for the finite element with that index, + // and finally determine the polynomial + // degree of that element. The result we + // put into a vector with one element per + // cell. The DataOut class requires this to + // be a vector of float or + // double, even though our + // values are all integers, so that it what + // we use: { - Vector fe_indices (triangulation.n_active_cells()); + Vector fe_degrees (triangulation.n_active_cells()); { typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (unsigned int index=0; cell!=endc; ++cell, ++index) - fe_indices(index) = cell->active_fe_index(); + fe_degrees(index) + = fe_collection[cell->active_fe_index()].degree; } - - const std::string filename = "solution-" + - Utilities::int_to_string (cycle, 2) + - ".gmv"; + + // With now all data vectors available -- + // solution, estimated errors and + // smoothness indicators, and finite + // element degrees --, we create a + // DataOut object for graphical output + // and attach all data. Note that the + // DataOut class has a second template + // argument (which defaults to + // DoFHandler@, which is why we + // have never seen it in previous + // tutorial programs) that indicates the + // type of DoF handler to be used. Here, + // we have to use the hp::DoFHandler + // class: DataOut > data_out; data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, "solution"); + data_out.add_data_vector (estimated_error_per_cell, "error"); data_out.add_data_vector (smoothness_indicators, "smoothness"); - data_out.add_data_vector (fe_indices, "fe_index"); + data_out.add_data_vector (fe_degrees, "fe_degree"); data_out.build_patches (); - + + // The final step in generating output is + // to determine a file name, open the + // file, and write the data into it: + const std::string filename = "solution-" + + Utilities::int_to_string (cycle, 2) + + ".gmv"; std::ofstream output (filename.c_str()); data_out.write_gmv (output); } - + + // After this, we would like to actually + // refine the mesh, in both $h$ and + // $p$. The way we are going to do this is + // as follows: first, we use the estimated + // error to flag those cells for refinement + // that have the largest error. This is + // what we have always done: { GridRefinement::refine_and_coarsen_fixed_number (triangulation, estimated_error_per_cell, 0.3, 0.03); - float max_smoothness = 0, - min_smoothness = smoothness_indicators.linfty_norm(); + // Next we would like to figure out which + // of the cells that have been flagged + // for refinement should actually have + // $p$ increased instead of $h$ + // decreased. The strategy we choose here + // is that we look at the smoothness + // indicators of those cells that are + // flagged for refinement, and increase + // $p$ for those with a smoothness larger + // than a certain threshold. For this, we + // first have to determine the maximal + // and minimal values of the smoothness + // indicators of all flagged cells, which + // we do using a loop over all cells and + // comparing current minimal and maximal + // values. (We start with the minimal and + // maximal values of all cells, a + // range within which the minimal and + // maximal values on cells flagged for + // refinement must surely lie.) Absent + // any better strategies, we will then + // set the threshold above which will + // increase $p$ instead of reducing $h$ + // as the mean value between minimal and + // maximal smoothness indicators on cells + // flagged for refinement: + float max_smoothness = *std::min_element (smoothness_indicators.begin(), + smoothness_indicators.end()), + min_smoothness = *std::max_element (smoothness_indicators.begin(), + smoothness_indicators.end()); { typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), @@ -440,8 +569,22 @@ void LaplaceProblem::postprocess (const unsigned int cycle) smoothness_indicators(index)); } } - - const float cutoff_smoothness = (max_smoothness + min_smoothness) / 2; + const float threshold_smoothness = (max_smoothness + min_smoothness) / 2; + + // With this, we can go back, loop over + // all cells again, and for those cells + // for which (i) the refinement flag is + // set, (ii) the smoothness indicator is + // larger than the threshold, and (iii) + // we still have a finite element with a + // polynomial degree higher than the + // current one in the finite element + // collection, we then increase the + // polynomial degree and in return remove + // the flag indicating that the cell + // should undergo bisection. For all + // other cells, the refinement flags + // remain untouched: { typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), @@ -449,23 +592,39 @@ void LaplaceProblem::postprocess (const unsigned int cycle) for (unsigned int index=0; cell!=endc; ++cell, ++index) if (cell->refine_flag_set() && - (smoothness_indicators(index) > cutoff_smoothness) + (smoothness_indicators(index) > threshold_smoothness) && - !(cell->active_fe_index() == fe_collection.size() - 1)) + (cell->active_fe_index()+1 < fe_collection.size())) { cell->clear_refine_flag(); - // REMOVE redundant std::min - cell->set_active_fe_index (std::min (cell->active_fe_index() + 1, - fe_collection.size() - 1)); + cell->set_active_fe_index (cell->active_fe_index() + 1); } } - + + // At the end of this procedure, we then + // refine the mesh. During this process, + // children of cells undergoing bisection + // inherit their mother cell's finite + // element index: triangulation.execute_coarsening_and_refinement (); } } // @sect4{LaplaceProblem::create_coarse_grid} + + // The following function is used when + // creating the initial grid. It is a + // specialization for the 2d case, i.e. a + // corresponding function needs to be + // implemented if the program is run in + // anything other then 2d. The function is + // actually stolen from step-14 and generates + // the same mesh used already there, i.e. the + // square domain with the square hole in the + // middle. The meaning of the different parts + // of this function are explained in the + // documentation of step-14: template <> void LaplaceProblem<2>::create_coarse_grid () { @@ -540,34 +699,45 @@ void LaplaceProblem<2>::create_coarse_grid () // @sect4{LaplaceProblem::run} + + // This function implements the logic of the + // program, as did the respective function in + // most of the previous programs already, see + // for example step-6. + // + // Basically, it contains the adaptive loop: + // in the first iteration create a coarse + // grid, and then set up the linear system, + // assemble it, solve, and postprocess the + // solution including mesh refinement. Then + // start over again. In the meantime, also + // output some information for those staring + // at the screen trying to figure out what + // the program does: template void LaplaceProblem::run () { - for (unsigned int cycle=0; cycle<5; ++cycle) + for (unsigned int cycle=0; cycle<6; ++cycle) { std::cout << "Cycle " << cycle << ':' << std::endl; if (cycle == 0) create_coarse_grid (); - std::cout << " Number of active cells: " - << triangulation.n_active_cells() - << std::endl; - setup_system (); - std::cout << " Number of degrees of freedom: " + std::cout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl + << " Number of degrees of freedom: " << dof_handler.n_dofs() - << std::endl; - std::cout << " Number of constraints : " + << std::endl + << " Number of constraints : " << hanging_node_constraints.n_constraints() << std::endl; assemble_system (); - - solve (); - postprocess (cycle); } } -- 2.39.5