[PPT] - IsoGeometric Analysis: B ezier techniques in Numerical Simulations PowerPoint Presentation

SLIDE 1

IsoGeometric Analysis: B´ ezier techniques in Numerical Simulations

Ahmed Ratnani

IPP, Garching, Germany July 30, 2015

A. Ratnani

IGA – CEMRACS 2015 1/1

1/1

SLIDE 2

Outline

Motivations Computer Aided Design (CAD) Zoology B´

ezier curves

Tensor Product surfaces B-Spline curves B´

ezier triangular surfaces

Splines/NURBS Finite Elements The IsoGeometric Approach Impact of the k-refinement strategy The discrete DeRham diagram Maxwell equations MultiGrid Methods Adaptive meshes B´

ezier techniques in Computational Plasma Physics

A. Ratnani

IGA – CEMRACS 2015 2/1

2/1

SLIDE 3

Motivations

Engineering Analysis Process

Finite Elements Analysis (FEA) models are created from CAD

representations

Fixing CAD geometry and creating FEA models accounts more

than 80% of overall analysis time and is a major engineering bottelneck

The geometry is approximated in the FEA mesh

A. Ratnani

IGA – CEMRACS 2015 3/1

3/1

SLIDE 4

Motivations

Engineering Analysis Process

Finite Elements Analysis (FEA) models are created from CAD

representations

Fixing CAD geometry and creating FEA models accounts more

than 80% of overall analysis time and is a major engineering bottelneck

The geometry is approximated in the FEA mesh

A. Ratnani

IGA – CEMRACS 2015 3/1

3/1

SLIDE 5

Motivations

A. Ratnani

IGA – CEMRACS 2015 4/1

4/1

SLIDE 6

Motivations

Even if the p-form: C(x(t)) = ∑n

i=0 tiPi, is a natural description for

curves, it presents some disadvantages :

the curve is not necessary regular everywhere. ➠ non-efficient

approximation,

the points (Pi)0≤i≤n do not have any geometric interpretation, unstable numerical evaluation.

A. Ratnani

IGA – CEMRACS 2015 5/1

5/1

SLIDE 7

Motivations

Even if the p-form: C(x(t)) = ∑n

i=0 tiPi, is a natural description for

curves, it presents some disadvantages :

the curve is not necessary regular everywhere. ➠ non-efficient

approximation,

the points (Pi)0≤i≤n do not have any geometric interpretation, unstable numerical evaluation.

Properties ➠ Bezier techniques provide a geometric-based method for describing and manipulating polynomial curves and surfaces.

brings sophisticated mathematical concepts into a highly

geometric and intuitive form.

this form facilitates the creative design process. Bezier techniques are an excellent choice in the context of

numerical stability of floating point operations[Farouki & Rajan]. ➠ Bezier techniques are at the core of 3D Modeling or Computer Aided Geometric Design (CAGD).

A. Ratnani

IGA – CEMRACS 2015 5/1

5/1

SLIDE 8

Computer Aided Design

B´ ezier curves

Rather than use {1, t, · · · , tn} as a basis of Π<n+1, we can take Bernstein polynomials; this leads to the B´ ezier-form. Therefore, it is equivalent to the p-form and writes : C(x(t)) =

n

∑

i=0

Bn

i (t)Pi,

0 ≤ t ≤ 1 (1) where Bn

i denote Bernstein polynomials :

Bn

i (t) =

n i

ti(1 − t)n−i =

n! i!(n − i)!ti(1 − t)n−i, 0 ≤ t ≤ 1 (2) The sequence (Pi)0≤i≤n is called control points. Conics can be descrived exactly using Non-Rational B´ ezier arcs.

A. Ratnani

IGA – CEMRACS 2015 6/1

6/1

SLIDE 9

Computer Aided Design

B´ ezier curves

A. Ratnani

IGA – CEMRACS 2015 7/1

7/1

SLIDE 10

Computer Aided Design

B´ ezier curves

A. Ratnani

IGA – CEMRACS 2015 7/1

7/1

SLIDE 11

Computer Aided Design

B´ ezier curves

A. Ratnani

IGA – CEMRACS 2015 7/1

7/1

SLIDE 12

Computer Aided Design

B´ ezier curves

A. Ratnani

IGA – CEMRACS 2015 7/1

7/1

SLIDE 13

Computer Aided Design

B´ ezier curves

A. Ratnani

IGA – CEMRACS 2015 7/1

7/1

SLIDE 14

Computer Aided Design

B´ ezier curves

Invariance under some transformations : rotation, translation, scaling; it is sufficient to

transform the control points,

Bn i (t) ≥ 0, ∀ 0 ≤ t ≤ 1 partition of unity : ∑n i=0 Bn i (t) = 1, ∀ 0 ≤ t ≤ 1 Bn 0 (0) = Bn n (1) = 1 each Bn i has exactly one maximum in [0, 1], at i n Bn i are symmetric with respect to 1 2 recursive property : Bn i (t) = (1 − t)Bn−1 i

(t) + tBn−1

i−1 (t), and

Bn

i (t) = 0,

if i < 0 or, i > n

deriving a curve : C′(t) = n{∑n−1 i=0 Bn−1 i

(t) (Pi+1 − Pi)}, then : C′(0) = n (P1 − P0) C′(1) = n (Pn − Pn−1) (3) C′′(0) = n(n − 1) (P0 − 2P1 + P2) C′′(1) = n (Pn − 2Pn−1 + Pn−2) (4)

DeCasteljau algorithm:

Cn(t; P0, · · · , Pn) = (1 − t)Cn−1(t; P0, · · · , Pn−1) + tCn−1(t; P1, · · · , Pn) (5)

A. Ratnani

IGA – CEMRACS 2015 8/1

8/1

SLIDE 15

Computer Aided Design

Tensor Product surfaces

Bezier patchs of arbitrary degrees

A B´ ezier patch of degrees (p, q) is defined as x(s, t) =

p,q

∑

i,j=0

xijBi(s)Bj(t), s, t ∈ [0, 1] (6) where (xij)0≤i≤p,0≤j≤q are called control points Properties

Endpoint interpolation: The patch passes through the four corner control points

{(s, t) = (0, 0), (0, 1), (1, 0), (1, 1)},

Each boundary corredponds to a Bezier curve Symmetry in the parametric domain Affine invariance (applied to the control points) Convex Hull, C 1 patchs can be easily created by solving (local/global) linear systems, using the

endpoints derivatives and moving some specific control points,

A. Ratnani

IGA – CEMRACS 2015 9/1

9/1

SLIDE 16

Computer Aided Design

Tensor Product surfaces

Geometric Operations

(Exact) Subdivision, (Exact) Degree Elevation, (Inexact) Patchs merge, (Exact if a Spline description is used) (Inexact) Degree Reduction

Figure : A mapping as a cubic Bezier patch (left) parametric domain with its domain

points, (right) the resulting physical domain

A. Ratnani

IGA – CEMRACS 2015 10/1

10/1

SLIDE 17

Computer Aided Design

B-Splines curves

For a fixed number of control points, we have a fixed number of

degrees of freedom to control the B´ ezier curve (= p + 1).

For a better control of the curve, one can subdivise it into a given

number of B´ ezier curves (using a refinement algorithm).

How to insure that the local regularity of the curve is preserved,

when controling these curves? Need the notion of Macro-Elements or Macro-Patchs, where given regularities are imposed between elements.

A. Ratnani

IGA – CEMRACS 2015 11/1

11/1

SLIDE 18

Computer Aided Design

B-Splines curves

(a) (b) (a) Original curve given as B´ ezier curves. (a) The quadratic B-spline curve and its control points. The knot vector is T = {000, 1

4 , 1 2 , 3 4 , 111}.

A. Ratnani

IGA – CEMRACS 2015 12/1

12/1

SLIDE 19

Computer Aided Design

B-Splines curves

Figure : (left) A quadratic B-Spline curve and its control points using the knot vector T = {000 1

2 3 4 3 4 111}, (right) the corresponding B-Splines.

A. Ratnani

IGA – CEMRACS 2015 13/1

13/1

SLIDE 20

Computer Aided Design

B-Splines

To create a family of B-splines, we need a non-decreasing sequence of knots T = (ti)1iN+k, also called knot vector, with k = p + 1. Each set of knots Tj = {tj, · · · , tj+p} will generate a B-spline Nj.

Definition (B-Spline serie)

The j-th B-Spline of order k is defined by the recurrence relation: Nk

j = wk j Nk−1 j

+ (1 − wk

j+1)Nk−1 j+1

where, wk

j (x) =

x − tj tj+k−1 − tj N1

j (x) = χ[tj,tj+1[(x)

for k ≥ 1 and 1 ≤ j ≤ N.

A. Ratnani

IGA – CEMRACS 2015 14/1

14/1

SLIDE 21

Computer Aided Design

B-Splines

0.5 1 1.5 2 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 N1 N2 N3 N4 N5 N6 N7 N8 0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 N1 N2 N3 N4 N5 N6 N7 N8 0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 N1 N2 N3 N4 N5 N6 N7 N8

Figure : B-splines functions associated to the knot vector T = {000 1 2 3 44 555}, of

rder k = 1, 2, 3

Figure : Quadratic B-Splines for T = {000, 111}, T = {000, 1

2 , 111} and

T = {000, 1

2 , 3 4 3 4 , 111}.

A. Ratnani

IGA – CEMRACS 2015 15/1

15/1

SLIDE 22

Computer Aided Design

B´ ezier triangular surfaces

Barycentric coordinates

Let T = {v1, v2, v3} be a non-degenerate triangle in the plance. Then for all point P in the plane, there exists τ = {τ1, τ2, τ3} such that P = ∑3

i=1 τivi. τ is unique if one add the normalization constraint ∑3 i=1 τi = 1 (will be assumed

during this talk)

P ∈ T if and only if τ ≥ 0 Affine invariance: if the triangle T together with the point P are transformed by an

affine transformation, the transformed point has unchanged barycentric coordinates

Bernstein polynomials on triangles

Let λ be a multi-index such that |λ| = n, T a triangle, and x a point in the plane, with τ as barycentric coordinates with respect to T . Bernstein polynomials are defined using the barycentric coordinates. Bn

λ(τ) =

n! λ1!λ2!λ3! τ1λ1τ2λ2τ3λ3 (7) Let ξijk = iv1+jv2+kv3

d

. The set Dd,T = {ξijk, i + j + k = d} is the set of domain-points.

A. Ratnani

IGA – CEMRACS 2015 16/1

16/1

SLIDE 23

Computer Aided Design

B´ ezier triangular surfaces

Properties of Bernstein polynomials

Partition of Unity: ∑|λ|=n Bn λ = 1, Positivity: Bn λ(τ) ≥ 0 if and only if τ ≥ 0, Bn λ has its maximum at τ = λ n (a domain point)

Bernstein triangular patchs

As for the rectangular case, the Bernstein patch can be defined as x(τ) = ∑

|λ|=n

xλBn

λ(τ),

∀τ ∈ [0, 1] (8) Properties of Bernstein Series

Endpoint interpolation: The patch passes through the three corner control points

that are x(1, 0, 0) = v1, x(0, 1, 0) = v2, x(0, 0, 1) = v3

Each boundary corredponds to a Bezier curve Symmetry in the parametric domain, Convex Hull Affine invariance (applied to the control points) C 1 patchs can be easily created by solving (local) linear systems, using the

endpoints derivatives and moving some specific control points,

A. Ratnani

IGA – CEMRACS 2015 17/1

17/1

SLIDE 24

Computer Aided Design

B´ ezier triangular surfaces

Geometric Operations

(Exact) Subdivision, (Exact) Degree Elevation, (Inexact) Patchs merge, (Exact if a Spline description is used) (Inexact) Degree Reduction ξ003 ξ102 ξ012 ξ201 ξ111 ξ021 ξ300 ξ210 ξ120 ξ030 x x x x x x x c003 c102 c012 c201 c111 c021 c300 c210 c120 c030 x x x x x x x

Figure : (left) Domain points and (right) B-coefficients for a cubic polynomial

A. Ratnani

IGA – CEMRACS 2015 18/1

18/1

SLIDE 25

Computer Aided Design

B´ ezier triangular surfaces

Figure : Examples of Bezier triangulations (left) quadratic (middle) cubic (right)

quadratic, curved triangles

A. Ratnani

IGA – CEMRACS 2015 19/1

19/1

SLIDE 26

Splines/NURBS Finite Elements

The IsoGeometric Approach

Q F

Patch Physical Domain

K

Q

F

Patch Physical Domain K

Grid generation: the use of h/p/k-refinement keeps the mapping F unchanged.

Compact support Partition of Unity Affine covariance IsoParametric concept Error estimates in Sobolev norms

A. Ratnani

IGA – CEMRACS 2015 20/1

20/1

SLIDE 27

Splines/NURBS Finite Elements

The IsoGeometric Approach

Refinement strategies

Refining the grid can be done in 3 different ways. This is the most interesting aspects of B-splines basis. h-refinement by inserting new knots. It is the equivalent of mesh refinement of the classical finite element method. p-refinement by elevating the B-spline degree. It is the equivalent of using higher finite element order in the classical FEM. k-refinement by increasing / decreasing the regularity of the basis functions (increasing / decreasing multiplicity of inserted knots). the use of k-refinement strategy is more efficient than the classical p-refinement, as it reduces the dimension of the basis.

A. Ratnani

IGA – CEMRACS 2015 21/1

21/1

SLIDE 28

Splines/NURBS Finite Elements

Impact of the k-refinement strategy

In the following table, we show the impact of the k-refinement on the resolution of the Poisson’s equation, on a square domain using:

B-splines of degree p and minimal regularity (i.e. C0) B-splines of degree p and maximal regularity (i.e. Cp−1)

number of d.o.f number of nnz cpu-SuperLU cpu-CG Cp−1 C0 Cp−1 C0 Cp−1 C0 Cp−1 C0 p=2 4’096 16’129 98’596 253’009 0.23 0.35 7 10−4 4.1 10−3 p=3 4’225 36’481 196’249 896’809 0.61 1.64 1.1 10−2 2 10−2 p=5 4’489 101’761 499’849 4’923’961 2.96 49.27 3.8 10−2 3.7 10−1

Table : Impact of the k-refinement on the resolution of the Poisson equation

n a grid 64 × 64 for quadratic, cubic and quintic B-splines.
A. Ratnani

IGA – CEMRACS 2015 22/1

22/1

SLIDE 29

Splines/NURBS Finite Elements

The discrete DeRham diagram

Vector fields transformations

grad curl div H1(Ω) − → H(curl, Ω) − → H(div, Ω) − → L2(Ω) ı0 ↑ ı1 ↑ ı2 ↑ ı3 ↑

grad
curl
div

H1(P) − → H(curl, P) − → H(div, P) − → L2(P)

Vector fields approximations

grad
curl
div

H1(P) − → H(curl, P) − → H(div, P) − → L2(P)

Π0↓
Π1 ↓
Π2 ↓
Π3 ↓
grad
curl
div

V − → Wcurl − → Wdiv − → X Sp,p

α,α

A. Ratnani

IGA – CEMRACS 2015 23/1

23/1

SLIDE 30

Splines/NURBS Finite Elements

The discrete DeRham diagram

Vector fields transformations

grad curl div H1(Ω) − → H(curl, Ω) − → H(div, Ω) − → L2(Ω) ı0 ↑ ı1 ↑ ı2 ↑ ı3 ↑

grad
curl
div

H1(P) − → H(curl, P) − → H(div, P) − → L2(P)

Vector fields approximations

grad
curl
div

H1(P) − → H(curl, P) − → H(div, P) − → L2(P)

Π0↓
Π1 ↓
Π2 ↓
Π3 ↓
grad
curl
div

V − → Wcurl − → Wdiv − → X Sp−1,p

α−1,α × Sp,p−1 α,α−1

A. Ratnani

IGA – CEMRACS 2015 23/1

23/1

SLIDE 31

Splines/NURBS Finite Elements

The discrete DeRham diagram

Vector fields transformations

grad curl div H1(Ω) − → H(curl, Ω) − → H(div, Ω) − → L2(Ω) ı0 ↑ ı1 ↑ ı2 ↑ ı3 ↑

grad
curl
div

H1(P) − → H(curl, P) − → H(div, P) − → L2(P)

Vector fields approximations

grad
curl
div

H1(P) − → H(curl, P) − → H(div, P) − → L2(P)

Π0↓
Π1 ↓
Π2 ↓
Π3 ↓
grad
curl
div

V − → Wcurl − → Wdiv − → X Sp,p−1

α,α−1 × Sp−1,p α−1,α

A. Ratnani

IGA – CEMRACS 2015 23/1

23/1

SLIDE 32

Splines/NURBS Finite Elements

The discrete DeRham diagram

Vector fields transformations

grad curl div H1(Ω) − → H(curl, Ω) − → H(div, Ω) − → L2(Ω) ı0 ↑ ı1 ↑ ı2 ↑ ı3 ↑

grad
curl
div

H1(P) − → H(curl, P) − → H(div, P) − → L2(P)

Vector fields approximations

grad
curl
div

H1(P) − → H(curl, P) − → H(div, P) − → L2(P)

Π0↓
Π1 ↓
Π2 ↓
Π3 ↓
grad
curl
div

V − → Wcurl − → Wdiv − → X Sp−1,p−1

α−1,α−1

A. Ratnani

IGA – CEMRACS 2015 23/1

23/1

SLIDE 33

Splines/NURBS Finite Elements

Maxwell Equations

We shall only consider in the sequel the TE mode which reads ∂E ∂t − rot H = −J, ∂H ∂t + rot E = 0, div E = ρ. (9)

1st variational formulation

Find (E, H) ∈ H0(curl, Ω) × L2(Ω) such that d dt

Ω E · ψ dX −
Ω H(rot ψ) dX = −
Ω J · ψ dX

∀ψ ∈ H0(curl, Ω), (10) d dt

Ω Hϕ dX +
Ω(rot E)ϕ dX = 0

∀ϕ ∈ L2(Ω). (11)

2nd variational formulation

Find (E, H) ∈ H(div, Ω) × H1(Ω) such that d dt

Ω E · ψ dX −
Ω(rot H) · ψ dX = −
Ω J · ψ dX

∀ψ ∈ H(div, Ω), (12) d dt

Ω Hϕ dX +
Ω E · (rot ϕ) dX = 0

∀ϕ ∈ H1(Ω). (13)

A. Ratnani

IGA – CEMRACS 2015 24/1

24/1

SLIDE 34

Splines/NURBS Finite Elements

Maxwell Equations

In the case of the 2nd formulation, we look for Eh ∈ Wdiv and Hh ∈ V . the linear system writes:

Without DeRham sequence

MW ˙

e = Kh MV ˙ h = −K Te

Using DeRham sequence

˙

e = Rh MV ˙ h = −K Te

Time discreatization using Leap-Frog scheme 2nd or 4th, Solving only one linear system at each time step,

A. Ratnani

IGA – CEMRACS 2015 25/1

25/1

SLIDE 35

Splines/NURBS Finite Elements

Maxwell Equations

In the case of the 2nd formulation, we look for Eh ∈ Wdiv and Hh ∈ V . the linear system writes:

Without DeRham sequence

MW ˙

e = Kh MV ˙ h = −K Te

Using DeRham sequence

˙

e = Rh MV ˙ h = −K Te

Time discreatization using Leap-Frog scheme 2nd or 4th, Solving only one linear system at each time step,

→ the use of regular elements leads to better CFL numbers.

A. Ratnani

IGA – CEMRACS 2015 25/1

25/1

SLIDE 36

Splines/NURBS Finite Elements

Maxwell Equations

Test on a square domain,

1e-13 1e-12 1e-11 1e-10 1e-09 1e-08 1e-07 1e-06 1e-05 0.0001 0.001 0.01 0.1 error (L2 norm) h quadratic Ch3 cubic Ch4 quartic Ch5 quintic Ch6 1e-12 1e-11 1e-10 1e-09 1e-08 1e-07 1e-06 1e-05 0.0001 0.001 0.001 0.01 0.1 error (L2 norm) h quadratic Ch2 cubic Ch3 quartic Ch4 quintic Ch5

Figure : Square test: the L2 norm error for (left) the magnetic field , (right) the electric field

A. Ratnani

IGA – CEMRACS 2015 26/1

26/1

SLIDE 37

Splines/NURBS Finite Elements

Maxwell Equations

Test on a square domain,

20000 40000 60000 80000 100000 120000 140000 0.01 0.02 0.03 0.04 0.05 0.06 0.07 dim W, m=1 dim W, m=2 dim V, m=1 dim V, m=2 1e-08 1e-07 1e-06 1e-05 0.0001 0.001 0.01 0.1 error (L2 norm) h m=1 m=2

Figure : Square test: (left) the dimension of the discrete spaces Wh and Vh , (right) the L2 norm error for the electric field, where the vector knots are multiplicity m = 1, 2 for quadratic B-splines

A. Ratnani

IGA – CEMRACS 2015 26/1

26/1

SLIDE 38

Splines/NURBS Finite Elements

Maxwell Equations

LF2Th LF2num LF4Th LF4num p = 2 0.3044 0.3056 0.8676 0.8720 p = 3 0.2058 0.1840 0.5866 0.5872 p = 4 0.1496 0.1520 0.4265 0.4272 p = 5 0.1151 0.1168 0.3281 0.3280

Table : Test case 1: CFL numbers (theoretical and numerical values), for splines of degree p = 2, · · · , 5

m=1 m=2 p = 2 0.8720 0.5200 p = 3 0.5872 0.4224 p = 4 0.4272 0.3056 p = 5 0.3280 0.2304

Table : Test case 1: CFL, using LF4, for splines of degee p = 2, · · · , 5 for singular knots (m = 1), and doubled knots (m = 2)

A. Ratnani

IGA – CEMRACS 2015 27/1

27/1

SLIDE 39

Splines/NURBS Finite Elements

MultiGrids Methods

One can insert a new knot t, where tj t < tj+1. Let N = N + 1,

k = k,
T = {t1, .., tj, t, tj+1, .., tN+k}

and αi =    1 1 i j − k + 1

t−ti ti+k−1−ti

j − k + 2 i j j + 1 i If Q are the new control points after the insertion of the knot t, then Qi = αiPi + (1 − αi)Pi−1 or equivalently Q = AP The basis transformation A is called the knot insertion matrix of degree k − 1 from T to

T.

Now let us consider a nested sequence of knot vectors T0 ⊂ T1 ⊂ ... ⊂ Tn, where #(Ti+1 − Ti) = 1. The knot insertion matrix from Ti to Ti+1 is denoted by Ai+1

i

. It is easy to see that the insertion matrix from T0 to Tn is simply: A := An

0 = A1 0 A2 1 ... An n−1

In the case of 2D, the interpolation operator can be constructed using the Kronecker product, as follows Aη ⊗ Aξ

A. Ratnani

IGA – CEMRACS 2015 28/1

28/1

SLIDE 40

Splines/NURBS Finite Elements

MultiGrids Methods: Numerical results

A. Ratnani

IGA – CEMRACS 2015 29/1

29/1

SLIDE 41

Splines/NURBS Finite Elements

MultiGrids Methods: Numerical results

A. Ratnani

IGA – CEMRACS 2015 30/1

30/1

SLIDE 42

Splines/NURBS Finite Elements

MultiGrids Methods: Numerical results

A. Ratnani

IGA – CEMRACS 2015 31/1

31/1

SLIDE 43

Adaptive meshes

The r-refinement strategy – Basic ideas

Idea: Move the control points in order to have a better resolution in high density regions.

Iteratively, by minimizing an estimation of the numerical error (see B. Mourrain

papers).

To ensure an equi-distribution property: needs a monitor function (works also with a

posteriori-estimates). Depending on the method, specific conditions (on the boundary) must be hold in order to keep the exact geometry. ➠ Injectivity property: the geometric transformation F must be a one-to-one map

Figure : Example of flux aligned mesh, used in MHD

A. Ratnani

IGA – CEMRACS 2015 32/1

32/1

SLIDE 44

Adaptive meshes

Optimal transport

Definition (L2 Monge-Kantorovich problem)

Let ρ0 and ρ1 be two given densities of equal masses, defined in Ω ⊂ Rd. Find a mapping x′ = Ψ(x), x, x′ ∈ Ω, that transfers the density ρ0 to ρ1 and minimizes the transport cost J [Ψ] =

Ω |Ψ(x) − x|2ρ0(x) dx

(14) The density transfert means that, ∀ω ⊂ Ω,

Ψ−1(ω) ρ0(x) dx =
ω ρ1(x′) dx′

(15) If the mapping Ψ is a smooth one-to-one map, then (Eq. ??) writes ρ1(Ψ(x)) det ∂Ψ ∂x

= ρ0(x)

(16)

A. Ratnani

IGA – CEMRACS 2015 33/1

33/1

SLIDE 45

Adaptive meshes

Optimal transport and MA eq.

Existence theorem [Brenier’91]

There exists a unique optimal mapping that satisfies the

equidistibution principle.

This mapping can be written as the gradient of a convex function

φ

φ is the solution of the Monge-Amp`

ere eq. det H(φ) = σ |Ωc|ρ(∇ξφ) (17) for the boundary conditions, we ensure that x(ξ) maps ∂Ωc to ∂Ω: ∇ξφ(∂Ωc) = ∂Ω (18)

A. Ratnani

IGA – CEMRACS 2015 34/1

34/1

SLIDE 46

Adaptive meshes

Monge-Amp` ere equation using Benamou-Froese-Oberman Method (BFO) [Benamou2010]

Let us define the operator T : H2(Ω) → H2(Ω) by T[u] = ∇2−1 (∇2u)2 + 2(f − det H(u)) (19) when u ∈ H2(Ω) is a solution of the Monge-Amp` ere equation, then it is a fixed point of the operator T.

Algorithm

Given an initial value u0, Compute un+1 as the solution of

∇2un+1 =

(∇2un)2 + 2(f − det H(un))
A. Ratnani

IGA – CEMRACS 2015 35/1

35/1

SLIDE 47

Adaptive meshes

Anisotropic case, on a square domain

A. Ratnani

IGA – CEMRACS 2015 36/1

36/1

SLIDE 48

Adaptive meshes

Anisotropic case, on a square domain

A. Ratnani

IGA – CEMRACS 2015 37/1

37/1

SLIDE 49

Adaptive meshes

Anisotropic case, on a square domain [Delzanno08] [Sulman11] Number of cells Error CPU-time Error CPU-time 16 × 16 4.22 × 10−1 0.2 6.8 × 10−2 0.01 32 × 32 2.17 × 10−1 0.5 3.28 × 10−2 0.09 64 × 64 9.45 × 10−2 1.8 1.66 × 10−2 0.6 128 × 128 2.88 × 10−2 6.9 8.7 × 10−3 5 256 × 256 7.16 × 10−3 27 4.5 × 10−3 33.6 512 × 512 1.76 × 10−3 109 p=2 p=3 p=5 Grid Error CPU Error CPU Error CPU 8 × 8 6.22 × 10−2 0.74 2.05 × 10−2 0.85 5.48 × 10−3 1.96 16 × 16 1.45 × 10−2 1.81 1.62 × 10−3 3.13 6.19 × 10−4 5.9 32 × 32 2.97 × 10−3 8.2 1.06 × 10−4 12.21 5.17 × 10−6 22.84 64 × 64 7.19 × 10−4 37.31 5.27 × 10−6 50.62 7.46 × 10−8 92.19 128 × 128 1.77 × 10−4 163.26 3.3 × 10−7 213.1 1.49 × 10−8 386.53

Table : Example 1: CPU-time needed to solve the Monge-Amp` ere equation and the adaptive Error Eadp using the two grids Picard algorithm using quadratic, cubic and quintic B-splines.

A. Ratnani

IGA – CEMRACS 2015 38/1

38/1

SLIDE 50

Adaptive meshes

Anisotropic case, on a square domain

A. Ratnani

IGA – CEMRACS 2015 39/1

39/1

SLIDE 51

Adaptive meshes

Anisotropic case, on a square domain [Delzanno08] [Sulman11] Number of cells Error CPU-time Error CPU-time 16 × 16 4.22 × 10−1 0.2 2.8 × 10−2 0.03 32 × 32 2.17 × 10−1 0.5 5.6 × 10−3 0.8 64 × 64 9.45 × 10−2 1.8 1.4 × 10−3 14.6 128 × 128 2.88 × 10−2 6.9 3.5 × 10−4 169 256 × 256 7.16 × 10−3 27 512 × 512 1.76 × 10−3 109 p=2 p=3 p=5 Grid Error CPU Error CPU Error CPU 8 × 8 2.37 × 10−1 0.46 2.12 × 10−1 1.29 1.78 × 10−1 2.59 16 × 16 1.38 × 10−1 1.65 1.31 × 10−1 2.96 1.16 × 10−1 10.65 32 × 32 6.35 × 10−2 9.59 1.75 × 10−2 17.81 8.06 × 10−3 38.23 64 × 64 1.31 × 10−2 43.21 7.9 × 10−4 76.58 3.95 × 10−4 166.74 128 × 128 2.3 × 10−3 188.61 2.14 × 10−5 360.3 1.57 × 10−5 768.91

Table : Example 2: CPU-time needed to solve the Monge-Amp` ere equation and the adaptive Error Eadp using the Picard algorithm with the two grids method.

A. Ratnani

IGA – CEMRACS 2015 40/1

40/1

SLIDE 52

Adaptive meshes

Anisotropic case, on a square domain

A. Ratnani

IGA – CEMRACS 2015 41/1

41/1

SLIDE 53

Adaptive meshes

Anisotropic case, on a square domain

A. Ratnani

IGA – CEMRACS 2015 42/1

42/1

SLIDE 54

Adaptive meshes

Anisotropic case, on a square domain

A. Ratnani

IGA – CEMRACS 2015 43/1

43/1

SLIDE 55

Adaptive meshes

Anisotropic case, on a square domain

A. Ratnani

IGA – CEMRACS 2015 44/1

44/1

SLIDE 56

Adaptive meshes

Anisotropic case, on a square domain

A. Ratnani

IGA – CEMRACS 2015 45/1

45/1

SLIDE 57

Adaptive meshes

Example on annulus domain

The annulus domain is contructed using NURBS. (NURBS curves model all conics)

A. Ratnani

IGA – CEMRACS 2015 46/1

46/1

SLIDE 58

B´ ezier techniques in Computational Plasma Physics

Applications

Mesh Generation Kinetic models Particle In Cell Semi-Lagrangian schemes MagnetoHydrodynamics (MHD) Equilibrium Non-linear Reduced MHD Full MHD Tomography

A. Ratnani

IGA – CEMRACS 2015 47/1