QGauss<3> quadrature_formula(fe.degree + 1);
QGauss<3 - 1> face_quadrature_formula(fe.degree + 1);
FEValues<3> fe_values(fe,
quadrature_formula,
update_values | update_gradients | update_JxW_values);
FEFaceValues<3> fe_face_values(fe,
face_quadrature_formula,
update_values | update_quadrature_points |
update_normal_vectors |
update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
const unsigned int n_face_q_points = face_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);
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
cell_matrix = 0;
cell_rhs = 0;
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) * // grad phi_i(x_q)
fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
fe_values.JxW(q_index)); // dx
for (unsigned int i = 0; i < dofs_per_cell; ++i)
cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
1 * // f(x_q)
fe_values.JxW(q_index)); // dx
}
cell->get_dof_indices(local_dof_indices);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
system_matrix.add(local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i, j));
for (unsigned int i = 0; i < dofs_per_cell; ++i)
system_rhs(local_dof_indices[i]) += cell_rhs(i);
// I know that the hull of the cylinder has boundary id 0, the top face has boundary id 1, and the bottom of the face has boundary id 2.
for(unsigned int face_number = 0; face_number < GeometryInfo<3>::faces_per_cell; ++face_number)
{
if(cell->face(face_number)->at_boundary() && (cell->face(face_number)->boundary_id() == 0))
{
fe_face_values.reinit(cell, face_number);
// If we come in here we have found a face that belongs to the boundary condtion of the hull
// I know that I am supposed to do something like the code in green below, but I don't know the exact solution.
// What I would like to do is set the derivative of my function to zero. My thinking is that it would entail
// taking the gradient of fe_face_values to ZeroFunction<3>(). I think that if I could understand how to apply
// the boundary condition to the hull of the cylinder, I could understand how to apply the the boundary condition
// to the top of the face just as easily. Here is the code:
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
{
const double neumann_value =
(exact_solution.gradient(
fe_face_values.quadrature_point(q_point)) *
fe_face_values.normal_vector(q_point));
for (unsigned int i = 0; i < dofs_per_cell; ++i)
cell_rhs(i) +=
(neumann_value * // g(x_q)
fe_face_values.shape_value(i, q_point) * // phi_i(x_q)
fe_face_values.JxW(q_point)); // dx
}