There was actually no bug here and I managed to fix the problem but thanks! On Wednesday, 2 June 2021 at 18:10:46 UTC-5 Wolfgang Bangerth wrote:
> > 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/ > > -- 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. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/17453237-ea2c-4b16-9eec-7ef487361c9dn%40googlegroups.com.
