Finite Element Method  3D Mesh Generator  Metfem3D
FEM Approach
Space  dependent term
After integration NK(Nu) by parts, one obtains a mixed set of linear and nonlinear equations
∫ 
(∇ N_{b})^{T} k_{+}
(∑_{a} N_{a}n^{+,a})

∇ (∑_{c} N_{c}φ^{c}) dΩ +

∫ 
(∇ N_{b})^{T} D_{+} ∇(∑_{a} N_{a} n^{+, a})dΩ

+ 1/(Δt)

∫ 
N^{b} ∑_{a} N_{a} ( n^{+, a}_{n}  n^{+, a}_{n1}) dΩ
= 0

 ∫ 
(∇ N_{b})^{T}k_{} (∑_{a} N_{a}n^{, a})
∇ (∑_{c} N_{c}φ^{c}) dΩ +

∫ 
(∇ N_{b})^{T}D_{} ∇(∑_{a} N_{a} n^{, a})dΩ

+ 1/(Δt)

∫ 
N^{b} ∑_{a} N_{a} (n^{,a}_{n}  n^{,a}_{n1}) dΩ = 0



∫ 
(∇ N_{b})^{T}ε_{0}ε∇(∑_{a} N_{a}φ^{a}) dΩ

+ z_{+}e

∫ 
N^{b} (∑_{a} N_{a} n^{+,a})dΩ

+ z_{}e

∫ 
N^{b} (∑_{a} N_{a} n^{,a})dΩ = 0

and a corresponding set of boundary terms for φ and n^{i} where i = {+, }
 k_{i}

∫ 
N^{b} ∂/(∂n) {(∑_{a} N_{a} n^{i,a})∇(∑_{c} N_{c}φ^{c})} dΓ

 D_{i}

∫ 
N^{b} ∂/(∂n){∇(∑_{a} N_{a}n^{i,a})} dΓ

+ εε_{0}

∫ 
N^{b} ∂/(∂n){∇(∑_{a} N_{a}φ^{a})} dΓ = 0

where ∂/(∂n) denotes derivative normal to Γ. Presentedabove spatially temporal discretization has been done for the case with β = 1 at each node. If the forced boundary conditions are imposed on Γ_{φ} and Γ_{n}, respectively, then all terms in equation can be neglected by restricting the choice of N_{b} functions to those which equal 0 on Γ. Let us denote integrals from equation as
K_{ac}^{b} =

∫ 
(∇ N_{b})^{T} N_{a}∇ N_{c} dΩ = ∑_{e}

∫ 
(∇ N_{b})^{T} N_{a}∇ N_{c} dΩ^{e} = ∑_{e} K_{ac}^{b, e}

K*_{a}^{b} =

∫ 
(∇ N_{b})^{T}∇ N_{a} dΩ = ∑_{e}

∫ 
(∇ N_{b})^{T}∇ N_{a} dΩ^{e} = ∑_{e} K*_{a}^{b, e}

K_{a}^{b} =

∫ 
N_{a} N^{b} dΩ = ∑_{e}

∫ 
N_{a} N^{b} dΩ^{e} = ∑_{e} K_{a}^{b,e}.

where ∑_{e} with e = 1, …, E denotes sum over elements. In 3D space tetrahedral elements seem to be a natural choice of finite volume elements. Then indices (a,b,c) take four values each (an element has four nodes) from the set of 1, …, M values.
References