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

Authors: Ouadie Koubaiti; Ahmed Elkhalfi; Jaouad El-mekkaoui
Journal: Abstract and Applied Analysis
 (2020)

---

## Abstract

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 calledCA,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.


---

## Body

## 1. 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 B-spline 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 conditionCA,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.


## 2. 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Ω⊂R2 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∈L2Ω, A,B∈L∞Γ2×2, two matricial functions that are invertible or zero. g∈H1Γ 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–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 pointx,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=u1,u2=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εI2+2με is the stress tensor, ε=1/2∇u+∇uT is the strain tensor, I2 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,
(1)∇·σ=λ∇trεI2+2μ∇·ε.Then, we have(2)∇·σ=λ∇trεI2+μ∇∇u+μ∇∇ut.With a simple calculation, we find that
(3)∇trεI2=∇∇ut=∇∇·u.Then, we get(4)ρa=μΔu+λ+μ∇∇·u.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
(5)f=−μΔu−λ+μ∇∇·u.We refer the reader to [15, 16] for more information of the elasticity problems.We create a new unknownψ=∇·u=∂u1/∂x+∂u2/∂y that is equal to the divergence of the displacement. The equation of Navier-Lamé becomes
(6)−μΔu−λ+μ∇ψ=finΩ,ψ−∇·u=0inΩ,Au+Bμ∂u∂n+λ∇·un=gonΓ.Our mathematical model is the Navier-Lamé system with the boundary condition notedCA,B such that A is called the Dirichlet matrix and B is the Neumann matrix.There are two strictly positive constantsα and β, such that
(7)αu·v<B−1Au·v<βu·v∀u,v∈R2.∣‖·∣‖is a matrix norm that will be defined below.If∣‖A∣‖≪∣‖B∣‖, then CA,B is the Neumann boundary; and if ∣‖B∣‖≪∣‖A∣‖, then CA,B is the Dirichlet boundary.We need functional spaces and norms(8)h1Ω=u:Ω⟶R\u,∂u∂x,∂u∂y∈L2Ω,(9)VΩ=H1Ω=h1Ω2,(10)MΩ=L02Ω=q∈L2Ω\∫Ωq=0,(11)∥v∥1,Ω=∫Ω∇v:∇vdΩ+∫Ωv.vdΩ1/2,(12)∥v∥0,Ω=∫Ωv.vdΩ1/2,(13)A=maxai,j,i=1.2,j=1.2.The variational formulation of the Navier-Lamé problem (6) is as follows.Findu,ψ∈VΩ×MΩ such that
(14)∫Ωμ∇u:∇vdΩ+∫ΓB−1Au·vdΓ,+∫Γμψn·vdΓ+∫Ωλ+μψ∇·vdΩ,=∫Ωf·vdΩ+∫ΓB−1g·vdΓ,∫Ωλ+μq∇·udΩ−∫Ωλ+μψqdΩ=0.The weak formulation (14) may be restated as follows.Findu,ψ∈VΩ×MΩ
(15)au,v+bΓv,ψ=Lv∀v∈VΩ,bu,q−dψ,q=0∀q∈MΩ.With the bilinear forms(16)au,v=∫Ωμ∇u:∇vdΩ+∫ΓB−1Au·vdΓ,bv,q=∫Ωλ+μq∇·vdΩ,bΓv,q=bv,q+∫Γμqn·vdΓ,dψ,q=∫Ωλ+μψqdΩ,Lv=∫Ωf·vdΩ+∫ΓB−1g·vdΓ.


## 3. Solvability of the Generalized Saddle Point System

In this section, we introduce some existing saddle point theory. LetV0Ω and M0Ω be two finite- or infinite-dimensional Hilbert spaces equipped with the inner products .,.V0Ω and .,.M0Ω and the induced norms ∥·∥1,Ω and ∥·∥0,Ω, respectively. Let au,v, b1v,q and b2v,q, dp,q be bilinear forms on V0Ω×V0Ω, V0Ω×M0Ω and V0Ω×M0Ω, M0Ω×M0Ω, respectively, which are bounded; i.e., there are positive constants ∥a∥, ∥b1∥, ∥b2∥, and ∥d∥ such that
(17)au,v≤∥a∥∥u∥1,Ω∥v∥1,Ω∀u,v∈V0Ω×V0Ω,b1u,q≤∥b1∥∥u∥1,Ω∥q∥0,Ω∀u,q∈V0Ω×M0Ω,b2u,q≤∥b2∥∥u∥1,Ω∥q∥0,Ω∀u,q∈V0Ω×M0Ω,dq,p≤∥d∥∥q∥0,Ω∥p∥0,Ω∀p,q∈M0Ω×M0Ω,Lv≤∥L∥∥v∥1,Ω∀v∈V0Ω,Gq≤∥G∥∥q∥0,Ω∀q∈M0Ω.Furthermore, we define the following Hilbert spaces:(18)Bi=u∈V0Ω\biu,q=0∀q∈M0Ω,fori=1,2,D=p∈M0Ω\dp,q=0∀q∈M0Ω.In addition to assumption (17), we assume that
(19)av,v≥δ∥v∥1,Ω2,forallv∈V0Ω,bi, for i=1,2, satisfying the inf-sup condition such that there exists a constant ρ>0
(20)sup∥q∥=1biv,q≥ρi∥v∥1,Ω∀v∈Bi⊥,supv∈Bi⊥biv,q>0∀q∈M0Ω.The bilinear formd satisfies the weak coerciveness such as there exists a constant ε>0 such that
(21)supq∈M0Ωdp,q≥ε∥p∥0,Ω2∀p∈M0Ω.Findu,ψ∈V0Ω×M0Ω
(22)au,v+b1v,ψ=Lv∀v∈V0Ω,b2u,q−dψ,q=Gq∀q∈M0Ω,where a, bi for i=1,2 and d are bounded bilinear forms, where f and g are bounded linear functionals on V0Ω and M0Ω, respectively. Reference [17] contains an existence and uniqueness, for the saddle point problem.Findu,ψ∈V0Ω×M0Ω such that
(23)au,v+b1v,ψ=Lv∀v∈V0Ω,b2u,q=Gq∀q∈M0Ω.Lemma 1.
LetD be a linear operator, defined by
(24)D:M0Ω⟶D⊥,ψ↦dψ,such that dψ,q=dψ,q for all q∈D⊥. Therefore, D is continuous and we have ImD=D⊥.Proof.
ImD is closed in D⊥. In fact, let ψjj be a sequence of M0Ω such as dψj⟶w on D⊥; then, dψjj is cauchy on D⊥.
By definition of the coercivity of the bilinear formd, we have
(25)∥ψj−ψk∥0,Ω≤1εsupq∈M0Ωdψj−ψk,q∥q∥0,Ω=1εsupq∈D⊥dψj−ψk,q∥q∥0,Ω=1ε∥dψj−ψk∥0,Ω=1ε∥dψj−dψk∥0,Ω.
We deduce thatψj is cauchy in M0Ω; then, there exist ψ∈M0Ω such as ψj⟶ψ on M0Ω.
By the continuity of the operatorD, 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.
Suppose thatImD⊥≠0; this implies that there exists p∈ImD⊥ nonzero. On the other hand, we have for all q∈M0Ωdq,p=dq,p=0 (because dq∈ImD); this implies that p∈D but p∈ImD⊥⊂D⊥; this contradicts the fact that p is nonzero. Finally, ImD=D⊥.Theorem 2.
With assumptions (17)–(21), the generalized variational problem (22) has a unique solution u,ψ∈V0Ω×M0Ω for any L∈V0Ω′ and G∈M0Ω′ as long as
(26)δ1=δε−1∥b1∥∥b2∥<1.
Further, the following stability estimates hold:(27)∥u∥1,Ω≤δ1+11−δ1∥u0∥1,Ω+δ−1∥b1∥1−δ1∥ψ0∥0,Ω,∥ψ∥0,Ω≤ε−1∥b2∥1−δ12∥u0∥1,Ω+δ−1∥b1∥∥ψ0∥0,Ω,where u0,ψ0 is the solution to (23) and thus has the bounds
(28)∥u0∥1,Ω≤c11∥L∥+c12∥G∥,∥ψ0∥1,Ω≤c21∥L∥+c22∥G∥,c11=δ−1,c12=ρ2−1δ−1∥a∥+1,c21=ρ2−1c11∥a∥+1,c22=ρ2−1∥a∥c12.Proof.
D is a closed set included in M0Ω, so D is a Hilbert space.
Problem (22) became
(29)au,v+b1v,ψ=Lv∀v∈V0Ω,b2u,q=Gq∀q∈D.
a, b1, and b2 satisfy the conditions of theorem 3.1 of the article in [17]; then, problem (29) has a unique solution u,ψ0∈VΩ×D (if D is reduced to 0, just take ψ0=0). We are now looking for ψ∈M0Ω such that
(30)dψ,q=−Gq+b2u,q∀q∈M0Ω.
For thisu, the displacement solution of (29) is fixed; we consider the following linear form defined by F:q↦−Gq+b2u,q. This form inherits its continuity from the continuity of b2 and G. We define the space D⊥=p∈M0Ω\p,q=0∀q∈D such that D⊥ is a Hilbert space with the norm and inner product inherited from M0Ω. Problem (30) is equivalent to the problem
(31)dψ,q=Fq∀q∈D⊥.
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 M0Ω, so we have M0Ω=D⊕D⊥; then, we obtain for all q∈M0Ω. Therefore, we have q=q1+q2 for q1∈D and q2∈D⊥; it implies that dψ,q=dψ,q1+dψ,q2=0+Fq2=F(q1+Fq2.
Conversely, letψ∈M0Ω be the solution of (30).
The linear formq↦dψ,q is continuous on D⊥ that is a Hilbert space with the scalar product on M0Ω.
From the Riesz theorem, there exists a uniquedψ∈D⊥ such that
(32)dψ,q=dψ,q∀q∈D⊥,∥dψ∥0,Ω=∥dψ,.∥D⊥′=supq∈D⊥dψ,q∥q∥0,Ω.
This defines an operatorD:ψ↦dψ which is linear and continuous on M0Ω such that ImD=dψ\ψ∈M0Ω. From Lemma (1), we have ImD=D⊥.
IfF∈D⊥′, Riesy implies that there exists w∈D⊥ such that Fq=w,qM0Ω for all q∈D⊥. Since ImD=D⊥, for this w∈D⊥, there exists ψ∈M0Ω such that w=dψ dψ,q=ddψ,q=w,q=Fq for all q∈D⊥.
Finally, it has been shown that problems (30) and (31) are equivalent.
Lax Milgram guarantees that problem (31) admits a unique solution ψ1∈D⊥. Just take ψ=ψ0+ψ1∈M0Ω=D⊕D⊥.Theorem 3.
With assumptions (17)–(21), the generalized variational problem (15) has a unique solution u,ψ∈V0Ω×M0Ω for any L∈V0Ω′ as long as
(33)δ1≔δε−1∥bΓ∥∥b∥<1.
Further, the following stability estimates hold:(34)∥u∥1,Ω≤δ1+11−δ1∥u0∥1,Ω+δ−1∥bΓ∥1−δ1∥ψ0∥0,Ω,∥ψ∥0,Ω≤ε−1∥b∥1−δ12∥u0∥1,Ω+δ−1∥bΓ∥∥ψ0∥0,Ω,where u0,ψ0 is the solution to (23) and thus has the bounds
(35)∥u0∥1,Ω≤c11∥L∥,∥ψ0∥1,Ω≤c21∥L∥,c11=δ−1,c21=ρ2−1c11∥a∥+1.Proof.
It is sufficient to show that the bilinear formsa, 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:V⟶V∗, B:V⟶M∗, BΓ:V→M∗, and D:M⟶M∗ by Au,v≔au,v,v∈VΩ, Bu,ψ≔bu,ψ,ψ∈MΩ, BΓu,ψ≔bΓu,ψ,ψ∈MΩ, and Dψ,q≔dψ,q,q∈MΩ so that (15) becomes
(36)ŁU≔ABΓBT−Duψ=L0≕F.The well posedness [17–19] of the above equation is written in the following form. For the variational problem (15), the mapping Ł:V×M→X∗×M∗ defines an isomorphism if and only if a, b, bΓ, and d satisfy assumptions (17)–(21).


## 4. 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 degreen is defined by the recursion [5]
(37)Bnx=xnBn−1x+n+1−xnBn−1x−1,starting from B0, 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=k0,k1∈ℤ2, and h>0 define the scaled translates:
(38)Bkx≔1hBxh−k0Byh−k1.Figure 1
Linear B-spline.Figure 2
Quadratic B-spline.Figure 3
Cubic B-spline.B is the univariate B-spline of degree n of support 0,n+1; the B-splines Bk are polynomials by pieces on the h-grid with vertices hℤ2 and scaled so that the L2-norm, ∥Bk∥=∥B0∥0,Ω, is independent of h.Knowing that the tensor product B-spline is the extension of B-spline to higher dimensions.In general, it can be defined as follows:(39)Bk,hn=Πν=1mBkν,hnν,such that the B-spline Bk,hn is an m-variant of degree nν in the νth variable, index k=k1,..,km and the width of the grid h.With the convention thatn1=⋯=nm unless otherwise indicated, for the sake of the problem, the element n can be taken as an integer, rather than an integer vector.ForI∈0,⋯,nm, we partition the grid cells QI=Ih+0,1mh with m=2 into interior, boundary, and exterior cells, depending on whether Q⊂Ω¯, the interior of Q intersects ∂Ω, or Q∩Ω=∅.Fork∈K≔I∈ℤ2:suppBI∩Ω≠∅ (the relevant index set for Ω), if suppBk has at least one grid cell completely inside Ω, then bk is named as inner B-spline; otherwise, it is outer.The corresponding subsets ofK are I and J (c.f Figure 4): K=I∪J.Figure 4
Relevant quadratic B-splines for domainA, marked at the center of their supports—inner B-splines with circles and outer ones with squares.More details about the construction of WEB-spline basis are given in [5].Now, it is tempting to useBh≔spanBk:k∈K 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 Bk by a smoothed distance function ωx~distx,∂Ω. Then, the weighted B-spline space (WB-space) Bωh≔spanωBk:k∈K spanned by weighted B-splines is a possible finite element subspace for Dirichlet boundary value problems yielding optimal order approximations. But, the condition number of the Galerkin matrix Gh can become extremely large. This is due to the outer B-splines which have only very small support inside Ω. One might think that these basis functions do not contribute much to the approximation 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]).
Fori∈I, WEB-spline Bi is defined as
(40)Bi=ωωxibi+∑j∈Jei,jbj,where xi denotes the center of a grid cell, which is completely inside the domain Ω. The coefficients ei,j satisfy ei,j≺¯1, ei,j=0 for ∥i−j∥±1 and are chosen so that all weighted polynomials ωp of order n are contained in the WEB-space Bh≔spanBi:i∈I.Theorem 6 (see [20]).
For an outer indexj∈J, let Ij=l+0,⋯,nm⊂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
(41)ei,j=∏ν=1m∏μ=0,lν+μ≠iνnjν−lν−μiν−lν−μare admissible for constructing WEB-splines according to Definition 5.


## 5.CA,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 conditionsCA,B that generalize all the cases that can be encountered on the edge of Ω.
(42)Au+Bμ∂u∂n+λ∇·un=gonΓ.WithA and B are two reversible square (in case 2D) matrices, A,B∈L∞Γ2×2 and g∈H1/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−1g,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 gi on each portion of boundary, u=A−1gi,onΓi, i=1,⋯N, and Γ=∪i=1NΓi, ∩iΓi=∅. If each Dirichlet boundary Γi is geometrically represented by implicit function ωi so that ωi=0 on Γi [23],
(43)u=A−1gi,∀ωi=0.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 ωix=0. As implicit function is not unique for the definition, weighting function of different forms can be used, as discussed in [24]. Here, two construction techniques are presented. First, a straightforward method is to construct the weighting function by means of the product ωx=∏i=1mωix. Clearly, ωx=0 if any ωix=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.
(44)ωx=1−x−x0R2−y−y0R2.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 especiallyCA,B in the case B=0. Let u be the displacement solution of problem (15), so we can write u=u0+uΓ with u0 being the solution of (15) in case of homogeneous Dirichlet boundary and uΓ=A−1g is the nonhomogeneous Dirichlet boundary. u is written as a linear combination of the Bi web-spline family cited in Definition (5)
(45)u=∑i=1NuiBi+A−1g.It can be rewritten again in the form(46)u=∑i=1NuiBi+∑i=1NbiA−1gi,where bix denotes the ith weighting coefficient associated with A−1gi in which basic properties of bix are as follows [23].bi=δ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 bi.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.
(47)bi=Πj=1,i≠jNωjθ∑k=1NΠj=1,i≠jNωjθ.The above equation shares a property of symmetry and similarity.ωjx is generally supposed to be positive in the physical domain to ensure the nonzero value of the denominator. Evidently, bi=1 only if ωi=0 and the partition of unity holds with ∑i=1Nbi=1. The value of θ is used to interpolate normal derivatives prescribed on Γ, which must be 1 greater than the order of the prescribed derivatives.Finally, the basis that we can propose for the nonhomogeneous boundary condition is as follows:(48)ϕii=1,..N=Bi+bii=1,..N.


## 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.


### 6.1. Linear Constant Element

More precisely, we giveVh and Mh, the displacement and div-displacement finite element spaces, respectively, as follows:
(49)Vh=u1,u2/uk=Σj∈Iuαjkϕj,αjk∈R,k=1,2,(50)Mh=ψ/ψ=Σi∈Iψβikφi,βi∈R.A sufficient condition forVh,Mh to satisfy the inf-sup condition is given in [21, 22]. A similar result has been proved in [26] for a different pair Vh,Mh.


### 6.2. Quadratic-Linear Element

Forj∈Iψ,
(51)φjx,y=1hφxh−j0φyh−j1,where φ is defined as
(52)φx≔xon0,1,2−xon1,2,2−xon2,3,x−4on3,4.For constructing the displacement approximation space, we do the following. Fori∈Iu, 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 ϕ
(53)ϕx≔x22on0,1,12+x−1−x−12on1,2,3−x22on2,3,Vb≔⊕K∈ThVK,where the quadrangulation Th 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 VK is the one dimensional subspace spanned by the function bK given by
(54)bK≔ωx,ybxh−k0byh−k1.The functionbx≔xx−1 is chosen so that the bubble function bK vanishes on the edges of the cell K. k0,k1 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δ−dx,0/δl, dx=distx,∂Ω, 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 beVh2≔Vh⊕Vb.This pairVh2,Mh of discretization spaces satisfies the discrete inf-sup condition.Lemma 7 (see [3]).
A web-basisBii∈I is a stable basis with respect to the L2-norm, C1∥aii∈I∥≤∥∑i∈IBi∥0≤C2∥aii∈I∥, where C1,C1 are constants.
We can also bound higher order Sobolev norms in terms of the 2-norm of the coefficients.Lemma 8 (see [3]).
A web-basisBii∈I satisfies ∥∑i∈IaiBi∥1≤C3h−1∥aii∈I∥, and C3>0 is a constant.Lemma 9 (inverse estimate).
LetVh be the finite element space considered as in (49). Then, there exists a constant C such that ∥∑i∈IaiBi∥1≤C6h−1∥∑i∈IaiBi∥0.
The proof follows from Lemmas7 and 8.Lemma 10 (see [21, 22]).
SupposeTh is a family of uniform quadrilateral mesh on the domain Ω and ∏h is the L2-projection operator onto Vh⊂H01Ω. Then
(55)∥∏hv∥1≤C∥v∥1,∀v∈H01Ω,where C is a constant which is h-independent.Lemma 11.
There exists an operatorPh:H1Ω⟶Vh2 which satisfies the following properties:
(56)bPhv,qh=bv,qh,∀qh∈Mh,(57)bΓPhv,qh=bΓv,qh,∀qh∈Mh,(58)∥Phv∥1≤C∥v∥1,C>0independentofh.Proof.
LetPh0:H1Ω⟶Vh be the usual L2-projection operator with
(59)∥Ph0v∥1≤c1∥v∥1,(60)∥v−Ph0v∥0≤c2h∥v∥1.
The first inequality (59) follows from Lemma (10) and the density of H01Ω in H1Ω H01Ω∩V¯∥·∥1,Ω⊂H1Ω, and the second (62) is an L2-error estimate. Now we fix a linear mapping Ph1:L2Ω⟶Vb such that
(61)∫kPh1v=∫kv,foreachK∈Th.
We may interpret the mapPh1 as a process with two steps. First, we apply the L2-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
(62)∥Ph1v∥0≤c3∥v∥0.
More precisely, we definePh1 as follows:
(63)Ph1v∣Kx,y=βKbKx,y,foreachK∈Th.
In view of condition (61), the constant βK is taken to be ∫kv/∫kbK. Observe that for ∫kbK=c4h2 and ∫kbK2=c5h2, where c4,c5>0 are constants,
(64)∥Ph1v∥02=∑K∈Th∫kPh1v2=∑K∈Th∫kβK2bK2=∑K∈Th∫kv∫kbK∫kbK2≤1c42h4∑K∈Th∫kv2∣K∣c5h2=c5c42∥v∥02.
The above inequality is due to Cauchy-Schwarz and∣K∣=h2 the area of the cell K. Let us set
(65)Phv=Ph0v+Ph1v−Ph0v.
Verify the two properties in the statement of the lemma. By the virtue of the construction ofPh and we impose that v−Ph0v=0 on ∂Ω
(66)∫Kv−Phv=∫KI−Ph1v−Ph0v=0∀K∈Th.
Becauseqh is piecewise linear, we apply Green’s theorem
(67)bv−Phv,qh=∫Ωdivv−Phvqh=∫∂Ωv−Phvnqh−∫Ωv−Phv∇qh=−∑K∈ThCK∫Ωv−Phv=0.
It is easy to prove (57) by applying Green’s theorem and using (56)
(68)∥Phv∥1≤∥Ph0v∥1+∥Ph1v−Ph0v∥1≤c1∥v∥1+c6h−1∥Ph1v−Ph0v∥0≤c1∥v∥1+c6h−1c3∥v−Ph0v∥0≤c1∥v∥1+c6h−1c3c2h∥v∥1≤c∥v∥1.
This proves the assertion (58) with c=c1+c6c3c2, such as c6 of Lemma (9).Theorem 12.
Let the space beBh=u∈Vh2\bPhu,qh=0,∀qh∈MhΩ such that Vh2=Bh⊕Bh⊥. The pair Vh2,Mh satisfies the discrete inf-sup condition:
(69)∃λ0>0,sup∥qh∥0=1bvh,qh∥vh∥Vh2≥λ0,∀vh∈Bh⊥,supvh∈Bh⊥bvh,qh>0,∀qh∈MhΩ.Proof.
LetPhv∈Bh⊥ be given. By the continuous form of the inf-sup condition (20) and (56) and (57) of Lemma 11, there exist c>0 and ρ>0 such that 1/c∥Phv∥Vh2≤∥v∥1,Ω≤1/ρsup∥qh∥=1bv,qh=1/ρsup∥qh∥=1bPhv,qh; then, we find sup∥qh∥=1bPhv,qh>c/ρ∥Phv∥Vh2,∀Phv∈Bh⊥. According to (56) and the second condition of (20), we have supvh∈Bh⊥bvh,qh=supv∈B⊥bv,qh>0,∀qh∈MhΩ.Remark 13.
With a similar demonstration as above, the bilinear formbΓ (16) on Vh2×Mh verifies the inf-sup condition cited in Theorem 12.


## 6.1. Linear Constant Element

More precisely, we giveVh and Mh, the displacement and div-displacement finite element spaces, respectively, as follows:
(49)Vh=u1,u2/uk=Σj∈Iuαjkϕj,αjk∈R,k=1,2,(50)Mh=ψ/ψ=Σi∈Iψβikφi,βi∈R.A sufficient condition forVh,Mh to satisfy the inf-sup condition is given in [21, 22]. A similar result has been proved in [26] for a different pair Vh,Mh.


## 6.2. Quadratic-Linear Element

Forj∈Iψ,
(51)φjx,y=1hφxh−j0φyh−j1,where φ is defined as
(52)φx≔xon0,1,2−xon1,2,2−xon2,3,x−4on3,4.For constructing the displacement approximation space, we do the following. Fori∈Iu, 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 ϕ
(53)ϕx≔x22on0,1,12+x−1−x−12on1,2,3−x22on2,3,Vb≔⊕K∈ThVK,where the quadrangulation Th 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 VK is the one dimensional subspace spanned by the function bK given by
(54)bK≔ωx,ybxh−k0byh−k1.The functionbx≔xx−1 is chosen so that the bubble function bK vanishes on the edges of the cell K. k0,k1 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δ−dx,0/δl, dx=distx,∂Ω, 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 beVh2≔Vh⊕Vb.This pairVh2,Mh of discretization spaces satisfies the discrete inf-sup condition.Lemma 7 (see [3]).
A web-basisBii∈I is a stable basis with respect to the L2-norm, C1∥aii∈I∥≤∥∑i∈IBi∥0≤C2∥aii∈I∥, where C1,C1 are constants.
We can also bound higher order Sobolev norms in terms of the 2-norm of the coefficients.Lemma 8 (see [3]).
A web-basisBii∈I satisfies ∥∑i∈IaiBi∥1≤C3h−1∥aii∈I∥, and C3>0 is a constant.Lemma 9 (inverse estimate).
LetVh be the finite element space considered as in (49). Then, there exists a constant C such that ∥∑i∈IaiBi∥1≤C6h−1∥∑i∈IaiBi∥0.
The proof follows from Lemmas7 and 8.Lemma 10 (see [21, 22]).
SupposeTh is a family of uniform quadrilateral mesh on the domain Ω and ∏h is the L2-projection operator onto Vh⊂H01Ω. Then
(55)∥∏hv∥1≤C∥v∥1,∀v∈H01Ω,where C is a constant which is h-independent.Lemma 11.
There exists an operatorPh:H1Ω⟶Vh2 which satisfies the following properties:
(56)bPhv,qh=bv,qh,∀qh∈Mh,(57)bΓPhv,qh=bΓv,qh,∀qh∈Mh,(58)∥Phv∥1≤C∥v∥1,C>0independentofh.Proof.
LetPh0:H1Ω⟶Vh be the usual L2-projection operator with
(59)∥Ph0v∥1≤c1∥v∥1,(60)∥v−Ph0v∥0≤c2h∥v∥1.
The first inequality (59) follows from Lemma (10) and the density of H01Ω in H1Ω H01Ω∩V¯∥·∥1,Ω⊂H1Ω, and the second (62) is an L2-error estimate. Now we fix a linear mapping Ph1:L2Ω⟶Vb such that
(61)∫kPh1v=∫kv,foreachK∈Th.
We may interpret the mapPh1 as a process with two steps. First, we apply the L2-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
(62)∥Ph1v∥0≤c3∥v∥0.
More precisely, we definePh1 as follows:
(63)Ph1v∣Kx,y=βKbKx,y,foreachK∈Th.
In view of condition (61), the constant βK is taken to be ∫kv/∫kbK. Observe that for ∫kbK=c4h2 and ∫kbK2=c5h2, where c4,c5>0 are constants,
(64)∥Ph1v∥02=∑K∈Th∫kPh1v2=∑K∈Th∫kβK2bK2=∑K∈Th∫kv∫kbK∫kbK2≤1c42h4∑K∈Th∫kv2∣K∣c5h2=c5c42∥v∥02.
The above inequality is due to Cauchy-Schwarz and∣K∣=h2 the area of the cell K. Let us set
(65)Phv=Ph0v+Ph1v−Ph0v.
Verify the two properties in the statement of the lemma. By the virtue of the construction ofPh and we impose that v−Ph0v=0 on ∂Ω
(66)∫Kv−Phv=∫KI−Ph1v−Ph0v=0∀K∈Th.
Becauseqh is piecewise linear, we apply Green’s theorem
(67)bv−Phv,qh=∫Ωdivv−Phvqh=∫∂Ωv−Phvnqh−∫Ωv−Phv∇qh=−∑K∈ThCK∫Ωv−Phv=0.
It is easy to prove (57) by applying Green’s theorem and using (56)
(68)∥Phv∥1≤∥Ph0v∥1+∥Ph1v−Ph0v∥1≤c1∥v∥1+c6h−1∥Ph1v−Ph0v∥0≤c1∥v∥1+c6h−1c3∥v−Ph0v∥0≤c1∥v∥1+c6h−1c3c2h∥v∥1≤c∥v∥1.
This proves the assertion (58) with c=c1+c6c3c2, such as c6 of Lemma (9).Theorem 12.
Let the space beBh=u∈Vh2\bPhu,qh=0,∀qh∈MhΩ such that Vh2=Bh⊕Bh⊥. The pair Vh2,Mh satisfies the discrete inf-sup condition:
(69)∃λ0>0,sup∥qh∥0=1bvh,qh∥vh∥Vh2≥λ0,∀vh∈Bh⊥,supvh∈Bh⊥bvh,qh>0,∀qh∈MhΩ.Proof.
LetPhv∈Bh⊥ be given. By the continuous form of the inf-sup condition (20) and (56) and (57) of Lemma 11, there exist c>0 and ρ>0 such that 1/c∥Phv∥Vh2≤∥v∥1,Ω≤1/ρsup∥qh∥=1bv,qh=1/ρsup∥qh∥=1bPhv,qh; then, we find sup∥qh∥=1bPhv,qh>c/ρ∥Phv∥Vh2,∀Phv∈Bh⊥. According to (56) and the second condition of (20), we have supvh∈Bh⊥bvh,qh=supv∈B⊥bv,qh>0,∀qh∈MhΩ.Remark 13.
With a similar demonstration as above, the bilinear formbΓ (16) on Vh2×Mh verifies the inf-sup condition cited in Theorem 12.


## 7. 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 asVh and Mh in general.) To find the discrete solution pair uh,ψh∈Vh×Mh, it is enough to find the coefficient vectors of displacement u¯≔u1¯,u2¯T, where uα¯=uiα¯i∈IuT, α=1,2, and div-displacement ψ¯≔ψj¯j∈IψT.Assuming thatB−1Aαβ=ξαβ, and B−1αβ=ζαβ for α,β=1,2. Let us define
(70)Ai,kα=∫Ωμ∇ϕi·∇ϕkdΩ+∫ΓξααϕiϕkdΓ,Bi,jα=∫Ωλ+μφj∂ϕi∂xαdΩ,BΓ,i,jα=∫Ωλ+μφj∂ϕi∂xαdΩ+∫ΓμφjnϕidΓ,Dm,n=∫Ωλ+μφmφndΩ,Ls1=∫Ωf1ϕsdΩ+∫Γζ11+ζ21g1ϕsdΓ,Ls2=∫Ωf2ϕsdΩ+∫Γζ22+ζ12g2ϕsdΓ.The matrix form of the discrete Navier-Lamé equations is written as(71)A10BΓ10A2BΓ2B1,TB2,T−Du1¯u2¯ψ¯=L1L20.By defining the assembled matricesA=DiagA1,A2, BΓ=BΓ1,BΓ2T B=B1,B2, F=L1,L2T, we rewrite the above matrix form as
(72)ABΓBT−Du¯ψ¯=F0.Decoupling of div-displacement and displacement: from the above matrix form, we write(73)Au¯+BΓψ¯=F,BTu¯−Dψ¯=0.Premultiplying the first equation above byA−1, followed by BT, and using the second equation, we get
(74)D+BTA−1BΓψ¯=BTA−1F.In the following, we show thatD+BTA−1BΓ is invertible, and hence, the above equation can be solved for ψ¯.The positive definiteness follows from the positiveD+BTA−1BΓ.Let us defineBΓ=B+bΓ, with bΓi,j=∫ΓμφjnϕidΓ,
(75)D+BTA−1BΓψ,ψRM=Dψ,ψRM+BTA−1BΓψ,ψRM=Dψ,ψRM+A−1B+bΓψ,BψRM=Dψ,ψRM+A−1Bψ,BψRM+A−1bΓψ,BψRM,where M=cardinality of Iψ.The positive definiteness ofBTA−1B 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).
(76)∥Bψ∥≥ρ∥ψ∥.Indeed, one has(77)A−1Bψ,BψRM≥β∥Bψ∥2≥βρ2∥ψ∥2.Condition (21) gives
(78)Dψ,ψRM≥ε∥ψ∥2.The inequality of Cauchy-Schwarz gives(79)A−1bΓψ,BψRM≥−∥A−1bΓψ∥∥Bψ∥.On the other hand,(80)∥A−1bΓψ∥=∑i∈IuA−1bΓψi1/2=∑i∈Iu∑m∈Iψ∑k∈IuAik−1bkmψm121/2≤∑i∈IucardIucardIψmaxk∈Iumaxm∈Iψ∣Aik−1bkm∣∑m∈Iψψm121/2≤cardIu2cardIψmaxk,i∈Iumaxm∈Iψ∣Aik−1bkm∣1/2∥ψ∥.By takingcardIψ=M, cardIu=N, and maxk,i∈Iumaxm∈Iψ∣Aik−1bkm∣=α. Then, we find that
(81)A−1bΓψ,BψRM≥−NMα1/2ρ∥ψ∥2.By combining (77), (78), and (81), we obtain
(82)D+BTA−1BΓψ,ψRM≥βρ2+ε−NρMα1/2∥ψ∥2,as soon as βρ2+ε−NρMα1/2>0.This proves the existence of the discrete solutionuh,ψh.


## 8. 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 Ω=x,y/x−1/22+y−1/22≤1/22\x,y/x−1/2+d2+y−1/22≤r2, with r=0.1, d=0.06, and the edge of Ω consists of two disjoint parts in the form of two circles ∂Ω=Γ=ΓD∪ΓN such as ΓD=x,y/x−1/2+d2+y−1/22−r2=0, ΓN=x−1/22+y−1/22−1/22=0.We give in (42) the Dirichlet matrix A=I2 on ΓD, 02 on ΓN, and we give the Neumann matrix B=I2 on ΓN, 02 on ΓD, the traction force g=0,0 on ∂Ω.Figures5 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.Figure 5
Numerical solution forH=8,16.Figure 6
Numerical solution forH=32,64.As mentioned in Table1, 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), errorpde=∥−μΔuh−λ+μ∇ψh−f∥0 and rate r=NaN−diffloge/log2. 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 H1Ω.Table 1
Different errors.
H
Error of pde
Rate
8
1.118e+000
NaN
16
3.644e−001
1.618e+000
32
1.831e−002
4.314e+000
64
7.003e−003
1.387e+000
128
1.775e−003
1.980e+000
256
4.462e−004
1.992e+000We show the relative residual for grid widths1/H=1/8,⋯,1/256. While the residuals decay, estimates of the rates exhibit a regular increasing behavior when H is 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–31].


## 9. Conclusion

Using the equations of linear elasticity as a model problem with the boundary conditionCA,B, we have described the implementation of finite element methods with weighted B-splines 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 order1/H, which is theoretically correct; these numerical results are attributed to the numerical stabilities of the solution. Moreover, the advantage of this problem with CA,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.



---

*Source: 4879723-2020-05-01.xml*