Deal all,
I am a beginner for deal.II, and have studied some examples in tutorial.
I tried to solve the laplace equation in 1D, but I could not get the correct
reuslt.
Especially, the right boundary value doesn`t conform to the value in boundary
function.
The source code is put in the attached file.
Thanks in advance.
Best
Feng Cui
#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/compressed_sparsity_pattern.h>
#include <lac/solver_cg.h>
#include <lac/precondition.h>
#include <numerics/data_out.h>
#include <fstream>
#include <iostream>
#include <base/logstream.h>
using namespace dealii;
template <int dim>
class LaplaceProblem
{
public:
LaplaceProblem ();
void run ();
private:
void make_grid_and_dofs ();
void assemble_system ();
void solve ();
void output_results () const;
Triangulation<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
template <int dim>
class RightHandSide : public Function<dim>
{
public:
RightHandSide () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
BoundaryValues () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double RightHandSide<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
{
// double return_value = 0;
// for (unsigned int i=0; i<dim; ++i)
// return_value += 4*std::pow(p(i), 4);
return 4*sin(2*p[0]);
}
template <int dim>
double BoundaryValues<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
{
std::cout << "===" << p[0] << std::endl;
return sin(2*p[0]);
}
template <int dim>
LaplaceProblem<dim>::LaplaceProblem () :
fe (1),
dof_handler (triangulation)
{}
template <int dim>
void LaplaceProblem<dim>::make_grid_and_dofs ()
{
GridGenerator::hyper_cube (triangulation, 0, 1);
triangulation.refine_global (8);
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl
<< " 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;
CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern (dof_handler, c_sparsity);
sparsity_pattern.copy_from(c_sparsity);
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 ()
{
// QGauss<dim> quadrature_formula(2);
//
// const RightHandSide<dim> right_hand_side;
//
// FEValues<dim> fe_values (fe, quadrature_formula,
// update_values | update_gradients |
// update_quadrature_points | 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);
//
// typename DoFHandler<dim>::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 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) += (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) *
// right_hand_side.value
(fe_values.quadrature_point (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));
//
// system_rhs(local_dof_indices[i]) += cell_rhs(i);
// }
// }
MatrixCreator::create_laplace_matrix(dof_handler, QGauss<dim>(3),
system_matrix);
// system_matrix *= -1;
RightHandSide<dim> rhs_function;
VectorTools::create_right_hand_side(dof_handler, QGauss<dim>(2),
rhs_function, system_rhs);
std::map<unsigned int,double> boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
0,
BoundaryValues<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);
cg.solve (system_matrix, solution, system_rhs,
PreconditionIdentity());
std::cout << " " << solver_control.last_step()
<< " CG iterations needed to obtain convergence."
<< std::endl;
}
template <int dim>
void LaplaceProblem<dim>::output_results () const
{
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, "solution");
data_out.build_patches ();
std::ofstream output (dim == 1 ?
"solution-1d.gnuplot" :
"solution-2d.gnuplot");
data_out.write_gnuplot (output);
}
template <int dim>
void LaplaceProblem<dim>::run ()
{
std::cout << "Solving problem in " << dim << " space dimensions." <<
std::endl;
make_grid_and_dofs();
assemble_system ();
solve ();
output_results ();
}
int main ()
{
deallog.depth_console (0);
{
LaplaceProblem<1> laplace_problem_1d;
laplace_problem_1d.run ();
}
// {
// LaplaceProblem<2> laplace_problem_2d;
// laplace_problem_2d.run ();
// }
//
// {
// LaplaceProblem<3> laplace_problem_3d;
// laplace_problem_3d.run ();
// }
return 0;
}
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
