From c48206771eec1f74c57fbd357d2b40efbfc7d9ef Mon Sep 17 00:00:00 2001 From: kanschat Date: Sat, 28 Nov 2009 04:42:33 +0000 Subject: [PATCH] MeshWorker::Assembler::MGMatrixSimple bug fixed git-svn-id: https://svn.dealii.org/trunk@20182 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/numerics/mesh_worker_assembler.h | 49 ++++++++---- deal.II/examples/step-39/step-39.cc | 78 +++++++++++++------ 2 files changed, 88 insertions(+), 39 deletions(-) diff --git a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h index e09a43f39b..33af1b0de0 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h @@ -579,8 +579,16 @@ namespace MeshWorker void assemble(MATRIX& G, const FullMatrix& M, const std::vector& i1, - const std::vector& i2, - const bool transpose = false); + const std::vector& i2); + + /** + * Assemble a single matrix + * into a global matrix. + */ + void assemble_transpose(MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2); /** * The global matrix being @@ -1178,11 +1186,11 @@ namespace MeshWorker template inline void - MGMatrixSimple::assemble(MATRIX& G, - const FullMatrix& M, - const std::vector& i1, - const std::vector& i2, - bool transpose) + MGMatrixSimple::assemble( + MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2) { AssertDimension(M.m(), i1.size()); AssertDimension(M.n(), i2.size()); @@ -1190,12 +1198,25 @@ namespace MeshWorker for (unsigned int j=0; j= threshold) - { - if (transpose) - G.add(i1[j], i2[k], M(j,k)); - else - G.add(i1[j], i2[k], M(j,k)); - } + G.add(i1[j], i2[k], M(j,k)); + } + + + template + inline void + MGMatrixSimple::assemble_transpose( + MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2) + { + AssertDimension(M.n(), i1.size()); + AssertDimension(M.m(), i2.size()); + + for (unsigned int j=0; j= threshold) + G.add(i1[j], i2[k], M(k,j)); } @@ -1232,7 +1253,7 @@ namespace MeshWorker // which is done by // the coarser cell assemble((*matrix)[level1], info1.M1[0].matrix, info1.indices, info1.indices); - assemble((*flux_up)[level1],info1.M2[0].matrix, info1.indices, info2.indices, true); + assemble_transpose((*flux_up)[level1],info1.M2[0].matrix, info2.indices, info1.indices); assemble((*flux_down)[level1], info2.M2[0].matrix, info2.indices, info1.indices); } } diff --git a/deal.II/examples/step-39/step-39.cc b/deal.II/examples/step-39/step-39.cc index 6a829c851e..f17955777c 100644 --- a/deal.II/examples/step-39/step-39.cc +++ b/deal.II/examples/step-39/step-39.cc @@ -28,6 +28,7 @@ // Include files for FiniteElement // classes and DoFHandler. #include +#include #include #include #include @@ -146,9 +147,16 @@ void MatrixIntegrator::face(typename MeshWorker::IntegrationWorker::Fa FullMatrix& matrix_v2u2 = info2.M1[0].matrix; const unsigned int deg = fe1.get_fe().tensor_degree(); - const double penalty1 = deg * (deg+1) * info1.face->measure() / info1.cell->measure(); - const double penalty2 = deg * (deg+1) * info2.face->measure() / info2.cell->measure(); - const double penalty = penalty1 + penalty2; + double penalty1 = deg * (deg+1) * info1.face->measure() / info1.cell->measure(); + double penalty2 = deg * (deg+1) * info2.face->measure() / info2.cell->measure(); + if (info1.cell->has_children() ^ info2.cell->has_children()) + { + Assert (info1.face == info2.face, ExcInternalError()); + Assert (info1.face->has_children(), ExcInternalError()); + Assert (info1.cell->has_children(), ExcInternalError()); + penalty1 *= 2; + } +const double penalty = penalty1 + penalty2; for (unsigned k=0;k triangulation; - const MappingQ1 mapping; + Triangulation triangulation; + const MappingQ1 mapping; const FiniteElement& fe; - MGDoFHandler dof_handler; + MGDoFHandler mg_dof_handler; + DoFHandler& dof_handler; SparsityPattern sparsity; SparseMatrix matrix; @@ -258,9 +267,12 @@ template Step39::Step39(const FiniteElement& fe) : fe(fe), - dof_handler(triangulation) + mg_dof_handler(triangulation), + dof_handler(mg_dof_handler) { GridGenerator::hyper_L(triangulation, -1, 1); + triangulation.begin()->set_refine_flag(); + triangulation.execute_coarsening_and_refinement(); } @@ -272,10 +284,8 @@ Step39::setup_system() unsigned int n_dofs = dof_handler.n_dofs(); CompressedSparsityPattern c_sparsity(n_dofs); - const DoFHandler& dof = dof_handler; - DoFTools::make_flux_sparsity_pattern(dof, c_sparsity); + DoFTools::make_flux_sparsity_pattern(dof_handler, c_sparsity); sparsity.copy_from(c_sparsity); - matrix.reinit(sparsity); solution.reinit(n_dofs); @@ -296,14 +306,14 @@ Step39::setup_system() for (unsigned int level=mg_sparsity.get_minlevel(); level<=mg_sparsity.get_maxlevel();++level) { - CompressedSparsityPattern c_sparsity(dof_handler.n_dofs(level)); + CompressedSparsityPattern c_sparsity(mg_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)); + ci_sparsity.reinit(mg_dof_handler.n_dofs(level-1), mg_dof_handler.n_dofs(level)); - MGTools::make_flux_sparsity_pattern(dof_handler, c_sparsity, level); + MGTools::make_flux_sparsity_pattern(mg_dof_handler, c_sparsity, level); if (level>0) - MGTools::make_flux_sparsity_pattern_edge(dof_handler, ci_sparsity, level); + MGTools::make_flux_sparsity_pattern_edge(mg_dof_handler, ci_sparsity, level); mg_sparsity[level].copy_from(c_sparsity); mg_matrix[level].reinit(mg_sparsity[level]); @@ -333,6 +343,8 @@ Step39::assemble_matrix() 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); + +// matrix.print_formatted(std::cout, 2, false, 5, "*"); } @@ -343,15 +355,29 @@ 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; + const unsigned int n_gauss_points = mg_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); + integrator.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down); + MeshWorker::IntegrationInfoBox info_box(mg_dof_handler); info_box.initialize(integrator, fe, mapping); - MeshWorker::integration_loop(dof_handler.begin(), dof_handler.end(), info_box, integrator); + MeshWorker::integration_loop(mg_dof_handler.begin(), mg_dof_handler.end(), info_box, integrator); + +// for (unsigned int l=0;l::solve() GrowingVectorMemory > mem; MGTransferPrebuilt > mg_transfer; - mg_transfer.build_matrices(dof_handler); + mg_transfer.build_matrices(mg_dof_handler); FullMatrix coarse_matrix; coarse_matrix.copy_from (mg_matrix[0]); MGCoarseGridHouseholder > mg_coarse; @@ -423,7 +449,7 @@ Step39::solve() // Now, we are ready to set up the // V-cycle operator and the // multilevel preconditioner. - Multigrid > mg(dof_handler, mgmatrix, + Multigrid > mg(mg_dof_handler, mgmatrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother); mg.set_edge_matrices(mgdown, mgup); @@ -432,7 +458,7 @@ Step39::solve() PreconditionMG, MGTransferPrebuilt > > - preconditioner(dof_handler, mg, mg_transfer); + preconditioner(mg_dof_handler, mg, mg_transfer); cg.solve(matrix, solution, right_hand_side, preconditioner); } @@ -480,6 +506,7 @@ template void Step39::run(unsigned int n_steps) { + deallog << "Element: " << fe.get_name() << std::endl; for (unsigned int s=0;s::run(unsigned int n_steps) << triangulation.n_levels() << " levels" << std::endl; setup_system(); - deallog << "DoFHandler " << dof_handler.n_dofs() << " dofs, levels"; + deallog << "DoFHandler " << dof_handler.n_dofs() << " dofs, level dofs"; for (unsigned int l=0;l::run(unsigned int n_steps) assemble_mg_matrix(); deallog << "Assemble right hand side" << std::endl; assemble_right_hand_side(); - + deallog << "Solve" << std::endl; solve(); + deallog << "Error" << std::endl; error(); output_results(s); } @@ -512,7 +540,7 @@ Step39::run(unsigned int n_steps) int main() { - FE_DGQ<2> dgq1(1); - Step39<2> test1(dgq1); + FE_DGQ<2> fe1(1); + Step39<2> test1(fe1); test1.run(7); } -- 2.39.5