feat(equations): added gravity force equation. This computes the point-wise gravitational force (which is the residual as generally it should be minimzed). Note that I let the gradient = the negative gravitational acceleration

src/include/serif/equations/gravity_force.hpp

@@ -0,0 +1,99 @@
+#pragma once
+
+/**
+ * @file gravity_force.hpp
+ * @author Emily M. Boudreaux
+ *
+ * @brief File containing the gravitational force contribution to the mechanical balance equation (displacement_force.hpp).
+ */
+
+#include <concepts>
+#include <string_view>
+
+#include "serif/discretization/blocks/fields.hpp"
+#include "serif/equations/contributions.hpp"
+#include "serif/equations/point.hpp"
+
+namespace serif::equations {
+
+    /**
+     * @class GravityForce
+     *
+     * @brief The gravitational force contribution to the mechanical balance equation.
+     *
+     * @contribution Residuals to the global displacement vector supported on the stellar domain.
+     *
+     * @note The gravitational force is given by \f$ \mathbf{F} = \rho \mathbf{g} \f$, where \f$ \rho \f$ is the mass density and \f$ \mathbf{g} \f$ is the gravitational acceleration. The residual is given by the difference between the model gravitational force and the gravitational force given by the equation of state. The partial derivatives implemented are:
+     * - \f$ \frac{\partial F}{\partial \rho} = g \f$
+     * - \f$ \frac{\partial F}{\partial g} = \rho I \f$
+     *
+     * @note For some reason which only past me can fathom I adopted \f$ \mathbf{g} = \nabla \Phi \f$ rather than \f$ \mathbf{g} = -\nabla \Phi \f$. The implication of this is that acceleration is \f$ -\mathbf{g} \f$. We may want to change this but we will need to be careful.
+     *
+     * @note The reason is actually completely fathomable: I just didn't think. Carry on now.
+     *
+     * @note You know I've written, over the last year about 60k lines of code for this project, much of it deleted and rewritten multiple time. I am
+     * going through now and adding comments to this version which I am finally happy with. In doing so I suspect that a reader can see
+     * the slow(ish?) deterioration of my mental state as these comments get more and more conversational. I wonder if anyone will have anything to say about this
+     * or if it will simply become some marginal curiosity buried in the docs. I suspect the former, though I am not so selfless as to not hope for the latter.
+     */
+    struct GravityForce : Contribution<discretization::forms::Operand<discretization::blocks::displacement::vector>, discretization::domain::StellarDomains> {
+        using Density = discretization::blocks::density::mass;                  ///< The density type is the model's mass density,
+        using Gradient = discretization::blocks::gravity::gradient;             ///< The gradient type is the model's negative gravitational acceleration (see class comment for why this is negative).
+
+        using DensityValue = discretization::forms::Operand<Density>;           ///< The density value type is a density Operand
+        using GradientValue = discretization::forms::Operand<Gradient>;         ///< The gradient value type is a gradient Operand
+
+        using Operands = utils::types::TypeList<DensityValue, GradientValue>;   ///< The operands for this equation are the density and gradient operands
+        using Result = discretization::forms::VectorValue;                      ///< This equation's results are vector valued.
+
+        static constexpr std::string_view name = "gravitational force";         ///< Compile time name for the gravitational force equation, used for debugging and logging.
+
+        /**
+         * @brief Evaluate the gravitational force equation at a given point.
+         *
+         * @param point The point at which to evaluate the gravitational force equation.
+         * @return A PointLinearization object containing the residual and field derivatives for the gravitational force equation at the given point.
+         *
+         * This function evaluates
+         * - \f$ F = \rho g \f$
+         * - \f$ \frac{\partial F}{\partial \rho} = g \f$
+         * - \f$ \frac{\partial F}{\partial g} = \rho I \f$
+         */
+        [[nodiscard]] static PointLinearization<Operands, Result> evaluate(
+            const EquationPoint<Operands>& point
+        ) {
+            const double density = point.fields.get<Density>();
+            const auto& gradient = point.fields.get<Gradient>();
+
+            PointLinearization<Operands, Result> result;
+
+            for (int axis = 0; axis < 3; ++axis) {
+                result.residual[axis] = density * gradient[axis];
+                result.field_derivatives.get<Density>()[axis] = gradient[axis];
+                result.field_derivatives.get<Gradient>()[axis][axis] = density;
+            }
+
+            return result;
+        }
+
+        /**
+         * @brief Return the additional quadrature order required for the gravitational force equation. We add one
+         * to the registered order of the gradient functional space to satisfy the InfSup condition for the density operand (which is in L2)
+         *
+         * @param operand The operand for which to determine the additional quadrature order.
+         * @param position_order The order of the position operand, which is not used in this implementation.
+         * @return The additional quadrature order required for the given operand.
+         *
+         * This function returns the additional quadrature order required for the gravitational force equation. For the density operand, it returns the order of the gradient functional space plus one. For other operands, it returns the order of the density functional space.
+         */
+        template <typename OperandT>
+        [[nodiscard]] static int additional_order(OperandT /*operand_type*/, int /*position_order*/) {
+            if constexpr (std::same_as<OperandT, DensityValue>) {
+                // The registered RT family of order p has vector components of degree up to p+1.
+                return Gradient::FunctionalSpace::order + 1;
+            } else {
+                return Density::FunctionalSpace::order;
+            }
+        }
+    };
+} // namespace serif::equations