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
13 // =============================================================================
14 //
15 // Base class for processing the interface between Chrono and fsi modules
16 // =============================================================================
17
18 #include "chrono_fsi/ChFsiInterface.h"
19 #include "chrono_fsi/utils/ChUtilsDevice.cuh"
20 #include "chrono_fsi/utils/ChUtilsTypeConvert.h"
21 #include "chrono/fea/ChElementCableANCF.h"
22 #include "chrono/fea/ChElementShellANCF_3423.h"
23 #include "chrono/fea/ChMesh.h"
24 #include "chrono/fea/ChNodeFEAxyzD.h"
25
26 namespace chrono {
27 namespace fsi {
28 //------------------------------------------------------------------------------------
ChFsiInterface(ChSystem & other_mphysicalSystem,std::shared_ptr<fea::ChMesh> other_fsiMesh,std::shared_ptr<SimParams> other_paramsH,std::shared_ptr<FsiBodiesDataH> other_fsiBodiesH,std::shared_ptr<FsiMeshDataH> other_fsiMeshH,std::vector<std::shared_ptr<ChBody>> & other_fsiBodies,std::vector<std::shared_ptr<fea::ChNodeFEAxyzD>> & other_fsiNodes,std::vector<std::shared_ptr<fea::ChElementCableANCF>> & other_fsiCables,std::vector<std::shared_ptr<fea::ChElementShellANCF_3423>> & other_fsiShells,thrust::host_vector<int2> & other_CableElementsNodesH,thrust::device_vector<int2> & other_CableElementsNodes,thrust::host_vector<int4> & other_ShellElementsNodesH,thrust::device_vector<int4> & other_ShellElementsNodes,thrust::device_vector<Real3> & other_rigid_FSI_ForcesD,thrust::device_vector<Real3> & other_rigid_FSI_TorquesD,thrust::device_vector<Real3> & other_Flex_FSI_ForcesD)29 ChFsiInterface::ChFsiInterface(ChSystem& other_mphysicalSystem,
30 std::shared_ptr<fea::ChMesh> other_fsiMesh,
31 std::shared_ptr<SimParams> other_paramsH,
32 std::shared_ptr<FsiBodiesDataH> other_fsiBodiesH,
33 std::shared_ptr<FsiMeshDataH> other_fsiMeshH,
34 std::vector<std::shared_ptr<ChBody>>& other_fsiBodies,
35 std::vector<std::shared_ptr<fea::ChNodeFEAxyzD>>& other_fsiNodes,
36 std::vector<std::shared_ptr<fea::ChElementCableANCF>>& other_fsiCables,
37 std::vector<std::shared_ptr<fea::ChElementShellANCF_3423>>& other_fsiShells,
38 thrust::host_vector<int2>& other_CableElementsNodesH,
39 thrust::device_vector<int2>& other_CableElementsNodes,
40 thrust::host_vector<int4>& other_ShellElementsNodesH,
41 thrust::device_vector<int4>& other_ShellElementsNodes,
42 thrust::device_vector<Real3>& other_rigid_FSI_ForcesD,
43 thrust::device_vector<Real3>& other_rigid_FSI_TorquesD,
44 thrust::device_vector<Real3>& other_Flex_FSI_ForcesD)
45 : mphysicalSystem(other_mphysicalSystem),
46 fsi_mesh(other_fsiMesh),
47 paramsH(other_paramsH),
48 fsiBodiesH(other_fsiBodiesH),
49 fsiMeshH(other_fsiMeshH),
50 fsiBodies(other_fsiBodies),
51 fsiNodes(other_fsiNodes),
52 fsiCables(other_fsiCables),
53 fsiShells(other_fsiShells),
54 CableElementsNodesH(other_CableElementsNodesH),
55 CableElementsNodes(other_CableElementsNodes),
56 ShellElementsNodesH(other_ShellElementsNodesH),
57 ShellElementsNodes(other_ShellElementsNodes),
58 rigid_FSI_ForcesD(other_rigid_FSI_ForcesD),
59 rigid_FSI_TorquesD(other_rigid_FSI_TorquesD),
60 Flex_FSI_ForcesD(other_Flex_FSI_ForcesD) {
61 size_t numBodies = mphysicalSystem.Get_bodylist().size();
62 chronoRigidBackup = chrono_types::make_shared<ChronoBodiesDataH>(numBodies);
63 chronoFlexMeshBackup = chrono_types::make_shared<ChronoMeshDataH>(0);
64 // int numShells = 0;
65 // int numCables = 0;
66 // if (mphysicalSystem.Get_otherphysicslist().size())
67 // numShells =
68 // std::dynamic_pointer_cast<fea::ChMesh>(mphysicalSystem.Get_otherphysicslist().at(0))->GetNelements();
69 int numNodes = 0;
70
71 if (mphysicalSystem.Get_otherphysicslist().size())
72 numNodes = std::dynamic_pointer_cast<fea::ChMesh>(mphysicalSystem.Get_otherphysicslist().at(0))->GetNnodes();
73 chronoFlexMeshBackup = chrono_types::make_shared<ChronoMeshDataH>(numNodes);
74 }
75 //------------------------------------------------------------------------------------
~ChFsiInterface()76 ChFsiInterface::~ChFsiInterface() {}
77 //------------------------------------------------------------------------------------
78
Add_Rigid_ForceTorques_To_ChSystem()79 void ChFsiInterface::Add_Rigid_ForceTorques_To_ChSystem() {
80 size_t numRigids = fsiBodies.size();
81 std::string delim = ",";
82 char filename[4096];
83 ChVector<> totalForce(0);
84 ChVector<> totalTorque(0);
85
86 for (size_t i = 0; i < numRigids; i++) {
87 ChVector<> mforce = ChUtilsTypeConvert::Real3ToChVector(ChUtilsDevice::FetchElement(rigid_FSI_ForcesD, i));
88 ChVector<> mtorque = ChUtilsTypeConvert::Real3ToChVector(ChUtilsDevice::FetchElement(rigid_FSI_TorquesD, i));
89
90 totalForce += mforce;
91 totalTorque + mtorque;
92 std::shared_ptr<ChBody> body = fsiBodies[i];
93 ChVector<> pos = body->GetPos();
94 ChVector<> vel = body->GetPos_dt();
95 ChQuaternion<> rot = body->GetRot();
96
97 sprintf(filename, "%s/FS_body%zd.csv", paramsH->demo_dir, i);
98 std::ofstream file;
99 if (mphysicalSystem.GetChTime() > 0)
100 file.open(filename, std::fstream::app);
101 else {
102 file.open(filename);
103 file << "Time" << delim << "x" << delim << "y" << delim << "z" << delim << "q0" << delim << "q1" << delim
104 << "q2" << delim << "q3" << delim << "Vx" << delim << "Vy" << delim << "Vz" << delim << "Fx" << delim
105 << "Fy" << delim << "Fz" << delim << "Tx" << delim << "Ty" << delim << "Tz" << std::endl;
106 }
107
108 if (1) // set as 1 if output to file
109 file << mphysicalSystem.GetChTime() << delim << pos.x() << delim << pos.y() << delim << pos.z() << delim
110 << rot.e0() << delim << rot.e1() << delim << rot.e2() << delim << rot.e3() << delim << vel.x() << delim
111 << vel.y() << delim << vel.z() << delim << mforce.x() << delim << mforce.y() << delim << mforce.z()
112 << delim << mtorque.x() << delim << mtorque.y() << delim << mtorque.z() << std::endl;
113
114 if (0) // set as 1 if output to screen
115 std::cout << mphysicalSystem.GetChTime() << delim << pos.x() << delim << pos.y() << delim << pos.z()
116 << delim << rot.e0() << delim << rot.e1() << delim << rot.e2() << delim << rot.e3() << delim
117 << vel.x() << delim << vel.y() << delim << vel.z() << delim << mforce.x() << delim << mforce.y()
118 << delim << mforce.z() << delim << mtorque.x() << delim << mtorque.y() << delim << mtorque.z()
119 << std::endl;
120
121 // note: when this FSI body goes back to Chrono system, the gravity
122 // will be automaticly added. Here only accumulate force from fluid
123 body->Empty_forces_accumulators();
124 body->Accumulate_force(mforce, body->GetPos(), false);
125 body->Accumulate_torque(mtorque, false);
126
127 file.close();
128 }
129 }
130
131 //------------------------------------------------------------------------------------
Copy_External_To_ChSystem()132 void ChFsiInterface::Copy_External_To_ChSystem() {
133 size_t numBodies = mphysicalSystem.Get_bodylist().size();
134 if (chronoRigidBackup->pos_ChSystemH.size() != numBodies) {
135 throw std::runtime_error(
136 "Size of the external data does not match the "
137 "ChSystem; thrown from Copy_External_To_ChSystem "
138 "!\n");
139 }
140
141 for (size_t i = 0; i < numBodies; i++) {
142 auto mBody = mphysicalSystem.Get_bodylist().at(i);
143 mBody->SetPos(ChUtilsTypeConvert::Real3ToChVector(chronoRigidBackup->pos_ChSystemH[i]));
144 mBody->SetPos_dt(ChUtilsTypeConvert::Real3ToChVector(chronoRigidBackup->vel_ChSystemH[i]));
145 mBody->SetPos_dtdt(ChUtilsTypeConvert::Real3ToChVector(chronoRigidBackup->acc_ChSystemH[i]));
146
147 mBody->SetRot(ChUtilsTypeConvert::Real4ToChQuaternion(chronoRigidBackup->quat_ChSystemH[i]));
148 mBody->SetWvel_par(ChUtilsTypeConvert::Real3ToChVector(chronoRigidBackup->omegaVelGRF_ChSystemH[i]));
149 ChVector<> acc = ChUtilsTypeConvert::Real3ToChVector(chronoRigidBackup->omegaAccGRF_ChSystemH[i]);
150 mBody->SetWacc_par(acc);
151 }
152 }
153 //------------------------------------------------------------------------------------
Copy_ChSystem_to_External()154 void ChFsiInterface::Copy_ChSystem_to_External() {
155 size_t numBodies = mphysicalSystem.Get_bodylist().size();
156 auto bodyList = mphysicalSystem.Get_bodylist();
157
158 if (chronoRigidBackup->pos_ChSystemH.size() != numBodies) {
159 throw std::runtime_error(
160 "Size of the external data does not match the "
161 "ChSystem; thrown from Copy_ChSystem_to_External "
162 "!\n");
163 }
164
165 chronoRigidBackup->resize(numBodies);
166 for (size_t i = 0; i < numBodies; i++) {
167 auto mBody = mphysicalSystem.Get_bodylist().at(i);
168 chronoRigidBackup->pos_ChSystemH[i] = ChUtilsTypeConvert::ChVectorToReal3(mBody->GetPos());
169 chronoRigidBackup->vel_ChSystemH[i] = ChUtilsTypeConvert::ChVectorToReal3(mBody->GetPos_dt());
170 chronoRigidBackup->acc_ChSystemH[i] = ChUtilsTypeConvert::ChVectorToReal3(mBody->GetPos_dtdt());
171
172 chronoRigidBackup->quat_ChSystemH[i] = ChUtilsTypeConvert::ChQuaternionToReal4(mBody->GetRot());
173 chronoRigidBackup->omegaVelGRF_ChSystemH[i] = ChUtilsTypeConvert::ChVectorToReal3(mBody->GetWvel_par());
174 chronoRigidBackup->omegaAccGRF_ChSystemH[i] = ChUtilsTypeConvert::ChVectorToReal3(mBody->GetWacc_par());
175 }
176
177 int numNodes = fsi_mesh->GetNnodes();
178 for (size_t i = 0; i < numNodes; i++) {
179 auto node = std::dynamic_pointer_cast<fea::ChNodeFEAxyz>(fsi_mesh->GetNode((unsigned int)i));
180 chronoFlexMeshBackup->posFlex_ChSystemH_H[i] = ChUtilsTypeConvert::ChVectorToReal3(node->GetPos());
181 chronoFlexMeshBackup->velFlex_ChSystemH_H[i] = ChUtilsTypeConvert::ChVectorToReal3(node->GetPos_dt());
182 chronoFlexMeshBackup->accFlex_ChSystemH_H[i] = ChUtilsTypeConvert::ChVectorToReal3(node->GetPos_dtdt());
183 }
184 }
185 //------------------------------------------------------------------------------------
Copy_fsiBodies_ChSystem_to_FluidSystem(std::shared_ptr<FsiBodiesDataD> fsiBodiesD)186 void ChFsiInterface::Copy_fsiBodies_ChSystem_to_FluidSystem(std::shared_ptr<FsiBodiesDataD> fsiBodiesD) {
187 size_t num_fsiBodies_Rigids = fsiBodies.size();
188 for (size_t i = 0; i < num_fsiBodies_Rigids; i++) {
189 std::shared_ptr<ChBody> bodyPtr = fsiBodies[i];
190 fsiBodiesH->posRigid_fsiBodies_H[i] = ChUtilsTypeConvert::ChVectorToReal3(bodyPtr->GetPos());
191 fsiBodiesH->velMassRigid_fsiBodies_H[i] =
192 ChUtilsTypeConvert::ChVectorToReal4(bodyPtr->GetPos_dt(), bodyPtr->GetMass());
193 fsiBodiesH->accRigid_fsiBodies_H[i] = ChUtilsTypeConvert::ChVectorToReal3(bodyPtr->GetPos_dtdt());
194 fsiBodiesH->q_fsiBodies_H[i] = ChUtilsTypeConvert::ChQuaternionToReal4(bodyPtr->GetRot());
195 fsiBodiesH->omegaVelLRF_fsiBodies_H[i] = ChUtilsTypeConvert::ChVectorToReal3(bodyPtr->GetWvel_loc());
196 fsiBodiesH->omegaAccLRF_fsiBodies_H[i] = ChUtilsTypeConvert::ChVectorToReal3(bodyPtr->GetWacc_loc());
197 }
198 fsiBodiesD->CopyFromH(*fsiBodiesH);
199 }
200
201 //------------------------------------------------------------------------------------
ResizeChronoBodiesData()202 void ChFsiInterface::ResizeChronoBodiesData() {
203 size_t numBodies = mphysicalSystem.Get_bodylist().size();
204 chronoRigidBackup->resize(numBodies);
205 }
206
207 //------------------------------------------------------------------------------------
208 //------------------------------------------------------------------------------------
209 //-----------------------Chrono FEA Specifics------------------------------------------
210 //------------------------------------------------------------------------------------
Add_Flex_Forces_To_ChSystem()211 void ChFsiInterface::Add_Flex_Forces_To_ChSystem() {
212 // int numShells = 0;
213 // int numNodes_ChSystem = 0;
214 std::string delim = ",";
215
216 size_t numNodes = fsiNodes.size();
217 ChVector<> total_force(0, 0, 0);
218 for (size_t i = 0; i < numNodes; i++) {
219 ChVector<> mforce = ChUtilsTypeConvert::Real3ToChVector(ChUtilsDevice::FetchElement(Flex_FSI_ForcesD, i));
220 // if (mforce.Length() != 0.0)
221 // printf("mforce= (%.3e,%.3e,%.3e)\n", mforce.x(), mforce.y(), mforce.z());
222 auto node = std::dynamic_pointer_cast<fea::ChNodeFEAxyzD>(fsi_mesh->GetNode((unsigned int)i));
223
224 ChVector<> pos = node->GetPos();
225 ChVector<> vel = node->GetPos_dt();
226 char filename[4096];
227 sprintf(filename, "%s/FS_node%zd.csv", paramsH->demo_dir, i);
228 std::ofstream file;
229 if (mphysicalSystem.GetChTime() > 0)
230 file.open(filename, std::fstream::app);
231 else {
232 file.open(filename);
233 file << "Time" << delim << "x" << delim << "y" << delim << "z" << delim << "Vx" << delim << "Vy" << delim
234 << "Vz" << delim << "Fx" << delim << "Fy" << delim << "Fz" << std::endl;
235 }
236
237 file << mphysicalSystem.GetChTime() << delim << pos.x() << delim << pos.y() << delim << pos.z() << delim
238 << vel.x() << delim << vel.y() << delim << vel.z() << delim << mforce.x() << delim << mforce.y() << delim
239 << mforce.z() << std::endl;
240 file.close();
241
242 // ChVector<> OldForce = node->GetForce();
243 node->SetForce(mforce);
244 }
245
246 // printf("Total Force from the fluid to the solid = (%.3e,%.3e,%.3e)\n", total_force.x(), total_force.y(),
247 // total_force.z());
248 }
249 //------------------------------------------------------------------------------------
250
Copy_fsiNodes_ChSystem_to_FluidSystem(std::shared_ptr<FsiMeshDataD> FsiMeshD)251 void ChFsiInterface::Copy_fsiNodes_ChSystem_to_FluidSystem(std::shared_ptr<FsiMeshDataD> FsiMeshD) {
252 size_t num_fsiNodes_Felx = fsiNodes.size();
253
254 for (size_t i = 0; i < num_fsiNodes_Felx; i++) {
255 auto NodePtr = fsiNodes[i];
256 fsiMeshH->pos_fsi_fea_H[i] = ChUtilsTypeConvert::ChVectorToReal3(NodePtr->GetPos());
257 fsiMeshH->vel_fsi_fea_H[i] = ChUtilsTypeConvert::ChVectorToReal3(NodePtr->GetPos_dt());
258 fsiMeshH->acc_fsi_fea_H[i] = ChUtilsTypeConvert::ChVectorToReal3(NodePtr->GetPos_dtdt());
259 }
260 FsiMeshD->CopyFromH(*fsiMeshH);
261 }
262 //------------------------------------------------------------------------------------
263
ResizeChronoNodesData()264 void ChFsiInterface::ResizeChronoNodesData() {
265 int numNodes = 0;
266 auto my_mesh = chrono_types::make_shared<fea::ChMesh>();
267 if (mphysicalSystem.Get_otherphysicslist().size()) {
268 printf("fsi_mesh.size in ResizeChronNodesData %d\n", numNodes);
269 }
270 numNodes = fsi_mesh->GetNnodes();
271 printf("numNodes in ResizeChronNodesData %d\n", numNodes);
272 chronoFlexMeshBackup->resize(numNodes);
273 }
274
275 //------------------------------------------------------------------------------------
ResizeChronoFEANodesData()276 void ChFsiInterface::ResizeChronoFEANodesData() {
277 int numNodes = 0;
278 auto my_mesh = chrono_types::make_shared<fea::ChMesh>();
279 if (mphysicalSystem.Get_otherphysicslist().size()) {
280 my_mesh = std::dynamic_pointer_cast<fea::ChMesh>(mphysicalSystem.Get_otherphysicslist().at(0));
281 }
282 numNodes = fsi_mesh->GetNnodes();
283 printf("fsi_mesh.size in ResizeChronoFEANodesData %d\n", numNodes);
284 printf("numNodes in ResizeChronoFEANodeData %d\n", numNodes);
285
286 chronoFlexMeshBackup->resize(numNodes);
287 }
288 //------------------------------------------------------------------------------------
289
ResizeChronoCablesData(std::vector<std::vector<int>> CableElementsNodesSTDVector,thrust::host_vector<int2> & CableElementsNodesH)290 void ChFsiInterface::ResizeChronoCablesData(std::vector<std::vector<int>> CableElementsNodesSTDVector,
291 thrust::host_vector<int2>& CableElementsNodesH) {
292 auto my_mesh = chrono_types::make_shared<fea::ChMesh>();
293 if (mphysicalSystem.Get_otherphysicslist().size()) {
294 my_mesh = std::dynamic_pointer_cast<fea::ChMesh>(mphysicalSystem.Get_otherphysicslist().at(0));
295 }
296
297 size_t numCables = 0;
298 for (size_t i = 0; i < fsi_mesh->GetNelements(); i++) {
299 if (std::dynamic_pointer_cast<fea::ChElementCableANCF>(fsi_mesh->GetElement((unsigned int)i)))
300 numCables++;
301 }
302 printf("fsi_mesh.size.shell in ResizeChronoCablesData %zd\n", numCables);
303 printf("numCables in ResizeChronoCablesData %zd\n", numCables);
304 printf("CableElementsNodesSTDVector.size() in ResizeChronoCablesData %zd\n", CableElementsNodesSTDVector.size());
305
306 if (CableElementsNodesSTDVector.size() != numCables) {
307 throw std::runtime_error(
308 "Size of the external data does not match the "
309 "ChSystem; thrown from ChFsiInterface::ResizeChronoCableData "
310 "!\n");
311 }
312
313 // ShellElementsNodesH is the elements connectivity
314 // Important: in CableElementsNodesH[i][j] j index starts from 1 not zero
315 // This is because of how the GMF files are read in Chrono
316 CableElementsNodesH.resize(numCables);
317 for (size_t i = 0; i < numCables; i++) {
318 CableElementsNodesH[i].x = CableElementsNodesSTDVector[i][0];
319 CableElementsNodesH[i].y = CableElementsNodesSTDVector[i][1];
320 }
321 }
322 //------------------------------------------------------------------------------------
323
ResizeChronoShellsData(std::vector<std::vector<int>> ShellElementsNodesSTDVector,thrust::host_vector<int4> & ShellElementsNodesH)324 void ChFsiInterface::ResizeChronoShellsData(std::vector<std::vector<int>> ShellElementsNodesSTDVector,
325 thrust::host_vector<int4>& ShellElementsNodesH) {
326 auto my_mesh = chrono_types::make_shared<fea::ChMesh>();
327 if (mphysicalSystem.Get_otherphysicslist().size()) {
328 my_mesh = std::dynamic_pointer_cast<fea::ChMesh>(mphysicalSystem.Get_otherphysicslist().at(0));
329 }
330
331 int numShells = 0;
332 for (unsigned int i = 0; i < fsi_mesh->GetNelements(); i++) {
333 if (std::dynamic_pointer_cast<fea::ChElementShellANCF_3423>(fsi_mesh->GetElement(i)))
334 numShells++;
335 }
336
337 printf("numShells in ResizeChronoShellsData %d\n", numShells);
338 printf("ShellElementsNodesSTDVector.size() in ResizeChronoShellsData %zd\n", ShellElementsNodesSTDVector.size());
339
340 if (ShellElementsNodesSTDVector.size() != numShells) {
341 throw std::runtime_error(
342 "Size of the external data does not match the "
343 "ChSystem; thrown from ChFsiInterface::ResizeChronoShellsData "
344 "!\n");
345 }
346
347 // ShellElementsNodesH is the elements connectivity
348 // Important: in ShellElementsNodesH[i][j] j index starts from 1 not zero
349 // This is because of how the GMF files are read in Chrono
350 ShellElementsNodesH.resize(numShells);
351 for (size_t i = 0; i < numShells; i++) {
352 ShellElementsNodesH[i].x = ShellElementsNodesSTDVector[i][0];
353 ShellElementsNodesH[i].y = ShellElementsNodesSTDVector[i][1];
354 ShellElementsNodesH[i].z = ShellElementsNodesSTDVector[i][2];
355 ShellElementsNodesH[i].w = ShellElementsNodesSTDVector[i][3];
356 }
357
358 // (*ShellelementsNodes).resize(numShells);
359 // thrust::copy(ShellelementsNodes_H.begin(), ShellelementsNodes_H.end(), (*ShellElementsNodes).begin());
360 }
361
362 } // end namespace fsi
363 } // end namespace chrono
364