Hi there Pietro and Andrew,
Two points about the code given below: the loop over vertices seems to
be unnecessary, and due to c style numbering fe.dofs_per_cell will
give you a value too large by one for cell_rhs. My two cents:
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));
}
}
if (cell->at_boundary(1) == true)
{
const double neumann_value = -0.5;
cell_rhs(1) += (neumann_value); // see comment below on
why 1 is used
}
} // end of loop on cells
// Copy local matrix into global matrix
In the code above I use the original test for a cell being at a
boundary, in the loop over cells which will iterate over all cells. If
we are at the end cell the if statement executes. There are further
complexities in choosing the correct dof to assign the neumann_value
to because the ordering of dofs in the FE is defined by the library
internals only and should not be relied on. For the FE_Q class the
documentation specifies that LHS dof gets 0 and RHS dof gets 1, the
remaining dofs are assigned to any internal points. Hence 1 seems to
be appropriate for FE_Q in this case but writing code that depends on
this ordering is not advised. I cannot think of the appropriate method
to use to avoid it in this case however.
Cheers,
Michael
On Tue, Sep 29, 2009 at 9:51 AM, Andrew McBride
<[email protected]> wrote:
> Hi,
> No, you can't use a face quadrature rule in one-dimension as the edge is
> simply a point and thus there is nothing to do. Thus, deal will throw an
> exception.
> This should work:
> for (unsigned int vertex = 0; vertex <
> GeometryInfo<1>::vertices_per_cell; ++vertex)
> if (cell->at_boundary(1) == true) {
> const double neumann_value = -0.5;
> cell_rhs(fe.dofs_per_cell) += (neumann_value);
> }
>
>
> Andrew
> On 29 Sep 2009, at 9:26 AM, Pietro Maximoff wrote:
>
>
> ________________________________
> From: Andrew McBride <[email protected]>
> To: Pietro Maximoff <[email protected]>
> Cc: [email protected]
> Sent: Tuesday, September 29, 2009 7:39:28 AM
> Subject: Re: [deal.II] 1D Boundary-Value Problem
>
> 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
>
> ============
>
> Do you mean using a face quadrature formula like so:
>
> QGauss<dim-1> face_quadrature_formula(3); ?
>
> Well, QGauss<0> generates an exception. So, doing it this way has been
> impossible so far.
>
> Incidentally, the exact solution is: 0.5ln(x) + 2/x.
>
> Please help.
>
> Thanks.
>
> Pietro
>
>
>
>
> _______________________________________________
> dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
>
>
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
