Euler's Elastica nonconvergence

variational formulation
1

I’m trying to solve Euler’s elastica equation (see this link) for a loop in the plane:

u: \mathbb{R} / \mathbb{Z} \to \mathbb{R}^2

with the bending energy:

E[u] = \int_0^1 (|u''|^2 + p(|u'|^2 - 1))dx

and p a Lagrange multiplier enforcing the inextensibility constraint |u'| = 1. I have expressed it in the weak form:

\int_0^1[u'' \cdot v'' + p u' \cdot v' - q(u'\cdot u' - 1)]dx = 0

Finally, here is my code implementing this:

import numpy as np
import pdb
from mpi4py import MPI
import dolfinx
import dolfinx_mpc
import ufl

from fenics_tools import *

# Create mesh & solution space
n_elem = 50
comm = MPI.COMM_WORLD
domain = dolfinx.mesh.create_unit_interval(comm=comm, nx=n_elem)
P1 = ufl.VectorElement('Lagrange', cell=domain.ufl_cell(), degree=2, dim=2)
P2 = ufl.FiniteElement('Lagrange', cell=domain.ufl_cell(), degree=1)
V = dolfinx.fem.FunctionSpace(domain, ufl.MixedElement([P1, P2]))

# Set periodic boundary conditions
mpc = apply_periodic_bcs(domain, V, 1)

# Create problem
state = dolfinx.fem.Function(V)
twopi = 2 * np.pi
state.sub(0).interpolate(lambda x: np.stack((np.sin(twopi * x[0]), np.cos(twopi * x[0]))) / twopi)
(u, p) = ufl.split(state)
(v, q) = ufl.TestFunctions(V)
dx = ufl.Measure('dx', domain=domain, metadata={'quadrature_degree': 4})
F = (
    ufl.inner(ufl.grad(ufl.grad(u)), ufl.grad(ufl.grad(v))) 
    + p * ufl.inner(ufl.grad(u), ufl.grad(v)) 
    - q * (ufl.inner(ufl.grad(u), ufl.grad(u)) - 1) 
) * dx
problem = dolfinx.fem.petsc.NonlinearProblem(F, state, bcs=[])

# Solve problem
solver = dolfinx.nls.petsc.NewtonSolver(comm, problem)
solver.atol = 1e-6
solver.rtol = 1e-6
solver.convergence_criterion = 'incremental'
solver.report = True
solver.solve(state)
u_sol, p_sol = state.split()
ux_sol, uy_sol = u_sol.split()
ux_sol, uy_sol, p_sol = ux_sol.x.array, uy_sol.x.array, p_sol.x.array

where apply_periodic_bcs is more or less exactly from the MPC example here.

However, it always returns:

RuntimeError: Newton solver did not converge because maximum number of iterations reached

I then tried interpolating an analytical solution, which should be a circle of circumference 1:

u(x) = (\sin(2\pi x)/2\pi, \cos(2\pi x)/2\pi)

But still obtain the same result. I’m new to the nonlinear solver and am wondering if I’m doing something incorrect in setting up this problem. Especially with the analytical solution interpolated, very few to no steps should be required in principle. Could I have some help with this?

2

You should not mix the assembly methods of dolfinx_mpc and dolfinx, as highlighted in all demos of dolfinx_mpc.

If you want to solve a non-linear problem, see
@nate’s excellent demo/test

3

I see, I’ve replaced the calls to NonlinearProblem and NewtonSolver above with the classes from the tests. However, upon calling solver.solve(state) I now see:

...
self.mpc.homogenize(self.u_mpc)
  File "/Users/asrvsn/mambaforge/envs/fenics/lib/python3.10/site-packages/dolfinx_mpc/multipointconstraint.py", line 387, in homogenize
    with vector.localForm() as vector_local:
AttributeError: 'Function' object has no attribute 'localForm'

where u_mpc is from the suggested NewtonSolverMPC class:

class NewtonSolverMPC(dolfinx.cpp.nls.petsc.NewtonSolver):
    def __init__(self, comm: MPI.Intracomm, problem: NonlinearMPCProblem,
                 mpc: dolfinx_mpc.MultiPointConstraint):
        """A Newton solver for non-linear MPC problems."""
        super().__init__(comm)
        self.mpc = mpc
        self.u_mpc = dolfinx.fem.Function(mpc.function_space)
...

It looks like the problem is that in my installed dolfinx_mpc (which is version 0.6.1), the homogenize definition:

    def homogenize(self, vector: _PETSc.Vec) -> None:
        """
        For a vector, homogenize (set to zero) the vector components at the multi-point
        constraint slave DoF indices. This is particularly useful for nonlinear problems.

        Args:
            vector: The input vector
        """
        with vector.localForm() as vector_local:
            self._cpp_object.homogenize(vector_local.array_w)
        vector.ghostUpdate(addv=_PETSc.InsertMode.INSERT, mode=_PETSc.ScatterMode.FORWARD)

is different from the one in the docs.

This looks like a simple installation issue, and I should be using code from this commit.

Should I simply change the homogenize() call to take u.vector instead? Or would it be more advisable to try and build the master branch locally? (Not sure if there would be any API conflicts with other fenics libraries, which were installed via conda.)

4

(I’ve just tried replacing u_mpc with u_mpc.vector in the calls to homogenize() and backsubstitution(), and still see the convergence issue. Is there anything else in my setup that could be causing this?)

5

To be able to give you further guidance, it would be helpful if you post your current attempt of solving the problem.

6

Here is my current code:

import numpy as np
import pdb
from mpi4py import MPI
import dolfinx
import dolfinx_mpc
import ufl

from fenics_tools import *

# Create mesh & solution space
n_elem = 50
comm = MPI.COMM_WORLD
domain = dolfinx.mesh.create_unit_interval(comm=comm, nx=n_elem)
P1 = ufl.VectorElement('Lagrange', cell=domain.ufl_cell(), degree=2, dim=2)
P2 = ufl.FiniteElement('Lagrange', cell=domain.ufl_cell(), degree=1)
V = dolfinx.fem.FunctionSpace(domain, ufl.MixedElement([P1, P2]))

# Set periodic boundary conditions
mpc = apply_periodic_bcs(domain, V, 1)

# Create problem
state = dolfinx.fem.Function(V)
twopi = 2 * np.pi
state.sub(0).interpolate(lambda x: np.stack((np.sin(twopi * x[0]), np.cos(twopi * x[0]))) / twopi)
(u, p) = ufl.split(state)
(v, q) = ufl.TestFunctions(V)
dx = ufl.Measure('dx', domain=domain, metadata={'quadrature_degree': 4})
F = (
    ufl.inner(ufl.grad(ufl.grad(u)), ufl.grad(ufl.grad(v))) 
    + p * ufl.inner(ufl.grad(u), ufl.grad(v)) 
    - q * (ufl.inner(ufl.grad(u), ufl.grad(u)) - 1) 
) * dx
problem = NonlinearMPCProblem(F, state, mpc, bcs=[])

# Solve problem
solver = NewtonSolverMPC(comm, problem, mpc)
solver.atol = 1e-6
solver.rtol = 1e-6
solver.convergence_criterion = 'incremental'
solver.report = True
solver.solve(state)
u_sol, p_sol = state.split()
ux_sol, uy_sol = u_sol.split()
ux_sol, uy_sol, p_sol = ux_sol.x.array, uy_sol.x.array, p_sol.x.array

The helper functions in fenics_tools can be seen here. I’m using dolfinx-0.6.0 and dolfinx-mpc-0.6.1.

7

Update, I downgraded to dolfin (2019) and think I have a working solution:

import matplotlib.pyplot as plt
import numpy as np
import pdb
import dolfin as df
from dolfin import (
    inner, grad, dx, DOLFIN_EPS
)

# Preconditions
if not df.has_linear_algebra_backend("PETSc"):
    df.info("DOLFIN has not been configured with PETSc. Exiting.")
    exit()

df.parameters["form_compiler"]["cpp_optimize"] = True

# Parameters
n_elem = 100

# Sub domain for Periodic boundary condition
class PeriodicBoundary(df.SubDomain):
    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        return bool( -DOLFIN_EPS < x[0] < DOLFIN_EPS and on_boundary)

    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
        y[0] = x[0] - 1.0

pbc = PeriodicBoundary()

# Create mesh and finite element
degree = 2
mesh = df.UnitIntervalMesh(n_elem)
W = df.VectorElement('CG', mesh.ufl_cell(), degree, dim=2)
Q = df.FiniteElement('CG', mesh.ufl_cell(), degree)
element = df.MixedElement([W, Q])
V = df.FunctionSpace(mesh, element, constrained_domain=pbc)

# Define variational problem
state = df.Function(V)

## Set initial guess (required to prevent singular residual at 0)
state.interpolate(df.Expression((
    "sin(2*pi*x[0])/(2*pi)", 
    "cos(2*pi*x[0])/(2*pi)",
    "1"
), degree=degree))

(u, p) = df.split(state) 
(v, q) = df.TestFunctions(V)
F = (
    inner(grad(grad(u)), grad(grad(v))) * dx
    + p * inner(grad(u), grad(v)) * dx
    + q * (inner(grad(u), grad(u)) - 1) * dx
)

df.solve(F == 0, state, solver_parameters={
    "newton_solver": {
        "linear_solver": "umfpack",
        "relative_tolerance": 1e-6,
        "absolute_tolerance": 1e-8,
    }
})
8

(For anyone revisiting this, I had to convert to a variational problem with first-order derivatives first, then use the umfpack solver with lowered Newton relaxation parameter).

9

Hey @asrvsn. I am trying to solve a similar problem and am unable to find a stable formulation even when I add an extra field to allow me avoid taking fourth order derivatives. Are you able to give me a bit more info about how you solved it?

10

Maybe this model can help. It is extensible and shearable.
https://bleyerj.github.io/comet-fenicsx/tours/beams/finite_rotation_nonlinear_beam/finite_rotation_nonlinear_beam.html
You may need to pay attention to locking issues if you want it inextensible and non-shearable.

11

Thanks I’ll have a look at this!

12

Hey @bleyerj , was wondering you would be able say about more about what you mean by locking issues. I do indeed want it inextnesible. I’ve been able to solve the problem by introducing an auxilliary function to represent the tangent, then everything converges nicely, however this work because when we are in euclidian space, everything can be written in terms of tangents and curvatures, without worrying about the actual configuration r (or u as called in the code above). In the end I am interested in embedding my elastica in a curved surface, which means I cannot work purely with tangents. So I am instead trying to redefine the auxilliary field in terms of the second derivative k (as is done in the Cahn-hilliard example in the fenics tutorial). But this seems to lead to a frightfully ill conditioned system, so the iterations only converge if your initial guess is essentially perfect. Happy to post code if it helps. I’m not sure if there is some big picture point I am missing here