An unfitted radial basis function generated finite difference method applied to thoracic diaphragm simulations

The thoracic diaphragm is the muscle that drives the respiratory cycle of a human being. Using a system of partial differential equations (PDEs) that models linear elasticity we compute displacements and stresses in a two-dimensional cross section of the diaphragm in its contracted state. The boundary data consists of a mix of displacement and traction conditions. If these are imposed as they are, and the conditions are not compatible, this leads to reduced smoothness of the solution. Therefore, the boundary data is first smoothed using the least-squares radial basis function generated finite difference (RBF-FD) framework. Then the boundary conditions are reformulated as a Robin boundary condition with smooth coefficients. The same framework is also used to approximate the boundary curve of the diaphragm cross section based on data obtained from a slice of a computed tomography (CT) scan. To solve the PDE we employ the unfitted least-squares RBF-FD method. This makes it easier to handle the geometry of the diaphragm, which is thin and non-convex. We show numerically that our solution converges with high-order towards a finite element solution evaluated on a fine grid. Through this simplified numerical model we also gain an insight into the challenges associated with the diaphragm geometry and the boundary conditions before approaching a more complex three-dimensional model.


Introduction
During the 2020 covid-19 pandemic we have all been made aware that intensive care units (ICU) have limited capacity with respect to the number of patients that can be cared for at the same time.The WHO report on covid-19 in China [1] indicates that covid-19 patients with severe symptoms need 3-6 weeks in ICU.Patients with severe respiratory symptoms are put under mechanical ventilation to save their lives.At the same time, the mechanical ventilation has adverse effects, such as ventilator induced diaphragmatic dysfunction (VIDD) [2], which prolongs the ICU time.Hence, improving mechanical ventilation can have a significant impact on ICU patient turnover.
This work is part of the INVIVE project 1 , where we aim to create a mechanically ventilated virtual patient [3] on whom we can perform tests with different ventilation strategies and counter measures against VIDD [4].With the virtual patient, we can also vary the gender, physiology, age and potential injuries that affect the response.Using simulations allows to create a controlled computer-based environment in a way that is not possible in a clinical setting.
The diaphragm is the main respiratory muscle, and the present focus of our study.The diaphragm is located between the thoracic and abdominal cavities.It has two domes, and is attached to the lower ribs, the spine, and the sternum.During mechanical ventilation, the normal action of the diaphragm, where inhalation follows the contraction of the muscle, is reversed, i.e., the muscles become passive and the air is pumped into the lungs by the ventilator.As air is entering the lungs with a positive pressure, the muscle fibres are instead extended.This sudden and extreme mechanical perturbation is the trigger of a chain of biological events causing the progression of VIDD.
Numerical simulation of the biomechanical action of the diaphragm during respiration or ventilation is a challenging problem.The shape of the diaphragm is non-trivial and there are gradual transitions between muscle and tendon with very different material response.Imperfect data can be extracted from medical images, and can then be converted into a geometry representation [5].The specific challenges of constructing a smooth geometry representation are addressed in a forthcoming paper [6].
In this paper, we model the diaphragm using a linear elastic PDE.This is a simplification and we plan to develop a more realistic tissue model as a part of our future work.
The boundary conditions for the elastic PDE system are given by a combination of traction boundary conditions, which contains first derivatives of the displacement, and (timedependent) Dirichlet boundary conditions for the displacement, where the diaphragm is attached.These are mixed boundary conditions, in general not fully compatible, leading to reduced regularity of the solution even when we expect a smooth solution from a physiological perspecitve.Therefore, we smooth the boundary data as well as the transition between the traction boundary data and Dirichlet boundary data before using it in the PDE solver.
We use the unfitted radial basis function generated finite difference method in the leastsquares setting (unfitted RBF-FD-LS) [7] to solve the diaphragm problem.A benefit of using the unfitted RBF-FD-LS method is that the PDE problem is solved on an extended domain (see Figure 3).This simplifies node generation as the node placement is then independent of the geometry.Another benefit of the unfitted setting is a smaller approximation error near the boundaries [7], where the stencils are typically highly skewed when using conventional RBF-FD methods such as the fitted RBF-FD-LS method [8] and the collocation RBF-FD method [9].
As model geometry we use a two-dimensional cross section of the diaphragm.To investigate the properties of the problem we use a combination of data from medical images and knowledge about the physiology of the diaphragm expressed in terms of boundary conditions.In particular, we construct two benchmark problems: (i) a pure Dirichlet case, (ii) a case with mixed boundary conditions (Dirichlet + traction), rephrased as a Robin condition with smooth coefficients.With the experiments performed for these benchmarks, we aim to answer the following questions: • Can we solve the benchmarks problems with unfitted RBF-FD-LS?Are there specific numerical challenges to be noted?
• Can we achieve high-order convergence to the solution of the elastic PDE when the imposed boundary data and the geometry is smooth?
• How is the performance and accuracy of the RBF-FD-LS solver affected by the type of boundary conditions that are imposed?
• How does the RBF-FD-LS solver compare with a basic FEM solver?Do the solvers give similar solutions?
Other authors have also modeled the diaphragm numerically.The most advanced diaphragm simulations in the literature can be found in a series of publications by Ladjal et al. [10,11,12,13,14].The application focus is to track the motion of lung tumours during respiration for radiotherapy purposes.The finite element models that are employed are highly elaborate, taking into account features such as patient-specific lung compliance in order to define the constitutive law.The simulation times reported are far from real-time capability.In these studies, the diaphragm is not the main target and it is not well-resolved in the thickness direction due to the high aspect ratio.Other relevant, but less detailed, diaphragm simulations can be found in [15], where the diaphragm is discretized using shell elements to compare healthy and pathological situations, and in [16], where the whole region under the lungs, including the diaphragm, forms one region in the simulation.
The outline of the paper is as follows: Section 2 describes the expected behavior of the diaphragm during respiration and the relation of the two-dimensional geometry to the real three-dimensional diaphragm geometry.Then the linear elasticity model problem is defined in Section 3.An RBF-FD algorithm for computing differentiation and evaluation matrices is described in Section 4. In Section 5 these matrices are then used in the least-squares unfitted RBF-FD setting to discretize the linear elasticity equations.Section 6 discusses how the boundary of the two-dimensional diaphragm cross-section, and the boundary conditions, are smoothed.In Section 7 and Section 8 we compute the solution to the linear elasticity equations for the two benchmark problems, and evaluate the convergence numerically.Finally, Section 9 contains the conclusions.

The expected behavior of the thoracic diaphragm
The diaphragm is a musculotendinous structure, approximately double-dome shaped, separating the thoracic and abdominal cavities.It is the main muscle of the physiological respiration, performing 70-80% of the work of breathing, although it also has non-ventilatory functions, e.g., coughing, hiccups, sneezing, vomiting, and postural functions.It is composed of three main parts, see Figure 1a.There are two muscular parts, one median and horizontal in upright position, separating the thoracic and the abdominal organs, and another lateral muscle part, which ends with the costal insertions (attachment to the lower ribs), and there is one central tendon, where the extremities of all the muscle fibers converge.When the diaphragm contracts during inspiration (inhalation) with a piston-like motion, the muscle zone thickens (inspiratory thickening), and the domes move caudally (downward) expanding the thorax.Therefore, the air enters the lungs under a pressure gradient, see Figure 1b.The two-dimensional geometry used for our simulations was extracted from a real patient medical image.We used a 3D CT scan image (resolution: 0.927 × 0.927 × 0.3 mm 3 ) that was acquired for medical reasons, see our previous work [5].The diaphragm was manually segmented in 3D following visual cues as explained in [17].The labeled voxels were exported to a triangle mesh with the Marching Cube algorithm and the mesh was simplified with a decimation algorithm [18].
The frontal plane slice we selected is in the middle of the body, corresponding roughly to the anatomic illustration of Figure 1.The raw data consists of a list of 2D vertices where a topology can easily be extracted.The raw geometry data contains noise from several sources.There is some CT scan device incertitude (noise, calibration etc) [19], the data is the result from a discretized process, and the labeling data comes from a manual segmentation that is prone to human error.Concerning this last point, the diaphragm is not entirely visible, sometimes part of its voxels also include other organs.The accuracy is linked to the imagination and the anatomical knowledge of the medical expert that performed the segmentation.
The expected displacement of the diaphragm between the relaxed and contracted states in the two-dimensional slice is is roughly illustrated in Figure 2.For the numerical simulations in Sections 7 and 8, we construct boundary conditions to replicate this motion qualitatively.

Equations of linear elasticity
The simplified model that we use is valid for studying small deformations of an isotropic and homogenuous diaphragm with a linear elasticity constitutive law.The deformation of a diaphragm Ω is described by applying a displacement field u = (u 1 , u 2 ) T , to an object location y = (y 1 , y 2 ) T ∈ Ω such that: where y * = (y * 1 , y * 2 ) T ∈ Ω * is then a deformed object.A field derived from the displacements is the stress tensor: which measures the internal forces in Ω as a consequence of deformation.It is related to the displacement field by: where ε ∈ R 2×2 is the strain tensor and Tr(ε) = (ε 11 + ε 22 ) I is its trace.The scalars λ and µ are the Lamé parameters, which are related to the Young modulus E and the Poisson ratio ν of the material through: For our computations we use E = 10 5 Pa and ν = 0.3.A special stress measure that is also of interest to us is the Von Mises stress: which provides a scalar measure of the total stress.The equations of linear elasticity are derived from the force equilibrium (Newton's second law) imposed on Ω.We have: where T is a field of internal forces in the horizontal and the vertical direction.The Dirichlet and traction boundary conditions are prescribed on two disjoint parts of the boundary, which together form the whole boundary ∂Ω = ∂Ω 0 ∪ ∂Ω 1 : The first boundary condition with the right hand side g = (g 1 (y), g 2 (y)) T is the displacement condition.The second boundary condition with the right hand side h = (h 1 (y), h 2 (y)) T is the traction condition.An equivalent form of ( 5) is the Robin boundary condition: where κ 0 and κ 1 correspond to two spatially dependent coefficients, in this case discontinuous: κ 0 (y) = 1 when y ∈ ∂Ω 0 and zero otherwise, and κ 1 (y) = 1 when y ∈ ∂Ω 1 and zero otherwise.
In Section 7 we first compute solutions using a Dirichlet condition for the whole boundary.(κ 0 (y) = 1, κ 1 (y) = 0, y ∈ ∂Ω).Then, in Section 8, we compute solutions using smoothed Dirichlet and traction conditions (smooth Robin coefficients), which approximates the imposition of these two conditions on two disjoint parts of the diaphragm.However, since physiology suggests a smooth solution, we see this as modeling rather than as an error.
For all of the computations in this paper we use the displacement formulation of ( 4), obtained by using the relation between stress and displacement given in (1).The force equilibrium (4) then expands to: and the traction boundary condition from ( 5) to: The displacement formulation of the linear elasticity equations with the boundary condition from ( 6) is then written as: For simplicity we rewrite the system above as: where: Here D 2 , D 1 , D 0 are the expanded operators that in (8) correspond to Ω, ∂Ω 1 and ∂Ω 0 respectively.If we let ∇ ij = ∂ 2 ∂yi∂yj and ∇ i = ∂ ∂yi , then the operator D 2 is: and the boundary operators D 1 and D 0 are:

Approximation of functions using the RBF-FD method
Here we describe the approximation framework which is used when representing the unknown displacement and the stress (which are the outputs of the simulation), and the known data vectors (the input data to the simulation) as continuous functions.The following discussion is general in the sense that the data vector can be known or unknown.Concrete examples of using the framework developed in this section are given in later sections.
Given a data vector ũ(X) = [ũ(x 1 ), ũ(x 2 ), ..., ũ(x N )] T where x i ∈ R d we construct a semi-discrete evaluation operator E = E(y, X), which, for any point y ∈ Ω ⊂ R d , returns the value of the interpolant of the data ũ(X) at that point, such that: Now we use (11) for each y ∈ Y to obtain a system of equations and form the matrix Note that in our notation this is equivalent to setting y = Y .The RBF-FD procedure constructs the matrix E(Y, X) using a sequence of local, stencil-based interpolation problems which are exact for a cubic or quintic polyharmonic spline basis (PHS) and a (multivariate) monomial basis of degree p.The matrix E(Y, X) is sparse, that is, it contains n N non-zero elements per row, where n is the stencil size defined by: Details about using PHS plus the monomial basis in RBF-FD approximations can be found in [20,21,22,23].In the same way, we construct a semi-discrete oprator for differentiation D L = D L (y, X) which locally evaluates any derivative L of ũ(Y ): Both matrices, E(Y, X) and D L (Y, X), can be formed using the MATLAB code available in [24].For completeness, we below provide the steps to compute the matrix elements (weights) and to assemble D L (Y, X) and E(Y, X) for the point sets Y and X in any dimension d. 3. Form a square interpolation matrix A (k) , where

Let x
j , i, j = 1, .., n belong to the stencil.4. Form a rectangular polynomial matrix P (k) , where ., m is a sampled d-dimensional monomial basis with m basis functions.When n is chosen as in ( 12) then m = n 2 and the size of the matrix P (k) is n × n 2 .5. Using A (k) and P (k) , form the augmented local interpolation matrix: 6. Repeat steps 1-5 for every x ∈ X in order to form all Ã(k) .
The local evaluation and differentiation weights can then be computed in the following way: 1. Take one evaluation point y l ∈ Y and find the first closest point from the X point set.We denote it by x When using these steps to compute the weights, it is possible to avoid storing the matrix Ã(k) by combining the two parts such that the evaluation/differentiation weights are computed for all y l which select x (k) 1 as its closest stencil center, immediately after step 5 of the first list.Once the matrix of local weights W L (Y ) of size M × n is computed, these weights are assembled into a global rectangular matrix, E(Y, X) for evaluation and D L (Y, X) for differentiation, both of size M × N .This can be done by using a matrix Γ ∈ Z N ×n , which is a list of indices of the local neighborhoods of points around every stencil point Additionally, a matrix κ ∈ Z M ×1 is needed, which is a list of indices of the stencil centers that are closest to each evaluation point y l .Γ and κ are found using the k-nearest neighbor method with k equal to n and 1, respectively.Using the MATLAB programming language, the two lists and the sparse matrix can be efficiently computed by invoking the following three commands: where the first two arguments to the function sparse() have the same shape as W L (Y ) and contain the row and column indices for inserting the locally computed weights into the global matrix.

Discretization of the linear elastic model using the uniftted RBF-FD method
In this section we employ the unfitted RBF-FD method [7] to discretize the linear elasticity equations ( 8), ( 9) and ( 10) over the diaphragm.The method relies on constructing rectangular differentiation matrices D L (Y, X) as described in Section 4. The matrices then replace the differential operators in the PDE, together with the corresponding discrete right hand sides F (Y ).
To form the differentiation matrices, we first decide the degree p of the monomial basis that we are appending to the PHS approximation, and compute the stencil size n using (12).Then we construct an interpolation point set X that extends over the diaphragm (see Figure 3), and an evaluation point set Y (see Figure 4) that conforms to the geometry of the diaphragm.Those two point sets are obtained in four simple steps: • The initial point set X 1 is a tilted Cartesian node layout with spacing h in a box that encloses the diaphragm.(the grey and blue points in Figure 3).• Then the point set Y 1 is generated by placing q points with average spacing h y in each Voronoi region inside the diaphragm, defined by the points in X 1 .(the red points in Figure 4).
• In addition, the point set Y b with the same average spacing h y is generated by placing points along the boundary of the diaphragm, see Figure 4.
• The final evaluation point set is given by Y = Y 1 ∪ Y b , see Figure 4.
• The final node set X (the blue points in Figure 3 and Figure 4) is formed by reducing X 1 by removing points that fall more than half stencil size outside the diaphragm.In Matlab this can be done by: X = X_1(unique(knnsearch(X_1,Y,'k',ceil(0.5*n)),:); The last step ensures that the columns of E(Y, X) and D L(Y,X) are non-zero [7].Note that the evaluation point set does not need to be constructed by placing precisely q points in every Voronoi region.It is possible to use any global point set, which is quasi-uniform by nature (e.g.Halton points), and as such on average samples every Voronoi region with approximately q points.Next, we use the method described in Section 4 to discretize the continuous operator D 2 in (9) using the interior evaluation points Y 1 , and to discretize the continuous operators D 1 and D 0 in (10) using the boundary evaluation points Y b .We let D ij k ( • , X) denote the differentiation matrix that approximates element i, j, of the operator D k , k = 0, 1, 2. If we let i, j = 1, 2, while i = j, then we can express the differentiation matrices as: We express the discrete Robin coefficients as In the numerical experiments, we use the following scale factors: where h and h y are the average internodal distances in the X and Y point sets, respectively.These are computed as: The choice of scale factors is based on the papers [7,8], where the scale factors h y for the interior and h 1/2 y for the boundary are used, such that the norms of the discrete least squares problem approximate the continuous L 2 -norm.The additional 1/h scaling increases the weight of the boundary conditions and improves convergence.The Lamé parameter µ is large and affects the scaling between different equations.Therefore, we also include the factor 10/µ in β 1 .
We solve the least-squares problem (15), for the nodal values ũ1 (X) and ũ2 (X), using backslash in MATLAB.After that the solution is evaluated at the Y point set, using the evaluation matrix E: The strains and the stresses (1), ( 3) are computed by applying appropriate differentiation matrices to the solution coefficients ũ1 (X) and ũ2 (X).

Smoothing of geometry and boundary data
If we view the two-dimensional diaphragm from the continuous perspective, the geometry can be described as a closed curve.We choose to parametrize this curve by t ∈ [0, 2π], where the starting point t = 0 is the same as the final point t = 2π.To benefit from the potential high-order convergence of the unfitted RBF-FD method, we need to approximate the boundary curve and the boundary data as functions of t with enough smoothness that the convergence of the PDE problem is not adversely affected.To avoid reduced accuracy due to boundary errors near the artificial end points of the interval, we extend the domain periodically to t ∈ [−2π, 4π] for the approximation.We discretize the extended domain using the uniformly spaced node points T = {t j } Ng j=1 .Given Mg data points, we replicate these periodically to get the extended data set , where M g = 3 Mg > N g .We use a one-dimensional RBF-FD approximation, based on a quintic PHS basis augmented with polynomials of degree p g = 6, to form an overdetermined linear system for the nodal values g(T ) To enforce continuity at t = 0, we add the following equality constraints and solve the constrained least squares problem where B g contains the p g constraints ( 20) and λ contains the corresponding Lagrange multipliers.
To find the smooth boundary curve from the initial vertex data xd i , i = 1, . . ., Mg , we first scale the data such that x d i = (p i , q i ) = s Ω xd i , i = 1, . . ., Mg , where s Ω = 156.92−1 mm −1 .The scaling was chosen such that all data points fall within [−1, 1] 2 .Then we compute an approximate arclength parametrization using the Euclidean distance between the scaled vertices, such that Then we replicate the data over the extended domain.Finally, system ( 21) is solved for each coordinate function p(t) and q(t).The resulting boundary curve is shown in Figure 5 and the individual coordinate functions are shown in Figure 6.Since the curve parametrization here is in the clockwise direction, the outward normals are computed as n(t) = (−q (t), p (t))/ (−q (t), p (t)) .For the boundary data functions g and h in (7), we manufacture data to mimic the expected physiological behaviour shown in Figure 2. A few data points are placed in the regions where we have some information (regions 1, 3, 5, 7, and 9 in Figure 2), and then the rest of the data is generated through linear interpolation.This results in gradual transitions in the regions where we lack information.The data points and the resulting curves are shown in Figures 7 and 15.

Benchmark I: Deformation of the diaphragm using the smoothed Dirichlet boundary conditions
In this section we solve the discretized linear elasticity equations (15).We are interested in whether the problem with smoothed Dirichlet boundary condition leads to a high-order convergence.A point of interest is also whether the resulting deformation is physiologically sensible, and whether the Von Mises stress is distributed as expected.
The boundary condition that we use is purely Dirichlet, which means that we set the Robin coefficients in (8) to κ 0 (y) = 1 and κ 1 (y) = 0. Then the only boundary data functions present in the system (15) are g 1 (y) and g 2 (y), which correspond to the imposition of displacements in the horizontal and vertical direction, respectively.

The imposition of boundary displacements
The displacements have been synthesized to reproduce the physiological behavior described in Section 2. Particularly, the translation of the horizontal part of the diaphragm (region 5) has been defined as a constant vertical displacement in the downward direction and the thickening of the appositional zone (regions 3 and 5) is also prescribed as a constant.As described in Section 6, we place a few data points based on this information, and then the rest of the data is generated by linear interpolation.From a physiological perspective all displacements should be smooth.We generate a smooth function by solving the constrained least squares problem (21).The results are shown in Figure 7. Imposing smooth boundary data makes it possible to obtain high-order convergence and provides a solution that is physically relevant.In Figure 8   The dashed lines show the location of the end points of the diaphragm (regions 2 and 8 in Figure 2).
and after application of the displacements from Figure 7.

Solution of Benchmark I
The solution is given in Figure 9, where we display the spatial distribution of the displacements and the Von Mises stress.Looking at the ũ1 distribution in the figure, we can observe thickening of the diaphragm in regions 2, 3, 7 and 8 according to the labels from Figure 2. Next, the distribution of ũ2 indicates that translation is largest in regions 4, 5 and 6.The Von Mises stress is largest at the interfaces between regions 2, 3 and regions 7, 8.This makes sense, since the change in the thickness is largest in these regions.We conclude that the behavior roughly follows the physiological cues described in Section 2.

Displacement ũ2
Von Mises stress The computed displacements and the corresponding von Mises stress over the diaphragm.This solution was obtained using the unfitted RBF-FD discretization with internodal distance h = 0.004, oversampling parameter q = 5, and an appended polynomial basis of degree p = 5.Due to the scaling applied to the geometry (see Section 6), the displayed results for displacement and stress should be multiplied with s −1 Ω = 0.15692 m to recover the results corresponding to the unscaled problem.

Convergence under node refinement
While the solution is roughly what we expect from a physiological perspective, we are yet to understand whether the simulation gives a correct answer from a numerical perspective.We investigate the convergence of the numerical solution under node refinement using several different polynomial degrees in the stencil-based approximation.Since we do not know what the true solution is, we measure convergence of the numerical solution ũ(Y ) towards a numerical reference solution ũ * (Y * ), where the node set is highly refined.We choose two numerical references: (i) computed using the unfitted RBF-FD-LS method with internodal distance h = 0.002, leading to N = 43 841 and polynomial degree p = 5, (ii) computed using the Galerkin finite element method with linear elements and 236 414 degrees of freedom (corresponding to h = 0.00085), using the GetDP solver [25].Every numerical solution ũ that we obtain is interpolated, consistent with the approximation order, to the point set of the reference solution Y * , where we then compute the approximation error: Table 1 shows the relation between the internodal distance h, computed according to (17),

Displacement ũ1 Displacement ũ2
Von Mises 100 200 300 p=2 p=3 p=4 p=5 Figure 10: Benchmark I: Convergence of the displacements ũ1 and ũ2 and the Von Mises stress for different polynomial degrees p, against a highly resolved numerical solution computed using the unfitted RBF-FD-LS method, with h = 0.002 and p = 5.
and the number of degrees of freedom N for the considered problem sizes.The results when the unfitted RBF-FD method is used as a reference are displayed in Figure 10.We observe that the error is small for all polynomial degrees.The convergence rate increases as p is increased.For all p the convergence rates of the displacements are close to p − 1.The convergence rate of the Von Mises stress is larger than p − 2 for every p.Stress is computed using the first derivatives of the displacements, which (in theory) lowers the convergence rate with one order.
The results when the finite element method is used as a reference are given in Figure 11.The convergence plots for the displacement are very similar to the results in Figure 10.The convergence rate of the stress for p = 5 is lower compared with the self-reference test provided in Figure 10.Furthermore, when p = 5, the convergence seems to be stalling at the last point of observation.Our speculative reasoning is that the derivatives in the finite element space do not approximate the stresses well enough there.This implies that the FEM mesh is too coarse to match the accuracy of the unfitted RBF-FD method for derivative approximation, when the polynomial degree is as high as p = 5, for the range of h used here.

Displacement ũ1
Displacement ũ2 Von Mises The spatial distribution of the error when using the self-reference solution and the finite element reference solution is shown in Figure 12 and Figure 13, respectively.The numerical solution was computed using h = 0.004, q = 5, p = 5.For all solution fields and both references we can observe that the error is larger in the regions with larger Von Mises stress.In addition, the error also tends to be larger in the regions where the boundary curve is concave.This implies that our problem could benefit from adaptive node refinement, which we are planning to use in our future work.

Benchmark II: Deformation of the diaphragm using the smoothed Robin boundary conditions
This benchmark includes a more difficult problem compared with Benchmark I from Section 7. As discussed in Section 2, we can measure the displacements over certain parts of the boundary, i.e., where the diaphragm is fixed near the spine, or where it moves together with Displacement ũ1

Displacement ũ2
Von Mises stress

Displacement ũ1
Displacement ũ2 Von Mises stress the sternum and some ribs.In the regions where the displacements are not known, we may instead have information about the thoracic or abdominal pressure.A pressure condition is a special case of a traction boundary condition.A straightforward way to handle these boundary conditions would be to impose the Dirichlet condition (the known displacements) and the traction condition (the known traction values) in disjoint regions.This implies a discontinuous imposition of boundary conditions, for which: (i) we can not obtain highorder convergence according to our preliminary tests, (ii) the resulting deformation might be discontinuous or have large local derivatives, which would not reflect the physicological behavior of the diaphragm.For this reason we introduce a smooth blending of the Dirichlet and the traction boundary conditions in Robin form (6).This setting implies that we solve the system of equations (15), where the supports of the Robin coefficients κ 0 (y) and κ 1 (y) are chosen to overlap slightly in the regions where the type of boundary data changes.In Figure 14 we display the Robin coefficients which we use in Benchmark II.The displayed coefficients are a function of the boundary parameter t.The coefficients were computed using a sum of the sigmoid functions 1 1+e ε(t−d i ) , where ε = ±20, t is the boundary parameter and d i , i = 1, 2, .. are the transition points (marked by dashed lines in Figure 14).The sign of ε depends on whether the coefficient is increasing (positive sign) or decreasing (negative sign).2).The parameter t corresponding to the boundary of the diaphragm is illustrated in Figure 5.
The smoothed boundary displacement values g 1 , g 2 and the boundary traction values h 1 , h 2 are given in Figure 15.Through the functions g 1 and g 2 we impose thickening in the regions 2, 3, 7 and 8 based on the distribution in Figure 2. Additionally, we also impose a gentle translation in the negative vertical direction in the regions 7 and 8. Through the functions h 1 and h 2 we impose the translation of zone 5 in the negative vertical direction.

Solution of Benchmark II
The solution is given in Figure 16, where we can see the spatial distribution of displacements and the Von Mises stress.To obtain this figure we used h = 0.006, p = 5, and q = 5.We observe that there is a slight thickening in regions 2, 3, 7 and 8, according to the labels from Figure 2, which corresponds to the imposed thickening up to some extent.In regions 4 and 6, we observe a slight bend towards the interior, and in region 5, a translation in the negative vertical direction.This behavior does not entirely mimic the physiological contraction of the diaphragm, however, constructing a more accurate model is planned as future work.This solution serves as a test to understand whether the unfitted RBF-FD method can handle the problem with the smoothed Robin boundary condition, as well as to explore how this type of boundary condition affects the behaviour of the solution.

Convergence under node refinement
We validate the numerical solution by studying convergence under node refinement.Our choice of reference solution for this study is a fine solution computed using the unfitted RBF-FD method with an internodal distance h = 0.002 and a polynomial degree p = 5 used for constructing the local approximations.The relative error against the reference solution is computed using (22).The relation between the considered h and the number of points N in X is given in Table 1.The errors for different choices of the polynomial degree p are displayed in Figure 17.We do not observe convergence for any p when h is too large.The reason for this is that the problem is not resolved yet.When h is sufficiently small, and when p = 4 or p = 5, the solution converges at least with order p − 1, which is desired.When p = 2, we do not see convergence, since we would need an even higher resolution for this (small) polynomial degree.When p = 3 we observe an approximately first order convergence for a sufficiently small h.This order is expected to increase to p − 1 = 2 if h is refined further.Low-order or no convergence when p is small advocates using a higher order method.
The spatial distribution of the error computed using h = 0.006 and p = 5 is given in Figure 18.Here the error is largest close to region 5 from Figure 2, where we have enforced the traction boundary condition, which includes derivatives.In our experience the error is normally larger in the regions where derivative boundary conditions are imposed compared with the regions where a Dirichlet condition is imposed, and we also observed this in [7,8].

Conclusions
The unfitted RBF-FD method in the least-squares setting [7] provided a robust framework for solving the linear elasticity system of PDEs over the simplified diaphragm geometry.
We could also use the unfitted RBF-FD method to smooth the geometry curve and the boundary data.Ensuring that all components of the model were smooth allowed us to

Displacement ũ2
Von Mises stress achieve high-order convergence, which reduces the number of unknowns needed for a given approximation error.In the experiments (not unexpectedly), Benchmark II, with Robin boundary conditions, proved more challenging to solve than Benchmark I with Dirichlet conditions.We needed higher resolution to achieve the same error level, and we did not see any convergence for the larger values of h.In the error plots, we could also see that the error is largest in the area where only the traction condition is active.In the present work, we manufactured the data for the boundary conditions, but the aim is to eventually use measured pressure values and the solid body rotation of the ribs as boundary input data.This is most similar to the more challenging Benchmark II.There are a large number of transition zones, where the type of boundary condition changes, while we still expect a smooth behaviour of the solution.We were able to achieve a smooth solution and highorder convergence by using the proposed smoothing approach, and we will use that for the real application.Further investigations are needed regarding how to choose smoothing parameters such as the size of the transition zone.
Expiration state with anatomical part (b) Inspiration state with muscle actions

Figure 1 :
Figure 1: Anatomy and physiology of the diaphragm (red and green parts) from expiration to inspiration

Figure 2 :
Figure 2: The expected displacement in different parts of the slice of the diaphragm.The red areas correspond to the real physiological behavior and the green areas are transition zones.In regions 3 and 7, there is both a thickening of the muscle and a downward motion, while in region 5, the downward motion dominates.

1 =
x i ∈ X (one stencil center) and find n closest neighbors x it using the Euclidean distance.2. Scale and shift the stencil points to a unit domain [−1, 1] d .Save the scaling as s(k) .

2 .
Scale y l to a unit domain [−1, 1] d using the previously computed scaling s (k) .3. Form a vector b 1 = L||y l − x (k) j || 3 2 , where {x (k) j } n j=1 is the local neighborhood of the center point x (k) 1 , where L = 1 for constructing evaluation weights, or a derivative operator for constructing differentiation weights.4. Form a vector b 2 = Lp j (y l ), j = 1, .., n 2 .5. Concatenate the two vectors into b(y l ) = [b 1 , b 2 ]. 6. Use the augmented interpolation matrix that belongs to x (k) 1 and compute the local weights by using the relation w L (y l ) = ( Ã(k) ) −1 b L (y l ). 7. Store w L (y l ) in the l-th row of the matrix W L (Y ). 8. Repeat steps 1-6 for every y l ∈ Y in order to form all local weights.

Figure 3 :
Figure 3: The initial point set X 1 (tilted Cartesian points) is distributed over a box that encloses the boundary of the diaphragm (black curve).Points that are more than half a stencil size away from the geometry (grey) are removed.The remaining points (blue) form the point set X.

Figure 4 :
Figure 4: The black curve represents a part of the boundary of the diaphragm.Node points in X (blue markers) and the corresponding Voronoi regions (grey lines) are shown together with interior evaluation points, Y 1 , and boundary evaluation points, Y b , (red markers).The same template of interior evaluation points is used in each Voronoi region inside the diaphragm geometry.
1, and introduce a scaling β i for equations connected with the operator D i .Finally, we form the rectangular system of size 2M × 2N that discretizes the Navier-Cauchy system with Robin boundary conditions (8):

Figure 5 :
Figure 5: The smoothed boundary geometry curve (solid line) is shown in both subfigures.The curve was computed using Ng = 133 node points, Mg = 177 • 3 = 531 data points, and stencil size n = 28.The markers show the Mg = 177 vertex data points (top) and uniform evaluation points (bottom).The normals computed from the approximation, as well as the values of the parameter t along the curve, are also shown in the bottom subfigure.

Figure 6 :
Figure6: The two smoothed coordinate functions approximating the geometry.The function p(t) (left) corresponds to the horizontal coordinate and q(t) (right) corresponds to the vertical coordinate.The markers show the initial data locations.
we display the boundary of the diaphragm before

Figure 7 :
Figure7: Benchmark I: The displacement in the horizontal direction (left) and the vertical direction (right) as a function of the boundary parametrization t (see Figure5, right image, for an illustration of t in relation to the boundary).The markers show the initially placed data points, which are then linearly interpolated into Mg = 80 data points.For the approximation, Ng = 120 node points and stencil size n = 28 were used.The dashed lines show the location of the end points of the diaphragm (regions 2 and 8 in Figure2).

Figure 8 :
Figure 8: Benchmark I: The diaphragm in its non-deformed state (dashed line) and after displacement of the boundary (solid line).

Figure 9 :
Figure 9: Benchmark I:The computed displacements and the corresponding von Mises stress over the diaphragm.This solution was obtained using the unfitted RBF-FD discretization with internodal distance h = 0.004, oversampling parameter q = 5, and an appended polynomial basis of degree p = 5.Due to the scaling applied to the geometry (see Section 6), the displayed results for displacement and stress should be multiplied with s −1 Ω = 0.15692 m to recover the results corresponding to the unscaled problem.

Figure 11 :
Figure 11: Benchmark I: Convergence of the displacements ũ1 and ũ2 and the Von Mises stress for different polynomial degrees p, against a highly resolved numerical solution computed using FEM (the GetDP solver[25]) with linear elements and 236 414 degrees of freedom (corresponding to h = 0.00085).

Figure 12 :
Figure 12: Benchmark I: Error distribution in logarithmic scale for h = 0.006 and p = 5 when a fine unfitted RBF-FD-LS solution (h = 0.002, p = 5) is taken as reference for computing the error.

Figure 13 :
Figure 13: Benchmark I: Error distribution in logarithmic scale for h = 0.006 and p = 5 when a fine unfitted RBF-FD-LS solution (h = 0.002, p = 5) is taken as reference for computing the error.

Figure 14 :
Figure 14: The Robin boundary coefficients displayed over the diaphragm show the impositions of the traction and the Dirichlet parts of the boundary condition.The orange line corresponds to the traction coefficient κ 1 , and the blue line corresponds to the Dirichlet coefficient κ 0 .The dashed lines show the location of the end points of the diaphragm (regions 2 and 8 in Figure2).The parameter t corresponding to the boundary of the diaphragm is illustrated in Figure5.

Figure 15 :
Figure15: Benchmark II: The smoothed boundary data over the diaphragm.For this benchmark we have both Dirichlet data (left column) and traction data (right column), each with a horizontal component (first row) and a vertical component (second row).The markers show the initially placed data points, which are then linearly interpolated into Mg = 80 data points.For the approximation, Ng = 120 node points and stencil size n = 28 were used.The dashed lines show the location of the end points of the diaphragm (regions 2 and 8 in Figure2).The parameter t corresponding to the boundary of the diaphragm is illustrated in Figure5(right image).

Figure 16 :
Figure16: Benchmark II: The computed displacements and the corresponding von Mises stress over the diaphragm.This solution was obtained using the unfitted RBF-FD discretization with internodal distance 0.006, an oversampling parameter q = 5 and an appended polynomial basis of degree p = 5.Due to the scaling applied to the geometry (see Section 6), the displayed results for displacement and stress should be multiplied with s −1 Ω = 0.15692 m to recover the results corresponding to the unscaled problem.

Table 1 :
(8)chmark I and II:The relation between the inverse internodal distance 1/h and (i) the internodal distance h and (ii) the number of interpolation points N used for discretizing the PDE problem(8).