Hi Pietro
Why are you note using a volume integral to evaluate the Neumann
contribution at a face (i.e. a point)? Surely you should be using the
equivalent of the face quadrature in 1D? I know this is covered in the
docs somewhere.
regards
Andrew
On 28 Sep 2009, at 7:21 PM, Pietro Maximoff wrote:
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
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
