MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
Macros | Functions | Variables
modal.cc File Reference
#include "mbconfig.h"
#include <cerrno>
#include <cstring>
#include <stdint.h>
#include <sys/stat.h>
#include <limits>
#include <algorithm>
#include "modal.h"
#include "dataman.h"
#include "Rot.hh"
#include "hint_impl.h"
Include dependency graph for modal.cc:

Go to the source code of this file.

Macros

#define MODAL_USE_GRAVITY
 

Functions

Joint * ReadModal (DataManager *pDM, MBDynParser &HP, const DofOwner *pDO, unsigned int uLabel)
 

Variables

static const char xyz [] = "xyz"
 
static const char * mdof []
 
static const char * rdof []
 
static const char * meq []
 
static const char * req []
 

Macro Definition Documentation

#define MODAL_USE_GRAVITY

Definition at line 109 of file modal.cc.

Function Documentation

Joint* ReadModal ( DataManager *  pDM,
MBDynParser &  HP,
const DofOwner *  pDO,
unsigned int  uLabel 
)

Definition at line 3189 of file modal.cc.

References Modal::a, AbsRefFrame, Mat3xN::AddMat3x3(), Mat3xN::AddVec(), ASSERT, c, Mat3xN::Copy(), Vec3::Cross(), DEBUGCOUT, Mat3xN::dGet(), MatNxN::dGet(), Modal::dMass, Eye3, grad::fabs(), DataManager::fReadOutput(), IncludeParser::GetFileName(), HighParser::GetInt(), WithLabel::GetLabel(), IncludeParser::GetLineData(), DataManager::GetLogFile(), MBDynParser::GetPosAbs(), MBDynParser::GetPosRel(), StructNode::GetRCurr(), HighParser::GetReal(), MBDynParser::GetRotAbs(), MBDynParser::GetRotRel(), HighParser::GetString(), HighParser::GetStringWithDelims(), StructNode::GetStructNodeType(), Mat3xN::GetVec(), Mat3x3::GetVec(), StructDispNode::GetXCurr(), HighParser::GetYesNoOrBool(), Modal::IdFEMNodes, Mat3xN::iGetNumCols(), VecN::iGetNumRows(), HighParser::IsKeyWord(), HighParser::IsStringWithDelims(), Elem::JOINT, Mat3xN::LeftMult(), MatCross, MatCrossCross, MBDYN_EXCEPT_ARGS, OutputHandler::MODAL, Modal::Modal(), StructNode::MODAL, Modal::NFEMNodes, Modal::NModes, Modal::NStrNodes, DataManager::OutputOpen(), DataManager::pFindNode(), StructDispNode::pGetRBK(), Modal::pInv10, Modal::pInv11, Modal::pInv3, Modal::pInv4, Modal::pInv5, Modal::pInv8, Modal::pInv9, Modal::pModalNode, Modal::pModeShapesr, Modal::pModeShapest, VecN::Put(), Mat3xN::Put(), MatNxN::Put(), Mat3xN::PutVec(), Modal::pXYZFEMNodes, Modal::R, SAFEDELETE, SAFENEWWITHCONSTRUCTOR, Modal::SND, grad::sqrt(), STRLENOF, Node::STRUCTURAL, Mat3xN::SubVec(), Modal::uModeNumber, Mat3x3::Write(), Zero3, and Zero3x3.

3193 {
3194  /* legge i dati d'ingresso e li passa al costruttore dell'elemento */
3195  Joint* pEl = 0;
3196 
3197  const ModalNode *pModalNode = 0;
3198  Vec3 X0(::Zero3);
3199  Mat3x3 R(::Eye3);
3200 
3201  /* If the modal element is clamped, a fixed position
3202  * and orientation of the reference point can be added */
3203  if (HP.IsKeyWord("clamped")) {
3204  if (HP.IsKeyWord("position")) {
3205  X0 = HP.GetPosAbs(::AbsRefFrame);
3206  }
3207 
3208  if (HP.IsKeyWord("orientation")) {
3209  R = HP.GetRotAbs(::AbsRefFrame);
3210  }
3211 
3212  /* otherwise a structural node of type "modal" must be given;
3213  * the "origin node" of the FEM model, if given, or the 0,0,0
3214  * coordinate point of the mesh will be put in the modal node */
3215  } else {
3216  unsigned int uNode = HP.GetInt();
3217 
3218  DEBUGCOUT("Linked to Modal Node: " << uNode << std::endl);
3219 
3220  /* verifica di esistenza del nodo */
3221  const StructNode* pTmpNode = pDM->pFindNode<const StructNode, Node::STRUCTURAL>(uNode);
3222 
3223  if (pTmpNode == 0) {
3224  silent_cerr("Modal(" << uLabel << "): "
3225  "StructuralNode(" << uNode << ") "
3226  "at line " << HP.GetLineData()
3227  << " not defined" << std::endl);
3228  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3229  }
3230 
3231  if (pTmpNode->GetStructNodeType() != StructNode::MODAL) {
3232  silent_cerr("Modal(" << uLabel << "): "
3233  "illegal type for "
3234  "StructuralNode(" << uNode << ") "
3235  "at line " << HP.GetLineData()
3236  << std::endl);
3237  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3238  }
3239  pModalNode = dynamic_cast<const ModalNode *>(pTmpNode);
3240  ASSERT(pModalNode != 0);
3241 
3242  if (pModalNode->pGetRBK() != 0) {
3243  silent_cerr("Modal(" << uLabel
3244  << "): StructuralNode(" << uNode
3245  << ") uses rigid body kinematics at line "
3246  << HP.GetLineData() << std::endl);
3247  throw ErrNotImplementedYet(MBDYN_EXCEPT_ARGS);
3248  }
3249  /* la posizione del nodo modale e' quella dell'origine del SdR
3250  * del corpo flessibile */
3251  X0 = pModalNode->GetXCurr();
3252 
3253  /* orientazione del corpo flessibile, data come orientazione
3254  * del nodo modale collegato */
3255  R = pModalNode->GetRCurr();
3256  }
3257 
3258  /* Legge i dati relativi al corpo flessibile */
3259  int tmpNModes = HP.GetInt(); /* numero di modi */
3260  if (tmpNModes <= 0) {
3261  silent_cerr("Modal(" << uLabel << "): "
3262  "illegal number of modes " << tmpNModes << " at line "
3263  << HP.GetLineData() << std::endl);
3264  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3265  }
3266  unsigned int NModes = (unsigned int)tmpNModes;
3267  unsigned int uLargestMode(0);
3268 
3269  std::vector<unsigned int> uModeNumber(NModes);
3270  if (HP.IsKeyWord("list")) {
3271  for (unsigned int iCnt = 0; iCnt < NModes; iCnt++) {
3272  int n = HP.GetInt();
3273 
3274  if (n <= 0) {
3275  silent_cerr("Modal(" << uLabel << "): "
3276  "illegal mode number "
3277  "for mode #" << iCnt + 1
3278  << " at line " << HP.GetLineData()
3279  << std::endl);
3280  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3281  }
3282 
3283  std::vector<unsigned int>::const_iterator
3284  iv = std::find(uModeNumber.begin(),
3285  uModeNumber.end(), (unsigned int)n);
3286  if (iv != uModeNumber.end()) {
3287  silent_cerr("Modal(" << uLabel << "): "
3288  "mode #" << iCnt + 1
3289  << " already defined "
3290  "at line " << HP.GetLineData()
3291  << std::endl);
3292  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3293  }
3294 
3295  if ((unsigned int)n > uLargestMode) {
3296  uLargestMode = (unsigned int)n;
3297  }
3298 
3299  uModeNumber[iCnt] = n;
3300  }
3301 
3302  } else {
3303  for (unsigned int iCnt = 1; iCnt <= NModes; iCnt++) {
3304  uModeNumber[iCnt - 1] = iCnt;
3305  }
3306  uLargestMode = NModes;
3307  }
3308 
3309  std::vector<doublereal> InitialValues(NModes, 0.);
3310  std::vector<doublereal> InitialValuesP(NModes, 0.);
3311  bool bInitialValues(false);
3312  if (HP.IsKeyWord("initial" "value")) {
3313  bInitialValues = true;
3314 
3315  if (HP.IsKeyWord("mode")) {
3316  int iCnt = 0;
3317  do {
3318  int n = HP.GetInt();
3319 
3320  if (n <= 0) {
3321  silent_cerr("Modal(" << uLabel << "): "
3322  "illegal mode number "
3323  "for initial value #" << iCnt + 1
3324  << " at line " << HP.GetLineData()
3325  << std::endl);
3326  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3327  }
3328 
3329  std::vector<unsigned int>::const_iterator
3330  iv = std::find(uModeNumber.begin(),
3331  uModeNumber.end(), (unsigned int)n);
3332  if (iv == uModeNumber.end()) {
3333  silent_cerr("Modal(" << uLabel << "): "
3334  "mode " << n << " not defined "
3335  "at line " << HP.GetLineData()
3336  << std::endl);
3337  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3338  }
3339 
3340  unsigned iIdx = iv - uModeNumber.begin();
3341  InitialValues[iIdx] = HP.GetReal();
3342  InitialValuesP[iIdx] = HP.GetReal();
3343 
3344  iCnt++;
3345  } while (HP.IsKeyWord("mode"));
3346 
3347  } else {
3348  for (unsigned int iCnt = 0; iCnt < NModes; iCnt++) {
3349  InitialValues[0] = HP.GetReal();
3350  InitialValuesP[0] = HP.GetReal();
3351  }
3352  }
3353  }
3354 
3355  /* numero di nodi FEM del modello */
3356  unsigned int NFEMNodes = unsigned(-1);
3357  if (!HP.IsKeyWord("from" "file")) {
3358  int tmpNFEMNodes = HP.GetInt();
3359  if (tmpNFEMNodes <= 0) {
3360  silent_cerr("Modal(" << uLabel << "): "
3361  "illegal number of FEM nodes " << tmpNFEMNodes
3362  << " at line " << HP.GetLineData() << std::endl);
3363  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3364  }
3365  NFEMNodes = unsigned(tmpNFEMNodes);
3366  }
3367 
3368 #ifdef MODAL_SCALE_DATA
3369  /* Legge gli eventuali fattori di scala per le masse nodali
3370  * (scale masses) e per i modi (scale modes) */
3371 
3372  /* NOTA: AL MOMENTO NON E' USATO */
3373  doublereal scalemasses(1.);
3374  doublereal scaleinertia(1.);
3375  doublereal scalemodes(1.);
3376 
3377  if (HP.IsKeyWord("scale" "masses")) {
3378  scalemasses = HP.GetReal();
3379  }
3380 
3381  if (HP.IsKeyWord("scale" "modes")) {
3382  scalemodes = HP.GetReal();
3383  }
3384 
3385  scaleinertia = scalemasses*(scalemodes*scalemodes);
3386 #endif /* MODAL_SCALE_DATA */
3387 
3388  /* Legge i coefficienti di smorzamento */
3389  enum {
3390  DAMPING_FROM_FILE = 0,
3391  DAMPING_NO,
3392  // DAMPING_PROPORTIONAL, // obsoleted
3393  DAMPING_RAYLEIGH,
3394  DAMPING_SINGLE_FACTOR,
3395  DAMPING_DIAG
3396  } eDamp(DAMPING_FROM_FILE);
3397  const char *sDamp[] = {
3398  "damping from file",
3399  "no damping",
3400  // "proportional damping", // obsoleted
3401  "Rayleigh damping",
3402  "single factor damping",
3403  "diag damping",
3404  NULL
3405  };
3406 
3407  doublereal damp_factor(0.), damp_coef_M(0.), damp_coef_K(0.);
3408  std::vector<doublereal> DampRatios;
3409 
3410  if (HP.IsKeyWord("no" "damping")) {
3411  eDamp = DAMPING_NO;
3412 
3413  } else if (HP.IsKeyWord("rayleigh" "damping")) {
3414  eDamp = DAMPING_RAYLEIGH;
3415  damp_coef_M = HP.GetReal();
3416  if (damp_coef_M < 0.) {
3417  silent_cerr("Modal(" << uLabel << "): "
3418  "warning, negative mass damping coefficient for \"Rayleigh damping\" "
3419  "at line " << HP.GetLineData() << std::endl);
3420  }
3421  damp_coef_K = HP.GetReal();
3422  if (damp_coef_K < 0.) {
3423  silent_cerr("Modal(" << uLabel << "): "
3424  "warning, negative stiffness damping coefficient for \"Rayleigh damping\" "
3425  "at line " << HP.GetLineData() << std::endl);
3426  }
3427 
3428  } else if (HP.IsKeyWord("single" "factor" "damping")) {
3429  eDamp = DAMPING_SINGLE_FACTOR;
3430  damp_factor = HP.GetReal();
3431  if (damp_factor < 0.) {
3432  silent_cerr("Modal(" << uLabel << "): "
3433  "warning, negative damping factor for \"single factor damping\" "
3434  "at line " << HP.GetLineData() << std::endl);
3435  }
3436 
3437  } else if (HP.IsKeyWord("proportional" "damping")) {
3438  silent_cerr("Modal(" << uLabel << "): "
3439  "warning, \"proportional damping\" is deprecated; "
3440  "use \"single factor damping\" with the desired damping factor instead "
3441  "at line " << HP.GetLineData() << std::endl);
3442  // NOTE: "proportional damping" deprecated, replaced by "single factor damping"
3443  eDamp = DAMPING_SINGLE_FACTOR;
3444  damp_factor = HP.GetReal();
3445  if (damp_factor < 0.) {
3446  silent_cerr("Modal(" << uLabel << "): "
3447  "warning, negative damping factor for \"proportional damping\" "
3448  "at line " << HP.GetLineData() << std::endl);
3449  }
3450 
3451  } else if (HP.IsKeyWord("diag" "damping")) {
3452  eDamp = DAMPING_DIAG;
3453  DampRatios.resize(NModes);
3454  fill(DampRatios.begin(), DampRatios.end(), 0.);
3455 
3456  if (HP.IsKeyWord("all")) {
3457  for (unsigned int iCnt = 0; iCnt < NModes; iCnt ++) {
3458  damp_factor = HP.GetReal();
3459  if (damp_factor < 0.) {
3460  silent_cerr("Modal(" << uLabel << "): "
3461  "warning, negative damping factor for \"diag damping\" "
3462  "of mode " << (iCnt + 1 ) << " "
3463  "at line " << HP.GetLineData() << std::endl);
3464  }
3465  DampRatios[iCnt] = damp_factor;
3466  }
3467 
3468  } else {
3469  int iDM = HP.GetInt();
3470  if (iDM <= 0 || unsigned(iDM) > NModes) {
3471  silent_cerr("invalid number of damped modes "
3472  "at line " << HP.GetLineData() <<std::endl);
3473  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3474  }
3475 
3476  std::vector<bool> gotit(NModes, false);
3477 
3478  for (unsigned int iCnt = 1; iCnt <= unsigned(iDM); iCnt ++) {
3479  integer iDampedMode = HP.GetInt();
3480  if (iDampedMode <= 0 || unsigned(iDampedMode) > NModes) {
3481  silent_cerr("invalid index " << iDampedMode
3482  << " for damped mode #" << iCnt
3483  << " of " << iDM
3484  << " at line " << HP.GetLineData() <<std::endl);
3485  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3486  }
3487 
3488  if (gotit[iDampedMode - 1]) {
3489  silent_cerr("damping already provided for mode " << iDampedMode
3490  << " (damped mode #" << iCnt
3491  << " of " << iDM
3492  << ") at line " << HP.GetLineData() <<std::endl);
3493  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3494 
3495  }
3496  gotit[iDampedMode - 1] = true;
3497  damp_factor = HP.GetReal();
3498  if (damp_factor < 0.) {
3499  silent_cerr("Modal(" << uLabel << "): "
3500  "warning, negative damping factor for \"diag damping\" "
3501  "of mode " << iDampedMode << " "
3502  "at line " << HP.GetLineData() << std::endl);
3503  }
3504  DampRatios[iDampedMode - 1] = damp_factor;
3505  }
3506  }
3507  } else if (HP.IsKeyWord("damping" "from" "file")) {
3508  eDamp = DAMPING_FROM_FILE;
3509  } else {
3510  silent_cout("Modal(" << uLabel << "): "
3511  "no damping is assumed at line "
3512  << HP.GetLineData() << " (deprecated)"
3513  << std::endl);
3514  }
3515 
3516  DEBUGCOUT("Number of Modes Imported : " << NModes << std::endl);
3517  DEBUGCOUT("Number of FEM Nodes Imported : " << NFEMNodes << std::endl);
3518  DEBUGCOUT("Origin of FEM Model : " << X0 << std::endl);
3519  DEBUGCOUT("Orientation of FEM Model : " << R << std::endl);
3520  /* DEBUGCOUT("Damping coefficient: "<< cdamp << std::endl); */
3521 
3522  doublereal dMass = 0; // total mass
3523  Vec3 XTmpIn(::Zero3); // center of mass location
3524  Vec3 STmp(::Zero3); // static moment
3525  Mat3x3 JTmp(::Zero3x3); // total inertia from input
3526  Mat3x3 JTmpIn(::Zero3x3); // total inertia after manipulation
3527  std::vector<doublereal> FEMMass; // nodal masses
3528  std::vector<Vec3> FEMJ; // nodal inertia (diagonal)
3529 
3530  Mat3xN *pModeShapest = 0; // displacement and rotation shapes
3531  Mat3xN *pModeShapesr = 0;
3532  Mat3xN PHIti(NModes, 0.); // i-th node shapes: 3*nmodes
3533  Mat3xN PHIri(NModes, 0.);
3534  Mat3xN *pXYZFEMNodes = 0; // pointer to nodal coords
3535  MatNxN *pGenMass = 0; // pointer to modal mass
3536  MatNxN *pGenStiff = 0; // pointer to modal stiffness
3537  MatNxN *pGenDamp = 0; // pointer to modal damping
3538 
3539  // inertia invariants
3540  Mat3xN *pInv3 = 0;
3541  Mat3xN *pInv4 = 0;
3542  Mat3xN *pInv5 = 0;
3543  Mat3xN *pInv8 = 0;
3544  Mat3xN *pInv9 = 0;
3545  Mat3xN *pInv10 = 0;
3546  Mat3xN *pInv11 = 0;
3547 
3548  // modal displacements and velocities
3549  VecN *a = 0;
3550  VecN *aP = 0;
3551 
3552  // FEM database file name
3553  const char *s = HP.GetFileName();
3554  if (s == 0) {
3555  silent_cerr("Modal(" << uLabel << "): unable to get "
3556  "modal data file name at line " << HP.GetLineData()
3557  << std::endl);
3558  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3559  }
3560 
3561  std::string sFileFEM(s);
3562 
3563  doublereal dMassThreshold = 0.;
3564  doublereal dDampingThreshold = 0.;
3565  doublereal dStiffnessThreshold = 0.;
3566  while (true) {
3567  if (HP.IsKeyWord("threshold")) {
3568  dStiffnessThreshold = HP.GetReal();
3569  if (dStiffnessThreshold < 0.) {
3570  silent_cerr("Modal(" << uLabel << "): "
3571  "invalid threshold " << dStiffnessThreshold
3572  << " at line " << HP.GetLineData()
3573  << std::endl);
3574  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3575  }
3576  dMassThreshold = dDampingThreshold = dStiffnessThreshold;
3577 
3578  } else if (HP.IsKeyWord("mass" "threshold")) {
3579  dMassThreshold = HP.GetReal();
3580  if (dMassThreshold < 0.) {
3581  silent_cerr("Modal(" << uLabel << "): "
3582  "invalid mass threshold " << dMassThreshold
3583  << " at line " << HP.GetLineData()
3584  << std::endl);
3585  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3586  }
3587 
3588  } else if (HP.IsKeyWord("damping" "threshold")) {
3589  dDampingThreshold = HP.GetReal();
3590  if (dDampingThreshold < 0.) {
3591  silent_cerr("Modal(" << uLabel << "): "
3592  "invalid damping threshold " << dDampingThreshold
3593  << " at line " << HP.GetLineData()
3594  << std::endl);
3595  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3596  }
3597 
3598  } else if (HP.IsKeyWord("stiffness" "threshold")) {
3599  dStiffnessThreshold = HP.GetReal();
3600  if (dStiffnessThreshold < 0.) {
3601  silent_cerr("Modal(" << uLabel << "): "
3602  "invalid stiffness threshold " << dStiffnessThreshold
3603  << " at line " << HP.GetLineData()
3604  << std::endl);
3605  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3606  }
3607 
3608  } else {
3609  break;
3610  }
3611  }
3612 
3613  /* loop for binary keywords */
3614  bool bCreateBinary = false;
3615  bool bUseBinary = false;
3616  bool bUpdateBinary = false;
3617  while (true) {
3618  if (HP.IsKeyWord("create" "binary")) {
3619  bCreateBinary = true;
3620 
3621  } else if (HP.IsKeyWord("use" "binary")) {
3622  bUseBinary = true;
3623 
3624  } else if (HP.IsKeyWord("update" "binary")) {
3625  bUpdateBinary = true;
3626 
3627  } else {
3628  break;
3629  }
3630  }
3631 
3632  std::string sEchoFileName;
3633  int iEchoPrecision(13);
3634  if (HP.IsKeyWord("echo")) {
3635  const char *s = HP.GetFileName();
3636  if (s == 0) {
3637  silent_cerr("Modal(" << uLabel << "): "
3638  "unable to parse echo file name at line " << HP.GetLineData()
3639  << std::endl);
3640  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3641  }
3642  sEchoFileName = s;
3643 
3644  if (HP.IsKeyWord("precision")) {
3645  iEchoPrecision = HP.GetInt();
3646  if (iEchoPrecision <= 0) {
3647  silent_cerr("Modal(" << uLabel << "): "
3648  "invalid precision at line " << HP.GetLineData()
3649  << std::endl);
3650  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3651  }
3652  }
3653  }
3654 
3655  /* stuff for binary FEM file handling */
3656  std::string sBinFileFEM;
3657  struct stat stFEM, stBIN;
3658  bool bReadFEM = true,
3659  bWriteBIN = false,
3660  bCheckBIN = false;
3661 
3662  if (stat(sFileFEM.c_str(), &stFEM) == -1) {
3663  int save_errno = errno;
3664  char *errmsg = strerror(save_errno);
3665 
3666  silent_cerr("Modal(" << uLabel << "): "
3667  "unable to stat(\"" << sFileFEM << "\") "
3668  "(" << save_errno << ": " << errmsg << ")"
3669  << std::endl);
3670  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3671  }
3672 
3673  if (bUseBinary || bCreateBinary || bUpdateBinary) {
3674  sBinFileFEM = sFileFEM + ".bin";
3675 
3676  if (stat(sBinFileFEM.c_str(), &stBIN) == -1) {
3677  int save_errno = errno;
3678  char *errmsg = strerror(save_errno);
3679 
3680  if (save_errno == ENOENT && (bCreateBinary || bUpdateBinary)) {
3681  silent_cerr("Modal(" << uLabel << "): "
3682  "creating binary file "
3683  "\"" << sBinFileFEM << "\" "
3684  "from file "
3685  "\"" << sFileFEM << "\""
3686  << std::endl);
3687 
3688  } else {
3689  silent_cerr("Modal(" << uLabel << "): "
3690  "unable to stat(\"" << sBinFileFEM << "\") "
3691  "(" << save_errno << ": " << errmsg << ")"
3692  << std::endl);
3693  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3694  }
3695 
3696  } else {
3697  /* can check bin */
3698  bCheckBIN = true;
3699  }
3700  }
3701 
3702  if (bUseBinary && bCheckBIN) {
3703  /* if timestamp of binary is later than of ASCII use binary */
3704  if (stBIN.st_mtime > stFEM.st_mtime) {
3705  bReadFEM = false;
3706 
3707  /* otherwise, if requested, update binary */
3708  } else {
3709  if (bUpdateBinary) {
3710  bWriteBIN = true;
3711 
3712  silent_cout("Modal(" << uLabel << "): "
3713  "binary database file \"" << sBinFileFEM << "\" "
3714  "older than text database file \"" << sFileFEM << "\"; "
3715  "updating"
3716  << std::endl);
3717 
3718  } else {
3719  silent_cerr("Modal(" << uLabel << "): "
3720  "warning, binary database file \"" << sBinFileFEM << "\" "
3721  "older than text database file \"" << sFileFEM << "\"; "
3722  "using text database file "
3723  "(enable \"update binary\" to refresh binary database file)"
3724  << std::endl);
3725  }
3726  }
3727 
3728  } else if (bCreateBinary || bUpdateBinary) {
3729  bWriteBIN = true;
3730  }
3731 
3732  /* alloca la memoria per le matrici necessarie a memorizzare i dati
3733  * relativi al corpo flessibile
3734  * nota: devo usare i puntatori perche' altrimenti non si riesce
3735  * a passarli al costruttore
3736  */
3737  SAFENEWWITHCONSTRUCTOR(pGenMass, MatNxN, MatNxN(NModes, 0.));
3738  SAFENEWWITHCONSTRUCTOR(pGenStiff, MatNxN, MatNxN(NModes, 0.));
3739  SAFENEWWITHCONSTRUCTOR(pGenDamp, MatNxN, MatNxN(NModes, 0.));
3740 
3741  SAFENEWWITHCONSTRUCTOR(a, VecN, VecN(NModes, 0.));
3742  SAFENEWWITHCONSTRUCTOR(aP, VecN, VecN(NModes, 0.));
3743 
3744  /* array contenente le label dei nodi FEM */
3745  std::vector<std::string> IdFEMNodes;
3746  std::vector<bool> bActiveModes;
3747 
3748  enum {
3749  MODAL_VERSION_UNKNOWN = 0,
3750 
3751  MODAL_VERSION_1 = 1,
3752  MODAL_VERSION_2 = 2,
3753  MODAL_VERSION_3 = 3, // after adding record 13 (generalized damping)
3754  MODAL_VERSION_4 = 4, // after making sizes portable
3755 
3756  MODAL_VERSION_LAST
3757  };
3758 
3759  /* NOTE: increment this each time the binary format changes! */
3760  const char magic[5] = "bmod";
3761  const uint32_t BinVersion = MODAL_VERSION_4;
3762 
3763  uint32_t currBinVersion = MODAL_VERSION_UNKNOWN;
3764  char checkPoint;
3765  uint32_t NFEMNodesFEM = 0, NModesFEM = 0;
3766 
3767  bool bBuildInvariants = false;
3768 
3769  enum {
3770  // note: we use 1-based array
3771  MODAL_RECORD_1 = 1,
3772  MODAL_RECORD_2 = 2,
3773  MODAL_RECORD_3 = 3,
3774  MODAL_RECORD_4 = 4,
3775  MODAL_RECORD_5 = 5,
3776  MODAL_RECORD_6 = 6,
3777  MODAL_RECORD_7 = 7,
3778  MODAL_RECORD_8 = 8,
3779  MODAL_RECORD_9 = 9,
3780  MODAL_RECORD_10 = 10,
3781  MODAL_RECORD_11 = 11,
3782  MODAL_RECORD_12 = 12,
3783  MODAL_RECORD_13 = 13,
3784 
3785  // NOTE: do not exceed 127
3786 
3787  MODAL_LAST_RECORD,
3788 
3789  MODAL_END_OF_FILE = char(-1)
3790  };
3791 
3792  /* Note: to keep it readable, we use a base 1 array */
3793  bool bRecordGroup[MODAL_LAST_RECORD] = { false };
3794 
3795  // record 1: header
3796  // record 2: FEM node labels
3797  // record 3: initial modal amplitudes
3798  // record 4: initial modal amplitude rates
3799  // record 5: FEM node X
3800  // record 6: FEM node Y
3801  // record 7: FEM node Z
3802  // record 8: shapes
3803  // record 9: generalized mass matrix
3804  // record 10: generalized stiffness matrix
3805  // record 11: (optional) diagonal mass matrix
3806  // record 12: (optional) rigid body inertia
3807  // record 13: (optional) damping matrix
3808 
3809  std::string fname;
3810  if (bReadFEM) {
3811  /* apre il file con i dati del modello FEM */
3812  std::ifstream fdat(sFileFEM.c_str());
3813  if (!fdat) {
3814  silent_cerr("Modal(" << uLabel << "): "
3815  "unable to open file \"" << sFileFEM << "\""
3816  "at line " << HP.GetLineData() << std::endl);
3817  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3818  }
3819 
3820  /* don't leave behind a corrupted .bin file */
3821  try {
3822 
3823  std::ofstream fbin;
3824  if (bWriteBIN) {
3825  fbin.open(sBinFileFEM.c_str(), std::ios::binary | std::ios::trunc);
3826  if (!fbin) {
3827  silent_cerr("Modal(" << uLabel << "): "
3828  "unable to open file \"" << sBinFileFEM << "\""
3829  "at line " << HP.GetLineData() << std::endl);
3830  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3831  }
3832 
3833  fbin.write(&magic[0], sizeof(4));
3834  currBinVersion = BinVersion;
3835  fbin.write((const char *)&currBinVersion, sizeof(currBinVersion));
3836  }
3837 
3838  /* carica i dati relativi a coordinate nodali, massa, momenti statici
3839  * e d'inerzia, massa e rigidezza generalizzate dal file nomefile.
3840  */
3841 
3842  doublereal d;
3843  unsigned int NRejModes = 0;
3844  char str[BUFSIZ];
3845 
3846  while (fdat && !fdat.eof()) { /* parsing del file */
3847  fdat.getline(str, sizeof(str));
3848 
3849  if (strncmp("** END OF FILE", str, STRLENOF("** END OF FILE")) == 0) {
3850  break;
3851  }
3852 
3853  /* legge il primo blocco (HEADER) */
3854  if (strncmp("** RECORD GROUP ", str,
3855  STRLENOF("** RECORD GROUP ")) != 0) {
3856  continue;
3857  }
3858 
3859  char *next;
3860  long rg = std::strtol(&str[STRLENOF("** RECORD GROUP ")], &next, 10);
3861  if (next == &str[STRLENOF("** RECORD GROUP ")]) {
3862  silent_cerr("file=\"" << sFileFEM << "\": "
3863  "unable to parse \"RECORD GROUP\" number <" << &str[STRLENOF("** RECORD GROUP ")] << ">"
3864  << std::endl);
3865  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3866  }
3867 
3868  if (!bRecordGroup[MODAL_RECORD_1] && rg != MODAL_RECORD_1) {
3869  silent_cerr("file=\"" << sFileFEM << "\": "
3870  "\"RECORD GROUP 1\" not parsed yet; cannot parse \"RECORD GROUP " << rg << "\""
3871  << std::endl);
3872  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3873  }
3874 
3875  switch (rg) {
3876  case MODAL_RECORD_1: {
3877  if (bRecordGroup[MODAL_RECORD_1]) {
3878  silent_cerr("file=\"" << sFileFEM << "\": "
3879  "\"RECORD GROUP 1\" already parsed"
3880  << std::endl);
3881  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3882  }
3883 
3884  fdat.getline(str, sizeof(str));
3885 
3886  /* ignore REV by now */
3887  fdat >> str;
3888 
3889  /* FEM nodes number */
3890  fdat >> NFEMNodesFEM;
3891 
3892  /* add to modes number */
3893  fdat >> NModesFEM;
3894 
3895  unsigned int i;
3896  fdat >> i;
3897  NModesFEM += i;
3898  fdat >> i;
3899  NModesFEM += i;
3900 
3901  /* "rejected" modes (subtract from modes number) */
3902  fdat >> NRejModes;
3903  NModesFEM -= NRejModes;
3904 
3905  /* consistency checks */
3906  if (NFEMNodes == unsigned(-1)) {
3907  NFEMNodes = NFEMNodesFEM;
3908 
3909  } else if (NFEMNodes != NFEMNodesFEM) {
3910  silent_cerr("Modal(" << uLabel << "), "
3911  "file \"" << sFileFEM << "\": "
3912  "FEM nodes number " << NFEMNodes
3913  << " does not match node number "
3914  << NFEMNodesFEM
3915  << std::endl);
3916  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3917  }
3918 
3919  if (NModes > NModesFEM) {
3920  silent_cerr("Modal(" << uLabel << "), "
3921  "file '" << sFileFEM << "': "
3922  "number of requested modes "
3923  "(" << NModes << ") "
3924  "exceeds number of available "
3925  "modes (" << NModesFEM << ")"
3926  << std::endl);
3927  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3928 
3929  }
3930 
3931  if (NModes < NModesFEM) {
3932  silent_cout("Modal(" << uLabel << "), "
3933  "file '" << sFileFEM << "': "
3934  "using " << NModes
3935  << " of " << NModesFEM
3936  << " available modes"
3937  << std::endl);
3938  }
3939 
3940  if (uLargestMode > NModesFEM) {
3941  silent_cerr("Modal(" << uLabel << "), "
3942  "file '" << sFileFEM << "': "
3943  "largest requested mode "
3944  "(" << uLargestMode << ") "
3945  "exceeds available modes "
3946  "(" << NModesFEM << ")"
3947  << std::endl);
3948  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
3949  }
3950 
3951  if (!bActiveModes.empty()) {
3952  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3953  }
3954 
3955  IdFEMNodes.resize(NFEMNodes);
3956 
3957  SAFENEWWITHCONSTRUCTOR(pXYZFEMNodes, Mat3xN, Mat3xN(NFEMNodes, 0.));
3958  SAFENEWWITHCONSTRUCTOR(pModeShapest, Mat3xN, Mat3xN(NFEMNodes*NModes, 0.));
3959  SAFENEWWITHCONSTRUCTOR(pModeShapesr, Mat3xN, Mat3xN(NFEMNodes*NModes, 0.));
3960 
3961  bActiveModes.resize(NModesFEM + 1, false);
3962 
3963  for (unsigned int iCnt = 0; iCnt < NModes; iCnt++) {
3964  if (uModeNumber[iCnt] > NModesFEM) {
3965  silent_cerr("Modal(" << uLabel << "): "
3966  "mode " << uModeNumber[iCnt]
3967  << " is not available (max = "
3968  << NModesFEM << ")"
3969  << std::endl);
3970  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3971  }
3972  bActiveModes[uModeNumber[iCnt]] = true;
3973  }
3974 
3975  if (bWriteBIN) {
3976  checkPoint = MODAL_RECORD_1;
3977  fbin.write(&checkPoint, sizeof(checkPoint));
3978  fbin.write((const char *)&NFEMNodesFEM, sizeof(NFEMNodesFEM));
3979  fbin.write((const char *)&NModesFEM, sizeof(NModesFEM));
3980  }
3981 
3982  bRecordGroup[MODAL_RECORD_1] = true;
3983  } break;
3984 
3985  case MODAL_RECORD_2: {
3986  /* legge il secondo blocco (Id.nodi) */
3987  if (bRecordGroup[MODAL_RECORD_2]) {
3988  silent_cerr("file=\"" << sFileFEM << "\": "
3989  "\"RECORD GROUP 2\" already parsed"
3990  << std::endl);
3991  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
3992  }
3993 
3994  for (unsigned int iNode = 0; iNode < NFEMNodes; iNode++) {
3995  // NOTE: this precludes the possibility
3996  // to have blanks in node labels
3997  fdat >> IdFEMNodes[iNode];
3998  }
3999 
4000  if (bWriteBIN) {
4001  checkPoint = MODAL_RECORD_2;
4002  fbin.write(&checkPoint, sizeof(checkPoint));
4003  for (unsigned int iNode = 0; iNode < NFEMNodes; iNode++) {
4004  uint32_t len = IdFEMNodes[iNode].size();
4005  fbin.write((const char *)&len, sizeof(len));
4006  fbin.write((const char *)IdFEMNodes[iNode].c_str(), len);
4007  }
4008  }
4009 
4010  bRecordGroup[2] = true;
4011  } break;
4012 
4013  case MODAL_RECORD_3: {
4014  /* deformate iniziali dei modi */
4015  if (bRecordGroup[MODAL_RECORD_3]) {
4016  silent_cerr("file=\"" << sFileFEM << "\": "
4017  "\"RECORD GROUP 3\" already parsed"
4018  << std::endl);
4019  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4020  }
4021 
4022  unsigned int iCnt = 1;
4023 
4024  if (bActiveModes.empty()) {
4025  silent_cerr("Modal(" << uLabel << "): "
4026  "input file \"" << sFileFEM << "\" "
4027  "is bogus (RECORD GROUP 3)"
4028  << std::endl);
4029  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4030  }
4031 
4032  if (bWriteBIN) {
4033  checkPoint = MODAL_RECORD_3;
4034  fbin.write(&checkPoint, sizeof(checkPoint));
4035  }
4036 
4037  for (unsigned int jMode = 1; jMode <= NModesFEM; jMode++) {
4038  fdat >> d;
4039 
4040  if (bWriteBIN) {
4041  fbin.write((const char *)&d, sizeof(d));
4042  }
4043 
4044  if (!bActiveModes[jMode]) {
4045  continue;
4046  }
4047  a->Put(iCnt, d);
4048  iCnt++;
4049  }
4050 
4051  bRecordGroup[MODAL_RECORD_3] = true;
4052  } break;
4053 
4054  case MODAL_RECORD_4: {
4055  /* velocita' iniziali dei modi */
4056  if (bRecordGroup[MODAL_RECORD_4]) {
4057  silent_cerr("file=\"" << sFileFEM << "\": "
4058  "\"RECORD GROUP 4\" already parsed"
4059  << std::endl);
4060  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4061  }
4062 
4063  unsigned int iCnt = 1;
4064 
4065  if (bActiveModes.empty()) {
4066  silent_cerr("Modal(" << uLabel << "): "
4067  "input file \"" << sFileFEM << "\" "
4068  "is bogus (RECORD GROUP 4)"
4069  << std::endl);
4070  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4071  }
4072 
4073  if (bWriteBIN) {
4074  checkPoint = MODAL_RECORD_4;
4075  fbin.write(&checkPoint, sizeof(checkPoint));
4076  }
4077 
4078  for (unsigned int jMode = 1; jMode <= NModesFEM; jMode++) {
4079  fdat >> d;
4080 
4081  if (bWriteBIN) {
4082  fbin.write((const char *)&d, sizeof(d));
4083  }
4084 
4085  if (!bActiveModes[jMode]) {
4086  continue;
4087  }
4088  aP->Put(iCnt, d);
4089  iCnt++;
4090  }
4091 
4092  bRecordGroup[MODAL_RECORD_4] = true;
4093  } break;
4094 
4095  case MODAL_RECORD_5: {
4096  /* Coordinate X dei nodi */
4097  if (bRecordGroup[MODAL_RECORD_5]) {
4098  silent_cerr("file=\"" << sFileFEM << "\": "
4099  "\"RECORD GROUP 5\" already parsed"
4100  << std::endl);
4101  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4102  }
4103 
4104  if (bWriteBIN) {
4105  checkPoint = MODAL_RECORD_5;
4106  fbin.write(&checkPoint, sizeof(checkPoint));
4107  }
4108 
4109  for (unsigned int iNode = 1; iNode <= NFEMNodes; iNode++) {
4110  fdat >> d;
4111 
4112  if (bWriteBIN) {
4113  fbin.write((const char *)&d, sizeof(d));
4114  }
4115 
4116 #ifdef MODAL_SCALE_DATA
4117  d *= scalemodes;
4118 #endif /* MODAL_SCALE_DATA */
4119  pXYZFEMNodes->Put(1, iNode, d);
4120  }
4121 
4122  bRecordGroup[MODAL_RECORD_5] = true;
4123  } break;
4124 
4125  case MODAL_RECORD_6: {
4126  /* Coordinate Y dei nodi*/
4127  if (bRecordGroup[MODAL_RECORD_6]) {
4128  silent_cerr("file=\"" << sFileFEM << "\": "
4129  "\"RECORD GROUP 6\" already parsed"
4130  << std::endl);
4131  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4132  }
4133 
4134  if (bWriteBIN) {
4135  checkPoint = MODAL_RECORD_6;
4136  fbin.write(&checkPoint, sizeof(checkPoint));
4137  }
4138 
4139  for (unsigned int iNode = 1; iNode <= NFEMNodes; iNode++) {
4140  fdat >> d;
4141 
4142  if (bWriteBIN) {
4143  fbin.write((const char *)&d, sizeof(d));
4144  }
4145 
4146 #ifdef MODAL_SCALE_DATA
4147  d *= scalemodes;
4148 #endif /* MODAL_SCALE_DATA */
4149  pXYZFEMNodes->Put(2, iNode, d);
4150  }
4151 
4152  bRecordGroup[MODAL_RECORD_6] = true;
4153  } break;
4154 
4155  case MODAL_RECORD_7: {
4156  /* Coordinate Z dei nodi*/
4157  if (bRecordGroup[MODAL_RECORD_7]) {
4158  silent_cerr("file=\"" << sFileFEM << "\": "
4159  "\"RECORD GROUP 7\" already parsed"
4160  << std::endl);
4161  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4162  }
4163 
4164  if (bWriteBIN) {
4165  checkPoint = MODAL_RECORD_7;
4166  fbin.write(&checkPoint, sizeof(checkPoint));
4167  }
4168 
4169  for (unsigned int iNode = 1; iNode <= NFEMNodes; iNode++) {
4170  fdat >> d;
4171 
4172  if (bWriteBIN) {
4173  fbin.write((const char *)&d, sizeof(d));
4174  }
4175 
4176 #ifdef MODAL_SCALE_DATA
4177  d *= scalemodes;
4178 #endif /* MODAL_SCALE_DATA */
4179  pXYZFEMNodes->Put(3, iNode, d);
4180  }
4181 
4182  bRecordGroup[MODAL_RECORD_7] = true;
4183  } break;
4184 
4185  case MODAL_RECORD_8: {
4186  /* Forme modali */
4187  if (bRecordGroup[MODAL_RECORD_8]) {
4188  silent_cerr("file=\"" << sFileFEM << "\": "
4189  "\"RECORD GROUP 8\" already parsed"
4190  << std::endl);
4191  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4192  }
4193 
4194  if (bWriteBIN) {
4195  checkPoint = MODAL_RECORD_8;
4196  fbin.write(&checkPoint, sizeof(checkPoint));
4197  }
4198 
4199  // we assume rejected modes come first
4200  for (unsigned int jMode = 1; jMode <= NRejModes; jMode++) {
4201  // header
4202  fdat.getline(str, sizeof(str));
4203 
4204  silent_cout("Modal(" << uLabel << "): rejecting mode #" << jMode
4205  << " (" << str << ")" << std::endl);
4206 
4207  // mode
4208  for (unsigned int iNode = 1; iNode <= NFEMNodes; iNode++) {
4209  doublereal n[6];
4210 
4211  fdat >> n[0] >> n[1] >> n[2]
4212  >> n[3] >> n[4] >> n[5];
4213  }
4214 
4215  /* FIXME: siamo sicuri di avere
4216  * raggiunto '\n'? */
4217  fdat.getline(str, sizeof(str));
4218  }
4219 
4220  if (bActiveModes.empty()) {
4221  silent_cerr("Modal(" << uLabel << "): "
4222  "input file \"" << sFileFEM << "\" "
4223  "is bogus (RECORD GROUP 8)"
4224  << std::endl);
4225  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4226  }
4227 
4228  unsigned int iCnt = 1;
4229  for (unsigned int jMode = 1; jMode <= NModesFEM; jMode++) {
4230  fdat.getline(str, sizeof(str));
4231  if (str[0] != '\0' && strncmp("**", str, STRLENOF("**")) != 0)
4232  {
4233  if (NRejModes) {
4234  silent_cerr("Modal(" << uLabel << "): "
4235  "input file \"" << sFileFEM << "\""
4236  "is bogus (\"**\" expected before mode #" << jMode
4237  << " of " << NModesFEM << "; "
4238  "#" << jMode + NRejModes
4239  << " of " << NModesFEM + NRejModes
4240  << " including rejected)"
4241  << std::endl);
4242  } else {
4243  silent_cerr("Modal(" << uLabel << "): "
4244  "input file \"" << sFileFEM << "\""
4245  "is bogus (\"**\" expected before mode #" << jMode
4246  << " of " << NModesFEM << ")"
4247  << std::endl);
4248  }
4249  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4250  }
4251 
4252  for (unsigned int iNode = 1; iNode <= NFEMNodes; iNode++) {
4253  doublereal n[6];
4254 
4255  fdat >> n[0] >> n[1] >> n[2]
4256  >> n[3] >> n[4] >> n[5];
4257 
4258  if (bWriteBIN) {
4259  fbin.write((const char *)n, sizeof(n));
4260  }
4261 
4262  if (!bActiveModes[jMode]) {
4263  continue;
4264  }
4265 
4266 #ifdef MODAL_SCALE_DATA
4267  n[0] *= scalemodes;
4268  n[1] *= scalemodes;
4269  n[2] *= scalemodes;
4270 #endif /* MODAL_SCALE_DATA */
4271  pModeShapest->PutVec((iCnt - 1)*NFEMNodes + iNode, Vec3(&n[0]));
4272  pModeShapesr->PutVec((iCnt - 1)*NFEMNodes + iNode, Vec3(&n[3]));
4273  }
4274 
4275  if (bActiveModes[jMode]) {
4276  iCnt++;
4277  }
4278  fdat.getline(str, sizeof(str));
4279  }
4280 
4281  bRecordGroup[MODAL_RECORD_8] = true;
4282  } break;
4283 
4284  case MODAL_RECORD_9: {
4285  /* Matrice di massa modale */
4286  if (bRecordGroup[MODAL_RECORD_9]) {
4287  silent_cerr("file=\"" << sFileFEM << "\": "
4288  "\"RECORD GROUP 9\" already parsed"
4289  << std::endl);
4290  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4291  }
4292 
4293  if (bActiveModes.empty()) {
4294  silent_cerr("Modal(" << uLabel << "): "
4295  "input file \"" << sFileFEM << "\" "
4296  "is bogus (RECORD GROUP 9)"
4297  << std::endl);
4298  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4299  }
4300 
4301  if (bWriteBIN) {
4302  checkPoint = MODAL_RECORD_9;
4303  fbin.write(&checkPoint, sizeof(checkPoint));
4304  }
4305 
4306  unsigned int iCnt = 1;
4307  for (unsigned int jMode = 1; jMode <= NModesFEM; jMode++) {
4308  unsigned int jCnt = 1;
4309 
4310  for (unsigned int kMode = 1; kMode <= NModesFEM; kMode++) {
4311  fdat >> d;
4312 
4313  if (dMassThreshold > 0.0 && fabs(d) < dMassThreshold) {
4314  d = 0.0;
4315  }
4316 
4317  if (bWriteBIN) {
4318  fbin.write((const char *)&d, sizeof(d));
4319  }
4320 
4321  if (!bActiveModes[jMode] || !bActiveModes[kMode]) {
4322  continue;
4323  }
4324  pGenMass->Put(iCnt, jCnt, d);
4325  jCnt++;
4326  }
4327 
4328  if (bActiveModes[jMode]) {
4329  iCnt++;
4330  }
4331  }
4332 
4333  bRecordGroup[MODAL_RECORD_9] = true;
4334  } break;
4335 
4336  case MODAL_RECORD_10: {
4337  /* Matrice di rigidezza modale */
4338  if (bRecordGroup[MODAL_RECORD_10]) {
4339  silent_cerr("file=\"" << sFileFEM << "\": "
4340  "\"RECORD GROUP 10\" already parsed"
4341  << std::endl);
4342  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4343  }
4344 
4345  if (bActiveModes.empty()) {
4346  silent_cerr("Modal(" << uLabel << "): "
4347  "input file \"" << sFileFEM << "\" "
4348  "is bogus (RECORD GROUP 10)"
4349  << std::endl);
4350  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4351  }
4352 
4353  if (bWriteBIN) {
4354  checkPoint = MODAL_RECORD_10;
4355  fbin.write(&checkPoint, sizeof(checkPoint));
4356  }
4357 
4358  unsigned int iCnt = 1;
4359  for (unsigned int jMode = 1; jMode <= NModesFEM; jMode++) {
4360  unsigned int jCnt = 1;
4361 
4362  for (unsigned int kMode = 1; kMode <= NModesFEM; kMode++) {
4363  fdat >> d;
4364 
4365  if (dStiffnessThreshold > 0.0 && fabs(d) < dStiffnessThreshold) {
4366  d = 0.0;
4367  }
4368 
4369  if (bWriteBIN) {
4370  fbin.write((const char *)&d, sizeof(d));
4371  }
4372 
4373  if (!bActiveModes[jMode] || !bActiveModes[kMode]) {
4374  continue;
4375  }
4376  pGenStiff->Put(iCnt, jCnt, d);
4377  jCnt++;
4378  }
4379 
4380  if (bActiveModes[jMode]) {
4381  iCnt++;
4382  }
4383  }
4384 
4385  bRecordGroup[MODAL_RECORD_10] = true;
4386  } break;
4387 
4388  case MODAL_RECORD_11: {
4389  /* Lumped Masses */
4390  if (bRecordGroup[MODAL_RECORD_11]) {
4391  silent_cerr("file=\"" << sFileFEM << "\": "
4392  "\"RECORD GROUP 11\" already parsed"
4393  << std::endl);
4394  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4395  }
4396 
4397  if (bWriteBIN) {
4398  checkPoint = MODAL_RECORD_11;
4399  fbin.write(&checkPoint, sizeof(checkPoint));
4400  }
4401 
4402  FEMMass.resize(NFEMNodes);
4403  FEMJ.resize(NFEMNodes);
4404 
4405  for (unsigned int iNode = 1; iNode <= NFEMNodes; iNode++) {
4406  doublereal tmpmass = 0.;
4407  for (unsigned int jCnt = 1; jCnt <= 6; jCnt++) {
4408  fdat >> d;
4409 
4410  if (bWriteBIN) {
4411  fbin.write((const char *)&d, sizeof(d));
4412  }
4413 
4414  switch (jCnt) {
4415  case 1:
4416  tmpmass = d;
4417 #ifdef MODAL_SCALE_DATA
4418  d *= scalemass;
4419 #endif /* MODAL_SCALE_DATA */
4420  FEMMass[iNode - 1] = d;
4421  break;
4422 
4423  case 2:
4424  case 3:
4425  if (d != tmpmass) {
4426  silent_cout("Warning: component(" << jCnt << ") "
4427  "of FEM node(" << iNode << ") mass "
4428  "differs from component(1)=" << tmpmass << std::endl);
4429  }
4430  break;
4431 
4432  case 4:
4433  case 5:
4434  case 6:
4435 #ifdef MODAL_SCALE_DATA
4436  d *= scaleinertia;
4437 #endif /* MODAL_SCALE_DATA */
4438  FEMJ[iNode - 1](jCnt - 3) = d;
4439  break;
4440  }
4441  }
4442  }
4443 
4444  bBuildInvariants = true;
4445 
4446  bRecordGroup[MODAL_RECORD_11] = true;
4447  } break;
4448 
4449  case MODAL_RECORD_12: {
4450  /* Rigid body inertia */
4451  if (bRecordGroup[MODAL_RECORD_12]) {
4452  silent_cerr("file=\"" << sFileFEM << "\": "
4453  "\"RECORD GROUP 12\" already parsed"
4454  << std::endl);
4455  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4456  }
4457 
4458  if (bWriteBIN) {
4459  checkPoint = MODAL_RECORD_12;
4460  fbin.write(&checkPoint, sizeof(checkPoint));
4461  }
4462 
4463  fdat >> d;
4464 
4465  if (bWriteBIN) {
4466  fbin.write((const char *)&d, sizeof(d));
4467  }
4468 
4469  dMass = d;
4470 
4471  // NOTE: the CM location is read and temporarily stored in STmp
4472  // later, it will be multiplied by the mass
4473  for (int iRow = 1; iRow <= 3; iRow++) {
4474  fdat >> d;
4475 
4476  if (bWriteBIN) {
4477  fbin.write((const char *)&d, sizeof(d));
4478  }
4479 
4480  XTmpIn(iRow) = d;
4481  }
4482 
4483  // here JTmp is with respect to the center of mass
4484  for (int iRow = 1; iRow <= 3; iRow++) {
4485  for (int iCol = 1; iCol <= 3; iCol++) {
4486  fdat >> d;
4487 
4488  if (bWriteBIN) {
4489  fbin.write((const char *)&d, sizeof(d));
4490  }
4491 
4492  JTmpIn(iRow, iCol) = d;
4493  }
4494  }
4495 
4496  bRecordGroup[MODAL_RECORD_12] = true;
4497  } break;
4498 
4499  case MODAL_RECORD_13: {
4500  /* Matrice di smorzamento modale */
4501  if (bRecordGroup[MODAL_RECORD_13]) {
4502  silent_cerr("file=\"" << sFileFEM << "\": "
4503  "\"RECORD GROUP 13\" already parsed"
4504  << std::endl);
4505  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4506  }
4507 
4508  if (bActiveModes.empty()) {
4509  silent_cerr("Modal(" << uLabel << "): "
4510  "input file \"" << sFileFEM << "\" "
4511  "is bogus (RECORD GROUP 13)"
4512  << std::endl);
4513  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4514  }
4515 
4516  if (bWriteBIN) {
4517  checkPoint = MODAL_RECORD_13;
4518  fbin.write(&checkPoint, sizeof(checkPoint));
4519  }
4520 
4521  unsigned int iCnt = 1;
4522  for (unsigned int jMode = 1; jMode <= NModesFEM; jMode++) {
4523  unsigned int jCnt = 1;
4524 
4525  for (unsigned int kMode = 1; kMode <= NModesFEM; kMode++) {
4526  fdat >> d;
4527 
4528  if (dDampingThreshold > 0.0 && fabs(d) < dDampingThreshold) {
4529  d = 0.0;
4530  }
4531 
4532  if (bWriteBIN) {
4533  fbin.write((const char *)&d, sizeof(d));
4534  }
4535 
4536  if (!bActiveModes[jMode] || !bActiveModes[kMode]) {
4537  continue;
4538  }
4539 
4540  if (eDamp == DAMPING_FROM_FILE) {
4541  pGenDamp->Put(iCnt, jCnt, d);
4542  }
4543 
4544  jCnt++;
4545  }
4546 
4547  if (bActiveModes[jMode]) {
4548  iCnt++;
4549  }
4550  }
4551 
4552  bRecordGroup[MODAL_RECORD_13] = true;
4553  } break;
4554 
4555  default:
4556  silent_cerr("file=\"" << sFileFEM << "\": "
4557  "\"RECORD GROUP " << rg << "\" unhandled"
4558  << std::endl);
4559  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4560  }
4561  }
4562 
4563  fdat.close();
4564  if (bWriteBIN) {
4565  checkPoint = MODAL_END_OF_FILE;
4566  fbin.write(&checkPoint, sizeof(checkPoint));
4567  fbin.close();
4568  }
4569 
4570  /* unlink binary file if create/update failed */
4571  } catch (...) {
4572  if (bWriteBIN) {
4573  // NOTE: might be worth leaving in place
4574  // for debugging purposes
4575  (void)unlink(sBinFileFEM.c_str());
4576  }
4577  throw;
4578  }
4579 
4580  fname = sFileFEM;
4581 
4582  } else {
4583  std::ifstream fbin(sBinFileFEM.c_str(), std::ios::binary);
4584  if (!fbin) {
4585  silent_cerr("Modal(" << uLabel << "): "
4586  "unable to open file \"" << sBinFileFEM << "\""
4587  "at line " << HP.GetLineData() << std::endl);
4588  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
4589  }
4590  silent_cout("Modal(" << uLabel << "): "
4591  "reading flexible body data from file "
4592  "\"" << sBinFileFEM << "\"" << std::endl);
4593 
4594  /* check version */
4595  // NOTE: any time the format of the binary file changes,
4596  // the BinVersion should be incremented.
4597  // When reading the binary file, the version of the file
4598  // is compared with that of the code. If they do not match:
4599  // if the file is newer than the code, abort.
4600  // if the file is older than the code, try to deal with it
4601  // as much as possible
4602  //
4603  // History:
4604  //
4605  // 1: first implementation
4606  //
4607  // 2: July 2008, at LaMSID
4608  // - string FEM labels
4609  // - 3 & 4 optional
4610  // - 11 & 12 optional and no longer mutually exclusive
4611  // - -1 as the last checkpoint
4612  // - version 1 tolerated
4613  // 3: November 2012
4614  // - record 13, generalized damping matrix
4615  // 4: December 2012; incompatible with previous ones
4616  // - bin file portable (fixed size types, magic signature)
4617  char currMagic[5] = { 0 };
4618  fbin.read((char *)&currMagic, sizeof(4));
4619  if (memcmp(magic, currMagic, 4) != 0) {
4620  silent_cerr("Modal(" << uLabel << "): "
4621  "file \"" << sBinFileFEM << "\" "
4622  "unexpected signature" << std::endl);
4623  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4624  }
4625 
4626  fbin.read((char *)&currBinVersion, sizeof(currBinVersion));
4627  if (currBinVersion > BinVersion) {
4628  silent_cerr("Modal(" << uLabel << "): "
4629  "file \"" << sBinFileFEM << "\" "
4630  "version " << currBinVersion << " unsupported" << std::endl);
4631  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4632 
4633  } else if (currBinVersion < BinVersion) {
4634  silent_cout("Modal(" << uLabel << "): "
4635  "file \"" << sBinFileFEM << "\" "
4636  "version " << currBinVersion << " is outdated; "
4637  "trying to read anyway..." << std::endl);
4638  }
4639 
4640  pedantic_cout("Modal(" << uLabel << "): "
4641  "binary version " << currBinVersion << std::endl);
4642 
4643  // NOTE: the binary file is expected to be in order;
4644  // however, if generated from a non-ordered textual file
4645  // it will be out of order... perhaps also the textual
4646  // file should always be ordered?
4647 
4648  /* legge il primo blocco (HEADER) */
4649  fbin.read(&checkPoint, sizeof(checkPoint));
4650  if (checkPoint != MODAL_RECORD_1) {
4651  silent_cerr("Modal(" << uLabel << "): "
4652  "file \"" << sBinFileFEM << "\" "
4653  "looks broken (expecting checkpoint 1)"
4654  << std::endl);
4655  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4656  }
4657 
4658  pedantic_cout("Modal(" << uLabel << "): "
4659  "reading block " << unsigned(checkPoint) << std::endl);
4660 
4661  fbin.read((char *)&NFEMNodesFEM, sizeof(NFEMNodesFEM));
4662  fbin.read((char *)&NModesFEM, sizeof(NModesFEM));
4663 
4664  /* consistency checks */
4665  if (NFEMNodes == unsigned(-1)) {
4666  NFEMNodes = NFEMNodesFEM;
4667 
4668  } else if (NFEMNodes != NFEMNodesFEM) {
4669  silent_cerr("Modal(" << uLabel << "), "
4670  "file \"" << sFileFEM << "\": "
4671  "FEM nodes number " << NFEMNodes
4672  << " does not match node number "
4673  << NFEMNodesFEM
4674  << std::endl);
4675  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
4676  }
4677 
4678  if (NModes != NModesFEM) {
4679  silent_cout("Modal(" << uLabel
4680  << "), file '" << sFileFEM
4681  << "': using " << NModes
4682  << " of " << NModesFEM
4683  << " modes" << std::endl);
4684  }
4685 
4686  if (!bActiveModes.empty()) {
4687  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4688  }
4689 
4690  IdFEMNodes.resize(NFEMNodes);
4691 
4692  SAFENEWWITHCONSTRUCTOR(pXYZFEMNodes, Mat3xN, Mat3xN(NFEMNodes, 0.));
4693  SAFENEWWITHCONSTRUCTOR(pModeShapest, Mat3xN, Mat3xN(NFEMNodes*NModes, 0.));
4694  SAFENEWWITHCONSTRUCTOR(pModeShapesr, Mat3xN, Mat3xN(NFEMNodes*NModes, 0.));
4695 
4696  bActiveModes.resize(NModesFEM + 1, false);
4697 
4698  for (unsigned int iCnt = 0; iCnt < NModes; iCnt++) {
4699  if (uModeNumber[iCnt] > NModesFEM) {
4700  silent_cerr("Modal(" << uLabel << "): "
4701  "mode " << uModeNumber[iCnt]
4702  << " is not available (max = "
4703  << NModesFEM << ")"
4704  << std::endl);
4705  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4706  }
4707  bActiveModes[uModeNumber[iCnt]] = true;
4708  }
4709 
4710  bRecordGroup[1] = true;
4711 
4712  for (;;) {
4713  fbin.read(&checkPoint, sizeof(checkPoint));
4714 
4715  if (checkPoint == MODAL_END_OF_FILE) {
4716  pedantic_cout("Modal(" << uLabel << "): "
4717  "end of binary file \"" << sBinFileFEM << "\""
4718  << std::endl);
4719  fbin.close();
4720  break;
4721  }
4722 
4723  if (checkPoint < MODAL_RECORD_1 || checkPoint >= MODAL_LAST_RECORD) {
4724  silent_cerr("Modal(" << uLabel << "): "
4725  "file \"" << sBinFileFEM << "\" "
4726  "looks broken (unknown block " << unsigned(checkPoint) << ")"
4727  << std::endl);
4728  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4729  }
4730 
4731  if (bRecordGroup[unsigned(checkPoint)]) {
4732  silent_cerr("Modal(" << uLabel << "): "
4733  "file \"" << sBinFileFEM << "\" "
4734  "looks broken (block " << unsigned(checkPoint) << " already parsed)"
4735  << std::endl);
4736  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4737  }
4738 
4739  pedantic_cout("Modal(" << uLabel << "): "
4740  "reading block " << unsigned(checkPoint)
4741  << " from file \"" << sBinFileFEM << "\"" << std::endl);
4742 
4743  switch (checkPoint) {
4744  case MODAL_RECORD_2:
4745  /* legge il secondo blocco (Id.nodi) */
4746  for (unsigned int iNode = 0; iNode < NFEMNodes; iNode++) {
4747  uint32_t len;
4748  fbin.read((char *)&len, sizeof(len));
4749  ASSERT(len > 0);
4750  IdFEMNodes[iNode].resize(len);
4751  fbin.read((char *)IdFEMNodes[iNode].c_str(), len);
4752  }
4753  break;
4754 
4755  case MODAL_RECORD_3:
4756  /* deformate iniziali dei modi */
4757  for (unsigned int iCnt = 1, jMode = 1; jMode <= NModesFEM; jMode++) {
4758  doublereal d;
4759 
4760  fbin.read((char *)&d, sizeof(d));
4761 
4762  if (!bActiveModes[jMode]) {
4763  continue;
4764  }
4765 
4766  a->Put(iCnt, d);
4767  iCnt++;
4768  }
4769  break;
4770 
4771  case MODAL_RECORD_4:
4772  /* velocita' iniziali dei modi */
4773  for (unsigned int iCnt = 1, jMode = 1; jMode <= NModesFEM; jMode++) {
4774  doublereal d;
4775 
4776  fbin.read((char *)&d, sizeof(d));
4777 
4778  if (!bActiveModes[jMode]) {
4779  continue;
4780  }
4781 
4782  aP->Put(iCnt, d);
4783  iCnt++;
4784  }
4785  break;
4786 
4787  case MODAL_RECORD_5:
4788  /* Coordinate X dei nodi */
4789  for (unsigned int iNode = 1; iNode <= NFEMNodes; iNode++) {
4790  doublereal d;
4791 
4792  fbin.read((char *)&d, sizeof(d));
4793 
4794 #ifdef MODAL_SCALE_DATA
4795  d *= scalemodes;
4796 #endif /* MODAL_SCALE_DATA */
4797 
4798  pXYZFEMNodes->Put(1, iNode, d);
4799  }
4800  break;
4801 
4802  case MODAL_RECORD_6:
4803  /* Coordinate Y dei nodi*/
4804  for (unsigned int iNode = 1; iNode <= NFEMNodes; iNode++) {
4805  doublereal d;
4806 
4807  fbin.read((char *)&d, sizeof(d));
4808 
4809 #ifdef MODAL_SCALE_DATA
4810  d *= scalemodes;
4811 #endif /* MODAL_SCALE_DATA */
4812 
4813  pXYZFEMNodes->Put(2, iNode, d);
4814  }
4815  break;
4816 
4817  case MODAL_RECORD_7:
4818  /* Coordinate Z dei nodi*/
4819  for (unsigned int iNode = 1; iNode <= NFEMNodes; iNode++) {
4820  doublereal d;
4821 
4822  fbin.read((char *)&d, sizeof(d));
4823 
4824 #ifdef MODAL_SCALE_DATA
4825  d *= scalemodes;
4826 #endif /* MODAL_SCALE_DATA */
4827 
4828  pXYZFEMNodes->Put(3, iNode, d);
4829  }
4830  break;
4831 
4832  case MODAL_RECORD_8:
4833  /* Forme modali */
4834  for (unsigned int iCnt = 1, jMode = 1; jMode <= NModesFEM; jMode++) {
4835  for (unsigned int iNode = 1; iNode <= NFEMNodes; iNode++) {
4836  doublereal n[6];
4837 
4838  fbin.read((char *)n, sizeof(n));
4839 
4840  if (!bActiveModes[jMode]) {
4841  continue;
4842  }
4843 
4844 #ifdef MODAL_SCALE_DATA
4845  n[0] *= scalemodes;
4846  n[1] *= scalemodes;
4847  n[2] *= scalemodes;
4848 #endif /* MODAL_SCALE_DATA */
4849  pModeShapest->PutVec((iCnt - 1)*NFEMNodes + iNode, Vec3(&n[0]));
4850  pModeShapesr->PutVec((iCnt - 1)*NFEMNodes + iNode, Vec3(&n[3]));
4851  }
4852 
4853  if (bActiveModes[jMode]) {
4854  iCnt++;
4855  }
4856  }
4857  break;
4858 
4859  case MODAL_RECORD_9:
4860  /* Matrice di massa modale */
4861  for (unsigned int iCnt = 1, jMode = 1; jMode <= NModesFEM; jMode++) {
4862  unsigned int jCnt = 1;
4863 
4864  for (unsigned int kMode = 1; kMode <= NModesFEM; kMode++) {
4865  doublereal d;
4866 
4867  fbin.read((char *)&d, sizeof(d));
4868 
4869  if (!bActiveModes[jMode] || !bActiveModes[kMode]) {
4870  continue;
4871  }
4872  pGenMass->Put(iCnt, jCnt, d);
4873  jCnt++;
4874  }
4875 
4876  if (bActiveModes[jMode]) {
4877  iCnt++;
4878  }
4879  }
4880  break;
4881 
4882  case MODAL_RECORD_10:
4883  /* Matrice di rigidezza modale */
4884  for (unsigned int iCnt = 1, jMode = 1; jMode <= NModesFEM; jMode++) {
4885  unsigned int jCnt = 1;
4886 
4887  for (unsigned int kMode = 1; kMode <= NModesFEM; kMode++) {
4888  doublereal d;
4889 
4890  fbin.read((char *)&d, sizeof(d));
4891 
4892  if (!bActiveModes[jMode] || !bActiveModes[kMode]) {
4893  continue;
4894  }
4895  pGenStiff->Put(iCnt, jCnt, d);
4896  jCnt++;
4897  }
4898 
4899  if (bActiveModes[jMode]) {
4900  iCnt++;
4901  }
4902  }
4903  break;
4904 
4905  case MODAL_RECORD_11:
4906  /* Lumped Masses */
4907  FEMMass.resize(NFEMNodes);
4908  FEMJ.resize(NFEMNodes);
4909 
4910  for (unsigned int iNode = 1; iNode <= NFEMNodes; iNode++) {
4911  for (unsigned int jCnt = 1; jCnt <= 6; jCnt++) {
4912  doublereal d;
4913 
4914  fbin.read((char *)&d, sizeof(d));
4915 
4916  switch (jCnt) {
4917  case 1:
4918 #ifdef MODAL_SCALE_DATA
4919  d *= scalemass;
4920 #endif /* MODAL_SCALE_DATA */
4921  FEMMass[iNode - 1] = d;
4922  break;
4923 
4924  case 4:
4925  case 5:
4926  case 6:
4927 #ifdef MODAL_SCALE_DATA
4928  d *= scaleinertia;
4929 #endif /* MODAL_SCALE_DATA */
4930  FEMJ[iNode - 1](jCnt - 3) = d;
4931  break;
4932  }
4933  }
4934  }
4935 
4936  bBuildInvariants = true;
4937  break;
4938 
4939  case MODAL_RECORD_12: {
4940  doublereal d;
4941  fbin.read((char *)&d, sizeof(d));
4942  dMass = d;
4943 
4944  // NOTE: the CM location is read and temporarily stored in XTmpIn
4945  // later, it will be multiplied by the mass
4946  for (int iRow = 1; iRow <= 3; iRow++) {
4947  fbin.read((char *)&d, sizeof(d));
4948 
4949  XTmpIn(iRow) = d;
4950  }
4951 
4952  // here JTmpIn is with respect to the center of mass
4953  for (int iRow = 1; iRow <= 3; iRow++) {
4954  for (int iCol = 1; iCol <= 3; iCol++) {
4955  fbin.read((char *)&d, sizeof(d));
4956 
4957  JTmpIn(iRow, iCol) = d;
4958  }
4959  }
4960  } break;
4961 
4962  case MODAL_RECORD_13:
4963  for (unsigned int iCnt = 1, jMode = 1; jMode <= NModesFEM; jMode++) {
4964  unsigned int jCnt = 1;
4965 
4966  for (unsigned int kMode = 1; kMode <= NModesFEM; kMode++) {
4967  doublereal d;
4968 
4969  fbin.read((char *)&d, sizeof(d));
4970 
4971  if (!bActiveModes[jMode] || !bActiveModes[kMode]) {
4972  continue;
4973  }
4974 
4975  if (eDamp == DAMPING_FROM_FILE) {
4976  pGenDamp->Put(iCnt, jCnt, d);
4977  }
4978 
4979  jCnt++;
4980  }
4981 
4982  if (bActiveModes[jMode]) {
4983  iCnt++;
4984  }
4985  }
4986  break;
4987 
4988  default:
4989  silent_cerr("Modal(" << uLabel << "): "
4990  "file \"" << sBinFileFEM << "\" "
4991  "looks broken (unknown block " << unsigned(checkPoint) << ")"
4992  << std::endl);
4993  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
4994  }
4995 
4996  bRecordGroup[unsigned(checkPoint)] = true;
4997  }
4998 
4999  fname = sBinFileFEM;
5000  }
5001 
5002  unsigned reqMR[] = {
5003  MODAL_RECORD_1,
5004  MODAL_RECORD_2,
5005  // 3 & 4 no longer required; explicit check is present when currBinVersion == 1
5006  MODAL_RECORD_5,
5007  MODAL_RECORD_6,
5008  MODAL_RECORD_7,
5009  MODAL_RECORD_8,
5010  MODAL_RECORD_9,
5011  MODAL_RECORD_10
5012  };
5013  bool bBailOut(false);
5014  for (unsigned iCnt = 0; iCnt < sizeof(reqMR)/sizeof(unsigned); iCnt++) {
5015  if (!bRecordGroup[reqMR[iCnt]]) {
5016  silent_cerr("Modal(" << uLabel << "): "
5017  "incomplete input file \"" << fname << "\", "
5018  "record group " << reqMR[iCnt] << " missing "
5019  "at line " << HP.GetLineData() << std::endl);
5020  bBailOut = true;
5021  }
5022  }
5023 
5024  if (bBailOut) {
5025  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
5026  }
5027 
5028  // JTmp is now referred to the origin
5029  JTmp = JTmpIn - Mat3x3(MatCrossCross, XTmpIn, XTmpIn*dMass);
5030 
5031  if (!sEchoFileName.empty()) {
5032  std::ofstream fecho(sEchoFileName.c_str());
5033 
5034  fecho.setf(std::ios::scientific);
5035  fecho.precision(iEchoPrecision);
5036 
5037  time_t t;
5038  time(&t);
5039 
5040  fecho
5041  << "** echo of file \"" << fname << "\" generated " << ctime(&t)
5042  << "** RECORD GROUP 1, HEADER" << std::endl
5043  << "** REVISION. NODES. NORMAL, ATTACHMENT, CONSTRAINT, REJECTED MODES." << std::endl
5044  << "VER" << currBinVersion << " " << NFEMNodes << " " << NModes << " 0 0 0" << std::endl
5045  << "**" << std::endl
5046  << "** RECORD GROUP 2, FINITE ELEMENT NODE LIST" << std::endl;
5047  for (unsigned r = 0; r <= (IdFEMNodes.size() - 1)/6; r++) {
5048  for (unsigned c = 0; c < std::min(6U, unsigned(IdFEMNodes.size() - 6*r)); c++) {
5049  fecho << IdFEMNodes[6*r + c] << " ";
5050  }
5051  fecho << std::endl;
5052  }
5053 
5054  fecho
5055  << "**" << std::endl
5056  << "** RECORD GROUP 3, INITIAL MODAL DISPLACEMENTS"<< std::endl;
5057  for (int r = 1; r <= a->iGetNumRows(); r++) {
5058  fecho << (*a)(r) << std::endl;
5059  }
5060 
5061  fecho
5062  << "**" << std::endl
5063  << "** RECORD GROUP 4, INITIAL MODAL VELOCITIES"<< std::endl;
5064  for (int r = 1; r <= aP->iGetNumRows(); r++) {
5065  fecho << (*aP)(r) << std::endl;
5066  }
5067 
5068  fecho
5069  << "**" << std::endl
5070  << "** RECORD GROUP 5, NODAL X COORDINATES" << std::endl;
5071  for (int r = 1; r <= pXYZFEMNodes->iGetNumCols(); r++) {
5072  fecho << (*pXYZFEMNodes)(1, r) << std::endl;
5073  }
5074 
5075  fecho
5076  << "**" << std::endl
5077  << "** RECORD GROUP 6, NODAL Y COORDINATES" << std::endl;
5078  for (int r = 1; r <= pXYZFEMNodes->iGetNumCols(); r++) {
5079  fecho << (*pXYZFEMNodes)(2, r) << std::endl;
5080  }
5081 
5082  fecho
5083  << "**" << std::endl
5084  << "** RECORD GROUP 7, NODAL Z COORDINATES" << std::endl;
5085  for (int r = 1; r <= pXYZFEMNodes->iGetNumCols(); r++) {
5086  fecho << (*pXYZFEMNodes)(3, r) << std::endl;
5087  }
5088 
5089  fecho
5090  << "**" << std::endl
5091  << "** RECORD GROUP 8, NON-ORTHOGONALIZED MODE SHAPES" << std::endl;
5092  for (unsigned m = 0; m < NModes; m++) {
5093  fecho
5094  << "** NORMAL MODE SHAPE # " << uModeNumber[m] << std::endl;
5095  for (unsigned n = 1; n <= NFEMNodes; n++) {
5096  fecho
5097  << (*pModeShapest)(1, m*NFEMNodes + n) << " "
5098  << (*pModeShapest)(2, m*NFEMNodes + n) << " "
5099  << (*pModeShapest)(3, m*NFEMNodes + n) << " "
5100  << (*pModeShapesr)(1, m*NFEMNodes + n) << " "
5101  << (*pModeShapesr)(2, m*NFEMNodes + n) << " "
5102  << (*pModeShapesr)(3, m*NFEMNodes + n) << std::endl;
5103  }
5104  }
5105 
5106  fecho
5107  << "**" << std::endl
5108  << "** RECORD GROUP 9, MODAL MASS MATRIX. COLUMN-MAJOR FORM" << std::endl;
5109  for (unsigned r = 1; r <= NModes; r++) {
5110  for (unsigned c = 1; c <= NModes; c++) {
5111  fecho << (*pGenMass)(r, c) << " ";
5112  }
5113  fecho << std::endl;
5114  }
5115 
5116  fecho
5117  << "**" << std::endl
5118  << "** RECORD GROUP 10, MODAL STIFFNESS MATRIX. COLUMN-MAJOR FORM" << std::endl;
5119  for (unsigned r = 1; r <= NModes; r++) {
5120  for (unsigned c = 1; c <= NModes; c++) {
5121  fecho << (*pGenStiff)(r, c) << " ";
5122  }
5123  fecho << std::endl;
5124  }
5125 
5126  if (bBuildInvariants) {
5127  fecho
5128  << "**" << std::endl
5129  << "** RECORD GROUP 11, DIAGONAL OF LUMPED MASS MATRIX" << std::endl;
5130  for (unsigned r = 0; r < FEMMass.size(); r++) {
5131  fecho
5132  << FEMMass[r] << " " << FEMMass[r] << " " << FEMMass[r] << " "
5133  << FEMJ[r] << std::endl;
5134  }
5135  }
5136 
5137  if (bRecordGroup[MODAL_RECORD_12]) {
5138  fecho
5139  << "**" << std::endl
5140  << "** RECORD GROUP 12, GLOBAL INERTIA" << std::endl
5141  << dMass << std::endl
5142  << XTmpIn << std::endl,
5143  JTmpIn.Write(fecho, " ", "\n") << std::endl;
5144  }
5145 
5146  if (bRecordGroup[MODAL_RECORD_13]) {
5147  fecho
5148  << "**" << std::endl
5149  << "** RECORD GROUP 13, MODAL DAMPING MATRIX. COLUMN-MAJOR FORM" << std::endl;
5150  for (unsigned r = 1; r <= NModes; r++) {
5151  for (unsigned c = 1; c <= NModes; c++) {
5152  fecho << (*pGenDamp)(r, c) << " ";
5153  }
5154  fecho << std::endl;
5155  }
5156  }
5157  }
5158 
5159  /* lettura dati di vincolo:
5160  * l'utente specifica la label del nodo FEM e del nodo rigido
5161  * d'interfaccia.
5162  * L'orientamento del nodo FEM e' quello del nodo modale, la
5163  * posizione e' la somma di quella modale e di quella FEM */
5164 
5165  /* array contenente le forme modali dei nodi d'interfaccia */
5166  Mat3xN* pPHItStrNode = 0;
5167  Mat3xN* pPHIrStrNode = 0;
5168 
5169  /* traslazione origine delle coordinate */
5170  Vec3 Origin(::Zero3);
5171  bool bOrigin(false);
5172 
5173  if (HP.IsKeyWord("origin" "node")) {
5174  /* numero di nodi FEM del modello */
5175  std::string FEMOriginNode;
5176  if (HP.IsStringWithDelims()) {
5177  FEMOriginNode = HP.GetStringWithDelims();
5178 
5179  } else {
5180  pedantic_cerr("Modal(" << uLabel << "): "
5181  "origin node expected as string with delimiters"
5182  << std::endl);
5183  FEMOriginNode = HP.GetString("");
5184  }
5185 
5186  unsigned int iNode;
5187  for (iNode = 0; iNode < NFEMNodes; iNode++) {
5188  if (FEMOriginNode == IdFEMNodes[iNode]) {
5189  break;
5190  }
5191  }
5192 
5193  if (iNode == NFEMNodes) {
5194  silent_cerr("Modal(" << uLabel << "): "
5195  "FEM node \"" << FEMOriginNode << "\""
5196  << " at line " << HP.GetLineData()
5197  << " not defined " << std::endl);
5198  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
5199  }
5200 
5201  iNode++;
5202  Origin = pXYZFEMNodes->GetVec(iNode);
5203 
5204  pedantic_cout("Modal(" << uLabel << "): origin x={" << Origin << "}" << std::endl);
5205  bOrigin = true;
5206 
5207  } else if (HP.IsKeyWord("origin" "position")) {
5208  Origin = HP.GetPosAbs(::AbsRefFrame);
5209  bOrigin = true;
5210  }
5211 
5212  if (bOrigin) {
5213  for (unsigned int iStrNode = 1; iStrNode <= NFEMNodes; iStrNode++) {
5214  pXYZFEMNodes->SubVec(iStrNode, Origin);
5215  }
5216 
5217  if (!bBuildInvariants) {
5218  XTmpIn -= Origin;
5219  JTmp = JTmpIn - Mat3x3(MatCrossCross, XTmpIn, XTmpIn*dMass);
5220  }
5221  }
5222 
5223  /* numero di nodi d'interfaccia */
5224  unsigned int NStrNodes = HP.GetInt();
5225  DEBUGCOUT("Number of Interface Nodes : " << NStrNodes << std::endl);
5226 
5227  SAFENEWWITHCONSTRUCTOR(pPHItStrNode, Mat3xN,
5228  Mat3xN(NStrNodes*NModes, 0.));
5229  SAFENEWWITHCONSTRUCTOR(pPHIrStrNode, Mat3xN,
5230  Mat3xN(NStrNodes*NModes, 0.));
5231 
5232  std::vector<Modal::StrNodeData> SND(NStrNodes);
5233 
5234  bool bOut = false;
5235  if (HP.IsKeyWord("output")) {
5236  bOut = HP.GetYesNoOrBool();
5237  }
5238 
5239  for (unsigned int iStrNode = 1; iStrNode <= NStrNodes; iStrNode++) {
5240  /* nodo collegato 1 (e' il nodo FEM) */
5241  std::string Node1 = HP.GetString("");
5242  // NOTE: we should check whether a label includes
5243  // whitespace. In any case, it wouldn't match
5244  // any of the labels read, so it is pointless,
5245  // but could help debugging...
5246  if (Node1.find(' ') != std::string::npos ) {
5247  silent_cout("Modal(" << uLabel << "): "
5248  "FEM node \"" << Node1 << "\""
5249  << " at line " << HP.GetLineData()
5250  << " contains a blank" << std::endl);
5251  }
5252 
5253  DEBUGCOUT("Linked to FEM Node " << Node1 << std::endl);
5254 
5255  unsigned int iNode;
5256  /* verifica di esistenza del nodo 1*/
5257  for (iNode = 0; iNode < NFEMNodes; iNode++) {
5258  if (Node1 == IdFEMNodes[iNode]) {
5259  break;
5260  }
5261  }
5262 
5263  if (iNode == NFEMNodes) {
5264  silent_cerr("Modal(" << uLabel << "): "
5265  "FEM node \"" << Node1 << "\""
5266  << " at line " << HP.GetLineData()
5267  << " not defined " << std::endl);
5268  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
5269  }
5270 
5271  iNode++;
5272 
5273  int iNodeCurr = iNode;
5274 
5275  /* recupera la posizione del nodo FEM, somma di posizione
5276  * e eventuale offset;
5277  *
5278  * HEADS UP: non piu' offset per i nodi FEM !!!!!!!!!
5279  *
5280  * nota: iNodeCurr contiene la posizione a cui si trova
5281  * il nodo FEM d'interfaccia nell'array pXYZNodes */
5282  SND[iStrNode-1].OffsetFEM = pXYZFEMNodes->GetVec(iNodeCurr);
5283 
5284  /* salva le forme modali del nodo d'interfaccia
5285  * nell'array pPHIStrNode */
5286  for (unsigned int jMode = 0; jMode < NModes; jMode++) {
5287  pPHItStrNode->PutVec(jMode*NStrNodes + iStrNode,
5288  pModeShapest->GetVec(jMode*NFEMNodes + iNodeCurr));
5289  pPHIrStrNode->PutVec(jMode*NStrNodes + iStrNode,
5290  pModeShapesr->GetVec(jMode*NFEMNodes + iNodeCurr));
5291  }
5292 
5293  /* nodo collegato 2 (e' il nodo multibody) */
5294  unsigned int uNode2 = (unsigned int)HP.GetInt();
5295  DEBUGCOUT("Linked to Multi-Body Node " << uNode2 << std::endl);
5296 
5297  /* verifica di esistenza del nodo 2 */
5298  SND[iStrNode - 1].pNode = pDM->pFindNode<const StructNode, Node::STRUCTURAL>(uNode2);
5299  if (SND[iStrNode - 1].pNode == 0) {
5300  silent_cerr("Modal(" << uLabel << "): "
5301  "StructuralNode(" << uNode2 << ") "
5302  "at line " << HP.GetLineData()
5303  << " not defined" << std::endl);
5304  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
5305  }
5306 
5307  /* offset del nodo Multi-Body */
5308  ReferenceFrame RF = ReferenceFrame(SND[iStrNode - 1].pNode);
5309  Vec3 d2(HP.GetPosRel(RF));
5310  Mat3x3 R2(::Eye3);
5311  if (HP.IsKeyWord("hinge") || HP.IsKeyWord("orientation")) {
5312  R2 = HP.GetRotRel(RF);
5313  }
5314 
5315  SND[iStrNode - 1].OffsetMB = d2;
5316  SND[iStrNode - 1].RotMB = R2; // FIXME: not used (so far)
5317 
5318  DEBUGCOUT("Multibody Node reference frame d2: " << std::endl
5319  << d2 << std::endl);
5320 
5321  /* salva le label dei nodi vincolati nell'array IntNodes;
5322  * puo' servire per il restart? */
5323  SND[iStrNode - 1].FEMNode = Node1;
5324 
5325  if (HP.IsKeyWord("output")) {
5326  SND[iStrNode - 1].bOut = HP.GetYesNoOrBool();
5327  }
5328 
5329  const Vec3& xMB(SND[iStrNode - 1].pNode->GetXCurr());
5330  pedantic_cout("Interface node " << iStrNode << ":" << std::endl
5331  << " MB node " << uNode2 << " x={" << xMB << "}" << std::endl);
5332  Vec3 xFEMRel(pXYZFEMNodes->GetVec(iNodeCurr));
5333  Vec3 xFEM(X0 + R*xFEMRel);
5334  pedantic_cout(" FEM node \"" << Node1 << "\" x={" << xFEM << "} "
5335  "xrel={" << xFEMRel << "}" << std::endl);
5336  pedantic_cout(" offset={" << xFEM - xMB << "}" << std::endl);
5337  } /* fine ciclo sui nodi d'interfaccia */
5338 
5339  if (bOut) {
5340  std::vector<Modal::StrNodeData>::iterator i = SND.begin();
5341  std::vector<Modal::StrNodeData>::iterator end = SND.end();
5342  for (; i != end; ++i) {
5343  i->bOut = bOut;
5344  }
5345  }
5346 
5347  /* fine ciclo caricamento dati */
5348 
5349  /*
5350  * calcola gli invarianti d'inerzia (massa, momenti statici
5351  * e d'inerzia, termini di accoppiamento nel SdR locale)
5352  */
5353 
5354  /* NOTE: right now, either they are all used (except for Inv9)
5355  * or none is used, and the global inertia of the body is expected */
5356  if (bBuildInvariants) {
5357  doublereal dMassInv = 0.;
5358  Vec3 STmpInv(::Zero3);
5359  Mat3x3 JTmpInv(::Zero3x3);
5360 
5361  MatNxN GenMass(NModes, 0.);
5362 
5363  /* TODO: select what terms are used */
5364 
5365  SAFENEWWITHCONSTRUCTOR(pInv3, Mat3xN, Mat3xN(NModes, 0.));
5366  SAFENEWWITHCONSTRUCTOR(pInv4, Mat3xN, Mat3xN(NModes, 0.));
5367 
5368  /* NOTE: we assume that Inv5 is used only if Inv4 is used as well */
5369  /* Inv5 e' un 3xMxM */
5370  SAFENEWWITHCONSTRUCTOR(pInv5, Mat3xN, Mat3xN(NModes*NModes, 0.));
5371 
5372  /* Inv8 e' un 3x3xM */
5373  SAFENEWWITHCONSTRUCTOR(pInv8, Mat3xN, Mat3xN(3*NModes, 0.));
5374 
5375  /* NOTE: we assume that Inv9 is used only if Inv8 is used as well */
5376  /* by default: no */
5377  if (HP.IsKeyWord("use" "invariant" "9")) {
5378  /* Inv9 e' un 3x3xMxM */
5379  SAFENEWWITHCONSTRUCTOR(pInv9, Mat3xN, Mat3xN(3*NModes*NModes, 0.));
5380  }
5381  /* Inv10 e' un 3x3xM */
5382  SAFENEWWITHCONSTRUCTOR(pInv10,Mat3xN, Mat3xN(3*NModes, 0.));
5383  SAFENEWWITHCONSTRUCTOR(pInv11,Mat3xN, Mat3xN(NModes, 0.));
5384 
5385  /* inizio ciclo scansione nodi */
5386  for (unsigned int iNode = 1; iNode <= NFEMNodes; iNode++) {
5387  doublereal mi = FEMMass[iNode - 1];
5388 
5389  /* massa totale (Inv 1) */
5390  dMassInv += mi;
5391 
5392  /* posizione nodi FEM */
5393  Vec3 ui = pXYZFEMNodes->GetVec(iNode);
5394 
5395  Mat3x3 uiWedge(MatCross, ui);
5396  Mat3x3 JiNodeTmp(::Zero3x3);
5397 
5398  JiNodeTmp(1, 1) = FEMJ[iNode - 1](1);
5399  JiNodeTmp(2, 2) = FEMJ[iNode - 1](2);
5400  JiNodeTmp(3, 3) = FEMJ[iNode - 1](3);
5401 
5402  JTmpInv += JiNodeTmp - Mat3x3(MatCrossCross, ui, ui*mi);
5403  STmpInv += ui*mi;
5404 
5405  /* estrae le forme modali del nodo i-esimo */
5406  for (unsigned int jMode = 1; jMode <= NModes; jMode++) {
5407  unsigned int iOffset = (jMode - 1)*NFEMNodes + iNode;
5408 
5409  PHIti.PutVec(jMode, pModeShapest->GetVec(iOffset));
5410  PHIri.PutVec(jMode, pModeShapesr->GetVec(iOffset));
5411  }
5412 
5413  /* TODO: only build what is required */
5414 
5415  Mat3xN Inv3Tmp(NModes, 0.);
5416  Mat3xN Inv4Tmp(NModes, 0.);
5417  Mat3xN Inv4JTmp(NModes, 0.);
5418  Inv3Tmp.Copy(PHIti);
5419 
5420  /* Inv3 = mi*PHIti, i = 1,...nnodi */
5421  Inv3Tmp *= mi;
5422 
5423  /* Inv4 = mi*ui/\*PHIti + Ji*PHIri, i = 1,...nnodi */
5424  Inv4Tmp.LeftMult(uiWedge*mi, PHIti);
5425  Inv4JTmp.LeftMult(JiNodeTmp, PHIri);
5426  Inv4Tmp += Inv4JTmp;
5427  *pInv3 += Inv3Tmp;
5428  *pInv4 += Inv4Tmp;
5429  *pInv11 += Inv4JTmp;
5430 
5431  /* inizio ciclo scansione modi */
5432  for (unsigned int jMode = 1; jMode <= NModes; jMode++) {
5433  Vec3 PHItij = PHIti.GetVec(jMode);
5434  Vec3 PHIrij = PHIri.GetVec(jMode);
5435 
5436  Mat3x3 PHItijvett_mi(MatCross, PHItij*mi);
5437  Mat3xN Inv5jTmp(NModes, 0);
5438 
5439  /* Inv5 = mi*PHItij/\*PHIti,
5440  * i = 1,...nnodi, j = 1,...nmodi */
5441  Inv5jTmp.LeftMult(PHItijvett_mi, PHIti);
5442  for (unsigned int kMode = 1; kMode <= NModes; kMode++) {
5443  pInv5->AddVec((jMode - 1)*NModes + kMode,
5444  Inv5jTmp.GetVec(kMode));
5445 
5446  /* compute the modal mass matrix
5447  * using the FEM inertia and the
5448  * mode shapes */
5449  GenMass(jMode, kMode) += (PHItij*PHIti.GetVec(kMode))*mi
5450  + PHIrij*(JiNodeTmp*PHIri.GetVec(kMode));
5451  }
5452 
5453  /* Inv8 = -mi*ui/\*PHItij/\,
5454  * i = 1,...nnodi, j = 1,...nmodi */
5455  Mat3x3 Inv8jTmp = -uiWedge*PHItijvett_mi;
5456  pInv8->AddMat3x3((jMode - 1)*3 + 1, Inv8jTmp);
5457 
5458  /* Inv9 = mi*PHItij/\*PHItik/\,
5459  * i = 1,...nnodi, j, k = 1...nmodi */
5460  if (pInv9 != 0) {
5461  for (unsigned int kMode = 1; kMode <= NModes; kMode++) {
5462  Mat3x3 PHItikvett(MatCross, PHIti.GetVec(kMode));
5463  Mat3x3 Inv9jkTmp = PHItijvett_mi*PHItikvett;
5464 
5465  pInv9->AddMat3x3((jMode - 1)*3*NModes + (kMode - 1)*3 + 1, Inv9jkTmp);
5466  }
5467  }
5468 
5469  /* Inv10 = [PHIrij/\][J0i],
5470  * i = 1,...nnodi, j = 1,...nmodi */
5471  Mat3x3 Inv10jTmp = PHIrij.Cross(JiNodeTmp);
5472  pInv10->AddMat3x3((jMode - 1)*3 + 1, Inv10jTmp);
5473  } /* fine ciclo scansione modi */
5474  } /* fine ciclo scansione nodi */
5475 
5476  if (bRecordGroup[12]) {
5477  Mat3x3 DJ = JTmp - JTmpInv;
5478  pedantic_cerr(" Rigid-body mass: "
5479  "input - computed" << std::endl
5480  << " " << dMass - dMassInv << std::endl);
5481  pedantic_cerr(" Rigid-body CM location: "
5482  "input - computed" << std::endl
5483  << " " << XTmpIn - STmpInv/dMassInv << std::endl);
5484  pedantic_cerr(" Rigid-body inertia: "
5485  "input - computed" << std::endl
5486  << " " << DJ.GetVec(1) << std::endl
5487  << " " << DJ.GetVec(2) << std::endl
5488  << " " << DJ.GetVec(3) << std::endl);
5489  }
5490 
5491  /*
5492  * TODO: check modal mass
5493  */
5494  pedantic_cerr(" Generalized mass: input - computed" << std:: endl);
5495  for (unsigned int jMode = 1; jMode <= NModes; jMode++) {
5496  pedantic_cerr(" ");
5497  for (unsigned int kMode = 1; kMode <= NModes; kMode++) {
5498  pedantic_cerr(" " << pGenMass->dGet(jMode, kMode) - GenMass(jMode, kMode));
5499  }
5500  pedantic_cerr(std::endl);
5501  }
5502 
5503  // if record 12 was read, leave its data in place,
5504  // otherwise replace rigid-body invariants with those from record 12
5505  if (!bRecordGroup[12]) {
5506  dMass = dMassInv;
5507  STmp = STmpInv;
5508  JTmp = JTmpInv;
5509  }
5510  }
5511 
5512  // if record 12 was read, fix static moment
5513  if (bRecordGroup[12]) {
5514  /* left over when reading XCG */
5515  STmp = XTmpIn*dMass;
5516  }
5517 
5518  /*
5519  * TODO: Check rank of modal stiffness matrix
5520  */
5521 
5522  // check if mass matrix is symmetric and diagonal
5523  bool bIsMSym = true;
5524  bool bIsMDiag = true;
5525  for (unsigned iRow = 2; iRow <= NModes; iRow++) {
5526  for (unsigned iCol = 1; iCol < iRow; iCol++) {
5527  doublereal mrc = pGenMass->dGet(iRow, iCol);
5528  doublereal mcr = pGenMass->dGet(iCol, iRow);
5529  if (mrc != mcr) {
5530  if (bIsMSym) {
5531  silent_cerr("Modal(" << uLabel << "): mass matrix is not symmetric: (at least) "
5532  "M(" << iRow << ", " << iCol << ")=" << mrc << " "
5533  "!= "
5534  "M(" << iCol << ", " << iRow << ")=" << mcr << " "
5535  << std::endl);
5536  }
5537  bIsMSym = false;
5538  }
5539 
5540  if (mrc != 0. || mcr != 0.) {
5541  if (bIsMDiag) {
5542  silent_cerr("Modal(" << uLabel << "): mass matrix is not diagonal: (at least) "
5543  "M(" << iRow << ", " << iCol << ")=" << mrc << " "
5544  "and/or "
5545  "M(" << iCol << ", " << iRow << ")=" << mcr << " "
5546  "!= 0.0"
5547  << std::endl);
5548  }
5549  bIsMDiag = false;
5550  }
5551  }
5552  }
5553 
5554  bool bIsKSym = true;
5555  bool bIsKDiag = true;
5556  for (unsigned iRow = 2; iRow <= NModes; iRow++) {
5557  for (unsigned iCol = 1; iCol < iRow; iCol++) {
5558  doublereal mrc = pGenStiff->dGet(iRow, iCol);
5559  doublereal mcr = pGenStiff->dGet(iCol, iRow);
5560  if (mrc != mcr) {
5561  if (bIsKSym) {
5562  silent_cerr("Modal(" << uLabel << "): stiffness matrix is not symmetric: (at least) "
5563  "K(" << iRow << ", " << iCol << ")=" << mrc << " "
5564  "!= "
5565  "K(" << iCol << ", " << iRow << ")=" << mcr << " "
5566  << std::endl);
5567  }
5568  bIsKSym = false;
5569  }
5570 
5571  if (mrc != 0. || mcr != 0.) {
5572  if (bIsKDiag) {
5573  silent_cerr("Modal(" << uLabel << "): stiffness matrix is not diagonal: (at least) "
5574  "K(" << iRow << ", " << iCol << ")=" << mrc << " "
5575  "and/or "
5576  "K(" << iCol << ", " << iRow << ")=" << mcr << " "
5577  "!= 0.0"
5578  << std::endl);
5579  }
5580  bIsKDiag = false;
5581  }
5582  }
5583  }
5584 
5585  if (eDamp == DAMPING_FROM_FILE) {
5586  bool bIsCSym = true;
5587  bool bIsCDiag = true;
5588  for (unsigned iRow = 2; iRow <= NModes; iRow++) {
5589  for (unsigned iCol = 1; iCol < iRow; iCol++) {
5590  doublereal mrc = pGenDamp->dGet(iRow, iCol);
5591  doublereal mcr = pGenDamp->dGet(iCol, iRow);
5592  if (mrc != mcr) {
5593  if (bIsCSym) {
5594  silent_cerr("Modal(" << uLabel << "): damping matrix is not symmetric: (at least) "
5595  "C(" << iRow << ", " << iCol << ") "
5596  "!= "
5597  "C(" << iCol << ", " << iRow << ") "
5598  << std::endl);
5599  }
5600  bIsCSym = false;
5601  }
5602 
5603  if (mrc != 0. || mcr != 0.) {
5604  if (bIsCDiag) {
5605  silent_cerr("Modal(" << uLabel << "): damping matrix is not diagonal: (at least) "
5606  "C(" << iRow << ", " << iCol << ") "
5607  "and/or "
5608  "C(" << iCol << ", " << iRow << ") "
5609  "!= 0.0"
5610  << std::endl);
5611  }
5612  bIsCDiag = false;
5613  }
5614  }
5615  }
5616 
5617  } else {
5618  /*
5619  * costruisce la matrice di smorzamento:
5620  * il termine diagonale i-esimo e' pari a
5621  * cii = 2*cdampi*(ki*mi)^.5
5622  */
5623 
5624  switch (eDamp) {
5625  case DAMPING_DIAG:
5626  case DAMPING_SINGLE_FACTOR:
5627  if (!bIsMDiag || !bIsKDiag) {
5628  silent_cerr("Modal(" << uLabel << "): "
5629  "warning, " << sDamp[eDamp]
5630  << " with non-diagonal mass and/or stiffness matrix" << std::endl);
5631  }
5632 
5633  for (unsigned int iCnt = 1; iCnt <= NModes; iCnt++) {
5634  doublereal k = pGenStiff->dGet(iCnt, iCnt);
5635  doublereal m = pGenMass->dGet(iCnt, iCnt);
5636  doublereal d = sqrt(k*m);
5637 
5638  if (eDamp == DAMPING_DIAG) {
5639  pGenDamp->Put(iCnt, iCnt, 2.*DampRatios[iCnt - 1]*d);
5640 
5641  } else if (eDamp == DAMPING_SINGLE_FACTOR) {
5642  pGenDamp->Put(iCnt, iCnt, 2.*damp_factor*d);
5643  }
5644  }
5645  break;
5646 
5647  case DAMPING_RAYLEIGH:
5648  for (unsigned int iRow = 1; iRow <= NModes; iRow++) {
5649  for (unsigned int iCol = 1; iCol <= NModes; iCol++) {
5650  doublereal m = pGenMass->dGet(iRow, iCol);
5651  doublereal k = pGenStiff->dGet(iRow, iCol);
5652  pGenDamp->Put(iRow, iCol, damp_coef_M*m + damp_coef_K*k);
5653  }
5654  }
5655  break;
5656 
5657  default:
5658  // nothing to do
5659  break;
5660  }
5661  }
5662 
5663 #if 1 // def DEBUG
5664  if (pedantic_output) {
5665  std::ostream &out = std::cout;
5666 
5667  out << " Total Mass: " << dMass << std::endl;
5668  out << " Inertia Matrix (referred to modal node): " << std::endl
5669  << " " << JTmp << std::endl;
5670  out << " Static Moment Vector: " << STmp << std::endl;
5671 
5672  out << " Generalized Stiffness: " << std::endl;
5673  for (unsigned int iCnt = 1; iCnt <= NModes; iCnt++) {
5674  out << " ";
5675  for (unsigned int jCnt = 1; jCnt <= NModes; jCnt++) {
5676  out << " " << pGenStiff->dGet(iCnt, jCnt);
5677  }
5678  out << std::endl;
5679  }
5680 
5681  out << " Generalized Mass: " << std::endl;
5682  for (unsigned int iCnt = 1; iCnt <= NModes; iCnt++) {
5683  out << " ";
5684  for (unsigned int jCnt = 1; jCnt <= NModes; jCnt++) {
5685  out << " " << pGenMass->dGet(iCnt, jCnt);
5686  }
5687  out << std::endl;
5688  }
5689 
5690  out << " Generalized Damping: " << std::endl;
5691  for (unsigned int iCnt = 1; iCnt <= NModes; iCnt++) {
5692  out << " ";
5693  for (unsigned int jCnt = 1; jCnt <= NModes; jCnt++) {
5694  out << " " << pGenDamp->dGet(iCnt, jCnt);
5695  }
5696  out << std::endl;
5697  }
5698 
5699  if (pInv3 != 0) {
5700  out << " Inv3: " << std::endl;
5701  for (unsigned int iCnt = 1; iCnt <= 3; iCnt++) {
5702  out << " ";
5703  for (unsigned int jCnt = 1; jCnt <= NModes; jCnt++) {
5704  out << " " << pInv3->dGet(iCnt, jCnt);
5705  }
5706  out << std::endl;
5707  }
5708  } else {
5709  out << " Inv3: unused" << std::endl;
5710  }
5711 
5712  if (pInv4 != 0) {
5713  out << " Inv4: " << std::endl;
5714  for (unsigned int iCnt = 1; iCnt <= 3; iCnt++) {
5715  out << " ";
5716  for (unsigned int jCnt = 1; jCnt <= NModes; jCnt++) {
5717  out << " " << pInv4->dGet(iCnt, jCnt);
5718  }
5719  out << std::endl;
5720  }
5721  } else {
5722  out << " Inv4: unused" << std::endl;
5723  }
5724 
5725  if (pInv5 != 0) {
5726  for (unsigned int jMode = 1; jMode <= NModes; jMode++) {
5727  out << " Inv5j(j=" << jMode << "): " << std::endl;
5728  for (unsigned int iCnt = 1; iCnt <= 3; iCnt++) {
5729  out << " ";
5730  for (unsigned int kMode = 1; kMode <= NModes; kMode++) {
5731  out << " " << pInv5->dGet(iCnt, (jMode - 1)*NModes + kMode);
5732  }
5733  out << std::endl;
5734  }
5735  }
5736  } else {
5737  out << " Inv5: unused" << std::endl;
5738  }
5739 
5740  if (pInv8 != 0) {
5741  for (unsigned int jMode = 1; jMode <= NModes; jMode++) {
5742  out << " Inv8j(j=" << jMode << "): " << std::endl;
5743  for (unsigned int iCnt = 1; iCnt <= 3; iCnt++) {
5744  out << " ";
5745  for (unsigned int jCnt = 1; jCnt <= 3; jCnt++) {
5746  out << " " << pInv8->dGet(iCnt, (jMode - 1)*3 + jCnt);
5747  }
5748  out << std::endl;
5749  }
5750  }
5751  } else {
5752  out << " Inv8: unused" << std::endl;
5753  }
5754 
5755  if (pInv9 != 0) {
5756  for (unsigned int jMode = 1; jMode <= NModes; jMode++) {
5757  for (unsigned int kMode = 1; kMode <= NModes; kMode++) {
5758  out << " Inv9jk(j=" << jMode << ",k=" << kMode << "): " << std::endl;
5759  for (unsigned int iCnt = 1; iCnt <= 3; iCnt++) {
5760  out << " ";
5761  for (unsigned int jCnt = 1; jCnt <= 3; jCnt++) {
5762  out << " " << pInv9->dGet(iCnt, (jMode - 1)*3*NModes + (kMode - 1)*3 + jCnt);
5763  }
5764  out << std::endl;
5765  }
5766  }
5767  }
5768  } else {
5769  out << " Inv9: unused" << std::endl;
5770  }
5771  }
5772 #endif /* DEBUG */
5773 
5774  if (HP.IsStringWithDelims()) {
5775  const char *sTmp = HP.GetFileName();
5776  silent_cout("Modal(" << uLabel << "): warning, the syntax changed "
5777  "since 1.2.7; the output now occurs to a common \".mod\" file, "
5778  "the per-element file \"" << sTmp << "\" is no longer required, "
5779  "and will actually be ignored." << std::endl);
5780  }
5781  flag fOut = pDM->fReadOutput(HP, Elem::JOINT);
5782 
5783  if (bInitialValues) {
5784  for (unsigned int iCnt = 0; iCnt < NModes; iCnt++) {
5785  a->Put(iCnt + 1, InitialValues[iCnt]);
5786  aP->Put(iCnt + 1, InitialValuesP[iCnt]);
5787  }
5788  }
5789 
5790  SAFENEWWITHCONSTRUCTOR(pEl,
5791  Modal,
5792  Modal(uLabel,
5793  pModalNode,
5794  X0,
5795  R,
5796  pDO,
5797  NModes,
5798  NStrNodes,
5799  NFEMNodes,
5800  dMass,
5801  STmp,
5802  JTmp,
5803  uModeNumber,
5804  pGenMass,
5805  pGenStiff,
5806  pGenDamp,
5807  IdFEMNodes,
5808  pXYZFEMNodes,
5809  SND,
5810  pPHItStrNode,
5811  pPHIrStrNode,
5812  pModeShapest,
5813  pModeShapesr,
5814  pInv3,
5815  pInv4,
5816  pInv5,
5817  pInv8,
5818  pInv9,
5819  pInv10,
5820  pInv11,
5821  a,
5822  aP,
5823  fOut));
5824 
5825  if (fOut) {
5826  pDM->OutputOpen(OutputHandler::MODAL);
5827  }
5828 
5829  SAFEDELETE(a);
5830  SAFEDELETE(aP);
5831 
5832  std::ostream& os = pDM->GetLogFile();
5833 
5834  // If we do not cast the label to integer, the output will be UINT_MAX if no modal node is defined
5835  os << "modal: " << pEl->GetLabel() << " "
5836  << static_cast<integer>(pModalNode ? pModalNode->GetLabel() : -1) << " "
5837  << dMass << " "
5838  << STmp / dMass << " "
5839  << JTmp << std::endl;
5840 
5841  return pEl;
5842 }
DataManager::fReadOutput
flag fReadOutput(MBDynParser &HP, const T &t) const
Definition: dataman.h:1064
MBDynParser::GetRotRel
Mat3x3 GetRotRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1795
Zero3
const Vec3 Zero3(0., 0., 0.)
Vec3::Cross
Vec3 Cross(const Vec3 &v) const
Definition: matvec3.h:218
MatNxN::Put
void Put(integer i, integer j, const doublereal &d)
Definition: matvec3n.h:578
flag
long int flag
Definition: mbdyn.h:43
Mat3xN::Put
void Put(int i, integer j, const doublereal &d)
Definition: matvec3n.h:306
MBDYN_EXCEPT_ARGS
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Vec3
Definition: matvec3.h:98
HighParser::GetInt
virtual integer GetInt(integer iDefval=0)
Definition: parser.cc:1050
MatCross
const MatCross_Manip MatCross
Definition: matvec3.cc:639
StructNode::GetRCurr
virtual const Mat3x3 & GetRCurr(void) const
Definition: strnode.h:1012
MBDynParser::GetRotAbs
Mat3x3 GetRotAbs(const ReferenceFrame &rf)
Definition: mbpar.cc:1857
Eye3
const Mat3x3 Eye3(1., 0., 0., 0., 1., 0., 0., 0., 1.)
IncludeParser::GetFileName
virtual const char * GetFileName(enum Delims Del=DEFAULTDELIM)
Definition: parsinc.cc:673
MatNxN::dGet
const doublereal & dGet(integer i, integer j) const
Definition: matvec3n.h:612
grad::fabs
GradientExpression< UnaryExpr< FuncFabs, Expr > > fabs(const GradientExpression< Expr > &u)
Definition: gradient.h:2973
Mat3x3::GetVec
Vec3 GetVec(unsigned short int i) const
Definition: matvec3.h:893
HighParser::GetYesNoOrBool
virtual bool GetYesNoOrBool(bool bDefval=false)
Definition: parser.cc:1038
Modal
Definition: modal.h:74
Zero3x3
const Mat3x3 Zero3x3(0., 0., 0., 0., 0., 0., 0., 0., 0.)
Mat3xN::GetVec
Vec3 GetVec(integer iCol) const
Definition: matvec3n.cc:553
AbsRefFrame
const ReferenceFrame AbsRefFrame(0, Vec3(0., 0., 0), Mat3x3(1., 0., 0., 0., 1., 0., 0., 0., 1.), Vec3(0., 0., 0), Vec3(0., 0., 0), EULER_123)
VecN::Put
void Put(integer i, const doublereal &d)
Definition: matvec3n.h:159
StructNode::GetStructNodeType
virtual StructNode::Type GetStructNodeType(void) const =0
MBDynParser::GetPosRel
Vec3 GetPosRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1331
HighParser::IsKeyWord
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
VecN::iGetNumRows
integer iGetNumRows(void) const
Definition: matvec3n.h:150
DEBUGCOUT
#define DEBUGCOUT(msg)
Definition: myassert.h:232
MBDynParser::GetPosAbs
Vec3 GetPosAbs(const ReferenceFrame &rf)
Definition: mbpar.cc:1401
HighParser::GetStringWithDelims
virtual const char * GetStringWithDelims(enum Delims Del=DEFAULTDELIM, bool escape=true)
Definition: parser.cc:1228
Mat3xN::SubVec
void SubVec(integer iCol, const Vec3 &v)
Definition: matvec3n.cc:593
Mat3xN::iGetNumCols
integer iGetNumCols(void) const
Definition: matvec3n.h:298
ASSERT
#define ASSERT(expression)
Definition: colamd.c:977
StructDispNode::pGetRBK
const RigidBodyKinematics * pGetRBK(void) const
Definition: strnode.cc:152
grad::sqrt
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
SAFENEWWITHCONSTRUCTOR
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
Mat3xN::dGet
const doublereal & dGet(int i, integer j) const
Definition: matvec3n.h:336
StructDispNode::GetXCurr
virtual const Vec3 & GetXCurr(void) const
Definition: strnode.h:310
Mat3xN::AddVec
void AddVec(integer iCol, const Vec3 &v)
Definition: matvec3n.cc:579
c
static std::stack< cleanup * > c
Definition: cleanup.cc:59
HighParser::GetString
virtual std::string GetString(const std::string &sDefVal)
Definition: parser.cc:1074
MatCrossCross
const MatCrossCross_Manip MatCrossCross
Definition: matvec3.cc:640
STRLENOF
#define STRLENOF(s)
Definition: mbdyn.h:166
HighParser::IsStringWithDelims
virtual bool IsStringWithDelims(enum Delims Del=DEFAULTDELIM)
Definition: parser.cc:1210
DataManager::GetLogFile
std::ostream & GetLogFile(void) const
Definition: dataman.h:326
Mat3xN::PutVec
void PutVec(integer iCol, const Vec3 &v)
Definition: matvec3n.cc:565
VecN
Definition: matvec3n.h:76
Joint
Definition: joint.h:50
a
static const doublereal a
Definition: hfluid_.h:289
DataManager::OutputOpen
void OutputOpen(const OutputHandler::OutFiles out)
Definition: dataman.cc:677
Mat3xN::AddMat3x3
void AddMat3x3(integer iCol, const Mat3x3 &m)
Definition: matvec3n.cc:652
integer
long int integer
Definition: colamd.c:51
IncludeParser::GetLineData
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
WithLabel::GetLabel
unsigned int GetLabel(void) const
Definition: withlab.cc:62
DataManager::pFindNode
Node * pFindNode(Node::Type Typ, unsigned int uNode) const
Definition: nodeman.cc:179
R
Mat3x3 R
Definition: mbdyn_util_oct.cc:74
SAFEDELETE
#define SAFEDELETE(pnt)
Definition: mynewmem.h:710
HighParser::GetReal
virtual doublereal GetReal(const doublereal &dDefval=0.0)
Definition: parser.cc:1056

Here is the call graph for this function:

Variable Documentation

const char* mdof[]
static
Initial value:
= {
"modal displacement q",
"modal velocity qP"
}

Definition at line 317 of file modal.cc.

Referenced by Modal::DescribeDof().

const char* meq[]
static
Initial value:
= {
"modal velocity definition q",
"modal equilibrium equation f"
}

Definition at line 327 of file modal.cc.

Referenced by Modal::DescribeEq().

const char* rdof[]
static
Initial value:
= {
"reaction force F",
"reaction couple M",
"reaction force derivative FP",
"reaction couple derivative MP"
}

Definition at line 321 of file modal.cc.

Referenced by Modal::DescribeDof().

const char* req[]
static
Initial value:
= {
"position constraint P",
"orientation constraint g",
"velocity constraint v",
"angular velocity constraint w"
}

Definition at line 331 of file modal.cc.

Referenced by Modal::DescribeEq().

const char xyz[] = "xyz"
static

Definition at line 316 of file modal.cc.

Referenced by Modal::DescribeDof(), and Modal::DescribeEq().