Hi Jane Lee,
I recently came across a similar problem.
On Thursday, August 30, 2018 at 12:35:47 PM UTC+2, Jane Lee wrote:
>
> I believe the Neumann conditions are strongly imposed.
>
> And yes - I realised that inhomogeneous Neumann bc is ambiguous phrasing.
>
> I mean that I have a conditions k grad p.n =g, or u.n = g equivalently, I
> think this is in point 3 in my notes in my original post but i think it was
> unclear what u was. I want to impose a nonzero condition on the normal
> derivative of the pressure.
>
> Yes, Neumann condition is a bit ambiguous here. It is often better to
talk about "essential boundary conditions". Essential means that the
boundary conditions need to be built in the function space. If you test
with v and do an integration by parts on the (non-mixed) problem -div
(K*grad p) = f you will see that the normal derivative of p along the
boundary appears in the variational formulation. A condition such as
K*(grad p) *n=g is then called "natural". A condition on p itself needs to
be hard coded into the function space you are solving in and is called
"essential".
Now, that you have a mixed form of the equation you do quite the same: you
test with a pair of test functions (tau and q) in suitable spaces (step-20:
H(div) x L^2) and then you integrate by parts. What you then see is that a
condition on p appears in the variational form (hence here it is natural)
whereas any condition on n*u~K*(grad p) *n needs to be hard coded into the
space H(div) (i.e., it is essential).
Now, if you use like in step-20 Raviart-Thomas elements of lowest order
(RT0) you need to keep in mind that the degrees of freedom are edge moments
(edge averages in case of RT0). It is possible to set them manually by
building them into a constraint matrix like that:
ConstraintMatrix 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);
This is how I did it. Hope that helps.
--
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.
