I am used to thinking of finite-differences as a special case of finite-elements, on a very constrained grid. So what are the conditions on how to choose between Finite Difference Method (FDM) and Finite Element Method (FEM) as a numerical method?

On the side of Finite Difference Method (FDM), one may count that they are conceptually simpler and easier to implement than Finite Element Method (FEM). FEM have the benefit of being very flexible, e.g., the grids may be very non-uniform and the domains may have arbitrary shape.

The only example I know where FDM has turned out superior to FEM is in
Celia, Bouloutas, Zarba, where the benefit is due to the FD method using a different discretization of time derivative, which, however, could be fixed for the finite element method.

It is possible to write most specific finite difference methods as Petrov-Galerkin finite element methods with some choice of local reconstruction and quadrature, and most finite element methods can also be shown to be algebraically equivalent to some finite difference method. Therefore, we should choose a method based on which analysis framework we want to use, which terminology we like, which system for extensibility we like, and how we would like to structure software. The following generalizations hold true in the vast majority of variations in practical use, but many points can be circumvented.

Finite Difference

Pros


efficient quadrature-free implementation
aspect ratio independence and local conservation for certain schemes (e.g. MAC for incompressible flow)
robust nonlinear methods for transport (e.g. ENO/WENO)
M-matrix for some problems
discrete maximum principle for some problems (e.g. mimetic finite differences)
diagonal (usually identity) mass matrix
inexpensive nodal residual permits efficient nonlinear multigrid (FAS)
cell-wise Vanka smoothers give efficient matrix-free smoothers for incompressible flow


Cons


more difficult to implement "physics"
staggered grids are sometimes quite technical
higher than second order on unstructured grids is difficult
no Galerkin orthogonality, so convergence may be more difficult to prove
not a Galerkin method, so discretization and adjoints do not commute (relevant to optimization and inverse problems)
self-adjoint continuum problems often yield non-symmetric matrices
solution is only defined pointwise, so reconstruction at arbitrary locations is not uniquely defined
boundary conditions tend to be complicated to implement
discontinuous coefficients usually make the methods first order
stencil grows if physics includes "cross terms"


Finite Element

Pros


Galerkin orthogonality (discrete solution to coercive problems is within a constant of the best solution in the space)
simple geometric flexibility
discontinuous Galerkin offers robust transport algorithm, arbitrary order on unstructured grids
cellwise entropy inequality guaranteeing $L^2$ stability holds independent of mesh, dimension, order of accuracy, and presence of discontinuous solutions, without needing nonlinear limiters
easy of implementing boundary conditions
can choose conservation statement by choosing test space
discretization and adjoints commute (for Galerkin methods)
elegant foundation in functional analysis
at high order, local kernels can exploit tensor product structure that is missing with FD
Lobatto quadrature can make methods energy-conserving (assuming a symplectic time integrator)
high order accuracy even with discontinuous coefficients, as long as you can align to boundaries
discontinuous coefficients inside elements can be accommodated with XFEM
easy to handle multiple inf-sup conditions


Cons


many elements have trouble at high aspect ratio
continuous FEM has trouble with transport (SUPG is diffusive and oscillatory)
DG usually has more degrees of freedom for same accuracy (though HDG is much better)
continuous FEM does not provide cheap nodal problems, so nonlinear smoothers have much poorer constants
usually more nonzeros in assembled matrices
have to choose between consistent mass matrix (some nice properties, but has full inverse, thus requiring an implicit solve per time step) and lumped mass matrix.

