1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/force.cc,v 1.64 2017/10/14 23:57:06 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 /* Forze */
33
34 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
35
36 #include <cfloat>
37
38 #include "dataman.h"
39 #include "force.h"
40 #include "strforce.h"
41 #include "strext.h"
42 #include "strmappingext.h"
43 #include "modalforce.h"
44 #include "modalext.h"
45 #include "modalmappingext.h"
46 #include "totalj.h"
47 #include "tpldrive_impl.h"
48
49 /* Force - begin */
50
51 std::ostream&
Restart(std::ostream & out) const52 Force::Restart(std::ostream& out) const
53 {
54 return out << " force: " << GetLabel();
55 }
56
57 /* Force - end */
58
59
60 /* AbstractForce - begin */
61
62 /* Costruttore non banale */
63
AbstractForce(unsigned int uL,const Node * pN,const DriveCaller * pDC,flag fOut)64 AbstractForce::AbstractForce(unsigned int uL, const Node* pN,
65 const DriveCaller* pDC, flag fOut)
66 : Elem(uL, fOut),
67 Force(uL, fOut),
68 DriveOwner(pDC),
69 pNode(pN)
70 {
71 NO_OP;
72 }
73
~AbstractForce(void)74 AbstractForce::~AbstractForce(void)
75 {
76 const Node2Scalar *pn2s = dynamic_cast<const Node2Scalar *>(pNode);
77 if (pn2s) {
78 SAFEDELETE(pn2s);
79 }
80 }
81
82 /* Contributo al file di restart */
83 std::ostream&
Restart(std::ostream & out) const84 AbstractForce::Restart(std::ostream& out) const
85 {
86 Force::Restart(out) << ", abstract, "
87 << pNode->GetLabel() << ", "
88 << psReadNodesNodes[pNode->GetNodeType()] << ", ";
89 return pGetDriveCaller()->Restart(out) << ';' << std::endl;
90 }
91
92
93 /* Assembla il residuo */
94 SubVectorHandler&
AssRes(SubVectorHandler & WorkVec,doublereal,const VectorHandler &,const VectorHandler &)95 AbstractForce::AssRes(SubVectorHandler& WorkVec,
96 doublereal /* dCoef */ ,
97 const VectorHandler& /* XCurr */ ,
98 const VectorHandler& /* XPrimeCurr */ )
99 {
100 DEBUGCOUT("Entering AbstractForce::AssRes()" << std::endl);
101
102 WorkVec.ResizeReset(1);
103
104 /* Dati */
105 doublereal dAmplitude = pGetDriveCaller()->dGet();
106
107 /* Indici delle incognite del nodo */
108 integer iFirstIndex = pNode->iGetFirstRowIndex();
109 WorkVec.PutRowIndex(1, iFirstIndex+1);
110
111 WorkVec.PutCoef(1, dAmplitude);
112
113 return WorkVec;
114 }
115
116 void
Output(OutputHandler & OH) const117 AbstractForce::Output(OutputHandler& OH) const
118 {
119 if (bToBeOutput()) {
120 if (OH.UseText(OutputHandler::FORCES)) {
121 OH.Forces()
122 << GetLabel()
123 << " " << pNode->GetLabel() << " " << dGet()
124 << std::endl;
125 }
126
127 /* TODO: NetCDF */
128 }
129 }
130
131 /* Contributo al residuo durante l'assemblaggio iniziale */
132 SubVectorHandler&
InitialAssRes(SubVectorHandler & WorkVec,const VectorHandler & XCurr)133 AbstractForce::InitialAssRes(SubVectorHandler& WorkVec,
134 const VectorHandler& XCurr)
135 {
136 DEBUGCOUT("Entering AbstractForce::InitialAssRes()" << std::endl);
137
138 return AssRes(WorkVec, 1., XCurr, XCurr);
139 }
140
141 /* AbstractForce - end */
142
143
144 /* AbstractInternalForce - begin */
145
146 /* Costruttore non banale */
147
AbstractInternalForce(unsigned int uL,const Node * pN1,const Node * pN2,const DriveCaller * pDC,flag fOut)148 AbstractInternalForce::AbstractInternalForce(unsigned int uL,
149 const Node* pN1, const Node *pN2,
150 const DriveCaller* pDC, flag fOut)
151 : Elem(uL, fOut),
152 Force(uL, fOut),
153 DriveOwner(pDC),
154 pNode1(pN1), pNode2(pN2)
155 {
156 NO_OP;
157 }
158
~AbstractInternalForce(void)159 AbstractInternalForce::~AbstractInternalForce(void)
160 {
161 const Node2Scalar *pn2s;
162
163 pn2s = dynamic_cast<const Node2Scalar *>(pNode1);
164 if (pn2s) {
165 SAFEDELETE(pn2s);
166 }
167
168 pn2s = dynamic_cast<const Node2Scalar *>(pNode2);
169 if (pn2s) {
170 SAFEDELETE(pn2s);
171 }
172 }
173
174 /* Contributo al file di restart */
175 std::ostream&
Restart(std::ostream & out) const176 AbstractInternalForce::Restart(std::ostream& out) const
177 {
178 Force::Restart(out) << ", abstract internal, "
179 << pNode1->GetLabel() << ", "
180 << psReadNodesNodes[pNode1->GetNodeType()] << ", "
181 << pNode1->GetLabel() << ", "
182 << psReadNodesNodes[pNode1->GetNodeType()] << ", ";
183 return pGetDriveCaller()->Restart(out) << ';' << std::endl;
184 }
185
186
187 /* Assembla il residuo */
188 SubVectorHandler&
AssRes(SubVectorHandler & WorkVec,doublereal,const VectorHandler &,const VectorHandler &)189 AbstractInternalForce::AssRes(SubVectorHandler& WorkVec,
190 doublereal /* dCoef */ ,
191 const VectorHandler& /* XCurr */ ,
192 const VectorHandler& /* XPrimeCurr */ )
193 {
194 DEBUGCOUT("Entering AbstractInternalForce::AssRes()" << std::endl);
195
196 WorkVec.ResizeReset(2);
197
198 /* Dati */
199 doublereal dAmplitude = pGetDriveCaller()->dGet();
200
201 /* Indici delle incognite del nodo */
202 integer iFirstIndex1 = pNode1->iGetFirstRowIndex();
203 integer iFirstIndex2 = pNode2->iGetFirstRowIndex();
204 WorkVec.PutRowIndex(1, iFirstIndex1 + 1);
205 WorkVec.PutRowIndex(2, iFirstIndex2 + 1);
206
207 WorkVec.PutCoef(1, dAmplitude);
208 WorkVec.PutCoef(2, -dAmplitude);
209
210 return WorkVec;
211 }
212
213 void
Output(OutputHandler & OH) const214 AbstractInternalForce::Output(OutputHandler& OH) const
215 {
216 if (bToBeOutput()) {
217 if (OH.UseText(OutputHandler::FORCES)) {
218 doublereal d = dGet();
219 OH.Forces()
220 << GetLabel()
221 << " " << pNode1->GetLabel() << " " << d
222 << " " << pNode2->GetLabel() << " " << -d
223 << std::endl;
224 }
225
226 /* TODO: NetCDF */
227 }
228 }
229
230 /* Contributo al residuo durante l'assemblaggio iniziale */
231 SubVectorHandler&
InitialAssRes(SubVectorHandler & WorkVec,const VectorHandler & XCurr)232 AbstractInternalForce::InitialAssRes(SubVectorHandler& WorkVec,
233 const VectorHandler& XCurr)
234 {
235 DEBUGCOUT("Entering AbstractInternalForce::InitialAssRes()" << std::endl);
236
237 return AssRes(WorkVec, 1., XCurr, XCurr);
238 }
239
240 /* AbstractInternalForce - end */
241
242
243 /* Legge una forza */
244
245 Elem *
ReadForce(DataManager * pDM,MBDynParser & HP,unsigned int uLabel,bool bCouple)246 ReadForce(DataManager* pDM,
247 MBDynParser& HP,
248 unsigned int uLabel,
249 bool bCouple)
250 {
251 const char* sKeyWords[] = {
252 "conservative", // deprecated
253 "absolute",
254 "absolute" "displacement",
255 "follower",
256
257 "conservative" "internal", // deprecated
258 "absolute" "internal",
259 "absolute" "internal" "displacement",
260 "follower" "internal",
261
262 "total", // not implented
263 "total" "internal",
264
265 "external" "structural",
266 "external" "structural" "mapping",
267
268 "modal",
269 "external" "modal",
270 "external" "modal" "mapping",
271
272 "abstract",
273 "abstract" "internal",
274
275 NULL
276 };
277
278 /* enum delle parole chiave */
279 enum KeyWords {
280 UNKNOWN = -1,
281
282 CONSERVATIVE, // deprecated
283 MB_ABSOLUTE, // NOTE: "ABSOLUTE" conflicts with gcc 4.6.1 on MinGW
284 ABSOLUTEDISPLACEMENT,
285 FOLLOWER,
286
287 CONSERVATIVEINTERNAL, // deprecated
288 ABSOLUTEINTERNAL,
289 ABSOLUTEINTERNALDISPLACEMENT,
290 FOLLOWERINTERNAL,
291
292 TOTAL,
293 TOTALINTERNAL,
294
295 EXTERNALSTRUCTURAL,
296 EXTERNALSTRUCTURALMAPPING,
297
298 MODALFORCE,
299 EXTERNALMODAL,
300 EXTERNALMODALMAPPING,
301
302 ABSTRACT,
303 ABSTRACTINTERNAL,
304
305 LASTKEYWORD
306 };
307
308 /* tabella delle parole chiave */
309 KeyTable K(HP, sKeyWords);
310
311 /* tipo di forza */
312 KeyWords CurrType = KeyWords(HP.GetWord());
313 if (CurrType == UNKNOWN) {
314 silent_cerr("Force(" << uLabel << "): unknown force type "
315 "at line " << HP.GetLineData() << std::endl);
316 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
317 }
318
319 switch (CurrType) {
320 case CONSERVATIVE:
321 case MB_ABSOLUTE:
322 case FOLLOWER:
323 case CONSERVATIVEINTERNAL:
324 case ABSOLUTEINTERNAL:
325 case FOLLOWERINTERNAL:
326 #if 0 /* not implemented yet */
327 case TOTAL:
328 #endif
329 case TOTALINTERNAL:
330 break;
331
332 default:
333 if (bCouple) {
334 silent_cerr("Force(" << uLabel << "): must be a \"force\"" << std::endl);
335 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
336 }
337 }
338
339 switch (CurrType) {
340 case CONSERVATIVE:
341 silent_cout("Force(" << uLabel << "): "
342 "deprecated \"conservative\" "
343 "at line " << HP.GetLineData() << "; "
344 "use \"absolute\" instead" << std::endl);
345 break;
346
347 case CONSERVATIVEINTERNAL:
348 silent_cout("Force(" << uLabel << "): "
349 "deprecated \"conservative internal\" "
350 "at line " << HP.GetLineData() << "; "
351 "use \"absolute internal\" instead" << std::endl);
352 break;
353
354 default:
355 break;
356 }
357
358 Elem* pEl = 0;
359
360 switch (CurrType) {
361 case ABSTRACT: {
362 /* tabella delle parole chiave */
363 KeyTable KDof(HP, psReadNodesNodes);
364 ScalarDof SD = ReadScalarDof(pDM, HP, true, false);
365 DriveCaller* pDC = HP.GetDriveCaller();
366 flag fOut = pDM->fReadOutput(HP, Elem::FORCE);
367
368 SAFENEWWITHCONSTRUCTOR(pEl,
369 AbstractForce,
370 AbstractForce(uLabel, SD.pNode, pDC, fOut));
371 } break;
372
373 case ABSTRACTINTERNAL: {
374 /* tabella delle parole chiave */
375 KeyTable KDof(HP, psReadNodesNodes);
376 ScalarDof SD1 = ReadScalarDof(pDM, HP, true, false);
377 ScalarDof SD2 = ReadScalarDof(pDM, HP, true, false);
378 DriveCaller* pDC = HP.GetDriveCaller();
379 flag fOut = pDM->fReadOutput(HP, Elem::FORCE);
380
381 SAFENEWWITHCONSTRUCTOR(pEl,
382 AbstractInternalForce,
383 AbstractInternalForce(uLabel, SD1.pNode, SD2.pNode, pDC, fOut));
384 } break;
385
386 case EXTERNALSTRUCTURAL:
387 pEl = ReadStructExtForce(pDM, HP, uLabel);
388 break;
389
390 case EXTERNALSTRUCTURALMAPPING:
391 pEl = ReadStructMappingExtForce(pDM, HP, uLabel);
392 break;
393
394 case MODALFORCE:
395 pEl = ReadModalForce(pDM, HP, uLabel);
396 break;
397
398 case EXTERNALMODAL:
399 pEl = ReadModalExtForce(pDM, HP, uLabel);
400 break;
401
402 case EXTERNALMODALMAPPING:
403 pEl = ReadModalMappingExtForce(pDM, HP, uLabel);
404 break;
405
406 case ABSOLUTEDISPLACEMENT:
407 pEl = ReadStructuralForce(pDM, HP, uLabel, true, bCouple, false, false);
408 break;
409
410 case CONSERVATIVE:
411 case MB_ABSOLUTE:
412 pEl = ReadStructuralForce(pDM, HP, uLabel, false, bCouple, false, false);
413 break;
414
415 case FOLLOWER:
416 pEl = ReadStructuralForce(pDM, HP, uLabel, false, bCouple, true, false);
417 break;
418
419 case ABSOLUTEINTERNALDISPLACEMENT:
420 pEl = ReadStructuralForce(pDM, HP, uLabel, true, bCouple, false, true);
421 break;
422
423 case CONSERVATIVEINTERNAL:
424 case ABSOLUTEINTERNAL:
425 pEl = ReadStructuralForce(pDM, HP, uLabel, false, bCouple, false, true);
426 break;
427
428 case FOLLOWERINTERNAL:
429 pEl = ReadStructuralForce(pDM, HP, uLabel, false, bCouple, true, true);
430 break;
431
432 case TOTAL:
433 case TOTALINTERNAL: {
434 /* nodo collegato 1 */
435 const StructNode* pNode1 = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
436
437 Vec3 f1(Zero3);
438 Mat3x3 R1h(Eye3);
439 Mat3x3 R1hr(Eye3);
440
441 ReferenceFrame RF(pNode1);
442
443 if (HP.IsKeyWord("position")) {
444 f1 = HP.GetPosRel(ReferenceFrame(pNode1));
445 }
446
447 if (HP.IsKeyWord("force" "orientation")) {
448 DEBUGCOUT("Force orientation matrix is supplied" << std::endl);
449 R1h = HP.GetRotRel(RF);
450 }
451
452 if (HP.IsKeyWord("moment" "orientation")) {
453 DEBUGCOUT("Moment orientation matrix is supplied" << std::endl);
454 R1hr = HP.GetRotRel(RF);
455 }
456
457 const StructNode* pNode2 = 0;
458 Vec3 f2(Zero3);
459 Mat3x3 R2h(Eye3);
460 Mat3x3 R2hr(Eye3);
461
462 if (CurrType == TOTALINTERNAL) {
463 /* nodo collegato 2 */
464 pNode2 = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
465
466 RF = ReferenceFrame(pNode2);
467
468 if (HP.IsKeyWord("position")) {
469 f2 = HP.GetPosRel(ReferenceFrame(pNode2));
470 }
471
472 if (HP.IsKeyWord("force" "orientation")) {
473 DEBUGCOUT("Force orientation matrix is supplied" << std::endl);
474 R2h = HP.GetRotRel(RF);
475 }
476
477 if (HP.IsKeyWord("moment" "orientation")) {
478 DEBUGCOUT("Moment orientation matrix is supplied" << std::endl);
479 R2hr = HP.GetRotRel(RF);
480 }
481 }
482
483 TplDriveCaller<Vec3>* pFDC = 0;
484 if (HP.IsKeyWord("force")) {
485 if (!HP.IsKeyWord("null")) {
486 pFDC = ReadDC3D(pDM, HP);
487 }
488 }
489
490 if (pFDC == 0) {
491 SAFENEW(pFDC, ZeroTplDriveCaller<Vec3>);
492 }
493
494 TplDriveCaller<Vec3>* pMDC = 0;
495 if (HP.IsKeyWord("moment")) {
496 if (!HP.IsKeyWord("null")) {
497 pMDC = ReadDC3D(pDM, HP);
498 }
499 }
500
501 if (pMDC == 0) {
502 SAFENEW(pMDC, ZeroTplDriveCaller<Vec3>);
503 }
504
505 flag fOut = pDM->fReadOutput(HP, Elem::FORCE);
506
507 switch (CurrType) {
508 case TOTALINTERNAL:
509 SAFENEWWITHCONSTRUCTOR(pEl,
510 TotalForce,
511 TotalForce(uLabel,
512 pFDC,
513 pMDC,
514 pNode1, f1, R1h, R1hr,
515 pNode2, f2, R2h, R2hr,
516 fOut));
517 break;
518
519 default:
520 silent_cerr("total force not implemented yet at line " << HP.GetLineData() << "; use a \"force\" and a \"couple\" element instead" << std::endl);
521 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
522 }
523 } break;
524
525 default:
526 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
527 }
528
529 /* Se non c'e' il punto e virgola finale */
530 if (HP.IsArg()) {
531 silent_cerr("Force(" << uLabel << "): "
532 "semicolon expected at line " << HP.GetLineData()
533 << std::endl);
534 throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
535 }
536
537 return pEl;
538 } /* End of ReadForce() */
539
540