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:
  1. Calculate the element matrix for Diffusion, in the traditional manner.
  2. 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]
  3. 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.
  4. 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];
    
  5. Determine the components (v.x, v.y, v.z) of the velocity vectors, in the middle of each line segment.
  6. Determine the inner products of the velocity vectors and the corresponding line segments:
          produkt := v.x*dx + v.y*dy + v.z*dz;
    
  7. Determine (absolute values of) the local Péclet numbers:
          Péclet := abs( ρ * C * produkt / λ );
    
  8. 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 ... )