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: [email protected]
www: http://www.math.colostate.edu/~bangerth/
--
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/8d261b04-0737-1ff3-9b9f-f601af352499%40colostate.edu.
