From: bangerth Date: Fri, 18 Feb 2011 14:11:54 +0000 (+0000) Subject: Convert to the hp framework. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=71f7f980add60ea32d47382f964d45f80a933459;p=dealii-svn.git Convert to the hp framework. git-svn-id: https://svn.dealii.org/trunk@23392 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-46/step-46.cc b/deal.II/examples/step-46/step-46.cc index b9bae96983..bb4bfc6ee5 100644 --- a/deal.II/examples/step-46/step-46.cc +++ b/deal.II/examples/step-46/step-46.cc @@ -31,15 +31,17 @@ #include #include -#include #include -#include #include +#include #include #include #include -#include + +#include +#include +#include #include #include @@ -56,7 +58,7 @@ template class StokesProblem { public: - StokesProblem (const unsigned int degree); + StokesProblem (const unsigned int stokes_degree); void run (); private: @@ -66,16 +68,17 @@ class StokesProblem void output_results (const unsigned int refinement_cycle) const; void refine_mesh (); - const unsigned int degree; + const unsigned int stokes_degree; - Triangulation triangulation; - FESystem fe; - DoFHandler dof_handler; + Triangulation triangulation; + FESystem stokes_fe; + hp::FECollection fe_collection; + hp::DoFHandler dof_handler; - ConstraintMatrix constraints; + ConstraintMatrix constraints; - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; Vector solution; Vector system_rhs; @@ -106,7 +109,7 @@ BoundaryValues::value (const Point &p, ExcIndexRange (component, 0, this->n_components)); if (component == 1) - return std::sin(2*numbers::PI*p[0]); + return std::sin(numbers::PI*p[0]); return 0; } @@ -163,14 +166,16 @@ RightHandSide::vector_value (const Point &p, template -StokesProblem::StokesProblem (const unsigned int degree) +StokesProblem::StokesProblem (const unsigned int stokes_degree) : - degree (degree), + stokes_degree (stokes_degree), triangulation (Triangulation::maximum_smoothing), - fe (FE_Q(degree+1), dim, - FE_Q(degree), 1), + stokes_fe (FE_Q(stokes_degree+1), dim, + FE_Q(stokes_degree), 1), dof_handler (triangulation) -{} +{ + fe_collection.push_back (stokes_fe); +} @@ -180,7 +185,7 @@ void StokesProblem::setup_dofs () { system_matrix.clear (); - dof_handler.distribute_dofs (fe); + dof_handler.distribute_dofs (fe_collection); DoFRenumbering::Cuthill_McKee (dof_handler); std::vector block_component (dim+1,0); @@ -200,19 +205,18 @@ void StokesProblem::setup_dofs () component_mask); } + for (unsigned int i=0; i dofs_per_block (2); - DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component); - const unsigned int n_u = dofs_per_block[0], - n_p = dofs_per_block[1]; - std::cout << " Number of active cells: " << triangulation.n_active_cells() << std::endl << " Number of degrees of freedom: " << dof_handler.n_dofs() - << " (" << n_u << '+' << n_p << ')' << std::endl; { @@ -237,17 +241,19 @@ void StokesProblem::assemble_system () system_matrix=0; system_rhs=0; - QGauss quadrature_formula(degree+2); - - FEValues fe_values (fe, quadrature_formula, - update_values | - update_quadrature_points | - update_JxW_values | - update_gradients); + QGauss stokes_quadrature(stokes_degree+2); + hp::QCollection q_collection; + q_collection.push_back (stokes_quadrature); + + hp::FEValues hp_fe_values (fe_collection, q_collection, + update_values | + update_quadrature_points | + update_JxW_values | + update_gradients); - const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int dofs_per_cell = stokes_fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); + const unsigned int n_q_points = stokes_quadrature.size(); FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); Vector local_rhs (dofs_per_cell); @@ -255,8 +261,6 @@ void StokesProblem::assemble_system () std::vector local_dof_indices (dofs_per_cell); const RightHandSide right_hand_side; - std::vector > rhs_values (n_q_points, - Vector(dim+1)); const FEValuesExtractors::Vector velocities (0); const FEValuesExtractors::Scalar pressure (dim); @@ -265,18 +269,18 @@ void StokesProblem::assemble_system () std::vector div_phi_u (dofs_per_cell); std::vector phi_p (dofs_per_cell); - typename DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) { - fe_values.reinit (cell); + hp_fe_values.reinit (cell); + + const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + local_matrix = 0; local_rhs = 0; - right_hand_side.vector_value_list(fe_values.get_quadrature_points(), - rhs_values); - for (unsigned int q=0; q::assemble_system () } for (unsigned int i=0; i::assemble_system () local_matrix(i,j) = local_matrix(j,i); cell->get_dof_indices (local_dof_indices); - constraints.distribute_local_to_global (local_matrix, local_rhs, + constraints.distribute_local_to_global (local_matrix, local_dof_indices, - system_matrix, system_rhs); + system_matrix); } } @@ -346,12 +339,13 @@ StokesProblem::output_results (const unsigned int refinement_cycle) const data_component_interpretation .push_back (DataComponentInterpretation::component_is_scalar); - DataOut data_out; + DataOut > data_out; data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, solution_names, - DataOut::type_dof_data, + DataOut >::type_dof_data, data_component_interpretation); - data_out.build_patches (); + data_out.build_patches (stokes_fe.degree); std::ostringstream filename; filename << "solution-" @@ -373,7 +367,7 @@ StokesProblem::refine_mesh () std::vector component_mask (dim+1, false); component_mask[dim] = true; KellyErrorEstimator::estimate (dof_handler, - QGauss(degree+1), + QGauss(stokes_degree+1), typename FunctionMap::type(), solution, estimated_error_per_cell, @@ -399,7 +393,7 @@ void StokesProblem::run () if (cell->face(f)->center()[dim-1] == 1) cell->face(f)->set_all_boundary_indicators(1); - triangulation.refine_global (4-dim); + triangulation.refine_global (3-dim); for (unsigned int refinement_cycle = 0; refinement_cycle<6; ++refinement_cycle)