MATH 676 Finite element methods in scientific computing Wolfgang - PowerPoint PPT Presentation

MATH 676 – Finite element methods in scientific computing Wolfgang Bangerth, Texas A&M University http://www.dealii.org/ Wolfgang Bangerth

Lecture 33.5: Which quadrature formula to use http://www.dealii.org/ Wolfgang Bangerth

The need for quadrature Recall from lecture 4 and many example programs: We compute A ij =(∇ φ i , ∇ φ j ) F i =(φ i , f ) by mapping back to the reference cell... A ij = (∇ φ i , ∇ φ j ) = ∑ K ∫ K ∇ φ i ( x )⋅∇ φ j ( x ) x ) ̂ x ) ̂ = ∑ K ∫ ̂ ∇ ̂ ∇ ̂ − 1 (̂ − 1 (̂ φ i (̂ φ j (̂ x ) ∣ det J K (̂ K J K x ) ⋅ J K x )∣ ...and quadrature: Q x q ) ̂ x q ) ̂ A ij ≈ ∑ K ∑ q = 1 ∇ ̂ ∇ ̂ − 1 (̂ − 1 ( ̂ ⏟ J K φ i ( ̂ x q ) ⋅ J K φ j ( ̂ x q ) ∣ det J K (̂ x q )∣ w q = :JxW Similarly for the right hand side F . http://www.dealii.org/ Wolfgang Bangerth

The need for quadrature Question: When approximating A ij =(∇ φ i , ∇ φ j ) by Q x q ) ̂ x q ) ̂ A ij ≈ ∑ K ∑ q = 1 ∇ ̂ ∇ ̂ − 1 (̂ − 1 ( ̂ J K φ i ( ̂ x q ) ⋅ J K φ j ( ̂ x q ) ∣ det J (̂ x q )∣ w q how should we choose the points and weights w q ? ̂ x q In other words: Which quadrature rule should we choose? http://www.dealii.org/ Wolfgang Bangerth

Considerations Question: Which quadrature rule should we choose? Q x q ) ̂ x q ) ̂ A ij ≈ ∑ K ∑ q = 1 ∇ ̂ ∇ ̂ − 1 (̂ − 1 ( ̂ J K φ i ( ̂ x q ) ⋅ J K φ j ( ̂ x q ) ∣ det J (̂ x q )∣ w q Goals: ● Efficient: Make Q as small as possible ● Accurate: Do not introduce unnecessary errors About accuracy: In particular, do nothing that affects the convergence order ! http://www.dealii.org/ Wolfgang Bangerth

1d: The matrix Question: Which quadrature rule should we choose? Q x q ) ̂ x q ) ̂ A ij ≈ ∑ K ∑ q = 1 ∇ ̂ ∇ ̂ − 1 (̂ − 1 ( ̂ J K φ i ( ̂ x q ) ⋅ J K φ j ( ̂ x q ) ∣ det J K (̂ x q )∣ w q Consider the 1d case: ● We use an element of polynomial degree k ● We use a linear mapping Then: ● are constant − 1 , det J K J K ,J K ● ̂ ∇ ̂ is a polynomial of degree k-1 φ j ( ̂ x q ) ● The integrand has polynomial degree 2(k-1) http://www.dealii.org/ Wolfgang Bangerth

1d: The matrix Question: Which quadrature rule should we choose? A ij = (∇ φ i , ∇ φ j ) Q x q ) ̂ x q ) ̂ ≈ ∑ K ∑ q = 1 ∇ ̂ ∇ ̂ − 1 (̂ − 1 (̂ φ i (̂ x q ) ⋅ J K φ j (̂ x q ) ∣ det J K ( ̂ x q )∣ w q J K Consider the 1d case: ● The integrand has polynomial degree 2(k-1) ● Gauss quadrature with n points is exact for polynomials up to degree 2n-1 Consequence: We can compute the integral A ij =(∇ φ i , ∇ φ j ) exactly via Gauss quadrature with n=k points! http://www.dealii.org/ Wolfgang Bangerth

1d: The right hand side Question: How about the right hand side? Q F i = (φ i , f ) ≈ ∑ K ∑ q = 1 ̂ φ i ( ̂ x q ) f ( x q ) ∣ det J K (̂ x q )∣ w q Consider the 1d case: ● We use an element of polynomial degree k ● We use a linear mapping Then: ● is constant det J K ● ̂ ∇ ̂ is a polynomial of degree k φ j ( ̂ x q ) ● is not in general a polynomial f ( x ) ● The integrand is not polynomial http://www.dealii.org/ Wolfgang Bangerth

1d: The right hand side Question: What to do here? Q F i = (φ i , f ) ≈ ∑ K ∑ q = 1 ̂ φ i ( ̂ x q ) f ( x q ) ∣ det J K (̂ x q )∣ w q Consider Gauss integration with n points: ● Integrates polynomials of degree 2n-1 exactly ● For general f(x) essentially integrates Q F i = (φ i , f ) ≈ ∑ K ∑ q = 1 ̂ φ i (̂ x q ) f ( x q ) ∣ det J K (̂ x q )∣ w q ≈ (φ i , I 2n − k f ) where I 2n-k f = f at the n quadrature points + n-k others http://www.dealii.org/ Wolfgang Bangerth

1d: The right hand side Consider Gauss integration with n points: ● Integrates polynomials of degree 2n-1 exactly ● For general f(x) essentially integrates Q F i = (φ i , f ) ≈ ∑ K ∑ q = 1 ̂ φ i ( ̂ x q ) f ( x q ) ∣ det J K (̂ x q )∣ w q = ∑ K ∫ K I 2n (φ i f ) ≈ (φ i , I 2n − k f ) where I 2n-k f = f at the n quadrature points + n-k others on every cell Consequence: Inexact integration is equivalent to approximating the solution of a slightly perturbed problem! http://www.dealii.org/ Wolfgang Bangerth

1d: The right hand side Consider the original and perturbed problems: u =̃ −Δ u = f in Ω −Δ ̃ f in Ω u = 0 on ∂Ω u = 0 on ∂Ω Consequence: ∥ u −̃ u h ∥ H 1 ≤ ∥ u −̃ u ∥ H +∥̃ u −̃ u h ∥ H ⏟ ⏟ 1 1 ≤ C 1 ∥ f −̃ 2n − k + 1 ∥ f ∥ H 2n − k k ∥̃ f ∥ H − 1 ≤ C 2 h ≤ C 3 h u ∥ H k We want the first term to be at least as good as the second. We need to choose n=k . http://www.dealii.org/ Wolfgang Bangerth

Higher dimensions: The matrix Question: Which quadrature rule should we choose? Q x q ) ̂ x q ) ̂ A ij ≈ ∑ K ∑ q = 1 ∇ ̂ ∇ ̂ − 1 (̂ − 1 ( ̂ J K φ i ( ̂ x q ) ⋅ J K φ j ( ̂ x q ) ∣ det J K (̂ x q )∣ w q Consider the higher dimensional case: ● Use an element of polynomial degree k in each direction ● Use a d- linear mapping Then: ● are polynomials of degree k, k d J K , det J K ● is a rational function − 1 J K ● ̂ ∇ ̂ is a polynomial of degree dk-1 φ j ( ̂ x q ) ● The integrand is rational ● For linear mappings, it is of degree 2(dk-1) http://www.dealii.org/ Wolfgang Bangerth

Higher dimensions: The matrix Question: Which quadrature rule should we choose? Q x q ) ̂ x q ) ̂ A ij ≈ ∑ K ∑ q = 1 ∇ ̂ ∇ ̂ − 1 (̂ − 1 ( ̂ J K φ i ( ̂ x q ) ⋅ J K φ j ( ̂ x q ) ∣ det J K (̂ x q )∣ w q Consider the higher dimensional case: ● The integrand is rational ● For linear mappings, it is of degree 2(dk-1) ● Gauss quadrature with n points per direction is exact for degree 2n-1 in each variable Nevertheless, using the tensor product structure: We need to use Gauss quadrature with n=k+1 points per direction. http://www.dealii.org/ Wolfgang Bangerth

Higher dimensions: The right hand side Question: Which quadrature rule should we choose? Q F i = (φ i , f ) ≈ ∑ K ∑ q = 1 ̂ φ i ( ̂ x q ) f ( x q ) ∣ det J K (̂ x q )∣ w q Similar considerations can be applied: We need to use Gauss quadrature with n=k+1 points per direction. http://www.dealii.org/ Wolfgang Bangerth

Summary As a general rule of thumb: ● Gauss quadrature with n=k+1 points per direction is sufficient – for the Laplace matrix A ij = (∇ φ i , ∇ φ j ) – for the mass matrix M ij = (φ i , φ j ) – for the right hand side F i = (φ i , f ) ● It is generally also sufficient with variable coefficients: A ij = ( a ( x ) ∇ φ i , ∇ φ j ) With n=k+1 , the quadrature error does not dominate the overall error (if a(x) is smooth). http://www.dealii.org/ Wolfgang Bangerth

Non-smooth coefficients What to with non-smooth terms? For example A ij = ( a ( x ) ∇ φ i , ∇ φ j ) F i = (φ i , f ) where a(x) or f(x) are discontinuous. Recall: Quadrature is equivalent to exact integration with an interpolated coefficient. For discontinuous functions, interpolation does not help very much: Quadrature produces large errors. http://www.dealii.org/ Wolfgang Bangerth

Non-smooth coefficients What to with non-smooth terms? For example A ij = ( a ( x ) ∇ φ i , ∇ φ j ) F i = (φ i , f ) where a(x) or f(x) are discontinuous. Before: ∥ u −̃ u h ∥ H 1 ≤ ∥ u −̃ u ∥ H +∥̃ u −̃ u h ∥ H ⏟ ⏟ 1 1 ≤ C 1 ∥ f −̃ 2n − k + 1 ∥ f ∥ H 2n − k k ∥̃ f ∥ H − 1 ≤ C 2 h ≤ C 3 h u ∥ H k Now: The interpolation step fails! We may only get ∥ u −̃ u h ∥ H 1 ≤ ∥ u −̃ u ∥ H +∥̃ u −̃ u h ∥ H ⏟ ⏟ 1 1 ≤ C 1 ∥ f −̃ s + 1 ∥ f ∥ H s k ∥̃ f ∥ H − 1 ≤ C 2 h ≤ C 3 h u ∥ H k http://www.dealii.org/ Wolfgang Bangerth

Non-smooth coefficients What to with non-smooth terms? For example A ij = ( a ( x ) ∇ φ i , ∇ φ j ) F i = (φ i , f ) where a(x) or f(x) are discontinuous. Now: ∥ u −̃ u h ∥ H 1 ≤ ∥ u −̃ u ∥ H +∥̃ u −̃ u h ∥ H ⏟ ⏟ 1 1 ≤ C 1 ∥ f −̃ s + 1 ∥ f ∥ H s k ∥̃ f ∥ H − 1 ≤ C 2 h ≤ C 3 h u ∥ H k Solution: Subdivide the cell into L pieces so that C 2 ( L ) s + 1 h s ≈ C 3 h k ∥̃ ∥ f ∥ H u ∥ H k This is what the QIterated class does. http://www.dealii.org/ Wolfgang Bangerth

Special purpose quadratures There are situations where we want quadrature rules other than Gauss: ● To affect stability properties of a discretization – Underintegration for nearly incompressible elasticity – Special purpose quadrature for mixed problems ● To improve sparsity of matrices – Make some terms zero – Make a matrix diagonal http://www.dealii.org/ Wolfgang Bangerth

Sparsifying matrices Using the trapezoidal rule for the Laplace matrix: Assume: ● Uniform mesh with square cells ● Q 1 element with shape functions φ 1 =( 1 −̂ x )( 1 − ̂ y ) , φ 2 =̂ x ( 1 − ̂ y ) , φ 3 =( 1 −̂ x ) ̂ y, φ 4 =̂ x ̂ y and gradients ∇ φ 1 = ( x ) ) , ̂ ∇ φ 2 = ( x ) −( 1 − ̂ y ) ( 1 − ̂ y ) ̂ −( 1 −̂ −̂ ∇ φ 3 = ( x ) , ̂ ∇ φ 4 = ( x ) − ̂ ̂ y y ̂ 1 −̂ ̂ ● Trapezoidal rule with integration points at the vertices http://www.dealii.org/ Wolfgang Bangerth