From: janssen Date: Wed, 21 Aug 2013 23:07:32 +0000 (+0000) Subject: one step further X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e46c6c742ca67936384604ade5da339101262371;p=dealii-svn.git one step further git-svn-id: https://svn.dealii.org/trunk@30395 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/multigrid/step-16-02.cc b/tests/multigrid/step-16-02.cc index 34a34c8565..6c1b72b426 100644 --- a/tests/multigrid/step-16-02.cc +++ b/tests/multigrid/step-16-02.cc @@ -106,8 +106,8 @@ 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)); +// Laplace::nitsche_matrix(dinfo.matrix(0,false).matrix, info.fe_values(0), +// Laplace::compute_penalty(dinfo, dinfo, deg, deg)); } @@ -117,10 +117,10 @@ void LaplaceMatrix::face( 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)); +// 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 @@ -135,7 +135,7 @@ private: void assemble_system (); void assemble_multigrid (const bool &use_mw); void solve (); - void refine_grid (); + void refine_grid (const std::string& reftype); void output_results (const unsigned int cycle) const; Triangulation triangulation; @@ -182,9 +182,9 @@ template double Coefficient::value (const Point &p, const unsigned int) const { - if (p.square() < 0.5*0.5) - return 20; - else +// if (p.square() < 0.5*0.5) +// return 20; +// else return 1; } @@ -354,16 +354,14 @@ void LaplaceProblem::assemble_multigrid (const bool& use_mw) 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()); + info_box.initialize(fe, mapping); - MeshWorker::DoFInfo dof_info(mg_dof_handler);//.block_info()); + MeshWorker::DoFInfo dof_info(mg_dof_handler); 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(), @@ -455,10 +453,11 @@ void LaplaceProblem::assemble_multigrid (const bool& use_mw) const unsigned int nlevels = triangulation.n_levels(); for (unsigned int level=1;level::solve () MGMatrix<> mg_matrix(&mg_matrices); MGMatrix<> mg_interface_up(&mg_interface_in); - MGMatrix<> mg_interface_down(&mg_interface_out); + MGMatrix<> mg_interface_down(&mg_interface_in); Multigrid > mg(mg_dof_handler, mg_matrix, @@ -566,19 +565,37 @@ void LaplaceProblem::solve () // today compared to those that were // available when step-6 was written. template -void LaplaceProblem::refine_grid () +void LaplaceProblem::refine_grid (const std::string& reftype) { - Vector estimated_error_per_cell (triangulation.n_active_cells()); - - KellyErrorEstimator::estimate (static_cast&>(mg_dof_handler), - QGauss(3), - typename FunctionMap::type(), - solution, - estimated_error_per_cell); - GridRefinement::refine_and_coarsen_fixed_number (triangulation, - estimated_error_per_cell, - 0.3, 0.03); - triangulation.execute_coarsening_and_refinement (); + bool cell_refined = false; + if (reftype == "center" || !cell_refined) + { + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + for (unsigned int vertex=0; + vertex < GeometryInfo::vertices_per_cell; + ++vertex) + { + { + const Point p = cell->vertex(vertex); + const Point origin = (dim == 2 ? + Point(0,0) : + Point(0,0,0)); + const double dist = p.distance(origin); + if(dist<0.25/M_PI) + { + cell->set_refine_flag (); + cell_refined = true; + break; + } + } + } + } + if (reftype=="global" || !cell_refined) + triangulation.refine_global(1); + else + triangulation.execute_coarsening_and_refinement (); } @@ -597,8 +614,8 @@ void LaplaceProblem::output_results (const unsigned int cycle) const << cycle << ".vtk"; -// std::ofstream output (filename.str().c_str()); -// data_out.write_vtk (output); + std::ofstream output (filename.str().c_str()); + data_out.write_vtk (output); } @@ -621,15 +638,15 @@ void LaplaceProblem::run () if (cycle == 0) { - GridGenerator::hyper_ball (triangulation); + GridGenerator::hyper_cube (triangulation, -1, 1); - static const HyperBallBoundary boundary; - triangulation.set_boundary (0, boundary); +// static const HyperBallBoundary boundary; +// triangulation.set_boundary (0, boundary); triangulation.refine_global (1); } else - refine_grid (); + refine_grid ("center"); deallog << " Number of active cells: "