OK, here are two small example programs - one with hp and one without.
Am 04.05.10 15:27, schrieb Wolfgang Bangerth:
Is there something special in 1D? Because in higher dimensions
everything is working fine.
Not that I would know of. It may be a bug, but it's hard to tell without
more information.
W.
-------------------------------------------------------------------------
Wolfgang Bangerth email: [email protected]
www: http://www.math.tamu.edu/~bangerth/
#include <grid/tria.h>
#include <hp/dof_handler.h>
#include <grid/grid_generator.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <dofs/dof_accessor.h>
#include <hp/fe_collection.h>
#include <fe/fe_q.h>
#include <dofs/dof_tools.h>
#include <hp/fe_values.h>
#include <hp/q_collection.h>
#include <base/quadrature_lib.h>
#include <base/function.h>
#include <numerics/vectors.h>
#include <numerics/matrices.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 <numerics/data_out.h>
#include <fstream>
#include <iostream>
using namespace dealii;
class ExactSolution: public Function<1> {
public:
ExactSolution ();
virtual double value (const Point<1>& p, const unsigned int component = 0)
const;
};
class LaplaceProblem
{
public:
LaplaceProblem ();
void run ();
private:
const ExactSolution exact_solution;
void make_grid_and_dofs ();
void assemble_system ();
void solve ();
void output_results () const;
Triangulation<1> triangulation;
hp::FECollection<1> fe;
hp::DoFHandler<1> dof_handler;
hp::QCollection<1> quadrature;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
double ExactSolution::value (const Point<1>& p, const unsigned int) const {
return p (0) * p (0);
}
ExactSolution::ExactSolution (): Function<1> () {
}
LaplaceProblem::LaplaceProblem () :
dof_handler (triangulation)
{
fe.push_back (FE_Q<1> (1));
quadrature.push_back (QGauss<1> (2));
}
void LaplaceProblem::make_grid_and_dofs ()
{
GridGenerator::hyper_cube (triangulation);
triangulation.refine_global (2);
std::cout << "Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl;
std::cout << "Total number of cells: "
<< triangulation.n_cells()
<< std::endl;
dof_handler.distribute_dofs (fe);
std::cout << "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());
}
void LaplaceProblem::assemble_system ()
{
hp::FEValues<1> fe_values (fe, quadrature,
update_values | update_gradients | update_JxW_values);
unsigned int dofs_per_cell;
unsigned int n_q_points;
FullMatrix<double> cell_matrix;
Vector<double> cell_rhs;
std::vector<unsigned int> local_dof_indices;
hp::DoFHandler<1>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
dofs_per_cell = cell->get_fe ().dofs_per_cell;
cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
cell_matrix = 0;
cell_rhs.reinit (dofs_per_cell);
cell_rhs = 0;
n_q_points = fe_values.get_present_fe_values ().n_quadrature_points;
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
cell_matrix(i,j) += (fe_values.get_present_fe_values ().shape_grad
(i, q_point) *
fe_values.get_present_fe_values ().shape_grad
(j, q_point) *
fe_values.get_present_fe_values ().JxW
(q_point));
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
cell_rhs(i) += (-2.0 * fe_values.get_present_fe_values ().shape_value
(i, q_point) *
fe_values.get_present_fe_values ().JxW (q_point));
local_dof_indices.resize (dofs_per_cell);
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));
for (unsigned int i=0; i<dofs_per_cell; ++i)
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
std::map<unsigned int,double> boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
0,
exact_solution,
boundary_values);
VectorTools::interpolate_boundary_values (dof_handler,
1,
exact_solution,
boundary_values);
MatrixTools::apply_boundary_values (boundary_values,
system_matrix,
solution,
system_rhs);
}
void LaplaceProblem::solve ()
{
SolverControl solver_control (1000, 1e-12);
SolverCG<> cg (solver_control);
cg.solve (system_matrix, solution, system_rhs,
PreconditionIdentity());
}
void LaplaceProblem::output_results () const
{
DataOut<1, hp::DoFHandler<1> > data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, "solution");
data_out.build_patches ();
std::ofstream output ("solution.gpl");
data_out.write_gnuplot (output);
for (unsigned int dof = 0; dof < dof_handler.n_dofs (); ++dof)
std::cout << solution (dof) << std::endl;
}
void LaplaceProblem::run ()
{
make_grid_and_dofs ();
assemble_system ();
solve ();
output_results ();
}
int main ()
{
LaplaceProblem laplace_problem;
laplace_problem.run ();
return 0;
}
#include <grid/tria.h>
#include <dofs/dof_handler.h>
#include <grid/grid_generator.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <dofs/dof_accessor.h>
#include <fe/fe_q.h>
#include <dofs/dof_tools.h>
#include <fe/fe_values.h>
#include <base/quadrature_lib.h>
#include <base/function.h>
#include <numerics/vectors.h>
#include <numerics/matrices.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 <numerics/data_out.h>
#include <fstream>
#include <iostream>
using namespace dealii;
class ExactSolution: public Function<1> {
public:
ExactSolution ();
virtual double value (const Point<1>& p, const unsigned int component = 0)
const;
};
class LaplaceProblem
{
public:
LaplaceProblem ();
void run ();
private:
const ExactSolution exact_solution;
void make_grid_and_dofs ();
void assemble_system ();
void solve ();
void output_results () const;
Triangulation<1> triangulation;
FE_Q<1> fe;
DoFHandler<1> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
double ExactSolution::value (const Point<1>& p, const unsigned int) const {
return p (0) * p (0);
}
ExactSolution::ExactSolution (): Function<1> () {
}
LaplaceProblem::LaplaceProblem () :
fe (1),
dof_handler (triangulation)
{}
void LaplaceProblem::make_grid_and_dofs ()
{
GridGenerator::hyper_cube (triangulation);
triangulation.refine_global (2);
std::cout << "Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl;
std::cout << "Total number of cells: "
<< triangulation.n_cells()
<< std::endl;
dof_handler.distribute_dofs (fe);
std::cout << "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());
}
void LaplaceProblem::assemble_system ()
{
QGauss<1> quadrature_formula(2);
FEValues<1> fe_values (fe, quadrature_formula,
update_values | update_gradients | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
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);
DoFHandler<1>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
cell_matrix = 0;
cell_rhs = 0;
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
cell_matrix(i,j) += (fe_values.shape_grad (i, q_point) *
fe_values.shape_grad (j, q_point) *
fe_values.JxW (q_point));
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
cell_rhs(i) += (-2.0 * fe_values.shape_value (i, q_point) *
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));
for (unsigned int i=0; i<dofs_per_cell; ++i)
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
std::map<unsigned int,double> boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
0,
exact_solution,
boundary_values);
VectorTools::interpolate_boundary_values (dof_handler,
1,
exact_solution,
boundary_values);
MatrixTools::apply_boundary_values (boundary_values,
system_matrix,
solution,
system_rhs);
}
void LaplaceProblem::solve ()
{
SolverControl solver_control (1000, 1e-12);
SolverCG<> cg (solver_control);
cg.solve (system_matrix, solution, system_rhs,
PreconditionIdentity());
}
void LaplaceProblem::output_results () const
{
DataOut<1> data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, "solution");
data_out.build_patches ();
std::ofstream output ("solution.gpl");
data_out.write_gnuplot (output);
}
void LaplaceProblem::run ()
{
make_grid_and_dofs ();
assemble_system ();
solve ();
output_results ();
}
int main ()
{
LaplaceProblem laplace_problem;
laplace_problem.run ();
return 0;
}
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
