1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/aero/indvel.cc,v 1.19 2017/01/12 14:45:58 masarati Exp $ */
2 /*
3 * MBDyn (C) is a multibody analysis code.
4 * http://www.mbdyn.org
5 *
6 * Copyright (C) 1996-2017
7 *
8 * Pierangelo Masarati <masarati@aero.polimi.it>
9 * Paolo Mantegazza <mantegazza@aero.polimi.it>
10 *
11 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12 * via La Masa, 34 - 20156 Milano, Italy
13 * http://www.aero.polimi.it
14 *
15 * Changing this copyright notice is forbidden.
16 *
17 * This program is free software; you can redistribute it and/or modify
18 * it under the terms of the GNU General Public License as published by
19 * the Free Software Foundation (version 2 of the License).
20 *
21 *
22 * This program is distributed in the hope that it will be useful,
23 * but WITHOUT ANY WARRANTY; without even the implied warranty of
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 * GNU General Public License for more details.
26 *
27 * You should have received a copy of the GNU General Public License
28 * along with this program; if not, write to the Free Software
29 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30 */
31
32 /* Rotore */
33
34 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
35
36 #include <limits>
37 #include <cmath>
38
39 #ifdef USE_MPI
40 #include "mysleep.h"
41 const int mysleeptime = 300;
42
43 #ifdef MPI_PROFILING
44 extern "C" {
45 #include <mpe.h>
46 }
47 #include <cstdio>
48 #endif /* MPI_PROFILING */
49 #include "mbcomm.h"
50 #endif /* USE_MPI */
51
52 #include "indvel.h"
53 #include "dataman.h"
54
55 /* InducedVelocity - begin */
56
InducedVelocity(unsigned int uL,const StructNode * pCraft,ResForceSet ** ppres,flag fOut)57 InducedVelocity::InducedVelocity(unsigned int uL,
58 const StructNode *pCraft,
59 ResForceSet **ppres, flag fOut)
60 : Elem(uL, fOut),
61 #ifdef USE_MPI
62 is_parallel(false),
63 pBlockLenght(0),
64 pDispl(0),
65 ReqV(MPI::REQUEST_NULL),
66 pIndVelDataType(0),
67 #endif /* USE_MPI */
68 pCraft(pCraft),
69 ppRes(ppres)
70 {
71 #ifdef USE_MPI
72 iForcesVecDim = 6;
73 for (int i = 0; ppRes && ppRes[i]; i++) {
74 iForcesVecDim += 6;
75 }
76 SAFENEWARR(pTmpVecR, doublereal, iForcesVecDim);
77 SAFENEWARR(pTmpVecS, doublereal, iForcesVecDim);
78 if (MPI::Is_initialized()) {
79 is_parallel = true;
80 IndVelComm = MBDynComm.Dup();
81 }
82 #endif /* USE_MPI */
83
84 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
85 pthread_mutex_init(&induced_velocity_mutex, NULL);
86 pthread_cond_init(&induced_velocity_cond, NULL);
87 pthread_mutex_init(&forces_mutex, NULL);
88 #endif /* USE_MULTITHREAD && MBDYN_X_MT_ASSRES */
89 }
90
~InducedVelocity(void)91 InducedVelocity::~InducedVelocity(void)
92 {
93 #ifdef USE_MPI
94 SAFEDELETEARR(pTmpVecR);
95 SAFEDELETEARR(pTmpVecS);
96 #endif /* USE_MPI */
97
98 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
99 pthread_mutex_destroy(&induced_velocity_mutex);
100 pthread_cond_destroy(&induced_velocity_cond);
101 pthread_mutex_destroy(&forces_mutex);
102 #endif /* USE_MULTITHREAD && MBDYN_X_MT_ASSRES */
103 }
104
105 bool
bSectionalForces(void) const106 InducedVelocity::bSectionalForces(void) const
107 {
108 return false;
109 }
110
111 unsigned int
iGetNumPrivData(void) const112 InducedVelocity::iGetNumPrivData(void) const {
113 return 6;
114 }
115
116 unsigned int
iGetPrivDataIdx(const char * s) const117 InducedVelocity::iGetPrivDataIdx(const char *s) const
118 {
119 ASSERT(s != 0);
120
121 unsigned int idx = 0;
122
123 // sanity check
124 if (s[0] == '\0' || s[1] == '\0' || s[2] != '\0' ) {
125 return 0;
126 }
127
128 switch (s[0]) {
129 case 'M':
130 idx += 3;
131 // fallthru
132
133 case 'T':
134 switch (s[1]) {
135 case 'x':
136 return idx + 1;
137
138 case 'y':
139 return idx + 2;
140
141 case 'z':
142 return idx + 3;
143 }
144 }
145
146 return 0;
147 }
148
149 doublereal
dGetPrivData(unsigned int i) const150 InducedVelocity::dGetPrivData(unsigned int i) const
151 {
152 ASSERT(i > 0 && i <= 6);
153
154 switch (i) {
155 case 1:
156 case 2:
157 case 3:
158 return Res.Force()(i);
159
160 case 4:
161 case 5:
162 case 6:
163 return Res.Moment()(i - 3);
164 }
165
166 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
167 }
168
169 void
AfterConvergence(const VectorHandler &,const VectorHandler &)170 InducedVelocity::AfterConvergence(const VectorHandler& /* X */ ,
171 const VectorHandler& /* XP */ )
172 {
173 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
174 ASSERT(bDone);
175 bDone = false;
176 #endif // USE_MULTITHREAD && MBDYN_X_MT_ASSRES
177 }
178
179 /* assemblaggio jacobiano (nullo per tutti tranne che per il DynamicInflow) */
180 VariableSubMatrixHandler&
AssJac(VariableSubMatrixHandler & WorkMat,doublereal,const VectorHandler &,const VectorHandler &)181 InducedVelocity::AssJac(VariableSubMatrixHandler& WorkMat,
182 doublereal /* dCoef */ ,
183 const VectorHandler& /* XCurr */ ,
184 const VectorHandler& /* XPrimeCurr */ )
185 {
186 DEBUGCOUT("Entering InducedVelocity::AssJac()" << std::endl);
187 WorkMat.SetNullMatrix();
188
189 return WorkMat;
190 }
191
192 #ifdef USE_MPI
193 void
ExchangeLoads(flag fWhat)194 InducedVelocity::ExchangeLoads(flag fWhat)
195 {
196 #ifdef MPI_PROFILING
197 MPE_Log_event(33, 0, "start Induced Velocity Loads Exchange ");
198 #endif /* MPI_PROFILING */
199 /* Se il rotore � connesso ad una sola macchina non � necessario scambiare messaggi */
200 if (is_parallel && IndVelComm.Get_size() > 1) {
201 if (fWhat) {
202 /* Scambia F e M */
203 Res.Force().PutTo(&pTmpVecS[0]);
204 Res.Moment().PutTo(&pTmpVecS[3]);
205 for (int i = 0; ppRes && ppRes[i]; i++) {
206 ppRes[i]->pRes->Force().PutTo(&pTmpVecS[6 + 6*i]);
207 ppRes[i]->pRes->Moment().PutTo(&pTmpVecS[9 + 6*i]);
208 }
209 IndVelComm.Allreduce(pTmpVecS, pTmpVecR, iForcesVecDim, MPI::DOUBLE, MPI::SUM);
210 Res.PutForces(Vec3(&pTmpVecR[0]), Vec3(&pTmpVecR[3]));
211 for (int i = 0; ppRes && ppRes[i]; i++) {
212 ppRes[i]->pRes->PutForces(Vec3(&pTmpVecR[6 + 6*i]), Vec3(&pTmpVecR[9 + 6*i]));
213 }
214
215 } else {
216 IndVelComm.Allreduce(Res.Force().pGetVec(), pTmpVecR, 3, MPI::DOUBLE, MPI::SUM);
217 Res.PutForce(Vec3(pTmpVecR));
218 }
219 }
220 #ifdef MPI_PROFILING
221 MPE_Log_event(34, 0, "end Induced Velocity Loads Exchange ");
222 #endif /* MPI_PROFILING */
223 }
224
225 void
InitializeIndVelComm(MPI::Intracomm * iv)226 InducedVelocity::InitializeIndVelComm(MPI::Intracomm* iv)
227 {
228 ASSERT(is_parallel);
229 IndVelComm = *iv;
230 }
231
232 void
ExchangeVelocity(void)233 InducedVelocity::ExchangeVelocity(void)
234 {
235 #define ROTDATATYPELABEL 100
236 if (is_parallel && IndVelComm.Get_size() > 1) {
237 if (IndVelComm.Get_rank() == 0) {
238 for (int i = 1; i < IndVelComm.Get_size(); i++) {
239 IndVelComm.Send(MPI::BOTTOM, 1, *pIndVelDataType,
240 i, ROTDATATYPELABEL);
241 }
242 } else {
243 ReqV = IndVelComm.Irecv((void *)MPI::BOTTOM, 1,
244 *pIndVelDataType, 0, ROTDATATYPELABEL);
245 }
246 }
247 }
248 #endif // USE_MPI
249
250 // Somma alla trazione il contributo di forza di un elemento generico
251 void
AddForce(const Elem * pEl,const StructNode * pNode,const Vec3 & F,const Vec3 & M,const Vec3 & X)252 InducedVelocity::AddForce(const Elem *pEl, const StructNode *pNode,
253 const Vec3& F, const Vec3& M, const Vec3& X)
254 {
255 for (int i = 0; ppRes && ppRes[i]; i++) {
256 if (ppRes[i]->is_in(pEl->GetLabel())) {
257 ppRes[i]->pRes->AddForces(F, M, X);
258 }
259 }
260 }
261
262 // Somma alla trazione il contributo di forza di un elemento generico
263 // usando la forza e il momento per unita' di apertura e il peso
264 void
AddSectionalForce(Elem::Type type,const Elem * pEl,unsigned uPnt,const Vec3 & F,const Vec3 & M,doublereal dW,const Vec3 & X,const Mat3x3 & R,const Vec3 & V,const Vec3 & W)265 InducedVelocity::AddSectionalForce(Elem::Type type,
266 const Elem *pEl, unsigned uPnt,
267 const Vec3& F, const Vec3& M, doublereal dW,
268 const Vec3& X, const Mat3x3& R,
269 const Vec3& V, const Vec3& W)
270 {
271 ASSERT(bSectionalForces() == true);
272 InducedVelocity::AddForce(pEl, 0, F*dW, M*dW, X);
273 }
274
275 void
ResetForce(void)276 InducedVelocity::ResetForce(void)
277 {
278 Res.Reset(pCraft->GetXCurr());
279 for (int i = 0; ppRes && ppRes[i]; i++) {
280 ppRes[i]->pRes->Reset();
281 }
282 }
283
284 #if defined(USE_MULTITHREAD) && defined(MBDYN_X_MT_ASSRES)
285 void
Wait(void) const286 InducedVelocity::Wait(void) const
287 {
288 pthread_mutex_lock(&induced_velocity_mutex);
289 if (!bDone) {
290 pthread_cond_wait(&induced_velocity_cond,
291 &induced_velocity_mutex);
292 }
293 pthread_mutex_unlock(&induced_velocity_mutex);
294 }
295
296 void
Done(void) const297 InducedVelocity::Done(void) const
298 {
299 pthread_mutex_lock(&induced_velocity_mutex);
300 ASSERT(!bDone);
301 bDone = true;
302 pthread_cond_broadcast(&induced_velocity_cond);
303 pthread_mutex_unlock(&induced_velocity_mutex);
304 }
305 #endif // USE_MULTITHREAD && MBDYN_X_MT_ASSRES
306
307 /* InducedVelocity - end */
308
309 /* InducedVelocityElem - begin */
310
InducedVelocityElem(unsigned int uL,const DofOwner * pDO,const StructNode * pCraft,ResForceSet ** ppres,flag fOut)311 InducedVelocityElem::InducedVelocityElem(unsigned int uL, const DofOwner* pDO,
312 const StructNode *pCraft,
313 ResForceSet **ppres, flag fOut)
314 : Elem(uL, fOut),
315 AerodynamicElem(uL, pDO, fOut),
316 InducedVelocity(uL, pCraft, ppres, fOut)
317 {
318 NO_OP;
319 }
320
~InducedVelocityElem(void)321 InducedVelocityElem::~InducedVelocityElem(void)
322 {
323 NO_OP;
324 }
325
326 /* Tipo dell'elemento (usato per debug ecc.) */
327 Elem::Type
GetElemType(void) const328 InducedVelocityElem::GetElemType(void) const
329 {
330 return Elem::INDUCEDVELOCITY;
331 }
332
333 /* Tipo dell'elemento (usato per debug ecc.) */
334 AerodynamicElem::Type
GetAerodynamicElemType(void) const335 InducedVelocityElem::GetAerodynamicElemType(void) const
336 {
337 return AerodynamicElem::INDUCEDVELOCITY;
338 }
339
340 /* InducedVelocity - end */
341
342