compute mixed derivatives

* refer #59

Signed-off-by: Lakshman Anumolu <[email protected]>

Author: acrlakshman

CMakeLists.txt

@@ -31,7 +31,7 @@ ENDIF()
 
 MARK_AS_ADVANCED(BUILD_COVERAGE)
 
-SET(CMAKE_CXX_STANDARD 14)
+SET(CMAKE_CXX_STANDARD 17)
 SET(CMAKE_CXX_STANDARD_REQUIRED ON)
 
 IF (BUILD_TESTS)

include/gals/cpu/levelset.h

@@ -84,6 +84,19 @@ class Levelset
    */
   ~Levelset();
 
+  /*! Compute mixed derivatives.
+   *
+   * Mixed derivatives of the levelset is computed using second order gradient approximation on
+   * the levelset gradient.
+   *
+   * \param psi gradient of levelset.
+   * \param phi_mixed_derivatives mixed derivatives of the levelset.
+   *
+   * \return Void.
+   */
+  void computeMixedDerivatives(const Array<T_GRID, Vec3<T>>& psi,
+                               Array<T_GRID, VecN<T, T_GRID::num_mixed_derivatives>>& phi_mixed_derivatives);
+
   /*! Return grid.
    *
    * \return grid.

src/cpu/levelset.cc

@@ -30,6 +30,7 @@
 ///////////////////////////////////////////////////////////////////////////////
 
 #include "gals/cpu/levelset.h"
+#include "gals/cpu/gradient.h"
 
 #include <iostream>
 
@@ -52,6 +53,32 @@ GALS::CPU::Levelset<T_GRID, T>::~Levelset()
 {
 }
 
+template <typename T_GRID, typename T>
+void GALS::CPU::Levelset<T_GRID, T>::computeMixedDerivatives(
+    const Array<T_GRID, Vec3<T>>& psi, Array<T_GRID, VecN<T, T_GRID::num_mixed_derivatives>>& phi_mixed_derivatives)
+{
+  // For now gradient computation is hardcoded here.
+  const Vec3<int> num_cells = psi.numCells();
+  const T_GRID& grid = psi.grid();
+  const Vec3<typename T_GRID::value_type> dx = grid.dX();
+  const auto& axis_vectors = GALS::CPU::Grid<typename T_GRID::value_type, T_GRID::dim>::axis_vectors;
+
+  if constexpr (T_GRID::dim == 1) {
+    for (int i = 0; i < num_cells[0]; ++i)
+      for (int j = 0; j < num_cells[1]; ++j)
+        for (int k = 0; k < num_cells[2]; ++k) {
+          for (int axis = 0; axis < T_GRID::dim; ++axis) {
+            typename T_GRID::value_type one_by_dx = static_cast<typename T_GRID::value_type>(1.) / dx[axis];
+
+            // grad_alpha(i, j, k)[axis] =
+            //(alpha(i + axis_vectors(axis, 0), j + axis_vectors(axis, 1), k + axis_vectors(axis, 2)) -
+            // alpha(i - axis_vectors(axis, 0), j - axis_vectors(axis, 1), k - axis_vectors(axis, 2))) *
+            // one_by_dx * static_cast<T>(0.5);
+          }
+        }
+  }
+}
+
 template <typename T_GRID, typename T>
 void GALS::CPU::Levelset<T_GRID, T>::print()
 {