From 2deba76953a6cf91593afff6dea66db6c106ddfa Mon Sep 17 00:00:00 2001 From: bangerth Date: Mon, 31 Jul 2006 19:41:45 +0000 Subject: [PATCH] Partly hp-ify these tests. Not all library calls work right now git-svn-id: https://svn.dealii.org/trunk@13548 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/hp/step-11.cc | 14 +-- tests/hp/step-12.cc | 212 +++++++++++++++++----------------- tests/hp/step-13.cc | 105 ++++++++--------- tests/hp/step-14.cc | 276 ++++++++++++++++++++++---------------------- tests/hp/step-15.cc | 106 ++++++++--------- 5 files changed, 360 insertions(+), 353 deletions(-) diff --git a/tests/hp/step-11.cc b/tests/hp/step-11.cc index c2d88a1ea4..c1247d8948 100644 --- a/tests/hp/step-11.cc +++ b/tests/hp/step-11.cc @@ -37,7 +37,7 @@ std::ofstream logfile("step-11/output"); #include #include #include -#include +#include #include #include #include @@ -63,9 +63,9 @@ class LaplaceProblem void solve (); Triangulation triangulation; - FE_Q fe; - DoFHandler dof_handler; - MappingQ mapping; + hp::FECollection fe; + hp::DoFHandler dof_handler; + hp::MappingCollection mapping; SparsityPattern sparsity_pattern; SparseMatrix system_matrix; @@ -81,9 +81,9 @@ class LaplaceProblem template LaplaceProblem::LaplaceProblem (const unsigned int mapping_degree) : - fe (1), + fe (FE_Q(1)), dof_handler (triangulation), - mapping (mapping_degree) + mapping (MappingQ(mapping_degree)) { deallog << "Using mapping with degree " << mapping_degree << ":" << std::endl @@ -134,7 +134,7 @@ void LaplaceProblem::assemble_and_solve () { const unsigned int gauss_degree - = std::max (static_cast(std::ceil(1.*(mapping.get_degree()+1)/2)), + = std::max (static_cast(std::ceil(1.*(static_cast&>(mapping[0]).get_degree()+1)/2)), 2U); MatrixTools::create_laplace_matrix (mapping, dof_handler, QGauss(gauss_degree), diff --git a/tests/hp/step-12.cc b/tests/hp/step-12.cc index bb4b20849b..42b1e0fbec 100644 --- a/tests/hp/step-12.cc +++ b/tests/hp/step-12.cc @@ -29,7 +29,7 @@ std::ofstream logfile("step-12/output"); #include #include #include -#include +#include #include #include #include @@ -132,21 +132,23 @@ class DGTransportEquation public: DGTransportEquation(); - void assemble_cell_term(const FEValues& fe_v, + void assemble_cell_term(const hp::FEValues& fe_v, FullMatrix &ui_vi_matrix, Vector &cell_vector) const; - void assemble_boundary_term(const FEFaceValues& fe_v, + void assemble_boundary_term(const hp::FEFaceValues& fe_v, FullMatrix &ui_vi_matrix, Vector &cell_vector) const; - - void assemble_face_term1(const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, + + template + void assemble_face_term1(const X& fe_v, + const Y& fe_v_neighbor, FullMatrix &ui_vi_matrix, FullMatrix &ue_vi_matrix) const; - void assemble_face_term2(const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, + template + void assemble_face_term2(const X& fe_v, + const Y& fe_v_neighbor, FullMatrix &ui_vi_matrix, FullMatrix &ue_vi_matrix, FullMatrix &ui_ve_matrix, @@ -169,149 +171,151 @@ DGTransportEquation::DGTransportEquation () template void DGTransportEquation::assemble_cell_term( - const FEValues &fe_v, + const hp::FEValues &fe_v, FullMatrix &ui_vi_matrix, Vector &cell_vector) const { - const std::vector &JxW = fe_v.get_JxW_values (); + const std::vector &JxW = fe_v.get_present_fe_values().get_JxW_values (); - std::vector > beta (fe_v.n_quadrature_points); - std::vector rhs (fe_v.n_quadrature_points); + std::vector > beta (fe_v.get_present_fe_values().n_quadrature_points); + std::vector rhs (fe_v.get_present_fe_values().n_quadrature_points); - beta_function.value_list (fe_v.get_quadrature_points(), beta); - rhs_function.value_list (fe_v.get_quadrature_points(), rhs); + beta_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), beta); + rhs_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), rhs); - for (unsigned int point=0; point void DGTransportEquation::assemble_boundary_term( - const FEFaceValues& fe_v, + const hp::FEFaceValues& fe_v, FullMatrix &ui_vi_matrix, Vector &cell_vector) const { - const std::vector &JxW = fe_v.get_JxW_values (); - const std::vector > &normals = fe_v.get_normal_vectors (); + const std::vector &JxW = fe_v.get_present_fe_values().get_JxW_values (); + const std::vector > &normals = fe_v.get_present_fe_values().get_normal_vectors (); - std::vector > beta (fe_v.n_quadrature_points); - std::vector g(fe_v.n_quadrature_points); + std::vector > beta (fe_v.get_present_fe_values().n_quadrature_points); + std::vector g(fe_v.get_present_fe_values().n_quadrature_points); - beta_function.value_list (fe_v.get_quadrature_points(), beta); - boundary_function.value_list (fe_v.get_quadrature_points(), g); + beta_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), beta); + boundary_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), g); - for (unsigned int point=0; point0) - for (unsigned int i=0; i +template void DGTransportEquation::assemble_face_term1( - const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, + const X& fe_v, + const Y& fe_v_neighbor, FullMatrix &ui_vi_matrix, FullMatrix &ue_vi_matrix) const { - const std::vector &JxW = fe_v.get_JxW_values (); - const std::vector > &normals = fe_v.get_normal_vectors (); + const std::vector &JxW = fe_v.get_present_fe_values().get_JxW_values (); + const std::vector > &normals = fe_v.get_present_fe_values().get_normal_vectors (); - std::vector > beta (fe_v.n_quadrature_points); - beta_function.value_list (fe_v.get_quadrature_points(), beta); + std::vector > beta (fe_v.get_present_fe_values().n_quadrature_points); + beta_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), beta); - for (unsigned int point=0; point0) - for (unsigned int i=0; i +template void DGTransportEquation::assemble_face_term2( - const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, + const X& fe_v, + const Y& fe_v_neighbor, FullMatrix &ui_vi_matrix, FullMatrix &ue_vi_matrix, FullMatrix &ui_ve_matrix, FullMatrix &ue_ve_matrix) const { - const std::vector &JxW = fe_v.get_JxW_values (); - const std::vector > &normals = fe_v.get_normal_vectors (); + const std::vector &JxW = fe_v.get_present_fe_values().get_JxW_values (); + const std::vector > &normals = fe_v.get_present_fe_values().get_normal_vectors (); - std::vector > beta (fe_v.n_quadrature_points); + std::vector > beta (fe_v.get_present_fe_values().n_quadrature_points); - beta_function.value_list (fe_v.get_quadrature_points(), beta); + beta_function.value_list (fe_v.get_present_fe_values().get_quadrature_points(), beta); - for (unsigned int point=0; point0) { - for (unsigned int i=0; i triangulation; - const MappingQ1 mapping; + const hp::MappingCollection mapping; - FE_DGQ fe; - DoFHandler dof_handler; + hp::FECollection fe; + hp::DoFHandler dof_handler; SparsityPattern sparsity_pattern; SparseMatrix system_matrix; - const QGauss quadrature; - const QGauss face_quadrature; + const hp::QCollection quadrature; + const hp::QCollection face_quadrature; Vector solution1; Vector solution2; @@ -358,11 +362,11 @@ class DGMethod template DGMethod::DGMethod () : - mapping (), - fe (1), + mapping (MappingQ1()), + fe (FE_DGQ(1)), dof_handler (triangulation), - quadrature (4), - face_quadrature (4), + quadrature (QGauss(4)), + face_quadrature (QGauss(4)), dg () {} @@ -382,7 +386,7 @@ void DGMethod::setup_system () sparsity_pattern.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), (GeometryInfo::faces_per_cell - *GeometryInfo::subfaces_per_face+1)*fe.dofs_per_cell); + *GeometryInfo::subfaces_per_face+1)*fe[0].dofs_per_cell); DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern); @@ -399,7 +403,7 @@ void DGMethod::setup_system () template void DGMethod::assemble_system1 () { - const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; + const unsigned int dofs_per_cell = dof_handler.get_fe()[0].dofs_per_cell; std::vector dofs (dofs_per_cell); std::vector dofs_neighbor (dofs_per_cell); @@ -415,16 +419,16 @@ void DGMethod::assemble_system1 () const UpdateFlags neighbor_face_update_flags = update_values; - FEValues fe_v ( + hp::FEValues fe_v ( mapping, fe, quadrature, update_flags); - FEFaceValues fe_v_face ( + hp::FEFaceValues fe_v_face ( mapping, fe, face_quadrature, face_update_flags); - FESubfaceValues fe_v_subface ( + hp::FESubfaceValues fe_v_subface ( mapping, fe, face_quadrature, face_update_flags); - FEFaceValues fe_v_face_neighbor ( + hp::FEFaceValues fe_v_face_neighbor ( mapping, fe, face_quadrature, neighbor_face_update_flags); - FESubfaceValues fe_v_subface_neighbor ( + hp::FESubfaceValues fe_v_subface_neighbor ( mapping, fe, face_quadrature, neighbor_face_update_flags); FullMatrix ui_vi_matrix (dofs_per_cell, dofs_per_cell); @@ -432,7 +436,7 @@ void DGMethod::assemble_system1 () Vector cell_vector (dofs_per_cell); - typename DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@ -451,7 +455,7 @@ void DGMethod::assemble_system1 () for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) { - typename DoFHandler::face_iterator face=cell->face(face_no); + typename hp::DoFHandler::face_iterator face=cell->face(face_no); ue_vi_matrix = 0; @@ -465,7 +469,7 @@ void DGMethod::assemble_system1 () } else { - typename DoFHandler::cell_iterator neighbor= + typename hp::DoFHandler::cell_iterator neighbor= cell->neighbor(face_no);; if (face->has_children()) @@ -477,7 +481,7 @@ void DGMethod::assemble_system1 () for (unsigned int subface_no=0; subface_non_children(); ++subface_no) { - typename DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator neighbor_child = cell->neighbor_child_on_subface (face_no, subface_no); @@ -565,7 +569,7 @@ void DGMethod::assemble_system1 () template void DGMethod::assemble_system2 () { - const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; + const unsigned int dofs_per_cell = dof_handler.get_fe()[0].dofs_per_cell; std::vector dofs (dofs_per_cell); std::vector dofs_neighbor (dofs_per_cell); @@ -581,13 +585,13 @@ void DGMethod::assemble_system2 () const UpdateFlags neighbor_face_update_flags = update_values; - FEValues fe_v ( + hp::FEValues fe_v ( mapping, fe, quadrature, update_flags); - FEFaceValues fe_v_face ( + hp::FEFaceValues fe_v_face ( mapping, fe, face_quadrature, face_update_flags); - FESubfaceValues fe_v_subface ( + hp::FESubfaceValues fe_v_subface ( mapping, fe, face_quadrature, face_update_flags); - FEFaceValues fe_v_face_neighbor ( + hp::FEFaceValues fe_v_face_neighbor ( mapping, fe, face_quadrature, neighbor_face_update_flags); @@ -599,7 +603,7 @@ void DGMethod::assemble_system2 () Vector cell_vector (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) @@ -617,7 +621,7 @@ void DGMethod::assemble_system2 () for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) { - typename DoFHandler::face_iterator face= + typename hp::DoFHandler::face_iterator face= cell->face(face_no); if (face->at_boundary()) @@ -632,7 +636,7 @@ void DGMethod::assemble_system2 () { Assert (cell->neighbor(face_no).state() == IteratorState::valid, ExcInternalError()); - typename DoFHandler::cell_iterator neighbor= + typename hp::DoFHandler::cell_iterator neighbor= cell->neighbor(face_no); if (face->has_children()) { @@ -642,7 +646,7 @@ void DGMethod::assemble_system2 () for (unsigned int subface_no=0; subface_non_children(); ++subface_no) { - typename DoFHandler::cell_iterator neighbor_child + typename hp::DoFHandler::cell_iterator neighbor_child = cell->neighbor_child_on_subface (face_no, subface_no); Assert (neighbor_child->face(neighbor2) == face->child(subface_no), ExcInternalError()); @@ -733,7 +737,7 @@ void DGMethod::solve (Vector &solution) PreconditionBlockSSOR > preconditioner; - preconditioner.initialize(system_matrix, fe.dofs_per_cell); + preconditioner.initialize(system_matrix, fe[0].dofs_per_cell); solver.solve (system_matrix, solution, right_hand_side, preconditioner); @@ -750,7 +754,7 @@ void DGMethod::refine_grid () solution2, gradient_indicator); - typename DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) @@ -785,7 +789,7 @@ void DGMethod::output_results (const unsigned int cycle) const deallog << "Writing solution to <" << filename << ">..." << std::endl << std::endl; - DataOut data_out; + DataOut > data_out; data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution2, "u"); diff --git a/tests/hp/step-13.cc b/tests/hp/step-13.cc index 2d80eb1603..8b8aea97c2 100644 --- a/tests/hp/step-13.cc +++ b/tests/hp/step-13.cc @@ -40,7 +40,7 @@ std::ofstream logfile("step-13/output"); #include #include #include -#include +#include #include #include #include @@ -64,7 +64,7 @@ namespace Evaluation void set_refinement_cycle (const unsigned int refinement_cycle); - virtual void operator () (const DoFHandler &dof_handler, + virtual void operator () (const hp::DoFHandler &dof_handler, const Vector &solution) const = 0; protected: unsigned int refinement_cycle; @@ -93,7 +93,7 @@ namespace Evaluation PointValueEvaluation (const Point &evaluation_point, TableHandler &results_table); - virtual void operator () (const DoFHandler &dof_handler, + virtual void operator () (const hp::DoFHandler &dof_handler, const Vector &solution) const; DeclException1 (ExcEvaluationPointNotFound, @@ -120,12 +120,12 @@ namespace Evaluation template void PointValueEvaluation:: - operator () (const DoFHandler &dof_handler, + operator () (const hp::DoFHandler &dof_handler, const Vector &solution) const { double point_value = 1e20; - typename DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); bool evaluation_point_found = false; @@ -135,7 +135,8 @@ namespace Evaluation ++vertex) if (cell->vertex(vertex) == evaluation_point) { - point_value = solution(cell->vertex_dof_index(vertex,0)); + point_value = solution(cell->vertex_dof_index(vertex,0, + cell->active_fe_index())); evaluation_point_found = true; break; @@ -159,7 +160,7 @@ namespace Evaluation SolutionOutput (const std::string &output_name_base, const typename DataOut::OutputFormat output_format); - virtual void operator () (const DoFHandler &dof_handler, + virtual void operator () (const hp::DoFHandler &dof_handler, const Vector &solution) const; private: const std::string output_name_base; @@ -179,10 +180,10 @@ namespace Evaluation template void - SolutionOutput::operator () (const DoFHandler &dof_handler, + SolutionOutput::operator () (const hp::DoFHandler &dof_handler, const Vector &solution) const { - DataOut data_out; + DataOut > data_out; data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, "solution"); data_out.build_patches (); @@ -241,8 +242,8 @@ namespace LaplaceSolver { public: Solver (Triangulation &triangulation, - const FiniteElement &fe, - const Quadrature &quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, const Function &boundary_values); virtual ~Solver (); @@ -260,9 +261,9 @@ namespace LaplaceSolver n_dofs () const; protected: - const SmartPointer > fe; - const SmartPointer > quadrature; - DoFHandler dof_handler; + const SmartPointer > fe; + const SmartPointer > quadrature; + hp::DoFHandler dof_handler; Vector solution; const SmartPointer > boundary_values; @@ -271,7 +272,7 @@ namespace LaplaceSolver private: struct LinearSystem { - LinearSystem (const DoFHandler &dof_handler); + LinearSystem (const hp::DoFHandler &dof_handler); void solve (Vector &solution) const; @@ -286,8 +287,8 @@ namespace LaplaceSolver void assemble_matrix (LinearSystem &linear_system, - const typename DoFHandler::active_cell_iterator &begin_cell, - const typename DoFHandler::active_cell_iterator &end_cell, + const typename hp::DoFHandler::active_cell_iterator &begin_cell, + const typename hp::DoFHandler::active_cell_iterator &end_cell, Threads::ThreadMutex &mutex) const; }; @@ -295,8 +296,8 @@ namespace LaplaceSolver template Solver::Solver (Triangulation &triangulation, - const FiniteElement &fe, - const Quadrature &quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, const Function &boundary_values) : Base (triangulation), @@ -349,7 +350,7 @@ namespace LaplaceSolver Solver::assemble_linear_system (LinearSystem &linear_system) { typedef - typename DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator active_cell_iterator; const unsigned int n_threads = multithread_info.n_default_threads; @@ -392,21 +393,21 @@ namespace LaplaceSolver template void Solver::assemble_matrix (LinearSystem &linear_system, - const typename DoFHandler::active_cell_iterator &begin_cell, - const typename DoFHandler::active_cell_iterator &end_cell, + const typename hp::DoFHandler::active_cell_iterator &begin_cell, + const typename hp::DoFHandler::active_cell_iterator &end_cell, Threads::ThreadMutex &mutex) const { - FEValues fe_values (*fe, *quadrature, + hp::FEValues fe_values (*fe, *quadrature, update_gradients | update_JxW_values); - const unsigned int dofs_per_cell = fe->dofs_per_cell; - const unsigned int n_q_points = quadrature->n_quadrature_points; + const unsigned int dofs_per_cell = (*fe)[0].dofs_per_cell; + const unsigned int n_q_points = (*quadrature)[0].n_quadrature_points; FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); std::vector local_dof_indices (dofs_per_cell); - for (typename DoFHandler::active_cell_iterator cell=begin_cell; + for (typename hp::DoFHandler::active_cell_iterator cell=begin_cell; cell!=end_cell; ++cell) { cell_matrix = 0; @@ -416,9 +417,9 @@ namespace LaplaceSolver for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); @@ -435,11 +436,11 @@ namespace LaplaceSolver template Solver::LinearSystem:: - LinearSystem (const DoFHandler &dof_handler) + LinearSystem (const hp::DoFHandler &dof_handler) { hanging_node_constraints.clear (); - void (*mhnc_p) (const DoFHandler &, + void (*mhnc_p) (const hp::DoFHandler &, ConstraintMatrix &) = &DoFTools::make_hanging_node_constraints; @@ -487,8 +488,8 @@ namespace LaplaceSolver { public: PrimalSolver (Triangulation &triangulation, - const FiniteElement &fe, - const Quadrature &quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, const Function &rhs_function, const Function &boundary_values); protected: @@ -500,8 +501,8 @@ namespace LaplaceSolver template PrimalSolver:: PrimalSolver (Triangulation &triangulation, - const FiniteElement &fe, - const Quadrature &quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, const Function &rhs_function, const Function &boundary_values) : @@ -518,32 +519,32 @@ namespace LaplaceSolver PrimalSolver:: assemble_rhs (Vector &rhs) const { - FEValues fe_values (*this->fe, *this->quadrature, + hp::FEValues fe_values (*this->fe, *this->quadrature, update_values | update_q_points | update_JxW_values); - const unsigned int dofs_per_cell = this->fe->dofs_per_cell; - const unsigned int n_q_points = this->quadrature->n_quadrature_points; + const unsigned int dofs_per_cell = (*this->fe)[0].dofs_per_cell; + const unsigned int n_q_points = (*this->quadrature)[0].n_quadrature_points; Vector cell_rhs (dofs_per_cell); std::vector rhs_values (n_q_points); std::vector local_dof_indices (dofs_per_cell); - typename DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator cell = this->dof_handler.begin_active(), endc = this->dof_handler.end(); for (; cell!=endc; ++cell) { cell_rhs = 0; fe_values.reinit (cell); - rhs_function->value_list (fe_values.get_quadrature_points(), + rhs_function->value_list (fe_values.get_present_fe_values().get_quadrature_points(), rhs_values); for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); for (unsigned int i=0; i &coarse_grid, - const FiniteElement &fe, - const Quadrature &quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, const Function &rhs_function, const Function &boundary_values); @@ -571,8 +572,8 @@ namespace LaplaceSolver template RefinementGlobal:: RefinementGlobal (Triangulation &coarse_grid, - const FiniteElement &fe, - const Quadrature &quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, const Function &rhs_function, const Function &boundary_values) : @@ -597,8 +598,8 @@ namespace LaplaceSolver { public: RefinementKelly (Triangulation &coarse_grid, - const FiniteElement &fe, - const Quadrature &quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, const Function &rhs_function, const Function &boundary_values); @@ -610,8 +611,8 @@ namespace LaplaceSolver template RefinementKelly:: RefinementKelly (Triangulation &coarse_grid, - const FiniteElement &fe, - const Quadrature &quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, const Function &rhs_function, const Function &boundary_values) : @@ -752,8 +753,8 @@ void solve_problem (const std::string &solver_name) Triangulation triangulation; GridGenerator::hyper_cube (triangulation, -1, 1); triangulation.refine_global (2); - const FE_Q fe(1); - const QGauss quadrature(4); + const hp::FECollection fe(FE_Q(1)); + const hp::QCollection quadrature(QGauss(4)); const RightHandSide rhs_function; const Solution boundary_values; diff --git a/tests/hp/step-14.cc b/tests/hp/step-14.cc index 725cb2cbe2..3c99840ca2 100644 --- a/tests/hp/step-14.cc +++ b/tests/hp/step-14.cc @@ -40,7 +40,7 @@ std::ofstream logfile("step-14/output"); #include #include #include -#include +#include #include #include #include @@ -65,7 +65,7 @@ namespace Evaluation void set_refinement_cycle (const unsigned int refinement_cycle); - virtual void operator () (const DoFHandler &dof_handler, + virtual void operator () (const hp::DoFHandler &dof_handler, const Vector &solution) const = 0; protected: unsigned int refinement_cycle; @@ -92,7 +92,7 @@ namespace Evaluation public: PointValueEvaluation (const Point &evaluation_point); - virtual void operator () (const DoFHandler &dof_handler, + virtual void operator () (const hp::DoFHandler &dof_handler, const Vector &solution) const; DeclException1 (ExcEvaluationPointNotFound, @@ -116,12 +116,12 @@ namespace Evaluation template void PointValueEvaluation:: - operator () (const DoFHandler &dof_handler, + operator () (const hp::DoFHandler &dof_handler, const Vector &solution) const { double point_value = 1e20; - typename DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); bool evaluation_point_found = false; @@ -154,7 +154,7 @@ namespace Evaluation public: PointXDerivativeEvaluation (const Point &evaluation_point); - virtual void operator () (const DoFHandler &dof_handler, + virtual void operator () (const hp::DoFHandler &dof_handler, const Vector &solution) const; DeclException1 (ExcEvaluationPointNotFound, @@ -177,19 +177,19 @@ namespace Evaluation template void PointXDerivativeEvaluation:: - operator () (const DoFHandler &dof_handler, + operator () (const hp::DoFHandler &dof_handler, const Vector &solution) const { double point_derivative = 0; QTrapez vertex_quadrature; - FEValues fe_values (dof_handler.get_fe(), + hp::FEValues fe_values (dof_handler.get_fe(), vertex_quadrature, update_gradients | update_q_points); std::vector > solution_gradients (vertex_quadrature.n_quadrature_points); - typename DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); unsigned int evaluation_point_hits = 0; @@ -200,12 +200,12 @@ namespace Evaluation if (cell->vertex(vertex) == evaluation_point) { fe_values.reinit (cell); - fe_values.get_function_grads (solution, + fe_values.get_present_fe_values().get_function_grads (solution, solution_gradients); unsigned int q_point = 0; for (; q_point &dof_handler, + virtual void operator () (const hp::DoFHandler &dof_handler, const Vector &solution) const; private: const std::string output_name_base; @@ -251,7 +251,7 @@ namespace Evaluation template void - GridOutput::operator () (const DoFHandler &dof_handler, + GridOutput::operator () (const hp::DoFHandler &dof_handler, const Vector &/*solution*/) const { std::ostringstream filename; @@ -323,9 +323,9 @@ namespace LaplaceSolver { public: Solver (Triangulation &triangulation, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, + const hp::QCollection &face_quadrature, const Function &boundary_values); virtual ~Solver (); @@ -343,10 +343,10 @@ namespace LaplaceSolver n_dofs () const; protected: - const SmartPointer > fe; - const SmartPointer > quadrature; - const SmartPointer > face_quadrature; - DoFHandler dof_handler; + const SmartPointer > fe; + const SmartPointer > quadrature; + const SmartPointer > face_quadrature; + hp::DoFHandler dof_handler; Vector solution; const SmartPointer > boundary_values; @@ -355,7 +355,7 @@ namespace LaplaceSolver private: struct LinearSystem { - LinearSystem (const DoFHandler &dof_handler); + LinearSystem (const hp::DoFHandler &dof_handler); void solve (Vector &solution) const; @@ -370,8 +370,8 @@ namespace LaplaceSolver void assemble_matrix (LinearSystem &linear_system, - const typename DoFHandler::active_cell_iterator &begin_cell, - const typename DoFHandler::active_cell_iterator &end_cell, + const typename hp::DoFHandler::active_cell_iterator &begin_cell, + const typename hp::DoFHandler::active_cell_iterator &end_cell, Threads::ThreadMutex &mutex) const; }; @@ -379,9 +379,9 @@ namespace LaplaceSolver template Solver::Solver (Triangulation &triangulation, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, + const hp::QCollection &face_quadrature, const Function &boundary_values) : Base (triangulation), @@ -435,7 +435,7 @@ namespace LaplaceSolver Solver::assemble_linear_system (LinearSystem &linear_system) { typedef - typename DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator active_cell_iterator; const unsigned int n_threads = multithread_info.n_default_threads; @@ -476,21 +476,21 @@ namespace LaplaceSolver template void Solver::assemble_matrix (LinearSystem &linear_system, - const typename DoFHandler::active_cell_iterator &begin_cell, - const typename DoFHandler::active_cell_iterator &end_cell, + const typename hp::DoFHandler::active_cell_iterator &begin_cell, + const typename hp::DoFHandler::active_cell_iterator &end_cell, Threads::ThreadMutex &mutex) const { - FEValues fe_values (*fe, *quadrature, + hp::FEValues fe_values (*fe, *quadrature, update_gradients | update_JxW_values); - const unsigned int dofs_per_cell = fe->dofs_per_cell; - const unsigned int n_q_points = quadrature->n_quadrature_points; + const unsigned int dofs_per_cell = (*fe)[0].dofs_per_cell; + const unsigned int n_q_points = (*quadrature)[0].n_quadrature_points; FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); std::vector local_dof_indices (dofs_per_cell); - for (typename DoFHandler::active_cell_iterator cell=begin_cell; + for (typename hp::DoFHandler::active_cell_iterator cell=begin_cell; cell!=end_cell; ++cell) { cell_matrix = 0; @@ -500,9 +500,9 @@ namespace LaplaceSolver for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); @@ -518,11 +518,11 @@ namespace LaplaceSolver template Solver::LinearSystem:: - LinearSystem (const DoFHandler &dof_handler) + LinearSystem (const hp::DoFHandler &dof_handler) { hanging_node_constraints.clear (); - void (*mhnc_p) (const DoFHandler &, + void (*mhnc_p) (const hp::DoFHandler &, ConstraintMatrix &) = &DoFTools::make_hanging_node_constraints; @@ -569,9 +569,9 @@ namespace LaplaceSolver { public: PrimalSolver (Triangulation &triangulation, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, + const hp::QCollection &face_quadrature, const Function &rhs_function, const Function &boundary_values); @@ -598,9 +598,9 @@ namespace LaplaceSolver template PrimalSolver:: PrimalSolver (Triangulation &triangulation, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, + const hp::QCollection &face_quadrature, const Function &rhs_function, const Function &boundary_values) : @@ -642,7 +642,7 @@ namespace LaplaceSolver void PrimalSolver::output_solution () const { - DataOut data_out; + DataOut > data_out; data_out.attach_dof_handler (this->dof_handler); data_out.add_data_vector (this->solution, "solution"); data_out.build_patches (); @@ -662,18 +662,18 @@ namespace LaplaceSolver PrimalSolver:: assemble_rhs (Vector &rhs) const { - FEValues fe_values (*this->fe, *this->quadrature, + hp::FEValues fe_values (*this->fe, *this->quadrature, update_values | update_q_points | update_JxW_values); - const unsigned int dofs_per_cell = this->fe->dofs_per_cell; - const unsigned int n_q_points = this->quadrature->n_quadrature_points; + const unsigned int dofs_per_cell = (*this->fe)[0].dofs_per_cell; + const unsigned int n_q_points = (*this->quadrature)[0].n_quadrature_points; Vector cell_rhs (dofs_per_cell); std::vector rhs_values (n_q_points); std::vector local_dof_indices (dofs_per_cell); - typename DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator cell = this->dof_handler.begin_active(), endc = this->dof_handler.end(); for (; cell!=endc; ++cell) @@ -682,14 +682,14 @@ namespace LaplaceSolver fe_values.reinit (cell); - rhs_function->value_list (fe_values.get_quadrature_points(), + rhs_function->value_list (fe_values.get_present_fe_values().get_quadrature_points(), rhs_values); for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); for (unsigned int i=0; i &coarse_grid, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, + const hp::QCollection &face_quadrature, const Function &rhs_function, const Function &boundary_values); @@ -718,9 +718,9 @@ namespace LaplaceSolver template RefinementGlobal:: RefinementGlobal (Triangulation &coarse_grid, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, + const hp::QCollection &face_quadrature, const Function &rhs_function, const Function &boundary_values) : @@ -746,9 +746,9 @@ namespace LaplaceSolver { public: RefinementKelly (Triangulation &coarse_grid, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, + const hp::QCollection &face_quadrature, const Function &rhs_function, const Function &boundary_values); @@ -760,9 +760,9 @@ namespace LaplaceSolver template RefinementKelly:: RefinementKelly (Triangulation &coarse_grid, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, + const hp::QCollection &face_quadrature, const Function &rhs_function, const Function &boundary_values) : @@ -798,9 +798,9 @@ namespace LaplaceSolver { public: RefinementWeightedKelly (Triangulation &coarse_grid, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, + const hp::QCollection &face_quadrature, const Function &rhs_function, const Function &boundary_values, const Function &weighting_function); @@ -816,9 +816,9 @@ namespace LaplaceSolver template RefinementWeightedKelly:: RefinementWeightedKelly (Triangulation &coarse_grid, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, + const hp::QCollection &face_quadrature, const Function &rhs_function, const Function &boundary_values, const Function &weighting_function) @@ -838,12 +838,12 @@ namespace LaplaceSolver { Vector estimated_error (this->triangulation->n_active_cells()); KellyErrorEstimator::estimate (this->dof_handler, - *this->face_quadrature, + (*this->face_quadrature)[0], typename FunctionMap::type(), this->solution, estimated_error); - typename DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator cell = this->dof_handler.begin_active(), endc = this->dof_handler.end(); for (unsigned int cell_index=0; cell!=endc; ++cell, ++cell_index) @@ -1112,7 +1112,7 @@ namespace DualFunctional public: virtual void - assemble_rhs (const DoFHandler &dof_handler, + assemble_rhs (const hp::DoFHandler &dof_handler, Vector &rhs) const = 0; }; @@ -1126,7 +1126,7 @@ namespace DualFunctional virtual void - assemble_rhs (const DoFHandler &dof_handler, + assemble_rhs (const hp::DoFHandler &dof_handler, Vector &rhs) const; DeclException1 (ExcEvaluationPointNotFound, @@ -1150,12 +1150,12 @@ namespace DualFunctional template void PointValueEvaluation:: - assemble_rhs (const DoFHandler &dof_handler, + assemble_rhs (const hp::DoFHandler &dof_handler, Vector &rhs) const { rhs.reinit (dof_handler.n_dofs()); - typename DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) @@ -1182,7 +1182,7 @@ namespace DualFunctional virtual void - assemble_rhs (const DoFHandler &dof_handler, + assemble_rhs (const hp::DoFHandler &dof_handler, Vector &rhs) const; DeclException1 (ExcEvaluationPointNotFound, @@ -1206,17 +1206,17 @@ namespace DualFunctional template void PointXDerivativeEvaluation:: - assemble_rhs (const DoFHandler &dof_handler, + assemble_rhs (const hp::DoFHandler &dof_handler, Vector &rhs) const { rhs.reinit (dof_handler.n_dofs()); QGauss quadrature(4); - FEValues fe_values (dof_handler.get_fe(), quadrature, + hp::FEValues fe_values (dof_handler.get_fe(), quadrature, update_gradients | update_q_points | update_JxW_values); - const unsigned int n_q_points = fe_values.n_quadrature_points; + const unsigned int n_q_points = fe_values.get_present_fe_values().n_quadrature_points; const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; Vector cell_rhs (dofs_per_cell); @@ -1224,7 +1224,7 @@ namespace DualFunctional double total_volume = 0; - typename DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) @@ -1237,9 +1237,9 @@ namespace DualFunctional for (unsigned int q=0; qget_dof_indices (local_dof_indices); @@ -1266,9 +1266,9 @@ namespace LaplaceSolver { public: DualSolver (Triangulation &triangulation, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, + const hp::QCollection &face_quadrature, const DualFunctional::DualFunctionalBase &dual_functional); virtual @@ -1298,9 +1298,9 @@ namespace LaplaceSolver template DualSolver:: DualSolver (Triangulation &triangulation, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, + const hp::FECollection &fe, + const hp::QCollection &quadrature, + const hp::QCollection &face_quadrature, const DualFunctional::DualFunctionalBase &dual_functional) : Base (triangulation), @@ -1354,10 +1354,10 @@ namespace LaplaceSolver { public: WeightedResidual (Triangulation &coarse_grid, - const FiniteElement &primal_fe, - const FiniteElement &dual_fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, + const hp::FECollection &primal_fe, + const hp::FECollection &dual_fe, + const hp::QCollection &quadrature, + const hp::QCollection &face_quadrature, const Function &rhs_function, const Function &boundary_values, const DualFunctional::DualFunctionalBase &dual_functional); @@ -1385,39 +1385,39 @@ namespace LaplaceSolver void solve_dual_problem (); typedef - typename DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator active_cell_iterator; typedef - typename std::map::face_iterator,double> + typename std::map::face_iterator,double> FaceIntegrals; struct CellData { - FEValues fe_values; + hp::FEValues fe_values; const SmartPointer > right_hand_side; std::vector cell_residual; std::vector rhs_values; std::vector dual_weights; typename std::vector > cell_grad_grads; - CellData (const FiniteElement &fe, - const Quadrature &quadrature, + CellData (const hp::FECollection &fe, + const hp::QCollection &quadrature, const Function &right_hand_side); }; struct FaceData { - FEFaceValues fe_face_values_cell; - FEFaceValues fe_face_values_neighbor; - FESubfaceValues fe_subface_values_cell; + hp::FEFaceValues fe_face_values_cell; + hp::FEFaceValues fe_face_values_neighbor; + hp::FESubfaceValues fe_subface_values_cell; std::vector jump_residual; std::vector dual_weights; typename std::vector > cell_grads; typename std::vector > neighbor_grads; - FaceData (const FiniteElement &fe, - const Quadrature &face_quadrature); + FaceData (const hp::FECollection &fe, + const hp::QCollection &face_quadrature); }; @@ -1459,8 +1459,8 @@ namespace LaplaceSolver template WeightedResidual::CellData:: - CellData (const FiniteElement &fe, - const Quadrature &quadrature, + CellData (const hp::FECollection &fe, + const hp::QCollection &quadrature, const Function &right_hand_side) : fe_values (fe, quadrature, @@ -1471,7 +1471,7 @@ namespace LaplaceSolver right_hand_side (&right_hand_side) { const unsigned int n_q_points - = quadrature.n_quadrature_points; + = quadrature[0].n_quadrature_points; cell_residual.resize(n_q_points); rhs_values.resize(n_q_points); @@ -1483,8 +1483,8 @@ namespace LaplaceSolver template WeightedResidual::FaceData:: - FaceData (const FiniteElement &fe, - const Quadrature &face_quadrature) + FaceData (const hp::FECollection &fe, + const hp::QCollection &face_quadrature) : fe_face_values_cell (fe, face_quadrature, update_values | @@ -1500,7 +1500,7 @@ namespace LaplaceSolver update_gradients) { const unsigned int n_face_q_points - = face_quadrature.n_quadrature_points; + = face_quadrature[0].n_quadrature_points; jump_residual.resize(n_face_q_points); dual_weights.resize(n_face_q_points); @@ -1514,10 +1514,10 @@ namespace LaplaceSolver template WeightedResidual:: WeightedResidual (Triangulation &coarse_grid, - const FiniteElement &primal_fe, - const FiniteElement &dual_fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, + const hp::FECollection &primal_fe, + const hp::FECollection &dual_fe, + const hp::QCollection &quadrature, + const hp::QCollection &face_quadrature, const Function &rhs_function, const Function &bv, const DualFunctional::DualFunctionalBase &dual_functional) @@ -1612,7 +1612,7 @@ namespace LaplaceSolver primal_hanging_node_constraints, dual_solution); - DataOut data_out; + DataOut > data_out; data_out.attach_dof_handler (primal_solver.dof_handler); data_out.add_data_vector (primal_solver.solution, @@ -1799,19 +1799,19 @@ namespace LaplaceSolver { cell_data.fe_values.reinit (cell); cell_data.right_hand_side - ->value_list (cell_data.fe_values.get_quadrature_points(), + ->value_list (cell_data.fe_values.get_present_fe_values().get_quadrature_points(), cell_data.rhs_values); - cell_data.fe_values.get_function_2nd_derivatives (primal_solution, + cell_data.fe_values.get_present_fe_values().get_function_2nd_derivatives (primal_solution, cell_data.cell_grad_grads); - cell_data.fe_values.get_function_values (dual_weights, + cell_data.fe_values.get_present_fe_values().get_function_values (dual_weights, cell_data.dual_weights); double sum = 0; - for (unsigned int p=0; pneighbor(face_no).state() == IteratorState::valid, @@ -1839,22 +1839,22 @@ namespace LaplaceSolver neighbor_neighbor = cell->neighbor_of_neighbor (face_no); const active_cell_iterator neighbor = cell->neighbor(face_no); face_data.fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor); - face_data.fe_face_values_neighbor.get_function_grads (primal_solution, + face_data.fe_face_values_neighbor.get_present_fe_values().get_function_grads (primal_solution, face_data.neighbor_grads); for (unsigned int p=0; pface(face_no)) != face_integrals.end(), ExcInternalError()); @@ -1876,11 +1876,11 @@ namespace LaplaceSolver FaceIntegrals &face_integrals) const { const unsigned int - n_q_points = face_data.fe_face_values_cell.n_quadrature_points; + n_q_points = face_data.fe_face_values_cell.get_present_fe_values().n_quadrature_points; - const typename DoFHandler::face_iterator + const typename hp::DoFHandler::face_iterator face = cell->face(face_no); - const typename DoFHandler::cell_iterator + const typename hp::DoFHandler::cell_iterator neighbor = cell->neighbor(face_no); Assert (neighbor.state() == IteratorState::valid, ExcInternalError()); @@ -1900,26 +1900,26 @@ namespace LaplaceSolver ExcInternalError()); face_data.fe_subface_values_cell.reinit (cell, face_no, subface_no); - face_data.fe_subface_values_cell.get_function_grads (primal_solution, + face_data.fe_subface_values_cell.get_present_fe_values().get_function_grads (primal_solution, face_data.cell_grads); face_data.fe_face_values_neighbor.reinit (neighbor_child, neighbor_neighbor); - face_data.fe_face_values_neighbor.get_function_grads (primal_solution, + face_data.fe_face_values_neighbor.get_present_fe_values().get_function_grads (primal_solution, face_data.neighbor_grads); for (unsigned int p=0; pface(neighbor_neighbor)] = face_integral; }; @@ -2000,10 +2000,10 @@ void Framework::run (const ProblemDescription &descriptor) triangulation (Triangulation::smoothing_on_refinement); descriptor.data->create_coarse_grid (triangulation); - const FE_Q primal_fe(descriptor.primal_fe_degree); - const FE_Q dual_fe(descriptor.dual_fe_degree); - const QGauss quadrature(descriptor.dual_fe_degree+1); - const QGauss face_quadrature(descriptor.dual_fe_degree+1); + const hp::FECollection primal_fe(FE_Q(descriptor.primal_fe_degree)); + const hp::FECollection dual_fe(FE_Q(descriptor.dual_fe_degree)); + const hp::QCollection quadrature(QGauss(descriptor.dual_fe_degree+1)); + const hp::QCollection face_quadrature(QGauss(descriptor.dual_fe_degree+1)); LaplaceSolver::Base * solver = 0; switch (descriptor.refinement_criterion) diff --git a/tests/hp/step-15.cc b/tests/hp/step-15.cc index ec7b7457f5..49c949c68b 100644 --- a/tests/hp/step-15.cc +++ b/tests/hp/step-15.cc @@ -38,7 +38,8 @@ std::ofstream logfile("step-15/output"); #include #include #include -#include +#include +#include #include #include #include @@ -101,7 +102,7 @@ class MinimizationProblem void output_results () const; void refine_grid (); - static double energy (const DoFHandler &dof_handler, + static double energy (const hp::DoFHandler &dof_handler, const Vector &function); @@ -109,8 +110,8 @@ class MinimizationProblem Triangulation triangulation; - FE_Q fe; - DoFHandler dof_handler; + hp::FECollection fe; + hp::DoFHandler dof_handler; ConstraintMatrix hanging_node_constraints; @@ -127,7 +128,7 @@ template MinimizationProblem::MinimizationProblem (const unsigned int run_number) : run_number (run_number), - fe (1), + fe (FE_Q(1)), dof_handler (triangulation) {} @@ -140,7 +141,7 @@ void MinimizationProblem<1>::initialize_solution () InitializationValues(), present_solution); - DoFHandler<1>::cell_iterator cell; + hp::DoFHandler<1>::cell_iterator cell; cell = dof_handler.begin(0); while (cell->at_boundary(0) == false) cell = cell->neighbor(0); @@ -183,13 +184,13 @@ void MinimizationProblem::assemble_step () matrix.reinit (sparsity_pattern); residual.reinit (dof_handler.n_dofs()); - QGauss quadrature_formula(4); - FEValues fe_values (fe, quadrature_formula, + hp::QCollection quadrature_formula(QGauss(4)); + hp::FEValues fe_values (fe, quadrature_formula, update_values | update_gradients | update_q_points | update_JxW_values); - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + const unsigned int dofs_per_cell = fe[0].dofs_per_cell; + const unsigned int n_q_points = quadrature_formula[0].n_quadrature_points; FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); Vector cell_rhs (dofs_per_cell); @@ -199,7 +200,7 @@ void MinimizationProblem::assemble_step () std::vector local_solution_values (n_q_points); std::vector > local_solution_grads (n_q_points); - typename DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) @@ -209,44 +210,44 @@ void MinimizationProblem::assemble_step () fe_values.reinit (cell); - fe_values.get_function_values (present_solution, + fe_values.get_present_fe_values().get_function_values (present_solution, local_solution_values); - fe_values.get_function_grads (present_solution, + fe_values.get_present_fe_values().get_function_grads (present_solution, local_solution_grads); for (unsigned int q_point=0; q_point u_prime = local_solution_grads[q_point]; for (unsigned int i=0; idiameter() * cell->diameter() + - fe_values.shape_value(i,q_point) * - fe_values.shape_value(j,q_point)) * - fe_values.JxW(q_point); + fe_values.get_present_fe_values().shape_value(i,q_point) * + fe_values.get_present_fe_values().shape_value(j,q_point)) * + fe_values.get_present_fe_values().JxW(q_point); for (unsigned int i=0; iget_dof_indices (local_dof_indices); @@ -366,7 +367,7 @@ template void MinimizationProblem::output_results () const { - DataOut data_out; + DataOut > data_out; data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (present_solution, "solution"); data_out.build_patches (); @@ -389,33 +390,34 @@ void MinimizationProblem<1>::refine_grid () Vector error_indicators (triangulation.n_active_cells()); - QTrapez quadrature; - FEValues fe_values (fe, quadrature, + QTrapez q; + hp::QCollection quadrature(q); + hp::FEValues fe_values (fe, quadrature, update_values | update_gradients | update_second_derivatives | update_q_points | update_JxW_values); - FEValues neighbor_fe_values (fe, quadrature, + hp::FEValues neighbor_fe_values (fe, quadrature, update_gradients); - std::vector local_values (quadrature.n_quadrature_points); - std::vector > local_gradients (quadrature.n_quadrature_points); - std::vector > local_2nd_derivs (quadrature.n_quadrature_points); + std::vector local_values (quadrature[0].n_quadrature_points); + std::vector > local_gradients (quadrature[0].n_quadrature_points); + std::vector > local_2nd_derivs (quadrature[0].n_quadrature_points); - DoFHandler::active_cell_iterator + hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active (), endc = dof_handler.end (); for (unsigned int cell_index = 0; cell!=endc; ++cell, ++cell_index) { fe_values.reinit (cell); - fe_values.get_function_values (present_solution, local_values); - fe_values.get_function_grads (present_solution, local_gradients); - fe_values.get_function_2nd_derivatives (present_solution, local_2nd_derivs); + fe_values.get_present_fe_values().get_function_values (present_solution, local_values); + fe_values.get_present_fe_values().get_function_grads (present_solution, local_gradients); + fe_values.get_present_fe_values().get_function_2nd_derivatives (present_solution, local_2nd_derivs); double cell_residual_norm = 0; - for (unsigned int q=0; q::refine_grid () 2*u_prime*(1-3*u*u*u_prime))); cell_residual_norm += (local_residual_value * local_residual_value * - fe_values.JxW(q)); + fe_values.get_present_fe_values().JxW(q)); } error_indicators(cell_index) = cell_residual_norm * cell->diameter() * cell->diameter(); - const double x_left = fe_values.quadrature_point(0)[0]; - const double x_right = fe_values.quadrature_point(1)[0]; + const double x_left = fe_values.get_present_fe_values().quadrature_point(0)[0]; + const double x_right = fe_values.get_present_fe_values().quadrature_point(1)[0]; Assert (x_left == cell->vertex(0)[0], ExcInternalError()); Assert (x_right == cell->vertex(1)[0], ExcInternalError()); @@ -447,12 +449,12 @@ void MinimizationProblem<1>::refine_grid () if (cell->at_boundary(0) == false) { - DoFHandler::cell_iterator left_neighbor = cell->neighbor(0); + hp::DoFHandler::cell_iterator left_neighbor = cell->neighbor(0); while (left_neighbor->has_children()) left_neighbor = left_neighbor->child(1); neighbor_fe_values.reinit (left_neighbor); - neighbor_fe_values.get_function_grads (present_solution, local_gradients); + neighbor_fe_values.get_present_fe_values().get_function_grads (present_solution, local_gradients); const double neighbor_u_prime_left = local_gradients[1][0]; @@ -465,12 +467,12 @@ void MinimizationProblem<1>::refine_grid () if (cell->at_boundary(1) == false) { - DoFHandler::cell_iterator right_neighbor = cell->neighbor(1); + hp::DoFHandler::cell_iterator right_neighbor = cell->neighbor(1); while (right_neighbor->has_children()) right_neighbor = right_neighbor->child(0); neighbor_fe_values.reinit (right_neighbor); - neighbor_fe_values.get_function_grads (present_solution, local_gradients); + neighbor_fe_values.get_present_fe_values().get_function_grads (present_solution, local_gradients); const double neighbor_u_prime_right = local_gradients[0][0]; @@ -508,39 +510,39 @@ void MinimizationProblem<1>::refine_grid () template double -MinimizationProblem::energy (const DoFHandler &dof_handler, +MinimizationProblem::energy (const hp::DoFHandler &dof_handler, const Vector &function) { - QGauss quadrature_formula(4); - FEValues fe_values (dof_handler.get_fe(), quadrature_formula, + hp::QCollection quadrature_formula(QGauss(4)); + hp::FEValues fe_values (dof_handler.get_fe(), quadrature_formula, update_values | update_gradients | update_q_points | update_JxW_values); - const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + const unsigned int n_q_points = quadrature_formula[0].n_quadrature_points; std::vector local_solution_values (n_q_points); std::vector > local_solution_grads (n_q_points); double energy = 0.; - 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); - fe_values.get_function_values (function, + fe_values.get_present_fe_values().get_function_values (function, local_solution_values); - fe_values.get_function_grads (function, + fe_values.get_present_fe_values().get_function_grads (function, local_solution_grads); for (unsigned int q_point=0; q_point