From 6e7452cbcda378b9e68e1ce7fc1e2b62079db831 Mon Sep 17 00:00:00 2001
From: Yu virtonomy <liao@virtonomy.io>
Date: Mon, 13 Jan 2025 16:31:08 +0100
Subject: [PATCH 1/7] add interaction_ck for ck version

---
 .../all_general_dynamics_ck.h                 |   1 +
 .../all_surface_indication_ck.h               |  37 ++++
 .../surface_indication_ck.h                   | 159 ++++++++++++++++
 .../surface_indication_ck.hpp                 | 169 ++++++++++++++++++
 .../test_2d_dambreak_ck/dambreak_ck.cpp       |   8 +-
 5 files changed, 373 insertions(+), 1 deletion(-)
 create mode 100644 src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/all_surface_indication_ck.h
 create mode 100644 src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.h
 create mode 100644 src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp

diff --git a/src/shared/shared_ck/particle_dynamics/general_dynamics/all_general_dynamics_ck.h b/src/shared/shared_ck/particle_dynamics/general_dynamics/all_general_dynamics_ck.h
index ab55f0a22..d703b1cdc 100644
--- a/src/shared/shared_ck/particle_dynamics/general_dynamics/all_general_dynamics_ck.h
+++ b/src/shared/shared_ck/particle_dynamics/general_dynamics/all_general_dynamics_ck.h
@@ -35,3 +35,4 @@
 #include "geometric_dynamics.hpp"
 #include "interpolation_dynamics.hpp"
 #include "kernel_correction_ck.hpp"
+#include "all_surface_indication_ck.h"
diff --git a/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/all_surface_indication_ck.h b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/all_surface_indication_ck.h
new file mode 100644
index 000000000..d41482a0b
--- /dev/null
+++ b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/all_surface_indication_ck.h
@@ -0,0 +1,37 @@
+/* ------------------------------------------------------------------------- *
+ *                                SPHinXsys                                  *
+ * ------------------------------------------------------------------------- *
+ * SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle *
+ * Hydrodynamics for industrial compleX systems. It provides C++ APIs for    *
+ * physical accurate simulation and aims to model coupled industrial dynamic *
+ * systems including fluid, solid, multi-body dynamics and beyond with SPH   *
+ * (smoothed particle hydrodynamics), a meshless computational method using  *
+ * particle discretization.                                                  *
+ *                                                                           *
+ * SPHinXsys is partially funded by German Research Foundation               *
+ * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1,            *
+ *  HU1527/12-1 and HU1527/12-4.                                             *
+ *                                                                           *
+ * Portions copyright (c) 2017-2023 Technical University of Munich and       *
+ * the authors' affiliations.                                                *
+ *                                                                           *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may   *
+ * not use this file except in compliance with the License. You may obtain a *
+ * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.        *
+ *                                                                           *
+ * ------------------------------------------------------------------------- */
+
+#ifndef B9F4AD2C_8E17_4C49_BCAA_4BD2C6A6951F
+#define B9F4AD2C_8E17_4C49_BCAA_4BD2C6A6951F
+/**
+ * @file    all_general_dynamics.h
+ * @brief   This is the header file that user code should include to pick up all
+ *          general dynamics used in SPHinXsys.
+ * @author	Chi Zhang and Xiangyu Hu
+ */
+
+
+#include "surface_indication_ck.hpp"
+
+
+#endif /* B9F4AD2C_8E17_4C49_BCAA_4BD2C6A6951F */
diff --git a/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.h b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.h
new file mode 100644
index 000000000..323391020
--- /dev/null
+++ b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.h
@@ -0,0 +1,159 @@
+#ifndef SURFACE_INDICATION_CK_H
+#define SURFACE_INDICATION_CK_H
+
+#include "base_general_dynamics.h"
+#include "base_fluid_dynamics.h"
+#include "interaction_ck.hpp"
+
+namespace SPH
+{
+namespace fluid_dynamics
+{
+
+template <typename...>
+class IndicatorCK;
+
+template <typename... RelationTypes>
+class FreeSurfaceIndicationCK;
+
+template <>
+class IndicatorCK<>
+{
+public:
+    IndicatorCK(BaseParticles *particles) {}
+
+    class ComputingKernel
+    {
+    public:
+        template <class ExecutionPolicy, class ComputingKernelType>
+        ComputingKernel(const ExecutionPolicy &ex_policy,
+                        IndicatorCK<> &encloser,
+                        ComputingKernelType &computing_kernel)
+            : threshold_by_dimensions_(0.75 * Dimensions) {}
+
+        Real operator()(Real &surface_indicator) { return surface_indicator; }
+
+    protected:
+        Real threshold_by_dimensions_;
+    };
+};
+
+template<>
+class IndicatorCK<Internal> : public IndicatorCK<> {
+public:
+    IndicatorCK(BaseParticles *particles) : IndicatorCK<>(particles) {}
+    using ComputingKernel = typename IndicatorCK<>::ComputingKernel;
+};
+
+template <template <typename...> class RelationType, typename... Parameters>
+class FreeSurfaceIndicationCK<Base, RelationType<Parameters...>>
+    : public Interaction<RelationType<Parameters...>>
+{
+public:
+    template <class BaseRelationType>
+    explicit FreeSurfaceIndicationCK(BaseRelationType &base_relation);
+    virtual ~FreeSurfaceIndicationCK() {}
+
+    class InteractKernel : public Interaction<RelationType<Parameters...>>::InteractKernel
+    {
+    public:
+        template <class ExecutionPolicy, typename... Args>
+        InteractKernel(const ExecutionPolicy &ex_policy,
+                       FreeSurfaceIndicationCK<Base, RelationType<Parameters...>> &encloser,
+                       Args &&...args);
+
+        virtual void interact(size_t index_i, Real dt = 0.0) {}
+
+    protected:
+        int *indicator_;
+        Real *pos_div_, *Vol_;
+        Real threshold_by_dimensions_;
+        Real smoothing_length_;
+    };
+
+protected:
+    DiscreteVariable<int> *dv_indicator_;
+    DiscreteVariable<Real> *dv_pos_div_, *dv_Vol_;
+    Real dv_threshold_by_dimensions_;
+    Real dv_smoothing_length_;
+};
+
+template <class FlowType, typename... Parameters>
+class FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>
+    : public FreeSurfaceIndicationCK<Base, Inner<Parameters...>>
+{
+    using IndicationKernel = typename IndicatorCK<FlowType>::ComputingKernel;
+
+public:
+    explicit FreeSurfaceIndicationCK(Relation<Inner<Parameters...>> &inner_relation);
+    virtual ~FreeSurfaceIndicationCK() {}
+
+    class InteractKernel : public FreeSurfaceIndicationCK<Base, Inner<Parameters...>>::InteractKernel
+    {
+    public:
+        template <class ExecutionPolicy>
+        InteractKernel(const ExecutionPolicy &ex_policy,
+                       FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>> &encloser);
+
+        void interact(size_t index_i, Real dt = 0.0) override;
+        int *previous_surface_indicator_;
+    };
+
+    class UpdateKernel : public FreeSurfaceIndicationCK<Base, Inner<Parameters...>>::InteractKernel
+    {
+    public:
+        template <class ExecutionPolicy>
+        UpdateKernel(const ExecutionPolicy &ex_policy,
+                     FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>> &encloser);
+
+        void update(size_t index_i, Real dt = 0.0);
+
+    protected:
+        IndicationKernel indication_;
+        int *previous_surface_indicator_;
+        FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>* outer_;
+    };
+
+protected:
+    IndicatorCK<FlowType> indicator_method_;
+    DiscreteVariable<int> *dv_previous_surface_indicator_;
+};
+
+using FreeSurfaceIndicationInnerCK = FreeSurfaceIndicationCK<Inner<WithUpdate, Internal>>;
+
+template <typename... Parameters>
+class FreeSurfaceIndicationCK<Contact<Parameters...>>
+    : public FreeSurfaceIndicationCK<Base, Contact<Parameters...>>
+{
+  public:
+    explicit FreeSurfaceIndicationCK(Relation<Contact<Parameters...>> &contact_relation);
+    virtual ~FreeSurfaceIndicationCK(){};
+
+    class InteractKernel
+        : public FreeSurfaceIndicationCK<Base, Contact<Parameters...>>::InteractKernel
+    {
+      public:
+        template <class ExecutionPolicy>
+        InteractKernel(const ExecutionPolicy &ex_policy,
+                       FreeSurfaceIndicationCK<Contact<Parameters...>> &encloser,
+                       size_t contact_index);
+        void interact(size_t index_i, Real dt = 0.0);
+
+      protected:
+
+        Real *contact_Vol_;
+    };
+
+  protected:
+
+        StdVec<DiscreteVariable<Real> *> dv_contact_Vol_;
+
+};
+
+using FreeSurfaceIndicationComplexCK = FreeSurfaceIndicationCK<Inner<WithUpdate, Internal>, Contact<>>;
+
+
+} // namespace fluid_dynamics
+} // namespace SPH
+
+#endif // SURFACE_INDICATION_CK_H
diff --git a/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp
new file mode 100644
index 000000000..1e31f71e7
--- /dev/null
+++ b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp
@@ -0,0 +1,169 @@
+#pragma once
+
+#include "surface_indication_ck.h"
+#include "base_particles.hpp"
+
+namespace SPH
+{
+namespace fluid_dynamics
+{
+
+template <template <typename...> class RelationType, typename... Parameters>
+template <class BaseRelationType>
+FreeSurfaceIndicationCK<Base, RelationType<Parameters...>>::FreeSurfaceIndicationCK(BaseRelationType &base_relation)
+    : Interaction<RelationType<Parameters...>>(base_relation),
+      dv_indicator_(this->particles_->template registerStateVariableOnly<int>("SurfaceIndicator")),
+      dv_pos_div_(this->particles_->template registerStateVariableOnly<Real>("PositionDivergence")),
+      dv_Vol_(this->particles_->template getVariableByName<Real>("VolumetricMeasure")),
+      dv_threshold_by_dimensions_(0.75 * Dimensions),
+      dv_smoothing_length_(this->sph_body_.sph_adaptation_->ReferenceSmoothingLength())
+{
+    std::cout << "surface_indication_ck" << std::endl;
+}
+
+template <template <typename...> class RelationType, typename... Parameters>
+template <class ExecutionPolicy, typename... Args>
+FreeSurfaceIndicationCK<Base, RelationType<Parameters...>>::InteractKernel::InteractKernel(
+    const ExecutionPolicy &ex_policy,
+    FreeSurfaceIndicationCK<Base, RelationType<Parameters...>> &encloser,
+    Args &&...args)
+    : Interaction<RelationType<Parameters...>>::InteractKernel(ex_policy, encloser, std::forward<Args>(args)...),
+      indicator_(encloser.dv_indicator_->DelegatedData(ex_policy)),
+      pos_div_(encloser.dv_pos_div_->DelegatedData(ex_policy)),
+      Vol_(encloser.dv_Vol_->DelegatedData(ex_policy)),
+      threshold_by_dimensions_(encloser.dv_threshold_by_dimensions_),
+      smoothing_length_(encloser.dv_smoothing_length_)
+{
+}
+
+template <class FlowType, typename... Parameters>
+FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::
+    FreeSurfaceIndicationCK(Relation<Inner<Parameters...>> &inner_relation)
+    : FreeSurfaceIndicationCK<Base, Inner<Parameters...>>(inner_relation),
+      indicator_method_(this->particles_),
+      dv_previous_surface_indicator_(this->particles_->template registerStateVariableOnly<int>("PreviousSurfaceIndicator"))
+{
+}
+
+template <class FlowType, typename... Parameters>
+template <class ExecutionPolicy>
+FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::InteractKernel::
+    InteractKernel(const ExecutionPolicy &ex_policy,
+                   FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>> &encloser)
+    : FreeSurfaceIndicationCK<Base, Inner<Parameters...>>::InteractKernel(ex_policy, encloser),
+      previous_surface_indicator_(encloser.dv_previous_surface_indicator_->DelegatedData(ex_policy))
+{
+}
+
+template <class FlowType, typename... Parameters>
+void FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::
+    InteractKernel::interact(size_t index_i, Real dt)
+{
+    Real pos_div = 0.0;
+    for (UnsignedInt n = this->FirstNeighbor(index_i); n != this->LastNeighbor(index_i); ++n)
+    {
+        UnsignedInt index_j = this->neighbor_index_[n];
+        Real r_ij = this->vec_r_ij(index_i, index_j).norm();
+        pos_div -= this->dW_ij(index_i, index_j) * this->Vol_[index_j] * r_ij;
+    }
+    this->pos_div_[index_i] = pos_div;
+}
+
+template <class FlowType, typename... Parameters>
+template <class ExecutionPolicy>
+FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::UpdateKernel::
+UpdateKernel(const ExecutionPolicy &ex_policy,
+             FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>> &encloser)
+    : FreeSurfaceIndicationCK<Base, Inner<Parameters...>>::InteractKernel(ex_policy, encloser),
+      indication_(ex_policy, encloser.indicator_method_, *this),
+      previous_surface_indicator_(encloser.dv_previous_surface_indicator_->DelegatedData(ex_policy)),
+      outer_(&encloser)
+{
+}
+
+template <class FlowType, typename... Parameters>
+void FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::
+UpdateKernel::update(size_t index_i, Real dt)
+{
+
+    int is_near_surface = 0;
+    for (UnsignedInt n = this->FirstNeighbor(index_i); n != this->LastNeighbor(index_i); ++n)
+    {
+        UnsignedInt index_j = this->neighbor_index_[n];
+        Real r_ij = this->vec_r_ij(index_i, index_j).norm();
+        if ( this->pos_div_[index_j]  < this->threshold_by_dimensions_ &&
+             r_ij < this->smoothing_length_)
+        {
+            is_near_surface = 1;
+            break;
+        }
+    }
+
+    int is_near_previous_surface = 0;
+    for (UnsignedInt n = this->FirstNeighbor(index_i); n != this->LastNeighbor(index_i); ++n)
+    {
+        if (this->previous_surface_indicator_[this->neighbor_index_[n]] == 1)
+        {
+            is_near_previous_surface = 1;
+            break;
+        }
+    }
+
+    if(this->pos_div_[index_i] < this->threshold_by_dimensions_ &&
+       is_near_surface != 1 &&
+       is_near_previous_surface != 1)
+    {
+        this->pos_div_[index_i] = 2.0 * this->threshold_by_dimensions_;
+    }
+    int indicator = 1;
+    if (this->pos_div_[index_i] > this->threshold_by_dimensions_ && is_near_surface !=1)
+        indicator = 0;
+
+
+    this->indicator_[index_i] = indicator;
+    this->previous_surface_indicator_[index_i] = indicator;
+}
+
+// Constructor for FreeSurfaceIndicationCK<Contact<Parameters...>>
+template <typename... Parameters>
+FreeSurfaceIndicationCK<Contact<Parameters...>>::
+FreeSurfaceIndicationCK(Relation<Contact<Parameters...>> &contact_relation)
+    : FreeSurfaceIndicationCK<Base, Contact<Parameters...>>(contact_relation)
+{
+    for (size_t k = 0; k != this->contact_particles_.size(); ++k)
+    {
+        dv_contact_Vol_.push_back(this->contact_particles_[k]->template getVariableByName<Real>("VolumetricMeasure"));
+    }
+}
+
+// Constructor for InteractKernel of FreeSurfaceIndicationCK<Contact<Parameters...>>
+template <typename... Parameters>
+template <class ExecutionPolicy>
+FreeSurfaceIndicationCK<Contact<Parameters...>>::InteractKernel::
+InteractKernel(const ExecutionPolicy &ex_policy,
+               FreeSurfaceIndicationCK<Contact<Parameters...>> &encloser,
+               size_t contact_index)
+    : FreeSurfaceIndicationCK<Base, Contact<Parameters...>>::InteractKernel(ex_policy, encloser, contact_index),
+      contact_Vol_(encloser.dv_contact_Vol_[contact_index]->DelegatedData(ex_policy))      
+{
+}
+
+// Implementation of InteractKernel::interact for FreeSurfaceIndicationCK<Contact<Parameters...>>
+template <typename... Parameters>
+void FreeSurfaceIndicationCK<Contact<Parameters...>>::InteractKernel::
+interact(size_t index_i, Real dt)
+{
+
+        Real pos_div = 0.0;
+
+    for (UnsignedInt n = this->FirstNeighbor(index_i); n != this->LastNeighbor(index_i); ++n)
+    {
+        UnsignedInt index_j = this->neighbor_index_[n];
+        Real r_ij = this->vec_r_ij(index_i, index_j).norm();
+
+        this->pos_div_[index_i] -= this->dW_ij(index_i, index_j) * this->contact_Vol_[index_j] * r_ij;
+    }
+    this->pos_div_[index_i] += pos_div;
+}
+} // namespace fluid_dynamics
+} // namespace SPH
diff --git a/tests/2d_examples/2d_examples_ck/test_2d_dambreak_ck/dambreak_ck.cpp b/tests/2d_examples/2d_examples_ck/test_2d_dambreak_ck/dambreak_ck.cpp
index e47e58b06..8df813330 100644
--- a/tests/2d_examples/2d_examples_ck/test_2d_dambreak_ck/dambreak_ck.cpp
+++ b/tests/2d_examples/2d_examples_ck/test_2d_dambreak_ck/dambreak_ck.cpp
@@ -113,7 +113,9 @@ int main(int ac, char *av[])
         fluid_acoustic_step_2nd_half(water_block_inner, water_wall_contact);
     InteractionDynamicsCK<MainExecutionPolicy, fluid_dynamics::DensityRegularizationComplexFreeSurface>
         fluid_density_regularization(water_block_inner, water_wall_contact);
-
+    InteractionDynamicsCK<MainExecutionPolicy, fluid_dynamics::FreeSurfaceIndicationComplexCK>
+        fluid_boundary_indicator(water_block_inner,water_wall_contact);
+        
     ReduceDynamicsCK<MainExecutionPolicy, fluid_dynamics::AdvectionTimeStepCK> fluid_advection_time_step(water_block, U_ref);
     ReduceDynamicsCK<MainExecutionPolicy, fluid_dynamics::AcousticTimeStepCK> fluid_acoustic_time_step(water_block);
     //----------------------------------------------------------------------
@@ -123,6 +125,9 @@ int main(int ac, char *av[])
     BodyStatesRecordingToVtp body_states_recording(sph_system);
     body_states_recording.addToWrite<Vecd>(wall_boundary, "NormalDirection");
     body_states_recording.addToWrite<Real>(water_block, "Density");
+    body_states_recording.addToWrite<int>(water_block, "SurfaceIndicator");
+    body_states_recording.addToWrite<Real>(water_block, "PositionDivergence");  
+
     RestartIO restart_io(sph_system);
 
     RegressionTestDynamicTimeWarping<ReducedQuantityRecording<MainExecutionPolicy, TotalMechanicalEnergyCK>>
@@ -187,6 +192,7 @@ int main(int ac, char *av[])
 
             fluid_density_regularization.exec();
             water_advection_step_setup.exec();
+            fluid_boundary_indicator.exec();
             Real advection_dt = fluid_advection_time_step.exec();
             fluid_linear_correction_matrix.exec();
             interval_computing_time_step += TickCount::now() - time_instance;

From 18dd206950b552985bc36d7def810fdd7fa53ac9 Mon Sep 17 00:00:00 2001
From: Yu virtonomy <liao@virtonomy.io>
Date: Tue, 14 Jan 2025 11:54:28 +0100
Subject: [PATCH 2/7] Refactor code structure and add comments for clarity

---
 .../surface_indication_ck.h                   |  12 +++
 .../surface_indication_ck.hpp                 | 100 ++++++++----------
 2 files changed, 59 insertions(+), 53 deletions(-)

diff --git a/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.h b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.h
index 323391020..5218c0770 100644
--- a/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.h
+++ b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.h
@@ -10,12 +10,24 @@ namespace SPH
 namespace fluid_dynamics
 {
 
+
 template <typename...>
 class IndicatorCK;
+/**
+ * @class FreeSurfaceIndicationCK
+ * @brief A free-surface indication class for specialized relation types.
+ */
 
 template <typename... RelationTypes>
 class FreeSurfaceIndicationCK;
 
+/**
+ * @class IndicatorCK
+ * @brief This class template can be specialized for different relation types. 
+ * Currently, this is not useful as this class was firstly construct for scenario with free surface + spatial.
+ * Can be extend further once other type of surface indication is required. 
+ * 
+ */
 template <>
 class IndicatorCK<>
 {
diff --git a/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp
index 1e31f71e7..65884a68f 100644
--- a/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp
+++ b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp
@@ -7,7 +7,7 @@ namespace SPH
 {
 namespace fluid_dynamics
 {
-
+//=================================================================================================//
 template <template <typename...> class RelationType, typename... Parameters>
 template <class BaseRelationType>
 FreeSurfaceIndicationCK<Base, RelationType<Parameters...>>::FreeSurfaceIndicationCK(BaseRelationType &base_relation)
@@ -16,11 +16,8 @@ FreeSurfaceIndicationCK<Base, RelationType<Parameters...>>::FreeSurfaceIndicatio
       dv_pos_div_(this->particles_->template registerStateVariableOnly<Real>("PositionDivergence")),
       dv_Vol_(this->particles_->template getVariableByName<Real>("VolumetricMeasure")),
       dv_threshold_by_dimensions_(0.75 * Dimensions),
-      dv_smoothing_length_(this->sph_body_.sph_adaptation_->ReferenceSmoothingLength())
-{
-    std::cout << "surface_indication_ck" << std::endl;
-}
-
+      dv_smoothing_length_(this->sph_body_.sph_adaptation_->ReferenceSmoothingLength()){}
+//=================================================================================================//
 template <template <typename...> class RelationType, typename... Parameters>
 template <class ExecutionPolicy, typename... Args>
 FreeSurfaceIndicationCK<Base, RelationType<Parameters...>>::InteractKernel::InteractKernel(
@@ -32,29 +29,23 @@ FreeSurfaceIndicationCK<Base, RelationType<Parameters...>>::InteractKernel::Inte
       pos_div_(encloser.dv_pos_div_->DelegatedData(ex_policy)),
       Vol_(encloser.dv_Vol_->DelegatedData(ex_policy)),
       threshold_by_dimensions_(encloser.dv_threshold_by_dimensions_),
-      smoothing_length_(encloser.dv_smoothing_length_)
-{
-}
-
+      smoothing_length_(encloser.dv_smoothing_length_){}
+//=================================================================================================//
 template <class FlowType, typename... Parameters>
 FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::
     FreeSurfaceIndicationCK(Relation<Inner<Parameters...>> &inner_relation)
     : FreeSurfaceIndicationCK<Base, Inner<Parameters...>>(inner_relation),
       indicator_method_(this->particles_),
-      dv_previous_surface_indicator_(this->particles_->template registerStateVariableOnly<int>("PreviousSurfaceIndicator"))
-{
-}
-
+      dv_previous_surface_indicator_(this->particles_->template registerStateVariableOnly<int>("PreviousSurfaceIndicator")){}
+//=================================================================================================//
 template <class FlowType, typename... Parameters>
 template <class ExecutionPolicy>
 FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::InteractKernel::
     InteractKernel(const ExecutionPolicy &ex_policy,
                    FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>> &encloser)
     : FreeSurfaceIndicationCK<Base, Inner<Parameters...>>::InteractKernel(ex_policy, encloser),
-      previous_surface_indicator_(encloser.dv_previous_surface_indicator_->DelegatedData(ex_policy))
-{
-}
-
+      previous_surface_indicator_(encloser.dv_previous_surface_indicator_->DelegatedData(ex_policy)){}
+//=================================================================================================//
 template <class FlowType, typename... Parameters>
 void FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::
     InteractKernel::interact(size_t index_i, Real dt)
@@ -68,7 +59,7 @@ void FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::
     }
     this->pos_div_[index_i] = pos_div;
 }
-
+//=================================================================================================//
 template <class FlowType, typename... Parameters>
 template <class ExecutionPolicy>
 FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::UpdateKernel::
@@ -77,54 +68,62 @@ UpdateKernel(const ExecutionPolicy &ex_policy,
     : FreeSurfaceIndicationCK<Base, Inner<Parameters...>>::InteractKernel(ex_policy, encloser),
       indication_(ex_policy, encloser.indicator_method_, *this),
       previous_surface_indicator_(encloser.dv_previous_surface_indicator_->DelegatedData(ex_policy)),
-      outer_(&encloser)
-{
-}
-
+      outer_(&encloser){}
+//=================================================================================================//
 template <class FlowType, typename... Parameters>
 void FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::
 UpdateKernel::update(size_t index_i, Real dt)
 {
-
-    int is_near_surface = 0;
+    // Check if the particle is near a surface
+    bool is_near_surface = false;
     for (UnsignedInt n = this->FirstNeighbor(index_i); n != this->LastNeighbor(index_i); ++n)
     {
-        UnsignedInt index_j = this->neighbor_index_[n];
-        Real r_ij = this->vec_r_ij(index_i, index_j).norm();
-        if ( this->pos_div_[index_j]  < this->threshold_by_dimensions_ &&
-             r_ij < this->smoothing_length_)
+        const UnsignedInt index_j = this->neighbor_index_[n];
+        const Real r_ij = this->vec_r_ij(index_i, index_j).norm();
+
+        // If a neighbor is under the threshold and within the smoothing length, mark near-surface
+        if ((this->pos_div_[index_j] < this->threshold_by_dimensions_) &&
+            (r_ij < this->smoothing_length_))
         {
-            is_near_surface = 1;
+            is_near_surface = true;
             break;
         }
     }
 
-    int is_near_previous_surface = 0;
+    // Check if the particle is near a previously marked surface
+    // Equvilent to SpatialTemporal
+    bool is_near_previous_surface = false;
     for (UnsignedInt n = this->FirstNeighbor(index_i); n != this->LastNeighbor(index_i); ++n)
     {
-        if (this->previous_surface_indicator_[this->neighbor_index_[n]] == 1)
+        const UnsignedInt index_j = this->neighbor_index_[n];
+        if (this->previous_surface_indicator_[index_j] == 1)
         {
-            is_near_previous_surface = 1;
+            is_near_previous_surface = true;
             break;
         }
     }
 
-    if(this->pos_div_[index_i] < this->threshold_by_dimensions_ &&
-       is_near_surface != 1 &&
-       is_near_previous_surface != 1)
+    // If the current particle is under the threshold but not near any surface (current or previous),
+    // artificially increase its pos_div to avoid misclassifying it as surface
+    if ((this->pos_div_[index_i] < this->threshold_by_dimensions_) &&
+        !is_near_surface && !is_near_previous_surface)
     {
         this->pos_div_[index_i] = 2.0 * this->threshold_by_dimensions_;
     }
-    int indicator = 1;
-    if (this->pos_div_[index_i] > this->threshold_by_dimensions_ && is_near_surface !=1)
-        indicator = 0;
 
+    // Determine the new indicator for this particle
+    int new_indicator = 1;  // Default: surface
+    // If the local pos_div is above the threshold and it is not near-surface, mark as not surface
+    if ((this->pos_div_[index_i] > this->threshold_by_dimensions_) && !is_near_surface)
+    {
+        new_indicator = 0;
+    }
 
-    this->indicator_[index_i] = indicator;
-    this->previous_surface_indicator_[index_i] = indicator;
+    // Update both current and previous indicators
+    this->indicator_[index_i] = new_indicator;
+    this->previous_surface_indicator_[index_i] = new_indicator;
 }
-
-// Constructor for FreeSurfaceIndicationCK<Contact<Parameters...>>
+//=================================================================================================//
 template <typename... Parameters>
 FreeSurfaceIndicationCK<Contact<Parameters...>>::
 FreeSurfaceIndicationCK(Relation<Contact<Parameters...>> &contact_relation)
@@ -135,8 +134,7 @@ FreeSurfaceIndicationCK(Relation<Contact<Parameters...>> &contact_relation)
         dv_contact_Vol_.push_back(this->contact_particles_[k]->template getVariableByName<Real>("VolumetricMeasure"));
     }
 }
-
-// Constructor for InteractKernel of FreeSurfaceIndicationCK<Contact<Parameters...>>
+//=================================================================================================//
 template <typename... Parameters>
 template <class ExecutionPolicy>
 FreeSurfaceIndicationCK<Contact<Parameters...>>::InteractKernel::
@@ -144,18 +142,13 @@ InteractKernel(const ExecutionPolicy &ex_policy,
                FreeSurfaceIndicationCK<Contact<Parameters...>> &encloser,
                size_t contact_index)
     : FreeSurfaceIndicationCK<Base, Contact<Parameters...>>::InteractKernel(ex_policy, encloser, contact_index),
-      contact_Vol_(encloser.dv_contact_Vol_[contact_index]->DelegatedData(ex_policy))      
-{
-}
-
-// Implementation of InteractKernel::interact for FreeSurfaceIndicationCK<Contact<Parameters...>>
+      contact_Vol_(encloser.dv_contact_Vol_[contact_index]->DelegatedData(ex_policy)){}
+//=================================================================================================//
 template <typename... Parameters>
 void FreeSurfaceIndicationCK<Contact<Parameters...>>::InteractKernel::
 interact(size_t index_i, Real dt)
 {
-
-        Real pos_div = 0.0;
-
+    Real pos_div = 0.0;
     for (UnsignedInt n = this->FirstNeighbor(index_i); n != this->LastNeighbor(index_i); ++n)
     {
         UnsignedInt index_j = this->neighbor_index_[n];
@@ -165,5 +158,6 @@ interact(size_t index_i, Real dt)
     }
     this->pos_div_[index_i] += pos_div;
 }
+//=================================================================================================//
 } // namespace fluid_dynamics
 } // namespace SPH

From 5fe08aa17e51a275170876c75939ac77cdb2541a Mon Sep 17 00:00:00 2001
From: Yu virtonomy <liao@virtonomy.io>
Date: Tue, 14 Jan 2025 14:20:28 +0100
Subject: [PATCH 3/7] Remove virtual keyword from interact method in
 surface_indication_ck.h for sycl

---
 .../surface_indication/surface_indication_ck.h                | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.h b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.h
index 5218c0770..9187f25c3 100644
--- a/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.h
+++ b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.h
@@ -74,7 +74,7 @@ class FreeSurfaceIndicationCK<Base, RelationType<Parameters...>>
                        FreeSurfaceIndicationCK<Base, RelationType<Parameters...>> &encloser,
                        Args &&...args);
 
-        virtual void interact(size_t index_i, Real dt = 0.0) {}
+        void interact(size_t index_i, Real dt = 0.0) {}
 
     protected:
         int *indicator_;
@@ -107,7 +107,7 @@ class FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>
         InteractKernel(const ExecutionPolicy &ex_policy,
                        FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>> &encloser);
 
-        void interact(size_t index_i, Real dt = 0.0) override;
+        void interact(size_t index_i, Real dt = 0.0);
         int *previous_surface_indicator_;
     };
 

From 68d15ddf220a9f7e9200e57cf57a16b05bcb15a7 Mon Sep 17 00:00:00 2001
From: Yu virtonomy <liao@virtonomy.io>
Date: Tue, 14 Jan 2025 16:22:00 +0100
Subject: [PATCH 4/7] Add fluid boundary indicator to dambreak_sycl test

---
 .../2d_examples/test_2d_dambreak_sycl/dambreak_sycl.cpp      | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/tests/tests_sycl/2d_examples/test_2d_dambreak_sycl/dambreak_sycl.cpp b/tests/tests_sycl/2d_examples/test_2d_dambreak_sycl/dambreak_sycl.cpp
index 2525defcb..fce6233f8 100644
--- a/tests/tests_sycl/2d_examples/test_2d_dambreak_sycl/dambreak_sycl.cpp
+++ b/tests/tests_sycl/2d_examples/test_2d_dambreak_sycl/dambreak_sycl.cpp
@@ -113,6 +113,8 @@ int main(int ac, char *av[])
         fluid_acoustic_step_2nd_half(water_block_inner, water_wall_contact);
     InteractionDynamicsCK<MainExecutionPolicy, fluid_dynamics::DensityRegularizationComplexFreeSurface>
         fluid_density_regularization(water_block_inner, water_wall_contact);
+    InteractionDynamicsCK<MainExecutionPolicy, fluid_dynamics::FreeSurfaceIndicationComplexCK>
+        fluid_boundary_indicator(water_block_inner,water_wall_contact);
 
     ReduceDynamicsCK<MainExecutionPolicy, fluid_dynamics::AdvectionTimeStepCK> fluid_advection_time_step(water_block, U_ref);
     ReduceDynamicsCK<MainExecutionPolicy, fluid_dynamics::AcousticTimeStepCK> fluid_acoustic_time_step(water_block);
@@ -123,6 +125,8 @@ int main(int ac, char *av[])
     BodyStatesRecordingToVtp body_states_recording(sph_system);
     body_states_recording.addToWrite<Vecd>(wall_boundary, "NormalDirection");
     body_states_recording.addToWrite<Real>(water_block, "Density");
+    body_states_recording.addToWrite<int>(water_block, "SurfaceIndicator");
+
     RestartIO restart_io(sph_system);
 
     RegressionTestDynamicTimeWarping<ReducedQuantityRecording<MainExecutionPolicy, TotalMechanicalEnergyCK>>
@@ -185,6 +189,7 @@ int main(int ac, char *av[])
             /** outer loop for dual-time criteria time-stepping. */
             time_instance = TickCount::now();
 
+            fluid_boundary_indicator.exec();
             fluid_density_regularization.exec();
             water_advection_step_setup.exec();
             Real advection_dt = fluid_advection_time_step.exec();

From 08cb82368f13e074ef7089292225f1699184241c Mon Sep 17 00:00:00 2001
From: Yu virtonomy <liao@virtonomy.io>
Date: Mon, 20 Jan 2025 10:46:18 +0100
Subject: [PATCH 5/7] Remove IndicatorCK; Refactor FreeSurfaceIndicationCK
 class for improved clarity and structure

In current structure, IndicatorCK for surface indicator is not needed
---
 .../surface_indication_ck.h                   | 132 ++++++++++--------
 .../surface_indication_ck.hpp                 | 101 ++++++++------
 2 files changed, 132 insertions(+), 101 deletions(-)

diff --git a/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.h b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.h
index 9187f25c3..d30d49a26 100644
--- a/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.h
+++ b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.h
@@ -9,54 +9,25 @@ namespace SPH
 {
 namespace fluid_dynamics
 {
-
-
-template <typename...>
-class IndicatorCK;
 /**
  * @class FreeSurfaceIndicationCK
- * @brief A free-surface indication class for specialized relation types.
+ * @brief Free-surface indication with configurable relationship types.
+ * 
+ * This template is specialized for different combinations of
+ * "Base" + "Inner" / "Contact" relations to handle free-surface
+ * detection and updating in SPH simulations.
  */
-
 template <typename... RelationTypes>
 class FreeSurfaceIndicationCK;
 
+//=================================================================================================//
+// Base relation version
+//=================================================================================================//
 /**
- * @class IndicatorCK
- * @brief This class template can be specialized for different relation types. 
- * Currently, this is not useful as this class was firstly construct for scenario with free surface + spatial.
- * Can be extend further once other type of surface indication is required. 
- * 
+ * @class FreeSurfaceIndicationCK<Base, RelationType<Parameters...>>
+ * @brief Basic free-surface indication logic, storing and computing
+ *        positional divergence, surface indicators, etc.
  */
-template <>
-class IndicatorCK<>
-{
-public:
-    IndicatorCK(BaseParticles *particles) {}
-
-    class ComputingKernel
-    {
-    public:
-        template <class ExecutionPolicy, class ComputingKernelType>
-        ComputingKernel(const ExecutionPolicy &ex_policy,
-                        IndicatorCK<> &encloser,
-                        ComputingKernelType &computing_kernel)
-            : threshold_by_dimensions_(0.75 * Dimensions) {}
-
-        Real operator()(Real &surface_indicator) { return surface_indicator; }
-
-    protected:
-        Real threshold_by_dimensions_;
-    };
-};
-
-template<>
-class IndicatorCK<Internal> : public IndicatorCK<> {
-public:
-    IndicatorCK(BaseParticles *particles) : IndicatorCK<>(particles) {}
-    using ComputingKernel = typename IndicatorCK<>::ComputingKernel;
-};
-
 template <template <typename...> class RelationType, typename... Parameters>
 class FreeSurfaceIndicationCK<Base, RelationType<Parameters...>>
     : public Interaction<RelationType<Parameters...>>
@@ -66,6 +37,12 @@ class FreeSurfaceIndicationCK<Base, RelationType<Parameters...>>
     explicit FreeSurfaceIndicationCK(BaseRelationType &base_relation);
     virtual ~FreeSurfaceIndicationCK() {}
 
+    //------------------------------------------------------------------------------------------//
+    /**
+     * @class InteractKernel
+     * @brief Base interaction kernel that calculates positional divergence,
+     *        used in detecting free-surface particles.
+     */
     class InteractKernel : public Interaction<RelationType<Parameters...>>::InteractKernel
     {
     public:
@@ -74,33 +51,47 @@ class FreeSurfaceIndicationCK<Base, RelationType<Parameters...>>
                        FreeSurfaceIndicationCK<Base, RelationType<Parameters...>> &encloser,
                        Args &&...args);
 
-        void interact(size_t index_i, Real dt = 0.0) {}
+        void interact(size_t index_i, Real dt = 0.0);
 
     protected:
         int *indicator_;
-        Real *pos_div_, *Vol_;
+        Real *pos_div_;
+        Real *Vol_;
         Real threshold_by_dimensions_;
         Real smoothing_length_;
     };
 
 protected:
     DiscreteVariable<int> *dv_indicator_;
-    DiscreteVariable<Real> *dv_pos_div_, *dv_Vol_;
+    DiscreteVariable<Real> *dv_pos_div_;
+    DiscreteVariable<Real> *dv_Vol_;
     Real dv_threshold_by_dimensions_;
     Real dv_smoothing_length_;
 };
 
+//=================================================================================================//
+// Inner relation version with "WithUpdate"
+//=================================================================================================//
+/**
+ * @class FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>
+ * @brief Extends the base free-surface indication for inner relations that also require
+ *        an updating step (e.g., "WithUpdate").
+ */
 template <class FlowType, typename... Parameters>
 class FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>
     : public FreeSurfaceIndicationCK<Base, Inner<Parameters...>>
 {
-    using IndicationKernel = typename IndicatorCK<FlowType>::ComputingKernel;
-
 public:
     explicit FreeSurfaceIndicationCK(Relation<Inner<Parameters...>> &inner_relation);
     virtual ~FreeSurfaceIndicationCK() {}
 
-    class InteractKernel : public FreeSurfaceIndicationCK<Base, Inner<Parameters...>>::InteractKernel
+    //------------------------------------------------------------------------------------------//
+    /**
+     * @class InteractKernel
+     * @brief Interaction kernel for detecting free surface in an inner relation.
+     */
+    class InteractKernel
+        : public FreeSurfaceIndicationCK<Base, Inner<Parameters...>>::InteractKernel
     {
     public:
         template <class ExecutionPolicy>
@@ -108,10 +99,19 @@ class FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>
                        FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>> &encloser);
 
         void interact(size_t index_i, Real dt = 0.0);
+
+        /// Pointer to the previously stored surface indicator.
         int *previous_surface_indicator_;
     };
 
-    class UpdateKernel : public FreeSurfaceIndicationCK<Base, Inner<Parameters...>>::InteractKernel
+    //------------------------------------------------------------------------------------------//
+    /**
+     * @class UpdateKernel
+     * @brief Post-processing/update kernel that modifies the surface indicator
+     *        based on the newly computed positional divergences and neighbor information.
+     */
+    class UpdateKernel
+        : public FreeSurfaceIndicationCK<Base, Inner<Parameters...>>::InteractKernel
     {
     public:
         template <class ExecutionPolicy>
@@ -121,49 +121,61 @@ class FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>
         void update(size_t index_i, Real dt = 0.0);
 
     protected:
-        IndicationKernel indication_;
         int *previous_surface_indicator_;
-        FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>* outer_;
+        FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>> *outer_;
     };
 
 protected:
-    IndicatorCK<FlowType> indicator_method_;
     DiscreteVariable<int> *dv_previous_surface_indicator_;
 };
 
+//------------------------------------------------------------------------------------------//
+// Common type alias for internal flow with free-surface indication (with update).
 using FreeSurfaceIndicationInnerCK = FreeSurfaceIndicationCK<Inner<WithUpdate, Internal>>;
 
+//=================================================================================================//
+// Contact relation version
+//=================================================================================================//
+/**
+ * @class FreeSurfaceIndicationCK<Contact<Parameters...>>
+ * @brief Extends the base free-surface indication for contact relations.
+ */
 template <typename... Parameters>
 class FreeSurfaceIndicationCK<Contact<Parameters...>>
     : public FreeSurfaceIndicationCK<Base, Contact<Parameters...>>
 {
-  public:
+public:
     explicit FreeSurfaceIndicationCK(Relation<Contact<Parameters...>> &contact_relation);
-    virtual ~FreeSurfaceIndicationCK(){};
+    virtual ~FreeSurfaceIndicationCK() {}
 
+    //------------------------------------------------------------------------------------------//
+    /**
+     * @class InteractKernel
+     * @brief Interaction kernel for a contact relation, combining data from multiple contact bodies.
+     */
     class InteractKernel
         : public FreeSurfaceIndicationCK<Base, Contact<Parameters...>>::InteractKernel
     {
-      public:
+    public:
         template <class ExecutionPolicy>
         InteractKernel(const ExecutionPolicy &ex_policy,
                        FreeSurfaceIndicationCK<Contact<Parameters...>> &encloser,
                        size_t contact_index);
-        void interact(size_t index_i, Real dt = 0.0);
 
-      protected:
+        void interact(size_t index_i, Real dt = 0.0);
 
+    protected:
         Real *contact_Vol_;
     };
 
-  protected:
-
-        StdVec<DiscreteVariable<Real> *> dv_contact_Vol_;
-
+protected:
+    StdVec<DiscreteVariable<Real> *> dv_contact_Vol_;
 };
 
+//------------------------------------------------------------------------------------------//
+// Common type alias for complex free-surface indication (inner + contact).
 using FreeSurfaceIndicationComplexCK = FreeSurfaceIndicationCK<Inner<WithUpdate, Internal>, Contact<>>;
-
+//------------------------------------------------------------------------------------------//
 
 } // namespace fluid_dynamics
 } // namespace SPH
diff --git a/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp
index 65884a68f..28404393c 100644
--- a/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp
+++ b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp
@@ -8,47 +8,62 @@ namespace SPH
 namespace fluid_dynamics
 {
 //=================================================================================================//
+// FreeSurfaceIndicationCK<Base, RelationType<Parameters...>>
+//=================================================================================================//
 template <template <typename...> class RelationType, typename... Parameters>
 template <class BaseRelationType>
-FreeSurfaceIndicationCK<Base, RelationType<Parameters...>>::FreeSurfaceIndicationCK(BaseRelationType &base_relation)
+FreeSurfaceIndicationCK<Base, RelationType<Parameters...>>::
+FreeSurfaceIndicationCK(BaseRelationType &base_relation)
     : Interaction<RelationType<Parameters...>>(base_relation),
       dv_indicator_(this->particles_->template registerStateVariableOnly<int>("SurfaceIndicator")),
       dv_pos_div_(this->particles_->template registerStateVariableOnly<Real>("PositionDivergence")),
       dv_Vol_(this->particles_->template getVariableByName<Real>("VolumetricMeasure")),
       dv_threshold_by_dimensions_(0.75 * Dimensions),
-      dv_smoothing_length_(this->sph_body_.sph_adaptation_->ReferenceSmoothingLength()){}
-//=================================================================================================//
+      dv_smoothing_length_(this->sph_body_.sph_adaptation_->ReferenceSmoothingLength())
+{
+}
+
 template <template <typename...> class RelationType, typename... Parameters>
 template <class ExecutionPolicy, typename... Args>
-FreeSurfaceIndicationCK<Base, RelationType<Parameters...>>::InteractKernel::InteractKernel(
-    const ExecutionPolicy &ex_policy,
-    FreeSurfaceIndicationCK<Base, RelationType<Parameters...>> &encloser,
-    Args &&...args)
+FreeSurfaceIndicationCK<Base, RelationType<Parameters...>>::InteractKernel::
+InteractKernel(const ExecutionPolicy &ex_policy,
+               FreeSurfaceIndicationCK<Base, RelationType<Parameters...>> &encloser,
+               Args &&...args)
     : Interaction<RelationType<Parameters...>>::InteractKernel(ex_policy, encloser, std::forward<Args>(args)...),
       indicator_(encloser.dv_indicator_->DelegatedData(ex_policy)),
       pos_div_(encloser.dv_pos_div_->DelegatedData(ex_policy)),
       Vol_(encloser.dv_Vol_->DelegatedData(ex_policy)),
       threshold_by_dimensions_(encloser.dv_threshold_by_dimensions_),
-      smoothing_length_(encloser.dv_smoothing_length_){}
+      smoothing_length_(encloser.dv_smoothing_length_)
+{
+}
+
+//=================================================================================================//
+// FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>
 //=================================================================================================//
 template <class FlowType, typename... Parameters>
 FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::
-    FreeSurfaceIndicationCK(Relation<Inner<Parameters...>> &inner_relation)
+FreeSurfaceIndicationCK(Relation<Inner<Parameters...>> &inner_relation)
     : FreeSurfaceIndicationCK<Base, Inner<Parameters...>>(inner_relation),
-      indicator_method_(this->particles_),
-      dv_previous_surface_indicator_(this->particles_->template registerStateVariableOnly<int>("PreviousSurfaceIndicator")){}
-//=================================================================================================//
+      dv_previous_surface_indicator_(
+          this->particles_->template registerStateVariableOnly<int>("PreviousSurfaceIndicator"))
+{
+}
+
+//------------------------------------------------------------------------------------------//
 template <class FlowType, typename... Parameters>
 template <class ExecutionPolicy>
 FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::InteractKernel::
-    InteractKernel(const ExecutionPolicy &ex_policy,
-                   FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>> &encloser)
+InteractKernel(const ExecutionPolicy &ex_policy,
+               FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>> &encloser)
     : FreeSurfaceIndicationCK<Base, Inner<Parameters...>>::InteractKernel(ex_policy, encloser),
-      previous_surface_indicator_(encloser.dv_previous_surface_indicator_->DelegatedData(ex_policy)){}
-//=================================================================================================//
+      previous_surface_indicator_(encloser.dv_previous_surface_indicator_->DelegatedData(ex_policy))
+{
+}
+
 template <class FlowType, typename... Parameters>
-void FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::
-    InteractKernel::interact(size_t index_i, Real dt)
+void FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::InteractKernel::
+interact(size_t index_i, Real dt)
 {
     Real pos_div = 0.0;
     for (UnsignedInt n = this->FirstNeighbor(index_i); n != this->LastNeighbor(index_i); ++n)
@@ -59,29 +74,30 @@ void FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::
     }
     this->pos_div_[index_i] = pos_div;
 }
-//=================================================================================================//
+
+//------------------------------------------------------------------------------------------//
 template <class FlowType, typename... Parameters>
 template <class ExecutionPolicy>
 FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::UpdateKernel::
 UpdateKernel(const ExecutionPolicy &ex_policy,
              FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>> &encloser)
     : FreeSurfaceIndicationCK<Base, Inner<Parameters...>>::InteractKernel(ex_policy, encloser),
-      indication_(ex_policy, encloser.indicator_method_, *this),
       previous_surface_indicator_(encloser.dv_previous_surface_indicator_->DelegatedData(ex_policy)),
-      outer_(&encloser){}
-//=================================================================================================//
+      outer_(&encloser)
+{
+}
+
 template <class FlowType, typename... Parameters>
-void FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::
-UpdateKernel::update(size_t index_i, Real dt)
+void FreeSurfaceIndicationCK<Inner<WithUpdate, FlowType, Parameters...>>::UpdateKernel::
+update(size_t index_i, Real dt)
 {
-    // Check if the particle is near a surface
+    // Detect if near surface based on neighbors
     bool is_near_surface = false;
     for (UnsignedInt n = this->FirstNeighbor(index_i); n != this->LastNeighbor(index_i); ++n)
     {
         const UnsignedInt index_j = this->neighbor_index_[n];
-        const Real r_ij = this->vec_r_ij(index_i, index_j).norm();
+        Real r_ij = this->vec_r_ij(index_i, index_j).norm();
 
-        // If a neighbor is under the threshold and within the smoothing length, mark near-surface
         if ((this->pos_div_[index_j] < this->threshold_by_dimensions_) &&
             (r_ij < this->smoothing_length_))
         {
@@ -90,8 +106,7 @@ UpdateKernel::update(size_t index_i, Real dt)
         }
     }
 
-    // Check if the particle is near a previously marked surface
-    // Equvilent to SpatialTemporal
+    // Check if near a previously marked surface
     bool is_near_previous_surface = false;
     for (UnsignedInt n = this->FirstNeighbor(index_i); n != this->LastNeighbor(index_i); ++n)
     {
@@ -103,17 +118,16 @@ UpdateKernel::update(size_t index_i, Real dt)
         }
     }
 
-    // If the current particle is under the threshold but not near any surface (current or previous),
-    // artificially increase its pos_div to avoid misclassifying it as surface
+    // If the particle’s pos_div is under threshold but isolated (no surface neighbors),
+    // push it above threshold to avoid false positives.
     if ((this->pos_div_[index_i] < this->threshold_by_dimensions_) &&
         !is_near_surface && !is_near_previous_surface)
     {
         this->pos_div_[index_i] = 2.0 * this->threshold_by_dimensions_;
     }
 
-    // Determine the new indicator for this particle
-    int new_indicator = 1;  // Default: surface
-    // If the local pos_div is above the threshold and it is not near-surface, mark as not surface
+    // Determine new indicator: surface (1) or not (0)
+    int new_indicator = 1;
     if ((this->pos_div_[index_i] > this->threshold_by_dimensions_) && !is_near_surface)
     {
         new_indicator = 0;
@@ -123,6 +137,9 @@ UpdateKernel::update(size_t index_i, Real dt)
     this->indicator_[index_i] = new_indicator;
     this->previous_surface_indicator_[index_i] = new_indicator;
 }
+
+//=================================================================================================//
+// FreeSurfaceIndicationCK<Contact<Parameters...>>
 //=================================================================================================//
 template <typename... Parameters>
 FreeSurfaceIndicationCK<Contact<Parameters...>>::
@@ -131,10 +148,11 @@ FreeSurfaceIndicationCK(Relation<Contact<Parameters...>> &contact_relation)
 {
     for (size_t k = 0; k != this->contact_particles_.size(); ++k)
     {
-        dv_contact_Vol_.push_back(this->contact_particles_[k]->template getVariableByName<Real>("VolumetricMeasure"));
+        dv_contact_Vol_.push_back(
+            this->contact_particles_[k]->template getVariableByName<Real>("VolumetricMeasure"));
     }
 }
-//=================================================================================================//
+
 template <typename... Parameters>
 template <class ExecutionPolicy>
 FreeSurfaceIndicationCK<Contact<Parameters...>>::InteractKernel::
@@ -142,8 +160,10 @@ InteractKernel(const ExecutionPolicy &ex_policy,
                FreeSurfaceIndicationCK<Contact<Parameters...>> &encloser,
                size_t contact_index)
     : FreeSurfaceIndicationCK<Base, Contact<Parameters...>>::InteractKernel(ex_policy, encloser, contact_index),
-      contact_Vol_(encloser.dv_contact_Vol_[contact_index]->DelegatedData(ex_policy)){}
-//=================================================================================================//
+      contact_Vol_(encloser.dv_contact_Vol_[contact_index]->DelegatedData(ex_policy))
+{
+}
+
 template <typename... Parameters>
 void FreeSurfaceIndicationCK<Contact<Parameters...>>::InteractKernel::
 interact(size_t index_i, Real dt)
@@ -153,11 +173,10 @@ interact(size_t index_i, Real dt)
     {
         UnsignedInt index_j = this->neighbor_index_[n];
         Real r_ij = this->vec_r_ij(index_i, index_j).norm();
-
-        this->pos_div_[index_i] -= this->dW_ij(index_i, index_j) * this->contact_Vol_[index_j] * r_ij;
+        pos_div -= this->dW_ij(index_i, index_j) * this->contact_Vol_[index_j] * r_ij;
     }
     this->pos_div_[index_i] += pos_div;
 }
-//=================================================================================================//
+
 } // namespace fluid_dynamics
 } // namespace SPH

From 7a33fe69d13e19d75d91c2dfd448c7b63a28a2c9 Mon Sep 17 00:00:00 2001
From: Yu virtonomy <liao@virtonomy.io>
Date: Mon, 20 Jan 2025 11:14:29 +0100
Subject: [PATCH 6/7] Add include guards to surface_indication_ck.hpp for
 improved header safety

---
 .../surface_indication/surface_indication_ck.hpp               | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp
index 28404393c..125c1063e 100644
--- a/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp
+++ b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp
@@ -1,4 +1,5 @@
-#pragma once
+#ifndef SURFACE_INDICATION_CK_HPP
+#define SURFACE_INDICATION_CK_HPP
 
 #include "surface_indication_ck.h"
 #include "base_particles.hpp"

From 4360200947be81c927fc95202ed9d4e620076793 Mon Sep 17 00:00:00 2001
From: Yu virtonomy <liao@virtonomy.io>
Date: Mon, 20 Jan 2025 11:16:33 +0100
Subject: [PATCH 7/7] Add include guard to surface_indication_ck.hpp for header
 safety

---
 .../surface_indication/surface_indication_ck.hpp                 | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp
index 125c1063e..61a8ef932 100644
--- a/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp
+++ b/src/shared/shared_ck/particle_dynamics/general_dynamics/surface_indication/surface_indication_ck.hpp
@@ -181,3 +181,4 @@ interact(size_t index_i, Real dt)
 
 } // namespace fluid_dynamics
 } // namespace SPH
+#endif // SURFACE_INDICATION_CK_HPP