> 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;
This looks like it could work, but probably only by coincidence. You are
evaluating the finite element shape function at point cell->vertex(0), i.e.
at a coordinate point in real space; on the other hand, the
FiniteElement::shape_value function wants as its argument a location on the
reference cell. In other words, this *may* work if you solve on the domain
[0,1] but it will definitely not otherwise.
The FEValues class is designed to eliminate exactly this mismatch. You should
use FEValues with a quadrature formula that has only two points, one at the
left end and one at the right (i.e. the QTrapez class) and then evaluate
these terms with FEValues just like you do for all the other parts of your
bilinear form.
Best
W.
-------------------------------------------------------------------------
Wolfgang Bangerth email: [email protected]
www: http://www.math.tamu.edu/~bangerth/
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
