Xiangyu/staggered surface tension

src/shared/particle_dynamics/base_local_dynamics.h

@@ -36,33 +36,77 @@
 
 namespace SPH
 {
-/** A Functor for Summation */
+//----------------------------------------------------------------------
+// Particle group scope functors
+//----------------------------------------------------------------------
+class AllParticles
+{
+  public:
+    AllParticles(BaseParticles *base_particles){};
+    bool operator()(size_t index_i)
+    {
+        return true;
+    };
+};
+
+template <int INDICATOR>
+class IndicatedParticles
+{
+    StdLargeVec<int> &indicator_;
+
+  public:
+    IndicatedParticles(BaseParticles *base_particles)
+        : indicator_(*base_particles->getVariableByName<int>("Indicator")){};
+    bool operator()(size_t index_i)
+    {
+        return indicator_[index_i] == INDICATOR;
+    };
+};
+using BulkParticles = IndicatedParticles<0>;
+
+template <int INDICATOR>
+class NotIndicatedParticles
+{
+    StdLargeVec<int> &indicator_;
+
+  public:
+    NotIndicatedParticles(BaseParticles *base_particles)
+        : indicator_(*base_particles->getVariableByName<int>("Indicator")){};
+    bool operator()(size_t index_i)
+    {
+        return indicator_[index_i] != INDICATOR;
+    };
+};
+
+//----------------------------------------------------------------------
+// Particle reduce functors
+//----------------------------------------------------------------------
 template <class ReturnType>
 struct ReduceSum
 {
     ReturnType operator()(const ReturnType &x, const ReturnType &y) const { return x + y; };
 };
-/** A Functor for Maximum */
+
 struct ReduceMax
 {
     Real operator()(Real x, Real y) const { return SMAX(x, y); };
 };
-/** A Functor for Minimum */
+
 struct ReduceMin
 {
     Real operator()(Real x, Real y) const { return SMIN(x, y); };
 };
-/** A Functor for OR operator */
+
 struct ReduceOR
 {
     bool operator()(bool x, bool y) const { return x || y; };
 };
-/** A Functor for AND operator */
+
 struct ReduceAND
 {
     bool operator()(bool x, bool y) const { return x && y; };
 };
-/** A Functor for lower bound */
+
 struct ReduceLowerBound
 {
     Vecd operator()(const Vecd &x, const Vecd &y) const
@@ -73,7 +117,6 @@ struct ReduceLowerBound
         return lower_bound;
     };
 };
-/** A Functor for upper bound */
 struct ReduceUpperBound
 {
     Vecd operator()(const Vecd &x, const Vecd &y) const

src/shared/particle_dynamics/fluid_dynamics/all_fluid_dynamics.h

@@ -40,3 +40,4 @@
 #include "fluid_dynamics_multi_phase.hpp"
 #include "fluid_surface_complex.h"
 #include "fluid_surface_inner.hpp"
+#include "transport_velocity_correction.h"

src/shared/particle_dynamics/fluid_dynamics/base_fluid_dynamics.h

@@ -0,0 +1,47 @@
+/* ------------------------------------------------------------------------- *
+ *                                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.        *
+ *                                                                           *
+ * ------------------------------------------------------------------------- */
+/**
+ * @file base_fluid_dynamics.h
+ * @brief Collection of headers and types used by all fluid dynamics classes.
+ * @author Xiangyu Hu
+ */
+
+#ifndef BASE_FLUID_DYNAMICS_H
+#define BASE_FLUID_DYNAMICS_H
+
+#include "all_body_relations.h"
+#include "all_particle_dynamics.h"
+#include "base_particles.hpp"
+#include "fluid_body.h"
+
+namespace SPH
+{
+namespace fluid_dynamics
+{
+typedef DataDelegateSimple<BaseParticles> FluidDataSimple;
+typedef DataDelegateInner<BaseParticles> FluidDataInner;
+typedef DataDelegateContact<BaseParticles, BaseParticles> FluidContactData;
+typedef DataDelegateContact<BaseParticles, BaseParticles, DataDelegateEmptyBase> FluidContactOnly;
+} // namespace fluid_dynamics
+} // namespace SPH
+#endif // BASE_FLUID_DYNAMICS_H

src/shared/particle_dynamics/fluid_dynamics/fluid_boundary.h

@@ -126,21 +126,21 @@ class FreeStreamVelocityCorrection : public LocalDynamics, public FluidDataSimpl
     Real rho0_;
     StdLargeVec<Real> &rho_sum_;
     StdLargeVec<Vecd> &pos_, &vel_;
-    StdLargeVec<int> &surface_indicator_;
+    StdLargeVec<int> &indicator_;
     TargetVelocity target_velocity;
 
   public:
     explicit FreeStreamVelocityCorrection(SPHBody &sph_body, const Transform &transform = Transform())
         : LocalDynamics(sph_body), FluidDataSimple(sph_body),
           transform_(transform), rho0_(DynamicCast<Fluid>(this, particles_->getBaseMaterial()).ReferenceDensity()),
           rho_sum_(*particles_->getVariableByName<Real>("DensitySummation")), pos_(particles_->pos_), vel_(particles_->vel_),
-          surface_indicator_(*particles_->getVariableByName<int>("SurfaceIndicator")),
+          indicator_(*particles_->getVariableByName<int>("Indicator")),
           target_velocity(*this){};
     virtual ~FreeStreamVelocityCorrection(){};
 
     void update(size_t index_i, Real dt = 0.0)
     {
-        if (surface_indicator_[index_i] == 1)
+        if (indicator_[index_i] == 1)
         {
             Vecd frame_position = transform_.shiftBaseStationToFrame(pos_[index_i]);
             Vecd frame_velocity = transform_.xformBaseVecToFrame(vel_[index_i]);

src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.h

@@ -31,6 +31,8 @@
 #ifndef FLUID_DYNAMICS_COMPLEX_H
 #define FLUID_DYNAMICS_COMPLEX_H
 
+#include "base_fluid_dynamics.h"
+
 #include "fluid_dynamics_inner.h"
 #include "fluid_dynamics_inner.hpp"
 #include "solid_body.h"
@@ -40,10 +42,7 @@ namespace SPH
 {
 namespace fluid_dynamics
 {
-typedef DataDelegateContact<BaseParticles, SolidParticles, DataDelegateEmptyBase>
-    FluidWallData;
-typedef DataDelegateContact<BaseParticles, BaseParticles, DataDelegateEmptyBase>
-    FluidContactData;
+typedef DataDelegateContact<BaseParticles, SolidParticles, DataDelegateEmptyBase> FluidWallData;
 typedef DataDelegateContact<BaseParticles, SolidParticles> FSIContactData;
 /**
  * @class InteractionWithWall
@@ -74,7 +73,7 @@ class InteractionWithWall : public BaseIntegrationType, public FluidWallData
  */
 template <class DensitySummationInnerType>
 class BaseDensitySummationComplex
-    : public BaseInteractionComplex<DensitySummationInnerType, FluidContactData>
+    : public BaseInteractionComplex<DensitySummationInnerType, FluidContactOnly>
 {
   public:
     template <typename... Args>
@@ -138,40 +137,6 @@ class BaseViscousAccelerationWithWall : public InteractionWithWall<ViscousAccele
 
 using ViscousAccelerationWithWall = BaseViscousAccelerationWithWall<ViscousAccelerationInner>;
 
-/**
- * @class TransportVelocityCorrectionComplex
- * @brief  transport velocity correction considering the contribution from contact bodies
- */
-class TransportVelocityCorrectionComplex
-    : public BaseInteractionComplex<TransportVelocityCorrectionInner, FluidContactData>
-{
-  public:
-    template <typename... Args>
-    TransportVelocityCorrectionComplex(Args &&...args)
-        : BaseInteractionComplex<TransportVelocityCorrectionInner, FluidContactData>(
-              std::forward<Args>(args)...){};
-    virtual ~TransportVelocityCorrectionComplex(){};
-
-    inline void interaction(size_t index_i, Real dt = 0.0);
-};
-
-/**
- * @class TransportVelocityCorrectionComplexAdaptive
- * @brief  transport velocity correction considering the contribution from contact bodies
- */
-class TransportVelocityCorrectionComplexAdaptive
-    : public BaseInteractionComplex<TransportVelocityCorrectionInnerAdaptive, FluidContactData>
-{
-  public:
-    template <typename... Args>
-    TransportVelocityCorrectionComplexAdaptive(Args &&...args)
-        : BaseInteractionComplex<TransportVelocityCorrectionInnerAdaptive, FluidContactData>(
-              std::forward<Args>(args)...){};
-    virtual ~TransportVelocityCorrectionComplexAdaptive(){};
-
-    inline void interaction(size_t index_i, Real dt = 0.0);
-};
-
 /**
  * @class BaseIntegration1stHalfWithWall
  * @brief  template class pressure relaxation scheme together with wall boundary

src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.hpp

@@ -32,55 +32,6 @@ void DensitySummationComplexAdaptive::
                          sph_adaptation_.NumberDensityScaleFactor(h_ratio_[index_i]);
 }
 //=================================================================================================//
-void TransportVelocityCorrectionComplex::
-    interaction(size_t index_i, Real dt)
-{
-    TransportVelocityCorrectionInner::interaction(index_i, dt);
-
-    Vecd acceleration_trans = Vecd::Zero();
-    for (size_t k = 0; k < contact_configuration_.size(); ++k)
-    {
-        Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i];
-        for (size_t n = 0; n != contact_neighborhood.current_size_; ++n)
-        {
-            Vecd nablaW_ijV_j = contact_neighborhood.dW_ijV_j_[n] * contact_neighborhood.e_ij_[n];
-
-            // acceleration for transport velocity
-            acceleration_trans -= 2.0 * nablaW_ijV_j;
-        }
-    }
-
-    /** correcting particle position */
-    if (surface_indicator_[index_i] == 0)
-        pos_[index_i] += coefficient_ * smoothing_length_sqr_ * acceleration_trans;
-}
-//=================================================================================================//
-void TransportVelocityCorrectionComplexAdaptive::
-    interaction(size_t index_i, Real dt)
-{
-    TransportVelocityCorrectionInnerAdaptive::interaction(index_i, dt);
-
-    Vecd acceleration_trans = Vecd::Zero();
-    for (size_t k = 0; k < contact_configuration_.size(); ++k)
-    {
-        Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i];
-        for (size_t n = 0; n != contact_neighborhood.current_size_; ++n)
-        {
-            Vecd nablaW_ijV_j = contact_neighborhood.dW_ijV_j_[n] * contact_neighborhood.e_ij_[n];
-
-            // acceleration for transport velocity
-            acceleration_trans -= 2.0 * nablaW_ijV_j;
-        }
-    }
-
-    /** correcting particle position */
-    if (surface_indicator_[index_i] == 0)
-    {
-        Real inv_h_ratio = 1.0 / sph_adaptation_.SmoothingLengthRatio(index_i);
-        pos_[index_i] += coefficient_ * smoothing_length_sqr_ * inv_h_ratio * inv_h_ratio * acceleration_trans;
-    }
-}
-//=================================================================================================//
 void Oldroyd_BIntegration1stHalfWithWall::
     interaction(size_t index_i, Real dt)
 {
@@ -160,7 +111,7 @@ template <class DensitySummationInnerType>
 template <typename... Args>
 BaseDensitySummationComplex<DensitySummationInnerType>::
     BaseDensitySummationComplex(Args &&...args)
-    : BaseInteractionComplex<DensitySummationInnerType, FluidContactData>(std::forward<Args>(args)...)
+    : BaseInteractionComplex<DensitySummationInnerType, FluidContactOnly>(std::forward<Args>(args)...)
 {
     for (size_t k = 0; k != this->contact_particles_.size(); ++k)
     {

src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.cpp

@@ -42,21 +42,6 @@ BaseViscousAccelerationInner::BaseViscousAccelerationInner(BaseInnerRelation &in
       mu_(DynamicCast<Fluid>(this, particles_->getBaseMaterial()).ReferenceViscosity()),
       smoothing_length_(sph_body_.sph_adaptation_->ReferenceSmoothingLength()) {}
 //=================================================================================================//
-TransportVelocityCorrectionInner::
-    TransportVelocityCorrectionInner(BaseInnerRelation &inner_relation, Real coefficient)
-    : LocalDynamics(inner_relation.getSPHBody()), FluidDataInner(inner_relation),
-      pos_(particles_->pos_), surface_indicator_(*particles_->getVariableByName<int>("SurfaceIndicator")),
-      smoothing_length_sqr_(pow(sph_body_.sph_adaptation_->ReferenceSmoothingLength(), 2)),
-      coefficient_(coefficient) {}
-//=================================================================================================//
-TransportVelocityCorrectionInnerAdaptive::
-    TransportVelocityCorrectionInnerAdaptive(BaseInnerRelation &inner_relation, Real coefficient)
-    : LocalDynamics(inner_relation.getSPHBody()), FluidDataInner(inner_relation),
-      sph_adaptation_(*sph_body_.sph_adaptation_),
-      pos_(particles_->pos_), surface_indicator_(*particles_->getVariableByName<int>("SurfaceIndicator")),
-      smoothing_length_sqr_(pow(sph_body_.sph_adaptation_->ReferenceSmoothingLength(), 2)),
-      coefficient_(coefficient) {}
-//=================================================================================================//
 AcousticTimeStepSize::AcousticTimeStepSize(SPHBody &sph_body, Real acousticCFL)
     : LocalDynamicsReduce<Real, ReduceMax>(sph_body, Real(0)),
       FluidDataSimple(sph_body), fluid_(DynamicCast<Fluid>(this, particles_->getBaseMaterial())),