Appendix A. Some basic computations between reference and real elements

Volume integral

One has

∫_T f(x) dx = ∫_T' f'(x') |vol((∂τ(x'))/(∂x'₀) ;(∂τ(x'))/(∂x'₁); ...; (∂τ(x'))/(∂x'_P-1 ))| dx'.

Denoting J_τ(x') the jacobian

J_τ(x') := |vol((∂τ(x'))/(∂x'₀) ;(∂τ(x'))/(∂x'₁); ...; (∂τ(x'))/(∂x'_P-1 ))| = (det(K(x')^T K(x')))^1/2,

one finally has

∫_T f(x) dx = ∫_T' f'(x') J_τ(x')dx'.

When P = N, the expression of the jacobian reduces to J_τ(x') = |det(K(x'))|.

Surface integral

With Γ a part of the boundary of T a real element and Γ' the corresponding boundary on the reference element T', one has

∫_Γ f(x) dσ= ∫_Γ' f'(x') ||B(x')n'|| J_τ(x') dσ',

where n' is the unit normal to T' on Γ'. In a same way

∫_Γ F(x).n dσ= ∫_Γ' F'(x').(B(x')n') J_τ(x') dσ'.

Second derivative computation

Denoting

∇² f = ((∂² f)/(∂x_i ∂x_j))_ij,

the N ×N matrix and

X'(x') = ∑_{k = 0}^N-1 ∇'² τ_k(x') (∂f)/(∂x_k)(x) = ∑_{k =
0}^N-1 ∑_{i = 0}^P-1 ∇'² τ_k(x') B_ki (∂f')/(∂x'_i)(x'),

the P ×P matrix, then

∇'² f'(x') = X'(x') + K(x')^T ∇² f(x) K(x'),

and thus

∇² f(x) = B(x') (∇'² f'(x') - X'(x')) B(x')^T.

In order to have uniform methods for the computation of elementary matrices, the Hessian is computed as a column vector H f whose components are (∂² f)/(∂x²₀), (∂² f)/(∂x₁ ∂x₀), ... (∂² f)/(∂x²_N-1). Then, with B₂ the P² ×P matrix defined as

(B₂(x'))_ij = ∑_{k =
0}^N-1 (∂² τ_k(x'))/(∂x'_{i / P} ∂x'_{i
mod P} ) B_kj(x'),

and B₃ the N² ×P² matrix defined as

(B₃(x'))_ij = B_{i / N, j /
P}(x') B_{i mod N, j mod
P}(x'),

one has

H f(x) = B₃(x') (H'f'(x') - B₂(x')∇'f'(x')).

Example of elementary matrix

Assume one needs to compute the elementary "matrix":

t(i₀, i₁, ..., i₇) = ∫_T φ_i₁^i₀ ∂_i₄ φ_i₃^i₂ ∂²_{i₇ / P, i₇
mod P} φ_i₆^i₅ dx,

The computations to be made on the reference elements are

t'₀(i₀, i₁, ..., i₇) = ∫_T' φ'_i₁^i₀ ∂_i₄ φ'_i₃^i₂ ∂²_{i₇ / P, i₇
mod P} φ'_i₆^i₅ J(x') dx',

and

t'₁(i₀, i₁, ..., i₇) = ∫_T' φ'_i₁^i₀ ∂_i₄ φ'_i₃^i₂ ∂_i₇ φ'_i₆^i₅ J(x') dx',

Those two tensor can be computed once on the whole reference element if the geometric transformation is linear (because J(x') is constant). If the geometric transformation is non-linear, what has to be stored is the value on each integration point. To compute the integral on the real element a certain number of reductions have to be made:

Concerning the first term (φ_i₁^i₀) nothing.
Concerning the second term (∂_i₄ φ_i₃^i₂) a reduction with respect to i₄ with the matrix B.
Concerning the third term (∂²_{i₇ / P, i₇
mod P} φ_i₆^i₅) a reduction of t'₀ with respect to i₇ with the matrix B₃ and a reduction of t'₁ with respect also to i₇ with the matrix B₃B₂

The reductions are to be made on each integration point if the geometric transformation is non-linear. Once those reductions are done, an addition of all the tensor resulting of those reductions is made (with a factor equal to the load of each integration point if the geometric transformation is non-linear).

If the finite element is non-τ-equivalent, a supplementary reduction of the resulting tensor with the matrix M has to be made.

Getfem++ project

Appendix A. Some basic computations between reference and real elements

Volume integral

Surface integral

Derivative computation

Second derivative computation

Example of elementary matrix