From: wolf Date: Wed, 3 May 2006 15:04:28 +0000 (+0000) Subject: Add more tests that all fail for various reasons at the moment. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f4e90ad7562e182ebbc8c674c807892bddf9c702;p=dealii-svn.git Add more tests that all fail for various reasons at the moment. git-svn-id: https://svn.dealii.org/trunk@13018 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/hp/step-10.cc b/tests/hp/step-10.cc new file mode 100644 index 0000000000..733b981986 --- /dev/null +++ b/tests/hp/step-10.cc @@ -0,0 +1,236 @@ +//---------------------------- step-10.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005, 2006 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- step-10.cc --------------------------- + + +// a hp-ified version of step-10 + + +#include +#include +std::ofstream logfile("step-10/output"); + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +const long double pi = 3.141592653589793238462643; + + + +template +void gnuplot_output() +{ + deallog << "Output of grids into gnuplot files:" << std::endl + << "===================================" << std::endl; + + Triangulation triangulation; + GridGenerator::hyper_ball (triangulation); + static const HyperBallBoundary boundary; + triangulation.set_boundary (0, boundary); + + for (unsigned int refinement=0; refinement<2; + ++refinement, triangulation.refine_global(1)) + { + deallog << "Refinement level: " << refinement << std::endl; + + std::string filename_base = "ball"; + filename_base += '0'+refinement; + + for (unsigned int degree=1; degree<4; ++degree) + { + deallog << "Degree = " << degree << std::endl; + + const MappingQ mapping (degree); + + + GridOut grid_out; + GridOutFlags::Gnuplot gnuplot_flags(false, 30); + grid_out.set_flags(gnuplot_flags); + + std::string filename = filename_base+"_mapping_q"; + filename += ('0'+degree); + filename += ".dat"; + std::ofstream gnuplot_file (filename.c_str()); + + grid_out.write_gnuplot (triangulation, gnuplot_file, &mapping); + } + deallog << std::endl; + } +} + +template +void compute_pi_by_area () +{ + deallog << "Computation of Pi by the area:" << std::endl + << "==============================" << std::endl; + + const hp::QCollection quadrature(QGauss(4)); + + for (unsigned int degree=1; degree<5; ++degree) + { + deallog << "Degree = " << degree << std::endl; + + Triangulation triangulation; + GridGenerator::hyper_ball (triangulation); + + static const HyperBallBoundary boundary; + triangulation.set_boundary (0, boundary); + + MappingQ m(degree); + const hp::MappingCollection mapping (m); + + const hp::FECollection dummy_fe (FE_Q(1)); + + hp::DoFHandler dof_handler (triangulation); + + hp::FEValues x_fe_values (mapping, dummy_fe, quadrature, + update_JxW_values); + + ConvergenceTable table; + + for (unsigned int refinement=0; refinement<6; + ++refinement, triangulation.refine_global (1)) + { + table.add_value("cells", triangulation.n_active_cells()); + + dof_handler.distribute_dofs (dummy_fe); + + long double area = 0; + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + x_fe_values.reinit (cell); + const FEValues &fe_values = x_fe_values.get_present_fe_values(); + for (unsigned int i=0; i (area)); + table.add_value("error", static_cast (std::fabs(area-pi))); + }; + + table.omit_column_from_convergence_rate_evaluation("cells"); + table.omit_column_from_convergence_rate_evaluation("eval.pi"); + table.evaluate_all_convergence_rates(ConvergenceTable::reduction_rate_log2); + + table.set_precision("eval.pi", 16); + table.set_scientific("error", true); + + table.write_text(deallog.get_file_stream()); + + deallog << std::endl; + }; +} + + +template +void compute_pi_by_perimeter () +{ + deallog << "Computation of Pi by the perimeter:" << std::endl + << "===================================" << std::endl; + + const hp::QCollection quadrature(4); + + for (unsigned int degree=1; degree<5; ++degree) + { + deallog << "Degree = " << degree << std::endl; + Triangulation triangulation; + GridGenerator::hyper_ball (triangulation); + + static const HyperBallBoundary boundary; + triangulation.set_boundary (0, boundary); + + const MappingQ m (degree); + const hp::MappingCollection mapping(m); + const hp::FECollection fe (FE_Q(1)); + + hp::DoFHandler dof_handler (triangulation); + + hp::FEFaceValues x_fe_face_values (mapping, fe, quadrature, + update_JxW_values); + ConvergenceTable table; + + for (unsigned int refinement=0; refinement<6; + ++refinement, triangulation.refine_global (1)) + { + table.add_value("cells", triangulation.n_active_cells()); + + dof_handler.distribute_dofs (fe); + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + long double perimeter = 0; + for (; cell!=endc; ++cell) + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + if (cell->face(face_no)->at_boundary()) + { + x_fe_face_values.reinit (cell, face_no); + const FEFaceValues &fe_face_values + = x_fe_face_values.get_present_fe_values (); + + for (unsigned int i=0; i (perimeter/2.)); + table.add_value("error", static_cast (std::fabs(perimeter/2.-pi))); + }; + + table.omit_column_from_convergence_rate_evaluation("cells"); + table.omit_column_from_convergence_rate_evaluation("eval.pi"); + table.evaluate_all_convergence_rates(ConvergenceTable::reduction_rate_log2); + + table.set_precision("eval.pi", 16); + table.set_scientific("error", true); + + table.write_text(deallog.get_file_stream()); + + deallog << std::endl; + }; +} + + +int main () +{ + logfile.precision(2); + + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + gnuplot_output<2>(); + + compute_pi_by_area<2> (); + compute_pi_by_perimeter<2> (); + + return 0; +} diff --git a/tests/hp/step-10/.cvsignore b/tests/hp/step-10/.cvsignore new file mode 100644 index 0000000000..d196dfb299 --- /dev/null +++ b/tests/hp/step-10/.cvsignore @@ -0,0 +1 @@ +obj.* exe OK output diff --git a/tests/hp/step-6.cc b/tests/hp/step-6.cc new file mode 100644 index 0000000000..2731346c4d --- /dev/null +++ b/tests/hp/step-6.cc @@ -0,0 +1,415 @@ +//---------------------------- step-6.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005, 2006 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- step-6.cc --------------------------- + + +// a hp-ified version of step-6 + + +#include +#include +std::ofstream logfile("step-6/output"); + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + + +#include + +#include + +#include + + + + +template +class LaplaceProblem +{ + public: + LaplaceProblem (); + ~LaplaceProblem (); + + void run (); + + private: + void setup_system (); + void assemble_system (); + void solve (); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + + hp::DoFHandler dof_handler; + hp::FECollection fe; + + ConstraintMatrix hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + + + +template +class Coefficient : public Function +{ + public: + Coefficient () : Function() {}; + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void value_list (const std::vector > &points, + std::vector &values, + const unsigned int component = 0) const; +}; + + + +template +double Coefficient::value (const Point &p, + const unsigned int) const +{ + if (p.square() < 0.5*0.5) + return 20; + else + return 1; +} + + + +template +void Coefficient::value_list (const std::vector > &points, + std::vector &values, + const unsigned int component) const +{ + const unsigned int n_points = points.size(); + + Assert (values.size() == n_points, + ExcDimensionMismatch (values.size(), n_points)); + + Assert (component == 0, + ExcIndexRange (component, 0, 1)); + + for (unsigned int i=0; i +LaplaceProblem::LaplaceProblem () : + dof_handler (triangulation), + fe (FE_Q(2)) +{} + + + +template +LaplaceProblem::~LaplaceProblem () +{ + dof_handler.clear (); +} + + + +template +void LaplaceProblem::setup_system () +{ + dof_handler.distribute_dofs (fe); + + sparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + + + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, + hanging_node_constraints); + + hanging_node_constraints.close (); + + hanging_node_constraints.condense (sparsity_pattern); + + sparsity_pattern.compress(); + + system_matrix.reinit (sparsity_pattern); +} + + +template +void LaplaceProblem::assemble_system () +{ + const hp::QCollection quadrature_formula(3); + + hp::FEValues x_fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_q_points | update_JxW_values); + + 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); + + std::vector local_dof_indices (dofs_per_cell); + + const Coefficient coefficient; + std::vector coefficient_values (n_q_points); + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_matrix = 0; + cell_rhs = 0; + + x_fe_values.reinit (cell); + const FEValues &fe_values = x_fe_values.get_present_fe_values(); + + coefficient.value_list (fe_values.get_quadrature_points(), + coefficient_values); + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + for (unsigned int i=0; i boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); +} + + + + + +template +void LaplaceProblem::solve () +{ + SolverControl solver_control (1000, 1e-12); + SolverCG<> cg (solver_control); + + PreconditionSSOR<> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + + hanging_node_constraints.distribute (solution); +} + + + +template +void LaplaceProblem::refine_grid () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + KellyErrorEstimator::estimate (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 (); +} + + + +template +void LaplaceProblem::output_results (const unsigned int cycle) const +{ + Assert (cycle < 10, ExcNotImplemented()); + + std::string filename = "grid-"; + filename += ('0' + cycle); + filename += ".eps"; + + std::ofstream output (filename.c_str()); + + GridOut grid_out; + grid_out.write_eps (triangulation, output); +} + + + + +template +void LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<8; ++cycle) + { + deallog << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_ball (triangulation); + + static const HyperBallBoundary boundary; + triangulation.set_boundary (0, boundary); + + triangulation.refine_global (1); + } + else + refine_grid (); + + + deallog << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + setup_system (); + + deallog << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + assemble_system (); + solve (); + output_results (cycle); + } + + DataOutBase::EpsFlags eps_flags; + eps_flags.z_scaling = 4; + + DataOut > data_out; + data_out.set_flags (eps_flags); + + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "solution"); + data_out.build_patches (); + + std::ofstream output ("final-solution.eps"); + data_out.write_eps (output); +} + + + +int main () +{ + logfile.precision(2); + + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + deallog.depth_console (0); + + LaplaceProblem<2> laplace_problem_2d; + laplace_problem_2d.run (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/hp/step-6/.cvsignore b/tests/hp/step-6/.cvsignore new file mode 100644 index 0000000000..d196dfb299 --- /dev/null +++ b/tests/hp/step-6/.cvsignore @@ -0,0 +1 @@ +obj.* exe OK output diff --git a/tests/hp/step-7.cc b/tests/hp/step-7.cc new file mode 100644 index 0000000000..e053f7c411 --- /dev/null +++ b/tests/hp/step-7.cc @@ -0,0 +1,742 @@ +//---------------------------- step-7.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005, 2006 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- step-7.cc --------------------------- + + +// a hp-ified version of step-7 + + +#include +#include +std::ofstream logfile("step-7/output"); + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + + + +template +class SolutionBase +{ + protected: + static const unsigned int n_source_centers = 3; + static const Point source_centers[n_source_centers]; + static const double width; +}; + + +template <> +const Point<1> +SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers] += { Point<1>(-1.0 / 3.0), + Point<1>(0.0), + Point<1>(+1.0 / 3.0) }; + +template <> +const Point<2> +SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers] += { Point<2>(-0.5, +0.5), + Point<2>(-0.5, -0.5), + Point<2>(+0.5, -0.5) }; + +template +const double SolutionBase::width = 1./3.; + + + +template +class Solution : public Function, + protected SolutionBase +{ + public: + Solution () : Function() {}; + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; +}; + + +template +double Solution::value (const Point &p, + const unsigned int) const +{ + double return_value = 0; + for (unsigned int i=0; in_source_centers; ++i) + { + const Point x_minus_xi = p - this->source_centers[i]; + return_value += std::exp(-x_minus_xi.square() / + (this->width * this->width)); + } + + return return_value; +} + + +template +Tensor<1,dim> Solution::gradient (const Point &p, + const unsigned int) const +{ + Tensor<1,dim> return_value; + + for (unsigned int i=0; in_source_centers; ++i) + { + const Point x_minus_xi = p - this->source_centers[i]; + + return_value += (-2 / (this->width * this->width) * + std::exp(-x_minus_xi.square() / + (this->width * this->width)) * + x_minus_xi); + } + + return return_value; +} + + + +template +class RightHandSide : public Function, + protected SolutionBase +{ + public: + RightHandSide () : Function() {}; + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + +template +double RightHandSide::value (const Point &p, + const unsigned int) const +{ + double return_value = 0; + for (unsigned int i=0; in_source_centers; ++i) + { + const Point x_minus_xi = p - this->source_centers[i]; + + return_value += ((2*dim - 4*x_minus_xi.square()/ + (this->width * this->width)) / + (this->width * this->width) * + std::exp(-x_minus_xi.square() / + (this->width * this->width))); + return_value += std::exp(-x_minus_xi.square() / + (this->width * this->width)); + } + + return return_value; +} + + + +template +class HelmholtzProblem +{ + public: + enum RefinementMode { + global_refinement, adaptive_refinement + }; + + HelmholtzProblem (const hp::FECollection &fe, + const RefinementMode refinement_mode); + + ~HelmholtzProblem (); + + void run (); + + private: + void setup_system (); + void assemble_system (); + void solve (); + void refine_grid (); + void process_solution (const unsigned int cycle); + + Triangulation triangulation; + hp::DoFHandler dof_handler; + + SmartPointer > fe; + + ConstraintMatrix hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + + const RefinementMode refinement_mode; + + ConvergenceTable convergence_table; +}; + + + + +template +HelmholtzProblem::HelmholtzProblem (const hp::FECollection &fe, + const RefinementMode refinement_mode) : + dof_handler (triangulation), + fe (&fe), + refinement_mode (refinement_mode) +{} + + + +template +HelmholtzProblem::~HelmholtzProblem () +{ + dof_handler.clear (); +} + + + +template +void HelmholtzProblem::setup_system () +{ + dof_handler.distribute_dofs (*fe); + DoFRenumbering::Cuthill_McKee (dof_handler); + + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, + hanging_node_constraints); + hanging_node_constraints.close (); + + sparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); + hanging_node_constraints.condense (sparsity_pattern); + sparsity_pattern.compress(); + + system_matrix.reinit (sparsity_pattern); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); +} + + + +template +void HelmholtzProblem::assemble_system () +{ + hp::QCollection quadrature_formula(3); + hp::QCollection face_quadrature_formula(3); + + const unsigned int n_q_points = quadrature_formula[0].n_quadrature_points; + const unsigned int n_face_q_points = face_quadrature_formula[0].n_quadrature_points; + + const unsigned int dofs_per_cell = (*fe)[0].dofs_per_cell; + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + hp::FEValues x_fe_values (*fe, quadrature_formula, + update_values | update_gradients | + update_q_points | update_JxW_values); + + hp::FEFaceValues x_fe_face_values (*fe, face_quadrature_formula, + update_values | update_q_points | + update_normal_vectors | update_JxW_values); + + const RightHandSide right_hand_side; + std::vector rhs_values (n_q_points); + + const Solution exact_solution; + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_matrix = 0; + cell_rhs = 0; + + x_fe_values.reinit (cell); + const FEValues &fe_values = x_fe_values.get_present_fe_values(); + + right_hand_side.value_list (fe_values.get_quadrature_points(), + rhs_values); + + for (unsigned int q_point=0; q_point::faces_per_cell; ++face) + if (cell->face(face)->at_boundary() + && + (cell->face(face)->boundary_indicator() == 1)) + { + x_fe_face_values.reinit (cell, face); + const FEFaceValues &fe_face_values + = x_fe_face_values.get_present_fe_values (); + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + for (unsigned int i=0; i boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + Solution(), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); +} + + + +template +void HelmholtzProblem::solve () +{ + SolverControl solver_control (1000, 1e-12); + SolverCG<> cg (solver_control); + + PreconditionSSOR<> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + + hanging_node_constraints.distribute (solution); +} + + + +template +void HelmholtzProblem::refine_grid () +{ + switch (refinement_mode) + { + case global_refinement: + { + triangulation.refine_global (1); + break; + } + + case adaptive_refinement: + { + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + typename FunctionMap::type neumann_boundary; + KellyErrorEstimator::estimate (dof_handler, + hp::QCollection(3), + neumann_boundary, + 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 (); + + break; + } + + default: + { + Assert (false, ExcNotImplemented()); + } + } +} + + + +template +void HelmholtzProblem::process_solution (const unsigned int cycle) +{ + Vector difference_per_cell (triangulation.n_active_cells()); + VectorTools::integrate_difference (dof_handler, + solution, + Solution(), + difference_per_cell, + hp::QCollection(3), + VectorTools::L2_norm); + const double L2_error = difference_per_cell.l2_norm(); + + VectorTools::integrate_difference (dof_handler, + solution, + Solution(), + difference_per_cell, + hp::QCollection(3), + VectorTools::H1_seminorm); + const double H1_error = difference_per_cell.l2_norm(); + + const QTrapez<1> q_trapez; + const QIterated q_iterated (q_trapez, 5); + VectorTools::integrate_difference (dof_handler, + solution, + Solution(), + difference_per_cell, + q_iterated, + VectorTools::Linfty_norm); + const double Linfty_error = difference_per_cell.linfty_norm(); + + const unsigned int n_active_cells=triangulation.n_active_cells(); + const unsigned int n_dofs=dof_handler.n_dofs(); + + deallog << "Cycle " << cycle << ':' + << std::endl + << " Number of active cells: " + << n_active_cells + << std::endl + << " Number of degrees of freedom: " + << n_dofs + << std::endl; + + convergence_table.add_value("cycle", cycle); + convergence_table.add_value("cells", n_active_cells); + convergence_table.add_value("dofs", n_dofs); + convergence_table.add_value("L2", L2_error); + convergence_table.add_value("H1", H1_error); + convergence_table.add_value("Linfty", Linfty_error); +} + + + +template +void HelmholtzProblem::run () +{ + for (unsigned int cycle=0; cycle<7; ++cycle) + { + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (1); + + typename Triangulation::cell_iterator + cell = triangulation.begin (), + endc = triangulation.end(); + for (; cell!=endc; ++cell) + for (unsigned int face=0; + face::faces_per_cell; + ++face) + if ((cell->face(face)->center()(0) == -1) + || + (cell->face(face)->center()(1) == -1)) + cell->face(face)->set_boundary_indicator (1); + } + else + refine_grid (); + + + setup_system (); + + assemble_system (); + solve (); + + process_solution (cycle); + } + + + std::string gmv_filename; + switch (refinement_mode) + { + case global_refinement: + gmv_filename = "solution-global"; + break; + case adaptive_refinement: + gmv_filename = "solution-adaptive"; + break; + default: + Assert (false, ExcNotImplemented()); + } + + switch ((*fe)[0].degree) + { + case 1: + gmv_filename += "-q1"; + break; + case 2: + gmv_filename += "-q2"; + break; + + default: + Assert (false, ExcNotImplemented()); + } + + gmv_filename += ".gmv"; + std::ofstream output (gmv_filename.c_str()); + + DataOut > data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "solution"); + + data_out.build_patches ((*fe)[0].degree); + data_out.write_gmv (output); + + + + convergence_table.set_precision("L2", 3); + convergence_table.set_precision("H1", 3); + convergence_table.set_precision("Linfty", 3); + + convergence_table.set_scientific("L2", true); + convergence_table.set_scientific("H1", true); + convergence_table.set_scientific("Linfty", true); + + convergence_table.set_tex_caption("cells", "\\# cells"); + convergence_table.set_tex_caption("dofs", "\\# dofs"); + convergence_table.set_tex_caption("L2", "$L^2$-error"); + convergence_table.set_tex_caption("H1", "$H^1$-error"); + convergence_table.set_tex_caption("Linfty", "$L^\\infty$-error"); + + convergence_table.set_tex_format("cells", "r"); + convergence_table.set_tex_format("dofs", "r"); + + deallog << std::endl; + convergence_table.write_text(deallog); + + std::string error_filename = "error"; + switch (refinement_mode) + { + case global_refinement: + error_filename += "-global"; + break; + case adaptive_refinement: + error_filename += "-adaptive"; + break; + default: + Assert (false, ExcNotImplemented()); + } + + switch ((*fe)[0].degree) + { + case 1: + error_filename += "-q1"; + break; + case 2: + error_filename += "-q2"; + break; + default: + Assert (false, ExcNotImplemented()); + } + + error_filename += ".tex"; + std::ofstream error_table_file(error_filename.c_str()); + + convergence_table.write_tex(error_table_file); + + + + if (refinement_mode==global_refinement) + { + convergence_table.add_column_to_supercolumn("cycle", "n cells"); + convergence_table.add_column_to_supercolumn("cells", "n cells"); + + std::vector new_order; + new_order.push_back("n cells"); + new_order.push_back("H1"); + new_order.push_back("L2"); + convergence_table.set_column_order (new_order); + + convergence_table + .evaluate_convergence_rates("L2", ConvergenceTable::reduction_rate); + convergence_table + .evaluate_convergence_rates("L2", ConvergenceTable::reduction_rate_log2); + convergence_table + .evaluate_convergence_rates("H1", ConvergenceTable::reduction_rate_log2); + + deallog << std::endl; + convergence_table.write_text(deallog); + + std::string conv_filename = "convergence"; + switch (refinement_mode) + { + case global_refinement: + conv_filename += "-global"; + break; + case adaptive_refinement: + conv_filename += "-adaptive"; + break; + default: + Assert (false, ExcNotImplemented()); + } + switch ((*fe)[0].degree) + { + case 1: + conv_filename += "-q1"; + break; + case 2: + conv_filename += "-q2"; + break; + default: + Assert (false, ExcNotImplemented()); + } + conv_filename += ".tex"; + + std::ofstream table_file(conv_filename.c_str()); + convergence_table.write_tex(table_file); + } +} + + +int main () +{ + logfile.precision(2); + + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + const unsigned int dim = 2; + + try + { + deallog.depth_console (0); + + + { + deallog << "Solving with Q1 elements, adaptive refinement" << std::endl + << "=============================================" << std::endl + << std::endl; + + hp::FECollection fe(FE_Q(1)); + HelmholtzProblem + helmholtz_problem_2d (fe, HelmholtzProblem::adaptive_refinement); + + helmholtz_problem_2d.run (); + + deallog << std::endl; + } + + { + deallog << "Solving with Q1 elements, global refinement" << std::endl + << "===========================================" << std::endl + << std::endl; + + hp::FECollection fe(FE_Q(1)); + HelmholtzProblem + helmholtz_problem_2d (fe, HelmholtzProblem::global_refinement); + + helmholtz_problem_2d.run (); + + deallog << std::endl; + } + + { + deallog << "Solving with Q2 elements, global refinement" << std::endl + << "===========================================" << std::endl + << std::endl; + + hp::FECollection fe(FE_Q(2)); + HelmholtzProblem + helmholtz_problem_2d (fe, HelmholtzProblem::global_refinement); + + helmholtz_problem_2d.run (); + + deallog << std::endl; + } + + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} + + +template const double SolutionBase<2>::width; diff --git a/tests/hp/step-7/.cvsignore b/tests/hp/step-7/.cvsignore new file mode 100644 index 0000000000..d196dfb299 --- /dev/null +++ b/tests/hp/step-7/.cvsignore @@ -0,0 +1 @@ +obj.* exe OK output diff --git a/tests/hp/step-8.cc b/tests/hp/step-8.cc new file mode 100644 index 0000000000..4744612bbb --- /dev/null +++ b/tests/hp/step-8.cc @@ -0,0 +1,468 @@ +//---------------------------- step-8.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005, 2006 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- step-8.cc --------------------------- + + +// a hp-ified version of step-8 + + +#include +#include +std::ofstream logfile("step-8/output"); + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + + + +template +class ElasticProblem +{ + public: + ElasticProblem (); + ~ElasticProblem (); + void run (); + + private: + void setup_system (); + void assemble_system (); + void solve (); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + hp::DoFHandler dof_handler; + + hp::FECollection fe; + + ConstraintMatrix hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + + +template +class RightHandSide : public Function +{ + public: + RightHandSide (); + + virtual void vector_value (const Point &p, + Vector &values) const; + + virtual void vector_value_list (const std::vector > &points, + std::vector > &value_list) const; +}; + + +template +RightHandSide::RightHandSide () + : + Function (dim) +{} + + +template +inline +void RightHandSide::vector_value (const Point &p, + Vector &values) const +{ + Assert (values.size() == dim, + ExcDimensionMismatch (values.size(), dim)); + Assert (dim >= 2, ExcNotImplemented()); + + Point point_1, point_2; + point_1(0) = 0.5; + point_2(0) = -0.5; + + if (((p-point_1).square() < 0.2*0.2) || + ((p-point_2).square() < 0.2*0.2)) + values(0) = 1; + else + values(0) = 0; + + if (p.square() < 0.2*0.2) + values(1) = 1; + else + values(1) = 0; +} + + + +template +void RightHandSide::vector_value_list (const std::vector > &points, + std::vector > &value_list) const +{ + Assert (value_list.size() == points.size(), + ExcDimensionMismatch (value_list.size(), points.size())); + + const unsigned int n_points = points.size(); + + for (unsigned int p=0; p::vector_value (points[p], + value_list[p]); +} + + + + + + +template +ElasticProblem::ElasticProblem () + : + dof_handler (triangulation), + fe (FESystem(FE_Q(1), dim)) +{} + + + + +template +ElasticProblem::~ElasticProblem () +{ + dof_handler.clear (); +} + + + +template +void ElasticProblem::setup_system () +{ + dof_handler.distribute_dofs (fe); + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, + hanging_node_constraints); + hanging_node_constraints.close (); + sparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); + + hanging_node_constraints.condense (sparsity_pattern); + + sparsity_pattern.compress(); + + system_matrix.reinit (sparsity_pattern); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); +} + + + +template +void ElasticProblem::assemble_system () +{ + hp::QCollection quadrature_formula(2); + + hp::FEValues x_fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_q_points | update_JxW_values); + + 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); + + std::vector local_dof_indices (dofs_per_cell); + + std::vector lambda_values (n_q_points); + std::vector mu_values (n_q_points); + + ConstantFunction lambda(1.), mu(1.); + + RightHandSide right_hand_side; + std::vector > rhs_values (n_q_points, + Vector(dim)); + + + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_matrix = 0; + cell_rhs = 0; + + x_fe_values.reinit (cell); + const FEValues &fe_values = x_fe_values.get_present_fe_values(); + + lambda.value_list (fe_values.get_quadrature_points(), lambda_values); + mu.value_list (fe_values.get_quadrature_points(), mu_values); + + right_hand_side.vector_value_list (fe_values.get_quadrature_points(), + rhs_values); + + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + for (unsigned int i=0; i boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(dim), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); +} + + + + +template +void ElasticProblem::solve () +{ + SolverControl solver_control (1000, 1e-12); + SolverCG<> cg (solver_control); + + PreconditionSSOR<> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + + hanging_node_constraints.distribute (solution); +} + + + +template +void ElasticProblem::refine_grid () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + typename FunctionMap::type neumann_boundary; + KellyErrorEstimator::estimate (dof_handler, + hp::QCollection(2), + neumann_boundary, + 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 (); +} + + + +template +void ElasticProblem::output_results (const unsigned int cycle) const +{ + std::string filename = "solution-"; + filename += ('0' + cycle); + Assert (cycle < 10, ExcInternalError()); + + filename += ".gmv"; + std::ofstream output (filename.c_str()); + + DataOut > data_out; + data_out.attach_dof_handler (dof_handler); + + + + std::vector solution_names; + switch (dim) + { + case 1: + solution_names.push_back ("displacement"); + break; + case 2: + solution_names.push_back ("x_displacement"); + solution_names.push_back ("y_displacement"); + break; + case 3: + solution_names.push_back ("x_displacement"); + solution_names.push_back ("y_displacement"); + solution_names.push_back ("z_displacement"); + break; + default: + Assert (false, ExcNotImplemented()); + } + + data_out.add_data_vector (solution, solution_names); + data_out.build_patches (); + data_out.write_gmv (output); +} + + + + +template +void ElasticProblem::run () +{ + for (unsigned int cycle=0; cycle<8; ++cycle) + { + deallog << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (2); + } + else + refine_grid (); + + deallog << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + setup_system (); + + deallog << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + assemble_system (); + solve (); + output_results (cycle); + } +} + + +int main () +{ + logfile.precision(2); + + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + deallog.depth_console (0); + + ElasticProblem<2> elastic_problem_2d; + elastic_problem_2d.run (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/hp/step-8/.cvsignore b/tests/hp/step-8/.cvsignore new file mode 100644 index 0000000000..d196dfb299 --- /dev/null +++ b/tests/hp/step-8/.cvsignore @@ -0,0 +1 @@ +obj.* exe OK output diff --git a/tests/hp/step-9.cc b/tests/hp/step-9.cc new file mode 100644 index 0000000000..c5db4b5e32 --- /dev/null +++ b/tests/hp/step-9.cc @@ -0,0 +1,761 @@ +//---------------------------- step-9.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005, 2006 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- step-9.cc --------------------------- + + +// a hp-ified version of step-9 + + +#include +#include +std::ofstream logfile("step-9/output"); + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include + +#include +#include + + + +#ifndef M_PI +# define M_PI 3.14159265358979323846 +#endif + + +template +class AdvectionProblem +{ + public: + AdvectionProblem (); + ~AdvectionProblem (); + void run (); + + private: + void setup_system (); + void assemble_system (); + void assemble_system_interval (const typename hp::DoFHandler::active_cell_iterator &begin, + const typename hp::DoFHandler::active_cell_iterator &end); + + void solve (); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + hp::DoFHandler dof_handler; + + hp::FECollection fe; + + ConstraintMatrix hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + + Threads::ThreadMutex assembler_lock; +}; + + + + +template +class AdvectionField : public TensorFunction<1,dim> +{ + public: + AdvectionField () : TensorFunction<1,dim> () {}; + + virtual Tensor<1,dim> value (const Point &p) const; + + virtual void value_list (const std::vector > &points, + std::vector > &values) const; + + DeclException2 (ExcDimensionMismatch, + unsigned int, unsigned int, + << "The vector has size " << arg1 << " but should have " + << arg2 << " elements."); +}; + + + +template +Tensor<1,dim> +AdvectionField::value (const Point &p) const +{ + Point value; + value[0] = 2; + for (unsigned int i=1; i +void +AdvectionField::value_list (const std::vector > &points, + std::vector > &values) const +{ + Assert (values.size() == points.size(), + ExcDimensionMismatch (values.size(), points.size())); + + for (unsigned int i=0; i::value (points[i]); +} + + + + +template +class RightHandSide : public Function +{ + public: + RightHandSide () : Function() {}; + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void value_list (const std::vector > &points, + std::vector &values, + const unsigned int component = 0) const; + + private: + static const Point center_point; +}; + + +template <> +const Point<1> RightHandSide<1>::center_point = Point<1> (-0.75); + +template <> +const Point<2> RightHandSide<2>::center_point = Point<2> (-0.75, -0.75); + +template <> +const Point<3> RightHandSide<3>::center_point = Point<3> (-0.75, -0.75, -0.75); + + + +template +double +RightHandSide::value (const Point &p, + const unsigned int component) const +{ + Assert (component == 0, ExcIndexRange (component, 0, 1)); + const double diameter = 0.1; + return ( (p-center_point).square() < diameter*diameter ? + .1/std::pow(diameter,dim) : + 0); +} + + + +template +void +RightHandSide::value_list (const std::vector > &points, + std::vector &values, + const unsigned int component) const +{ + Assert (values.size() == points.size(), + ExcDimensionMismatch (values.size(), points.size())); + + for (unsigned int i=0; i::value (points[i], component); +} + + + +template +class BoundaryValues : public Function +{ + public: + BoundaryValues () : Function() {}; + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void value_list (const std::vector > &points, + std::vector &values, + const unsigned int component = 0) const; +}; + + + +template +double +BoundaryValues::value (const Point &p, + const unsigned int component) const +{ + Assert (component == 0, ExcIndexRange (component, 0, 1)); + + const double sine_term = std::sin(16*M_PI*std::sqrt(p.square())); + const double weight = std::exp(-5*p.square()) / std::exp(-5.); + return sine_term * weight; +} + + + +template +void +BoundaryValues::value_list (const std::vector > &points, + std::vector &values, + const unsigned int component) const +{ + Assert (values.size() == points.size(), + ExcDimensionMismatch (values.size(), points.size())); + + for (unsigned int i=0; i::value (points[i], component); +} + + + + +class GradientEstimation +{ + public: + template + static void estimate (const hp::DoFHandler &dof, + const Vector &solution, + Vector &error_per_cell); + + DeclException2 (ExcInvalidVectorLength, + int, int, + << "Vector has length " << arg1 << ", but should have " + << arg2); + DeclException0 (ExcInsufficientDirections); + + private: + typedef std::pair IndexInterval; + + template + static void estimate_interval (const hp::DoFHandler &dof, + const Vector &solution, + const IndexInterval &index_interval, + Vector &error_per_cell); +}; + + + + + +template +AdvectionProblem::AdvectionProblem () : + dof_handler (triangulation), + fe(FE_Q(1)) +{} + + + +template +AdvectionProblem::~AdvectionProblem () +{ + dof_handler.clear (); +} + + + +template +void AdvectionProblem::setup_system () +{ + dof_handler.distribute_dofs (fe); + + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, + hanging_node_constraints); + hanging_node_constraints.close (); + + sparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); + + hanging_node_constraints.condense (sparsity_pattern); + + sparsity_pattern.compress(); + + system_matrix.reinit (sparsity_pattern); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); +} + + + +template +void AdvectionProblem::assemble_system () +{ + const unsigned int n_threads = multithread_info.n_default_threads; + Threads::ThreadGroup<> threads; + + typedef typename hp::DoFHandler::active_cell_iterator active_cell_iterator; + std::vector > + thread_ranges + = Threads::split_range (dof_handler.begin_active (), + dof_handler.end (), + n_threads); + + for (unsigned int thread=0; thread::assemble_system_interval) + (thread_ranges[thread].first, + thread_ranges[thread].second); + + threads.join_all (); + + + hanging_node_constraints.condense (system_matrix); + hanging_node_constraints.condense (system_rhs); +} + + + +template +void +AdvectionProblem:: +assemble_system_interval (const typename hp::DoFHandler::active_cell_iterator &begin, + const typename hp::DoFHandler::active_cell_iterator &end) +{ + const AdvectionField advection_field; + const RightHandSide right_hand_side; + const BoundaryValues boundary_values; + + hp::QCollection quadrature_formula(QGauss(2)); + hp::QCollection face_quadrature_formula(QGauss(2)); + + hp::FEValues x_fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_q_points | update_JxW_values); + hp::FEFaceValues x_fe_face_values (fe, face_quadrature_formula, + update_values | update_q_points | + update_JxW_values | update_normal_vectors); + + const unsigned int dofs_per_cell = fe[0].dofs_per_cell; + const unsigned int n_q_points = quadrature_formula[0].n_quadrature_points; + const unsigned int n_face_q_points = face_quadrature_formula[0].n_quadrature_points; + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + std::vector rhs_values (n_q_points); + std::vector > advection_directions (n_q_points); + std::vector face_boundary_values (n_face_q_points); + std::vector > face_advection_directions (n_face_q_points); + + typename hp::DoFHandler::active_cell_iterator cell; + for (cell=begin; cell!=end; ++cell) + { + cell_matrix = 0; + cell_rhs = 0; + + x_fe_values.reinit (cell); + const FEValues &fe_values = x_fe_values.get_present_fe_values(); + + advection_field.value_list (fe_values.get_quadrature_points(), + advection_directions); + right_hand_side.value_list (fe_values.get_quadrature_points(), + rhs_values); + + const double delta = 0.1 * cell->diameter (); + + for (unsigned int q_point=0; q_point::faces_per_cell; ++face) + if (cell->face(face)->at_boundary()) + { + x_fe_face_values.reinit (cell, face); + const FEFaceValues &fe_face_values + = x_fe_face_values.get_present_fe_values (); + + boundary_values.value_list (fe_face_values.get_quadrature_points(), + face_boundary_values); + advection_field.value_list (fe_face_values.get_quadrature_points(), + face_advection_directions); + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + + assembler_lock.acquire (); + for (unsigned int i=0; i +void AdvectionProblem::solve () +{ + SolverControl solver_control (1000, 1e-12); + SolverBicgstab<> bicgstab (solver_control); + + PreconditionJacobi<> preconditioner; + preconditioner.initialize(system_matrix, 1.0); + + bicgstab.solve (system_matrix, solution, system_rhs, + preconditioner); + + hanging_node_constraints.distribute (solution); +} + + +template +void AdvectionProblem::refine_grid () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + GradientEstimation::estimate (dof_handler, + solution, + estimated_error_per_cell); + + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.5, 0.03); + + triangulation.execute_coarsening_and_refinement (); +} + + + +template +void AdvectionProblem::output_results (const unsigned int cycle) const +{ + std::string filename = "grid-"; + filename += ('0' + cycle); + Assert (cycle < 10, ExcInternalError()); + + filename += ".eps"; + std::ofstream output (filename.c_str()); + + GridOut grid_out; + grid_out.write_eps (triangulation, output); +} + + +template +void AdvectionProblem::run () +{ + for (unsigned int cycle=0; cycle<6; ++cycle) + { + deallog << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (4); + } + else + { + refine_grid (); + }; + + + deallog << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + setup_system (); + + deallog << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + assemble_system (); + solve (); + output_results (cycle); + }; + + DataOut > data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "solution"); + data_out.build_patches (); + + std::ofstream output ("final-solution.gmv"); + data_out.write_gmv (output); +} + + + + +template +void +GradientEstimation::estimate (const hp::DoFHandler &dof_handler, + const Vector &solution, + Vector &error_per_cell) +{ + Assert (error_per_cell.size() == dof_handler.get_tria().n_active_cells(), + ExcInvalidVectorLength (error_per_cell.size(), + dof_handler.get_tria().n_active_cells())); + + const unsigned int n_threads = multithread_info.n_default_threads; + std::vector index_intervals + = Threads::split_interval (0, dof_handler.get_tria().n_active_cells(), + n_threads); + + Threads::ThreadGroup<> threads; + void (*estimate_interval_ptr) (const hp::DoFHandler &, + const Vector &, + const IndexInterval &, + Vector &) + = &GradientEstimation::template estimate_interval; + for (unsigned int i=0; i +void +GradientEstimation::estimate_interval (const hp::DoFHandler &dof_handler, + const Vector &solution, + const IndexInterval &index_interval, + Vector &error_per_cell) +{ + QMidpoint xmidpoint; + hp::QCollection midpoint_rule (xmidpoint); + hp::FEValues x_fe_midpoint_value (dof_handler.get_fe(), + midpoint_rule, + update_values | update_q_points); + + Tensor<2,dim> Y; + + typename hp::DoFHandler::active_cell_iterator cell, endc; + + cell = dof_handler.begin_active(); + advance (cell, static_cast(index_interval.first)); + + endc = dof_handler.begin_active(); + advance (endc, static_cast(index_interval.second)); + + Vector::iterator + error_on_this_cell = error_per_cell.begin() + index_interval.first; + + + std::vector::active_cell_iterator> active_neighbors; + active_neighbors.reserve (GeometryInfo::faces_per_cell * + GeometryInfo::subfaces_per_face); + + for (; cell!=endc; ++cell, ++error_on_this_cell) + { + x_fe_midpoint_value.reinit (cell); + const FEValues &fe_midpoint_value = x_fe_midpoint_value.get_present_fe_values (); + + Y.clear (); + + Tensor<1,dim> projected_gradient; + + + active_neighbors.clear (); + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + if (! cell->at_boundary(face_no)) + { + const typename hp::DoFHandler::face_iterator + face = cell->face(face_no); + const typename hp::DoFHandler::cell_iterator + neighbor = cell->neighbor(face_no); + + if (neighbor->active()) + active_neighbors.push_back (neighbor); + else + { + if (dim == 1) + { + typename hp::DoFHandler::cell_iterator + neighbor_child = neighbor; + while (neighbor_child->has_children()) + neighbor_child = neighbor_child->child (face_no==0 ? 1 : 0); + + Assert (neighbor_child->neighbor(face_no==0 ? 1 : 0)==cell, + ExcInternalError()); + + active_neighbors.push_back (neighbor_child); + } + else + for (unsigned int subface_no=0; subface_non_children(); ++subface_no) + active_neighbors.push_back ( + cell->neighbor_child_on_subface(face_no, subface_no)); + }; + }; + + const Point this_center = fe_midpoint_value.quadrature_point(0); + + std::vector this_midpoint_value(1); + fe_midpoint_value.get_function_values (solution, this_midpoint_value); + + + std::vector neighbor_midpoint_value(1); + typename std::vector::active_cell_iterator>::const_iterator + neighbor_ptr = active_neighbors.begin(); + for (; neighbor_ptr!=active_neighbors.end(); ++neighbor_ptr) + { + const typename hp::DoFHandler::active_cell_iterator + neighbor = *neighbor_ptr; + + x_fe_midpoint_value.reinit (neighbor); + const FEValues &fe_midpoint_value = x_fe_midpoint_value.get_present_fe_values (); + const Point neighbor_center = fe_midpoint_value.quadrature_point(0); + + fe_midpoint_value.get_function_values (solution, + neighbor_midpoint_value); + + Point y = neighbor_center - this_center; + const double distance = std::sqrt(y.square()); + y /= distance; + + for (unsigned int i=0; i Y_inverse = invert(Y); + + Point gradient; + contract (gradient, Y_inverse, projected_gradient); + + *error_on_this_cell = (std::pow(cell->diameter(), + 1+1.0*dim/2) * + std::sqrt(gradient.square())); + }; +} + + + +int main () +{ + logfile.precision(2); + + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + deallog.depth_console (0); + + AdvectionProblem<2> advection_problem_2d; + advection_problem_2d.run (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; + + return 0; +} diff --git a/tests/hp/step-9/.cvsignore b/tests/hp/step-9/.cvsignore new file mode 100644 index 0000000000..d196dfb299 --- /dev/null +++ b/tests/hp/step-9/.cvsignore @@ -0,0 +1 @@ +obj.* exe OK output