title: A C0 interior penalty finite element method for flexoelectricity
authors: Jordi Ventura David Codony Sonia Fernández-Méndez
text: 
  Efficient closed-form estimation of large spatial autoregressionsI am grateful to the associate editor, three referees and Peter Robinson for very helpful comments that improved the paper significantly, and to Christos Kolympiris for permission to use the data in the empirical illustration.
      Abhimanyu GuptaDepartment of Economics, University of Essex, Wivenhoe Park, Colchester CO4 3SQ, UK. E-mail: a.gupta@essex.ac.uk. Research supported by ESRC grant ES/R006032/1. 
      December 30, 2023
  ===================================================================================================================================================================================================================================================================================================
  
  
   
  We propose a 𝒞^0 Interior Penalty Method (C0-IPM) for the computational modelling of flexoelectricity, with application also to strain gradient elasticity, as a simplified case. Standard high-order 𝒞^0 finite element approximations, with nodal basis, are considered. The proposed C0-IPM formulation involves second derivatives in the interior of the elements, plus integrals on the mesh faces (sides in 2D), that impose 𝒞^1 continuity of the displacement in weak form.  The formulation is stable for large enough interior penalty parameter, which can be estimated solving an eigenvalue problem. The applicability and convergence of the method is demonstrated with 2D and 3D numerical examples.
  
  
  Keywords: 
  4th order PDE 𝒞^0 finite elements interior penalty method strain gradient elasticity flexoelectricity
  
  
  
  
  § INTRODUCTION
  
  
  
  The rising interest on microtechnology  evidences the need for mathematical and computational models suitable for small scales, often giving rise to 4^th order Partial Differential Equations (PDEs). In particular flexoelectric effects become relevant, and may be crucial, in the design of small electromechanical devices or for the understanding of physical phenomena <cit.>.
  
  The modelling of flexoelectricity involves a two-way coupling between strain gradient and electric field; strain gradient elasticity is frequently also included in the model to regularise the problem, leading to a system of 4^th order PDEs. 
  
  Several numerical strategies have been recently proposed for the solution of  flexoelectricity problems, based on the use of 𝒞^1 approximation spaces or on mixed formulations.  Mixed formulations split the PDE in two 2^nd order PDEs, allowing the use of 𝒞^0 Finite Element (FE) approximations <cit.>. The approximation spaces, for the primal unknown and for the additional unknowns, must fulfill some conditions for stability that lead to approximation spaces with cumbersome definitions, and difficult extension to 3D or high-order approximations. However, the main drawback of mixed formulations in the high computational cost due to the additional unknowns. 
  
  On other hand, 𝒞^1 approximations can be directly used for the discretization of the weak form in ℋ^2, involving 2^nd order derivatives, without additional unknowns. The first successful attempt in this direction considered a meshless method: the maximum entropy method <cit.>. 
  
  Unfortunately, the computational cost of accurate meshless methods is high, mainly due to the excessive number of integration points for accurate solutions, and the large stencils in the discrete matrices. Aiming to improve the efficiency, a solution based on Isogeometric Analysis (IGA) is proposed in <cit.>. The solution is approximated by means of Non-Uniform Rational B-Splines (NURBS). An interesting critical comparison of IGA, meshless and mixed methods for flexoelectricity can be found at <cit.>, although numerical examples are restricted to simple geometries. The conclusion is that IGA is very efficient on regular grids, corresponding to a transformation of a rectangle grid, where plain B-Spline approximations can be easily defined. 
  However, in a more general context, defining a NURBS approximation with 𝒞^1 continuity in a whole domain with complex shape may not be straightforward, and the numerical integration of the resulting NURBS may be, again, very expensive <cit.>. 
  An efficient alternative for complex domains is the immersed B-Spline method proposed in <cit.>. It considers B-Spline approximations based on a background regular grid, with an embedded domain. The applicability of the proposal is demonstrated with 2D and 3D complex geometries. The weak points of immersed B-Splines are the usual ones in the context of embedded domains: the robust definition of numerical integration in cropped cells (intersected by the domain boundary), which is specially challenging in 3D, and the ill-conditioning problems in the presence of cells with a small portion in the domain, that can be alleviated with specific techniques <cit.>.
  
  
  Here we propose a 𝒞^0 Interior Penalty Method (C0-IPM) for the solution of flexoelectricity. A standard 𝒞^0 FE approximation is considered, and 𝒞^1 continuity between elements is imposed in weak form by means of the Interior Penalty Method (IPM). The procedure for the derivation of the C0-IPM weak form is analogous to the derivation of the IPM in the context of DG methods for 2^nd order PDEs <cit.>, or Nitsche's method for weak imposition of Dirichlet boundary conditions <cit.>, but now applied to the continuity of normal derivatives on element boundaries.
  The resulting weak form involves second derivatives in the interior of the elements, plus integrals on the faces (sides in 2D) of the mesh, that impose 𝒞^1 continuity of the displacement in weak form. 
  C0-IPM formulations overcome the disadvantages of other methods, because they allow the use of standard 𝒞^0 FE approximations. Namely, (i) the computational mesh can be adapted to any geometry, with localised refinement were needed, (ii) there is no need to use embedded discretizations, avoiding the consequent ill-conditioning problems and the definition of special numerical integration for cropped elements, (iii) there are no additional unknowns and (iv) they handle material interfaces in a natural way. In summary, C0-IPM retains the computational efficiency and the versatility that make standard FEs the preferred method for many practitioners in the computational mechanics community.
  
  In <cit.>, C0-IPM formulations, there referred to as continuous/discontinuous finite elements, are applied to several problems modelled by 4^th order PDEs, including Kirchhoff plates and 1D strain gradient elasticity.  Numerical experiments show the applicability of the formulation in both applications, but convergence studies are limited to 1D examples. The C0-IPM formulation is then analysed in <cit.> for the 2D biharmonic equation, with first and second Dirichlet conditions, including a convergence analysis that shows that the method is convergent for p≥ 2, but may have suboptimal convergence depending on the degree and the penalty parameter. An experimental convergence study for Kirchhoff plates can be found at <cit.>. The numerical results demonstrate the applicability of the method for degree greater or equal to 3, and also show slightly suboptimal convergence, that slowly deteriorates for larger penalty parameter, in agreement with the analysis in <cit.>. Variations of C0-IPM have also been applied to strain gradient dependent damage models in <cit.> and to the Cahn-Hilliard equation in <cit.>. 
  
  This paper develops the C0-IPM method for flexoelectricity and, as a simplified case, for strain gradient elasticity, for 2D and 3D computations. Section <ref> presents the problem statement and recalls the weak form in ℋ^2. The derivation of the C0-IPM method for 𝒞^0 approximations, not in ℋ^2, is presented in section <ref>. An eigenvalue problem to determine a large enough penalty parameter, ensuring coercivity of the strain gradient bilinear form, is derived in section <ref>. Finally, in section <ref>, 2D and 3D numerical experiments  demonstrate the applicability of the method and show, as expected, slightly suboptimal convergence under uniform mesh refinement, but still with a robust high-order convergence for p≥ 3.
  
  
  
  
  
  § PROBLEM STATEMENT
  
  
  We consider the model in <cit.>, where flexoelectricity is ruled by the following set of PDEs and boundary conditions:
  
   
  
      · ((,ϕ)-·(,ϕ)) + b = 0    in Ω
      ∇·(,ϕ) - q = 0    in Ω
       = g_1     on Γ^_D_1
      t (,ϕ)  = t_n     on Γ^_N_1
      ∂∂n = g_2     on Γ^_D_2
      r (,ϕ) = r_n     on Γ^_N_2
      j(,ϕ) = j^ext     on  C_N^∂Ω
      ϕ = g_3     on Γ_D^ϕ
      
          w(u,ϕ)=w_n     on Γ_N^ϕ
  
  
  
  where Ω⊂ℝ^ is the domain, the displacement  and the electric potential  ϕ are the unknowns,  and  are the local and double stress tensors and  is the electric displacement tensor, that is, 
  
  
      = ℂ : ε - E·e     ≡  σ_ij = C_ijkℓ ε_kℓ -  E_ℓ e_ℓ ij
       = h ⋮ ∇ε-E·μ     ≡  σ_ijk = h_ijkℓ mn ∂ε_ℓ m∂ x_n - E_ℓ μ_ℓ ijk
      D = κ·E  +e: ε +μ⋮∇ε     ≡  D_ℓ = κ_ℓ m E_m + e_ℓ ijε_ij+ μ_ℓ ijk∂ε_ij∂ x_k
  
  ε is the strain tensor, that is ε_ij=(∂ u_i / ∂ x_j + ∂ u_j / ∂ x_i)/2, =-ϕ is the electric field, 
   C is the elasticity tensor, that depends on the Young modulus E and the Poisson ratio ν, h is the strain-gradient tensor, defined as
  h_ijkℓ mn=l^2 C_ijℓ mδ_kn 
   with the internal length scale parameter l,
  e and  μ are the tensors of piezoelectric and flexoelectric coefficients  and κ contains the dielectricity constants, see appendix B in <cit.> for detailed definitions.
  
  In the previous equations, and in the rest of the document,
   Einstein's notation is assumed. That is, repeated indexes sum over the spatial dimensions. 
   
  The boundary of the domain is split in Dirichlet and Neumann boundaries, for the first and second conditions of the mechanical problem and for the electric problem, that is
  
      ∂Ω = Γ^_D_1∪Γ^_N_1 =  Γ^_D_2∪Γ^_N_2  =  Γ^ϕ_D∪Γ^ϕ_N.
  
  Note that all volume and boundary domains are assumed to be open domains, not including their boundaries.
  
  The first mechanical boundary condition, (<ref>)  or (<ref>), sets the displacement  or the traction
  
      t_i (u,ϕ)= (σ_ij -∂σ_ijk∂ x_k - ∇^S_k σ_ikj) n_j + σ_ijk N_jk,
  
  where ∇^S_k σ_ikj is the surface divergence of σ_ikj, and Ñ is the second order geometry tensor, see <cit.> for details. The second mechanical boundary condition, (<ref>)  or (<ref>), sets the normal derivative of the displacement ∂/∂n or the double traction 
  
      r_i(u,ϕ)=σ_ijkn_jn_k.
  
  
  
  
  The condition (<ref>) sets forces on the Neumann boundary edges. That is, the domain boundary  is assumed to be composed of smooth surfaces (curves in 2D) that are joined on sharp boundary edges (corners in 2D). C_N^∂Ω denotes the union of the boundary edges that are shared by two surfaces with first Neumann conditions, i.e. the edges in the interior of Γ^_N_1. At the edges shared by at least one Dirichlet surface the value on the edge is assumed to be the one set on the surface, i.e. =g_1 for all edges in Γ^_D_1. Line forces (punctual forces in 2D) are defined on boundary edges as
  
      j_i(u,ϕ) = τ_j^L σ_ijℓ^Ln_ℓ^L + τ_j^R σ_ijℓ^Rn_ℓ^R,
  
  being n^L and n^R  the unitary exterior normals on the left and right surfaces sharing the boundary edge, and τ^L and τ^R the tangent vectors on each surface pointing outward and perpendicular to the edge, see an example in figure <ref> left. In 2D, τ^L and τ^R at a corner are just the tangent vectors on each curve sharing the corner and pointing outward, as depicted in figure <ref> right.
  
  Finally, the electric boundary condition,  (<ref>)  or (<ref>), sets the electric potential ϕ or the surface charge density
  
      w(u,ϕ) =-D_ℓ(u,ϕ)  n_ℓ.
  
  
  
  For the sake of simplicity, we initially restrict to the case with second Neumann boundary conditions in the whole boundary, i.e. Γ^_N_2=∂Ω and Γ^_D_2=∅. The treatment of second Dirichlet conditions (<ref>) is commented in Remark <ref>.
  It is also worth mentioning that the rationales in this work can also be applied to other models, including converse flexoelectricity or expressed in terms of polarization, see for instance <cit.>.  
  
  
  If an approximation in ℋ^2(Ω) can be considered, multiplying  (<ref>) by a weighting vector , applying integration by parts twice, and using the symmetries of the stress tensors, leads to
  
      ∫_Ω·b dΩ =    ∫_Ωε() : (,ϕ)  dΩ + ∫_Ωε() ⋮ (,ϕ) dΩ
         - ∫_∂Ωv_i (σ_ij-∂σ_ijk∂ x_k)  n_j  dS - ∫_∂Ω∂ v_i∂ x_j σ_ijk  n_k  dS,
  
   where ε⋮=∂ε_ij/∂ x_k  σ_ijk.
   
  Now, to properly treat boundary conditions, the derivative ∂ v_i/∂ x_j on the boundary is split in normal and tangential derivatives, and the surface divergence theorem is applied to the term with the tangential derivative, leading to 
  
      ∫_Ω·b dΩ = ∫_Ωε() : (,ϕ)  dΩ + ∫_Ωε() ⋮ (,ϕ)  dΩ
         - ∫_∂Ω·t( u,ϕ)  dS- ∫_∂Ω∂∂n·r(u,ϕ) dS - ∫_C^∂Ω·j(u,ϕ)   dl,
  
  where C^∂Ω is the union of all sharp edges of the domain, and the integral on it reduces to a sum evaluating at the boundary corners in 2D.
  
  Thus, applying boundary conditions (<ref>)-(<ref>), under the assumption Γ^_D_2=∅, and adding the weighted residual of the electric potential problem (<ref>) with (<ref>) and (<ref>), the weak form of  (<ref>) in  ℋ^2(Ω) is: find ∈ [ℋ^2(Ω)]^ and ϕ∈ℋ^1(Ω) such that (<ref>) and  (<ref>) hold and
  
      ∫_Ωε() : (,ϕ)  dΩ   + ∫_Ωε() ⋮ (,ϕ)  dΩ  + ∫_Ωψ·(,ϕ) dΩ
         = s()
          - ∫_Ωψ q  dΩ -∫_Γ^ϕ_Nψ  w_n dS
  
  for all ∈ [ℋ^2(Ω)]^ and ψ∈ℋ^1(Ω) such that =0 on Γ^_D_1 and ψ=0 on Γ^ϕ_D where
  
      s()= ∫_Γ^_N_1·t_n  dS+∫_Γ^_N_2∂∂n·r_n dS + ∫_C_N^∂Ω·j^ext dl  + ∫_Ω·b dΩ.
  
  
  This weak form is not suitable when considering 𝒞^0 FE approximations, but the same derivation can be applied in the interior of each element of the mesh, as detailed next.
  
  
  
  
  § 𝒞^0 INTERIOR PENALTY FINITE ELEMENT METHOD
  
  
  
  The domain Ω is now split in FEs {Ω_e}_e=1^, and a 𝒞^0 piece-wise polynomial approximation is considered. That is, the approximation space for the components of the displacement and for the potential is
  
      𝒱^h={v∈ℋ^1(Ω)  such that φ_e^-1(v|_Ω_e )∈𝒫^p  for  e=1,…,}⊄ℋ^2(Ω),
  
  where φ_e is the isoparametric transformation from the reference element to the physical element Ω_e, and 𝒫^p is the space of polynomials of degree less or equal to p for simplexes, and less or equal to  p  in each direction for quadrilaterals and hexahedra.  
  
  Since the approximation space in not in ℋ^2(Ω), we can not consider the weak form (<ref>). However, the approximation is ℋ^2(Ω_e); thus, considering 
  (<ref>) in each element we have
  
  
      ∫_Ω_e·b dΩ = ∫_Ω_eε() : (,ϕ)  dΩ + ∫_Ω_eε() ⋮ (,ϕ)  dΩ
         
          - ∫_∂Ω_e·t^e( u,ϕ)  dS - ∫_∂Ω_e∂∂n^e·r^e(u,ϕ) dS - ∫_C^∂Ω_e·j^e(u,ϕ)  dl,
  
  where C^∂Ω_e is the union of the edges (corners in 2D) of the element Ω_e and j^e is the line force (punctual force in 2D) on C^∂Ω_e; see a representation in figure <ref>a. The superscripts highlight that  the surface and line forces, and the normal vector, are from the element Ω_e.
  
  Summing  for all elements, and noting that  is continuous but  ∂/∂n is not, we get
  
      ∫_Ω·b dΩ =    ∫_Ωε() : (,ϕ)  dΩ + ∫_Ωε() ⋮ (,ϕ)  dΩ
         - ∫_·(t^L( u,ϕ)+t^R( u,ϕ))  dS
         -∫_( ∂^L ∂n^L·r^L(u,ϕ)+ ∂^R ∂n^R·r^R(u,ϕ)) dS
         -∑_k=1^∫_C_k·(∑_e∈ E(k)j^e(u,ϕ) ) dl
         -∫_∂Ω·t( u,ϕ)  dS 
          - ∫_∂Ω∂∂n·r(u,ϕ) dS,
  
  where  is the union of all the internal element faces and Ω is the union of the interior of the elements, where second derivatives are well-defined, i.e. 
  
      Ω= ⋃_e=1^Ω_e,   =⋃_e=1^∂Ω_e \∂Ω,
  
  see figure <ref>b.
  The supercripts R and L now denote the evaluation from the elements to the left and right side of the face in  (see figure <ref>c), and 
  {C_k}_k=1^ are all the  edges (corners in 2D) in the mesh, being  E(k) the set of indexes of the elements sharing the edge C_k.
  
  
  
  
  Now, let us recall that 
  the conditions for interfaces in the domain (also in the case of discontinuous material parameters) are the ones corresponding to both continuity for the Dirichlet values and equilibrium of Neumann forces. That is, 
   
  
      ⊗n = 0
         t (,ϕ)  = 0
         ∂∂n= 0
         r (,ϕ)⊗n  =  0
  
  
  on the faces in , where
  
  the jump operator is defined as
  
  
      a = a^L+a^R
  
  and it is used always involving a change of sign due to an odd appearance of the normal vector.
  
  In addition, we have to impose equilibrium of forces on the mesh edges.
  That is, each element Ω_e contributes with a force j^e on its edges (corners in 2D); and for each edge C_k, the sum of the forces for all elements sharing the edge (i.e. for Ω_e with e∈ E(k)) must be zero, or in internal equilibrium with the external forces. That is,
  
  
      ∑_e∈ E(k)j^e(u,ϕ) =  
      {[                0    on  C_k ⊄C^∂Ω;            j^ext on  C_k ⊂ C_N^∂Ω ].
   
  
  where j^ext is the force set in  (<ref>).
  Note that {C_k ⊄C^∂Ω} includes interior edges and also element edges on ∂Ω, just excluding the ones in the domain sharp edges.
  For the edges in Γ_D_1^ no value is set, and the sum of the forces will be in equilibrium with the reaction forces associated to the prescribed displacement (<ref>).
  
  On other hand, using the algebraic identity (a^Ln^L)b^L+(a^Rn^R)b^R = abn + anb, and the equilibrium condition (<ref>), we can rewrite 
  
      ∂^L ∂n^L·r^L(u,ϕ)   + ∂^R ∂n^R·r^R(u,ϕ) 
          
         
          = ∇ v : r(u,ϕ)⊗n + ∂v∂n·r(u,ϕ)
          
          = ∂v∂n·r(u,ϕ),
  
  with the mean operator {a} = (a^L+a^R)/2.
  
  Now, replacing in (<ref>) the identity (<ref>), the Neumann boundary conditions (<ref>) and (<ref>), the homogeneous Dirichlet condition =0 on Γ_D_1^ related to (<ref>), the first interface equilibrium condition (<ref>) and the equilibrium at interior edges (<ref>), and under the assumption Γ^_D_2=∅, (<ref>) simplifies to
  
      ∫_Ωε() :  (,ϕ)  dΩ   + ∫_Ωε() ⋮ (,ϕ)  dΩ
         
         -∫_∂v∂n·r(u,ϕ)dS
          = s(),
  
  with s() defined in (<ref>).
  
  
  The first two integrals in (<ref>) are symmetric and coercive bilinear forms in  and , as expected for the weak form of a strain gradient elasticity operator. However, it is not the case for the integral on the interior faces .
  
  The idea of IPM is adding terms that are analytically zero, thanks to the continuity interface condition (<ref>), to recover symmetry and coercivity of the strain gradient bilinear form. 
  
  The resulting weak form for flexoelectricity, under the assumption Γ^_D_2=∅, is: find ∈ [ℋ^1(Ω)∩ℋ^2(Ω)]^ and ϕ∈ℋ^1(Ω) such that (<ref>) and  (<ref>) hold and
  
      ∫_Ωε() : (,ϕ)    dΩ +∫_Ωε() ⋮ (,ϕ)  dΩ 
          + ∫_Ωψ·(,ϕ)  dΩ
         -∫_∂v∂n·r(u,ϕ)dS
           -∫_r(v,ψ)·∂u∂ndS
         +∫_β∂v∂n·∂u∂ndS
         = s()- ∫_Ωψ q  dΩ -∫_Γ^ϕ_Nψ  w_n dS
  
  for all ∈ [ℋ^1(Ω)∩ℋ^2(Ω)]^ and ψ∈ℋ^1(Ω) such that =0 on Γ^_D_1 and ψ=0 on Γ^ϕ_D.
  
  The parameter β is a stabilization parameter that must be taken large enough to ensure coercivity of the strain gradient bilinear form, to get a well-defined saddle point problem <cit.>. Although it is usually called penalty parameter, thanks to the consistency of the formulation, high-order convergence can be achieved with β or order h^-1. In practice, not very large values are needed for accurate solutions, avoiding the unaccuracy or ill-conditioning that typically suffer non-consistent penalty methods <cit.>.
  
  The minimum value of the stabilization parameter β can be estimated solving an eigenvalue problem, as commented in section <ref>.
  
   
  If second Dirichlet boundary conditions (<ref>) are imposed (i.e. Γ_D_2^≠∅), an additional term ∫_Γ_D_2^∂v/∂n·r(u,ϕ)  dS appears in (<ref>) and, consequently, in the weak form (<ref>). Following the same IPM rationale, two new terms, that are null thanks to (<ref>), are also added in (<ref>) to recover again symmetry and coercivity, namely 
  ∫_Γ_D_2^r(v,ψ) ·(∂u/∂n- g_2 )  dS 
  +
  β_D ∫_Γ_D_2^∂v/∂n·(∂u/∂n- g_2 )  dS, where β_D is a new stabilization parameter that can be taken equal to β or tuned separately.
  
  
  
  
  It is interesting to note that the C0-IPM weak form (<ref>) reduces to the one for ℋ^2(Ω), i.e.  (<ref>), when a  𝒞^1(Ω) approximation is considered. The C0-IPM formulation keeps the consistency and is valid for standard 𝒞^0 FE approximations, just introducing the proper integrals on the faces . Also note that the second integral in (<ref>) is in the interior of the elements, Ω, to account for the fact that second derivatives are not defined on . 
  
  
  
  
  
   §.§ Estimate of the interior penalty parameter β
  
  
  
  In this section we derive an eigenvalue problem to estimate a lower bound for β. The derivation is the usual one in IPM and Nitsche's formulations <cit.>.
  
  The bilinear form of the strain-gradient elasticity operator is 
  
  
      𝒜(,) = a(,) 
             -∫_∂v∂n·r^sg(u)dS
           -∫_r^sg(v)·∂u∂ndS
         
          + β∫_∂v∂n·∂u∂ndS
  
  
  with 
  
      a(,) = 
      ∫_Ωε() ⋮^sg ()  dΩ
  
  where r^sg_i=σ^sg_ijkn_jn_k
  and ^sg()=h ⋮ ∇ε(), that is, the mechanical part of the second traction and the double stress tensor.
  
  The bilinear form a is semicoercive (i.e. a(,)>0 for any  such that ε () ≠0, and a(,)=0 otherwise), leading to well-posed strain gradient elasticity  and flexoelectricity  problems for any value of the internal length scale parameter l. However, the addition of the integrals on the faces , leads to a bilinear form 𝒜 that retains semicoercivity only for large enough β. 
  
  Thus, to ensure well-posedness of the discrete problem for any value of l, we want β such that
  𝒜(,)>0   ∀ ∈𝒰^h={∈ [𝒱^h]^ such that ε () ≠0} with 
  
  
      𝒜(,) = a(,) 
          -2∫_∂u∂n·r^sg(u)dS
          +β∫_∂u∂n·∂u∂ndS.
  
  
  
  Using the Cauchy-Schwarz and Young's inequalities we can bound the interface terms as
  
      -2∫_   ∂u∂n·r^sg(u)dS
          +β∫_∂u∂n·∂u∂ndS
         ≥ - 2r^sg(u)∂u∂n 
         + β∂u∂n^2
       
         ≥   - 1/ϵr^sg(u) ^2
         + (β- ϵ) ∂u∂n^2,
  
  for any positive ϵ.
  
  Thus, considering a positive constant K such that
  
      r^sg(u) ^2 ≤ K a(,)
        ∀∈𝒰^h,
  
  we have
  
      𝒜(,) ≥(1-K/ϵ) a(,)+ (β- ϵ) ∂u∂n^2
  
  and the bilinear form is then positive definite if 1-K/ϵ>0 and β- ϵ>0, for any positive ϵ.
  
  In conclusion, the strain gradient  bilinear form 𝒜 is positive definite in the reduced discrete space if β>K, where K is the constant satisfying (<ref>). 
  This constant can be computed as the largest eigenvalue of the generalised problem
  
      Bx = λAx
  
  where B and A are the discrete matrices corresponding to 
  the bilinear forms b(,) and a(,) in the reduced discrete space 𝒰^h, with
  
      b(,)= ∫_r^sg(v)·r^sg(u)dS.
   
  
  The computation of the maximum eigenvalue in the reduced space 𝒰^h can be done from the problem stated in the complete discrete space [𝒱^h]^ setting nodal values to reduce the space or using the so-called eigenvalue problem deflation <cit.>.
  
  
  Matrices B and A scale as 𝒪(E^2 l^4 h^-5) and 𝒪(E l^2 h^-4), respectively, with characteristic element size h. Thus, the maximum eigenvalue of (<ref>) scales as 𝒪(E l^2/h). Consequently, we can consider
  
  
      β = α El^2/h,
  
  
  with a large enough constant α, that can be computed solving the eigenvalue problem, or simply tuned, in a coarse mesh with any value of E and l. 
  
  It is important noting that, differently to non-consistent penalty methods, IPM and Nitsche's methods provide accurate solutions and high-order convergence  with moderate values of β of order  𝒪(h^-1) for any degree of approximation, as shown in the numerical examples in section <ref> and in <cit.>. 
  
  
  An alternative sufficient condition to have a well-posed discrete problem can be stated  including also the elasticity term in the bilinear form, that is, with 
  
      a(,) = 
      ∫_Ωε() : ^sg() dΩ +
      ∫_Ωε() ⋮^sg ()  dΩ ,
  
  where ^sg()= ℂ : ε()  is the mechanical part of the local stress tensor.
  This option leads to a smaller (sharper) bound for β, specially for small l or large h. However, since the matrix corresponding to the first elasticity term scales as 𝒪(Eh^-2),  the dependency on the mesh size and material parameters is not so obvious.
  
  
  
   §.§ Implementation aspects
  
  
  
  The current implementation considers high-order Lagrange nodal basis, with Fekete nodes in the reference element to minimise the condition number of elemental matrices. For degree p≥ 3, special attention must be paid to the position of interior nodes in curved physical elements to keep high-order convergence, see <cit.>. High-order mesh generators, see for instance <cit.>, produce curved meshes taking care of this important aspect.
  
  The computation of the system involves two separated loops: in elements for volume integrals, and in faces for the computation of integrals on .
  To do so, the standard 𝒞^0 reference element is extended including second derivatives of the basis functions at the element integration points, the value of element basis functions and their derivatives at the integration points of the reference element faces, a list of the nodes corresponding to each face in the reference element and permutations for the integration points of the reference face for flipping. 
  
  The so-called flipping is a permutation (usually for the nodes in DG methods, but for integration points in our IPM implementation) that has to be applied to the face when seen from the second element, to match the orientation  of the corresponding face in the first element. In 2D the flipping is the same for any side of the mesh, just using a reverse ordering for the second element sharing the side. In 3D the possible rotations of the face have to be taken into account to choose the proper permutation for the integration points. 
  
  A variable storing, for each face, the number of the elements sharing the face, the local numbering of the face in each one of the two elements and the rotation to be applied for the second element, is also computed from the mesh as a preprocess.
  
  
  
  Dirichlet conditions (<ref>) are imposed in strong form, just setting the corresponding nodal values, as usual in standard FE computations. Second Dirichlet conditions can be imposed in weak form as commented in Remark <ref>.
  
  
  
  
  § NUMERICAL EXAMPLES
  
  
  
  Several numerical examples are included in this section to study the convergence of the C0-IPM formulation in 2D and 3D, and to validate the computational tool by comparison with previous works. 
  Homogeneous first, second and corner Neumann boundary conditions are assumed where no boundary condition is specified.
  
  
  
  
   §.§ 2D convergence test
  
  
  
  The convergence of the method for the solution of 
  problem (<ref>) is studied in this section. To test the method with non-regular meshes and curved boundaries, the problem is solved in a square with a hole, Ω=(0,1)^2\ B((0.5,0.5),0.2). Figure <ref> shows the coarsest mesh for nested refinement, with degree p=4.
  
  
  
  
  
  First Dirichlet  and second Neumann conditions are imposed on all the boundary. The body force b, the free charge q and the boundary data are set so that the solution is 
  
      [ u = [sin(2π(x_1+x_2)), cos(2π(x_1+x_2))]^T,;     ϕ = sin(2π(x_1+x_2))+ cos(2π(x_1+x_2)), ]
  
  
  the material parameters are
  
      [         E = 2.5,         ν =0.25,                 ;           l=1.1,      κ_L = 1.21,; e_L = 7.2
           , e_T =1.33
           ,      e_S = 1.73,;       μ_L = 1.5,      μ_T = 1.34,     μ_S = 5.47 , ]
  
  and the piezoelectric principal direction is x_1. The definition of the material tensors in terms of these parameters can be found, for instance, in appendix B of <cit.>.
  
  
  First we consider the uncoupled problem (with e=0 and μ=0), that is, an uncoupled solution of a strain gradient elasticity problem and an electric potential problem. The convergence plots are shown in figure <ref>, for penalty parameter β=100El^2/h, and degree p=1… 4. For strain gradient elasticity, the displacement error behaves in agreement with the results for Kirchhoff plates in <cit.>. 
  
  With degree p=1, the approximation space is not rich enough to impose 𝒞^1 continuity. Moreover, the second derivatives of the displacement in the strain gradient elasticity terms and the flexoelectricity terms cancel out, or are almost zero for curved elements. Thus, the method does not converge for linear approximation. For degree p= 2 much finer meshes would be necessary to reach assymptotic convergence, reducing its practical applicability. Accurate results with high-order convergence are obtained for degree p≥ 3, with slightly suboptimal convergence for p=3, in agreement with the analysis in <cit.> for the biharmonic equation. In this particular example, p=4 behaves better than expected, exhibiting slightly superoptimal convergence. The expected optimal convergence is observed for the uncoupled electric potential problem for any degree.  
  
  Figure <ref> shows the convergence plots for the flexoelectricity problem, with piezoelectric and flexoelectric coupling. The coupling leads to a reduction in the convergence rate, not relevant for the displacement, but around one for the potential, for p≥ 3. This is probably due to the, small but still present, discontinuity of the displacement derivative across element sides, affecting the potential through the flexoelectricity coupling.
  
  The conclusion is then that, even though convergence is suboptimal, the method is able to reach high accuracy with high-order convergence for degree p≥ 3. The C0-IPM method is therefore promising for an efficient solution of flexoelectricity.
  
  Similar results can be observed with quadrilateral meshes, with better behaviour for the p=2 approximation thanks to the richer approximation space and the presence of interior nodes in the element.
  
  
  
  
   §.§ Robustness with respect to the interior penalty parameter β 
  
  
  The effect of the interior penalty parameter in the accuracy of the numerical solution is studied next, with the 2D example and meshes of the previous section. Following Remark <ref>, the parameter is taken as (<ref>),
  with different orders of magnitude for α, independent of h.
  
  
  
  Figure <ref> shows the convergence plots for the flexoelectricity coupled problem, for the displacement  (left) and for the potential ϕ (right), for degree p=3 (top) and p=4 (bottom). The slopes of the segments are shown for the plots with α=10 for p=3, and with α=100 for p=4. We can observe the poor performance of the method for α=1, due to the fact that it is not large enough for a coercive mechanical bilinear form. 
  
  For degree p=3, α=10 is large enough and provides the best results. Larger values of α, several orders of magnitude larger, also lead to high-order convergence, proving the robustness of the method; but, in agreement with the analysis in <cit.> the convergence rate slowly decreases for increasing α.
  
  Looking to the results for p=4 we can observe that, for α=10, the bilinear form is coercive for the first meshes, because the elasticity part dominates in the coefficients of the matrix. This is not the case for the last mesh, where higher order terms become more relevant. With α≥ 100 the condition in Remark <ref> is satisfied and convergence is close to p+1=5 for  and around p=4 for ϕ, with almost no loss in the accuracy for increasing β.
  
  Thus, from this experiment we conclude that C0-IPM with degree p=4 provides excellent results, with convergence rates close to p+1=5 for  and around p=4 for ϕ,  and with little dependency on the particular value of β, for β≥ 100 El^2/h.
  
  
  
  The same analysis is performed now for strain gradient elasticity. Figure <ref> shows the convergence plots for the displacement  for degree p=3 (left) and degree p=4 (right), with the same conclusions.
  
  
  
  
   §.§ Cantilever beam
  
  
  The cantilever beam depicted in figure <ref> is considered. The aspect ratio is 20, and the width a varies to show the size-dependent nature of flexoelectricity.
  The beam is fixed to a wall and grounded on its left end, and it undergoes a punctual force F at the top-right corner. 
  
  
  
  The boundary conditions are thus 
  
      [            u=             0    at   x_1=0;      j_2(,ϕ)=            -F at  x=(L,a/2);            ϕ=             0   at   x_1=L, ]
  
  where L=20a is the beam length.
  
  To reproduce the results obtained in <cit.> with B-splines, the material parameters are
  
      [ E = 100 ,  
      	, κ_11 = κ_22 = 11  □,;             e_T = -4.4 □, μ_T = 1 ,;    l=ν = μ_L = μ_S = e_L = e_S = 0, ]
  
  and the piezoelectric principal direction is x_2. A uniform discretization with 2× 2 × 40 triangular elements (with characteristic element size h=0.5a) of degree p=4, and with β=100, is considered. Since l=0, any positive value of β provides good results.
  
  
  
  
  Figure <ref> left shows the normalised effective piezoelectric constant, e', versus the normalised beam thickness, a', defined as
  
      a'=-ae_Tμ_T^-1,   e':=k_effk_eff|_μ= 0,   	k_eff:=√(∫_ΩE·κ·E dΩ∫_Ωε·C·ε dΩ),
  
  where k_eff|_μ= 0 is the effective piezoelectric constant in the absence of flexoelectric effects, i.e. with μ=0.
  
  The results are in perfect agreement with the B-spline results in <cit.>, and with the analytical approximation in <cit.>:
  
      e'|_flexo(a')≃√(12a' ^2),   e'|_flexo-piezo(a')≃√(1+12a'^2).
  
  The plots in figure <ref>  also illustrate how flexoelectricity is a size dependent phenomenon, with relevant and even crucial effect for very small scales.  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Open and Closed circuit in the cantilever beam
  
  For further validation of the C0-IPM computational model, we now consider the open and closed circuit example in <cit.>, where maximum-entropy approximations (LME) were used. The problem is solved on the same beam with the same FE mesh.
  The material parameters are  now
  
  
      [            E = 100 , ν=0.37,; κ_11 = 11 □, κ_22 = 12.48 □,; e_T =-4.4 □, μ_T = μ_L = 1 ,;     l=  μ_S = e_S =  e_L= 0, ]
  
  
  and the the piezoelectric principal direction is again x_2.
  The mechanical boundary conditions are the same as in the previous case.
  
  For the electrostatic boundary conditions, two different cases are considered: open and closed circuit. The open circuit is the one considered in the previous example, with grounded right end, that is ϕ=0 at x_1=L, as shown in figure <ref>. In the closed circuit, the upper side is grounded and an electrode is placed on the bottom side, that is
  
      ϕ={[               0     for x_2=a/2;               V for  x_2=-a/2 , ].
  
  where V is a free constant value, see figure <ref>. The electrode condition is enforced setting all potential nodal values on the bottom boundary to be equal to the first one, with Lagrange multipliers in our implementation. 
  
  
  Figure <ref> shows the normalised effective piezoelectric constant e' as a function of the normalised thickness a'. Again, we observe that flexoelectricy becomes relevant for small scales.  For the open circuit, comparing to the previous results in figure <ref>, where ν=μ_L = 0, this more general model gives lower values for the normalised effective piezoelectric constant. On other hand, the open circuit setting leads to larger values of the effective piezoelectric constant. The numerical results are  in perfect agreement with the LME results in <cit.> demonstrating again the applicability of C0-IPM for the study and design of flexoelectric devices.
  
  
  
  
  
  
   §.§ Actuator example
  
  
  In this section we consider an actuator beam also from <cit.>. The displacement is fixed on the left boundary, and a potential difference is applied at the top and bottom sides, leading to a bending of the beam. That is,
  
      u(0,x_2)=0,  ϕ( x_1,a/2)=0,   ϕ( x_1,-a/2)=V,
  
  on the same beam, i.e. Ω=(0,20a)×(-a/2,a/2). 
  
  
  
  The material parameters are (<ref>) and the applied voltage is V=-8a.
  
  Figure <ref> shows the potential on the deformed beam  for width a=2.5. The potential seems to be smooth, but the section along x_2=0 in figure <ref> reveals a sharp variation close to the right end. Consequently, the electric field also presents sharp variations close to the right end, as shown in figure <ref>. 
  
  
  
  
  
  
  
  
  
  
  These results have  been computed on the adapted quadrilateral mesh in figure <ref> (top), with degree p=4 and β=1. The mesh has been refined to capture the sharp variations in the solution; otherwise, numerical oscillations spoil the solution in the whole domain. It is also worth mentioning that the plot in figure <ref> coincides in magnitude and shape with the results in <cit.> with LME, but getting rid of the smooth oscillations.
  
  Sharp variations along the boundary in the solution of flexoelectricity problems can be even more pronounced, as can be observed in figure <ref>. In this case the problem is solved with strain gradient elasticity, with l=0.1a,  on the p=4 adapted mesh in figure <ref> with min(h)=a/2^5; with smaller element size along the boundary to capture the high curvatures in the electric field. The stabilization parameter is again taken as β=100El^2/min(h), providing stable results.
  
  
  
  
   §.§ Periodicity
  
  
  The implementation of periodicity boundary conditions in the C0-IPM method is straightforward, by simply considering the periodicity faces as interior faces and imposing the periodicity constraints on the boundary nodal values. 
  
  Considering the periodicity faces as interior faces, that is in , ensures that 𝒞^1 continuity is enforced in weak form and that internal forces are equilibrated also on the periodicity boundary. The periodicity conditions for the nodal values can be implemented, for instance, by means of Lagrange multipliers, or reducing the system to the periodic space.
  
  
  
  As a verification example, figure <ref> shows the evolution of the error under nested refinement for the solution of the flexoelectricity coupled problem (<ref>) in a square domain Ω=(0,1)^2 with a regular triangular mesh. First Dirichlet and second Neumann conditions, (<ref>) and (<ref>),  are set on the top and bottom boundaries, and periodicity is imposed in the x_1 direction. That is, (<ref>) is imposed identifying the left and right boundary as the same boundary and including it in . The body force b, the free charge q, and the data for the boundary conditions on the top and bottom boundaries, are set so that the analytical solution is (<ref>). The stabilization parameter is (<ref>) with α=100.
  
  The errors exhibit the same behavior as in the convergence analysis in section <ref>.
  
  
  
  
   §.§ 3D convergence test
  
  
  The flexoelectricity equations are now solved in a cube, Ω=(0,0.5)^3, to show the applicability of the method also in 3D. The mesh for degree p=2 and the second level of refinement is shown in figure <ref>. First Dirichlet and second Neumann boundary conditions are considered in the whole boundary, and the material parameters are (<ref>). The data is set so that the solution is 
  
      =[ cos(2 π(x_1+2x_2-x_3)),sin(2 π(x_1+2x_2-x_3)),cos(2 π(x_1+2x_2-x_3))]^T, 
         ϕ=sin(2 π(x_1+2x_2-x_3)).
  
  
  
  
  
  
  
  
  Figures <ref> and <ref> show the convergence plots with β=100El^2/h, for strain gradient elasticity (solving the decoupled problem) and for flexoelectricity, respectively. As in 2D, the method does not converge for degree p=1; thus, we show the results for p=2,3,4. 
  
  Robust high-order convergence is observed in all cases, providing accurate results. Again, in agreement with the analysis in <cit.>, the convergence is suboptimal; but still with order close to p+1 for the displacement in the strain gradient elasticity problem for p≥ 3. Again, we also observe that the flexoelectricy coupling provokes a loss in the convergence rate and the accuracy of the solution; with order close to p for the displacement and the potential in this example.
  
  Further numerical experiments show that, with these regular hexahedra meshes, the error has very little dependency on the particular value of β≥ 100El^2/h. 
  
  
  
  § CONCLUSIONS
  
  A novel C0-IPM formulation for strain gradient elasticity and flexoelectricity is proposed. The weak form involves second derivatives of the displacement in the interior of the elements, plus integrals on the element faces, weakly imposing continuity of the displacement derivatives, as well as equilibrium of internal forces across element faces and on interior edges (vertexes in 2D). 
  
  The formulation is stable, with a symmetric and positive definite matrix for the strain gradient elasticity operator, for large enough interior parameter β. An eigenvalue problem is stated to determine a bound for β, which leads to a general formula for the parameter: β=α El^2/h, with constant α independent of the element size. Thus, differently to non-consistent penalty methods, and as usual in interior penalty methods, moderate values for β provide stable and accurate results.
  
  Standard 𝒞^0 FE approximations are considered, retaining the advantages and computational efficiency of high-order FE. The implementation is based on  assembly of elemental matrices, with standard FE numerical integration and nodal approximation, the discretization can be adapted to the geometry and locally refined where needed, no additional unknowns are needed, and material interfaces can be directly considered just adapting the mesh, as usual in FE computations.  
  
  The application of C0-IPM to problems with periodicity boundary conditions is straightforward, just considering the periodicity faces as interior faces (thus, imposing 𝒞^1 continuity and equilibrium of forces in weak form) and setting the periodicity conditions on the nodal values.
  
  Convergence tests, on 2D non-uniform curved triangular meshes and on 3D hexahedra regular meshes, show high-order convergence of the method for degree p≥ 3. A slow continuous loss in the convergence rate for increasing β is observed for p=3, which is in agreement with the analysis for the biharmonic equation in <cit.> and the results for Kirchhoff plates in <cit.>. Fortunately, for p=4 the convergence shows little dependency on β. In any case, in all examples, the convergence rates are at least close to p for both variables,  demonstrating the good behaviour of the method for p≥ 3.
  
  The computational tool is also validated by comparison with previous works solving realistic actuator and sensor problems on a beam, with perfect agreement.
  
  
  
  
  § ACKNOWLEDGEMENTS
  
  This work was supported by the European Research Council (StG-679451 to Irene Arias) and Generalitat de Catalunya (2017-SGR-1278)
  
  unsrt
  
  
  