{'fraction_non_alphanumeric': 0.05044798864920255, 'fraction_numerical': 0.011114322499277438, 'mean_word_length': 4.383309759547383, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Meshfree finite differences for vector Poisson and pressure Poisson equations with electric boundary conditions

We demonstrate how meshfree finite difference methods can be applied to solve vector Poisson problems with electric boundary conditions. In these, the tangential velocity and the incompressibility of the vector field are prescribed at the boundary. Even on irregular domains with only convex corners, canonical nodal-based finite elements may converge to the wrong solution due to a version of the Babuska paradox. In turn, straightforward meshfree finite differences converge to the true solution, and even high-order accuracy can be achieved in a simple fashion. The methodology is then extended to a specific pressure Poisson equation reformulation of the Navier-Stokes equations that possesses the same type of boundary conditions. The resulting numerical approach is second order accurate and allows for a simple switching between an explicit and implicit treatment of the viscosity terms.


Introduction
The numerical approximation of vector fields that are incompressible is often times challenging because incompressibility, ∇ · u = 0, is a global constraint that may not fit within the framework of simple discretization approaches of the complete problem. The instationary incompressible Navier-Stokes equations (NSE) represent a prime example, in which the timeevolution of the velocity field is given, but not of the pressure (which is a Lagrange multiplier associated with ∇ · u = 0). As a consequence, there is no single canonical way to advance the NSE forward in time. Similarly, in electrostatics the electric field and the magnetic potential are solutions to vector Poisson equations with divergence constraints.
One methodology to circumvent the incompressibility constraint inside the computational domain is to formulate a different problem that imposes ∇ · u = 0 as a boundary condition instead. Under certain circumstances, this new problem has the same solution as the original problem, while at the same time giving rise to new numerical approximation methods. For electrostatic problems (see §2) this approach is employed and analyzed in [23,16], and in the context of incompressible fluid flows (see §3) it has been proposed in [30]. The specific boundary conditions for these problems consist of enforcing the tangential component(s) of the solution, together with the condition ∇ · u = 0. Due to their occurrence in electrostatics, they are often called (perfect) electric boundary conditions (EBC), a terminology that we follow in this paper.
A fundamental question is what are simple numerical approaches to approximate the solutions to vector-valued problems with EBC. In [30] an immersed boundary staggered grid approach has been proposed. While convergent, this approach is definitely not a simple or canonical method. A seemingly more natural approach for problems on irregular domains is the standard nodal finite element method (FEM). Interestingly, for the class of problems at hand, nodal FEM can exhibit a Babuska paradox, i.e., for domains with curved boundaries the sequence of FEM approximations can converge to a wrong solution (see §2.2). Domains with re-entrant corners (see §2.1) pose additional challenges; however, this last aspect is not the focus of this paper. The convergence of the FEM can be recovered by converting to a mixed FEM formulation, as conducted in [5]. However, the simplicity of nodal FEM does not carry over to mixed FEM.
In this paper, we present yet another alternative for problems with EBC, namely meshfree finite differences (FD). These are generalizations of traditional grid-based FD that apply to clouds of points without any connectivity between them. The meshfree FD methodology is presented in §2.3 and §2.4, and its application to the vector Poisson equation with EBC is shown in §2.5. Structurally, the approximation of a general PDE boundary value problem via meshfree FD is very straightforward: any differential operator, whether in the domain's interior or on the boundary, is approximated via a meshfree FD stencil. Thus, the problem is directly transformed into a finite dimensional system, in which each individual equation corresponds to the governing condition that holds at a particular point of the cloud.
We demonstrate (in §2.5) that for the vector Poisson equation with EBC, meshfree finite differences do not exhibit the Babuska paradox, and furthermore that there is no conceptual problem to obtain higher-order accuracy (we test the method up to third order convergence).
We then move on to time-dependent problems. First, the meshfree FD method is extended to the vector heat equation with EBC (see §3.2), where an explicit or an implicit time-stepping can be conducted. Then, by adding nonlinear convective terms and a pressure, the approach is further extended to the incompressible Navier-Stokes equations (NSE). Specifically, we consider a pressure Poisson equation (PPE) reformulation of the Navier-Stokes equations. The idea of PPE reformulations (see §3) is that an operator function p = P (u) is formulated that yields (via the solution of a Poisson equation) the pressure p to any given velocity field u that solves the NSE. Here, we focus on a specific PPE reformulation, proposed in [30], which prescribes EBC for the fluid velocity, the motivation for which is outlined in §3.1. We demonstrate how a meshfree FD approximation for the full PPE reformulation can be constructed (see §3.3), and show computational results for a resulting numerical scheme that is second order accurate in space and time, and that allows for a choice of an explicit or an implicit treatment of the viscosity (see §3.4).

Vector Poisson Equation
The vector Poisson equation (VPE) arises, for instance, in problems in electrostatics. The electric field satisfies ∇ · E = ρ, where ρ = ρ(x) is the (normalized) charge density. Using the fact that ∇ × E = 0, this implies the VPE ∆E = ∇ρ. Moreover, if the boundaries of the domain are perfect conductors, then the vector field is perpendicular to the boundary, i.e., n × E = 0, where n is the outer surface normal vector. Another example is the magnetic potential, which satisfies the VPE ∆A = −J, where J = J(x) is the (normalized) electric current density. The Coulomb gauge yields ∇ · A = 0, and the boundary condition n × A = 0 represents a zero magnetic magnetic field (see [9] for more details). Motivated by the structure of these examples, we here consider the VPE on Ω n × u = n × g on ∂Ω . (1) As motivated in §1, it can be desirable to remove the divergence condition that holds in the whole domain. In the following, we outline how this can be achieved.

Electric Boundary Conditions.
Let Ω be a bounded, simply connected domain with Lipschitz boundary ∂Ω. Moreover, in this paper we restrict to domains with boundaries ∂Ω which are piecewise C 2 and convex (see Remark 1). We denote by n the outward unit normal vector along the boundary (that is defined almost everywhere). The vector Poisson equation (VPE) with electric boundary conditions (EBC) takes the form where the source is incompressible, i.e., ∇ · f = 0. Note that in contrast to problem (1), problem (2) possesses no source-free condition in the domain's interior. Instead, ∇ · u = 0 is specified as an additional boundary condition. Clearly, any solution of (1) is also a solution of (2). Moreover. . . (2) is in H 2 (Ω), then it is also a solution to (1).

Lemma 1. If the solution to
Proof. Define φ = ∇ · u. Then φ is a (weak) solution of the problem which has the unique solution φ ≡ 0. Hence ∇ · u = 0 in Ω.
Remark 1. As shown in [18], the assumption of Lemma 1 is satisfied for the domains considered in this paper. However, it is not satisfied if the domain Ω has re-entrant (i.e., non-convex) corners. In such a case, the physically relevant (i.e., source-free) solution to (1) is not in H 1 , while problem (2) possesses a solution in H 1 , however, one that does not satisfy ∇ · u = 0 inside Ω. In this paper we exclude this possibility, and for the domains considered here (see above) the problems (1) and (2) are in fact equivalent (see [18] for a proof).  (2), Moreover, let H 1 0t (Ω) N denote the associated homogeneous (i.e., g = 0) Hilbert space. There are then two equivalent weak formulations of (2). To obtain the first formulation, we use the identity −∆u = ∇ × (∇ × u) − ∇(∇ · u) and follow the standard procedure of multiplying the first equation in (2) by a test function v ∈ H 0t (Ω) N , integrating by parts, and applying the boundary conditions to the boundary integral to obtain where the bilinear form is Based on this, the first variational formulation of (2) reads as: Note that due to (3), the condition (∇ · u) ∂Ω = 0 arises as a natural boundary condition. It is this formulation (VP1) that we implement in the numerical test in §2.2.3.
In obtaining the second variational formulation, we restrict the derivation to the case g = 0, because this case is enough to show that the Babuska paradox can arise. Assume for a moment that u ∈ H 2 (Ω) N . Then, using the fact that ∆u = ∇ · (∇u), we multiply the left hand side of (2) by v ∈ H 0t (Ω) N and integrate by parts to obtain Combining (3) and (4), and using that n × v = 0 on ∂Ω, we can rewrite a(u, v) as a new bilinear form Moreover, since n × u = 0 on ∂Ω, we can expand the divergence on the boundary as where κ is the local curvature which is defined almost everywhere. Using (5), we can write thus giving rise to a different variational formulation: Clearly, by construction the bilinear forms are equal, a(u, v) = b(u, v), for functions in H 2 (Ω) N . In fact, as shown in [18], the equality also holds if the functions are in H 1 (Ω) N .
Moreover, the standard theory shows that the bilinear forms are coercive and continuous on H 0t (Ω) N (since Ω is simply connected) so that by the Lax-Milgram theorem there is a unique solution to (VP1) and (VP2). Moreover the variational problems (VP1) and (VP2) have the same solution.  2.2.2. Babuska paradox. Using the just derived weak formulations, we prove the possible occurrence of the Babuska paradox. Note that other proofs have been provided before, such as in [32].
Theorem 2. When solving the vector Poisson equation (2) using a nodal FEM implementation of (VP1), one may encounter the Babuska paradox.
Proof. Suppose that u h solves (VP1) using nodal elements, a triangular mesh, and a regular polygonal domain Ω h . Here h denotes the diameter of the largest mesh element, so that Ω h → Ω (in the appropriate sup-norm sense) as h → 0. Now, under the current assumptions on the domain Ω, the problems (VP1) and (VP2) have the same weak solution. Then by the equivalence of the two problems, u h also solves (VP2). For any given mesh, however, the boundary of Ω h has flat sides with κ = 0. Consequently, the weak solution u h solves (VP2) with κ = 0, i.e., for each test function v h . Therefore the solutions u h converge to the function u * which solves the limit problem In other words, the nodal FEM solutions u h converge to a solution where κ is artificially set to zero, or equivalently to a problem where one replaces the boundary condition ∇ · u = 0 with du dn = 0. Hence, for an arbitrary domain (with boundaries that are at least partially curved), generally u * does not equal the true solution of the problem (VP1).
We consider two triangulations (of different resolutions) that approximate the domain Ω via polygons with straight edges. To define normal vectors at the boundary vertices, we adopt the method introduced in [10], which obtains normal vectors as suitable averages of normal vectors at the edges connecting to the boundary vertex. Using standard nodal-based finite elements, or Lagrange FE [7], we implement the tangential boundary condition in an essential fashion (by choosing the solution space H 1 gt (Ω) N ), and leave the normal component unprescribed, with the idea that ∇ · u = 0 follows naturally (see above). We use quadratic C 0 elements.
The results of the two FEM approximations are shown in Fig. 1. The true vector field (6) is shown by red arrows, and the approximate FEM solution is given by black arrows. In each case, the mesh that is used to conduct the computation is shown in the background. The right panel shows a computation with a mesh that is twice as fine as the one in the left panel. The two numerical solutions are almost identical, and hence they can be interpreted as converged (in the eye-norm) to the limit (h → 0) solution u * of the nodal FEM. Clearly the FEM solution u * is different from the true solution u, thus confirming the occurrence of the Babuska paradox.
Remark 2. Within the framework of FEM, the Babuska paradox can be overcome by moving to a mixed FEM formulation (cf. [27,32,2,1]). The idea of mixed FEM for the VPE (2) is introduce an additional variable, σ = ∇ × u, and transform −∆u = f into two equations: σ = ∇ × u and ∇ × σ − ∇(∇ · u) = f . This framework allows one to use Raviart-Thomas elements [27] for the approximate vector field u h (and standard nodal elements (2D) or Nédélec elements [24] (3D) for σ h ), and to incorporate the tangential velocity boundary condition n × (u − g) = 0 as boundary integrals into the weak formulation, rather than into the solution space. In a companion paper [5], we apply high-order mixed FEM to the pressure Poisson equation reformulations of the Navier-Stokes equations devised in [30]. While mixed FEM overcome the Babuska paradox, this framework is clearly not as simple as nodal FEM, or as meshfree FD, described below.
2.3. Meshfree Finite Difference Method. Meshfree finite differences (FD) generalize classical FD that are defined on regular grids: at a given point, a differential operator of a smooth function is approximated via a combination of function values at nearby points. The selection of points and the corresponding weights are called the stencil. In the same way as grid-based FD, meshfree FD can be derived in two ways: as derivatives of suitable local interpolants of the data (cf. [22,8,21]), or via Taylor expansion of the solution (cf. [29]). Here we outline the second methodology.
Consider a point cloud that consists of interior points (inside Ω) and boundary points (on ∂Ω); see Fig. 2 for an example. For a point x i , let a neighborhood B i be defined. Here we employ circular neighborhoods, i.e., B i = {j : x j − x i ≤ r}, where r is an appropriately chosen radius (see below). However, many other types of neighborhoods are possible [28]. Now define the relative coordinatesx ij = x j − x i and Taylor-expand the solution u(x) around Note that in the quadratic term, the matrix scalar product A : C = i,j A ij C ij and the outer product ofx ij with itself are used. While here we stop at the quadratic term, the expansion can of course be carried out further (or less far). A linear combination (with weights a ij ) of nearby solution values yields If (7) is supposed to approximate a given differential operator applied to the solution, then the stencil weights a ij must satisfy certain constraints. For instance, for (7) to approximate ∆u(x i ), it is required that which in 2D gives rise to the linear system of constraints Here, the stencil vector a i does not contain the diagonal entry a ii . Its value is obtained as a ii = − j∈B i \{i} a ij . Moreover, m i = |B i | − 1 is the number of neighbors of x i . If the radius r is chosen large enough that m i ≥ 5 ∀ i, and if the point cloud generation (see §2.4) ensures that no pathological point configurations arise (see [28] for examples), then system (8) always has a solution, and the resulting approximation is (at least) first order accurate.
If m i > 5, system (8) in general has infinitely many solutions. One way (employed here) to single out a unique solution is via a weighted least-squares (WLSQ) minimization problem where the weights are decreasing with the distance, . The solution of (9) is Note that an alternative approach (not employed here) would be to solve system (8) in an 1 sense, i.e., which would generate optimally sparse stencils [29].
Other differential operators are approximated in an analogous fashion. For instance, a first order approximation to ∂ x u is obtained by setting and a second order approximation to ∂ x u is obtained by setting Applying this procedure to each equation (at all interior points) and boundary condition (at all boundary points) of the vector Poisson equation (2) leads (here is 2D) to the linear system * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * In this system, u = (u x , u y ), f = (f x , f y ), and the function g = n × g at the boundary, and g i = n i × g i . Moreover, the total number of points is N , and the number of interior points is N i . In the sparse block matrix, the first two rows of blocks correspond to the two components of the Poisson equation at the interior points; the third block row encodes the ∇ · u = 0 boundary condition; and the fourth block row represents the n × (u − g) = 0 condition. The first two block columns corresponds to the two vector field components at the interior points x 1 , . . . , x N i ; and likewise the last two block columns correspond to u = (u x , u y ) at the boundary points x N i , . . . , x N . Each empty block is devoid of nonzero entries, and each block with stars possesses some nonzero entries, corresponding to the neighboring points of each central point.  The generation of a cloud of N points is initiated by placing N points randomly inside Ω. After that, the points are moved according to repulsive forces experienced from nearby points, i.e., the point x i moves according tȯ where B i is a circular neighborhood around x i , and v max is some upper bound on the repulsion. The law of motion (11) is further constrained by x i ∈Ω, i.e., no point can ever move beyond the domain boundary to leave the domain. A simple way to implement this constraint is to allow points to slightly leaveΩ, but then to immediately project them back onto ∂Ω, using the normal ∇φ/|∇φ| that is defined also outside ofΩ. Finally, points that are inside Ω but too close to the boundary are also projected onto ∂Ω, thus preventing interior points from being too close to the boundary (see [29] for why this would be undesirable). The law of motion (11) is then applied to all points until the amount of motion has fallen below a given threshold. The resulting point clouds are unstructured, and tend to be quite uniform (i.e., the ratio between the minimum distance between points and the radius of the largest ball that contains no points (cf. [21]) is quite large). An example of a point cloud associated with the domain defined in §2.2.3 is shown in Fig. 2.
For the definition of a "mesh" resolution h of a given point cloud, there is a variety of possible choices [21]. Here we use a simple averaged concept of resolution, defined as follows. One type of configuration in which nearby points are equidistant is an optimal sphere packing, which in 2D is a hexagonal lattice, i.e., adjacent points form equilateral triangles. Each point is a corner of six triangles, and each triangle is shared by three points. One can therefore associate to each point 1/3 of each of the six triangles, resulting in an area per point of A =     (2) is discretized via the procedure described in §2.3. We conduct the meshfree FD approximation for three different orders: first order accuracy, i.e., the Taylor expansion in (7) is carried out up to the quadratic term for ∆u and up to the linear term for ∇ · u; and second and third order accuracies, for which the Taylor expansions in (7) are carried out further accordingly.
The numerical approximations obtained for the different orders and mesh resolutions are then compared to the true solution in the maximum norm, taken over all points. We consider the errors in the vector field u itself, as well as its Jacobian ∇u (which is important for calculating stresses at the boundary when u represents a velocity field) and its divergence ∇ · u (which by the equivalence of problems (1) and (2) should be close to zero). All derivative quantities are obtained from the vector field u via meshfree FD stencils of fourth order. Hence, if an order of less than 4 is observed, we know that this is the true accuracy of the numerical result.
The error convergence of these quantities is shown in Fig. 3, for the approximation orders 1 (left panel), 2 (middle panel), and 3 (right panel). The results show that all approaches converge as h → 0, and the convergence orders equal the local approximation orders. In particular, the convergence orders of the derivative quantities are the same as those of the vector field itself. This is an important property that finite difference methods commonly exhibit, and that is in contrast to finite element methods that frequently lose an order of accuracy when a derivative quantity is evaluated (in the sense of functions).

Pressure Poisson Equation Reformulation of the Navier-Stokes Equations
We consider the time-dependent incompressible Navier-Stokes equations (NSE) with compatibility conditions u = g on ∂Ω × {t = 0} (continuity between i.c. and b.c.) Due to the lack of a time evolution of the pressure, there is no single canonical way to numerically advance (12) forward in time. One class of approaches to do so is based on approximating the time derivative in the momentum equation, and solving for u and p in a fully coupled fashion. This methodology is accurate, but also costly, because a large system must be solved that possesses a saddle point structure. An alternative class of approaches decouples the pressure solve from the velocity update. This methodology was first proposed in the form of projection methods [6,31], and later employed in approaches based on pressure Poisson equation (PPE) reformulations of the NSE [13,14,17,15,30].
A fundamental difference between projection methods (see [12] for an overview) and PPE reformulations is that projection methods are based on a fractional step approach in which the time-evolution of the velocity field and its projection onto the space of incompressible fields are alternated. In contrast, PPE reformulations derive an equation for the pressure that replaces the incompressibility constraint in (12) by a global pressure function p = P (u) that is designed so that the solutions of the PPE reformulation are identical to the solutions of the original NSE. As a consequence, numerical methods based on PPE reformulations are structurally easy to extend to high order accuracy in time. In addition, they do not suffer from poor spatial accuracy near boundaries. A difficulty of PPE reformulations is that the Poisson equation for the pressure can involve complicated expressions, whose interaction with the velocity field equation is not always easy to understand and analyze. Another important property of PPE reformulations is that, unlike the original NSE, they are also defined if the initial conditions are not incompressible, see §3.1.

PPE Reformulation with Electric Boundary Conditions.
In this paper, we are concerned with the particular PPE reformulation of the NSE proposed in [30]. Its fundamental difference from previously proposed PPE reformulations [13,14,17] is that the velocity field satisfies electric boundary conditions, i.e., incompressibility and the tangential flow are prescribed at the boundary. In turn, the normal velocity is enforced via a relaxation term in the pressure equation (see [30] for a discussion on the choice of λ). The PPE reformulation consists of the momentum equation where P (u) is a solution of the associated pressure Poisson equation While a variety of modifications and additions can be applied to these equations (cf. [30,5]), here we study the equations exactly in the given form (with one small caveat regarding the compatibility of the pressure boundary conditions, see below).
If the initial conditions are not incompressible, i.e., ∇ ·ů = 0, then the solution of (13) relaxes towards a solutions of (12), for the following reason. Let φ = ∇ · u. Then the application of ∇· to the momentum equation in (12), and the use of the first equation in (14) yields that φ satisfies the heat equation with homogeneous Dirichlet boundary conditions This property is of great relevance. It means that in PPE reformulations, there is no need to impose a discrete incompressibility principle. If, due to numerical approximation errors, the numerical solution starts to drift away from the ∇ · u = 0 manifold, equation (15) ensures that it is pulled back towards incompressibility.
At the same time, the fact that the numerical solution may not be exactly incompressible, implies that the compatibility condition in the pressure Poisson equation (14) may be violated. Specifically, (14) has a solution if and due to numerical approximation errors (or because a problem with ∇·ů = 0 is considered) this condition may be violated. However, whenever this occurs, the solution of the augmented system (22), described in §3.3 projects the right hand side of (14) in a way that the solvability condition is satisfied.
Below, we first generalize the meshfree finite difference methods developed in §2 for the vector Poisson equation to the vector heat equation ( §3.2). Then, we extend the methodology to the PPE reformulation ( §3.3).
The schemes that treat ∆u implicitly lead to relatively simple modifications of the linear system (10) of the VPE (17). For instance, the first-order scheme, amounts to the following modifications of system (10): (i) The vector of unknowns becomes the new vector field at time t + ∆t. (ii) In the right hand side vector, the function g is evaluated at time t + ∆t. (iii) In the top two block rows of the system matrix, multiply all entries by ν and add 1/∆t to the diagonal entries. (iv) Add 1/∆t times the solution at time t to the top two blocks of the right hand side vector. Everything else remains unchanged.
In the forward Euler case, the update at the interior points becomes explicit, while the boundary conditions are still implicit. Hence, to advance the solution from time t to t + ∆t, one first updates at each interior point explicitly according to (7) and (9), and after that solves a small linear system for the boundary points that results from the two bottom block rows of system (10), where the left two block columns are brought to the right hand side (using the just updated interior point values).
To study these numerical schemes, we generalize the manufactured solution from §2.2.3 to the time-dependent case. We set ν = 1, prescribe the incompressible solution u(x, y, t) = π cos(t) sin(2πy) sin 2 (πx) −π cos(t) sin(2πx) sin 2 (πy) , and calculate the forcing f = ∂ t u − ν∆u, the boundary velocity g = u, and the initial conditionsů = u(t = 0) accordingly. Using this test case, we first determine, via numerical experiments, the maximum time step that the forward Euler scheme admits to be stable. We find ∆t ≤ Ch 2 ν , where for the given point clouds, C always lies between 0.2 and 0.3. Then, we study the convergence orders of the numerical schemes. The results are shown in Fig. 4. We use a second order meshfree FD approximation in space, and conduct five kinds of time stepping: forward Euler and backward Euler with ∆t = 0.2h 2 (left panel; the plots of forward vs. backward Euler are indistinguishable); backward Euler with ∆t = h (not shown; due to very small temporal errors, the convergence looks like second order); backward Euler with ∆t = 100h (middle panel; the temporal error is visible and yields the expected drop to first order); and the second-order ImEx scheme (19) (right panel). From these results,  we can see that the solution, and its derivatives, are in fact second-order accurate in space. Moreover, explicit and implicit time stepping (of first and second order) can be conducted without problems; and the temporal errors are relatively small. (13) is the same as the vector heat equation (17), "just" with the nonlinear term N(u) = (u · ∇)u and the pressure term ∇P (u) added to the time evolution. We treat both of these terms, as well as the forcing, explicitly. The first-order ImEx time-stepping (18) gives rise to the update rule 1 ∆t (u(t + ∆t) − u(t)) = −N(u(t)) − ∇P (u(t)) + ν∆u(t + θ∆t) + f (t) ,

Meshfree Finite Differences for the PPE Reformulation. Structurally the PPE reformulation
where θ ∈ {0, 1} allows to switch between an explicit/implicit treatment of viscosity. In the forward Euler case (θ = 0), stability requires ∆t = O(h 2 ), and thus the scheme's accuracy is O(h 2 ) overall. In the semi-implicit case (θ = 1), one can choose larger time steps, and consequently the temporal accuracy is not sufficient. We therefore use instead the secondorder ImEx time-stepping (19) with R(u) = ν∆u and Q(u) = f − N(u) − ∇P (u), which allows for time steps ∆t = O(h) and yields an O(h 2 ) accurate scheme.
The Jacobi matrix ∇u needed for the nonlinear terms is approximated very simply via point-centered meshfree finite differences, via the methodology described in §2.3. Clearly, such a centered treatment of advection is not the most effective choice for high Reynolds numbers (i.e., ν 1). And in fact, meshfree FD are quite easily amenable to an upwind treatment (e.g., by centering the approximation around a position x i − βu(x i ), where β is a suitably chosen parameter). However, for the purpose of demonstrating the convergence of meshfree FD methods for the PPE reformulation (13), the central treatment of N(u) is sufficient.
The pressure P (u) results from the solution of the pressure Poisson equation (14). We discretize this problem via the same meshfree FD method described in §2.3, with one important deviation from the standard procedure. The right hand side of the boundary conditions in (14) requires the evaluation of ∆u at the boundary ∂Ω. While straightforward meshfree FD  for ∆u yield accurate approximations inside the domain, it turns out that low accuracy (in the form of bounded but noticeable spatial oscillations along ∂Ω) is achieved when using the same procedure at a boundary point. The reason is that the Laplacian is an operator that naturally "likes" to use data around the approximation point (cf. [29]); however, at the boundary, data in such a configuration is not accessible. We therefore employ a different approach that remedies the problem: we use the meshfree approximation of w = ∆u at the interior points (as calculated for the viscosity), and extrapolate this field w to the boundary points, using moving least squares (MLS) interpolation [19]. This aspect is visualized in Fig. 5: the black dots are the approximation errors when approximating ∆u at ∂Ω via meshfree FD; the red dots are the errors obtained when using MLS interpolation.
Since (14) is a Neumann problem, its discretization leads to a linear system A · p = r, in which the Poisson matrix A has corank 1. In fact, because in the meshfree FD expansion (7) the first term must vanish for any differential operator, the kernel of A is e = (1, . . . , 1) T . In line with the approach described in [13], we solve the augmented system whose unique solution is the one satisfying A · p = r − αe, where the new right hand side is the projection of r onto the range of A. Moreover, a unique solution is singled out by the condition e T · p = 0. This approach in particular addresses the possibility that the PPE compatibility condition (16) may be not satisfied exactly. The gradient of the resulting pressure is then approximated via standard meshfree FD at all interior points.
The maximum admissible time step of the numerical scheme is determined by the viscosity term in the explicit case, i.e., ∆t = O(h 2 ), and by the nonlinear advection term in the semiimplicit case, i.e., ∆t = O(h). These stability time step restrictions are in line with those observed for the numerical method presented in [30]. However, they are different from the phenomenon observed and analyzed in [26,15] for a different PPE reformulation and a different numerical discretization. In that study, the parabolic scaling ∆t = O(h 2 ) is observed to be required for stability, even if viscosity is treated implicitly. 3.4. Numerical Results. In order to investigate the convergence of the numerical scheme developed in §3.3, we use the same manufactured solution (20) as for the VHE, set ν = 1, and calculate the pressure p(x, y, t) = − cos(t) cos(πx) sin(πy), the forcing f = ∂ t u + (u · ∇)u + ∇P (u) − ν∆u, and the boundary velocity g = u accordingly. We use a spatially second order meshfree FD scheme (with the special treatment of ∆u| ∂Ω , see §3.3), and two types of time stepping: a) forward Euler with ∆t = 0.2h 2 ; and second-order ImEx with ∆t = 0.2h. In all cases the boundary relaxation value is chosen λ = 30. The numerical results, shown in Fig. 6, demonstrate that we obtain an overall second order convergence rate for all quantities of interest: the velocity field, its gradient, the divergence, the pressure, and the pressure gradient.
Moreover, to demonstrate the applicability of the numerical methodology, we conduct the standard benchmark lid-driven cavity test [4] for Reynolds number 100, i.e., ν = 0.01. On the domain Ω = (0, 1) 2 , the velocity field is zero at the boundaries, except for the tangential velocity at y = 1, which is 1. The initial velocity field is zero everywhere except for y = 1, where it equals the boundary condition. The numerical approach used here is the same forward Euler-based scheme as in the manufactured solution test. The results of the steadystate profile (t = 20) are shown in Fig. 7. The velocity field (left panel) is depicted in the form of normalized vectors at the approximation points. The large center vortex and the two vortices in the bottom corners are captured. Moreover, a comparison (right panel) of the velocity through the centerlines of the cavity with reference data [11] shows a good agreement even on this not very highly resolved point cloud. Note that one particularity of the PPE reformulation (13) is that the flow through the boundary need not necessarily be exactly zero due to numerical approximation errors. This is why a flow through the boundaries is visible in the scaled quiver plot. However, the actual flow through the boundary is very small.

Conclusions and Outlook
The results in this paper show that meshfree finite differences (FD) provide a relatively straightforward methodology to approximate the solutions of vector-valued elliptic, parabolic, and fluid flow problems with electric boundary conditions (EBC), on domains without reentrant corners. This is in contrast to finite element methods (FEM), whose simplest version, nodal-based FEM on triangular elements, fails at generating the correct solution. Instead, a Babuška paradox arises, which is shown to arise from the fact that nodal-based FEM do not capture the domain boundary's curvature.
For the vector Poisson equation, meshfree FD lead to a linear system that discretizes the Laplace operator at interior points, and the divergence operator at boundary points, in a natural and very systematic fashion. The same methodology is shown to yield first, second, and third order convergent numerical schemes. Analogous statements hold for the vectorheat equation. Implicit time stepping is a straightforward extension of the vector Poisson case; explicit time stepping is a bit more interesting (because boundary conditions remain implicit), but poses no conceptual complication.
The extension of the methods to a PPE reformulation of the Navier-Stokes equation with EBC is, again, conceptually not complicated. There is one challenge that must be overcome, namely the approximation of the Laplacian of the velocity field at the domain boundary. Once this issue has been addressed, a second-order accurate numerical scheme is obtained in which the pressure solve and the viscosity solve are decoupled. As a consequence, one can choose between an explicit and an implicit treatment of viscosity.
Being FD approaches, the numerical schemes yield the values of the velocity field u at the approximation points only. However, meshfree stencils can be employed to also calculate spatial derivatives of the velocity field. For instance, the velocity gradient ∇u is crucial in computing forces and stresses acting on the boundary. Our investigation of the accuracy of these derivative quantities reveals that, for all studied problems, they show no degradation in order: a k th order scheme yields k th order convergence in u, and also in ∇u.
While the results demonstrate the potential of meshfree FD for these types of problems, they also give rise to further questions. One important question is whether domains with reentrant corners can also be treated. At first glance, one would think "no", because of the lack of smoothness of the solutions on such domains. However, FD methods are known to be able to yield correct answers even for certain problems that lack smoothness (such as hyperbolic conservation laws [20]). This work also gives rise to a number of questions regarding the high-order accurate meshfree FD approximations. First, extensions of the numerical schemes for the PPE reformulation to convergence orders higher than two are of interest. Second, to avoid excessively large stencils for higher approximation orders, it is of interest whether the numerical schemes would work equally well (or better) if the meshfree FD approximations were obtained in a different fashion, such as via radial basis functions, compact FD, or deferred correction.
NSF through grant DMS-1318942, and B. Seibold would like to acknowledge partial support by NSF through grants DMS-1318641 and DMS-1318709. D. Shirokoff acknowledges partial funding by NSERC.