implicit BCs and Laplacian operator

include/NeoFOAM/core/primitives/vector.hpp

@@ -145,9 +145,9 @@ Vector operator*(const scalar& sclr, Vector rhs)
 }
 
 KOKKOS_INLINE_FUNCTION
-Vector operator&(const Vector& lhs, Vector rhs)
+scalar operator&(const Vector& lhs, Vector rhs)
 {
-    return Vector(rhs[0] * lhs[0], rhs[1] * lhs[1], rhs[2] * lhs[2]);
+    return lhs[0] * rhs[0] + lhs[1] * rhs[1] + lhs[2] * rhs[2];
 }
 
 KOKKOS_INLINE_FUNCTION

include/NeoFOAM/core/tokenList.hpp

@@ -11,6 +11,11 @@ namespace NeoFOAM
 {
 
 
+void logOutRange(
+    const std::out_of_range& e, const std::size_t& key, const std::vector<std::any>& data
+);
+
+
 /**
  * @class TokenList
  * @brief A class representing a list of tokens.
@@ -22,13 +27,15 @@ class TokenList
 {
 public:
 
-    TokenList() = default;
+    TokenList();
+
+    TokenList(const TokenList&);
 
     /**
      * @brief Construct a TokenList object from a vector of std::any.
      * @param data A vector of std::any.
      */
-    TokenList(const std::vector<std::any>& data);
+    TokenList(const std::vector<std::any>& data, size_t nextIndex = 0);
 
     /**
      * @brief Construct a TokenList object from an initializer list of std::any.
@@ -94,8 +101,27 @@ class TokenList
             logBadAnyCast<ReturnType>(e, idx, data_);
             throw e;
         }
+        catch (const std::out_of_range& e)
+        {
+            logOutRange(e, idx, data_);
+            throw e;
+        }
     }
 
+    /**
+     * @brief Retrieves the value associated with the nextIndex, casting it to
+     * the specified type.
+     * @tparam T The type to cast the value to.
+     * @return A reference to the value associated with the index, casted to
+     * type T.
+     */
+    template<typename ReturnType>
+    ReturnType& next()
+    {
+        ReturnType& retValue = get<ReturnType&>(nextIndex_);
+        nextIndex_++;
+        return retValue;
+    }
 
     /**
      * @brief Retrieves the value associated with the given index, casting it to
@@ -117,6 +143,26 @@ class TokenList
             logBadAnyCast<ReturnType>(e, idx, data_);
             throw e;
         }
+        catch (const std::out_of_range& e)
+        {
+            logOutRange(e, idx, data_);
+            throw e;
+        }
+    }
+
+    /**
+     * @brief Retrieves the value associated with the nextIndex, casting it to
+     * the specified type.
+     * @tparam T The type to cast the value to.
+     * @return A const reference to the value associated with the index, casted to
+     * type T.
+     */
+    template<typename ReturnType>
+    const ReturnType& next() const
+    {
+        const ReturnType& retValue = get<ReturnType>(nextIndex_);
+        nextIndex_++;
+        return retValue;
     }
 
     /**
@@ -132,6 +178,7 @@ class TokenList
 private:
 
     std::vector<std::any> data_;
+    mutable size_t nextIndex_;
 };
 
 } // namespace NeoFOAM

include/NeoFOAM/dsl/expression.hpp

@@ -195,6 +195,15 @@ template<typename leftOperator, typename rightOperator>
     return lhs;
 }
 
+template<typename leftOperator, typename rightOperator>
+[[nodiscard]] inline Expression operator-(leftOperator lhs, rightOperator rhs)
+{
+    Expression expr(lhs.exec());
+    expr.addOperator(lhs);
+    expr.addOperator(Coeff(-1) * rhs);
+    return expr;
+}
+
 [[nodiscard]] inline Expression operator-(const SpatialOperator& lhs, const SpatialOperator& rhs)
 {
     Expression expr(lhs.exec());

include/NeoFOAM/dsl/implicit.hpp

@@ -35,4 +35,10 @@ div(fvcc::SurfaceField<NeoFOAM::scalar>& faceFlux, fvcc::VolumeField<NeoFOAM::sc
     return SpatialOperator(fvcc::DivOperator(dsl::Operator::Type::Implicit, faceFlux, phi));
 }
 
+SpatialOperator
+laplacian(fvcc::SurfaceField<NeoFOAM::scalar>& gamma, fvcc::VolumeField<NeoFOAM::scalar>& phi)
+{
+    return SpatialOperator(fvcc::LaplacianOperator(dsl::Operator::Type::Implicit, gamma, phi));
+}
+
 } // namespace NeoFOAM

include/NeoFOAM/dsl/operator.hpp

@@ -16,4 +16,49 @@ class Operator
     };
 };
 
+
+/* @class OperatorMixin
+ * @brief A mixin class to simplify implementations of concrete operators
+ * in NeoFOAMs dsl
+ *
+ * @ingroup dsl
+ */
+template<typename FieldType>
+class OperatorMixin
+{
+
+public:
+
+    OperatorMixin(const Executor exec, const Coeff& coeffs, FieldType& field, Operator::Type type)
+        : exec_(exec), coeffs_(coeffs), field_(field), type_(type) {};
+
+
+    Operator::Type getType() const { return type_; }
+
+    virtual ~OperatorMixin() = default;
+
+    virtual const Executor& exec() const final { return exec_; }
+
+    Coeff& getCoefficient() { return coeffs_; }
+
+    const Coeff& getCoefficient() const { return coeffs_; }
+
+    FieldType& getField() { return field_; }
+
+    const FieldType& getField() const { return field_; }
+
+    /* @brief Given an input this function reads required coeffs */
+    void build([[maybe_unused]] const Input& input) {}
+
+protected:
+
+    const Executor exec_; //!< Executor associated with the field. (CPU, GPU, openMP, etc.)
+
+    Coeff coeffs_;
+
+    FieldType& field_;
+
+    Operator::Type type_;
+};
+
 } // namespace NeoFOAM::dsl

include/NeoFOAM/dsl/spatialOperator.hpp

@@ -58,6 +58,8 @@ class SpatialOperator
 
     SpatialOperator(SpatialOperator&& eqnOperator);
 
+    SpatialOperator& operator=(const SpatialOperator& eqnOperator);
+
     void explicitOperation(Field<scalar>& source);
 
     void implicitOperation(la::LinearSystem<scalar, localIdx>& ls);
@@ -213,46 +215,4 @@ SpatialOperator operator*([[maybe_unused]] CoeffFunction coeffFunc, const Spatia
     return result;
 }
 
-/* @class OperatorMixin
- * @brief A mixin class to simplify implementations of concrete operators
- * in NeoFOAMs dsl
- *
- * @ingroup dsl
- */
-template<typename FieldType>
-class OperatorMixin
-{
-
-public:
-
-    OperatorMixin(const Executor exec, FieldType& field, Operator::Type type)
-        : exec_(exec), coeffs_(), field_(field), type_(type) {};
-
-    Operator::Type getType() const { return type_; }
-
-    virtual ~OperatorMixin() = default;
-
-    virtual const Executor& exec() const final { return exec_; }
-
-    Coeff& getCoefficient() { return coeffs_; }
-
-    const Coeff& getCoefficient() const { return coeffs_; }
-
-    FieldType& getField() { return field_; }
-
-    const FieldType& getField() const { return field_; }
-
-    /* @brief Given an input this function reads required coeffs */
-    void build([[maybe_unused]] const Input& input) {}
-
-protected:
-
-    const Executor exec_; //!< Executor associated with the field. (CPU, GPU, openMP, etc.)
-
-    Coeff coeffs_;
-
-    FieldType& field_;
-
-    Operator::Type type_;
-};
 } // namespace NeoFOAM::dsl

include/NeoFOAM/fields/boundaryFields.hpp

@@ -90,7 +90,7 @@ class BoundaryFields
      * @brief Get the view storing the fraction of the boundary value.
      * @return The view storing the fraction of the boundary value.
      */
-    NeoFOAM::Field<scalar>& valueFraction() { return refValue_; }
+    NeoFOAM::Field<scalar>& valueFraction() { return valueFraction_; }
 
     /** @copydoc BoundaryFields::refGrad()*/
     const NeoFOAM::Field<T>& refGrad() const { return refGrad_; }

include/NeoFOAM/finiteVolume/cellCentred.hpp

@@ -12,10 +12,18 @@
 #include "cellCentred/fields/surfaceField.hpp"
 #include "cellCentred/fields/volumeField.hpp"
 
+#include "cellCentred/operators/linearSystem.hpp"
 #include "cellCentred/operators/divOperator.hpp"
 #include "cellCentred/operators/gaussGreenDiv.hpp"
+
 #include "cellCentred/operators/ddtOperator.hpp"
+
 #include "cellCentred/operators/sourceTerm.hpp"
 
+#include "cellCentred/operators/laplacianOperator.hpp"
+#include "cellCentred/operators/gaussGreenLaplacian.hpp"
+
 #include "cellCentred/interpolation/linear.hpp"
 #include "cellCentred/interpolation/upwind.hpp"
+
+#include "cellCentred/faceNormalGradient/uncorrected.hpp"

include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/fixedGradient.hpp

@@ -26,10 +26,16 @@ void setGradientValue(
 )
 {
     const auto iField = domainField.internalField().span();
-    auto refGradient = domainField.boundaryField().refGrad().span();
-    auto value = domainField.boundaryField().value().span();
-    auto faceCells = mesh.boundaryMesh().faceCells(static_cast<localIdx>(patchID));
-    auto deltaCoeffs = mesh.boundaryMesh().deltaCoeffs(static_cast<localIdx>(patchID));
+
+    auto [refGradient, value, valueFraction, refValue] = spans(
+        domainField.boundaryField().refGrad(),
+        domainField.boundaryField().value(),
+        domainField.boundaryField().valueFraction(),
+        domainField.boundaryField().refValue()
+    );
+
+    auto faceCells = mesh.boundaryMesh().faceCells();
+    auto deltaCoeffs = mesh.boundaryMesh().deltaCoeffs();
 
     NeoFOAM::parallelFor(
         domainField.exec(),
@@ -39,6 +45,8 @@ void setGradientValue(
             // operator / is not defined for all ValueTypes
             value[i] =
                 iField[static_cast<size_t>(faceCells[i])] + fixedGradient * (1 / deltaCoeffs[i]);
+            valueFraction[i] = 0.0;      // only used refGrad
+            refValue[i] = fixedGradient; // not used
         }
     );
 }

include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/fixedValue.hpp

@@ -22,15 +22,21 @@ void setFixedValue(
     DomainField<ValueType>& domainField, std::pair<size_t, size_t> range, ValueType fixedValue
 )
 {
-    auto refValue = domainField.boundaryField().refValue().span();
-    auto value = domainField.boundaryField().value().span();
+    auto [refGradient, value, valueFraction, refValue] = spans(
+        domainField.boundaryField().refGrad(),
+        domainField.boundaryField().value(),
+        domainField.boundaryField().valueFraction(),
+        domainField.boundaryField().refValue()
+    );
 
     NeoFOAM::parallelFor(
         domainField.exec(),
         range,
         KOKKOS_LAMBDA(const size_t i) {
             refValue[i] = fixedValue;
             value[i] = fixedValue;
+            valueFraction[i] = 1.0;      // only used refValue
+            refGradient[i] = fixedValue; // not used
         }
     );
 }