1 // =============================================================================
2 // PROJECT CHRONO - http://projectchrono.org
3 //
4 // Copyright (c) 2014 projectchrono.org
5 // All rights reserved.
6 //
7 // Use of this source code is governed by a BSD-style license that can be found
8 // in the LICENSE file at the top level of the distribution and at
9 // http://projectchrono.org/license-chrono.txt.
10 //
11 // =============================================================================
12 // Author: Milad Rakhsha, Arman Pazouki, Wei Hu
13 // =============================================================================
14 //
15 // Implementation of FSI system that includes all subclasses for proximity and
16 // force calculation, and time integration.
17 //
18 // =============================================================================
19 
20 #ifndef CH_SYSTEMFSI_IMPL_H_
21 #define CH_SYSTEMFSI_IMPL_H_
22 
23 #include "chrono/ChConfig.h"
24 #include <thrust/device_vector.h>
25 #include <thrust/host_vector.h>
26 #include <thrust/iterator/detail/normal_iterator.h>
27 #include <thrust/iterator/transform_iterator.h>
28 #include <thrust/iterator/zip_iterator.h>
29 #include <thrust/tuple.h>
30 
31 #include "chrono_fsi/physics/ChFsiGeneral.h"
32 #include "chrono_fsi/math/custom_math.h"
33 #include "chrono_fsi/utils/ChUtilsDevice.cuh"
34 
35 namespace chrono {
36 namespace fsi {
37 
38 /// @addtogroup fsi_physics
39 /// @{
40 
41     /// typedef device iterators for shorthand SPH operation of thrust vectors of Real3
42     typedef thrust::device_vector<Real3>::iterator r3IterD;
43     /// typedef device iterators for shorthand SPH operation of thrust vectors of Real4
44     typedef thrust::device_vector<Real4>::iterator r4IterD;
45     /// typedef device tuple for holding SPH data pos,vel,[rho,pressure,mu,type]
46     typedef thrust::tuple<r4IterD, r3IterD, r4IterD, r3IterD, r3IterD> iterTupleSphD;
47     typedef thrust::zip_iterator<iterTupleSphD> zipIterSphD;
48 
49     /// typedef host iterators for shorthand SPH operation of thrust vectors of Real3
50     typedef thrust::host_vector<Real3>::iterator r3IterH;
51     /// typedef host iterators for shorthand SPH operation of thrust vectors of Real4
52     typedef thrust::host_vector<Real4>::iterator r4IterH;
53     /// typedef host tuple for holding SPH data pos,vel,[rho,pressure,mu,type]
54     typedef thrust::tuple<r4IterH, r3IterH, r4IterH, r3IterH, r3IterH> iterTupleH;
55     typedef thrust::zip_iterator<iterTupleH> zipIterSphH;
56 
57     /// typedef device iterators for shorthand rigid body states:
58     /// pos,orientation in position, velocity and acceleration level
59     typedef thrust::tuple<r3IterD, r4IterD, r3IterD, r4IterD, r3IterD, r3IterD> iterTupleRigidD;
60     typedef thrust::zip_iterator<iterTupleRigidD> zipIterRigidD;
61 
62     /// typedef host iterators for shorthand rigid body states:
63     /// pos,orientation in position, velocity and acceleration level
64     typedef thrust::tuple<r3IterH, r4IterH, r3IterH, r4IterH, r3IterH, r3IterH> iterTupleRigidH;
65     typedef thrust::zip_iterator<iterTupleRigidH> zipIterRigidH;
66 
67     /// typedef device iterators for shorthand chrono bodies operations
68     typedef thrust::tuple<r3IterH, r3IterH, r3IterH, r4IterH, r3IterH, r3IterH> iterTupleChronoBodiesH;
69     typedef thrust::zip_iterator<iterTupleChronoBodiesH> zipIterChronoBodiesH;
70 
71 
72     /// Struct to store the information of SPH particles on the device
73     struct SphMarkerDataD {
74         thrust::device_vector<Real4> posRadD;     ///< Vector of the positions of particles + characteristic radius
75         thrust::device_vector<Real3> velMasD;     ///< Vector of the velocities of particles
76         thrust::device_vector<Real4> rhoPresMuD;  ///< Vector of the rho+pressure+mu+type of particles
77         thrust::device_vector<Real3> tauXxYyZzD;  ///< Vector of the shear stress (diagonal) of particles
78         thrust::device_vector<Real3> tauXyXzYzD;  ///< Vector of the shear stress (off-diagonal) of particles
79 
80         zipIterSphD iterator();
81         void resize(size_t s);
82     };
83 
84     /// Struct to store the information of SPH particles on the host
85     struct SphMarkerDataH {
86         thrust::host_vector<Real4> posRadH;     ///< Vector of the positions of particles
87         thrust::host_vector<Real3> velMasH;     ///< Vector of the velocities of particles
88         thrust::host_vector<Real4> rhoPresMuH;  ///< Vector of the rho+pressure+mu+type of particles
89         thrust::host_vector<Real3> tauXxYyZzH;  ///< Vector of the shear stress (diagonal) of particles
90         thrust::host_vector<Real3> tauXyXzYzH;  ///< Vector of the shear stress (off-diagonal) of particles
91 
92         zipIterSphH iterator();
93         void resize(size_t s);
94     };
95 
96     /// Struct to store the information of rigid bodies on the host
97     struct FsiBodiesDataH {
98         thrust::host_vector<Real3> posRigid_fsiBodies_H;      ///< Vector of the linear positions of rigid bodies
99         thrust::host_vector<Real4> velMassRigid_fsiBodies_H;  ///< Vector of the linear velocities of rigid bodies
100         thrust::host_vector<Real3> accRigid_fsiBodies_H;      ///< Vector of the linear acceleration of rigid bodies
101         thrust::host_vector<Real4> q_fsiBodies_H;             ///< Vector of the orientations (Euler parameters as Quaternion) of rigid bodies
102         thrust::host_vector<Real3> omegaVelLRF_fsiBodies_H;   ///< Vector of the angular velocities of rigid bodies
103         thrust::host_vector<Real3> omegaAccLRF_fsiBodies_H;   ///< Vector of the angular acceleration of rigid bodies
104 
105         zipIterRigidH iterator();
106         void resize(size_t s);
107     };
108 
109     /// Struct to store the information of rigid bodies on the device
110     struct FsiBodiesDataD {
111         thrust::device_vector<Real3> posRigid_fsiBodies_D;      ///< Vector of the linear positions of rigid bodies
112         thrust::device_vector<Real4> velMassRigid_fsiBodies_D;  ///< Vector of the linear velocities of rigid bodies
113         thrust::device_vector<Real3> accRigid_fsiBodies_D;      ///< Vector of the linear acceleration of rigid bodies
114         thrust::device_vector<Real4> q_fsiBodies_D;             ///< Vector of the orientations (Euler parameters as Quaternion) of rigid bodies
115         thrust::device_vector<Real3> omegaVelLRF_fsiBodies_D;   ///< Vector of the angular velocities of rigid bodies
116         thrust::device_vector<Real3> omegaAccLRF_fsiBodies_D;   ///< Vector of the angular acceleration of rigid bodies
117 
118         zipIterRigidD iterator();
119         void CopyFromH(const FsiBodiesDataH& other);
120         FsiBodiesDataD& operator=(const FsiBodiesDataD& other);
121         void resize(size_t s);
122     };
123 
124     /// Struct to store the information of mesh on the host
125     struct FsiMeshDataH {
126         thrust::host_vector<Real3> pos_fsi_fea_H;   ///< Vector of the mesh position
127         thrust::host_vector<Real3> vel_fsi_fea_H;   ///< Vector of the mesh velocity
128         thrust::host_vector<Real3> acc_fsi_fea_H;   ///< Vector of the mesh acceleration
129 
130         // zipIterFlexH iterator();
131         void resize(size_t s);
sizechrono::fsi::FsiMeshDataH132         size_t size() { return pos_fsi_fea_H.size(); };
133     };
134 
135     /// Struct to store the information of mesh on the device
136     struct FsiMeshDataD {
137         thrust::device_vector<Real3> pos_fsi_fea_D; ///< Vector of the mesh position
138         thrust::device_vector<Real3> vel_fsi_fea_D; ///< Vector of the mesh velocity
139         thrust::device_vector<Real3> acc_fsi_fea_D; ///< Vector of the mesh acceleration
140 
141         // zipIterFlexD iterator();
142         void CopyFromH(const FsiMeshDataH& other);
143         FsiMeshDataD& operator=(const FsiMeshDataD& other);
144         void resize(size_t s);
145     };
146 
147     /// Struct to store the information of shell elements on the host
148     struct FsiShellsDataH {
149         thrust::host_vector<Real3> posFlex_fsiBodies_nA_H; ///< Vector of the node A position
150         thrust::host_vector<Real3> posFlex_fsiBodies_nB_H; ///< Vector of the node B position
151         thrust::host_vector<Real3> posFlex_fsiBodies_nC_H; ///< Vector of the node B position
152         thrust::host_vector<Real3> posFlex_fsiBodies_nD_H; ///< Vector of the node D position
153 
154         thrust::host_vector<Real3> velFlex_fsiBodies_nA_H; ///< Vector of the node A velocity
155         thrust::host_vector<Real3> velFlex_fsiBodies_nB_H; ///< Vector of the node B velocity
156         thrust::host_vector<Real3> velFlex_fsiBodies_nC_H; ///< Vector of the node C velocity
157         thrust::host_vector<Real3> velFlex_fsiBodies_nD_H; ///< Vector of the node D velocity
158 
159         thrust::host_vector<Real3> accFlex_fsiBodies_nA_H; ///< Vector of the node A acceleration
160         thrust::host_vector<Real3> accFlex_fsiBodies_nB_H; ///< Vector of the node B acceleration
161         thrust::host_vector<Real3> accFlex_fsiBodies_nC_H; ///< Vector of the node C acceleration
162         thrust::host_vector<Real3> accFlex_fsiBodies_nD_H; ///< Vector of the node D acceleration
163 
164         // zipIterFlexH iterator();
165         void resize(size_t s);
166     };
167 
168     /// Struct to store the information of shell elements on the device
169     struct FsiShellsDataD {
170         thrust::device_vector<Real3> posFlex_fsiBodies_nA_D; ///< Vector of the node A position
171         thrust::device_vector<Real3> posFlex_fsiBodies_nB_D; ///< Vector of the node B position
172         thrust::device_vector<Real3> posFlex_fsiBodies_nC_D; ///< Vector of the node C position
173         thrust::device_vector<Real3> posFlex_fsiBodies_nD_D; ///< Vector of the node D position
174 
175         thrust::device_vector<Real3> velFlex_fsiBodies_nA_D; ///< Vector of the node A velocity
176         thrust::device_vector<Real3> velFlex_fsiBodies_nB_D; ///< Vector of the node B velocity
177         thrust::device_vector<Real3> velFlex_fsiBodies_nC_D; ///< Vector of the node C velocity
178         thrust::device_vector<Real3> velFlex_fsiBodies_nD_D; ///< Vector of the node D velocity
179 
180         thrust::device_vector<Real3> accFlex_fsiBodies_nA_D; ///< Vector of the node A acceleration
181         thrust::device_vector<Real3> accFlex_fsiBodies_nB_D; ///< Vector of the node B acceleration
182         thrust::device_vector<Real3> accFlex_fsiBodies_nC_D; ///< Vector of the node C acceleration
183         thrust::device_vector<Real3> accFlex_fsiBodies_nD_D; ///< Vector of the node D acceleration
184 
185         // zipIterFlexD iterator();
186         void CopyFromH(const FsiShellsDataH& other);
187         FsiShellsDataD& operator=(const FsiShellsDataD& other);
188         void resize(size_t s);
189     };
190 
191     /// Struct to store neighbor search information on the device
192     struct ProximityDataD {
193         thrust::device_vector<uint> gridMarkerHashD;   ///< gridMarkerHash=s(i,j,k)= k*n_x*n_y + j*n_x + i (numAllMarkers);
194         thrust::device_vector<uint> gridMarkerIndexD;  ///< (numAllMarkers);
195         thrust::device_vector<uint> cellStartD;        ///< Index of the particle starts a cell in sorted list (m_numGridCells)
196         thrust::device_vector<uint> cellEndD;          ///< Index of the particle ends a cell in sorted list (m_numGridCells)
197         thrust::device_vector<uint> mapOriginalToSorted; ///< Index mapping from the original to the sorted
198 
199         void resize(size_t numAllMarkers);
200     };
201 
202     /// Struct to store Chrono rigid bodies information on the host
203     struct ChronoBodiesDataH {
ChronoBodiesDataHchrono::fsi::ChronoBodiesDataH204         ChronoBodiesDataH() {}
205         ChronoBodiesDataH(size_t s);
206         thrust::host_vector<Real3> pos_ChSystemH;          ///< Vector of the linear positions of rigid bodies
207         thrust::host_vector<Real3> vel_ChSystemH;          ///< Vector of the linear velocities of rigid bodies
208         thrust::host_vector<Real3> acc_ChSystemH;          ///< Vector of the linear accelerations of rigid bodies
209         thrust::host_vector<Real4> quat_ChSystemH;         ///< Vector of the orientations (Euler parameters as Quaternion) of rigid bodies
210         thrust::host_vector<Real3> omegaVelGRF_ChSystemH;  ///< Vector of the angular velocities of rigid bodies
211         thrust::host_vector<Real3> omegaAccGRF_ChSystemH;  ///< Vector of the angular acceleraion of rigid bodies
212 
213         zipIterChronoBodiesH iterator();
214         void resize(size_t s);
215     };
216 
217     /// Struct to store Chrono shell elements information on the host
218     struct ChronoShellsDataH {
ChronoShellsDataHchrono::fsi::ChronoShellsDataH219         ChronoShellsDataH() {}
220         ChronoShellsDataH(size_t s);
221 
222         // zipIterChronoShellsH iterator();
223 
224         thrust::host_vector<Real3> posFlex_ChSystemH_nA_H; ///< Vector of the node A position
225         thrust::host_vector<Real3> posFlex_ChSystemH_nB_H; ///< Vector of the node B position
226         thrust::host_vector<Real3> posFlex_ChSystemH_nC_H; ///< Vector of the node C position
227         thrust::host_vector<Real3> posFlex_ChSystemH_nD_H; ///< Vector of the node D position
228 
229         thrust::host_vector<Real3> velFlex_ChSystemH_nA_H; ///< Vector of the node A velocity
230         thrust::host_vector<Real3> velFlex_ChSystemH_nB_H; ///< Vector of the node B velocity
231         thrust::host_vector<Real3> velFlex_ChSystemH_nC_H; ///< Vector of the node C velocity
232         thrust::host_vector<Real3> velFlex_ChSystemH_nD_H; ///< Vector of the node D velocity
233 
234         thrust::host_vector<Real3> accFlex_ChSystemH_nA_H; ///< Vector of the node A acceleration
235         thrust::host_vector<Real3> accFlex_ChSystemH_nB_H; ///< Vector of the node B acceleration
236         thrust::host_vector<Real3> accFlex_ChSystemH_nC_H; ///< Vector of the node C acceleration
237         thrust::host_vector<Real3> accFlex_ChSystemH_nD_H; ///< Vector of the node D acceleration
238 
239         void resize(size_t s);
240     };
241 
242     /// Struct to store Chrono mesh information on the host
243     struct ChronoMeshDataH {
ChronoMeshDataHchrono::fsi::ChronoMeshDataH244         ChronoMeshDataH() {}
245         ChronoMeshDataH(size_t s);
246 
247         thrust::host_vector<Real3> posFlex_ChSystemH_H; ///< Vector of the mesh position
248         thrust::host_vector<Real3> velFlex_ChSystemH_H; ///< Vector of the mesh velocity
249         thrust::host_vector<Real3> accFlex_ChSystemH_H; ///< Vector of the mesh acceleration
250 
251         void resize(size_t s);
252     };
253 
254     /// Struct to store fluid/granular system information that need to be passed to Chrono
255     struct FsiGeneralData {
256         // ----------------
257         //  host
258         // ----------------
259         // fluidfsiBodiesIndex
260         thrust::host_vector<int4> referenceArray;        ///< Holds information of each phase in the array of SPH particles
261         thrust::host_vector<int4> referenceArray_FEA;    ///< Holds information of each phase in the array of SPH particles for Flexible elements
262         // ----------------
263         //  device
264         // ----------------
265         // fluid
266         thrust::device_vector<Real4> derivVelRhoD;       ///< dv/dt and d(rho)/dt for particles
267         thrust::device_vector<Real4> derivVelRhoD_old;   ///< dv/dt and d(rho)/dt for particles
268         thrust::device_vector<Real3> derivTauXxYyZzD;    ///< d(tau)/dt for particles
269         thrust::device_vector<Real3> derivTauXyXzYzD;    ///< d(tau)/dt for particles
270         thrust::device_vector<Real3> vel_XSPH_D;         ///< XSPH velocity for particles
271         thrust::device_vector<Real3> vis_vel_SPH_D;      ///< IISPH velocity for particles
272         thrust::device_vector<Real4> sr_tau_I_mu_i;      ///< I2SPH strain-rate, stress, inertia number, friction
273 
274         // BCE
275         thrust::device_vector<Real3> rigidSPH_MeshPos_LRF_D;  ///< Position of a particle attached to a rigid body in a local
276         thrust::device_vector<Real3> FlexSPH_MeshPos_LRF_D;   ///< Position of a particle attached to a mesh in a local on device
277         thrust::host_vector<Real3> FlexSPH_MeshPos_LRF_H;     ///< Position of a particle attached to a mesh in a local on host
278 
279         thrust::device_vector<uint> rigidIdentifierD;         ///< Identifies which rigid body a particle belongs to
280         thrust::device_vector<uint> FlexIdentifierD;          ///< Identifies which flexible body a particle belongs to
281 
282         // fsi bodies
283         thrust::device_vector<Real3> rigid_FSI_ForcesD;       ///< Vector of the surface-integrated forces to rigid bodies
284         thrust::device_vector<Real3> rigid_FSI_TorquesD;      ///< Vector of the surface-integrated torques to rigid bodies
285         thrust::device_vector<Real3> Flex_FSI_ForcesD;        ///< Vector of the surface-integrated force on FEA nodes
286 
287         thrust::host_vector<int2> CableElementsNodesH;        ///< Vector of the cable elements dodes on host
288         thrust::device_vector<int2> CableElementsNodes;       ///< Vector of the cable elements dodes on device
289 
290         thrust::host_vector<int4> ShellElementsNodesH;        ///< Vector of the shell elements dodes on host
291         thrust::device_vector<int4> ShellElementsNodes;       ///< Vector of the shell elements dodes on device
292     };
293 
294 /// @brief Data related function implementations for FSI system
295 class ChSystemFsi_impl : public ChFsiGeneral {
296   public:
297     ChSystemFsi_impl();
298     virtual ~ChSystemFsi_impl();
299 
300     /// Add an SPH particle given its position, physical properties, velocity, and stress
301     void AddSphMarker(Real4 pos, Real4 rhoPresMu, Real3 vel = mR3(0.0), Real3 tauXxYyZz = mR3(0.0), Real3 tauXyXzYz = mR3(0.0));
302 
303     /// Resize the simulation data based on the FSI system constructed.
304     void ResizeDataManager(int numNode = 0);
305 
306     std::shared_ptr<NumberOfObjects> numObjects;        ///< Number of objects (SPH particles, BCE particles, rigid bodies and such).
307 
308     std::shared_ptr<SphMarkerDataD> sphMarkersD1;       ///< Information of SPH particles at state 1 on device
309     std::shared_ptr<SphMarkerDataD> sphMarkersD2;       ///< Information of SPH particles at state 2 on device
310     std::shared_ptr<SphMarkerDataD> sortedSphMarkersD;  ///< Sorted information of SPH particles at state 1 on device
311     std::shared_ptr<SphMarkerDataH> sphMarkersH;        ///< Information of SPH particles on host
312 
313     std::shared_ptr<FsiBodiesDataD> fsiBodiesD1;  ///< Information of rigid bodies at state 1 on device
314     std::shared_ptr<FsiBodiesDataD> fsiBodiesD2;  ///< Information of rigid bodies at state 2 on device
315     std::shared_ptr<FsiBodiesDataH> fsiBodiesH;   ///< Information of rigid bodies at state 1 on host
316     std::shared_ptr<FsiMeshDataD> fsiMeshD;       ///< Information of mesh on device
317     std::shared_ptr<FsiMeshDataH> fsiMeshH;       ///< Information of mesh on host
318 
319     std::shared_ptr<FsiGeneralData> fsiGeneralData;  ///< General FSI data needed in the simulation
320 
321     std::shared_ptr<ProximityDataD> markersProximityD; ///< Information of neighbor search on the device
322 
323   private:
324     void ArrangeDataManager();
325     void ConstructReferenceArray();
326     void InitNumObjects();
327     void CalcNumObjects();  ///< Calculates the number of rigid bodies, flexible bodies, etc. based on the type of particles
328 
329   public:
330     /// Base class of FSI system.
331     friend class ChSystemFsi;
332 
333 };
334 
335 /// @} fsi_physics
336 
337 }  // end namespace fsi
338 }  // end namespace chrono
339 
340 #endif
341