WEB-Spline Finite Elements for the Approximation of Navier-Lamé System with CA,B Boundary Condition

The objective of this article is to discuss the existence and the uniqueness of a weighted extended B-spline(WEB-spline-) based discrete solution for the 2D Navier-Lamé equation of linear elasticity with a different type of mixed boundary condition called CA,B boundary condition. Along with the usual weak mixed formulation, we give existence and uniqueness results for weak solution. Then, we illustrate the performance of Ritz–Galerkin schemes for a model problem and applications in linear elasticity. Finally, we discuss several implementation aspects. The numerical tests confirm that, due to the new integration routines, the weighted B-spline solvers have become considerably more efficient.


Introduction
The finite element method has become the method of choice for solving many types of partial differential equations in engineering and physical sciences. Important applications include structural mechanics, fluid flow, thermodynamics, and electromagnetic fields [1], using the approximation of Lagrange [2]. This type of approximation has experienced a great restriction in the level of the geometric domain, especially in the case of complicated boundaries such as that in the form of curvilinear graphs. A new way of approximation came a few years ago called B-spline that will solve this problem with excellence.
B-splines have become standard tools in approximation, computer graphics, design, and manufacturing. Recently, they have also been used to construct basis functions for finite element methods. The resulting techniques combine the computational efficiency of regular grids with the geometric flexibility of classical finite elements. Some key advantages are the free choice of order and smoothness, a simple data structure with one parameter per grid point, and the exact representation of boundary conditions. A type of Bspline is called weighted extended B-spline (WEB-spline) [3]. As will be described in the next section, a weight function, ω, is being multiplied to the B-spline functions to ensure that all WEB-splines satisfy the boundary conditions exactly. The books [4,5] provide a comprehensive treatment of the relevant theory. Which technique is best suited for a particular problem depends to some extent on the topological form of the simulation domain and its representation. Weighted methods are a good choice for problems with natural (Neumann) boundary conditions or if the part of the boundary, where essential (Dirichlet) boundary conditions are prescribed, has a convenient implicit description. Isogeometric methods can handle domains well which are parametrized over rectangles or cuboids or which can be expressed as union of few such parametrizations, e.g., NURBS representations for CAD/CAM applications. There are also problems where a combination of both techniques might be useful [6].
In this project, we use the implementation of weighted finite element methods to solve the Navier-Lamé system with a different type of mixed boundary condition C A,B [7] that generalizes the well-known basis conditions, especially the Dirichlet and the Neumann conditions. This boundary condition helps us to implement a single MATLAB code that summarizes any kind of boundary conditions that we can meet. At the programming level of the method, we can reap several programs in one. The outline of this article is as follows. In the next section, we introduce the modeling of the Navier-Lamé equation and we will give the variational problem that corresponds to this equation. We review the preliminary aspects on WEB-spline-based methods for Navier-Lamé equations. In Section 3, we prove inf-sup conditions for the variational problem in two-dimensional case using WEB-spline-based mesh-free method, which are necessary for the well posedness of the problem.

Governing Equation
Linear elasticity is the mathematical study of how each point of the solid object is displaced. The latter produces a deformation while the object becomes subjected to internal stresses due to the prescribed loading conditions, knowing that the linear elasticity models materials in a continuous state. Linear elasticity is a simple case of the nonlinear elasticity theory and is a branch of continuous domain mechanics. The basic theory of linear elasticity is based on the following: the strains (or stress) are infinitesimal or small and the relationships between the components of the stress and the strain are linear, and the linear elasticity is valid only for the states of stress which do not produce a yield. These assumptions cited above are reasonable for the analysis of many engineering materials and in technical design. It is often that this analysis is done using the finite element method.
Let us consider Ω ⊂ R 2 to be a bounded Lipschitz domain with boundary condition Γ which will be presented in a new form that generalizes the Neumann and Dirichlet boundary conditions. Given f ∈ L 2 ðΩÞ, A, B ∈ L ∞ ðΓÞ 2×2 , two matricial functions that are invertible or zero. g ∈ H 1 ðΓÞ as well as the positive parameters λ and μ. When μ the shear modulus is small, then this kind of problems is called singularly perturbed problems, where the uniform mesh or L2 norm-based error analysis does not converge to the original problem, and hence, one must use adaptive mesh with uniform error analysis. This numerical analysis is given in the papers [8][9][10][11][12][13][14].
When solid objects are subjected to external or internal loads, they deform and lead to stress. If the deformation of the solid is relatively small, linear relationships between the components of stress and strain are maintained. Consequently, linear elasticity theory is valid. In practice, linear elasticity theory is applicable to a wide range of natural and engineering materials and thus extensively used in structural analysis and engineering design. The equation of Navier-Lamé below is governed as follows.
A solid object is deformed under the action of forces applied. Any fixed point ðx, yÞ of the domain and when surface forces are applied to the domain, this point moves to another point ðX, YÞ. Then, the vector u = ðu 1 , u 2 Þ = ðX − x, Y − yÞ is called displacement. When the movement is small and the solid is elastic, then Hooke's law gives a relationship between the stress tensor and the strain tensor. σ = λtrðεÞ I 2 + 2με is the stress tensor, ε = 1/2ð∇u + ð∇uÞ T Þ is the strain tensor, I 2 is the identity matrix, and μ is the shear modulus (or rigidity), where λ is Lam's first parameter. The Navier-Lamé equation is given by the law of conservation moment ρa = ∇ ·σ where a is the acceleration and ρ is the density of the material. On the other hand, Then, we have With a simple calculation, we find that Then, we get If the solid is in dynamic equilibrium, then we have ρa + f = 0; f are the external forces applied to the solid. Finally, we find out the equation We refer the reader to [15,16] for more information of the elasticity problems.
We create a new unknown ψ = ∇ · u = ð∂u 1 /∂xÞ + ð∂u 2 /∂ yÞ that is equal to the divergence of the displacement. The equation of Navier-Lamé becomes Our mathematical model is the Navier-Lamé system with the boundary condition noted C A,B such that A is called the Dirichlet matrix and B is the Neumann matrix.
There are two strictly positive constants α and β, such that |‖·|‖is a matrix norm that will be defined below. If |‖A|‖≪|‖B|‖, then C A,B is the Neumann boundary; and if |‖B|‖≪|‖A|‖, then C A,B is the Dirichlet boundary.
We need functional spaces and norms V 2 Abstract and Applied Analysis A j j j j j j= max a i,j , i = 1:2, j = 1:2: The variational formulation of the Navier-Lamé problem (6) is as follows.

Solvability of the Generalized Saddle Point System
In this section, we introduce some existing saddle point theory. Let Furthermore, we define the following Hilbert spaces: In addition to assumption (17), we assume that b i , for i = 1, 2, satisfying the inf-sup condition such that there exists a constant ρ > 0 The bilinear form d satisfies the weak coerciveness such as there exists a constant ε > 0 such that Find ðu, ψÞ where a, b i for i = 1, 2 and d are bounded bilinear forms, where f and g are bounded linear functionals on V 0 ðΩÞ and M 0 ðΩÞ, respectively. Reference [17] contains an existence and uniqueness, for the saddle point problem.
Find ðu, ψÞ ∈ V 0 ðΩÞ × M 0 ðΩÞ such that Abstract and Applied Analysis Lemma 1. Let D be a linear operator, defined by such that dðψ, qÞ = ðd ψ , qÞ for all q ∈ D ⊥ . Therefore, D is continuous and we have ImD By definition of the coercivity of the bilinear form d, we have We deduce that ψ j is cauchy in M 0 ðΩÞ; then, there exist ψ ∈ M 0 ðΩÞ such as ψ j ⟶ ψ on M 0 ðΩÞ.
By the continuity of the operator D, we have d ψ j ⟶ d ψ . The uniqueness of the limit gives w = d ψ ∈ D ⊥ , so we obtain that ImD is a closed set of D ⊥ ; it gives that D ⊥ = ImD ⊕ ðImDÞ ⊥ and ImD ∩ ðImDÞ ⊥ = 0.

Theorem 2. With assumptions
Further, the following stability estimates hold: where ðu 0 , ψ 0 Þ is the solution to (23) and thus has the bounds Proof. D is a closed set included in M 0 ðΩÞ, so D is a Hilbert space.
( ð29Þ a, b 1 , and b 2 satisfy the conditions of theorem 3.1 of the article in [17]; then, problem (29) has a unique solution ðu, We are now looking for ψ ∈ M 0 ðΩÞ such that For this u, the displacement solution of (29) is fixed; we consider the following linear form defined by F : q↦−G ðqÞ + b 2 ðu, qÞ. This form inherits its continuity from the continuity of b 2 and G. We define the space D ⊥ = fp ∈ M 0 ðΩÞ \ ðp, qÞ = 0∀q ∈ Dg such that D ⊥ is a Hilbert space with the norm and inner product inherited from M 0 ðΩÞ.
According to the Lax Milgram theorem, problem (31) admits a unique solution ψ ∈ D ⊥ . Assume that (31) is verified, let us show that (30) is verified.
D is closed in M 0 ðΩÞ, so we have M 0 ðΩÞ = D ⊕ D ⊥ ; then, we obtain for all q ∈ M 0 ðΩÞ. Therefore, we have q = q 1 + q 2 for q 1 ∈ D and q 2 ∈ D ⊥ ; it implies that dðψ, qÞ = dðψ, Conversely, let ψ ∈ M 0 ðΩÞ be the solution of (30). The linear form q ↦ dðψ, qÞ is continuous on D ⊥ that is a Hilbert space with the scalar product on M 0 ðΩÞ.
From the Riesz theorem, there exists a unique d ψ ∈ D ⊥ such that This defines an operator D : ψ ↦ d ψ which is linear and continuous on M 0 ðΩÞ such that ImD = fd ψ \ ψ ∈ M 0 ðΩÞg. From Lemma (1), we have ImD = D ⊥ .

Theorem 3. With assumptions
Further, the following stability estimates hold: where ðu 0 , ψ 0 Þ is the solution to (23) and thus has the bounds Proof. It is sufficient to show that the bilinear forms a, b, b Γ , and d have checked the conditions (17)- (21), applying Theorem 2 in the case of G = 0.
Remark 4. We rewrite the system (15) in a single equation. For that, let us define the operators A :  (15) becomes The well posedness [17][18][19] of the above equation is written in the following form.

WEB-Spline Process
Several types of mesh-less approximations were proposed for various applications. A central problem of the mesh-less Galerkin method is to incorporate boundary conditions of Dirichlet type. The weighted extended B-spline approximation takes care of not only the boundary constraints but also the issue of well conditioning of the Galerkin systems [3,20]. This essential new feature of constructing well-conditioned basis is of cardinal importance in this work. We are interested in applying the approximation properties of the WEB-spaces for discretizing the Navier-Lamé equations as it have been done for Stokes and Navier-Stokes problems [21,22].
The standard uniform B-spline of degree n is defined by the recursion [5] starting from B 0 , the characteristic function of the unit interval between zero and one. Figures 1-3 show the uniform B-splines of degrees one, two, and three. These are also known as linear, quadratic, and cubic B-splines, respectively.
In this paper, we use the following notational conventions [5]. For functions f and g, we write f ≺g if f ≤ cg with a constant c which does not depend on the grid width h, indices, or arguments of functions. The symbols ≻ and ≈. are defined analogously. We describe the procedure of constructing the WEB-splines and discuss the approximating properties of the web-space as a finite element space. For k = ðk 0 , k 1 Þ ∈ ℤ 2 , and h > 0 define the scaled translates: B is the univariate B-spline of degree n of support ½0, n + 1Þ; the B-splines B k are polynomials by pieces on the h-grid with vertices hℤ 2 and scaled so that the L 2 -norm, Knowing that the tensor product B-spline is the extension of B-spline to higher dimensions.
In general, it can be defined as follows: such that the B-spline B n k,h is an m-variant of degree n ν in the νth variable, index k = k 1 , ::, k m and the width of the grid h.
With the convention that n 1 = ⋯ = n m unless otherwise indicated, for the sake of the problem, the element n can be taken as an integer, rather than an integer vector.
For I ∈ f0, ⋯, ng m , we partition the grid cells Q I = I h + ½0, 1 m h with m = 2 into interior, boundary, and exterior cells, depending on whether Q ⊂ Ω, the interior of Q intersects ∂Ω, For k ∈ K ≔ fI ∈ ℤ 2 : supp ðB I Þ ∩ Ω ≠ ∅g (the relevant index set for Ω), if supp ðB k Þ has at least one grid cell completely inside Ω, then bk is named as inner B-spline; otherwise, it is outer.
The corresponding subsets of K are I and J (c.f Figure 4): More details about the construction of WEB-spline basis are given in [5]. Now, it is tempting to use B h ≔ spanfB k : k ∈ Kg as a finite element approximation space. At first sight, this does not seem feasible since B-splines do not conform to the boundary conditions. But this difficulty can be resolved easily by multiplying B k by a smoothed distance function ωðxÞ~di   Abstract and Applied Analysis power and can simply be omitted. Unfortunately, this is not the case. A suitable solution to the problem of controlling the unstable outer B-splines is provided by adjoining them appropriately with the inner B-splines. It has been done in such a way that the approximation power of the finite element subspace is retained.
Definition 5 (see [20]). For i ∈ I, WEB-spline B i is defined as 7 Abstract and Applied Analysis where x i denotes the center of a grid cell, which is completely inside the domain Ω. The coefficients e i,j satisfy je i,j j≺1, e i,j = 0 for ∥i − j∥±1 and are chosen so that all weighted polynomials ðωpÞ of order n are contained in the WEB-space B h ≔ spanfB i : i ∈ Ig.
Theorem 6 (see [20]). For an outer index j ∈ J, let IðjÞ = l + f0, ⋯, ng m ⊂ I be an m-dimensional array of inner indices closest to j assuming that h is small enough so that such an array exists. Then, the coefficients are admissible for constructing WEB-splines according to Definition 5.

C A,B Boundary Conditions with Resolution Adapted
The imposition of inhomogeneous Dirichlet boundary conditions is essential in numerical analysis of a structure. It is especially difficult and no longer straightforward whenever nonconformal mesh is used to discretize a structure. One of the contributions of this paper is to develop a weighted extended spline basis with high computing accuracy and appropriate with the mixed formulation of boundary conditions C A,B that generalize all the cases that can be encountered on the edge of Ω.
With A and B are two reversible square (in case 2D) matrices, A, B ∈ L ∞ ðΓÞ 2×2 and g ∈ H 1/2 ðΓÞ where the function g is a priori known function of spatial coordinates. The inhomogeneous Dirichlet condition is expressed in the case of B = 0; then, we obtain u = A −1 g, on Γ. In general, a boundary condition function g may not be known explicitly or is not prescribed globally. More typically, the boundary conditions are prescribed in a piecewise manner, i.e., the boundary conditions are specified by individual functions g i on each portion of boundary, In fact, as the weighting function and boundary value function are expressed in the form of implicit functions, no unique expressions exist but typical forms are given below.
Weighting function ωðxÞ acts as a multiplier to modify the original interpolation functions and has to vanish on any Dirichlet boundary ðω i ðxÞ = 0Þ. As implicit function is not unique for the definition, weighting function of different forms can be used, as discussed in [24]. Here, two construc-tion techniques are presented. First, a straightforward method is to construct the weighting function by means of the product ωðxÞ = Q m i=1 ω i ðxÞ. Clearly, ωðxÞ = 0 if any ω i ðxÞ = 0 on Γ i . Nevertheless, a product of all the implicit functions might lead to a surge of function values which will cause a numerical overflow and then result in poor computing accuracy and robustness. Therefore, each implicit function ω i is preferred to be normalized in advance according to the size of the physical domain. For instance, the implicit function of a circle can be normalized in terms of its radius.
Now, we will choose a basis function web-spline that will be called the web-spline basis function adapted. This basis will respond to any type of nonhomogeneous boundary conditions especially C A,B in the case B = 0. Let u be the displacement solution of problem (15), so we can write u = u 0 + u Γ with u 0 being the solution of (15) in case of homogeneous Dirichlet boundary and u Γ = A −1 g is the nonhomogeneous Dirichlet boundary. u is written as a linear combination of the B i web-spline family cited in Definition (5) It can be rewritten again in the form where b i ðxÞ denotes the ith weighting coefficient associated with A −1 g i in which basic properties of b i ðxÞ are as follows [23]. b i = δ ij , i, j = 1, 2 ⋯ N, with δ ij denoting the Kronecker delta function. With this condition, different forms are proposed below for the definition of b i .
Transfinite interpolation form: the weighting coefficients can be defined by extending the transfinite interpolation that was studied by Rvachev et al. [25] in the CAD community.
The above equation shares a property of symmetry and similarity. ω j ðxÞ is generally supposed to be positive in the physical domain to ensure the nonzero value of the denominator. Evidently, b i = 1 only if ω i = 0 and the partition of unity holds with ∑ N i=1 b i = 1. The value of θ is used to interpolate normal derivatives prescribed on Γ, which must be 1 greater than the order of the prescribed derivatives. 8 Abstract and Applied Analysis Finally, the basis that we can propose for the nonhomogeneous boundary condition is as follows: 6. Discretization of the Navier-Lamé Problem Using WEB-Spline Basis In the Navier-Lamé equation, Δu and ∇ψ are the terms with derivatives of the highest order for the displacement and div-displacement ðψÞ, respectively. Thus, the orders of the differential operators differ by 1. This suggests the rule of thumb. The degree of the basis functions used to approximate the displacement should be one larger than the approximation of the div-displacement. Also, to satisfy the Dirichlet boundary conditions, the displacement basis functions are multiplied by a suitable weight function ω. Hence, in this article, we choose ϕ j -linear weighted extended B-spline for displacement approximation along with φ i -Haar wavelet basis function for div-displacement approximation as our linear-constant element and for quadratic-linear element we denote ϕ j as quadratic WEB-spline and φ i the mean zero linear function (see Figures 1 and 2) as defined above. In the following, the inf-sup condition and, therefore, the well posedness of the discrete Navier-Lamé problem are settled first for the linear-constant element and then for the quadratic-linear element.

Linear Constant Element.
More precisely, we give V h and M h , the displacement and div-displacement finite element spaces, respectively, as follows: A sufficient condition for ðV h , M h Þ to satisfy the inf-sup condition is given in [21,22]. A similar result has been proved in [26] for a different pair ðV h , M h Þ.

Quadratic-Linear Element. For
where φ is defined as x − 4 on 3, 4 ½ Þ: For constructing the displacement approximation space, we do the following. For i ∈ I u , let ϕ i be the WEB-spline of order 3, which is given by (40) and (48) where ϕ i is defined as in (51), the tensor product of the scaled translate of the function ϕ ϕ x ð Þ ≔ x 2 2 on 0, 1 ½ Þ, where the quadrangulation T h is the collection of all cells K ⊂ Ω such that K ∩ ∂Ω = ∅ (those which are fully inside the domain Ω) and K ∂Ω ≔ K ∩ Ω ≠ ∅ such that K ∩ ∂Ω ≠ ∅ (those portions which intersect the boundary ∂Ω) and V K is the one dimensional subspace spanned by the function b K given by The function bðxÞ ≔ xðx − 1Þ is chosen so that the bubble function b K vanishes on the edges of the cell K. ðk 0 , k 1 Þ is the node corresponding to the cell K. The weight function w is multiplied to ensure that the bubble function also vanishes on the boundary ∂Ω. The distance function is used near the boundary, where it is free of singularities and blended smoothly with a plateau inside the domain. More precisely, ω is defined as ω = 1 − ðmax ðδ − dðxÞ, 0Þ/δÞ l , dðxÞ = distðx, ∂ΩÞ, where δ controls the height of the plateau and l the smoothness of the weight function. The plateau facilitates the use of precomputed values while assembling the Galerkin matrix and also avoids the use of high-order quadratures for the integration of the bubble functions supported on the cells which are fully inside the domain.
The displacement approximation space is taken to be This pair ðV 2 h , M h Þ of discretization spaces satisfies the discrete inf-sup condition.
Lemma 7 (see [3]). A web-basis fB i g i∈I is a stable basis with respect to the L 2 -norm, We can also bound higher order Sobolev norms in terms of the 2-norm of the coefficients.
Lemma 8 (see [3]). A web-basis fB i g i∈I satisfies ∥∑ i∈I a i B i ∥ 1 ≤ C 3 h −1 ∥fa i g i∈I ∥, and C 3 > 0 is a constant.
Lemma 9 (inverse estimate). Let V h be the finite element space considered as in (49). Then, there exists a constant C such that ∥∑ i∈I a i B i ∥ 1 ≤ C 6 h −1 ∥∑ i∈I a i B i ∥ 0 .
The proof follows from Lemmas 7 and 8.
Lemma 10 (see [21,22]). Suppose T h is a family of uniform quadrilateral mesh on the domain Ω and Q h is the L 2 -projection operator onto V h ⊂ H 1 0 ðΩÞ. Then 9 Abstract and Applied Analysis where C is a constant which is h-independent.
Lemma 11. There exists an operator P h : H 1 ðΩÞ ⟶ V 2 h which satisfies the following properties: Proof. Let P 0 h : H 1 ðΩÞ ⟶ V h be the usual L 2 -projection operator with The first inequality (59) follows from Lemma (10) We may interpret the map P 1 h as a process with two steps. First, we apply the L 2 -projection onto the space of piecewise constant functions. Afterwards, in each cell K the constant is replaced by a bubble function with the same integral. In this way, we get More precisely, we define P 1 h as follows: In view of condition (61), the constant βðKÞ is taken to be where c 4 , c 5 > 0 are constants, The above inequality is due to Cauchy-Schwarz and |K| = h 2 the area of the cell K. Let us set Verify the two properties in the statement of the lemma. By the virtue of the construction of P h and we impose that Because q h is piecewise linear, we apply Green's theorem It is easy to prove (57) by applying Green's theorem and using (56) This proves the assertion (58) with c = c 1 + c 6 c 3 c 2 , such as c 6 of Lemma (9).

Theorem 12. Let the space be
The pair ðV 2 h , M h Þ satisfies the discrete inf-sup condition: Proof. Let P h v ∈ B ⊥ h be given. By the continuous form of the inf-sup condition (20) and (56)

Matrix Form of the Navier-Lamé Problem
After establishing the discrete inf-sup condition, we are ready to prove the existence and uniqueness of the discrete solution. (This part of the section is common to both the linear-constant and quadratic-linear elements, so we call the finite element spaces as V h and M h in general.) To find the discrete solution pair ðu h , ψ h Þ ∈ V h × M h , it is enough to 10 Abstract and Applied Analysis where u α = ðð u α i Þ i∈I u Þ T , α = 1, 2, and div-displacement ψ ≔ ðð ψ j Þ j∈I ψ Þ T . Assuming that ðB −1 AÞ αβ = ξ αβ , and ðB −1 Þ αβ = ζ αβ for α, β = 1, 2. Let us define The matrix form of the discrete Navier-Lamé equations is written as By defining the assembled matrices A = DiagðA 1 , A 2 Þ, Decoupling of div-displacement and displacement: from the above matrix form, we write Premultiplying the first equation above by A −1 , followed by B T , and using the second equation, we get In the following, we show that D + B T A −1 B Γ is invertible, and hence, the above equation can be solved for ψ. The positive definiteness follows from the positive D + B T A −1 B Γ : where M = cardinality of I ψ .
The positive definiteness of B T A −1 B follows from the positive definiteness of matrix A (which is guaranteed by the coercivity condition) and the following inequality which is an equivalent condition to the inf-sup condition (20).

Numerical Results
We tested the performance of the WEB-spline method for the Navier-Lamé problem (6). We take the centrifugal force ∥f ∥ = r and the domain Ω is a rotating steel disc defined by Ω = fðx, yÞ/ðx − ð1/2ÞÞ 2 + ðy − ð1/2ÞÞ 2 ≤ ð1/2Þ 2 g \ fðx, yÞ/ ðx − ð1/2Þ + dÞ 2 + ðy − ð1/2ÞÞ 2 ≤ r 2 g, with r = 0:1, d = 0:06, and the edge of Ω consists of two disjoint parts in the form of two circles We give in (42) the Dirichlet matrix A = I 2 on Γ D , 0 2 on Γ N , and we give the Neumann matrix B = I 2 on Γ N , 0 2 on Γ D , the traction force g = ð0, 0Þ on ∂Ω. Figures 5 and 6 show the numerical solution for rotating steel disc with ðE = 1, ν = 0:28Þ that is fixed at a slightly eccentric axis with radius r. We show the domain as well as direction and size (colored) of displacement for each H = ð8,16,32,64Þ number of cells per coordinate direction.
As mentioned in Table 1, we check for this example the accuracy of the numerical approximations. Since in view of the smoothness of the weighted B-spline basis functions that are continuously differentiable for n ≥ 2 we can compute the pointwise residuals [27] for the Navier-Lamé boundary value problem (6), error pde = ∥−μΔu h − ðλ + μÞ∇ψ h − f ∥ 0 and rate r = ½NaN − diff ðlog ðeÞÞ/log ð2Þ. We can substitute the finite element approximation into the partial differential equation Navier-Lamé where g = 0 for the rotating disc. This would not be possible for standard finite element approximation which usually merely belong to H 1 ðΩÞ.
We show the relative residual for grid widths 1/H = 1/8, ⋯, 1/256. While the residuals decay, estimates of the rates exhibit a regular increasing behavior when H is

12
Abstract and Applied Analysis changed from 64 to 256 (e.g., the error pde norm [27] has dropped from 1:118e + 000 to 4:462e − 004 when H is changed from 8 to 256) which we attribute to numerical stabilities. Note that in the case when μ is quite small, we will refer the readers to see [28][29][30][31].

Conclusion
Using the equations of linear elasticity as a model problem with the boundary condition C A,B , we have described the implementation of finite element methods with weighted Bsplines using mixed finite element. In this article, we have shown that our WEB-spline-based quadratic-linear finite elements satisfy the inf-sup condition, which is necessary for the existence and uniqueness of the solution, and we proved the existence of the discrete solution. We have shown the convergence of the numerical solution for the quadratic case. Because of limited regularity, the Navier-Lamé problem will not change by increasing the degree of the WEB-spline.
We have computed the relative errors and rates and have shown that it is of order 1/H, which is theoretically correct; these numerical results are attributed to the numerical stabilities of the solution. Moreover, the advantage of this problem with C A,B boundary condition is the program level MATLAB; it is enough to make a single program MATLAB and can be reduced to ordinary problems as Dirichlet and Neumann.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that they have no conflicts of interest.