ENH: overhaul of the adjoint optimisation library

Parts of the adjoint optimisation library were re-designed to generalise
the way sensitivity derivatives (SDs) are computed and to allow easier
extension to primal problems other than the ones governed by
incompressible flows. In specific:
- the adjoint solver now holds virtual functions returning the part of
  SDs that depends only on the primal and the adjoint fields.
- a new class named designVariables was introduced which, apart from
  defining the design variables of the optimisation problem and
  providing hooks for updating them in an optimisation loop, provides
  the part of the SDs that affects directly the flow residuals (e.g.
  geometric variations in shape optimisation, derivatives of source
  terms in topology optimisation, etc). The final assembly of the SDs
  happens here, with the updated sensitivity class acting as an
  intermediate.

With the new structure, when the primal problem changes (for instance,
passive scalars are included), the same design variables and sensitivity
classes can be re-used for all physics, with additional contributions to
the SDs being limited (and contained) to the new adjoint solver to be
implemented. The old code structure would require new SD classes for
each additional primal problem.

As a side-effect, setting up a case has arguably become a bit easier and
more intuitive.

Additional changes include:
---------------------------

- Changes in the formulation and computation of shape sensitivity derivatives
  using the E-SI approach. The latter is now derived directly from the
  FI approach, with proper discretization for the terms and boundary
  conditions that emerge from applying the Gauss divergence theorem used
  to transition from FI to E-SI. When E-SI and FI are based on the same
  Laplace grid displacement model, they are now numerically equivalent
  (the previous formulation proved the theoretical equivalence of the
  two approaches but numerical results could differ, depending on the
  case).
- Sensitivity maps at faces are now computed based (and are deriving
  from) sensitivity maps at points, with a constistent point-to-face
  interpolation (requires the differentiation of volPointInterpolation).
- The objective class now allocates only the member pointers that
  correspond to the non-zero derivatives of the objective w.r.t. the
  flow and geometric quantities, leading to a reduced memory footprint.
  Additionally, contributions from volume-based objectives to the
  adjoint equations have been re-worked, removing the need for
  objectiveManager to be virtual.
- In constrained optimisation, an adjoint solver needs to be present for
  each constraint function. For geometric constraints though, no adjoint
  equations need to solved. This is now accounted for through the null
  adjoint solver and the geometric objectives which do not allocate
  adjoint fields for this kind of constraints, reducing memory
  requirements and file clutter.
- Refactoring of the updateMethod to collaborate with the new
  designVariables. Additionally, all updateMethods can now read and
  write restart data in binary, facilitating exact continuation.
  Furthermore, code shared by various quasi-Newton methods (BFGS, DBFGS,
  LBFGS, SR1) has been organised in the namesake class. Over and above,
  an SQP variant capable of tackling inequality constraints has been
  added (ISQP, with I indicating that the QP problem in the presence of
  inequality constraints is solved through an interior point method).
  Inequality constraints can be one-sided (constraint < upper-value)
  or double-sided (lower-value < constraint < upper-value).
- Bounds can now be defined for the design variables.
  For volumetricBSplines in specific, these can be computed as the
  mid-points of the control points and their neighbouring ones. This
  usually leads to better-defined optimisation problems and reduces the
  chances of an invalid mesh during optimisation.
- Convergence criteria can now be defined for the optimisation loop
  which will stop if the relative objective function reduction over
  the last objective value is lower than a given threshold and
  constraints are satisfied within a give tolerance. If no criteria are
  defined, the optimisation will run for the max. given number of cycles
  provided in controlDict.
- Added a new grid displacement method based on the p-Laplacian
  equation, which seems to outperform other PDE-based approaches.

TUT: updated the shape optimisation tutorials and added a new one
showcasing the use of double-sided constraints, ISQP, applying
no-overlapping constraints to volumetric B-Splines control points
and defining convergence criteria for the optimisation loop.