If cell is a DoFHandler::active_cell_iterator, will cell->at_boundary(0) and
cell->vertex(0) always mean the left vertex of the cell? I didn't find any
information from the documentation.
I'm programing on a 1D complex valued helmholtz equation with mixed boundary
conditions and I can't use quadrature points or FEFaceValues to evaluate the
boundary terms. Instead, I use cell->at_boundary(0) to identify the left
boundary and FiniteElement->shape_value(i, cell->vertex(0)) to compute the
shape value at the bounary. The program runs but the solution is not correct
compared with the exact solution. I don't know whether my mathematical
formulation or the code is wrong.
This is the major part of my assemble_system function:
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
if (fe.system_to_component_index(i).first ==
fe.system_to_component_index(j).first)
for (unsigned int q=0; q<n_q_points; ++q)
cell_matrix(i,j) += ( fe_values.shape_value(i,q)
* fe_values.shape_value(j,q)
* coefficient_values[q]
-
fe_values.shape_grad(i,q)
* fe_values.shape_grad(j,q)
/ mu_values[q]
) * fe_values.JxW(q);
if ( cell -> at_boundary(0) )
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
if ( (fe.system_to_component_index(i).first !=
fe.system_to_component_index(j).first) )
cell_matrix(i,j) += ( (
fe.system_to_component_index(i).first ) ? -1 : 1 )
* fe.shape_value(i, cell ->
vertex(0) )
* fe.shape_value(j, cell ->
vertex(0) )
* alpha_inc;
if ( fe.system_to_component_index(i).first )
cell_rhs(i) += fe.shape_value( i, cell -> vertex(0) ) * 2.0 *
alpha_inc;
}
if ( cell -> at_boundary(1) )
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
if ( (fe.system_to_component_index(i).first !=
fe.system_to_component_index(j).first) )
cell_matrix(i,j) += ( (
fe.system_to_component_index(i).first ) ? -1 : 1 )
* fe.shape_value(i, cell ->
vertex(1) )
* fe.shape_value(j, cell ->
vertex(1) )
* alpha_sub;
--
*Yuliang*
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
