From: kanschat Date: Wed, 21 Aug 2013 19:07:25 +0000 (+0000) Subject: test DG with constraints and remove text X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7d149344ccf935c6aed86fd427eb541895056bfd;p=dealii-svn.git test DG with constraints and remove text git-svn-id: https://svn.dealii.org/trunk@30380 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/multigrid/step-16.cc b/tests/multigrid/step-16.cc index 1d881e62b5..01bd6949bb 100644 --- a/tests/multigrid/step-16.cc +++ b/tests/multigrid/step-16.cc @@ -23,25 +23,6 @@ #include "../tests.h" #include -// As discussed in the introduction, most of -// this program is copied almost verbatim -// from step-6, which itself is only a slight -// modification of step-5. Consequently, a -// significant part of this program is not -// new if you've read all the material up to -// step-6, and we won't comment on that part -// of the functionality that is -// unchanged. Rather, we will focus on those -// aspects of the program that have to do -// with the multigrid functionality which -// forms the new aspect of this tutorial -// program. - -// @sect3{Include files} - -// Again, the first few include files -// are already known, so we won't -// comment on them: #include #include #include @@ -71,19 +52,6 @@ #include #include -// These, now, are the include necessary for -// the multi-level methods. The first two -// declare classes that allow us to enumerate -// degrees of freedom not only on the finest -// mesh level, but also on intermediate -// levels (that's what the MGDoFHandler class -// does) as well as allow to access this -// information (iterators and accessors over -// these cells). -// -// The rest of the include files deals with -// the mechanics of multigrid as a linear -// operator (solver or preconditioner). #include #include #include @@ -92,24 +60,11 @@ #include #include -// This is C++: #include #include -// The last step is as in all -// previous programs: using namespace dealii; - -// @sect3{The LaplaceProblem class template} - -// This main class is basically the same -// class as in step-6. As far as member -// functions is concerned, the only addition -// is the assemble_multigrid -// function that assembles the matrices that -// correspond to the discrete operators on -// intermediate levels: template class LaplaceProblem { @@ -140,37 +95,6 @@ private: const unsigned int degree; - // The following three objects are the - // only additional member variables, - // compared to step-6. They represent the - // operators that act on individual - // levels of the multilevel hierarchy, - // rather than on the finest mesh as do - // the objects above. - // - // To facilitate having objects on each - // level of a multilevel hierarchy, - // deal.II has the MGLevelObject class - // template that provides storage for - // objects on each level. What we need - // here are matrices on each level, which - // implies that we also need sparsity - // patterns on each level. As outlined in - // the @ref mg_paper, the operators - // (matrices) that we need are actually - // twofold: one on the interior of each - // level, and one at the interface - // between each level and that part of - // the domain where the mesh is - // coarser. In fact, we will need the - // latter in two versions: for the - // direction from coarse to fine mesh and - // from fine to coarse. Fortunately, - // however, we here have a self-adjoint - // problem for which one of these is the - // transpose of the other, and so we only - // have to build one; we choose the one - // from coarse to fine. MGLevelObject mg_sparsity_patterns; MGLevelObject > mg_matrices; MGLevelObject > mg_interface_matrices; @@ -178,13 +102,6 @@ private: }; - -// @sect3{Nonconstant coefficients} - -// The implementation of nonconstant -// coefficients is copied verbatim -// from step-5 and step-6: - template class Coefficient : public Function { @@ -231,33 +148,6 @@ void Coefficient::value_list (const std::vector > &points, } -// @sect3{The LaplaceProblem class implementation} - -// @sect4{LaplaceProblem::LaplaceProblem} - -// The constructor is left mostly -// unchanged. We take the polynomial degree -// of the finite elements to be used as a -// constructor argument and store it in a -// member variable. -// -// By convention, all adaptively refined -// triangulations in deal.II never change by -// more than one level across a face between -// cells. For our multigrid algorithms, -// however, we need a slightly stricter -// guarantee, namely that the mesh also does -// not change by more than refinement level -// across vertices that might connect two -// cells. In other words, we must prevent the -// following situation: -// -// @image html limit_level_difference_at_vertices.png "" -// -// This is achieved by passing the -// Triangulation::limit_level_difference_at_vertices -// flag to the constructor of the -// triangulation class. template LaplaceProblem::LaplaceProblem (const unsigned int degree) : @@ -269,22 +159,10 @@ LaplaceProblem::LaplaceProblem (const unsigned int degree) {} - -// @sect4{LaplaceProblem::setup_system} - -// The following function extends what the -// corresponding one in step-6 did. The top -// part, apart from the additional output, -// does the same: template void LaplaceProblem::setup_system () { mg_dof_handler.distribute_dofs (fe); - - // Here we output not only the - // degrees of freedom on the finest - // level, but also in the - // multilevel structure deallog << "Number of degrees of freedom: " << mg_dof_handler.n_dofs(); @@ -303,28 +181,6 @@ void LaplaceProblem::setup_system () solution.reinit (mg_dof_handler.n_dofs()); system_rhs.reinit (mg_dof_handler.n_dofs()); - // But it starts to be a wee bit different - // here, although this still doesn't have - // anything to do with multigrid - // methods. step-6 took care of boundary - // values and hanging nodes in a separate - // step after assembling the global matrix - // from local contributions. This works, - // but the same can be done in a slightly - // simpler way if we already take care of - // these constraints at the time of copying - // local contributions into the global - // matrix. To this end, we here do not just - // compute the constraints do to hanging - // nodes, but also due to zero boundary - // conditions. Both kinds of constraints - // can be put into the same object - // (constraints), and we will - // use this set of constraints later on to - // help us copy local contributions - // correctly into the global linear system - // right away, without the need for a later - // clean-up stage: constraints.clear (); hanging_node_constraints.clear (); DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints); @@ -345,21 +201,6 @@ void LaplaceProblem::setup_system () mg_constrained_dofs.clear(); mg_constrained_dofs.initialize(mg_dof_handler, dirichlet_boundary); - // Now for the things that concern the - // multigrid data structures. First, we - // resize the multi-level objects to hold - // matrices and sparsity patterns for every - // level. The coarse level is zero (this is - // mandatory right now but may change in a - // future revision). Note that these - // functions take a complete, inclusive - // range here (not a starting index and - // size), so the finest level is - // n_levels-1. We first have - // to resize the container holding the - // SparseMatrix classes, since they have to - // release their SparsityPattern before the - // can be destroyed upon resizing. const unsigned int n_levels = triangulation.n_levels(); mg_interface_matrices.resize(0, n_levels-1); @@ -368,32 +209,6 @@ void LaplaceProblem::setup_system () mg_matrices.clear (); mg_sparsity_patterns.resize(0, n_levels-1); - // Now, we have to provide a matrix on each - // level. To this end, we first use the - // MGTools::make_sparsity_pattern function - // to first generate a preliminary - // compressed sparsity pattern on each - // level (see the @ref Sparsity module for - // more information on this topic) and then - // copy it over to the one we really - // want. The next step is to initialize - // both kinds of level matrices with these - // sparsity patterns. - // - // It may be worth pointing out that the - // interface matrices only have entries for - // degrees of freedom that sit at or next - // to the interface between coarser and - // finer levels of the mesh. They are - // therefore even sparser than the matrices - // on the individual levels of our - // multigrid hierarchy. If we were more - // concerned about memory usage (and - // possibly the speed with which we can - // multiply with these matrices), we should - // use separate and different sparsity - // patterns for these two kinds of - // matrices. for (unsigned int level=0; level::setup_system () } -// @sect4{LaplaceProblem::assemble_system} - -// The following function assembles the -// linear system on the finesh level of the -// mesh. It is almost exactly the same as in -// step-6, with the exception that we don't -// eliminate hanging nodes and boundary -// values after assembling, but while copying -// local contributions into the global -// matrix. This is not only simpler but also -// more efficient for large problems. template void LaplaceProblem::assemble_system () { @@ -475,23 +279,6 @@ void LaplaceProblem::assemble_system () } -// @sect4{LaplaceProblem::assemble_multigrid} - -// The next function is the one that builds -// the linear operators (matrices) that -// define the multigrid method on each level -// of the mesh. The integration core is the -// same as above, but the loop below will go -// over all existing cells instead of just -// the active ones, and the results must be -// entered into the correct matrix. Note also -// that since we only do multi-level -// preconditioning, no right-hand side needs -// to be assembled here. -// -// Before we go there, however, we have to -// take care of a significant amount of book -// keeping: template void LaplaceProblem::assemble_multigrid () { @@ -511,70 +298,11 @@ void LaplaceProblem::assemble_multigrid () const Coefficient coefficient; std::vector coefficient_values (n_q_points); - // Next a few things that are specific to - // building the multigrid data structures - // (since we only need them in the current - // function, rather than also elsewhere, we - // build them here instead of the - // setup_system - // function). Some of the following may be - // a bit obscure if you're not familiar - // with the algorithm actually implemented - // in deal.II to support multilevel - // algorithms on adaptive meshes; if some - // of the things below seem strange, take a - // look at the @ref mg_paper. - // - // Our first job is to identify those - // degrees of freedom on each level that - // are located on interfaces between - // adaptively refined levels, and those - // that lie on the interface but also on - // the exterior boundary of the domain. As - // in many other parts of the library, we - // do this by using boolean masks, - // i.e. vectors of booleans each element of - // which indicates whether the - // corresponding degree of freedom index is - // an interface DoF or not: std::vector > interface_dofs = mg_constrained_dofs.get_refinement_edge_indices (); std::vector > boundary_interface_dofs = mg_constrained_dofs.get_refinement_edge_boundary_indices (); - - // The indices just identified will later - // be used to impose zero boundary - // conditions for the operator that we will - // apply on each level. On the other hand, - // we also have to impose zero boundary - // conditions on the external boundary of - // each level. So let's identify these - // nodes as well (this time as a set of - // degrees of freedom, rather than a - // boolean mask; the reason for this being - // that we will not need fast tests whether - // a certain degree of freedom is in the - // boundary list, though we will need such - // access for the interface degrees of - // freedom further down below): - - // The third step is to construct - // constraints on all those degrees of - // freedom: their value should be zero - // after each application of the level - // operators. To this end, we construct - // ConstraintMatrix objects for each level, - // and add to each of these constraints for - // each degree of freedom. Due to the way - // the ConstraintMatrix stores its data, - // the function to add a constraint on a - // single degree of freedom and force it to - // be zero is called - // Constraintmatrix::add_line(); doing so - // for several degrees of freedom at once - // can be done using - // Constraintmatrix::add_lines(): std::vector boundary_constraints (triangulation.n_levels()); std::vector boundary_interface_constraints (triangulation.n_levels()); for (unsigned int level=0; level::assemble_multigrid () boundary_interface_constraints[level].close (); } - // Now that we're done with most of our - // preliminaries, let's start the - // integration loop. It looks mostly like - // the loop in - // assemble_system, with two - // exceptions: (i) we don't need a right - // han side, and more significantly (ii) we - // don't just loop over all active cells, - // but in fact all cells, active or - // not. Consequently, the correct iterator - // to use is MGDoFHandler::cell_iterator - // rather than - // MGDoFHandler::active_cell_iterator. Let's - // go about it: typename MGDoFHandler::cell_iterator cell = mg_dof_handler.begin(), endc = mg_dof_handler.end(); @@ -621,35 +335,8 @@ void LaplaceProblem::assemble_multigrid () fe_values.shape_grad(j,q_point) * fe_values.JxW(q_point)); - // The rest of the assembly is again - // slightly different. This starts with - // a gotcha that is easily forgotten: - // The indices of global degrees of - // freedom we want here are the ones - // for current level, not for the - // global matrix. We therefore need the - // function - // MGDoFAccessorLLget_mg_dof_indices, - // not MGDoFAccessor::get_dof_indices - // as used in the assembly of the - // global system: cell->get_mg_dof_indices (local_dof_indices); - // Next, we need to copy local - // contributions into the level - // objects. We can do this in the same - // way as in the global assembly, using - // a constraint object that takes care - // of constrained degrees (which here - // are only boundary nodes, as the - // individual levels have no hanging - // node constraints). Note that the - // boundary_constraints - // object makes sure that the level - // matrices contains no contributions - // from degrees of freedom at the - // interface between cells of different - // refinement level. boundary_constraints[cell->level()] .distribute_local_to_global (cell_matrix, local_dof_indices, diff --git a/tests/multigrid/step-39-02a.cc b/tests/multigrid/step-39-02a.cc new file mode 100644 index 0000000000..6eb5dafa25 --- /dev/null +++ b/tests/multigrid/step-39-02a.cc @@ -0,0 +1,767 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2013 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. +// +// --------------------------------------------------------------------- + + +// Add edge matrices and MGConstraints to make sure they are empty and do not mess things up + +#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 +#include + +#include +#include +#include +#include + +#include +#include + +namespace Step39 +{ + using namespace dealii; + + Functions::SlitSingularityFunction<2> exact_solution; + + + + + template + class MatrixIntegrator : public MeshWorker::LocalIntegrator + { + public: + void cell(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; + }; + + + template + void MatrixIntegrator::cell( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const + { + LocalIntegrators::Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values()); + } + + + template + void MatrixIntegrator::boundary( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const + { + const unsigned int deg = info.fe_values(0).get_fe().tensor_degree(); + LocalIntegrators::Laplace::nitsche_matrix( + dinfo.matrix(0,false).matrix, info.fe_values(0), + LocalIntegrators::Laplace::compute_penalty(dinfo, dinfo, deg, deg)); + } + + template + void MatrixIntegrator::face( + MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const + { + const unsigned int deg = info1.fe_values(0).get_fe().tensor_degree(); + LocalIntegrators::Laplace::ip_matrix( + dinfo1.matrix(0,false).matrix, dinfo1.matrix(0,true).matrix, + dinfo2.matrix(0,true).matrix, dinfo2.matrix(0,false).matrix, + info1.fe_values(0), info2.fe_values(0), + LocalIntegrators::Laplace::compute_penalty(dinfo1, dinfo2, deg, deg)); + } + + template + class RHSIntegrator : public MeshWorker::LocalIntegrator + { + public: + void cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; + }; + + + template + void RHSIntegrator::cell(MeshWorker::DoFInfo &, typename MeshWorker::IntegrationInfo &) const + {} + + + template + void RHSIntegrator::boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + Vector &local_vector = dinfo.vector(0).block(0); + + std::vector boundary_values(fe.n_quadrature_points); + exact_solution.value_list(fe.get_quadrature_points(), boundary_values); + + const unsigned int deg = fe.get_fe().tensor_degree(); + const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure(); + + for (unsigned k=0; k + void RHSIntegrator::face(MeshWorker::DoFInfo &, + MeshWorker::DoFInfo &, + typename MeshWorker::IntegrationInfo &, + typename MeshWorker::IntegrationInfo &) const + {} + + + template + class Estimator : public MeshWorker::LocalIntegrator + { + public: + void cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; + }; + + + template + void Estimator::cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + + const std::vector > &DDuh = info.hessians[0][0]; + for (unsigned k=0; kdiameter() * trace(DDuh[k]); + dinfo.value(0) += t*t * fe.JxW(k); + } + dinfo.value(0) = std::sqrt(dinfo.value(0)); + } + + template + void Estimator::boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + + std::vector boundary_values(fe.n_quadrature_points); + exact_solution.value_list(fe.get_quadrature_points(), boundary_values); + + const std::vector &uh = info.values[0][0]; + + const unsigned int deg = fe.get_fe().tensor_degree(); + const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure(); + + for (unsigned k=0; k + void Estimator::face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const + { + const FEValuesBase &fe = info1.fe_values(); + const std::vector &uh1 = info1.values[0][0]; + const std::vector &uh2 = info2.values[0][0]; + const std::vector > &Duh1 = info1.gradients[0][0]; + const std::vector > &Duh2 = info2.gradients[0][0]; + + const unsigned int deg = fe.get_fe().tensor_degree(); + const double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure(); + const double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure(); + const double penalty = penalty1 + penalty2; + const double h = dinfo1.face->measure(); + + for (unsigned k=0; k + class ErrorIntegrator : public MeshWorker::LocalIntegrator + { + public: + void cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; + }; + + + template + void ErrorIntegrator::cell( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + std::vector > exact_gradients(fe.n_quadrature_points); + std::vector exact_values(fe.n_quadrature_points); + + exact_solution.gradient_list(fe.get_quadrature_points(), exact_gradients); + exact_solution.value_list(fe.get_quadrature_points(), exact_values); + + const std::vector > &Duh = info.gradients[0][0]; + const std::vector &uh = info.values[0][0]; + + for (unsigned k=0; k + void ErrorIntegrator::boundary( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + + std::vector exact_values(fe.n_quadrature_points); + exact_solution.value_list(fe.get_quadrature_points(), exact_values); + + const std::vector &uh = info.values[0][0]; + + const unsigned int deg = fe.get_fe().tensor_degree(); + const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure(); + + for (unsigned k=0; k + void ErrorIntegrator::face( + MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const + { + const FEValuesBase &fe = info1.fe_values(); + const std::vector &uh1 = info1.values[0][0]; + const std::vector &uh2 = info2.values[0][0]; + + const unsigned int deg = fe.get_fe().tensor_degree(); + const double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure(); + const double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure(); + const double penalty = penalty1 + penalty2; + + for (unsigned k=0; k + class InteriorPenaltyProblem + { + public: + typedef MeshWorker::IntegrationInfo CellInfo; + + InteriorPenaltyProblem(const FiniteElement &fe); + + void run(unsigned int n_steps); + + private: + void setup_system (); + void assemble_matrix (); + void assemble_mg_matrix (); + void assemble_right_hand_side (); + void error (); + double estimate (); + void solve (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + const MappingQ1 mapping; + const FiniteElement &fe; + MGDoFHandler mg_dof_handler; + DoFHandler &dof_handler; + MGConstrainedDoFs mg_constraints; + + SparsityPattern sparsity; + SparseMatrix matrix; + Vector solution; + Vector right_hand_side; + BlockVector estimates; + + MGLevelObject mg_sparsity; + MGLevelObject > mg_matrix; + + MGLevelObject mg_sparsity_dg_interface; + MGLevelObject > mg_matrix_dg_down; + MGLevelObject > mg_matrix_dg_up; + MGLevelObject > mg_matrix_in_out; + }; + + + template + InteriorPenaltyProblem::InteriorPenaltyProblem(const FiniteElement &fe) + : + mapping(), + fe(fe), + mg_dof_handler(triangulation), + dof_handler(mg_dof_handler), + estimates(1) + { + GridGenerator::hyper_cube_slit(triangulation, -1, 1); + } + + + template + void + InteriorPenaltyProblem::setup_system() + { + dof_handler.distribute_dofs(fe); + unsigned int n_dofs = dof_handler.n_dofs(); + solution.reinit(n_dofs); + right_hand_side.reinit(n_dofs); + + mg_constraints.clear(); + mg_constraints.initialize(dof_handler); + + CompressedSparsityPattern c_sparsity(n_dofs); + DoFTools::make_flux_sparsity_pattern(dof_handler, c_sparsity); + sparsity.copy_from(c_sparsity); + matrix.reinit(sparsity); + + const unsigned int n_levels = triangulation.n_levels(); + mg_matrix.resize(0, n_levels-1); + mg_matrix.clear(); + mg_matrix_dg_up.resize(0, n_levels-1); + mg_matrix_dg_up.clear(); + mg_matrix_dg_down.resize(0, n_levels-1); + mg_matrix_dg_down.clear(); + mg_matrix_in_out.resize(0, n_levels-1); + mg_matrix_in_out.clear(); + mg_sparsity.resize(0, n_levels-1); + mg_sparsity_dg_interface.resize(0, n_levels-1); + + for (unsigned int level=mg_sparsity.min_level(); + level<=mg_sparsity.max_level(); ++level) + { + CompressedSparsityPattern c_sparsity(mg_dof_handler.n_dofs(level)); + MGTools::make_flux_sparsity_pattern(mg_dof_handler, c_sparsity, level); + mg_sparsity[level].copy_from(c_sparsity); + mg_matrix[level].reinit(mg_sparsity[level]); + mg_matrix_in_out[level].reinit(mg_sparsity[level]); + + if (level>0) + { + CompressedSparsityPattern ci_sparsity; + ci_sparsity.reinit(mg_dof_handler.n_dofs(level-1), mg_dof_handler.n_dofs(level)); + MGTools::make_flux_sparsity_pattern_edge(mg_dof_handler, ci_sparsity, level); + mg_sparsity_dg_interface[level].copy_from(ci_sparsity); + mg_matrix_dg_up[level].reinit(mg_sparsity_dg_interface[level]); + mg_matrix_dg_down[level].reinit(mg_sparsity_dg_interface[level]); + } + } + } + + + template + void + InteriorPenaltyProblem::assemble_matrix() + { + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = update_values | update_gradients; + info_box.add_update_flags_all(update_flags); + info_box.initialize(fe, mapping); + + MeshWorker::DoFInfo dof_info(dof_handler); + + MeshWorker::Assembler::MatrixSimple > assembler; + assembler.initialize(matrix); + + MatrixIntegrator integrator; + MeshWorker::integration_loop( + dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, + integrator, assembler); + } + + + template + void + InteriorPenaltyProblem::assemble_mg_matrix() + { + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = update_values | update_gradients; + info_box.add_update_flags_all(update_flags); + info_box.initialize(fe, mapping); + + MeshWorker::DoFInfo dof_info(mg_dof_handler); + + MeshWorker::Assembler::MGMatrixSimple > assembler; + assembler.initialize(mg_matrix); + assembler.initialize(mg_constraints); + assembler.initialize_interfaces(mg_matrix_in_out, mg_matrix_in_out); + assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down); + + MatrixIntegrator integrator; + MeshWorker::integration_loop ( + mg_dof_handler.begin(), mg_dof_handler.end(), + dof_info, info_box, + integrator, assembler); + + for (unsigned int level=mg_matrix_in_out.min_level(); + level<=mg_matrix_in_out.min_level(); ++level) + if (mg_matrix_in_out[level].frobenius_norm() != 0.) + deallog << "Oops!" << std::endl; + } + + + template + void + InteriorPenaltyProblem::assemble_right_hand_side() + { + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients; + info_box.add_update_flags_all(update_flags); + info_box.initialize(fe, mapping); + + MeshWorker::DoFInfo dof_info(dof_handler); + + MeshWorker::Assembler::ResidualSimple > assembler; + NamedData* > data; + Vector *rhs = &right_hand_side; + data.add(rhs, "RHS"); + assembler.initialize(data); + + RHSIntegrator integrator; + MeshWorker::integration_loop( + dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, + integrator, assembler); + + right_hand_side *= -1.; + } + + + template + void + InteriorPenaltyProblem::solve() + { + SolverControl control(1000, 1.e-12); + SolverCG > solver(control); + + MGTransferPrebuilt > mg_transfer; + mg_transfer.build_matrices(mg_dof_handler); + + FullMatrix coarse_matrix; + coarse_matrix.copy_from (mg_matrix[0]); + MGCoarseGridHouseholder > mg_coarse; + mg_coarse.initialize(coarse_matrix); + + GrowingVectorMemory > mem; + typedef PreconditionSOR > RELAXATION; + mg::SmootherRelaxation > + mg_smoother; + RELAXATION::AdditionalData smoother_data(1.); + mg_smoother.initialize(mg_matrix, smoother_data); + + mg_smoother.set_steps(2); + mg_smoother.set_symmetric(true); + mg_smoother.set_variable(false); + + MGMatrix, Vector > mgmatrix(&mg_matrix); + MGMatrix, Vector > mgdown(&mg_matrix_dg_down); + MGMatrix, Vector > mgup(&mg_matrix_dg_up); + MGMatrix, Vector > mgedge(&mg_matrix_in_out); + + Multigrid > mg(mg_dof_handler, mgmatrix, + mg_coarse, mg_transfer, + mg_smoother, mg_smoother); + mg.set_edge_flux_matrices(mgdown, mgup); + mg.set_edge_matrices(mgedge, mgedge); + + PreconditionMG, + MGTransferPrebuilt > > + preconditioner(mg_dof_handler, mg, mg_transfer); + solver.solve(matrix, solution, right_hand_side, preconditioner); + } + + + template + double + InteriorPenaltyProblem::estimate() + { + std::vector old_user_indices; + triangulation.save_user_indices(old_user_indices); + + estimates.block(0).reinit(triangulation.n_active_cells()); + unsigned int i=0; + for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell,++i) + cell->set_user_index(i); + + MeshWorker::IntegrationInfoBox info_box; + const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; + info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points); + + NamedData* > solution_data; + solution_data.add(&solution, "solution"); + + info_box.cell_selector.add("solution", false, false, true); + info_box.boundary_selector.add("solution", true, true, false); + info_box.face_selector.add("solution", true, true, false); + + info_box.add_update_flags_boundary(update_quadrature_points); + info_box.initialize(fe, mapping, solution_data); + + MeshWorker::DoFInfo dof_info(dof_handler); + + MeshWorker::Assembler::CellsAndFaces assembler; + NamedData* > out_data; + BlockVector *est = &estimates; + out_data.add(est, "cells"); + assembler.initialize(out_data, false); + + Estimator integrator; + MeshWorker::integration_loop ( + dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, + integrator, assembler); + + triangulation.load_user_indices(old_user_indices); + return estimates.block(0).l2_norm(); + } + + + template + void + InteriorPenaltyProblem::error() + { + BlockVector errors(2); + errors.block(0).reinit(triangulation.n_active_cells()); + errors.block(1).reinit(triangulation.n_active_cells()); + unsigned int i=0; + for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell,++i) + cell->set_user_index(i); + + MeshWorker::IntegrationInfoBox info_box; + const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; + info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points); + + NamedData* > solution_data; + solution_data.add(&solution, "solution"); + + info_box.cell_selector.add("solution", true, true, false); + info_box.boundary_selector.add("solution", true, false, false); + info_box.face_selector.add("solution", true, false, false); + + info_box.add_update_flags_cell(update_quadrature_points); + info_box.add_update_flags_boundary(update_quadrature_points); + info_box.initialize(fe, mapping, solution_data); + + MeshWorker::DoFInfo dof_info(dof_handler); + + MeshWorker::Assembler::CellsAndFaces assembler; + NamedData* > out_data; + BlockVector *est = &errors; + out_data.add(est, "cells"); + assembler.initialize(out_data, false); + + ErrorIntegrator integrator; + MeshWorker::integration_loop ( + dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, + integrator, assembler); + + deallog << "energy-error: " << errors.block(0).l2_norm() << std::endl; + deallog << "L2-error: " << errors.block(1).l2_norm() << std::endl; + } + + + template + void InteriorPenaltyProblem::output_results (const unsigned int cycle) const + { + char *fn = new char[100]; + sprintf(fn, "step-39-02/sol-%02d", cycle); + + std::string filename(fn); + filename += ".gnuplot"; + deallog << "Writing solution to <" << filename << ">..." + << std::endl << std::endl; + std::ofstream gnuplot_output (filename.c_str()); + + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "u"); + data_out.add_data_vector (estimates.block(0), "est"); + + data_out.build_patches (); + + data_out.write_gnuplot(gnuplot_output); + } + + template + void + InteriorPenaltyProblem::run(unsigned int n_steps) + { + deallog << "Element: " << fe.get_name() << std::endl; + for (unsigned int s=0; s fe1(2); + InteriorPenaltyProblem<2> test1(fe1); + test1.run(6); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/multigrid/step-39-02a/cmp/generic b/tests/multigrid/step-39-02a/cmp/generic new file mode 100644 index 0000000000..0af72150d3 --- /dev/null +++ b/tests/multigrid/step-39-02a/cmp/generic @@ -0,0 +1,74 @@ + +DEAL::Element: FE_DGQ<2>(2) +DEAL::Step 0 +DEAL::Triangulation 16 cells, 2 levels +DEAL::DoFHandler 144 dofs, level dofs 36 144 +DEAL::Assemble matrix +DEAL::Assemble multilevel matrix +DEAL::Assemble right hand side +DEAL::Solve +DEAL:cg::Starting value 23.0730 +DEAL:cg::Convergence step 11 value 5.78162e-14 +DEAL::energy-error: 0.439211 +DEAL::L2-error: 0.0109342 +DEAL::Estimate 0.979555 +DEAL::Step 1 +DEAL::Triangulation 25 cells, 3 levels +DEAL::DoFHandler 225 dofs, level dofs 36 144 108 +DEAL::Assemble matrix +DEAL::Assemble multilevel matrix +DEAL::Assemble right hand side +DEAL::Solve +DEAL:cg::Starting value 23.0730 +DEAL:cg::Convergence step 11 value 9.37499e-13 +DEAL::energy-error: 0.332582 +DEAL::L2-error: 0.00548996 +DEAL::Estimate 0.838060 +DEAL::Step 2 +DEAL::Triangulation 37 cells, 4 levels +DEAL::DoFHandler 333 dofs, level dofs 36 144 144 108 +DEAL::Assemble matrix +DEAL::Assemble multilevel matrix +DEAL::Assemble right hand side +DEAL::Solve +DEAL:cg::Starting value 23.0730 +DEAL:cg::Convergence step 12 value 3.41649e-13 +DEAL::energy-error: 0.237705 +DEAL::L2-error: 0.00250869 +DEAL::Estimate 0.611449 +DEAL::Step 3 +DEAL::Triangulation 58 cells, 5 levels +DEAL::DoFHandler 522 dofs, level dofs 36 144 180 180 144 +DEAL::Assemble matrix +DEAL::Assemble multilevel matrix +DEAL::Assemble right hand side +DEAL::Solve +DEAL:cg::Starting value 23.4898 +DEAL:cg::Convergence step 12 value 5.58516e-13 +DEAL::energy-error: 0.170860 +DEAL::L2-error: 0.00120805 +DEAL::Estimate 0.452418 +DEAL::Step 4 +DEAL::Triangulation 85 cells, 6 levels +DEAL::DoFHandler 765 dofs, level dofs 36 144 360 180 180 108 +DEAL::Assemble matrix +DEAL::Assemble multilevel matrix +DEAL::Assemble right hand side +DEAL::Solve +DEAL:cg::Starting value 25.9490 +DEAL:cg::Convergence step 13 value 2.71720e-13 +DEAL::energy-error: 0.122755 +DEAL::L2-error: 0.000602816 +DEAL::Estimate 0.332525 +DEAL::Step 5 +DEAL::Triangulation 130 cells, 7 levels +DEAL::DoFHandler 1170 dofs, level dofs 36 144 432 432 180 180 144 +DEAL::Assemble matrix +DEAL::Assemble multilevel matrix +DEAL::Assemble right hand side +DEAL::Solve +DEAL:cg::Starting value 27.1421 +DEAL:cg::Convergence step 13 value 2.50045e-13 +DEAL::energy-error: 0.0869445 +DEAL::L2-error: 0.000292783 +DEAL::Estimate 0.236647