Hi
My main problem is the application of the Neumann B.C in 1D. Initially, I
thought this was all fine but compared to the exact solution, it's far from it.
Below is the assemble_system function I'm using. I'm not applying the Neumann
B.C properly. So, could you please help look over the code and point out what
I'm doing wrong.
The governing differential equation is:
d/dx (x * du/dx) = 2/x^2 on the domain 1 < x < 2
subject to these boundary conditions:
u(1) = 2
(-x * du/dx) = 0.5 @ x = 2
The assemble function:
template <int dim>
void LaplaceProblem<dim>::assemble_system ()
{
QGauss<dim> quadrature_formula(3);
//QGauss<dim-1> face_quadrature_formula(3);
const RightHandSide<dim> right_hand_side;
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
// FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
// update_values | update_quadrature_points |
// update_normal_vectors | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
// const unsigned int n_face_q_points = face_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)) *
fe_values.quadrature_point(q_point)[0];
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));
}
for (unsigned int vertex=0;
vertex < GeometryInfo<2>::vertices_per_cell;
++vertex)
for (unsigned int vertex = 0; vertex <
GeometryInfo<1>::vertices_per_cell; ++vertex
if (cell->at_boundary(1) == true) {
fe_values.reinit(cell);
for (unsigned int q_point=0; q_point<n_q_points; ++q_point){
const double neumann_value = -0.5;
for (unsigned int i=0; i<dofs_per_cell; ++i)
cell_rhs(i) += (neumann_value *
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));
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
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);
}
Thank you.
Pietro
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
