[Study] Shell-to-solid coupling

src/shared/body_relations/contact_body_relation.cpp

@@ -250,4 +250,30 @@ void SurfaceContactRelation::updateConfiguration()
     }
 }
 //=================================================================================================//
+MaxSmoothingLengthContactRelation::MaxSmoothingLengthContactRelation(SPHBody &sph_body,
+                                                                     RealBodyVector contact_bodies,
+                                                                     const std::vector<Real> &factors)
+    : ContactRelationCrossResolution(sph_body, contact_bodies)
+{
+    for (size_t k = 0; k != contact_bodies_.size(); ++k)
+    {
+        Real factor = factors.empty() ? 1 : factors[k];
+        get_contact_neighbors_.push_back(
+            neighbor_builder_contact_ptrs_keeper_.createPtr<MaxSmoothingLengthNeighborBuilder>(
+                sph_body_, *contact_bodies_[k], factor));
+    }
+}
+//=================================================================================================//
+void MaxSmoothingLengthContactRelation::updateConfiguration()
+{
+    resetNeighborhoodCurrentSize();
+    for (size_t k = 0; k != contact_bodies_.size(); ++k)
+    {
+        Mesh &mesh = target_cell_linked_lists_[k]->getMesh();
+        target_cell_linked_lists_[k]->searchNeighborsByMesh(
+            mesh, 0, sph_body_, contact_configuration_[k],
+            *get_search_depths_[k], *get_contact_neighbors_[k]);
+    }
+}
+//=================================================================================================//
 } // namespace SPH

src/shared/body_relations/contact_body_relation.h

@@ -60,7 +60,7 @@ class ContactRelationCrossResolution : public BaseContactRelation
                     sph_body_, target_cell_linked_list->getMesh()));
         }
     };
-    virtual ~ContactRelationCrossResolution(){};
+    virtual ~ContactRelationCrossResolution() {};
     StdVec<CellLinkedList *> getContactCellLinkedList() { return target_cell_linked_lists_; }
 
   protected:
@@ -79,7 +79,7 @@ class ContactRelation : public ContactRelationCrossResolution
 
   public:
     ContactRelation(SPHBody &sph_body, RealBodyVector contact_bodies);
-    virtual ~ContactRelation(){};
+    virtual ~ContactRelation() {};
     virtual void updateConfiguration() override;
 
   protected:
@@ -103,8 +103,8 @@ class ShellSurfaceContactRelation : public ContactRelationCrossResolution
     ShellSurfaceContactRelation(SPHBody &sph_body, RealBodyVector contact_bodies);
     ShellSurfaceContactRelation(SelfSurfaceContactRelation &solid_body_relation_self_contact,
                                 RealBodyVector contact_bodies)
-        : ShellSurfaceContactRelation(*solid_body_relation_self_contact.real_body_, contact_bodies){};
-    virtual ~ShellSurfaceContactRelation(){};
+        : ShellSurfaceContactRelation(*solid_body_relation_self_contact.real_body_, contact_bodies) {};
+    virtual ~ShellSurfaceContactRelation() {};
     virtual void updateConfiguration() override;
 
   protected:
@@ -127,7 +127,7 @@ class ContactRelationToBodyPart : public ContactRelationCrossResolution
     StdVec<NeighborBuilderContactBodyPart *> get_part_contact_neighbors_;
 
     ContactRelationToBodyPart(SPHBody &sph_body, BodyPartVector contact_body_parts_);
-    virtual ~ContactRelationToBodyPart(){};
+    virtual ~ContactRelationToBodyPart() {};
 
     virtual void updateConfiguration() override;
 };
@@ -149,7 +149,7 @@ class AdaptiveContactRelation : public BaseContactRelation
 
   public:
     AdaptiveContactRelation(SPHBody &body, RealBodyVector contact_bodies);
-    virtual ~AdaptiveContactRelation(){};
+    virtual ~AdaptiveContactRelation() {};
 
     virtual void updateConfiguration() override;
 };
@@ -211,7 +211,7 @@ class SurfaceContactRelation : public ContactRelationCrossResolution
                            StdVec<bool> normal_corrections = {})
         : SurfaceContactRelation(*solid_body_relation_self_contact.real_body_,
                                  std::move(contact_bodies),
-                                 std::move(normal_corrections)){};
+                                 std::move(normal_corrections)) {};
     ~SurfaceContactRelation() override = default;
     void updateConfiguration() override;
     BodySurfaceLayer *get_body_surface_layer() { return body_surface_layer_; }
@@ -223,5 +223,22 @@ class SurfaceContactRelation : public ContactRelationCrossResolution
 
     void resetNeighborhoodCurrentSize() override;
 };
+
+/**
+ * @class MaxSmoothingLengthContactRelation
+ * @brief The relation between a SPH body and its contact SPH bodies, where the cut-off radius is set as factor * 2h_max
+ */
+class MaxSmoothingLengthContactRelation : public ContactRelationCrossResolution
+{
+  private:
+    UniquePtrsKeeper<MaxSmoothingLengthNeighborBuilder> neighbor_builder_contact_ptrs_keeper_;
+
+  public:
+    MaxSmoothingLengthContactRelation(SPHBody &sph_body, RealBodyVector contact_bodies, const std::vector<Real> &factors = {});
+    void updateConfiguration() override;
+
+  private:
+    StdVec<MaxSmoothingLengthNeighborBuilder *> get_contact_neighbors_;
+};
 } // namespace SPH
 #endif // CONTACT_BODY_RELATION_H

src/shared/particle_dynamics/solid_dynamics/all_solid_dynamics.h

@@ -37,5 +37,7 @@
 #include "inelastic_dynamics.h"
 #include "loading_dynamics.h"
 #include "solid_dynamics_variable.h"
+#include "solid_to_shell_coupling.h"
+#include "solid_to_shell_coupling.hpp"
 #include "thin_structure_dynamics.h"
 #include "thin_structure_math.h"

src/shared/particle_dynamics/solid_dynamics/solid_to_shell_coupling.h

@@ -0,0 +1,104 @@
+/* ------------------------------------------------------------------------- *
+ *                                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 solid_to_shell_coupling.h
+ * @brief Constraints for solid to shell coupling.
+ * @author	Xiangyu Hu, Weiyi Kong (Virtonomy GmbH)
+ */
+
+#ifndef SOLID_TO_SHELL_CONSTRAINT_H
+#define SOLID_TO_SHELL_CONSTRAINT_H
+
+#include "force_prior.h"
+#include "general_constraint.h"
+
+namespace SPH
+{
+namespace solid_dynamics
+{
+/**@class ComputeWeight
+ * @brief Compute the total weight for each solid particle from the coupled shell particles.
+ * This weight is used for interpolation of velocity and force
+ */
+template <class DynamicsIdentifier>
+class TotalWeightComputation : public BaseLocalDynamics<DynamicsIdentifier>,
+                               public DataDelegateContact
+{
+  private:
+    Real *total_weight_; // \sum_i W_ij * Vol_j for j of all contact bodies
+    StdVec<Real *> contact_Vol_;
+
+  public:
+    TotalWeightComputation(DynamicsIdentifier &identifier, BaseContactRelation &contact_relation);
+    void update(size_t index_i, Real dt = 0.0);
+};
+
+/**@class ConsistentMapping
+ * @brief Constrain the variable by the interpolation from coupled particle velocities with consistent mapping.
+ * Usually used for velocity.
+ * When there is more than one contact body, interpolate from all contact bodies.
+ */
+template <class DynamicsIdentifier, typename DataType>
+class ConsistentMapping : public BaseLocalDynamics<DynamicsIdentifier>,
+                          public DataDelegateContact
+{
+  private:
+    DataType *interpolated_quantities_;
+    Real *total_weight_;
+    StdVec<Real *> contact_Vol_;
+    StdVec<DataType *> contact_data_;
+
+  public:
+    ConsistentMapping(DynamicsIdentifier &identifier,
+                      BaseContactRelation &contact_relation,
+                      const std::string &variable_name,          // variable name of the source body, e.g., "Velocity"
+                      const std::string &contact_variable_name); // variable name of the contact body, e.g., "Velocity"
+    void update(size_t index_i, Real dt = 0.0);
+};
+
+/**@class ConservativeMapping
+ * @brief Constrain the variable by the interpolation from coupled particle with conservative mapping.
+ * Usually used for force. Note that for prior force, the force needs to be registered as a new prior force.
+ * When there is more than one contact body, interpolate force from each contact body and sum up.
+ */
+template <class DynamicsIdentifier, typename DataType>
+class ConservativeMapping : public BaseLocalDynamics<DynamicsIdentifier>,
+                            public DataDelegateContact
+{
+  private:
+    DataType *interpolated_quantities_;
+    Real *Vol_;
+    StdVec<Real *> contact_total_weight_;
+    StdVec<DataType *> contact_data_;
+
+  public:
+    ConservativeMapping(DynamicsIdentifier &identifier,
+                        BaseContactRelation &contact_relation,
+                        const std::string &variable_name,          // variable name of the source body, e.g., "CouplingForce"
+                        const std::string &contact_variable_name); // variable name of the contact body, e.g., "Force"
+    void update(size_t index_i, Real dt = 0.0);
+};
+} // namespace solid_dynamics
+} // namespace SPH
+
+#endif // SOLID_TO_SHELL_CONSTRAINT_H

src/shared/particle_dynamics/solid_dynamics/solid_to_shell_coupling.hpp

@@ -0,0 +1,121 @@
+#ifndef SOLID_TO_SHELL_CONSTRAINT_HPP
+#define SOLID_TO_SHELL_CONSTRAINT_HPP
+
+#include "solid_to_shell_coupling.h"
+
+namespace SPH
+{
+namespace solid_dynamics
+{
+template <class DynamicsIdentifier>
+TotalWeightComputation<DynamicsIdentifier>::
+    TotalWeightComputation(DynamicsIdentifier &identifier, BaseContactRelation &contact_relation)
+    : BaseLocalDynamics<DynamicsIdentifier>(identifier),
+      DataDelegateContact(contact_relation),
+      total_weight_(this->particles_->template registerStateVariable<Real>("TotalWeight"))
+{
+    for (auto *contact_particle : contact_particles_)
+        contact_Vol_.emplace_back(contact_particle->template getVariableDataByName<Real>("VolumetricMeasure"));
+}
+
+template <class DynamicsIdentifier>
+void TotalWeightComputation<DynamicsIdentifier>::update(size_t index_i, Real dt)
+{
+    Real weight_ttl = 0;
+    for (size_t k = 0; k != contact_configuration_.size(); ++k)
+    {
+        const Real *Vol_k = contact_Vol_[k];
+        const Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i];
+        for (size_t n = 0; n != contact_neighborhood.current_size_; ++n)
+        {
+            size_t index_j = contact_neighborhood.j_[n];
+            Real weight_j = contact_neighborhood.W_ij_[n] * Vol_k[index_j];
+            weight_ttl += weight_j;
+        }
+    }
+    total_weight_[index_i] = weight_ttl;
+}
+//=================================================================================================//
+template <class DynamicsIdentifier, typename DataType>
+ConsistentMapping<DynamicsIdentifier, DataType>::
+    ConsistentMapping(DynamicsIdentifier &identifier,
+                      BaseContactRelation &contact_relation,
+                      const std::string &variable_name,
+                      const std::string &contact_variable_name)
+    : BaseLocalDynamics<DynamicsIdentifier>(identifier),
+      DataDelegateContact(contact_relation),
+      interpolated_quantities_(this->particles_->template registerStateVariable<DataType>(variable_name)),
+      total_weight_(this->particles_->template getVariableDataByName<Real>("TotalWeight"))
+{
+    for (auto *contact_particle : contact_particles_)
+    {
+        contact_Vol_.emplace_back(contact_particle->template getVariableDataByName<Real>("VolumetricMeasure"));
+        contact_data_.emplace_back(contact_particle->template getVariableDataByName<DataType>(contact_variable_name));
+    }
+};
+
+template <class DynamicsIdentifier, typename DataType>
+void ConsistentMapping<DynamicsIdentifier, DataType>::update(size_t index_i, Real)
+{
+    // only consider particles with contact neighbors
+    if (total_weight_[index_i] < TinyReal)
+        return;
+
+    DataType interpolated_quantities = ZeroData<DataType>::value;
+    for (size_t k = 0; k != contact_configuration_.size(); ++k)
+    {
+        const Real *Vol_k = contact_Vol_[k];
+        const DataType *data_k = contact_data_[k];
+        const Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i];
+        for (size_t n = 0; n != contact_neighborhood.current_size_; ++n)
+        {
+            size_t index_j = contact_neighborhood.j_[n];
+            Real weight_j = contact_neighborhood.W_ij_[n] * Vol_k[index_j];
+            interpolated_quantities += weight_j * data_k[index_j];
+        }
+    }
+    interpolated_quantities_[index_i] = interpolated_quantities / total_weight_[index_i];
+}
+//=================================================================================================//
+template <class DynamicsIdentifier, typename DataType>
+ConservativeMapping<DynamicsIdentifier, DataType>::
+    ConservativeMapping(DynamicsIdentifier &identifier,
+                        BaseContactRelation &contact_relation,
+                        const std::string &variable_name,
+                        const std::string &contact_variable_name)
+    : BaseLocalDynamics<DynamicsIdentifier>(identifier),
+      DataDelegateContact(contact_relation),
+      interpolated_quantities_(this->particles_->template registerStateVariable<DataType>(variable_name)),
+      Vol_(this->particles_->template getVariableDataByName<Real>("VolumetricMeasure"))
+{
+    for (auto *contact_particle : contact_particles_)
+    {
+        contact_total_weight_.emplace_back(contact_particle->template getVariableDataByName<Real>("TotalWeight"));
+        contact_data_.emplace_back(contact_particle->template getVariableDataByName<DataType>(contact_variable_name));
+    }
+};
+
+template <class DynamicsIdentifier, typename DataType>
+void ConservativeMapping<DynamicsIdentifier, DataType>::update(size_t index_i, Real)
+{
+    DataType interpolated_quantities = ZeroData<DataType>::value;
+    for (size_t k = 0; k != contact_configuration_.size(); ++k)
+    {
+        const DataType *data_k = contact_data_[k];
+        const Real *total_weight_k = contact_total_weight_[k];
+        const Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i];
+        for (size_t n = 0; n != contact_neighborhood.current_size_; ++n)
+        {
+            size_t index_j = contact_neighborhood.j_[n];
+            if (total_weight_k[index_j] < TinyReal)
+                continue;
+            Real weight_j = contact_neighborhood.W_ij_[n] * Vol_[index_i] / total_weight_k[index_j];
+            interpolated_quantities += weight_j * data_k[index_j];
+        }
+    }
+    interpolated_quantities_[index_i] = interpolated_quantities;
+}
+} // namespace solid_dynamics
+} // namespace SPH
+
+#endif // SOLID_TO_SHELL_CONSTRAINT_HPP

src/shared/particle_neighborhood/neighborhood.cpp

@@ -558,5 +558,15 @@ operator()(Neighborhood &neighborhood, const Vecd &pos_i, size_t index_i, const
     }
 };
 //=================================================================================================//
+MaxSmoothingLengthNeighborBuilder::MaxSmoothingLengthNeighborBuilder(SPHBody &body, SPHBody &contact_body, Real factor)
+    : NeighborBuilderContact(body, contact_body)
+{
+    // create a kernel with a cut-off radius of factor * dp_max
+    Real h_max = SMAX(body.getSPHAdaptation().ReferenceSmoothingLength(), contact_body.getSPHAdaptation().ReferenceSmoothingLength());
+    Real h = factor * h_max;
+    // smoothing length: cut_off_radius/2.0
+    kernel_ = kernel_keeper_.createPtr<KernelWendlandC2>(h);
+}
+//=================================================================================================//
 } // namespace SPH
 //=================================================================================================//