Cut Finite Element Methods for Linear Elasticity Problems

TL;DR

This paper introduces a cut finite element method for linear elasticity problems, incorporating stabilization techniques, error analysis, and applications to various structural and material models, including fiber-reinforced composites.

Contribution

The paper develops a novel cut finite element framework for linear elasticity with higher order elements, stabilization, and applications to complex material models.

Findings

Error estimates and condition number bounds established.

Numerical examples demonstrate method effectiveness.

Application to fiber-reinforced materials with stabilized beam modeling.

Abstract

We formulate a cut finite element method for linear elasticity based on higher order elements on a fixed background mesh. Key to the method is a stabilization term which provides control of the jumps in the derivatives of the finite element functions across faces in the vicinity of the boundary. We then develop the basic theoretical results including error estimates and estimates of the condition number of the mass and stiffness matrices. We apply the method to the standard displacement problem, the frequency response problem, and the eigenvalue problem. We present several numerical examples including studies of thin bending dominated structures relevant for engineering applications. Finally, we develop a cut finite element method for fibre reinforced materials where the fibres are modeled as a superposition of a truss and a Euler-Bernoulli beam. The beam model leads to a fourth order problem which we discretize using the restriction of the bulk finite element space to the fibre together with a continuous/discontinuous finite element formulation. Here the bulk material stabilizes the problem and it is not necessary to add additional stabilization terms.

Figures

Tables

Table 1. Table 1 : Numerical condition number estimates for stiffness matrix 𝐀 𝐀 \mathbf{A} and mass matrix 𝐌 𝐌 \mathbf{M} in the perfectly fitted situation in Figure  11 . The preconditioning here is a simple diagonal scaling of the matrix.

$𝐀$	$p$	Plain	Precond.		$𝐌$	$p$	Plain	Precond.
	1	$5.2187\cdot {10}^5$	$3.4096\cdot {10}^3$			1	$2.2392\cdot {10}^1$	$1.1646\cdot {10}^1$
	2	$5.0805\cdot {10}^6$	$2.7717\cdot {10}^4$			2	$6.1285\cdot {10}^1$	$1.3824\cdot {10}^1$
	3	$1.9148\cdot {10}^7$	$1.1609\cdot {10}^5$			3	$1.2169\cdot {10}^2$	$2.0935\cdot {10}^1$
	4	$6.0930\cdot {10}^7$	$4.6446\cdot {10}^5$			4	$3.1309\cdot {10}^2$	$3.4586\cdot {10}^1$
	5	$1.7828\cdot {10}^8$	$1.7832\cdot {10}^6$			5	$9.6345\cdot {10}^2$	$6.2514\cdot {10}^1$

Table 2. Table 2 : Numerical condition number estimates for the stiffness matrix 𝐀 𝐀 \mathbf{A} and the mass matrix 𝐌 𝐌 \mathbf{M} for the ‘worst case’ cut situation in Figure  11 where only 1/1000 of boundary elements are inside the domain. The preconditioning here is a simple diagonal scaling of the matrix.

$𝐀$	$p$	Plain	Precond.	Stabilized	Stab. + Precond.
	1	$1.6675\cdot {10}^{16}$	$2.8270\cdot {10}^3$	$4.6725\cdot {10}^6$	$2.9021\cdot {10}^3$
	2	$6.0692\cdot {10}^{30}$	$1.4264\cdot {10}^{16}$	$1.0239\cdot {10}^9$	$2.3952\cdot {10}^4$
	3	$1.6325\cdot {10}^{32}$	$4.0341\cdot {10}^{19}$	$1.3580\cdot {10}^{11}$	$4.7715\cdot {10}^5$
	4	$3.0219\cdot {10}^{32}$	$1.0862\cdot {10}^{20}$	$2.2118\cdot {10}^{13}$	$1.1224\cdot {10}^8$
	5	$4.1334\cdot {10}^{33}$	$2.0851\cdot {10}^{21}$	$4.4914\cdot {10}^{15}$	$3.5444\cdot {10}^{10}$
$𝐌$	$p$	Plain	Precond.	Stabilized	Stab. + Precond.
	1	$2.6963\cdot {10}^{19}$	$1.1848\cdot {10}^1$	$3.4553\cdot {10}^3$	$1.5677\cdot {10}^1$
	2	$3.8825\cdot {10}^{34}$	$5.1576\cdot {10}^{16}$	$1.8307\cdot {10}^5$	$1.3623\cdot {10}^3$
	3	$1.2091\cdot {10}^{35}$	$5.7858\cdot {10}^{18}$	$1.0203\cdot {10}^7$	$2.6414\cdot {10}^5$
	4	$3.6643\cdot {10}^{35}$	$6.0559\cdot {10}^{19}$	$1.4240\cdot {10}^9$	$6.7336\cdot {10}^7$
	5	$4.7194\cdot {10}^{35}$	$5.1345\cdot {10}^{19}$	$5.3499\cdot {10}^{11}$	$2.4344\cdot {10}^{10}$

Equations

$$-\mathbf{\sigma }(\mathbf{u})\cdot \nabla$$

$$\mathbf{\sigma }(\mathbf{u})\cdot \mathbf{n}$$

$$\mathbf{u}$$

$$\bm{\sigma}(\bm{u})=2\mu\bm{\epsilon}(\bm{u})+\lambda\text{tr}(\bm{\epsilon}(\bm{u})),\qquad\bm{\epsilon}(\bm{u})=\frac{1}{2}\Big{(}\bm{u}\otimes\nabla+\nabla\otimes\bm{u}\Big{)}$$

$$-\mathbf{\sigma }(\mathbf{u})\cdot \nabla -{\omega }^2\mathbf{u}$$

$$\mathbf{\sigma }(\mathbf{u})\cdot \mathbf{n}$$

$$\mathbf{u}$$

$$E(\omega )={\int }_{\Omega }\mathbf{\sigma }(\mathbf{u}(\omega )):\mathbf{ϵ}(\mathbf{u}(\omega ))$$

$$-\mathbf{\sigma }(\mathbf{u})\cdot \nabla -\lambda \mathbf{u}$$

$$\mathbf{\sigma }(\mathbf{u})\cdot \mathbf{n}$$

$$\mathbf{u}$$

$$m(\mathbf{v},\mathbf{w})$$

$$a(\mathbf{v},\mathbf{w})$$

$$a(\mathbf{u},\mathbf{v})=l(\mathbf{v}),\forall \mathbf{v}\in {\mathbf{V}}_0$$

$$l(\mathbf{v})$$

$$a(\mathbf{u},\mathbf{v})-{\omega }^{2}m(\mathbf{u},\mathbf{v})=l(\mathbf{v}),\forall \mathbf{v}\in {\mathbf{V}}_0$$

$$a(\mathbf{u},\mathbf{v})-\lambda m(\mathbf{u},\mathbf{v})=0,\forall \mathbf{v}\in {\mathbf{V}}_0$$

$$l_h(\mathbf{v})=0,\mathbf{v}\in \mathbf{RM}$$

$${\mathcal{K}}_h=\{K\in {\mathcal{K}}_{0,h}:K\cap \Omega \ne \varnothing \}$$

$$N_h(\Omega )={\cup }_{K\in {\mathcal{K}}_h}$$

$${\mathcal{K}}_h(\Gamma )$$

$${\mathcal{F}}_h(\Gamma )$$

$${\mathbf{V}}_h={\mathbf{V}}_{0,h}{|}_{N_h(\Omega )}$$

$$j_h(\mathbf{v},\mathbf{w})$$

$$m_h(\mathbf{v},\mathbf{w})$$

$$a_h(\mathbf{v},\mathbf{w})$$

$$A_h(\mathbf{v},\mathbf{w})$$

$$b_h(\mathbf{v},\mathbf{w})=2\mu (\mathbf{v},\mathbf{w}{)}_{\partial {\Omega }_D}+\lambda (\mathbf{v}\cdot \mathbf{n},\mathbf{w}\cdot \mathbf{n}{)}_{\partial {\Omega }_D}$$

$$j_{h,D}(\mathbf{v},\mathbf{w})$$

$$j_{h,N}(\mathbf{v},\mathbf{w})$$

$$a_h(v,w)=a(v,w)+{\gamma }_a\left(j_{h,N}(v,w)+h^{-2}{j}_{h,D}(v,w)\right)$$

$$A_h({\mathbf{u}}_h,\mathbf{v})=L_h(\mathbf{v}),\forall \mathbf{v}\in {\mathbf{V}}_h$$

$$L_h(\mathbf{v})$$

$$-({\mathbf{g}}_D,\mathbf{\sigma }(\mathbf{v})\cdot \mathbf{n}{)}_{\partial {\Omega }_D}+\beta h^{-1}{b}_h({\mathbf{g}}_D,\mathbf{v}{)}_{\partial {\Omega }_D}$$

$$A_h({\mathbf{u}}_h,\mathbf{v})-{\omega }^2{m}_h({\mathbf{u}}_h,\mathbf{v})=L_h(\mathbf{v}),\forall \mathbf{v}\in {\mathbf{V}}_h$$

$$A_h({\mathbf{u}}_h,\mathbf{v})-{\lambda }_h{m}_h({\mathbf{u}}_h,\mathbf{v})=0,\forall \mathbf{v}\in {\mathbf{V}}_h$$

Full text

1. 1 Introduction

2. 2 Linear Elasticity and the Cut Finite Element Method

3. 2.1 Linear Elasticity and Model Problems

4. 2.2 Weak Formulations

5. 2.3 The Mesh and Finite Element Spaces

6. 2.4 Bilinear Forms

7. 2.5 The Finite Element Method

8. 3 Implementational Aspects

9. 3.1 Boundary Representation

10. 3.2 Quadrature

11. 4 Theoretical Results

12. 4.1 Norms

13. 4.2 Inverse Estimates

14. 4.3 Properties of the Forms

15. 4.4 Conditioning of the Mass and Stiffness Matrices

16. 4.5 Interpolation

17. 4.6 A Priori Error Estimates

18. 5 Examples

19. 5.1 Benchmarks

20. 5.2 Condition Numbers

21. 5.3 Thin Geometries

22. 5.4 Compound Bodies

23. 5.5 Two-Grid Eigenvalue Estimation

24. 6 Modeling of Embedded Structures with Application to Fibre Reinforced Materials

25. 6.1 Modeling of Trusses and Beams

26. 6.2 Finite Element Discretization

27. 6.3 Numerical Examples

1 Introduction

The CutFEM Paradigm

In this contribution we develop a cut finite element method (CutFEM) for linear elasticity. The CutFEM paradigm is based on the following main concepts

• •

The computational domain is considered as a subset of polygonal domain which is equipped with a so called background mesh. The computational domain is represented on the background mesh and its boundary and interior interfaces are allowed to cut through the elements in an arbitrary fashion. On the background mesh there is a finite element space.

• •

The active mesh is the set of elements that intersects the computational domain. The finite element space on the active mesh is taken to be the restriction of the finite element space on the background mesh.

• •

The finite element method is based on a variational formulation on the computational domain, which includes weak enforcement of all interface and boundary conditions.

• •

Stabilization terms are added to handle the loss of control emanating from the presence of cut elements. The stabilization terms essentially restores the properties of standard finite element formulations.

The CutFEM is designed in such a way that we can establish the following theoretical results

• •

Optimal order a priori error estimates.

• •

Optimal order bounds on the condition number of the stiffness and mass matrices.

These results hold also for higher order polynomials provided the geometric representation of the domain and quadrature computations are sufficiently accurate. See the overview article **[3]** and the references therein for further details.

Contributions

In this contribution we focus on applications of CutFEM in solid mechanics. In particular we consider:

• •

Formulation of CutFEM for linear elasticity including an easy to access account of the theoretical results. We also provide extensive numerical results illustrating the performance of the method for test problems of practical interest. We study the displacement problem, the frequency response problem, and the eigenvalue problem. We also study the behavior of the method for thin structures which are common in engineering applications.

• •

Employing the CutFEM paradigm to the study of reinforced bulk materials, where the reinforcement consists of fibres, which are modeled as a superposition of beams and trusses. The beam model leads to a fourth order problem which we discretize using techniques from continuous/discontinuous Galerkin methods and the finite element space is the restriction of the bulk finite element space to the fibre. This approach requires a continuous finite element space consisting of at least second order piecewise polynomials. In its simplest form, we simply superimpose the additional stiffness, from the fibres to the bulk stiffness, by summing the bilinear forms for the fibres to the bulk material bilinear form and using the same finite element space. Note that here the the bulk material provides sufficient stability and no further stabilization is necessary. We provide illustrating numerical examples.

Previous Work

The possibility of discretizing partial differential equations on the cut part of elements was first proposed by Olshanskii, Reusken, and Grande **[20]**, and has been further developed by this group in, e.g., **[21, 12, 11]**. The work presented here builds on the CutFEM developments by the authors and their coworkers concerning stabilization of cut elements in **[1, 4, 5, 18, 16]** and concerning beam, plate, and membrane formulations in **[15, 6, 14]**. Related work has been presented by Fries et al. **[10, 9]** and by Rank and coworkers **[22, 23]**, see also **[24, 7]**.

Outline

In Section 2 we formulate the governing equations, the different model problems, and the finite element method, in Section 3 we present the theoretical results including error estimates and bounds for the condition number of the mass and stiffness matrix, in Section 4 we present extensive numerical results, in Section 5 we develop a cut finite element method for fibre reinforced materials.

2 Linear Elasticity and the Cut Finite Element Method

2.1 Linear Elasticity and Model Problems

Let $\Omega$ be a domain in $\mathbb{R}^{d}$, $d=2$ or $3$, with boundary $\partial\Omega={\partial\Omega_{D}}\cup{\partial\Omega_{N}}$, ${\partial\Omega_{D}}\cap{\partial\Omega_{N}}=\emptyset$, and exterior unit normal $\bm{n}$. We consider the following problems:

• •

(The Displacement Problem)  Find the displacement $\bm{u}:\Omega\rightarrow\mathbb{R}^{d}$ such that

[TABLE]

where the stress and strain tensors are defined by

[TABLE]

with Lamé parameters $\lambda$ and $\mu$, $\bm{f}$, $\bm{g}_{N}$, $\bm{g}_{D}$ are given data, and $\bm{a}\otimes\bm{b}$ is the tensor product of vectors $\bm{a}$ and $\bm{b}$ with elements $(\bm{a}\otimes\bm{b})_{ij}=a_{i}b_{j}$.

• •

(The Frequency Response Problem)  Given a frequency $\omega$ find the displacement $\bm{u}:\Omega\rightarrow\mathbb{R}^{d}$, such that

[TABLE]

Here we often compute the energy

[TABLE]

for a range of values of the parameter $\omega$.

• •

(The Eigenvalue Problem)  Find the displacement eigenvector $\bm{u}:\Omega\rightarrow\mathbb{R}^{d}$ and eigenvalue $\lambda\in\mathbb{R}$, such that

[TABLE]

This problem has an infinite sequence of solutions $(\bm{u}_{i},\lambda_{i})$ where $0\leq\lambda_{1}<\lambda_{2}<\dots$, such that $\lambda_{i}\rightarrow\infty$ and there are no accumulation points. There may be multiple eigenvalues, which is common in applications due to symmetries of the domain.

2.2 Weak Formulations

Let \bm{V}_{\bm{g}}=\{\bm{v}\in\bm{H}^{1}(\Omega):\text{\bm{v}=\bm{g}$on$\partial\Omega_{D}}\}, and define the forms

[TABLE]

Here and below we use the standard notation $H^{s}(\omega)$ for the $\mathbb{R}^{d}$ valued Sobolev space of order $s$ on the set $\omega$, with norm $\|\cdot\|_{H^{s}(\omega)}$, for $s=0$ we use the notation $L^{2}(\omega)$ and denote the scalar product by $(\cdot,\cdot)_{\omega}$, and the norm by $\|\cdot\|_{\omega}$.

The weak formulations of problems (2.1), (2.3), and (2.5) take the form:

• •

(The Displacement Problem)  Find $\bm{u}\in\bm{V}_{\bm{g}_{D}}$ such that

[TABLE]

where the linear form on the right hand side is defined by

[TABLE]

• •

(The Frequency Response Problem)  Given $\omega\in\mathbb{R}$ find $\bm{u}\in\bm{V}_{0}$ such that

[TABLE]

• •

(The Eigenvalue Problem)  Find $(\bm{u},\lambda)\in\bm{V}_{0}\times\mathbb{R}$ such that

[TABLE]

Let $\bm{RM}=\text{ker}(\bm{\epsilon})$ be the space of linearized rigid body motions, which has dimension $3$ in $\mathbb{R}^{2}$ and $6$ in $\mathbb{R}^{3}$, see **[2]** Section 11. We assume for simplicity that the Dirichlet boundary condition completely determines the linearized rigid body part of the solution, i.e., if $v\in\bm{RM}$ and $v=0$ on $\partial\Omega_{D}$ then $v=0$, then problem (2.8) is well posed. More general situations may be considered by working in a quotient space. For instance, an important case in practice is a free elastic body, $\partial\Omega_{D}=\emptyset$. We then seek the solution in the quotient space $H^{1}(\Omega)/\bm{RM}$ and problem (2.8) is well posed if the right hand side satisfies

[TABLE]

The frequency response problem (2.10) is well posed if $\omega$ is not an eigenvalue.

2.3 The Mesh and Finite Element Spaces

• •

Let $\Omega_{0}$ be polygonal domain such that $\Omega\subset\Omega_{0}\subset\mathbb{R}^{d}$ and let $\{\mathcal{K}_{0,h},h\in(0,h_{0}]\}$ be a family of quasiuniform partitions, with mesh parameter $h$, of $\Omega_{0}$ into shape regular elements $K$. We denote the set of faces in $\mathcal{K}_{h,0}$ by $\mathcal{F}_{h,0}$. We refer to $\mathcal{K}_{h,0}$ as the background mesh.

• •

Let

[TABLE]

be the submesh of $\mathcal{K}_{h,0}$ consisting of elements that intersect $\Omega$. We denote the set of faces in $\mathcal{K}_{h}$ by $\mathcal{F}_{h}$ and we refer to $\mathcal{K}_{h}$ as the active mesh. The domain of the active mesh is

[TABLE]

and thus $N_{h}(\Omega)$ covers $\Omega$ and we obviously have $\Omega\subset N_{h}(\Omega)$.

• •

Given some subset of the boundary $\Gamma\subset\partial\Omega$ let

[TABLE]

Thus $\mathcal{F}_{h}(\Gamma)$ is the set of interior faces belonging to elements in $\mathcal{K}_{h}$ that intersects $\Gamma$.

• •

Let $\bm{V}_{0,h}$ be the space of piecewise continuous $\mathbb{R}^{d}$ valued polynomials of order $p$ defined on $\mathcal{K}_{0,h}$. Let the finite element space on $N_{h}(\Omega)$ be defined by

[TABLE]

We illustrate the computational domain $\Omega$ and the background mesh $\mathcal{K}_{h,0}$ in Figure 1, the domain of active elements $N_{h}(\Omega)$ in Figure 2, and the sets of faces $\mathcal{K}_{h}(\Gamma_{N})$ and $\mathcal{K}_{h}(\Gamma_{D})$ in Figure 3.

2.4 Bilinear Forms

• •

Let $\mathcal{F}_{h}(\partial\Omega)$ in be the set of interior faces that belongs to an element $K$ such that $K\cap\partial\Omega\neq\emptyset$, and define the stabilization form

[TABLE]

where $D_{\bm{n}_{F}}^{l}$ is the $l$:th partial derivative in the direction of the normal $\bm{n}_{F}$ to the face $F\in\mathcal{F}_{h}(\partial\Omega)$ and $[\cdot]$ denotes the jump in a discontinuous function at $F$.

• •

Define the stabilized forms

[TABLE]

where $\gamma_{m},\gamma_{a}>0$ are parameters. We will see that the stabilized forms have the same general properties as the corresponding forms on standard (not cut) elements.

• •

Define the stabilized Nitsche form

[TABLE]

where $\beta>0$ is a parameter and

[TABLE]

Note that the stabilization is included in the form $a_{h}$, see (2.20).

Remark 2.1

We will see in the analysis that we may use an alternative form of the stabilization term in $a_{h}$ defined as follows. Define

[TABLE]

with the obvious notation, and let

[TABLE]

Thus we use a stronger stabilization only on the Dirichlet boundary, which is needed in the proof of the coercivity of $A_{h}$. The weaker control on the Neumann part of the boundary is sufficient to establish the bounds on the condition number. See Figure 3. The stabilization of the form $m_{h}$ is always done using $j_{h}$ as defined in (2.18). As a rule of thumb weaker stabilization yields more accurate numerical results and therefore (2.25) may be prefered in practice.

2.5 The Finite Element Method

With $m_{h}$ and $A_{h}$ defined in (2.19) and (2.21) we have the finite element methods

• •

(Static Load) Find $\bm{u}_{h}\in\bm{V}_{h}$ such that

[TABLE]

where the right hand side is given by

[TABLE]

• •

(Frequency Response) Given $\omega\in\mathbb{R}$ find $\bm{u}_{h}\in\bm{V}_{h}$ such that

[TABLE]

• •

(Eigenvalue Problem) Find $(\bm{u}_{h},\lambda_{h})\in\bm{V}_{h}\times\mathbb{R}$ such that

[TABLE]

3 Implementational Aspects

3.1 Boundary Representation

While it is common for fictitious domain methods to be based on implicit representations of the geometry, such as level-sets or signed distance functions, mechanical components typically are parametrically described using CAD. We therefore in the present work focus on parametric boundary representations.

For simplicity we in our implementation use a piecewise linear boundary representation albeit with a resolution which is independent of the finite element mesh. Thus, within a single element in the mesh we can represent geometric features of arbitrary complexity, see for example Figure 4.

3.2 Quadrature

To assemble our matrices we need to correctly integrate over the intersection between the computational domain $\Omega$ and each active element $K$. As the geometry of $\Omega\cap K$ is allowed to be arbitrarily complex we use the 2D quadrature rule described in **[16]** for accurate integration of both full and tensor product polynomials of higher order over this intersection. This rule assumes the domain of integration is described as a number of closed loops of piecewise polynomial curves and by using the divergence theorem the integral is first posed as a boundary integral. Using the fundamental theorem of calculus it is then reformulated as a sequence of one dimensional integrals, one in each spatial dimension, which are evaluated using appropriate Gauss quadrature rules chosen based on the polynomial structure of the integrand and the polynomial order of the boundary description. Example quadrature points for integration of a tensor product polynomial of order $2$ are displayed in Figure 4. Note that each boundary segment produces a number of quadrature points.

4 Theoretical Results

In this section we develop the basic theoretical results concerning stability, conditioning of the stiffness and mass matrices, and error estimates for the displacement problem (2.1a-2.1c). The main results are:

• •

The stabilized forms $m_{h}$ and $a_{h}$ enjoy the same coercivity and continuity properties with respect to the proper norms on $N_{h}(\Omega)$ as corresponding standard forms on $\Omega$ equipped with a fitted mesh.

• •

Optimal order approximation holds in the relevant norms since there is a stable interpolation operator obtained by composing a standard Scott-Zhang interpolation operator with an extension operator that in a $H^{s}$ stable way extends a function from $\Omega$ to a neighborhood of $\Omega$ containing $N_{h}(\Omega)$.

• •

Optimal order error estimates follows from the stability and continuity properties combined with the approximation properties.

• •

The mass and stiffness matrices have condition numbers that scale with the mesh size in the same way as for standard fitted meshes.

4.1 Norms

• •

Define the norms

[TABLE]

We recall that $m_{h}$ and $a_{h}$ are the stabilized forms, see (2.19) and (2.20), so that

[TABLE]

• •

Define the discontinuous Galerkin norm associated with the form $A_{h}$,

[TABLE]

4.2 Inverse Estimates

The form $j_{h}$ provides additional control of the variation of functions in $\bm{V}_{h}$ close to the boundary. This additional control implies that certain inverse inequalities that hold for standard elements also hold for the stabilized forms. We have the following results.

Lemma 4.1

The following inverse estimates hold

[TABLE]

Proof.* (4.4). We first consider two neighboring elements $K_{1}$ and $K_{2}$ sharing a face $F$. We shall prove that*

[TABLE]

Iterating (4.7) a uniformly bounded number of times we reach an element in the interior of $\Omega$ and thus we can control the elements on the boundary in terms of elements in the interior and the stabilization term. To prove (4.7) we consider $\bm{v}\in V_{h}$ and let $\bm{v}_{i}$ be the polynomial on $K_{1}\cup K_{2}$ such that $\bm{v}_{i}=v|_{K_{i}}$, $i=1,2$. We then note that

[TABLE]

where we used the inverse inequality

[TABLE]

which holds for all $\bm{w}\in\mathcal{P}_{p}(K_{1}\cup K_{2})$ since $K_{1}$ and $K_{2}$ are shape regular.

Now for each $\bm{x}\in K_{1}$ we have the identity

[TABLE]

where $\bm{x}=\bm{x}_{F}+x_{n}\bm{n}_{F}$, with $\bm{n}_{F}$ the unit normal to $F$ exterior to $K_{2}$, and $x_{n}\in\mathbb{R}$ the distance from $\bm{x}$ to the hyperplane containing $F$. Let $P_{F}(K_{1})=\{\bm{x}_{F}=\bm{x}-x_{n}\bm{n}_{F}:\bm{x}\in K_{1}\}$ be the projection of $K_{1}$ onto the hyperplane containing $F$ and $h_{n}\lesssim h$ be the maximal value of $x_{n}$ over $K_{1}$. Then we have the estimate

[TABLE]

where we used the inverse estimate

[TABLE]

*which holds with a uniform constant due to shaperegularity and the fact that $[D^{j}_{n}\bm{v}]$ is a polynomial on $P_{F}(K_{1})$.

---

(4.5).* Here we have*

[TABLE]

where we used the inverse inequality

[TABLE]

for all $w\in\mathcal{P}_{p}(K_{1}\cup K_{2})$, which holds since $\nabla w\in\mathcal{P}_{p-1}(K_{1}\cup K_{2})$ is also a polynomial. Computing the gradient of (4.10) we have the identity

[TABLE]

where $\nabla_{F}=(1-\bm{n}_{F}\otimes\bm{n}_{F})\nabla$ is the gradient tangent to $F$. We then have the estimate

[TABLE]

*where we used an inverse estimate to remove the tangent gradient in the first term on the right hand side and finally we used the inverse estimate (4.14).

---

(4.6).* Proceeding as above we have the estimate*

[TABLE]

where we used the inverse inequalities

[TABLE]

for all $w\in\mathcal{P}_{k}(K_{1}\cup K_{2})$, which hold since $\bm{\epsilon}(w)\in\mathcal{P}_{p-1}(K_{1}\cup K_{2})$ and $\text{tr}\bm{\epsilon}(w)\in\mathcal{P}_{p-1}(K_{1}\cup K_{2})$ are both polynomials. Next, in view of the identity

[TABLE]

*we may prove (4.6) using the same approach as for the gradient estimate (4.5).

---

**

4.3 Properties of the Forms

Properties of $m_{h}$

Lemma 4.2

The bilinear form $m_{h}$, defined by (2.19), is continuous and coercive

[TABLE]

Proof.* Coercivity follows directly from the inverse estimate (4.4) and to show continuity we use an inverse estimate to conclude that*

[TABLE]

---

**

Properties of $A_{h}$

First we note that we have the following Poincaré inequality which shows that $|\kern-1.07639pt|\kern-1.07639pt|\cdot|\kern-1.07639pt|\kern-1.07639pt|_{h}$ is indeed a norm on $\bm{V}_{h}$ and plays an important role in the estimate of the condition number of the stiffness matrix associated with $A_{h}$.

Lemma 4.3

The following Poincaré estimate holds

[TABLE]

Next we turn to continuity and coercivity of $A_{h}$.

Lemma 4.4

The bilinear form $A_{h}$, defined by (2.20), is continuous

[TABLE]

and, if the stabilization parameter $\beta$ is large enough, coercive

[TABLE]

Proof.* We show continuity of $A_{h}$ on $\bm{V}+\bm{V}_{h}$, directly using the Cauchy-Schwarz inequality. To show coercivity we first note that using the Cauchy-Schwarz inequality*

[TABLE]

for each $\delta>0$. Using the elementwise inverse inequality

[TABLE]

where the hidden constant is independent of the position of the boundary $\partial\Omega_{D}$, see **[13]**, followed by the inverse estimate (4.5) we conclude that

[TABLE]

This is the crucial inverse inequality required in the standard proof, see **[17*]** Section 14.2, of coercivity of $A_{h}$ on $\bm{V}_{h}$, which holds if $\beta$ is large enough. Note that the size of the penalty parameter is completely independent of the actual position of the domain $\Omega$ in the background mesh.

---

**

Properties of $L_{h}$

The form $L_{h}$ is continuous on $\bm{V}_{h}$,

[TABLE]

Remark 4.1

For fixed $h\in(0,h_{0}]$ we may then use the coercivity and continuity of $A_{h}$, see Lemma 4.4, the continuity (4.36), and apply the Lax-Milgram lemma to conclude that there exists a unique solution to the finite element problem (2.26).

4.4 Conditioning of the Mass and Stiffness Matrices

Let $\{\bm{\varphi}_{j}\}_{j=1}^{N}$ be the standard Lagrange basis in $\bm{V}_{h}$ and $N$ the dimension of $\bm{V}_{h}$. Let $\bm{v}=\sum_{j=1}^{N}\widehat{v}_{j}\bm{\varphi}_{j}$ be the expansion of $\bm{v}$ in the Lagrange basis. We define the mass and stiffness matrices, ${\bf M}$ and ${\bf A}$ with elements

[TABLE]

which is equivalent to the identities

[TABLE]

where $(\cdot,\cdot)_{\mathbb{R}^{N}}$ is the inner product in $\mathbb{R}^{N}$.

Recall that the condition number of an $N\times N$ matrix $\bf B$ is defined by $\kappa({\bf B})=\|{\bf B}\|_{\mathbb{R}^{N}}\|{\bf B}^{-1}\|_{\mathbb{R}^{N}}$, where $\|\cdot\|_{\mathbb{R}^{N}}$ is the Euclidian norm in $\mathbb{R}^{N}$. We then have the following result.

Theorem 4.1

The condition numbers of mass and stiffness matrices, ${\bf M}$ and ${\bf A}$, defined by (4.37), satisfy the estimates

[TABLE]

The proof uses the approach in **[8]** and builds on the standard equivalence of norms

[TABLE]

where $\widehat{\bm{v}}\in\mathbb{R}^{N}$ is the coefficient vector with elements $\widehat{v}_{i}$.

Proof.* We present the proof for the bound on the condition number of the stiffness matrix. The bound for the condition number of the mass matrix follows in the same way but is slightly simpler since only the equivalence (4.40) and the coercivity and continuity properties in Lemma 4.2 are used, while in the the case of the stiffness matrix we also employ the Poincaré inequality and the inverse bound*

[TABLE]

We have

[TABLE]

Starting with the estimate of $\|{\bf A}\|_{\mathbb{R}^{N}}$ we obtain, using the definition of the stiffness matrix and the inverse bound (4.41),

[TABLE]

We conclude that $\|{\bf A}\widehat{\bm{v}}\|_{\mathbb{R}^{N}}\lesssim h^{d-2}\|\widehat{\bm{v}}\|_{\mathbb{R}^{N}}$, and thus

[TABLE]

Next to estimate $\|{\bf A}^{-1}\|_{\mathbb{R}^{N}}$ we note that

[TABLE]

where in the last step we used th eequivalence (4.40) and the Poincaré estimate (4.28) to conclude that

[TABLE]

Thus we arrive at

[TABLE]

and setting $\widehat{\bm{v}}={\bf A}^{-1}\widehat{\bm{w}}$ we find that

[TABLE]

and thus

[TABLE]

*Combining (4.45) and (4.51) the estimate of the condition number for the stiffness matrix follows.

---

**

4.5 Interpolation

To define the interpolation operator $\bm{\pi}_{h}$ we use a bounded extension operator $E:H^{s}(\Omega)\rightarrow H^{s}(U(\Omega))$, where $U(\Omega)$ is a neighborhood of $\Omega$ such that $N_{h}(\Omega)\subset U(\Omega)$, $h\in(0,h_{0}]$, and then using the Scott-Zhang interpolation operator $\bm{\pi}_{h,SZ}:\bm{H}^{1}(N_{h}(\Omega))\rightarrow\bm{V}_{h}$, we may define

[TABLE]

Using the simplified notation $\bm{u}=E\bm{u}$ on $U(\Omega)$ we have have the interpolation error estimate

[TABLE]

Summing over the elements $K\in\mathcal{K}_{h}$ and using the stability of the extension operator we obtain

[TABLE]

Furthermore, we have the following estimate for the interpolation error in the energy norm

[TABLE]

To verify (4.54) we note that using the trace inequality

[TABLE]

where the hidden constant is independent of the position of the boundary $\partial\Omega_{D}$ in element $K$ see **[13]** for further details, we may estimate the boundary terms as follows

[TABLE]

[TABLE]

where $\mathcal{K}_{h}(\partial\Omega_{D})=\{K\in\mathcal{K}_{h}:K\cap\partial\Omega_{D}\neq\emptyset\}$. Using the standard trace inequality

[TABLE]

where $F$ is a face associated with element $K$ we obtain

[TABLE]

where $\mathcal{K}_{h}(\mathcal{F}_{h}(\partial\Omega))$ is the set of elements that have a face in $\mathcal{F}_{h}(\partial\Omega)$, i.e. an element with a stabilized face. Thus we conclude that we have the estimate

[TABLE]

Finally, setting $\bm{v}=\bm{u}-\bm{\pi}_{h}\bm{u}$ and using (4.53) and the stability of the extension operator (4.55) follows.

4.6 A Priori Error Estimates

Theorem 4.2

Let $\bm{u}$ be the solution to problem (2.1a-2.1c) and $\bm{u}_{h}$ the corresponding finite element approximation defined by (2.26), then the following a priori error estimates hold

[TABLE]

Proof.* Using the coercivity and continuity properties of $A_{h}$, Galerkin orthogonality, and the interpolation estimate (4.55) we prove the error estimates using the standard approach, see [17].*

Adding and subtracting an interpolant and using the triangle inequality we have

[TABLE]

For the second term we use coercivity of $A_{h}$,

[TABLE]

and thus

[TABLE]

The $L^{2}$ estimate is proved using a standard duality argument, see **[17*]**.

---

**

5 Examples

The numerical results presented below include benchmarks of CutFEM for standard problems in linear elasticity, a study of the condition numbers of mass and stiffness matrices, and also various examples which demonstrate different modeling possibilities which are naturally combined with CutFEM.

Unless otherwise stated we in the numerical results below assume a linear elastic isotropic material with the material constants of steel, i.e., Young’s modulus $E=200\cdot 10^{9}$ Pa, Poisson’s ratio $\nu=0.3$, and density $\rho=7850$ kg/m3*. Default parameter values used in the CutFEM method are $\gamma_{D}=1000\cdot p^{2}$ for the Nitche penalty parameter, $\gamma_{m}=\rho\cdot 10^{-4}$ for the mass matrix stabilization and $\gamma_{a}=(2\mu+\lambda)\cdot 10^{-4}$ for the stiffness matrix stabilization. The finite elements used in are Lagrange elements with evenly distributed nodes; either full polynomials of order $p$ on triangles or tensor product polynomials of order $p$ on quadrilaterals.*

5.1 Benchmarks

We begin our numerical examples by studying the performance of CutFEM in the three standard problems of linear elasticity described in Section 2.1.

• •

(The Displacement Problem)  In the stationary load problem we benchmark CutFEM using a manufactured problem on the unit square. The geometry and the solution is given by

[TABLE]

and from this we deduce expressions for the input data $f$, $g_{N}$ and $g_{D}$. Visually inspecting numerical solutions using $p=\{1,2\}$ elements and different rotations of the background grid we in Figure 5 note that $p=1$ elements are sensitive to the grid rotation with regards to the quality of the stresses while $p=2$ elements give results which are visually invariant to the grid rotation. In Figure 6 we present convergence results in $L^{2}(\Omega)$ norm and we see that the performance of CutFEM using $p=\{1,2,3,4\}$ elements in this problem is equivalent to that of conforming FEM using the same elements.

• •

(The Frequency Response Problem)  To evaluate the performance of CutFEM in the frequency response problem we consider the cantilever beam with holes used to illustrate the cut meshes in Figures 1–3. This steel beam is fixated along left side and is under the influence of an oscillating gravitational load with an oscillatory frequency $\omega$. As reference we use $p=4$ elements on a conforming triangle mesh and we compare this to CutFEM calculations using $p=2$ elements on two structured meshes, see Figure 7. We are interested in accurately estimating the energy $E(\omega)$ defined in (2.4) and in Figure 8 we present the results for CutFEM compared to the higher order conforming method. It seems the coarse grid CutFEM estimation fails to capture some of the higher frequency details while the higher resolution CutFEM estimation succeeds. This is quite natural as the oscillatory modes corresponding to higher frequency loads most likely have more local details.

• •

(The Eigenvalue Problem)  We benchmark CutFEM for the eigenvalue problem on a free steel beam of length 3 m and height 0.3 m, i.e.

[TABLE]

In Figure 9 we give a visual comparison of an eigenmode computed using $p=2$ elements with conforming FEM and CutFEM in two different cut situations and we note no noticeable differences. To investigate the convergence in eigenvalues we as reference use conforming FEM with $p=4$ elements on a fine grid. The convergence results for an eigenvalue is presented in Figure 10 and we see that CutFEM in this problem performs equivalently to conforming FEM.

5.2 Condition Numbers

As the domain is allowed to cut the elements in the background grid in an arbitrary fashion we may end up with cut situations in which the method becomes ill-conditioned unless properly stabilized or preconditioned. While we focus our work on the stabilization side an alternative, or complement, to stabilization is preconditioning techniques which can yield good results, see for example **[7]**. In the numerical results below we include preconditioning in the form of simple diagonal scaling, demonstrated for the stationary load problem as follows. The stationary load problem in matrix form reads

[TABLE]

where $\mathbf{A}$ is the stiffness matrix, $\mathbf{U}$ is the vector of unknowns, and $\mathbf{L}$ is the load vector. By the substitution $\mathbf{U}=\mathbf{D}\mathbf{Y}$ where $\mathbf{D}$ is the diagonal matrix given by $\mathbf{D}_{ii}=\mathbf{A}_{ii}^{-1/2}$ we can instead consider the equivalent matrix problem

[TABLE]

where the (proconditioned) matrix $\mathbf{D}^{T}\mathbf{A}\mathbf{D}$ have diagonal elements which are all one and in which the symmetry of $\mathbf{A}$ is preserved.

We investigate the effect of preconditioning and stabilization using the model static load problem from Section 5.1. First we construct two extreme cut situations as illustrated in Figure 11; a perfectly fitted (conforming) mesh and a mesh where all boundary elements have a very small intersection with the domain. Comparing the condition number results for a perfectly fitted mesh in Table 1 to the ‘worst case’ mesh results in Table 2 we note that a combination of stabilization and preconditioning by simple diagonal scaling gives the best results. However, we also see that even when applying both stabilization and preconditioning condition numbers when using $p\geq 3$ elements are orders of magnitude greater in the worst case scenario. We see the same behavior when investigating how the condition numbers scale with the mesh size $h$. By rotating the background grid $\theta=\{0,\frac{\pi}{9},\frac{\pi}{7},\frac{\pi}{5}\}$ radians we construct both conforming meshes and various cut situations for which we compare condition numbers from preconditioning by diagonal scaling alone to condition numbers when we also add stabilization, i.e. CutFEM. In Figure 12 we for $p=\{1,2\}$ elements see good results, both regarding scaling and size of condition numbers. Note that for $p=1$ elements preconditioning by diagonal scaling alone is sufficient. For $p=\{3,4,5\}$ elements we in Figure 13 see that the condition numbers scale in the correct manner but their size in cut situations are orders of magnitude higher compared to those of conforming FEM. Reviewing the analytical results above we attribute this effect to the constant in the inverse estimates (4.9), (4.16) and (4.23) which can become quite large for higher order polynomials.

5.3 Thin Geometries

It is well known that low order elements ($p=1$) in elasticity problems suffer from locking on thin geometries when the number of elements in the thickness direction approaches one. This effect stems from the boundary condition $\bm{\sigma}(\bm{u})\cdot\bm{n}=0$ which in such situations constrains the non-zero components of the stress to the tangential plane. In CutFEM the situation may become even more extreme as the number of elements in the thickness direction can be a small fraction of an element. This effect is illustrated in Figure 14 where a thin cantilever beam under gravitational load using large elements exhibits extreme locking when using $p=1$ elements while there is no such effects when using higher order elements.

Furthermore, in CutFEM it is possible for the geometry to be curved inside an element, yielding a tangential plane with non-zero curvature. To avoid locking in such situations when using coarse meshes we must use even higher order elements or increase resolution until the curvature is small compared to the element size. As an illustration of this we here consider a free ring under a centrifugal load. In Figure 15 we note that while $p=2$ elements yield fairly good results for this mesh resolution we need $p=3$ elements for visually perfect results with regards to rotational symmetry for the stress distribution.

5.4 Compound Bodies

Compound bodies consist of several subdomains where the discrete solution on each subdomain is described on a separate mesh. Such a construction gives a number of powerful modeling possibilities, for example

• •

adding or changing small geometric details, for example fillets, without needing to reassemble all matrices

• •

in each subdomain using different material properties, mesh resolution or elements

• •

constructing objects from a number of predefined and preassembled building blocks

The weak enforcement of Dirichlet boundary conditions by Nitsche’s method **[19]** in CutFEM makes it easy to construct compound bodies by a simple adaptation of Nitsche’s method to interface conditions.

Consider the schematic compound body illustrated in Figure 16 for which we have the following interface conditions on $\Gamma_{12}$

[TABLE]

Summing the standard CutFEM formulations for all subdomains, with discrete functions on each subdomain described on its own background grid (and possibly its own type of elements), we for this example have the following interface term

[TABLE]

From this term and the interface conditions (5.5) we derive Nitsche’s method for the interface conditions via the calculation

[TABLE]

As an illustration of compound bodies we in Figures 17 and 18 consider an L-shape domain where the inside corner has been drilled to avoid a stress singularity in the solution. We have constructed this domain out of three subdomains described on three different background grids and we use $p=2$ elements to represent the solution on each mesh. Note that there is no need for the various meshes to match at the interface as the interface conditions are enforced weakly.

5.5 Two-Grid Eigenvalue Estimation

The use of structured background grids in CutFEM is very convenient when employing multi-grid techniques. We here illustrate this with a two-grid method for estimation of eigenvalues **[25]**. Consider two background grids, one coarse grid with mesh size $H$ and one fine with mesh size $h$, constructed such that the fine grid is a refinement of the coarse grid. The two-grid algorithm reads

Solve the eigenvalue problem on the coarse background grid: Find $(\bm{u}_{H},\lambda_{H})\in V_{H}\times\mathbb{R}$ such that

[TABLE] 2. 2.

Solve a single linear problem on the fine background grid: Find $\bm{u}_{h}\in V_{h}$

[TABLE] 3. 3.

Compute the Rayleigh quotient

[TABLE]

The analysis in **[25]** for a conforming two-grid method implies that the optimal choice of mesh size $H$ for the coarse grid is given by the relationship

[TABLE]

However in a cut situation, especially in light of the results in Section 5.3, it is very reasonable to assume that the coarse grid mesh size must be small enough to give an acceptable approximation of the eigenvalues and eigenmodes for this procedure to be efficient.

As a simple numerical example we consider the eigenvalue problem of a free ring with coarse and fine grids illustrated in Figure 19. As this is a free eigenvalue problem we in step 2 of the algorithm must seek the solution $\bm{u}_{h}\in\bm{V}_{h}/\bm{RM}$. The resulting approximations of the 10:th eigenvalue and the corresponding eigenmode are presented in Figure 20 together with a reference solution computed using high order parametric conforming FEM.

6 Modeling of Embedded Structures with Application to Fibre Reinforced Materials

In this Section we consider the embedding of membranes, thin plates, truss and beam elements into a higher dimensional elastic body. Simulations of this type are suitable for reinforced structures, e.g., reinforced concrete, fiber reinforced polymers, and composite materials in general.

The approach is based on:

• •

Given a continuous finite element space, based on at least second second order polynomials, we define the finite element space for the thin structure as the restriction of the bulk finite element space to the thin structure which is geometrically modeled by an embedded curve or surface.

• •

To formulate a finite element method on the restricted or trace finite element space we employ continuous/discontinuous Galerkin approximations of plate and beam models.

The thin structures are then modelled using the CutFEM paradigm and the stiffness of the embedded structure is in the most basic version, which we consider here, simply added to the bulk stiffness. This in turn means that the bulk stabilizes the CutFEM contributions and no additional stabilization is needed. In this initial report, we limit ourselves to two-dimensional bulk problems and thus the embedded structures will be limited to trusses and beams. In future work, the handling of Kirchhoff plates and beams and trusses of codimension 2 will be considered. The bulk problem may also be viewed as an interface problem on the case of an embedded surface in three dimensions or curve in two dimensions in order to more accurately approximate the bulk problem in the vicinity of the thin structure.

The work presented here is and extension of earlier work **[6]** where only membrane structures were considered, in which case a linear approximation in the bulk suffices.

6.1 Modeling of Trusses and Beams

Consider the elasticity equations (2.1) in $\Omega\subset\mathbb{R}^{2}$. Embedded in $\Omega$ we have straight truss and beam elements. Curved elements can be modeled using piecewise affine approximations of the curve, following **[14]**.

We consider modeling the truss and beam using tangential differential calculus. To this end we follow the exposition in **[15]** and then make some simplifying assumptions. Let then $\Sigma$ denote line embedded in $\mathbb{R}^{2}$, with tangent vector $\bm{t}$. We let $\bm{p}:\mathbb{R}^{2}\rightarrow\Sigma$ be the closest point mapping, i.e. $\bm{p}(\bm{x})=\bm{y}$ where $\bm{y}\in\Sigma$ minimizes the Euclidean norm $|\bm{x}-\bm{y}|_{\mathbb{R}^{3}}$. We define $\zeta$ as the signed distance function $\zeta(\bm{x}):=\pm|\bm{x}-\bm{p}(\bm{x})|$, positive on one side of $\Sigma$ and negative on the other. The line $\Sigma$ is assumed to be the center line of a truss element with thickness $t$, which we for simplicity assume is constant.

The linear projector $\bm{P}_{\Sigma}=\bm{P}_{\Sigma}(\bm{x})$, onto the tangent line of $\Sigma$ at $\bm{x}\in\Sigma$, is given by

[TABLE]

Based on the assumption that planar cross sections orthogonal to the midline remain plane after deformation we assume that the displacement takes the form

[TABLE]

where $\bm{u}_{0}:\Sigma\rightarrow\mathbb{R}^{2}$ is the deformation of the midline, decomposed as

[TABLE]

where $\bm{n}\perp\bm{t}$, and $\theta:\Sigma\rightarrow\mathbb{R}$ is an angle representing an infinitesimal rotation. Both $\bm{u}_{0}$ and $\theta$ are assumed constant in the normal plane.

In Euler–Bernoulli beam theory the beam cross-section is assumed plane and orthogonal to the beam midline after deformation and no shear deformations occur. This means that we have

[TABLE]

This definition for $\theta$ in combination with (6.2) constitutes the Euler–Bernoulli kinematic assumption.

We assume the usual Hooke’s law for one dimensional structural members

[TABLE]

where $\bm{\epsilon}_{\Sigma}(\bm{u}):=\bm{P}_{\Sigma}\bm{\epsilon}(\bm{u})\bm{P}_{\Sigma}$. We then notice that the strain energy density can be written

[TABLE]

Assuming computations are being done per unit length in the third dimension, the second moment of inertia $I={t^{3}}/{12}$ and the area $A=t$ and the total energy of the truss and beam structures are postulated as

[TABLE]

where the truss energy

[TABLE]

and the beam energy

[TABLE]

where $\partial^{2}_{tt}:=\partial_{t}\partial_{t}$. Assuming zero displacements and rotations at the end points of $\Sigma$, we thus seek $(u_{t},u_{n})\in H^{1}_{0}(\Sigma)\times H^{2}_{0}(\Sigma)$ such that

[TABLE]

for all $(v_{t},v_{n})\in H^{1}_{0}(\Sigma)\times H^{2}_{0}(\Sigma)$.

6.2 Finite Element Discretization

In analogy with the preceding discussion, we now let $\mathcal{K}_{h}$ be a quasiuniform partition, with mesh parameter $h$, of the bulk domain $\Omega$ into shape regular elements $K$, and denote the set of faces in $\mathcal{K}_{h}$ by $\mathcal{F}_{h}$. The elements of the bulk mesh cut by $\Sigma$, and the corresponding element faces, are denoted

[TABLE]

The intersection points between $\Sigma$ and element faces in $\mathcal{K}_{h}(\Sigma)$ is denoted

[TABLE]

and we assume that this is a discrete set of points (thus excluding the case where any $F\in\mathcal{F}_{h}$ coincides with a part of $\Sigma$ of nonzero measure). Setting as above

[TABLE]

and introducing

[TABLE]

and

[TABLE]

and defining $a_{\text{tot}}(\bm{v},\bm{w}):=a(\bm{v},\bm{w})+b(\bm{v},\bm{w})+c(\bm{v},\bm{w})$ we propose the following continuous/discontinuous Galerkin method. With

[TABLE]

we seek $\bm{u}\in V_{h}$ such that

[TABLE]

for all $\bm{v}\in V_{h}$. Here

[TABLE]

The terms on the discrete set $\mathcal{P}_{h}(\Sigma)$ are associated with the work of the end moments on the end rotation which occur due to the lack of $C^{1}(\Omega)$ continuity of the approximation, cf. **[14]**. In this formulation, stabilization terms of the kind discussed above can be added to ensure coercivity. Due to the presence of the bulk equations, the need for these terms is however mitigiated.

6.3 Numerical Examples

We consider a bulk material on the domain $\Omega=(0,4)\times(0,1)$ with $E=300$, $\nu=1/3$ in plane strain. The domain is fixed at the left end and unloaded on the rest of the boundary. It is first loaded with a volume load $\bm{f}=(0,-1)$ and the displacements using only the bulk material are shown in Fig. 21 with corresponding (Frobenius) norm of the stresses shown in Fig. 22. We next introduce two truss elements as shown in Fig 23 at $y=0.249$ and at $y=0.751$. The thickness of the trusses was set to $t=0.1$ and Young’s modulus $E=10^{4}$. In Figs. 24 and 25 we show the corresponding deformations and norm of the stresses. Finally we remove the trusses and add a beam at $y=0.501$, with thickness $t=0.1$ and $E=10^{6}$. In Figs. 26 and 27 we show the corresponding deformations and norm of the stresses.

In our second example we use the same boundary conditions but pull the material using $\bm{f}=(100,0)$. The corresponding displacements using only the bulk material are shown in Fig. 28 with associated norm of the stresses shown in Fig. 29. We then insert the truss elements from the previous example and apply the load again to obtain the deformation shown in Fig. 30 with associated norm of stresses in Fig. 31.

In all these examples we see that the added lower dimensional structural elements have a profound impact on the deformation and stress distribution in the bulk material. The possibility of modeling lower dimensional structural elements without meshing is clearly beneficial and can be used for example in optimization algorithms where the structural members are required to move around to find their optimal shape and place.