From 04b375716f5681bd0ca8b1cd4ce307760fb92e0c Mon Sep 17 00:00:00 2001 From: janssen <janssen@0785d39b-7218-0410-832d-ea1e28bc413d> Date: Wed, 17 Mar 2010 15:50:39 +0000 Subject: [PATCH] boundary_indices are not set correctly for systems git-svn-id: https://svn.dealii.org/trunk@20840 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/multigrid/mg_renumbered_02.cc | 49 +++++++++++++++++++++++------ tests/multigrid/mg_renumbered_03.cc | 47 ++++++++++++++++++++------- 2 files changed, 75 insertions(+), 21 deletions(-) diff --git a/tests/multigrid/mg_renumbered_02.cc b/tests/multigrid/mg_renumbered_02.cc index 0d0e367662..7da2d3b221 100644 --- a/tests/multigrid/mg_renumbered_02.cc +++ b/tests/multigrid/mg_renumbered_02.cc @@ -87,7 +87,6 @@ void initialize (const MGDoFHandler<dim> &dof, } - template <int dim> void initialize (const MGDoFHandler<dim> &dof, MGLevelObject<Vector<double> > &u) @@ -193,6 +192,7 @@ class LaplaceProblem private: void setup_system (); void test (); + void test_boundary (); void output_gpl(const MGDoFHandler<dim> &dof, MGLevelObject<Vector<double> > &v, const bool renumbered); @@ -205,12 +205,9 @@ class LaplaceProblem MGDoFHandler<dim> mg_dof_handler_renumbered; const unsigned int degree; + std::vector<std::set<unsigned int> > + boundary_indices_renumbered; - //std::vector<std::set<unsigned int> > - //boundary_indices; - - //std::vector<std::set<unsigned int> > - //boundary_indices_renumbered; }; @@ -245,6 +242,8 @@ void LaplaceProblem<dim>::setup_system () //DoFRenumbering::Cuthill_McKee (mg_dof_handler_renumbered, level); } + boundary_indices_renumbered.resize(nlevels); + deallog << "Number of degrees of freedom: " << mg_dof_handler_renumbered.n_dofs(); @@ -302,13 +301,16 @@ LaplaceProblem<dim>::output_gpl(const MGDoFHandler<dim> &dof, } - - template <int dim> void LaplaceProblem<dim>::test () { + typename FunctionMap<dim>::type dirichlet_boundary; + ZeroFunction<dim> dirichlet_bc(fe.n_components()); + dirichlet_boundary[0] = &dirichlet_bc; + MGTools::make_boundary_list (mg_dof_handler_renumbered, dirichlet_boundary, + boundary_indices_renumbered); MGTransferPrebuilt<Vector<double> > mg_transfer_renumbered; - mg_transfer_renumbered.build_matrices(mg_dof_handler_renumbered); + mg_transfer_renumbered.build_matrices(mg_dof_handler_renumbered, boundary_indices_renumbered); Vector<double> test; test.reinit(mg_dof_handler_renumbered.n_dofs()); @@ -332,6 +334,33 @@ void LaplaceProblem<dim>::test () } + template <int dim> +void LaplaceProblem<dim>::test_boundary () +{ + typename FunctionMap<dim>::type dirichlet_boundary; + ZeroFunction<dim> dirichlet_bc(fe.n_components()); + dirichlet_boundary[0] = &dirichlet_bc; + MGTools::make_boundary_list (mg_dof_handler_renumbered, dirichlet_boundary, + boundary_indices_renumbered); + + MGLevelObject<Vector<double> > u(0, triangulation.n_levels()-1); + MGLevelObject<Vector<double> > d(0, triangulation.n_levels()-1); + + for(unsigned int l=0; l<triangulation.n_levels(); ++l) + { + u[l].reinit(mg_dof_handler_renumbered.n_dofs(l)); + u[l] = 1.; + + for(std::set<unsigned int>::const_iterator i = boundary_indices_renumbered[l].begin(); + i!=boundary_indices_renumbered[l].end(); ++i) + { + u[l](*i) = 0; + } + } + output_gpl(mg_dof_handler_renumbered, u, false); +} + + template <int dim> void LaplaceProblem<dim>::refine_local () @@ -384,7 +413,7 @@ void LaplaceProblem<dim>::run () //triangulation.refine_global (1); refine_local (); setup_system (); - test(); + test_boundary(); }; } diff --git a/tests/multigrid/mg_renumbered_03.cc b/tests/multigrid/mg_renumbered_03.cc index a640fe1935..bd3d0bd4b6 100644 --- a/tests/multigrid/mg_renumbered_03.cc +++ b/tests/multigrid/mg_renumbered_03.cc @@ -101,6 +101,32 @@ void initialize (const MGDoFHandler<dim> &dof, } +template <int dim> +void diff (Vector<double> &diff, const MGDoFHandler<dim> &dof, + const Vector<double> &v, const unsigned int level) +{ + diff.reinit (v); + const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell; + std::vector<unsigned int> mg_dof_indices(dofs_per_cell); + const unsigned int n_comp = dof.get_fe().n_components(); + for (typename MGDoFHandler<dim>::cell_iterator + cell= dof.begin(level); + cell != dof.end(level); ++cell) + { + cell->get_mg_dof_indices(mg_dof_indices); + for(unsigned int i=0; i<dofs_per_cell/n_comp; ++i) + { + const unsigned int ni = n_comp* i; + const unsigned int i0 = mg_dof_indices[ni]; + const unsigned int i1 = mg_dof_indices[ni+1]; + const unsigned int i2 = mg_dof_indices[ni+2]; + diff(i0) = 2*v(i0) - v(i1); + diff(i1) = 3*v(i1) - 2*v(i2); + diff(i2) = v(i2) - 3*v(i0); + } + } +} + template <int dim> void diff (Vector<double> &diff, const MGDoFHandler<dim> &dof_1, const MGDoFHandler<dim> &dof_2, const Vector<double> &u, const Vector<double> &v, const unsigned int level) @@ -197,12 +223,9 @@ class LaplaceProblem MGDoFHandler<dim> mg_dof_handler_renumbered; const unsigned int degree; - std::vector<std::set<unsigned int> > - boundary_indices; + boundary_indices, boundary_indices_renumbered; - std::vector<std::set<unsigned int> > - boundary_indices_renumbered; }; @@ -211,7 +234,7 @@ LaplaceProblem<dim>::LaplaceProblem (const unsigned int deg) : triangulation (Triangulation<dim>::limit_level_difference_at_vertices), //fe (deg), //fe (FE_Q<dim> (deg),2, FE_Q<dim> (deg),2), - fe (FE_Q<dim> (deg),2), + fe (FE_Q<dim> (deg),3), mg_dof_handler (triangulation), mg_dof_handler_renumbered (triangulation), degree(deg) @@ -239,6 +262,9 @@ void LaplaceProblem<dim>::setup_system () //DoFRenumbering::Cuthill_McKee (mg_dof_handler_renumbered, level); } + boundary_indices.resize(nlevels); + boundary_indices_renumbered.resize(nlevels); + deallog << "Number of degrees of freedom: " << mg_dof_handler.n_dofs(); @@ -247,8 +273,6 @@ void LaplaceProblem<dim>::setup_system () << mg_dof_handler.n_dofs(l); deallog << std::endl; - boundary_indices.resize(triangulation.n_levels()); - boundary_indices_renumbered.resize(triangulation.n_levels()); } template <int dim> @@ -260,7 +284,7 @@ LaplaceProblem<dim>::output_gpl(const MGDoFHandler<dim> &dof, MeshWorker::IntegrationWorker<dim> integration_worker; MeshWorker::Assembler::GnuplotPatch assembler; - const unsigned int n_gauss_points = dof.get_fe().tensor_degree()+1; + const unsigned int n_gauss_points = dof.get_fe().tensor_degree(); QTrapez<1> trapez; QIterated<dim> quadrature(trapez, n_gauss_points); integration_worker.cell_quadrature = quadrature; @@ -273,7 +297,7 @@ LaplaceProblem<dim>::output_gpl(const MGDoFHandler<dim> &dof, cs.add("mg_vector"); integration_worker.cell_selector = cs; - assembler.initialize(dim, quadrature.size(), dim+mg_dof_handler.get_fe().n_components()); + assembler.initialize(dim, quadrature.size(), dim+dof.get_fe().n_components()); MeshWorker::IntegrationInfoBox<dim> info_box; MeshWorker::DoFInfo<dim> dof_info(dof); info_box.initialize(integration_worker, fe, mapping, data); @@ -332,14 +356,15 @@ void LaplaceProblem<dim>::test () //output_gpl(mg_dof_handler_renumbered, u, true); for(unsigned int l=0; l<triangulation.n_levels(); ++l) { - diff(d[l], mg_dof_handler, mg_dof_handler_renumbered, v[l],u[l],l); + diff(d[l], mg_dof_handler_renumbered, u[l],l); + //diff(d[l], mg_dof_handler, mg_dof_handler_renumbered, v[l],u[l],l); deallog << l << " " << u[l].l2_norm() << '\t' << v[l].l2_norm() << '\t' << d[l].l2_norm()<< std::endl; for(unsigned int i=0; i<d[l].size(); ++i) if(d[l](i)!=0) deallog << i << " " << d[l](i) << std::endl; } - output_gpl(mg_dof_handler, d, false); + output_gpl(mg_dof_handler_renumbered, d, false); } -- 2.39.5