--- /dev/null
+//---------------------------- step-5.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-5.cc ---------------------------
+
+
+// a hp-ified version of step-5
+
+
+#include <base/logstream.h>
+#include <fstream>
+std::ofstream logfile("step-5/output");
+
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <dofs/hp_dof_handler.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_q.h>
+#include <fe/hp_fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+
+#include <grid/grid_in.h>
+
+#include <grid/tria_boundary_lib.h>
+
+#include <fstream>
+#include <sstream>
+
+
+
+template <int dim>
+class LaplaceProblem
+{
+ public:
+ LaplaceProblem ();
+ void run ();
+
+ private:
+ void setup_system ();
+ void assemble_system ();
+ void solve ();
+ void output_results (const unsigned int cycle) const;
+
+ Triangulation<dim> triangulation;
+ hp::FECollection<dim> fe;
+ hp::DoFHandler<dim> dof_handler;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+};
+
+
+
+template <int dim>
+class Coefficient : public Function<dim>
+{
+ public:
+ Coefficient () : Function<dim>() {};
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ virtual void value_list (const std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+};
+
+
+
+template <int dim>
+double Coefficient<dim>::value (const Point<dim> &p,
+ const unsigned int /*component*/) const
+{
+ if (p.square() < 0.5*0.5)
+ return 20;
+ else
+ return 1;
+}
+
+
+
+
+template <int dim>
+void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component) const
+{
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch (values.size(), points.size()));
+
+ Assert (component == 0,
+ ExcIndexRange (component, 0, 1));
+
+ const unsigned int n_points = points.size();
+
+ for (unsigned int i=0; i<n_points; ++i)
+ {
+ if (points[i].square() < 0.5*0.5)
+ values[i] = 20;
+ else
+ values[i] = 1;
+ }
+}
+
+
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem () :
+ fe (FE_Q<dim>(1)),
+ dof_handler (triangulation)
+{}
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+ dof_handler.distribute_dofs (fe);
+
+ deallog << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << std::endl;
+
+ 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);
+ sparsity_pattern.compress();
+
+ system_matrix.reinit (sparsity_pattern);
+
+ solution.reinit (dof_handler.n_dofs());
+ system_rhs.reinit (dof_handler.n_dofs());
+}
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_system ()
+{
+ hp::QCollection<dim> quadrature_formula(QGauss<dim>(2));
+
+ hp::FEValues<dim> 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<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs (dofs_per_cell);
+
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+ const Coefficient<dim> coefficient;
+ std::vector<double> coefficient_values (n_q_points);
+
+ typename hp::DoFHandler<dim>::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<2> &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_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += (coefficient_values[q_point] *
+ fe_values.shape_grad(i,q_point) *
+ fe_values.shape_grad(j,q_point) *
+ fe_values.JxW(q_point));
+
+ cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+ 1.0 *
+ fe_values.JxW(q_point));
+ }
+
+
+ cell->get_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ system_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i,j));
+
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ }
+ }
+
+ std::map<unsigned int,double> 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 <int dim>
+void LaplaceProblem<dim>::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);
+
+ deallog << " " << solver_control.last_step()
+ << " CG iterations needed to obtain convergence."
+ << std::endl;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
+{
+ DataOut<dim,hp::DoFHandler<dim> > data_out;
+
+ data_out.attach_dof_handler (dof_handler);
+ data_out.add_data_vector (solution, "solution");
+
+ data_out.build_patches ();
+
+ DataOutBase::EpsFlags eps_flags;
+ eps_flags.z_scaling = 4;
+ eps_flags.azimut_angle = 40;
+ eps_flags.turn_angle = 10;
+ data_out.set_flags (eps_flags);
+
+ data_out.write_eps (deallog.get_file_stream());
+}
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run ()
+{
+ for (unsigned int cycle=0; cycle<6; ++cycle)
+ {
+ deallog << "Cycle " << cycle << ':' << std::endl;
+
+ if (cycle != 0)
+ triangulation.refine_global (1);
+ else
+ {
+ GridIn<dim> grid_in;
+ grid_in.attach_triangulation (triangulation);
+ std::ifstream input_file("step-5/circle-grid.inp");
+ Assert (dim==2, ExcInternalError());
+
+ grid_in.read_ucd (input_file);
+
+ static const HyperBallBoundary<dim> boundary;
+ triangulation.set_boundary (0, boundary);
+ }
+
+ deallog << " Number of active cells: "
+ << triangulation.n_active_cells()
+ << std::endl
+ << " Total number of cells: "
+ << triangulation.n_cells()
+ << std::endl;
+
+ setup_system ();
+ assemble_system ();
+ solve ();
+ output_results (cycle);
+ }
+}
+
+
+
+int main ()
+{
+ logfile.precision(2);
+
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ deallog.depth_console (0);
+
+ LaplaceProblem<2> laplace_problem_2d;
+ laplace_problem_2d.run ();
+
+/*
+ Coefficient<2> coefficient;
+ std::vector<Point<2> > points (2);
+ std::vector<double> coefficient_values (1);
+ coefficient.value_list (points, coefficient_values);
+*/
+
+ return 0;
+}
--- /dev/null
+
+DEAL::Cycle 0:
+DEAL:: Number of active cells: 20
+DEAL:: Total number of cells: 20
+DEAL:: Number of degrees of freedom: 25
+DEAL:cg::Starting value 0.51
+DEAL:cg::Convergence step 13 value 0
+DEAL:: 13 CG iterations needed to obtain convergence.
+DEAL::Cycle 1:
+DEAL:: Number of active cells: 80
+DEAL:: Total number of cells: 100
+DEAL:: Number of degrees of freedom: 89
+DEAL:cg::Starting value 0.31
+DEAL:cg::Convergence step 18 value 0
+DEAL:: 18 CG iterations needed to obtain convergence.
+DEAL::Cycle 2:
+DEAL:: Number of active cells: 320
+DEAL:: Total number of cells: 420
+DEAL:: Number of degrees of freedom: 337
+DEAL:cg::Starting value 0.17
+DEAL:cg::Convergence step 29 value 0
+DEAL:: 29 CG iterations needed to obtain convergence.
+DEAL::Cycle 3:
+DEAL:: Number of active cells: 1280
+DEAL:: Total number of cells: 1700
+DEAL:: Number of degrees of freedom: 1313
+DEAL:cg::Starting value 0.092
+DEAL:cg::Convergence step 52 value 0
+DEAL:: 52 CG iterations needed to obtain convergence.
+DEAL::Cycle 4:
+DEAL:: Number of active cells: 5120
+DEAL:: Total number of cells: 6820
+DEAL:: Number of degrees of freedom: 5185
+DEAL:cg::Starting value 0.047
+DEAL:cg::Convergence step 95 value 0
+DEAL:: 95 CG iterations needed to obtain convergence.
+DEAL::Cycle 5:
+DEAL:: Number of active cells: 20480
+DEAL:: Total number of cells: 27300
+DEAL:: Number of degrees of freedom: 20609
+DEAL:cg::Starting value 0.024
+DEAL:cg::Convergence step 182 value 0
+DEAL:: 182 CG iterations needed to obtain convergence.
projects / dealii.git / commitdiff