A short FE implementation for a 2d homogeneous Dirichlet problem of a Fractional Laplacian

• Numerical methods for the nonlocal problem remain scarce, with singular kernels and unbounded integration posing major challenges.The cited examples of proposed methods concern one-dimensional settings.Source: §1 · p. 2
• The short MATLAB finite-element code complements the complete n-dimensional analysis in.Source: §1 · p. 2
• Duffy-type transforms split singular element integrals into an explicitly integrable singular part and a smooth part treated numerically, enabling arbitrary precision.Source: §1 · p. 3
• An extended secondary mesh reduces the exterior-integral computation to a case with a circular boundary, including for complex domains.The authors state that this idea extends to arbitrary space dimensions.Source: §1 · p. 3
• The code balances readability and efficiency through vectorized computations, and small modifications may adapt it to other problems.Source: §1 · p. 3
• The code has been successfully used for eigenvalue computations and time-dependent problems, including semi- and full-fractional settings.Source: §1 · p. 4

Fractional spaces

• The section introduces fractional Sobolev spaces and their characterizations, including the distinction at s = 1/2 and the Lions-Magenes space.Source: §2 · pp. 4–5

Weak formulation

• The weak formulation seeks u in the zero-boundary fractional space, with its inner product involving integrals over R^n.Source: §2 · p. 5

Regularity of solutions

• Regardless of the right-hand side regularity, solutions should not be expected to have more than s + 1/2 derivatives in L2(Ω).Near the boundary, u(x)/d(x, ∂Ω)^s remains finite; establishes 1+s−ε derivatives in a suitable weighted space when f belongs to C1−s(Ω).Source: §2 · p. 6

Finite-element discretization

• Continuous piecewise-linear elements on an admissible triangulation define the discrete space, and the resulting method is conforming for 0 < s < 1.The mesh family is assumed shape-regular; the authors choose continuous P1 elements for a unified conforming approach over this range.Source: §3 · p. 6

Linear system

• The discrete problem becomes a linear system with a symmetric positive-definite fractional stiffness matrix and a unique solution.Source: §3 · p. 6

Auxiliary triangulation

• An auxiliary mesh triangulates the region between Ω and a containing ball B, making the complete triangulation admissible.Nodes in the auxiliary domain are treated as Dirichlet nodes; Figure 1 depicts the domain and the two triangulations.Source: §3 · p. 7, Figure 1 · p. 7

Assembly integrals

• Element-pair integrals pose separate challenges: singular integrands for intersecting elements and an unbounded-domain contribution requiring accurate computation of h(x).Because h(x) is radial, it can be computed quickly on the fly or precomputed to arbitrary precision.Source: §3 · pp. 7–8

4. Data structure and auxiliary variables

• The input mesh is supplied through variables for node coordinates, triangle connectivity, boundary and free nodes, auxiliary-element count, and the radius of the containing ball.Triangles in the auxiliary region must be listed last, and the data must be available in the MATLAB workspace before execution.Source: §4 · p. 8
• The supplied data.mat file stores nodes and weights for quadrature rules used throughout the vectorized computations.Its contents are listed in Table 2, with further details in Appendix C.Source: §4 · p. 9, Table 2 · p. 9
• Auxiliary elements extend the original mesh to a ball containing Ω, and nodes outside Ω are treated as Dirichlet nodes.The code initializes the discrete solution, stiffness matrix, and right-hand side before assembling element measures.Source: §4 · pp. 9–10
• The patches structure maps each node to the indices of its neighboring elements, allowing element-pair intersections to be classified efficiently.It is built using a linear number of operations.Source: §4 · p. 10

5. Main loop

• Because disjoint basis supports still interact, the fractional stiffness matrix is full and can require quadratic memory and operations; the code uses vectorized operations and computes only element pairs with ℓ ≤ m.The latter restriction uses stiffness-matrix symmetry to avoid redundant pair computations.Source: §5 · pp. 10–11
• For each domain element, the loop classifies unvisited element pairs as disjoint, vertex-sharing, or edge-sharing, while the right-hand side is assembled using quadrature.The patches structure supports the classification in a linear number of operations; the force is interpolated at element-edge midpoints for the right-hand-side quadrature.Source: §5 · pp. 11–12

5.1. Identical elements.

• For identical elements, triangle_quad and comp_quad numerically compute the two integral contributions, using prebuilt quadrature data.Their outputs are 3 × 3 matrices.Source: §5.1 · p. 12

5.2. Non-touching elements.

• For disjoint triangles, the code computes distances between quadrature points, raises them to −(2 + 2s), and combines them with prebuilt matrices to form local interactions.It assembles the resulting 6 × 6 local matrix for each element pair.Source: §5.2 · pp. 13–14

5.3. Vertex-touching elements.

• For triangles sharing a vertex, vertex_quad computes interactions involving the five nodal basis functions in the pair.The routine uses prebuilt quadrature variables and returns a 5 × 5 array.Source: §5.3 · pp. 14–15

5.4. Edge-touching elements.

• For triangles sharing an edge, edge_quad computes the interaction using four nodal basis functions, with the shared-edge endpoints first in the local ordering.Its output is a 4 × 4 array.Source: §5.4 · p. 15

5.5. Discrete solution.

• After assembly, the code solves the finite-element system; only free-node entries are needed because boundary and auxiliary-domain nodes are set to zero.The displayed solution omits the auxiliary domain.Source: §5.5 · p. 15

6. Numerical Experiments

• The experiments include explicit ball-domain solutions, including k = 0 with constant forcing, and vary s while fixing the auxiliary-ball radius at 1.1.The paper also displays an FE solution for s = 0.75 and k = 2 on a mesh of about 14,000 triangles.Source: §6 · pp. 16, 18, Figure 3 · p. 18
• Computational convergence rates for problem (6.1), measured in the L2(Ω) and energy norms, are reported as being in accordance with theory.The rates are shown with respect to mesh size.Source: §6 · p. 16, Table 3 · p. 17
• Table 4 shows that changing the auxiliary-ball radius does not substantially affect error, while CPU time grows linearly with the number of auxiliary-domain elements.The comparison uses a fixed regular triangulation with 4,228 elements and s = 0.5.Source: §6 · p. 17, Table 4 · p. 17
• For k = 2, the L2(Ω) rate is ≈ (#DOFs)^-3/8 at s = 0.25 and ≈ (#DOFs)^-1/2 at s = 0.75, in agreement with theory.These rates are shown for the Theorem 6.1 solution.Source: §6 · p. 17, Figure 2 · p. 18
• Over 99% of CPU time is spent assembling K, whose expected assembly complexity is quadratic in the number of elements.The authors' tests appear consistent with the expected complexity.Source: §6 · p. 18

A.1. Non-touching elements.

• For non-touching elements, a nonsingular integrand is estimated with 36-point nested Gaussian quadrature, while vectorized basis-product and distance data support efficient block assembly.Quadrature-dependent matrices are precomputed; distances are calculated for element pairs and can be batched across multiple elements.Source: §A.1 · pp. 20, 22, 23

A.2. Vertex-touching elements.

• For vertex-touching elements, a change of variables explicitly computes the singularity, leaving a regular integral estimated using 27 nested quadrature nodes.The resulting quadrature matrices are precomputed and supplied in data.mat.Source: §A.2 · pp. 24–25

A.3. Edge-touching elements.

• For edge-touching elements, the integrals are reduced to the unit cube and estimated with 27 quadrature points, using precomputed matrices for the five transformed terms.The matrices Ψ1 through Ψ5 are stored in data.mat.Source: §A.3 · p. 28

A.4. Identical elements.

• For identical elements, the integration region is decomposed, and changes of variables transform its first two regions to the same domain.The supplied passage begins a further decomposition using a four-dimensional simplex.Source: §A.4 · p. 30
• The domains are transformed into [0, 1]^4 by the respective transformations Th and the Duffy transformation (A.7); integrals over the unit interval are estimated using a 9-point Gaussian quadrature rule.The first three integrals are calculated by hand, while the last has a regular integrand and is estimated by Gaussian quadrature.Source: §A.4 · pp. 31–32

A.5. Complement.

• Because ψ is radial, it suffices to estimate it at x = (x1, 0), where x1 > 0; computing Jℓ uses nested quadrature with 12 points over the reference triangle and 9 over [0, 2π].The computation is expressed as a precomputed matrix multiplied by a vector that depends on the elements under consideration.Source: §A.5 · pp. 33–34

Appendix B. Two auxiliary functions

• The auxiliary functions provide a precondition-based set difference and compute the three-entry right-hand-side vector using edge-midpoint evaluations weighted by |Tℓ|/6.setdiff_ assumes A is consecutive positive integers in ascending order and B satisfies the stated length and maximum constraints.Source: Appendix B · pp. 35–36

Appendix C. Auxiliary data

• The appendix documents auxiliary data stored in data.mat; the displayed code is illustrative and is not part of the program itself.Source: Appendix C · p. 36

C.1. Quadrature points and weights:

• The cube rule uses 27 quadrature points over [0, 1]^3 as input to vertex_quad and edge_quad.Source: §C.1 · p. 36
• For the non-touching-element case, the code uses six points on the reference triangle; comp_quad uses a 12-point rule.Source: §C.1 · p. 37
• The interval rule supplies nine quadrature points and their weights to comp_quad, and its points are also used in triangle_quad.Source: §C.1 · p. 37

C.2. Auxiliary variables to compute non-touching elements case: phiA, phiB and phiD.

• phiA, phiB, and phiD correspond to ΦA, ΦB, and ΦD, respectively, and are set up using six-point quadrature data on the reference triangle.Source: §C.2 · pp. 37–38
• Six local functions encode signed basis-function values on the reference triangle.Source: §C.2 · p. 38
• Nested loops accumulate weighted products into three 18-by-18 auxiliary matrices, whose blocks are reshaped into phiB, phiA, and phiD.Source: §C.2 · pp. 38–39

C.3. Auxiliary variables to compute vertex-touching elements case: vpsi1 and vpsi2.

• vpsi1 and vpsi2 serve as vertex_quad arguments and correspond to Ψ1 and Ψ2 from Appendix A.2.Source: §C.3 · p. 39
• The code weights evaluations of each function at cube quadrature points and stores them in vpsi1 and vpsi2.Source: §C.3 · pp. 39–40
• The setup defines two groups of five functions, psi_D1 and psi_D2, for evaluating the basis-function expressions.Source: §C.3 · pp. 39–40

C.4. Auxiliary variables to compute edge-touching elements case: epsi1,

• The epsi1–epsi5 arrays supply edge_quad with the Ψ1–Ψ5 quantities defined in Appendix A.3.Their entries are formed by evaluating products of the corresponding ψ_D functions at p_cube quadrature points and multiplying by w_cube.Source: §C.4 · pp. 40–41
• Each epsi array is initialized with 16 rows and 27 quadrature-point columns.The code fills each row indexed by an ordered pair of local indices, using sub2ind([4 4], i, j).Source: §C.4 · p. 41
• The five integrands use products of paired ψ_D functions, with x^2 and, for epsi2–epsi5, an additional y factor.These integrands are evaluated at the points in p_cube before their values are weighted by w_cube.Source: §C.4 · p. 41

C.5. Auxiliary variables to compute identical elements case: tpsi1, tpsi2 and tpsi3.

• The tpsi1, tpsi2, and tpsi3 matrices serve as the Ψ1–Ψ3 inputs to triangle_quad, as defined in Appendix A.4.The code evaluates products of paired λ_D functions at p_I and multiplies by the quadrature weights w_I.Source: §C.5 · pp. 41–42
• The three λ_D function families provide the component functions used to form the respective products for tpsi1–tpsi3.Each product pairs the same family’s functions indexed by i and j.Source: §C.5 · p. 42

C.6. Auxiliary variable to compute quadrature over complement: cphi.

• The cphi matrix provides the Φ input to comp_quad, using products of local functions evaluated at p_T_12 and weighted by w_T_12.The local functions are 1−x, x−y, and y.Source: §C.6 · pp. 42–43

Appendix D. Main Code

• The main program computes triangle areas, builds vertex patches, and preallocates arrays before its main loop over elements.Within the loop it forms the right-hand side and adds triangle, complement, vertex-touching, edge-touching, and disjoint-element contributions to K.Source: Appendix D · pp. 43–45
• After assembly, the program solves the system on the free nodes, scales the result by cns, and plots uh over the mesh.The solve is given by uh(nf) = (K(nf,nf)\b(nf))./cns.Source: Appendix D · p. 45