Dear All,
I'm currently adapting step-21, but have come up with a problem. I want
to write the routine in a way that loops over each cell, then each face,
only visiting a face when the index of the neighbor is higher. This way
each face is only visited once (compare to the two alternatives in step 12).
I know the current value of u and p, so want to construct the right hand
side F_3 - \Delta t H (in the notation of step-21). However, now that I
am only visiting each face once, I need to calculate for both sides of
the edge in one visit.
In the code below I want to calculate the term \int_e jump{dh} \cdot
average{DD(uh)\nabla ch} where jump(a) = a1 n1 + a2 n2 for cells 1 and
2, n1 being the normal from cell 1, and average(b) = 0.5 (b1 + b2).
The code compiles and runs, but does not give a correct (or even
sensible) solution. I don't know where it is going wrong, but suspect it
has something to do with nbr_rhs not being put into the correct position
in system_rhs. I should only be looking at the jump between the ith dof
and the corresponding dof on the neighbor, but I don't think I am doing
this.
I would be thankful if someone could pass an expert eye over this and
tell me where I'm going wrong.
Many thanks,
John
for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
++face_no)
{
fe_face_values.reinit (cell, face_no);
// We are only interested in edges on the interior
// Our no flow boundary conditions are weakly enforced through the
// bilinear forms
if (cell->at_boundary(face_no)) continue;
const typename DoFHandler<dim>::active_cell_iterator nbr =
cell->neighbor(face_no);
const unsigned int nbr_face = cell->neighbor_of_neighbor(face_no);
// We only want to visit each face once, so we arbitarialy chose to sum
// from the lower numbered cell
// This will be more involved for adaptive meshes, so we put an assert
// in now to remined us
Assert (nbr->level() == cell->level(), ExcNotImplemented());
// Only integrate face when nbr is higher
if (nbr->index() > cell->index()) continue;
// Now we know we will certainly use them, we can get the function values
// on this face and the neighbours
nbr_rhs = 0;
fe_face_values.get_function_values
(old_solution,old_solution_values_face);
fe_face_values.get_function_values
(solution,present_solution_values_face);
const std::vector<Point<dim> > &normals =
fe_face_values.get_normal_vectors();
fe_face_values[concentration].get_function_gradients
(old_solution,old_solution_grads_face);
fe_face_values_nbr.reinit (nbr, nbr_face);
fe_face_values_nbr.get_function_values
(old_solution,old_solution_values_face_nbr);
fe_face_values_nbr.get_function_values
(solution,present_solution_values_face_nbr);
fe_face_values_nbr[concentration].get_function_gradients
(old_solution,old_solution_grads_face_nbr);
phi.value_list (fe_face_values.get_quadrature_points(),phi_face_values);
// pick up the concentration at the neighbouring quadrature points
for (unsigned int q=0; q<n_face_q_points; ++q)
{
nbr_face_concentration[q] = old_solution_values_face_nbr[q](dim+1);
}
// Now we are ready to integrate this face
for (unsigned int q=0; q<n_face_q_points; ++q)
{
// We are working with the concentration at the nth time step
const double old_c = old_solution_values_face[q](dim+1);
const double old_c_nbr = old_solution_values_face_nbr[q](dim+1);
const Tensor<1,dim> old_c_grad = old_solution_grads_face[q];
const Tensor<1,dim> old_c_grad_nbr = old_solution_grads_face_nbr[q];
const double JxW = fe_face_values.JxW(q);
Tensor<1,dim> present_u_face;
Tensor<1,dim> present_u_face_nbr;
// Pick up the velocity at n+1.
for (unsigned int d=0; d<dim; ++d)
{
present_u_face[d] = present_solution_values_face[q](d);
present_u_face_nbr[d] = present_solution_values_face_nbr[q](d);
}
// loop over each dof on the face and construct those for this
// dof and the neighbor
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const double dh_i = fe_face_values[concentration].value (i,q);
const double dh_i_nbr
= fe_face_values_nbr[concentration].value(i,q);
const Tensor<1,dim> grad_dh_i
= fe_face_values[concentration].gradient(i,q);
const Tensor<1,dim> grad_dh_i_nbr
= fe_face_values_nbr[concentration].gradient(i,q);
const double Dtphi = time_step/phi_face_values[q];
const Tensor<2,dim> DD
= diffusion_dispersion.tensor_value
(phi_face_values[q],present_u_face);
const Tensor<2,dim> DD_nbr
= diffusion_dispersion.tensor_value
(phi_face_values[q],present_u_face_nbr);
const int theta = 1;
const Point<dim> n = normals[q];
// $-\jump{dh}*\average{DD(uh^{n+1}} \nabla ch^n}
local_rhs(i) += theta*(0.5*Dtphi)*
(dh_i*((DD*old_c_grad+DD_nbr*old_c_grad_nbr)*n))*JxW;
nbr_rhs(i) -= theta*0.5*Dtphi*
(dh_i_nbr*((DD*old_c_grad+DD_nbr*old_c_grad_nbr)*n))*JxW;
}
}
// We have sum the neighbor dofs here before nbr goes out of scope
nbr->get_dof_indices (nbr_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
system_rhs(nbr_dof_indices[i]) += nbr_rhs(i);
}
} // End of loop over faces
// Add local contributions to system, including the cell terms already
calculated
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
system_rhs(local_dof_indices[i]) += local_rhs(i);
}
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
