On 03/19/2017 09:21 PM, [email protected] wrote:
As you can see above, what I change here is....
1. const unsigned int dofs_per_cell = fe.dofs_per_cell; -> const
unsigned int dofs_per_face = fe.dofs_per_face;
2. Vector<double> local_rhs(dofs_per_cell); -> Vector<double>
local_rhs (dofs_per_face);
3.std::vector<types::global_dof_index> local_dof_indices
(dofs_per_cell); -> std::vector<types::global_dof_index>
local_dof_indices (dofs_per_face);
4.
for (unsigned int i=0; i<dofs_per_cell; ++i)
local_rhs(i) += fe_face_values.shape_value(i,q)*
(p(0)*p(0)+p(1)*p(1));
->
for (unsigned int i=0; i<dofs_per_face; ++i)
local_rhs(i) += (p(0)*p(0)+p(1)*p(1));
5.
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
->
cell->face(face_no)->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_face; ++i)
But the result is still strange...
I don't know how I can get the indices of G on the boundary.
Could you please teach me how to do that?
It is not correct because (i) you really do need to loop over all
dofs_per_cell, and (ii) think of a vertex that is shared by two faces.
The way you do things here you get a contribution from each face, but
because you add things up, you get twice the value you want in the
result. If you output what you compute on every cell and what ends up in
the G vector, you can see this -- in fact, you should have tried that
yourself already, printf and debuggers are your friends, as are problems
on small meshes for which you can write things down on a piece of paper
and compare with what your program computes.
The approach you use, making things look like an assembly step, is not
appropriate for what you want to do. I would just call
VectorTools::interpolate_boundary_values to get the pairs of dof_index
and value, and then convert the result of that function into a vector.
Best
W.
--
------------------------------------------------------------------------
Wolfgang Bangerth email: [email protected]
www: http://www.math.colostate.edu/~bangerth/
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
For more options, visit https://groups.google.com/d/optout.
