Virtual Element Methods F. Brezzi IMATI-C.N.R., Pavia, Italy - - PowerPoint PPT Presentation

▶
virtual element methods f brezzi
SMART_READER_LITE
LIVE PREVIEW

Virtual Element Methods F. Brezzi IMATI-C.N.R., Pavia, Italy - - PowerPoint PPT Presentation

Apr 09, 2023 228 likes •976 views

Virtual Element Methods F. Brezzi IMATI-C.N.R., Pavia, Italy IUSS-Istituto Universitario di Studi Superiori, Pavia, Italy KAU-King Abdulaziz University, Jeddah, Saudi Arabia From joint works with: L. Beira o da Veiga, A. Cangiani G. Manzini,

slide-1
SLIDE 1

Virtual Element Methods

  • F. Brezzi

IMATI-C.N.R., Pavia, Italy IUSS-Istituto Universitario di Studi Superiori, Pavia, Italy KAU-King Abdulaziz University, Jeddah, Saudi Arabia From joint works with: L. Beira˜

  • da Veiga, A. Cangiani G. Manzini,

L.D. Marini and A. Russo Milano, 17-19 September, 2012 Workshop on Polygonal and Polyhedral Meshes

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 1 / 71

slide-2
SLIDE 2

Outline

1 Traditional FEM 2 Basic Virtual Elements 3 A more precise version 4 Robustness 5 Higher order VEM 6 Linear elasticity 7 Nearly incompressible elasticity 8 Plate bending K-L 9 Conclusions

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 2 / 71

slide-3
SLIDE 3

Generalities - Reminders on Classic FEM In FEM the degrees of freedom are used to reconstruct polynomials (or isoparametric images of polynomials) in each element. Ingredients: the geometry of the element (e.g.: triangles) the degrees of freedom; say, N d.o.f. per element in each element, a space of polynomials of dim. N. The ingredients must match Unisolvence N numbers ↔ one and only one polynomial Continuity

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 3 / 71

slide-4
SLIDE 4

Virtual Elements on pentagons Assume now that we want to decompose the domain Ω, for instance, in pentagons, obviously not necessarily regular. How to take a polynomial space of dimension 5 (e.g., to be associated to the nodal values)?

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 4 / 71

slide-5
SLIDE 5

Virtual Elements for Laplace Equation VEM: We take, as unknowns, the values at the vertices. Then, for every pentagon E and for every fixed set of 5 vertex-values we define the corresponding function (say, ϕh) first on ∂E, by: ϕh is linear on each edge.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 5 / 71

slide-6
SLIDE 6

Virtual Elements for Laplace Equation VEM: We take, as unknowns, the values at the vertices. Then, for every pentagon E and for every fixed set of 5 vertex-values we define the corresponding function (say, ϕh) first on ∂E, by: ϕh is linear on each edge. Hence, from the vertex values we have ϕh on the whole ∂E. At this point we decide that all our functions, inside E, should be harmonic. Hence, from its 5 nodal values ϕh is first defined on ∂E (by linear interpolation), and from its value on ∂E it is defined inside (by harmonic extension). In summary, the value of ϕh on E is uniquely determined by its nodal values.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 5 / 71

slide-7
SLIDE 7

Virtual Elements for Laplace Equation In E we have therefore a local space V E, of dimension 5, that contains as a subspace the space P1 of polynomials

  • f degree ≤ 1, plus two additional functions that are not

polynomials (and that, unfortunately, we cannot compute,

  • r, at least, not in a cheap-enough way).

E

P

1 1 x y

V

Classical option: use some numerical integration formula. In VEM, however, we proceed differently.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 6 / 71

slide-8
SLIDE 8

Reminder: implementation of Classic FEM Find u ∈ V ≡ H1

0(Ω) s. t. −∆u = f . That is:

  • Ω

∇u · ∇v dΩ =

  • Ω

f v dΩ ∀ v ∈ V . Setting Vh = continuous piecewise linear functions vanishing at the boundary, we look for uh in Vh such that

  • Ω

∇uh · ∇vh dΩ =

  • Ω

f vh dΩ ∀ vh ∈ Vh. The final matrix is then computed as the sum of the contributions of the single elements: Ai,j ≡

  • Ω

∇v j · ∇v i dΩ =

  • E
  • E

∇v j · ∇v i dE.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 7 / 71

slide-9
SLIDE 9

How to use VEM

E

E

P

1 1 x y

V

Vh := {v ∈ V : v linear on each edge, −∆v = 0 in E ∀E} aE(p1, v) =

  • E

∇p1 · ∇v dE =

  • ∂E

∂p1 ∂n v dℓ=: aE

h (p1, v)

If u is in P1(E), then aE(u, v) can be computed exactly. We have then to decide how to chose ahE(u, v) when both u and v are not in P1(E).

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 8 / 71

slide-10
SLIDE 10

Continuous and discretized problem We consider the continuous problem Find u ∈ V ≡ H1

0(Ω)

such that a(u, v) ≡

  • Ω

∇u · ∇v dΩ =

  • Ω

f v dΩ ∀ v ∈ V , and its discretized version: Find uh ∈ Vh such that ah(uh, vh) = (fh, vh) ∀ vh ∈ Vh, and we look for sufficient conditions on ah that ensure all the good properties that you would have with standard Finite Elements

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 9 / 71

slide-11
SLIDE 11

The two basic properties H1 (Consistency) aE

h (p1, v) = aE(p1, v)

∀E, ∀v ∈ V E, ∀p1 ∈ P1(E). H2 (Stability) ∃ α∗, α∗ > 0 such that: α∗ aE(v, v) ≤ aE

h (v, v) ≤ α∗ aE(v, v)

∀E, ∀ v ∈ V E. Under Assumptions H1 and H2 the discrete problem has a unique solution. Moreover the Patch Test of order 1 is satisfied:on any patch of elements, if the exact solution is a global polynomial of degree 1, then the exact solution and the approximate solution coincide. Incidentally: u − uh1 = O(h).

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 10 / 71

slide-12
SLIDE 12

Convergence Theorem

Theorem

Under the above assumptions H1 and H2, for every approximation uI of u in Vh and for every approximation up of u that is piecewise in P1, we have u − uhV ≤ C

  • u − uIV + u − uph,V + f − fh(Vh)′
  • where

f − fh(Vh)′ := sup

vh∈Vh

(f , vh) − (fh, vh) vhV

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 11 / 71

slide-13
SLIDE 13

Proof of convergence

Set δh := uh − uI α∗ αδh2

V ≤ α∗ a(δh, δh) ≤ ah(δh, δh) = ah(uh, δh) − ah(uI, δh)

= (fh, δh) −

  • E

aE

h (uI, δh)

= (fh, δh) −

  • E
  • aE

h (uI − up, δh) + aE h (up, δh)

  • = (fh, δh) −
  • E
  • aE

h (uI − up, δh) + aE(up, δh)

  • = (fh, δh) −
  • E
  • aE

h (uI − up, δh) + aE(up − u, δh)

  • − a(u, δh)

= (fh, δh) −

  • E
  • aE

h (uI − up, δh) + aE(up − u, δh)

  • − (f , δh)

= (fh − f , δh) −

  • E
  • aE

h (uI − up, δh) + aE(up − u, δh)

  • .

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 12 / 71

slide-14
SLIDE 14

How to satisfy H1 and H2 We construct first on each E an operator Π∇

1 from V E into P1(E) defined by

  • Vi=vertex of E

(v−Π∇

1 v)(Vi) = 0

aE(v−Π∇

1 v, p1) = 0 ∀ p1

Note that Π∇

1 p1 = p1 for all p1 in P1(E).

Then we set, for all u and v in V E aE

h (u, v) := aE(Π∇ 1 u, Π∇ 1 v) + S(u − Π∇ 1 u, v − Π∇ 1 v)

where the stabilizing bilinear form S is (for instance) the Euclidean inner product in R5.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 13 / 71

slide-15
SLIDE 15

Proof of H1 and H2 Consistency: aE

h (pk, vh) = aE(pk, Πkvh) = aE(Πkvh, pk) = aE(pk, vh)

Stability (upper bound): aE

h (vh, vh) ≤ aE(Πkvh, Πkvh)+c1aE(vh −Πkvh, vh −Πkvh)

= aE(vh, Πkvh) + c1aE(vh − Πkvh, vh) ≤ α∗aE(vh, vh) Stability (lower bound): aE

h (vh, vh) ≥ aE(Πkvh, Πkvh)+c0aE(vh −Πkvh, vh −Πkvh)

≥ α∗(aE(vh, Πkvh) + aE(vh − Πkvh, vh)) = α∗aE(vh, vh)

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 14 / 71

slide-16
SLIDE 16

Robustness of the method (by A. Russo)

0.2 0.4 0.6 0.8 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 15 / 71

slide-17
SLIDE 17

Robustness of the method (by A. Russo)

0.2 0.4 0.6 0.8 1 0.5 1 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 16 / 71

slide-18
SLIDE 18

Robustness of the method (by A. Russo)

0.2 0.4 0.6 0.8 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 17 / 71

slide-19
SLIDE 19

Robustness of the method (by A. Russo)

0.2 0.4 0.6 0.8 1 0.5 1 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 18 / 71

slide-20
SLIDE 20

Robustness of the method (by A. Russo)

0.2 0.4 0.6 0.8 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 19 / 71

slide-21
SLIDE 21

Robustness of the method (by A. Russo)

0.2 0.4 0.6 0.8 1 0.5 1 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 20 / 71

slide-22
SLIDE 22

Robustness of the method (by A. Russo)

0.2 0.4 0.6 0.8 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 21 / 71

slide-23
SLIDE 23

Robustness of the method (by A. Russo)

0.2 0.4 0.6 0.8 1 0.5 1 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 22 / 71

slide-24
SLIDE 24

Robustness of the method (by A. Russo)

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 512 polygons, 2849 vertices Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 23 / 71

slide-25
SLIDE 25

Robustness of the method Note that the pink element is a polygon with 9 edges, while the blue element is a polygon (not simply connected) with 13 edges. We are exact on linears...

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 24 / 71

slide-26
SLIDE 26

Robustness of the method (by A. Russo)

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 −2 −1 1 2 max |u−uh| = 0.008783

For reasons of ”glastnost”, we take as exact solution w = x(x − 0.3)3(2 − y)2 sin(2πx) sin(2πy) + sin(10xy)

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 25 / 71

slide-27
SLIDE 27

Robustness of the method (by A. Russo)

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 −2 −1 1 2 max |u−uh| = 0.074424

This is on a mesh of 512 (16 × 16 little squares) elements.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 26 / 71

slide-28
SLIDE 28

Robustness of the method (by A. Russo)

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 −2 −1 1 2 max |u−uh| = 0.019380

This is on a mesh of 2048 (32 × 32 little squares) elements.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 27 / 71

slide-29
SLIDE 29

Robustness of the method (by A. Russo)

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 −2 −1 1 2 max |u−uh| = 0.005035

And this is on a mesh of 8192 (64 × 64 little squares) elements.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 28 / 71

slide-30
SLIDE 30

The next steps? (by M.C. Escher)

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 29 / 71

slide-31
SLIDE 31

The next steps? (by M.C. Escher)

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 30 / 71

slide-32
SLIDE 32

Going berserk (by A. Russo)

−0.5 0.5 1 1.5 −0.2 0.2 0.4 0.6 0.8 1 1.2 1.4 1 polygons, 82 vertices

The first step: a pegasus-shaped polygon

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 31 / 71

slide-33
SLIDE 33

Going berserk (by A. Russo)

−0.5 0.5 1 1.5 −0.2 0.2 0.4 0.6 0.8 1 1.2 1.4 1 polygons, 82 vertices 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82

The second step: local numbering of nodes.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 32 / 71

slide-34
SLIDE 34

Going berserk (by A. Russo)

−0.4 −0.2 0.2 0.4 0.6 0.8 1 1.2 1.4 −0.2 0.2 0.4 0.6 0.8 1 1.2 4 polygons, 243 vertices

The third step: a mesh of 2 × 2 pegasus.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 33 / 71

slide-35
SLIDE 35

Going totally berserk (by A. Russo)

0.2 0.4 0.6 0.8 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 400 cells, 16821 vertices

A mesh of 20 × 20 pegasus.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 34 / 71

slide-36
SLIDE 36

Going totally berserk (by A. Russo)

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 −2 −1 1 2 max |u−uh| = 0.077167

Solution on a 20 × 20-pegasus mesh.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 35 / 71

slide-37
SLIDE 37

Going totally berserk (by A. Russo)

0.2 0.4 0.6 0.8 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1600 cells, 65641 vertices

A mesh of 40 × 40 pegasus.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 36 / 71

slide-38
SLIDE 38

Going totally berserk (by A. Russo)

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 −2 −1 1 2 max |u−uh| = 0.026436

Solution on a 40 × 40-pegasus mesh.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 37 / 71

slide-39
SLIDE 39

Higher order VEM . Vh,k := {v ∈ V : v of degree k on each edge, with ∆v ∈ Pk−2 in E ∀E} . In each E the functions in V E are identified by their value at ∂E, (for k > 1)the moments up to the order k − 2 in E One can prove that these d.o.f. are unisolvent. aE(pk, v) =

  • E

∇pk · ∇v dE = −

  • E

∆pk v dE +

  • ∂E

∂pk ∂n v dℓ=: aE

h (pk, v)

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 38 / 71

slide-40
SLIDE 40

Higher order VEM . On an element with n edges, the local V E space of degree k will have: nk d.o.f. at the boundary and k(k − 1)/2 internal d.o.f.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 39 / 71

slide-41
SLIDE 41

Higher order VEM . On an element with n edges, the local V E space of degree k will have: nk d.o.f. at the boundary and k(k − 1)/2 internal d.o.f.

k=3 k=1 k=2 k=4

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 39 / 71

slide-42
SLIDE 42

How to satisfy H1 and H2 We construct first on each E an operator Π∇

k

from V E into Pk(E) defined by

  • E

(v − Π∇

k v) dE = 0

aE(v − Π∇

k v, pk) = 0 ∀ pk ∈ P1

Then we set, for all u and v in V E aE

h (u, v) := aE(Π∇ k u, Π∇ k v) + S(u − Π∇ k u, v − Π∇ k v)

where, here too, the stabilizing bilinear form S is (for instance) the Euclidean inner product in RN. The resulting scheme will satisfy a Patch Test of order k.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 40 / 71

slide-43
SLIDE 43

Linear elasticity problems Consider now a (toy) 2d linear elasticity problem, with (unrealistic) homogeneous kinematic boundary conditions all over ∂Ω . The internal energy, in terms of the displacements u is Internal Energy = µ

  • Ω

ε(u) : ε(v)dx + λ 2

  • Ω

divu divvdx where ε(u) (classical symmetric gradient) is the strain tensor and µ and λ are the classical Lam´ e coefficients. The bilinear form associated with the internal energy is a(u, v) = 2µ

  • Ω

ε(u) : ε(v)dx + λ

  • Ω

divu divvdx.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 41 / 71

slide-44
SLIDE 44

Linear elasticity problems The corresponding Green formula is a(u, v) = −2µ

  • Ω

(div(ε(u)) · vdx − λ

  • Ω

(∇divu) · vdx + 2µ

  • ∂Ω

Mn(u) · vds + λ

  • ∂Ω

divu (v · n)ds where Mn is a first order trace operator. We also set Aµ := −divε, Aλ := −∇div and Aλ,µ := 2µ Aµ + λ Aλ.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 42 / 71

slide-45
SLIDE 45

Towards VEM for linear elasticity aE(u, v) = −2µ

  • E

div(ε(u)) · vdx − λ

  • E

(∇divu) · vdx + 2µ

  • ∂E

Mn(u) · vds + λ

  • ∂E

divu (v · n)ds. Assume that our degrees of freedom allow the reconstruction, at the boundary, of vector valued polynomials of degree k. Assume further that we have, as internal degrees of freedom, the moments up to the order k − 2. Then for every element E, for every u ∈ Pk, and for every v in the discretized subspace the local matrix aE(u, v) is uniquely computable.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 43 / 71

slide-46
SLIDE 46

VEM for elasticity Vh,k := {v ∈ (H1

0(Ω))2 : such that

v ∈ (Pk(e))2 ∀ e, Aλ,µv ∈ (Pk−2(E))2 ∀E} In each E the vectors v in V E are identified by their values at ∂E, (for k > 1) the moments up to the order k − 2 in E One can prove that these d.o.f. are unisolvent. Moreover aE(pk, v) = −2µ

  • E

div(ε(pk))·vdx−λ

  • E

(∇divpk)·vdx + 2µ

  • ∂E

Mn(pk) · vds + λ

  • ∂E

divpk (v · n)ds =: 2µ aE

h,µ(pk, v) + λ aE h,λ(pk, v).

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 44 / 71

slide-47
SLIDE 47

Linear Elasticity Elements

2

d=24+6 d=16+2 d=8 k=1 k=2 k=3

(P )

1 + 2

(P )

2

2+ 6

(P )

3

2+ 10

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 45 / 71

slide-48
SLIDE 48

How to satisfy H1 and H2 Again, we construct on each E a projection

  • perator Πel

k : V E → (Pk(E))2 defined by

  • vertices

(v − Πel

k v) = 0, and aE(v − Πel k v, pk) = 0 ∀ pk,

and then set, for all u and v in V E aE

h (u, v) := aE(Πel k u, Πel k v) + S(u − Πel k u, v − Πel k v)

where the stabilizing bilinear form S, once more, is for instance the Euclidean inner product in RN. The resulting scheme will again satisfy a Patch Test of order k

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 46 / 71

slide-49
SLIDE 49

Nearly incompressible elasticity To understand what to do for the nearly incompressible case (λ >> µ) it is convenient (as usual) to consider the (u, p) formulation, introducing p := λu as an additional unknown. The simplest approach consists in introducing the space Qh of discontinuous local polynomials of degree k − 1 (with zero global mean value), and define an operator Div from Vh to Qh by

  • E

Divvhqh dx = −

  • E

vh · ∇qh dx +

  • ∂E

vh · nqh ds Note that ∇qh ∈ (Pk−2)2 and hence Divvh is computable.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 47 / 71

slide-50
SLIDE 50

Nearly incompressible elasticity The discrete problem will then be        Find (uh, ph) ∈ Vh × Qh such that 2µ ah,µ(uh, vh) + (ph, Divvh) = (f, vh) ∀ vh ∈ Vh (qh, Divuh) − 1 λ(ph, qh) = 0 ∀ qh ∈ Qh It is not difficult to check that, for k ≥ 2, the inf-sup condition is satisfied for the operator Div, and that this implies optimal order convergence. Then we learn that, in the original formulation, we can write 2µ ah,µ(uh, vh) + λ(Divuh, Divvh) = (f, vh) ∀ vh ∈ Vh

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 48 / 71

slide-51
SLIDE 51

Nearly incompressible Elasticity and Stokes

k=3

u u p p

k=2

These elements could be seen as a generalization of the (Pk + Bk+1)2-(discontinuous Pk−1) elements for Stokes.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 49 / 71

slide-52
SLIDE 52

Plate bending - Kirchhoff-Love We consider now the Kirchhoff-Love model for the bending of thin plates. At rest the midsection of the plate occupies the region Ω. After deformation the (scaled) total energy of the plate could be written as 1 2

  • Ω

D

  • (1 − ν) w/ijw/ij + νw/iiw/jj
  • dx −
  • Ω

f w dx where w=transversal displacement, f =external load density, ν=Poisson ratio, D=elastic coefficient. It is well known that the approximation of plate problems (in the K-L formulation) requires the use of C 1 elements.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 50 / 71

slide-53
SLIDE 53

C 1 Finite Elements There are relatively few C 1 Finite Elements on the

  • market. Here are some:

Bell HCT reduced HCT Argyris

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 51 / 71

slide-54
SLIDE 54

Programming C 1 elements Cod liver oil (Dorschleber¨

  • l)

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 52 / 71

slide-55
SLIDE 55

General idea of VEM for plates Step 1. We first choose the degree, r,of our shape functions, and the degree, s, of their normal derivative on each edge. Step 2. We choose d.o.f. such that: i) they determine uniquely v and vn on ∂E, ii) their restriction on each edge determine uniquely the value of v and vn on that edge. Step 3. Finally we fix a polynomial degree, k such that k ≤ r and k − 1 ≤ s. This will be our order of accuracy. Step 4. If k ≥ 4 we then add to the boundary degrees of freedom the moments

  • E

v q dx, q ∈ Pk−4(E).

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 53 / 71

slide-56
SLIDE 56

Constructing VEM for plates For r, s, k, with r ≥ k and s ≥ k − 1 we set Vh := {v ∈ V : v ∈ Pr(e), vn ∈ Ps(e) ∀ edge e and ∆2v ∈ Pk−4(E) ∀ element E} In each E the functions in V E are identified by their value and the value of their derivatives on ∂E, (for k > 3) the moments up to the order k − 4 in E These d.o.f. are unisolvent. Moreover, for p ∈ Pk aE(v, p) :=

  • E

D

  • (1 − ν) v/ijp/ij + νv/iip/jj
  • dx

= D

  • E

v ∆2p dx+

  • ∂E

∇v·Mn(p) −v Qn(p) dℓ =: aE

h (v, p)

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 54 / 71

slide-57
SLIDE 57

Example 1

r=3 d=12 k=2 s=1

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 55 / 71

slide-58
SLIDE 58

Example 2

s=2 d=16 r=3 k=3

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 56 / 71

slide-59
SLIDE 59

Example 3

d=24+1 r=4 s=3 k=4

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 57 / 71

slide-60
SLIDE 60

Example 4

r=5 s=3 k=4 d=24+1

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 58 / 71

slide-61
SLIDE 61

Example 5

d=28+3 r=5 s=4 k=5

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 59 / 71

slide-62
SLIDE 62

Example 6 For the same r, s, and k we might use different degrees of freedom.

d=28+1 r=5 s=3 k=4 d=24+1 r=5 s=3 k=4

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 60 / 71

slide-63
SLIDE 63

Example 7 For the same r, s, and k we might use different degrees of freedom.

d=32+3 d=28+3 r=5 k=5 s=4 r=5 k=5 s=4

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 61 / 71

slide-64
SLIDE 64

How to satisfy H1 and H2 We construct on each E a projection

  • perator ΠM

k from V E into Pk(E) defined by

  • vertices

(v − ΠM

k v) = 0,

  • vertices

∇(v − ΠM

k v) = 0 and

aE(v − ΠM

k v, pk) = 0 ∀ pk,

and then set, for all u and v in V E aE

h (u, v) := aE(ΠM k u, ΠM k v) + S(u − ΠM k u, v − ΠM k v)

where the stabilizing bilinear form is again (for instance) the Euclidean inner product in RN. The resulting scheme will again satisfy a Patch Test of order k

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 62 / 71

slide-65
SLIDE 65

Numerical experiments on the 3-2 element

s=2 d=16 r=3 k=3

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 63 / 71

slide-66
SLIDE 66

Exact and approximate solution

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 −1 1 2 3 4 x 10

−3

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 −1 1 2 3 4 x 10

−3

Exact solution (left): w = x2(1 − x)2y 2(1 − y)2 on the unit square ]0, 1[×]0, 1[. The approximate solution is computed with the r = 3, s = 2, k = 3 element on a grid

  • f uniform 32 × 32 square (BLUSH!) elements.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 64 / 71

slide-67
SLIDE 67

Behaviour of the L2 error

10

−1

10

−6

10

−5

10

−4

10

−3

10

−2

L2 −error log h ||wI−wh||0

error h4

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 65 / 71

slide-68
SLIDE 68

Behaviour of the H1 error

10

−1

10

−6

10

−5

10

−4

10

−3

10

−2

H1 −error log h |wI−wh|1

error h4 h3

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 66 / 71

slide-69
SLIDE 69

Behaviour of the H2 error

10

−1

10

−8

10

−6

10

−4

10

−2

H2 −error log h |wI−wh|2

error h2 h3 h4

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 67 / 71

slide-70
SLIDE 70

VEM on Triangles

+

v, Dv v, Dv, D v vn v, v , v

n nt 2 Internal moments

2+ 3

P3 + 2 P P

4 + 4

P4+ 4 P

5 3

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 68 / 71

slide-71
SLIDE 71

VEM on general shapes

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 69 / 71

slide-72
SLIDE 72

VEM on general shapes

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 70 / 71

slide-73
SLIDE 73

Conclusions

  • VEM preserve all the good features of MFD, but are

much more simple and more elegant.

  • In some cases, they could be seen ”just” as a

generalization to more general geometries of classical FEM

  • In other cases, they offer additional possibilities. This is

for instance the case for general quadrilaterals, hanging nodes, and C k elements with k ≥ 1.

  • They are almost newborn. To assess their true interest

for engineering computations still requires a long long way to go.

Franco Brezzi (IMATI-CNR & IUSS, Pavia) Virtual Elements Milano, September 17-19, 2012 71 / 71