## Recipe for Convection & Diffusion in 3-D

A conventional Finite Element method with **linear triangles** (2-D) or
**linear tetrahedra** (3-D)

is assumed in the first place. Then do the following:
- Calculate the element matrix for Diffusion, in the traditional manner.
- Any finite element matrix for diffusion can always be considered as a
superposition of 1-D resistor elements,

according to
[ 1 ] for 2-D and
[ 2 ] for 3-D.
Split the element into these superpositions and call, each of them: e[ie,je]

- Eventually divide these 2 x 2 matrix entries by the conductivity:
e[ie,je] := e[ie,je] / λ ;

Then the resistor matrices have become dependent on Geometry only.
- Determine the coordinates (x,y,z) of the end points of all straight line
segments

joining any two vertices of the finite element shape.

These line segments are interpreted as vectors with components (dx, dy, dz).

The directions in which the line-segments are pointed should be calculated
whith care:
dx := x[je] - x[ie] ; dy := y[je] - y[ie] ; dz := z[je] - z[ie];

- Determine the components (v.x, v.y, v.z) of the velocity vectors, in the
middle of each line segment.
- Determine the inner products of the velocity vectors and the corresponding
line segments:
produkt := v.x*dx + v.y*dy + v.z*dz;

- Determine (absolute values of) the local Péclet numbers:
Péclet := abs( ρ * C * produkt / λ );

- At last, account for Convection And Diffusion contributions in the finite
element matrices:
e[ie,je] := e[ie,je] * (General(Péclet) + max(0,- Péclet));

Here *General* is the function as given by A(|P|) in table 5.2 from the
book by S.V. Patankar: *Numerical Heat Transfer and Fluid Flow*.

The assembly is such that the case of pure Diffusion is recovered for
(Péclet = 0). Upwind schemes for Convection will be recovered as well,

as soon as local Péclet numbers assume high absolute values.
In case of **Anisotropic** Conduction, it's almost impossible to implement
something else than the Upwind scheme.

The reason is that the Conductive
Tensor ( Λ ) being singular is far from an exception to the rule.

e[ie,je] := e[ie,je] * ( Λ + max(0,- abs(ρ * C * produkt) )

A live example is found in the *Physics unit* of the 3-D
ZonWind program.
### ( Disclaimer: the above has not been put into practice, not yet ... )