Dear All,
I don't know if it is common to answer ones own posts on the newsgroup,
but I don't want anyone wasting any time trying to solve my problems
when I have worked them out myself.
The solution is quite obvious - two separate loops are required on each
face for the dofs on the current cell and those on the neighbors. There
is no interaction as the other terms are entirely known at the
quadrature points. The loops should look something like this
loop faces
loop quadrature points
determine known values at quadrature points
loop dofs on this cell
assemble for this face/cell
end loop dofs on this cell
loop dofs on neighbor
assemble for this face/neighbor
end loop dofs on neighbor
end
end
All the best,
John
On 22/03/11 14:39, John Chapman wrote:
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
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
