From 21157015eb765518aca83d5b53c29f9cc49b0892 Mon Sep 17 00:00:00 2001 From: kanschat Date: Wed, 11 Nov 2009 13:23:09 +0000 Subject: [PATCH] working on inconsistency git-svn-id: https://svn.dealii.org/trunk@20089 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-39/step-39.cc | 198 ++++++++++++++++++++++++++-- 1 file changed, 185 insertions(+), 13 deletions(-) diff --git a/deal.II/examples/step-39/step-39.cc b/deal.II/examples/step-39/step-39.cc index 638037436f..1ab4b2f563 100644 --- a/deal.II/examples/step-39/step-39.cc +++ b/deal.II/examples/step-39/step-39.cc @@ -19,13 +19,15 @@ #include #include #include +#include // Include files for setting up the // mesh #include - // Include files for FinitElement + // Include files for FiniteElement // classes and DoFHandler. +#include #include #include #include @@ -35,9 +37,21 @@ #include #include + // Support for multigrid methods +#include +#include +#include +#include +#include +#include + // Finally, we take our exact - // solution from the library. + // solution from the library as well + // as quadrature and additional tools. #include +#include +#include +#include #include #include @@ -49,6 +63,11 @@ // as well. using namespace dealii; + // This is the function we use to set + // the boundary values and also the + // exact solution we compare to. +//Functions::LSingularityFunction exact_solution; +Functions::Q1WedgeFunction<2> exact_solution; // @sect3{The local integrators} @@ -177,16 +196,15 @@ void RHSIntegrator::bdry(typename MeshWorker::IntegrationWorker::FaceI const FEFaceValuesBase& fe = info.fe(); Vector& local_vector = info.R[0].block(0); - static Functions::LSingularityFunction lshaped; std::vector boundary_values(fe.n_quadrature_points); - lshaped.value_list(fe.get_quadrature_points(), boundary_values); + 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) * info.face->measure() / info.cell->measure(); for (unsigned k=0;k& fe; MGDoFHandler dof_handler; - SparsityPattern sparsity_pattern; + SparsityPattern sparsity; SparseMatrix matrix; Vector solution; Vector right_hand_side; + + MGLevelObject mg_sparsity; + MGLevelObject mg_sparsity_dg_interface; + MGLevelObject > mg_matrix; + MGLevelObject > mg_matrix_dg_up; + MGLevelObject > mg_matrix_dg_down; }; @@ -234,7 +260,7 @@ Step39::Step39(const FiniteElement& fe) fe(fe), dof_handler(triangulation) { - GridGenerator::hyper_L(triangulation); + GridGenerator::hyper_cube(triangulation, -1, 1); } @@ -248,12 +274,46 @@ Step39::setup_system() CompressedSparsityPattern c_sparsity(n_dofs); const DoFHandler& dof = dof_handler; DoFTools::make_flux_sparsity_pattern(dof, c_sparsity); - sparsity_pattern.copy_from(c_sparsity); + sparsity.copy_from(c_sparsity); - matrix.reinit(sparsity_pattern); + matrix.reinit(sparsity); solution.reinit(n_dofs); right_hand_side.reinit(n_dofs); + + 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_sparsity.resize(0, n_levels-1); + mg_sparsity_dg_interface.resize(0, n_levels-1); + + for (unsigned int level=mg_sparsity.get_minlevel(); + level<=mg_sparsity.get_maxlevel();++level) + { + CompressedSparsityPattern c_sparsity(dof_handler.n_dofs(level)); + CompressedSparsityPattern ci_sparsity; + if (level>0) + ci_sparsity.reinit(dof_handler.n_dofs(level-1), dof_handler.n_dofs(level)); + + MGTools::make_flux_sparsity_pattern(dof_handler, c_sparsity, level); + if (level>0) + MGTools::make_flux_sparsity_pattern_edge(dof_handler, ci_sparsity, level); + + mg_sparsity[level].copy_from(c_sparsity); + mg_matrix[level].reinit(mg_sparsity[level]); + if (level>0) + { + 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]); + } + } } @@ -276,6 +336,25 @@ Step39::assemble_matrix() } +template +void +Step39::assemble_mg_matrix() +{ + const MatrixIntegrator local; + MeshWorker::AssemblingIntegrator >, MatrixIntegrator > + integrator(local); + const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; + integrator.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points); + UpdateFlags update_flags = update_values | update_gradients; + integrator.add_update_flags(update_flags, true, true, true, true); + + integrator.initialize(mg_matrix); + MeshWorker::IntegrationInfoBox info_box(dof_handler); + info_box.initialize(integrator, fe, mapping); + MeshWorker::integration_loop(dof_handler.begin(), dof_handler.end(), info_box, integrator); +} + + template void Step39::assemble_right_hand_side() @@ -284,7 +363,7 @@ Step39::assemble_right_hand_side() MeshWorker::AssemblingIntegrator >, RHSIntegrator > integrator(local); const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; - integrator.initialize_gauss_quadrature(1, n_gauss_points, n_gauss_points); + integrator.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points); UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients; integrator.add_update_flags(update_flags, true, true, true, true); @@ -295,6 +374,8 @@ Step39::assemble_right_hand_side() MeshWorker::IntegrationInfoBox info_box(dof_handler); info_box.initialize(integrator, fe, mapping); MeshWorker::integration_loop(dof_handler.begin_active(), dof_handler.end(), info_box, integrator); + + right_hand_side *= -1.; } @@ -304,7 +385,94 @@ Step39::solve() { SolverControl control(1000, 1.e-12); SolverCG > cg(control); - cg.solve(matrix, solution, right_hand_side, PreconditionIdentity()); + + GrowingVectorMemory > mem; + MGTransferPrebuilt > mg_transfer; + mg_transfer.build_matrices(dof_handler); + FullMatrix coarse_matrix; + coarse_matrix.copy_from (mg_matrix[0]); + MGCoarseGridHouseholder > mg_coarse; + mg_coarse.initialize(coarse_matrix); + typedef PreconditionSOR > RELAXATION; + MGSmootherRelaxation, RELAXATION, Vector > + mg_smoother(mem); + RELAXATION::AdditionalData smoother_data(1.); +// RELAXATION::AdditionalData smoother_data(fe.dofs_per_cell, 1.); + mg_smoother.initialize(mg_matrix, smoother_data); + + // Do two smoothing steps per level + mg_smoother.set_steps(2); + // Since the SOR method is not + // symmetric, but we use conjugate + // gradient iteration below, here + // is a trick to make the + // multilevel preconditioner a + // symmetric operator even for + // nonsymmetric smoothers. + mg_smoother.set_symmetric(true); + mg_smoother.set_variable(false); + + // We must wrap our matrices in an + // object having the required + // multiplication functions. + MGMatrix, Vector > mgmatrix(&mg_matrix); + MGMatrix, Vector > mgdown(&mg_matrix_dg_down); + MGMatrix, Vector > mgup(&mg_matrix_dg_up); + + + // Now, we are ready to set up the + // V-cycle operator and the + // multilevel preconditioner. + Multigrid > mg(dof_handler, mgmatrix, + mg_coarse, mg_transfer, + mg_smoother, mg_smoother); + mg.set_edge_matrices(mgdown, mgup); + mg.set_debug(0); + mg_smoother.set_debug(0); + + PreconditionMG, + MGTransferPrebuilt > > + preconditioner(dof_handler, mg, mg_transfer); + + cg.solve(matrix, solution, right_hand_side, preconditioner); +} + + +template +void +Step39::error() +{ + const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+2; + Vector cell_errors(triangulation.n_active_cells()); + + QGauss quadrature(n_gauss_points); + VectorTools::integrate_difference(mapping, dof_handler, solution, exact_solution, + cell_errors, quadrature, VectorTools::L2_norm); + deallog << "Error " << cell_errors.l2_norm() << std::endl; +} + + +template +void Step39::output_results (const unsigned int cycle) const +{ + // Output of the solution in + // gnuplot format. + std::string filename = "sol-"; + filename += ('0' + cycle); + Assert (cycle < 10, ExcInternalError()); + + filename += ".gnuplot"; + std::cout << "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.build_patches (); + + data_out.write_gnuplot(gnuplot_output); } @@ -329,15 +497,19 @@ Step39::run(unsigned int n_steps) deallog << std::endl; assemble_matrix(); + assemble_mg_matrix(); assemble_right_hand_side(); + solve(); + error(); + output_results(s); } } int main() { - FE_DGQ<2> dgq1(1); + FE_Q<2> dgq1(1); Step39<2> test1(dgq1); test1.run(7); } -- 2.39.5