TY - JOUR AU - Bartels, Sören AB - Abstract We consider variational problems that model the bending behavior of curves that are constrained to belong to given hypersurfaces. Finite element discretizations of corresponding functionals are justified rigorously via |$\varGamma $|-convergence. The stability of semi-implicit discretizations of gradient flows is investigated, which provide a practical method to determine stationary configurations. A particular application of the considered models arises in the description of conical sheet deformations. 1. Introduction The elastic flow of curves has attracted considerable attention among applied and numerical analysts within the last decades, cf., e.g., Langer & Singer (1985); Dziuk et al. (2002); Dall’Acqua et al. (2014) for analytical results and Deckelnick et al. (2005); Barrett et al. (2008, 2010, 2011, 2012, 2019); Deckelnick & Dziuk (2009); Bartels (2013); Pozzi & Stinner (2017); Bartels et al. (2018) for results concerning the discretization. Corresponding applications occur in the modeling of phase transitions, the description of large deformations of elastic rods and ribbons (Audoly & Pomeau, 2010), and prediction of prefered shapes of molecules (Coleman & Swigon, 2000; Chouaieb et al., 2006). For the class of inextensible curves, which arise naturally as dimensionally reduced descriptions in nonlinear elasticity (Mora & Müller, 2003; Antman, 2005), recent developments concerning the numerical treatment of partial differential equations with holonomic constraints such as harmonic maps turned out be useful for their efficient approximation, cf. Bartels (2005, 2016); Barrett et al. (2007). In this article we consider curves that are restricted to belong to given surfaces and whose behavior is determined by appropriate bending energies. To model their relaxation dynamics and find stationary configurations of low energy, we adapt techniques developed in Bartels (2013) to develop convergent finite element discretizations and stable iterative numerical schemes. Our approach provides an alternative to the methods developed in Barrett et al. (2012, 2019, 2020). Here, motivated by applications in nonlinear elasticity, we consider curves in Euclidean space that are parametrized by arclength, which allows for an efficient numerical treatment. Related analytical contributions are contained in Koiso (1996), Linnér (1991); Gerlach & von der Mosel (2011). 1.1 Constrained nonlinear bending We first consider relaxation processes of curves |$u$| on a given surface |$S$| whose bending behavior is determined by the functional $$\begin{align*} & I[u] = \frac12 \int_0^L |u^{\prime\prime}|^2 \,{\textrm d}x. \end{align*}$$ Here, we require |$u:(0,L)\to{\mathbb{R}}^3$| to be an arclength-parametrized curve, i.e., that |$|u^{\prime}(x)|=1$| for all |$x\in (0,L)$|⁠, so that |$|u^{\prime\prime}|^2$| is the squared curvature of the curve parametrized by the function |$u$|⁠. The functional arises from a rigorous dimension reduction for circular elastic rods from three-dimensional hyperelasticity, cf. Mora & Müller (2003). The constraint $$\begin{align*} & u(x) \in S \end{align*}$$ for all |$x\in (0,L)$| restricts the curve to belong to the regular hypersurface |$S\subset{\mathbb{R}}^3$|⁠. We also incorporate boundary conditions modeled by a bounded and linear functional |$L_{\textrm bc}:H^2(0,L;{\mathbb{R}}^3)\to{\mathbb{R}}^\ell $| and a vector |$\ell _{\textrm bc}\in{\mathbb{R}}^\ell $|⁠. The setting may describe the behavior of a wire on a magnetic surface neglecting effects related to twist. Corresponding torsion contributions can however be directly included, cf. Bartels & Reiter (2020). We thus consider the following constrained minimization problem. $$\begin{align*} \left\{\begin{array}{l} \text{Find a minimizing curve {$u \in H^2(0,L;{\mathbb{R}}^3)$} for} \\ \ \displaystyle{ I[u] = \frac12 \int_0^L |u^{\prime\prime}|^2 \,{\textrm d}x} \\[3mm] \textrm{subject to}\ \displaystyle{ u(x)\in S, \, |u^{\prime}(x)|^2 =1} \text{ for all {$x\in [0,L]$}} \\[3mm] \textrm{and}\ L_{\textrm bc}[u] = \ell_{\textrm bc}. \end{array}\right. \end{align*}$$($$\textrm {P}_{\textrm {bend}}$$) For an initial configuration described by a function |$u_0$| and for given boundary conditions, e.g., that the wire is clamped at one end, the relaxation of the bending energy is modeled by the formal gradient flow evolution $$\begin{align*} & \partial_t u = - I^{\prime}[u] + (\lambda u^{\prime})^{\prime} + \mu \varPhi_S^{\prime}(u) \end{align*}$$ for a family of curves |$(u(t))_{t\in [0,T]}$| satisfying the the initial, holonomic and boundary conditions $$\begin{align*} & u(0) = u_0, \quad |u^{\prime}|^2 = 1, \quad \varPhi_S(u) = 0, \quad L_{\textrm bc}[u] = \ell_{\textrm bc}. \end{align*}$$ The functions |$\lambda $| and |$\mu $| are Lagrange multipliers related to the arclength and surface constraints, respectively, where we assume that the surface |$S$| is given as the zero level set of the function |$\varPhi _S$|⁠. With the backward difference quotient operator $$\begin{align*} & d_t a^k = \frac1\tau(a^k-a^{k-1}) \end{align*}$$ we use a time-stepping scheme that linearizes the constraints at a previous approximation. By restricting to test functions that belong to the intersection of the kernels of the linearized constraints this eliminates the explicit occurence of the Lagrange multipliers. Since the time-derivative obeys the same linear constraints, we obtain for an appropriate inner product |$(\cdot ,\cdot )_*$| and the |$L^2$| inner product |$(\cdot ,\cdot )$| the time-stepping scheme $$\begin{equation*} (d_t u^k,v)_\ast + ([u^k]^{\ast\ast},v^{\ast\ast}) = 0 \end{equation*}$$ subject to the inclusions $$\begin{equation*} d_t u^k, \, v \in{\mathcal{F}}[u^{k-1}], \end{equation*}$$ where the set |${\mathcal{F}}[u^{k-1}]$| contains the linearized constraints, i.e., for a given curve |$\widehat{u}$|⁠, we have $$\begin{align*} & {\mathcal{F}}[\widehat{u}] = \big\{ v\in H^2(0,L;{\mathbb{R}}^3): \widehat{u}^{\prime}\cdot v = 0, \, \varPhi_S^{\prime}(\widehat{u})\cdot v = 0, \, L_{\textrm bc}[v] = 0 \big\}. \end{align*}$$ The time-stepping scheme thus requires solving linearly constrained linear systems of equations, where the constraints are pointwise. We show that the scheme is unconditionally energy decreasing and that the violation of the constraints is controlled by the step size independently of the number of iterations. Our spatial discretization uses an |$H^2$|-conforming ansatz, and imposes the constraints at the nodes of a partitioning of the reference interval |$(0,L)$|⁠. We justify the spatial discretization by proving its |$\varGamma $|-convergence to the continuous minimization problem. 1.2 Geodesic curvature An intrinsic variant of the constrained variational problem arises, e.g., in the description of phase separation processes on surfaces, cf. Adkins & Zhou (2017) for a model used to describe the formation of microdomains in lipid bilayer sheets. It replaces the curvature |$\kappa = |u^{\prime\prime}|$| by the geodesic curvature |$\kappa _g$|⁠. For an arclength parametrized curve |$u:(0,L)\to S$| it is defined as $$\begin{align*} & \kappa_g^2 = |u^{\prime\prime}|^2 - |u^{\prime\prime}\cdot n_S(u)|^2 = |u^{\prime\prime} \times n_S(u)|^2, \end{align*}$$ where |$n_S= \varPhi _S^{\prime}/|\varPhi _S^{\prime}|$| is a unit normal field on |$S$| and where we used that |$u^{\prime\prime}\cdot u^{\prime} = 0$|⁠. The corresponding energy functional $$\begin{align*} & I[u] = \frac12 \int_0^L \kappa_g^2 \,{\textrm d}s \end{align*}$$ still controls the |$H^2$| norm of |$u$| since the normal part of the curvature is bounded by the curvature of |$S$|⁠, i.e., we have $$\begin{align*} & |u^{\prime\prime}|^2 \le \kappa_g^2 + c_S^2, \end{align*}$$ where |$c_S$| is the maximum of the principal curvatures of |$S$|⁠. This estimate is not available when only nodal values of a piecewise polynomial curve |$u_h$| belong to |$S$|⁠. To cope with this aspect we introduce a stabilization via a damping parameter |$\gamma \le 1$| in the energy functional. $$\begin{align*} & \left\{\begin{array}{l} \text{Find a minimizing curve {$u \in H^2(0,L;{\mathbb{R}}^3)$} for} \\ \ \displaystyle{ I_\gamma[u] = \frac12 \int_0^L |u^{\prime\prime}|^2 - \gamma |u^{\prime\prime}\cdot n_S(u)|^2 \,{\textrm d}x} \\[3mm] \textrm{subject to}\ \displaystyle{ u(x)\in S, \, |u^{\prime}(x)|^2 =1} \text{ for all {$x\in [0,L]$}} \\[3mm] \textrm{and}\ L_{\textrm bc}[u] = \ell_{\textrm bc}. \end{array}\right. \end{align*}$$($${P_geo}^gamma$$) We prove that the stabilized problems converge in a variational sense to the unstabilized original problem as |$\gamma \to 1$|⁠. The stabilization allows us to prove convergence of discretizations. As an alternative to or in combination with stabilizations additional constraints may be imposed to ensure that discrete curves remain sufficiently close to the surface |$S$|⁠, so that their second derivative in normal direction is controlled by the curvature of the surface. This approach however leads to difficulties in the iterative solution. For the stabilized problem we follow the ideas described above with an explicit treatment of the nonlinear term. Hence, we compute a sequence |$(u^k)_{k=0,1,\dots }$| via the recursion $$\begin{align*} (d_t u^k&,v)_* + ([u^k]^{\prime\prime},v^{\prime\prime}) \\ & = \gamma \big([u^{k-1}]^{\prime\prime}\cdot n_S(u^{k-1}), v ^{\prime\prime}\cdot n_S(u^{k-1}) + [u^{k-1}]^{\prime\prime}\cdot n_S^{\prime}(u^{k-1}) v \big) \end{align*}$$ subject to |$d_t u^k,v\in{\mathcal{F}}[u^{k-1}]$|⁠. Under moderate conditions on the step size |$\tau $| in terms of |$\gamma $|⁠, we obtain a monotonicity property for the iteration. Fig. 1. Open in new tabDownload slide A point |$C$| of an initially flat elastic sheet (gray line representing cross section) is displaced by a distance |$\delta $|⁠. The resulting deformation is constrained by a circular obstacle at distance |$r$| to |$C$|⁠. Points on the deformed sheet (black lines) at distance |$(r^2+\delta ^2)^{1/2}$| to the center |$C$| either touch the obstacle (right end point) or are above it (left end point). Fig. 1. Open in new tabDownload slide A point |$C$| of an initially flat elastic sheet (gray line representing cross section) is displaced by a distance |$\delta $|⁠. The resulting deformation is constrained by a circular obstacle at distance |$r$| to |$C$|⁠. Points on the deformed sheet (black lines) at distance |$(r^2+\delta ^2)^{1/2}$| to the center |$C$| either touch the obstacle (right end point) or are above it (left end point). 1.3 Conical sheets Motivated by the problem of understanding folding and crumpling deformations of thin elastic sheets, the articles (Cerda & Mahadevan, 2005; Brandman et al., 2013; Müller & Olbermann, 2014; Olbermann, 2016; Figalli & Mooney, 2018) address the situation in which an elastic plate is placed on a circular obstacle of radius |$r$|⁠, and then indented by an amount |$\delta $| at the center |$C$|⁠. The resulting deformation is homogeneous along rays starting from the center, points at a distance |$(r^2+\delta ^2)^{1/2}$| from the center are either in contact with the obstacle or above it. The displacement of these points entirely determines the full deformation of the sheet, and it therefore suffices to compute the deformation of the points belonging to this circle. The displaced points belong to a sphere and are constrained by the obstacle. By an appropriate rescaling we may assume that |$r^2 + \delta ^2 = 1$|⁠. A cross section of the rotationally symmetric setting through the center |$C$| is depicted in Fig. 1. The solutions of the two-dimensional problem and its one-dimensional reduction cannot be rotationally symmetric unless the indentation depth |$\delta $| is trivial. The corresponding reduced description has been rigorously identified in Figalli & Mooney (2018) and characterizes the deformation |$u:S^1\to{\mathbb{R}}^3$| of the unit circle |$S^1 \subset{\mathbb{R}}^2$| via a minimization of the functional $$\begin{align*} & I[u] = \frac12 \int_{S^1} \kappa_g^2 \,{\textrm d}x \end{align*}$$ in the set of periodic curves |$u\in H^2(S^1;{\mathbb{R}}^3)$| subject to the constraints that |$u$| attains its values on the unit sphere |$S=S^2\subset{\mathbb{R}}^3$| and is inextensible, i.e., $$\begin{align*} & |u(x)|^2 = 1, \quad |u^{\prime}(x)|^2 = 1, \end{align*}$$ and that the curve does not penetrate the obstacle, i.e., for the vertical component |$u_3$| of |$u$|⁠, we have $$\begin{align*} & u_3(x) \ge \delta, \end{align*}$$ for all |$x\in S^1$|⁠. Because of the unit-length constraints on |$u$| and |$u^{\prime}$| we have that the normal curvature |$\kappa _n$| of |$u$| is given by $$\begin{align*} & \kappa_n = u^{\prime\prime}\cdot u = (u^{\prime}\cdot u)^{\prime} - |u^{\prime}|^2 = -1, \end{align*}$$ so that for the geodesic part we have $$\begin{align*} & \kappa_g^2 = \kappa^2 -\kappa_n^2 = |u^{\prime\prime}|^2 -1. \end{align*}$$ The reduced indentation problem thus leads to the following minimization problem for a given indentation depth |$\delta \ge 0$|⁠. $$\begin{align*} & \left\{ \begin{array}{l} \text{Find a minimizing curve {$u \in H^2(S^1;{\mathbb{R}}^3)$} for} \\ \ \displaystyle{ I[u] = \frac12 \int_{S^1} |u^{\prime\prime}|^2 \,{\textrm d}s - \pi} \\[3mm] \textrm{subject to}\ \displaystyle{ |u(x)|^2 = 1, \ |u^{\prime}(x)|^2 =1, \ u_3(x) \ge \delta} \text{ for all {$x\in S^1$}.} \end{array}\right. \end{align*}$$($$\textrm {P}_{\textrm {ind}}$$) Various features of minimizers have been characterized in Figalli & Mooney (2018), e.g., that the non-contact zone |$\{s\in S^1: u_3(s)> \delta \}$| is an interval. Via less rigorous arguments it has been stated in Cerda & Mahadevan (2005) that minimizers have, in a certain projection, a unique maximum, i.e., that single folds of the indented sheet are preferred over double folds, as is observed in reality. To investigate such questions via numerical experiments we approximate the problem by imposing the inequality constraint using a penalty approximation, i.e., we consider $$\begin{align*} & I_\varepsilon[u] = \frac12 \int_{S^1} |u^{\prime\prime}|^2 \,{\textrm d}x + \frac{1}{2\varepsilon} \int_{S^1} (u_3-\delta)_-^2 \,{\textrm d}x - \pi. \end{align*}$$ The minimization of |$I_\varepsilon $| is done with a gradient flow that linearizes the constraints, and which uses an implicit-explicit treatment of the penalty term defined via the convex-concave splitting $$\begin{align*} & (s-\delta)_-^2 = (s-\delta)^2 - (s-\delta)_+^2 \end{align*}$$ i.e., we compute a sequence |$(u^k)_{k=0,1,\dots }$| via $$\begin{align*} & (d_t u^k,v)_* + ([u^k]^{\prime\prime},v^{\prime\prime}) + \varepsilon^{-1} (u_3^k-\delta,v_3) = \varepsilon^{-1} ((u_3^{k-1}-\delta)_+,v_3) \end{align*}$$ for all |$v\in H^2(S^1;{\mathbb{R}}^3)$| subject to the linearized unit-length constraints and periodicity conditions contained in the space |${\mathcal{F}}[u^{k-1}]$| $$\begin{align*} & d_t u^k,v\in{\mathcal{F}}[u^{k-1}]. \end{align*}$$ The resulting iterative method is unconditionally energy monotone and converges to stationary configurations of low bending energy. To avoid case distinctions we will occasionally model closed curves in |$H^2(S^1;{\mathbb{R}}^3)$| as curves |$u\in H^2(0,L;{\mathbb{R}}^3)$| with |$L=2\pi $| satisfying the boundary conditions |$u(0)=u(L)$| and |$u^{\prime}(0)=u^{\prime}(L)$|⁠. 1.4 Outline The article is organized as follows. In Section 2 we introduce the finite element spaces used to approximate |$H^2$| curves and prove |$\varGamma $|-convergence results for the model problems. Section 3 is devoted to the development of stable gradient flow discretizations used to compute stationary configurations. In Section 4 we illustrate the theoretical findings by numerical experiments. 2. Discretization and |$\varGamma $|-convergence In this section we define suitable finite element spaces to approximate curves, devise discretizations of the constrained minimization problems and prove their variational convergence as discretization parameters tend to zero. 2.1 Finite element spaces We discretize the constrained minimization problems using |$H^2$| conforming finite element spaces for partitions $$\begin{align*} & 0 = z_0 < z_1 < \dots 1$|⁠. With the interpolation operator |${\mathcal{I}}_h$| we define discrete inner products and norms via $$\begin{align*} & (v,w)_h = \int_0^L {\mathcal{I}}_h [vw] \,{\textrm d}x, \quad \|v\|_{L^p_h(0,L)}^p = \int_0^L {\mathcal{I}}_h [|v|^p] \,{\textrm d}x \end{align*}$$ for |$v,w\in C([0,L];{\mathbb{R}}^\ell )$| and |$1\le p\le \infty $|⁠, where |$\|v\|_{L^\infty _h(0,L)} = \max _{j=0,\dots ,J} |v(z_h)|$|⁠. 2.2 Discrete minimization problems The pointwise constraints and the nonlinearities require making certain approximations which lead to inconsistency terms. We impose the arclength condition and the surface constraints at the nodes of a partitioning, i.e., we impose that $$\begin{align*} & {\mathcal{I}}_h |v_h^{\prime}|^2 = 1, \quad{\mathcal{I}}_h \varPhi_S(v_h) = 0, \end{align*}$$ which is equivalent to the nodal constraints $$\begin{align*} & |v_h^{\prime}(z_j)| = 1, \quad v_h(z_j) \in S \end{align*}$$ for |$j=0,1,\dots ,J$|⁠. The discrete set of admissible curves is then given by $$\begin{align*} & {\mathcal{A}}_h = \big\{ v_h \in V_h: {\mathcal{I}}_h|v_h^{\prime}|^2 = 1, \, {\mathcal{I}}_h \varPhi_S(v_h) = 0, \, L_{\textrm bc}[v_h] = \ell_{\textrm bc} \big\}. \end{align*}$$ It provides an approximation of the continuous set of admissible curves defined as $$\begin{align*} & {\mathcal{A}} = \big\{ v \in H^2(0,L;{\mathbb{R}}^3): |v^{\prime}|^2=1, \, \varPhi_S(v) =0, \, L_{\textrm bc}[v] = \ell_{\textrm bc} \big\}. \end{align*}$$ We note that if the continuous admissible set is nonempty then also the discrete admissible set is nonempty, i.e., we have the implication $$\begin{align*} & v\in{\mathcal{A}} \, \implies \, {\mathcal{I}}_h^{3,1} v \in{\mathcal{A}}_h, \end{align*}$$ where we assume that |$L_{\textrm bc}[v]$| only depends on the boundary values of |$v$| and |$v^{\prime}$|⁠. Our convergence result considers the minimization of $$\begin{align*} & I_\gamma[u] = \begin{cases} \displaystyle{\frac12 \int_0^L |u^{\prime\prime}|^2 - \gamma |u^{\prime\prime}\cdot n_S(u)|^2 \,{\textrm d}x} & \mbox{for } u \in{\mathcal{A}}, \\ + \infty & \mbox{for } H^2(0,L;{\mathbb{R}}^3)\setminus{\mathcal{A}}, \end{cases} \end{align*}$$ with a parameters |$\gamma \in [0,1)$|⁠. The approximating discrete functionals are given by $$\begin{align*} & I_{\gamma,h}[u_h] = \begin{cases} \displaystyle{\frac12 \int_0^L |u_h^{\prime\prime}|^2 - \gamma |u_h^{\prime\prime} \cdot n_S (u_h)|^2 \,{\textrm d}x} & \mbox{for } u_h \in{\mathcal{A}}_h, \\ + \infty & \mbox{for } H^2(0,L;{\mathbb{R}}^3) \setminus{\mathcal{A}}_h. \end{cases} \end{align*}$$ To prove the convergence |$I_{\gamma ,h}\to I_\gamma $| we impose a definiteness property on the boundary condition operator |$L_{\textrm bc}$| and an approximability condition on |${\mathcal{A}}$|⁠. Assumption 2.1 (Definiteness). The seminorm |$v\mapsto \|v^{\prime\prime}\|$| is a norm on the kernel of the operator |$L_{\textrm bc}:H^2(0,L;{\mathbb{R}}^3)\to{\mathbb{R}}^\ell $|⁠. The assumption is satisfied for clamped boundary conditions, e.g., |$L_{\textrm bc}[v] = (v(0),v^{\prime}(0))$| and boundary conditions that fix both end points, i.e., |$L_{\textrm bc}[v] = (v(0),v(L))$|⁠. We always assume that the boundary conditions lead to a nonempty set |${\mathcal{A}}$|⁠. Assumption 2.2 (Density of smooth curves). The subset of smooth curves |${\mathcal{A}}\cap H^3(0,L;{\mathbb{R}}^3)$| is dense in |${\mathcal{A}}$| with respect to strong convergence in |$H^2$|⁠. A relaxation of the assumption is discussed below in Remark 2.4. The assumption can be justified by regularizing curves in |${\mathcal{A}}$|⁠, projecting regular curves on |$S$|⁠, adjusting the boundary conditions and carrying out a suitable reparametrization. We refer the reader to Bartels & Reiter (2020) for related ideas. Proposition 2.3 (⁠|$\varGamma $|-convergence). If |$0\le \gamma <1$|⁠, |$\varPhi _S\in C^1({\mathbb{R}}^3)$| and Assumptions 2.1 and 2.2 are satisfied, then we have |$I_{\gamma ,h}\to I_\gamma $| in the sense of |$\varGamma $|-convergence with respect to weak convergence in |$H^2$|⁠, i.e., we have the following: (i) If |$(u_h)_{h>0} \subset H^2(0,L;{\mathbb{R}}^3)$| such that |$u_h\in{\mathcal{A}}_h$| for every |$h>0$| and |$I_{\gamma ,h}[u_h]\le c$|⁠, then there exists |$u\in{\mathcal{A}}$| such that |$u_h \rightharpoonup u$| in |$H^2$| and $$\begin{align*} & I_\gamma[u] \le \liminf_{h\to 0} I_{\gamma,h}[u_h]. \end{align*}$$ (ii) For every |$u\in{\mathcal{A}}$| there exists a sequence |$(u_h)_{h>0}\subset H^2(0,L;{\mathbb{R}}^3)$| such that |$u_h \to u$| in |$H^2$| and $$\begin{align*} & I_\gamma[u] = \lim_{h\to 0} I_{\gamma,h}[u_h]. \end{align*}$$ (iii) Weak accumulation points of sequences of quasiminimizers |$(u_h)_{h>0}$| for the functionals |$I_{\gamma ,h}$| in |$H^2$| are minimizers for |$I_\gamma $|⁠. Proof. (i) If |$I_{\gamma ,h}[u_h] \le c$| for a sequence |$(u_h)_{h>0}$| then, since |$\gamma <1$| and since $$\begin{align} |u_h^{\prime\prime}|^2 - \gamma |u_h^{\prime\prime}\cdot & n_S(u_h)|^2 \\ \nonumber & = (1-\gamma) |u_h^{\prime\prime}|^2 + \gamma \big|\big(I_3 - n_S(u_h)\otimes n_S(u_h)\big) u_h^{\prime\prime}\big|^2, \end{align}$$(1) we have that the sequence is bounded in |$H^2(0,L;{\mathbb{R}}^3)$| and there exists a weak limit |$u\in H^2(0,L;{\mathbb{R}}^3)$| of an appropriate subsequence which is not relabeled. The boundedness of the linear operator |$L_{\textrm bc}[v]$| shows that we have |$L_{\textrm bc}[u]=\ell _{\textrm bc}$|⁠. The compactness of the embedding |$H^2(0,L) \to W^{1,\infty }(0,L)$| implies that the sequence |$(u_h^{\prime})_{h>0}$| is strongly convergent in |$L^\infty (0,L;{\mathbb{R}}^3)$|⁠. Using that |${\mathcal{I}}_{\gamma ,h} |u_h^{\prime}|^2 = 1$| we thus deduce that $$\begin{align*} \big\||u_h^{\prime}|^2-1\big\|_{L^2(0,L)} &= \big\||u_h^{\prime}|^2-{\mathcal{I}}_h|u_h^{\prime}|^2\big\|_{L^2(0,L)} \\ &\le 2 ch \|u_h^{\prime}\|_{L^\infty(0,L)} \| u_h^{\prime\prime}\|_{L^2(0,L)}^2, \end{align*}$$ which implies that |$|u_h^{\prime}|^2 \to 1$| in |$L^2(0,L)$|⁠. We have that $$\begin{align*} \|\varPhi_S(u_h)\|_{L^\infty(0,L)} &= \|\varPhi_S(u_h) - {\mathcal{I}}_h \varPhi_S(u_h)\|_{L^\infty(0,L)} \\ &\le c h \|\varPhi_S^{\prime}(u_h) u_h^{\prime}\|_{L^\infty(0,L)}. \end{align*}$$ The pointwise convergence |$u_h\to u$| and continuity of |$\varPhi _S$| imply that |$\varPhi _S(u)=0$| in |$(0,L)$|⁠. Hence, we have that |$u\in{\mathcal{A}}$|⁠. Since $$\begin{align*} & P_{u_h} = I_3 - n_S(u_h)\otimes n_S(u_h) \to P_u = I_3 - n_S(u) \times n_S(u) \end{align*}$$ strongly in |$L^\infty (0,L;{\mathbb{R}}^{3\times 3})$| it follows that |$P_{u_h} u_h^{\prime\prime} \rightharpoonup P_u u^{\prime\prime}$| in |$L^2(0,L;{\mathbb{R}}^3)$| and the weak lower semicontinuity of the |$L^2$| norm in combination with the identity (1) shows that $$\begin{align*} & \int_0^L |u^{\prime\prime}|^2 - \gamma |u^{\prime\prime}\cdot n_S(u)|^2 \,{\textrm d}x \le \liminf_{h\to 0} \int_0^L |u_h^{\prime\prime}|^2 - \gamma |u_h^{\prime\prime}\cdot n_S(u_h)|^2 \,{\textrm d}x, \end{align*}$$ i.e., that |$I_\gamma [u]\le \liminf _{h\to 0} I_{\gamma ,h}[u_h]$|⁠. (ii) Since |$I_\gamma $| is continuous on |${\mathcal{A}}$| with respect to strong convergence in |$H^2$| and because of Assumption 2.2, we may assume that |$u\in{\mathcal{A}} \cap H^3(0,L;{\mathbb{R}}^3)$|⁠. Letting |$u_h = {\mathcal{I}}_h^{3,1}u$| we have that |$u_h \in{\mathcal{A}}_h$|⁠, |$u_h \to u$| in |$H^2$| and |$I_\gamma [u]= \lim _{h\to 0} I_{\gamma ,h}[u_h]$|⁠. (iii) The convergence of quasi-minimizers is an immediate consequence of the equicoercivity of the functionals |$I_{\gamma ,h}$| owing to the condition |$\gamma <1$| and assertions (i) and (ii). Remark 2.4 To avoid Assumption 2.2 one may impose the arclength and surface constraints in a relaxed sense in defining |${\mathcal{A}}_h$|⁠, i.e., using $$\begin{align*} \widetilde{{\mathcal{A}}}_h & = \big\{ v_h \in V_h:, \, L_{\textrm bc}[v_h] = \ell_{\textrm bc}, \\ & \qquad \|{\mathcal{I}}_h|v_h^{\prime}|^2 - 1\|_{L^\infty(0,L)} \le \alpha_h, \, \| {\mathcal{I}}_h \varPhi_S(v_h)\|_{L^\infty(0,L)} \le \beta_h \big\}, \end{align*}$$ with |$h$|-dependent parameters |$\alpha _h,\beta _h>0$|⁠. In this case, one may construct a recovery sequence |$u_h$| in part (ii) of the Proposition by letting |$\widetilde{u}\in C^\infty (0,L;{\mathbb{R}}^3)$| be a regularization of |$u\in{\mathcal{A}}$| which obeys the boundary conditions and define |$u_h = {\mathcal{I}}_h^{3,1} u$|⁠. If |$\alpha _h,\beta _h$| are appropriately chosen we have |$u_h\in \widetilde{{\mathcal{A}}}_h$| and |$u_h \to u$| in |$H^2$|⁠. 2.3 Application to model problems We next apply the abstract |$\varGamma $|-convergence result to the model problems defined by the variational problems (Pbend), (PgeodΓ) and (Pind). We assume throughout the following that |$\varPhi _S\in C^1({\mathbb{R}}^3)$| and that Assumptions 2.1 and 2.2 are satisfied, and always consider weak convergence in |$H^2$|⁠. The discretization of the constrained nonlinear bending problem (Pbend) is defined as follows: $$\begin{align*} & \left\{\begin{array}{l} \text{Find a minimizing curve {$u_h\in{\mathcal{A}}_h$} for} \\ \ \displaystyle{I[u_h] = \frac12 \int_0^L |u_h^{\prime\prime}|^2 \,{\textrm d}x.} \end{array}\right. \end{align*}$$($$\textrm {P}_{\textrm {bend}}^h$$) A convergence result is obtained from choosing |$\gamma =0$| in Proposition 2.3. Corollary 2.5 (Constrained nonlinear bending). The minimization problems (Pbendh) approximate the problem (Pbend) as |$h\to 0$|⁠. A discretization of the geodesic curvature minimization problem (PgeodΓ) is defined as $$\begin{align*} & \left\{\begin{array}{l} \text{Find a minimizing curve {$u_h\in{\mathcal{A}}_h$} for} \\ \ \displaystyle{I_{\gamma,h} [u_h] = \frac12 \int_0^L |u_h^{\prime\prime}|^2 - \gamma |u_h^{\prime\prime} \cdot n_S(u_h)|^2 \,{\textrm d}x.} \end{array}\right. \end{align*}$$($$\textrm {P}_{\textrm {geod}}^{\gamma ,h}$$) This problem approximates for fixed |$0<\gamma <1$| the stabilized problem (PgeodΓ) which is a direct consequence of Proposition 2.3. We also have that the regularized minimization problems converge for |$\gamma \to 1$| to the original, unstabilized problem defined with |$\gamma =1$|⁠. Corollary 2.6 (Geodesic curvature minimization). The minimization problems (PgeodΓ ,h) approximate problem (PgeodΓ) as |$h\to 0$|⁠. For |$\gamma \to 1$| problems (PgeodΓ) approximate problem (PgeodΓ) with |$\gamma =1$|⁠. Proof. The first part follows from Proposition 2.3. To prove the second part one uses that second derivatives of arclength-parametrized curves on |$S$| are bounded by their geodesic curvature. Remark 2.7 (i) For an efficient numerical realization it is helpful to replace the function |$u_h^{\prime\prime}\cdot n_S(u_h)$| by |$u_h^{\prime\prime}\cdot n_S(\overline{u}_h)$|⁠, where |$\overline{u}_h$| is a piecewise constant approximation of |$u_h$|⁠. The approximation result remains valid if |$u_h - \overline{u}_h \to 0$| in |$L^\infty (0,L;{\mathbb{R}}^3)$| for every bounded sequence |$(u_h)_{h>0}$| in |$H^2(0,L;{\mathbb{R}}^3)$|⁠, e.g., if |$\overline{u}_h$| is defined via the midpoint values of |$u_h$|⁠. (ii) A modification of the method is necessary to justify a joint limit passage |$(h,\gamma ) \to (0,1)$|⁠. In particular, control on the normal part of |$u_h^{\prime\prime}$| is needed, e.g., via requiring that |$u_h^{\prime}(z_j)$| is a tangent vector at every node |$z_j$|⁠, |$j=0,1,\dots ,J$|⁠. A discretization of the sheet indentation problem (Pind) is defined as $$\begin{align*} & \left\{\begin{array}{l} \text{Find a minimizing curve {$u_h\in{\mathcal{A}}_h$} for} \\ \ \displaystyle{I_{h,\varepsilon}[u_h] = \frac12 \int_0^L |u_h^{\prime\prime}|^2 + \frac1{2\varepsilon} \int_0^L {\mathcal{I}}_h (u_{3,h}-\delta)_-^2\,{\textrm d}x.} \end{array}\right. \end{align*}$$($$\textrm {P}_{\textrm {ind}}^{h,\varepsilon }$$) A convergence result is obtained from choosing |$\gamma =0$| in Proposition 2.3 and showing that the penalty term turns into a rigid constraint as |$(h,\varepsilon )\to 0$|⁠. Corollary 2.8 (Constrained nonlinear bending). Assume that Assumption 2.2 holds with |${\mathcal{A}}$| replaced by the set of functions |$u\in{\mathcal{A}}$| with |$u_3 \ge \delta $|⁠. Then the minimization problems (Pindh, Γ) approximate the problem (Pind) as |$(h,\varepsilon )\to 0$|⁠. Proof. Certain modifications of the proof of Proposition 2.3 are required. If the sequence |$(u_h)_{h>0}$| is such that |$I_{h,\varepsilon }[u_h] \le c$| then we have |$\|(u_{3,h}-\delta )_-\|_{L^2_h(0,L)}^2 \le 2 c \varepsilon $|⁠, and every weak accumulation point |$u\in H^2(0,L;{\mathbb{R}}^3)$| satisfies |$u_3\ge \delta $|⁠. Since the penalty term is non-negative we have that |$\liminf _{(h,\varepsilon )\to 0} I_{h,\varepsilon }[u_h] \ge I[u]$|⁠. For a curve |$u\in{\mathcal{A}}\cap C^\infty (0,L;{\mathbb{R}}^3)$| obeying the constraint |$u_3\ge \delta $| we have that the interpolants |$u_h = {\mathcal{I}}_h^{3,1}u$| also satisfy |${\mathcal{I}}_h u_{3,h} \ge \delta $| so that the penalty term in the functional disappears and the second part of the proof of Proposition 2.3 applies verbatimly. 3. Discrete gradient flows on surfaces We investigate in this section the stability of gradient flows for curvature energies defined on classes of arclength parametrized curves that belong to a given surface. The first model uses the full bending energy, the second one is defined by the geodesic curvature, while the third problem involves an obstacle constraint. 3.1 Constrained elastic flow of curves Minimizing the bending energy of curves restricted to a surface |$S$| subject to inextensibility and boundary conditions as formulated in problem (Pbend) leads to gradient flows such as $$\begin{align*} & \partial_t u = - u^{(4)} + (\lambda u^{\prime})^{\prime} + \mu \varPhi_S^{\prime} (u), \end{align*}$$ where |$\lambda $| and |$\mu $| are Lagrange multipliers related to inextensibility and surface constraints. More generally, given a metric |$(\cdot ,\cdot )_*$| defined on |$L^2(0,L;{\mathbb{R}}^3)$|⁠, we consider the evolution problem $$\begin{align*} & (\partial_t u,v)_* + (u^{\prime\prime},v^{\prime\prime}) = 0 \end{align*}$$ that determines a family |$u:[0,T] \to H^2(0,L;{\mathbb{R}}^3)$| of curves satisfying $$\begin{align*} & u(0)=u_0, \quad u(t)\in{\mathcal{A}} \end{align*}$$ for all |$t\in [0,T]$|⁠. We require the test functions |$v\in H^2(0,L;{\mathbb{R}}^3)$| to belong to the linearization of |${\mathcal{A}}$| at |$u(t)$|⁠, i.e., that |$v\in{\mathcal{F}}[u(t)]$|⁠, where $$\begin{align*} & {\mathcal{F}}[\widehat{u}] = \big\{ v\in H^2(0,L;{\mathbb{R}}^3): \varPhi_S^{\prime}(\widehat{u}) \cdot v = 0, \quad \widehat{u}^{\prime}\cdot v = 0, \quad L_{\textrm bc}[v]= 0 \big\}. \end{align*}$$ Note that also |$\partial _t u(t) \in{\mathcal{F}}[u(t)]$|⁠. To discretize the evolution equation we use a step size |$\tau>0$| and the backward difference operator $$\begin{align*} & d_t u^k = \frac{1}{\tau} (u^k-u^{k-1}). \end{align*}$$ For a partition |$z_00$|⁠. Proposition 3.2 (i) Algorithm 3.1 defines a sequence |$(u_h^k)_{k=0,1,...} \subset V_h$| such that for every |$K\ge 0$|⁠, we have $$\begin{align*} & I[u_h^K] + \tau \sum_{k=1}^K \|d_t u_h^k\|_*^2 \le I[u_h^0] = e_{0,h}. \end{align*}$$ (ii) Assume that |$u_h^0\in{\mathcal{A}}_h$| and that the inner product |$(\cdot ,\cdot )_*$| induces a norm |$\|\cdot \|_*$| with $$\begin{align*} & \|v_h^{\prime}\|_{L^\infty_h(0,L)}^2 = \|{\mathcal{I}}_h v_h^{\prime}\|_{L^\infty(0,L)}^2 \le c_* \|v_h\|_*^2 \end{align*}$$ for all |$v_h \in V_h$| and |$|\varPhi _S^{\prime\prime}(s)| \le c_{S,2}(1+ |s|^r)$| for all |$s\in{\mathbb{R}}^3$|⁠. Then, we have for every |$K \ge 0$| that $$\begin{align*} & \max_{k=0,1,\dots,K} \||[u_h^k]^{\prime}|^2-1\|_{L_h^\infty(0,L)} \le c_* \tau e_{0,h}, \end{align*}$$ and $$\begin{align*} & \max_{k=0,1,\dots,K} \|\varPhi_S(u_h^k)\|_{L_h^\infty(0,L)} \le c_* c_{S,2} c^{\prime} \tau e_{0,h}^{r+1}. \end{align*}$$ Proof. We test the formulation of Step (1) of Algorithm 3.1 with |$v_h = d_t u_h^k$| to deduce with a binomial formula that $$\begin{align*} & \|d_t u_h^k\|_*^2 + d_t \frac12 \|[u_h^k]^{\prime\prime}\|^2 + \frac{\tau}{2} \|[d_t u_h^k]^{\prime\prime}\|^2 = 0. \end{align*}$$ A summation over |$k=1,2,\dots ,K$| yields the asserted energy estimate. The nodewise orthogonality |$[d_t u_h^k]^{\prime} \cdot [u_h^{k-1}]^{\prime} = 0$| and the relation |$u_h^k = u_h^{k-1}+\tau d_t u_h^k$| imply that at every node |$z\in{\mathcal{N}}_h$|⁠, we have $$\begin{align*} & |[u_h^k]^{\prime}|^2 = |[u_h^{k-1}]^{\prime}|^2 + \tau^2 |[d_t u_h^k]^{\prime}|^2 = \dots = 1 + \tau^2 \sum_{\ell = 1}^k |[d_t u_h^k]^{\prime}|^2. \end{align*}$$ The energy bound and the assumed inequality for |$\|\cdot \|_*$| imply the bound for the arclength violation. For the surface constraint we note that the application of a Taylor formula and, the fact that |$d_t u_h^k \in{\mathcal{F}}_h[u_h^{k-1}]$| yield that at every node, we have $$\begin{align*} & \varPhi_S(u_h^k) = \varPhi_S(u_h^{k-1}) + \frac12 \tau^2 \varPhi_S^{\prime\prime}(\xi_h^k) [d_t u_h^k,d_t u_h^k]. \end{align*}$$ Repeating this argument and using |$\varPhi _S(u_h^0)=0$| at the nodes we infer with the assumed estimate for |$\varPhi _S^{\prime\prime}$| that $$\begin{align*} & \|{\mathcal{I}}_h \varPhi_S(u_h^k)\|_{L^\infty(0,L)} \le \frac12 \tau c_{S,2} \big(1+\|{\mathcal{I}}_h \xi_h^k\|_{L^\infty(0,L)}^r\big) \tau \sum_{\ell=1}^k \|{\mathcal{I}}_h d_t u_h^\ell\|_{L^\infty(0,L)}^2. \end{align*}$$ Since the nodal values of |$\xi _h^k$| belongs to the line segment connecting |$u_h^k$| and |$u_h^{k-1}$|⁠, we may incorporate the discrete |$L^\infty $| estimates to deduce the estimate for the nodewise surface constraint violation. 3.2 Geodesic curvature flow To develop an iterative scheme for the approximate solution of the geodesic curvature problem (PgeodΓ), we follow the ideas used for the constrained bending problem and use that $$\begin{align*} & \kappa_g^2 = |u^{\prime\prime}|^2 - |u^{\prime\prime}\cdot n_S(u)|^2. \end{align*}$$ To control the nonlinear second term by the first one, we introduce a stabilization via a damping factor |$\gamma _\varepsilon = (1-\varepsilon ^2)$|⁠. This leads to the functional $$\begin{align*} & I_\varepsilon[u] = \frac12 \int_0^L |u^{\prime\prime}|^2 - \gamma_\varepsilon |u^{\prime\prime}\cdot n_S(u)|^2 \,{\textrm d}x. \end{align*}$$ Because of the stabilization we have the implication $$\begin{align*} & I_\varepsilon[u] \le c_0 \quad \implies \quad \|u^{\prime\prime}\|^2 \le 2 c_0 \varepsilon^{-2}. \end{align*}$$ While on the continuous level the geodesic curvature of a curve on the surface |$S$| controls the full curvature, this is not the case for the discretization and hence necessitates the stabilization. We assume that $$\begin{align*} & n_S: {\mathbb{R}}^3 \to{\mathbb{R}}^3 \end{align*}$$ is a |$C^2$| vector field which coincides with the normal field on |$S$|⁠, i.e., we have |$n_S|_S = \frac{\varPhi _S^{\prime}(u)}{|\varPhi _S^{\prime}(u)|}$|⁠. We further assume that |$n_S$| has bounded derivatives. To simplify notation we use the mapping $$\begin{align*} & G_\varepsilon[u] = \frac{\gamma_\varepsilon}{2} \int_0^L |u^{\prime\prime}\cdot n_S(u)|^2 \,{\textrm d}x. \end{align*}$$ The constrained gradient flow for |$I_\varepsilon $| can thus be represented as $$\begin{align*} & (\partial_t u,v)_* + (u^{\prime\prime},v^{\prime\prime}) = G_\varepsilon^{\prime}[u;v], \end{align*}$$ where $$\begin{align*} & G_\varepsilon^{\prime}[u;v] = \gamma_\varepsilon \int_0^L u^{\prime\prime} \cdot n_S(u) \big(v^{\prime\prime} \cdot n_S(u) + u^{\prime\prime}\cdot n_S^{\prime}(u)v \big)\,{\textrm d}x. \end{align*}$$ We note that we have $$\begin{align*} & G_\varepsilon^{\prime}[u;v] \le \gamma_\varepsilon \big(\|u^{\prime\prime}\|\|v^{\prime\prime}\| + c_{n_S} \|u^{\prime\prime}\|^2 \|v\|_{L^\infty(0,L)}\big). \end{align*}$$ For ease of presentation we consider a semi-discrete setting. All arguments carry over to the case of a spatially discrete scheme. Algorithm 3.3 (Constrained geodesic curvature flow). Choose |$u^0\in V$| such that |$\varPhi _S(u^0) = 0$| and |$|[u^0]^{\prime}|^2 =1$| and |$L_{\textrm bc}[u^0] = \ell _{\textrm bc}$|⁠. Set |$k=0$|⁠. (1) Compute |$d_t u^k \in V$| such that $$\begin{align*} & (d_t u^k, v)_* + ([u^{k-1}+\tau d_t u^k]^{\prime\prime},v^{\prime\prime}) = G_\varepsilon^{\prime}[u^{k-1};v] \end{align*}$$ for all |$v \in V$| subject to the constraints $$\begin{align*} & d_t u^k, \, v \in{\mathcal{F}}[u^{k-1}]. \end{align*}$$ (2) Define |$u^k = u^{k-1} + \tau d_t u^k$|⁠; set |$k\to k+1$|⁠, and continue with (1). We have the following stability properties for Algorithm 3.3. Proposition 3.4 Assume that there exists |$c_*>0$| such that $$\begin{align*} & \|v\|_{L^\infty(0,L)} + \|v^{\prime\prime}\| \le c_* \|v\|_* \end{align*}$$ for all |$v\in V$|⁠. (i) There exists |$c_3 \ge 0$| such that if |$c_3 \tau \varepsilon ^{-1} \le 1/2$| then the iterates of Algorithm 3.3 satisfy for all |$K\ge 0$| $$\begin{align*} & I_\varepsilon[u^K] + (1-c_3 \tau \varepsilon^{-1}) \tau \sum_{k=1}^K \|d_t u^k\|_*^2 \le I_\varepsilon[u^0]. \end{align*}$$ (ii) Under the above condition the bounds on the constraint violation errors apply as in Proposition 3.2 (ii). Proof. We argue by induction and assume that the energy estimate and the constraint violation bounds have been established up to some number |$k-1\ge 0$| so that $$\begin{align*} & I_\varepsilon[u^{k-1}] + \frac\tau2 \sum_{\ell=1}^{k-1} \|d_t u^\ell\|_*^2 \le I_\varepsilon[u^0] =e_0. \end{align*}$$ This implies that $$\begin{align*} & \|[u^{k-1}]^{\prime\prime}\|\le \sqrt{2} e_0^{1/2} \varepsilon^{-1}. \end{align*}$$ By the assumption on the boundary data we thus have that |$\|u^{k-1}\|_{H^2(0,L)} \le c_1 \varepsilon ^{-1}$|⁠. Moreover, we have that |$\|u^{k-1}\|_{L^\infty (0,L)} \le c$|⁠. To derive an auxiliary bound we choose |$v=d_t u^k$| in Step (1) of Algorithm 3.3. Incorporating the bound for |$G_\varepsilon ^{\prime}$| and noting |$\gamma _\varepsilon \le 1$| this leads to $$\begin{align*} \|d_t u^k\|_*^2 & + d_t \frac12 \|[u^k]^{\prime\prime}\|^2 + \frac{\tau}{2} \|[d_t u^k]^{\prime\prime}\|^2 \\ & \le \|[u^{k-1}]^{\prime\prime}\| \|[d_t u^k]^{\prime\prime}\| + \|[u^{k-1}]^{\prime\prime}\|^2 c_{n_S} \|d_t u^k\|_{L^\infty(0,L)}. \end{align*}$$ By the assumption on the inner product |$(\cdot ,\cdot )_*$| we have that $$\begin{align*} & \|d_t u^k\|_{L^\infty(0,L)} + \|[d_t u^k]^{\prime\prime}\| \le c_* \|d_t u^k\|_* \end{align*}$$ and we deduce that $$\begin{align*} & \frac12 \|d_t u^k\|_*^2 + d_t \frac12 \|[u^k]^{\prime\prime}\|^2 \le c_1 \varepsilon^{-2}. \end{align*}$$ Hence, by choosing |$\tau $| sufficiently small, we have that $$\begin{align*} & \|[u^k]^{\prime\prime}\|^2 \le \|[u^{k-1}]^{\prime\prime}\|^2 + 2\tau c_1 \varepsilon^{-2} \le 5 e_0 \varepsilon^{-2}. \end{align*}$$ We next improve the latter bound by choosing again |$v=d_t u^k$| and using $$\begin{align*} & G_\varepsilon [u^k] - G_\varepsilon[u^{k-1}] = \tau G_\varepsilon^{\prime}[u^{k-1};d_t u^k] + \tau^2 G_\varepsilon^{\prime\prime}[\xi^k;d_t u^k,d_t u^k], \end{align*}$$ where |$G_\varepsilon ^{\prime\prime}[\xi ^k;d_t u^k,d_t u^k]$| is a formal representation of the Taylor remainder term $$\begin{align*} & {\mathcal{R}}_{G_\varepsilon}[u^{k-1},u^k;d_t u^k,d_t u^k] = \int_0^1 (1-s)G_\varepsilon^{\prime\prime}[u^{k-1}+s(u^k-u^{k-1});d_t u^k,d_t u^k] \,{\textrm d}s. \end{align*}$$ With the bounds for |$u^k$| and |$u^{k-1}$|⁠, we obtain that $$\begin{align*} & \big|G_\varepsilon^{\prime\prime}[\xi^k;d_t u^k,d_t u^k]\big| \le c_2 (1+\|[\xi^k]^{\prime\prime}\|) \|[d_t u^k]^{\prime\prime}\|^2 \le c_2^{\prime} \varepsilon^{-1} \|[d_t u^k]^{\prime\prime}\|^2. \end{align*}$$ We thus obtain that $$\begin{align*} \|d_t u^k\|_*^2 + d_t & \frac12 \|[u^k]^{\prime\prime}\|^2 + \frac12 \|[d_t u^k]^{\prime\prime}\|^2 = G_\varepsilon^{\prime}[u^{k-1};d_t u^k] \\ &= d_t G_\varepsilon[u^k] - \tau G_\varepsilon^{\prime\prime}[\xi^k;d_t u^k,d_t u^k] \le d_t G_\varepsilon[u^k] + \tau c_2^{\prime} \varepsilon^{-1} \|d_t u^k\|_*^2. \end{align*}$$ This proves the energy monotonicity and hence part (i) of the proposition. Part (ii) follows as in the proof of Proposition 3.2. 3.3 Conical sheet indentation flow To iteratively solve the reduced conical sheet indentation problem (Pind), we include the obstacle condition |$u_s(x)\ge \delta $| via a penalty term in the energy functional, i.e., $$\begin{align*} & I_\varepsilon[u] = \frac12 \int_{S^1} |u^{\prime\prime}|^2 \,{\textrm d}x + \frac{1}{2\varepsilon} \int_{S^1} (u_3-\delta)_-^2 \,{\textrm d}x, \end{align*}$$ where |$(s)_-=\min \{s,0\}$|⁠. The discretization of the related gradient flow $$\begin{align*} & (\partial_t u,v)_* + (u^{\prime\prime},v^{\prime\prime}) + \varepsilon^{-1} ((u_3-\delta)_-,v_3) = 0 \end{align*}$$ uses the convex–concave splitting $$\begin{align*} & (u_3-\delta)_-^2 = (u_3-\delta)^2 - (u_3-\delta)_+^2 \end{align*}$$ and an implicit treatment of the corresponding monotone and an explicit treatment of the corresponding antimonotone terms, i.e., we use the time-stepping scheme $$\begin{align*} & (d_t u^k,v)_* + ([u^k]^{\prime\prime},v^{\prime\prime}) + \varepsilon^{-1} (u_3^k-\delta,v_3) = \varepsilon^{-1} \big((u_3^{k-1}-\delta)_+,v_3\big). \end{align*}$$ A spatial discretization leads to the following algorithm where periodicity is guaranteed via an appropriate definition of the operator |$L_{\textrm bc}$|⁠. Algorithm 3.5 (Conical sheet flow). Choose |$u_h^0\in V_h$| such that |${\mathcal{I}}_h \varPhi _S(u_h^0) = 0$| and |${\mathcal{I}}_h |[u_h^0]^{\prime}|^2 =1$|⁠, and |${\mathcal{I}}_h u_h^0 \ge \delta $| and |$L_{\textrm bc}[u_h^0] = \ell _{\textrm bc}$|⁠. Set |$k=0$|⁠. (1) Compute |$d_t u_h^k \in V_h$| such that $$\begin{align*} &\begin{split} (d_t u_h^k, v_h)_* + ([u_h^{k-1}+\tau d_t u_h^k]^{\prime\prime},v_h^{\prime\prime}) + \varepsilon^{-1} & (u_{h,3}^k-\delta,v_{h,3})_h \\ & = \varepsilon^{-1} \big((u_{h,3}^{k-1}-\delta)_+,v_{h,3}\big)_h \end{split}\end{align*}$$ for all |$v_h \in V_h$| subject to the constraints $$\begin{align*} & d_t u_h^k, \, v_h \in{\mathcal{F}}_h[u_h^{k-1}]. \end{align*}$$ (2) Define |$u_h^k = u_h^{k-1} + \tau d_t u_h^k$|⁠; set |$k\to k+1$|⁠, and continue with (1). The iteration of Algorithm 3.5 has the same features as that of Algorithm 3.1. We use the discrete penalized energy functional $$\begin{align*} & I_{h,\varepsilon}[u_h] = \frac12 \int_{S^1} |u_h^{\prime\prime}|^2 \,{\textrm d}x + \frac{1}{2\varepsilon} \int_{S^1} {\mathcal{I}}_h (u_{h,3} -\delta)_-^2 \,{\textrm d}x. \end{align*}$$ Proposition 3.6 (i) Assume that |$u_h^0\in{\mathcal{A}}_h$| with |$u_{h,3}^0 \ge \delta $|⁠. Algorithm 3.5 defines a sequence |$(u_h^k)_{k=0,1,...} \subset V_h$| such that for every |$K\ge 0$|⁠, we have $$\begin{align*} & I_{h,\varepsilon}[u_h^K] + \tau \sum_{k=1}^K \|d_t u_h^k\|_*^2 \le I_{h,\varepsilon}[u_h^0] = e_{0,h}. \end{align*}$$ (ii) Under the above conditions the bounds on the constraint violation errors apply as in Proposition 3.2 (ii). Proof. We follow the steps of the proof of Proposition 3.2 and use |$v_h = d_t u_h^k$| in Step (1) of Algorithm 3.5. Defining the convex and concave functions |$p_{cx}$| and |$p_{cv}$|⁠, suitably embedded into |${\mathbb{R}}^3$|⁠, via $$\begin{align*} & p_{cx}(s) = (s_3-\delta)^2 e_3, \quad p_{cv}(s) = -(s_3-\delta)_+^2 e_3, \end{align*}$$ with the canonical basis vector |$e_3\in{\mathbb{R}}^3$|⁠, we thus have $$\begin{align*} \|d_t u^k\|_*^2 + d_t & \frac12 \|[u_h^k]^{\prime\prime}\|^2 + \frac{\tau}{2} \|[d_t u_h^k]^{\prime\prime}\|^2 \\ & = - \varepsilon^{-1} (p_{cx}^{\prime}(u_h^k),d_t u_h^k) - \varepsilon^{-1} (p_{cv}^{\prime}(u_h^{k-1}),d_tu_h^k). \end{align*}$$ The convexity of |$p_{cx}$| and |$-p_{cv}$| imply that we have $$\begin{align*} p_{cx}^{\prime}(u_h^k) \cdot (u_h^{k-1}-u_h^k) + p_{cx}(u_h^k) &\le p_{cx}(u_h^{k-1}), \\ -p_{cv}^{\prime}(u_h^{k-1}) \cdot (u_h^k-u_h^{k-1}) - p_{cv}(u_h^{k-1}) &\le -p_{cv}(u_h^k). \end{align*}$$ By adding the inequalities and dividing by |$\tau $|⁠, we find that $$\begin{align*} & - \big(p_{cx}^{\prime}(u_h^k) + p_{cv}^{\prime}(u_h^{k-1})\big) \cdot d_t u_h^k \le - d_t \big( p_{cx}(u_h^k) + p_{cv}(u_h^k)\big). \end{align*}$$ Combining the estimates implies the unconditional energy decay property. The remaining part (ii) is derived as in the proof of Proposition 3.2. Remark 3.7 To obtain a consistency property for the discrete gradient flow as an approximation of a corresponding continuous gradient flow, a condition relating the step-size |$\tau $| and the penalty parameter |$\varepsilon $| is required. 4. Numerical experiments We illustrate the performance of the numerical methods devised in the previous sections by various numerical experiments which are specified in the following subsections. The implementation of the algorithms was realized in Matlab with a direct solution of the linear systems of equations. The evolution metric |$(\cdot ,\cdot )_*$| was always chosen to coincide with the |$L^2$| inner product which leads to a mesh-dependent constant |$c_*$| in Propositions 3.2, 3.4 and 3.6. We observe however good stability properties for the resulting discrete |$L^2$| flow. Fig. 2. Open in new tabDownload slide Snapshots of the discrete gradient flow evolutions after |$n=0,20,40,\dots ,160$| iterations for the bending (solid curves) and geodesic curvature (dashed curves) energies from the same initial curve. Fig. 2. Open in new tabDownload slide Snapshots of the discrete gradient flow evolutions after |$n=0,20,40,\dots ,160$| iterations for the bending (solid curves) and geodesic curvature (dashed curves) energies from the same initial curve. Fig. 3. Open in new tabDownload slide Decay of the bending energy and geodesic curvature functional for the evolution of clamped curves on a torus illustrated in Fig. 2. Fig. 3. Open in new tabDownload slide Decay of the bending energy and geodesic curvature functional for the evolution of clamped curves on a torus illustrated in Fig. 2. Fig. 4. Open in new tabDownload slide Initial closed curve (dashed) and corresponding relaxed curves for bending energy and geodesic curvature with stabilization (nearly coinciding solid curves) and configuration for geodesic curvature flow without stabilization (dotted irregular curve). Fig. 4. Open in new tabDownload slide Initial closed curve (dashed) and corresponding relaxed curves for bending energy and geodesic curvature with stabilization (nearly coinciding solid curves) and configuration for geodesic curvature flow without stabilization (dotted irregular curve). Fig. 5. Open in new tabDownload slide Snapshots of the discrete gradient flow for the sheet indentation problem after |$n=0,10,20,30,40,70,190,430$| iterations. The family of curves (solid lines) relax their bending energies until only one fold is present which is stationary and energy minimizing. Fig. 5. Open in new tabDownload slide Snapshots of the discrete gradient flow for the sheet indentation problem after |$n=0,10,20,30,40,70,190,430$| iterations. The family of curves (solid lines) relax their bending energies until only one fold is present which is stationary and energy minimizing. Fig. 6. Open in new tabDownload slide Penetration of the obstacle at height |$\delta = 0.25$| (straight line) by the third component |$u_{3,h}^n$| of the iterates in the sheet indentation problem after a fixed number |$n=160$| of iterations for different choices of penalty parameters |$\varepsilon $|⁠. Fig. 6. Open in new tabDownload slide Penetration of the obstacle at height |$\delta = 0.25$| (straight line) by the third component |$u_{3,h}^n$| of the iterates in the sheet indentation problem after a fixed number |$n=160$| of iterations for different choices of penalty parameters |$\varepsilon $|⁠. 4.1 Elastic and geodesic flows on a torus We compare discrete relaxation dynamics for curves on a torus that are determined by the elastic bending energy and by the geodesic curvature functional. The torus |$T_{r,R}$| has radii |$R=2$| and |$r=1$|⁠, and is described by the zero level set of the function $$\begin{align*} & \varPhi_S(s) = (|s|^2 + R^2 - r^2)^2-4R^2(|s|^2-s_3^2). \end{align*}$$ The following example defines an open curve on |$T_{r,R}$|⁠. Example 4.1 Let |$\widetilde{L} = (7/4) \pi $|⁠, |$a=1/2$|⁠, |$b=1$| and for |$x\in (0,\widetilde{L})$| define |$\widetilde{u}^0:(0,\widetilde{L})\to T_{r,R}$| via $$\begin{align*} & \widetilde{u}^0(x) = \left[ \begin{array}{ccc} \sin(ax)\big(R+\sin(bx)r\big) \\ \cos(ax) \big(R+\sin(bx)r\big) \\ \cos(bx) r \end{array}\right]. \end{align*}$$ The curve |$u^0:(0,L) \to T_{r,R}$| is obtained from a re-parametrization of |$\widetilde{u}^0$|⁠. We use clamped boundary conditions at |$x=0$| that fix the initial position and tangent, i.e., we have $$\begin{align*} & L_{\textrm bc}[u] = (u(0),u^{\prime}(0)), \quad \ell_{\textrm bc} = (u^0(0),[u^0]^{\prime}(0)). \end{align*}$$ For a partition of the interval |$(0,L)$| we ran Algorithms 3.1 and 3.3 with the parameters $$\begin{align*} & J = 80, \quad h = L /J, \quad \tau = h, \quad \gamma = 1-h, \end{align*}$$ where |$L\approx 2.55\, \pi $|⁠. Figure 2 shows snapshots of the iterations. We observe that the curve changes quicker initially in the case of the bending energy and slightly slower for the geodesic curvature functional. This behavior is also seen in the energy plot shown in Fig. 3 where we plotted the energies in dependence of the iteration numbers. In Fig. 4 we illustrate for the evolution of closed curves on the torus the necessity of a stabilizing damping factor for the geodesic curvature flow. When no stabilization is used, i.e., in case that |$\gamma =1$|⁠, then energy monotonicity fails and the discrete curves fail to belong to a small neighborhood of the given surface. The initial closed curve used here was obtained by combining circular pieces with pieces of the curve defined in Example 4.1. 4.2 Conical sheet indentation We consider the following specification of the conical sheet indentation problem (Pind). Fig. 7. Open in new tabDownload slide Energy decay |$n\mapsto I_{h,\varepsilon }[u_h^n]$| in the sheet indentation problem for different spatial resolutions for mesh-dependent randomly generated initial configurations of large bending energy. Fig. 7. Open in new tabDownload slide Energy decay |$n\mapsto I_{h,\varepsilon }[u_h^n]$| in the sheet indentation problem for different spatial resolutions for mesh-dependent randomly generated initial configurations of large bending energy. Example 4.2 (Conical sheet). Let |$\delta = 1/4$| and |$r= (15/16)^2$|⁠. We used uniform partitions with mesh size |$h>0$| and nodes |$0=z_0<\dots < z_J=2\pi $| of the circle |$S^1$|⁠, where |$z_0$| and |$z_J$| are identified in the sense that we impose the periodic boundary condition |$L_{\textrm bc}[u]=0$| with $$\begin{align*} & L_{\textrm bc}[u] = \big(u(z_J)-u(z_0),u^{\prime}(z_J)-u^{\prime}(z_0)\big), \quad \ell_{\textrm bc} = 0. \end{align*}$$ Figure 5 shows snapshots of the discrete evolution computed with Algorithm 3.5 for the discretization parameters $$\begin{align*} & J = 80, \quad h = 2\pi/J, \quad \varepsilon = h^2, \quad \tau = h. \end{align*}$$ The visualization displays the two-dimensional deformation of the elastic sheet by linearly connecting the origin with points on the curve. The initial configuration |$u_h^0$| is a randomly generated function with corrected values to satisfy the condition |$u_h^0\in{\mathcal{A}}_h$|⁠. The discrete evolution shows a rapid change to a smooth curve approximately obeying the obstacle constraint. In the following iterations the number of local maxima decreases until finally only one fold can be observed, while the remaining part of the curve is in contact with the obstacle. Only a small penetration error occurs as can be seen in Fig. 6, where we plotted the third component of the iterates |$u_h^n$| with |$n$| such that |$t_n = n\tau = 2$|⁠, i.e., |$n=160$|⁠, for the choices $$\begin{align*} & \textrm{(i)} \quad \varepsilon = h, \quad \textrm{(ii)} \quad \varepsilon = h^2, \quad \textrm{(iii)} \quad \varepsilon = h^3. \end{align*}$$ For |$\varepsilon = h$| we observe a strong penetration of the obstacle. Our energy monotonicity property implies the estimate $$\begin{align*} & \|(u_{3,h}^n - \delta)_-\| \le (2e_{0,h})^{1/2} \varepsilon^{1/2} \end{align*}$$ and from the experimental results we infer that |$\varepsilon =h^2$| leads to the best results. It is also interesting to see how smaller penalization terms decrease the speed of the relaxation process. For |$\varepsilon =h$| only one fold is present indicating stationarity, while for |$\varepsilon = h^2$| and |$\varepsilon = h^3$| a larger number of local maxima can be observed after 160 iteration steps. Figure 7 shows the decay of the bending energy for different resolutions, and confirms the energy monotonicity and convergence to a stationary configuration. The large values of the energies are related to a strong dependence of minimal energies on the indentation depth |$\delta $|⁠. For the significantly smaller choice |$\delta = 0.05$| we obtained the stationary energy values |$I_{h,\varepsilon }[u_h^{n^*}] = 24.104, 17.828, 21.822, 21.569, 21.565$| for discretizations with |$J=40,80,\dots ,640$| grid points. These values also confirm convergence of the discrete minimal energies as |$h\to 0$|⁠. Acknowledgements The author is grateful to Rebecca Kromer for providing first versions of the numerical experiments. Dedicated to the memory of John W. Barrett References Adkins , M. R. & Zhou , Y. C. ( 2017 ) Geodesic curvature driven surface microdomain formation . J. Comput. Phys. , 345 , 260 – 274 . Google Scholar Crossref Search ADS PubMed WorldCat Antman , S. S. ( 2005 ) Nonlinear Problems of Elasticity , 2nd edn. Volume 107 of Applied Mathematical Sciences . New York : Springer . Audoly , B. & Pomeau , Y. (2010) Elasticity and Geometry . From hair curls to the non-linear response of shells, with a foreword by John W. Hutchinson . Oxford : Oxford University Press . Barrett , J. W. , Bartels , S., Feng , X. & Prohl , A. ( 2007 ) A convergent and constraint-preserving finite element method for the |$p$|-harmonic flow into spheres . SIAM J. Numer. Anal. , 45 , 905 – 927 . Google Scholar Crossref Search ADS WorldCat Barrett , J. W. , Garcke , H. & Nürnberg , R. ( 2008 ) Numerical approximation of anisotropic geometric evolution equations in the plane . IMA J. Numer. Anal. , 28 , 292 – 330 . Google Scholar Crossref Search ADS WorldCat Barrett , J. W. , Garcke , H. & Nürnberg , R. ( 2010 ) Numerical approximation of gradient flows for closed curves in |${\mathbb{R}}^d$| . IMA J. Numer. Anal., 30 , 4 – 60 . Barrett , J. W. , Garcke , H. & Nürnberg , R. ( 2011 ) The approximation of planar curve evolutions by stable fully implicit finite element schemes that equidistribute . Numer. Methods Partial Differential Equations , 27 , 1 – 30 . Google Scholar Crossref Search ADS WorldCat Barrett , J. W. , Garcke , H. & Nürnberg , R. ( 2012 ) Parametric approximation of isotropic and anisotropic elastic flow for closed and open curves . Numer. Math. , 120 , 489 – 542 . Google Scholar Crossref Search ADS WorldCat Barrett , J. W. , Garcke , H. & Nürnberg , R. ( 2019 ) Stable discretizations of elastic flow in Riemannian manifolds . SIAM J. Numer. Anal. , 57 , 1987 – 2018 . Google Scholar Crossref Search ADS WorldCat Barrett , J. W. , Garcke , H. & Nürnberg , R. ( 2020 ) Numerical approximation of curve evolutions in Riemannian manifolds . IMA J. Numer. Anal. , 40 , 1601 – 1651 . Google Scholar Crossref Search ADS WorldCat Bartels , S. ( 2005 ) Stability and convergence of finite-element approximation schemes for harmonic maps . SIAM J. Numer. Anal. , 43 , 220 – 238 . Google Scholar Crossref Search ADS WorldCat Bartels , S. ( 2013 ) A simple scheme for the approximation of the elastic flow of inextensible curves . IMA J. Numer. Anal. , 33 , 1115 – 1125 . Google Scholar Crossref Search ADS WorldCat Bartels , S. ( 2016 ) Projection-free approximation of geometrically constrained partial differential equations . Math. Comp. , 85 , 1033 – 1049 . Google Scholar Crossref Search ADS WorldCat Bartels , S. & Reiter , P. ( 2020 ) Numerical solution of a bending-torsion model for elastic rods . Numer. Math. , 146 , 661 – 697 . Bartels , S. , Reiter , P. & Riege , J. ( 2018 ) A simple scheme for the approximation of self-avoiding inextensible curves . IMA J. Numer. Anal. , 38 , 543 – 565 . Google Scholar Crossref Search ADS WorldCat Brandman , J. , Kohn , R. V. & Nguyen , H.-M. ( 2013 ) Energy scaling laws for conically constrained thin elastic sheets . J. Elasticity , 113 , 251 – 264 . Google Scholar Crossref Search ADS WorldCat Brenner , S. C. & Scott , L. R. ( 2008 ). The Mathematical Theory of Finite Element Methods , 3rd edn. Volume 15 of Texts in Applied Mathematics . New York : Springer . Cerda , E. & Mahadevan , L. ( 2005 ). Confined developable elastic surfaces: cylinders, cones and the elastica . Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. , 461 , 671 – 700 . Google Scholar OpenURL Placeholder Text WorldCat Chouaieb , N. , Goriely , A. & Maddocks , J. H. ( 2006 ) Helices . Proc. Nat. Acad. Sci. U.S.A. , 103 , 9398 – 9403 . Google Scholar Crossref Search ADS WorldCat Coleman , B. D. and Swigon , D. ( 2000 ) Theory of supercoiled elastic rings with self-contact and its application to DNA plasmids . J. Elasticity , 60 , 173 – 221 . Dall’Acqua , A. , Lin , C.-C. & Pozzi , P. (2014) Evolution of open elastic curves in |${\mathbb{R}}^n$| subject to fixed length and natural boundary conditions . Analysis (Berlin) , 34 , 209 – 222 . Deckelnick , K. & Dziuk , G. ( 2009 ) Error analysis for the elastic flow of parametrized curves . Math. Comp. , 78 , 645 – 671 . Google Scholar Crossref Search ADS WorldCat Deckelnick , K. , Dziuk , G. & Elliott , C. M. ( 2005 ) Computation of geometric partial differential equations and mean curvature flow . Acta Numerica , 14 , 139 – 232 . Google Scholar Crossref Search ADS WorldCat Dziuk , G. , Kuwert , E. & Schätzle , R. ( 2002 ) Evolution of elastic curves in |${\mathbb{R}}^n$|⁠: existence and computation . SIAM J. Math. Anal. , 33 , 1228 – 1245 . Figalli , A. & Mooney , C. ( 2018 ) An obstacle problem for conical deformations of thin elastic sheets . Arch. Rational Mech. Anal. , 228 , 401 – 429 . Google Scholar Crossref Search ADS WorldCat Gerlach , H. & von der Mosel , H. ( 2011 ) What are the longest ropes on the unit sphere? Arch. Rational Mech. Anal. , 201 , 303 – 342 . Google Scholar Crossref Search ADS WorldCat Koiso , N. ( 1996 ) On the motion of a curve towards elastica . Actes de la Table Ronde de Géométrie Différentielle (Luminy, 1992), volume 1 of Sémin. Congr . Paris : Soc. Math. France , pp. 403 – 436 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Langer , J. & Singer , D. A. ( 1985 ) Curve straightening and a minimax argument for closed elastic curves . Topology , 24 , 75 – 88 . Google Scholar Crossref Search ADS WorldCat Linnér , A. ( 1991 ) Curve-straightening in closed Euclidean submanifolds . Comm. Math. Phys. , 138 , 33 – 49 . Google Scholar Crossref Search ADS WorldCat Mora , M. G. & Müller , S. ( 2003 ) Derivation of the nonlinear bending-torsion theory for inextensible rods by |$\varGamma $|-convergence . Calc. Var. Partial Differential Equations , 18 , 287 – 305 . Google Scholar Crossref Search ADS WorldCat Müller , S. & Olbermann , H. ( 2014 ). Conical singularities in thin elastic sheets . Calc. Var. Partial Differential Equations , 49 , 1177 – 1186 . Google Scholar Crossref Search ADS WorldCat Olbermann , H. ( 2016 ) The one-dimensional model for d-cones revisited . Adv. Calc. Var. , 9 , 201 – 215 . Google Scholar Crossref Search ADS WorldCat Pozzi , P. & Stinner , B. ( 2017 ) Curve shortening flow coupled to lateral diffusion . Numer. Math. , 135 , 1171 – 1205 . Google Scholar Crossref Search ADS WorldCat © The Author(s) 2021. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) TI - Simulation of constrained elastic curves and application to a conical sheet indentation problem JO - IMA Journal of Numerical Analysis DO - 10.1093/imanum/drab008 DA - 2021-02-24 UR - https://www.deepdyve.com/lp/oxford-university-press/simulation-of-constrained-elastic-curves-and-application-to-a-conical-7O0oIXRFST SP - 1 EP - 1 VL - Advance Article IS - DP - DeepDyve ER -