Re: [deal.II] Odd behaviour in VectorTools::interpolate_boundary_values

[Wolfgang Bangerth]

Stephen,
apologies for taking 2 weeks to get to your question here:

> I have two different codes to solve the same problem but one works with
> VectorTools::interpolate_boundary_values whereas the other doesn't and I'm
> trying to figure out if I've done something wrong or there's a bug in the
> library. In both codes, I have an FESystem made up of two smaller FESystems:
>
> fe_trial_interior (FE_DGQ<dim>(degree), 1, FE_DGQ<dim>(degree), dim),
> fe_trial_trace (FE_TraceQ<dim>(degree + 1), 1, FE_FaceQ<dim>(degree), 1),
> fe_trial (fe_trial_interior, 1, fe_trial_trace, 1)
>
> In one code, I apply the non-zero dirichlet boundary conditions to the trace
> variables directly via
>
> const FEValuesExtractors::Scalar dirichlet_index (0); const ComponentMask
> dirichlet_comp_mask = fe_trial_trace.component_mask (dirichlet_index);
> const FEValuesExtractors::Scalar robin_index (1); const ComponentMask
> robin_comp_mask = fe_trial_trace.component_mask (robin_index);
>
> trace_constraints.clear ();
> VectorTools::interpolate_boundary_values (dof_handler_trial_trace, 0,
> DirichletBoundaryValues<dim>(), trace_constraints, dirichlet_comp_mask); //
> Dirichlet boundary constraints.
> VectorTools::interpolate_boundary_values (dof_handler_trial_trace, 1,
> RobinBoundaryValues<dim>(), trace_constraints, robin_comp_mask); // Robin
> boundary constraints.
> DoFTools::make_hanging_node_constraints (dof_handler_trial_trace,
> trace_constraints); // Hanging node constraints.
> trace_constraints.close ();
>
> which works perfectly (error confirmed to go to zero for this code). In the
> second code, I'm trying to apply the same non-zero boundary conditions but to
> the full system via
>
> const FEValuesExtractors::Scalar dirichlet_index (dim + 1); const
> ComponentMask dirichlet_comp_mask = fe_trial.component_mask (dirichlet_index);
> const FEValuesExtractors::Scalar robin_index (dim + 2); const ComponentMask
> robin_comp_mask = fe_trial.component_mask (robin_index);
>
> constraints.clear ();
> VectorTools::interpolate_boundary_values (dof_handler_trial, 0,
> DirichletBoundaryValues<dim>(), constraints, dirichlet_comp_mask); //
> Dirichlet boundary constraints.
> VectorTools::interpolate_boundary_values (dof_handler_trial, 1,
> RobinBoundaryValues<dim>(), constraints, robin_comp_mask); // Robin boundary
> constraints.
> DoFTools::make_hanging_node_constraints (dof_handler_trial, constraints); //
> Hanging node constraints.
> constraints.close ();
>
> and this is giving me issues. In particular,
> VectorTools::interpolate_boundary_values does seem to be constraining the
> correct nodes but it is not correctly setting the inhomogeneity (thus it is
> setting them to zero). Would welcome any thoughts as to if I'm doing anything
> wrong here or if it's a bug.

Interesting problem. In general, these sorts of things are most easily
debugged by coming up with a simple, minimal working example (MWE). Here, this
would involve a mesh with just a single cell, where you create the constraints
as you do above and output what you get. That should make it relatively
straightforward to reason about the correctness of what you get. Having such a
MWE would be excellent, if you could try to get one given that you already
have a good idea of what is going wrong!

In this particular case, I would not be greatly surprised if the root of the
problem is due to the fact that you have a nested FESystem going on:

fe_trial_interior (FE_DGQ<dim>(degree), 1, FE_DGQ<dim>(degree), dim),
fe_trial_trace (FE_TraceQ<dim>(degree + 1), 1, FE_FaceQ<dim>(degree), 1),
fe_trial (fe_trial_interior, 1, fe_trial_trace, 1)

It would be interesting to see whether the problem disappears if you unnest
things:

fe_trial (FE_DGQ<dim>(degree), 1, FE_DGQ<dim>(degree), dim,
FE_TraceQ<dim>(degree + 1), 1, FE_FaceQ<dim>(degree), 1)

Of course, even if this now produces correct results, it would be interesting
to trace the origin of the bug and see whether we can't fix it anyway!

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Stephen]

There was actually no bug here and I managed to fix the problem but thanks!