constraints.distribute does not seem to be enforcing Dirichlet conditions

[Lucas Myers]

Hi everyone,

I'm trying to solve a nonlinear ODE with Dirichlet conditions using Newton's method. For this I give an initial guess for the solution ( y = (1/2) x ) and then update each Newton step with an update that ought to be zero at the boundaries (hence preserving the boundary conditions of the initial guess). However, for some reason my update is taking on nonzero values at the boundaries.

To enforce these boundary conditions, I initialize constraints with:

```

constraints.clear();
dealii::DoFTools::make_hanging_node_constraints(dof_handler, constraints);
 
dealii::VectorTools::
    interpolate_boundary_values(dof_handler,
                                                         0,
                                                         dealii::Functions::ZeroFunction<dim>(1),
                                                         constraints);
constraints.close();

```

I use the `distribute_local_to_global` in the assembly function, and then once I've solved the linear system (directly, with UMFPACK), I use `constraints.distribute`. Not sure what else I need to do to enforce those boundary conditions.

I'm attaching the code (<300 lines including white space) in case that is helpful. Any help is appreciated!

- Lucas

[Wolfgang Bangerth]

On 8/30/22 13:56, Lucas Myers wrote:
>
> To enforce these boundary conditions, I initialize constraints with:
> ```
> constraints.clear();
> dealii::DoFTools::make_hanging_node_constraints(dof_handler, constraints);
>
> dealii::VectorTools::
>     interpolate_boundary_values(dof_handler,
>                                                          0,
>
> dealii::Functions::ZeroFunction<dim>(1),
>                                                          constraints);
> constraints.close();
> ```
> I use the `distribute_local_to_global` in the assembly function, and then once
> I've solved the linear system (directly, with UMFPACK), I use
> `constraints.distribute`. Not sure what else I need to do to enforce those
> boundary conditions.
>
> I'm attaching the code (<300 lines including white space) in case that is
> helpful. Any help is appreciated!

Rather than trying to debug this myself, let me ask some questions:
* Right after the call you show above, are the necessary constraints in your
'constraint' object?
* Right before you call 'distribute_local_to_blobal()' are the necessary
constraints still in the 'constraints' object?
* Right before you call 'distribute()' after you solve the linear system, are
the necessary constraints still in the 'constraints' object?
* Right after you call 'distribute()', are the entries in the solution vector
that correspond to boundary nodes correct?

If the answer to any of the questions above is 'no', then you've got a spot to
further debug. (All of this might be easiest to check on a very coarse mesh,
say a 2x2 mesh.) Right now, all you see is that the *end result* is wrong.
You've got to narrow things down to a smaller part of the code, and the
questions above should help you with it!

Best
W.

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

[Lucas Myers]

Thanks for the help! I printed out the number of constraints after each of your suggestions, and after setting the hanging node constraints the number of constraints is 0 (as it should be, I'm working with a 1D mesh that has been globally refined some number of times), but after interpolating the boundary values the number of constraints is only 1. Given that dim = 1 and I'm trying to enforce Dirichlet boundaries, I reckon there should be 2 constraints. Running the std::map version of interpolate_boundary_values, I still only get 1 constraint.

So my question is: is this a bug? I'm trying to step through it in GDB, but it hangs when I try to step into this function, so I'm a little bit at a loss as to how to continue debugging.

- Lucas

[Wolfgang Bangerth]

On 8/31/22 11:09, Lucas Myers wrote:
> Thanks for the help! I printed out the number of constraints after each
> of your suggestions, and after setting the hanging node constraints the
> number of constraints is 0 (as it should be, I'm working with a 1D mesh
> that has been globally refined some number of times), but after
> interpolating the boundary values the number of constraints is only 1.
> Given that dim = 1 and I'm trying to enforce Dirichlet boundaries, I
> reckon there should be 2 constraints. Running the std::map version of
> interpolate_boundary_values, I still only get 1 constraint.

Ah, that's a commonly asked question. In 2d and 3d, the default is for
the entire boundary to have boundary ID equal to zero. In 1d, it is zero
for the left end point and one for the right end point. So you have to
call interpolate_b_v() twice, once for each of the two boundary IDs.