MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
dataman2.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/dataman2.cc,v 1.194 2017/05/12 17:28:14 morandini 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 /* DataManager -
33  * continua qui perche' il file dataman.cc sta diventando lungo */
34 
35 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
36 #include <algorithm>
37 #include <set>
38 #include <cmath>
39 #include <sstream>
40 // #include <typeinfo>
41 
42 #include "dataman.h"
43 #include "dataman_.h"
44 
45 #include "gravity.h"
46 #include "aerodyn.h"
47 #include "solver.h"
48 #include "ls.h"
49 
50 #include "Rot.hh"
51 #include "naivemh.h"
52 #include "spmapmh.h"
53 
54 #include "bufferstream_out_elem.h"
55 #include "bufferstreamdrive.h"
56 
57 const LoadableCalls *
58 DataManager::GetLoadableElemModule(std::string name) const
59 {
60  for (int j = 0; name[j]; j++) {
61  name[j] = tolower(name[j]);
62  }
63 
64  typedef std::map<std::string,const LoadableCalls *> mleh;
65  mleh::const_iterator i = MapOfLoadableElemHandlers.find(name);
66  if (i == MapOfLoadableElemHandlers.end()) {
67  return 0;
68  }
69  return i->second;
70 }
71 
72 void
73 DataManager::SetLoadableElemModule(std::string name,
74  const LoadableCalls *calls, ModuleInsertMode mode)
75 {
76  for (int j = 0; name[j]; j++) {
77  name[j] = tolower(name[j]);
78  }
79 
80  const LoadableCalls *tmp = GetLoadableElemModule(name);
81 
82  if (tmp != 0) {
83  switch (mode) {
84  case MIM_FAIL:
85  default:
86  silent_cerr("DataManager::SetLoadableElemModule(): "
87  "loadable element handler \"" << name
88  << "\" already defined" << std::endl);
89  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
90 
91  case MIM_IGNORE:
92  silent_cout("DataManager::SetLoadableElemModule(): "
93  "loadable element handler \"" << name
94  << "\" already defined; "
95  "new definition ignored" << std::endl);
96  return;
97 
98  case MIM_REPLACE:
99  silent_cout("DataManager::SetLoadableElemModule(): "
100  "loadable element handler \"" << name
101  << "\" already defined; "
102  "replaced by new definition" << std::endl);
103  break;
104  }
105  }
106 
107  MapOfLoadableElemHandlers[name] = calls;
108 }
109 
110 const doublereal&
111 DataManager::dGetInitialPositionStiffness(void) const
112 {
113  return dInitialPositionStiffness;
114 }
115 
116 const doublereal&
117 DataManager::dGetInitialVelocityStiffness(void) const
118 {
119  return dInitialVelocityStiffness;
120 }
121 
122 bool
123 DataManager::bDoesOmegaRotate(void) const
124 {
125  return bOmegaRotates;
126 }
127 
128 void
129 DataManager::IncElemCount(Elem::Type type)
130 {
131  /* FIXME: assert the data structure has not been allocated yet */
132  ElemData[type].iExpectedNum++;
133 }
134 
135 /* Setta il valore della variabile Time nel DataManager
136  * usato dal metodo numerico all'inizio di ogni step temporale */
137 
138 void
139 DataManager::SetTime(const doublereal& dTime, const doublereal& dTimeStep,
140  const integer& iStep, bool bServePending)
141 {
142  /* Setta il tempo nel DriveHandler */
143  DrvHdl.SetTime(dTime, dTimeStep, iStep);
144 
145  DEBUGLCOUT(MYDEBUG_INIT|MYDEBUG_ASSEMBLY,
146  "Global symbol table:" << std::endl
147  << MathPar.GetSymbolTable() << std::endl);
148 
149  /* serve i drive pending */
150  if (bServePending) {
151  for (int iType = 0; iType < Drive::LASTDRIVETYPE; iType++) {
152  for (unsigned int iCnt = 0; iCnt < DriveData[iType].iNum; iCnt++) {
153  DriveData[iType].ppFirstDrive[iCnt]->ServePending(dTime);
154  }
155  }
156  }
157 
158  // updates rigid body kinematics, if any
159  if (pRBK) {
160  pRBK->Update();
161  }
162 } /* End of DataManager::SetTime() */
163 
164 doublereal
165 DataManager::dGetTime(void) const
166 {
167  return DrvHdl.dGetTime();
168 } /* End of DataManager::dGetTime() */
169 
170 /* Collega il DataManager ed il DriveHandler alla soluzione */
171 void
172 DataManager::LinkToSolution(VectorHandler& XCurr,
173  VectorHandler& XPrimeCurr)
174 {
175  pXCurr = &XCurr;
176  pXPrimeCurr = &XPrimeCurr;
177  DrvHdl.LinkToSolution(XCurr, XPrimeCurr);
178 }
179 
180 /* Inizializzatore dei dof di ogni elemento */
181 void
182 DataManager::DofOwnerInit(void)
183 {
184  DEBUGCOUTFNAME("DataManager::DofOwnerInit");
185  ASSERT(!Dofs.empty());
186  ASSERT(!Nodes.empty());
187 
188  if ( uPrintFlags & PRINT_TO_FILE ) {
189  OutHdl.Open(OutputHandler::DOFSTATS);
190  }
191 
192  std::ostream& out_ds = (uPrintFlags & PRINT_TO_FILE)
193  ? OutHdl.DofStats()
194  : std::cout;
195 
196  bool pds =
197 #ifdef DEBUG
198  DEBUG_LEVEL_MATCH(MYDEBUG_INIT|MYDEBUG_ASSEMBLY) ||
199 #endif /* DEBUG */
200  (!silent_output && (uPrintFlags & PRINT_DOF_STATS));
201 
202  /* NOTE: further direct use of std::cout instead
203  * of silent_cout() macro because silent_cout is
204  * tested in "pds".
205  */
206  if (pds) {
207  out_ds << "Regular steps dof stats" << std::endl;
208  }
209 
210  /* per ogni nodo strutturale */
211  if (!NodeData[Node::STRUCTURAL].NodeContainer.empty()) {
212 
213  /*
214  * used by POD stuff: if any, output
215  * the list of the first dof (minus 1)
216  * of each structural node, so it's easy
217  * to get the struct node values
218  * in MATLAB: given a vector "X" with all
219  * the states, and a vector
220  * "v" with the first dof of each
221  * structural node, then the x coordinate
222  * is X(v+1) and so forth
223  */
224 
225  OutHdl.Log() << "struct node dofs:";
226  for (NodeContainerType::const_iterator i = NodeData[Node::STRUCTURAL].NodeContainer.begin();
227  i != NodeData[Node::STRUCTURAL].NodeContainer.end(); ++i)
228  {
229  const StructNode *pNode = dynamic_cast<const StructNode *>(i->second);
230  if (pNode) {
231  if (pNode->GetStructNodeType() == StructNode::DUMMY) {
232  continue;
233  }
234  OutHdl.Log() << " " << pNode->iGetFirstPositionIndex();
235  }
236  }
237  OutHdl.Log() << std::endl;
238 
239  OutHdl.Log() << "struct node eqs:";
240  for (NodeContainerType::const_iterator i = NodeData[Node::STRUCTURAL].NodeContainer.begin();
241  i != NodeData[Node::STRUCTURAL].NodeContainer.end(); ++i)
242  {
243  const StructNode *pNode = dynamic_cast<const StructNode *>(i->second);
244  if (pNode) {
245  if (pNode->GetStructNodeType() == StructNode::DUMMY) {
246  continue;
247  }
248  OutHdl.Log() << " " << pNode->iGetFirstMomentumIndex();
249  }
250  }
251  OutHdl.Log() << std::endl;
252 
253  OutHdl.Log() << "struct node momentum dofs:";
254  for (NodeContainerType::const_iterator i = NodeData[Node::STRUCTURAL].NodeContainer.begin();
255  i != NodeData[Node::STRUCTURAL].NodeContainer.end(); ++i)
256  {
257  const StructNode *pNode = dynamic_cast<const StructNode *>(i->second);
258  if (pNode) {
259  switch (pNode->GetStructNodeType()) {
260  case StructNode::STATIC:
261  case StructNode::DUMMY:
262  continue;
263 
264  default:
265  break;
266  }
267  OutHdl.Log() << " " << pNode->iGetFirstMomentumIndex();
268  }
269  }
270  OutHdl.Log() << std::endl;
271 
272  OutHdl.Log() << "struct node momentum eqs:";
273  for (NodeContainerType::const_iterator i = NodeData[Node::STRUCTURAL].NodeContainer.begin();
274  i != NodeData[Node::STRUCTURAL].NodeContainer.end(); ++i)
275  {
276  const StructNode *pNode = dynamic_cast<const StructNode *>(i->second);
277  if (pNode) {
278  switch (pNode->GetStructNodeType()) {
279  case StructNode::STATIC:
280  case StructNode::DUMMY:
281  continue;
282 
283  default:
284  break;
285  }
286  OutHdl.Log() << " " << pNode->iGetFirstIndex();
287  }
288  }
289  OutHdl.Log() << std::endl;
290 
291  OutHdl.Log() << "struct displacement node dofs:";
292  for (NodeContainerType::const_iterator i = NodeData[Node::STRUCTURAL].NodeContainer.begin();
293  i != NodeData[Node::STRUCTURAL].NodeContainer.end(); ++i)
294  {
295  const StructDispNode *pDispNode = dynamic_cast<const StructDispNode *>(i->second);
296  if (pDispNode) {
297  const StructNode* pNode = dynamic_cast<const StructNode*>(pDispNode);
298  if (pNode && pNode->GetStructNodeType() == StructNode::DUMMY) {
299  continue;
300  }
301  OutHdl.Log() << " " << pDispNode->iGetFirstPositionIndex();
302  }
303  }
304  OutHdl.Log() << std::endl;
305 
306  OutHdl.Log() << "struct displacement node eqs:";
307  for (NodeContainerType::const_iterator i = NodeData[Node::STRUCTURAL].NodeContainer.begin();
308  i != NodeData[Node::STRUCTURAL].NodeContainer.end(); ++i)
309  {
310  const StructDispNode *pDispNode = dynamic_cast<const StructDispNode *>(i->second);
311  if (pDispNode) {
312  const StructNode* pNode = dynamic_cast<const StructNode*>(pDispNode);
313  if (pNode && pNode->GetStructNodeType() == StructNode::DUMMY) {
314  continue;
315  }
316  OutHdl.Log() << " " << pDispNode->iGetFirstMomentumIndex();
317  }
318  }
319  OutHdl.Log() << std::endl;
320 
321  OutHdl.Log() << "struct displacement node momentum dofs:";
322  for (NodeContainerType::const_iterator i = NodeData[Node::STRUCTURAL].NodeContainer.begin();
323  i != NodeData[Node::STRUCTURAL].NodeContainer.end(); ++i)
324  {
325  const StructDispNode *pDispNode = dynamic_cast<const StructDispNode *>(i->second);
326  if (pDispNode) {
327  const StructNode* pNode = dynamic_cast<const StructNode*>(pDispNode);
328  if (pNode && pNode->GetStructNodeType() == StructNode::DUMMY) {
329  continue;
330  }
331  switch (pDispNode->GetStructDispNodeType()) {
332  case StructDispNode::STATIC:
333  continue;
334 
335  default:
336  break;
337  }
338  OutHdl.Log() << " " << pDispNode->iGetFirstMomentumIndex();
339  }
340  }
341  OutHdl.Log() << std::endl;
342 
343  OutHdl.Log() << "struct displacement node momentum eqs:";
344  for (NodeContainerType::const_iterator i = NodeData[Node::STRUCTURAL].NodeContainer.begin();
345  i != NodeData[Node::STRUCTURAL].NodeContainer.end(); ++i)
346  {
347  const StructDispNode *pDispNode = dynamic_cast<const StructDispNode *>(i->second);
348  if (pDispNode) {
349  const StructNode* pNode = dynamic_cast<const StructNode*>(pDispNode);
350  if (pNode && pNode->GetStructNodeType() == StructNode::DUMMY) {
351  continue;
352  }
353  switch (pDispNode->GetStructDispNodeType()) {
354  case StructNode::STATIC:
355  continue;
356 
357  default:
358  break;
359  }
360  OutHdl.Log() << " " << pDispNode->iGetFirstIndex();
361  }
362  }
363  OutHdl.Log() << std::endl;
364 
365  OutHdl.Log() << "struct node labels:";
366  for (NodeContainerType::const_iterator i = NodeData[Node::STRUCTURAL].NodeContainer.begin();
367  i != NodeData[Node::STRUCTURAL].NodeContainer.end(); ++i)
368  {
369  const StructNode *pNode = dynamic_cast<const StructNode *>(i->second);
370  if (pNode) {
371  if (pNode->GetStructNodeType() == StructNode::DUMMY) {
372  continue;
373  }
374  OutHdl.Log() << " " << pNode->GetLabel();
375  }
376  }
377  OutHdl.Log() << std::endl;
378 
379  OutHdl.Log() << "struct displacement node labels:";
380  for (NodeContainerType::const_iterator i = NodeData[Node::STRUCTURAL].NodeContainer.begin();
381  i != NodeData[Node::STRUCTURAL].NodeContainer.end(); ++i)
382  {
383  const StructDispNode *pDispNode = dynamic_cast<const StructDispNode *>(i->second);
384  if (pDispNode) {
385  const StructNode* pNode = dynamic_cast<const StructNode*>(pDispNode);
386  if (pNode && pNode->GetStructNodeType() == StructNode::DUMMY) {
387  continue;
388  }
389  OutHdl.Log() << " " << pDispNode->GetLabel();
390  }
391  }
392  OutHdl.Log() << std::endl;
393  }
394 
395  /* per ogni nodo */
396  for (NodeVecType::const_iterator i = Nodes.begin(); i != Nodes.end(); ++i) {
397  DEBUGLCOUT(MYDEBUG_INIT|MYDEBUG_ASSEMBLY,
398  psNodeNames[(*i)->GetNodeType()]
399  << "(" << (*i)->GetLabel() << ")"
400  << std::endl);
401 
402  /* chiede al nodo quanti dof possiede */
403  unsigned int iNumDof = (*i)->iGetNumDof();
404  if (iNumDof > 0) {
405  ASSERT((*i)->iGetFirstIndex() >= 0);
406 
407  /* si fa passare il primo Dof */
408  Dof* pDf = &Dofs[(*i)->iGetFirstIndex()];
409 
410 #ifdef DEBUG
411  DEBUGLCOUT(MYDEBUG_INIT|MYDEBUG_ASSEMBLY,
412  psNodeNames[(*i)->GetNodeType()]
413  << "(" << (*i)->GetLabel() << "): "
414  "first dof = " << pDf->iIndex + 1
415  << std::endl);
416 #endif /* DEBUG */
417 
418  if (pds) {
419  unsigned int nd = (*i)->iGetNumDof();
420  integer fd = pDf->iIndex;
421 
422  out_ds << psNodeNames[(*i)->GetNodeType()]
423  << "(" << (*i)->GetLabel() << "): "
424  << nd << " " << fd + 1;
425  if (nd > 1) {
426  out_ds << "->" << fd + nd;
427  }
428  out_ds << std::endl;
429  if (uPrintFlags & PRINT_DOF_DESCRIPTION) {
430  (*i)->DescribeDof(out_ds,
431  " ");
432  }
433 
434  if (uPrintFlags & PRINT_EQ_DESCRIPTION) {
435  (*i)->DescribeEq(out_ds,
436  " ");
437  }
438  }
439 
440  /* per ogni Dof, chiede al nodo di che tipo e' e lo
441  * setta nel DofOwner */
442  std::vector<std::string> DofDesc;
443  (*i)->DescribeDof(DofDesc);
444  if (DofDesc.size() == iNumDof) {
445  for (unsigned int iCnt = 0; iCnt < iNumDof; iCnt++) {
446  pDf[iCnt].Description = DofDesc[iCnt];
447  }
448 
449  } else {
450  std::ostringstream os;
451  os << psNodeNames[(*i)->GetNodeType()]
452  << "(" << (*i)->GetLabel() << ")";
453  std::string name(os.str());
454 
455  for (unsigned int iCnt = 0; iCnt < iNumDof; iCnt++) {
456  os.str(name);
457  os.seekp(0, std::ios_base::end);
458  os << ": dof(" << iCnt + 1 << ")";
459  pDf[iCnt].Description = os.str();
460  }
461  }
462 
463  std::vector<std::string> EqDesc;
464  (*i)->DescribeEq(EqDesc);
465  if (EqDesc.size() == iNumDof) {
466  for (unsigned int iCnt = 0; iCnt < iNumDof; iCnt++) {
467  pDf[iCnt].EqDescription = EqDesc[iCnt];
468  }
469 
470  } else {
471  std::ostringstream os;
472  os << psNodeNames[(*i)->GetNodeType()]
473  << "(" << (*i)->GetLabel() << ")";
474  std::string name(os.str());
475 
476  for (unsigned int iCnt = 0; iCnt < iNumDof; iCnt++) {
477  os.str(name);
478  os.seekp(0, std::ios_base::end);
479  os << ": equation(" << iCnt + 1 << ")";
480  pDf[iCnt].EqDescription = os.str();
481  }
482  }
483 
484  for (unsigned int iCnt = 0; iCnt < iNumDof; iCnt++) {
485  pDf[iCnt].Order = (*i)->GetDofType(iCnt);
486  pDf[iCnt].EqOrder = (*i)->GetEqType(iCnt);
487  }
488  }
489  }
490 
491  /* per ogni elemento */
492  Elem* pEl = NULL;
493  if (ElemIter.bGetFirst(pEl)) {
494  do {
495  ASSERT(pEl != NULL);
496  DEBUGLCOUT(MYDEBUG_INIT|MYDEBUG_ASSEMBLY,
497  "Elem type " << pEl->GetElemType()
498  << " (" << psElemNames[pEl->GetElemType()]
499  << "(" << pEl->GetLabel() << "))" << std::endl);
500 
501  /* chiede all'elemento quanti dof possiede */
502  unsigned int iNumDof = pEl->iGetNumDof();
503  if (iNumDof > 0) {
504  ElemWithDofs* pEWD = Cast<ElemWithDofs>(pEl);
505 
506  ASSERT(pEWD->iGetFirstIndex() >= 0);
507 
508  /* si fa passare il DofOwner */
509  Dof* pDf = &Dofs[pEWD->iGetFirstIndex()];
510 
511 #ifdef DEBUG
512  DEBUGLCOUT(MYDEBUG_INIT|MYDEBUG_ASSEMBLY,
513  psElemNames[pEl->GetElemType()]
514  << "(" << pEWD->GetLabel() << "): "
515  "first dof = " << pDf->iIndex + 1
516  << std::endl);
517 #endif /* DEBUG */
518 
519  if (pds) {
520  unsigned int nd = pEWD->iGetNumDof();
521  integer fd = pDf->iIndex;
522 
523  out_ds << psElemNames[pEWD->GetElemType()]
524  << "(" << pEWD->GetLabel() << "): "
525  << nd << " " << fd + 1;
526  if (nd > 1) {
527  out_ds << "->" << fd + nd;
528  }
529  out_ds << std::endl;
530  if (uPrintFlags & PRINT_DOF_DESCRIPTION) {
531  pEWD->DescribeDof(out_ds,
532  " ");
533  }
534 
535  if (uPrintFlags & PRINT_EQ_DESCRIPTION) {
536  pEWD->DescribeEq(out_ds,
537  " ");
538  }
539  }
540 
541 
542  /* per ogni Dof, chiede all'elemento
543  * di che tipo e' e lo setta
544  * nel DofOwner */
545  std::vector<std::string> DofDesc;
546  pEWD->DescribeDof(DofDesc);
547  if (DofDesc.size() == iNumDof) {
548  for (unsigned int iCnt = 0; iCnt < iNumDof; iCnt++) {
549  pDf[iCnt].Description = DofDesc[iCnt];
550  }
551 
552  } else {
553  std::ostringstream os;
554  os << psElemNames[pEWD->GetElemType()]
555  << "(" << pEWD->GetLabel() << ")";
556  std::string name(os.str());
557 
558  for (unsigned int iCnt = 0; iCnt < iNumDof; iCnt++) {
559  os.str(name);
560  os.seekp(0, std::ios_base::end);
561  os << ": dof(" << iCnt + 1 << ")";
562  pDf[iCnt].Description = os.str();
563  }
564  }
565 
566  std::vector<std::string> EqDesc;
567  pEWD->DescribeEq(EqDesc);
568  if (EqDesc.size() == iNumDof) {
569  for (unsigned int iCnt = 0; iCnt < iNumDof; iCnt++) {
570  pDf[iCnt].EqDescription = EqDesc[iCnt];
571  }
572 
573  } else {
574  std::ostringstream os;
575  os << psElemNames[pEWD->GetElemType()]
576  << "(" << pEWD->GetLabel() << ")";
577  std::string name(os.str());
578 
579  for (unsigned int iCnt = 0; iCnt < iNumDof; iCnt++) {
580  os.str(name);
581  os.seekp(0, std::ios_base::end);
582  os << ": equation(" << iCnt + 1 << ")";
583  pDf[iCnt].EqDescription = os.str();
584  }
585  }
586 
587  for (unsigned int iCnt = 0; iCnt < iNumDof; iCnt++) {
588  pDf[iCnt].Order = pEWD->GetDofType(iCnt);
589  pDf[iCnt].EqOrder = pEWD->GetEqType(iCnt);
590  }
591  }
592  } while (ElemIter.bGetNext(pEl));
593  }
594 
595  /* FIXME: this should rather go before initial assembly */
596  /* NOTE: we run code anyway, but only print if requested/allowed
597  * for consistency checking purposes */
598 
599  /* per ogni elemento */
600  if (ElemIter.bGetFirst(pEl)) {
601  /* create node connection structure */
602  typedef std::set<const Elem *> elmap;
603  typedef std::map<const Node *, elmap *> nodemap;
604  std::vector<nodemap> connectedElems(Node::LASTNODETYPE);
605 
606  /* element connections get populated directly by elements */
607  std::vector<const Node *> connectedNodes;
608 
609  if (uPrintFlags & PRINT_EL_CONNECTION) {
610  out_ds << "Element connections" << std::endl;
611  }
612 
613  do {
614  pEl->GetConnectedNodes(connectedNodes);
615 
616  if (connectedNodes.size() > 0) {
617  if (uPrintFlags & PRINT_EL_CONNECTION) {
618  out_ds << psElemNames[pEl->GetElemType()]
619  << "(" << pEl->GetLabel() << ") connecting" << std::endl;
620  }
621  for (std::vector<const Node *>::const_iterator i = connectedNodes.begin();
622  i != connectedNodes.end(); ++i)
623  {
624  const Node *real_i = (*i)->GetNode();
625  if (uPrintFlags & PRINT_EL_CONNECTION) {
626  out_ds << " "
627  << psNodeNames[real_i->GetNodeType()]
628  << "(" << real_i->GetLabel() << ")" << std::endl;
629  }
630 
631  nodemap::iterator n = connectedElems[real_i->GetNodeType()].find(real_i);
632  if (n == connectedElems[real_i->GetNodeType()].end()) {
633  connectedElems[real_i->GetNodeType()][real_i] = new elmap;
634  }
635  connectedElems[real_i->GetNodeType()][real_i]->insert(pEl);
636  }
637 
638  } else {
639  if (uPrintFlags & PRINT_EL_CONNECTION) {
640  out_ds << psElemNames[pEl->GetElemType()]
641  << "(" << pEl->GetLabel() << ") not connected" << std::endl;
642  }
643  }
644 
645  } while (ElemIter.bGetNext(pEl));
646 
647 
648  if (uPrintFlags & PRINT_NODE_CONNECTION) {
649  out_ds << "Node connections" << std::endl;
650  }
651  for (unsigned t = 0; t < Node::LASTNODETYPE; t++) {
652  for (nodemap::iterator n = connectedElems[t].begin();
653  n != connectedElems[t].end(); ++n)
654  {
655  if (uPrintFlags & PRINT_NODE_CONNECTION) {
656  out_ds << psNodeNames[n->first->GetNodeType()]
657  << "(" << n->first->GetLabel() << ") connected to" << std::endl;
658  }
659  for (elmap::const_iterator e = n->second->begin();
660  e != n->second->end(); ++e)
661  {
662  if (uPrintFlags & PRINT_NODE_CONNECTION) {
663  out_ds << " "
664  << psElemNames[(*e)->GetElemType()]
665  << "(" << (*e)->GetLabel() << ")" << std::endl;
666  }
667  }
668 
669  delete n->second;
670  n->second = 0;
671  }
672  }
673  }
674 } /* End of DataManager::DofOwnerInit() */
675 
676 /* Inizializzazione della struttura dei dof
677  * per l'assemblaggio iniziale dei vincoli */
678 void
679 DataManager::InitialJointAssembly(void)
680 {
681  if ( uPrintFlags & PRINT_TO_FILE ) {
682  OutHdl.Open(OutputHandler::DOFSTATS);
683  }
684 
685  std::ostream& out_ds = (uPrintFlags & PRINT_TO_FILE)
686  ? OutHdl.DofStats()
687  : std::cout;
688 
689  /* Costruisce la struttura temporanea dei Dof */
690 
691  ASSERTMSG(DofData[DofOwner::JOINT].iNum > 0,
692  "Warning, no joints are defined; "
693  "You shouldn't have reached this point");
694  ASSERT(DofData[DofOwner::STRUCTURALNODE].iNum > 0);
695 
696  /* Nodi strutturali: mette gli indici ai DofOwner */
697  bool pds =
698 #ifdef DEBUG
699  DEBUG_LEVEL_MATCH(MYDEBUG_INIT|MYDEBUG_ASSEMBLY) ||
700 #endif /* DEBUG */
701  (!silent_output && (uPrintFlags & PRINT_DOF_STATS));
702 
703  if (pds) {
704  out_ds << "Initial assembly dof stats" << std::endl;
705  }
706 
707  /* Numero totale di Dof durante l'assemblaggio iniziale */
708  integer iInitialNumDofs = 0;
709  for (NodeContainerType::const_iterator i = NodeData[Node::STRUCTURAL].NodeContainer.begin();
710  i != NodeData[Node::STRUCTURAL].NodeContainer.end(); ++i)
711  {
712  iInitialNumDofs += dynamic_cast<const StructDispNode*>(i->second)->iGetInitialNumDof();
713  }
714 
715  /* Elementi: mette gli indici agli eventuali DofOwner */
716  for (int iCnt1 = 0; iCnt1 < Elem::LASTELEMTYPE; iCnt1++) {
717  /* Per ogni tipo di elemento */
718  if (ElemData[iCnt1].bToBeUsedInAssembly() && !ElemData[iCnt1].ElemContainer.empty()) {
719  /* Se deve essere usato nell'assemblaggio e ne sono definiti */
720 
721  /* Tipo di dof dell'elemento corrente */
722  DofOwner::Type CurrDofType =
723  ElemData[iCnt1].DofOwnerType;
724 
725  if (CurrDofType != DofOwner::UNKNOWN) {
726  ASSERT((unsigned)DofData[CurrDofType].iNum == ElemData[iCnt1].ElemContainer.size());
727 
728  /* Iterazione sugli Elem */
729  for (ElemContainerType::const_iterator e = ElemData[iCnt1].ElemContainer.begin();
730  e != ElemData[iCnt1].ElemContainer.end(); ++e)
731  {
732  InitialAssemblyElem *pEl = dynamic_cast<InitialAssemblyElem *>(e->second);
733  if (pEl == 0) {
734  /* Ignore elements
735  * not subjected
736  * to initial assembly */
737  continue;
738  }
739 
740  ElemWithDofs *pDOEl = dynamic_cast<ElemWithDofs *>(e->second);
741  if (pDOEl == 0) {
742  /* Ignore elements subjected
743  * to initial assembly
744  * but without dofs */
745  continue;
746  }
747 
748  iInitialNumDofs += pEl->iGetInitialNumDof();
749  }
750  }
751  }
752  }
753 
754  /*
755  * Alla fine, i DofOwner di nodi e joint contengono gli indici giusti per
756  * l'assemblaggio iniziale. Corrispondono a:
757  * - per ogni nodo:
758  * - posizione x
759  * - parametri di rotazione g
760  * - velocita' xP
761  * - velocita' angolare omega
762  * - per ogni joint:
763  * - se vincolo in posizione, reazione e sua derivata
764  * - se vincolo in velocita', reazione.
765  * - per vincoli misti si hanno reazioni ed eventualmente loro derivate
766  * in base al tipo
767  */
768 
769  /* Creazione e costruzione array Dof */
770  Dofs.resize(iInitialNumDofs);
771 
772  // just to make sure nothing strange occurs...
773  iTotDofs = iInitialNumDofs;
774 
775  integer iIndex; /* Indice dei gradi di liberta' */
776  for (iIndex = 0; iIndex < iInitialNumDofs; iIndex++) {
777  Dofs[iIndex].iIndex = iIndex;
778  }
779 
780  /* mette a posto i dof */
781  iIndex = 0;
782  for (NodeContainerType::const_iterator n = NodeData[Node::STRUCTURAL].NodeContainer.begin();
783  n != NodeData[Node::STRUCTURAL].NodeContainer.end(); ++n)
784  {
785  // numero di dof di un owner
786  const StructDispNode *pNode = dynamic_cast<const StructDispNode *>(n->second);
787  unsigned int iNumDofs = pNode->iGetInitialNumDof();
788  if (iNumDofs > 0) {
789  DofOwner* pFDO = const_cast<DofOwner *>(pNode->pGetDofOwner());
790  pFDO->iNumDofs = iNumDofs;
791  pFDO->iFirstIndex = iIndex;
792  if (pds) {
793  unsigned int nd = iNumDofs;
794  integer fd = iIndex;
795 
796  out_ds << psNodeNames[pNode->GetNodeType()]
797  << "(" << pNode->GetLabel()
798  << "): " << nd << " " << fd + 1;
799  if (nd > 1) {
800  out_ds << "->" << fd + nd;
801  }
802  out_ds << std::endl;
803  if (uPrintFlags & PRINT_DOF_DESCRIPTION) {
804  pNode->DescribeDof(out_ds,
805  " ", true);
806  }
807 
808  if (uPrintFlags & PRINT_EQ_DESCRIPTION) {
809  pNode->DescribeEq(out_ds,
810  " ", true);
811  }
812  }
813 
814  std::vector<std::string> DofDesc;
815  pNode->DescribeDof(DofDesc, true);
816  if (DofDesc.size() == iNumDofs) {
817  for (unsigned iCnt = 0; iCnt < iNumDofs; iCnt++) {
818  Dofs[iIndex + iCnt].Description = DofDesc[iCnt];
819  }
820 
821  } else {
822  std::ostringstream os;
823  os << psNodeNames[pNode->GetNodeType()]
824  << "(" << pNode->GetLabel() << ")";
825  std::string name(os.str());
826 
827  for (unsigned int iCnt = 0; iCnt < iNumDofs; iCnt++) {
828  os.str(name);
829  os.seekp(0, std::ios_base::end);
830  os << ": dof(" << iCnt + 1 << ")";
831  Dofs[iIndex + iCnt].Description = os.str();
832  }
833  }
834 
835  std::vector<std::string> EqDesc;
836  pNode->DescribeEq(EqDesc, true);
837  if (EqDesc.size() == iNumDofs) {
838  for (unsigned iCnt = 0; iCnt < iNumDofs; iCnt++) {
839  Dofs[iIndex + iCnt].EqDescription = EqDesc[iCnt];
840  }
841 
842  } else {
843  std::ostringstream os;
844  os << psNodeNames[pNode->GetNodeType()]
845  << "(" << pNode->GetLabel() << ")";
846  std::string name(os.str());
847 
848  for (unsigned int iCnt = 0; iCnt < iNumDofs; iCnt++) {
849  os.str(name);
850  os.seekp(0, std::ios_base::end);
851  os << ": equation(" << iCnt + 1 << ")";
852  Dofs[iIndex + iCnt].EqDescription = os.str();
853  }
854  }
855 
856  iIndex += iNumDofs;
857 
858  } else {
859  pedantic_cerr(psNodeNames[pNode->GetNodeType()]
860  << "(" << n->second->GetLabel() << ") has 0 dofs"
861  << std::endl);
862  }
863  }
864 
865  /* Elementi: mette gli indici agli eventuali DofOwner */
866  for (int iCnt1 = 0; iCnt1 < Elem::LASTELEMTYPE; iCnt1++) {
867  /* Pre ogni tipo di elemento */
868  if (ElemData[iCnt1].bToBeUsedInAssembly() && !ElemData[iCnt1].ElemContainer.empty()) {
869  /* Se deve essere usato nell'assemblaggio e ne sono definiti */
870 
871  /* Tipo di dof dell'elemento corrente */
872  DofOwner::Type CurrDofType =
873  ElemData[iCnt1].DofOwnerType;
874 
875  if (CurrDofType != DofOwner::UNKNOWN) {
876  ASSERT((unsigned)DofData[CurrDofType].iNum == ElemData[iCnt1].ElemContainer.size());
877 
878  /* Iterazione sugli Elem */
879  for (ElemContainerType::const_iterator p = ElemData[iCnt1].ElemContainer.begin();
880  p != ElemData[iCnt1].ElemContainer.end();
881  ++p)
882  {
883  InitialAssemblyElem *pEl = dynamic_cast<InitialAssemblyElem *>(p->second);
884  if (pEl == 0) {
885  /* Ignore elements
886  * not subjected
887  * to initial assembly */
888  continue;
889  }
890 
891  ElemWithDofs *pDOEl = dynamic_cast<ElemWithDofs *>(p->second);
892  if (pDOEl == 0) {
893  /* Ignore elements subjected
894  * to initial assembly
895  * but without dofs */
896  continue;
897  }
898 
899  DofOwner *pDO = const_cast<DofOwner *>(pDOEl->pGetDofOwner());
900  ASSERT(pDO != 0);
901  // numero di dof di un owner
902  unsigned int iNumDofs = pEl->iGetInitialNumDof();
903  pDO->iNumDofs = iNumDofs;
904  if (iNumDofs > 0) {
905  pDO->iFirstIndex = iIndex;
906  if (pds) {
907  unsigned int nd = iNumDofs;
908  integer fd = iIndex;
909  ElemWithDofs* pEWD = Cast<ElemWithDofs>(p->second);
910 
911  out_ds << psElemNames[pEl->GetElemType()]
912  << "(" << pEl->GetLabel()
913  << "): " << nd << " " << fd + 1;
914  if (nd > 1) {
915  out_ds << "->" << fd + nd;
916  }
917  out_ds << std::endl;
918  if (uPrintFlags & PRINT_DOF_DESCRIPTION) {
919  pEWD->DescribeDof(out_ds,
920  " ", true);
921  }
922 
923  if (uPrintFlags & PRINT_EQ_DESCRIPTION) {
924  pEWD->DescribeEq(out_ds,
925  " ", true);
926  }
927  }
928 
929  std::vector<std::string> DofDesc;
930  pEl->DescribeDof(DofDesc, true);
931  if (DofDesc.size() == iNumDofs) {
932  for (unsigned iCnt = 0; iCnt < iNumDofs; iCnt++) {
933  Dofs[iIndex + iCnt].Description = DofDesc[iCnt];
934  }
935 
936  } else {
937  std::ostringstream os;
938  os << psElemNames[pEl->GetElemType()]
939  << "(" << pEl->GetLabel() << ")";
940  std::string name(os.str());
941 
942  for (unsigned int iCnt = 0; iCnt < iNumDofs; iCnt++) {
943  os.str(name);
944  os.seekp(0, std::ios_base::end);
945  os << ": dof(" << iCnt + 1 << ")";
946  Dofs[iIndex + iCnt].Description = os.str();
947  }
948  }
949 
950  std::vector<std::string> EqDesc;
951  pEl->DescribeEq(EqDesc, true);
952  if (EqDesc.size() == iNumDofs) {
953  for (unsigned iCnt = 0; iCnt < iNumDofs; iCnt++) {
954  Dofs[iIndex + iCnt].EqDescription = EqDesc[iCnt];
955  }
956 
957  } else {
958  std::ostringstream os;
959  os << psElemNames[pEl->GetElemType()]
960  << "(" << pEl->GetLabel() << ")";
961  std::string name(os.str());
962 
963  for (unsigned int iCnt = 0; iCnt < iNumDofs; iCnt++) {
964  os.str(name);
965  os.seekp(0, std::ios_base::end);
966  os << ": equation(" << iCnt + 1 << ")";
967  Dofs[iIndex + iCnt].EqDescription = os.str();
968  }
969  }
970 
971  iIndex += iNumDofs;
972 
973  } else {
974  pedantic_cerr(psElemNames[iCnt1]
975  << "(" << pEl->GetLabel() << ") "
976  "has 0 dofs" << std::endl);
977  }
978  }
979  }
980  }
981  }
982 
983  ASSERT(iIndex == iInitialNumDofs);
984 
985  /* Trova le massime dimensioni del workspace
986  * per l'assemblaggio iniziale */
987  integer iMaxRowsRes = 0;
988  integer iMaxRowsJac = 0;
989  integer iMaxColsJac = 0;
990  integer iMaxItemsJac = 0;
991 
992  InitialAssemblyIterator IAIter(&ElemData);
993  InitialAssemblyElem* pEl = IAIter.GetFirst();
994  while (pEl != NULL) {
995  integer iCurrRows = 0;
996  integer iCurrCols = 0;
997  pEl->InitialWorkSpaceDim(&iCurrRows, &iCurrCols);
998 
999  if (iCurrRows >= 0) {
1000  // Assume a full Jacobian matrix
1001  iMaxRowsJac = std::max(iMaxRowsJac, iCurrRows);
1002  iMaxColsJac = std::max(iMaxColsJac, iCurrCols);
1003  } else {
1004  // Assume a sparse Jacobian matrix
1005  iCurrRows = std::abs(iCurrRows);
1006  }
1007 
1008  iMaxRowsRes = std::max(iMaxRowsRes, iCurrRows);
1009  iMaxItemsJac = std::max(iMaxItemsJac, iCurrRows * iCurrCols);
1010 
1011  pEl = IAIter.GetNext();
1012  }
1013 
1014  /* Ciclo di iterazioni fino a convergenza */
1015 
1016  /* Crea il solutore lineare, tenendo conto dei tipi
1017  * supportati, di quanto scelto nel file di configurazione
1018  * e di eventuali paraametri extra */
1019  SolutionManager* pSM = CurrSolver.GetSolutionManager(iInitialNumDofs);
1020 
1021 #ifdef DEBUG_MEMMANAGER
1022  DEBUGLCOUT(MYDEBUG_MEM|MYDEBUG_ASSEMBLY,
1023  "After initialisation in InitialJointAssembly" << std::endl
1024  << defaultMemoryManager << std::endl);
1025 #endif /* DEBUG_MEMMANAGER */
1026 
1027  MyVectorHandler X(iInitialNumDofs);
1028  X.Reset();
1029 
1030  /* Linka il DriveHandler al vettore soluzione */
1031  LinkToSolution(X, X);
1032 
1033  /* Setta i valori iniziali dei gradi di liberta' dei nodi strutturali
1034  * durante l'assemblaggio iniziale */
1035  for (NodeContainerType::iterator i = NodeData[Node::STRUCTURAL].NodeContainer.begin();
1036  i != NodeData[Node::STRUCTURAL].NodeContainer.end(); ++i)
1037  {
1038  dynamic_cast<StructDispNode *>(i->second)->SetInitialValue(X);
1039  }
1040 
1041  /* Setta i valori iniziali dei gradi di liberta' dei vincoli
1042  * durante l'assemblaggio iniziale */
1043  for (int iCnt1 = 0; iCnt1 < Elem::LASTELEMTYPE; iCnt1++) {
1044  /* Pre ogni tipo di elemento */
1045  if (ElemData[iCnt1].DofOwnerType != DofOwner::UNKNOWN &&
1046  ElemData[iCnt1].bToBeUsedInAssembly() &&
1047  !ElemData[iCnt1].ElemContainer.empty())
1048  {
1049  for (ElemContainerType::const_iterator p = ElemData[iCnt1].ElemContainer.begin();
1050  p != ElemData[iCnt1].ElemContainer.end(); ++p)
1051  {
1052  ElemWithDofs *pEWD = Cast<ElemWithDofs>(p->second);
1053  pEWD->SetInitialValue(X);
1054  }
1055  }
1056  }
1057 
1058  /* Vettore di lavoro */
1059  VectorHandler* pResHdl = pSM->pResHdl();
1060  MySubVectorHandler WorkVec(iMaxRowsRes);
1061 
1062  /* Matrice di lavoro */
1063  MatrixHandler* pMatHdl = pSM->pMatHdl();
1064  VariableSubMatrixHandler WorkMat(iMaxRowsJac, iMaxColsJac, iMaxItemsJac);
1065 
1066  /* Soluzione */
1067  VectorHandler* pSolHdl = pSM->pSolHdl();
1068 
1069  if (
1070 #ifdef DEBUG
1071  DEBUG_LEVEL_MATCH(MYDEBUG_ASSEMBLY) ||
1072 #endif /* DEBUG */
1073  outputIters())
1074  {
1075  silent_cout("Assembly Tol=" << dInitialAssemblyTol << std::endl);
1076  }
1077 
1078  /* Ciclo di assemblaggio */
1079  for (integer iNumIter = 0; ; iNumIter++) {
1080  /* Assemblo il residuo */
1081  pResHdl->Reset();
1082 
1083  /* Contributo dei nodi */
1084  for (NodeContainerType::iterator i = NodeData[Node::STRUCTURAL].NodeContainer.begin();
1085  i != NodeData[Node::STRUCTURAL].NodeContainer.end(); ++i)
1086  {
1087  const StructDispNode *pDispNode = dynamic_cast<const StructDispNode *>(i->second);
1088  const StructNode *pNode = dynamic_cast<const StructNode *>(i->second);
1089  if (pNode && pNode->GetStructNodeType() == StructNode::DUMMY) {
1090  continue;
1091  }
1092 
1093  int iOffset = 3;
1094  if (pNode) {
1095  iOffset = 6;
1096  }
1097 
1098  integer iFirstIndex = pDispNode->iGetFirstPositionIndex();
1099 
1100  /* Nuova feature: ogni nodo ha la sua stiffness */
1101  doublereal dPosStiff = pDispNode->dGetPositionStiffness();
1102  doublereal dVelStiff = pDispNode->dGetVelocityStiffness();
1103 
1104  Vec3 TmpVec;
1105 
1106  /* Posizione: k*Delta_x = k(x_0-x) + F */
1107  TmpVec = pDispNode->GetXPrev() - pDispNode->GetXCurr();
1108  pResHdl->Add(iFirstIndex + 1, TmpVec*dPosStiff);
1109 
1110  if (pNode) {
1111  /* Rotazione: k*Delta_g = -k*g(R_Delta) + M */
1112  Mat3x3 RDelta = pNode->GetRPrev().MulMT(pNode->GetRCurr());
1113  TmpVec = Vec3(CGR_Rot::Param, RDelta);
1114  pResHdl->Add(iFirstIndex + 4, TmpVec*dPosStiff);
1115 
1116  /* Velocita' angolare: k*(Delta_w+(R_Delta*w0)/\Delta_g) =
1117  * k*(R_Delta*w0-w) + M */
1118  const Vec3& wPrev(pNode->GetWPrev());
1119  const Vec3& wCurr(pNode->GetWCurr());
1120 
1121  if (pNode->bOmegaRotates()) {
1122  /* con questa la velocita' angolare e' solidale
1123  * con il nodo */
1124  TmpVec = RDelta*wPrev - wCurr;
1125  } else {
1126  /* con questa la velocita' angolare e' solidale
1127  * col riferimento assoluto */
1128  TmpVec = wPrev - wCurr;
1129  }
1130 
1131  pResHdl->Add(iFirstIndex + iOffset + 4, TmpVec*dVelStiff);
1132  }
1133 
1134  /* Velocita': k*Delta_v = k*(v0-Delta_v) + F */
1135  TmpVec = pDispNode->GetVPrev() - pDispNode->GetVCurr();
1136  pResHdl->Add(iFirstIndex + iOffset + 1, TmpVec*dVelStiff);
1137 
1138  }
1139 
1140  /* Elementi (con iteratore): */
1141  pEl = IAIter.GetFirst();
1142  while (pEl != NULL) {
1143  try {
1144  *pResHdl += pEl->InitialAssRes(WorkVec, X);
1145  }
1146  catch (Elem::ChangedEquationStructure) {
1147  // do nothing: Jacobian matrix
1148  // is always recomputed anyway...
1149  }
1150  pEl = IAIter.GetNext();
1151  }
1152 
1153  if (
1154 #ifdef DEBUG
1155  DEBUG_LEVEL_MATCH(MYDEBUG_ASSEMBLY|MYDEBUG_RESIDUAL) ||
1156 #endif /* DEBUG */
1157  outputRes())
1158  {
1159  /* Output del residuo */
1160  PrintResidual(*pResHdl, iNumIter);
1161  }
1162 
1163  /* Eseguo il test di convergenza; se e' positivo, esco */
1164  /* FIXME: why /(1.+X.Dot()) ??? */
1165  doublereal dTest = pResHdl->Dot()/(1. + X.Dot());
1166  if (!std::isfinite(dTest)) {
1167  silent_cerr("Assembly diverged; aborting..." << std::endl);
1168  throw DataManager::ErrAssemblyDiverged(MBDYN_EXCEPT_ARGS);
1169  }
1170  dTest = sqrt(dTest);
1171 
1172  if ((dTest <= dInitialAssemblyTol ||
1173  iNumIter >= iMaxInitialIterations) && (
1174 #ifdef DEBUG
1175  DEBUG_LEVEL_MATCH(MYDEBUG_ASSEMBLY) ||
1176 #endif /* DEBUG */
1177  outputIters()))
1178  {
1179  silent_cout("\tIteration(" << iNumIter << ") "
1180  "" << dTest << std::endl);
1181  }
1182 
1183  /* Se la tolleranza e' raggiunta, esce dal ciclo */
1184  if (dTest <= dInitialAssemblyTol) {
1185  DEBUGLCOUT(MYDEBUG_ASSEMBLY, "Initial assembly "
1186  "performed successfully in "
1187  << iNumIter << " iterations"
1188  << std::endl);
1189  goto endofcycle;
1190  }
1191 
1192  /* Se ho raggiunto il numero massimo di iterazioni */
1193  if (iNumIter >= iMaxInitialIterations) {
1194  silent_cerr("Initial assembly iterations "
1195  "reached maximum number "
1196  << iMaxInitialIterations << "; aborting..."
1197  << std::endl);
1198  throw DataManager::ErrAssemblyMaxIters(MBDYN_EXCEPT_ARGS);
1199  }
1200 
1201  /* Assemblo lo jacobiano e risolvo */
1202  pSM->MatrInitialize();
1203 
1204  /* Contributo dei nodi */
1205  for (NodeContainerType::iterator i = NodeData[Node::STRUCTURAL].NodeContainer.begin();
1206  i != NodeData[Node::STRUCTURAL].NodeContainer.end(); ++i)
1207  {
1208  const StructDispNode *pDispNode = dynamic_cast<const StructDispNode *>(i->second);
1209  const StructNode *pNode = dynamic_cast<const StructNode *>(i->second);
1210  if (pNode && pNode->GetStructNodeType() == StructNode::DUMMY) {
1211  continue;
1212  }
1213 
1214  int iOffset = 3;
1215  if (pNode) {
1216  iOffset = 6;
1217  }
1218 
1219  integer iFirstIndex = pDispNode->iGetFirstPositionIndex();
1220 
1221  /* Nuova feature: ogni nodo ha la sua stiffness */
1222  doublereal dPosStiff = pDispNode->dGetPositionStiffness();
1223  doublereal dVelStiff = pDispNode->dGetVelocityStiffness();
1224 
1225  /* NOTE: iOffset equal to number of position/velocity equations */
1226  for (int iCnt = 1; iCnt <= iOffset; iCnt++) {
1227  /* Posizione, rotazione */
1228  integer iTmp = iFirstIndex + iCnt;
1229  pMatHdl->PutCoef(iTmp, iTmp, dPosStiff);
1230 
1231  /* Velocita', velocita' angolare */
1232  iTmp += iOffset;
1233  pMatHdl->PutCoef(iTmp, iTmp, dVelStiff);
1234  }
1235 
1236  if (pNode && pNode->bOmegaRotates()) {
1237  /* con questi la velocita' angolare e' solidale con il nodo */
1238 
1239  /* Velocita' angolare - termine di rotazione: R_Delta*w0/\ */
1240  const Mat3x3& R0 = pNode->GetRPrev();
1241  const Mat3x3& R = pNode->GetRCurr();
1242  const Vec3& W0 = pNode->GetWPrev();
1243  Vec3 TmpVec = R*(R0.MulTV(W0*dVelStiff));
1244 
1245  /* W1 in m(3, 2), -W1 in m(2, 3) */
1246  doublereal d = TmpVec(1);
1247  pMatHdl->PutCoef(iFirstIndex + iOffset + 6, iFirstIndex + 5, d);
1248  pMatHdl->PutCoef(iFirstIndex + iOffset + 5, iFirstIndex + 6, -d);
1249 
1250  /* W2 in m(1, 3), -W2 in m(3, 1) */
1251  d = TmpVec(2);
1252  pMatHdl->PutCoef(iFirstIndex + iOffset + 4, iFirstIndex + 6, d);
1253  pMatHdl->PutCoef(iFirstIndex + iOffset + 6, iFirstIndex + 4, -d);
1254 
1255  /* W3 in m(2, 1), -W3 in m(1, 2) */
1256  d = TmpVec(3);
1257  pMatHdl->PutCoef(iFirstIndex + iOffset + 5, iFirstIndex + 4, d);
1258  pMatHdl->PutCoef(iFirstIndex + iOffset + 4, iFirstIndex + 5, -d);
1259  } /* altrimenti la velocita' angolare e' solidale con il nodo */
1260  }
1261 
1262  /* Contributo degli elementi */
1263  pEl = IAIter.GetFirst();
1264  while (pEl != NULL) {
1265  *pMatHdl += pEl->InitialAssJac(WorkMat, X);
1266  pEl = IAIter.GetNext();
1267  }
1268 
1269  if (
1270 #ifdef DEBUG
1271  DEBUG_LEVEL_MATCH(MYDEBUG_ASSEMBLY|MYDEBUG_JAC) ||
1272 #endif /* DEBUG */
1273  outputJac())
1274  {
1275  silent_cout("Jacobian:" << std::endl
1276  << *pMatHdl);
1277  }
1278 
1279  /* Fattorizza e risolve con jacobiano e residuo appena calcolati */
1280  try {
1281  pSM->Solve();
1282  }
1283  catch (LinearSolver::ErrFactor& err) {
1284  silent_cerr("Initial assembly failed because no pivot element "
1285  "could be found for column " << err.iCol
1286  << " (" << GetDofDescription(err.iCol) << "); "
1287  "aborting..." << std::endl);
1288  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1289  }
1290 
1291  if (
1292 #ifdef DEBUG
1293  DEBUG_LEVEL_MATCH(MYDEBUG_ASSEMBLY|MYDEBUG_RESIDUAL) ||
1294 #endif /* DEBUG */
1295  outputSol())
1296  {
1297  /* Output della soluzione */
1298  PrintSolution(*pSolHdl, iNumIter);
1299  }
1300 
1301  if (outputIters() || outputSolverConditionNumber()) {
1302  if (outputIters()) {
1303  silent_cout("\tIteration(" << iNumIter << ") " << std::setw(12) << dTest << " J");
1304  }
1305 
1306  if (outputSolverConditionNumber()) {
1307  silent_cout(" cond=");
1308  doublereal dCond;
1309  if (pSM->bGetConditionNumber(dCond)) {
1310  silent_cout(dCond);
1311  } else {
1312  silent_cout("NA");
1313  }
1314  }
1315 
1316  silent_cout(std::endl);
1317  }
1318 
1319  /* Aggiorno la soluzione */
1320  if (dEpsilon != 1.) {
1321  *pSolHdl *= dEpsilon;
1322  }
1323  X += *pSolHdl;
1324 
1325  /* Correggo i nodi */
1326  for (NodeContainerType::iterator i = NodeData[Node::STRUCTURAL].NodeContainer.begin();
1327  i != NodeData[Node::STRUCTURAL].NodeContainer.end(); ++i)
1328  {
1329  dynamic_cast<StructDispNode *>(i->second)->InitialUpdate(X);
1330  }
1331  }
1332 
1333 endofcycle:
1334  /* Resetta e distrugge la struttura temporanea dei Dof */
1335 
1336  /* Elementi: rimette a posto il numero di Dof propri dei vincoli */
1337  for (int iCnt1 = 0; iCnt1 < Elem::LASTELEMTYPE; iCnt1++) {
1338  /* Per ogni tipo di elemento */
1339  if (ElemData[iCnt1].DofOwnerType != DofOwner::UNKNOWN &&
1340  ElemData[iCnt1].bToBeUsedInAssembly() &&
1341  !ElemData[iCnt1].ElemContainer.empty())
1342  {
1343  /* Se possiede dofs, se deve essere usato nell'assemblaggio
1344  * e se ne sono presenti */
1345  for (ElemContainerType::const_iterator p = ElemData[iCnt1].ElemContainer.begin();
1346  p != ElemData[iCnt1].ElemContainer.end();
1347  ++p)
1348  {
1349  ElemWithDofs *pEWD = Cast<ElemWithDofs>(p->second);
1350  DofOwner *pDO = const_cast<DofOwner *>(pEWD->pGetDofOwner());
1351  pDO->iNumDofs = p->second->iGetNumDof();
1352  }
1353  }
1354  }
1355 
1356  /* Dealloca il vettore dei Dof */
1357  ASSERT(!Dofs.empty());
1358 
1359  // restore
1360  iTotDofs = 0;
1361 
1362  SAFEDELETE(pSM);
1363 } /* End of InitialJointAssembly */
1364 
1365 /* Aggiorna i DofOwner con il numero di dofs dell'elemento */
1366 
1367 void
1368 DataManager::DofOwnerSet(void)
1369 {
1370  DEBUGCOUTFNAME("DataManager::DofOwnerSet");
1371 
1372  /* Setta i DofOwner dei nodi */
1373  for (NodeVecType::iterator i = Nodes.begin(); i != Nodes.end(); ++i) {
1374  DofOwner* pDO = const_cast<DofOwner *>((*i)->pGetDofOwner());
1375  pDO->iNumDofs = (*i)->iGetNumDof();
1376  }
1377 
1378  /* Setta i DofOwner degli elementi (chi li possiede) */
1379  for (int iCnt = 0; iCnt < Elem::LASTELEMTYPE; iCnt++) {
1380  DofOwner::Type DT = ElemData[iCnt].DofOwnerType;
1381  if (DT != DofOwner::UNKNOWN) {
1382  DEBUGLCOUT(MYDEBUG_INIT, "Elem type " << iCnt
1383  << " (" << psElemNames[iCnt] << ")"
1384  << std::endl);
1385 
1386  for (ElemContainerType::const_iterator p = ElemData[iCnt].ElemContainer.begin();
1387  p != ElemData[iCnt].ElemContainer.end(); ++p)
1388  {
1389  ElemWithDofs* pEWD = Cast<ElemWithDofs>(p->second);
1390 
1391  DEBUGLCOUT(MYDEBUG_INIT, " " << psElemNames[pEWD->GetElemType()]
1392  << "(" << pEWD->GetLabel() << ")" << std::endl);
1393 
1394  DofOwner* pDO = const_cast<DofOwner *>(pEWD->pGetDofOwner());
1395  pDO->iNumDofs = pEWD->iGetNumDof();
1396  DEBUGLCOUT(MYDEBUG_INIT, " num dofs: " << pDO->iNumDofs << std::endl);
1397  }
1398  }
1399  }
1400 } /* end of DofOwnerSet() */
1401 
1402 
1403 void
1404 DataManager::SetValue(VectorHandler& X, VectorHandler& XP)
1405 {
1406  /* Nodi */
1407  for (NodeVecType::iterator i = Nodes.begin(); i != Nodes.end(); ++i) {
1408  (*i)->SetValue(this, X, XP);
1409  }
1410 
1411  /* Elementi */
1412  /* Versione con iteratore: */
1413  Elem* pEl = NULL;
1414  if (ElemIter.bGetFirst(pEl)) {
1415  do {
1416  pEl->SetValue(this, X, XP);
1417  } while (ElemIter.bGetNext(pEl));
1418  }
1419  if (solArrFileName != NULL) {
1420  std::ifstream fp(solArrFileName);
1421 #ifdef HAVE_ISOPEN
1422  if (!fp.is_open()) {
1423  silent_cerr("DataManager::SetValue(): "
1424  "Cannot open file \"" << solArrFileName << "\""
1425  << std::endl);
1426  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1427  }
1428 #endif /* HAVE_ISOPEN */
1429  fp.read((char*)X.pdGetVec() , X.iGetSize()*sizeof(double));
1430  if (fp.gcount() != std::streamsize(X.iGetSize()*sizeof(double))) {
1431  silent_cerr("DataManager::SetValue(): "
1432  "File(" << solArrFileName << ") too short!"
1433  << std::endl);
1434  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1435  }
1436  fp.read((char*)XP.pdGetVec() , XP.iGetSize()*sizeof(double));
1437  if (fp.gcount() != std::streamsize(XP.iGetSize()*sizeof(double))) {
1438  silent_cerr("DataManager::SetValue(): "
1439  "File(" << solArrFileName << ") too short!"
1440  << std::endl);
1441  throw DataManager::ErrGeneric(MBDYN_EXCEPT_ARGS);
1442  }
1443  SAFEDELETEARR(solArrFileName);
1444  fp.close();
1445  }
1446 } /* End of SetValue */
1447 
1448 
1449 /* Output dati */
1450 void
1451 DataManager::OutputPrepare(void)
1452 {
1453 #ifdef USE_NETCDF
1454  /* Set up NetCDF stuff if required */
1455  if (OutHdl.UseNetCDF(OutputHandler::NETCDF)) {
1456  OutHdl.Open(OutputHandler::NETCDF);
1457  ASSERT(OutHdl.IsOpen(OutputHandler::NETCDF));
1458 
1459  Var_Step = OutHdl.CreateVar<integer>("run.step", "-", "time step index");
1460  Var_Time = OutHdl.CreateVar<doublereal>("time", "s", "simulation time");
1461  Var_TimeStep = OutHdl.CreateVar<doublereal>("run.timestep", "s", "integration time step");
1462  }
1463 #endif /* USE_NETCDF */
1464 
1465  /* Dati dei nodi */
1466  NodeOutputPrepare(OutHdl);
1467 
1468  /* Dati degli elementi */
1469  ElemOutputPrepare(OutHdl);
1470 }
1471 
1472 /* Output setup for Eigenanalysis parameters */
1473 void
1474 DataManager::OutputEigPrepare(const integer iNumAnalyses, const integer iSize)
1475 {
1476 #ifdef USE_NETCDF
1477  /* Set up additional NetCDF stuff for eigenanalysis output */
1478  if (OutHdl.UseNetCDF(OutputHandler::NETCDF)) {
1479 
1480  OutputHandler::NcDimVec dim(1);
1481  dim[0] = OutHdl.CreateDim("eigensolutions", iNumAnalyses);
1482 
1483  m_Dim_Eig_iSize = OutHdl.CreateDim("eig_iSize", iSize);
1484  m_Dim_Eig_iComplex = OutHdl.CreateDim("complex_var_dim", 2);
1485 
1486  OutputHandler::AttrValVec attrs2(2);
1487  attrs2[0] = OutputHandler::AttrVal("type", "integer");
1488  attrs2[1] = OutputHandler::AttrVal("description",
1489  "timestep index of eigensolution");
1490 
1491  Var_Eig_lStep = OutHdl.CreateVar("eig.step", ncInt, attrs2, dim);
1492 
1493  OutputHandler::AttrValVec attrs3(3);
1494  attrs3[0] = OutputHandler::AttrVal("units", "s");
1495  attrs3[1] = OutputHandler::AttrVal("type", "doublereal");
1496  attrs3[2] = OutputHandler::AttrVal("description",
1497  "simulation time at which the eigensolution was computed");
1498 
1499  Var_Eig_dTime = OutHdl.CreateVar("eig.time", ncDouble, attrs3, dim);
1500 
1501  attrs3[0] = OutputHandler::AttrVal("units", "-");
1502  attrs3[1] = OutputHandler::AttrVal("type", "doublereal");
1503  attrs3[2] = OutputHandler::AttrVal("description",
1504  "coefficient used to build Aplus and Aminus matrices");
1505 
1506  Var_Eig_dCoef = OutHdl.CreateVar("eig.dCoef", ncDouble, attrs3, dim);
1507 
1508  OutputHandler::NcDimVec dim2(2);
1509  integer iNumNodes = NodeData[Node::STRUCTURAL].NodeContainer.size();
1510 
1511  dim2[0] = dim[0];
1512  dim2[1] = OutHdl.CreateDim("eig_iIdxSize", iNumNodes);
1513 
1514  attrs2[0] = OutputHandler::AttrVal("type", "integer");
1515  attrs2[1] = OutputHandler::AttrVal("description",
1516  "structural nodes base index");
1517 
1518  Var_Eig_Idx = OutHdl.CreateVar("eig.idx", ncInt, attrs2, dim2);
1519  }
1520 #endif /* USE_NETCDF */
1521 }
1522 
1523 /* Output of Eigenanalysis parameters */
1524 void
1525 DataManager::OutputEigParams(const doublereal& dTime,
1526  const doublereal& dCoef,
1527  const unsigned uCurrEigSol,
1528  const int iResultsPrecision)
1529 {
1530  if (OutHdl.UseText(OutputHandler::EIGENANALYSIS)) {
1531  std::ostream& out = OutHdl.Eigenanalysis();
1532 
1533  if (iResultsPrecision) {
1534  // 7 = number of characters requested by scientific notation
1535  int iNewWidth = iResultsPrecision + 7;
1536  out.width(iNewWidth);
1537  out.precision(iResultsPrecision);
1538  }
1539 
1540  // header
1541  out
1542  << "% time: " << dTime << std::endl;
1543  out
1544  << "dTime = " << dTime << ';' << std::endl;
1545 
1546  // coefficient
1547  out
1548  << "% coefficient" << std::endl
1549  << "dCoef = " << dCoef << ";" << std::endl;
1550  }
1551 #ifdef USE_NETCDF
1552  if (OutHdl.UseNetCDF(OutputHandler::NETCDF)) {
1553 
1554  long lStep = OutHdl.GetCurrentStep();
1555 
1556  Var_Eig_dTime->put_rec(&dTime, uCurrEigSol);
1557  Var_Eig_lStep->put_rec(&lStep, uCurrEigSol);
1558  Var_Eig_dCoef->put_rec(&dCoef, uCurrEigSol);
1559 
1560  }
1561 #endif /* USE_NETCDF */
1562 }
1563 
1564 void
1565 DataManager::OutputEigFullMatrices(const MatrixHandler* pMatA,
1566  const MatrixHandler* pMatB,
1567  const unsigned uCurrEigSol,
1568  const int iMatrixPrecision)
1569 {
1570  const FullMatrixHandler& MatB = dynamic_cast<const FullMatrixHandler &>(*pMatB);
1571  const FullMatrixHandler& MatA = dynamic_cast<const FullMatrixHandler &>(*pMatA);
1572  integer nrows = MatB.iGetNumRows();
1573  integer ncols = MatB.iGetNumCols();
1574 
1575 #ifdef USE_NETCDF
1576  if (OutHdl.UseNetCDF(OutputHandler::NETCDF)) {
1577  OutputHandler::NcDimVec dim2(2);
1578  dim2[0] = m_Dim_Eig_iSize;
1579  dim2[1] = m_Dim_Eig_iSize;
1580 
1581  OutputHandler::AttrValVec attrs3(3);
1582  attrs3[0] = OutputHandler::AttrVal("units", "-");
1583  attrs3[1] = OutputHandler::AttrVal("type", "doublereal");
1584  attrs3[2] = OutputHandler::AttrVal("description", "F/xPrime + dCoef * F/x");
1585 
1586  std::stringstream varname_ss;
1587  varname_ss << "eig." << uCurrEigSol << ".Aplus";
1588  Var_Eig_dAplus = OutHdl.CreateVar(varname_ss.str(),ncDouble, attrs3, dim2);
1589  attrs3[2] = OutputHandler::AttrVal("description", "F/xPrime - dCoef * F/x");
1590 
1591  varname_ss.str("");
1592  varname_ss.clear();
1593  varname_ss << "eig." << uCurrEigSol << ".Aminus";
1594  Var_Eig_dAminus = OutHdl.CreateVar(varname_ss.str(),ncDouble, attrs3, dim2);
1595 
1596  Var_Eig_dAplus->put(MatB.pdGetMat(), nrows, ncols);
1597  Var_Eig_dAminus->put(MatA.pdGetMat(), nrows, ncols);
1598 
1599  }
1600 #endif /* USE_NETCDF */
1601 
1602  if (OutHdl.UseText(OutputHandler::EIGENANALYSIS)) {
1603  std::ostream& out = OutHdl.Eigenanalysis();
1604 
1605  if (iMatrixPrecision) {
1606  // 7 = number of characters requested by scientific notation
1607  int iNewWidth = iMatrixPrecision + 7;
1608  out.width(iNewWidth);
1609  out.precision(iMatrixPrecision);
1610  }
1611 
1612  out
1613  << "% F/xPrime + dCoef * F/x" << std::endl
1614  << "Aplus" << " = [";
1615 
1616 
1617  for (integer r = 1; r <= nrows; r++) {
1618  for (integer c = 1; c <= ncols; c++) {
1619  out << MatB(r, c) << ' ';
1620  }
1621 
1622  if (r == nrows) {
1623  out << "];" << std::endl;
1624 
1625  } else {
1626  out << ";" << std::endl;
1627  }
1628  }
1629 
1630  out
1631  << "% F/xPrime - dCoef * F/x" << std::endl
1632  << "Aminus" << " = [";
1633 
1634  for (integer r = 1; r <= nrows; r++) {
1635  for (integer c = 1; c <= ncols; c++) {
1636  out << MatA(r, c) << ' ';
1637  }
1638 
1639  if (r == nrows) {
1640  out << "];" << std::endl;
1641 
1642  } else {
1643  out << ";" << std::endl;
1644  }
1645  }
1646 
1647  }
1648 }
1649 
1650 void
1651 DataManager::OutputEigSparseMatrices(const MatrixHandler* pMatA,
1652  const MatrixHandler* pMatB,
1653  const unsigned uCurrEigSol,
1654  const int iMatrixPrecision)
1655 {
1656  const SpMapMatrixHandler& MatB = dynamic_cast<const SpMapMatrixHandler &>(*pMatB);
1657  const SpMapMatrixHandler& MatA = dynamic_cast<const SpMapMatrixHandler &>(*pMatA);
1658 
1659  if (OutHdl.UseText(OutputHandler::EIGENANALYSIS)) {
1660  std::ostream& out = OutHdl.Eigenanalysis();
1661 
1662  if (iMatrixPrecision) {
1663  // 7 = number of characters requested by scientific notation
1664  int iNewWidth = iMatrixPrecision + 7;
1665  out.width(iNewWidth);
1666  out.precision(iMatrixPrecision);
1667  }
1668 
1669  out
1670  << "% F/xPrime + dCoef *F/x" << std::endl
1671  << "Aplus" << " = [";
1672 
1673  for (SpMapMatrixHandler::const_iterator i = MatB.begin();
1674  i != MatB.end(); ++i)
1675  {
1676  if (i->dCoef != 0.) {
1677  out << i->iRow + 1 << " " << i->iCol + 1 << " " << i->dCoef << ";" << std::endl;
1678  }
1679  }
1680 
1681  out << "];" << std::endl
1682  << "Aplus = spconvert(Aplus);" << std::endl;
1683 
1684  out
1685  << "% F/xPrime - dCoef *F/x" << std::endl
1686  << "Aminus" << " = [";
1687 
1688  for (SpMapMatrixHandler::const_iterator i = MatA.begin();
1689  i != MatA.end(); ++i)
1690  {
1691  if (i->dCoef != 0.) {
1692  out << i->iRow + 1 << " " << i->iCol + 1 << " " << i->dCoef << ";" << std::endl;
1693  }
1694  }
1695 
1696  out << "];" << std::endl
1697  << "Aminus = spconvert(Aminus);" << std::endl;
1698  }
1699 #ifdef USE_NETCDF
1700  if (OutHdl.UseNetCDF(OutputHandler::NETCDF)) {
1701  OutputHandler::NcDimVec dim2(2);
1702  std::stringstream dimname_ss;
1703  dimname_ss << "eig_" << uCurrEigSol << "_Aplus_sp_iSize";
1704  dim2[0] = OutHdl.CreateDim(dimname_ss.str(), MatB.Nz());
1705  dim2[1] = OutHdl.DimV3();
1706 
1707  OutputHandler::AttrValVec attrs3(3);
1708  attrs3[0] = OutputHandler::AttrVal("units", "-");
1709  attrs3[1] = OutputHandler::AttrVal("type", "doublereal");
1710  attrs3[2] = OutputHandler::AttrVal("description", "F/xPrime + dCoef * F/x");
1711 
1712  dimname_ss.str("");
1713  dimname_ss.clear();
1714  dimname_ss << "eig_" << uCurrEigSol << "_Aminus_sp_iSize";
1715  dim2[0] = OutHdl.CreateDim(dimname_ss.str(), MatA.Nz());
1716 
1717  std::stringstream varname_ss;
1718  varname_ss << "eig." << uCurrEigSol << ".Aplus";
1719  Var_Eig_dAplus = OutHdl.CreateVar(varname_ss.str(), ncDouble, attrs3, dim2);
1720  attrs3[2] = OutputHandler::AttrVal("description", "F/xPrime - dCoef * F/x");
1721 
1722  varname_ss.str("");
1723  varname_ss.clear();
1724  varname_ss << "eig." << uCurrEigSol << ".Aminus";
1725  Var_Eig_dAminus = OutHdl.CreateVar(varname_ss.str(), ncDouble, attrs3, dim2);
1726 
1727  int iCnt = 0;
1728  Vec3 v;
1729  for (SpMapMatrixHandler::const_iterator i = MatB.begin();
1730  i != MatB.end(); ++i)
1731  {
1732  if (i->dCoef != 0.) {
1733  v = Vec3(i->iRow, i->iCol, i->dCoef);
1734  Var_Eig_dAplus->put_rec(v.pGetVec(), iCnt);
1735  iCnt++;
1736  }
1737  }
1738 
1739  iCnt = 0;
1740  for (SpMapMatrixHandler::const_iterator j = MatA.begin();
1741  j != MatA.end(); ++j)
1742  {
1743  if (j->dCoef != 0.) {
1744  v = Vec3(j->iRow, j->iCol, j->dCoef);
1745  Var_Eig_dAminus->put_rec(v.pGetVec(), iCnt);
1746  }
1747  }
1748 
1749  }
1750 #endif
1751 }
1752 
1753 void
1754 DataManager::OutputEigNaiveMatrices(const MatrixHandler* pMatA,
1755  const MatrixHandler* pMatB,
1756  const unsigned uCurrEigSol,
1757  const int iMatrixPrecision)
1758 {
1759  const NaiveMatrixHandler& MatB = dynamic_cast<const NaiveMatrixHandler &>(*pMatB);
1760  const NaiveMatrixHandler& MatA = dynamic_cast<const NaiveMatrixHandler &>(*pMatA);
1761 
1762  if (OutHdl.UseText(OutputHandler::EIGENANALYSIS)) {
1763  std::ostream& out = OutHdl.Eigenanalysis();
1764 
1765  if (iMatrixPrecision) {
1766  // 7 = number of characters requested by scientific notation
1767  int iNewWidth = iMatrixPrecision + 7;
1768  out.width(iNewWidth);
1769  out.precision(iMatrixPrecision);
1770  }
1771 
1772  out
1773  << "% F/xPrime + dCoef *F/x" << std::endl
1774  << "Aplus" << " = [";
1775 
1776  for (NaiveMatrixHandler::const_iterator i = MatB.begin();
1777  i != MatB.end(); ++i)
1778  {
1779  if (i->dCoef != 0.) {
1780  out << i->iRow + 1 << " " << i->iCol + 1 << " " << i->dCoef << ";" << std::endl;
1781  }
1782  }
1783 
1784  out << "];" << std::endl
1785  << "Aplus = spconvert(Aplus);" << std::endl;
1786 
1787  out
1788  << "% F/xPrime - dCoef *F/x" << std::endl
1789  << "Aminus" << " = [";
1790 
1791  for (NaiveMatrixHandler::const_iterator j = MatA.begin();
1792  j != MatA.end(); ++j)
1793  {
1794  if (j->dCoef != 0.) {
1795  out << j->iRow + 1 << " " << j->iCol + 1 << " " << j->dCoef << ";" << std::endl;
1796  }
1797  }
1798 
1799  out << "];" << std::endl
1800  << "Aminus = spconvert(Aminus);" << std::endl;
1801  }
1802 #ifdef USE_NETCDF
1803  if (OutHdl.UseNetCDF(OutputHandler::NETCDF)) {
1804 
1805  OutputHandler::NcDimVec dim2(2);
1806 
1807  std::stringstream dimname_ss;
1808  dimname_ss << "eig_" << uCurrEigSol << "_Aplus_sp_iSize";
1809 
1810  // FIXME: Is there a more efficient way of doing this??
1811  integer iMatBNz = 0;
1812  for (NaiveMatrixHandler::const_iterator i = MatB.begin();
1813  i != MatB.end(); ++i)
1814  {
1815  if (i->dCoef != 0.) {
1816  iMatBNz++;
1817  }
1818  }
1819 
1820  integer iMatANz = 0;
1821  for (NaiveMatrixHandler::const_iterator j = MatA.begin();
1822  j != MatA.end(); ++j)
1823  {
1824  if(j->dCoef != 0.) {
1825  iMatANz++;
1826  }
1827  }
1828 
1829  dim2[0] = OutHdl.CreateDim(dimname_ss.str(), iMatBNz);
1830  dim2[1] = OutHdl.DimV3();
1831 
1832  OutputHandler::AttrValVec attrs3(3);
1833  attrs3[0] = OutputHandler::AttrVal("units", "-");
1834  attrs3[1] = OutputHandler::AttrVal("type", "doublereal");
1835  attrs3[2] = OutputHandler::AttrVal("description", "F/xPrime + dCoef * F/x");
1836 
1837  std::stringstream varname_ss;
1838  varname_ss << "eig." << uCurrEigSol << ".Aplus";
1839  Var_Eig_dAplus = OutHdl.CreateVar(varname_ss.str(), ncDouble, attrs3, dim2);
1840 
1841  dimname_ss.str("");
1842  dimname_ss.clear();
1843  dimname_ss << "eig_" << uCurrEigSol << "_Aminus_sp_iSize";
1844  dim2[0] = OutHdl.CreateDim(dimname_ss.str(), iMatANz);
1845 
1846  attrs3[2] = OutputHandler::AttrVal("description", "F/xPrime - dCoef * F/x");
1847 
1848  varname_ss.str("");
1849  varname_ss.clear();
1850  varname_ss << "eig." << uCurrEigSol << ".Aminus";
1851  Var_Eig_dAminus = OutHdl.CreateVar(varname_ss.str(), ncDouble, attrs3, dim2);
1852 
1853  int iCnt = 0;
1854  Vec3 v;
1855  for (NaiveMatrixHandler::const_iterator i = MatB.begin();
1856  i != MatB.end(); ++i)
1857  {
1858  if (i->dCoef != 0.) {
1859  v = Vec3(i->iRow, i->iCol, i->dCoef);
1860  Var_Eig_dAplus->put_rec(v.pGetVec(), iCnt);
1861  iCnt++;
1862  }
1863  }
1864 
1865  iCnt = 0;
1866  for (NaiveMatrixHandler::const_iterator i = MatA.begin();
1867  i != MatA.end(); ++i)
1868  {
1869  if (i->dCoef != 0.) {
1870  v = Vec3(i->iRow, i->iCol, i->dCoef);
1871  Var_Eig_dAminus->put_rec(v.pGetVec(), iCnt);
1872  }
1873  }
1874 
1875  }
1876 #endif
1877 }
1878 
1879 void
1880 DataManager::OutputEigGeometry(const unsigned uCurrEigSol, const int iResultsPrecision)
1881 {
1882  NodeContainerType::const_iterator i = NodeData[Node::STRUCTURAL].NodeContainer.begin();
1883  NodeContainerType::const_iterator e = NodeData[Node::STRUCTURAL].NodeContainer.end();
1884 
1885  // no structural nodes!
1886  if (i == e) {
1887  return;
1888  }
1889 
1890  if (OutHdl.UseText(OutputHandler::EIGENANALYSIS)) {
1891  std::ostream& out = OutHdl.Eigenanalysis();
1892 
1893  if (iResultsPrecision) {
1894  // 7 = number of characters requested by scientific notation
1895  int iNewWidth = iResultsPrecision + 7;
1896  out.width(iNewWidth);
1897  out.precision(iResultsPrecision);
1898  }
1899 
1900  out
1901  << "% structural nodes labels" << std::endl
1902  << "labels = [" << std::endl;
1903 
1904  for (; i != e; ++i) {
1905  const StructDispNode *pN = dynamic_cast<const StructDispNode *>(i->second);
1906  const StructNode *pSN = dynamic_cast<const StructNode *>(pN);
1907  ASSERT(pN != 0);
1908 
1909  if (pSN && pSN->GetStructNodeType() == StructNode::DUMMY) {
1910  continue;
1911  }
1912 
1913  out << pN->GetLabel() << ";" << std::endl;
1914  }
1915 
1916  out << "];" << std::endl;
1917 
1918  out
1919  << "% structural nodes base index" << std::endl
1920  << "idx = [" << std::endl;
1921 
1922  for (i = NodeData[Node::STRUCTURAL].NodeContainer.begin(); i != e; ++i) {
1923  const StructDispNode *pN = dynamic_cast<const StructDispNode *>(i->second);
1924  const StructNode *pSN = dynamic_cast<const StructNode *>(pN);
1925  ASSERT(pN != 0);
1926 
1927  if (pSN && pSN->GetStructNodeType() == StructNode::DUMMY) {
1928  continue;
1929  }
1930 
1931  out << pN->iGetFirstIndex() << ";" << std::endl;
1932  }
1933 
1934  out << "];" << std::endl;
1935 
1936  out
1937  << "% structural nodes reference configuration (X, Phi)" << std::endl
1938  << "X0 = [" << std::endl;
1939 
1940  for (i = NodeData[Node::STRUCTURAL].NodeContainer.begin(); i != e; ++i) {
1941  const StructDispNode *pN = dynamic_cast<const StructDispNode *>(i->second);
1942  const StructNode *pSN = dynamic_cast<const StructNode *>(pN);
1943  ASSERT(pN != 0);
1944 
1945  if (pSN && pSN->GetStructNodeType() == StructNode::DUMMY) {
1946  continue;
1947  }
1948 
1949  const Vec3& X(pN->GetX());
1950  Vec3 Phi(mb_zero<Vec3>());
1951  if (pSN) {
1952  Phi = RotManip::VecRot(pSN->GetR());
1953  }
1954 
1955  out
1956  << X(1) << ";" << std::endl
1957  << X(2) << ";" << std::endl
1958  << X(3) << ";" << std::endl
1959  << Phi(1) << ";" << std::endl
1960  << Phi(2) << ";" << std::endl
1961  << Phi(3) << ";" << std::endl;
1962  }
1963 
1964  out << "];" << std::endl;
1965  }
1966 #ifdef USE_NETCDF
1967  if (OutHdl.UseNetCDF(OutputHandler::NETCDF)) {
1968 
1969  /* start corner and count vector for NetCDF matrix output.
1970  * Since we are writing a matrix element-by-element, count will
1971  * always be (1, 1) and start will move to the desired place in the
1972  * matrix */
1973  std::vector<long> start (2, 0);
1974  start[0] = uCurrEigSol;
1975  std::vector<long> count (2, 1);
1976 
1977  for (i = NodeData[Node::STRUCTURAL].NodeContainer.begin(); i != e; ++i) {
1978  const StructDispNode *pN = dynamic_cast<const StructDispNode *>(i->second);
1979  const StructNode *pSN = dynamic_cast<const StructNode *>(pN);
1980  integer iNodeIndex;
1981  ASSERT(pN != 0);
1982 
1983  if (pSN && pSN->GetStructNodeType() == StructNode::DUMMY) {
1984  continue;
1985  }
1986 
1987  iNodeIndex = pSN->iGetFirstIndex();
1988 
1989  Var_Eig_Idx->set_cur(&start[0]);
1990  Var_Eig_Idx->put(&iNodeIndex, &count[0]);
1991  start[1]++;
1992  }
1993 
1994  }
1995 #endif // USE_NETCDF
1996 }
1997 
1998 void
1999 DataManager::OutputEigenvectors(const VectorHandler *pBeta,
2000  const VectorHandler& R, const VectorHandler& I,
2001  const doublereal& dShiftR,
2002  const MatrixHandler *pVL, const MatrixHandler& VR,
2003  const std::vector<bool>& vOut,
2004  const unsigned uCurrEigSol,
2005  const int iResultsPrecision)
2006 {
2007  const char signs[] = {'-', '+'};
2008  int iSign;
2009 
2010  integer iSize = VR.iGetNumRows();
2011  integer iNVec = VR.iGetNumCols();
2012 
2013  integer iEigenValues = 0;
2014 
2015  for (integer r = 1; r <= iNVec; r++) {
2016  if (!vOut[r - 1]) {
2017  continue;
2018  }
2019  ++iEigenValues;
2020  }
2021 
2022  if (iEigenValues == 0)
2023  return; // this allows to load the .m file into Matlab/Octave even
2024  // if no eigenvalues have converged
2025 
2026  // alphar, alphai, beta
2027  if (OutHdl.UseText(OutputHandler::EIGENANALYSIS)) {
2028  std::ostream& out = OutHdl.Eigenanalysis();
2029 
2030  if (iResultsPrecision) {
2031  // 7 = number of characters requested by scientific notation
2032  int iNewWidth = iResultsPrecision + 7;
2033  out.width(iNewWidth);
2034  out.precision(iResultsPrecision);
2035  }
2036 
2037  out
2038  << "% alphar, alphai, beta" << std::endl
2039  << "alpha = [";
2040 
2041  for (integer r = 1; r <= iNVec; r++) {
2042  if (!vOut[r - 1]) {
2043  continue;
2044  }
2045 
2046  out
2047  << R(r) + dShiftR << ' '
2048  << I(r) << ' '
2049  << (pBeta ? (*pBeta)(r) : 1.)
2050  << ";" << std::endl;
2051  }
2052 
2053  out << "];" << std::endl;
2054 
2055  if (pVL) {
2056  // VL
2057  out
2058  << "% left eigenvectors" << std::endl
2059  << "VL = [" << std::endl;
2060  for (integer r = 1; r <= iSize; r++) {
2061  for (integer c = 1; c <= iNVec; c++) {
2062  if (!vOut[c - 1]) {
2063  continue;
2064  }
2065 
2066  if (I(c) != 0.) {
2067  ASSERTMSG(c < iNVec, "partial eigenanalysis output: complex eigenvalue with real part of left eigenvector only");
2068  ASSERT(I(c) > 0.);
2069 
2070  doublereal re = (*pVL)(r, c);
2071  // NOTE: we cannot assume that if c == iNVec
2072  // it corresponds to a real-valued eigenvalue;
2073  // "im" will be wrong, but at least we no not sigsegv
2074  doublereal im = (c < iNVec) ? (*pVL)(r, c + 1) : 0.;
2075  if (im < 0) {
2076  iSign = 0;
2077  im = -im;
2078  } else {
2079  iSign = 1;
2080  }
2081 
2082  out
2083  << re << signs[iSign] << "i*" << im << ' ';
2084  if (vOut[c]) {
2085  out
2086  << re << signs[1 - iSign] << "i*" << im << ' ';
2087  }
2088  c++;
2089  } else {
2090  out
2091  << (*pVL)(r, c) << ' ';
2092  }
2093  }
2094 
2095  if (r < iSize) {
2096  out << ";" << std::endl;
2097  } else {
2098  out << "];" << std::endl;
2099  }
2100  }
2101  }
2102 
2103  // VR
2104  out
2105  << "% right eigenvectors" << std::endl
2106  << "VR = [" << std::endl;
2107  for (integer r = 1; r <= iSize; r++) {
2108  for (integer c = 1; c <= iNVec; c++) {
2109  if (!vOut[c - 1]) {
2110  continue;
2111  }
2112 
2113  if(I(c) != 0.) {
2114  ASSERTMSG(c < iNVec, "partial eigenanalysis output: complex eigenvalue with real part of right eigenvector only");
2115  ASSERT(I(c) > 0.);
2116 
2117  doublereal re = VR(r, c);
2118  // NOTE: we cannote assume that if c == iNVec
2119  // it corresponds to a real-valued eigenvalue;
2120  // "im" will be wrong, but at least we do not sigsev
2121  doublereal im = (c < iNVec) ? VR(r, c + 1) : 0.;
2122  if (im < 0.) {
2123  iSign = 0;
2124  im = -im;
2125  } else {
2126  iSign = 1;
2127  }
2128  out
2129  << re << signs[iSign] << "i*" << im << ' ';
2130  if (vOut[c]) {
2131  out
2132  << re << signs[1 - iSign] << "i*" << im << ' ';
2133  }
2134  c++;
2135  } else {
2136  out
2137  << VR(r, c) << ' ';
2138  }
2139  }
2140 
2141  if (r < iSize) {
2142  out << ";" << std::endl;
2143  } else {
2144  out << "];" << std::endl;
2145  }
2146  }
2147  }
2148 
2149 #ifdef USE_NETCDF
2150  if (OutHdl.UseNetCDF(OutputHandler::NETCDF)) {
2151  OutputHandler::NcDimVec dim_alpha(2);
2152 
2153  std::stringstream dimname_ss;
2154  dimname_ss << "eig_" << uCurrEigSol << "_iNVec_out";
2155 
2156  integer iNVecOut = 0;
2157  for (integer r = 1; r <= iNVec; r++)
2158  {
2159  if(vOut[r -1]){
2160  iNVecOut++;
2161  }
2162  }
2163 
2164  dim_alpha[0] = OutHdl.CreateDim(dimname_ss.str(), iNVecOut);
2165  dim_alpha[1] = OutHdl.DimV3();
2166 
2167  OutputHandler::AttrValVec attrs3(3);
2168  attrs3[0] = OutputHandler::AttrVal("units", "-");
2169  attrs3[1] = OutputHandler::AttrVal("type", "doublereal");
2170  attrs3[2] = OutputHandler::AttrVal("description", "alpha matrix");
2171 
2172  std::stringstream varname_ss;
2173  varname_ss << "eig." << uCurrEigSol << ".alpha";
2174  Var_Eig_dAlpha = OutHdl.CreateVar(varname_ss.str(), ncDouble, attrs3, dim_alpha);
2175 
2176  Vec3 v;
2177  unsigned uNRec = 0;
2178  for (integer r = 1; r <= iNVec; r++) {
2179  if (!vOut[r - 1]) {
2180  continue;
2181  }
2182 
2183  v(1) = R(r) + dShiftR;
2184  v(2) = I(r);
2185  v(3) = (pBeta ? (*pBeta)(r) : 1.);
2186  Var_Eig_dAlpha->put_rec(v.pGetVec(), uNRec);
2187  uNRec++;
2188  }
2189 
2190  OutputHandler::NcDimVec dim_v(3);
2191  dim_v[0] = m_Dim_Eig_iComplex;
2192  dim_v[1] = dim_alpha[0];
2193  dim_v[2] = m_Dim_Eig_iSize;
2194 
2195  /* start corner and count vector for NetCDF matrix output.
2196  * Since we are writing a matrix element-by-element, count will
2197  * always be (1, 1, 1) and start will move to the desired place in the
2198  * matrix */
2199  std::vector<long> start (3, 0);
2200  const std::vector<long> count (3, 1);
2201 
2202  if (pVL) {
2203  // VL
2204  OutputHandler::AttrValVec attrs3(3);
2205  attrs3[0] = OutputHandler::AttrVal("units", "-");
2206  attrs3[1] = OutputHandler::AttrVal("type", "doublereal");
2207  attrs3[2] = OutputHandler::AttrVal("description", "VL - Left eigenvectors matrix");
2208 
2209  std::stringstream varname_ss;
2210  varname_ss << "eig." << uCurrEigSol << ".VL";
2211  Var_Eig_dVL = OutHdl.CreateVar(varname_ss.str(), ncDouble, attrs3, dim_v);
2212 
2213  doublereal re;
2214  doublereal im;
2215 
2216 
2217  for (integer r = 1; r <= iSize; r++) {
2218  start[1] = 0;
2219  for (integer c = 1; c <= iNVec; c++) {
2220  if (!vOut[c - 1]) {
2221  continue;
2222  }
2223 
2224  start[2] = r - 1;
2225  if (I(c) != 0.) {
2226  ASSERTMSG(c < iNVec, "partial eigenanalysis output: complex eigenvalue with real part of left eigenvector only");
2227  ASSERT(I(c) > 0.);
2228 
2229  re = (*pVL)(r, c); // see above comment
2230  im = (c < iNVec) ? (*pVL)(r, c + 1) : 0.;
2231 
2232  // NetCDF indexing is zero-based!
2233  start[0] = 0; // real part in first "page" of VL
2234  Var_Eig_dVL->set_cur(&start[0]);
2235  Var_Eig_dVL->put(&re, &count[0]);
2236 
2237  start[0] = 1; // imaginary part in second "page"
2238  Var_Eig_dVL->set_cur(&start[0]);
2239  Var_Eig_dVL->put(&im, &count[0]);
2240 
2241  start[1]++;
2242 
2243  if (vOut[c]) {
2244  im = -im;
2245 
2246  start[0] = 0;
2247  Var_Eig_dVL->set_cur(&start[0]);
2248  Var_Eig_dVL->put(&re, &count[0]);
2249 
2250  start[0] = 1;
2251  Var_Eig_dVL->set_cur(&start[0]);
2252  Var_Eig_dVL->put(&im, &count[0]);
2253 
2254  start[1]++;
2255  }
2256  c++;
2257  } else {
2258  re = (*pVL)(r, c);
2259  im = 0.;
2260 
2261  start[0] = 0;
2262  Var_Eig_dVL->set_cur(&start[0]);
2263  Var_Eig_dVL->put(&re, &count[0]);
2264 
2265  start[0] = 1;
2266  Var_Eig_dVL->set_cur(&start[0]);
2267  Var_Eig_dVL->put(&im, &count[0]);
2268 
2269  start[1]++;
2270  }
2271  }
2272  }
2273  }
2274 
2275  // VR
2276  attrs3[0] = OutputHandler::AttrVal("units", "-");
2277  attrs3[1] = OutputHandler::AttrVal("type", "doublereal");
2278  attrs3[2] = OutputHandler::AttrVal("description", "VR - Left eigenvectors matrix");
2279 
2280  varname_ss.str("");
2281  varname_ss.clear();
2282  varname_ss << "eig." << uCurrEigSol << ".VR";
2283  Var_Eig_dVR = OutHdl.CreateVar(varname_ss.str(), ncDouble, attrs3, dim_v);
2284 
2285  doublereal re;
2286  doublereal im;
2287  for (integer r = 1; r <= iSize; r++) {
2288  start[1] = 0;
2289  for (integer c = 1; c <= iNVec; c++) {
2290  if (!vOut[c - 1]) {
2291  continue;
2292  }
2293 
2294  start[2] = r - 1;
2295  if (I(c) != 0.) {
2296  ASSERTMSG(c < iNVec, "partial eigenanalysis output: complex eigenvalue with real part of right eigenvector only");
2297  ASSERT(I(c) > 0.);
2298 
2299  re = VR(r, c); // see above comments
2300  im = (c < iNVec) ? VR(r, c + 1) : 0.;
2301 
2302  // NetCDF indexing is zero-based!
2303  start[0] = 0; // real part in first "page" of VL
2304  Var_Eig_dVR->set_cur(&start[0]);
2305  Var_Eig_dVR->put(&re, &count[0]);
2306 
2307  start[0] = 1; // imaginary part in second "page"
2308  Var_Eig_dVR->set_cur(&start[0]);
2309  Var_Eig_dVR->put(&im, &count[0]);
2310 
2311  start[1]++;
2312 
2313  if (vOut[c]) {
2314  im = -im;
2315 
2316  start[0] = 0;
2317  Var_Eig_dVR->set_cur(&start[0]);
2318  Var_Eig_dVR->put(&re, &count[0]);
2319 
2320  start[0] = 1;
2321  Var_Eig_dVR->set_cur(&start[0]);
2322  Var_Eig_dVR->put(&im, &count[0]);
2323  start[1]++;
2324  }
2325  c++;
2326  } else {
2327  re = VR(r, c);
2328  im = 0.;
2329 
2330  start[0] = 0;
2331  Var_Eig_dVR->set_cur(&start[0]);
2332  Var_Eig_dVR->put(&re, &count[0]);
2333 
2334  start[0] = 1;
2335  Var_Eig_dVR->set_cur(&start[0]);
2336  Var_Eig_dVR->put(&im, &count[0]);
2337 
2338  start[1]++;
2339  }
2340  }
2341  }
2342 
2343  }
2344 #endif /* USE_NETCDF */
2345 }
2346 
2347 bool
2348 DataManager::OutputEigClose(void)
2349 {
2350  return OutHdl.Close(OutputHandler::EIGENANALYSIS);
2351 }
2352 
2353 /* Output dati */
2354 bool
2355 DataManager::Output(long lStep,
2356  const doublereal& dTime,
2357  const doublereal& dTimeStep,
2358  bool force) const
2359 {
2360  /* Nota: il casting di OutHdl e' necessario in quanto la funzione propria
2361  * <void DataManager::Output(void) const> e' dichiarata, appunto, <const>.
2362  * Questo fa si' che un oggetto proprio della classe DataManager sia
2363  * implicitamente definito come <const> agli occhi della funzione.
2364  * Dal momento che le funzioni
2365  * <void NodeManager::Output(OutputHandler&) const> e
2366  * <void ElemManager::Output(OutputHandler&) const> ricevono come argomento
2367  * un oggetto di tipo <OutputHandler&> che non e' <const> in quanto su di
2368  * esso si scrive, il casting e' necessario per spiegare alla funzione
2369  * <void DataManager::Output(void) const> che le funzioni invocate
2370  * modificano si' l'<OutputHandler> passato loro, ma solo nel modo
2371  * consentito e quindi la sua dichiarazione come funzione <const> e'
2372  * dovuta al fatto che i dati propri non vengono modificati in modo
2373  * incontrollabile */
2374 
2375  DriveTrace(OutHdl); // trace output will be written for every time step
2376 
2377  /* output only when allowed by the output meter */
2378  if (!force && !pOutputMeter->dGet()) {
2379  return false;
2380  }
2381 
2382  /*
2383  * Write general simulation data to binary NetCDF file
2384  * the current time step index
2385  * the current simulatin time
2386  * the current integration time step
2387  */
2388 #ifdef USE_NETCDF
2389  if (OutHdl.UseNetCDF(OutputHandler::NETCDF)) {
2390  Var_Step->put_rec(&lStep, OutHdl.GetCurrentStep());
2391  Var_Time->put_rec(&dTime, OutHdl.GetCurrentStep());
2392  Var_TimeStep->put_rec(&dTimeStep, OutHdl.GetCurrentStep());
2393  }
2394 #endif /* USE_NETCDF */
2395 
2396  /* Dati dei nodi */
2397  NodeOutput(OutHdl);
2398 
2399  /* Dati degli elementi */
2400  ElemOutput(OutHdl);
2401 
2402  DriveOutput(OutHdl);
2403 
2404  OutHdl.IncCurrentStep();
2405 #ifdef USE_NETCDF
2406  if (bNetCDFsync) {
2407  OutHdl.pGetBinFile()->sync();
2408  }
2409 #endif /* USE_NETCDF */
2410 
2411  return true;
2412 }
2413 
2414 void
2415 DataManager::DriveTrace(OutputHandler& OH) const
2416 {
2417  if (!OH.IsOpen(OutputHandler::TRACES)) {
2418  return;
2419  }
2420 
2421  const MBDynParser::DCType& DC = MBPar.GetDriveCallerContainer();
2422 
2423  for (MBDynParser::DCType::const_iterator i = DC.begin(); i != DC.end(); ++i) {
2424  i->second->Trace(OH);
2425  }
2426 }
2427 
2428 void
2429 DataManager::DriveOutput(OutputHandler& OH) const
2430 {
2431  if (!OH.IsOpen(OutputHandler::DRIVECALLERS)) {
2432  return;
2433  }
2434 
2435  const MBDynParser::DCType& DC = MBPar.GetDriveCallerContainer();
2436 
2437  for (MBDynParser::DCType::const_iterator i = DC.begin(); i != DC.end(); ++i) {
2438  i->second->Output(OH);
2439  }
2440 }
2441 
2442 /* Output dati */
2443 void
2444 DataManager::Output(const VectorHandler& X, const VectorHandler& XP) const
2445 {
2446  /* Dati dei nodi */
2447  NodeOutput(OutHdl, X, XP);
2448 
2449  /* Dati degli elementi */
2450  ElemOutput(OutHdl, X, XP);
2451 }
2452 
2453 void
2454 DataManager::BeforePredict(VectorHandler& X, VectorHandler& XP,
2455  VectorHandler& XPrev,
2456  VectorHandler& XPPrev) const
2457 {
2458  for (NodeVecType::const_iterator i = Nodes.begin(); i != Nodes.end(); ++i) {
2459  (*i)->BeforePredict(X, XP, XPrev, XPPrev);
2460  }
2461 
2462  /* Versione con iteratore: */
2463  Elem* pEl = NULL;
2464  if (ElemIter.bGetFirst(pEl)) {
2465  do {
2466  pEl->BeforePredict(X, XP, XPrev, XPPrev);
2467  } while (ElemIter.bGetNext(pEl));
2468  }
2469 }
2470 
2471 void
2472 DataManager::AfterPredict(void) const
2473 {
2474  /* reset any external convergence requirement before starting
2475  * a new step */
2476  for (Converged_t::iterator i = m_IsConverged.begin();
2477  i != m_IsConverged.end(); ++i)
2478  {
2479  *i = Converged::NOT_CONVERGED;
2480  }
2481 
2482  for (NodeVecType::const_iterator i = Nodes.begin(); i != Nodes.end(); ++i) {
2483  try {
2484  (*i)->AfterPredict(*pXCurr, *pXPrimeCurr);
2485  }
2486  catch (Elem::ChangedEquationStructure& e) {
2487  // ignore by now
2488  silent_cerr("DataManager::AfterPredict: "
2489  "warning, caught Elem::ChangedEquationStructure while processing "
2490  << psNodeNames[(*i)->GetNodeType()] << "(" << (*i)->GetLabel() << ")" << std::endl);
2491  }
2492  }
2493 
2494  Elem* pEl = 0;
2495  if (ElemIter.bGetFirst(pEl)) {
2496  do {
2497  try {
2498  pEl->AfterPredict(*pXCurr, *pXPrimeCurr);
2499  }
2500  catch (Elem::ChangedEquationStructure) {
2501  // ignore by now
2502  silent_cerr("DataManager::AfterPredict: "
2503  "warning, caught Elem::ChangedEquationStructure while processing "
2504  << psElemNames[pEl->GetElemType()] << "(" << pEl->GetLabel() << ")" << std::endl);
2505  }
2506  } while (ElemIter.bGetNext(pEl));
2507  }
2508 }
2509 
2510 void
2511 DataManager::Update(void) const
2512 {
2513  for (NodeVecType::const_iterator i = Nodes.begin(); i != Nodes.end(); ++i) {
2514  (*i)->Update(*pXCurr, *pXPrimeCurr);
2515  }
2516 
2517  /* Versione con iteratore: */
2518  Elem* pEl = NULL;
2519  if (ElemIter.bGetFirst(pEl)) {
2520  do {
2521  pEl->Update(*pXCurr, *pXPrimeCurr);
2522  } while (ElemIter.bGetNext(pEl));
2523  }
2524 }
2525 
2526 void
2527 DataManager::AfterConvergence(void) const
2528 {
2529  for (NodeVecType::const_iterator i = Nodes.begin(); i != Nodes.end(); ++i) {
2530  (*i)->AfterConvergence(*pXCurr, *pXPrimeCurr);
2531  }
2532 
2533  /* Versione con iteratore: */
2534  Elem* pEl = NULL;
2535  if (ElemIter.bGetFirst(pEl)) {
2536  do {
2537  pEl->AfterConvergence(*pXCurr,
2538  *pXPrimeCurr);
2539  } while (ElemIter.bGetNext(pEl));
2540  }
2541 
2542  /* Restart condizionato */
2543  switch (RestartEvery) {
2544  case NEVER:
2545  break;
2546 
2547  case ITERATIONS:
2548  if (++iCurrRestartIter == iRestartIterations) {
2549  iCurrRestartIter = 0;
2550  const_cast<DataManager *>(this)->MakeRestart();
2551  }
2552  break;
2553 
2554  case TIME: {
2555  doublereal dT = DrvHdl.dGetTime();
2556  if (dT - dLastRestartTime >= dRestartTime) {
2557  dLastRestartTime = dT;
2558  const_cast<DataManager *>(this)->MakeRestart();
2559  }
2560  break;
2561  }
2562 
2563  case TIMES: {
2564  doublereal dT = DrvHdl.dGetTime()
2565  + pSolver->dGetInitialTimeStep()/100.;
2566  if (iCurrRestartTime == iNumRestartTimes) {
2567  break;
2568  }
2569 
2570  ASSERT(iCurrRestartTime < iNumRestartTimes);
2571 
2572  if (dT >= pdRestartTimes[iCurrRestartTime]) {
2573  iCurrRestartTime++;
2574  const_cast<DataManager *>(this)->MakeRestart();
2575  }
2576  break;
2577  }
2578 
2579  default:
2580  ASSERT(0);
2581  break;
2582  }
2583 }
2584 
2585 
2586 void
2587 DataManager::DerivativesUpdate(void) const
2588 {
2589  for (NodeVecType::const_iterator i = Nodes.begin(); i != Nodes.end(); ++i) {
2590  (*i)->DerivativesUpdate(*pXCurr, *pXPrimeCurr);
2591  }
2592 
2593  /* Versione con iteratore: */
2594  Elem* pEl = NULL;
2595  if (ElemIter.bGetFirst(pEl)) {
2596  do {
2597  pEl->DerivativesUpdate(*pXCurr, *pXPrimeCurr);
2598  } while (ElemIter.bGetNext(pEl));
2599  }
2600 }
2601 
2602 void
2603 DataManager::PrintResidual(const VectorHandler& Res, integer iIterCnt) const
2604 {
2605  silent_cout("Residual(" << DrvHdl.iGetStep() << ":" << iIterCnt << ") "
2606  "t=" << DrvHdl.dGetTime()
2607  << " dt=" << DrvHdl.dGetTimeStep() << std::endl);
2608  integer iSize = Res.iGetSize();
2609  for (int iTmpCnt = 1; iTmpCnt <= iSize; iTmpCnt++) {
2610  silent_cout("Eq " << std::setw(8)
2611  << iTmpCnt << ": "
2612  << std::setw(20) << Res(iTmpCnt)
2613  << " " << Dofs[iTmpCnt - 1].EqDescription
2614  << std::endl);
2615  }
2616 }
2617 
2618 void
2619 DataManager::PrintSolution(const VectorHandler& Sol, integer iIterCnt) const
2620 {
2621  silent_cout("Solution(" << DrvHdl.iGetStep() << ":" << iIterCnt << ") "
2622  "t=" << DrvHdl.dGetTime()
2623  << " dt=" << DrvHdl.dGetTimeStep() << std::endl);
2624  integer iSize = Sol.iGetSize();
2625  for (integer iTmpCnt = 1; iTmpCnt <= iSize; iTmpCnt++) {
2626  silent_cout("Dof " << std::setw(8)
2627  << iTmpCnt << ": "
2628  << std::setw(20) << Sol(iTmpCnt)
2629  << " " << Dofs[iTmpCnt - 1].Description
2630  << std::endl);
2631  }
2632 }
2633 
2634 const std::string&
2635 DataManager::GetDofDescription(int i) const
2636 {
2637  ASSERT(i > 0 && i <= iTotDofs);
2638  return Dofs[i - 1].Description;
2639 }
2640 
2641 const std::string&
2642 DataManager::GetEqDescription(int i) const
2643 {
2644  ASSERT(i > 0 && i <= iTotDofs);
2645  return Dofs[i - 1].EqDescription;
2646 }
2647 
2648 DofOrder::Order
2649 DataManager::GetDofType(int i) const
2650 {
2651  ASSERT(i > 0 && i <= iTotDofs);
2652  return Dofs[i - 1].Order;
2653 }
2654 
2655 DofOrder::Order
2656 DataManager::GetEqType(int i) const
2657 {
2658  ASSERT(i > 0 && i <= iTotDofs);
2659  return Dofs[i - 1].EqOrder;
2660 }
2661 
2662 #ifdef MBDYN_FDJAC
2663 bool
2664 DataManager::bFDJac(void) const
2665 {
2666  if (pFDJacMeter) {
2667  return (pFDJacMeter->dGet() != 0.);
2668  }
2669 
2670  return true;
2671 }
2672 #endif // MBDYN_FDJAC
2673 
2674 unsigned
2675 DataManager::ConvergedRegister(void)
2676 {
2677  unsigned idx = m_IsConverged.size();
2678  m_IsConverged.resize(idx + 1);
2679  m_IsConverged[idx] = Converged::CONVERGED;
2680  return idx;
2681 }
2682 
2683 void
2684 DataManager::ConvergedSet(unsigned idx, Converged::State s)
2685 {
2686  ASSERT(idx < m_IsConverged.size());
2687 
2688  m_IsConverged[idx] = s;
2689 }
2690 
2691 bool
2692 DataManager::IsConverged(void) const
2693 {
2694  for (Converged_t::const_iterator i = m_IsConverged.begin();
2695  i != m_IsConverged.end(); ++i)
2696  {
2697  if (*i == Converged::NOT_CONVERGED) {
2698  return false;
2699  }
2700  }
2701 
2702  return true;
2703 }
2704 
2705 bool
2706 DataManager::EndOfSimulation(void) const
2707 {
2708  for (Converged_t::const_iterator i = m_IsConverged.begin();
2709  i != m_IsConverged.end(); ++i)
2710  {
2711  if (*i == Converged::END_OF_SIMULATION) {
2712  return true;
2713  }
2714  }
2715 
2716  return false;
2717 }
2718 
2719 std::vector<doublereal>&
2720 DataManager::GetBufIn(unsigned uL)
2721 {
2722  Drive* pD = pFindDrive(Drive::FILEDRIVE, uL);
2723  if (pD == 0) {
2724  silent_cerr("unable to find FileDrive(" << uL << ")" << std::endl);
2725  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2726  }
2727 
2728  BufferStreamDrive *pBSD = dynamic_cast<BufferStreamDrive *>(pD);
2729  if (pBSD == 0) {
2730  silent_cerr("FileDrive(" << uL << ") is not a BufferStreamDrive" << std::endl);
2731  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2732  }
2733 
2734  return pBSD->GetBuf();
2735 }
2736 
2737 const std::vector<doublereal>&
2738 DataManager::GetBufOut(unsigned uL) const
2739 {
2740  BufferStreamElem *pBSE = pFindElem<BufferStreamElem, StreamOutElem, Elem::SOCKETSTREAM_OUTPUT>(uL);
2741  if (pBSE == 0) {
2742  silent_cerr("unable to find StreamOutElem(" << uL << "), or not a BufferStreamElem" << std::endl);
2743  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2744  }
2745 
2746  return pBSE->GetBuf();
2747 }
2748 
2749 doublereal *
2750 DataManager::GetBufInRaw(unsigned uL)
2751 {
2752  Drive* pD = pFindDrive(Drive::FILEDRIVE, uL);
2753  if (pD == 0) {
2754  silent_cerr("unable to find FileDrive(" << uL << ")" << std::endl);
2755  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2756  }
2757 
2758  BufferStreamDrive_base *pBSD = dynamic_cast<BufferStreamDrive_base *>(pD);
2759  if (pBSD == 0) {
2760  silent_cerr("FileDrive(" << uL << ") is not a BufferStreamDrive_base" << std::endl);
2761  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2762  }
2763 
2764  // NOTE: this cast is needed because the input buf must be writeable.
2765  return (doublereal *)pBSD->GetBufRaw();
2766 }
2767 
2768 void
2769 DataManager::SetBufInRaw(unsigned uL, integer n, const doublereal *p)
2770 {
2771  Drive* pD = pFindDrive(Drive::FILEDRIVE, uL);
2772  if (pD == 0) {
2773  silent_cerr("unable to find FileDrive(" << uL << ")" << std::endl);
2774  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2775  }
2776 
2777  BufferStreamDriveRaw *pBSD = dynamic_cast<BufferStreamDriveRaw *>(pD);
2778  if (pBSD == 0) {
2779  silent_cerr("FileDrive(" << uL << ") is not a BufferStreamDriveRaw" << std::endl);
2780  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2781  }
2782 
2783  if (pBSD->bOwnsMemory()) {
2784  silent_cerr("FileDrive(" << uL << ") owns its memory, unable to set buffer" << std::endl);
2785  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2786  }
2787 
2788  pBSD->SetBufRaw(n, p);
2789 }
2790 
2791 const doublereal *
2792 DataManager::GetBufOutRaw(unsigned uL) const
2793 {
2794  BufferStreamElem_base *pBSE = pFindElem<BufferStreamElem_base, StreamOutElem, Elem::SOCKETSTREAM_OUTPUT>(uL);
2795  if (pBSE == 0) {
2796  silent_cerr("unable to find StreamOutElem(" << uL << "), or not a BufferStreamElem_base" << std::endl);
2797  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2798  }
2799 
2800  return pBSE->GetBufRaw();
2801 }
2802 
2803 void
2804 DataManager::SetBufOutRaw(unsigned uL, integer n, const doublereal *p)
2805 {
2806  BufferStreamElemRaw *pBSE = pFindElem<BufferStreamElemRaw, StreamOutElem, Elem::SOCKETSTREAM_OUTPUT>(uL);
2807  if (pBSE == 0) {
2808  silent_cerr("unable to find StreamOutElem(" << uL << "), or not a BufferStreamElemRaw" << std::endl);
2809  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2810  }
2811 
2812  if (pBSE->bOwnsMemory()) {
2813  silent_cerr("StreamOutElem(" << uL << ") owns its memory, unable to set buffer" << std::endl);
2814  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
2815  }
2816 
2817  pBSE->SetBufRaw(n, p);
2818 }
2819 
2820 /* DataManager - end */
SimulationEntity::GetEqType
virtual DofOrder::Order GetEqType(unsigned int i) const
Definition: simentity.h:138
DataManager::IsConverged
bool IsConverged(void) const
Definition: dataman2.cc:2692
StructNode::GetWPrev
virtual const Vec3 & GetWPrev(void) const
Definition: strnode.h:1018
DataManager::RestartEvery
eRestart RestartEvery
Definition: dataman.h:164
VectorHandler::Reset
virtual void Reset(void)=0
SolutionManager::pResHdl
virtual VectorHandler * pResHdl(void) const =0
DEBUG_LEVEL_MATCH
#define DEBUG_LEVEL_MATCH(level)
Definition: myassert.h:241
DataManager::pFindDrive
Drive * pFindDrive(Drive::Type Typ, unsigned int uL) const
Definition: elman.cc:705
calls
void * calls
Definition: module-template.cc:425
DataManager::ElemDataStructure::ElemContainer
ElemContainerType ElemContainer
Definition: dataman.h:591
DataManager::OutputEigParams
void OutputEigParams(const doublereal &dTime, const doublereal &dCoef, const unsigned uCurrEigSol, const int iResultsPrecision)
Definition: dataman2.cc:1525
DataManager::GetBufIn
std::vector< doublereal > & GetBufIn(unsigned uL)
Definition: dataman2.cc:2720
MatrixHandler::iGetNumCols
virtual integer iGetNumCols(void) const =0
CGR_Rot::Param
const Param_Manip Param
Definition: matvec3.cc:644
SolverDiagnostics::outputIters
bool outputIters(void) const
Definition: solverdiagnostics.h:84
DataManager::dEpsilon
doublereal dEpsilon
Definition: dataman.h:126
StructDispNode::GetStructDispNodeType
virtual StructDispNode::Type GetStructDispNodeType(void) const =0
SubjectToInitialAssembly::iGetInitialNumDof
virtual unsigned int iGetInitialNumDof(void) const =0
MBDYN_EXCEPT_ARGS
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
DataManager::GetBufOutRaw
const doublereal * GetBufOutRaw(unsigned uL) const
Definition: dataman2.cc:2792
DataManager::NodeOutput
void NodeOutput(OutputHandler &OH) const
Definition: nodeman.cc:152
Vec3
Definition: matvec3.h:98
OutputHandler::DofStats
std::ostream & DofStats(void) const
Definition: output.h:594
Node
Definition: node.h:67
DataManager::MathPar
MathParser & MathPar
Definition: dataman.h:98
DataManager::ElemDataStructure::DofOwnerType
DofOwner::Type DofOwnerType
Definition: dataman.h:565
VecIter::bGetNext
bool bGetNext(T &TReturn) const
Definition: veciter.h:118
DEBUGCOUTFNAME
#define DEBUGCOUTFNAME(fname)
Definition: myassert.h:256
VectorHandler::pdGetVec
virtual doublereal * pdGetVec(void) const =0
DataManager::dGetInitialVelocityStiffness
const doublereal & dGetInitialVelocityStiffness(void) const
Definition: dataman2.cc:117
StructDispNode::DescribeDof
virtual std::ostream & DescribeDof(std::ostream &out, const char *prefix="", bool bInitial=false) const
Definition: strnode.cc:200
DofOwner::iNumDofs
unsigned int iNumDofs
Definition: dofown.h:95
OutputHandler::UseNetCDF
bool UseNetCDF(int out) const
Definition: output.cc:491
MBDynParser::GetDriveCallerContainer
const DCType & GetDriveCallerContainer(void) const
Definition: mbpar.cc:854
StructNode::GetRCurr
virtual const Mat3x3 & GetRCurr(void) const
Definition: strnode.h:1012
StructDispNode::GetNodeType
virtual Node::Type GetNodeType(void) const
Definition: strnode.cc:145
DataManager::Dofs
DofVecType Dofs
Definition: dataman.h:813
DataManager::ConvergedRegister
unsigned ConvergedRegister(void)
Definition: dataman2.cc:2675
DataManager::GetEqType
virtual DofOrder::Order GetEqType(int i) const
Definition: dataman2.cc:2656
Drive
Definition: drive.h:89
DataManager::DriveData
struct DataManager::@29 DriveData[Drive::LASTDRIVETYPE]
DataManager::OutputEigClose
bool OutputEigClose(void)
Definition: dataman2.cc:2348
DataManager::Nodes
NodeVecType Nodes
Definition: dataman.h:743
DataManager::dInitialPositionStiffness
doublereal dInitialPositionStiffness
Definition: dataman.h:121
StructDispNode::GetX
const Vec3 & GetX(void) const
Definition: strnode.cc:158
NaiveMatrixHandler::begin
NaiveMatrixHandler::const_iterator begin(void) const
Definition: naivemh.h:97
DataManager::EndOfSimulation
bool EndOfSimulation(void) const
Definition: dataman2.cc:2706
DataManager::bOmegaRotates
bool bOmegaRotates
Definition: dataman.h:123
StructDispNode::GetXPrev
virtual const Vec3 & GetXPrev(void) const
Definition: strnode.h:304
DataManager::dGetTime
doublereal dGetTime(void) const
Definition: dataman2.cc:165
DataManager::iCurrRestartTime
integer iCurrRestartTime
Definition: dataman.h:170
DataManager::OutHdl
OutputHandler OutHdl
Definition: dataman.h:105
DataManager::GetDofDescription
virtual const std::string & GetDofDescription(int i) const
Definition: dataman2.cc:2635
LinSol::GetSolutionManager
SolutionManager *const GetSolutionManager(integer iNLD, integer iLWS=0) const
Definition: linsol.cc:455
DataManager::Output
virtual bool Output(long lStep, const doublereal &dTime, const doublereal &dTimeStep, bool force=false) const
Definition: dataman2.cc:2355
DofOwnerOwner::SetInitialValue
virtual void SetInitialValue(VectorHandler &X)
Definition: dofown.cc:56
DataManager::BeforePredict
virtual void BeforePredict(VectorHandler &X, VectorHandler &XP, VectorHandler &XPrev, VectorHandler &XPPrev) const
Definition: dataman2.cc:2454
DataManager::m_IsConverged
Converged_t m_IsConverged
Definition: dataman.h:399
OutputHandler::Close
bool Close(const OutputHandler::OutFiles out)
Definition: output.cc:536
SAFEDELETEARR
#define SAFEDELETEARR(pnt)
Definition: mynewmem.h:713
mb_zero< Vec3 >
const Vec3 & mb_zero< Vec3 >(void)
Definition: matvec3.h:1560
ASSERTMSG
#define ASSERTMSG(expr, msg)
Definition: myassert.h:219
NaiveMatrixHandler::end
const NaiveMatrixHandler::const_iterator & end(void) const
Definition: naivemh.h:101
DataManager::ConvergedSet
void ConvergedSet(unsigned idx, Converged::State s)
Definition: dataman2.cc:2684
DataManager::DriveTrace
void DriveTrace(OutputHandler &OH) const
Definition: dataman2.cc:2415
DataManager::DofOwnerSet
void DofOwnerSet(void)
Definition: dataman2.cc:1368
SolverDiagnostics::outputRes
bool outputRes(void) const
Definition: solverdiagnostics.h:88
DataManager::GetBufInRaw
doublereal * GetBufInRaw(unsigned uL)
Definition: dataman2.cc:2750
SolverDiagnostics::outputJac
bool outputJac(void) const
Definition: solverdiagnostics.h:96
SparseMatrixHandler::Nz
const integer Nz() const
Definition: spmh.h:104
StructDispNode::dGetVelocityStiffness
virtual const doublereal & dGetVelocityStiffness(void) const
Definition: strnode.h:427
DataManager::iMaxInitialIterations
integer iMaxInitialIterations
Definition: dataman.h:125
SubjectToInitialAssembly::InitialWorkSpaceDim
virtual void InitialWorkSpaceDim(integer *piNumRows, integer *piNumCols) const =0
DataManager::SetLoadableElemModule
void SetLoadableElemModule(std::string, const LoadableCalls *, ModuleInsertMode=MIM_FAIL)
Definition: dataman2.cc:73
Elem::GetElemType
virtual Elem::Type GetElemType(void) const =0
DataManager::OutputEigFullMatrices
void OutputEigFullMatrices(const MatrixHandler *pmMatA, const MatrixHandler *pmMatB, const unsigned uCurrEigSol, const int iMatrixPrecision)
Definition: dataman2.cc:1565
DofOwner::iFirstIndex
integer iFirstIndex
Definition: dofown.h:94
DataManager::OutputEigGeometry
void OutputEigGeometry(const unsigned uCurrSol, const int iResultsPrecision)
Definition: dataman2.cc:1880
DataManager::DerivativesUpdate
virtual void DerivativesUpdate(void) const
Definition: dataman2.cc:2587
DataManager::MapOfLoadableElemHandlers
std::map< std::string, const LoadableCalls * > MapOfLoadableElemHandlers
Definition: dataman.h:102
VectorHandler::Dot
virtual doublereal Dot(void) const
Definition: vh.cc:244
DataManager::GetBufOut
const std::vector< doublereal > & GetBufOut(unsigned uL) const
Definition: dataman2.cc:2738
DataManager::ElemOutput
void ElemOutput(OutputHandler &OH) const
Definition: elman.cc:583
DataManager::DofOwnerInit
void DofOwnerInit(void)
Definition: dataman2.cc:182
DataManager::InitialJointAssembly
void InitialJointAssembly(void)
Definition: dataman2.cc:679
VectorHandler::iGetSize
virtual integer iGetSize(void) const =0
BufferStreamDriveRaw::bOwnsMemory
bool bOwnsMemory(void) const
Definition: bufferstreamdrive.cc:172
DataManager::GetDofType
virtual DofOrder::Order GetDofType(int i) const
Definition: dataman2.cc:2649
DataManager::GetLoadableElemModule
const LoadableCalls * GetLoadableElemModule(std::string) const
Definition: dataman2.cc:58
StructNode::bOmegaRotates
virtual bool bOmegaRotates(void) const
Definition: strnode.h:1236
SpMapMatrixHandler::end
const SpMapMatrixHandler::const_iterator & end(void) const
Definition: spmapmh.h:120
SolverDiagnostics::outputSolverConditionNumber
bool outputSolverConditionNumber(void) const
Definition: solverdiagnostics.h:117
DataManager::NodeDataStructure::NodeContainer
NodeContainerType NodeContainer
Definition: dataman.h:731
DataManager::AfterConvergence
virtual void AfterConvergence(void) const
Definition: dataman2.cc:2527
RotManip::VecRot
Vec3 VecRot(const Mat3x3 &Phi)
Definition: Rot.cc:136
DataManager::iTotDofs
integer iTotDofs
Definition: dataman.h:809
OutputHandler::Open
bool Open(const OutputHandler::OutFiles out)
Definition: output.cc:298
Mat3x3::MulTV
Vec3 MulTV(const Vec3 &v) const
Definition: matvec3.cc:482
BufferStreamDriveRaw::SetBufRaw
void SetBufRaw(integer n, const doublereal *p)
Definition: bufferstreamdrive.cc:178
DofOwnerOwner::pGetDofOwner
virtual const DofOwner * pGetDofOwner(void) const
Definition: dofown.h:113
DataManager::AfterPredict
virtual void AfterPredict(void) const
Definition: dataman2.cc:2472
DataManager::OutputEigPrepare
virtual void OutputEigPrepare(const integer iNumAnalyses, const integer iSize)
Definition: dataman2.cc:1474
DataManager::SetValue
void SetValue(VectorHandler &X, VectorHandler &XP)
Definition: dataman2.cc:1404
SolverDiagnostics::outputSol
bool outputSol(void) const
Definition: solverdiagnostics.h:92
DataManager::PrintSolution
virtual void PrintSolution(const VectorHandler &Sol, integer iIterCnt) const
Definition: dataman2.cc:2619
DataManager::LinkToSolution
void LinkToSolution(VectorHandler &XCurr, VectorHandler &XPrimeCurr)
Definition: dataman2.cc:172
SpMapMatrixHandler::const_iterator
Definition: spmapmh.h:91
DataManager::dRestartTime
doublereal dRestartTime
Definition: dataman.h:166
MBDynParser::DCType
std::map< unsigned, const DriveCaller * > DCType
Definition: mbpar.h:207
StructNode::GetStructNodeType
virtual StructNode::Type GetStructNodeType(void) const =0
DataManager::begin
DataManager::ElemContainerType::const_iterator begin(Elem::Type t) const
Definition: dataman.cc:902
DataManager::IncElemCount
void IncElemCount(Elem::Type type)
Definition: dataman2.cc:129
DriveHandler::iGetStep
integer iGetStep(void) const
Definition: drive.h:400
MyVectorHandler::Dot
doublereal Dot(void) const
Definition: vh.cc:713
DataManager::pSolver
Solver * pSolver
Definition: dataman.h:99
SolutionManager::pMatHdl
virtual MatrixHandler * pMatHdl(void) const =0
SubjectToInitialAssembly::InitialAssJac
virtual VariableSubMatrixHandler & InitialAssJac(VariableSubMatrixHandler &WorkMat, const VectorHandler &XCurr)=0
DataManager::pOutputMeter
DriveCaller * pOutputMeter
Definition: dataman.h:179
BufferStreamElemRaw::bOwnsMemory
bool bOwnsMemory(void) const
Definition: bufferstream_out_elem.cc:186
DataManager::MakeRestart
virtual void MakeRestart(void)
Definition: dataman.cc:699
OutputHandler::IncCurrentStep
void IncCurrentStep(void)
Definition: output.h:113
DataManager::dLastRestartTime
doublereal dLastRestartTime
Definition: dataman.h:173
DataManager::pRBK
RigidBodyKinematics * pRBK
Definition: dataman.h:129
StructNode::GetRPrev
virtual const Mat3x3 & GetRPrev(void) const
Definition: strnode.h:1000
RigidBodyKinematics::Update
virtual void Update(void)
Definition: rbk.cc:45
OutputHandler::GetCurrentStep
long GetCurrentStep(void) const
Definition: output.h:116
DriveHandler::dGetTimeStep
doublereal dGetTimeStep(void) const
Definition: drive.h:393
DataManager::OutputEigSparseMatrices
void OutputEigSparseMatrices(const MatrixHandler *pmMatA, const MatrixHandler *pmMatB, const unsigned uCurrEigSol, const int iMatrixPrecision)
Definition: dataman2.cc:1651
DataManager::NodeData
struct DataManager::NodeDataStructure NodeData[Node::LASTNODETYPE]
LinearSolver::ErrFactor::iCol
integer iCol
Definition: ls.h:60
DataManager::bDoesOmegaRotate
bool bDoesOmegaRotate(void) const
Definition: dataman2.cc:123
StructDispNode::GetVPrev
virtual const Vec3 & GetVPrev(void) const
Definition: strnode.h:316
count
static int count
Definition: modules.cc:41
DataManager::OutputPrepare
virtual void OutputPrepare(void)
Definition: dataman2.cc:1451
SubjectToInitialAssembly::InitialAssRes
virtual SubVectorHandler & InitialAssRes(SubVectorHandler &WorkVec, const VectorHandler &XCurr)=0
StructNode::GetR
const Mat3x3 & GetR(void) const
Definition: strnode.cc:1481
StructDispNode::dGetPositionStiffness
virtual const doublereal & dGetPositionStiffness(void) const
Definition: strnode.h:421
SpMapMatrixHandler::begin
SpMapMatrixHandler::const_iterator begin(void) const
Definition: spmapmh.h:116
MatrixHandler::PutCoef
virtual void PutCoef(integer iRow, integer iCol, const doublereal &dCoef)
Definition: mh.cc:384
Elem::DescribeEq
virtual std::ostream & DescribeEq(std::ostream &out, const char *prefix="", bool bInitial=false) const
Definition: elem.cc:137
DataManager::pXPrimeCurr
VectorHandler * pXPrimeCurr
Definition: dataman.h:109
StructDispNode::iGetFirstMomentumIndex
virtual integer iGetFirstMomentumIndex(void) const =0
StructDispNode::iGetFirstPositionIndex
virtual integer iGetFirstPositionIndex(void) const
Definition: strnode.h:452
StructNode::GetWCurr
virtual const Vec3 & GetWCurr(void) const
Definition: strnode.h:1030
SolutionManager::Solve
virtual void Solve(void)=0
DataManager::dInitialAssemblyTol
doublereal dInitialAssemblyTol
Definition: dataman.h:124
DataManager::iRestartIterations
integer iRestartIterations
Definition: dataman.h:165
StructDispNode::iGetInitialNumDof
virtual unsigned int iGetInitialNumDof(void) const
Definition: strnode.h:298
OutputHandler::IsOpen
bool IsOpen(int out) const
Definition: output.cc:395
ASSERT
#define ASSERT(expression)
Definition: colamd.c:977
DataManager::ElemIter
VecIter< Elem * > ElemIter
Definition: dataman.h:612
grad::sqrt
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
SolutionManager::MatrInitialize
virtual void MatrInitialize(void)
Definition: solman.cc:72
DataManager::dInitialVelocityStiffness
doublereal dInitialVelocityStiffness
Definition: dataman.h:122
MyVectorHandler::Reset
virtual void Reset(void)
Definition: vh.cc:459
SimulationEntity::BeforePredict
virtual void BeforePredict(VectorHandler &, VectorHandler &, VectorHandler &, VectorHandler &) const
Definition: simentity.cc:82
StructDispNode::GetXCurr
virtual const Vec3 & GetXCurr(void) const
Definition: strnode.h:310
Elem::iGetNumDof
virtual unsigned int iGetNumDof(void) const
Definition: elem.cc:118
VectorHandler::Add
virtual void Add(integer iRow, const Vec3 &v)
Definition: vh.cc:63
bufferstream_out_elem.h
DataManager::iSize
integer iSize
Definition: dataman.h:788
DataManager::SetTime
void SetTime(const doublereal &dTime, const doublereal &dTimeStep=-1., const integer &iStep=-1, bool bServePending=true)
Definition: dataman2.cc:139
defaultMemoryManager
#define defaultMemoryManager
Definition: mynewmem.h:691
DataManager::iCurrRestartIter
integer iCurrRestartIter
Definition: dataman.h:172
DataManager::iNum
unsigned int iNum
Definition: dataman.h:618
InitialAssemblyIterator::GetNext
InitialAssemblyElem * GetNext(void) const
Definition: elman.cc:791
DataManager::solArrFileName
char * solArrFileName
Definition: dataman.h:176
FullMatrixHandler::iGetNumCols
virtual integer iGetNumCols(void) const
Definition: fullmh.h:229
BufferStreamDrive::GetBuf
std::vector< doublereal > & GetBuf(void)
Definition: bufferstreamdrive.cc:130
DriveCaller::dGet
virtual doublereal dGet(const doublereal &dVar) const =0
DataManager::ElemOutputPrepare
void ElemOutputPrepare(OutputHandler &OH)
Definition: elman.cc:545
BufferStreamElem_base
Definition: bufferstream_out_elem.h:43
DataManager::OutputEigenvectors
void OutputEigenvectors(const VectorHandler *pBeta, const VectorHandler &R, const VectorHandler &I, const doublereal &dShiftR, const MatrixHandler *pVL, const MatrixHandler &VR, const std::vector< bool > &vOut, const unsigned uCurrEigSol, const int iResultsPrecision)
Definition: dataman2.cc:1999
Elem::GetConnectedNodes
virtual void GetConnectedNodes(std::vector< const Node * > &connectedNodes) const
Definition: elem.h:243
c
static std::stack< cleanup * > c
Definition: cleanup.cc:59
FullMatrixHandler::pdGetMat
const doublereal * pdGetMat(void) const
Definition: fullmh.h:144
DataManager::PrintResidual
virtual void PrintResidual(const VectorHandler &Res, integer iIterCnt) const
Definition: dataman2.cc:2603
BufferStreamElem_base::GetBufRaw
virtual const doublereal * GetBufRaw(void) const =0
psNodeNames
const char * psNodeNames[]
Definition: enums.cc:372
DriveHandler::LinkToSolution
void LinkToSolution(const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: drive.cc:260
Elem
Definition: elem.h:75
Node::GetNodeType
virtual Node::Type GetNodeType(void) const =0
Elem::Type
Type
Definition: elem.h:91
DataManager::pdRestartTimes
doublereal * pdRestartTimes
Definition: dataman.h:168
psElemNames
const char * psElemNames[]
Definition: enums.cc:39
DataManager::dGetInitialPositionStiffness
const doublereal & dGetInitialPositionStiffness(void) const
Definition: dataman2.cc:111
SimulationEntity::AfterConvergence
virtual void AfterConvergence(const VectorHandler &X, const VectorHandler &XP)
Definition: simentity.cc:120
SolutionManager::bGetConditionNumber
bool bGetConditionNumber(doublereal &dCond) const
Definition: solman.cc:107
DataManager::SetBufOutRaw
void SetBufOutRaw(unsigned uL, integer n, const doublereal *p)
Definition: dataman2.cc:2804
Vec3::pGetVec
const doublereal * pGetVec(void) const
Definition: matvec3.h:192
DataManager::MBPar
MBDynParser & MBPar
Definition: dataman.h:97
Elem::GetDofType
virtual DofOrder::Order GetDofType(unsigned int) const
Definition: elem.cc:150
DriveHandler::dGetTime
doublereal dGetTime(void) const
Definition: drive.h:386
DataManager::CurrSolver
LinSol CurrSolver
Definition: dataman.h:127
SimulationEntity::DerivativesUpdate
virtual void DerivativesUpdate(const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: simentity.cc:105
DriveHandler::SetTime
void SetTime(const doublereal &dt, const doublereal &dts, const integer &s)
Definition: drive.cc:208
StructDispNode::GetVCurr
virtual const Vec3 & GetVCurr(void) const
Definition: strnode.h:322
DataManager::Update
virtual void Update(void) const
Definition: dataman2.cc:2511
DataManager::DriveOutput
void DriveOutput(OutputHandler &OH) const
Definition: dataman2.cc:2429
DataManager::iNumRestartTimes
integer iNumRestartTimes
Definition: dataman.h:169
OutputHandler::Log
std::ostream & Log(void) const
Definition: output.h:541
BufferStreamElemRaw::SetBufRaw
virtual void SetBufRaw(integer n, const doublereal *p)
Definition: bufferstream_out_elem.cc:192
DataManager::ModuleInsertMode
ModuleInsertMode
Definition: dataman.h:513
VecIter::bGetFirst
bool bGetFirst(T &TReturn) const
Definition: veciter.h:88
Mat3x3::MulMT
Mat3x3 MulMT(const Mat3x3 &m) const
Definition: matvec3.cc:444
InitialAssemblyIterator::GetFirst
InitialAssemblyElem * GetFirst(void) const
Definition: elman.cc:763
SimulationEntity::AfterPredict
virtual void AfterPredict(VectorHandler &X, VectorHandler &XP)
Definition: simentity.cc:91
DataManager::DrvHdl
DriveHandler DrvHdl
Definition: dataman.h:104
DataManager::SetBufInRaw
void SetBufInRaw(unsigned uL, integer n, const doublereal *p)
Definition: dataman2.cc:2769
SolutionManager::pSolHdl
virtual VectorHandler * pSolHdl(void) const =0
BufferStreamDrive_base::GetBufRaw
virtual const doublereal * GetBufRaw(void)=0
SimulationEntity::SetValue
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *h=0)
Definition: simentity.cc:63
DataManager::DofData
struct DataManager::@30 DofData[DofOwner::LASTDOFTYPE]
DofOwnerOwner::iGetFirstIndex
virtual integer iGetFirstIndex(void) const
Definition: dofown.h:127
doublereal
double doublereal
Definition: colamd.c:52
BufferStreamElem::GetBuf
const std::vector< doublereal > & GetBuf(void) const
Definition: bufferstream_out_elem.cc:145
dataman_.h
OutputHandler::Eigenanalysis
std::ostream & Eigenanalysis(void) const
Definition: output.h:615
integer
long int integer
Definition: colamd.c:51
DataManager::NodeOutputPrepare
void NodeOutputPrepare(OutputHandler &OH)
Definition: nodeman.cc:119
MathParser::GetSymbolTable
Table & GetSymbolTable(void) const
Definition: mathp.cc:1931
Dof
Definition: dofown.h:59
Node::GetNode
const Node * GetNode(void) const
Definition: node.h:111
WithLabel::GetLabel
unsigned int GetLabel(void) const
Definition: withlab.cc:62
DataManager::ElemData
struct DataManager::ElemDataStructure ElemData[Elem::LASTELEMTYPE]
DataManager::OutputEigNaiveMatrices
void OutputEigNaiveMatrices(const MatrixHandler *pmMatA, const MatrixHandler *pmMatB, const unsigned uCurrEigSol, const int iMatrixPrecision)
Definition: dataman2.cc:1754
Solver::dGetInitialTimeStep
virtual doublereal dGetInitialTimeStep(void) const
Definition: solver.h:417
DEBUGLCOUT
#define DEBUGLCOUT(level, msg)
Definition: myassert.h:244
DataManager::GetEqDescription
virtual const std::string & GetEqDescription(int i) const
Definition: dataman2.cc:2642
DataManager::uPrintFlags
unsigned uPrintFlags
Definition: dataman.h:156
FullMatrixHandler::iGetNumRows
virtual integer iGetNumRows(void) const
Definition: fullmh.h:225
Elem::DescribeDof
virtual std::ostream & DescribeDof(std::ostream &out, const char *prefix="", bool bInitial=false) const
Definition: elem.cc:124
MatrixHandler::iGetNumRows
virtual integer iGetNumRows(void) const =0
NaiveMatrixHandler::const_iterator
Definition: naivemh.h:72
R
Mat3x3 R
Definition: mbdyn_util_oct.cc:74
DataManager::pXCurr
VectorHandler * pXCurr
Definition: dataman.h:108
StructDispNode::DescribeEq
virtual std::ostream & DescribeEq(std::ostream &out, const char *prefix="", bool bInitial=false) const
Definition: strnode.cc:259
OutputHandler::UseText
bool UseText(int out) const
Definition: output.cc:446
BufferStreamDrive_base
Definition: bufferstreamdrive.h:41
SAFEDELETE
#define SAFEDELETE(pnt)
Definition: mynewmem.h:710
SimulationEntity::Update
virtual void Update(const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: simentity.cc:98