Divergence‐free tangential finite element methods for incompressible flows on surfaces

Summary In this work we consider the numerical solution of incompressible flows on two‐dimensional manifolds. Whereas the compatibility demands of the velocity and the pressure spaces are known from the flat case one further has to deal with the approximation of a velocity field that lies only in the tangential space of the given geometry. Abandoning H 1‐conformity allows us to construct finite elements which are—due to an application of the Piola transformation—exactly tangential. To reintroduce continuity (in a weak sense) we make use of (hybrid) discontinuous Galerkin techniques. To further improve this approach, H(divΓ)‐conforming finite elements can be used to obtain exactly divergence‐free velocity solutions. We present several new finite element discretizations. On a number of numerical examples we examine and compare their qualitative properties and accuracy.


INTRODUCTION
Partial Differential Equations (PDEs) that are posed on curved surfaces play an important role in several applications in engineering, physics and mathematics.Surface PDEs describing flows on surfaces appear for instance in the modeling of emulsions, foams and biological membranes, cf.Slattery et al. 1 or Brenner 2 , or liquid crystals, cf.de Gennes and Prost 3 or Napoli and Vergori 4 .The numerical treatment of these PDEs gained an increasing amount of attention in the field of numerical simulations and numerical analysis in the last two decades.In this work we consider vector valued PDEs for viscous incompressible flows on surfaces that are immersed in the three dimensional space.
A main source of difficulty for vector valued PDEs on surfaces is the fact that the unknown vector field is typically tangential, i.e. for a two dimensional manifold Γ embedded in ℝ 3 the unknown field is only two-dimensional.In recent works almost exclusively [ 1 (Γ)] 3 -conforming finite elements have been used to approximate the unknown tangential vector field.Tangential vector fields are only weakly imposed through the variational formulation.In this work we follow a different approach: we abandon continuity of the finite elements.This loss of conformity however allows us to construct exactly tangential vector fields.This is achieved by mapping finite element functions from the two-dimensional reference element by a straight-forward generalization of the well-known Piola transformation.This guarantees that the resulting (possibly higher order) vectorial basis functions are exactly tangential to the surface.This specifically means that no additional enforcement of the tangentiality condition is needed to be enforced through the variational formulation.One could say that we trade one structure property (continuity) for the other (tangential vector fields).To deal with the missing continuity we apply well established techniques from the flat case: discontinuous Galerkin (DG) methods and variants such as the hybrid DG (HDG) methods.
It turns out that we do not have to abandon continuity completely, but can preserve continuity at least for the co-normal component, i.e. the in-plane normal component across element interfaces, resulting in (div Γ )-conforming finite elements.These finite elements in conjuction with suitable (hybrid) DG techniques have been proven to be excellent discretizations for incompressible fluid flows in the flat case due to benefitial properties such as exactly divergence-free solutions, pressure robustness and energy stability.These properties can easily be transfered to the case of surface Stokes and surface Navier-Stokes equations as we will explain in the sequel of this paper.

State of the art
Let us briefly give an overview on the state of the art in the literature.Initially surface finite element methods (SFEM) for scalar PDEs have been introduced in the seminal work by Dziuk 5 , we also refer to the survey papers by Dziuk and Elliott 6 and Bonito et al. 7 .The use and analysis of 1 (Γ)-conforming surface finite elements has been extended to higher order discretizations and adaptive schemes by Demlow et al. 8,9 , including the analysis of geometry errors 10 .In the last decade the extension to non-( 1 (Γ)-)conforming surface finite elements has been done in several works, cf. for instance Refs. 11,12and eventually also HDG formulations have been considered by Cockburn and Demlow 13 .In all these works an explicit surface mesh is used.A different approach is based on a mesh of the surrounding 3D space and a level set description of the surface.In the "TraceFEM" the trace on the level set surface of finite element functions of this background mesh are used for the approximation of the solution.The method was introduced by Olshanskii et al. 14 .Sometimes the method is also known as "CutFEM", cf.e.g.Ref. 15 .We also refer to the overview paper on TraceFEM by Olshanskii and Reusken 16 .
Mixed formulations of the surface Poisson problem can be seen as an intermediate step towards vector valued PDEs as they involve the vector valued surface flux that is approximated seperately from the primal unknown, resulting in a system of first order surface PDEs.These mixed formulations have been considered for instance by Rognes et al. 17 and -as part of their mixed formulations of DG and HDG methods -by Antonietti et al. 12 and Cockburn and Demlow 13 .Here, tangential surface finite elements are constructed for the flux.Primal DG formulations in the context of TraceFEM/CutFEM have also been considered in 18 .In the works by Rognes et al. 17 and Cockburn et al. 13 the construction of the spaces is based on the Piola transformation, resulting in (broken) surface Raviart-Thomas spaces as considered already by Nedelec 19 and Bendali 20 .Let us note that the analysis of mixed Poisson formulations including variational crimes in a general framework has been considered in Holst and Stern 10 .In all these works, where the primal unknown is scalar, an isoparametric geometry approximation, i.e. using order for approximating the scalar unknown and order for the geometry approximation is sufficient to obtain order + 1 error estimates in the 2 -norm.Only for the superconvergence property in the HDG method by Cockburn and Demlow 13 an increased geometry order is necessary.
While scalar problems on surfaces and their numerical treatment seem to be well understood vector valued problems have drawn an increased interest recently.Relevant models of viscous fluidic surfaces based on 3D Cartesian differential operators are described in Refs. 21,22.In the context of finite element methods these PDEs require vector valued finite element spaces on surfaces that are tangential.Starting with the Vector Laplacian surface finite elements that are [ 1 (Γ)] 3 -conforming, i.e. three dimensional, and impose the tangential condition through the variational formulation have been proposed by Hansbo et al. 23 .In that paper a penalty formulation and a Lagrange multiplier based formulation are considered to drive the normal component of the discrete approximation to zero.Further, it is already observed that -in contrast to the scalar problem -an isoparametric discretization is not sufficient to preserve optimal order 2 -errors.In Refs. 24,25similar approaches are considered and analysed in the context of TraceFEM discretizations.Hansbo et al. 26 extended their approach to Darcy problems on surfaces using [ 1 (Γ)] 3 -conforming (low order) surface FEM.The surface Stokes problem based on a velocity-pressure formulation has been discretized using TraceFEM based on a stabilized P1-P1 discretization in Ref. 27 .In the very recent paper by Olshanskii et al. 28 the TraceFEM with a P2-P1 discretization has been considered.Also only very recently Bonito et al. 29 presented a low order (div Γ )-conforming discretization for the surface Stokes problem on 4 smooth closed surfaces with a focus on the numerical treatment of killing fields.Their approach to discretize the surface Stokes problem is similar to the discretizations that we treat for the surface Navier-Stokes equations in this paper.Similar methods to Hansbo et al. 23 have been proposed in Nestler et al. 30,31 , where vector-and tensor-valued surface PDE models are considered.
A vorticity formulation has been considered for the surface Stokes problem in Refs. 32,33,34where the explicit construction of tangential vector fields is circumvented.Vorticity formulations have also been used to solve the surface Navier-Stokes equations in Refs. 35,36.Velocity-pressure formulations for the surface Navier-Stokes problem have been recently considered using higher order [ 1 (Γ)] 3 -conforming surface FEM by Reuther and Voigt 37 (low order) and by Fries 38 (higher order) and based on a low order TraceFEM discretization with penalty by Olshanskii and Yushutin 39 .We also mention the numerical approaches based on discrete exterior calculus in 33 and spectral methods in 40 .
In all the previous works either closed smooth surfaces or at least smooth surfaces with boundaries, e.g. in Fries 38 , have been considered.The case of only piecewise smooth geometries has not been addressed in the finite element literature so far to the best of our knowledge.
As we will base our discretization for the surface Navier-Stokes equations on (div Γ )-conforming elements, let us also briefly mention previous works on (div Γ )-conforming methods in the plane.In the context of DG discretizations Cockburn et al. 41 were the first to realize that energy stability and local mass conservation for DG methods is only achieved for (div Γ )-conforming finite elements which result in pointwise divergence-free solutions.We extended this idea to HDG methods and considered several extensions and improvements and evaluated the computational efficiency of the resulting methods in Refs. 42,43,44,45,46,47,48, cf. also the discussion in Section 3.5 below.

Content and structure
In this paper we introduce non-( 1 (Γ)-)conforming finite elements for incompressible surface flows, starting from the Vector Laplacian to the unsteady surface Navier-Stokes equations.We present different DG and HDG discretizations -most of which are new -and compare them to [ 1 (Γ)] 3 -conforming methods.The use of (div Γ )-conforming finite element methods results in pointwise divergence-free and exactly tangential solutions.Our methods are high order accurate but allow for surfaces which are only piecewise smooth.Moreover, arbitrary surface meshes can be dealed with, i.e. an exact geometry description does not need to be known.
Model problems are presented in Section 2 before several numerical schemes with tangential finite elements are introduced in Section 3. Several numerical examples for the different problems and discretizations are presented and discussed in Section 4 before we conclude the manuscript.

Notation and surface differential operators
Let Γ be a sufficiently smooth connected two-dimensional stationary and oriented surface embedded in ℝ 3 .At every point ∈ Γ we denote by ( ) a uniquely oriented unit normal vector, and by ( ) = − ( ) ( ) , with the identity matrix , the corresponding orthogonal projection onto the tangential plane of Γ at .In this work we assume that there exists a sufficiently smooth extension into the neighborhood (Γ) of Γ which induces a projection ∶ (Γ) → Γ. Assume a given scalar-valued, sufficiently smooth function ∶ Γ → ℝ, and let = • ∶ (Γ) → ℝ denote its extensions to the neighborhood (Γ) of Γ.Then, we define the scalar surface gradient through the gradient of in the embedding space and the projection onto the tangential plane.Hence, for any ∈ Γ we have where ∇ is the column vector consisting of all partial derivatives.As a direct consequence we realize that the scalar surface gradient lies in the tangential plane of Γ. Further, let us note that the surface gradient is independent of the concrete choice of the projection .In the same manner let = ( 1 , 2 , 3 ) ∶ Γ → ℝ 3 be a given vector-valued and sufficiently smooth function and denote by ∶ (Γ) → ℝ 3 its extension to the neighborhood.According to the above definition we can define the componentwise surface gradient through ∇ , where ∇ = (∇ 1 , ∇ 2 , ∇ 3 ) is the standard Jacobian matrix of .Hence, ∇ is the matrix where each row gives the scalar surface gradient of the components of .We define another operator called the tangential surface gradient by applying an additional projection from the left Note, that in the literature this operator is usually known as the covariant derivative on Γ.We are now able to state the symmetric surface strain tensor which is -following Gurtin and Murdoch 49 -defined as and the surface divergence operator div Γ ∶= tr(∇ Γ ).
As usual the divergence operator div Γ applied on a matrix valued function reads as the row-wise divergence.So far we assumed that and are sufficiently smooth so that the above differential operators exist in a point-wise sense on Γ.We can generalize to the notion of weak derivatives in the usual sense.For instance we define the weak gradient ∈ [ 2 (Γ)] 3 (if it exists) of a given function ∈ 2 (Γ), as the function that fulfills ∫ Γ ⋅ = − ∫ Γ div Γ for all ∈ [ ∞ 0 (Γ)] 3 .In the next three subsections we introduce the surface PDE problems considered in this work.To this end let ∈ [ 2 (Γ)] 3 such that ⋅ = 0 on Γ be a given force vector.

Vector-valued elliptic problem on Γ:
We seek for a solution ∶ Γ → ℝ 3 with ⋅ = 0 on Γ of the second order partial differential equation given by = 0 on Γ.
For the ease of presentation we only consider homogeneous Dirichlet boundary conditions in this part but the introduced methods can all be extended to more general boundary conditions, as demonstrated in the numerical examples.Further note that in the case of a closed surface ( Γ = ∅) no boundary conditions are prescribed.
Following Ref. 24 , the variational formulation of (2) is given by: Find ∈ such that where While finite element discretization of such variational problems are well understood on a flat surface, the tangential constraint ⋅ = 0 on Γ makes the construction of an appropriate numerical scheme on surfaces much more difficult.In Section 3 we discuss a natural approach how to deal with this challenge.
Remark 1.Note that the above variational formulation can also be defined on piecewise smooth connected surfaces.Then the strong form of the partial differential equation is given by equation (2) defined on each (smooth) sub domain and continuity conditions of the trace and the normal fluxes at the common interfaces.In Section 4.1.3we demonstrate that our methods are applicable for such problems.

Stationary Stokes equations on Γ:
We consider a Newtonian fluid on Γ, see Refs. 21,50, and assume for the ease of representation that Γ ≠ ∅, see Remark 2. Adding the incompressibility constraint and the pressure as Lagrange multiplier we now seek for a solution ∶ Γ → ℝ 3 with ⋅ = 0 on Γ and ∶ Γ → ℝ such that = 0 on Γ.
Here, (4a) can also be read as − div Γ ( ) = with the stress tensor ( , ) = − + 2 Γ ( ) where is the kinematic viscosity.Following the derivation in Ref. 38 the weak formulation is given by with ( , ) = ∫ Γ div Γ ( ) d , and the pressure space ∶= { ∈ 2 (Γ) ∶ ∫ Γ = 0}.Beside the crucial constraint that is in the tangential space of Γ we now further have to deal with the divergence constraint.In particular this plays a key role in the discretization as one has to find compatible function spaces for the discrete velocity and the pressure spaces.
Remark 2. As the variational formulation of the standard stationary (flat) Stokes equations is defined without a mass bilinear form ( , ) it demands further constraints to filter out the possibly non-trivial kernel of the composite bilinear form involving (⋅, ⋅) and (⋅, ⋅).On surfaces with a boundary sufficient constraints can be obtained from suitable boundary conditions.However, for the surface Stokes equations on closed surfaces the non-trivial kernel, the so-called killing fields, need to be taken care of to guarantee uniqueness.For a more detailed discussion we refer to Refs. 29,37.

Unsteady Navier-Stokes equations on Γ:
We extend the Stokes model to the fully unsteady Navier-Stokes case, i.e. we seek for a solution ∶ Γ × (0, ] → ℝ 3 with ⋅ = 0 on Γ and ∶ Γ × (0, ] → ℝ such that with > 0, a given boundary differential operator , boundary values and initial values 0 .Beside the difficulties already discussed in the Stokes setting, the challenging part now is to deal with the nonlinear convection and the time derivative.

CONSTRUCTION OF TANGENTIAL FINITE ELEMENT METHODS
In this section we focus on the derivation of new tangential DG schemes for the problems introduced in Section 2. After introducing some basic notation in Section 3.1, we concentrate on rather standard DG versions first in Sections 3.2-3.4.Benefits of choosing a (div Γ )-conforming formulations are discussed in Section 3.5.In Section 3.6 we introduce similar HDG formulations and explain their computational advantage over the standard DG formulations.

Basic ingredients and motivation
In this section we discuss a natural approach for the construction of tangential vector fields.To this end we summarize the findings of the works 17,51 .Let  ℎ be a consistent triangulation of the smooth manifold Γ into curved triangles such that for every element ∈  ℎ there exists a not degenerated polynomial mapping Φ of order from the reference element to the physical element , i.e.Φ ∶ ̂ → .With respect to this triangulation we write Γ ℎ ∶= ∪ ∈ ℎ for the corresponding locally smooth discrete manifold.In the following we will use the notation Γ , and ∇ Γ also for operations w.r.t. to Γ ℎ (instead of Γ).Further, we define the set  ℎ as the union of all element interfaces.We assume that  ℎ is shape regular and quasi uniform, thus there exists a mesh size ℎ with ℎ ≃ diam( ) for all ∈  ℎ .We denote by = Φ ′ ∈ ℝ 3×2 the Jacobian of the finite element mapping and remind the reader that the columns of span the tangential space for each point in .In the following, we will drop the subscript in the transformation Φ and the derived quantities such as the Jacobian unless the association to the element shall be highlighted.Next, we write for the Moore-Penrose pseudo inverse of the Jacobian , and set ∶= √ det( ) as the functional determinant.Now let ̂ be a differentiable function defined on the reference element ̂ and let = Φ( ̂ ) ∈ for all points ̂ ∈ ̂ .Using the classical pull back we define a function on by Following the ideas of Rognes et al. 17 , this pullback allows us to calculate the surface gradient of the function by As − = ( ) − , the above equation shows that the gradient ∇ is a linear combination of the two tangent vectors 2 and hence lies in the tangent space as expected from differential geometry, see for example 6 .The above construction allows us to define an 1 -conforming finite element space of order on the surface Γ ℎ by where ℙ ( ̂ ) is the space of polynomials up to degree on ̂ .Whereas the above technique allows an easy construction of a scalar approximation space, the problems ( 2), ( 4) and ( 5) demand for vector valued approximation spaces.In particular we aim for a space that can handle the constraint ⋅ = 0 in a proper way.Obviously, the simple product space ℎ × ℎ × ℎ is not convenient and we need a different construction.The solution to this is given by using the Piola transformation instead of the classical pull back.Originally, thus on flat surfaces, this mapping is used for the construction of (div)-conforming finite element spaces as it preserves the normal components on element interfaces.To this end let ̂ be a vector valued function on ̂ , then we define on the function Due to the multiplication with , we again see that the constructed vector field lies in the tangential plane of .This finding is the key argument and motivation for the construction of the numerical schemes in this work.As we will explain below, cf.Lemma 1, the factor 1∕ is important for the construction of (div Γ )-conforming finite element spaces, i.e. finite element spaces with continuous in-plane normal components.We want to mention that in Ref. 17 and Ref. 51 the authors already realized that the Piola mapping can be used on surfaces to construct appropriate finite element spaces for the approximation of the spaces (div Γ ) and (curl Γ ) on Γ ℎ .There, the according finite elements on the surface triangulation  ℎ are based on the Raviart-Thomas and Brezzi-Douglas-Marini finite element families on the reference element ̂ , cf.Refs. 52,53,54.

A discontinuous Galerkin discretization for vector valued elliptic problems
Based on the findings from the last section we now motivate a new discontinuous Galerkin (DG) method for the approximation of problem (2).However, we want to mention that the discretization can also be adjusted to other elliptic problems like the Vector Laplacian without a mass term.To this end we need several operators motivated by their definitions on a flat surface, see Arnold et al. 55 .Let 1 , 2 ∈  ℎ be two arbitrary elements with common edge = 1 ∩ 2 .Let ℎ be the oriented unit normal vector on Γ ℎ , and let be a uniquely oriented tangential vector on such that 1 , 2 are located on the left and on the right side respectively with respect to the direction of , see Figure 1.Using these vectors we define the in-plane unit (outer) normal vectors Let us stress here that ℎ will in general be discontinuous across element interfaces so that 1 and 2 will not be parallel.
Let and be vector and matrix valued functions respectively.We define the vector valued average and jump on by with the scalar valued normal and tangential jump and the averaged in-plane normal In the case of ⊂ Γ ≠ ∅ we set [[ ]] = ⋅ and the operators { {⋅} } and [[⋅]] are replaced by the identity.Based on the Piola mapping, see equation ( 6), we define a finite element space of order by

FIGURE 1
Two elements 1 , 2 ∈  ℎ with a common edge = 1 ∩ 2 and the corresponding in-plane normal vectors 1 , 2 , the common tangential vector and the oriented normal vector ℎ .On the left with a linear geometry approximation, on the right side with a high order approximation.
Here, ℙ ( ̂ ) is the scalar-valued polynomial space of order on ̂ .By construction, we have that thus discrete functions in ℎ are exactly in the tangent plane of Γ ℎ .Assuming that , originally defined on Γ, has a smooth extension on Γ ℎ , we follow Arnold et al. 55 to define the DG method (based on an symmetric interior penalty formulation): Find where with the constant chosen big enough, see Ref. 55 , and ℎ ( ℎ , ℎ ) = ∫ Γ ℎ ℎ ℎ d .Note that the projection in the definition of the bilinear form ℎ (⋅, ⋅) is redundant for functions in ℎ .However, we will use ℎ (⋅, ⋅) for non-tangential functions as well later on.
Remark 3. In the case of nonhomogeneous Dirichlet or other types of boundary conditions one uses the standard techniques employed for DG methods on flat surfaces, see Ref. 55 .
Remark 4. The application of the Piola mapping in the definition of functions in ℎ results in exactly tangential fields.Furthermore, it results in a discretization that is -as the continuous problem -invariant under isometries, i.e. if ℎ solves (9b) on Γ ℎ , ( ℎ ) solves the corresponding problem on Γ ′ ℎ = Φ(Γ ℎ ) if Φ is an isometry and  the corresponding Piola mapping.We discuss this in more detail below in the numerical examples, cf.Sections 4.1.3and4.3.1 below.

Discretization of the surface Stokes equations
In this section we focus on the construction of a numerical scheme to solve the Stokes problem on Γ, see equation ( 4).Note that we assumed Γ ≠ ∅ and thus also Γ ℎ ≠ ∅.As known from the literature, see for example Refs. 54,56, discrete stability demands to find a compatible pair of the discrete velocity and pressure space.We aim to construct a method based on the works 57,41,43 for the flat case.In the latter works, the discrete velocity space is chosen to be the (div)-conforming BDM space of order , and the pressure is approximated by discontinuous polynomials of order − 1.These two spaces fulfill the property that the divergence of the velocity space is a subspace of the discrete pressure space.This ensures stability in the sense of Babuška-Brezzi, see for example Ref. 44 , and leads to major benefits such as exactly divergence-free discrete velocities, see also Section 3.5.Following the ideas of Rognes et al. 17 and Castro et al. 51 , we now map this method to the surface using the previously introduced Piola mapping, see equation (6), to define an (div Γ )-conforming velocity space on Γ ℎ .Then for an arbitrary order we set where 0 (div Γ ) is the subspace of (div Γ ) with zero normal trace at the boundary Γ.Regarding the implementation of ℎ note, that the Piola mapping not only helps to incorporate property (8), but can also be used to incorporate the zero normal jump.This is done by choosing the BDM basis on the reference element ̂ and map it with the Piola transformation to the physical elements ∈  ℎ .As the BDM basis functions are associated to the scalar normal moments on the edges of the reference element, the mapped functions are then associated to the according (in-plane) normal moments on the edges of the mapped element .This automatically results in a normal continuous function, i.e. [[ ]] = 0 on all edges.A detailed discussion on this topic is given in 51 .Next we define the discrete pressure space as For a fixed velocity approximation order the (div Γ )-conforming DG method then read as: with Note that the appearing jumps in ℎ ( ℎ , ℎ ) now only read as the tangential jump since functions in ℎ are normal continuous.
As a consequence we have the following lemma.
where div( ̂ ℎ, ) reads as the divergence on the flat reference element which shows that div( ̂ ℎ, ) = 0, thus we conclude by equation 12.
From Lemma 1 we conclude that the solution ℎ of ( 11) is exactly divergence-free.Note that this statement is completely independent of the geometry approximation.This further leads to pressure robustness (see Sections 3.5 and 4.2) and shows that the function spaces ℎ and −1 ℎ are compatible, thus the stability proof of (11) follows the same steps as in the flat case, see Refs. 42,43.

Discretization of the surface Navier-Stokes equations
The discretization of the Navier-Stokes problem (5) follows the ideas of Refs. 43,42which we briefly summarize in the following.We aim to use a semi-discrete method to decouple the discretization in space and time.This is then further combined with an (high order) operator splitting method to efficiently deal with the nonlinear convection term.For the latter one we use a standard upwinding scheme in space which guarantees energy stability.To this end let ℎ , ℎ , ℎ ∈ ℎ , then we define the form where the upwinding value is chosen in upstream direction with respect to the convection velocity ℎ : Let ∈ and ′ be the neighboring element, s.t.∈ , ′ , then we define We notice that the upwinding only affects the discontinuous tangential component.This type of formulation is typically derived by applying partial integration in an element-by-element fashion and introducing numerical fluxes, cf.e.g.Ref. 58 .Let us further note that the integration by parts formula for vector fields that are not exactly tangential involves an additional term including the mean curvature, cf.Ref. 38 Eq. ( 3).
In combination with the bilinear forms defined in the previous section and again assuming that the initial velocity 0 has a smooth extension onto Γ ℎ , we derive the following spatially discrete DAE problem: find where (14c) resembles the 2 -projection of the initial velocity onto the discrete velocity space ℎ .Now let ℎ ( ) and ℎ ( ) be the finite element coefficient vectors of the functions ℎ , ℎ respectively, and let with = , … , where 0 = 0, = be an equidistant mesh for [0, ] with time step size Δ .To get the solution at time +1 we solve a step of the fully discrete lowest order implicit explicit (IMEX) splitting scheme Here the matrices , , are the corresponding matrices of the bilinear forms ℎ , ℎ and ℎ respectively, and ( ℎ ( )) ℎ ( ) reads as the evaluation of the convection trilinear form ℎ with the wind ℎ = ℎ ( ).Above system shows that we treat the stiffness and divergence constraint implicitly which guarantees exactly divergence-free velocity solutions at each point in time.The convection however is treated explicitly.For more details regarding the efficiency and accuracy of above splitting methods for the Navier-Stokes equations we refer to Refs. 43,42for the flat case.The simplest variant to generalize (15) to higher order order in time are IMEX schemes, cf. the works by Ascher et al. 59,60 .Below in the numerical examples we use a second order IMEX based on two compatible explicit and implicit Runge-Kutta schemes of second order.

On the benefits of (div Γ )-conformity
Additionally to the fact that the (div Γ )-conforming space ℎ is tangential there are several advantages over other Stokes discretizations (see also the HDG Stokes discretization introduced in Section 4.2).Below, we focus on two important ones.For other aspects -which transfer directly from the flat case -such as the ability to reduce the set of basis functions, convection stability (beyond energy stability) and good dissipation properties we refer to the literature, e.g.Refs. 42,43,61.

Pressure robustness
As discussed in Section 3.3, Lemma 1 shows that solutions ℎ of ( 11) are exactly divergence-free independently of the geometry approximation.Further, the combination of the velocity space ℎ and the pressure space −1 ℎ allows to test equation (11)  with exactly divergence-free velocity test functions.This is a crucial advantage as it enables us to derive velocity error estimates that are independent of the pressure approximation (and the viscosity ).This property is known in the literature as pressure robustness and was first introduced by Linke 62 .In the following we briefly sketch the main idea (hence the occurring problem): let the r.h.s. of the Stokes problem be decomposed as 3 such that ⋅ = 0. Testing the continuous problem (4) with a test function ∈ 0 ∶= { ∈ (div Γ , Γ) ∶ div Γ ( ) = 0} we see that 2 ( , ) = ∫ Γ ⋅ d , hence the velocity is not steered by the gradient ∇ Γ .If the same observation can be made in the discrete setting the method is called pressure robust and one can derive velocity error estimates that are independent of the pressure approximation (see Section 4.2 for a numerical observation of this phenomenon).A rigorous analysis of such (Helmholtz) decompositions on smooth manifolds is given in Ref. 34 .For the above findings in the continuous setting it was crucial that one tests with exactly divergence-free test functions ∈ 0 .Lemma 1 shows that the same can be done in the discrete setting which leads to pressure robustness of the discretization given by equation (11).Note, that also standard methods can be made pressure robust, see Ref. 63,64,65,66 .Further, it also plays an important role in the Navier-Stokes case 67,68,69,70 .

Energy stability
The standard DG upwind formulation of the convection bilinear form (13) allows to show the following stability result (up to geometry errors).Let Γ ∶= { ∈ Γ ∶ ℎ ⋅ < 0} then there holds Note that the crucial assumption here is that the transport velocity ℎ is exactly (surface) divergence-free.Considering the time discretization (15) of the Navier-Stokes equations, the transport velocity was chosen as the discrete velocity of the old time step, i.e. ℎ = ℎ ( ).Thus, as div Γ ( ℎ ( )) = 0 (exactly) we can apply above stability estimate showing energy stability of the Navier-Stokes discretization.

Hybrid discontinuous Galerkin formulations
Although the formulation of Section 3.2 fulfills property (8), the computational overhead introduced by the DG formulation is a major drawback.A well known technique to overcome this circumstance is to use hybrid discontinuous Galerkin (HDG) formulations and a static condensation strategy.The idea is to introduce additional unknowns on the skeleton  ℎ which circumvents the direct coupling of element (interior) unknowns with neighboring elements.We give two examples how this can be achieved and comment on further efficiency improvements.Below, we will use HDG discretizations for the numerical examples.

An HDG method for the Vector Laplacian
Let Ψ be the unique polynomial map from the reference edge ̂ = [0, 1] to an edge ∈  ℎ , with = 1 ∩ 2 , 1 , 2 ∈  ℎ as before.We define the space of (discontinuous) piecewise polynomials on the skeleton as Here, ℙ ( ̂ ) is the scalar-valued polynomial space of order on ̂ .Then we define the following two vector valued functions 1 ∶= + 1 and 2 ∶= + 2 for , ∈ Λ ℎ , where and 1 , 2 are defined as above.Here is introduced to decouple the weak (in the DG sense) tangential continuity whereas is introduced for the weak normal continuity.Note, that 1 might not be equal to 2 as the in-plane normal vector may be different.The HDG method then reads as: Find ℎ , where The HDG method has more unknowns compared to the DG version, but the element unknowns ℎ ∈ ℎ depend only on the facet function .Hence, in a linear system, we can form the Schur complement with respect to those unknowns which is known as static condensation.This can be done in an element-by-element fashion already during the setup of the linear systems.The Schur complement system is typically much smaller, especially for higher order elements (note that the facet unknowns scale with whereas the element unknowns scale with 2 ).An example for the impact of static condensation with HDG methods in comparison to DG methods will be given in Section 4.1.2.

An (div Γ )-conforming HDG method for the Stokes and the Navier-Stokes equations
Similarly to the standard DG method of Section 3.2, we can also use a hybridization technique for the (div Γ )-conforming DG methods of Sections 3.3 and 3.4.As normal continuity is already incorporated in the space ℎ , we do not need the facet variable . Hence the facet variable on the skeleton  ℎ reduces to the single valued quantity as the tangential vector is the same on both sides.With the definitions ( , ℎ ) ∶= + ( ℎ ⋅ ) and ( , ℎ ) ∶= + ( ℎ ⋅ ) , we can use the same HDG bilinear form HDG ℎ for the treatment of viscosity.For the convection operator we have to replace the upwind flux up ℎ , cf. (13), to avoid communication between interior element unknowns in ℎ .Hence, we make the following choice: Further, to make sure that takes the value of the upwind neighbor also in the convection dominated regime we add a consistent stabilization (the "downwind glue") as introduced in Ref. 71 and define where out denotes the outflow boundary of an element, i.e. out ∶= ∩ { ℎ ⋅ > 0}.Let us stress that in the convective limit → 0, the HDG solution converges to the DG solution discussed above.
The semi-discretization of the new hybrid (div Γ )-conforming DG method for the Navier-Stokes equations then reads: Find where we implicitly made use of = ( , ℎ ) and = ( , ℎ ).Using this spatial discretization we can then again use the implicit-explicit splitting scheme for the discretization of the Navier-Stokes equations on Γ ℎ , see equation (15).

Efficiency aspects and superconvergence
One advantage of HDG methods that we did not address so far is the ability to achieve superconvergence in diffusion dominated problems.If the (trace) unknowns on the skeleton are approximated with polynomials of order , one can (for example) apply a local postprocessing strategy (see Ref. 72 ) to define a local element-wise approximation of order + 1 which has order + 2 accuracy in the 2 norm.Alternatively, one can consider the previous HDG formulation and reduce the facet degree by one order while preserving the order of accuracy by introducing a slight modification in the formulation, cf.Ref. 43 Section 2.2.In the context of (div)-conforming methods, the reduction of the polynomial degree for the normal component requires a bit more care in order to preserve the benefitial properties of (div)-conforming methods, cf.Ref. 45,46 .

NUMERICAL EXAMPLES
In this section we consider several numerical examples.We start with a comparison of [ 1 (Γ)] 3 -conforming and the previously introduced non-conforming finite element methods for the Vector Laplacian for smooth and piecewise smooth manifolds in Subsection 4.1.In the subsequent Subsection 4.2 we consider and compare the non-conforming methods for a surface Stokes problem.In the remainder we fix the method to the (div Γ )-conforming HDG method and consider several surface versions of a well-known benchmark problem in 2D in Subsection 4.3 and the Kelvin-Helmholtz instability problem on different rotationally symmetric surfaces in Subsection 4.4.Finally, in Subsection 4.5 we consider a self-organisation process on the Stanford bunny geometry.All numerical examples were implemented within the finite element library Netgen/NGSolve, see Refs. 73,74and www.ngsolve.org.The data, i.e. time series, mesh refinement series, etc. as well as scripts for reproduction for all numerical examples can be found at DOI: 10.5281/zenodo.3406173.‖ FIGURE 2 Exact and discrete geometry and solution for the Vector Laplacian in Section 4.1.2.Coloring corresponds to velocity magnitude.From left to right: exact geometry and solution, (uncurved) computational mesh on coarsest level = 0, mesh and discrete solution for Hdiv-HDG with = 1, = 1, and discrete solution for Hdiv-HDG with = 1, = 2.

Vector Laplacian
First, we consider two Vector Laplace problems with different discretizations, compare and discuss them.Additionally to those discretizations introduced before in Section 3.2, we consider two [ 1 (Γ ℎ )] 3 -conforming discretizations from the literature that we briefly summarize in Section 4.1.1below and which we denote as H1-L (weak enforcement of tangential condition through Lagrangian multipliers) and H1-P (weak enforcement of tangential condition through penalties).We introduce labels for the aforementioned methods.The method in (9) will be denoted as DG whereas the hybrid version of it is denoted as HDG.For the (div, Γ ℎ )-conforming discretization, i.e. using ( 9) with ℎ replaced by ℎ , cf. ( 10), we use the label Hdiv-DG and denote the HDG version as Hdiv-HDG.For the numerical computations we only consider the HDG versions.Note however, that the differences between DG and HDG are primarily in a computational aspect, see also the explanations in Section 3.6.Hence, the DG versions are only considered for a comparison of this observation.

1 -conforming discretizations for the Vector Laplacian
In the literature mostly [ 1 (Γ ℎ )] 3 -conforming surface FEM discretizations for vector-valued problems are considered by either applying Lagrangian multipliers for the tangential constraint or using a penalty formulation, see e.g.Refs. 38,24,25,26,23,37.As we will use two instances of these methods for comparison to the methods discussed in Section 3, we briefly state these two formulations and discuss the choices involved.The methods are defined as follows: The [ 1 (Γ ℎ )] 3 -conforming method H1-L in (17a) uses Langrange multipliers to enforce the tangential constraint whereas (17b) enforces the tangential constraint through penalization with a penalty parameter (H1-P).For H1-L the choice of the Lagrange multiplier space is = as in Ref. 23 or Ref. 25 (in a TraceFEM context).We note that in Ref. 38 an ℎ −6 scaling of the condition number for Navier-Stokes problems has been observed for = and = − 1 has been used instead.To drive the normal component sufficiently small in the H1-P method, we make the simple choice = 10 ⋅ ℎ −( +1) and accept the severe ill-conditioning of arising linear systems here.In Refs. 23,25choosing such a large penalty is circumvented by replacing Γ ( ℎ ) (and Γ ( ℎ )) by Γ ( ℎ ) (and Γ ( ℎ )) which is a consistent manipulation.Thereby the normal and tangential parts of the discrete solution decouple and a soft penalty with a penalty parameter ∼ ℎ −2 suffices.Both choices taken here, = for H1-L and = 10 ⋅ ℎ −( +1) for H1-P will eventually result in severe ill-conditioning of arising linear systems.However, for the test cases under consideration with H1-L and H1-P we used sparse direct solvers and obtained results of sufficient meaningfulness to serve as candidates for comparison.

Vector Laplacian on the sphere
We consider the Vector Laplace problem with the (extended) exact solution = (− 2 , , ) on the unit sphere Γ = 1 (0), i.e. we choose the r.h.s.so that = | Γ solves (2).For the geometry approximation we considered two choices, = and = + 1 and for in the SIP formulations we take = 10.On five successively and uniformly refined meshes we evaluate 2 -and 1 -errors.The initial mesh (mesh level = 0) consists of 124 triangles, see Fig. 2. accurate geometry approximation is not surprising and has already been discussed in e.g.Refs. 23,25.In Fig. 2 a comparison of the pictures on the right illustrates the difference between = 1 and = 2 for = 1.In the rows two to four of Fig. 3 we fix = + 1, i.e. a superparametric geometry approximation, and choose ∈ {1, 2, 4}.For all methods we now observe optimal order of convergence, i.e. (ℎ ) in the 1 -semi-norms and (ℎ +1 ) in the 2 -norms.For = 4, the results of H1-P are not behaving optimally on the finest meshes.This is probably due to the ill-conditioning stemming from the penalty parameter of magnitude 10ℎ −5 .Let us further mention that we tried (for all methods) = and obtained essentially the same results, i.e. in this example it seems that the superparametric approximation is not necessary as long as > 1.
Overall, we can conclude that all methods perform similarly well for this example.The (H)DG methods are obviouslyby construction -perfect in approximating the tangential constraint, but for the other error measures there is no significant difference.

Computational aspects
In this paragraph we want to discuss and compare computational aspects of the different discretizaton methods under consideration for the Vector Laplacian.Starting with the 1 -conforming methods, one immediately notices that the advantage of these methods is their simplicity.Lagrangian finite elements are available in most finite element packages and hence, a realization of these methods is comparably simple.Then again, the two-dimensional vector field is approximated with a three-dimensional vector field which can be seen as undesirable in terms of the computational overhead.To overcome this issue we introduced tangential finite elements.These however came at the price of abandoning 1 -conformity which requires the implementation of weak continuity (at least in tangential direction) through the discrete variational formulation.This approach results in DG methods which come at the disadvantage of introducing more unknowns and more couplings.To alleviate these costs we also discussed the use of hybrid versions of DG methods.In Table 1 we compare the six different methods introduced before for polynomial orders 1 to 5 on the finest mesh = 4 of the previous example (Vector Laplacian on the sphere).While the 1conforming and the HDG methods are exactly those investigated in the previous paragraph, the methods DG and Hdiv-DG are only considered with respect to their computational costs, here.For the HDG methods and the 1 -conforming methods we apply static condensation in an element-by-element fashion.As measures for the computational costs we take the number of degrees of freedom (dof), the number of global dof that remain in the Schur complement after static condensation (gdof) and the number of non-zero entries (nze) in the Schur complement.
Let us first take a look at the dof measure.Here, the H1-P method has the smallest number of dof in the low order case ≤ 3. The additional dof for approximating a three-dimensional vector field are still less than those obtained from approximating a two-dimensional vector field with discontinuous piecewise polynomials.Only for ≥ 4 the Hdiv-DG method with normalcontinuity results in less dof.When normal-continuity is also broken up, i.e. when going to DG it requires at least order ≥ 6 to beat H1-P in terms of unknowns.When further going to HDG methods there is obviously no advantage over DG methods in terms of dof as only additional unknowns are introduced.As we can apply static condensation with HDG methods it is worth taking a look at those dof that remain in the Schur complement.The number of gdof can be reduced for DG for all polynomial degrees whereas the step from Hdiv-DG to Hdiv-HDG pays off in terms of gdof for ≥ 3. We observe that gdof is the same for HDG and Hdiv-HDG.When considering nze the picture shifts even further towards HDG and Hdiv-HDG.For = 2 already the methods HDG and Hdiv-HDG outperform the other two DG methods and for ≥ 3 the method generates even less unknowns than the 1 -conforming methods.We conclude that the hybrid formulations with tangential fields are benefitial not only because of their additional structure properties, but also computationally advantageous when going for higher order discretizations.Remark 5. HDG superconvergence.Let us note that we did not consider further tweaks of the HDG methods related to the aspect of superconvergence, cf. the paragraph on efficiency aspects and superconvergence in Section 3.6.With a corresponding modification the polynomial degree on the facets can be reduced by one order while keeping the same order of accuracy (for diffusion dominated problems), i.e. the costs in terms of gdof and nze for a method with accuracy (in the 1norm) are the same of those HDG methods without this modification with order − 1.In the lowest order case we obtatin gdof(HDG)=gdof(Hdiv-HDG)=96.8K and nze(HDG)=nze(Hdiv-HDG)=967.7K and hence, the HDG schemes are level with H1-P in terms of gdof, but already more efficient in terms of nze.

Vector Laplacian in the plane and on a house of cards
In this example we consider two configurations.First, we consider a flat domain case with Γ ∶= (0, 2) × (0, 1) where we pose the Vector-Laplacian problem (2) and replace (2b) with inhomogeneous Dirichlet boundary conditions = on Γ.As a second example we fold the geometry in the middle and lift it up to a house of two cards: where is the height and is the width of one card with 2 + 2 = 1 so that ∶= Φ ′ ∈ ℝ 3×2 is a length-preserving map with = ∈ ℝ 2×2 , cf. Figure 4 for a sketch.We define with ̂ 1 = 1 = (1, 0) and ̂ 2 = 2 = (0, 1) the tangent unit vectors of Γ.Then, = ̂ , = 1, 2 are the tangent unit vectors for Γ a.e., i.e. we can interprete -and due to = 1 also the Piola transformation -as the operator that realizes a basis transformation from one tangent space to another.
On Γ we pose the Vector Laplacian with r.h.s.̂ and boundary data ̂ such that the exact solution to (2) on Γ is For the same problem on Γ (with properly transformed data = ̂ Next, we discuss how this characterization of the solution translates to the discrete level.We note that Φ is piecewise affine, so that the Jacobian is piecewise constant.Let ̂ ℎ be the discrete solution to the Vector Laplace problem on Γ and ℎ =  ̂ ℎ .Then, we can expand .
and obtain the following relation for the surface gradients: One may ask, if after applying the Piola transformation also on the test function, ℎ solves the discrete Vector Laplace problem on Γ for the DG methods.One can easily check that this is true for the DG methods discussed here as there holds for instance Hence, ℎ solves the discrete Vector Laplace problem on Γ for the DG discretizations if ̂ ℎ solves the discrete problem on Γ.For the 1 -conforming methods this does not apply as -for ̂ ℎ being the discrete 1 -conforming solution to the Vector Laplacian  on Γ -the mapped function ℎ =  ̂ ℎ will not be in [ ℎ ] 3 for > 0. By construction ̂ ℎ is a continuous (2D) vector field s.t. after the transformation with Φ , which has a discontinuous Jacobian, ℎ is discontinuous and hence not included in [ ℎ ] 3 .
To illustrate this effect we consider the aforementioned methods with = 3 and ∈ {0, √ 3∕4}.In Figure 5 we display the error behavior for the two situations on 5 consecutive meshes, starting from the mesh displayed in Figure 4.Note that the geometry is piecewise planar, so that Γ ℎ = Γ.For Γ we observe that all considered methods behave essentially the same.The deviation between the errors is only marginal.All methods yield an exactly tangential field ℎ ⋅ ℎ = 0 as the equations for the -component of the 1 -conforming method completely seperates from those in and -direction.When going to Γ we observe that the 1 -conforming methods fail to converge.This is easily explained by the fact that the solution is discontinuous when represented in the embedding space ℝ 3 .Considering the representation of the solution (18), we notice that the Piola mapping incorporates the tangential field automatically and we observe that the convergence behaviour of the DG methods stay the same when going from Γ to Γ. Actually, the numbers used in the plots are identical (up to round-off errors).

Stokes
Whereas the last section demonstrated the accuracy and advantages of the (H)DG methods introduced in this work for second order elliptic problems on surfaces, we now aim to solve the stationary incompressible Stokes equations, see problem 4. Similarly as in the last section we define a Stokes problem on the flat reference domain Γ = (0, 1) × (0, 1) and map it isometrically onto a smooth manifold given by half of the surface of an open cylinder with radius 1∕ Again, Φ 1∕ preserves the length, i.e.
, and 1 , 2 as before.We set the corresponding right hand side to ̂ = −2 ( ̂ ) + ∇ ̂ .Defining the tangential vectors on Γ by = 1∕ , = 1, 2, the exact solutions are given by (see Figure 6) with ( ̂ , ̂ ) = Φ −1 1∕ ( , , ) = ( , arcsin( − 1)∕ + 1∕2), and the right hand side In the following we compare two different discretizations.The first one is the (div Γ )-conforming HDG discretization, see Section 3.6.1 and Section 3.3.The second one is based on the HDG formulation for the Vector Laplacian of section 3.6.2.For the divergence constraint we now further introduce the bilinear form For a fixed velocity approximation order the HDG Stokes discretization then read as: Find ℎ , In Figure 7 the error behavior for the above Stokes problem (with right-hand side (20)) is given for both discussed discretizations with a fixed viscosity = 1 and different polynomial orders = 2, 3.As before we use the labels HDG and Hdiv-HDG.For in the SIP stabilization we take 10.Note, that in contrast to the previous section we used the geometry approximation order = + 2. A numerical investigation showed that we have to apply this enhanced geometry approximation (only needed for this particular example) due to the ill-conditioned inverse of the sine function at 1 and −1 (needed for the calculation of the right-hand side).The first, third and fourth plot show the 1 -semi-norms error and the 2 -norm error of the velocity and the 2 -norm error of the pressure for both methods where we used an initial triangulation with 100 elements and 3 refinement levels.As we can see, all errors convergence with optimal order and the accuracy of both methods is approximately the same.The second plot shows the 2 -norm of the surface divergence.As proven by Lemma 1 the solution of the (div Γ )-conforming method is exactly divergence-free, whereas the standard HDG method is only weakly divergence-free.However the divergence error of the HDG method stills shows the optimal convergence order.
To emphasize the importance of exactly divergence-free velocity solutions we now focus on pressure robustness.Beside the observations discussed in Section 3.5 we want to discuss pressure robustness with respect to the arising error estimates.A standard a priori error estimate of inf-sup stable Stokes discretizations usually reads as (ℎ 3 )

FIGURE 7
Error behavior for the Stokes problem of section 4.2 for the (div Γ )-conforming HDG method and the HDG method given by ( 21) on 4 successively refined meshes (uniform refinements) for a fixed viscosity = 1, two polynomial orders = 2, 3 and an geometry approximation order = + 1.
where ‖ ⋅ ‖ 1 ,ℎ is an appropriate (H)DG-version of an 1 -norm, and is a constant independent on the mesh size ℎ and the viscosity .Above estimate shows that the velocity error may depend on the best approximation of the pressure including the factor 1∕ , hence the velocity error can blow up in the case of vanishing viscosity ≪ 1.As discussed in the literature, see for example 45,75 , methods that provide exactly divergence-free velocities allow to derive an error estimate that reads as where ( ) is a function that only depends on the exact solution (and not ) and shows optimal convergence properties (for methods yielding exactly divergence-free solutions we have ( ) = 0).Hence, these methods show no bad behavior for small values of the viscosity .Note that above estimate assumes an exact geometry representation Γ ℎ = Γ, compare the discussion below.
In Figure 8 the 1 -seminorm error of the solution of the same problem as above for varying viscosities = 10 −6 , … , 1, two different polynomial orders = 2, 3 and different geometry approximation orders = + 2 for = 1, 2, 3 is given.We can make several observations.First, the errors of the weakly divergence-free HDG method, see equation (21), show the expected behavior: Reducing the viscosity leads to a blow up of the 1 -semi norm error of the velocity independently of the geometry approximation.We can clearly see the scaling 1∕ for all combination of approximation orders and geometry approximations .Next, note that the (div Γ )-conforming method is expected to be pressure robust as the discrete velocity is exactly divergence-free.As we can see, this is indeed true if the geometry approximation is accurate enough.In the case = 3 + 2 the error is constant and independent of the choice of .Reducing the approximation order shows that the error is again affected by a change of the viscosity, however the error is still much better compared to the standard HDG method.This behavior is also known from the flat case, and the problem comes from an inexact evaluation of the right-hand side integral As all integrals are computed using a quadrature rule the right-hand side integral might not vanish in the case where = ∇ Γ and ℎ is divergence-free, i.e. div Γ ( ℎ ) = 0, compare Section 3.5.Beside choosing a high order quadrature rule, we now also have to increase the geometry approximation to ensure that a gradient field = ∇ Γ is also a gradient field on Γ ℎ , otherwise again the integral might not vanish even if a high order quadrature rule is used.

Generalization of the Schäfer Turek benchmark for surfaces
We consider the setup of a standard benchmark test case, cf.Ref. 76 Case 2D-2.First we recall the setup in the plane and afterwards apply different mappings to obtain geometries in 3D, some of them as in Ref. 38 .
Here, ̄ = 1 is the average inflow velocity and the viscosity is fixed to = 10 −3 which results in a Reynolds number = 100.This setup results in a time-periodic flow with vortex shedding behind the obstacle, cf. Figure 9.
The (time dependent) quantities of interest in this example are the forces that act on the disc ̂ • = 0.05 ((0.2, 0.2)) and the pressure difference before and behind the obstacle: This benchmark problem is well studied in the literature and reference values are available in Ref. 76 .Below, we consider similar setups on geometries Γ = Φ ( Γ) for mappings Φ , = 1, .., 4 specified below.For ∈ {in, out, , •} we define the boundary segments accordingly, i.e. ∶= Φ (̂ ) and prescribe natural outflow conditions on out , homogeneuous Dirichlet conditions on and prescribe inflow velocities on in .For the inflow velocities we take (for sake of comparibility) the choice of Ref. 38 : Let us stress that this results in a tangential velocity, but may result in an average inflow velocity that may deviate from ̄ = 1.We note that the degrees of freedoms of our (div Γ )-conforming HDG formulation naturally fit the tangential boundary conditions.The inflow which is in co-normal direction corresponds to the degrees of freedom of the ℎ whereas tangential flow conditions (here: zero) correspond to the boundary degrees of freedom of Λ ℎ .

Computational setup
We solve the surface Navier-Stokes problem for four different mappings, but use -for the most part -the same computational setup that we want to describe here first.As initial conditions for the Navier-Stokes problem, we take the solution of a Stokes problem with the same boundary data.To make sure that the simulation reached the time where the periodicity is established we simulate the problem for 30 time units.For the time stepping we use a second order implicit-explicit (IMEX) time stepping method and consider up to three different time step sizes: Δ = 2 1− ⋅ 10 −3 , ∈ {1, 2, 3}.We consider up to three mesh levels with a characteristic mesh size (w.r.t. the reference configuration) ĥ = 2 1− ⋅ 0.05, ∈ {1, 2, 3} corresponding to mesh levels = 1, 2, 3.For the non-isometric mapping in Subsection 4.3.2we use the same unstructured meshes as in the flat case, resulting in 857 ( = 1), 3179 ( = 2) and 13081 ( = 3) triangular elements, cf.Fig. 9 for the first two levels.The meshes are not specifically adapted to improve the approximation of potential boundary layers.For the isometric mappings in Subsection 4.3.1 we consider only one mesh with 865 elements.In all examples in this section we fix = 4 and = 5.

Isometric mappings
First, we consider two mappings Φ , = 1, 2 that are isometric, i.e. there holds ̂ ( ̂ , ̂ ) = Γ (Φ ( ̂ ), Φ ( ̂ )), ̂ , ̂ ∈ Γ where ̂ is the two-dimensional Euclidean distance and Γ is the geodesic distance on the surface Γ = Φ ( Γ) or equivalently det( ⋅ ) = 1 with = Φ ′ .Similar to the setup in Section 4.1.3we consider as a first case To consider the kink in the geometry properly we use a slightly different mesh in this setup than in all others of this subsection by making sure that the line ̂ = 1.1 corresponds to a mesh line, c.f. Fig. 10.Thereby the mapping Φ can be represented exactly in a finite element space so that the mapped mesh introduces no additional geometrical error.We note that already the mesh for Γ includes (small) approximation errors due to the approximation of the circle.
As a second example we consider Φ 2 ( ̂ , ̂ ) = (− 0.41 cos( ̂ 0.41 ), 0.41 (1 − sin( ̂ 0.41 ))), which corresponds to a bending of the flat geometry around the -axis.This time, when projecting into a finite element space, we can not represent Φ 2 exactly, hence the isometry property will only be fulfilled approximately.In Fig. 10 we show the geometries, the vorticity at a fixes time ( = 26) and plots that compare Δ ( ) in a small range of time ( ∈ [26, 26.2] on the coarsest mesh level = 1 with time level = 2.We observe that there is hardly any difference between the results.The differences between the results for Γ and Γ 1 are only due to round-off errors and not visible even after heavily zooming in.For Γ 2 we also only observe differences in the 6th digit.

Discussion of results
We observe that the energy dissipation is similar for all cases (up to a constant shift due to different initial kinetic energies).Especially, the kinetic energy is monotonely decreasing.The same holds for the enstrophy for ≤ 100.For later times the results deviate significantly.The reason for those deviations can be best explained by the palinstrophy evolution in Fig. 14.Until ̄ = 100 the simulations of the four cases agree very much, at least qualitatively.Initially four vortices form which eventually merge to two vortices around ̄ ∈ (30, 60).Until that point of the evolution the rigid body rotations in the most upper and most lower part are essentially not influenced by the interactions in the center.This changes at around ̄ = 90 where the rotations of the top and bottom are perturbed for Γ 2 and Γ 3 and somewhat later (around ̄ = 120) also for Γ 1 .For Γ 1 and Γ 2 this results in the merging of the latest two vortices to one vortex which is also reflected by an increase in the palinstrophy.For Γ 3 the perturbation of the rotations seems to stay confined and the decreased height seems to surpress the interaction of the vortices so that even at ̄ = 200 the latest two vortices did not merge yet.
Let us note that we also observe a deviation of the evolution computed for Γ 0 and the reference solution, especially after ̄ = 130 which is in agreement with the extreme sensitivity of the problem to (numerical) perturbations observed in Ref. 48.The final merge is typically observed sooner the higher the perturbations in the simulation are.In contrast to the 2D reference solution we consider a much coarser, not structured and not symmetric mesh with an additional geometry error due to the curved representation.

Smooth manifold: the sphere
Now, we consider the Kelvin-Helmholtz problem on a smooth manifold, the unit sphere, Γ 4 = 1 (0).The upper and lower halfs rotate in different directions with a perturbation similar as before in (22).However, we change the perturbation magnitude to = 2 ⋅ 10 −3 (due to the increase of the equador length from 1 to 2 ) and choose a different initial perturbation mode with = 16, = 1, = 20 and = 0.1.We use = 5, = 6 on a unstructured mesh consisting of 5442 triangles.A few snapshots of the flow and the evolution of the palinstrophy are depicted in Fig. 15.We observe that initially eight vortices form from the initial perturbation at around ̄ = 50.Diffusion takes its time until it drives the interaction of two neighboring vortices which pair up to 4 vortices at around ̄ = 200.These vortices take even longer to eventually pair up to two vortices at around ̄ = 500.The two vortices are positioned on opposite sides rotating in opposite directions.The overall evolution of the palinstrophy and the times where vortices pair up is very similar to those from the cylindrical setup in the previous section.However, due to the larger length scales the interaction between the vortices takes more time.

Stanford bunny
As a final example we want to demonstrate that the aforementioned methods can be applied on essentially arbitrary geometries.As a prototype of a complex geometry we take the famous Stanford bunny which fits approximately in a bounding box of size 80 × 80 × 60 and initially put 60 vortices on the surface.Let us note that in Reuther and Voigt 32 this geometry has also been considered (with different initial data) with the aim to investigate defect configurations which is not the aim here.The initial velocity is taken as the surface curl of

CONCLUSION AND OUTLOOK
In this work we introduced new numerical methods for the discretization of incompressible flows and vector valued elliptic problems on two-dimensional manifolds.Abandoning the 1 -conformity (originally demanded by the considered problems) we applied the Piola transformation to construct finite elements which are exactly tangential.Based on these findings we presented non-( 1 (Γ)-)conforming (hybrid) discontinuous Galerkin discretizations which performed extremely well in several considered numerical examples.Among other benefits, it was shown that the resulting methods can outperform 1 -conforming discretizations in the aspect of computational costs, and that they can deal with piecewise smooth manifolds which is a unique property so far.
For incompressible flow problems the (hybrid) discontinuous Galerkin methods were modified to be (div Γ )-conforming.The resulting methods provide important properties such as exactly (surface) divergence-free velocity fields, energy stability and pressure robustness.We studied certain numerical examples showing that the methods are highly accurate and applicable to deal with complex geometries.In particular we also showed a robustness of the new methods with respect to isometric mappings.
The extension of the discretizations to non stationary (evolving) manifolds is left for future research.Further, due to the increasing interest on foldable LC displays and flows on cell membranes, the coupling of the presented discretizations with other PDEs on surfaces is an interesting task for the future.

FIGURE 5
FIGURE 5 Error behavior for Vector Laplace problem on the flat (top row) and the house of cards geometry (bottom row) from Section 4.1.3for four different discretization methods on 5 successively refined meshes (uniform refinements) for fixed discretization order = 3 and an exact geometry approximation.

FIGURE 9
FIGURE 9 Flat Schäfer-Turek benchmark problem 2D-2.The first two meshes used in the simulations (left), discrete solution of the Stokes problem (center) and discrete solution of the unsteady Navier-Stokes problem at a fixed time (right).

20 ‖
− ‖ 2 )where the , = 1, .., 60 are some randomly located but sufficiently separated (‖ − ‖ 2 ≥ 3.5, ≠ ) vertices of the mesh.For the viscosity we choose =1 50 .The flow is again only driven by its initial condition.Pictures of the numerical solution on a mesh with 6054 triangles, = = 4, time step size Δ = 1 200 and = 80 are shown in Fig.16.We observe that the flow undergoes a process of self-organization as it is well-known from 2D flows.At = = 500 the vortices merged successively into two remaining vortices: One around the ears of the Standford bunny, one on the breast.In contrast to the studies in Reuther and Voigt 32 , the flow is not yet in a quasi stationary state at = 500, i.e. these two vortices will continue moving and possibly merge while the magnitude of the velocity will decay over time.

TABLE 1
Comparison of different computational quantities (dof: degrees of freedom, gdof: global degrees of freedom that remain after static condensation (if applicable), nze: non-zero entries in system (Schur complement) matrix) for the different schemes on finest level.Numbers in green indicate that the corresponding method yields the best values in the current column.