Dr. Bangerth,
Thank you very much for your kind reply.
> What is the content of 'constraints' in your case?
>
> DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
VectorTools::interpolate_boundary_values (dof_handler,
0,
ZeroFunction<dim>(),
constraints);
constraints.close ();
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler,
dsp,
constraints,
/*keep_constrained_dofs = */ false);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit (sparsity_pattern);
}
As you said in the lecture, I apply the ZeroFunction<dim>() for
A*U_0=F-A*G_tilt
>
> This is not correct. You are computing the G vector as
> G_i = \int_{\partial\Omega} \varphi_i(x) g(x) ds
> i.e., as a boundary integral. But G should be the vector that
> *interpolates*
> g(x) on the boundary.
>
> template <int dim>
void Step6<dim>::assemble_G ()
{
const QGauss<dim-1> face_quadrature_formula(3);
FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
update_values | update_gradients |
update_quadrature_points |
update_JxW_values);
const unsigned int dofs_per_face = fe.dofs_per_face;
const unsigned int n_face_q_points = face_quadrature_formula.size();
Vector<double> local_rhs (dofs_per_face);
G_tilt.reinit (dof_handler.n_dofs());
std::vector<types::global_dof_index> local_dof_indices (dofs_per_face);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
local_rhs = 0;
for (unsigned int face_no=0;
face_no<GeometryInfo<dim>::faces_per_cell;
++face_no)
if (cell->at_boundary(face_no))
{
fe_face_values.reinit (cell, face_no);
for (unsigned int q=0; q<n_face_q_points; ++q)
{
const Point<dim> &p=fe_face_values.quadrature_point(q);
for (unsigned int i=0; i<dofs_per_face; ++i)
local_rhs(i) += (p(0)*p(0)+p(1)*p(1));
}
cell->face(face_no)->get_dof_indices (local_dof_indices);
}
for (unsigned int i=0; i<dofs_per_face; ++i)
G_tilt(local_dof_indices[i]) += local_rhs(i);
}
}
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?
Thank you.
Kyusik.
--
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
