Hi Konrad,
I may add as a comment (without going into details): Your problem is - as
you correctly describe above - a mixed Poisson problem. You essentially
have two options to solve it and the way you choose influences the way you
need to treat the boundary conditions: in primal form and in mixed form. In
your formulation you have natural boundary conditions (Dirichlet) on Phi
which means that you need to prescribe them in the variational form. This
will lead to a boundary integral on the relevant part of the boundary and
result in an additional term on the right-hand side. The boundary condition
for E is a flux condition and in your mixed form this (Neumann) condition
is essential. This means that you need to enforce them on the
Raviart-Thomas space as a hard constraints. There are essentially two
options for this:
1. If you have given a vector field that describes E on the relevant part
of the boundary you can use the function
*VectorTools::project_boundary_values_div_conforming. *This will project
the given vector field onto its normal component (since this is the one you
must describe with RT elements in H(div)) and use the projected value as
boundary condition for the normal flux of E
2. If you only are given the scalar value of the normal flux of E then you
may have to set the relevant dofs in the RT-space by hand. This you can do
through constraints that way:
AffineConstraints<double> constraints;
constraints.clear ();
DoFTools::make_hanging_node_constraints (dof_handler, constraints);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int dofs_per_face = fe.dofs_per_face;
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_face_indices
(dofs_per_face);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
cell->get_dof_indices(local_dof_indices);
for (unsigned int face_n=0;
face_n<GeometryInfo<dim>::faces_per_cell;
++face_n)
{
cell->face(face_n)->get_dof_indices(local_dof_face_indices);
if (cell->at_boundary(face_n))
{
if (cell->face(face_n)->boundary_id()==id_bottom_boundary)
{
for (unsigned int i = 0; i<dofs_per_face; ++i)
{
// see documentation for these functions
constraints.add_line (local_dof_face_indices.at(i));
constraints.set_inhomogeneity (
local_dof_face_indices.at(i),
value_of_dof);
}
}
else // add some constraints on other boundaries. Here:
zero-flux
{
for (unsigned int i = 0; i<dofs_per_face; ++i)
{
constraints.add_line (local_dof_face_indices.at(i));
}
}
}
}
++cell_n;
}
constraints.close();
Once you use the DoFTools::make_sparsity_pattern function you need to make
sure to keep the inhomogenities:
DoFTools::make_sparsity_pattern(dof_handler,
dsp,
constraints,
/*keep_constrained_dofs = */
true);
The same in the assembly routine:
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(local_matrix,
local_rhs_local,
local_dof_indices,
system_matrix,
system_rhs,
/* use
inhomogeneities for rhs */ true);
Keep in mind that the dofs for RT elements are edge moments in 2D and
face moments in 3D (and not nodal values).
Best,
Konrad
--
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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/48e77a20-69a8-4108-a5bc-d61adf3fba95%40googlegroups.com.
