From 091d6271cf831489f7b9241427454ecf574d2ff5 Mon Sep 17 00:00:00 2001 From: janssen Date: Wed, 17 Mar 2010 00:05:41 +0000 Subject: [PATCH] small modifications git-svn-id: https://svn.dealii.org/trunk@20838 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/multigrid/mg_renumbered_02.cc | 118 ++++-------- tests/multigrid/mg_renumbered_03.cc | 273 +++++++++++++--------------- 2 files changed, 161 insertions(+), 230 deletions(-) diff --git a/tests/multigrid/mg_renumbered_02.cc b/tests/multigrid/mg_renumbered_02.cc index f65994abb2..0d0e367662 100644 --- a/tests/multigrid/mg_renumbered_02.cc +++ b/tests/multigrid/mg_renumbered_02.cc @@ -1,4 +1,4 @@ -#define baerbel_mg_test 1 +#define BAERBEL_MG_TEST 1 #include #include @@ -104,21 +104,28 @@ void initialize (const MGDoFHandler &dof, template -void diff (Vector &diff, const MGDoFHandler &dof_1, const MGDoFHandler &dof_2, - const Vector &u, const Vector &v, const unsigned int level) +void diff (Vector &diff, const MGDoFHandler &dof, + const Vector &v, const unsigned int level) { - diff.reinit (u); - const unsigned int dofs_per_cell = dof_1.get_fe().dofs_per_cell; - std::vector dof_indices_1(dofs_per_cell); - std::vector dof_indices_2(dofs_per_cell); + diff.reinit (v); + const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell; + std::vector mg_dof_indices(dofs_per_cell); + const unsigned int n_comp = dof.get_fe().n_components(); for (typename MGDoFHandler::cell_iterator - cell_1 = dof_1.begin(level), cell_2 = dof_2.begin(level); - cell_1 != dof_1.end(level); ++cell_1, ++cell_2) + cell= dof.begin(level); + cell != dof.end(level); ++cell) { - cell_1->get_mg_dof_indices(dof_indices_1); - cell_2->get_mg_dof_indices(dof_indices_2); - for(unsigned int i=0; iget_mg_dof_indices(mg_dof_indices); + for(unsigned int i=0; i::cell( const FEValuesBase& fe = info.fe_values(); const std::vector >& uh = info.values[0]; - if(dinfo.cell->level() == 3) - { - for(unsigned int i=0; i mapping; FESystem fe; //FE_Q fe; - MGDoFHandler mg_dof_handler; MGDoFHandler mg_dof_handler_renumbered; - MGLevelObject mg_sparsity_renumbered; - MGLevelObject > mg_matrices_renumbered; - MGLevelObject mg_sparsity; - MGLevelObject > mg_matrices; - const unsigned int degree; - std::vector > - boundary_indices; + //std::vector > + //boundary_indices; - std::vector > - boundary_indices_renumbered; + //std::vector > + //boundary_indices_renumbered; }; @@ -224,8 +219,7 @@ LaplaceProblem::LaplaceProblem (const unsigned int deg) : triangulation (Triangulation::limit_level_difference_at_vertices), //fe (deg), //fe (FE_Q (deg),2, FE_Q (deg),2), - fe (FE_Q (deg),4), - mg_dof_handler (triangulation), + fe (FE_Q (deg),3), mg_dof_handler_renumbered (triangulation), degree(deg) {} @@ -234,7 +228,6 @@ LaplaceProblem::LaplaceProblem (const unsigned int deg) : template void LaplaceProblem::setup_system () { - mg_dof_handler.distribute_dofs (fe); mg_dof_handler_renumbered.distribute_dofs (fe); std::vector block_component (2*dim,0); @@ -253,39 +246,12 @@ void LaplaceProblem::setup_system () } deallog << "Number of degrees of freedom: " - << mg_dof_handler.n_dofs(); + << mg_dof_handler_renumbered.n_dofs(); for (unsigned int l=0;l @@ -297,7 +263,7 @@ LaplaceProblem::output_gpl(const MGDoFHandler &dof, MeshWorker::IntegrationWorker 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 quadrature(trapez, n_gauss_points); integration_worker.cell_quadrature = quadrature; @@ -310,7 +276,7 @@ LaplaceProblem::output_gpl(const MGDoFHandler &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 info_box; MeshWorker::DoFInfo dof_info(dof); info_box.initialize(integration_worker, fe, mapping, data); @@ -341,42 +307,28 @@ LaplaceProblem::output_gpl(const MGDoFHandler &dof, template void LaplaceProblem::test () { - typename FunctionMap::type dirichlet_boundary; - ZeroFunction dirichlet_bc(fe.n_components()); - dirichlet_boundary[0] = &dirichlet_bc; - MGTools::make_boundary_list (mg_dof_handler, dirichlet_boundary, - boundary_indices); - MGTools::make_boundary_list (mg_dof_handler_renumbered, dirichlet_boundary, - boundary_indices_renumbered); - MGTransferPrebuilt > mg_transfer; - mg_transfer.build_matrices(mg_dof_handler, boundary_indices); MGTransferPrebuilt > mg_transfer_renumbered; - mg_transfer_renumbered.build_matrices(mg_dof_handler_renumbered, boundary_indices_renumbered); + mg_transfer_renumbered.build_matrices(mg_dof_handler_renumbered); Vector test; - test.reinit(mg_dof_handler.n_dofs()); + test.reinit(mg_dof_handler_renumbered.n_dofs()); - MGLevelObject > v(0, triangulation.n_levels()-1); MGLevelObject > u(0, triangulation.n_levels()-1); MGLevelObject > d(0, triangulation.n_levels()-1); - initialize(mg_dof_handler, test); - mg_transfer.copy_to_mg(mg_dof_handler, v, test); - //output_gpl(mg_dof_handler, v, false); - initialize(mg_dof_handler_renumbered, test); mg_transfer_renumbered.copy_to_mg(mg_dof_handler_renumbered, u, test); - //output_gpl(mg_dof_handler_renumbered, u, true); + for(unsigned int l=0; l #include #include @@ -33,6 +31,9 @@ #include #include #include +#include +#include +#include #include #include @@ -46,6 +47,7 @@ #include #include +#include using namespace dealii; @@ -67,7 +69,6 @@ template void initialize (const MGDoFHandler &dof, Vector &u) { - unsigned int counter=0; const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell; std::vector dof_indices(dofs_per_cell); for (typename MGDoFHandler::active_cell_iterator @@ -76,7 +77,10 @@ void initialize (const MGDoFHandler &dof, { cell->get_dof_indices(dof_indices); for(unsigned int i=0; i &dof, template -void print_diff (const MGDoFHandler &dof_1, const MGDoFHandler &dof_2, - const Vector &u, const Vector &v) +void diff (Vector &diff, const MGDoFHandler &dof_1, const MGDoFHandler &dof_2, + const Vector &u, const Vector &v, const unsigned int level) { - Vector diff; diff.reinit (u); const unsigned int dofs_per_cell = dof_1.get_fe().dofs_per_cell; std::vector dof_indices_1(dofs_per_cell); std::vector dof_indices_2(dofs_per_cell); - for (typename MGDoFHandler::active_cell_iterator - cell_1 = dof_1.begin_active(), cell_2 = dof_2.begin_active(); - cell_1 != dof_1.end(); ++cell_1, ++cell_2) + for (typename MGDoFHandler::cell_iterator + cell_1 = dof_1.begin(level), cell_2 = dof_2.begin(level); + cell_1 != dof_1.end(level); ++cell_1, ++cell_2) { - cell_1->get_dof_indices(dof_indices_1); - cell_2->get_dof_indices(dof_indices_2); + cell_1->get_mg_dof_indices(dof_indices_1); + cell_2->get_mg_dof_indices(dof_indices_2); for(unsigned int i=0; i @@ -139,37 +140,68 @@ void print(const MGDoFHandler &dof, std::vector > &interf } } +template +class OutputCreator : public Subscriptor +{ + public: + void cell(MeshWorker::DoFInfo& dinfo, + typename MeshWorker::IntegrationWorker::CellInfo& info); + +}; + +template +void OutputCreator::cell( + MeshWorker::DoFInfo& dinfo, + typename MeshWorker::IntegrationWorker::CellInfo& info) +{ + const FEValuesBase& fe = info.fe_values(); + const std::vector >& uh = info.values[0]; + + const unsigned int square_root = std::pow(fe.n_quadrature_points, 1./dim)+.5; + for (unsigned int k1=0; k1 class LaplaceProblem { public: + typedef typename MeshWorker::IntegrationWorker::CellInfo CellInfo; + typedef typename MeshWorker::IntegrationWorker::FaceInfo FaceInfo; + LaplaceProblem (const unsigned int deg); void run (); private: void setup_system (); void test (); - void test_renumber (); - void test_interface_dofs (); + void output_gpl(const MGDoFHandler &dof, + MGLevelObject > &v, + const bool renumbered); void refine_local (); Triangulation triangulation; + const MappingQ1 mapping; FESystem fe; + //FE_Q fe; MGDoFHandler mg_dof_handler; MGDoFHandler mg_dof_handler_renumbered; - MGLevelObject mg_sparsity_renumbered; - MGLevelObject > mg_matrices_renumbered; - MGLevelObject mg_sparsity; - MGLevelObject > mg_matrices; - const unsigned int degree; - std::vector > + std::vector > boundary_indices; - std::vector > + std::vector > boundary_indices_renumbered; }; @@ -177,7 +209,9 @@ class LaplaceProblem template LaplaceProblem::LaplaceProblem (const unsigned int deg) : triangulation (Triangulation::limit_level_difference_at_vertices), - fe (FE_Q (deg),2, FE_Q (deg),2), + //fe (deg), + //fe (FE_Q (deg),2, FE_Q (deg),2), + fe (FE_Q (deg),2), mg_dof_handler (triangulation), mg_dof_handler_renumbered (triangulation), degree(deg) @@ -197,11 +231,11 @@ void LaplaceProblem::setup_system () const unsigned int nlevels = triangulation.n_levels(); DoFHandler &dof = mg_dof_handler_renumbered; - DoFRenumbering::component_wise (dof, block_component); + DoFRenumbering::component_wise (dof);//, block_component); //DoFRenumbering::Cuthill_McKee (dof); for (unsigned int level=0;level::setup_system () boundary_indices.resize(triangulation.n_levels()); boundary_indices_renumbered.resize(triangulation.n_levels()); - - mg_matrices.resize(0, nlevels-1); - mg_matrices.clear (); - mg_matrices_renumbered.resize(0, nlevels-1); - mg_matrices_renumbered.clear (); - mg_sparsity.resize(0, nlevels-1); - mg_sparsity_renumbered.resize(0, nlevels-1); - - for (unsigned int level=0;level -void LaplaceProblem::test_interface_dofs () -{ - std::vector > interface_dofs; - std::vector > boundary_interface_dofs; - for (unsigned int level = 0; level tmp (mg_dof_handler.n_dofs(level)); - interface_dofs.push_back (tmp); - boundary_interface_dofs.push_back (tmp); - } - - MGTools::extract_inner_interface_dofs(mg_dof_handler_renumbered, - interface_dofs, boundary_interface_dofs); - deallog << "1. Test" << std::endl; - print(mg_dof_handler_renumbered, interface_dofs); - - deallog << std::endl; - - MGTools::extract_inner_interface_dofs(mg_dof_handler, - interface_dofs, boundary_interface_dofs); - - deallog << "2. Test" << std::endl; - print(mg_dof_handler, interface_dofs); -} - - template -void LaplaceProblem::test_renumber () +template +void +LaplaceProblem::output_gpl(const MGDoFHandler &dof, + MGLevelObject > &v, + const bool renumbered) { - std::vector > v(triangulation.n_levels()); - for(unsigned int l=0; l integration_worker; + MeshWorker::Assembler::GnuplotPatch assembler; + + const unsigned int n_gauss_points = dof.get_fe().tensor_degree()+1; + QTrapez<1> trapez; + QIterated quadrature(trapez, n_gauss_points); + integration_worker.cell_quadrature = quadrature; + UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients; + integration_worker.add_update_flags(update_flags, true, true, true, true); + + NamedData >* > data; + data.add(&v, "mg_vector"); + MeshWorker::VectorSelector cs; + cs.add("mg_vector"); + integration_worker.cell_selector = cs; + + assembler.initialize(dim, quadrature.size(), dim+mg_dof_handler.get_fe().n_components()); + MeshWorker::IntegrationInfoBox info_box; + MeshWorker::DoFInfo dof_info(dof); + info_box.initialize(integration_worker, fe, mapping, data); - const unsigned int dofs_per_cell = - mg_dof_handler_renumbered.get_fe().dofs_per_cell; - std::vector dof_indices(dofs_per_cell); - for (typename MGDoFHandler::cell_iterator - cell = mg_dof_handler_renumbered.begin(); - cell != mg_dof_handler_renumbered.end(); ++cell) - { - cell->get_mg_dof_indices(dof_indices); - for(unsigned int i=0; ilevel()][dof_indices[i]] = true; - } for(unsigned int l=0; l matrix_integrator; + MeshWorker::loop, dim, dim> ( + dof.begin(l), dof.end(l), + dof_info, info_box, + std_cxx1x::bind(&OutputCreator::cell, &matrix_integrator, _1, _2), + 0, + 0, + assembler); } } + + template -void LaplaceProblem::test () +void LaplaceProblem::test () { typename FunctionMap::type dirichlet_boundary; ZeroFunction dirichlet_bc(fe.n_components()); @@ -310,60 +311,39 @@ void LaplaceProblem::test () boundary_indices); MGTools::make_boundary_list (mg_dof_handler_renumbered, dirichlet_boundary, boundary_indices_renumbered); - const unsigned int min_l = mg_matrices.get_minlevel(); - const unsigned int max_l = mg_matrices.get_maxlevel(); - for(unsigned int l=min_l; l vector_memory; - GrowingVectorMemory<> vector_memory_renumbered; - MGTransferPrebuilt > mg_transfer; mg_transfer.build_matrices(mg_dof_handler, boundary_indices); MGTransferPrebuilt > mg_transfer_renumbered; mg_transfer_renumbered.build_matrices(mg_dof_handler_renumbered, boundary_indices_renumbered); - Vector u(mg_dof_handler.n_dofs()); - initialize(mg_dof_handler,u); - deallog << "initialize " << u.l2_norm() << std::endl; + Vector test; + test.reinit(mg_dof_handler.n_dofs()); MGLevelObject > v(0, triangulation.n_levels()-1); - reinit_vector(mg_dof_handler, v); - - mg_transfer.copy_to_mg(mg_dof_handler, v, u); - //print(mg_dof_handler, v); - - u=0; - initialize(mg_dof_handler, v); - for(unsigned int l=0; l > u(0, triangulation.n_levels()-1); + MGLevelObject > d(0, triangulation.n_levels()-1); - mg_transfer.copy_from_mg(mg_dof_handler, u, v); - Vector diff = u; + initialize(mg_dof_handler, test); + mg_transfer.copy_to_mg(mg_dof_handler, v, test); + //output_gpl(mg_dof_handler, v, false); - initialize(mg_dof_handler_renumbered,u); - deallog << "initialize " << u.l2_norm() << std::endl; - - reinit_vector(mg_dof_handler_renumbered, v); - mg_transfer_renumbered.copy_to_mg(mg_dof_handler_renumbered, v, u); - deallog << "copy_to_mg" << std::endl; - //print(mg_dof_handler_renumbered, v); - - u=0; - initialize(mg_dof_handler_renumbered, v); - for(unsigned int l=0; l void LaplaceProblem::refine_local () { @@ -401,9 +381,9 @@ void LaplaceProblem::refine_local () } template -void LaplaceProblem::run () +void LaplaceProblem::run () { - for (unsigned int cycle=0; cycle<7; ++cycle) + for (unsigned int cycle=0; cycle<6; ++cycle) { deallog << "Cycle " << cycle << std::endl; @@ -419,7 +399,7 @@ void LaplaceProblem::run () }; } -int main () +int main () { std::ofstream logfile("mg_renumbered_03/output"); deallog << std::setprecision(4); @@ -429,5 +409,4 @@ int main () LaplaceProblem<2> laplace_problem_2d(1); laplace_problem_2d.run (); - } -- 2.39.5