# Discretization error estimates for penalty formulations of a linearized Canham–Helfrich-type energy

Discretization error estimates for penalty formulations of a linearized Canham–Helfrich-type... Abstract This article is concerned with minimization of a fourth-order linearized Canham–Helfrich energy subject to Dirichlet boundary conditions on curves inside the domain. Such problems arise in the modeling of the mechanical interaction of biomembranes with embedded particles. There, the curve conditions result from the imposed particle–membrane coupling. We prove almost-$$H^{\frac{5}{2}}$$ regularity of the solution and then consider two possible penalty formulations. For the combination of these penalty formulations with a Bogner–Fox–Schmit finite element discretization, we prove discretization error estimates that are optimal in view of the solution’s reduced regularity. The error estimates are based on a general estimate for linear penalty problems in Hilbert spaces. Finally, we illustrate the theoretical results by numerical computations. An important feature of the presented discretization is that it does not require the particle boundary to be resolved. This is crucial to avoid re-meshing if the presented problem arises as a subproblem in a model where particles are allowed to move or rotate. 1. Introduction A standard model for the behavior of biomembranes on a macroscale is the Canham–Helfrich model, which describes a biological membrane mathematically as a hypersurface $$\mathcal{M}\subseteq \mathbb{R}^3$$ minimizing the Canham–Helfrich energy   \begin{align*} \mathcal{J}_{\rm CHS}(\mathcal{M}) = \int_{\mathcal{M}} \frac{1}{2} \kappa H^2 + \kappa_{\rm G} K + \sigma \, \mathrm{d}\mathcal{H}^2\text{.} \end{align*} Here $$H$$ and $$K$$ are the mean and Gaussian curvature of $$\mathcal{M}$$, the coefficients $$\kappa>0$$ and $$\kappa_{\rm G}\geq 0$$ are the corresponding bending rigidities and $$\sigma \geq 0$$ is the membrane’s surface tension. Under the assumption that the membrane is ‘almost flat’, one can justify a geometric linearization of this functional for a membrane patch, leading to the so-called Monge–gauge approximation   \begin{align*} \mathcal{J}_{{\it{\Omega}}}(u) = \int_{\it{\Omega}} \frac{1}{2} \kappa ({\it{\Delta}} u)^2 + \frac{1}{2} \sigma \vert\nabla u\vert^2 \,\mathrm{{\rm d}}x \end{align*} of the Canham–Helfrich energy. Here the membrane patch is considered to be the graph $$\mathcal{M} = \{(x,u(x)) \mathrel\vert x \in {\it{\Omega}} \}$$ of a function $$u\colon {\it{\Omega}} \to \mathbb{R}$$ over some reference domain $${\it{\Omega}}\subseteq \mathbb{R}^2$$. A variety of hybrid models for the coupling of embedded particles to the membrane have been considered. These hybrid models are based on a continuous surface description of the membrane while particles are described by discrete entities (see, e.g., Helfrich & Jakobsson, 1990; Weikl et al., 1998; Dommersnes & Fournier, 1999, 2002; Bahrami et al., 2014). For an overview on hybrid models, we refer to Elliott et al. (2016) and the references cited therein. In this article we consider coupling conditions imposed on the particle boundaries and follow the notation introduced in Elliott et al. (2016). There, the reference domain $${\it{\Omega}}$$ is split into the membrane’s domain $${\it{\Omega}}_B \subseteq {\it{\Omega}}$$ and the particles’ domain $$B = {\it{\Omega}} \setminus {\it{\Omega}}_B$$. The model introduces membrane–particle interactions by functions that prescribe the membrane’s height profile and slope on the interface $${\it{\Gamma}} := \partial {\it{\Omega}}_B \cap \partial B$$. Altogether this yields an energy minimization problem subject to Dirichlet boundary conditions on $${\it{\Gamma}}$$. This article introduces and analyses a discretization based on a penalty formulation for the corresponding boundary value problem that avoids the resolution of particle boundaries. The reason for considering penalized boundary conditions is the following: since a biological membrane behaves like a fluid in tangential directions, particles can in principle move and rotate in plane. Any model simulating moving particles and any method computing optimal particle positions will thus have to solve multiple problems with varying particle positions. If the particle boundaries had to be resolved, this would require mesh deformation or re-meshing which can be computationally quite expensive. An alternative is to replace the strict boundary conditions by adequate penalty terms that can be formulated without having to resolve the boundary. One such penalty approach was introduced in Elliott et al. (2016) and is called soft curve formulation. A novel second formulation in this article is the soft bulk formulation. Both penalty problems will be defined later on. Our goal is to derive discretization error estimates for those penalty problems. For the second-order equations such estimates are well known (see, e.g., Nitsche, 1971; Babuška, 1973; Barrett & Elliott, 1986). It turns out that the fourth-order problems that we are interested in can be treated sufficiently well using a simple general framework for penalty problems in Hilbert spaces. This is mostly due to the fact that the regularity of our solutions is rather limited in the first place and so we can use rather simple estimates to still obtain optimal rates of convergence. Thanks to the abstract formulations, we get as a by-product a general error theory for finite element penalty problems with low regularity. The article is structured as follows. In Section 2, we introduce the notation and problem formulations that are used throughout this article. Section 3 is devoted to an abstract error result for linear penalty approximations on Hilbert spaces. As a foundation for the application of this result to the problem at hand, we then discuss regularity of solutions in Section 4. There we prove that a solution of the original problem lies in $$H^{\frac{5}{2}-\delta}({\it{\Omega}})$$ for all $$\delta>0$$. In Section 5, we combine the regularity with the developed abstract results to show the optimal (in view of the restricted regularity) convergence rate $${\mathcal O}(h^{1/2-\delta})$$ for a discretization with Bogner–Fox–Schmit finite elements. Finally, Section 6 illustrates our results by numerical examples that reproduce our theoretical findings. 2. Notation and problem formulations We consider a membrane with $$k$$ embedded particles. To this end, let $${\it{\Omega}} \subset \mathbb{R}^2$$ be a bounded reference domain with Lipschitz boundary $$\partial {\it{\Omega}}$$ and $$B_i \subset {\it{\Omega}}$$, $$i=1,\dots, k$$ the area occupied by the $$i$$th particle. For simplicity, we assume that the $$B_i$$ are closed, nonempty, connected and pairwise disjoint. Furthermore we require that each particle $$B_i$$ has a $$C^{1,1}$$-boundary. By $$B = \bigcup_{i=1}^k B_i$$ and $${\it{\Omega}}_B = {\it{\Omega}} \setminus B,$$ we denote the areas occupied by the particles and the membrane, respectively. Since the $$B_i$$ are closed and disjoint, the total membrane–particle interface is given by   \begin{align*} {\it{\Gamma}} = \bigcup_{i=1}^k \partial B_i = \partial {\it{\Omega}}_B \setminus \partial {\it{\Omega}} = \partial {\it{\Omega}}_B \cap B. \end{align*} For simplicity, we consider the space   \begin{align*} H = \{v \in H^2({\it{\Omega}}_B) \mathrel\vert v|_{\partial {\it{\Omega}}} = \partial_\nu v|_{\partial{\it{\Omega}}} = 0 \} \end{align*} with homogeneous Dirichlet boundary conditions on the boundary of $$\partial {\it{\Omega}}$$. Notice that we can treat Navier and periodic boundary conditions on $$\partial {\it{\Omega}}$$ using the same techniques (cf. Elliott et al., 2016). We assume that the membrane–particle interaction is governed by boundary conditions on the interface. More precisely, each particle $$B_i$$ enforces a height profile given by $$f^i_1\colon \partial B_i \to \mathbb{R}$$ and a slope given by $$f^i_2\colon \partial B_i \to \mathbb{R}$$ to the membrane on its boundary $$\partial B_i$$. That is, we consider boundary values   \begin{align} u|_{\partial B_i} &= f^i_1 + \gamma^i, & \partial_\nu u|_{\partial B_i} &= f^i_2, \end{align} (2.1) where $$\nu$$ is the unit outward normal to $$\partial {\it{\Omega}}_B$$, and the parameter $$\gamma^i \in \mathbb{R}$$ is allowed to vary freely to factor out the average height on $$\partial B_i$$. This is necessary because we want to prescribe the height profile only, while the absolute or average height is not fixed. For the following, we collect all such boundary data in a function $$f = (f_1,f_2) = \sum_{i=1}^{k} f^i\colon {\it{\Gamma}} \to \mathbb{R}^2$$, where $$f^i = (f^i_1,f^i_2)\colon \partial B_i \to \mathbb{R}^2$$ is extended by zero to the whole of $${\it{\Gamma}}$$. Using this notation, we consider the minimization problem: Problem 2.1 Find $$u\in H$$ and $$\gamma \in \mathbb{R}^k$$ minimizing $$\mathcal{J}_{{\it{\Omega}}_B}(u)$$ subject to   \begin{align} u|_{{\it{\Gamma}}} = f_1 + \sum_{i=1}^k \gamma_{i} \eta_{i}, \qquad \partial_\nu u|_{{\it{\Gamma}}} = f_2\text{ .} \end{align} (2.2) Notice that the effect of the free parameter $$\gamma_i$$ is localized to $$\partial B_i$$ via the use of the indicator function $$\eta_i := \chi_{\partial B_i} \in L^2({\it{\Gamma}}),$$ which is 1 on $$\partial B_i$$ but vanishes on all $$\partial B_j$$, $$j \neq i$$. In view of the trace theorem, we from now on assume that $$f = (f_1,f_2) \in H^{\frac{3}{2}}({\it{\Gamma}}) \times H^{\frac{1}{2}}({\it{\Gamma}})$$. Under this assumption it is known that Problem 2.1 admits a unique solution by application of Lax–Milgram’s theorem (Elliott et al., 2016). It is in fact possible to simplify this problem formulation. For this purpose, we make use of the trace operator $$T_{{\it{\Gamma}}} = (T_{{\it{\Gamma}}}^1, T_{{\it{\Gamma}}}^2)$$,   \begin{align*} T_{{\it{\Gamma}}}\colon H^2({\it{\Omega}}) & \rightarrow H^{\frac{3}{2}}({\it{\Gamma}}) \times H^{\frac{1}{2}}({\it{\Gamma}}), & \qquad T_{{\it{\Gamma}}}(v) = (v|_{{\it{\Gamma}}}, \ \partial_\nu v|_{{\it{\Gamma}}} ). \end{align*} This operator is well defined, continuous, surjective and admits a continuous right inverse (see e.g. Grisvard, 1985, Theorem 1.5.1.2). Note that we can also view $$T_{{\it{\Gamma}}}$$ as a trace operator for $$H^2({\it{\Omega}}_B)$$ due to our regularity assumptions on $${\it{\Gamma}}$$. To simplify the notation for boundary values up to the average height in (2.2), we furthermore introduce the projection operator   \begin{align*} P_{{\it{\Gamma}}} \colon H^{\frac{3}{2}}({\it{\Gamma}}) \times H^{\frac{1}{2}}({\it{\Gamma}}) &\longrightarrow H^{\frac{3}{2}}({\it{\Gamma}}) \times H^{\frac{1}{2}}({\it{\Gamma}}) \end{align*} defined by $$P_{\it{\Gamma}}(v_1,v_2) = \left(P_{\it{\Gamma}}^1(v_1),P_{\it{\Gamma}}^2(v_2)\right),$$ where   \begin{align*} P_{\it{\Gamma}}^1(v) = v - \sum_{i=1}^k |\partial B_i|^{-1}\int_{\partial B_i} v \,\mathrm{d}\sigma = v - \sum_{i=1}^k \frac{(v, \eta_i)_{L^2({\it{\Gamma}})}}{(\eta_i,\eta_i)_{L^2({\it{\Gamma}})}} \eta_i, \qquad P_{\it{\Gamma}}^2(v) = v. \end{align*} Using this operator, we can formulate the boundary conditions (2.2) as $$P_{\it{\Gamma}} (T_{\it{\Gamma}} u - f ) = 0,$$ which gives rise to the space of admissible functions defined by   \begin{align*} V_f := \left\{v \in H_0^2({\it{\Omega}}) \mid P_{\it{\Gamma}} (T_{\it{\Gamma}} v - f ) = 0 \right \} \end{align*} and the corresponding minimization problem: Problem 2.2 Find $$u\in V_f$$ minimizing $$\mathcal{J}_{{\it{\Omega}}}(u)$$. Notice that Problem 2.2 differs from Problem 2.1 not only due to the different notation for the boundary conditions on $${\it{\Gamma}}$$ but also because it minimizes $$J_{\it{\Omega}}$$ for functions defined on the whole of $${\it{\Omega}}$$ and not just on $${\it{\Omega}}_B$$. However, the following result from Elliott et al. (2016) shows that a solution of Problem 2.2 also immediately yields a solution to Problem 2.1. Proposition 2.3 Let $$u \in H_0^2({\it{\Omega}})$$ be the solution of Problem 2.2. Define $$\gamma_i = \vert \partial B_i\vert^{-1} (u-f_1, \eta_i )_{L^2({\it{\Gamma}})}$$ for $$i = 1,\dots,k$$. Then $$(u|_{{\it{\Omega}}_B},\gamma)$$ is the solution of Problem 2.1. While Problem 2.2 allows for variable height of particles, the full problem considered in Elliott et al. (2016) has additional degrees of freedom. This is due to the fact that particles can move and rotate in the plane of the fluid membrane. To avoid mesh deformation or re-meshing whenever particle positions change, we will in the following drop the hard constraints at the particle boundaries in favor of a penalized approach. Replacing the hard curve constraints by penalty terms in the energy functional leads to the following problem. Problem 2.4 Find $$u_\varepsilon \in H^2_0({\it{\Omega}})$$ minimizing   \begin{align*} \mathcal{J}_{{\it{\Omega}}}(u_\varepsilon ) + \sum_{i=1}^2 \frac{1}{\varepsilon_i} \left\Vert P_{\it{\Gamma}}^i ( T_{{\it{\Gamma}}}^i u_\varepsilon - f_i ) \right\Vert_{L^2({\it{\Gamma}})}^2\text{.} \end{align*} This formulation is more favorable than Problem 2.2 insofar as it admits a straightforward conforming finite element discretization without re-meshing in the case of variable particle positions. At this point, we want to mention an alternative penalty formulation that is based on the idea that, for known shapes $$B_i$$, the solution $$u|_B$$ could be computed a priori up to a constant per component $$B_i$$. We define the restriction operator   \begin{align*} T_B\colon H^2({\it{\Omega}}) & \longrightarrow H^{2}(B), & \qquad & T_B(v)= v|_B. \end{align*} Analogously to the curve constraint formulation, we introduce the associated projection operator   \begin{align*} P_B\colon H^2(B) & \longrightarrow H^2(B), & P_B(v) &= v - \sum_{i=1}^k \frac{(v,\psi_i)_{H^s(B)}}{(\psi_i,\psi_i)_{H^s(B)}} \psi_i, \end{align*} for $$\psi_i := \chi_{B_i}$$ and some fixed $$s \in [0,2]$$. Using this notation, the bulk constrained problem reads as follows. Problem 2.5 Find $$u\in H_0^2({\it{\Omega}})$$ minimizing $$\mathcal{J}_{{\it{\Omega}}}(u)$$ subject to   \begin{align*} P_B T_B u = P_B u|_B \text{.} \end{align*} One quickly verifies that the solutions of Problems 2.2 and 2.5 coincide. Analogously to Problem 2.4, a penalty formulation of Problem 2.5 is given by the following. Problem 2.6 Find $$u_\varepsilon \in H_0^2({\it{\Omega}})$$ minimizing   \begin{align*} \mathcal{J}_{{\it{\Omega}}}(u_\varepsilon) + \frac{1}{\varepsilon} \Vert P_B(T_Bu_\varepsilon - u|_B) \Vert_{H^s(B)}^2\text{.} \end{align*} While the penalized formulations are more flexible in the sense that the domain $${\it{\Omega}}_B$$ does not have to be resolved for discretization, one will be faced with the problem of balancing the penalty parameter and discretization errors. To this end, we will first develop an abstract error estimate for penalized problems and then analyse the regularity of the solution of the hard constrained Problem 2.2. Using a suitable regularity result later allows us to derive optimal $$h$$-dependent values of $$\varepsilon_i$$ for a discretization with mesh size $$h$$. 3. An error estimate for linear penalty problems on Hilbert spaces In this section, we derive an abstract energy error estimate for penalized discretizations of linearly constrained problems in a Hilbert space setting. Let $$H$$ be a Hilbert space and $$U_0 \subset H$$ a closed subspace. We consider the affine closed subspace $$U = U_0 + u_0 \subseteq H$$ of $$H$$ for some given $$u_0 \in H$$. In addition, let $$a\colon H \times H \rightarrow \mathbb{R}$$ be a symmetric positive semidefinite bounded bilinear form and $$\ell \colon H \rightarrow \mathbb{R}$$ a bounded linear form on $$H$$. We furthermore assume that $$a(\cdot,\cdot)$$ is coercive on $$U_0$$. In this setting, we consider the affine constrained minimization problem: Problem 3.1 Find $$u \in U$$ minimizing   \begin{align*} \tfrac{1}{2} a(u,u) - \ell(u)\text{.} \end{align*} Application of Lax–Milgram’s theorem yields the existence of a unique solution $$u \in U$$ of Problem 3.1 characterized by the variational equation   \begin{align} a(u, v) = \ell(v) \qquad \forall\,{v \in U_0}. \end{align} (3.1) From now on, $$u \in U$$ denotes this unique solution. Now let us assume that we want to approximate this problem by a penalty formulation over some closed linear subspace $$X \subseteq H$$ with $$m\in\mathbb{N}$$ penalty terms. Those penalty terms shall be given by symmetric positive semidefinite bounded bilinear forms $$b_i\colon H \times H \rightarrow \mathbb{R}$$ and penalty parameters $$\varepsilon_i >0$$. Denoting by $$\Vert v \Vert_{c} := \sqrt{ c(v,v) }$$ the seminorm induced by a symmetric positive semidefinite bilinear form $$c(\cdot,\cdot)$$ on $$H$$ the penalized problem reads as follows. Problem 3.2 Find $$u^X_\varepsilon \in X$$ minimizing   \begin{align*} \frac{1}{2} a\left(u^X_\varepsilon,u^X_\varepsilon\right) - \ell\left(u^X_\varepsilon\right) + \sum_{i=1}^m \frac{1}{2\varepsilon_i} \left\Vert u^X_\varepsilon-u\right\Vert_{b_i}^2. \end{align*} For the sake of convenience, we write   \begin{align} a_\varepsilon(\cdot,\cdot) = a(\cdot,\cdot) + \sum_{i=1}^m \frac{1}{\varepsilon_i} b_i(\cdot,\cdot), \qquad \ell_\varepsilon(\cdot) = \ell(\cdot) + \sum_{i=1}^m \frac{1}{\varepsilon_i} b_i(u,\cdot) \end{align} (3.2) for $$\varepsilon = (\varepsilon_i)_{i=1,\dots,m} \in \mathbb{R}_{+}^m = (0,\infty)^m$$. We require that $$a_1(\cdot,\cdot)$$ is coercive on $$X$$, such that $$a_\varepsilon(\cdot,\cdot)$$ is also coercive for all $$\varepsilon \in (0,1]^m$$. Under these assumptions, Lax–Milgram’s theorem implies existence of a unique solution $$u^X_\varepsilon \in X$$ for any $$\varepsilon \in (0,1]^m$$ characterized by the variational equation   \begin{align} a_\varepsilon\left(u^X_\varepsilon, v\right) = \ell_\varepsilon(v) \qquad \forall\, v \in X. \end{align} (3.3) The following result states a Céa-type estimate for the error $$u-u^X_\varepsilon$$ resulting from penalization and discretization in $$X$$. Theorem 3.3 For the fixed $$u \in U$$ used in the definition of $$\ell_\varepsilon$$, suppose there exist constants $$c_i >0$$ for $$i=1,\dots,m$$ such that   \begin{align} \vert a(u, v) - \ell(v) \vert \leq \sum_{i=1}^m c_i \Vert v\Vert_{b_i} \qquad \forall\,{v\in X}. \end{align} (3.4) Then the error for the solution $$u^X_\varepsilon$$ of Problem 3.2 can be bounded by   \begin{align} \Vert u- u^X_\varepsilon \Vert_{a_\varepsilon}^2 \leq 3 \inf_{v\in X} \left( \Vert u-v \Vert_{a_\varepsilon}^2 + \sum_{i=1}^m \varepsilon_i c_i^2 \right)\!\text{.} \end{align} (3.5) Proof. Let $$e = u - u^X_\varepsilon$$. Then we get   \begin{align} \Vert e \Vert_{a_\varepsilon}^2 = a_\varepsilon(e,e) = a_\varepsilon(e,u-v) + a_\varepsilon\left(e,v-u^X_\varepsilon\right) \end{align} (3.6) for all $$v\in X$$. Using Young’s inequality, we can bound the first term by   \begin{align} \begin{aligned} a_\varepsilon(e,u-v) &\leq \Vert e\Vert_{a_\varepsilon} \, \Vert u - v\Vert_{a_\varepsilon} \leq \tfrac{1}{4} \Vert e\Vert_{a_\varepsilon}^2 + \Vert u - v\Vert_{a_\varepsilon}^2. \end{aligned} \end{align} (3.7) The definition (3.2) of the $$u$$-dependent penalized functional $$\ell_\varepsilon$$ yields   \begin{align*} a_\varepsilon(u,w) - a(u,w) & = \ell_\varepsilon(w) - \ell(w) \qquad \forall\,{w \in X}. \end{align*} Combining this identity with (3.3) we get   \begin{align} a_\varepsilon(e,w) = a(u,w) - \ell(w) \qquad \forall\,{w \in X}. \end{align} (3.8) Because of $$v - u^X_\varepsilon \in X,$$ this implies   \begin{align} \begin{aligned} a_\varepsilon\left(e,v-u^X_\varepsilon\right) & = a\left(u,v-u^X_\varepsilon\right) - \ell\left(v-u^X_\varepsilon\right) \\ & \leq \sum_{i=1}^m c_i \left\Vert v - u^X_\varepsilon\right\Vert_{b_i} \\ & \leq \sum_{i=1}^m c_i \left( \Vert e\Vert_{b_i} + \Vert u - v \Vert_{b_i} \right) \\ & \leq \sum_{i=1}^m \left( \frac{1}{4\varepsilon_i} \Vert e\Vert_{b_i}^2 + \frac{1}{2\varepsilon_i} \Vert u - v\Vert_{b_i}^2 + \frac{3}{2}\varepsilon_i c_i^2\right) \\ & \leq \frac{1}{4}\Vert e\Vert_{a_\varepsilon}^2 + \frac{1}{2}\Vert u-v\Vert_{a_\varepsilon}^2 + \sum_{i=1}^m \frac{3}{2}\varepsilon_i c_i^2, \end{aligned} \end{align} (3.9) where the second to last inequality follows from Young’s inequality. Inserting the estimates (3.7) and (3.9) into (3.6), we obtain   \begin{align*} \Vert e\Vert_{a_\varepsilon}^2 \leq \frac{1}{2}\Vert e\Vert_{a_\varepsilon}^2 + \frac{3}{2}\Vert u-v\Vert_{a_\varepsilon}^2 + \sum_{i=1}^m \frac{3}{2}\varepsilon_i c_i^2, \end{align*} which proves the assertion. □ Next we consider the case where the evaluation of $$a(\cdot,\cdot)$$, $$\ell(\cdot)$$ and $$b_i(\cdot,\cdot)$$ is not performed exactly but approximated by some $$\widetilde{a}$$, $$\widetilde{\ell}$$ and $$\widetilde{b}_i$$, respectively. We use a notation analogous to $$a_\varepsilon$$ and $$\ell_\varepsilon$$:   \begin{align*} \widetilde{a}_\varepsilon(\cdot,\cdot) = \widetilde{a}(\cdot,\cdot) + \sum_{i=1}^m \frac{1}{\varepsilon_i} \widetilde{b}_i(\cdot,\cdot) , \qquad \widetilde{\ell}_\varepsilon(\cdot) = \widetilde{\ell}(\cdot) + \sum_{i=1}^m \frac{1}{\varepsilon_i} \widetilde{b}_i(u,\cdot). \end{align*} Instead of solving (3.3) for $$u^X_\varepsilon$$ directly, one now computes a $$\widetilde{u}^X_\varepsilon$$ by solving   \begin{align} \widetilde{a}_\varepsilon(\widetilde{u}^X_\varepsilon, v ) = \widetilde{\ell}_\varepsilon(v) \qquad \forall\,{v \in X} \text{.} \end{align} (3.10) In this setting, we can prove the following Strang-type result: Proposition 3.4 Let $$\left\Vert{a}\right\Vert$$ be the continuity constant of $$a$$ with respect to the $$H$$-norm and let $$a_1 = a+\sum_{i=1}^m b_i$$ and $$\widetilde{a}_1 = \widetilde{a} +\sum_{i=1}^m \widetilde{b}_i$$ be coercive with respect to the constants $$\alpha$$ and $$\widetilde{\alpha}$$, respectively. In addition to the assumptions of Theorem 3.3, suppose that for all $$i\in\{1,\dots,m\}$$ there exists $$c_i \in \mathbb{R}$$ such that for all $$v \in X$$,   \begin{align} \left\vert{ b_i(v,v) - \widetilde{b}_i(v,v) }\right\vert \leq c_i \Vert v \Vert^2\text{.} \end{align} (3.11) Then there exists a constant $$C = C(\alpha,\widetilde{\alpha},\left\Vert{a}\right\Vert) >0$$ such that   $$\matrix{ {{{\left\| {u - {\rm{ }}\tilde u_\varepsilon ^X} \right\|}_{{a_\varepsilon }}} \le C\left( {1 + \sum\limits_{i = 1}^m {{{{c_i}} \over {{\varepsilon _i}}}} } \right)\mathop {\inf }\limits_{v \in X} \mathop {\sup }\limits_{w \in X} {\rm{ }}\left[ {{{\left\| {u - v} \right\|}_{{a_\varepsilon }}} + {1 \over {\left\| w \right\|}}\left| {(a - \tilde a)(v,w)} \right|} \right.} \hfill \cr \hspace{20mm} {\quad \left. { + {1 \over {\left\| w \right\|}}\sum\limits_{i = 1}^m {{1 \over {{\varepsilon _i}}}} \left| {({b_i} - {{\tilde b}_i})(u - v,w)} \right| + {1 \over {\left\| w \right\|}}\left| {(\ell - \tilde \ell )(w)} \right|} \right].} \hfill \cr }$$ Proof. See the Appendix. □ 4. Regularity of the hard curve constraint problem To apply the results of the previous section, it will be crucial to prove (3.4). For this purpose and to prove convergence rates by bounding the best approximation errors in Theorem 3.3, we will now derive regularity results for the solution of Problem 2.2. Let $$u \in V_f$$ be the solution of the hard curve minimization formulation, Problem 2.2. Knowing the regularity of $$u$$ is central to proving discretization errors since the maximal regularity of the solution immediately reveals the optimal rates of convergence that one would expect for the corresponding discretization errors. Our strategy is to rewrite the hard curve minimization problem as a system of elliptic partial differential equations to which we apply standard regularity theory for elliptic interface problems. To this end, we define the fourth-order elliptic differential operator $$L := \kappa {\it{\Delta}}^2 - \sigma {\it{\Delta}}$$ associated with the energy functional $$\mathcal{J}$$. Proposition 4.1 Suppose $$(u_1,u_2) \in H^2({\it{\Omega}}_B) \times H^2(B)$$ is a weak solution of the system of partial differential equations   \begin{align} \begin{aligned} Lu_1 & = 0 \text{ on}\,{\it{\Omega}}_B, & T_{\it{\Gamma}} u_1 & = T_{\it{\Gamma}} u \text{ on}\, {\it{\Gamma}}, & u_1 = \partial_\nu u_1 = 0 \text{ on}\, \partial{\it{\Omega}}, \\ Lu_2 & = 0 \text{ on}\, B, & T_{\it{\Gamma}} u_2 & = T_{\it{\Gamma}} u \text{ on}\,{\it{\Gamma}}. \end{aligned} \end{align} (4.1) Then $$u_1 = u|_{{\it{\Omega}}_B}$$ and $$u_2 = u|_B$$ where $$u$$ is the solution of Problem 2.2. Conversely, if $$u$$ is the solution of Problem 2.2, then $$(u|_{{\it{\Omega}}_B}, u|_B)$$ solves (4.1). Proof. A proof for the equivalence of Problem 2.2 with a weak formulation of (4.1) is given in Elliott et al. (2016, Proposition 4.1). □ Lemma 4.2 Let $${\it{\Omega}}$$ be a piecewise polygonal domain, let $${\it{\Gamma}}$$ be smooth and suppose $$(f_1,f_2) \in H^{\frac{7}{2}}({\it{\Gamma}})\times H^{\frac{5}{2}}({\it{\Gamma}})$$. If all corners of $${\it{\Omega}}$$ have an inner angle $$\omega$$ with $$\omega \leq 126^\circ$$ then $$(u|_{{\it{\Omega}}_B},u|_B) \in H^{4}({\it{\Omega}}_B)\times H^{4}(B)$$. Proof Let $$g := \frac{\sigma}{\kappa} {\it{\Delta}} u|_{{\it{\Omega}}_B} \in L^2({\it{\Omega}}_B)$$. From Proposition 4.1, we know that $$u|_{{\it{\Omega}}_B}$$ is a weak solution of   \begin{align*} {\it{\Delta}}^2 u = g \text{ on}\, {\it{\Omega}}_B, \qquad T_{{\it{\Gamma}}}u = T_{{\it{\Gamma}}} u|_{{\it{\Omega}}_B} \text{ on}\, {\it{\Gamma}}, \qquad u = \partial_\nu u = 0 \text{ on}\, \partial{\it{\Omega}}. \end{align*} Since $$g \in L^2({\it{\Omega}}_B)$$ and $$T_{{\it{\Gamma}}}^1 u|_{{\it{\Omega}}_B} = f_1 + \sum_{i=1}^k \vert \partial B_i\vert^{-1} ( u|_{{\it{\Omega}}_B} - f_1, \eta_i)_{L^2({\it{\Gamma}})}\eta_i \in H^{\frac{7}{2}}({\it{\Gamma}})$$ as well as $$T_{{\it{\Gamma}}}^2 u|_{{\it{\Omega}}_B} = f_2 \in H^{\frac{5}{2}}({\it{\Gamma}}),$$ it follows from Grisvard (1985, Theorem 7.2.2.3) and the computations in Blum & Rannacher (1980) on the associated characteristic equation that $$u|_{{\it{\Omega}}_B} \in H^4({\it{\Omega}}_B)$$. An analogous argument on $$B$$ yields $$u|_B \in H^4(B)$$ and thus proves the assertion. □ Having the intrinsic regularity $$u \in H^2({\it{\Omega}})$$ and the piecewise regularity from Lemma 4.2 allows us to show an improved global regularity result. The key ingredient is the following technical lemma. Lemma 4.3 Let $$v \in H^2({\it{\Omega}})$$ such that $$v|_{{\it{\Omega}}_B} \in H^4({\it{\Omega}}_B)$$ and $$v|_{B} \in H^4(B)$$. Then $$v \in H^{2+\frac{1}{2}-\delta}({\it{\Omega}})$$ for all $$\delta > 0$$. Proof. See the Appendix. □ Combining Lemmas 4.2 and 4.3 now gives a corollary. Corollary 4.4 Let the assumptions from Lemma 4.2 hold. Then $$u \in H^{\frac{5}{2}-\delta}({\it{\Omega}})$$ for every $$\delta >0$$. 5. Discretization errors In this section, we apply the results from the previous sections to derive discretization errors for penalized finite element approximations of the form of Problems 2.4 and 2.6 but with $$h$$-dependent penalty parameters. We suppose from now on that $${\it{\Omega}}$$ is a rectangular domain. In view of Lemma 4.2, this implies that $$u|_{{\it{\Omega}}_B} \in H^4({\it{\Omega}}_B)$$, $$u|_B \in H^4(B)$$ and $$u \in H^{\frac{5}{2}-\delta}({\it{\Omega}})$$ for any $$\delta >0$$. On the domain $${\it{\Omega}}$$ we establish a quadrilateral grid equipped with Bogner–Fox–Schmit finite elements (see, e.g., Ciarlet, 1978). Given the set of grid nodes $$N$$ and the set of multiindices $${\it{\Lambda}} = \{(0,0), (1,0), (0,1), (1,1)\}$$, the resulting finite element space is spanned by a basis $$(\psi_{p,\alpha})_{p\in N, \alpha \in {\it{\Lambda}}}$$ of piecewise bicubic polynomials such that each $$\psi_{p,\alpha}$$ satisfies   \begin{align*} \partial^\beta \psi_{p,\alpha}(q) = \delta_{\alpha,\beta} \delta_{p,q} \qquad \forall\, q \in N, \beta \in {\it{\Lambda}}, \end{align*} where $$\partial^\beta$$ is the partial derivative associated with the multiindex $$\beta,$$ and $$\delta_{a,b}$$ is the Kronecker delta. That is, the degrees of freedom are the values, first-order partial derivatives and mixed second-order partial derivatives at the vertices (see Fig. 1). Fig. 1. View largeDownload slide Degrees of freedom for the Bogner–Fox–Schmit element. Fig. 1. View largeDownload slide Degrees of freedom for the Bogner–Fox–Schmit element. By setting the degrees of freedom on $$\partial{\it{\Omega}}$$ to zero this yields a conforming subspace of $$C^1(\overline{{\it{\Omega}}}) \cap H_0^2({\it{\Omega}})$$. Given a family $$(S_h)_{h\in I}$$ of such finite element discretizations over $${\it{\Omega}}$$ based on quasi-uniform grids with mesh size $$h,$$ we will use the following well-known result: there exists a constant $$c>0$$ such that for all $$K \in [0,4],$$ the approximation estimate   \begin{align} \forall\,{v \in H^K({\it{\Omega}}) \cap H_0^2({\it{\Omega}})},\,\exists{\overline{v} \in S_h}, \, \forall\,{k\in [0,\min\{2,K\}]}, \Vert v - \overline{v} \Vert_{H^k({\it{\Omega}})} \leq c h^{K-k} \Vert v \Vert_{H^K({\it{\Omega}})} \end{align} (5.1) holds true for all $$h \in I$$. (This is a direct consequence of classical approximation estimates (see, e.g., Ciarlet, 1978, Theorem 3.1.4), and interpolation theory of function spaces, (see e.g., Lunardi, 2009, Theorem 1.1.6).) To put our minimization problems into the variational framework used in Section 3, we introduce the bilinear form   \begin{align*} a \colon H^2({\it{\Omega}}) \times H^2({\it{\Omega}}) &\longrightarrow \mathbb{R}, & a(w,v) &= \kappa ( {\it{\Delta}} w, {\it{\Delta}} v )_{L^2({\it{\Omega}})} + \sigma ( \nabla w, \nabla v )_{L^2({\it{\Omega}})}, \end{align*} together with the curve penalty terms   \begin{align*} b_{{\it{\Gamma}}}^1 \colon H^2({\it{\Omega}}) \times H^2({\it{\Omega}}) & \longrightarrow \mathbb{R}, & b_{{\it{\Gamma}}}^1(w,v) &= ( P_{\it{\Gamma}}^1T_{\it{\Gamma}}^1 w, P_{\it{\Gamma}}^1T_{\it{\Gamma}}^1 v )_{L^2({\it{\Gamma}})},\\ b_{{\it{\Gamma}}}^2 \colon H^2({\it{\Omega}}) \times H^2({\it{\Omega}}) & \longrightarrow \mathbb{R}, & b_{{\it{\Gamma}}}^2(w,v) &= ( P_{\it{\Gamma}}^2T_{\it{\Gamma}}^2 w, P_{\it{\Gamma}}^2T_{\it{\Gamma}}^2 v )_{L^2({\it{\Gamma}})} \end{align*} and the bulk penalty term   \begin{align*} b_{B}\colon H^2({\it{\Omega}}) \times H^2({\it{\Omega}}) & \longrightarrow \mathbb{R}, & b_{B}(w,v) &= ( P_BT_B w, P_BT_B v )_{H^s(B)}\text{.} \end{align*} Then the solution $$u \in V_f \subset H_0^2({\it{\Omega}})$$ of Problem 2.2 satisfies the variational equation   \begin{align} a(u,v) = 0 \qquad \forall\,{v \in V_0}. \end{align} (5.2) Discretization of the penalized curve constrained Problem 2.4 in $$S_h$$ leads to the soft curve problem: find $$u_\varepsilon^h \in S_h$$ minimizing   \begin{align} \mathcal{J}_{{\it{\Omega}}}(u_\varepsilon^h) + \sum_{i=1}^2 \frac{1}{\varepsilon_i} \left\Vert P_{\it{\Gamma}}^i (T_{{\it{\Gamma}}}^i u_\varepsilon^h - f_i)\right \Vert_{L^2({\it{\Gamma}})}^2\text{,} \end{align} (5.3) which is equivalently characterized by the variational equation   \begin{align} u_\varepsilon^h \in S_h: \qquad a\left(u^h_\varepsilon,v\right) + \sum_{i=1}^2 \frac{1}{\varepsilon_i} b_{\it{\Gamma}}^i\left(u^h_\varepsilon-u,v\right) = 0 \qquad \forall\,{v \in S_h}. \end{align} (5.4) Analogously, discretization of the penalized bulk constrained Problem 2.6 in $$S_h$$ leads to the soft bulk problem: find $$u_\varepsilon^h \in S_h$$ minimizing   \begin{align} \mathcal{J}_{{\it{\Omega}}}\left(u_\varepsilon^h\right) + \frac{1}{\varepsilon} \left\Vert P_BT_B(u_\varepsilon^h - u)\right \Vert_{H^s(B)}^2 \end{align} (5.5) with the equivalent variational equation   \begin{align} u_\varepsilon^h \in S_h: \qquad a\left(u^h_\varepsilon,v\right) + \frac{1}{\varepsilon} b_B\left(u^h_\varepsilon-u,v\right) = 0 \qquad \forall\,{v \in S_h}. \end{align} (5.6) At this point, we want to emphasize that the solutions of (5.4) and (5.6) are in general different. We nevertheless refer in a slight abuse of notation to both solutions as $$u^h_\varepsilon$$. In the following it will, however, be clear from the context whether $$u^h_\varepsilon$$ denotes the solution of the soft curve problem or the soft bulk problem. In the next result, we show that the soft curve formulation meets the requirement from Theorem 3.3, which we will use afterwards to get an asymptotic error estimate. Lemma 5.1 For all $$v \in H_0^2({\it{\Omega}})$$ the solution $$u$$ of Problem 2.2 satisfies   \begin{align*} a(u,v) \leq \kappa \left(\Vert [{\it{\Delta}} u] \Vert_{L^2({\it{\Gamma}})} + \Vert [\partial_\nu {\it{\Delta}} u]\Vert_{L^2({\it{\Gamma}})} \right) \left( \Vert v\Vert_{b^1_{\it{\Gamma}}} + \Vert v \Vert_{b^2_{\it{\Gamma}}} \right)\text{.} \end{align*} Here $$[w] = w|_{{\it{\Omega}}_B} - w|_B$$ denotes the jump of the function $$w$$ on $${\it{\Gamma}}$$. Proof. By integration by parts and as $$u \in H^2({\it{\Omega}})$$, $$u|_{{\it{\Omega}}_B} \in H^4({\it{\Omega}}_B)$$, $$u|_B \in H^4(B)$$ and $$Lu = 0$$, we are able to state for all $$v \in H_0^2({\it{\Omega}})$$,   \begin{align*} \frac{1}{\kappa} a(u,v) & = \int_{{\it{\Omega}}} {\it{\Delta}} u {\it{\Delta}} v + \frac{\sigma}{\kappa} \nabla u \cdot \nabla v \\ & = \int_{{\it{\Omega}}_B} \frac{1}{\kappa} Lu|_{{\it{\Omega}}_B} \, v + \int_B \frac{1}{\kappa} Lu|_B \, v + \int_{{\it{\Gamma}}} \frac{\sigma}{\kappa} \left( \partial_\nu u|_{{\it{\Omega}}_B} - \partial_\nu u|_B \right) v \\ & \quad + \int_{{\it{\Gamma}}} \left( - \partial_\nu {\it{\Delta}} u|_{{\it{\Omega}}_B} + \partial_\nu {\it{\Delta}} u|_B \right) v + \int_{{\it{\Gamma}}} \left( {\it{\Delta}} u|_{{\it{\Omega}}_B} - {\it{\Delta}} u|_B \right ) \partial_\nu v \\ & = \left( -[\partial_\nu {\it{\Delta}} u], T_{\it{\Gamma}}^1 v\right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u], T_{\it{\Gamma}}^2 v\right)_{L^2({\it{\Gamma}})} \\ & = \left( -[\partial_\nu {\it{\Delta}} u], P_{\it{\Gamma}}^1 T_{\it{\Gamma}}^1 v\right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u], P_{\it{\Gamma}}^2 T_{\it{\Gamma}}^2 v\right)_{L^2({\it{\Gamma}})} \\ & \quad + \left( -[\partial_\nu {\it{\Delta}} u], \left(\operatorname{id}_{L^2({\it{\Gamma}})}-P_{\it{\Gamma}}^1\right)T_{\it{\Gamma}}^1 v\right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u], \left(\operatorname{id}_{L^2({\it{\Gamma}})}-P_{\it{\Gamma}}^2\right)T_{\it{\Gamma}}^2 v\right)_{L^2({\it{\Gamma}})}\text{.} \end{align*} Since $$\widetilde{v}_i := \left(\operatorname{id}_{L^2({\it{\Gamma}})} - P_{{\it{\Gamma}}}^i\right) T_{\it{\Gamma}}^i v \in H^{\frac{5}{2}-i}({\it{\Gamma}})$$ and because $${\it{\Gamma}}$$ is sufficiently smooth, there exists a function $$w \in H_0^2({\it{\Omega}})$$ such that $$T_{\it{\Gamma}}^i w = \widetilde{v}_i$$Lions & Magenes, 1972, Theorem 9.4). In particular, we have $$P_{\it{\Gamma}}^i T_{\it{\Gamma}}^i w = 0,$$ which yields $$w \in V_0$$ and   \begin{align*} 0 & = \frac{1}{\kappa} a(u,w) = \left( -[\partial_\nu {\it{\Delta}} u], T_{\it{\Gamma}}^1 w\right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u], T_{\it{\Gamma}}^2 w\right)_{L^2({\it{\Gamma}})} \\ & = \left( -[\partial_\nu {\it{\Delta}} u], (\operatorname{id}_{L^2({\it{\Gamma}})}-P_{\it{\Gamma}}^1)T_{\it{\Gamma}}^1 v\right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u], (\operatorname{id}_{L^2({\it{\Gamma}})}-P_{\it{\Gamma}}^2)T_{\it{\Gamma}}^2 v\right)_{L^2({\it{\Gamma}})} \end{align*} by applying the variational equation (5.2) for $$u$$. Combined with the Cauchy–Schwarz inequality, this implies   \begin{align*} a(u,v) \leq \kappa \left(\Vert [{\it{\Delta}} u] \Vert_{L^2({\it{\Gamma}})} + \Vert [\partial_\nu {\it{\Delta}} u]\Vert_{L^2({\it{\Gamma}})} \right) \left( \Vert P_{\it{\Gamma}}^1 T_{\it{\Gamma}}^1 v\Vert_{L^2({\it{\Gamma}})} + \Vert P_{\it{\Gamma}}^2 T_{\it{\Gamma}}^2 v\Vert_{L^2({\it{\Gamma}})} \right)\text{,} \end{align*} which was to be shown. □ Now we are in the situation to show the main results, namely the discretization error estimate for the discretizations (5.3) and (5.5). Theorem 5.2 Let $$u^h_\varepsilon \in S_h$$ be the solution of the discrete soft curve problem (5.4) and assume that $$\varepsilon_1 = c_1 h^{\lambda_1}$$ and $$\varepsilon_2 = c_2 h^{\lambda_2}$$. For any $$\delta >0$$ define   \begin{align*} \gamma = \min_{i\in\{1,2\}}\left( \frac{1}{2}-\delta, 3-i-2\delta-\frac{\lambda_i}{2}, \frac{\lambda_i}{2} \right)\text{.} \end{align*} Then there exists a constant $$c>0$$ independent of $$h$$ such that   \begin{align*} \left\Vert u - u^h_\varepsilon \right\Vert_{H^2({\it{\Omega}})} \leq c h^{\gamma}\text{.} \end{align*} In particular, $$\left\Vert u - u^h_\varepsilon \right\Vert_{H^2({\it{\Omega}})} \in {\mathcal O}(h^{\frac{1}{2}-\delta})$$ for $$\lambda_1 \in [1-2\delta,3-2\delta]$$ and $$\lambda_2 = 1-2\delta$$. Proof. In this proof, we use the notation $$A \lesssim B$$ whenever $$A \leq c B$$ holds with a constant $$c$$ that is independent of $$\varepsilon_1$$, $$\varepsilon_2$$ and $$h$$. As of Lemma 5.1, we can apply Theorem 3.3. Using this together with the coercivity of $$a(\cdot,\cdot)$$ on $$H_0^2({\it{\Omega}}),$$ we conclude for all $$v\in S_h$$,   \begin{align*} \left\Vert u - u^h_\varepsilon \right\Vert_{H^2({\it{\Omega}})} & \lesssim \left\Vert u - u^h_\varepsilon \right\Vert_{a} \\ & \lesssim \left\Vert u - v \right\Vert_a + \sum_{i=1}^2 \frac{1}{\sqrt{\varepsilon_i}} \Vert u - v \Vert_{b_i} + \sum_{i=1}^2 \sqrt{\varepsilon_i}. \end{align*} Note that $$a(\cdot,\cdot)$$ is continuous on $$H^2({\it{\Omega}})$$ and $$b_i(\cdot,\cdot)$$ is continuous on $$H^{i-\frac{1}{2}+\delta}({\it{\Omega}})$$. Furthermore, $$u \in H^{\frac{5}{2}-\delta}({\it{\Omega}})$$ according to Corollary 4.4. Choosing the interpolant $$v = \overline{u}$$ and applying the interpolation estimate (5.1) yields   \begin{align*} \left\Vert u - u^h_\varepsilon \right\Vert_{H^2({\it{\Omega}})} & \lesssim \Vert u - \overline{u} \Vert_{H^2({\it{\Omega}})} + \sum_{i=1}^2 \frac{1}{\sqrt{\varepsilon_i}} \Vert u - \overline{u} \Vert_{H^{i-\frac{1}{2}+\delta}({\it{\Omega}})} + \sum_{i=1}^2 \sqrt{\varepsilon_i} \\ & \lesssim h^{\frac{1}{2}-\delta} \Vert u \Vert_{H^{\frac{5}{2}-\delta}({\it{\Omega}})} + \sum_{i=1}^2 \frac{1}{\sqrt{\varepsilon_i}} h^{3-i-2\delta} \Vert u \Vert_{H^{\frac{5}{2}-\delta}({\it{\Omega}})} + \sum_{i=1}^2 \sqrt{\varepsilon_i} \\ & \lesssim h^{\frac{1}{2}-\delta} + \sum_{i=1}^2 h^{3-i-2\delta-\frac{\lambda_i}{2}} + \sum_{i=1}^2 h^{\frac{\lambda_i}{2}} \\ & \lesssim h^\gamma\text{.} \end{align*} This proves the assertion. □ We can proceed similarly to prove convergence rates for the soft bulk formulation. Lemma 5.3 Suppose that there is a constant $$c_0 >0$$ such that for every $$h \in I$$, $$K \in [0,2]$$ and $$k\in [0,K]$$ the inverse estimate   \begin{align} \forall\,{v \in S_h}, \ \vert v \vert_{H^K(B)} \leq c_0 h^{k-K} \vert v \vert_{H^k(B)} \end{align} (5.7) is fulfilled. Let $$u$$ be the solution of Problem 2.2 and $$\delta>0$$. Then there exists a constant $$c>0$$ such that for all $$h \in I$$ and $$v \in S_h$$,   \begin{align*} a(u,v) \leq c h^{\min\left(0,s-\frac{3}{2}-\delta\right)} \Vert v \Vert_{b_B}\text{.} \end{align*} Proof. From the proof of Lemma 5.1, we know for all $$v \in S_h$$,   \begin{align*} \frac{1}{\kappa} a(u,v) & = \left( -[\partial_\nu {\it{\Delta}} u], v \right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u],\partial_\nu v\right)_{L^2({\it{\Gamma}})} \\ & = \left( -[\partial_\nu {\it{\Delta}} u], T_{\it{\Gamma}}^1 P_B T_Bv \right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u], T_{\it{\Gamma}}^2 P_BT_Bv\right)_{L^2({\it{\Gamma}})} \\ & \quad + \left( -[\partial_\nu {\it{\Delta}} u], T_{\it{\Gamma}}^1 (\operatorname{id}_B - P_B) T_B v \right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u], T_{\it{\Gamma}}^2 (\operatorname{id}_B - P_B) T_Bv\right)_{L^2({\it{\Gamma}})}\text{.} \end{align*} Also as in the proof of Lemma 5.1, there exists a function $$w \in V_0 \subseteq H_0^2({\it{\Omega}})$$ such that $$T_{\it{\Gamma}}^i w = T_{\it{\Gamma}}^i\left(\operatorname{id}_B - P_B\right) T_B v$$ for $$i\in\{1,2\}$$. From this we conclude   \begin{align*} 0 = \frac{1}{\kappa} a(u,w) = \left( -[\partial_\nu {\it{\Delta}} u], T_{\it{\Gamma}}^1 (\operatorname{id}_B - P_B) T_B v \right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u], T_{\it{\Gamma}}^2 (\operatorname{id}_B - P_B) T_Bv\right)_{L^2({\it{\Gamma}})} \end{align*} as in the proof of Lemma 5.1. Furthermore, note that the operators $$T_{\it{\Gamma}}^i$$ are continuous and linear over $$H^{\frac{3}{2}+\delta}(B)$$ and thus the corresponding operator norms $${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {T_{{\it{\Gamma}}}^i} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_\delta$$ are bounded. It follows that   \begin{align*} \frac{1}{\kappa} a(u,v) & = \left( -[\partial_\nu {\it{\Delta}} u], T_{\it{\Gamma}}^1 P_B T_Bv \right)_{L^2({\it{\Gamma}})} +\left ( [{\it{\Delta}} u], T_{\it{\Gamma}}^2 P_BT_Bv\right)_{L^2({\it{\Gamma}})} \\ & \leq \Vert [\partial_\nu {\it{\Delta}} u] \Vert_{L^2({\it{\Gamma}})}{\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {T_{\it{\Gamma}}^1} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert} _{\delta} \Vert P_B T_B v \Vert_{H^{\frac{3}{2}+\delta}(B)} \\ & \quad + \Vert [{\it{\Delta}} u] \Vert_{L^2({\it{\Gamma}})} {\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {T_{\it{\Gamma}}^2}\right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{\delta} \Vert P_B T_B v \Vert_{H^{\frac{3}{2}+\delta}(B)}\text{.} \end{align*} Together with (5.7) this implies   \begin{align*} a(u,v) & \leq c h^{\min(0,s-\frac{3}{2}-\delta)}\Vert P_BT_Bv \Vert_{H^s(B),} \end{align*} which proves the assertion. □ Note that (5.7) is indeed an assumption. We cannot apply the standard finite element inverse estimates here as the grid is in general not matched to the domain $$B$$. A closer look to the proof of the standard inverse estimate reveals that the constant in (5.7) would go to infinity as $$\left\vert{E_h\cap B}\right\vert / \left\vert{E_h}\right\vert \to 0$$ for $$I \ni h \to 0$$ and some element $$E_h$$ in the grid associated to $$S_h$$. This means that the above inverse estimate over $$B$$ is only valid for grids that resolve $$B$$ sufficiently well. Theorem 5.4 Let the assumptions from Lemma 5.3 hold, $$u^h_\varepsilon \in S_h$$ the solution of the discrete soft bulk problem (5.6) and $$\varepsilon = c_0 h^\lambda$$ for some $$\delta > 0$$. Define   \begin{align*} \gamma = \min\left( \frac{1}{2}-\delta, \frac{5}{2}-\delta-s-\frac{\lambda}{2}, \min\left(0,s-\frac{3}{2}-\delta\right)+\frac{\lambda}{2}\right)\text{.} \end{align*} Then there exists a constant $$c>0$$, independent of $$h$$ such that   \begin{align*} \left\Vert u - u^h_\varepsilon \right\Vert_{H^2({\it{\Omega}})} & \leq c h^\gamma\text{.} \end{align*} In particular, $$\left\Vert u-u^h_\varepsilon\right\Vert_{H^2({\it{\Omega}})} \in \mathcal{O}\left(h^{\frac{1}{2}-\delta}\right)$$ for $$s \leq \frac{3}{2}+\delta$$ and $$\lambda = 4-2s$$. Proof. As of Lemma 5.3, we can apply Theorem 3.3. Noting the continuity of $$b_B(\cdot,\cdot)$$ on $$H^s({\it{\Omega}}),$$ we obtain as in the proof of Theorem 5.2 that   \begin{align*} \left\Vert u - u^h_\varepsilon \right\Vert_{H^2({\it{\Omega}})} & \lesssim \Vert u - v \Vert_{H^2({\it{\Omega}})}({\it{\Omega}}) + \frac{1}{\sqrt{\varepsilon}} \Vert u - v\Vert_{H^s({\it{\Omega}})} + h^{\min\left(0,s-\frac{3}{2}-\delta\right)}\sqrt{\varepsilon} \text{.} \end{align*} Again choosing the interpolant $$v = \overline{u}$$ and applying the interpolation estimate (5.1) finally gives   \begin{align*} \left\Vert u - u^h_\varepsilon \right\Vert_{H^2({\it{\Omega}})} & \lesssim h^{\frac{1}{2}-\delta} + h^{\frac{5}{2}-\delta-s-\frac{\lambda}{2}} + h^{\min\left(0,s-\frac{3}{2}-\delta\right)+\frac{\lambda}{2}} \lesssim h^\gamma \end{align*} and thus proves the assertion. □ 6. Numerical example 6.1. A symmetric example problem with known exact solution We select a problem with a known, analytically computable solution. We consider the circular domain $$\hat{{\it{\Omega}}} = \{x \in \mathbb{R}^2 \mathrel\vert |x| < r_2\}$$ with one embedded particle $$\hat{B} = \{x \in \mathbb{R}^2 \mathrel\vert |x| < r_1\}$$ for $$0 < r_1 < r_2 < 1$$. On the boundary of $$\hat{{\it{\Omega}}}$$ we consider homogeneous Dirichlet boundary conditions, whereas on the particle boundary $${\it{\Gamma}}_1 = \partial B_1$$, we impose variable-height coupling conditions of the form (2.1), i.e.,   \begin{align} u|_{\partial B_1} &= f^1_1 + \gamma^1, & \partial_\nu u|_{\partial B_1} &= f^1_2 \end{align} (6.1) for variable $$\gamma^1 \in \mathbb{R}$$ with the given functions   \begin{align*} f_1^1(re^{{\rm i}\theta}) &= \cos(n \theta), & f_2^1(x) &= 0. \end{align*} We select the parameters   \begin{align*} r_2 &= \tfrac23, & r_1 &= \tfrac13, & \kappa &= 1, & \sigma &= 0, & n &=4. \end{align*} Then the exact solution of Problem 2.1 on $$\hat{{\it{\Omega}}}$$ is the fourfold symmetric function $$u \in H^2(\hat{{\it{\Omega}}})$$ given by   \begin{align*} u(re^{{\rm i}\theta}) = \begin{cases} \cos(4\theta) \, (c_1 r^4 + c_2 r^6) & \text{ for}\, r \in (0,r_1], \\ \cos(4\theta) \, ( c_3 r^{-2} + c_4 r^{-4} + c_5 r^4 + c_6 r^6 ) & \text{ for}\, r \in [r_1,r_2], \\ 0 & \text{ otherwise}, \end{cases} \end{align*} with the constants   \begin{align*} c_1 & = 243, \qquad & c_2 & = -1458, \qquad & c_3 & =-6909/689, \\ c_4 & = -7936/502281, \qquad & c_5 & = 11502/689, \qquad & c_6 & = 44288/167427 \text{.} \end{align*} To discretize the circular domain $$\hat{{\it{\Omega}}}$$ with the presented Bogner–Fox–Schmit discretization, we embed $$\hat{{\it{\Omega}}}$$ into the larger rectangular domain $${\it{\Omega}} = [-1,1]^2$$ and treat the condition on $${\it{\Gamma}}_2 = \partial \hat{{\it{\Omega}}}$$ like a second particle boundary condition. It is easily seen that the solution of Problem 2.1 on $${\it{\Omega}}$$ is obtained by extending the solution given on $$\hat{{\it{\Omega}}}$$ by zero. The solution of this problem is depicted in Fig. 2. Fig. 2. View largeDownload slide Exact solution of the symmetric example problem. Left: top view. Right: rendered three-dimensional view. Fig. 2. View largeDownload slide Exact solution of the symmetric example problem. Left: top view. Right: rendered three-dimensional view. 6.2. Discretization, quadrature and error measurement For the finite element discretization $$S_h \subseteq H_0^2({\it{\Omega}}),$$ we divide $${\it{\Omega}}$$ into uniform squares with edges of length $$h$$ and equip the resulting grid with a Bogner–Fox–Schmit finite element basis. The boundary conditions over $$\partial{\it{\Omega}}$$ are enforced by setting the corresponding degrees of freedom to zero while the boundary conditions over the $${\it{\Gamma}}_i$$ are replaced by penalized constraints. The resulting soft curve problem with penalized curve constraints is then given by (5.4) while the soft bulk problem with penalized bulk constraints is given by (5.6), where the area of the virtual second particle is now given by $$B_2 = {\it{\Omega}} \setminus \hat{{\it{\Omega}}}$$. We note that we drop the projection for the constraints on $${\it{\Gamma}}_2$$ and $$B_2$$ because we do not allow variable height there. While our analysis is formulated only in terms of the variable-height constraints incorporated using projections, this is in fact no limitation, since all arguments directly carry over to the case without these projections. Since the assembly of the bilinear forms $$b^i_{\it{\Gamma}}(\cdot,\cdot)$$ and $$b_B(\cdot,\cdot)$$ involves integrals over the curves $${\it{\Gamma}}_i$$ or the nontrivial domains $$B_i$$, respectively, we briefly mention how these can be approximated. Integrals over full grid cells are evaluated exactly using standard tensor Gaussian quadrature rules. For quadrature over the $${\it{\Gamma}}_i$$, we use a trigonometric Gaussian quadrature as described in Da Fies & Vianello (2012) which in our case is both exact for integration of finite element functions and for integration of $$u$$. Integrals over the nonrectangular domains $$B_i$$ are approximated using the local parameterization quadrature method introduced in Olshanskii & Safin (2016), which is no longer exact. The discretization error is measured in terms of the norm $$\| {\it{\Delta}} (u^h_\varepsilon - u)\|_{L^2({\it{\Omega}})}$$ on $$H_0^2({\it{\Omega}})$$. It is important to note that the error cannot be computed accurately by simple quadrature over the grid elements, because $$u$$ is not smooth across the $${\it{\Gamma}}_i$$. Instead we split the integration into the subdomains $$B_1$$, $$B_2 = {\it{\Omega}} \setminus \hat{{\it{\Omega}}}$$ and $${\it{\Omega}} \setminus (B_1 \cup B_2),$$ where we apply the quadrature rules described above. 6.3. Results for the soft curve formulation In case of the soft curve formulation, we expect convergence with order $${\mathcal O}(h^{\frac12})$$ due to Theorem 5.2 if the penalty parameters are chosen suitably. For the numerical examples, we selected penalty parameters $$\varepsilon = ( c h^{\lambda_1}, c h )$$ with $$\lambda_1 \in [1,3]$$ according to Theorem 5.2. The constant was fixed to be $$c=10^{-3}$$. Figure 3 shows the behavior of the $$H^2$$-discretization error over the mesh size $$h$$ for uniform grids and penalty parameters with $$\lambda_1 \in \{1, 2, 3\}$$. In accordance with Theorem 5.2, we observe convergence of at least order $${\mathcal O}(h^{1/2})$$. While the observed rate seems to behave like $${\mathcal O}(h^{1/2})$$ for $$\lambda_1 \in \{2,3\}$$ it seems that the rate is slightly better for $$\lambda_1 = 1$$. A possible explanation for this observation is that the error contribution for the term $$\|\cdot\|_{b_1}$$ in the proof of Theorem 5.2 can be improved by smaller values of $$\lambda_1$$. However, this does not lead to a better theoretical bound due to other dominating terms. Fig. 3. View largeDownload slide $$H^2$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 3. View largeDownload slide $$H^2$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Figures 4 and 5 show the $$H^1$$- and $$L^2$$-discretization errors, respectively, for the same set of example problems. In both cases, we observe convergence with order $${\mathcal O}(h)$$. While improved convergence order for weaker norms is a well-known property of many discretizations, we emphasize that this is not covered by the theory presented here and may be considered in the future with a refined analysis. Fig. 4. View largeDownload slide $$H^1$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 4. View largeDownload slide $$H^1$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 5. View largeDownload slide $$L^2$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 5. View largeDownload slide $$L^2$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. 6.4. Results for the soft bulk formulation Because we use uniform grids and do not refine with regard to the geometry of the $${\it{\Gamma}}_i$$, we lose control over the inverse estimate assumption (5.7). Thus, we cannot prove condition (3.4) in general. Assuming that (3.4) still holds true, we could use Theorem 5.4 in the case of exact integration. However, as discussed above, we approximate integrals using only the quadrature from Olshanskii & Safin (2016). If this is accurate enough, we can in view of the Strang-type result in Proposition 3.4 expect convergence of order $${\mathcal O}(h^{1/2})$$ when using $$\varepsilon = c h^{4-2s}$$ and the penalty norm $$\|\cdot\|_{H^s(B_i)}$$ for $$s\in [0,1]$$. For our numerical experiments we considered $$c = 10^{-3}$$ and $$s \in \{0,1\}$$. Figure 6 shows the behavior of the discretization error over the mesh size $$h$$ for uniform grids and $$s\in \{0, 1\}$$. Again the observed order is in accordance with the expected order $${\mathcal O}(h^{1/2})$$. However, especially for $$s=1,$$ this is perturbed by strong oscillations. This is maybe due to the fact that the constant resulting from (3.4) is not uniformly bounded and can strongly vary depending on how the curve $${\it{\Gamma}}$$ intersects mesh elements. Again the $$H^1$$- and $$L^2$$-errors depicted in Figs 7 and 8, respectively, exhibit similar behavior but the improved convergence order $${\mathcal O}(h)$$. Fig. 6. View largeDownload slide $$H^2$$-errors for the soft bulk formulation for $$s=0$$ and $$s=1$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 6. View largeDownload slide $$H^2$$-errors for the soft bulk formulation for $$s=0$$ and $$s=1$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 7. View largeDownload slide $$H^1$$-errors for the soft bulk formulation for $$s=0$$ and $$s=1$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 7. View largeDownload slide $$H^1$$-errors for the soft bulk formulation for $$s=0$$ and $$s=1$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 8. View largeDownload slide $$L^2$$-errors for the soft bulk formulation for $$s=0$$ and $$s=1$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 8. View largeDownload slide $$L^2$$-errors for the soft bulk formulation for $$s=0$$ and $$s=1$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. 6.5. A nonsymmetric example problem While the first example allowed us to compute the exact solution due to its symmetry properties, we will now assess the discretization for a nonsymmetric example with several particles but without known exact solution. We consider the domain $${\it{\Omega}} = [-1,1]^2$$ with homogeneous Dirichlet boundary conditions on $$\partial {\it{\Omega}}$$ and we again select the parameters   \begin{align*} \kappa &= 1, & \sigma &= 0. \end{align*} This time we consider four elliptical particles embedded nonsymmetrically into the membrane domain $${\it{\Omega}}$$. The ellipses’ major and minor axes have lengths $$0.4$$ and $$0.2$$, respectively. Their positions and orientations are depicted in Fig. 9. Fig. 9. View largeDownload slide Solution of the nonsymmetric example problem. Left: top view. Right: rendered three-dimensional view. Fig. 9. View largeDownload slide Solution of the nonsymmetric example problem. Left: top view. Right: rendered three-dimensional view. Regarding the particle–membrane coupling conditions, we assume that all particles have a constant height profile $$f^i_1 = 0$$ and a constant slope $$f^i_2$$ varying for different particles. More precisely, we use   \begin{align*} f^1_2 &= 1, & f^2_2 &= 1, & f^3_2 &= -1, & f^4_2 &= -1. \end{align*} The solution resulting from that setup is depicted in Fig. 9. Since we do not know the exact solution for this problem, we estimate the error by comparing our discrete solutions to a reference solution computed on a finer grid. We computed the solutions for the different problem formulations on uniform grids with mesh size $$h=2^{-k}$$ for $$k= 3,\dots,6,$$ while the reference solution was computed on a uniform grid with mesh size $$h=2^{-8}$$ using the soft curve formulation with $$\varepsilon = 10^{-3}(h^3,h)$$. The (approximate) discretization errors for the soft curve and soft bulk formulations are shown in Figs 10 and 11, respectively. As before, we used $$\varepsilon = (ch^\lambda, ch)$$ with $$\lambda \in \{1,2,3\}$$ for the soft curve formulation and $$\varepsilon = ch^{4-2s}$$ with $$s \in \{0,1\}$$ for the soft bulk formulation. In both cases, we selected $$c = 10^{-3}$$. Fig. 10. View largeDownload slide $$H^2$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$. Fig. 10. View largeDownload slide $$H^2$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$. Fig. 11. View largeDownload slide $$H^2$$-errors for the soft bulk formulation for $$s=0$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,4,5$$ and $$s=1$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$ (from left to right). Fig. 11. View largeDownload slide $$H^2$$-errors for the soft bulk formulation for $$s=0$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,4,5$$ and $$s=1$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$ (from left to right). For the soft bulk formulation, the projection $$P_{B}$$ may lead to a densely populated matrix for those degrees of freedom located inside the particles. To avoid the resulting computational effort for $$s=1,$$ we used the $$H^1(B_i)$$-norm in the form   \begin{align*} \|v\|_{H^1(B_i)}^2 = \|\nabla v\|_{L^2(B_i)}^2 + \Bigl(\int_{B_i} v \, \,{\rm d}x\Bigr)^2. \end{align*} Then the penalty term that incorporates the $$H^1(B_i)$$-projections takes the form   \begin{align*} b_{B}(w,v) = (\nabla w, \nabla v)_{L^2(B)}, \end{align*} which results in a sparse matrix again. For both formulations, soft curve and soft bulk, we observe a rate that is slightly better than the expected order of convergence, namely $${\mathcal O}\left(h^{1/2}\right)$$. Note that we restricted the mesh size $$h=2^{-k}$$ to powers of 2 in this example to simplify comparison with the reference solution. Furthermore, in the case of soft bulk constraints with $$s=0$$, we were not able to carry out the computations for $$h=2^{-6}$$ due to hardware restrictions. This explains the much smaller number of data points in the plots. As before, we also computed the $$H^1$$- and $$L^2$$-errors for the same set of example problems. For the soft curve formulation, we again observe that the order of the $$H^1$$- and $$L^2$$-errors depicted in Figs 12 and 13, respectively, is essentially squared in comparison with the $$H^2$$-error leading to a convergence order which is approximately $${\mathcal O}(h)$$. The situation is different for the soft bulk formulation. Here the $$H^1$$-error depicted in Fig. 14 is of order $${\mathcal O}(h^{3/2})$$ and the $$L^2$$-error depicted in Fig. 15 is of order $${\mathcal O}(h^{5/2})$$. While the improved $$H^1$$- and $$L^2$$-order is not covered by the presented theory anyway, we also cannot explain the surprising difference in the observed order for the symmetric and nonsymmetric example problems. Fig. 12. View largeDownload slide $$H^1$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$. Fig. 12. View largeDownload slide $$H^1$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$. Fig. 13. View largeDownload slide $$L^2$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$. Fig. 13. View largeDownload slide $$L^2$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$. Fig. 14. View largeDownload slide $$H^1$$-errors for the soft bulk formulation for $$s=0$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,4,5$$ and $$s=1$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$ (from left to right). Fig. 14. View largeDownload slide $$H^1$$-errors for the soft bulk formulation for $$s=0$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,4,5$$ and $$s=1$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$ (from left to right). Fig. 15. View largeDownload slide $$L^2$$-errors for the soft bulk formulation for $$s=0$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,4,5$$ and $$s=1$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$ (from left to right). Fig. 15. View largeDownload slide $$L^2$$-errors for the soft bulk formulation for $$s=0$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,4,5$$ and $$s=1$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$ (from left to right). Funding This research has been partially supported by Deutsche Forschungsgemeinschaft (DFG) through grant CRC 1114 ‘Scaling Cascades in Complex Systems’, Project AP01 ‘Particles in lipid bilayers’. Appendix Proof of Lemma 4.3. Let $$u_1 := u|_{{\it{\Omega}}_B}$$ and $$u_2 := u|_B$$ as well as $${\it{\Omega}}_1 := {\it{\Omega}}_B$$ and $${\it{\Omega}}_2 := B$$. Define $$\varphi_{\alpha}(x,y) := {(\vert D^\alpha u(x) - D^\alpha u(y) \vert^2)}/{(\Vert x-y\Vert^{2s+2})}$$. For $$s \in (0,\frac{1}{2}),$$ the Sobolev–Slobodeckij norm of $$u$$ is given by   \begin{align*} \Vert u\Vert_{H^{2+s}({\it{\Omega}})}^2 & = \Vert u\Vert_{H^2({\it{\Omega}})}^2 + \max_{\vert\alpha\vert = 2} \Vert D^\alpha u\Vert_{H^s({\it{\Omega}})}^2 \\ & = \Vert u\Vert_{H^2({\it{\Omega}})}^2 + \max_{\vert\alpha\vert = 2} \int_{{\it{\Omega}} \times {\it{\Omega}}} \varphi_{\alpha}(x,y) \, \mathrm{d}(x,y)\text{.} \end{align*} Splitting the domain of integration yields   \begin{align*} \int_{{\it{\Omega}} \times {\it{\Omega}}} \varphi_{\alpha}(x,y) \, \mathrm{d}(x,y) & = \left\Vert D^\alpha u_1 \right\Vert_{H^s({\it{\Omega}}_1)}^2 + \left\Vert D^\alpha u_2\right\Vert_{H^s({\it{\Omega}}_2)}^2 + 2\int_{{\it{\Omega}}_1 \times {\it{\Omega}}_2} \varphi_{\alpha}(x,y) \, \mathrm{d}(x,y)\text{.} \end{align*} Using the layer-cake principle, we write   \begin{align*} \int_{{\it{\Omega}}_1 \times {\it{\Omega}}_2} \varphi_{\alpha}(x,y) = \int_0^\infty \left \vert \left\{ (x,y) \in {\it{\Omega}}_1\times {\it{\Omega}}_2 \mid \varphi_{\alpha}(x,y) \geq t \right\} \right\vert \, \mathrm{d}t\text{.} \end{align*} By the assumptions from Section 2 it holds that $${\it{\Omega}}_1$$ and $${\it{\Omega}}_2$$ are Lipschitz domains. By virtue of the Sobolev embedding theorems, we conclude $$u_i \in C^2(\overline{{\it{\Omega}}_i})$$ (see, e.g., Adams, 1975, Theorem 5.4). Consequently, the constant $$C := \max_{\vert\alpha\vert = 2} \sup_{{\it{\Omega}}_1 \times {\it{\Omega}}_2} \vert D^\alpha u_1(x) - D^\alpha u_2(y)\vert^2 \in \mathbb{R}_{>0}$$ exists and thus the implication   \begin{align*} \frac{\vert D^\alpha u(x) - D^\alpha u(y)\vert^2}{\Vert x-y\Vert^{2s+2}} \geq t \ \Longrightarrow \ \Vert x-y\Vert \leq C^{\frac{1}{2s+2}} t^{-\frac{1}{2s+2}} =: \widetilde{C} t^{-\frac{1}{2s+2}} \end{align*} holds for all $$(x,y) \in {\it{\Omega}}_1 \times {\it{\Omega}}_2$$. Let $${\it{\Gamma}} := \partial{\it{\Omega}}_1 \cap \partial {\it{\Omega}}_2$$. From   \begin{align*} \left\vert \left\{ (x,y) \in {\it{\Omega}}_1\times {\it{\Omega}}_2 \mid \Vert x-y\Vert \leq \widetilde{C} t^{-\frac{1}{2s+2}} \right\} \right\vert \leq 4\vert{\it{\Gamma}}\vert (\widetilde{C} t^{-\frac{1}{2s+2}})(\widetilde{C} t^{-\frac{1}{2s+2}})^2 =: \hat{C} t^{\frac{-3}{2s+2}} \end{align*} we infer   \begin{align*} & \int_0^\infty \left \vert \left\{ (x,y) \in {\it{\Omega}}_1\times {\it{\Omega}}_2 \mid \varphi_{\alpha}(x,y) \geq t \right\} \right\vert \, \mathrm{d}t \\ & \leq \vert {\it{\Omega}} \vert^2 + \int_1^\infty \left \vert \left\{ (x,y) \in {\it{\Omega}}_1\times {\it{\Omega}}_2 \mid \varphi_{\alpha}(x,y) \geq t \right\} \right\vert \, \mathrm{d}t \\ & \leq \vert {\it{\Omega}} \vert^2 + \int_1^\infty \hat{C} t^{\frac{-3}{2s+2}} \, \mathrm{d}t\text{.} \end{align*} This expression is finite for $$s < \frac{1}{2}$$. This implies $$u \in H^{2+\frac{1}{2}-\delta}({\it{\Omega}})$$ for all $$\delta >0$$ as was to be shown. □ Proof of Proposition 3.4. For brevity we define $$\widetilde{u} := \widetilde{u}^X_\varepsilon$$ and assume without loss of generality that $$\varepsilon = (\varepsilon_i)_{i=1,\dots,m} \leq 1$$ holds componentwise. Let $$v \in X$$ and define $$\widetilde{w}:= v -\widetilde{u}$$. We have by (3.11),   \begin{align*} \left\Vert{\widetilde{w}}\right\Vert_{a_\varepsilon}^2 & = \left\Vert{\widetilde{w}}\right\Vert_a^2 + \sum_{i=1}^m \frac{1}{\varepsilon_i} \left\Vert{\widetilde{w}}\right\Vert_{b_i}^2 \\ & \leq \left\Vert{\widetilde{w}}\right\Vert_a^2 + \sum_{i=1}^m \frac{1}{\varepsilon_i} \left\Vert{\widetilde{w}}\right\Vert_{\widetilde{b}_i}^2 + \sum_{i=1}^m \frac{c_i}{\varepsilon_i} \left\Vert{\widetilde{w}}\right\Vert^2\text{.} \end{align*} Using the continuity of $$a$$, the coercivity of $$\widetilde{a}_\varepsilon$$ and $$\varepsilon \leq 1,$$ this leads to the inequality   \begin{align} \left\Vert{\widetilde{w}}\right\Vert_{a_\varepsilon}^2 & \leq \left(\frac{\left\Vert{a}\right\Vert}{\widetilde{\alpha}} + 1 + \sum_{i=1}^m \frac{c_i}{\varepsilon_i \widetilde{\alpha}} \right) \left\Vert{\widetilde{w}}\right\Vert_{\widetilde{a}_\varepsilon}^2\text{.} \end{align} (A.1) On the other hand, we have the equality   \begin{align*} \left\Vert{\widetilde{w}}\right\Vert_{\widetilde{a}_\varepsilon}^2 & = \widetilde{a}_\varepsilon(v,\widetilde{w}) - \widetilde{a}_\varepsilon(\widetilde{u},\widetilde{w}) + \ell_\varepsilon(\widetilde{w}) - a_\varepsilon(u,\widetilde{w}) + a_\varepsilon(v,\widetilde{w}) - a_\varepsilon(v,\widetilde{w}), \end{align*} which by application of (3.10), the Cauchy–Schwarz inequality, coercivity of $$a_\varepsilon$$, definition of $$\widetilde{a}_\varepsilon$$ and $$\widetilde{\ell}_\varepsilon$$ and $$\varepsilon \leq 1$$ yields   \begin{align} \begin{aligned} \left\Vert{\widetilde{w}}\right\Vert_{\widetilde{a}_\varepsilon}^2 & = a_\varepsilon(v-u,\widetilde{w}) + (\widetilde{a}_\varepsilon-a_\varepsilon)(v,\widetilde{w}) + (\ell_\varepsilon-\widetilde{\ell}_\varepsilon)(\widetilde{w}) \\ & \leq \left\Vert{v-u}\right\Vert_{a_\varepsilon} \left\Vert{\widetilde{w}}\right\Vert_{a_\varepsilon} \\ & \quad + \frac{\left\vert{(\widetilde{a}-a)(v,\widetilde{w})}\right\vert + \sum_{i=1}^m \frac{1}{\varepsilon_i} \left\vert{(\widetilde{b}_i-b_i)(u-v,\widetilde{w})}\right\vert + \left\vert{(\ell-\widetilde{\ell})(\widetilde{w})}\right\vert}{\sqrt{\alpha} \left\Vert{\widetilde{w}}\right\Vert} \left\Vert{\widetilde{w}}\right\Vert_{a_\varepsilon}. \end{aligned} \end{align} (A.2) Then from the triangle inequality and (A.2) inserted into (A.1), we get   \begin{align*} \left\Vert{u - \widetilde{u}}\right\Vert_{a_\varepsilon} & \leq \left\Vert{u - v}\right\Vert_{a_\varepsilon} + \left\Vert{v - \widetilde{u}}\right\Vert_{a_\varepsilon} \\ & \leq \left\Vert{u - v}\right\Vert_{a_\varepsilon} + \left(\frac{\left\Vert{a}\right\Vert}{\widetilde{\alpha}} + 1 + \sum_{i=1}^m \frac{c_i}{\varepsilon_i \widetilde{\alpha}} \right) \cdot \left( \left\Vert{v-u}\right\Vert_{a_\varepsilon} + \phantom{\frac{\sum_{i=1}^m \frac{1}{\varepsilon_i}(\widetilde{b}}{\left\Vert{\widetilde{w}}\right\Vert}}\right. \\ & \qquad + \left. \frac{\left\vert{(\widetilde{a}-a)(v,\widetilde{w})}\right\vert + \sum_{i=1}^m \frac{1}{\varepsilon_i} \left\vert{(\widetilde{b}_i-b_i)(u-v,\widetilde{w})}\right\vert + \left\vert{(\ell-\widetilde{\ell})(\widetilde{w})}\right\vert}{\sqrt{\alpha} \left\Vert{\widetilde{w}}\right\Vert} \right), \end{align*} which proves the statement by taking the infimum over all $$v \in X$$ and replacing $$\widetilde{w} = v - \widetilde{u}$$ by the supremum over all $$w \in X$$. □ References Adams R. ( 1975) Sobolev Spaces . Pure and Applied Mathematics. New York-London: Academic Press. Babuška I. ( 1973) The finite element method with penalty. Math. Comp. , 27, 221– 228. Google Scholar CrossRef Search ADS   Bahrami A. H., Raatz M., Agudo-Canalejo J., Michel R., Curtis E. M., Hall C. K., Gradzielski M., Lipowsky R. & Weikl T. R. ( 2014) Wrapping of nanoparticles by membranes. Adv. Colloid Interface Sci. , 208, 214– 224. Special issue in honour of Wolfgang Helfrich. Google Scholar CrossRef Search ADS PubMed  Barrett J. W. & Elliott C. M. ( 1986) Finite element approximation of the Dirichlet problem using the boundary penalty method. Numer. Math. , 49, 343– 366. Google Scholar CrossRef Search ADS   Blum H. & Rannacher R. ( 1980) On the Boundary Value Problem of the Biharmonic Operator on Domains with Angular Corners . Preprint: Sonderforschungsbereich Approximation und Mathematische Optimierung in einer Anwendungsbezogenen Mathematik. Mathematical Methods in the Applied Sciences, John Wiley & Sons, Ltd., SFB 72. Ciarlet P. G. ( 1978) The Finite Element Method for Elliptic Problems . Amsterdam-New York-Oxford: North-Holland Publishing Company. Da Fies G. & Vianello M. ( 2012) Trigonometric Gaussian quadrature on subintervals of the period. Electron. Trans. Numer. Anal. , 39, 102– 112. Dommersnes P. G. & Fournier J.-B. ( 1999) Casimir and mean-field interactions between membrane inclusions subject to external torques. Europhys. Lett. , 46, 256. Google Scholar CrossRef Search ADS   Dommersnes P. G. & Fournier J.-B. ( 2002) The many-body problem for anisotropic membrane inclusions and the self-assembly of ‘saddle’ defects into an ‘egg carton’. Biophys. J. , 83, 2898– 2905. Google Scholar CrossRef Search ADS PubMed  Elliott C. M., Gräser C., Hobbs G., Kornhuber R. & Wolf M.-W. ( 2016) A variational approach to particles in lipid membranes. Arch. Ration. Mech. Anal. , 222, 1011– 1075. Google Scholar CrossRef Search ADS   Grisvard P. ( 1985) Elliptic Problems in Nonsmooth Domains . Monographs and Studies in Mathematics. Boston, London, Melbourne: Massachusetts, USA-London: Pitman Publishing. Helfrich P. & Jakobsson E. ( 1990) Calculation of deformation energies and conformations in lipid membranes containing gramicidin channels. Biophys. J. , 57, 1075– 1084. Google Scholar CrossRef Search ADS PubMed  Lions J. & Magenes E. ( 1972) Non-homogeneous Boundary Value Problems and Applications . Massachusetts, USA-London: Springer. Lunardi A. ( 2009) Interpolation Theory  ( Publications of the Scuola Normale Superiore). Scuola Normale Superiore, Pisa PI, Italy: Edizioni della Normale. Nitsche J. ( 1971) Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind. Abh. Math. Semin. Univ. Hambg. , 36, 9– 15. Google Scholar CrossRef Search ADS   Olshanskii M. A. & Safin D. ( 2016) Numerical integration over implicitly defined domains for higher order unfitted finite element methods. Lobachevskii J. Math. , 37, 582– 596. Google Scholar CrossRef Search ADS   Weikl T. R., Kozlov M. M. & Helfrich W. ( 1998) Interaction of conical membrane inclusions: effect of lateral tension. Phys. Rev. E , 57, 6988. Google Scholar CrossRef Search ADS   © The authors 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

# Discretization error estimates for penalty formulations of a linearized Canham–Helfrich-type energy

, Volume Advance Article – Jan 4, 2018
24 pages

/lp/ou_press/discretization-error-estimates-for-penalty-formulations-of-a-J0cBYHdXfu
Publisher
Oxford University Press
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/drx071
Publisher site
See Article on Publisher Site

### Abstract

Abstract This article is concerned with minimization of a fourth-order linearized Canham–Helfrich energy subject to Dirichlet boundary conditions on curves inside the domain. Such problems arise in the modeling of the mechanical interaction of biomembranes with embedded particles. There, the curve conditions result from the imposed particle–membrane coupling. We prove almost-$$H^{\frac{5}{2}}$$ regularity of the solution and then consider two possible penalty formulations. For the combination of these penalty formulations with a Bogner–Fox–Schmit finite element discretization, we prove discretization error estimates that are optimal in view of the solution’s reduced regularity. The error estimates are based on a general estimate for linear penalty problems in Hilbert spaces. Finally, we illustrate the theoretical results by numerical computations. An important feature of the presented discretization is that it does not require the particle boundary to be resolved. This is crucial to avoid re-meshing if the presented problem arises as a subproblem in a model where particles are allowed to move or rotate. 1. Introduction A standard model for the behavior of biomembranes on a macroscale is the Canham–Helfrich model, which describes a biological membrane mathematically as a hypersurface $$\mathcal{M}\subseteq \mathbb{R}^3$$ minimizing the Canham–Helfrich energy   \begin{align*} \mathcal{J}_{\rm CHS}(\mathcal{M}) = \int_{\mathcal{M}} \frac{1}{2} \kappa H^2 + \kappa_{\rm G} K + \sigma \, \mathrm{d}\mathcal{H}^2\text{.} \end{align*} Here $$H$$ and $$K$$ are the mean and Gaussian curvature of $$\mathcal{M}$$, the coefficients $$\kappa>0$$ and $$\kappa_{\rm G}\geq 0$$ are the corresponding bending rigidities and $$\sigma \geq 0$$ is the membrane’s surface tension. Under the assumption that the membrane is ‘almost flat’, one can justify a geometric linearization of this functional for a membrane patch, leading to the so-called Monge–gauge approximation   \begin{align*} \mathcal{J}_{{\it{\Omega}}}(u) = \int_{\it{\Omega}} \frac{1}{2} \kappa ({\it{\Delta}} u)^2 + \frac{1}{2} \sigma \vert\nabla u\vert^2 \,\mathrm{{\rm d}}x \end{align*} of the Canham–Helfrich energy. Here the membrane patch is considered to be the graph $$\mathcal{M} = \{(x,u(x)) \mathrel\vert x \in {\it{\Omega}} \}$$ of a function $$u\colon {\it{\Omega}} \to \mathbb{R}$$ over some reference domain $${\it{\Omega}}\subseteq \mathbb{R}^2$$. A variety of hybrid models for the coupling of embedded particles to the membrane have been considered. These hybrid models are based on a continuous surface description of the membrane while particles are described by discrete entities (see, e.g., Helfrich & Jakobsson, 1990; Weikl et al., 1998; Dommersnes & Fournier, 1999, 2002; Bahrami et al., 2014). For an overview on hybrid models, we refer to Elliott et al. (2016) and the references cited therein. In this article we consider coupling conditions imposed on the particle boundaries and follow the notation introduced in Elliott et al. (2016). There, the reference domain $${\it{\Omega}}$$ is split into the membrane’s domain $${\it{\Omega}}_B \subseteq {\it{\Omega}}$$ and the particles’ domain $$B = {\it{\Omega}} \setminus {\it{\Omega}}_B$$. The model introduces membrane–particle interactions by functions that prescribe the membrane’s height profile and slope on the interface $${\it{\Gamma}} := \partial {\it{\Omega}}_B \cap \partial B$$. Altogether this yields an energy minimization problem subject to Dirichlet boundary conditions on $${\it{\Gamma}}$$. This article introduces and analyses a discretization based on a penalty formulation for the corresponding boundary value problem that avoids the resolution of particle boundaries. The reason for considering penalized boundary conditions is the following: since a biological membrane behaves like a fluid in tangential directions, particles can in principle move and rotate in plane. Any model simulating moving particles and any method computing optimal particle positions will thus have to solve multiple problems with varying particle positions. If the particle boundaries had to be resolved, this would require mesh deformation or re-meshing which can be computationally quite expensive. An alternative is to replace the strict boundary conditions by adequate penalty terms that can be formulated without having to resolve the boundary. One such penalty approach was introduced in Elliott et al. (2016) and is called soft curve formulation. A novel second formulation in this article is the soft bulk formulation. Both penalty problems will be defined later on. Our goal is to derive discretization error estimates for those penalty problems. For the second-order equations such estimates are well known (see, e.g., Nitsche, 1971; Babuška, 1973; Barrett & Elliott, 1986). It turns out that the fourth-order problems that we are interested in can be treated sufficiently well using a simple general framework for penalty problems in Hilbert spaces. This is mostly due to the fact that the regularity of our solutions is rather limited in the first place and so we can use rather simple estimates to still obtain optimal rates of convergence. Thanks to the abstract formulations, we get as a by-product a general error theory for finite element penalty problems with low regularity. The article is structured as follows. In Section 2, we introduce the notation and problem formulations that are used throughout this article. Section 3 is devoted to an abstract error result for linear penalty approximations on Hilbert spaces. As a foundation for the application of this result to the problem at hand, we then discuss regularity of solutions in Section 4. There we prove that a solution of the original problem lies in $$H^{\frac{5}{2}-\delta}({\it{\Omega}})$$ for all $$\delta>0$$. In Section 5, we combine the regularity with the developed abstract results to show the optimal (in view of the restricted regularity) convergence rate $${\mathcal O}(h^{1/2-\delta})$$ for a discretization with Bogner–Fox–Schmit finite elements. Finally, Section 6 illustrates our results by numerical examples that reproduce our theoretical findings. 2. Notation and problem formulations We consider a membrane with $$k$$ embedded particles. To this end, let $${\it{\Omega}} \subset \mathbb{R}^2$$ be a bounded reference domain with Lipschitz boundary $$\partial {\it{\Omega}}$$ and $$B_i \subset {\it{\Omega}}$$, $$i=1,\dots, k$$ the area occupied by the $$i$$th particle. For simplicity, we assume that the $$B_i$$ are closed, nonempty, connected and pairwise disjoint. Furthermore we require that each particle $$B_i$$ has a $$C^{1,1}$$-boundary. By $$B = \bigcup_{i=1}^k B_i$$ and $${\it{\Omega}}_B = {\it{\Omega}} \setminus B,$$ we denote the areas occupied by the particles and the membrane, respectively. Since the $$B_i$$ are closed and disjoint, the total membrane–particle interface is given by   \begin{align*} {\it{\Gamma}} = \bigcup_{i=1}^k \partial B_i = \partial {\it{\Omega}}_B \setminus \partial {\it{\Omega}} = \partial {\it{\Omega}}_B \cap B. \end{align*} For simplicity, we consider the space   \begin{align*} H = \{v \in H^2({\it{\Omega}}_B) \mathrel\vert v|_{\partial {\it{\Omega}}} = \partial_\nu v|_{\partial{\it{\Omega}}} = 0 \} \end{align*} with homogeneous Dirichlet boundary conditions on the boundary of $$\partial {\it{\Omega}}$$. Notice that we can treat Navier and periodic boundary conditions on $$\partial {\it{\Omega}}$$ using the same techniques (cf. Elliott et al., 2016). We assume that the membrane–particle interaction is governed by boundary conditions on the interface. More precisely, each particle $$B_i$$ enforces a height profile given by $$f^i_1\colon \partial B_i \to \mathbb{R}$$ and a slope given by $$f^i_2\colon \partial B_i \to \mathbb{R}$$ to the membrane on its boundary $$\partial B_i$$. That is, we consider boundary values   \begin{align} u|_{\partial B_i} &= f^i_1 + \gamma^i, & \partial_\nu u|_{\partial B_i} &= f^i_2, \end{align} (2.1) where $$\nu$$ is the unit outward normal to $$\partial {\it{\Omega}}_B$$, and the parameter $$\gamma^i \in \mathbb{R}$$ is allowed to vary freely to factor out the average height on $$\partial B_i$$. This is necessary because we want to prescribe the height profile only, while the absolute or average height is not fixed. For the following, we collect all such boundary data in a function $$f = (f_1,f_2) = \sum_{i=1}^{k} f^i\colon {\it{\Gamma}} \to \mathbb{R}^2$$, where $$f^i = (f^i_1,f^i_2)\colon \partial B_i \to \mathbb{R}^2$$ is extended by zero to the whole of $${\it{\Gamma}}$$. Using this notation, we consider the minimization problem: Problem 2.1 Find $$u\in H$$ and $$\gamma \in \mathbb{R}^k$$ minimizing $$\mathcal{J}_{{\it{\Omega}}_B}(u)$$ subject to   \begin{align} u|_{{\it{\Gamma}}} = f_1 + \sum_{i=1}^k \gamma_{i} \eta_{i}, \qquad \partial_\nu u|_{{\it{\Gamma}}} = f_2\text{ .} \end{align} (2.2) Notice that the effect of the free parameter $$\gamma_i$$ is localized to $$\partial B_i$$ via the use of the indicator function $$\eta_i := \chi_{\partial B_i} \in L^2({\it{\Gamma}}),$$ which is 1 on $$\partial B_i$$ but vanishes on all $$\partial B_j$$, $$j \neq i$$. In view of the trace theorem, we from now on assume that $$f = (f_1,f_2) \in H^{\frac{3}{2}}({\it{\Gamma}}) \times H^{\frac{1}{2}}({\it{\Gamma}})$$. Under this assumption it is known that Problem 2.1 admits a unique solution by application of Lax–Milgram’s theorem (Elliott et al., 2016). It is in fact possible to simplify this problem formulation. For this purpose, we make use of the trace operator $$T_{{\it{\Gamma}}} = (T_{{\it{\Gamma}}}^1, T_{{\it{\Gamma}}}^2)$$,   \begin{align*} T_{{\it{\Gamma}}}\colon H^2({\it{\Omega}}) & \rightarrow H^{\frac{3}{2}}({\it{\Gamma}}) \times H^{\frac{1}{2}}({\it{\Gamma}}), & \qquad T_{{\it{\Gamma}}}(v) = (v|_{{\it{\Gamma}}}, \ \partial_\nu v|_{{\it{\Gamma}}} ). \end{align*} This operator is well defined, continuous, surjective and admits a continuous right inverse (see e.g. Grisvard, 1985, Theorem 1.5.1.2). Note that we can also view $$T_{{\it{\Gamma}}}$$ as a trace operator for $$H^2({\it{\Omega}}_B)$$ due to our regularity assumptions on $${\it{\Gamma}}$$. To simplify the notation for boundary values up to the average height in (2.2), we furthermore introduce the projection operator   \begin{align*} P_{{\it{\Gamma}}} \colon H^{\frac{3}{2}}({\it{\Gamma}}) \times H^{\frac{1}{2}}({\it{\Gamma}}) &\longrightarrow H^{\frac{3}{2}}({\it{\Gamma}}) \times H^{\frac{1}{2}}({\it{\Gamma}}) \end{align*} defined by $$P_{\it{\Gamma}}(v_1,v_2) = \left(P_{\it{\Gamma}}^1(v_1),P_{\it{\Gamma}}^2(v_2)\right),$$ where   \begin{align*} P_{\it{\Gamma}}^1(v) = v - \sum_{i=1}^k |\partial B_i|^{-1}\int_{\partial B_i} v \,\mathrm{d}\sigma = v - \sum_{i=1}^k \frac{(v, \eta_i)_{L^2({\it{\Gamma}})}}{(\eta_i,\eta_i)_{L^2({\it{\Gamma}})}} \eta_i, \qquad P_{\it{\Gamma}}^2(v) = v. \end{align*} Using this operator, we can formulate the boundary conditions (2.2) as $$P_{\it{\Gamma}} (T_{\it{\Gamma}} u - f ) = 0,$$ which gives rise to the space of admissible functions defined by   \begin{align*} V_f := \left\{v \in H_0^2({\it{\Omega}}) \mid P_{\it{\Gamma}} (T_{\it{\Gamma}} v - f ) = 0 \right \} \end{align*} and the corresponding minimization problem: Problem 2.2 Find $$u\in V_f$$ minimizing $$\mathcal{J}_{{\it{\Omega}}}(u)$$. Notice that Problem 2.2 differs from Problem 2.1 not only due to the different notation for the boundary conditions on $${\it{\Gamma}}$$ but also because it minimizes $$J_{\it{\Omega}}$$ for functions defined on the whole of $${\it{\Omega}}$$ and not just on $${\it{\Omega}}_B$$. However, the following result from Elliott et al. (2016) shows that a solution of Problem 2.2 also immediately yields a solution to Problem 2.1. Proposition 2.3 Let $$u \in H_0^2({\it{\Omega}})$$ be the solution of Problem 2.2. Define $$\gamma_i = \vert \partial B_i\vert^{-1} (u-f_1, \eta_i )_{L^2({\it{\Gamma}})}$$ for $$i = 1,\dots,k$$. Then $$(u|_{{\it{\Omega}}_B},\gamma)$$ is the solution of Problem 2.1. While Problem 2.2 allows for variable height of particles, the full problem considered in Elliott et al. (2016) has additional degrees of freedom. This is due to the fact that particles can move and rotate in the plane of the fluid membrane. To avoid mesh deformation or re-meshing whenever particle positions change, we will in the following drop the hard constraints at the particle boundaries in favor of a penalized approach. Replacing the hard curve constraints by penalty terms in the energy functional leads to the following problem. Problem 2.4 Find $$u_\varepsilon \in H^2_0({\it{\Omega}})$$ minimizing   \begin{align*} \mathcal{J}_{{\it{\Omega}}}(u_\varepsilon ) + \sum_{i=1}^2 \frac{1}{\varepsilon_i} \left\Vert P_{\it{\Gamma}}^i ( T_{{\it{\Gamma}}}^i u_\varepsilon - f_i ) \right\Vert_{L^2({\it{\Gamma}})}^2\text{.} \end{align*} This formulation is more favorable than Problem 2.2 insofar as it admits a straightforward conforming finite element discretization without re-meshing in the case of variable particle positions. At this point, we want to mention an alternative penalty formulation that is based on the idea that, for known shapes $$B_i$$, the solution $$u|_B$$ could be computed a priori up to a constant per component $$B_i$$. We define the restriction operator   \begin{align*} T_B\colon H^2({\it{\Omega}}) & \longrightarrow H^{2}(B), & \qquad & T_B(v)= v|_B. \end{align*} Analogously to the curve constraint formulation, we introduce the associated projection operator   \begin{align*} P_B\colon H^2(B) & \longrightarrow H^2(B), & P_B(v) &= v - \sum_{i=1}^k \frac{(v,\psi_i)_{H^s(B)}}{(\psi_i,\psi_i)_{H^s(B)}} \psi_i, \end{align*} for $$\psi_i := \chi_{B_i}$$ and some fixed $$s \in [0,2]$$. Using this notation, the bulk constrained problem reads as follows. Problem 2.5 Find $$u\in H_0^2({\it{\Omega}})$$ minimizing $$\mathcal{J}_{{\it{\Omega}}}(u)$$ subject to   \begin{align*} P_B T_B u = P_B u|_B \text{.} \end{align*} One quickly verifies that the solutions of Problems 2.2 and 2.5 coincide. Analogously to Problem 2.4, a penalty formulation of Problem 2.5 is given by the following. Problem 2.6 Find $$u_\varepsilon \in H_0^2({\it{\Omega}})$$ minimizing   \begin{align*} \mathcal{J}_{{\it{\Omega}}}(u_\varepsilon) + \frac{1}{\varepsilon} \Vert P_B(T_Bu_\varepsilon - u|_B) \Vert_{H^s(B)}^2\text{.} \end{align*} While the penalized formulations are more flexible in the sense that the domain $${\it{\Omega}}_B$$ does not have to be resolved for discretization, one will be faced with the problem of balancing the penalty parameter and discretization errors. To this end, we will first develop an abstract error estimate for penalized problems and then analyse the regularity of the solution of the hard constrained Problem 2.2. Using a suitable regularity result later allows us to derive optimal $$h$$-dependent values of $$\varepsilon_i$$ for a discretization with mesh size $$h$$. 3. An error estimate for linear penalty problems on Hilbert spaces In this section, we derive an abstract energy error estimate for penalized discretizations of linearly constrained problems in a Hilbert space setting. Let $$H$$ be a Hilbert space and $$U_0 \subset H$$ a closed subspace. We consider the affine closed subspace $$U = U_0 + u_0 \subseteq H$$ of $$H$$ for some given $$u_0 \in H$$. In addition, let $$a\colon H \times H \rightarrow \mathbb{R}$$ be a symmetric positive semidefinite bounded bilinear form and $$\ell \colon H \rightarrow \mathbb{R}$$ a bounded linear form on $$H$$. We furthermore assume that $$a(\cdot,\cdot)$$ is coercive on $$U_0$$. In this setting, we consider the affine constrained minimization problem: Problem 3.1 Find $$u \in U$$ minimizing   \begin{align*} \tfrac{1}{2} a(u,u) - \ell(u)\text{.} \end{align*} Application of Lax–Milgram’s theorem yields the existence of a unique solution $$u \in U$$ of Problem 3.1 characterized by the variational equation   \begin{align} a(u, v) = \ell(v) \qquad \forall\,{v \in U_0}. \end{align} (3.1) From now on, $$u \in U$$ denotes this unique solution. Now let us assume that we want to approximate this problem by a penalty formulation over some closed linear subspace $$X \subseteq H$$ with $$m\in\mathbb{N}$$ penalty terms. Those penalty terms shall be given by symmetric positive semidefinite bounded bilinear forms $$b_i\colon H \times H \rightarrow \mathbb{R}$$ and penalty parameters $$\varepsilon_i >0$$. Denoting by $$\Vert v \Vert_{c} := \sqrt{ c(v,v) }$$ the seminorm induced by a symmetric positive semidefinite bilinear form $$c(\cdot,\cdot)$$ on $$H$$ the penalized problem reads as follows. Problem 3.2 Find $$u^X_\varepsilon \in X$$ minimizing   \begin{align*} \frac{1}{2} a\left(u^X_\varepsilon,u^X_\varepsilon\right) - \ell\left(u^X_\varepsilon\right) + \sum_{i=1}^m \frac{1}{2\varepsilon_i} \left\Vert u^X_\varepsilon-u\right\Vert_{b_i}^2. \end{align*} For the sake of convenience, we write   \begin{align} a_\varepsilon(\cdot,\cdot) = a(\cdot,\cdot) + \sum_{i=1}^m \frac{1}{\varepsilon_i} b_i(\cdot,\cdot), \qquad \ell_\varepsilon(\cdot) = \ell(\cdot) + \sum_{i=1}^m \frac{1}{\varepsilon_i} b_i(u,\cdot) \end{align} (3.2) for $$\varepsilon = (\varepsilon_i)_{i=1,\dots,m} \in \mathbb{R}_{+}^m = (0,\infty)^m$$. We require that $$a_1(\cdot,\cdot)$$ is coercive on $$X$$, such that $$a_\varepsilon(\cdot,\cdot)$$ is also coercive for all $$\varepsilon \in (0,1]^m$$. Under these assumptions, Lax–Milgram’s theorem implies existence of a unique solution $$u^X_\varepsilon \in X$$ for any $$\varepsilon \in (0,1]^m$$ characterized by the variational equation   \begin{align} a_\varepsilon\left(u^X_\varepsilon, v\right) = \ell_\varepsilon(v) \qquad \forall\, v \in X. \end{align} (3.3) The following result states a Céa-type estimate for the error $$u-u^X_\varepsilon$$ resulting from penalization and discretization in $$X$$. Theorem 3.3 For the fixed $$u \in U$$ used in the definition of $$\ell_\varepsilon$$, suppose there exist constants $$c_i >0$$ for $$i=1,\dots,m$$ such that   \begin{align} \vert a(u, v) - \ell(v) \vert \leq \sum_{i=1}^m c_i \Vert v\Vert_{b_i} \qquad \forall\,{v\in X}. \end{align} (3.4) Then the error for the solution $$u^X_\varepsilon$$ of Problem 3.2 can be bounded by   \begin{align} \Vert u- u^X_\varepsilon \Vert_{a_\varepsilon}^2 \leq 3 \inf_{v\in X} \left( \Vert u-v \Vert_{a_\varepsilon}^2 + \sum_{i=1}^m \varepsilon_i c_i^2 \right)\!\text{.} \end{align} (3.5) Proof. Let $$e = u - u^X_\varepsilon$$. Then we get   \begin{align} \Vert e \Vert_{a_\varepsilon}^2 = a_\varepsilon(e,e) = a_\varepsilon(e,u-v) + a_\varepsilon\left(e,v-u^X_\varepsilon\right) \end{align} (3.6) for all $$v\in X$$. Using Young’s inequality, we can bound the first term by   \begin{align} \begin{aligned} a_\varepsilon(e,u-v) &\leq \Vert e\Vert_{a_\varepsilon} \, \Vert u - v\Vert_{a_\varepsilon} \leq \tfrac{1}{4} \Vert e\Vert_{a_\varepsilon}^2 + \Vert u - v\Vert_{a_\varepsilon}^2. \end{aligned} \end{align} (3.7) The definition (3.2) of the $$u$$-dependent penalized functional $$\ell_\varepsilon$$ yields   \begin{align*} a_\varepsilon(u,w) - a(u,w) & = \ell_\varepsilon(w) - \ell(w) \qquad \forall\,{w \in X}. \end{align*} Combining this identity with (3.3) we get   \begin{align} a_\varepsilon(e,w) = a(u,w) - \ell(w) \qquad \forall\,{w \in X}. \end{align} (3.8) Because of $$v - u^X_\varepsilon \in X,$$ this implies   \begin{align} \begin{aligned} a_\varepsilon\left(e,v-u^X_\varepsilon\right) & = a\left(u,v-u^X_\varepsilon\right) - \ell\left(v-u^X_\varepsilon\right) \\ & \leq \sum_{i=1}^m c_i \left\Vert v - u^X_\varepsilon\right\Vert_{b_i} \\ & \leq \sum_{i=1}^m c_i \left( \Vert e\Vert_{b_i} + \Vert u - v \Vert_{b_i} \right) \\ & \leq \sum_{i=1}^m \left( \frac{1}{4\varepsilon_i} \Vert e\Vert_{b_i}^2 + \frac{1}{2\varepsilon_i} \Vert u - v\Vert_{b_i}^2 + \frac{3}{2}\varepsilon_i c_i^2\right) \\ & \leq \frac{1}{4}\Vert e\Vert_{a_\varepsilon}^2 + \frac{1}{2}\Vert u-v\Vert_{a_\varepsilon}^2 + \sum_{i=1}^m \frac{3}{2}\varepsilon_i c_i^2, \end{aligned} \end{align} (3.9) where the second to last inequality follows from Young’s inequality. Inserting the estimates (3.7) and (3.9) into (3.6), we obtain   \begin{align*} \Vert e\Vert_{a_\varepsilon}^2 \leq \frac{1}{2}\Vert e\Vert_{a_\varepsilon}^2 + \frac{3}{2}\Vert u-v\Vert_{a_\varepsilon}^2 + \sum_{i=1}^m \frac{3}{2}\varepsilon_i c_i^2, \end{align*} which proves the assertion. □ Next we consider the case where the evaluation of $$a(\cdot,\cdot)$$, $$\ell(\cdot)$$ and $$b_i(\cdot,\cdot)$$ is not performed exactly but approximated by some $$\widetilde{a}$$, $$\widetilde{\ell}$$ and $$\widetilde{b}_i$$, respectively. We use a notation analogous to $$a_\varepsilon$$ and $$\ell_\varepsilon$$:   \begin{align*} \widetilde{a}_\varepsilon(\cdot,\cdot) = \widetilde{a}(\cdot,\cdot) + \sum_{i=1}^m \frac{1}{\varepsilon_i} \widetilde{b}_i(\cdot,\cdot) , \qquad \widetilde{\ell}_\varepsilon(\cdot) = \widetilde{\ell}(\cdot) + \sum_{i=1}^m \frac{1}{\varepsilon_i} \widetilde{b}_i(u,\cdot). \end{align*} Instead of solving (3.3) for $$u^X_\varepsilon$$ directly, one now computes a $$\widetilde{u}^X_\varepsilon$$ by solving   \begin{align} \widetilde{a}_\varepsilon(\widetilde{u}^X_\varepsilon, v ) = \widetilde{\ell}_\varepsilon(v) \qquad \forall\,{v \in X} \text{.} \end{align} (3.10) In this setting, we can prove the following Strang-type result: Proposition 3.4 Let $$\left\Vert{a}\right\Vert$$ be the continuity constant of $$a$$ with respect to the $$H$$-norm and let $$a_1 = a+\sum_{i=1}^m b_i$$ and $$\widetilde{a}_1 = \widetilde{a} +\sum_{i=1}^m \widetilde{b}_i$$ be coercive with respect to the constants $$\alpha$$ and $$\widetilde{\alpha}$$, respectively. In addition to the assumptions of Theorem 3.3, suppose that for all $$i\in\{1,\dots,m\}$$ there exists $$c_i \in \mathbb{R}$$ such that for all $$v \in X$$,   \begin{align} \left\vert{ b_i(v,v) - \widetilde{b}_i(v,v) }\right\vert \leq c_i \Vert v \Vert^2\text{.} \end{align} (3.11) Then there exists a constant $$C = C(\alpha,\widetilde{\alpha},\left\Vert{a}\right\Vert) >0$$ such that   $$\matrix{ {{{\left\| {u - {\rm{ }}\tilde u_\varepsilon ^X} \right\|}_{{a_\varepsilon }}} \le C\left( {1 + \sum\limits_{i = 1}^m {{{{c_i}} \over {{\varepsilon _i}}}} } \right)\mathop {\inf }\limits_{v \in X} \mathop {\sup }\limits_{w \in X} {\rm{ }}\left[ {{{\left\| {u - v} \right\|}_{{a_\varepsilon }}} + {1 \over {\left\| w \right\|}}\left| {(a - \tilde a)(v,w)} \right|} \right.} \hfill \cr \hspace{20mm} {\quad \left. { + {1 \over {\left\| w \right\|}}\sum\limits_{i = 1}^m {{1 \over {{\varepsilon _i}}}} \left| {({b_i} - {{\tilde b}_i})(u - v,w)} \right| + {1 \over {\left\| w \right\|}}\left| {(\ell - \tilde \ell )(w)} \right|} \right].} \hfill \cr }$$ Proof. See the Appendix. □ 4. Regularity of the hard curve constraint problem To apply the results of the previous section, it will be crucial to prove (3.4). For this purpose and to prove convergence rates by bounding the best approximation errors in Theorem 3.3, we will now derive regularity results for the solution of Problem 2.2. Let $$u \in V_f$$ be the solution of the hard curve minimization formulation, Problem 2.2. Knowing the regularity of $$u$$ is central to proving discretization errors since the maximal regularity of the solution immediately reveals the optimal rates of convergence that one would expect for the corresponding discretization errors. Our strategy is to rewrite the hard curve minimization problem as a system of elliptic partial differential equations to which we apply standard regularity theory for elliptic interface problems. To this end, we define the fourth-order elliptic differential operator $$L := \kappa {\it{\Delta}}^2 - \sigma {\it{\Delta}}$$ associated with the energy functional $$\mathcal{J}$$. Proposition 4.1 Suppose $$(u_1,u_2) \in H^2({\it{\Omega}}_B) \times H^2(B)$$ is a weak solution of the system of partial differential equations   \begin{align} \begin{aligned} Lu_1 & = 0 \text{ on}\,{\it{\Omega}}_B, & T_{\it{\Gamma}} u_1 & = T_{\it{\Gamma}} u \text{ on}\, {\it{\Gamma}}, & u_1 = \partial_\nu u_1 = 0 \text{ on}\, \partial{\it{\Omega}}, \\ Lu_2 & = 0 \text{ on}\, B, & T_{\it{\Gamma}} u_2 & = T_{\it{\Gamma}} u \text{ on}\,{\it{\Gamma}}. \end{aligned} \end{align} (4.1) Then $$u_1 = u|_{{\it{\Omega}}_B}$$ and $$u_2 = u|_B$$ where $$u$$ is the solution of Problem 2.2. Conversely, if $$u$$ is the solution of Problem 2.2, then $$(u|_{{\it{\Omega}}_B}, u|_B)$$ solves (4.1). Proof. A proof for the equivalence of Problem 2.2 with a weak formulation of (4.1) is given in Elliott et al. (2016, Proposition 4.1). □ Lemma 4.2 Let $${\it{\Omega}}$$ be a piecewise polygonal domain, let $${\it{\Gamma}}$$ be smooth and suppose $$(f_1,f_2) \in H^{\frac{7}{2}}({\it{\Gamma}})\times H^{\frac{5}{2}}({\it{\Gamma}})$$. If all corners of $${\it{\Omega}}$$ have an inner angle $$\omega$$ with $$\omega \leq 126^\circ$$ then $$(u|_{{\it{\Omega}}_B},u|_B) \in H^{4}({\it{\Omega}}_B)\times H^{4}(B)$$. Proof Let $$g := \frac{\sigma}{\kappa} {\it{\Delta}} u|_{{\it{\Omega}}_B} \in L^2({\it{\Omega}}_B)$$. From Proposition 4.1, we know that $$u|_{{\it{\Omega}}_B}$$ is a weak solution of   \begin{align*} {\it{\Delta}}^2 u = g \text{ on}\, {\it{\Omega}}_B, \qquad T_{{\it{\Gamma}}}u = T_{{\it{\Gamma}}} u|_{{\it{\Omega}}_B} \text{ on}\, {\it{\Gamma}}, \qquad u = \partial_\nu u = 0 \text{ on}\, \partial{\it{\Omega}}. \end{align*} Since $$g \in L^2({\it{\Omega}}_B)$$ and $$T_{{\it{\Gamma}}}^1 u|_{{\it{\Omega}}_B} = f_1 + \sum_{i=1}^k \vert \partial B_i\vert^{-1} ( u|_{{\it{\Omega}}_B} - f_1, \eta_i)_{L^2({\it{\Gamma}})}\eta_i \in H^{\frac{7}{2}}({\it{\Gamma}})$$ as well as $$T_{{\it{\Gamma}}}^2 u|_{{\it{\Omega}}_B} = f_2 \in H^{\frac{5}{2}}({\it{\Gamma}}),$$ it follows from Grisvard (1985, Theorem 7.2.2.3) and the computations in Blum & Rannacher (1980) on the associated characteristic equation that $$u|_{{\it{\Omega}}_B} \in H^4({\it{\Omega}}_B)$$. An analogous argument on $$B$$ yields $$u|_B \in H^4(B)$$ and thus proves the assertion. □ Having the intrinsic regularity $$u \in H^2({\it{\Omega}})$$ and the piecewise regularity from Lemma 4.2 allows us to show an improved global regularity result. The key ingredient is the following technical lemma. Lemma 4.3 Let $$v \in H^2({\it{\Omega}})$$ such that $$v|_{{\it{\Omega}}_B} \in H^4({\it{\Omega}}_B)$$ and $$v|_{B} \in H^4(B)$$. Then $$v \in H^{2+\frac{1}{2}-\delta}({\it{\Omega}})$$ for all $$\delta > 0$$. Proof. See the Appendix. □ Combining Lemmas 4.2 and 4.3 now gives a corollary. Corollary 4.4 Let the assumptions from Lemma 4.2 hold. Then $$u \in H^{\frac{5}{2}-\delta}({\it{\Omega}})$$ for every $$\delta >0$$. 5. Discretization errors In this section, we apply the results from the previous sections to derive discretization errors for penalized finite element approximations of the form of Problems 2.4 and 2.6 but with $$h$$-dependent penalty parameters. We suppose from now on that $${\it{\Omega}}$$ is a rectangular domain. In view of Lemma 4.2, this implies that $$u|_{{\it{\Omega}}_B} \in H^4({\it{\Omega}}_B)$$, $$u|_B \in H^4(B)$$ and $$u \in H^{\frac{5}{2}-\delta}({\it{\Omega}})$$ for any $$\delta >0$$. On the domain $${\it{\Omega}}$$ we establish a quadrilateral grid equipped with Bogner–Fox–Schmit finite elements (see, e.g., Ciarlet, 1978). Given the set of grid nodes $$N$$ and the set of multiindices $${\it{\Lambda}} = \{(0,0), (1,0), (0,1), (1,1)\}$$, the resulting finite element space is spanned by a basis $$(\psi_{p,\alpha})_{p\in N, \alpha \in {\it{\Lambda}}}$$ of piecewise bicubic polynomials such that each $$\psi_{p,\alpha}$$ satisfies   \begin{align*} \partial^\beta \psi_{p,\alpha}(q) = \delta_{\alpha,\beta} \delta_{p,q} \qquad \forall\, q \in N, \beta \in {\it{\Lambda}}, \end{align*} where $$\partial^\beta$$ is the partial derivative associated with the multiindex $$\beta,$$ and $$\delta_{a,b}$$ is the Kronecker delta. That is, the degrees of freedom are the values, first-order partial derivatives and mixed second-order partial derivatives at the vertices (see Fig. 1). Fig. 1. View largeDownload slide Degrees of freedom for the Bogner–Fox–Schmit element. Fig. 1. View largeDownload slide Degrees of freedom for the Bogner–Fox–Schmit element. By setting the degrees of freedom on $$\partial{\it{\Omega}}$$ to zero this yields a conforming subspace of $$C^1(\overline{{\it{\Omega}}}) \cap H_0^2({\it{\Omega}})$$. Given a family $$(S_h)_{h\in I}$$ of such finite element discretizations over $${\it{\Omega}}$$ based on quasi-uniform grids with mesh size $$h,$$ we will use the following well-known result: there exists a constant $$c>0$$ such that for all $$K \in [0,4],$$ the approximation estimate   \begin{align} \forall\,{v \in H^K({\it{\Omega}}) \cap H_0^2({\it{\Omega}})},\,\exists{\overline{v} \in S_h}, \, \forall\,{k\in [0,\min\{2,K\}]}, \Vert v - \overline{v} \Vert_{H^k({\it{\Omega}})} \leq c h^{K-k} \Vert v \Vert_{H^K({\it{\Omega}})} \end{align} (5.1) holds true for all $$h \in I$$. (This is a direct consequence of classical approximation estimates (see, e.g., Ciarlet, 1978, Theorem 3.1.4), and interpolation theory of function spaces, (see e.g., Lunardi, 2009, Theorem 1.1.6).) To put our minimization problems into the variational framework used in Section 3, we introduce the bilinear form   \begin{align*} a \colon H^2({\it{\Omega}}) \times H^2({\it{\Omega}}) &\longrightarrow \mathbb{R}, & a(w,v) &= \kappa ( {\it{\Delta}} w, {\it{\Delta}} v )_{L^2({\it{\Omega}})} + \sigma ( \nabla w, \nabla v )_{L^2({\it{\Omega}})}, \end{align*} together with the curve penalty terms   \begin{align*} b_{{\it{\Gamma}}}^1 \colon H^2({\it{\Omega}}) \times H^2({\it{\Omega}}) & \longrightarrow \mathbb{R}, & b_{{\it{\Gamma}}}^1(w,v) &= ( P_{\it{\Gamma}}^1T_{\it{\Gamma}}^1 w, P_{\it{\Gamma}}^1T_{\it{\Gamma}}^1 v )_{L^2({\it{\Gamma}})},\\ b_{{\it{\Gamma}}}^2 \colon H^2({\it{\Omega}}) \times H^2({\it{\Omega}}) & \longrightarrow \mathbb{R}, & b_{{\it{\Gamma}}}^2(w,v) &= ( P_{\it{\Gamma}}^2T_{\it{\Gamma}}^2 w, P_{\it{\Gamma}}^2T_{\it{\Gamma}}^2 v )_{L^2({\it{\Gamma}})} \end{align*} and the bulk penalty term   \begin{align*} b_{B}\colon H^2({\it{\Omega}}) \times H^2({\it{\Omega}}) & \longrightarrow \mathbb{R}, & b_{B}(w,v) &= ( P_BT_B w, P_BT_B v )_{H^s(B)}\text{.} \end{align*} Then the solution $$u \in V_f \subset H_0^2({\it{\Omega}})$$ of Problem 2.2 satisfies the variational equation   \begin{align} a(u,v) = 0 \qquad \forall\,{v \in V_0}. \end{align} (5.2) Discretization of the penalized curve constrained Problem 2.4 in $$S_h$$ leads to the soft curve problem: find $$u_\varepsilon^h \in S_h$$ minimizing   \begin{align} \mathcal{J}_{{\it{\Omega}}}(u_\varepsilon^h) + \sum_{i=1}^2 \frac{1}{\varepsilon_i} \left\Vert P_{\it{\Gamma}}^i (T_{{\it{\Gamma}}}^i u_\varepsilon^h - f_i)\right \Vert_{L^2({\it{\Gamma}})}^2\text{,} \end{align} (5.3) which is equivalently characterized by the variational equation   \begin{align} u_\varepsilon^h \in S_h: \qquad a\left(u^h_\varepsilon,v\right) + \sum_{i=1}^2 \frac{1}{\varepsilon_i} b_{\it{\Gamma}}^i\left(u^h_\varepsilon-u,v\right) = 0 \qquad \forall\,{v \in S_h}. \end{align} (5.4) Analogously, discretization of the penalized bulk constrained Problem 2.6 in $$S_h$$ leads to the soft bulk problem: find $$u_\varepsilon^h \in S_h$$ minimizing   \begin{align} \mathcal{J}_{{\it{\Omega}}}\left(u_\varepsilon^h\right) + \frac{1}{\varepsilon} \left\Vert P_BT_B(u_\varepsilon^h - u)\right \Vert_{H^s(B)}^2 \end{align} (5.5) with the equivalent variational equation   \begin{align} u_\varepsilon^h \in S_h: \qquad a\left(u^h_\varepsilon,v\right) + \frac{1}{\varepsilon} b_B\left(u^h_\varepsilon-u,v\right) = 0 \qquad \forall\,{v \in S_h}. \end{align} (5.6) At this point, we want to emphasize that the solutions of (5.4) and (5.6) are in general different. We nevertheless refer in a slight abuse of notation to both solutions as $$u^h_\varepsilon$$. In the following it will, however, be clear from the context whether $$u^h_\varepsilon$$ denotes the solution of the soft curve problem or the soft bulk problem. In the next result, we show that the soft curve formulation meets the requirement from Theorem 3.3, which we will use afterwards to get an asymptotic error estimate. Lemma 5.1 For all $$v \in H_0^2({\it{\Omega}})$$ the solution $$u$$ of Problem 2.2 satisfies   \begin{align*} a(u,v) \leq \kappa \left(\Vert [{\it{\Delta}} u] \Vert_{L^2({\it{\Gamma}})} + \Vert [\partial_\nu {\it{\Delta}} u]\Vert_{L^2({\it{\Gamma}})} \right) \left( \Vert v\Vert_{b^1_{\it{\Gamma}}} + \Vert v \Vert_{b^2_{\it{\Gamma}}} \right)\text{.} \end{align*} Here $$[w] = w|_{{\it{\Omega}}_B} - w|_B$$ denotes the jump of the function $$w$$ on $${\it{\Gamma}}$$. Proof. By integration by parts and as $$u \in H^2({\it{\Omega}})$$, $$u|_{{\it{\Omega}}_B} \in H^4({\it{\Omega}}_B)$$, $$u|_B \in H^4(B)$$ and $$Lu = 0$$, we are able to state for all $$v \in H_0^2({\it{\Omega}})$$,   \begin{align*} \frac{1}{\kappa} a(u,v) & = \int_{{\it{\Omega}}} {\it{\Delta}} u {\it{\Delta}} v + \frac{\sigma}{\kappa} \nabla u \cdot \nabla v \\ & = \int_{{\it{\Omega}}_B} \frac{1}{\kappa} Lu|_{{\it{\Omega}}_B} \, v + \int_B \frac{1}{\kappa} Lu|_B \, v + \int_{{\it{\Gamma}}} \frac{\sigma}{\kappa} \left( \partial_\nu u|_{{\it{\Omega}}_B} - \partial_\nu u|_B \right) v \\ & \quad + \int_{{\it{\Gamma}}} \left( - \partial_\nu {\it{\Delta}} u|_{{\it{\Omega}}_B} + \partial_\nu {\it{\Delta}} u|_B \right) v + \int_{{\it{\Gamma}}} \left( {\it{\Delta}} u|_{{\it{\Omega}}_B} - {\it{\Delta}} u|_B \right ) \partial_\nu v \\ & = \left( -[\partial_\nu {\it{\Delta}} u], T_{\it{\Gamma}}^1 v\right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u], T_{\it{\Gamma}}^2 v\right)_{L^2({\it{\Gamma}})} \\ & = \left( -[\partial_\nu {\it{\Delta}} u], P_{\it{\Gamma}}^1 T_{\it{\Gamma}}^1 v\right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u], P_{\it{\Gamma}}^2 T_{\it{\Gamma}}^2 v\right)_{L^2({\it{\Gamma}})} \\ & \quad + \left( -[\partial_\nu {\it{\Delta}} u], \left(\operatorname{id}_{L^2({\it{\Gamma}})}-P_{\it{\Gamma}}^1\right)T_{\it{\Gamma}}^1 v\right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u], \left(\operatorname{id}_{L^2({\it{\Gamma}})}-P_{\it{\Gamma}}^2\right)T_{\it{\Gamma}}^2 v\right)_{L^2({\it{\Gamma}})}\text{.} \end{align*} Since $$\widetilde{v}_i := \left(\operatorname{id}_{L^2({\it{\Gamma}})} - P_{{\it{\Gamma}}}^i\right) T_{\it{\Gamma}}^i v \in H^{\frac{5}{2}-i}({\it{\Gamma}})$$ and because $${\it{\Gamma}}$$ is sufficiently smooth, there exists a function $$w \in H_0^2({\it{\Omega}})$$ such that $$T_{\it{\Gamma}}^i w = \widetilde{v}_i$$Lions & Magenes, 1972, Theorem 9.4). In particular, we have $$P_{\it{\Gamma}}^i T_{\it{\Gamma}}^i w = 0,$$ which yields $$w \in V_0$$ and   \begin{align*} 0 & = \frac{1}{\kappa} a(u,w) = \left( -[\partial_\nu {\it{\Delta}} u], T_{\it{\Gamma}}^1 w\right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u], T_{\it{\Gamma}}^2 w\right)_{L^2({\it{\Gamma}})} \\ & = \left( -[\partial_\nu {\it{\Delta}} u], (\operatorname{id}_{L^2({\it{\Gamma}})}-P_{\it{\Gamma}}^1)T_{\it{\Gamma}}^1 v\right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u], (\operatorname{id}_{L^2({\it{\Gamma}})}-P_{\it{\Gamma}}^2)T_{\it{\Gamma}}^2 v\right)_{L^2({\it{\Gamma}})} \end{align*} by applying the variational equation (5.2) for $$u$$. Combined with the Cauchy–Schwarz inequality, this implies   \begin{align*} a(u,v) \leq \kappa \left(\Vert [{\it{\Delta}} u] \Vert_{L^2({\it{\Gamma}})} + \Vert [\partial_\nu {\it{\Delta}} u]\Vert_{L^2({\it{\Gamma}})} \right) \left( \Vert P_{\it{\Gamma}}^1 T_{\it{\Gamma}}^1 v\Vert_{L^2({\it{\Gamma}})} + \Vert P_{\it{\Gamma}}^2 T_{\it{\Gamma}}^2 v\Vert_{L^2({\it{\Gamma}})} \right)\text{,} \end{align*} which was to be shown. □ Now we are in the situation to show the main results, namely the discretization error estimate for the discretizations (5.3) and (5.5). Theorem 5.2 Let $$u^h_\varepsilon \in S_h$$ be the solution of the discrete soft curve problem (5.4) and assume that $$\varepsilon_1 = c_1 h^{\lambda_1}$$ and $$\varepsilon_2 = c_2 h^{\lambda_2}$$. For any $$\delta >0$$ define   \begin{align*} \gamma = \min_{i\in\{1,2\}}\left( \frac{1}{2}-\delta, 3-i-2\delta-\frac{\lambda_i}{2}, \frac{\lambda_i}{2} \right)\text{.} \end{align*} Then there exists a constant $$c>0$$ independent of $$h$$ such that   \begin{align*} \left\Vert u - u^h_\varepsilon \right\Vert_{H^2({\it{\Omega}})} \leq c h^{\gamma}\text{.} \end{align*} In particular, $$\left\Vert u - u^h_\varepsilon \right\Vert_{H^2({\it{\Omega}})} \in {\mathcal O}(h^{\frac{1}{2}-\delta})$$ for $$\lambda_1 \in [1-2\delta,3-2\delta]$$ and $$\lambda_2 = 1-2\delta$$. Proof. In this proof, we use the notation $$A \lesssim B$$ whenever $$A \leq c B$$ holds with a constant $$c$$ that is independent of $$\varepsilon_1$$, $$\varepsilon_2$$ and $$h$$. As of Lemma 5.1, we can apply Theorem 3.3. Using this together with the coercivity of $$a(\cdot,\cdot)$$ on $$H_0^2({\it{\Omega}}),$$ we conclude for all $$v\in S_h$$,   \begin{align*} \left\Vert u - u^h_\varepsilon \right\Vert_{H^2({\it{\Omega}})} & \lesssim \left\Vert u - u^h_\varepsilon \right\Vert_{a} \\ & \lesssim \left\Vert u - v \right\Vert_a + \sum_{i=1}^2 \frac{1}{\sqrt{\varepsilon_i}} \Vert u - v \Vert_{b_i} + \sum_{i=1}^2 \sqrt{\varepsilon_i}. \end{align*} Note that $$a(\cdot,\cdot)$$ is continuous on $$H^2({\it{\Omega}})$$ and $$b_i(\cdot,\cdot)$$ is continuous on $$H^{i-\frac{1}{2}+\delta}({\it{\Omega}})$$. Furthermore, $$u \in H^{\frac{5}{2}-\delta}({\it{\Omega}})$$ according to Corollary 4.4. Choosing the interpolant $$v = \overline{u}$$ and applying the interpolation estimate (5.1) yields   \begin{align*} \left\Vert u - u^h_\varepsilon \right\Vert_{H^2({\it{\Omega}})} & \lesssim \Vert u - \overline{u} \Vert_{H^2({\it{\Omega}})} + \sum_{i=1}^2 \frac{1}{\sqrt{\varepsilon_i}} \Vert u - \overline{u} \Vert_{H^{i-\frac{1}{2}+\delta}({\it{\Omega}})} + \sum_{i=1}^2 \sqrt{\varepsilon_i} \\ & \lesssim h^{\frac{1}{2}-\delta} \Vert u \Vert_{H^{\frac{5}{2}-\delta}({\it{\Omega}})} + \sum_{i=1}^2 \frac{1}{\sqrt{\varepsilon_i}} h^{3-i-2\delta} \Vert u \Vert_{H^{\frac{5}{2}-\delta}({\it{\Omega}})} + \sum_{i=1}^2 \sqrt{\varepsilon_i} \\ & \lesssim h^{\frac{1}{2}-\delta} + \sum_{i=1}^2 h^{3-i-2\delta-\frac{\lambda_i}{2}} + \sum_{i=1}^2 h^{\frac{\lambda_i}{2}} \\ & \lesssim h^\gamma\text{.} \end{align*} This proves the assertion. □ We can proceed similarly to prove convergence rates for the soft bulk formulation. Lemma 5.3 Suppose that there is a constant $$c_0 >0$$ such that for every $$h \in I$$, $$K \in [0,2]$$ and $$k\in [0,K]$$ the inverse estimate   \begin{align} \forall\,{v \in S_h}, \ \vert v \vert_{H^K(B)} \leq c_0 h^{k-K} \vert v \vert_{H^k(B)} \end{align} (5.7) is fulfilled. Let $$u$$ be the solution of Problem 2.2 and $$\delta>0$$. Then there exists a constant $$c>0$$ such that for all $$h \in I$$ and $$v \in S_h$$,   \begin{align*} a(u,v) \leq c h^{\min\left(0,s-\frac{3}{2}-\delta\right)} \Vert v \Vert_{b_B}\text{.} \end{align*} Proof. From the proof of Lemma 5.1, we know for all $$v \in S_h$$,   \begin{align*} \frac{1}{\kappa} a(u,v) & = \left( -[\partial_\nu {\it{\Delta}} u], v \right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u],\partial_\nu v\right)_{L^2({\it{\Gamma}})} \\ & = \left( -[\partial_\nu {\it{\Delta}} u], T_{\it{\Gamma}}^1 P_B T_Bv \right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u], T_{\it{\Gamma}}^2 P_BT_Bv\right)_{L^2({\it{\Gamma}})} \\ & \quad + \left( -[\partial_\nu {\it{\Delta}} u], T_{\it{\Gamma}}^1 (\operatorname{id}_B - P_B) T_B v \right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u], T_{\it{\Gamma}}^2 (\operatorname{id}_B - P_B) T_Bv\right)_{L^2({\it{\Gamma}})}\text{.} \end{align*} Also as in the proof of Lemma 5.1, there exists a function $$w \in V_0 \subseteq H_0^2({\it{\Omega}})$$ such that $$T_{\it{\Gamma}}^i w = T_{\it{\Gamma}}^i\left(\operatorname{id}_B - P_B\right) T_B v$$ for $$i\in\{1,2\}$$. From this we conclude   \begin{align*} 0 = \frac{1}{\kappa} a(u,w) = \left( -[\partial_\nu {\it{\Delta}} u], T_{\it{\Gamma}}^1 (\operatorname{id}_B - P_B) T_B v \right)_{L^2({\it{\Gamma}})} + \left( [{\it{\Delta}} u], T_{\it{\Gamma}}^2 (\operatorname{id}_B - P_B) T_Bv\right)_{L^2({\it{\Gamma}})} \end{align*} as in the proof of Lemma 5.1. Furthermore, note that the operators $$T_{\it{\Gamma}}^i$$ are continuous and linear over $$H^{\frac{3}{2}+\delta}(B)$$ and thus the corresponding operator norms $${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {T_{{\it{\Gamma}}}^i} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_\delta$$ are bounded. It follows that   \begin{align*} \frac{1}{\kappa} a(u,v) & = \left( -[\partial_\nu {\it{\Delta}} u], T_{\it{\Gamma}}^1 P_B T_Bv \right)_{L^2({\it{\Gamma}})} +\left ( [{\it{\Delta}} u], T_{\it{\Gamma}}^2 P_BT_Bv\right)_{L^2({\it{\Gamma}})} \\ & \leq \Vert [\partial_\nu {\it{\Delta}} u] \Vert_{L^2({\it{\Gamma}})}{\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {T_{\it{\Gamma}}^1} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert} _{\delta} \Vert P_B T_B v \Vert_{H^{\frac{3}{2}+\delta}(B)} \\ & \quad + \Vert [{\it{\Delta}} u] \Vert_{L^2({\it{\Gamma}})} {\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {T_{\it{\Gamma}}^2}\right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{\delta} \Vert P_B T_B v \Vert_{H^{\frac{3}{2}+\delta}(B)}\text{.} \end{align*} Together with (5.7) this implies   \begin{align*} a(u,v) & \leq c h^{\min(0,s-\frac{3}{2}-\delta)}\Vert P_BT_Bv \Vert_{H^s(B),} \end{align*} which proves the assertion. □ Note that (5.7) is indeed an assumption. We cannot apply the standard finite element inverse estimates here as the grid is in general not matched to the domain $$B$$. A closer look to the proof of the standard inverse estimate reveals that the constant in (5.7) would go to infinity as $$\left\vert{E_h\cap B}\right\vert / \left\vert{E_h}\right\vert \to 0$$ for $$I \ni h \to 0$$ and some element $$E_h$$ in the grid associated to $$S_h$$. This means that the above inverse estimate over $$B$$ is only valid for grids that resolve $$B$$ sufficiently well. Theorem 5.4 Let the assumptions from Lemma 5.3 hold, $$u^h_\varepsilon \in S_h$$ the solution of the discrete soft bulk problem (5.6) and $$\varepsilon = c_0 h^\lambda$$ for some $$\delta > 0$$. Define   \begin{align*} \gamma = \min\left( \frac{1}{2}-\delta, \frac{5}{2}-\delta-s-\frac{\lambda}{2}, \min\left(0,s-\frac{3}{2}-\delta\right)+\frac{\lambda}{2}\right)\text{.} \end{align*} Then there exists a constant $$c>0$$, independent of $$h$$ such that   \begin{align*} \left\Vert u - u^h_\varepsilon \right\Vert_{H^2({\it{\Omega}})} & \leq c h^\gamma\text{.} \end{align*} In particular, $$\left\Vert u-u^h_\varepsilon\right\Vert_{H^2({\it{\Omega}})} \in \mathcal{O}\left(h^{\frac{1}{2}-\delta}\right)$$ for $$s \leq \frac{3}{2}+\delta$$ and $$\lambda = 4-2s$$. Proof. As of Lemma 5.3, we can apply Theorem 3.3. Noting the continuity of $$b_B(\cdot,\cdot)$$ on $$H^s({\it{\Omega}}),$$ we obtain as in the proof of Theorem 5.2 that   \begin{align*} \left\Vert u - u^h_\varepsilon \right\Vert_{H^2({\it{\Omega}})} & \lesssim \Vert u - v \Vert_{H^2({\it{\Omega}})}({\it{\Omega}}) + \frac{1}{\sqrt{\varepsilon}} \Vert u - v\Vert_{H^s({\it{\Omega}})} + h^{\min\left(0,s-\frac{3}{2}-\delta\right)}\sqrt{\varepsilon} \text{.} \end{align*} Again choosing the interpolant $$v = \overline{u}$$ and applying the interpolation estimate (5.1) finally gives   \begin{align*} \left\Vert u - u^h_\varepsilon \right\Vert_{H^2({\it{\Omega}})} & \lesssim h^{\frac{1}{2}-\delta} + h^{\frac{5}{2}-\delta-s-\frac{\lambda}{2}} + h^{\min\left(0,s-\frac{3}{2}-\delta\right)+\frac{\lambda}{2}} \lesssim h^\gamma \end{align*} and thus proves the assertion. □ 6. Numerical example 6.1. A symmetric example problem with known exact solution We select a problem with a known, analytically computable solution. We consider the circular domain $$\hat{{\it{\Omega}}} = \{x \in \mathbb{R}^2 \mathrel\vert |x| < r_2\}$$ with one embedded particle $$\hat{B} = \{x \in \mathbb{R}^2 \mathrel\vert |x| < r_1\}$$ for $$0 < r_1 < r_2 < 1$$. On the boundary of $$\hat{{\it{\Omega}}}$$ we consider homogeneous Dirichlet boundary conditions, whereas on the particle boundary $${\it{\Gamma}}_1 = \partial B_1$$, we impose variable-height coupling conditions of the form (2.1), i.e.,   \begin{align} u|_{\partial B_1} &= f^1_1 + \gamma^1, & \partial_\nu u|_{\partial B_1} &= f^1_2 \end{align} (6.1) for variable $$\gamma^1 \in \mathbb{R}$$ with the given functions   \begin{align*} f_1^1(re^{{\rm i}\theta}) &= \cos(n \theta), & f_2^1(x) &= 0. \end{align*} We select the parameters   \begin{align*} r_2 &= \tfrac23, & r_1 &= \tfrac13, & \kappa &= 1, & \sigma &= 0, & n &=4. \end{align*} Then the exact solution of Problem 2.1 on $$\hat{{\it{\Omega}}}$$ is the fourfold symmetric function $$u \in H^2(\hat{{\it{\Omega}}})$$ given by   \begin{align*} u(re^{{\rm i}\theta}) = \begin{cases} \cos(4\theta) \, (c_1 r^4 + c_2 r^6) & \text{ for}\, r \in (0,r_1], \\ \cos(4\theta) \, ( c_3 r^{-2} + c_4 r^{-4} + c_5 r^4 + c_6 r^6 ) & \text{ for}\, r \in [r_1,r_2], \\ 0 & \text{ otherwise}, \end{cases} \end{align*} with the constants   \begin{align*} c_1 & = 243, \qquad & c_2 & = -1458, \qquad & c_3 & =-6909/689, \\ c_4 & = -7936/502281, \qquad & c_5 & = 11502/689, \qquad & c_6 & = 44288/167427 \text{.} \end{align*} To discretize the circular domain $$\hat{{\it{\Omega}}}$$ with the presented Bogner–Fox–Schmit discretization, we embed $$\hat{{\it{\Omega}}}$$ into the larger rectangular domain $${\it{\Omega}} = [-1,1]^2$$ and treat the condition on $${\it{\Gamma}}_2 = \partial \hat{{\it{\Omega}}}$$ like a second particle boundary condition. It is easily seen that the solution of Problem 2.1 on $${\it{\Omega}}$$ is obtained by extending the solution given on $$\hat{{\it{\Omega}}}$$ by zero. The solution of this problem is depicted in Fig. 2. Fig. 2. View largeDownload slide Exact solution of the symmetric example problem. Left: top view. Right: rendered three-dimensional view. Fig. 2. View largeDownload slide Exact solution of the symmetric example problem. Left: top view. Right: rendered three-dimensional view. 6.2. Discretization, quadrature and error measurement For the finite element discretization $$S_h \subseteq H_0^2({\it{\Omega}}),$$ we divide $${\it{\Omega}}$$ into uniform squares with edges of length $$h$$ and equip the resulting grid with a Bogner–Fox–Schmit finite element basis. The boundary conditions over $$\partial{\it{\Omega}}$$ are enforced by setting the corresponding degrees of freedom to zero while the boundary conditions over the $${\it{\Gamma}}_i$$ are replaced by penalized constraints. The resulting soft curve problem with penalized curve constraints is then given by (5.4) while the soft bulk problem with penalized bulk constraints is given by (5.6), where the area of the virtual second particle is now given by $$B_2 = {\it{\Omega}} \setminus \hat{{\it{\Omega}}}$$. We note that we drop the projection for the constraints on $${\it{\Gamma}}_2$$ and $$B_2$$ because we do not allow variable height there. While our analysis is formulated only in terms of the variable-height constraints incorporated using projections, this is in fact no limitation, since all arguments directly carry over to the case without these projections. Since the assembly of the bilinear forms $$b^i_{\it{\Gamma}}(\cdot,\cdot)$$ and $$b_B(\cdot,\cdot)$$ involves integrals over the curves $${\it{\Gamma}}_i$$ or the nontrivial domains $$B_i$$, respectively, we briefly mention how these can be approximated. Integrals over full grid cells are evaluated exactly using standard tensor Gaussian quadrature rules. For quadrature over the $${\it{\Gamma}}_i$$, we use a trigonometric Gaussian quadrature as described in Da Fies & Vianello (2012) which in our case is both exact for integration of finite element functions and for integration of $$u$$. Integrals over the nonrectangular domains $$B_i$$ are approximated using the local parameterization quadrature method introduced in Olshanskii & Safin (2016), which is no longer exact. The discretization error is measured in terms of the norm $$\| {\it{\Delta}} (u^h_\varepsilon - u)\|_{L^2({\it{\Omega}})}$$ on $$H_0^2({\it{\Omega}})$$. It is important to note that the error cannot be computed accurately by simple quadrature over the grid elements, because $$u$$ is not smooth across the $${\it{\Gamma}}_i$$. Instead we split the integration into the subdomains $$B_1$$, $$B_2 = {\it{\Omega}} \setminus \hat{{\it{\Omega}}}$$ and $${\it{\Omega}} \setminus (B_1 \cup B_2),$$ where we apply the quadrature rules described above. 6.3. Results for the soft curve formulation In case of the soft curve formulation, we expect convergence with order $${\mathcal O}(h^{\frac12})$$ due to Theorem 5.2 if the penalty parameters are chosen suitably. For the numerical examples, we selected penalty parameters $$\varepsilon = ( c h^{\lambda_1}, c h )$$ with $$\lambda_1 \in [1,3]$$ according to Theorem 5.2. The constant was fixed to be $$c=10^{-3}$$. Figure 3 shows the behavior of the $$H^2$$-discretization error over the mesh size $$h$$ for uniform grids and penalty parameters with $$\lambda_1 \in \{1, 2, 3\}$$. In accordance with Theorem 5.2, we observe convergence of at least order $${\mathcal O}(h^{1/2})$$. While the observed rate seems to behave like $${\mathcal O}(h^{1/2})$$ for $$\lambda_1 \in \{2,3\}$$ it seems that the rate is slightly better for $$\lambda_1 = 1$$. A possible explanation for this observation is that the error contribution for the term $$\|\cdot\|_{b_1}$$ in the proof of Theorem 5.2 can be improved by smaller values of $$\lambda_1$$. However, this does not lead to a better theoretical bound due to other dominating terms. Fig. 3. View largeDownload slide $$H^2$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 3. View largeDownload slide $$H^2$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Figures 4 and 5 show the $$H^1$$- and $$L^2$$-discretization errors, respectively, for the same set of example problems. In both cases, we observe convergence with order $${\mathcal O}(h)$$. While improved convergence order for weaker norms is a well-known property of many discretizations, we emphasize that this is not covered by the theory presented here and may be considered in the future with a refined analysis. Fig. 4. View largeDownload slide $$H^1$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 4. View largeDownload slide $$H^1$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 5. View largeDownload slide $$L^2$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 5. View largeDownload slide $$L^2$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. 6.4. Results for the soft bulk formulation Because we use uniform grids and do not refine with regard to the geometry of the $${\it{\Gamma}}_i$$, we lose control over the inverse estimate assumption (5.7). Thus, we cannot prove condition (3.4) in general. Assuming that (3.4) still holds true, we could use Theorem 5.4 in the case of exact integration. However, as discussed above, we approximate integrals using only the quadrature from Olshanskii & Safin (2016). If this is accurate enough, we can in view of the Strang-type result in Proposition 3.4 expect convergence of order $${\mathcal O}(h^{1/2})$$ when using $$\varepsilon = c h^{4-2s}$$ and the penalty norm $$\|\cdot\|_{H^s(B_i)}$$ for $$s\in [0,1]$$. For our numerical experiments we considered $$c = 10^{-3}$$ and $$s \in \{0,1\}$$. Figure 6 shows the behavior of the discretization error over the mesh size $$h$$ for uniform grids and $$s\in \{0, 1\}$$. Again the observed order is in accordance with the expected order $${\mathcal O}(h^{1/2})$$. However, especially for $$s=1,$$ this is perturbed by strong oscillations. This is maybe due to the fact that the constant resulting from (3.4) is not uniformly bounded and can strongly vary depending on how the curve $${\it{\Gamma}}$$ intersects mesh elements. Again the $$H^1$$- and $$L^2$$-errors depicted in Figs 7 and 8, respectively, exhibit similar behavior but the improved convergence order $${\mathcal O}(h)$$. Fig. 6. View largeDownload slide $$H^2$$-errors for the soft bulk formulation for $$s=0$$ and $$s=1$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 6. View largeDownload slide $$H^2$$-errors for the soft bulk formulation for $$s=0$$ and $$s=1$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 7. View largeDownload slide $$H^1$$-errors for the soft bulk formulation for $$s=0$$ and $$s=1$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 7. View largeDownload slide $$H^1$$-errors for the soft bulk formulation for $$s=0$$ and $$s=1$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 8. View largeDownload slide $$L^2$$-errors for the soft bulk formulation for $$s=0$$ and $$s=1$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. Fig. 8. View largeDownload slide $$L^2$$-errors for the soft bulk formulation for $$s=0$$ and $$s=1$$ (from left to right) over the grid size $$h = 1/N$$ for $$N = 8,\dots,75$$. 6.5. A nonsymmetric example problem While the first example allowed us to compute the exact solution due to its symmetry properties, we will now assess the discretization for a nonsymmetric example with several particles but without known exact solution. We consider the domain $${\it{\Omega}} = [-1,1]^2$$ with homogeneous Dirichlet boundary conditions on $$\partial {\it{\Omega}}$$ and we again select the parameters   \begin{align*} \kappa &= 1, & \sigma &= 0. \end{align*} This time we consider four elliptical particles embedded nonsymmetrically into the membrane domain $${\it{\Omega}}$$. The ellipses’ major and minor axes have lengths $$0.4$$ and $$0.2$$, respectively. Their positions and orientations are depicted in Fig. 9. Fig. 9. View largeDownload slide Solution of the nonsymmetric example problem. Left: top view. Right: rendered three-dimensional view. Fig. 9. View largeDownload slide Solution of the nonsymmetric example problem. Left: top view. Right: rendered three-dimensional view. Regarding the particle–membrane coupling conditions, we assume that all particles have a constant height profile $$f^i_1 = 0$$ and a constant slope $$f^i_2$$ varying for different particles. More precisely, we use   \begin{align*} f^1_2 &= 1, & f^2_2 &= 1, & f^3_2 &= -1, & f^4_2 &= -1. \end{align*} The solution resulting from that setup is depicted in Fig. 9. Since we do not know the exact solution for this problem, we estimate the error by comparing our discrete solutions to a reference solution computed on a finer grid. We computed the solutions for the different problem formulations on uniform grids with mesh size $$h=2^{-k}$$ for $$k= 3,\dots,6,$$ while the reference solution was computed on a uniform grid with mesh size $$h=2^{-8}$$ using the soft curve formulation with $$\varepsilon = 10^{-3}(h^3,h)$$. The (approximate) discretization errors for the soft curve and soft bulk formulations are shown in Figs 10 and 11, respectively. As before, we used $$\varepsilon = (ch^\lambda, ch)$$ with $$\lambda \in \{1,2,3\}$$ for the soft curve formulation and $$\varepsilon = ch^{4-2s}$$ with $$s \in \{0,1\}$$ for the soft bulk formulation. In both cases, we selected $$c = 10^{-3}$$. Fig. 10. View largeDownload slide $$H^2$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$. Fig. 10. View largeDownload slide $$H^2$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$. Fig. 11. View largeDownload slide $$H^2$$-errors for the soft bulk formulation for $$s=0$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,4,5$$ and $$s=1$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$ (from left to right). Fig. 11. View largeDownload slide $$H^2$$-errors for the soft bulk formulation for $$s=0$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,4,5$$ and $$s=1$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$ (from left to right). For the soft bulk formulation, the projection $$P_{B}$$ may lead to a densely populated matrix for those degrees of freedom located inside the particles. To avoid the resulting computational effort for $$s=1,$$ we used the $$H^1(B_i)$$-norm in the form   \begin{align*} \|v\|_{H^1(B_i)}^2 = \|\nabla v\|_{L^2(B_i)}^2 + \Bigl(\int_{B_i} v \, \,{\rm d}x\Bigr)^2. \end{align*} Then the penalty term that incorporates the $$H^1(B_i)$$-projections takes the form   \begin{align*} b_{B}(w,v) = (\nabla w, \nabla v)_{L^2(B)}, \end{align*} which results in a sparse matrix again. For both formulations, soft curve and soft bulk, we observe a rate that is slightly better than the expected order of convergence, namely $${\mathcal O}\left(h^{1/2}\right)$$. Note that we restricted the mesh size $$h=2^{-k}$$ to powers of 2 in this example to simplify comparison with the reference solution. Furthermore, in the case of soft bulk constraints with $$s=0$$, we were not able to carry out the computations for $$h=2^{-6}$$ due to hardware restrictions. This explains the much smaller number of data points in the plots. As before, we also computed the $$H^1$$- and $$L^2$$-errors for the same set of example problems. For the soft curve formulation, we again observe that the order of the $$H^1$$- and $$L^2$$-errors depicted in Figs 12 and 13, respectively, is essentially squared in comparison with the $$H^2$$-error leading to a convergence order which is approximately $${\mathcal O}(h)$$. The situation is different for the soft bulk formulation. Here the $$H^1$$-error depicted in Fig. 14 is of order $${\mathcal O}(h^{3/2})$$ and the $$L^2$$-error depicted in Fig. 15 is of order $${\mathcal O}(h^{5/2})$$. While the improved $$H^1$$- and $$L^2$$-order is not covered by the presented theory anyway, we also cannot explain the surprising difference in the observed order for the symmetric and nonsymmetric example problems. Fig. 12. View largeDownload slide $$H^1$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$. Fig. 12. View largeDownload slide $$H^1$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$. Fig. 13. View largeDownload slide $$L^2$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$. Fig. 13. View largeDownload slide $$L^2$$-errors for the soft curve formulation for $$\lambda=1$$, $$\lambda=2$$ and $$\lambda=3$$ (from left to right) over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$. Fig. 14. View largeDownload slide $$H^1$$-errors for the soft bulk formulation for $$s=0$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,4,5$$ and $$s=1$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$ (from left to right). Fig. 14. View largeDownload slide $$H^1$$-errors for the soft bulk formulation for $$s=0$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,4,5$$ and $$s=1$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$ (from left to right). Fig. 15. View largeDownload slide $$L^2$$-errors for the soft bulk formulation for $$s=0$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,4,5$$ and $$s=1$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$ (from left to right). Fig. 15. View largeDownload slide $$L^2$$-errors for the soft bulk formulation for $$s=0$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,4,5$$ and $$s=1$$ over the grid sizes $$h=2^{-k}$$, $$k= 3,\dots,6$$ (from left to right). Funding This research has been partially supported by Deutsche Forschungsgemeinschaft (DFG) through grant CRC 1114 ‘Scaling Cascades in Complex Systems’, Project AP01 ‘Particles in lipid bilayers’. Appendix Proof of Lemma 4.3. Let $$u_1 := u|_{{\it{\Omega}}_B}$$ and $$u_2 := u|_B$$ as well as $${\it{\Omega}}_1 := {\it{\Omega}}_B$$ and $${\it{\Omega}}_2 := B$$. Define $$\varphi_{\alpha}(x,y) := {(\vert D^\alpha u(x) - D^\alpha u(y) \vert^2)}/{(\Vert x-y\Vert^{2s+2})}$$. For $$s \in (0,\frac{1}{2}),$$ the Sobolev–Slobodeckij norm of $$u$$ is given by   \begin{align*} \Vert u\Vert_{H^{2+s}({\it{\Omega}})}^2 & = \Vert u\Vert_{H^2({\it{\Omega}})}^2 + \max_{\vert\alpha\vert = 2} \Vert D^\alpha u\Vert_{H^s({\it{\Omega}})}^2 \\ & = \Vert u\Vert_{H^2({\it{\Omega}})}^2 + \max_{\vert\alpha\vert = 2} \int_{{\it{\Omega}} \times {\it{\Omega}}} \varphi_{\alpha}(x,y) \, \mathrm{d}(x,y)\text{.} \end{align*} Splitting the domain of integration yields   \begin{align*} \int_{{\it{\Omega}} \times {\it{\Omega}}} \varphi_{\alpha}(x,y) \, \mathrm{d}(x,y) & = \left\Vert D^\alpha u_1 \right\Vert_{H^s({\it{\Omega}}_1)}^2 + \left\Vert D^\alpha u_2\right\Vert_{H^s({\it{\Omega}}_2)}^2 + 2\int_{{\it{\Omega}}_1 \times {\it{\Omega}}_2} \varphi_{\alpha}(x,y) \, \mathrm{d}(x,y)\text{.} \end{align*} Using the layer-cake principle, we write   \begin{align*} \int_{{\it{\Omega}}_1 \times {\it{\Omega}}_2} \varphi_{\alpha}(x,y) = \int_0^\infty \left \vert \left\{ (x,y) \in {\it{\Omega}}_1\times {\it{\Omega}}_2 \mid \varphi_{\alpha}(x,y) \geq t \right\} \right\vert \, \mathrm{d}t\text{.} \end{align*} By the assumptions from Section 2 it holds that $${\it{\Omega}}_1$$ and $${\it{\Omega}}_2$$ are Lipschitz domains. By virtue of the Sobolev embedding theorems, we conclude $$u_i \in C^2(\overline{{\it{\Omega}}_i})$$ (see, e.g., Adams, 1975, Theorem 5.4). Consequently, the constant $$C := \max_{\vert\alpha\vert = 2} \sup_{{\it{\Omega}}_1 \times {\it{\Omega}}_2} \vert D^\alpha u_1(x) - D^\alpha u_2(y)\vert^2 \in \mathbb{R}_{>0}$$ exists and thus the implication   \begin{align*} \frac{\vert D^\alpha u(x) - D^\alpha u(y)\vert^2}{\Vert x-y\Vert^{2s+2}} \geq t \ \Longrightarrow \ \Vert x-y\Vert \leq C^{\frac{1}{2s+2}} t^{-\frac{1}{2s+2}} =: \widetilde{C} t^{-\frac{1}{2s+2}} \end{align*} holds for all $$(x,y) \in {\it{\Omega}}_1 \times {\it{\Omega}}_2$$. Let $${\it{\Gamma}} := \partial{\it{\Omega}}_1 \cap \partial {\it{\Omega}}_2$$. From   \begin{align*} \left\vert \left\{ (x,y) \in {\it{\Omega}}_1\times {\it{\Omega}}_2 \mid \Vert x-y\Vert \leq \widetilde{C} t^{-\frac{1}{2s+2}} \right\} \right\vert \leq 4\vert{\it{\Gamma}}\vert (\widetilde{C} t^{-\frac{1}{2s+2}})(\widetilde{C} t^{-\frac{1}{2s+2}})^2 =: \hat{C} t^{\frac{-3}{2s+2}} \end{align*} we infer   \begin{align*} & \int_0^\infty \left \vert \left\{ (x,y) \in {\it{\Omega}}_1\times {\it{\Omega}}_2 \mid \varphi_{\alpha}(x,y) \geq t \right\} \right\vert \, \mathrm{d}t \\ & \leq \vert {\it{\Omega}} \vert^2 + \int_1^\infty \left \vert \left\{ (x,y) \in {\it{\Omega}}_1\times {\it{\Omega}}_2 \mid \varphi_{\alpha}(x,y) \geq t \right\} \right\vert \, \mathrm{d}t \\ & \leq \vert {\it{\Omega}} \vert^2 + \int_1^\infty \hat{C} t^{\frac{-3}{2s+2}} \, \mathrm{d}t\text{.} \end{align*} This expression is finite for $$s < \frac{1}{2}$$. This implies $$u \in H^{2+\frac{1}{2}-\delta}({\it{\Omega}})$$ for all $$\delta >0$$ as was to be shown. □ Proof of Proposition 3.4. For brevity we define $$\widetilde{u} := \widetilde{u}^X_\varepsilon$$ and assume without loss of generality that $$\varepsilon = (\varepsilon_i)_{i=1,\dots,m} \leq 1$$ holds componentwise. Let $$v \in X$$ and define $$\widetilde{w}:= v -\widetilde{u}$$. We have by (3.11),   \begin{align*} \left\Vert{\widetilde{w}}\right\Vert_{a_\varepsilon}^2 & = \left\Vert{\widetilde{w}}\right\Vert_a^2 + \sum_{i=1}^m \frac{1}{\varepsilon_i} \left\Vert{\widetilde{w}}\right\Vert_{b_i}^2 \\ & \leq \left\Vert{\widetilde{w}}\right\Vert_a^2 + \sum_{i=1}^m \frac{1}{\varepsilon_i} \left\Vert{\widetilde{w}}\right\Vert_{\widetilde{b}_i}^2 + \sum_{i=1}^m \frac{c_i}{\varepsilon_i} \left\Vert{\widetilde{w}}\right\Vert^2\text{.} \end{align*} Using the continuity of $$a$$, the coercivity of $$\widetilde{a}_\varepsilon$$ and $$\varepsilon \leq 1,$$ this leads to the inequality   \begin{align} \left\Vert{\widetilde{w}}\right\Vert_{a_\varepsilon}^2 & \leq \left(\frac{\left\Vert{a}\right\Vert}{\widetilde{\alpha}} + 1 + \sum_{i=1}^m \frac{c_i}{\varepsilon_i \widetilde{\alpha}} \right) \left\Vert{\widetilde{w}}\right\Vert_{\widetilde{a}_\varepsilon}^2\text{.} \end{align} (A.1) On the other hand, we have the equality   \begin{align*} \left\Vert{\widetilde{w}}\right\Vert_{\widetilde{a}_\varepsilon}^2 & = \widetilde{a}_\varepsilon(v,\widetilde{w}) - \widetilde{a}_\varepsilon(\widetilde{u},\widetilde{w}) + \ell_\varepsilon(\widetilde{w}) - a_\varepsilon(u,\widetilde{w}) + a_\varepsilon(v,\widetilde{w}) - a_\varepsilon(v,\widetilde{w}), \end{align*} which by application of (3.10), the Cauchy–Schwarz inequality, coercivity of $$a_\varepsilon$$, definition of $$\widetilde{a}_\varepsilon$$ and $$\widetilde{\ell}_\varepsilon$$ and $$\varepsilon \leq 1$$ yields   \begin{align} \begin{aligned} \left\Vert{\widetilde{w}}\right\Vert_{\widetilde{a}_\varepsilon}^2 & = a_\varepsilon(v-u,\widetilde{w}) + (\widetilde{a}_\varepsilon-a_\varepsilon)(v,\widetilde{w}) + (\ell_\varepsilon-\widetilde{\ell}_\varepsilon)(\widetilde{w}) \\ & \leq \left\Vert{v-u}\right\Vert_{a_\varepsilon} \left\Vert{\widetilde{w}}\right\Vert_{a_\varepsilon} \\ & \quad + \frac{\left\vert{(\widetilde{a}-a)(v,\widetilde{w})}\right\vert + \sum_{i=1}^m \frac{1}{\varepsilon_i} \left\vert{(\widetilde{b}_i-b_i)(u-v,\widetilde{w})}\right\vert + \left\vert{(\ell-\widetilde{\ell})(\widetilde{w})}\right\vert}{\sqrt{\alpha} \left\Vert{\widetilde{w}}\right\Vert} \left\Vert{\widetilde{w}}\right\Vert_{a_\varepsilon}. \end{aligned} \end{align} (A.2) Then from the triangle inequality and (A.2) inserted into (A.1), we get   \begin{align*} \left\Vert{u - \widetilde{u}}\right\Vert_{a_\varepsilon} & \leq \left\Vert{u - v}\right\Vert_{a_\varepsilon} + \left\Vert{v - \widetilde{u}}\right\Vert_{a_\varepsilon} \\ & \leq \left\Vert{u - v}\right\Vert_{a_\varepsilon} + \left(\frac{\left\Vert{a}\right\Vert}{\widetilde{\alpha}} + 1 + \sum_{i=1}^m \frac{c_i}{\varepsilon_i \widetilde{\alpha}} \right) \cdot \left( \left\Vert{v-u}\right\Vert_{a_\varepsilon} + \phantom{\frac{\sum_{i=1}^m \frac{1}{\varepsilon_i}(\widetilde{b}}{\left\Vert{\widetilde{w}}\right\Vert}}\right. \\ & \qquad + \left. \frac{\left\vert{(\widetilde{a}-a)(v,\widetilde{w})}\right\vert + \sum_{i=1}^m \frac{1}{\varepsilon_i} \left\vert{(\widetilde{b}_i-b_i)(u-v,\widetilde{w})}\right\vert + \left\vert{(\ell-\widetilde{\ell})(\widetilde{w})}\right\vert}{\sqrt{\alpha} \left\Vert{\widetilde{w}}\right\Vert} \right), \end{align*} which proves the statement by taking the infimum over all $$v \in X$$ and replacing $$\widetilde{w} = v - \widetilde{u}$$ by the supremum over all $$w \in X$$. □ References Adams R. ( 1975) Sobolev Spaces . Pure and Applied Mathematics. New York-London: Academic Press. Babuška I. ( 1973) The finite element method with penalty. Math. Comp. , 27, 221– 228. Google Scholar CrossRef Search ADS   Bahrami A. H., Raatz M., Agudo-Canalejo J., Michel R., Curtis E. M., Hall C. K., Gradzielski M., Lipowsky R. & Weikl T. R. ( 2014) Wrapping of nanoparticles by membranes. Adv. Colloid Interface Sci. , 208, 214– 224. Special issue in honour of Wolfgang Helfrich. Google Scholar CrossRef Search ADS PubMed  Barrett J. W. & Elliott C. M. ( 1986) Finite element approximation of the Dirichlet problem using the boundary penalty method. Numer. Math. , 49, 343– 366. Google Scholar CrossRef Search ADS   Blum H. & Rannacher R. ( 1980) On the Boundary Value Problem of the Biharmonic Operator on Domains with Angular Corners . Preprint: Sonderforschungsbereich Approximation und Mathematische Optimierung in einer Anwendungsbezogenen Mathematik. Mathematical Methods in the Applied Sciences, John Wiley & Sons, Ltd., SFB 72. Ciarlet P. G. ( 1978) The Finite Element Method for Elliptic Problems . Amsterdam-New York-Oxford: North-Holland Publishing Company. Da Fies G. & Vianello M. ( 2012) Trigonometric Gaussian quadrature on subintervals of the period. Electron. Trans. Numer. Anal. , 39, 102– 112. Dommersnes P. G. & Fournier J.-B. ( 1999) Casimir and mean-field interactions between membrane inclusions subject to external torques. Europhys. Lett. , 46, 256. Google Scholar CrossRef Search ADS   Dommersnes P. G. & Fournier J.-B. ( 2002) The many-body problem for anisotropic membrane inclusions and the self-assembly of ‘saddle’ defects into an ‘egg carton’. Biophys. J. , 83, 2898– 2905. Google Scholar CrossRef Search ADS PubMed  Elliott C. M., Gräser C., Hobbs G., Kornhuber R. & Wolf M.-W. ( 2016) A variational approach to particles in lipid membranes. Arch. Ration. Mech. Anal. , 222, 1011– 1075. Google Scholar CrossRef Search ADS   Grisvard P. ( 1985) Elliptic Problems in Nonsmooth Domains . Monographs and Studies in Mathematics. Boston, London, Melbourne: Massachusetts, USA-London: Pitman Publishing. Helfrich P. & Jakobsson E. ( 1990) Calculation of deformation energies and conformations in lipid membranes containing gramicidin channels. Biophys. J. , 57, 1075– 1084. Google Scholar CrossRef Search ADS PubMed  Lions J. & Magenes E. ( 1972) Non-homogeneous Boundary Value Problems and Applications . Massachusetts, USA-London: Springer. Lunardi A. ( 2009) Interpolation Theory  ( Publications of the Scuola Normale Superiore). Scuola Normale Superiore, Pisa PI, Italy: Edizioni della Normale. Nitsche J. ( 1971) Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind. Abh. Math. Semin. Univ. Hambg. , 36, 9– 15. Google Scholar CrossRef Search ADS   Olshanskii M. A. & Safin D. ( 2016) Numerical integration over implicitly defined domains for higher order unfitted finite element methods. Lobachevskii J. Math. , 37, 582– 596. Google Scholar CrossRef Search ADS   Weikl T. R., Kozlov M. M. & Helfrich W. ( 1998) Interaction of conical membrane inclusions: effect of lateral tension. Phys. Rev. E , 57, 6988. Google Scholar CrossRef Search ADS   © The authors 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

### Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Jan 4, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations