From: Timo Heister Date: Fri, 25 Aug 2017 18:56:57 +0000 (-0400) Subject: add mesh_loop() test: step-11 X-Git-Tag: v9.0.0-rc1~1166^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4aa95116c9d119a674f1274cadc54d17150e362d;p=dealii.git add mesh_loop() test: step-11 --- diff --git a/tests/meshworker/step-11-mesh_loop.cc b/tests/meshworker/step-11-mesh_loop.cc new file mode 100644 index 0000000000..0a4dd4e47a --- /dev/null +++ b/tests/meshworker/step-11-mesh_loop.cc @@ -0,0 +1,404 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// step-11 but rewritten using mesh_loop() + + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Just this one is new: it declares a class +// DynamicSparsityPattern, which we will use and explain +// further down below. +#include + +// We will make use of the std::find algorithm of the C++ standard library, so +// we have to include the following file for its declaration: +#include +#include +#include +#include + +// The last step is as in all previous programs: +namespace Step11 +{ + using namespace dealii; + + // Then we declare a class which represents the solution of a Laplace + // problem. As this example program is based on step-5, the class looks + // rather the same, with the sole structural difference that the functions + // assemble_system now calls solve itself, and is + // thus called assemble_and_solve, and that the output function + // was dropped since the solution function is so boring that it is not worth + // being viewed. + // + // The only other noteworthy change is that the constructor takes a value + // representing the polynomial degree of the mapping to be used later on, + // and that it has another member variable representing exactly this + // mapping. In general, this variable will occur in real applications at the + // same places where the finite element is declared or used. + template + class LaplaceProblem + { + public: + LaplaceProblem (const unsigned int mapping_degree); + void run (); + + private: + void setup_system (); + void assemble_and_solve (); + void solve (); + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + MappingQ mapping; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + ConstraintMatrix constraints; + + Vector solution; + Vector system_rhs; + + TableHandler output_table; + }; + + + + // Construct such an object, by initializing the variables. Here, we use + // linear finite elements (the argument to the fe variable + // denotes the polynomial degree), and mappings of given order. Print to + // screen what we are about to do. + template + LaplaceProblem::LaplaceProblem (const unsigned int mapping_degree) : + fe (1), + dof_handler (triangulation), + mapping (mapping_degree) + { + deallog << "Using mapping with degree " << mapping_degree << ":\n" + << "============================" + << std::endl; + } + + + + // The first task is to set up the variables for this problem. This includes + // generating a valid DoFHandler object, as well as the + // sparsity patterns for the matrix, and the object representing the + // constraints that the mean value of the degrees of freedom on the boundary + // be zero. + template + void LaplaceProblem::setup_system () + { + dof_handler.distribute_dofs (fe); + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + + std::vector boundary_dofs (dof_handler.n_dofs(), false); + DoFTools::extract_boundary_dofs (dof_handler, + ComponentMask(), + boundary_dofs); + + const unsigned int first_boundary_dof + = std::distance (boundary_dofs.begin(), + std::find (boundary_dofs.begin(), + boundary_dofs.end(), + true)); + + // Then generate a constraints object with just this one constraint. First + // clear all previous content (which might reside there from the previous + // computation on a once coarser grid), then add this one line + // constraining the first_boundary_dof to the sum of other + // boundary DoFs each with weight -1. Finally, close the constraints + // object, i.e. do some internal bookkeeping on it for faster processing + // of what is to come later: + constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.add_line (first_boundary_dof); + for (unsigned int i=first_boundary_dof+1; i + struct ScratchData + { + ScratchData (const Mapping &mapping, + const FiniteElement &fe, + const unsigned int quadrature_degree) + : + fe_values (mapping, + fe, + QGauss(quadrature_degree), + update_values | update_gradients | + update_quadrature_points | update_JxW_values), + fe_face_values (mapping, + fe, + QGauss(quadrature_degree), + update_values | update_quadrature_points | + update_JxW_values | update_normal_vectors) + {} + + ScratchData (const ScratchData &scratch_data) + : + fe_values (scratch_data.fe_values.get_mapping(), + scratch_data.fe_values.get_fe(), + scratch_data.fe_values.get_quadrature(), + update_values | update_gradients | + update_quadrature_points | update_JxW_values), + fe_face_values (scratch_data.fe_face_values.get_mapping(), + scratch_data.fe_face_values.get_fe(), + scratch_data.fe_face_values.get_quadrature(), + update_values | update_quadrature_points | + update_JxW_values | update_normal_vectors) + {} + + FEValues fe_values; + FEFaceValues fe_face_values; + + }; + struct CopyData + { + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + }; + + // The next function then assembles the linear system of equations, solves + // it, and evaluates the solution. + template + void LaplaceProblem::assemble_and_solve () + { + + + typedef decltype(dof_handler.begin_active()) Iterator; + + auto cell_worker = [] (const Iterator &cell, ScratchData &scratch_data, CopyData ©_data) + { + const unsigned int dofs_per_cell = scratch_data.fe_values.get_fe().dofs_per_cell; + const unsigned int n_q_points = scratch_data.fe_values.get_quadrature().size(); + + copy_data.cell_matrix.reinit (dofs_per_cell, dofs_per_cell); + copy_data.cell_rhs.reinit (dofs_per_cell); + + copy_data.local_dof_indices.resize(dofs_per_cell); + cell->get_dof_indices (copy_data.local_dof_indices); + + scratch_data.fe_values.reinit (cell); + + std::vector rhs_values (n_q_points); + ConstantFunction right_hand_side (-2.0); + right_hand_side.value_list (scratch_data.fe_values.get_quadrature_points(), + rhs_values); + + // ... and assemble the local contributions to the system matrix and + // right hand side as also discussed above: + for (unsigned int q_point=0; q_point &scratch_data, CopyData ©_data) + { + const unsigned int dofs_per_cell = scratch_data.fe_values.get_fe().dofs_per_cell; + const unsigned int n_face_q_points = scratch_data.fe_face_values.get_quadrature().size(); + + std::vector face_boundary_values (n_face_q_points); + ConstantFunction boundary_values (1.0); + + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->at_boundary()) + { + scratch_data.fe_face_values.reinit (cell, face); + boundary_values.value_list (scratch_data.fe_face_values.get_quadrature_points(), + face_boundary_values); + + for (unsigned int q_point=0; q_point(std::ceil(1.*(mapping.get_degree()+1)/2)), + 2U); + + MeshWorker::mesh_loop(dof_handler.begin_active(), + dof_handler.end(), + cell_worker, + copier, + ScratchData(mapping, fe, gauss_degree), + CopyData(), + MeshWorker::assemble_own_cells | MeshWorker::assemble_boundary_faces | MeshWorker::cells_first, + boundary_worker); + + + solve (); + + + // Finally, evaluate what we got as solution. As stated in the + // introduction, we are interested in the H1 semi-norm of the + // solution. Here, as well, we have a function in the library that does + // this, although in a slightly non-obvious way: the + // VectorTools::integrate_difference function integrates the + // norm of the difference between a finite element function and a + // continuous function. If we therefore want the norm of a finite element + // field, we just put the continuous function to zero. Note that this + // function, just as so many other ones in the library as well, has at + // least two versions, one which takes a mapping as argument (which we + // make us of here), and the one which we have used in previous examples + // which implicitly uses MappingQ1. Also note that we take a + // quadrature formula of one degree higher, in order to avoid + // superconvergence effects where the solution happens to be especially + // close to the exact solution at certain points (we don't know whether + // this might be the case here, but there are cases known of this, and we + // just want to make sure): + Vector norm_per_cell (triangulation.n_active_cells()); + VectorTools::integrate_difference (mapping, dof_handler, + solution, + ZeroFunction(), + norm_per_cell, + QGauss(gauss_degree+1), + VectorTools::H1_seminorm); + // Then, the function just called returns its results as a vector of + // values each of which denotes the norm on one cell. To get the global + // norm, we do the following: + const double norm = VectorTools::compute_global_error(triangulation, + norm_per_cell, + VectorTools::H1_seminorm); + + // Last task -- generate output: + output_table.add_value ("cells", triangulation.n_active_cells()); + output_table.add_value ("|u|_1", norm); + output_table.add_value ("error", std::fabs(norm-std::sqrt(3.14159265358/2))); + } + + + + // The following function solving the linear system of equations is copied + // from step-5 and is explained there in some detail: + template + void LaplaceProblem::solve () + { + SolverControl solver_control (1000, 1e-12, false, false); + SolverCG<> cg (solver_control); + + PreconditionSSOR<> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + + constraints.distribute (solution); + } + + + + // Finally the main function controlling the different steps to be + // performed. + template + void LaplaceProblem::run () + { + GridGenerator::hyper_ball (triangulation); + static const SphericalManifold boundary; + triangulation.set_all_manifold_ids_on_boundary(0); + triangulation.set_manifold (0, boundary); + + for (unsigned int cycle=0; cycle<5; ++cycle) + { + if (cycle>0) + triangulation.refine_global(1); + + setup_system (); + assemble_and_solve (); + }; + + // After all the data is generated, write a table of results to the + // screen: + output_table.set_precision("|u|_1", 6); + output_table.set_precision("error", 6); + output_table.write_text (deallog.get_file_stream()); + deallog << std::endl; + } +} + +int main () +{ + initlog(); + + for (unsigned int mapping_degree=1; mapping_degree<=3; ++mapping_degree) + Step11::LaplaceProblem<2>(mapping_degree).run (); + + return 0; +} diff --git a/tests/meshworker/step-11-mesh_loop.output b/tests/meshworker/step-11-mesh_loop.output new file mode 100644 index 0000000000..b1c061c30e --- /dev/null +++ b/tests/meshworker/step-11-mesh_loop.output @@ -0,0 +1,28 @@ + +DEAL::Using mapping with degree 1: +============================ +cells |u|_1 error +5 0.680402 0.572912 +20 1.088141 0.165173 +80 1.210142 0.043172 +320 1.242375 0.010939 +1280 1.250569 0.002745 +DEAL:: +DEAL::Using mapping with degree 2: +============================ +cells |u|_1 error +5 1.050963 0.202351 +20 1.203014 0.050300 +80 1.241120 0.012194 +320 1.250319 0.002996 +1280 1.252572 0.000742 +DEAL:: +DEAL::Using mapping with degree 3: +============================ +cells |u|_1 error +5 1.069381 0.183934 +20 1.206606 0.046708 +80 1.241611 0.011703 +320 1.250382 0.002932 +1280 1.252580 0.000734 +DEAL::