Kyusik,
I'm studying Inhomogenous Dirichlet boundary conditions watching Lecture 21.65.
To practice this boundary condition, I'm trying to reproduce the result of
step-4(in which Inhomogenous Dirichlet boundary condition is used) in step-6.
As I learned in the lecture...
First, I assembled A and F with zero boundary condition by....
template <int dim>
void Step6<dim>::assemble_system ()
{
const QGauss<dim> quadrature_formula(3);
const RightHandSide<dim> right_hand_side;
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit (cell);
for (unsigned int q_index=0; q_index<n_q_points; ++q_index)
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_index) *
fe_values.shape_grad(j,q_index) *
fe_values.JxW(q_index));
cell_rhs(i) += (fe_values.shape_value(i,q_index) *
right_hand_side.value(fe_values.quadrature_point(q_index)) *
fe_values.JxW(q_index));
}
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (cell_matrix,
cell_rhs,
local_dof_indices,
system_matrix,
system_rhs);
}
}
What is the content of 'constraints' in your case?
And then, I made G_tilt by...
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_cell = fe.dofs_per_cell;
const unsigned int n_face_q_points = face_quadrature_formula.size();
Vector<double> local_rhs (dofs_per_cell);
G_tilt.reinit (dof_handler.n_dofs());
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
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_cell; ++i)
local_rhs(i) += fe_face_values.shape_value(i,q)*
(p(0)*p(0)+p(1)*p(1));
}
}
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
G_tilt(local_dof_indices[i]) += local_rhs(i);
}
}
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.
And I solved AU_0=F-AG_tilt by...
Vector<double> tmp;
tmp.reinit(solution.size());
system_matrix.vmult(tmp,G_tilt);
system_rhs -= tmp;
solve ();
solution += G_tilt;
This looks correct.
Best
W.
--
------------------------------------------------------------------------
Wolfgang Bangerth email: bange...@colostate.edu
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
