add aortic_valve.cpp

.github/workflows/ci.yml

@@ -38,7 +38,9 @@ jobs:
             pkg-config `# for installing libraries with vcpkg`\
             git \
             cmake \
-            ninja-build
+            libx11-dev \
+            freeglut3-dev \
+            ninja-build 
 
       - uses: hendrikmuhs/ccache-action@v1.2
         with:
@@ -62,7 +64,9 @@ jobs:
             simbody \
             gtest \
             xsimd \
-            pybind11
+            pybind11 \
+            opencascade
+
 
       - name: Generate buildsystem
         run: |
@@ -149,7 +153,8 @@ jobs:
             simbody `
             gtest `
             xsimd `
-            pybind11
+            pybind11 `
+            opencascade
 
       - uses: ilammy/msvc-dev-cmd@v1
 
@@ -245,7 +250,8 @@ jobs:
             gtest \
             simbody \
             xsimd \
-            pybind11
+            pybind11 \
+            opencascade
 
       - name: Generate buildsystem
         run: |

CMakeLists.txt

@@ -63,6 +63,10 @@ target_compile_definitions(sphinxsys_core INTERFACE _SILENCE_CXX17_ITERATOR_BASE
 find_package(Eigen3 CONFIG REQUIRED)
 target_link_libraries(sphinxsys_core INTERFACE Eigen3::Eigen)
 
+# ## OCCT
+find_package(OpenCASCADE CONFIG REQUIRED)
+target_link_libraries(sphinxsys_core INTERFACE TKSTEP  ) 
+
 # ## TBB
 find_package(TBB CONFIG REQUIRED)
 target_link_libraries(sphinxsys_core INTERFACE TBB::tbb TBB::tbbmalloc $<$<PLATFORM_ID:Windows>:TBB::tbbmalloc_proxy>)

src/shared/common/vector_functions.cpp

@@ -3,6 +3,21 @@
 namespace SPH
 {
 	//=================================================================================================//
+     Vec3d OcctToEigen(const gp_Pnt &occt_vector)
+   {
+    return Vec3d(occt_vector.X(), occt_vector.Y(), occt_vector.Z());
+   }
+   //=================================================================================================//
+     gp_Pnt EigenToOcct(const Vec3d &eigen_vector)
+    {
+    return gp_Pnt(eigen_vector[0], eigen_vector[1], eigen_vector[2]);
+    }
+   //=================================================================================================//
+   Vec3d OcctVecToEigen(const gp_Vec &occt_vector)
+   {
+    return Vec3d(occt_vector.X(), occt_vector.Y(), occt_vector.Z());
+   }
+	//=================================================================================================//
 	SimTK::Vec2 EigenToSimTK(const Vec2d &eigen_vector)
 	{
 		return SimTK::Vec2(eigen_vector[0], eigen_vector[1]);

src/shared/common/vector_functions.h

@@ -29,10 +29,13 @@
 #define SMALL_VECTORS_H
 
 #include "base_data_type.h"
+#include "opencascade/gp_Pnt.hxx"
 
 namespace SPH
 {
-
+	 Vec3d OcctToEigen(const gp_Pnt &occt_vector);
+    gp_Pnt EigenToOcct(const Vec3d &eigen_vector);
+    Vec3d OcctVecToEigen(const gp_Vec &occt_vector);
 	SimTK::Vec2 EigenToSimTK(const Vec2d &eigen_vector);
 	SimTK::Vec3 EigenToSimTK(const Vec3d &eigen_vector);
 	Vec2d SimTKToEigen(const SimTK::Vec2 &simTK_vector);

tests/user_examples/ExtraSources/extra_sources_share/relax_dynamics_surface.cpp

@@ -0,0 +1,105 @@
+#include "relax_dynamics.h"
+#include "base_particles.hpp"
+#include "particle_generator_lattice.h"
+#include "level_set_shape.h"
+#include "surface_shape.h"
+#include "relax_dynamics_surface.h"
+//========================================================================================================//
+namespace SPH
+{
+	//=====================================================================================================//
+	namespace relax_dynamics
+	{
+        //=================================================================================================//
+        ShapeSurfaceBounding2::ShapeSurfaceBounding2(RealBody &real_body_)
+            : LocalDynamics(real_body_), RelaxDataDelegateSimple(real_body_), pos_(particles_->pos_) /*, surface_shape_(DynamicCast<SurfaceShape>(this, sph_body_.body_shape_))*/
+        {
+            shape_ = real_body_.body_shape_;
+
+        }
+        //=================================================================================================//
+        void ShapeSurfaceBounding2::update(size_t index_i, Real dt)
+        {
+            pos_[index_i] = shape_->findClosestPoint(pos_[index_i]);
+
+        }
+        //=================================================================================================//
+        RelaxationStepInnerFirstHalf::
+            RelaxationStepInnerFirstHalf(BaseInnerRelation &inner_relation, bool level_set_correction)
+            : BaseDynamics<void>(inner_relation.getSPHBody()), real_body_(inner_relation.real_body_),
+              inner_relation_(inner_relation)
+        {
+            relaxation_acceleration_inner_ =
+                makeUnique<InteractionDynamics<RelaxationAccelerationInner>>(inner_relation);
+        }
+        //=================================================================================================//
+        void RelaxationStepInnerFirstHalf::exec(Real dt)
+        {
+            real_body_->updateCellLinkedList();
+            inner_relation_.updateConfiguration();
+            relaxation_acceleration_inner_->exec();
+        }
+
+        //=================================================================================================//
+        RelaxationStepInnerSecondHalf::
+            RelaxationStepInnerSecondHalf(BaseInnerRelation &inner_relation, bool level_set_correction)
+            : BaseDynamics<void>(inner_relation.getSPHBody()), real_body_(inner_relation.real_body_),
+              get_time_step_square_(*real_body_), update_particle_position_(*real_body_),
+              surface_bounding_(*real_body_)
+        {
+        }
+        //=================================================================================================//
+        void RelaxationStepInnerSecondHalf::exec(Real dt)
+        {
+            Real dt_square = get_time_step_square_.exec();
+            update_particle_position_.exec(dt_square);
+            surface_bounding_.exec();
+        }
+
+		//=================================================================================================//
+        SurfaceNormalDirection::SurfaceNormalDirection(SPHBody &sph_body)
+            : RelaxDataDelegateSimple(sph_body), LocalDynamics(sph_body), 
+            surface_shape_(DynamicCast<SurfaceShape>(this, sph_body.body_shape_)), 
+            pos_(particles_->pos_), n_(*particles_->getVariableByName<Vecd>("NormalDirection")) {}
+
+        //=================================================================================================//
+        void SurfaceNormalDirection::update(size_t index_i, Real dt)
+        {
+            Extrema_ExtAlgo Algo = Extrema_ExtAlgo_Tree;
+            gp_Vec tangent_u, tangent_v;
+            gp_Vec norm;
+            Vecd normal_direction_;
+
+
+            gp_Pnt point1 = EigenToOcct(pos_[index_i]);
+            GeomAPI_ProjectPointOnSurf pointonsurf(point1, surface_shape_->surface_, Algo);
+
+
+            Standard_Real u;
+            Standard_Real v;
+            gp_Pnt Point;
+
+            pointonsurf.Parameters(1, u, v);
+
+            surface_shape_->surface_->D1(u, v, Point, tangent_u, tangent_v);
+            norm = tangent_u.Crossed(tangent_v);
+            normal_direction_ = OcctVecToEigen(norm);
+            n_[index_i] = normal_direction_.normalized();
+        }
+        //=================================================================================================//
+        ConstrainSurfaceBodyRegion::
+            ConstrainSurfaceBodyRegion(BodyPartByParticle &body_part)
+            : BaseLocalDynamics<BodyPartByParticle>(body_part), RelaxDataDelegateSimple(sph_body_),
+              acc_(particles_->acc_)
+        {
+        }
+        //=================================================================================================//
+        void ConstrainSurfaceBodyRegion::update(size_t index_i, Real dt)
+        {
+            acc_[index_i] = Vecd::Zero();
+        }
+        //=================================================================================================//
+    }
+	//=================================================================================================//
+}
+//=================================================================================================//

tests/user_examples/ExtraSources/extra_sources_share/relax_dynamics_surface.h

@@ -0,0 +1,127 @@
+/* -------------------------------------------------------------------------*
+ *								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-2022 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 	relax_dynamics.h
+ * @brief 	This is the classes of particle relaxation in order to produce body fitted
+ * 			initial particle distribution.
+ * @author	Chi Zhang and Xiangyu Hu
+ */
+
+#ifndef RELAX_DYNAMICS_SURFACE_H
+#define RELAX_DYNAMICS_SURFACE_H
+
+#include "all_particle_dynamics.h"
+#include "base_kernel.h"
+#include "cell_linked_list.h"
+#include "all_body_relations.h"
+#include "general_dynamics.h"
+#include "relax_dynamics.h"
+
+namespace SPH
+{
+
+    class SurfaceShape;
+
+	namespace relax_dynamics
+    {
+     class ShapeSurfaceBounding2 : public LocalDynamics,
+                                  public RelaxDataDelegateSimple
+    {
+      public:
+        ShapeSurfaceBounding2(RealBody &real_body_);
+        virtual ~ShapeSurfaceBounding2(){};
+        void update(size_t index_i, Real dt = 0.0);
+
+      protected:
+        StdLargeVec<Vecd> &pos_;
+        Shape *shape_;
+    };
+
+
+      class RelaxationStepInnerFirstHalf : public BaseDynamics<void>
+    {
+      public:
+        explicit RelaxationStepInnerFirstHalf(BaseInnerRelation &inner_relation,
+                                              bool level_set_correction = false);
+        virtual ~RelaxationStepInnerFirstHalf(){};
+        virtual void exec(Real dt = 0.0) override;
+
+
+      protected:
+        RealBody *real_body_;
+        BaseInnerRelation &inner_relation_;
+        UniquePtr<BaseDynamics<void>> relaxation_acceleration_inner_;
+    };
+
+    class RelaxationStepInnerSecondHalf : public BaseDynamics<void>
+    {
+      public:
+        explicit RelaxationStepInnerSecondHalf(BaseInnerRelation &inner_relation,
+                                               bool level_set_correction = false);
+        virtual ~RelaxationStepInnerSecondHalf(){};
+        SimpleDynamics<ShapeSurfaceBounding2> &SurfaceBounding() { return surface_bounding_; };
+        virtual void exec(Real dt = 0.0) override;
+
+
+      protected:
+        RealBody *real_body_;
+        ReduceDynamics<GetTimeStepSizeSquare> get_time_step_square_;
+        SimpleDynamics<UpdateParticlePosition> update_particle_position_;
+        SimpleDynamics<ShapeSurfaceBounding2> surface_bounding_;
+    };
+
+    /**
+     * @class SurfaceNormalDirection
+     * @brief get the normal direction of surface particles.
+     */
+    class SurfaceNormalDirection : public RelaxDataDelegateSimple, public LocalDynamics
+    {
+      public:
+        explicit SurfaceNormalDirection(SPHBody &sph_body);
+        virtual ~SurfaceNormalDirection(){};
+        void update(size_t index_i, Real dt = 0.0);
+
+      protected:
+        SurfaceShape *surface_shape_;
+        StdLargeVec<Vecd> &pos_, &n_;
+    };
+
+
+    /**@class ConstrainSuefaceBodyRegion
+     * @brief Fix the position surafce body part.
+     */
+    class ConstrainSurfaceBodyRegion : public BaseLocalDynamics<BodyPartByParticle>, public RelaxDataDelegateSimple
+    {
+      public:
+        ConstrainSurfaceBodyRegion(BodyPartByParticle &body_part);
+        virtual ~ConstrainSurfaceBodyRegion(){};
+        void update(size_t index_i, Real dt = 0.0);
+
+      protected:
+        StdLargeVec<Vecd> &acc_;
+    };
+
+
+    }
+  }
+#endif // RELAX_DYNAMICS_H

tests/user_examples/ExtraSources/extra_sources_share/surface_shape.cpp

@@ -0,0 +1,58 @@
+#include "surface_shape.h"
+
+namespace SPH
+{
+	//=================================================================================================//
+Vecd SurfaceShape::findClosestPoint(const Vecd &input_pnt)
+{
+
+    Vecd closest_pnt;
+    gp_Pnt point1 = EigenToOcct(input_pnt);
+    gp_Pnt point2;
+    Extrema_ExtAlgo Algo = Extrema_ExtAlgo_Tree;
+    point2 = GeomAPI_ProjectPointOnSurf(point1, surface_, Algo);
+    closest_pnt[0] = point2.X();
+    closest_pnt[1] = point2.Y();
+    closest_pnt[2] = point2.Z();
+
+    return closest_pnt;
+}
+
+////=================================================================================================//
+bool SurfaceShape::checkContain(const Vecd &pnt, bool BOUNDARY_INCLUDED) { return 0; }
+BoundingBox SurfaceShape::findBounds() { return BoundingBox(); }
+//=================================================================================================//
+
+Vecd SurfaceShape::findActualPoint(Standard_Real u, Standard_Real v)
+{
+    gp_Pnt point;
+    Vec3d actual_pnt;
+    point = surface_->Value(u, v);
+    actual_pnt[0] = point.X();
+    actual_pnt[1] = point.Y();
+    actual_pnt[2] = point.Z();
+    return Vecd(actual_pnt[0], actual_pnt[1], actual_pnt[2]);
+}
+
+//=================================================================================================//
+SurfaceShapeSTEP::
+    SurfaceShapeSTEP(Standard_CString &filepathname, const std::string &shape_name)
+    : SurfaceShape(shape_name)
+{
+    STEPControl_Reader step_reader;
+    step_reader.ReadFile(filepathname);
+
+    for (Standard_Integer i = 1; i <= step_reader.NbRootsForTransfer(); i++)
+        step_reader.TransferRoot(i);
+
+    TopoDS_Shape step_shape;
+    for (Standard_Integer i = 1; i <= step_reader.NbShapes(); i++)
+        step_shape = step_reader.Shape(i);
+
+    TopExp_Explorer explorer(step_shape, TopAbs_FACE);
+    TopoDS_Face face = TopoDS::Face(explorer.Current());
+    surface_ = BRep_Tool::Surface(face);
+}
+
+	//=================================================================================================//
+}