Catenary Problem

2

Keep in mind that the Lagrange multiplier in the arc length constraint is a single scalar to enforce one integral constraint over the whole domain, so it should be in a FunctionSpace of type "R", not a DG0 space where it can have a different value on each element. The following code snippet seems to solve the problem correctly:

from dolfin import *
N = 100

# Length of interval
Lhat = 2 

# Length of curve; needs to be larger than Lhat to make sense
L = Constant(10)

mesh = IntervalMesh(N,-0.5*Lhat,0.5*Lhat)
cell = mesh.ufl_cell()
Ye = FiniteElement("CG",cell,1)
# The "R" type space refers to a space with 1 DoF, for a global
# constant value.  
Le = FiniteElement("R",cell,0)
V = FunctionSpace(mesh,MixedElement([Ye,Le]))

# Set up the Lagrangian for the problem:
yl = Function(V)
y,lam = split(yl)
ds_dx = sqrt(1+(y.dx(0))**2)
constraint = (ds_dx-L/Constant(Lhat))
Lagrangian = (y*ds_dx - lam*constraint)*dx
dLagrangian = derivative(Lagrangian,yl)
bc = DirichletBC(V.sub(0),Constant(0),"on_boundary")

# Solver needs some kind of guess to set the sign of the solution.
# If you change the sign of the guess, it converges to a cantenary facing
# the other direction.  May need to adjust amplitude for different values
# of L.
initGuessAmp = 1.0
yl.interpolate(Expression(("amp*(x[0]+0.5*Lhat)*(x[0]-0.5*Lhat)","0"),
                          degree=2,Lhat=Lhat,amp=initGuessAmp))

# Need a more robust nonlinear solver than Newton iteration; try
# PETSc SNES with back-tracking ("bt") line search:
problem = NonlinearVariationalProblem(dLagrangian,yl,bc,
                                      derivative(dLagrangian,yl))
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['nonlinear_solver'] = 'snes'
prm['snes_solver']['line_search'] = 'bt'
prm['snes_solver']['maximum_iterations'] = 1000
solver.solve()

# Sanity check; should be zero
print(assemble(ds_dx*dx) - float(L))

# Plot result:
from matplotlib import pyplot as plt
plot(y)
plt.show()
1 Like