Defining evolving non homogenous Dirichlet boundary conditions in a NonlinearForm

[thelfer]

Dear all,

I have defined a NonlinearForm on which I want to impose evolving non homogenous Dirichlet boundary conditions. By evolving, I mean that the NonlinearForm will be used for several time steps and I want to change the values of the imposed Dirichlet boundary conditions at each time step.

It seems that there is no example of using SetEssentialBC with a right hand side.

To make things clearer, here is a (non functional, of course) code snipset which tries to impose the boundary conditions on a unit cube. The faces of the cube are xy0 (normal along Oz), xy1, etc...

I came to the following code which I call at each time step:

    auto fixed_dirichlet_dofs = mfem::Array<int>{};
    auto fixed_dirichlet_values = std::vector<double>{};
    auto append_dof = [&fixed_dirichlet_dofs, &fixed_dirichlet_values, &mesh,
                       &fespace](const auto boundary, const auto normal,
                                 const auto value) {
      auto tmp = mfem::Array<int>{};
      auto ess_bdr = mfem::Array<int>(mesh->bdr_attributes.Max());
      ess_bdr = 0;
      ess_bdr[boundary] = 1;
      fespace->GetEssentialTrueDofs(ess_bdr, tmp, normal);
      fixed_dirichlet_dofs.Append(tmp);
      fixed_dirichlet_values.resize(fixed_dirichlet_values.size() + tmp.Size(),
                                    value);
    };
    append_dof(0, 1, 0.);     // xz0
    append_dof(1, 2, 0.);     // xy0
    append_dof(4, 0, 0.);     // yz0
    append_dof(2, 0, 1e-2 * (t + dt) / 10);  // yz1
    auto rhs = mfem::Vector{fixed_dirichlet_values.data(),
                            fixed_dirichlet_values.size()};
    p.SetEssentialBC(fixed_dirichlet_dofs, &rhs);

where p is my non linear form.

I don't know is this is the right way to do this. Indeed, there are many thing questionable here:

• I selected SetEssentialBC but its implementation shows that I am probably totally wrong.
• I call SetEssentialBC at each time step. I don't know if it is valid.
• I shall probably precomputed fixed boundary conditions (i.e. boundary conditions whose values don't evolve)...

Currently, my residuals are 0, so probably no BC is applied and my code finishes with a memory fault (double free)...

Any help would be appreciated.

[jandrej]

You don't have to call SetEssentialBC after every time step with a NonlinearForm.

This can be made more clear by looking at how the NonlinearForm is used in computing the residual. If you define F(x) = b where F(x) is the nonlinear part and b is a linear term, then SetEssentialBC marks the dofs in x to automatically evaluate to zero if the form is applied to x.

Pseudocode example

NonlinearForm N;

Array<int> all_boundaries_dofs = .. mark all dofs 

N.AddDomainIntegrator(new NonlinearIntegrator);
N.SetEssentialBC(all_boundaries_dofs);

N.Mult(x);

Then, no matter what value is in x at all_boundaries_dofs, the residual is overwritten with zero values at those dofs.

If you want to change the essential dofs for the next timestep, you have to set the values in your solution vector before the computation. The rhs, contrary to a linear computation, does not matter in case of essential dofs.

[thelfer]

Thanks @jandrej
I understand things a bit better. Sorry for this newbie question.
What would be the simpliest way to assign the value of a component of the unknown on a given boundary ?

[thelfer]

I have tried the following:

  auto fixed_dirichlet_dofs = mfem::Array<int>{};
  auto append_dof = [&](const auto boundary, const auto normal) {
    auto tmp = mfem::Array<int>{};
    auto boundaries_markers = mfem::Array<int>(mesh->bdr_attributes.Max());
    boundaries_markers = 0;
    boundaries_markers[boundary] = 1;
    fespace->GetEssentialTrueDofs(boundaries_markers, tmp, normal);
    fixed_dirichlet_dofs.Append(tmp);
    return boundaries_markers;
  };
  append_dof(0, 1);  // xz0
  append_dof(1, 2);  // xy0
  append_dof(4, 0);  // yz0
  yz1_markers = append_dof(2, 0);  // yz1
  p.SetEssentialBC(fixed_dirichlet_dofs);

And at each time step:

   mfem::GridFunction x(fespace.get());
    x.MakeTRef(fespace.get(), u1, 0);
    x.SetFromTrueVector();
    mfem::Vector u_values(3);
    u_values(0) = 1e-2 * (t + dt) / 10;
    u_values(1) = u_values(2) = 0;
    mfem::VectorConstantCoefficient u_coefficients(u_values);
    x.ProjectBdrCoefficient(u_coefficients, yz1_markers);

where u1 is an mfem::Vector which is used to initialize the Newton algorithm.

   newton_solver.Solve(u1);

However:

• It does not work as my residual is still zero
• I assign values to components that are not affected by the boundary condition (u_values(1) and u_values(2) which is probably bad...

[jandrej]

Without a RHS you call

Vector zero;
newton_solver.Mult(zero, x);

where x has to be your solution vector, initialized with the proper essential dofs on your marked boundary points.

Can you try what you did above without making x a reference but rather use u1 directly in the newton solver?

[thelfer]

Hi,
I just changed my code to a much simplier version. I now call setEssentialTrueDofs Still no effect...

  auto fixed_dirichlet_dofs = mfem::Array<mfem_mgis::size_type>{};
  auto append_dof = [&](const auto boundary, const auto normal) {
    auto tmp = mfem::Array<int>{};
    auto boundaries_markers = mfem::Array<int>(mesh->bdr_attributes.Max());
    boundaries_markers = 0;
    boundaries_markers[boundary] = 1;
    fespace->GetEssentialTrueDofs(boundaries_markers, tmp, normal);
    fixed_dirichlet_dofs.Append(tmp);
    return tmp;
  };
  append_dof(0, 1);  // xz0
  append_dof(1, 2);  // xy0
  append_dof(4, 0);  // yz0
  auto yz1_ux_dofs = append_dof(2, 0);  // yz1
  p.SetEssentialTrueDofs(fixed_dirichlet_dofs);

and set the imposed values as follows:

    const auto u = 1e-2 * (t + dt) / 10;
    for (mfem_mgis::size_type idx = 0; idx != yz1_ux_dofs.Size(); ++idx) {
      u1[yz1_ux_dofs[idx]] = u;
    }

[thelfer]

@jandrej With the previous solution, u1 values seems properly set. However, the local displacements passed to the AssembleElementVector and AssembleElementGrad are always zero. Hence my strain are zero, so are my stresses and in fine my residuals...
Any idea ? Could this be a bug in the NonlinearForm ?

[jandrej]

This will not work with a NonlinearForm.

It would help if you can provide a minimal failing example where someone else can also chime in to help with setting a certain component on a GridFunction.

[thelfer]

@jandrej I reopened Issue #1946 with a minimal example. Altough it does not use a GridFunction, it may help, don't you think ?

[thelfer]

Ok. The culprit has been found. I copy/pasted code from ex19 where the iterative_mode was set to false as an option to the Newton solver. When this flag is set to true, it seems to works as expected.... Thanks for your help.