Hi all,
I was following tutorial 22 and 34 to implement inhomogeneous Dirichlet
boundary conditions using the ConsraintMatrix class but couldn't get it
right. I think the problem mainly lies in the assemblage part (i was
solving the basic Poisson's eq.)
if(constraints.is_inhomogeneously_constrained
(local_dof_indices[i])){
for(unsigned int j=0; j <dofs_per_cell; ++j){
matrix_for_bc(j,i) +=
grad_basis_phi[i] * grad_basis_phi[j] *
fe_values.JxW(q);
}//j-loop
}//if
...
constraints.distribute_local_to_global(local_rhs,
local_dof_indices,
rhs,
matrix_for_bc);
However, when I comment out the if statement, the program seems to work
correctly. The code is attached in this email.
Your help will be greatly appreciated.
Thanks!
Huan
#include <base/quadrature_lib.h>
#include <base/logstream.h>
#include <base/utilities.h>
#include <base/timer.h>
#include <lac/full_matrix.h>
#include <lac/sparse_matrix.h>
#include <lac/solver_gmres.h>
#include <lac/solver_cg.h>
#include <lac/block_sparsity_pattern.h>
#include <lac/constraint_matrix.h>
#include <lac/precondition.h>
#include <grid/tria.h>
#include <grid/grid_generator.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/grid_tools.h>
#include <grid/grid_refinement.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_renumbering.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
#include <fe/fe_q.h>
#include <fe/fe_system.h>
#include <fe/fe_values.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
#include <numerics/data_out.h>
#include <numerics/error_estimator.h>
#include <numerics/solution_transfer.h>
#include <lac/trilinos_sparse_matrix.h>
#include <lac/trilinos_block_sparse_matrix.h>
#include <lac/trilinos_vector.h>
#include <lac/trilinos_block_vector.h>
#include <lac/trilinos_precondition.h>
#include <fstream>
#include <sstream>
#include<cstdlib>
using namespace dealii;
namespace EquationData{
//geometry info
const double xl = 0.0;
const double xr = numbers::PI;
const double yb = 0.0;
const double yt = numbers::PI;
// 2d version
const Point<2> left_bottom(xl, yb);
const Point<2> right_top(xr, yt);
//solution
template<int dim>
class ExactSolution : public Function<dim>{
public:
ExactSolution() : Function<dim>(1){}
virtual double value(const Point<dim> &p,
const unsigned int component = 0) const;
};
template<int dim>
double ExactSolution<dim>::value(const Point<dim> &p,
const unsigned int component)const{
return std::sin(p[0]);
}//value()
template<int dim>
class RightHandSide : public Function<dim>{
public:
RightHandSide() : Function<dim>(1){}
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{
return std::sin(p[0]);
}//value()
}//namespace EquationData
template<int dim>
class PoissonProblem{
public:
PoissonProblem();
void run();
private:
void setup_dofs();
void assemble_system();
void assemble_stiffness_matrix();
void solve();
void compute_errors() const;
Triangulation<dim> triangulation;
const unsigned int degree;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
ConstraintMatrix constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> stiffness_matrix;
Vector<double> solution;
Vector<double> rhs;
};//class PoissonProblem
template<int dim>
PoissonProblem<dim>::PoissonProblem()
:
triangulation(Triangulation<dim>::maximum_smoothing),
degree(1),
fe(degree),
dof_handler(triangulation)
{}
template<int dim>
void PoissonProblem<dim>::setup_dofs(){
{
dof_handler.distribute_dofs(fe);
constraints.clear();
//let's assume dirichlet boundary condition.
VectorTools::interpolate_boundary_values
(dof_handler,
0,
EquationData::ExactSolution<dim>(),
constraints);
DoFTools::make_hanging_node_constraints(dof_handler,
constraints);
constraints.close();
}//localize
std::cout << "Number of active cells: "
<< triangulation.n_active_cells()
<< " ( on "
<< triangulation.n_levels()
<< " levels )"
<< std::endl
<< "Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl << std::endl;
{
stiffness_matrix.clear();
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,
constraints, false);
sparsity_pattern.compress();
stiffness_matrix.reinit(sparsity_pattern);
}//localize
solution.reinit(dof_handler.n_dofs());
rhs.reinit(dof_handler.n_dofs());
}//setup_dofs()
template<int dim>
void
PoissonProblem<dim>::assemble_system(){
stiffness_matrix = 0;
rhs = 0;
const QGauss<dim> quadrature_formula(degree + 2);
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> local_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> local_rhs(dofs_per_cell);
std::vector<unsigned int> local_dof_indices(dofs_per_cell);
std::vector<double> basis_phi(dofs_per_cell);
std::vector<Tensor<1,dim> > grad_basis_phi(dofs_per_cell);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
EquationData::RightHandSide<dim> rhs_function;
FullMatrix<double> matrix_for_bc(dofs_per_cell,
dofs_per_cell);
for(; cell != endc; ++cell){
local_matrix = 0;
local_rhs = 0;
matrix_for_bc = 0;
fe_values.reinit(cell);
for(unsigned int q=0; q<n_q_points;++q){
for(unsigned int k=0; k<dofs_per_cell; ++k){
basis_phi[k] = fe_values.shape_value(k,q);
grad_basis_phi[k] = fe_values.shape_grad(k,q);
}//k-loop
for(unsigned int i=0; i<dofs_per_cell;++i){
for(unsigned int j =0; j < dofs_per_cell; ++j){
local_matrix(i,j) +=
grad_basis_phi[i] * grad_basis_phi[j] * fe_values.JxW(q);
}//j-loop
local_rhs(i) +=
rhs_function.value(fe_values.quadrature_point(q))
* basis_phi[i]* fe_values.JxW(q);
if(constraints.is_inhomogeneously_constrained
(local_dof_indices[i])){
for(unsigned int j=0; j <dofs_per_cell; ++j){
matrix_for_bc(j,i) +=
grad_basis_phi[i] * grad_basis_phi[j] *
fe_values.JxW(q);
}//j-loop
}//if
}//i-loop
}//q-loop
cell->get_dof_indices(local_dof_indices);
// constraints.distribute_local_to_global(local_matrix,
// local_rhs,
// local_dof_indices,
// stiffness_matrix,
// rhs);
constraints.distribute_local_to_global(local_matrix,
local_dof_indices,
stiffness_matrix);
constraints.distribute_local_to_global(local_rhs,
local_dof_indices,
rhs,
matrix_for_bc);
}//cell-loop
}//assemble_system()
template<int dim>
void PoissonProblem<dim>::solve(){
//one update in Newton's method
std::cout << " Solving ..." << std::endl;
assemble_system();
{
SolverControl solver_control(stiffness_matrix.m(),
1e-12*rhs.l2_norm());
SolverCG<> cg(solver_control);
// for(unsigned int i = 0; i < solution.size(); ++i){
// if(constraints.is_constrained(i))
// solution(i) = 0;
// }//i-loop
cg.solve(stiffness_matrix, solution,
rhs, PreconditionIdentity());
constraints.distribute(solution);
std::cout << " "
<< solver_control.last_step()
<< " CG iterations for solution."
<< std::endl;
double min_phi = solution(0),
max_phi = solution(0);
for(unsigned int i=0; i < solution.size(); ++i){
min_phi = std::min<double>(min_phi, solution(i));
max_phi = std::max<double>(max_phi, solution(i));
}//i-loop
std::cout << " solution range: "
<< min_phi << ' ' << max_phi
<< std::endl;
}//localize
}//solve()
template<int dim>
void PoissonProblem<dim>::compute_errors()const{
EquationData::ExactSolution<dim> exact_solution;
Vector<double> cellwise_errors(triangulation.n_active_cells());
QTrapez<1> q_trapez;
QIterated<dim> quadrature(q_trapez, degree+2);
VectorTools::integrate_difference(dof_handler,
solution,
exact_solution,
cellwise_errors,
quadrature,
VectorTools::L2_norm);
const double L2_error = cellwise_errors.l2_norm();
std::cout << "Errors: " << L2_error << std::endl;
}//compute_errors()
template<int dim>
void PoissonProblem<dim>::run(){
const unsigned int initial_refinement = 2;
GridGenerator::hyper_rectangle(triangulation,
EquationData::left_bottom,
EquationData::right_top);
start_time_iteration:
unsigned int step_number = 0;
triangulation.refine_global(initial_refinement);
double time = 0;
do{//time marching
triangulation.refine_global(1);
setup_dofs();
solve();
compute_errors();
std::cout << std::endl;
++step_number;
}while(step_number < 4);
std::cout << "job done!" << std::endl;
}//run()
int main(int argc, char* argv[]){
try{
deallog.depth_console(0);
PoissonProblem<2> problem;
problem.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;
}
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
