From: janssen Date: Wed, 21 Aug 2013 20:14:22 +0000 (+0000) Subject: run into trouble X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=517b993103200ed1822b4186b9774cecd3dfc49e;p=dealii-svn.git run into trouble git-svn-id: https://svn.dealii.org/trunk@30385 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/multigrid/step-16-02.cc b/tests/multigrid/step-16-02.cc index 01bd6949bb..2c3283c5ae 100644 --- a/tests/multigrid/step-16-02.cc +++ b/tests/multigrid/step-16-02.cc @@ -42,6 +42,8 @@ #include #include +#include + #include #include @@ -60,10 +62,64 @@ #include #include +#include +#include +#include +#include +#include + #include #include using namespace dealii; +using namespace LocalIntegrators; + +template +class LaplaceMatrix : public MeshWorker::LocalIntegrator +{ +public: + LaplaceMatrix(); + virtual void cell(MeshWorker::DoFInfo& dinfo, MeshWorker::IntegrationInfo& info) const; + virtual void boundary(MeshWorker::DoFInfo& dinfo, MeshWorker::IntegrationInfo& info) const; + virtual void face(MeshWorker::DoFInfo& dinfo1, MeshWorker::DoFInfo& dinfo2, + MeshWorker::IntegrationInfo& info1, MeshWorker::IntegrationInfo& info2) const; +}; + + +template +LaplaceMatrix::LaplaceMatrix() +{} + + +template +void LaplaceMatrix::cell(MeshWorker::DoFInfo& dinfo, MeshWorker::IntegrationInfo& info) const +{ + AssertDimension (dinfo.n_matrices(), 1); + Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values(0)); +} + + +template +void LaplaceMatrix::boundary(MeshWorker::DoFInfo& dinfo, + typename MeshWorker::IntegrationInfo& info) const +{ + const unsigned int deg = info.fe_values(0).get_fe().tensor_degree(); + Laplace::nitsche_matrix(dinfo.matrix(0,false).matrix, info.fe_values(0), + Laplace::compute_penalty(dinfo, dinfo, deg, deg)); +} + + +template +void LaplaceMatrix::face( + MeshWorker::DoFInfo& dinfo1, MeshWorker::DoFInfo& dinfo2, + MeshWorker::IntegrationInfo& info1, MeshWorker::IntegrationInfo& info2) const +{ + const unsigned int deg = info1.fe_values(0).get_fe().tensor_degree(); + 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), + Laplace::compute_penalty(dinfo1, dinfo2, deg, deg)); +} template class LaplaceProblem @@ -75,7 +131,7 @@ public: private: void setup_system (); void assemble_system (); - void assemble_multigrid (); + void assemble_multigrid (const bool &use_mw); void solve (); void refine_grid (); void output_results (const unsigned int cycle) const; @@ -94,10 +150,12 @@ private: Vector system_rhs; const unsigned int degree; + LaplaceMatrix matrix_integrator; MGLevelObject mg_sparsity_patterns; MGLevelObject > mg_matrices; - MGLevelObject > mg_interface_matrices; + MGLevelObject > mg_interface_in; + MGLevelObject > mg_interface_out; MGConstrainedDoFs mg_constrained_dofs; }; @@ -155,7 +213,8 @@ LaplaceProblem::LaplaceProblem (const unsigned int degree) limit_level_difference_at_vertices), fe (degree), mg_dof_handler (triangulation), - degree(degree) + degree(degree), + matrix_integrator() {} @@ -203,8 +262,10 @@ void LaplaceProblem::setup_system () mg_constrained_dofs.initialize(mg_dof_handler, dirichlet_boundary); const unsigned int n_levels = triangulation.n_levels(); - mg_interface_matrices.resize(0, n_levels-1); - mg_interface_matrices.clear (); + mg_interface_in.resize(0, n_levels-1); + mg_interface_in.clear (); + mg_interface_out.resize(0, n_levels-1); + mg_interface_out.clear (); mg_matrices.resize(0, n_levels-1); mg_matrices.clear (); mg_sparsity_patterns.resize(0, n_levels-1); @@ -219,7 +280,8 @@ void LaplaceProblem::setup_system () mg_sparsity_patterns[level].copy_from (csp); mg_matrices[level].reinit(mg_sparsity_patterns[level]); - mg_interface_matrices[level].reinit(mg_sparsity_patterns[level]); + mg_interface_in[level].reinit(mg_sparsity_patterns[level]); + mg_interface_out[level].reinit(mg_sparsity_patterns[level]); } } @@ -280,124 +342,122 @@ void LaplaceProblem::assemble_system () template -void LaplaceProblem::assemble_multigrid () +void LaplaceProblem::assemble_multigrid (const bool& use_mw) { - QGauss quadrature_formula(1+degree); + if(use_mw == true) + { + mg_matrices = 0.; + + MappingQ1 mapping; + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = update_values | update_gradients | update_hessians; + info_box.add_update_flags_all(update_flags); + info_box.initialize(fe, mapping);//, &dof_handler.block_info()); + + MeshWorker::DoFInfo dof_info(mg_dof_handler);//.block_info()); + + MeshWorker::Assembler::MGMatrixSimple > assembler; + assembler.initialize(mg_constrained_dofs); + // assembler.initialize_local_blocks(dof_handler.block_info().local()); + assembler.initialize(mg_matrices); + assembler.initialize_interfaces(mg_interface_in, mg_interface_out); +// assembler.initialize_fluxes(mg_matrix_up, mg_matrix_down); + + MeshWorker::integration_loop ( + mg_dof_handler.begin(), mg_dof_handler.end(), + dof_info, info_box, matrix_integrator, assembler); + + const unsigned int nlevels = triangulation.n_levels(); + for (unsigned int level=0;level quadrature_formula(1+degree); - FEValues fe_values (fe, quadrature_formula, - update_values | update_gradients | - update_quadrature_points | update_JxW_values); + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); - const Coefficient coefficient; - std::vector coefficient_values (n_q_points); + const Coefficient coefficient; + std::vector coefficient_values (n_q_points); - std::vector > interface_dofs - = mg_constrained_dofs.get_refinement_edge_indices (); - std::vector > boundary_interface_dofs - = mg_constrained_dofs.get_refinement_edge_boundary_indices (); + std::vector > interface_dofs + = mg_constrained_dofs.get_refinement_edge_indices (); + std::vector > boundary_interface_dofs + = mg_constrained_dofs.get_refinement_edge_boundary_indices (); - std::vector boundary_constraints (triangulation.n_levels()); - std::vector boundary_interface_constraints (triangulation.n_levels()); - for (unsigned int level=0; level boundary_constraints (triangulation.n_levels()); + std::vector boundary_interface_constraints (triangulation.n_levels()); + for (unsigned int level=0; level::cell_iterator cell = mg_dof_handler.begin(), - endc = mg_dof_handler.end(); + typename MGDoFHandler::cell_iterator cell = mg_dof_handler.begin(), + endc = mg_dof_handler.end(); - for (; cell!=endc; ++cell) + for (; cell!=endc; ++cell) { cell_matrix = 0; fe_values.reinit (cell); coefficient.value_list (fe_values.get_quadrature_points(), - coefficient_values); + coefficient_values); for (unsigned int q_point=0; q_pointget_mg_dof_indices (local_dof_indices); boundary_constraints[cell->level()] - .distribute_local_to_global (cell_matrix, - local_dof_indices, - mg_matrices[cell->level()]); - - // The next step is again slightly more - // obscure (but explained in the @ref - // mg_paper): We need the remainder of - // the operator that we just copied - // into the mg_matrices - // object, namely the part on the - // interface between cells at the - // current level and cells one level - // coarser. This matrix exists in two - // directions: for interior DoFs (index - // $i$) of the current level to those - // sitting on the interface (index - // $j$), and the other way around. Of - // course, since we have a symmetric - // operator, one of these matrices is - // the transpose of the other. - // - // The way we assemble these matrices - // is as follows: since the are formed - // from parts of the local - // contributions, we first delete all - // those parts of the local - // contributions that we are not - // interested in, namely all those - // elements of the local matrix for - // which not $i$ is an interface DoF - // and $j$ is not. The result is one of - // the two matrices that we are - // interested in, and we then copy it - // into the - // mg_interface_matrices - // object. The - // boundary_interface_constraints - // object at the same time makes sure - // that we delete contributions from - // all degrees of freedom that are not - // only on the interface but also on - // the external boundary of the domain. - // - // The last part to remember is how to - // get the other matrix. Since it is - // only the transpose, we will later - // (in the solve() - // function) be able to just pass the - // transpose matrix where necessary. + .distribute_local_to_global (cell_matrix, + local_dof_indices, + mg_matrices[cell->level()]); + for (unsigned int i=0; ilevel()][local_dof_indices[i]]==true && - interface_dofs[cell->level()][local_dof_indices[j]]==false)) + interface_dofs[cell->level()][local_dof_indices[j]]==false)) cell_matrix(i,j) = 0; boundary_interface_constraints[cell->level()] - .distribute_local_to_global (cell_matrix, - local_dof_indices, - mg_interface_matrices[cell->level()]); + .distribute_local_to_global (cell_matrix, + local_dof_indices, + mg_interface_in[cell->level()]); } + } + + const unsigned int nlevels = triangulation.n_levels(); + for (unsigned int level=1;level::solve () mg_smoother.set_symmetric(true); MGMatrix<> mg_matrix(&mg_matrices); - MGMatrix<> mg_interface_up(&mg_interface_matrices); - MGMatrix<> mg_interface_down(&mg_interface_matrices); + MGMatrix<> mg_interface_up(&mg_interface_in); + MGMatrix<> mg_interface_down(&mg_interface_out); Multigrid > mg(mg_dof_handler, mg_matrix, @@ -586,7 +646,7 @@ void LaplaceProblem::run () deallog << std::endl; assemble_system (); - assemble_multigrid (); + assemble_multigrid (true); solve (); // output_results (cycle);