A conservative, consistent, and scalable meshfree mimetic method

TL;DR

This paper introduces a novel meshfree mimetic divergence operator for point clouds that ensures local conservation, is scalable, and can be instantiated with or without a background mesh, enabling robust discretization of conservation laws.

Contribution

The work defines a new meshfree mimetic divergence operator that satisfies a discrete divergence theorem and can be implemented with or without a background mesh, improving scalability and robustness.

Findings

The proposed MMD operator is first-order accurate.

It satisfies a discrete divergence theorem and local conservation.

Numerical results show mimetic properties similar to finite volume methods.

Abstract

Mimetic methods discretize divergence by restricting the Gauss theorem to mesh cells. Because point clouds lack such geometric entities, construction of a compatible meshfree divergence remains a challenge. In this work, we define an abstract Meshfree Mimetic Divergence (MMD) operator on point clouds by contraction of field and virtual face moments. This MMD satisfies a discrete divergence theorem, provides a discrete local conservation principle, and is first-order accurate. We consider two MMD instantiations. The first one assumes a background mesh and uses generalized moving least squares (GMLS) to obtain the necessary field and face moments. This MMD instance is appropriate for settings where a mesh is available but its quality is insufficient for a robust and accurate mesh-based discretization. The second MMD operator retains the GMLS field moments but defines virtual face moments using computationally efficient weighted graph-Laplacian equations. This MMD instance does not require a background grid and is appropriate for applications where mesh generation creates a computational bottleneck. It allows one to trade an expensive mesh generation problem for a scalable algebraic one, without sacrificing compatibility with the divergence operator. We demonstrate the approach by using the MMD operator to obtain a virtual finite-volume discretization of conservation laws on point clouds. Numerical results in the paper confirm the mimetic properties of the method and show that it behaves similarly to standard finite volume methods.

Figures

Tables

Table 1. Table 1: Average solution time and average number of the AMG iterations for the weighted graph Laplacian equations ( 38 ) as function of the point cloud size. Solution time and AMG iterations are averaged over the set of d ​ ( d + 1 ) 𝑑 𝑑 1 d(d+1) solves. Results indicate optimal O ​ ( p Ω ) 𝑂 subscript 𝑝 Ω O(p_{\Omega}) performance of the AMG solver.

$p_{\omega }$	Solution time (sec)	AMG iterations
${16}^2$	0.0015	5.750
${32}^2$	0.0038	6.125
${64}^2$	0.0094	5.000
${128}^2$	0.0368	6.500
${256}^2$	0.1839	6.625
${512}^2$	0.5782	7.125

Table 2. Table 2: Convergence rates of the virtual finite volume scheme ( 42 ) for the manufactured solution u e ​ x ​ ( x , y ) = sin ⁡ 2 ​ π ​ x ​ sin ⁡ 2 ​ π ​ y subscript 𝑢 𝑒 𝑥 𝑥 𝑦 2 𝜋 𝑥 2 𝜋 𝑦 u_{ex}(x,y)=\sin 2\pi x\sin 2\pi y .

N	${\ell }_2$-Dirichlet	$h_1$-Dirichlet	${\ell }_2$-Neumann	$h_1$-Neumann
$16$	0.0983022	0.193076	0.494616	0.353165
$32$	0.0516438	0.0943233	0.307722	0.200986
$64$	0.0284723	0.0484218	0.165076	0.103267
$128$	0.0184388	0.0288341	0.0771193	0.0512336
$256$	0.00831139	0.0132393	0.0405609	0.0257531
Rate	1.14958	1.12295	0.927002	0.992344

Table 3. Table 3: Convergence rates in ℓ 2 subscript ℓ 2 \ell_{2} -norm of the virtual finite volume approximation to u e ​ x ​ ( x , y ) = sin ⁡ 2 ​ π ​ x ​ sin ⁡ 2 ​ π ​ y subscript 𝑢 𝑒 𝑥 𝑥 𝑦 2 𝜋 𝑥 2 𝜋 𝑦 u_{ex}(x,y)=\sin 2\pi x\sin 2\pi y using the centered advective flux reconstruction 𝝈 ~ a ​ ( u h ) subscript ~ 𝝈 𝑎 superscript 𝑢 ℎ \widetilde{\bm{\sigma}}_{a}({u}^{h}) . An entry of n . c . formulae-sequence 𝑛 𝑐 n.c. denotes that the linear solver failed to converge due to ill-conditioning.

$N$	$Pe=1$	$Pe=10$	$Pe=100$	$Pe=1000$
$16$	0.0535674	0.0959928	n.c.	n.c.
$32$	0.0301618	0.060635	n.c.	n.c.
$64$	0.0152526	0.0333928	0.0590015	n.c.
$128$	0.00756073	0.0175604	0.0314838	n.c.
$256$	0.00376267	0.00901208	0.0161104	n.c.
Rate	1.006	0.957	0.963	n.c.

Table 4. Table 4: Convergence rates in ℓ 2 subscript ℓ 2 \ell_{2} -norm of the virtual finite volume approximation to u e ​ x ​ ( x , y ) = sin ⁡ 2 ​ π ​ x ​ sin ⁡ 2 ​ π ​ y subscript 𝑢 𝑒 𝑥 𝑥 𝑦 2 𝜋 𝑥 2 𝜋 𝑦 u_{ex}(x,y)=\sin 2\pi x\sin 2\pi y using the upwind advective flux reconstruction 𝝈 ~ → a ​ ( u h ) subscript → ~ 𝝈 𝑎 superscript 𝑢 ℎ \maxvec{\widetilde{\bm{\sigma}}}_{a}({u}^{h}) .

$N$	$Pe=1$	$Pe=10$	$Pe=100$	$Pe=1000$	$Pe=10000$
$16$	0.061851	0.11673	0.146463	0.153675	0.15467
$32$	0.031884	0.0633952	0.0911815	0.0971759	0.0979677
$64$	0.0155568	0.0339481	0.0550493	0.0595179	0.0600997
$128$	0.00762096	0.0177066	0.0304351	0.0331474	0.0337972
$256$	0.00377766	0.009056	0.0158462	0.0172527	0.0174534
Rate	1.02	0.967	0.944	0.944	0.954

Table 5. Table 5: Truncation errors in ℓ 2 subscript ℓ 2 \ell_{2} -norm and convergence rate of the MMD operator ( 39 ) implemented using the uniform ( 33 ) and non-uniform ( 34 ) virtual volume definitions.

$N$	${\mu }_i^{\text{U}}$	${\mu }_i^{\text{MR}}$
$16$	0.312637	0.267592
$32$	0.188097	0.145183
$64$	0.102245	0.076769
$128$	0.054184	0.040412
$256$	0.027495	0.020119
Rate	0.978711	1.00622

Equations

$${\int }_{\omega }d\mathbf{u}={\int }_{\partial \omega }\mathbf{u}.$$

$$d\mathbf{u}(\mathbf{x})=\underset{\mu (\omega )\to 0}{\lim }\frac{{\int }_{\partial \omega }\mathbf{u}}{\mu (\omega )},$$

$$\chi_{i}(\gamma)=\left\{\begin{array}[]{rl}\Gamma_{i}&\mbox{if $\bm{v}_{i}\in\gamma$}\\[2.15277pt] \emptyset&\mbox{otherwise}\end{array}\right.\quad\mbox{and}\quad\chi_{i}(\gamma^{\prime})=\left\{\begin{array}[]{rl}\Gamma_{i}&\mbox{if $\bm{v}_{i}\in\gamma^{\prime}$}\\[2.15277pt] \emptyset&\mbox{otherwise}\end{array}\right.,$$

$$h_X=\underset{\mathbf{x}\in \Omega }{\sup }\underset{{\mathbf{x}}_i\in X}{\min }|\mathbf{x}-{\mathbf{x}}_i|,$$

$$q_X=\frac{1}{2}\underset{i\ne j}{\min }|{\mathbf{x}}_i-{\mathbf{x}}_j|,$$

$$q_X\le h_X\le c_{qu}{q}_X.$$

$$S^h\ni u^h$$

$$U^h\ni {\mathbf{u}}^h$$

$$S_{\gamma }^h\ni u^h$$

$$U_{\gamma }^h\ni {\mathbf{u}}^h$$

$$\Vert u^h{\Vert }_{{\ell }^{\infty },X_{\omega }}:=\underset{{\mathbf{x}}_i\in X_{\omega }}{\max }|u_i^h|.$$

$$\Vert u^h{\Vert }_{S^h}:=\Vert u^h{\Vert }_{{\ell }^{\infty },X}\text{mbox}and\Vert u^h{\Vert }_{S_{\gamma }^h}:=\Vert u^h{\Vert }_{{\ell }^{\infty },X\setminus X_{\gamma }}$$

$$\{\varphi \in \Phi |{\lambda }_i(\varphi )=0,i=1,\dots ,n\}=\{0\}.$$

$$(\mathbf{\tau }(\varphi ){)}_i=\tau ({\varphi }_i),i=1,\dots ,q,$$

$$B_{ij}={\lambda }_i({\varphi }_j);i=1,\dots ,n;j=1,\dots ,q.$$

$$|\mathbf{b}{|}_{W(\tau )}^2={\mathbf{b}}^{⊺}W(\tau )\mathbf{b}\forall \mathbf{b}\in {\mathbb{R}}^n.$$

$$\tilde{\tau }(u):=\mathbf{c}(u)\cdot \mathbf{\tau }(\varphi ),$$

$$\mathbf{c}(u)=\underset{\mathbf{c}\in {\mathbb{R}}^q}{\mathrm{argmin}}\frac{1}{2}{\left|B\mathbf{c}-\mathbf{\lambda }(u)\right|}_{W(\tau )}^2$$

$$\mathbf{c}(u)={\left(B^{⊺}W(\tau )B\right)}^{-1}(B^{⊺}W(\tau ))\mathbf{\lambda }(u).$$

$${\mathbf{u}}^h=\{u_{ij}\in \mathbb{R}|{\mathbf{f}}_{ij}\in F\}\cup \{u_{{\Gamma }_i}\in \mathbb{R}|{\Gamma }_i\in {\gamma }^{\prime }\}\text{mbox}and{s}^h=\{{\varphi }_i\in \mathbb{R}|{\omega }_i\in C\},$$

$$V_i=\{{\mathbf{v}}_j\in V|{\mathbf{e}}_{ij}=({\mathbf{v}}_i,{\mathbf{v}}_j)\in E\}$$

$$\partial\omega_{i}=\Big{(}\bigcup_{j\in V_{i}}\bm{f}_{ij}\Big{)}\cup\chi_{i}(\gamma^{\prime})\cup\chi_{i}(\gamma)\,.$$

$$(DIV{\mathbf{u}}^h{)}_i:=\frac{1}{{\mu }_i}\left[\sum _{{\mathbf{f}}_{ij}\in \partial {\omega }_i}{u}_{ij}{\mu }_{ij}+u_{{\chi }_i({\gamma }^{\prime })}{\mu }_{{\chi }_i({\gamma }^{\prime })}+{\int }_{{\chi }_i(\gamma )}\mathbf{u}\cdot \mathbf{n}dS\right]\forall {\omega }_i\in C.$$

$$(DIV\mathcal{R}(\mathbf{u}){)}_i=\nabla \cdot \mathbf{u}({\mathbf{x}}_i)+O(h).$$

$$\omega =\bigcup _{i\in C_{\omega }}{\omega }_i.$$

$$\sum _{i\in C_{\omega }}{\mu }_i(DIV{\mathbf{u}}^h{)}_i=\sum _{{\mathbf{f}}_{ij}\in \partial {\omega }_0}{u}_{ij}{\mu }_{ij}+\sum _{{\Gamma }_k\in \partial {\omega }_{{\gamma }^{\prime }}}{u}_{{\Gamma }_k}{\mu }_{{\Gamma }_k}+\sum _{{\Gamma }_k\in \partial {\omega }_{\gamma }}{\int }_{{\Gamma }_k}\mathbf{u}\cdot {\mathbf{n}}_{k}dS.$$

$$\begin{array}[]{rcl}\displaystyle\sum_{i\in C_{\omega}}\mu_{i}(DIV\bm{u}^{h})_{i}&=&\displaystyle\sum_{i\in C_{\omega}}\left[\sum_{\bm{f}_{ij}\in\partial\omega_{i}}u_{ij}\mu_{ij}+u_{\chi_{i}(\gamma^{\prime})}\mu_{\chi_{i}(\gamma^{\prime})}+\int_{\chi_{i}(\gamma)}\bm{u}\cdot\bm{n}dS\right]\\[12.91663pt] &=&\displaystyle\sum_{i\in C_{\omega}}\sum_{\bm{f}_{ij}\in\partial\omega_{i}}u_{ij}\mu_{ij}+\sum_{\Gamma_{k}\in\partial\omega_{\gamma^{\prime}}}u_{\Gamma_{k}}\mu_{\Gamma_{k}}+\sum_{\Gamma_{k}\in\partial\omega_{\gamma}}\int_{\Gamma_{k}}\bm{u}\cdot\bm{n}_{k}dS\,.\end{array}$$

$$\sum _{{\omega }_i\in C}{\mu }_i(DIV{\mathbf{u}}^h{)}_i={\int }_{\Gamma }\mathbf{u}\cdot \mathbf{n}dS.$$

$$\bm{e}_{ij}=\left\{(ij)\,\big{|}\,(ij)\in V\times V\right\}\,.$$

$$\partial\omega_{i}=\left\{\begin{array}[]{ll}\{\bm{f}_{ij}\},\ j\in V_{i}&\mbox{if $i\in V_{0}$}\\[8.61108pt] \{\bm{f}_{ij}\}\cup\Gamma_{i},\ j\in V_{i}&\mbox{if $i\in V_{\Gamma}$}\end{array}\right.;\quad\forall\omega_{i}\in C.$$

$$\partial\omega_{i}=\Big{(}\bigcup_{j\in V_{i}}\bm{f}_{ij}\Big{)}\cup\chi_{i}(\gamma^{\prime})\cup\chi_{i}(\gamma)\,.$$

Full text

††journal: Journalsand-blurbsand-blurbfootnotetext: Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA-0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. CCRCCRfootnotetext: Center for Computing Research, Sandia National Laboratories, Albuquerque, NM 87185

url]http://www.sandia.gov/ natrask/

url]http://www.sandia.gov/ pbboche/

url]http://www.sandia.gov/ mperego/

1 Introduction

The vast majority of numerical methods for PDEs rely on partitioning of the spatial domain into a mesh to both represent the unknown fields and to define the discrete operators acting on these fields. The quality of the resulting mesh-based schemes usually depends on the quality of the underlying mesh and may suffer when the latter deteriorates. For example, shape functions defined on “sliver” elements can result in badly conditioned or even singular stiffness matrices; see [1] and the more recent survey [2].

Accordingly, automated generation of high-quality grids is a task of tremendous practical importance for mesh-based PDE discretizations. Yet, for problems with complex geometries this task remains a major challenge, and creates a significant performance bottleneck in the simulation workflow [3]. As future applications continue to evolve beyond forward simulation, overcoming this bottleneck will be of essence for enabling tasks such as shape optimization and uncertainty quantification. For applications where the domain evolves in time, maintaining a high-quality mesh is even more challenging, typically requiring the introduction of remap/remeshing algorithms as elements collapse and tangle under Lagrangian motion. This situation is typical of continuum mechanics problems characterized by large deformation and topology change, such as free-surface flows, hypervelocity impact, and ductile fracture.

Meshfree methods can significantly alleviate or even completely eliminate the mesh generation bottleneck. For example, the domain integration in Galerkin meshfree methods can be performed on substandard meshes, by overlaying the domain with a mesh that is independent of the point cloud [4], or even by a non-conforming set of cells as in the SNNI family of schemes [5]. Collocation meshfree methods on the other hand only require a formation of well-distributed point clouds, which is much easier than generation of high-quality grids, especially for complex geometric domains.

However, since their inception, meshfree methods have struggled to simultaneously maintain rigorous notions of consistency, accuracy, stability, and compatibility that are now commonplace in mesh-based approaches. While significant improvements have been made in recent years addressing the approximation theory of meshfree methods, a meshfree framework mirroring the conservation and consistency properties of mesh-based compatible discretizations has remained elusive. The root cause is that the topological and geometric data employed by typical compatible mesh-based discretizations such as mimetic finite differences (MFD) [6], Whitney elements [7], and finite volumes [8, 9, 10] cannot be assumed readily available in the meshfree setting, especially in the absence of a conforming background mesh. To successfully transplant the properties of these methods to such a setting it is important to recognize that their compatibility is a consequence of discrete differentiation based on the generalized Stokes theorem

[TABLE]

In (1) $d$ is the exterior derivative, $\bm{u}$ is a $k$-differential form and $\omega$ is a $k+1$ dimensional manifold. This theorem implies that for a sufficiently smooth $\bm{u}$ its derivative at a point $\bm{x}$ can be characterized as

[TABLE]

where $\omega$ is a $k+1$-dimensional region with measure $\mu(\omega)$ containing the point $\bm{x}$. Restriction of the right hand side in (2) to a $k+1$-dimensional mesh entity $\omega_{h}$ then yields a notion of a discrete derivative that satisfies a discrete Stokes theorem by construction. In particular, the resulting discrete grad, curl and div operators mimic the vector calculus identities $curl\,grad\,\bm{u}=0$ and $div\,curl\,\bm{u}=0$; see, e.g., [11].

If one wishes to use this procedure to define the discrete derivatives, then the underlying discretization infrastructure must provide the following three key pieces of information:

field data representing the integrand $\bm{u}$;

• 2.

topological data describing the action of the boundary operator on a geometric entity $\omega$, and

• 3.

metric data providing the measures of the geometric entity $\omega$ and its boundary $\partial\omega$.

The first one, i.e., the field data, can be presumed available irregardless of the type of the underlying discretization as both mesh-based and meshfree methods should at a minimum be able to represent the approximate PDE solutions. The second and third information pieces are virtually always available in mesh-based methods. The nodes, edges, faces and cells in most grids form a chain complex, which provides the necessary topological data for compatible discrete grad, curl and div operators; see, e.g., [12, 13]. Likewise, the necessary metric data can be calculated from the mesh description without much difficulty.

This is not the case when the spatial domain is represented by a point cloud. The amount of metric and topological information that can be extracted from the cloud at a reasonable cost is quite limited. One can build, at a linear cost, the $\epsilon$-ball graph of the cloud, construct its vertex-to-edge incidence matrix and compute the lengths of its edges. This would provide (2) with just enough information to build a compatible discrete gradient [12] on the point cloud but not much else.

Here we develop a computationally efficient and scalable approach to generate the topological and metric data for (2) needed to define a compatible divergence operator on point clouds. Using as a template a mimetic divergence on a primal-dual mesh, we employ a virtual primal-dual mesh to construct an abstract meshfree version of this operator by contraction of field moments, associated with the virtual faces, and face moments, which characterize the latter algebraically.

We consider two instantiations of this abstract mimetic meshfree divergence (MMD) operator. The first one assumes a background mesh and uses generalized moving least squares (GMLS) [14] to obtain the necessary field and face moments. This MMD instance is appropriate for settings where a background mesh exists but its quality is insufficient for a robust and accurate mesh-based discretization. In this case the mesh is used only to define a boundary operator and to integrate the local GMLS polynomial basis over the cell faces. Both of these tasks can be performed reliably even on substandard grids.

The second MMD operator does not assume a background mesh and uses instead the $\epsilon$-ball graph of the point cloud and its formal dual as surrogates for a primal and a virtual dual mesh, respectively. The face-to-cell incidence matrix of the latter provides the necessary topological information for the divergence operator. This MMD instance retains the GMLS field moments but defines the virtual face moments in terms of virtual area potentials, leading to a graph Laplacian problem that can be solved in an efficient and scalable manner by standard multigrid preconditioners. In so doing we trade a challenging mesh generation problem for a scalable algebraic one, without sacrificing compatibility with the divergence operator.

In [15] we used the local connectivity graph of each particle and its formal dual to mimic the staggered arrangement in div-grad stencils on primal-dual grids. The resulting “staggered” meshfree scheme behaves similarly to mesh-based div-compatible discretizations but is not locally conservative because the primal-dual mesh surrogate is defined independently for each particle. This paper is a logical extension of [15] aiming to deliver a meshfree scheme mirroring the conservation properties of traditional compatible discretizations.

While there is an abundant literature on compatible mesh-based discretizations, this is not the case for meshfree methods. To the best of our knowledge, the extant work comprises the meshless volume scheme [16], the Uncertain Grid Method (UGM) [17], the meshfree framework for conservation laws in [18], and a few related references. Of particular interest to us are [17] and [18] because, as in the second MMD instance, they define the necessary metric data in a purely algebraic fashion without assuming a background grid.

The UGM appears to be the first example of a locally conservative mesfree scheme, which uses virtual geometric entities characterized algebraically by numbers representing their “measures”. In UGM, the term “uncertain” refers to an “uncertain face” between two neighboring particles, which is similar to our virtual faces. The meshfree framework in [18] also seeks sets of numbers that can be interpreted as measures of virtual volumes, centered at each particle, and the areas of their virtual faces, respectively. The single most important difference between these works and our approach is in the type of algebraic problems that generate the necessary metric data. Both [17] and [18] find this data by solving global constrained optimization problems. In UGM this problem is a linear program solved by a primal-dual logarithmic barrier method, whereas in [18] it is a quadratic program (QP) requiring a QP-specific solver. In contrast, our approach requires solution of several graph Laplacian problems, which can be accomplished in linear time.

We have organized the paper as follows. Section 2 introduces notation and reviews the necessary technical background. Section 3 formulates the abstract MMD operator and states general conditions for its consistency and accuracy, while Section 4 presents its two instances. Efficient and scalable computation of the virtual face moments necessary for the second MMD instance is discussed in Section 4.2.3. In Section 5 we use the MMD operator to discretize a model conservation law problem and then in Section 6 we provide some representative numerical examples. We summarize our findings and discuss future work in Section 7.

2 Preliminaries

2.1 Notation

Throughout the paper upper case fonts are reserved for function spaces, operators and sets of various entities, while lower case fonts stand for scalar fields, linear functionals, indices, etc. We denote the standard Euclidean norm on $\mathbb{R}^{n}$ by $|\cdot|$ and use bold face fonts to denote vector quantities, e.g., $\bm{x}=(x_{1},\ldots,x_{d})$ is a point in the Euclidean space $\mathbb{R}^{d}$, $\bm{e}_{ij}=(\bm{x}_{i},\bm{x}_{j})$ is an edge connecting two such points, $\bm{u}=(u_{1},\ldots,u_{d})$ is a vector field in $\mathbb{R}^{d}$, and $\bm{n}=(n_{1},\ldots,n_{d})$ is a unit normal vector, i.e., $|\bm{n}|=1$.

As usual, $C^{k}(\Omega)$ is the space of all $k$-continuously differentiable functions, and $P_{m}(\mathbb{R}^{d})$, or simply $P_{m}$ is the space of all multivariate polynomials of degree less than or equal to $m$. We denote the standard norm on $C^{k}(\Omega)$ and its restriction to $\omega\subset\Omega$ by $\|\cdot\|_{C^{k}(\Omega)}$ and $\|\cdot\|_{C^{k}(\omega)},$ respectively.

Let $\Omega\subset\mathbb{R}^{d}$, $d=2,3$ be a bounded connected region with a Lipschitz continuous boundary $\Gamma=\partial\Omega$. We assume that $\Gamma$ may be partitioned into a collection of $p_{\Gamma}$ faces $\Gamma=\cup\Gamma_{j}$, each having maximal diameter $h$, measure $\mu_{\Gamma_{j}}$, centroid $\bm{x}_{j}$ and outward unit normal $\bm{n}_{j}$. The collection of these centroids forms the set ${X}_{\Gamma}$ of all boundary particles. In the interior of the domain, we assume a set ${X}_{0}$ containing $p_{0}$ interior particles. The union ${X}={X}_{0}\cup{X}_{\Gamma}$ has $p_{\Omega}=p_{0}+p_{\Gamma}$ particles and defines the meshfree discretization of $\Omega$. Since by assumption $\Omega$ is bounded we can choose the coordinate system in such a way that the point cloud is strictly contained in the first quadrant, i.e., ${x}_{k}>0$, $k=1,\ldots,d$ for all $\bm{x}\in{X}$. We denote the subset of $X$ contained in an entity $\omega$ by $X_{\omega}:=\omega\cap X$. Whenever required, boundary conditions will be imposed on simply connected non-empty subsets $\gamma\subseteq\Gamma$. Without loss of generality we shall assume that $\gamma$ and its complement $\gamma^{\prime}:=\Gamma\setminus\gamma$ are exact unions of boundary segments $\Gamma_{j}$ and define their indicator functions as

[TABLE]

respectively. The quality of the point cloud ${X}$ can be characterized by its fill distance [14, p.14] given by

[TABLE]

and its separation distance [14, p.41] defined as

[TABLE]

We shall assume that ${X}$ is quasi-uniform, namely that there exists $c_{qu}>0$ such that

[TABLE]

We denote the spaces of all discrete scalar and vector fields on $X$ by $S^{h}$ and $U^{h}$, respectively. These fields represent functions by their point samples on $X$, which we group as

[TABLE]

respectively. Note that $\bm{u}^{h}=(u^{h}_{1},\ldots,u^{h}_{d})$, where $u^{h}_{k}\in S^{h}$, $k=1,2,\ldots,d$. We denote by $S^{h}_{\gamma}$ and $U^{h}_{\gamma}$ the sets of all discrete scalar and vector fields defined by their point values in $X_{0}\cup X_{\gamma^{\prime}}$, i.e.,

[TABLE]

We use $S^{h}_{\gamma}$ or $U^{h}_{\gamma}$ whenever a scalar field $u(\bm{x})$ or a vector field $\bm{u}(\bm{x})$ is prescribed on $\gamma$. Let

[TABLE]

We equip $S^{h}$ and $S^{h}_{\gamma}$ with the norms

[TABLE]

respectively, and similarly for $U^{h}$ and $U^{h}_{\gamma}$.

2.2 GMLS essentials

Generalized moving least squares (GMLS) is a non-parametric regression technique for the approximation of bounded linear functionals from scattered samples of their arguments [14]. GMLS has been used in meshfree collocation schemes to approximate point values of derivatives; see [19], [20], [15] and the references therein. Here we shall apply GMLS to approximate integrals as our goal is to develop a meshfree divergence operator based on the definition (2). We provide a brief summary of the GMLS framework, specialized to our needs. The abstract GMLS setting comprises the following key ingredients; see [14, Section 4.3]:

a function space $U$ with a dual $U^{*}$;

• 2.

a finite dimensional space $\Phi=\mbox{span}\{\phi_{1},\ldots,\phi_{q}\}\subset U$;

• 3.

a finite set of linear functionals $\Lambda=\{\lambda_{1},\ldots,\lambda_{n}\}\subset U^{*}$; and

• 4.

a correlation (kernel) function $w:U^{*}\times U^{*}\mapsto\mathbb{R}^{+}\cup\{0\}$.

The set $\Lambda$ is assumed to be $\Phi$-unisolvent, that is

[TABLE]

Given a target functional $\tau\in U^{*}$ and an arbitrary $u\in U$, GMLS seeks an approximation $\widetilde{\tau}(u)$ of $\tau(u)$ in terms of the sample set $\bm{\lambda}(u)=(\lambda_{1}(u),\ldots,\lambda_{n}(u))\in\mathbb{R}^{n}$, such that $\widetilde{\tau}(\phi)=\tau(\phi)$ for all $\phi\in\Phi$. We call such a GMLS approximation $\Phi$-reproducing. To describe the GMLS solution to this problem we introduce the vector $\bm{\tau}(\phi)\in\mathbb{R}^{q}$ with elements

[TABLE]

the diagonal weight matrix $W(\tau)\in\mathbb{R}^{n\times n}$ with element $W_{ii}(\tau)=w(\tau;\lambda_{i})$, and the basis sample matrix $B\in\mathbb{R}^{n\times q}$ with element

[TABLE]

In what follows, $|\cdot|_{W(\tau)}$ denotes the Euclidean norm on $\mathbb{R}^{n}$ weighted by $W(\tau)$, i.e.,

[TABLE]

One can then show that the $\Phi$-reproducing GMLS approximation of the target functional is given by

[TABLE]

where the GMLS coefficient vector $\bm{c}(u)\in\mathbb{R}^{q}$ solves the weighted least-squares problem

[TABLE]

It is not hard to see that

[TABLE]

We refer to previous work [15] for details regarding the efficient and stable calculation of $\bm{c}(u)$.

2.3 Mimetic divergence operator

In Section 3 we define an abstract meshfree divergence operator by mimicking a mimetic divergence operator on a primal-dual mesh [21, 10, 22]. Remark 2 explains the appropriateness of this choice. To describe the template operator assume that $\Omega$ is equipped with discretization infrastructure comprising

a primal mesh $\tau_{h}(V,E)$ with vertices $V=\{\bm{v}_{i}\}_{i=1}^{p_{V}}$ and edges $E=\{\bm{e}_{ij}\}$; $\bm{e}_{ij}=(\bm{v}_{i},\bm{v}_{j})\in V\times V$,

• 2.

a topologically dual mesh $\tau_{h}^{\prime}(C,F)$ with cells $C=\{\omega_{i}\}_{i=1}^{p_{V}}$ and affine faces $F=\{\bm{f}_{ij}\}$, and

• 3.

a pair of discrete spaces $\widetilde{U}^{h}_{\gamma}$ and $\widetilde{S}^{h}$ with elements

[TABLE]

respectively, representing the domain and the range of the discrete divergence, respectively.

Every dual cell $\omega_{i}$ corresponds to a primal vertex $\bm{v}_{i}$ and has a measure $\mu_{i}$. We assume that every dual face $\bm{f}_{ij}$ corresponds to a primal edge $\bm{e}_{ij}$, has a unit normal $\bm{n}_{ij}=\bm{e}_{ij}/|\bm{e}_{ij}|$, and oriented measure $\mu_{ij}=\int_{\bm{f}_{ij}}dS$. If $\omega_{i}$ and $\omega_{j}$ are the dual cells corresponding to the endpoints of this edge, then $\bm{f}_{ij}=\partial\omega_{i}\cap\partial\omega_{j}$. The space $\widetilde{U}^{h}_{\gamma}$ approximates vector fields in H^{0}_{\gamma}(div,\Omega)=\{\bm{u}\in H(div,\Omega)\cap(C^{0}(\Omega))^{d}\,|\,\mbox{\bm{u}=\bm{u}_{\gamma}$on$\gamma}\} by their scalar projections onto $\bm{n}_{ij}$ or $\bm{n}_{ji}$, whereas $\widetilde{S}^{h}$ approximates $L^{2}(\Omega)$ functions by their dual cell averages. Without a loss of generality we shall assume that the restriction $\mathcal{R}(\bm{u})\in\widetilde{U}^{h}_{\gamma}$ of a vector field $\bm{u}$ is given by $u_{ij}=u_{ji}=\bm{u}\cdot\bm{n}_{ij}$ for $i<j$. The neighborhood of $\bm{v}_{i}$ is defined as

[TABLE]

and contains all mesh vertices connected to $\bm{v}_{i}$ by an edge. If a vertex $\bm{v}_{i}\in\Gamma$ then the boundary of its dual cell $\omega_{i}$ has a non-empty intersection with $\Gamma$, which we denote as $\Gamma_{i}:=\partial\omega_{i}\cap\Gamma$. As in Section 2.1 we shall always assume that $\gamma$ and $\gamma^{\prime}$ are exact unions of such segments. It follows that

[TABLE]

Restriction of (2) to dual cells $\omega_{i}$, followed by approximation of the integrals on $\partial\omega_{i}$ by single point quadrature yields a discrete divergence operator $DIV:\widetilde{U}^{h}_{\gamma}\rightarrow\widetilde{S}^{h}$, whose action on $\bm{u}^{h}\in\widetilde{U}^{h}_{\gamma}$ is given by

[TABLE]

In this operator

the field data is provided by the degrees-of-freedom in $\widetilde{U}^{h}_{\gamma}$, i.e., the sets of numbers $\{u_{ij}\}$ and $\{u_{\Gamma_{i}}\}$, such that $u_{ij}=u_{ji}$;

• 2.

the metric data is provided by the cell volumes and the oriented face areas in the dual mesh, i.e., the sets of numbers $\{\mu_{i}\}$, $\{\mu_{ij}\}$ and $\{\mu_{\Gamma_{i}}\}$, respectively, such that $\mu_{ij}=-\mu_{ji}$;

• 3.

the topological data is provided by formula (10).

Remark 1**.**

In mixed discretizations [23] of, e.g., diffusion problems , $\{u_{ij}\}$ and $\{u_{\Gamma_{i}}\}$ are degrees-of-freedom describing an $H(div)$-conforming discretization of a flux field. In contrast, in finite volume discretizations of conservation laws [24] the sets $\{u_{ij}\}$ and $\{u_{\Gamma_{i}}\}$ represent values derived from another set of degrees-of-freedom, usually describing a scalar dependent variable such as a pressure field.

One can show that (11) is first-order accurate, i.e., for all sufficiently regular vector fields $\bm{u}$ there holds

[TABLE]

The following result confirms that (11) satisifes a discrete divergence theorem, i.e., that it is compatible.

Lemma 1**.**

Consider a collection $C_{\omega}$ of dual cell indices and define $\omega\subseteq\Omega$ as

[TABLE]

Let $\partial\omega_{0}=\partial\omega\cap\Omega$, $\partial\omega_{\gamma}=\partial\omega\cap\gamma$, and $\partial\omega_{\gamma^{\prime}}=\partial\omega\cap\gamma^{\prime}$. Then,

[TABLE]

Proof.

Summing (11) over the cells in $\omega$ gives

[TABLE]

Since $\mu_{ij}=-\mu_{ji}$ all terms in the double sum above corresponding to pairs of adjacent cells cancel each other, leaving only the terms for which $\bm{f}_{ij}\in\partial\omega_{0}$. This completes the proof. ∎

Global conservation in $\Omega$ a direct consequence of Lemma 1. For simplicity we state it with $\gamma=\Gamma$.

Corollary 1**.**

Assume that $\bm{u}^{h}\in\widetilde{U}^{h}_{\Gamma}$ and let $\bm{u}\in H_{\Gamma}(div,\Omega)$ be given. Then,

[TABLE]

Proof.

The result follows by setting $\omega=\Omega$ in Lemma 1 and noting that $\partial\omega_{0}=\emptyset$, $\partial\omega_{\gamma^{\prime}}=\emptyset$, and $\partial\omega_{\gamma}=\Gamma$. ∎

3 An abstract meshfree mimetic divergence operator

We consider a meshfree discretization infrastructure comprising

a point cloud $X=X_{0}\cup X_{\Gamma}$, representing the computational domain $\Omega$, and

• 2.

discrete spaces $S^{h}$ and $U^{h}_{\gamma}$, representing scalar and vector fields by their point samples on the cloud.

Our goal is to endow this infrastructure with an abstract meshfree mimetic divergence (MMD) operator whose properties mirror those of its mesh-based prototype (11). Thus, we seek a mapping ${DIV}:U^{h}_{\gamma}\rightarrow S^{h}$ that is first-order accurate and satisfies a discrete divergence theorem. We shall define this mapping by mimicking the action of (11) which assigns to every dual cell $\omega_{i}$ a contraction of the field and metric data living on the faces in $\partial\omega_{i}$, and weighted by the reciprocal of the cell’s measure. This task requires appropriate abstractions for the field and metric data, as well as a suitable notion of a boundary operator.

We shall define this operator through a virtual primal-dual mesh complex. To reduce notational clutter we reuse the nomenclature from Section 2.3. We call a collection of integer indices $V=\{i\}_{i=1}^{p_{V}}$ a virtual vertex set and refer to its elements as $\bm{v}_{i}$. A virtual edge is simply an ordered pair of virtual vertices, i.e.,

[TABLE]

A collection of $p_{E}$ virtual edges forms a set $E$, which together with $V$ defines the virtual primal mesh $\tau_{h}(V,E)$. To every $\bm{v}_{i}\in V$ and $\bm{e}_{ij}\in E$ we assign a virtual dual cell $\omega_{i}$ and a virtual dual face $\bm{f}_{ij}$, respectively. The sets of these dual entities are denoted by $C$ and $F$, respectively, and they form a virtual topologically dual mesh $\tau^{\prime}_{h}(C,F)$. A practically useful virtual boundary operator $\partial C\rightarrow F$ should be able to recover the physical domain boundary $\Gamma$. This requires some degree of association between the virtual mesh structures and the physical domain $\Omega$. A simple and effective way to establish such an association is to tie every virtual vertex to a point in the cloud $X$. This induces a partition of $V$ into interior vertices $V_{0}=X_{0}$ and boundary vertices $V_{\Gamma}=X_{\Gamma}$, and prompts the following definition:

[TABLE]

This definition allows us to write the virtual cell boundary in the same form as in (10), i.e.,

[TABLE]

Note that while $\omega_{i}$ is a purely virtual entity, its boundary may contain “real” geometric entities from $\Gamma$.

Consider next the field data abstraction. To mimic (11) the domain of the MMD operator should be associated with the virtual faces $F$, whereas its range should contain fields living on the virtual cells $C$. To ensure first-order accuracy the latter should be piecewise constant with respect to $C$. Thus, we can reuse the definition of $\widetilde{S}^{h}$ from (9) with the understanding that each degree-of-freedom is now assigned to a virtual, rather than to a physical, dual cell. Likewise, the domain space should be able to reproduce exactly the preimage of $\widetilde{S}^{h}$, i.e., any vector field $\bm{u}$ such that $\nabla\cdot\bm{u}\in\widetilde{S}^{h}$. In the mesh-based case this can be accomplished by a single111For example, the lowest-order Raviart-Thomas elements on $d$-simplices [25], also known as Whitney 2-forms [7], are incomplete linear polynomial fields of the form $\bm{a}+b\bm{x}$ which have $d+1$ degrees-of-freedom. degree-of-freedom per face. However, this may not be enough in the meshfree setting and so we shall allow for $n_{F}\geq 1$ field moments per virtual face $\bm{f}_{ij}\in F$ and boundary segment $\Gamma_{i}\in\gamma$. In sum, we define the domain $\widetilde{U}^{h}_{\gamma}$ and the range $\widetilde{S}^{h}$ for the abstract MMD operator by modifying (9) to

[TABLE]

respectively, with the understanding that all dual mesh entities in (17) are virtual.

Finally, let us consider the abstraction for the metric data. We seek this abstraction in the form of virtual face and cell moments, i.e., sets of numbers assigned to every element of $F$ and $C$. To allow contraction of these moments with the field moments in (17) we chose the former to be the algebraic duals of the latter. Thus, we assign a real number $\mu_{i}$ to every virtual dual cell $\omega_{i}\in C$, a real vector $\bm{\mu}_{ij}\in\mathbb{R}^{n_{F}}$ to every virtual dual face $\bm{f}_{ij}\in F$ and a real vector $\bm{\mu}_{\Gamma_{i}}\in\mathbb{R}^{n_{F}}$ to every boundary segment $\Gamma_{i}\in\gamma^{\prime}$. Following (11) we then define an abstract operator $DIV:\widetilde{U}^{h}_{\gamma}\rightarrow\widetilde{S}^{h}$ by contracting the metric and the field data:

[TABLE]

This operator acts on the field moments $\{\bm{u}_{ij}\}$ and $\{\bm{u}_{\Gamma_{i}}\}$ rather than on the point samples $U^{h}_{\gamma}$ that are the “native” representation of vector fields on the point cloud. Thus, deployment of (18) in our meshfree infrastructure requires a particle-to-virtual face data transfer operator222Formally we also need a data transfer operator $T_{F}:\widetilde{S}^{h}\rightarrow S^{h}$ to map the discrete divergence back to the cloud points. However, since in the present context $S^{h}\equiv\widetilde{S}^{h}$, this operator is the identity and is omitted for simplicity.

[TABLE]

which maps point samples into field moments. Combining the action of this operator with (18) yields the desired abstract MMD operator $DIV:U^{h}_{\gamma}\rightarrow S^{h}$, where

[TABLE]

The following section studies the properties of (20).

3.1 Analysis of the abstract MMD operator

To ensure that (20) has the same mimetic properties as its mesh-based prototype we require that

T.1

The virtual face volumes satisfy $\mu_{i}>0$, $\mu_{i}=O(h^{d}_{X})$, and $\sum_{i}\mu_{i}=\mu(\Omega)$.

T.2

The virtual face moments $\{\bm{\mu}_{ij}\}$ are antisymmetric: $\bm{\mu}_{ij}=-\bm{\mu}_{ji}$.

T.3

The data transfer operator (19) is symmetric: $\bm{t}_{ij}(\bm{u}^{h})=\bm{t}_{ji}(\bm{u}^{h}).$

T.1 ensures that the virtual volumes behave like physical volumes. The “topological” requirements T.2-T.3 are prompted by Lemma 1, which reveals that mimetic properties of (11) hinge on the symmetry of the face degrees-of-freedom $u_{ij}$ and the antisymmetry of the oriented face areas $\mu_{ij}$. The following result shows that T.1-T.3 are indeed sufficient for the abstract MMD operator to be locally and globally conservative.

Theorem 1**.**

Assume that T.1-T.3 hold. Then the abstract MMD operator (20) satisfies (13) (local conservation) and (14) (global conservation).

Proof.

Owing to T.2 and T.3 the proof of Lemma 1 applies verbatim to the abstract MMD operator as well. The global conservation (14) then follows from T.1 and (13) as in Corollary 1. ∎

While conservation is a purely “topological” property that can be secured by independent conditions on the field and face moments, this is not the case with the consistency of the MMD operator, which requires these moments to work well together. The following two definitions introduce conditions that coordinate the properties of the metric and the field data to ensure the desired first-order accuracy.

Definition 1**.**

A set of virtual cell volumes and virtual face moments $\{\mu_{i},\bm{\mu}_{ij},\bm{\mu}_{\Gamma_{i}}\}$ and a data transfer operator $T_{F}:U^{h}_{\gamma}\rightarrow\widetilde{U}^{h}_{\gamma}$ are called a $P_{1}$-reproducing pair if for every $\bm{p}\in(P_{1})^{d}$ and its point sample $\bm{p}^{h}\in U^{h}_{\gamma}$ there holds

[TABLE]

In other words, $\{T_{F};\{\mu_{i},\bm{\mu}_{ij},\bm{\mu}_{\Gamma_{i}}\}\}$ is a $P_{1}$-reproducing pair iff for every $\bm{p}\in(P_{1})^{d}$ there holds

[TABLE]

Definition 2**.**

A set of virtual face moments $\{\bm{\mu}_{ij},\bm{\mu}_{\Gamma_{i}}\}$ and a data transfer operator $T_{F}:U^{h}_{\gamma}\rightarrow\widetilde{U}^{h}_{\gamma}$ are called a locally Lipschitz-continuous pair on $X$ if there exist balls $\{B_{\bm{f}}^{h}\}$ of radius $O(h_{X})$, centered at the barycenters of the virtual faces $\bm{f}\in F$ and a constant $C>0$, such that for every $\bm{u},\bm{v}\in(C^{2}(\Omega))^{d}$ with point samples $\bm{u}^{h},\bm{v}^{h}\in U^{h}_{\gamma}$ there holds

[TABLE]

These two properties are sufficient for the abstract MMD operator to be first-order accurate.

Theorem 2**.**

Assume that $\{T_{F};\{\mu_{i},\bm{\mu}_{ij},\bm{\mu}_{\Gamma_{i}}\}\}$ is a $P_{1}$-reproducing pair, $\{T_{F},\{\bm{\mu}_{ij},\bm{\mu}_{\Gamma_{i}}\}\}$ is a locally Lipschitz continuous pair on $X$, and that T.1 holds. Then, there exists a constant $C$, independent of the point cloud and such that for every $\bm{u}\in(C^{2}(\Omega))^{d}$ and its point sample $\bm{u}^{h}\in U^{h}_{\gamma}$ there holds the error bound

[TABLE]

Proof.

Consider a linear vector field $\bm{p}\in(P_{1})^{d}$ with a point sample $\bm{p}^{h}\in U^{h}_{\gamma}$, and an arbitrary point $\bm{x}_{i}\in X$. From the triangle inequality and the $P_{1}$-reproducing pair assumption it follows that

[TABLE]

Consider first the case $\chi_{i}(\gamma^{\prime})=\emptyset$. Using (22) then yields the following bound for the second term:

[TABLE]

where $\mu(\chi_{i}(\gamma))$ is the measure of $\chi_{i}(\gamma)$. Assumption T.1 implies that $\mu_{i}=O(h_{X}^{d})$, while $\mu(\chi_{i}(\gamma))$ is either zero or $O(h_{X}^{d-1})$. As a result, we have that

[TABLE]

Taking $\bm{p}$ to be the linear Taylor polynomial of $\bm{u}$ at the point $\bm{x}_{i}$ then yields the bound

[TABLE]

for the second term and the bound

[TABLE]

for the first term. This completes the proof for $\chi_{i}(\gamma^{\prime})=\emptyset$. When $\chi_{i}(\gamma^{\prime})=\Gamma_{i}$ the extra term $\|\bm{p}^{h}-\bm{u}^{h}\|_{\ell^{\infty},X_{B^{h}_{i}}}$ can be estimated by $\|\bm{p}-\bm{u}\|_{C^{0}(B^{h}_{i})}$. The proof follows by combining this and the above bounds. ∎

4 Meshfree mimetic divergence (MMD) instantiations

In this section we present two instances of the abstract MMD operator. Section 4.1 implements (20) using a background mesh. This hybrid formulation is useful for applications where such a grid is available but its quality is poor. For example, the mesh may contain valid, but almost degenerate, “sliver” elements that render the mesh-based shape functions nearly singular. In this case, a meshfree approach such as (20) can provide a more robust and accurate alternative to traditional mesh-based discretizations by using the mesh only to supply the metric and topological information but not to define the field data and the operators acting on it. Section 4.2 shows how to implement (20) without a background mesh. This MMD instance targets applications where mesh generation creates a computational bottleneck and should be avoided.

4.1 MMD with a background mesh

We assume that the virtual primal-dual mesh complex can be associated with a conforming physical grid on $\Omega$. Existence of a physical dual mesh implies that a set $\{\mu_{i}\}$ satisfying T.1 is readily available. Thus, to instantiate (20) it remains to find a data transfer operator $T_{F}$ and virtual face moments $\{\bm{\mu}_{ij},\bm{\mu}_{\Gamma_{i}}\}$ such that assumptions T.2–T.3 hold, $\{T_{F};\{\mu_{i},\bm{\mu}_{ij},\bm{\mu}_{\Gamma_{i}}\}\}$ is a $P_{1}$-reproducing pair, and $\{T_{F},\{\bm{\mu}_{ij},\bm{\mu}_{\Gamma_{i}}\}\}$ is a locally Lipschitz continuous pair on $X$. To that end we shall use GMLS with targets

[TABLE]

where $\bm{f}$ is a dual face or a segment on $\gamma^{\prime}$, i.e., $\bm{f}\in F\cup\{\Gamma_{i}\,|\,\bm{x}_{i}\in\gamma^{\prime}\}$. We apply GMLS with

a function space $U=(C^{1}(\Omega))^{d}$;

• 2.

a reproduction space $\Phi=(P_{1})^{d}$ with dimension $n_{1}=d(d+1)$ and basis $\{\bm{\phi}_{s}\}_{s=1}^{n_{1}}$ with basis functions

[TABLE]

where $\mathbf{e}_{k}$ is the $k$th canonical basis vector in $\mathbb{R}^{d}$ and $\{\phi_{r}(\bm{x})\}_{r=1}^{d+1}$ is a basis for $P_{1}$;

• 3.

sampling functionals333We refer to [14, Chapter 3] for conditions on the sampling locations $\{\bm{x}_{i}\}$ that are sufficient to ensure polynomial reproduction and rigorous approximation error estimates. $\lambda_{i}=\delta_{\bm{x}_{i}}$ mapping vector fields $\bm{u}\in U$ into point samples $\bm{u}^{h}\in U^{h}_{\gamma}$ on the cloud;

• 4.

kernels $w(\tau_{\bm{f}};\lambda_{i})=\kappa(|\bm{x}_{\bm{f}}-\bm{x}_{i}|)$ where $\bm{x}_{\bm{f}}$ is the centroid of $\bm{f}\in F\cup\{\Gamma_{j}\,|\,\bm{x}_{j}\in\gamma^{\prime}\}$ and $\kappa(\rho)=\left(1-{\rho}/{\epsilon_{\kappa}}\right)^{4}_{+}$, with $\epsilon_{\kappa}$ is large enough to obtain $P_{1}$-unisolvency. The face centroids $\{\bm{x}_{\bm{f}}\}$ are not required to coincide with the midpoints $\{\bm{x}_{ij}\}$ of the corresponding primal edges.

With these choices the GMLS approximation of (24) is exact for linear polynomials and is given by

[TABLE]

The GMLS coefficients $\bm{c}_{\bm{f}}(\bm{u}^{h})\in\mathbb{R}^{n_{1}}$ solve the local weighted least-squares problems

[TABLE]

which transform point samples into “field moments” and define a mapping

[TABLE]

We chose (28) as the data transfer operator for (20). The second term in (26), i.e., the vector $\bm{\tau}_{\bm{f}}(\bm{\phi})\in\mathbb{R}^{n_{1}}$, contains the integrals of the polynomial basis functions on $\bm{f}$ and defines the face moments for (20). Setting $\bm{\mu}_{ij}=\bm{\tau}_{ij}(\bm{\phi})$ and $\bm{\mu}_{\Gamma_{j}}=\bm{\tau}_{\Gamma_{j}}(\phi)$ in (20) we obtain the following GMLS instance of the MMD operator:

[TABLE]

4.1.1 Analysis

Since $\mu_{i}$ are volumes of physical dual cells on a conforming mesh, assumption T.1 readily holds. We now check the rest of the assumptions of the abstract theory in §3.1. To avoid non-essential technicalities we shall assume that all dual faces are affine.

Lemma 2**.**

The pair $\{T_{F},\{\mu_{i},\bm{\mu}_{ij},\bm{\mu}_{\Gamma_{i}}\}\}$ is $P_{1}$-reproducing.

Proof.

Let $\bm{p}\in(P_{1})^{d}$ and $\bm{p}^{h}\in U^{h}_{\gamma}$ denote its point sample on the cloud. By construction, the GMLS approximation (26) is exact for linear polynomials and so,

[TABLE]

Inserting these identities in (29), using the divergence theorem and the fact that $\nabla\cdot\bm{p}$ is constant, yields

[TABLE]

∎

Lemma 3**.**

The face moments $\{\bm{\mu}_{ij}\}$ and the operator $T_{F}$ satisfy assumptions T.2 and T.3, respectively.

Proof.

To show that T.3 holds note that the kernel $w(\tau_{ij};\lambda_{i})$ is anchored at the centroid $\bm{x}_{\bm{f}_{ij}}$ of $\bm{f}_{ij}$. As a result, (8) implies that $\bm{c}_{ij}(\bm{u}^{h})=\bm{c}_{ji}(\bm{u}^{h})$, i.e., $T_{F}$ is symmetric. For T.2 we use that $\bm{n}_{ij}=-\bm{n}_{ji}$, and so,

[TABLE]

∎

Lemma 4**.**

The pair $\{T_{F},\{\bm{\mu}_{ij},\bm{\mu}_{\Gamma_{i}}\}\}$ is locally Lipschitz continuous on $X$.

Proof.

Let $\bm{u},\bm{v}\in(C^{2}(\Omega))^{d}$ with point samples $\bm{u}^{h},\bm{v}^{h}\in U^{h}_{\gamma}$. Given $\bm{f}\in F\cup\{\Gamma_{j}\,|\,\bm{x}_{j}\in\gamma^{\prime}\}$ we have that

[TABLE]

By assumption, all dual faces are affine and so their unit normals are constant. As a result, $\bm{\phi}_{s}(\bm{x})\cdot\bm{n}$ is linear on $\bm{f}$ and

[TABLE]

It follows that

[TABLE]

The expression in the square brackets is the Moving Least Squares (MLS) approximant $\widetilde{\bm{u}}(\bm{x}_{\bm{f}})$ of $\bm{u}(\bm{x}_{\bm{f}})$; see [14, Chapter 4]. Using the fact that $\widetilde{\bm{u}}(\bm{x}_{\bm{f}})$ is uniformly bounded by the data we obtain

[TABLE]

The proof follows by noting that for physical faces $|\bm{f}|=O(h^{d-1}_{X})$. ∎

Lemmata (2)–(4) verify all requirements of the abstract theory in §3.1. Thus, we have the following result.

Theorem 3**.**

The operator (29) satisfies (13) (local conservation), (14) (global conservation), and is first-order accurate.

4.2 MMD without a background mesh

In this section we drop the assumption of a background mesh and only require that the virtual vertices $V$ are tied to the point cloud $X$, see §3. To implement (20) under these conditions one has to construct the necessary topological, field, and metric data without referencing physical cells and faces, except for the boundary segments $\Gamma_{i}\in\gamma^{\prime}$ where field and face moments can be defined as in §4.1.

4.2.1 Topological data

The virtual boundary in (15) is completely determined by the list of edges connected to the primal vertex $\bm{v}_{i}$. Any graph $G(V,E)$ built on the point cloud can provide this information. In particular, here we shall use the $\varepsilon_{g}$-ball graph $G_{\varepsilon_{g}}(V,E)$ of $X$ with vertices $V:={X}$ and edges

[TABLE]

where $\varepsilon_{g}$ is a given positive real number. This graph may be trivially constructed with $O(p_{\Omega})$ computational complexity using standard binning algorithms. We view $G_{\varepsilon_{g}}(V,E)$ as a surrogate for a primal mesh and associate it with the virtual primal grid $\tau_{h}(V,E)$. The latter induces a virtual dual grid $\tau^{\prime}_{h}(C,F)$ that does not possess any physical entities, except for the boundary segments $\Gamma_{i}$.

Remark 2**.**

We use a (virtual) primal-dual setting to define the abstract MMD operator precisely because any graph on $X$ induces a boundary operator on $\tau^{\prime}_{h}(C,F)$. Building such a graph only requires defining its edges, which can be done efficiently on any point cloud. In contrast, if one were to adopt a single virtual mesh, one would have to construct suitable surrogates for its virtual vertices, edges, faces and cells in order to obtain a correct notion of a virtual boundary. This is a substantially more difficult task.

4.2.2 Field data

In Section 4.1 we used the GMLS coefficients $\bm{c}_{ij}(\bm{u}^{h})\in\mathbb{R}^{n_{1}}$ to define the data transfer operator $T_{F}:U^{h}_{\gamma}\rightarrow\widetilde{U}^{h}_{\gamma}$.

The only physical mesh entities needed to set up the least-squares problem (27) for these coefficients were the dual face centroids $\bm{x}_{\bm{f}_{ij}}$, involved in the kernel $\omega(\tau_{ij},\lambda_{i})$. In the present context the dual faces and their centroids are purely virtual. These virtual centroids can be tied to any reasonable physical location. Here we choose the midpoint $\bm{x}_{ij}=(\bm{x}_{i}+\bm{x}_{j})/{2}$ of the primal edge corresponding to $\bm{f}_{ij}$ as a physical location of the face’s virtual centroid, i.e., we set $\bm{x}_{\bm{f}_{ij}}:=\bm{x}_{ij}$.

Efficient computation and storage of the field moments requires some care though. The number of virtual faces in $\partial\omega_{i}$ equals the valence of the primal vertex $\bm{x}_{i}$. Since in the absence of a background mesh there are no restrictions on the number of edges connected to this vertex, computing and storing the GMLS coefficients $\bm{c}_{ij}(\bm{u}^{h})$ for every virtual dual face may become ineffective. To improve efficiency we compute and store these coefficients on the point cloud. This can be accomplished by changing the kernel from $w(\tau_{\bm{f}};\lambda_{i})=\kappa(|\bm{x}_{\bm{f}}-\bm{x}_{i}|)$ to $\omega(\lambda_{j},\lambda_{i})=\kappa(|\bm{x}_{j}-\bm{x}_{i}|)$, which replaces the matrix $W(\tau_{\bm{f}})$ in (27) by a new matrix $W(\bm{x}_{j})$ with element $W_{ii}(\bm{x}_{j})=\kappa(|\bm{x}_{j}-\bm{x}_{i}|)$. Solution of the modified problem (27) yields the “nodal” GMLS coefficients $\bm{c}_{i}(\bm{u}^{h})$, which are then ”blended” into the field moments according to

[TABLE]

This blending preserves polynomial reproduction property and the blended coefficients (31) satisfy T.3.

4.2.3 Metric data

In Section 4.1 we used the GMLS target vector $\bm{\tau}_{\bm{f}}(\bm{\phi})\in\mathbb{R}^{n_{1}}$ to define the virtual face moments $\bm{\mu}_{\bm{f}}$. In contrast to the GMLS coefficient vector $\bm{c}_{ij}(\bm{u}^{h})$, which only requires the centroids $\bm{x}_{\bm{f}_{ij}}:=\bm{x}_{ij}$, computation of $\bm{\tau}_{ij}(\bm{\phi})$ is impossible without physical dual faces. Likewise, in the absence of physical dual cells we also need a method for assigning a volume $\mu_{i}$ to every virtual dual cell.

In this section we formulate a scalable and computationally efficient algebraic procedure for computing virtual cell volumes and virtual face moments. The procedure starts with construction of virtual volumes $\{\mu_{i}\}$ satisfying T.1. We then seek the components of the face moments $\bm{\mu}_{ij}$ as gradients of scalar potentials. Inserting this ansatz into the $P_{1}$-reproduction condition (21) results in $n_{1}$ independent graph Laplacian problems for the face moments that can be solved in an efficient and scalable way by using algebraic multigrid preconditioners.

Definition of virtual volumes

To construct the virtual volumes we consider the quadratic program

[TABLE]

where $\rho_{i}>0$ are “densities” assigned to each virtual dual cell. In this work we use two different density distributions. The first one is $\rho^{\mbox{\tiny\sf U}}_{i}=1$ in which case solution of (32) yields the uniform virtual volumes

[TABLE]

However, the uniform volumes (33) may not be the best choice for multiresolution problems that call for non-quasiuniform point clouds. For such point clouds one can use the non-uniform density distribution

[TABLE]

where $\kappa(\cdot)$ is the kernel used in the GMLS problem. In this case solution of (32) results in

[TABLE]

Both (33) and (34) satisfy T.1 and numerical results in §6.4 reveal that either one provides first-order convergence on quasi-uniform point clouds. This suggests that the actual volume definition is not critical for the accuracy of the MMD operator on such clouds. Nonetheless, unless otherwise noted, we use (34) because it is simple to implement, while still providing a sense of adaptivity for multiresolution problems.

Definition of virtual face moments

For clarity we discuss construction of the face moments assuming that $\gamma=\Gamma$ and $\gamma^{\prime}=\emptyset$. For the same reason it will be convenient to relabel the polynomial basis functions in (25) and their point samples on the cloud using double indices as

[TABLE]

and $\bm{\phi}^{h}_{k,r}\in U^{h}_{\Gamma}$, respectively. We also relabel the components of the field moments and the virtual face moments to match the indexing of the basis functions, i.e., we write $(\bm{c}_{ij})_{k,r}$ and $(\bm{\mu}_{ij})_{k,r}$ for the components of $\bm{c}_{ij}$ and $\bm{\mu}_{ij}$, respectively, corresponding to $\bm{\phi}_{k,r}$.

Consistency of the MMD operator requires $\{T_{F},\{\mu_{i},\bm{\mu}_{ij}\}\}$ to be a $P_{1}$-reproducing pair. A sufficient condition for this property is that (21) holds for all basis functions, i.e., we require that for all $\bm{x}_{i}\in X_{0}$

[TABLE]

Note that $\nabla\cdot\bm{\phi}_{k,r}=\partial_{k}\phi_{r}$ and $\bm{\phi}_{k,r}\cdot\bm{n}=\bm{n}_{k}\phi_{r}$, while the $P_{1}$-reproduction property of GMLS implies that

[TABLE]

Therefore, (35) is equivalent to the statement that for all $\bm{x}_{i}\in X_{0}$ there holds

[TABLE]

The sum over $\bm{f}_{ij}$ in (36) is the graph divergence of a vector containing the $(k,r)$th components of $\bm{\mu}_{ij}$. This prompts us to seek each component as a scaled graph gradient of a virtual area potential $\bm{\psi}^{h}_{k,r}\in S^{h}$, i.e.,

[TABLE]

where $\bm{x}_{ij}$ is the physical centroid of the virtual face $\bm{f}_{ij}$. Inserting the ansatz (37) into (36) yields the following set of linear equations for the scalar potentials: for $r=1,\ldots,d+1$ solve

[TABLE]

for $k=1,\ldots,d$. Equation (38) is a weighted graph Laplacian problem. Thus, to determine the virtual face moments we need to solve $d+1$ graph Laplacians with different weights, each with $d$ different right-hand-sides. We write these linear systems as

[TABLE]

The standard $n\times n$ graph Laplacian is an M-matrix, and therefore, it may be solved with $O(n)$ complexity using standard algebraic multigrid. For the weighted graph Laplacian $\mathcal{L}_{r}$ we must ensure that the weights $\phi_{r}(\bm{x}_{ij})$ preserve the M-matrix structure. With the assumption that the coordinate system is such that all points in the cloud have positive coordinates, this can be easily accomplished by, e.g., choosing the scalar basis set as $\{\phi_{r}\}_{r=1}^{d}=\{1,x_{1},\ldots,x_{d}\}$.

The operator $\mathcal{L}_{r}$ also has a one-dimensional null space comprising the constant vector $\mathbf{1}$, i.e., $\mathcal{L}_{r}\mathbf{1}=0$. As a result, the $d$ right-hand sides $\bm{b}_{k,r}$, for every one of the $d+1$ graph Laplacian equations in (38) are subject to the compatibility condition $\mathbf{1}^{\intercal}\bm{b}_{k,r}=0$ for $k=1,\ldots,d$. It is easy to see that the volume assumption in T.1 is sufficient for this compatibility condition to hold. Indeed, from the fact that $\partial_{k}\phi_{r}=const$ and $\sum_{i}\mu_{i}=|\Omega|$ it follows that for all $k=1,\ldots,d$ there holds

[TABLE]

4.2.4 MMD without a background mesh: setup and properties

Assuming that $\gamma=\Gamma$ we have the following instance of the abstract MMD operator (20)

[TABLE]

where the field moments $\bm{c}_{ij}$ are given by (31), the virtual volumes $\mu_{i}$ are given by (33) or (34), and the virtual face moments $\bm{\mu}_{ij}$ are defined by solving the graph Laplacian systems (38) for the scalar potentials $\psi^{h}_{k,r}$ and then using the ansatz (37). Thus, although (29) and (39) have the same structure, the metric information in the latter is determined in a purely algebraic way without requiring a physical mesh structure. It is, therefore, of some interest to examine the comparative costs of setting up these operators. To make this comparison more equitable we do not count the mesh generation towards the setup cost of (29), as the purpose of this operator is to improve discretization robustness and accuracy on an already existing mesh.

The setup of (29) involves computing the GMLS coefficients and the face moments for every dual face. Assuming that the valence of every primal node in the mesh is bounded by a constant $C_{\bm{e}}$ the total number of these faces is444Recall that the background mesh nodes are assumed to coincide with the point cloud $X$ and so their number equals $p_{\Omega}$. $O(p_{\Omega})$. The work to compute each $\bm{c}_{ij}$ is proportional to the average number $n_{\kappa}$ of mesh nodes in the support of the kernel function, while the cost to compute every face moment $\bm{\mu}_{ij}$ is proportional to $d+1$. As a result, the setup cost for (29) is $O(p_{\Omega}(n_{\kappa}+d))$.

The setup cost for (39) involves computing the nodal GMLS coefficients $\bm{c}_{i}$ and the scalar potentials $\psi^{h}_{k,r}$ on the point cloud. The work for every $\bm{c}_{i}$ is again proportional to $n_{\kappa}$, while the potentials require the solution of $d+1$ linear systems with $d$ different right hand sides per point. Assuming that these systems are solved using algebraic multigrid the cost will scale as $O(p_{\Omega}d^{2})$, where $p_{\Omega}$ is the number of points in the cloud $X$. Thus, the total setup work for (39) scales as $O(p_{\Omega}(n_{\kappa}+d^{2}))$.

Note that $n_{\kappa}$ can be assumed fixed as long as the polynomial degree of the GMLS reproducing space $\Phi$ remains the same. Therefore, the setup cost for both (29) and (39) scales linearly with the number of points in the cloud $X$.

Remark 3**.**

In the absence of a background mesh the valence of each point in $X$ may grow with the size of the point cloud, leading to an increase of the number of virtual edges in $\partial\omega_{i}$. Using (31) to compute the field moments and the scalar potentials $\psi^{h}_{k,r}$ to compute the virtual face moments allows us to store all field and metric data, required for the evaluation of (39), on the point cloud $X$. This is essential for ensuring that the setup cost and the storage requirements of (39) scale linearly with the size of the point cloud.

Remark 4**.**

Although a detailed discussion of a performant implementation of (39) is beyond the scope of this work we include in Table 1 the number of iterations and total wall-clock time when using preconditioned conjugate gradient method with a classic algebraic multigrid preconditioner [26] to solve (38). The information in the table confirms that these standard techniques work well for the weighted graph Laplacian operators in (38) and achieve linear computational complexity.

Properties

By construction the metric and field data for (39) satisfies the topological requirements T.1-T.3 and $\{T_{F}\{\mu_{i},\bm{\mu}_{ij}\}\}$ is a $P_{1}$-reproducing pair. This is sufficient for Theorem 1 to hold for (39), i.e., this MMD instance satisfies (13) (local conservation) and (14) (global conservation). On the other hand, to show that (39) is also first-order accurate requires $\{T_{F},\{\bm{\mu}_{ij}\}\}$ to be a locally Lipschitz continuous pair on $X$, i.e., that (22) holds. The proof of this property is highly non-trivial and will be addressed in a forthcoming paper. However, preliminary numerical convergence studies in Section 6 suggest that (39) is indeed first-order accurate.

5 A virtual finite volume scheme for a scalar advection-diffusion equation

In this section we apply the abstract MMD operator (20) to construct a meshfree mimetic discretization of the model boundary value problem

[TABLE]

where $u$ is a scalar variable, $f$, $g$ and $h$ are prescribed data, and $\sigma(\cdot)$ is a flux function. To illustrate the application of the MMD operator it suffices to consider a simple advective-diffusive flux given by

[TABLE]

where $\varepsilon>0$ is a picewise $C^{0}$ diffusion coefficient and $\mathbf{a}\in(C^{0}(\Omega))^{d}$ is a given velocity field.

We highlight the ability of our approach to deliver a truly conservative meshfree scheme without a physical background mesh by implementing the second instance (39) of the abstract MMD operator. Thus, we consider the discretization infrastructure of Section 4.2, comprising a quasi-uniform point cloud $X$, an $\epsilon_{g}$-ball graph $G_{\epsilon_{g}}(V,E)$ serving as a surrogate primal mesh $\tau_{h}(V,E)$, and the associated virtual dual mesh $\tau^{\prime}_{h}(C,F)$. The virtual volumes $\mu_{i}$ are set by (34), while the virtual face moments $\bm{\mu}_{ij}$ are defined by (37) with potentials determined by solving the graph Laplacian systems (38). The face moments on the Dirichlet boundary segments are computed using the target functionals in (26), i.e.,

[TABLE]

where $\{\bm{\phi}_{s}\}$ is a basis of $\Phi=(P_{1})^{d}$. To discretize (40) we represent the unknown field $u(\bm{x})$ by its point samples on $X_{0}\cup X_{\gamma}$, the source term $f$ by its point samples on $X$ and assume that the boundary data $g$ and $h$ are given on the Dirichlet and Neumann parts of the boundary, respectively. We thus obtain the following meshfree mimetic discretization of (40): given $f^{h}\in S^{h}$ find $u^{h}\in S^{h}_{\gamma^{\prime}}$ such that

[TABLE]

where $\widetilde{\bm{\sigma}}(u^{h})\in\widetilde{U}^{h}_{\gamma}$ is a “numerical flux”. Note that the Dirichlet boundary condition is enforced in (42) through the definition of the space $S^{h}_{\gamma^{\prime}}$, which implies that $u^{h}_{i}=g(\bm{x}_{i})$ for all $\bm{x}_{i}\in\gamma^{\prime}$. In contrast, the Neumann boundary condition is enforced through the definition of the MMD operator (39) by applying the specified boundary data on all segments $\Gamma_{j}\subset\gamma$. To complete the formulation of (42) it remains to construct a numerical flux $\widetilde{\bm{\sigma}}^{h}$ from the field data $u^{h}\in S^{h}_{\gamma}$. Therefore, our application of the MMD operator to (40) is similar to the second use case described in Remark 1, and that is typical of finite volume schemes. Thus, we shall refer to (42) as a virtual finite volume discretization of the model problem. We discuss construction of the numerical flux in the next section.

5.1 Numerical flux

We will define the numerical flux as a sum

[TABLE]

where $\widetilde{\bm{\sigma}}_{d}(u^{h})\in\widetilde{U}^{h}_{\gamma}$ and $\widetilde{\bm{\sigma}}_{a}(u^{h})\in\widetilde{U}^{h}_{\gamma}$ are diffusive and advective flux moments, respectively. To compute these moments we follow the procedure outlined in Section 4.2.2 and first generate the necessary GMLS coefficients on the point cloud. Then we use (31) to define $\widetilde{\bm{\sigma}}_{d}(u^{h})$ and $\widetilde{\bm{\sigma}}_{a}(u^{h})$ on the virtual faces. In both cases we assume that the diffusion coefficient and the velocity field are represented on the point cloud by their point samples $\varepsilon^{h}\in S^{h}$ and $\mathbf{a}^{h}\in U^{h}_{\gamma}$, respectively.

5.1.1 Advective flux moments

The first step in the construction of $\widetilde{\bm{\sigma}}_{a}(u^{h})$ involves computing a set of “nodal” GMLS coefficients at all points $\bm{x}_{i}\in X_{0}\cup X_{\gamma}$ from a point sample $\bm{\sigma}^{h}_{a}(u)\in U^{h}_{\gamma}$ of the advective flux. We consider two different types of such “nodal” coefficients. The first one, denoted by $\bm{c}^{a}_{i}(u^{h})$, solves the weighted least-squares problem (27) with weight matrix $W(\bm{x}_{i})$ and a point sample555Here $\circ$ denotes the Hadamar (pointwise) product of two vectors. $\bm{\sigma}^{h}_{a}(u)=\mathbf{a}^{h}\circ u^{h}\in U^{h}_{\gamma}$, i.e.,

[TABLE]

The second “nodal” GMLS coefficient, denoted by $\maxvec{\bm{c}}^{\,a}_{i}(u^{h})$, is an upwind version of $\bm{c}^{a}_{i}(u^{h})$ intended for advection-dominated problems. To incorporate upwinding we modify (43) as follows. Given a point $\bm{x}_{i}\in X_{0}\cup X_{\gamma}$ we define the upwind part of the point cloud at this point as

[TABLE]

Thus, $\maxvec{X}_{i}$ is the subset of $X$ containing the point $\bm{x}_{i}$ and all cloud points strictly upwind of that point. Next, we define the upwind point sample of the advective flux as

[TABLE]

Finally, we redefine the basis sample matrix $B$ and the weight matrix $W(\bm{x}_{i})$ in (43) to include only the rows corresponding to the points in $\maxvec{X}_{i}$ and denote the new matrices as $\maxvec{B}_{i}$ and $\maxvec{W}(\bm{x}_{i})$, respectively. The upwind “nodal” GMLS coefficient is then obtained by solving the modified weighted least-squares problem

[TABLE]

We use (31) with $\theta_{ij}=1/2$ and the first type of “nodal” GMLS coefficients to define the standard advective flux moments $\widetilde{\bm{\sigma}}_{a}({u}^{h})$. Setting $\theta_{ij}=1$ if $\mathbf{a}(\bm{x}_{ij})\cdot\left(\bm{x}_{j}-\bm{x}_{i}\right)\geq 0$ and $\theta_{ij}=0$ otherwise, and using the second type of GMLS coefficients defines the upwind moments $\maxvec{\widetilde{\bm{\sigma}}}_{a}({u}^{h})$. Succinctly,

[TABLE]

5.1.2 Diffusive flux moments

As in the advective flux case we obtain $\widetilde{\bm{\sigma}}_{d}(u^{h})$ by blending “nodal” GMLS coefficients. However, generating these coefficients requires a different approach because the point values of $\nabla u$ are not specified on the cloud and so, a point sample of $\bm{\sigma}_{d}(u)$ is not readily available on $X$. This means that we cannot use the setting of (43), which requires such a sample. Instead, we shall obtain the “nodal” coefficients through an auxiliary GMLS problem for the scalar field $u$ that is exact for quadratic polynomials. Specifically, we consider the GMLS framework in Section 2.2 with $U=C^{1}(\Omega)$, $\Phi=P_{2}$, sampling functionals $\lambda_{i}=\delta_{\bm{x}_{i}}$, a target $\tau_{i}(u)=\bm{\sigma}_{d}(u)|_{\bm{x}_{i}}=\varepsilon\nabla u(\bm{x}_{i})$, and kernel $w(\tau_{i},\lambda_{j})=\kappa(|\bm{x}_{i}-\bm{x}_{j}|)$. The GMLS approximation of $\tau_{i}(u)$ at $\bm{x}_{i}\in X$ is then given by

[TABLE]

where $\phi=\{\phi_{q}\}_{i=1}^{n_{2}}$ is basis of $P_{2}$, $\nabla\phi(\bm{x}_{i})=(\nabla\phi_{1}(\bm{x}_{i}),\ldots,\nabla\phi_{n_{2}}(\bm{x}_{i}))$, and

[TABLE]

Since $\nabla{\phi}_{q}\in(P_{1})^{d}$ for all $q=1,\ldots,n_{2}$, it follows that there exists a matrix $G$ such that

[TABLE]

where $\bm{\phi}=\{\bm{\phi}_{k}\}_{k=1}^{n_{1}}$ is basis of $(P_{1})^{d}$. As a result, we can express the target approximation (45) in terms of this basis as

[TABLE]

where $\bm{g}_{i}({u}^{h})=G^{\intercal}\bm{c}_{i}({u}^{h})$ and $\varepsilon_{i}=\varepsilon|_{\bm{x}_{i}}$. We refer to $\bm{g}_{i}({u}^{h})$ as the derived “nodal” GMLS coefficients because they are not based on a point sample of the diffusive flux but rather derived from a point sample of the scalar field $u$. To define the diffusive flux moments on the virtual faces we then use (31) with $\theta_{ij}=1/2$, the derived “nodal” coefficients, and an averaged diffusion coefficient $\varepsilon_{ij}$:

[TABLE]

In Section 6.2, we will explore the use of both the arithmetic mean $\varepsilon^{A}_{ij}=(\varepsilon_{i}+\varepsilon_{j})/2$ and the harmonic mean $\varepsilon^{H}_{ij}=2\varepsilon_{i}\varepsilon_{j}/(\varepsilon_{i}+\varepsilon_{j})$, which is more appropriate when $\varepsilon$ has large jumps.

6 Numerical results

In the following sections, we study numerically the properties of the virtual finite volume scheme (42) for the model problem (40). We consider two versions of this problem corresponding to $\mathbf{a}=\mathbf{0}$ and $\mathbf{a}\neq\mathbf{0}$, respectively. We refer to the first version as the Darcy problem and to the second one as the advection-diffusion problem. For simplicity, we restrict attention to spatial domains in $\mathbb{R}^{2}$.

Unless otherwise noted, we generate the point clouds for our study using the following procedure. Given a desired fill distance $h_{X}>0$ we first construct a uniform partition of the boundary comprising segments $\Gamma_{i}$ of length $h_{X}$. The barycenters of these segments define ${X}_{\Gamma}$. Then, we obtain ${X}_{0}$ by generating a uniform lattice $\mathcal{L}_{h_{X}}$ of spacing $O(h_{X})$ over a square containing $\Omega$, and finally define ${X}_{0}=\mathcal{L}_{h_{X}}\cap\Omega$.

To remove any symmetries that might prompt superconvergence, we then perturb each point in ${X}_{0}$ by a uniformly distributed random variable with magnitude $h_{X}/5$. In so doing we obtain point clouds satisfying the assumption of quasi-uniformity (4). Figure 2 shows a point cloud on the unit disk generated by this procedure. To facilitate post-processing of results, it will often be convenient to generate the Delaunay triangulation corresponding to $X$; it will be made explicit if results are presented on the Delaunay triangulation rather than natively on the point cloud.

6.1 Darcy problem: convergence study

We perform a numerical convergence study of (42) using (40) with $\varepsilon=1$, $\mathbf{a}=\mathbf{0}$, pure Neumann and pure Dirichlet boundary conditions, the smooth manufactured solution $u_{ex}(x,y)=\sin 2\pi x\sin 2\pi y$, and the domain shown the in Figure 3. Substitution of $u_{ex}(x,y)$ into the governing equations defines the necessary data fields for the study. In the pure Neumann case the solution is determined up to an arbitrary constant. To ensure uniqueness we use a Lagrange multiplier to enforce a zero mean condition on the solution. We refer to [27] and the references therein for further information on handling pure Neumann conditions.

To estimate convergence rates we use successively refined point clouds constructed from a sequence of $N\times N$ uniform lattices for $N=16,32,64,128,256$. Rates are estimated with respect to the error measures

[TABLE]

respectively, where $u^{h}_{ex}$ is a point sample of the manufactured solution and $\bm{\sigma}_{d}(u^{h})$ is the GMLS reconstruction of the diffusive flux defined in (45). We refer to these measures as the $\ell_{2}$-norm and the $h_{1}$-seminorm, respectively. Table 2 shows convergence results for the virtual finite volume solution with pure Dirichlet and pure Neumann boundary conditions in these norms. We observe roughly the same convergence rate in both norms, and for both types of boundary conditions, which confirms numerically that (42) is first-order accurate.

6.2 Darcy problem: $H(div)$-compatibility study

In this section we consider Darcy problems with discontinuous coefficients $\varepsilon$. In this case, for any interface $\eta$ corresponding to a discontinuity in $\varepsilon$, solutions of (40) must satisfy the transmission conditions

[TABLE]

where $\bm{n}_{\eta}$ is a unit normal to $\eta$ and $+$ and $-$ denote entities associated with the opposite sides of the interface. Conditions (48) imply that the normal component of the flux $\varepsilon\nabla u$ is always continuous across the interface but its tangential component and the gradient $\nabla u$ may be discontinuous on $\eta$. A hallmark of any $H(div)$-compatible discretization for such Darcy problems is the ability to accurately represent this flux behavior, which is essential for, e.g., accurate simulations of subsurface flow problems.

To evaluate the extent to which the virtual finite volume scheme (42) may be deemed $H(div)$-compatible, we consider three benchmark examples adapted from [28, 29] and driven by pure Neumann conditions. In each one of these examples $\varepsilon$ is a piecewise constant field defined with respect to a disjoint partition $\Omega={\cup}_{i}\Omega_{i}$ of the unit square. The value of $\varepsilon$ on $\Omega_{i}$ is denoted by $\varepsilon_{i}$. Figure 4 shows the partitions for each one of the three examples. We refer to these examples as the two strip, five strip and five spot problems, respectively.

Two strip problem

In this example $\varepsilon$ is discontinuous along the vertical line $x=1/2$ and

[TABLE]

The exact solution satisfies the interface condition (48) and the pure Neumann condition

[TABLE]

where $\bm{n}$ is the outer unit normal to the boundary. Note that the gradient of the exact solution has a discontinuity proportional to the diffusivity ratio $R=\varepsilon_{1}/\varepsilon_{2}$. We use this fact to compare and contrast the arithmetic $\varepsilon^{A}_{ij}$ and geometric $\varepsilon^{M}_{ij}$ mean reconstructions of the diffusivity defined in §5.1.2. Figure 5 compares approximations of $\partial_{x}u_{ex}$ computed by the virtual finite volume scheme using these two reconstructions. The solution obtained with $\varepsilon^{A}_{ij}$ exhibits oscillations in the vicinity of the discontinuity, while using the harmonic mean allows the scheme to accurately represent the discontinuity in the solution. Figure 6 compares solution profiles along $x=0.5$ for increasing diffusivity ratios $R$. We see that the results are comparable for $\varepsilon_{1}/\varepsilon_{2}=O(1)$, but as the ratio is increased the harmonic mean is better able to preserve the location of the jump. Motivated by these results, we use the harmonic mean for all results in the remainder of this paper.

The five-strip problem

In this example $\varepsilon$ has discontinuities along the lines $y=0.2$, $y=0.4$, $y=0.6$ and $y=0.8$, which divide the unit square into five horizontal strips. We set the diffusivity value on each strip to $\varepsilon_{1}=16$, $\varepsilon_{2}=6$, $\varepsilon_{2}=1$, $\varepsilon_{4}=10$, and $\varepsilon_{5}=2$, respectively and let

[TABLE]

This exact solution satisfies the interface conditions (48), the pure Neumann condition

[TABLE]

has a continuous vertical flux component $\varepsilon\partial_{y}u_{ex}=0$, and a piecewise constant horizontal flux component $\varepsilon\partial_{x}u_{ex}=\varepsilon_{i}$ on $\Omega_{i}$. The five strip benchmark is designed to test how well a scheme can represent this flux behavior. Figure 7 shows the approximation of $\varepsilon\partial_{x}u_{ex}$ for both uniform and non-uniform point clouds. The profiles of the horizontal flux approximation along $x=0.5$ for increasing point cloud resolutions are shown in Figures 8. These results show that the virtual finite volume scheme (42) can accurately represent the flux components on both uniform and non-uniform point clouds. In particular, approximations of the discontinuous, horizontal flux component on uniform point clouds are indistinguishable from those obtained by an $H(div)$-conforming, mesh-based scheme using interface-fitted grids. The small wiggles in the non-uniform case are caused by the non-symmetric distribution of points across the interfaces, which is the meshless analogue of an unfitted grid. We note that an $H(div)$ scheme on an unfitted grid would have a similar behavior.

The five-spot problem

This example models an injection-extraction application in which the flow is driven from an injection well in the bottom left corner of the domain to an extraction well in the top right. The diffusivity $\varepsilon$ is a piecewise constant which assigned the same value on the bottom left and top right subdomains and another value on the bottom right and top left subdomains; see the right plot in Fig. 4.

To drive the flow we use the same Neumann conditions as in [29, Figure 26]. Specifically, we apply a flux of $-\frac{1}{8}$ to both faces coincident to the bottom left corner, and $\frac{1}{8}$ to those coincident with the top right (Figure 9). In the case where $\varepsilon_{1}=\varepsilon_{2}$, this problem may be solved analytically via superposition of Greens function solutions, assuming a periodic domain. Figure 10 demonstrates convergence of the virtual finite volume scheme to the analytic solution for both the pressure and the flux.

To highlight the robustness of the virtual finite volume scheme we consider the five spot problem with different diffusivity ratios $R=\varepsilon_{1}/\varepsilon_{2}$. The left plot in Figure 11 compares the pressure computed by (42) with a solution obtained by a mixed method implemented with the RT0-P0 element pair. We see that for all diffusivity ratios the two pressure profiles are nearly identical. Finally, the right plot in the same figure shows that the virtual finite volume scheme continues to deliver stable solutions even for very large material contrasts.

The results in this section suggest that for a range of representative example problems the virtual finite volume scheme does indeed behave in a manner that mimics the behavior of mesh-based $H(div)$-conforming methods, such as mixed finite elements implemented with the RT0-P0 element pair.

6.3 Advection-diffusion problem

In this section we examine the performance of the virtual finite volume scheme (42) for the advection-diffusion problem using a range of Péclet numbers

[TABLE]

spanning both the diffusion-dominated, i.e., $Pe=O(1)$, and the advection dominated, i.e., $Pe\gg O(1)$, regimes. Our objectives are twofold. First, using manufactured solutions, we aim to demonstrate the consistency and the accuracy of the virtual finite volume scheme for both the centered $\widetilde{\bm{\sigma}}_{a}({u}^{h})$ and the upwind $\maxvec{\widetilde{\bm{\sigma}}}_{a}({u}^{h})$ reconstructions of the advective flux defined in (44). Our second goal is to confirm that the construction of $\maxvec{\widetilde{\bm{\sigma}}}_{a}({u}^{h})$ does provide an appropriate notion of upwinding in the advection-dominated regime.

In both cases we consider the unit square with a quasi-uniform point cloud $X$ defined according to the procedure described at the beginning of §6. To study the accuracy of the virtual finite volume scheme we use the same manufactured solution as in §6.1, i.e., $u_{ex}(x,y)=\sin 2\pi x\sin 2\pi y$, and apply Dirichlet boundary conditions at all particles in $X_{\Gamma}$.

Tables 3 and 4 present convergence rates for increasing Péclet numbers using the centered $\widetilde{\bm{\sigma}}_{a}({u}^{h})$ and the upwind $\maxvec{\widetilde{\bm{\sigma}}}_{a}({u}^{h})$ flux reconstructions, respectively. The results in Table 3 show that virtual finite volume scheme with the centered advective flux reconstruction is able to maintain $O(h)$ convergence for small to moderate Péclet numbers. However, as the problem becomes strongly advection-dominated, using the centered flux leads to ill-conditioned discrete equations and failure of the iterative solver to converge. On the other hand, Table 4 reveals that with the upwind reconstruction $\maxvec{\widetilde{\bm{\sigma}}}_{a}({u}^{h})$ the virtual finite volume scheme remains first-order accurate over a wide range of Péclet numbers.

To further study the appropriateness of $\maxvec{\widetilde{\bm{\sigma}}}_{a}({u}^{h})$ in the advection-dominated case we consider an example with discontinuous boundary data and velocity field $\mathbf{a}=(1,2)$ similar to the skew advection test in [30]. To describe the problem setup let $\Gamma_{t}$, $\Gamma_{r}$, $\Gamma_{b}$ and $\Gamma_{l}$ denote the top, right, bottom, and left sides of the unit square domain, respectively. On the bottom and left boundaries we impose Dirichlet boundary conditions given by

[TABLE]

respectively. On the remaining parts of the boundary we impose either the Dirichlet conditions

[TABLE]

or the following outflow condition, derived by setting $\bm{n}\cdot\bm{\sigma}_{d}(u)=0$ in the flux definition.

[TABLE]

Figure 12 shows plots of the solution for both the pure Dirichlet and Dirichet/outflow cases computed by the virtual finite volume scheme with the upwind advective flux for $Pe\in\left\{1,10,100\right\}$. A solution plot for the Dirichet/outflow boundary conditions and $Pe=1000$ is shown on Figure 13. We see that in this case the solution develops moderate crosswind oscillations that remain localized near the internal layer. This behavior is similar to the one observed in, e.g., streamline upwind finite element methods [31] and can be corrected by including an appropriate discontinuity capturing term [32]. However, design of an advective flux reconstruction incorporating such a term is beyond the scope of this paper.

6.4 Sensitivity of MMD with respect to virtual volume definitions

We conclude by examining the impact of the virtual volume definition on the truncation error of the MMD operator (39). For this study we use the smooth vector field $\bm{u}=2\pi\left(\cos 2\pi x\sin 2\pi y,\sin 2\pi x\cos 2\pi y\right)$, the non-trivial geometry illustrated in Figure 3, and quasi-uniform point clouds generated according to the procedure described earlier in §6. Table 5 shows the truncation errors and the corresponding convergence rates for \overset{{\mbox{\tiny\sf GMLS{}^{\star}}}}{DIV} implemented with the uniform volumes $\mu_{i}^{\mbox{\tiny\sf U}}$ and the non-uniform volumes $\mu_{i}^{\mbox{\tiny\sf MR}}$, respectively. The results in this table show that in both cases the MMD operator remains first-order accurate.

7 Conclusion

The principal contribution of this work is the formulation of a consistent and conservative mimetic meshfree divergence operator that can be used with or without a background grid while remaining computationally efficient in the latter case. Specifically, implementation of the operator in the absence of background mesh involves solution of graph Laplacian problems, which can be accomplished in a scalable and efficient manner by using standard multigrid techniques. We demonstrate the accuracy and robustness of our operator by using it to define a virtual finite volume scheme for a scalar advection-diffusion equation. Our numerical results show that the scheme is first-order accurate, exhibits the same traits as an $H(div)$-conforming mesh-based discretization and can be equipped with an appropriate notion of upwinding for advection-dominated problems.

Acknowledgments

This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research under Award Number DE-SC-0000230927, and the Laboratory Directed Research and Development program at Sandia National Laboratories. The first author acknowledges support through the NSF-MSPRF program.