MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
strmappingext.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/struct/strmappingext.cc,v 1.36 2017/01/12 14:46:44 masarati Exp $ */
2 /*
3  * MBDyn (C) is a multibody analysis code.
4  * http://www.mbdyn.org
5  *
6  * Copyright (C) 2007-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 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
33 
34 #include "dataman.h"
35 #include "strmappingext.h"
36 #include "modalmappingext.h"
37 #include "Rot.hh"
38 
39 #include <fstream>
40 #include <cerrno>
41 #include <cstdlib>
42 #include <ctime>
43 
44 /* StructMappingExtForce - begin */
45 
46 /* Costruttore */
48  DataManager *pDM,
49  const StructNode *pRefNode,
50  bool bUseReferenceNodeForces,
51  bool bRotateReferenceNodeForces,
52  std::vector<const StructDispNode *>& nodes,
53  std::vector<Vec3>& offsets,
54  std::vector<unsigned>& labels,
56  std::vector<uint32_t>& mappedlabels,
57  bool bLabels,
58  bool bOutputAccelerations,
59  unsigned uRRot,
60  ExtFileHandlerBase *pEFH,
61  bool bSendAfterPredict,
62  int iCoupling,
63  flag fOut)
64 : Elem(uL, fOut),
65 ExtForce(uL, pDM, pEFH, bSendAfterPredict, iCoupling, fOut),
66 pRefNode(pRefNode),
67 bUseReferenceNodeForces(bUseReferenceNodeForces),
68 bRotateReferenceNodeForces(bRotateReferenceNodeForces),
69 F0(Zero3), M0(Zero3),
70 F1(Zero3), M1(Zero3),
71 F2(Zero3), M2(Zero3),
72 pH(pH),
73 m_uResSize(0),
74 uPoints(nodes.size()),
75 uMappedPoints(pH ? unsigned(pH->iGetNumRows())/3 : 0),
76 bLabels(bLabels),
77 bOutputAccelerations(bOutputAccelerations),
78 uRRot(uRRot),
79 m_qlabels(pH ? mappedlabels : labels),
80 m_x(3*uPoints),
81 m_xP(3*uPoints),
82 m_xPP(0),
83 m_q(3*uMappedPoints),
84 m_qP(3*uMappedPoints),
85 m_qPP(0),
86 m_f(3*uPoints),
87 m_p(3*uMappedPoints)
88 {
89  ASSERT(nodes.size() == offsets.size());
90  if (pH) {
91  ASSERT(3*uPoints == unsigned(pH->iGetNumCols()));
92  ASSERT(3*uMappedPoints == unsigned(pH->iGetNumRows()));
93  }
94 
95  if (pRefNode) {
96  switch (uRRot) {
97  case MBC_ROT_THETA:
98  case MBC_ROT_MAT:
99  case MBC_ROT_EULER_123:
100  break;
101 
102  default:
103  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
104  "invalid reference node rotation type " << uRRot << std::endl);
106  }
107 
108  if (bUseReferenceNodeForces) {
109  m_uResSize += 6;
110  }
111  }
112 
113  const StructDispNode *pNode = 0;
114  unsigned uNodes = 0;
115  std::vector<const StructDispNode *>::const_iterator p;
116  for (p = nodes.begin(); p != nodes.end(); ++p) {
117  if (*p != pNode) {
118  pNode = *p;
119  uNodes++;
120  }
121  }
122 
123  Nodes.resize(uNodes);
124  p = nodes.begin();
125  std::vector<const StructDispNode *>::const_iterator pPrev = p;
126  std::vector<NodeData>::iterator n = Nodes.begin();
127  while (true) {
128  ++p;
129  if ((p == nodes.end()) || (*p != *pPrev)) {
130  n->pNode = *pPrev;
131  n->Offsets.resize(p - pPrev);
132  for (std::vector<OffsetData>::iterator i = n->Offsets.begin(); i != n->Offsets.end(); ++i) {
133  i->Offset = ::Zero3;
134  }
135  n->F = Zero3;
136  n->M = Zero3;
137 
138  if (dynamic_cast<const StructNode *>(n->pNode) == 0) {
139  switch (n->Offsets.size()) {
140  case 1: // regular rotationless node
141  case 2: // special rotationless node for "membrane"
142  break;
143 
144  default:
145  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
146  "too many offsets for StructDispNode(" << n->pNode->GetLabel() << ")"
147  << std::endl);
149  }
150 
151  m_uResSize += 3;
152 
153  } else {
154  m_uResSize += 6;
155  }
156 
157  if (p == nodes.end()) {
158  break;
159  }
160 
161  ++n;
162  pPrev = p;
163  }
164  }
165 
166  unsigned uPts = 0;
167  n = Nodes.begin();
168  std::vector<Vec3>::const_iterator o = offsets.begin();
169  std::vector<uint32_t>::iterator l = labels.begin();
170  for (; o != offsets.end(); ++o, uPts++) {
171  if (uPts == n->Offsets.size()) {
172  ++n;
173  uPts = 0;
174 
175  if (dynamic_cast<const StructNode *>(n->pNode) == 0) {
176  for (std::vector<StructMappingExtForce::OffsetData>::const_iterator i = n->Offsets.begin();
177  i != n->Offsets.end(); i++)
178  {
179  if (!i->Offset.IsNull()) {
180  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
181  "offset #" << (i - n->Offsets.begin())
182  << " (" << i->Offset << ") "
183  "for StructDispNode(" << n->pNode->GetLabel() << ") is not zero"
184  << std::endl);
186  }
187  }
188  }
189  }
190 
191  // FIXME: pass labels
192  n->Offsets[uPts].uLabel = unsigned(-1);
193  n->Offsets[uPts].Offset = *o;
194  n->Offsets[uPts].F = Zero3;
195 
196  if (bLabels) {
197  n->Offsets[uPts].uLabel = *l;
198  ++l;
199  }
200  }
201 
202  if (bOutputAccelerations) {
203  for (unsigned i = 0; i < nodes.size(); i++) {
204  const DynamicStructDispNode *pDSN = dynamic_cast<const DynamicStructDispNode *>(Nodes[i].pNode);
205  if (pDSN == 0) {
206  silent_cerr("StructMappingExtForce"
207  "(" << GetLabel() << "): "
208  "StructNode(" << Nodes[i].pNode->GetLabel() << ") "
209  "is not dynamic"
210  << std::endl);
212  }
213 
214  const_cast<DynamicStructDispNode *>(pDSN)->ComputeAccelerations(true);
215  }
216 
217  m_xPP.resize(3*uPoints);
218  m_qPP.resize(3*uMappedPoints);
219 
220  }
221 }
222 
224 {
225  NO_OP;
226 }
227 
228 void
230 {
231  if (iCoupling == COUPLING_NONE) {
232  *piNumRows = 0;
233  *piNumCols = 0;
234 
235  } else {
236  *piNumRows = (pRefNode ? 6 : 0) + 6*Nodes.size();
237  *piNumCols = 1;
238  }
239 }
240 
241 bool
243 {
244  bool bResult = true;
245 
246  switch (pEFH->NegotiateRequest()) {
248  break;
249 
251  std::ostream *outfp = pEFH->GetOutStream();
252  if (outfp) {
253 
254 #ifdef USE_SOCKET
255  } else {
256  char buf[sizeof(uint32_t) + sizeof(uint32_t)];
257  uint32_t *uint32_ptr;
258 
259  uint32_ptr = (uint32_t *)&buf[0];
260  uint32_ptr[0] = MBC_NODAL | MBC_U_ROT_2_REF_NODE_ROT(uRRot);
261  if (pRefNode != 0) {
262  uint32_ptr[0] |= MBC_REF_NODE;
263  }
264 
265  if (bLabels) {
266  uint32_ptr[0] |= MBC_LABELS;
267  }
268 
269  if (bOutputAccelerations) {
270  uint32_ptr[0] |= MBC_ACCELS;
271  }
272 
273  uint32_ptr[1] = uPoints;
274 
275  ssize_t rc = send(pEFH->GetOutFileDes(),
276  (const void *)buf, sizeof(buf),
277  pEFH->GetSendFlags());
278  if (rc == -1) {
279  int save_errno = errno;
280  char *err_msg = strerror(save_errno);
281  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
282  "negotiation request send() failed "
283  "(" << save_errno << ": " << err_msg << ")"
284  << std::endl);
286 
287  } else if (rc != sizeof(buf)) {
288  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
289  "negotiation request send() failed "
290  "(sent " << rc << " of " << sizeof(buf) << " bytes)"
291  << std::endl);
293  }
294 #endif // USE_SOCKET
295  }
296  } break;
297 
299  unsigned uN;
300  unsigned uNodal;
301  bool bRef;
302  unsigned uRR;
303  unsigned uR;
304  bool bA;
305  bool bL;
306 
307  std::istream *infp = pEFH->GetInStream();
308  if (infp) {
309  // TODO: stream negotiation?
310 
311 #ifdef USE_SOCKET
312  } else {
313  char buf[sizeof(uint32_t) + sizeof(uint32_t)];
314  uint32_t *uint32_ptr;
315 
316  ssize_t rc = recv(pEFH->GetInFileDes(),
317  (void *)buf, sizeof(buf),
318  pEFH->GetRecvFlags());
319  if (rc == -1) {
320  int save_errno = errno;
321  char *err_msg = strerror(save_errno);
322  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
323  "negotiation response recv() failed "
324  "(" << save_errno << ": " << err_msg << ")"
325  << std::endl);
327 
328  } else if (rc != sizeof(buf)) {
329  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
330  "negotiation response recv() failed "
331  "(got " << rc << " of " << sizeof(buf) << " bytes)"
332  << std::endl);
334  }
335 
336  uint32_ptr = (uint32_t *)&buf[0];
337  uNodal = (uint32_ptr[0] & MBC_MODAL_NODAL_MASK);
338  bRef = (uint32_ptr[0] & MBC_REF_NODE);
339  uRR = (uint32_ptr[0] & MBC_REF_NODE_ROT_MASK);
340  uR = (uint32_ptr[0] & MBC_ROT_MASK);
341  bL = (uint32_ptr[0] & MBC_LABELS);
342  bA = (uint32_ptr[0] & MBC_ACCELS);
343 
344  uN = uint32_ptr[1];
345 #endif // USE_SOCKET
346  }
347 
348  if (uNodal != MBC_NODAL) {
349  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
350  "negotiation response failed: expecting MBC_NODAL "
351  "(=" << MBC_MODAL << "), got " << uNodal
352  << std::endl);
353  bResult = false;
354  }
355 
356  if ((pRefNode != 0 && !bRef) || (pRefNode == 0 && bRef)) {
357  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
358  "negotiation response failed: reference node configuration mismatch "
359  "(local=" << (pRefNode != 0 ? "yes" : "no") << ", remote=" << (bRef ? "yes" : "no") << ")"
360  << std::endl);
361  bResult = false;
362  }
363 
364  if (uR != MBC_ROT_NONE) {
365  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
366  "negotiation response failed: orientation output mismatch "
367  "(local=" << MBC_ROT_NONE << ", remote=" << uR << ")"
368  << std::endl);
369  bResult = false;
370  }
371 
372  if (uRR != MBC_U_ROT_2_REF_NODE_ROT(uRRot)) {
373  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
374  "negotiation response failed: reference node orientation output mismatch "
375  "(local=" << MBC_U_ROT_2_REF_NODE_ROT(uRRot) << ", remote=" << uR << ")"
376  << std::endl);
377  bResult = false;
378  }
379 
380  if (bL != bLabels) {
381  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
382  "negotiation response failed: labels output mismatch "
383  "(local=" << (bLabels ? "yes" : "no") << ", remote=" << (bL ? "yes" : "no") << ")"
384  << std::endl);
385  bResult = false;
386  }
387 
388  if (bA != bOutputAccelerations) {
389  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
390  "negotiation response failed: acceleration output mismatch "
391  "(local=" << (bOutputAccelerations ? "yes" : "no") << ", remote=" << (bA ? "yes" : "no") << ")"
392  << std::endl);
393  bResult = false;
394  }
395 
396  unsigned uMP = pH ? uMappedPoints : uPoints;
397  if (uN != uMP) {
398  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
399  "negotiation response failed: node number mismatch "
400  "(local=" << uMP << ", remote=" << uN << ")"
401  << std::endl);
402  bResult = false;
403  }
404  } break;
405 
406  default:
408  }
409 
410  return bResult;
411 }
412 
413 /*
414  * Send output to companion software
415  */
416 void
418 {
419  std::ostream *outfp = pEFH->GetOutStream();
420  if (outfp) {
421  SendToStream(*outfp, when);
422 
423  } else {
424  SendToFileDes(pEFH->GetOutFileDes(), when);
425  }
426 }
427 
428 void
430 {
431  if (pRefNode) {
432  const Vec3& xRef = pRefNode->GetXCurr();
433  const Mat3x3& RRef = pRefNode->GetRCurr();
434  const Vec3& xpRef = pRefNode->GetVCurr();
435  const Vec3& wRef = pRefNode->GetWCurr();
436  const Vec3& xppRef = pRefNode->GetXPPCurr();
437  const Vec3& wpRef = pRefNode->GetWPCurr();
438 
439  if (bLabels) {
440  outf
441  << pRefNode->GetLabel()
442  << " ";
443  }
444 
445  outf
446  << xRef
447  << " " << RRef
448  << " " << xpRef
449  << " " << wRef;
450 
451  if (bOutputAccelerations) {
452  outf
453  << " " << xppRef
454  << " " << wpRef;
455  }
456  outf << std::endl;
457 
458  for (unsigned i = 0; i < Nodes.size(); i++) {
459 #if 0
460  Vec3 f(Nodes[i]->GetRCurr()*Offsets[i]);
461  Vec3 x(Nodes[i]->GetXCurr() + f);
462  Vec3 Dx(x - xRef);
463  Mat3x3 DR(RRef.MulTM(Nodes[i]->GetRCurr()));
464  Vec3 v(Nodes[i]->GetVCurr() + Nodes[i]->GetWCurr().Cross(f));
465  Vec3 Dv(v - xpRef - wRef.Cross(Dx));
466  const Vec3& w(Nodes[i]->GetWCurr());
467 
468  // manipulate
469 
470  if (bLabels) {
471  outf
472  << Nodes[i]->GetLabel()
473  << " ";
474  }
475 
476  outf
477  << RRef.MulTV(Dx);
478 
479  switch (uRot) {
480  case MBC_ROT_NONE:
481  break;
482 
483  case MBC_ROT_MAT:
484  outf
485  << " " << DR;
486  break;
487 
488  case MBC_ROT_THETA:
489  outf
490  << " " << RotManip::VecRot(DR);
491  break;
492 
493  case MBC_ROT_EULER_123:
494  outf
495  << " " << MatR2EulerAngles123(DR)*dRaDegr;
496  break;
497  }
498 
499  outf
500  << " " << RRef.MulTV(Dv);
501 
502  if (uRot != MBC_ROT_NONE) {
503  outf
504  << " " << RRef.MulTV(w - wRef);
505  }
506 
507  if (bOutputAccelerations) {
508  const Vec3& xpp(Nodes[i]->GetXPPCurr());
509 
510  outf
511  << " " << RRef.MulTV(xpp - xppRef - wpRef.Cross(Dx)
512  - wRef.Cross(wRef.Cross(Dx) + Dv*2));
513  if (uRot != MBC_ROT_NONE) {
514  const Vec3& wp(Nodes[i]->GetWPCurr());
515 
516  outf
517  << " " << RRef.MulTV(wp - wpRef - wRef.Cross(w));
518  }
519  }
520  outf << std::endl;
521 #endif
522  }
523 
524  } else {
525  for (unsigned i = 0; i < Nodes.size(); i++) {
526 #if 0
527  /*
528  p = x + f
529  R = R
530  v = xp + w cross f
531  w = w
532  a = xpp + wp cross f + w cross w cross f
533  wp = wp
534  */
535 
536  // Optimization of the above formulas
537  const Mat3x3& R = Nodes[i]->GetRCurr();
538  Vec3 f = R*Offsets[i];
539  Vec3 x = Nodes[i]->GetXCurr() + f;
540  const Vec3& w = Nodes[i]->GetWCurr();
541  Vec3 wCrossf = w.Cross(f);
542  Vec3 v = Nodes[i]->GetVCurr() + wCrossf;
543 
544  if (bLabels) {
545  outf
546  << Nodes[i]->GetLabel()
547  << " ";
548  }
549 
550  outf
551  << x;
552 
553  switch (uRot) {
554  case MBC_ROT_NONE:
555  break;
556 
557  case MBC_ROT_MAT:
558  outf
559  << " " << R;
560  break;
561 
562  case MBC_ROT_THETA:
563  outf
564  << " " << RotManip::VecRot(R);
565  break;
566 
567  case MBC_ROT_EULER_123:
568  outf
569  << " " << MatR2EulerAngles123(R)*dRaDegr;
570  break;
571  }
572  outf
573  << " " << v;
574 
575  if (uRot != MBC_ROT_NONE) {
576  outf
577  << " " << w;
578  }
579 
580  if (bOutputAccelerations) {
581  const Vec3& wp = Nodes[i]->GetWPCurr();
582  Vec3 a = Nodes[i]->GetXPPCurr() + wp.Cross(f) + w.Cross(wCrossf);
583 
584  outf
585  << " " << a;
586 
587  if (uRot != MBC_ROT_NONE) {
588  outf
589  << " " << wp;
590  }
591  }
592 
593  outf << std::endl;
594 #endif
595  }
596  }
597 }
598 
599 void
601 {
602 #ifdef USE_SOCKET
603  if (pRefNode) {
604  const Vec3& xRef = pRefNode->GetXCurr();
605  const Mat3x3& RRef = pRefNode->GetRCurr();
606  const Vec3& xpRef = pRefNode->GetVCurr();
607  const Vec3& wRef = pRefNode->GetWCurr();
608  const Vec3& xppRef = pRefNode->GetXPPCurr();
609  const Vec3& wpRef = pRefNode->GetWPCurr();
610 
611  if (bLabels) {
612  uint32_t l = pRefNode->GetLabel();
613  send(outfd, (void *)&l, sizeof(l), 0);
614  }
615 
616  send(outfd, (void *)xRef.pGetVec(), 3*sizeof(doublereal), 0);
617  switch (uRRot) {
618  case MBC_ROT_MAT:
619  send(outfd, (void *)RRef.pGetMat(), 9*sizeof(doublereal), 0);
620  break;
621 
622  case MBC_ROT_THETA: {
623  Vec3 Theta(RotManip::VecRot(RRef));
624  send(outfd, (void *)Theta.pGetVec(), 3*sizeof(doublereal), 0);
625  } break;
626 
627  case MBC_ROT_EULER_123: {
629  send(outfd, (void *)E.pGetVec(), 3*sizeof(doublereal), 0);
630  } break;
631  }
632  send(outfd, (void *)xpRef.pGetVec(), 3*sizeof(doublereal), 0);
633  send(outfd, (void *)wRef.pGetVec(), 3*sizeof(doublereal), 0);
634  if (bOutputAccelerations) {
635  send(outfd, (void *)xppRef.pGetVec(), 3*sizeof(doublereal), 0);
636  send(outfd, (void *)wpRef.pGetVec(), 3*sizeof(doublereal), 0);
637  }
638 
639  for (unsigned p3 = 0, n = 0; n < Nodes.size(); n++) {
640  const StructNode *pNode(dynamic_cast<const StructNode *>(Nodes[n].pNode));
641  if (pNode != 0) {
642  for (unsigned o = 0; o < Nodes[n].Offsets.size(); o++, p3 += 3) {
643  Vec3 f(pNode->GetRCurr()*Nodes[n].Offsets[o].Offset);
644  Vec3 x(pNode->GetXCurr() + f);
645  Vec3 Dx(x - xRef);
646  Vec3 v(pNode->GetVCurr() + pNode->GetWCurr().Cross(f));
647  Vec3 Dv(v - xpRef - wRef.Cross(Dx));
648  const Vec3& w(pNode->GetWCurr());
649 
650  Vec3 xTilde(RRef.MulTV(Dx));
651  m_x.Put(p3 + 1, xTilde);
652 
653  Vec3 vTilde(RRef.MulTV(Dv));
654  m_xP.Put(p3 + 1, vTilde);
655 
656  if (bOutputAccelerations) {
657  const Vec3& xpp = pNode->GetXPPCurr();
658  const Vec3& wp = pNode->GetWPCurr();
659 
660  Vec3 xppTilde(RRef.MulTV(xpp - xppRef - wpRef.Cross(Dx)
661  - wRef.Cross(wRef.Cross(Dx) + Dv*2)
662  + wp.Cross(f) + w.Cross(w.Cross(f))));
663  m_xPP.Put(p3 + 1, xppTilde);
664  }
665  }
666 
667  } else {
668  Vec3 Dx(Nodes[n].pNode->GetXCurr() - xRef);
669  Vec3 Dv(Nodes[n].pNode->GetVCurr() - xpRef - wRef.Cross(Dx));
670 
671  Vec3 xTilde(RRef.MulTV(Dx));
672  m_x.Put(p3 + 1, xTilde);
673 
674  Vec3 vTilde(RRef.MulTV(Dv));
675  m_xP.Put(p3 + 1, vTilde);
676 
677  if (bOutputAccelerations) {
678  Vec3 xppTilde(RRef.MulTV(Nodes[n].pNode->GetXPPCurr()
679  - xppRef - wpRef.Cross(Dx)
680  - wRef.Cross(wRef.Cross(Dx) + Dv*2)));
681  m_xPP.Put(p3 + 1, xppTilde);
682  }
683 
684  p3 += 3;
685  }
686  }
687 
688  } else {
689  for (unsigned p3 = 0, n = 0; n < Nodes.size(); n++) {
690  const StructNode *pNode(dynamic_cast<const StructNode *>(Nodes[n].pNode));
691  if (pNode != 0) {
692  for (unsigned o = 0; o < Nodes[n].Offsets.size(); o++, p3 +=3 ) {
693  /*
694  p = x + f
695  R = R
696  v = xp + w cross f
697  w = w
698  a = xpp + wp cross f + w cross w cross f
699  wp = wp
700  */
701 
702  // Optimization of the above formulas
703  const Mat3x3& R = pNode->GetRCurr();
704  Vec3 f = R*Nodes[n].Offsets[o].Offset;
705  Vec3 x = pNode->GetXCurr() + f;
706  const Vec3& w = pNode->GetWCurr();
707  Vec3 wCrossf = w.Cross(f);
708  Vec3 v = pNode->GetVCurr() + wCrossf;
709 
710  m_x.Put(p3 + 1, x);
711  m_xP.Put(p3 + 1, v);
712 
713  if (bOutputAccelerations) {
714  const Vec3& wp = pNode->GetWPCurr();
715  Vec3 xpp = pNode->GetXPPCurr() + wp.Cross(f) + w.Cross(wCrossf);
716 
717  m_xPP.Put(p3 + 1, xpp);
718  }
719  }
720 
721  } else {
722  m_x.Put(p3 + 1, Nodes[n].pNode->GetXCurr());
723  m_xP.Put(p3 + 1, Nodes[n].pNode->GetVCurr());
724 
725  if (bOutputAccelerations) {
726  m_xPP.Put(p3 + 1, Nodes[n].pNode->GetXPPCurr());
727  }
728 
729  p3 += 3;
730  }
731  }
732  }
733 
734  if (bLabels) {
735  send(outfd, &m_qlabels[0], sizeof(uint32_t)*m_qlabels.size(), 0);
736  }
737 
738  if (pH) {
739  pH->MatVecMul(m_q, m_x);
740  pH->MatVecMul(m_qP, m_xP);
741 
742  send(outfd, &m_q[0], sizeof(double)*m_q.size(), 0);
743  send(outfd, &m_qP[0], sizeof(double)*m_qP.size(), 0);
744 
745  if (bOutputAccelerations) {
746  pH->MatVecMul(m_qPP, m_xPP);
747  send(outfd, &m_qPP[0], sizeof(double)*m_qPP.size(), 0);
748  }
749 
750  } else {
751  send(outfd, &m_x[0], sizeof(double)*m_x.size(), 0);
752  send(outfd, &m_xP[0], sizeof(double)*m_xP.size(), 0);
753 
754  if (bOutputAccelerations) {
755  send(outfd, &m_xPP[0], sizeof(double)*m_xPP.size(), 0);
756  }
757  }
758 
759 #else // ! USE_SOCKET
761 #endif // ! USE_SOCKET
762 }
763 
764 void
766 {
767  std::istream *infp = pEFH->GetInStream();
768  if (infp) {
769  RecvFromStream(*infp);
770 
771  } else {
772  RecvFromFileDes(pEFH->GetInFileDes());
773  }
774 }
775 
776 void
778 {
779 #if 0
780  if (pRefNode) {
781  unsigned l;
782  doublereal *f = F0.pGetVec(), *m = M0.pGetVec();
783 
784  if (bLabels) {
785  inf >> l;
786  }
787 
788  inf >> f[0] >> f[1] >> f[2];
789  if (uRRot != MBC_ROT_NONE) {
790  inf >> m[0] >> m[1] >> m[2];
791  }
792  }
793 
794  for (unsigned i = 0; i < Nodes.size(); i++) {
795  /* assume unsigned int label */
796  unsigned l;
797  doublereal f[3], m[3];
798 
799  if (bLabels) {
800  inf >> l;
801 
802  if (Nodes[i]->GetLabel() != l) {
803  silent_cerr("StructMappingExtForce"
804  "(" << GetLabel() << "): "
805  "invalid " << i << "-th label " << l
806  << std::endl);
808  }
809  }
810 
811  inf
812  >> f[0] >> f[1] >> f[2];
813  if (uRot != MBC_ROT_NONE) {
814  inf >> m[0] >> m[1] >> m[2];
815  }
816 
817  if (!inf) {
818  break;
819  }
820 
821  F[i] = Vec3(f);
822  if (uRot != MBC_ROT_NONE) {
823  M[i] = Vec3(m);
824  }
825  }
826 #endif
827 }
828 
829 void
831 {
832 #ifdef USE_SOCKET
833  if (pRefNode) {
834  size_t ulen = 0;
835  char buf[sizeof(uint32_t) + 6*sizeof(doublereal)];
836  doublereal *f;
837  ssize_t len;
838 
839  if (bLabels) {
840  ulen = sizeof(uint32_t);
841  }
842 
843  ulen += 6*sizeof(doublereal);
844 
845  len = recv(infd, (void *)buf, ulen, pEFH->GetRecvFlags());
846  if (len == -1) {
847  int save_errno = errno;
848  char *err_msg = strerror(save_errno);
849  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
850  "recv() failed (" << save_errno << ": "
851  << err_msg << ")" << std::endl);
853 
854  } else if (unsigned(len) != ulen) {
855  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
856  "recv() failed " "(got " << len << " of "
857  << ulen << " bytes)" << std::endl);
859  }
860 
861  if (bLabels) {
862  uint32_t *uint32_ptr = (uint32_t *)buf;
863  unsigned l = uint32_ptr[0];
864  if (l != pRefNode->GetLabel()) {
865  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
866  "invalid reference node label "
867  "(wanted " << pRefNode->GetLabel() << ", got " << l << ")"
868  << std::endl);
870  }
871  f = (doublereal *)&uint32_ptr[1];
872 
873  } else {
874  f = (doublereal *)buf;
875  }
876 
877  F0 = Vec3(&f[0]);
878  M0 = Vec3(&f[3]);
879  }
880 
881  if (bLabels) {
882  // Hack!
883  ssize_t len = recv(infd, (void *)&m_p[0], sizeof(uint32_t)*m_p.size(),
884  pEFH->GetRecvFlags());
885  if (len == -1) {
886  int save_errno = errno;
887  char *err_msg = strerror(save_errno);
888  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
889  "recv() failed (" << save_errno << ": "
890  << err_msg << ")" << std::endl);
892 
893  } else if (unsigned(len) != sizeof(double)*m_p.size()) {
894  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
895  "recv() failed " "(got " << len << " of "
896  << sizeof(uint32_t)*m_p.size() << " bytes)" << std::endl);
898  }
899 
900  uint32_t *labels = (uint32_t *)&m_p[0];
901  for (unsigned l = 0; l < m_qlabels.size(); l++) {
902  if (labels[l] != m_qlabels[l]) {
903  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
904  "label mismatch for point #" << l << "/" << m_qlabels.size()
905  << " local=" << m_qlabels[l] << " remote=" << labels[l] << std::endl);
907  }
908  }
909  }
910 
911  size_t fsize;
912  double *fp;
913 
914  if (pH) {
915  fp = &m_p[0];
916  fsize = sizeof(double)*m_p.size();
917 
918  } else {
919  fp = &m_f[0];
920  fsize = sizeof(double)*m_f.size();
921  }
922 
923  ssize_t len = recv(infd, (void *)fp, fsize, pEFH->GetRecvFlags());
924  if (len == -1) {
925  int save_errno = errno;
926  char *err_msg = strerror(save_errno);
927  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
928  "recv() failed (" << save_errno << ": "
929  << err_msg << ")" << std::endl);
931 
932  } else if (unsigned(len) != fsize) {
933  silent_cerr("StructMappingExtForce(" << GetLabel() << "): "
934  "recv() failed " "(got " << len << " of "
935  << fsize << " bytes)" << std::endl);
937  }
938 
939  if (pH) {
940  pH->MatTVecMul(m_f, m_p);
941  }
942 
943  if (pRefNode) {
944  for (unsigned p3 = 0, n = 0; n < Nodes.size(); n++) {
945  Nodes[n].F = Zero3;
946  Nodes[n].M = Zero3;
947  const StructNode *pNode(dynamic_cast<const StructNode *>(Nodes[n].pNode));
948  if (pNode != 0) {
949  for (unsigned o = 0; o < Nodes[n].Offsets.size(); o++, p3 += 3) {
950  Nodes[n].Offsets[o].F = pRefNode->GetRCurr()*Vec3(&m_f[p3]);
951  Nodes[n].F += Nodes[n].Offsets[o].F;
952  Vec3 f(pNode->GetRCurr()*Nodes[n].Offsets[o].Offset);
953  Nodes[n].M += f.Cross(Nodes[n].Offsets[o].F);
954  }
955 
956  } else {
957  Nodes[n].Offsets[0].F = pRefNode->GetRCurr()*Vec3(&m_f[p3]);
958  Nodes[n].F += Nodes[n].Offsets[0].F;
959  }
960  }
961 
962  } else {
963  for (unsigned p3 = 0, n = 0; n < Nodes.size(); n++) {
964  Nodes[n].F = Zero3;
965  Nodes[n].M = Zero3;
966  const StructNode *pNode(dynamic_cast<const StructNode *>(Nodes[n].pNode));
967  if (pNode != 0) {
968  for (unsigned o = 0; o < Nodes[n].Offsets.size(); o++, p3 += 3) {
969  Nodes[n].Offsets[o].F = Vec3(&m_f[p3]);
970  Nodes[n].F += Nodes[n].Offsets[o].F;
971  Vec3 f(pNode->GetRCurr()*Nodes[n].Offsets[o].Offset);
972  Nodes[n].M += f.Cross(Nodes[n].Offsets[o].F);
973  }
974 
975  } else {
976  Nodes[n].Offsets[0].F = Vec3(&m_f[p3]);
977  Nodes[n].F += Nodes[n].Offsets[0].F;
978  }
979  }
980  }
981 #else // ! USE_SOCKET
983 #endif // ! USE_SOCKET
984 }
985 
988  doublereal dCoef,
989  const VectorHandler& XCurr,
990  const VectorHandler& XPrimeCurr)
991 {
992  ExtForce::Recv();
993 
994  if (iCoupling == COUPLING_NONE) {
995  WorkVec.Resize(0);
996  return WorkVec;
997  }
998 
999 
1000  if (pRefNode) {
1001  integer iSize(0);
1002 
1003  WorkVec.ResizeReset(m_uResSize);
1004 
1005  const Vec3& xRef = pRefNode->GetXCurr();
1006  const Mat3x3& RRef = pRefNode->GetRCurr();
1007 
1008  // manipulate
1011  F1 = RRef*F0;
1012  M1 = RRef*M0;
1013 
1014  } else {
1015  F1 = F0;
1016  M1 = M0;
1017  }
1018  }
1019 
1020  F2 = Zero3;
1021  M2 = Zero3;
1022  for (unsigned n = 0; n < Nodes.size(); n++) {
1023  integer iFirstIndex = Nodes[n].pNode->iGetFirstMomentumIndex();
1024  integer iDim;
1025 
1026  WorkVec.Add(iSize + 1, Nodes[n].F);
1027  if (dynamic_cast<const StructNode *>(Nodes[n].pNode)) {
1028  WorkVec.Add(iSize + 4, Nodes[n].M);
1029  iDim = 6;
1030 
1031  } else {
1032  iDim = 3;
1033  }
1034 
1035  for (int r = 1; r <= iDim; r++) {
1036  WorkVec.PutRowIndex(iSize + r, iFirstIndex + r);
1037  }
1038 
1040  // compute Global Reference Node Forces, even if they are not used, for output only :)
1041  F2 += Nodes[n].F;
1042  M2 += Nodes[n].M + (Nodes[n].pNode->GetXCurr() - xRef).Cross(Nodes[n].F);
1043  }
1044 
1045  iSize += iDim;
1046  }
1047 
1049  integer iFirstIndex = pRefNode->iGetFirstMomentumIndex();
1050  for (int r = 1; r <= 6; r++) {
1051  WorkVec.PutRowIndex(iSize + r, iFirstIndex + r);
1052  }
1053 
1054  F1 -= F2;
1055  M1 -= M2;
1056  WorkVec.Add(iSize + 1, F1);
1057  WorkVec.Add(iSize + 4, M1);
1058 
1059  iSize += 6;
1060  }
1061 
1062  ASSERT(iSize == m_uResSize);
1063 
1064  } else {
1065  integer iSize(0);
1066 
1067  WorkVec.ResizeReset(m_uResSize);
1068 
1069  for (unsigned n = 0; n < Nodes.size(); n++) {
1070  integer iFirstIndex = Nodes[n].pNode->iGetFirstMomentumIndex();
1071  integer iDim;
1072 
1073  WorkVec.Add(iSize + 1, Nodes[n].F);
1074  if (dynamic_cast<const StructNode *>(Nodes[n].pNode)) {
1075  WorkVec.Add(iSize + 4, Nodes[n].M);
1076  iDim = 6;
1077 
1078  } else {
1079  iDim = 3;
1080  }
1081 
1082  for (int r = 1; r <= iDim; r++) {
1083  WorkVec.PutRowIndex(iSize + r, iFirstIndex + r);
1084  }
1085 
1086  iSize += iDim;
1087  }
1088 
1089  ASSERT(iSize == m_uResSize);
1090  }
1091 
1092  return WorkVec;
1093 }
1094 
1095 void
1097 {
1098  if (bToBeOutput()) {
1099  if (OH.UseText(OutputHandler::FORCES)) {
1100  std::ostream& out = OH.Forces();
1101 
1102  if (pRefNode) {
1103  out << GetLabel() << "#" << pRefNode->GetLabel()
1104  << " " << F0
1105  << " " << M0
1106  << " " << F1
1107  << " " << M1
1108  << " " << F2
1109  << " " << M2
1110  << std::endl;
1111  }
1112 
1113  for (unsigned n = 0; n < Nodes.size(); n++) {
1114  out << GetLabel() << "@" << Nodes[n].pNode->GetLabel()
1115  << " " << Nodes[n].F
1116  << " " << Nodes[n].M
1117  << std::endl;
1118  }
1119  }
1120 
1121  /* TODO: NetCDF */
1122  }
1123 }
1124 
1125 void
1126 StructMappingExtForce::GetConnectedNodes(std::vector<const Node *>& connectedNodes) const
1127 {
1128  unsigned n = Nodes.size();
1129  if (pRefNode) {
1130  n++;
1131  }
1132  connectedNodes.resize(n);
1133 
1134  for (n = 0; n < Nodes.size(); n++) {
1135  connectedNodes[n] = Nodes[n].pNode;
1136  }
1137 
1138  if (pRefNode) {
1139  connectedNodes[n] = pRefNode;
1140  }
1141 }
1142 
1143 /* StructMappingExtForce - end */
1144 
1145 
1146 /* StructMembraneMappingExtForce - begin */
1147 
1148 /* Costruttore */
1150  DataManager *pDM,
1151  const StructNode *pRefNode,
1152  bool bUseReferenceNodeForces,
1153  bool bRotateReferenceNodeForces,
1154  std::vector<const StructDispNode *>& nodes,
1155  std::vector<Vec3>& offsets,
1156  std::vector<unsigned>& labels,
1157  std::vector<NodeConnData>& nodesConn,
1158  SpMapMatrixHandler *pH,
1159  std::vector<uint32_t>& mappedlabels,
1160  bool bLabels,
1161  bool bOutputAccelerations,
1162  unsigned uRRot,
1163  ExtFileHandlerBase *pEFH,
1164  bool bSendAfterPredict,
1165  int iCoupling,
1166  flag fOut)
1167 : Elem(uL, fOut),
1168 StructMappingExtForce(uL, pDM,
1169  pRefNode, bUseReferenceNodeForces, bRotateReferenceNodeForces,
1170  nodes, offsets, labels, pH, mappedlabels,
1171  bLabels, bOutputAccelerations, uRRot,
1172  pEFH, bSendAfterPredict, iCoupling,
1173  fOut)
1174 {
1175  NodesConn.resize(Nodes.size());
1176  for (unsigned n = 0; n < NodesConn.size(); n++) {
1177  NodesConn[n] = nodesConn[n];
1178  }
1179 }
1180 
1182 {
1183  NO_OP;
1184 }
1185 
1186 /*
1187  * Send output to companion software
1188  */
1189 void
1191 {
1192  if (pRefNode) {
1193  const Vec3& xRef = pRefNode->GetXCurr();
1194  const Mat3x3& RRef = pRefNode->GetRCurr();
1195  const Vec3& xpRef = pRefNode->GetVCurr();
1196  const Vec3& wRef = pRefNode->GetWCurr();
1197  const Vec3& xppRef = pRefNode->GetXPPCurr();
1198  const Vec3& wpRef = pRefNode->GetWPCurr();
1199 
1200  if (bLabels) {
1201  outf
1202  << pRefNode->GetLabel()
1203  << " ";
1204  }
1205 
1206  outf
1207  << xRef
1208  << " " << RRef
1209  << " " << xpRef
1210  << " " << wRef;
1211 
1212  if (bOutputAccelerations) {
1213  outf
1214  << " " << xppRef
1215  << " " << wpRef;
1216  }
1217  outf << std::endl;
1218 
1219  for (unsigned i = 0; i < Nodes.size(); i++) {
1220 #if 0
1221  Vec3 f(Nodes[i]->GetRCurr()*Offsets[i]);
1222  Vec3 x(Nodes[i]->GetXCurr() + f);
1223  Vec3 Dx(x - xRef);
1224  Mat3x3 DR(RRef.MulTM(Nodes[i]->GetRCurr()));
1225  Vec3 v(Nodes[i]->GetVCurr() + Nodes[i]->GetWCurr().Cross(f));
1226  Vec3 Dv(v - xpRef - wRef.Cross(Dx));
1227  const Vec3& w(Nodes[i]->GetWCurr());
1228 
1229  // manipulate
1230 
1231  if (bLabels) {
1232  outf
1233  << Nodes[i]->GetLabel()
1234  << " ";
1235  }
1236 
1237  outf
1238  << RRef.MulTV(Dx);
1239 
1240  switch (uRot) {
1241  case MBC_ROT_NONE:
1242  break;
1243 
1244  case MBC_ROT_MAT:
1245  outf
1246  << " " << DR;
1247  break;
1248 
1249  case MBC_ROT_THETA:
1250  outf
1251  << " " << RotManip::VecRot(DR);
1252  break;
1253 
1254  case MBC_ROT_EULER_123:
1255  outf
1256  << " " << MatR2EulerAngles123(DR)*dRaDegr;
1257  break;
1258  }
1259 
1260  outf
1261  << " " << RRef.MulTV(Dv);
1262 
1263  if (uRot != MBC_ROT_NONE) {
1264  outf
1265  << " " << RRef.MulTV(w - wRef);
1266  }
1267 
1268  if (bOutputAccelerations) {
1269  const Vec3& xpp(Nodes[i]->GetXPPCurr());
1270 
1271  outf
1272  << " " << RRef.MulTV(xpp - xppRef - wpRef.Cross(Dx)
1273  - wRef.Cross(wRef.Cross(Dx) + Dv*2));
1274  if (uRot != MBC_ROT_NONE) {
1275  const Vec3& wp(Nodes[i]->GetWPCurr());
1276 
1277  outf
1278  << " " << RRef.MulTV(wp - wpRef - wRef.Cross(w));
1279  }
1280  }
1281  outf << std::endl;
1282 #endif
1283  }
1284 
1285  } else {
1286  for (unsigned i = 0; i < Nodes.size(); i++) {
1287 #if 0
1288  /*
1289  p = x + f
1290  R = R
1291  v = xp + w cross f
1292  w = w
1293  a = xpp + wp cross f + w cross w cross f
1294  wp = wp
1295  */
1296 
1297  // Optimization of the above formulas
1298  const Mat3x3& R = Nodes[i]->GetRCurr();
1299  Vec3 f = R*Offsets[i];
1300  Vec3 x = Nodes[i]->GetXCurr() + f;
1301  const Vec3& w = Nodes[i]->GetWCurr();
1302  Vec3 wCrossf = w.Cross(f);
1303  Vec3 v = Nodes[i]->GetVCurr() + wCrossf;
1304 
1305  if (bLabels) {
1306  outf
1307  << Nodes[i]->GetLabel()
1308  << " ";
1309  }
1310 
1311  outf
1312  << x;
1313 
1314  switch (uRot) {
1315  case MBC_ROT_NONE:
1316  break;
1317 
1318  case MBC_ROT_MAT:
1319  outf
1320  << " " << R;
1321  break;
1322 
1323  case MBC_ROT_THETA:
1324  outf
1325  << " " << RotManip::VecRot(R);
1326  break;
1327 
1328  case MBC_ROT_EULER_123:
1329  outf
1330  << " " << MatR2EulerAngles123(R)*dRaDegr;
1331  break;
1332  }
1333  outf
1334  << " " << v;
1335 
1336  if (uRot != MBC_ROT_NONE) {
1337  outf
1338  << " " << w;
1339  }
1340 
1341  if (bOutputAccelerations) {
1342  const Vec3& wp = Nodes[i]->GetWPCurr();
1343  Vec3 a = Nodes[i]->GetXPPCurr() + wp.Cross(f) + w.Cross(wCrossf);
1344 
1345  outf
1346  << " " << a;
1347 
1348  if (uRot != MBC_ROT_NONE) {
1349  outf
1350  << " " << wp;
1351  }
1352  }
1353 
1354  outf << std::endl;
1355 #endif
1356  }
1357  }
1358 }
1359 
1360 void
1362 {
1363 #ifdef USE_SOCKET
1364  if (pRefNode) {
1365  const Vec3& xRef = pRefNode->GetXCurr();
1366  const Mat3x3& RRef = pRefNode->GetRCurr();
1367  const Vec3& xpRef = pRefNode->GetVCurr();
1368  const Vec3& wRef = pRefNode->GetWCurr();
1369  const Vec3& xppRef = pRefNode->GetXPPCurr();
1370  const Vec3& wpRef = pRefNode->GetWPCurr();
1371 
1372  if (bLabels) {
1373  uint32_t l = pRefNode->GetLabel();
1374  send(outfd, (void *)&l, sizeof(l), 0);
1375  }
1376 
1377  send(outfd, (void *)xRef.pGetVec(), 3*sizeof(doublereal), 0);
1378  switch (uRRot) {
1379  case MBC_ROT_MAT:
1380  send(outfd, (void *)RRef.pGetMat(), 9*sizeof(doublereal), 0);
1381  break;
1382 
1383  case MBC_ROT_THETA: {
1384  Vec3 Theta(RotManip::VecRot(RRef));
1385  send(outfd, (void *)Theta.pGetVec(), 3*sizeof(doublereal), 0);
1386  } break;
1387 
1388  case MBC_ROT_EULER_123: {
1389  Vec3 E(MatR2EulerAngles123(RRef)*dRaDegr);
1390  send(outfd, (void *)E.pGetVec(), 3*sizeof(doublereal), 0);
1391  } break;
1392  }
1393  send(outfd, (void *)xpRef.pGetVec(), 3*sizeof(doublereal), 0);
1394  send(outfd, (void *)wRef.pGetVec(), 3*sizeof(doublereal), 0);
1395  if (bOutputAccelerations) {
1396  send(outfd, (void *)xppRef.pGetVec(), 3*sizeof(doublereal), 0);
1397  send(outfd, (void *)wpRef.pGetVec(), 3*sizeof(doublereal), 0);
1398  }
1399 
1400  for (unsigned p3 = 0, n = 0; n < Nodes.size(); n++) {
1401  const StructNode *pNode(dynamic_cast<const StructNode *>(Nodes[n].pNode));
1402  if (pNode != 0) {
1403  for (unsigned o = 0; o < Nodes[n].Offsets.size(); o++, p3 += 3) {
1404  Vec3 f(pNode->GetRCurr()*Nodes[n].Offsets[o].Offset);
1405  Vec3 x(pNode->GetXCurr() + f);
1406  Vec3 Dx(x - xRef);
1407  Vec3 v(pNode->GetVCurr() + pNode->GetWCurr().Cross(f));
1408  Vec3 Dv(v - xpRef - wRef.Cross(Dx));
1409  const Vec3& w(pNode->GetWCurr());
1410 
1411  Vec3 xTilde(RRef.MulTV(Dx));
1412  m_x.Put(p3 + 1, xTilde);
1413 
1414  Vec3 vTilde(RRef.MulTV(Dv));
1415  m_xP.Put(p3 + 1, vTilde);
1416 
1417  if (bOutputAccelerations) {
1418  const Vec3& xpp = pNode->GetXPPCurr();
1419  const Vec3& wp = pNode->GetWPCurr();
1420 
1421  Vec3 xppTilde(RRef.MulTV(xpp - xppRef - wpRef.Cross(Dx)
1422  - wRef.Cross(wRef.Cross(Dx) + Dv*2)
1423  + wp.Cross(f) + w.Cross(w.Cross(f))));
1424  m_xPP.Put(p3 + 1, xppTilde);
1425  }
1426  }
1427 
1428  } else {
1429  if (NodesConn[n].pNode[0] != 0) {
1430  Vec3 e1(NodesConn[n].pNode[1]->GetXCurr() - NodesConn[n].pNode[0]->GetXCurr());
1431  Vec3 e2(NodesConn[n].pNode[3]->GetXCurr() - NodesConn[n].pNode[2]->GetXCurr());
1432  Vec3 e3(e1.Cross(e2));
1433  e3 /= e3.Norm();
1434 
1435  Nodes[n].Offsets[0].Offset = e3*NodesConn[n].h1;
1436  Nodes[n].Offsets[1].Offset = e3*NodesConn[n].h2;
1437 
1438  for (unsigned o = 0; o < 2; o++, p3 += 3) {
1439  Vec3 Dx(Nodes[n].pNode->GetXCurr() + Nodes[n].Offsets[o].Offset - xRef);
1440  Vec3 Dv(Nodes[n].pNode->GetVCurr() - xpRef - wRef.Cross(Dx));
1441 
1442  Vec3 xTilde(RRef.MulTV(Dx));
1443  m_x.Put(p3 + 1, xTilde);
1444 
1445  Vec3 vTilde(RRef.MulTV(Dv));
1446  m_xP.Put(p3 + 1, vTilde);
1447 
1448  if (bOutputAccelerations) {
1449  Vec3 xppTilde(RRef.MulTV(Nodes[n].pNode->GetXPPCurr()
1450  - xppRef - wpRef.Cross(Dx)
1451  - wRef.Cross(wRef.Cross(Dx) + Dv*2)));
1452  m_xPP.Put(p3 + 1, xppTilde);
1453  }
1454  }
1455 
1456  } else {
1457  Vec3 Dx(Nodes[n].pNode->GetXCurr() - xRef);
1458  Vec3 Dv(Nodes[n].pNode->GetVCurr() - xpRef - wRef.Cross(Dx));
1459 
1460  Vec3 xTilde(RRef.MulTV(Dx));
1461  m_x.Put(p3 + 1, xTilde);
1462 
1463  Vec3 vTilde(RRef.MulTV(Dv));
1464  m_xP.Put(p3 + 1, vTilde);
1465 
1466  if (bOutputAccelerations) {
1467  Vec3 xppTilde(RRef.MulTV(Nodes[n].pNode->GetXPPCurr()
1468  - xppRef - wpRef.Cross(Dx)
1469  - wRef.Cross(wRef.Cross(Dx) + Dv*2)));
1470  m_xPP.Put(p3 + 1, xppTilde);
1471  }
1472 
1473  p3 += 3;
1474  }
1475  }
1476  }
1477 
1478  } else {
1479  for (unsigned p3 = 0, n = 0; n < Nodes.size(); n++) {
1480  const StructNode *pNode(dynamic_cast<const StructNode *>(Nodes[n].pNode));
1481  if (pNode != 0) {
1482  for (unsigned o = 0; o < Nodes[n].Offsets.size(); o++, p3 +=3 ) {
1483  /*
1484  p = x + f
1485  R = R
1486  v = xp + w cross f
1487  w = w
1488  a = xpp + wp cross f + w cross w cross f
1489  wp = wp
1490  */
1491 
1492  // Optimization of the above formulas
1493  const Mat3x3& R = pNode->GetRCurr();
1494  Vec3 f = R*Nodes[n].Offsets[o].Offset;
1495  Vec3 x = pNode->GetXCurr() + f;
1496  const Vec3& w = pNode->GetWCurr();
1497  Vec3 wCrossf = w.Cross(f);
1498  Vec3 v = pNode->GetVCurr() + wCrossf;
1499 
1500  m_x.Put(p3 + 1, x);
1501  m_xP.Put(p3 + 1, v);
1502 
1503  if (bOutputAccelerations) {
1504  const Vec3& wp = pNode->GetWPCurr();
1505  Vec3 xpp = pNode->GetXPPCurr() + wp.Cross(f) + w.Cross(wCrossf);
1506 
1507  m_xPP.Put(p3 + 1, xpp);
1508  }
1509  }
1510 
1511  } else {
1512  if (NodesConn[n].pNode[0] != 0) {
1513  Vec3 e1(NodesConn[n].pNode[1]->GetXCurr() - NodesConn[n].pNode[0]->GetXCurr());
1514  Vec3 e2(NodesConn[n].pNode[3]->GetXCurr() - NodesConn[n].pNode[2]->GetXCurr());
1515  Vec3 e3(e1.Cross(e2));
1516  e3 /= e3.Norm();
1517 
1518  Nodes[n].Offsets[0].Offset = e3*NodesConn[n].h1;
1519  Nodes[n].Offsets[1].Offset = e3*NodesConn[n].h2;
1520 
1521  for (unsigned o = 0; o < 2; o++, p3 += 3) {
1522  m_x.Put(p3 + 1, Nodes[n].pNode->GetXCurr() + Nodes[n].Offsets[o].Offset);
1523  m_xP.Put(p3 + 1, Nodes[n].pNode->GetVCurr());
1524 
1525  if (bOutputAccelerations) {
1526  m_xPP.Put(p3 + 1, Nodes[n].pNode->GetXPPCurr());
1527  }
1528  }
1529 
1530  } else {
1531  m_x.Put(p3 + 1, Nodes[n].pNode->GetXCurr());
1532  m_xP.Put(p3 + 1, Nodes[n].pNode->GetVCurr());
1533 
1534  if (bOutputAccelerations) {
1535  m_xPP.Put(p3 + 1, Nodes[n].pNode->GetXPPCurr());
1536  }
1537 
1538  p3 += 3;
1539  }
1540  }
1541  }
1542  }
1543 
1544  if (bLabels) {
1545  send(outfd, &m_qlabels[0], sizeof(uint32_t)*m_qlabels.size(), 0);
1546  }
1547 
1548  if (pH) {
1549  pH->MatVecMul(m_q, m_x);
1550  pH->MatVecMul(m_qP, m_xP);
1551 
1552  send(outfd, &m_q[0], sizeof(double)*m_q.size(), 0);
1553  send(outfd, &m_qP[0], sizeof(double)*m_qP.size(), 0);
1554 
1555  if (bOutputAccelerations) {
1556  pH->MatVecMul(m_qPP, m_xPP);
1557  send(outfd, &m_qPP[0], sizeof(double)*m_qPP.size(), 0);
1558  }
1559 
1560  } else {
1561  send(outfd, &m_x[0], sizeof(double)*m_x.size(), 0);
1562  send(outfd, &m_xP[0], sizeof(double)*m_xP.size(), 0);
1563 
1564  if (bOutputAccelerations) {
1565  send(outfd, &m_xPP[0], sizeof(double)*m_xPP.size(), 0);
1566  }
1567  }
1568 
1569 #else // ! USE_SOCKET
1571 #endif // ! USE_SOCKET
1572 }
1573 
1574 void
1576 {
1577 #if 0
1578  if (pRefNode) {
1579  unsigned l;
1580  doublereal *f = F0.pGetVec(), *m = M0.pGetVec();
1581 
1582  if (bLabels) {
1583  inf >> l;
1584  }
1585 
1586  inf >> f[0] >> f[1] >> f[2];
1587  if (uRRot != MBC_ROT_NONE) {
1588  inf >> m[0] >> m[1] >> m[2];
1589  }
1590  }
1591 
1592  for (unsigned i = 0; i < Nodes.size(); i++) {
1593  /* assume unsigned int label */
1594  unsigned l;
1595  doublereal f[3], m[3];
1596 
1597  if (bLabels) {
1598  inf >> l;
1599 
1600  if (Nodes[i]->GetLabel() != l) {
1601  silent_cerr("StructMembraneMappingExtForce"
1602  "(" << GetLabel() << "): "
1603  "invalid " << i << "-th label " << l
1604  << std::endl);
1606  }
1607  }
1608 
1609  inf
1610  >> f[0] >> f[1] >> f[2];
1611  if (uRot != MBC_ROT_NONE) {
1612  inf >> m[0] >> m[1] >> m[2];
1613  }
1614 
1615  if (!inf) {
1616  break;
1617  }
1618 
1619  F[i] = Vec3(f);
1620  if (uRot != MBC_ROT_NONE) {
1621  M[i] = Vec3(m);
1622  }
1623  }
1624 #endif
1625 }
1626 
1627 void
1629 {
1630 #ifdef USE_SOCKET
1631  if (pRefNode) {
1632  size_t ulen = 0;
1633  char buf[sizeof(uint32_t) + 6*sizeof(doublereal)];
1634  doublereal *f;
1635  ssize_t len;
1636 
1637  if (bLabels) {
1638  ulen = sizeof(uint32_t);
1639  }
1640 
1641  ulen += 6*sizeof(doublereal);
1642 
1643  len = recv(infd, (void *)buf, ulen, pEFH->GetRecvFlags());
1644  if (len == -1) {
1645  int save_errno = errno;
1646  char *err_msg = strerror(save_errno);
1647  silent_cerr("StructMembraneMappingExtForce(" << GetLabel() << "): "
1648  "recv() failed (" << save_errno << ": "
1649  << err_msg << ")" << std::endl);
1651 
1652  } else if (unsigned(len) != ulen) {
1653  silent_cerr("StructMembraneMappingExtForce(" << GetLabel() << "): "
1654  "recv() failed " "(got " << len << " of "
1655  << ulen << " bytes)" << std::endl);
1657  }
1658 
1659  if (bLabels) {
1660  uint32_t *uint32_ptr = (uint32_t *)buf;
1661  unsigned l = uint32_ptr[0];
1662  if (l != pRefNode->GetLabel()) {
1663  silent_cerr("StructMembraneMappingExtForce(" << GetLabel() << "): "
1664  "invalid reference node label "
1665  "(wanted " << pRefNode->GetLabel() << ", got " << l << ")"
1666  << std::endl);
1668  }
1669  f = (doublereal *)&uint32_ptr[1];
1670 
1671  } else {
1672  f = (doublereal *)buf;
1673  }
1674 
1675  F0 = Vec3(&f[0]);
1676  M0 = Vec3(&f[3]);
1677  }
1678 
1679  if (bLabels) {
1680  // Hack!
1681  ssize_t len = recv(infd, (void *)&m_p[0], sizeof(uint32_t)*m_p.size(),
1682  pEFH->GetRecvFlags());
1683  if (len == -1) {
1684  int save_errno = errno;
1685  char *err_msg = strerror(save_errno);
1686  silent_cerr("StructMembraneMappingExtForce(" << GetLabel() << "): "
1687  "recv() failed (" << save_errno << ": "
1688  << err_msg << ")" << std::endl);
1690 
1691  } else if (unsigned(len) != sizeof(double)*m_p.size()) {
1692  silent_cerr("StructMembraneMappingExtForce(" << GetLabel() << "): "
1693  "recv() failed " "(got " << len << " of "
1694  << sizeof(uint32_t)*m_p.size() << " bytes)" << std::endl);
1696  }
1697 
1698  uint32_t *labels = (uint32_t *)&m_p[0];
1699  for (unsigned l = 0; l < m_qlabels.size(); l++) {
1700  if (labels[l] != m_qlabels[l]) {
1701  silent_cerr("StructMembraneMappingExtForce(" << GetLabel() << "): "
1702  "label mismatch for point #" << l << "/" << m_qlabels.size()
1703  << " local=" << m_qlabels[l] << " remote=" << labels[l] << std::endl);
1705  }
1706  }
1707  }
1708 
1709  size_t fsize;
1710  double *fp;
1711 
1712  if (pH) {
1713  fp = &m_p[0];
1714  fsize = sizeof(double)*m_p.size();
1715 
1716  } else {
1717  fp = &m_f[0];
1718  fsize = sizeof(double)*m_f.size();
1719  }
1720 
1721  ssize_t len = recv(infd, (void *)fp, fsize, pEFH->GetRecvFlags());
1722  if (len == -1) {
1723  int save_errno = errno;
1724  char *err_msg = strerror(save_errno);
1725  silent_cerr("StructMembraneMappingExtForce(" << GetLabel() << "): "
1726  "recv() failed (" << save_errno << ": "
1727  << err_msg << ")" << std::endl);
1729 
1730  } else if (unsigned(len) != fsize) {
1731  silent_cerr("StructMembraneMappingExtForce(" << GetLabel() << "): "
1732  "recv() failed " "(got " << len << " of "
1733  << fsize << " bytes)" << std::endl);
1735  }
1736 
1737  if (pH) {
1738  pH->MatTVecMul(m_f, m_p);
1739  }
1740 
1741  if (pRefNode) {
1742  for (unsigned p3 = 0, n = 0; n < Nodes.size(); n++) {
1743  Nodes[n].F = Zero3;
1744  Nodes[n].M = Zero3;
1745  const StructNode *pNode(dynamic_cast<const StructNode *>(Nodes[n].pNode));
1746  if (pNode != 0) {
1747  for (unsigned o = 0; o < Nodes[n].Offsets.size(); o++, p3 += 3) {
1748  Nodes[n].Offsets[o].F = pRefNode->GetRCurr()*Vec3(&m_f[p3]);
1749  Nodes[n].F += Nodes[n].Offsets[o].F;
1750  Vec3 f(pNode->GetRCurr()*Nodes[n].Offsets[o].Offset);
1751  Nodes[n].M += f.Cross(Nodes[n].Offsets[o].F);
1752  }
1753 
1754  } else {
1755  if (NodesConn[n].pNode[0] != 0) {
1756  for (unsigned o = 0; o < 2; o++, p3 += 3) {
1757  Nodes[n].Offsets[o].F = pRefNode->GetRCurr()*Vec3(&m_f[p3]);
1758  Nodes[n].F += Nodes[n].Offsets[o].F;
1759  }
1760 
1761  } else {
1762  Nodes[n].Offsets[0].F = pRefNode->GetRCurr()*Vec3(&m_f[p3]);
1763  Nodes[n].F += Nodes[n].Offsets[0].F;
1764  }
1765  }
1766  }
1767 
1768  } else {
1769  for (unsigned p3 = 0, n = 0; n < Nodes.size(); n++) {
1770  Nodes[n].F = Zero3;
1771  Nodes[n].M = Zero3;
1772  const StructNode *pNode(dynamic_cast<const StructNode *>(Nodes[n].pNode));
1773  if (pNode != 0) {
1774  for (unsigned o = 0; o < Nodes[n].Offsets.size(); o++, p3 += 3) {
1775  Nodes[n].Offsets[o].F = Vec3(&m_f[p3]);
1776  Nodes[n].F += Nodes[n].Offsets[o].F;
1777  Vec3 f(pNode->GetRCurr()*Nodes[n].Offsets[o].Offset);
1778  Nodes[n].M += f.Cross(Nodes[n].Offsets[o].F);
1779  }
1780 
1781  } else {
1782  if (NodesConn[n].pNode[0] != 0) {
1783  for (unsigned o = 0; o < 2; o++, p3 += 3) {
1784  Nodes[n].Offsets[o].F = Vec3(&m_f[p3]);
1785  Nodes[n].F += Nodes[n].Offsets[o].F;
1786  }
1787 
1788  } else {
1789  Nodes[n].Offsets[0].F = Vec3(&m_f[p3]);
1790  Nodes[n].F += Nodes[n].Offsets[0].F;
1791  }
1792  }
1793  }
1794  }
1795 #else // ! USE_SOCKET
1797 #endif // ! USE_SOCKET
1798 }
1799 
1800 /* StructMembraneMappingExtForce - end */
1801 
1802 
1803 Elem*
1805  MBDynParser& HP,
1806  unsigned int uLabel)
1807 {
1808  ExtFileHandlerBase *pEFH = 0;
1809  int iCoupling;
1810 
1811  bool bSendAfterPredict;
1812  ReadExtForce(pDM, HP, uLabel, pEFH, bSendAfterPredict, iCoupling);
1813 
1814  const StructNode *pRefNode(0);
1815  if (HP.IsKeyWord("reference" "node")) {
1816  pRefNode = pDM->ReadNode<const StructNode, Node::STRUCTURAL>(HP);
1817  }
1818 
1819  bool bLabels(false);
1820  unsigned uRRot = pRefNode ? MBC_ROT_MAT : MBC_ROT_NONE;
1821  bool bOutputAccelerations(false);
1822  bool bUseReferenceNodeForces(true);
1823  bool bRotateReferenceNodeForces(true);
1824 
1825  bool bGotLabels(false);
1826  bool bGotRot(false);
1827  bool bGotAccels(false);
1828  bool bGotUseRefForces(false);
1829 
1830  while (HP.IsArg()) {
1831  if (HP.IsKeyWord("no" "labels")) {
1832  silent_cerr("StructMappingExtForce(" << uLabel << "): "
1833  "use of \"no labels\" deprecated in favor of \"labels, { yes | no }\" at line "
1834  << HP.GetLineData() << std::endl);
1835 
1836  if (bGotLabels) {
1837  silent_cerr("StructMappingExtForce(" << uLabel << "): "
1838  "\"no labels\" already specified at line "
1839  << HP.GetLineData() << std::endl);
1841  }
1842 
1843  bLabels = false;
1844  bGotLabels = true;
1845 
1846  } else if (HP.IsKeyWord("labels")) {
1847  if (bGotLabels) {
1848  silent_cerr("StructMappingExtForce(" << uLabel << "): "
1849  "\"labels\" already specified at line "
1850  << HP.GetLineData() << std::endl);
1852  }
1853 
1854  if (!HP.GetYesNo(bLabels)) {
1855  silent_cerr("StructMappingExtForce(" << uLabel << "): "
1856  "\"labels\" must be either \"yes\" or \"no\" at line "
1857  << HP.GetLineData() << std::endl);
1859  }
1860  bGotLabels = true;
1861 
1862  } else if (HP.IsKeyWord("orientation")) {
1863  if (bGotRot) {
1864  silent_cerr("StructMappingExtForce(" << uLabel << "): "
1865  "\"orientation\" already specified at line "
1866  << HP.GetLineData() << std::endl);
1868  }
1869 
1870  if (HP.IsKeyWord("none")) {
1871  uRRot = MBC_ROT_NONE;
1872 
1873  } else if (HP.IsKeyWord("orientation" "vector")) {
1874  uRRot = MBC_ROT_THETA;
1875 
1876  } else if (HP.IsKeyWord("orientation" "matrix")) {
1877  uRRot = MBC_ROT_MAT;
1878 
1879  } else if (HP.IsKeyWord("euler" "123")) {
1880  uRRot = MBC_ROT_EULER_123;
1881 
1882  } else {
1883  silent_cerr("StructMappingExtForce(" << uLabel << "): "
1884  "unknown \"orientation\" format at line "
1885  << HP.GetLineData() << std::endl);
1887  }
1888 
1889  bGotRot = true;
1890 
1891  } else if (HP.IsKeyWord("accelerations")) {
1892  if (bGotAccels) {
1893  silent_cerr("StructMappingExtForce(" << uLabel << "): "
1894  "\"accelerations\" already specified at line "
1895  << HP.GetLineData() << std::endl);
1897  }
1898 
1899  if (!HP.GetYesNo(bOutputAccelerations)) {
1900  silent_cerr("StructMappingExtForce(" << uLabel << "): "
1901  "\"accelerations\" must be either \"yes\" or \"no\" at line "
1902  << HP.GetLineData() << std::endl);
1904  }
1905  bGotAccels = true;
1906 
1907  } else if (HP.IsKeyWord("use" "reference" "node" "forces")) {
1908  if (pRefNode == 0) {
1909  silent_cerr("StructMappingExtForce(" << uLabel << "): "
1910  "\"use reference node forces\" only meaningful when reference node is used at line "
1911  << HP.GetLineData() << std::endl);
1913  }
1914 
1915  if (bGotUseRefForces) {
1916  silent_cerr("StructMappingExtForce(" << uLabel << "): "
1917  "\"use reference node forces\" already specified at line "
1918  << HP.GetLineData() << std::endl);
1920  }
1921 
1922  if (!HP.GetYesNo(bUseReferenceNodeForces)) {
1923  silent_cerr("StructMappingExtForce(" << uLabel << "): "
1924  "\"use reference node forces\" must be either \"yes\" or \"no\" at line "
1925  << HP.GetLineData() << std::endl);
1927  }
1928  bGotUseRefForces = true;
1929 
1930  if (bUseReferenceNodeForces && HP.IsKeyWord("rotate" "reference" "node" "forces")) {
1931  if (!HP.GetYesNo(bRotateReferenceNodeForces)) {
1932  silent_cerr("StructMappingExtForce(" << uLabel << "): "
1933  "\"rotate reference node forces\" must be either \"yes\" or \"no\" at line "
1934  << HP.GetLineData() << std::endl);
1936  }
1937  }
1938 
1939  } else {
1940  break;
1941  }
1942  }
1943 
1944  if (!HP.IsKeyWord("points" "number")) {
1945  silent_cerr("StructMappingExtForce(" << uLabel << "): "
1946  "\"points number\" keyword expected "
1947  "at line " << HP.GetLineData() << std::endl);
1949  }
1950 
1951  int nPoints = HP.GetInt();
1952  if (nPoints <= 0) {
1953  silent_cerr("StructMappingExtForce(" << uLabel << "): illegal points number " << nPoints <<
1954  " at line " << HP.GetLineData() << std::endl);
1956  }
1957 
1958  std::vector<const StructDispNode *> Nodes(nPoints);
1959  std::vector<Vec3> Offsets(nPoints, ::Zero3);
1960  std::vector<uint32_t> Labels;
1961  if (bLabels) {
1962  Labels.resize(nPoints);
1963  }
1964  bool bMembrane(false);
1965  std::vector<StructMembraneMappingExtForce::NodeConnData> NodesConn;
1966 
1967  std::map<unsigned, bool> Got;
1968 
1969  for (unsigned n = 0, p = 0; p < unsigned(nPoints); n++) {
1970  const StructDispNode *pNode = pDM->ReadNode<const StructDispNode, Node::STRUCTURAL>(HP);
1971  unsigned uL(pNode->GetLabel());
1972  if (Got[uL]) {
1973  silent_cerr("StructMappingExtForce(" << uLabel << "): "
1974  "StructNode(" << uL << ") out of order "
1975  "at line " << HP.GetLineData() << std::endl);
1977  }
1978  Got[uL] = true;
1979 
1981  ncd.pNode[0] = 0;
1982 
1983  ReferenceFrame RF(pNode);
1984  if (HP.IsKeyWord("membrane")) {
1985  bMembrane = true;
1986 
1987  for (unsigned n = 0; n < 4; n++) {
1988  const StructDispNode *pNn = pDM->ReadNode<const StructDispNode, Node::STRUCTURAL>(HP);
1989  if ((n%2) && pNn == ncd.pNode[n - 1]) {
1990  silent_cerr("StructMappingExtForce(" << uLabel << "): "
1991  "nodes #" << n << " and #" << n - 1 << " are the same in \"membrane\" mapping for StructNode(" << uL << ") at line " << HP.GetLineData() << std::endl);
1993  }
1994  ncd.pNode[n] = pNn;
1995  }
1996 
1997  // first node
1998  Nodes[p] = pNode;
1999  Offsets[p] = ::Zero3;
2000 
2001  doublereal h1 = HP.GetReal();
2002  if (h1 == 0.) {
2003  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2004  "invalid h1 in \"membrane\" mapping for StructNode(" << uL << ") at line " << HP.GetLineData() << std::endl);
2006  }
2007  ncd.h1 = h1;
2008 
2009  if (bLabels) {
2010  int l = HP.GetInt();
2011  if (l < 0) {
2012  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2013  "invalid (negative) point label " << l
2014  << " at line " << HP.GetLineData() << std::endl);
2016  }
2017  Labels[p] = l;
2018  }
2019 
2020  p++;
2021 
2022  if (p == unsigned(nPoints)) {
2023  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2024  "second point in \"membrane\" mapping for StructNode(" << uL << ") exceeds points number " << nPoints << " at line " << HP.GetLineData() << std::endl);
2026  }
2027 
2028  // second node
2029  Nodes[p] = pNode;
2030  Offsets[p] = ::Zero3;
2031 
2032  doublereal h2 = HP.GetReal();
2033  if (h2 == 0. || h1*h2 > 0.) {
2034  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2035  "invalid h2 in \"membrane\" mapping for StructNode(" << uL << ") at line " << HP.GetLineData() << std::endl);
2037  }
2038  ncd.h2 = h2;
2039 
2040  if (bLabels) {
2041  int l = HP.GetInt();
2042  if (l < 0) {
2043  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2044  "invalid (negative) point label " << l
2045  << " at line " << HP.GetLineData() << std::endl);
2047  }
2048  Labels[p] = l;
2049  }
2050 
2051  p++;
2052 
2053  } else {
2054  while (HP.IsKeyWord("offset")) {
2055  if (p == unsigned(nPoints)) {
2056  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2057  "point " << p << " offset from StructNode(" << pNode->GetLabel() << ") exceeds expected value " << nPoints
2058  << " at line " << HP.GetLineData() << std::endl);
2060  }
2061 
2062  if (bLabels) {
2063  int l = HP.GetInt();
2064  if (l < 0) {
2065  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2066  "invalid (negative) point label " << l
2067  << " at line " << HP.GetLineData() << std::endl);
2069  }
2070  Labels[p] = l;
2071  }
2072 
2073  Nodes[p] = pNode;
2074  Offsets[p] = HP.GetPosRel(RF);
2075  p++;
2076  }
2077  }
2078 
2079  NodesConn.push_back(ncd);
2080  }
2081 
2082  if (HP.IsKeyWord("echo")) {
2083  const char *s = HP.GetFileName();
2084  if (s == NULL) {
2085  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2086  "unable to parse echo file name "
2087  "at line " << HP.GetLineData()
2088  << std::endl);
2090  }
2091 
2092  std::ofstream out(s);
2093 
2094  out.setf(std::ios::scientific);
2095 
2096  bool bGotSurface(false);
2097  bool bGotOutput(false);
2098  bool bGotOrder(false);
2099  bool bGotBaseNode(false);
2100  bool bGotWeight(false);
2101 
2102  bool bGotPrecision(false);
2103 
2104  /*
2105  surface: basicgrid.dat
2106  output: blade1H.dat
2107  order: 2
2108  basenode: 12
2109  weight: 2
2110 
2111  */
2112 
2113  std::string surface;
2114  std::string output;
2115  int order;
2116  int basenode;
2117  int weight;
2118  bool bWeightInf(false);
2119 
2120  while (true) {
2121  if (HP.IsKeyWord("precision")) {
2122  if (bGotPrecision) {
2123  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2124  "\"precision\" already specified "
2125  "at line " << HP.GetLineData() << std::endl);
2127  }
2128  bGotPrecision = true;
2129 
2130  int iPrecision = HP.GetInt();
2131  if (iPrecision <= 0) {
2132  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2133  "invalid echo precision " << iPrecision
2134  << " at line " << HP.GetLineData()
2135  << std::endl);
2137  }
2138  out.precision(iPrecision);
2139 
2140  } else if (HP.IsKeyWord("surface")) {
2141  if (bGotSurface) {
2142  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2143  "\"surface\" already specified "
2144  "at line " << HP.GetLineData() << std::endl);
2146  }
2147  bGotSurface = true;
2148 
2149  surface = HP.GetFileName();
2150  if (surface.empty()) {
2151  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2152  "invalid \"surface\" "
2153  "at line " << HP.GetLineData() << std::endl);
2155  }
2156 
2157  } else if (HP.IsKeyWord("output")) {
2158  if (bGotOutput) {
2159  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2160  "\"output\" already specified "
2161  "at line " << HP.GetLineData() << std::endl);
2163  }
2164  bGotOutput = true;
2165 
2166  output = HP.GetFileName();
2167  if (output.empty()) {
2168  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2169  "invalid \"output\" "
2170  "at line " << HP.GetLineData() << std::endl);
2172  }
2173 
2174  } else if (HP.IsKeyWord("order")) {
2175  if (bGotOrder) {
2176  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2177  "\"order\" already specified "
2178  "at line " << HP.GetLineData() << std::endl);
2180  }
2181  bGotOrder = true;
2182 
2183  order = HP.GetInt();
2184  if (order < 1 || order > 3) {
2185  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2186  "invalid order=" << order << ", must be 1 <= order <= 3 "
2187  "at line " << HP.GetLineData() << std::endl);
2189  }
2190 
2191  } else if (HP.IsKeyWord("basenode")) {
2192  if (bGotBaseNode) {
2193  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2194  "\"basenode\" already specified "
2195  "at line " << HP.GetLineData() << std::endl);
2197  }
2198  bGotBaseNode = true;
2199 
2200  basenode = HP.GetInt();
2201  if (basenode <= 0) {
2202  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2203  "invalid basenode=" << basenode << ", must be > 0 "
2204  "at line " << HP.GetLineData() << std::endl);
2206  }
2207 
2208  } else if (HP.IsKeyWord("weight")) {
2209  if (bGotWeight) {
2210  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2211  "\"weight\" already specified "
2212  "at line " << HP.GetLineData() << std::endl);
2214  }
2215  bGotWeight = true;
2216 
2217  if (HP.IsKeyWord("inf")) {
2218  bWeightInf = true;
2219 
2220  } else {
2221  weight = HP.GetInt();
2222  if (weight < -2) {
2223  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2224  "invalid weight=" << weight << ", must be >= -2 or \"inf\" "
2225  "at line " << HP.GetLineData() << std::endl);
2227  }
2228  }
2229 
2230  } else {
2231  break;
2232  }
2233  }
2234 
2235  time_t t = std::time(NULL);
2236  const char *user = std::getenv("USER");
2237  const char *host = std::getenv("HOSTNAME");
2238  if (user == 0) user = "nobody";
2239  if (host == 0) host = "unknown";
2240 
2241  out
2242  << "# Generated by MBDyn StructMappingExtForce(" << uLabel << ")" << std::endl
2243  << "# " << user << "@" << host << std::endl
2244  << "# " << std::ctime(&t)
2245  << "# labels: " << (bLabels ? "on" : "off") << std::endl;
2246  Vec3 xRef(::Zero3);
2247  Mat3x3 RRef(::Eye3);
2248  if (pRefNode) {
2249  xRef = pRefNode->GetXCurr();
2250  RRef = pRefNode->GetRCurr();
2251 
2252  out
2253  << "# reference: " << pRefNode->GetLabel() << std::endl
2254  << "# position: " << xRef << std::endl
2255  << "# orientation: " << MatR2EulerAngles123(RRef)*dRaDegr << std::endl;
2256  }
2257 
2258  out
2259  << "# points: " << nPoints << std::endl;
2260 
2261  if (bGotSurface) {
2262  out << "# surface: " << surface << std::endl;
2263  }
2264 
2265  if (bGotOutput) {
2266  out << "# output: " << output << std::endl;
2267  }
2268 
2269  if (bGotOrder) {
2270  out << "# order: " << order << std::endl;
2271  }
2272 
2273  if (bGotBaseNode) {
2274  out << "# basenode: " << basenode << std::endl;
2275  }
2276 
2277  if (bGotWeight) {
2278  if (bWeightInf) {
2279  out << "# weight: inf" << std::endl;
2280  } else {
2281  out << "# weight: " << weight << std::endl;
2282  }
2283  }
2284 
2285  for (unsigned n = 0, p = 0; p < unsigned(nPoints); p++) {
2286  if (p > 0 && Nodes[p] != Nodes[p - 1]) {
2287  n++;
2288  }
2289 
2290  if (bLabels) {
2291  out << Labels[p] << " ";
2292  }
2293 
2294  Vec3 x;
2295  const StructNode *pNode(dynamic_cast<const StructNode *>(Nodes[p]));
2296  if (pNode != 0) {
2297  x = RRef.MulTV(pNode->GetXCurr() + pNode->GetRCurr()*Offsets[p] - xRef);
2298 
2299  } else {
2300  Vec3 h(::Zero3);
2301  if (NodesConn[n].pNode[0] != 0) {
2302  Vec3 e1(NodesConn[n].pNode[1]->GetXCurr() - NodesConn[n].pNode[0]->GetXCurr());
2303  Vec3 e2(NodesConn[n].pNode[3]->GetXCurr() - NodesConn[n].pNode[2]->GetXCurr());
2304  Vec3 e3(e1.Cross(e2));
2305  e3 /= e3.Norm();
2306 
2307  if (Nodes[p + 1] == Nodes[p]) {
2308  h = e3*NodesConn[n].h1;
2309 
2310  } else {
2311  h = e3*NodesConn[n].h2;
2312  }
2313  }
2314 
2315  x = RRef.MulTV(Nodes[p]->GetXCurr() + h - xRef);
2316  }
2317  out << x << std::endl;
2318  }
2319 
2320  if (HP.IsKeyWord("stop")) {
2321  throw NoErr(MBDYN_EXCEPT_ARGS);
2322  }
2323  }
2324 
2325  SpMapMatrixHandler *pH = 0;
2326  std::vector<uint32_t> MappedLabels;
2327  if (HP.IsKeyWord("mapped" "points" "number")) {
2328  int nMappedPoints = 0;
2329  if (HP.IsKeyWord("from" "file")) {
2330  nMappedPoints = -1;
2331 
2332  } else {
2333  nMappedPoints = HP.GetInt();
2334  if (nMappedPoints <= 0) {
2335  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2336  "invalid mapped points number " << nMappedPoints
2337  << " at line " << HP.GetLineData() << std::endl);
2339  }
2340 
2341  nMappedPoints *= 3;
2342  }
2343 
2344  integer nCols = 3*nPoints;
2345  pH = ReadSparseMappingMatrix(HP, nMappedPoints, nCols);
2346  ASSERT((nMappedPoints%3) == 0);
2347  ASSERT(nCols == 3*nPoints);
2348  nMappedPoints /= 3;
2349 
2350  if (bLabels) {
2351  MappedLabels.resize(nMappedPoints);
2352  if (HP.IsKeyWord("mapped" "labels" "file")) {
2353  const char *sFileName = HP.GetFileName();
2354  if (sFileName == 0) {
2355  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2356  "unable to read mapped labels file name "
2357  "at line " << HP.GetLineData() << std::endl);
2359  }
2360 
2361  std::ifstream in(sFileName);
2362  if (!in) {
2363  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2364  "unable to open mapped labels file "
2365  "\"" << sFileName << "\" "
2366  "at line " << HP.GetLineData() << std::endl);
2368  }
2369 
2370  char c = in.get();
2371  while (c == '#') {
2372  do {
2373  c = in.get();
2374  } while (c != '\n');
2375  c = in.get();
2376  }
2377  in.putback(c);
2378 
2379  for (unsigned l = 0; l < unsigned(nMappedPoints); l++) {
2380  int i;
2381  in >> i;
2382  if (!in) {
2383  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2384  "unable to read mapped label #" << l << "/" << nMappedPoints
2385  << " from mapped labels file \"" << sFileName << "\"" << std::endl);
2387  }
2388  if (i < 0) {
2389  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2390  "invalid (negative) mapped label #" << l << "/" << nMappedPoints
2391  << " from mapped labels file \"" << sFileName << "\"" << std::endl);
2393  }
2394  MappedLabels[l] = i;
2395  }
2396 
2397  } else {
2398  for (unsigned l = 0; l < unsigned(nMappedPoints); l++) {
2399  int i = HP.GetInt();
2400  if (i < 0) {
2401  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2402  "invalid (negative) mapped label #" << l << "/" << nMappedPoints
2403  << std::endl);
2405  }
2406  MappedLabels[l] = i;
2407  }
2408  }
2409 
2410  int duplicate = 0;
2411  for (unsigned l = 1; l < unsigned(nMappedPoints); l++) {
2412  for (unsigned c = 0; c < l; c++) {
2413  if (MappedLabels[l] == MappedLabels[c]) {
2414  duplicate++;
2415  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2416  "duplicate mapped label " << MappedLabels[l] << ": "
2417  "#" << l << "==#" << c << std::endl);
2418  }
2419  }
2420  }
2421 
2422  if (duplicate) {
2423  silent_cerr("StructMappingExtForce(" << uLabel << "): "
2424  << duplicate << " duplicate mapped labels"
2425  << std::endl);
2427  }
2428  }
2429  }
2430 
2431  flag fOut = pDM->fReadOutput(HP, Elem::FORCE);
2432  Elem *pEl = 0;
2433  if (bMembrane) {
2435  StructMembraneMappingExtForce(uLabel, pDM, pRefNode,
2436  bUseReferenceNodeForces, bRotateReferenceNodeForces,
2437  Nodes, Offsets, Labels, NodesConn,
2438  pH,
2439  MappedLabels,
2440  bLabels, bOutputAccelerations, uRRot,
2441  pEFH, bSendAfterPredict, iCoupling, fOut));
2442 
2443  } else {
2445  StructMappingExtForce(uLabel, pDM, pRefNode,
2446  bUseReferenceNodeForces, bRotateReferenceNodeForces,
2447  Nodes, Offsets, Labels,
2448  pH,
2449  MappedLabels,
2450  bLabels, bOutputAccelerations, uRRot,
2451  pEFH, bSendAfterPredict, iCoupling, fOut));
2452  }
2453 
2454  return pEl;
2455 }
2456 
flag fReadOutput(MBDynParser &HP, const T &t) const
Definition: dataman.h:1064
virtual std::istream * GetInStream(void)
Definition: extforce.cc:78
const Vec3 Zero3(0., 0., 0.)
Vec3 Cross(const Vec3 &v) const
Definition: matvec3.h:218
STLVectorHandler m_xPP
Definition: strmappingext.h:86
long int flag
Definition: mbdyn.h:43
STLVectorHandler m_qPP
Definition: strmappingext.h:89
virtual bool bToBeOutput(void) const
Definition: output.cc:890
static void output(const LoadableElem *pEl, OutputHandler &OH)
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Definition: matvec3.h:98
virtual int GetInFileDes(void)
Definition: extforce.cc:99
virtual int GetOutFileDes(void)
Definition: extforce.cc:85
virtual void ResizeReset(integer)
Definition: vh.cc:55
void Recv(void)
Definition: extforce.cc:798
virtual integer GetInt(integer iDefval=0)
Definition: parser.cc:1050
virtual const Mat3x3 & GetRCurr(void) const
Definition: strnode.h:1012
doublereal Norm(void) const
Definition: matvec3.h:263
virtual void RecvFromStream(std::istream &inf)
virtual ~StructMappingExtForce(void)
virtual void GetConnectedNodes(std::vector< const Node * > &connectedNodes) const
Definition: mbc.h:65
integer iGetNumRows(void) const
Definition: spmh.h:113
const Mat3x3 Eye3(1., 0., 0., 0., 1., 0., 0., 0., 1.)
virtual const char * GetFileName(enum Delims Del=DEFAULTDELIM)
Definition: parsinc.cc:673
StructMappingExtForce(unsigned int uL, DataManager *pDM, const StructNode *pRefNode, bool bUseReferenceNodeForces, bool bRotateReferenceNodeForces, std::vector< const StructDispNode * > &Nodes, std::vector< Vec3 > &Offsets, std::vector< unsigned > &Labels, SpMapMatrixHandler *pH, std::vector< uint32_t > &MappedLabels, bool bLabels, bool bOutputAccelerations, unsigned uRRot, ExtFileHandlerBase *pEFH, bool bSendAfterPredict, int iCoupling, flag fOut)
std::vector< uint32_t > m_qlabels
Definition: strmappingext.h:82
void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
#define NO_OP
Definition: myassert.h:74
Definition: mbc.h:71
STLVectorHandler m_x
Definition: strmappingext.h:84
std::vector< NodeData > Nodes
Definition: strmappingext.h:73
Vec3 VecRot(const Mat3x3 &Phi)
Definition: Rot.cc:136
virtual void SendToStream(std::ostream &outf, ExtFileHandlerBase::SendWhen when)
virtual void SendToFileDes(int outfd, ExtFileHandlerBase::SendWhen when)
Vec3 MulTV(const Vec3 &v) const
Definition: matvec3.cc:482
virtual Negotiate NegotiateRequest(void) const
Definition: extforce.cc:63
virtual void PutRowIndex(integer iSubRow, integer iRow)=0
bool Prepare(ExtFileHandlerBase *pEFH)
enum @55 order
integer iGetNumCols(void) const
Definition: spmh.h:117
static int labels
virtual void RecvFromFileDes(int infd)
Definition: mbc.h:66
Vec3 GetPosRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1331
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
virtual void Put(integer iRow, const Vec3 &v)
Definition: stlvh.cc:181
const char * host
Definition: autopilot.c:142
virtual int GetRecvFlags(void) const
Definition: extforce.cc:106
ExtFileHandlerBase * pEFH
Definition: extforce.h:193
virtual void Output(OutputHandler &OH) const
Vec3 MatR2EulerAngles123(const Mat3x3 &R)
Definition: matvec3.cc:893
virtual integer iGetFirstMomentumIndex(void) const =0
virtual void SendToStream(std::ostream &outf, ExtFileHandlerBase::SendWhen when)
virtual const Vec3 & GetWCurr(void) const
Definition: strnode.h:1030
const doublereal dRaDegr
Definition: matvec3.cc:884
virtual const Vec3 & GetWPCurr(void) const
Definition: strnode.h:1042
#define ASSERT(expression)
Definition: colamd.c:977
virtual std::ostream * GetOutStream(void)
Definition: extforce.cc:72
Mat3x3 MulTM(const Mat3x3 &m) const
Definition: matvec3.cc:500
VectorExpression< VectorCrossExpr< VectorLhsExpr, VectorRhsExpr >, 3 > Cross(const VectorExpression< VectorLhsExpr, 3 > &u, const VectorExpression< VectorRhsExpr, 3 > &v)
Definition: matvec.h:3248
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
virtual const Vec3 & GetXCurr(void) const
Definition: strnode.h:310
virtual void Add(integer iRow, const Vec3 &v)
Definition: vh.cc:63
std::vector< NodeConnData > NodesConn
Definition: except.h:79
SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
static int nodes
static std::stack< cleanup * > c
Definition: cleanup.cc:59
const doublereal * pGetMat(void) const
Definition: matvec3.h:743
STLVectorHandler m_q
Definition: strmappingext.h:87
virtual void RecvFromStream(std::istream &inf)
virtual bool GetYesNo(bool &bRet)
Definition: parser.cc:1022
virtual bool IsArg(void)
Definition: parser.cc:807
Definition: elem.h:75
#define MBC_U_ROT_2_REF_NODE_ROT(u)
Definition: mbc.h:248
virtual ~StructMembraneMappingExtForce(void)
const doublereal * pGetVec(void) const
Definition: matvec3.h:192
void Send(ExtFileHandlerBase *pEFH, ExtFileHandlerBase::SendWhen when)
virtual const Vec3 & GetVCurr(void) const
Definition: strnode.h:322
STLVectorHandler m_p
Definition: strmappingext.h:92
StructMembraneMappingExtForce(unsigned int uL, DataManager *pDM, const StructNode *pRefNode, bool bUseReferenceNodeForces, bool bRotateReferenceNodeForces, std::vector< const StructDispNode * > &Nodes, std::vector< Vec3 > &Offsets, std::vector< unsigned > &Labels, std::vector< NodeConnData > &NodesConn, SpMapMatrixHandler *pH, std::vector< uint32_t > &MappedLabels, bool bLabels, bool bOutputAccelerations, unsigned uRRot, ExtFileHandlerBase *pEFH, bool bSendAfterPredict, int iCoupling, flag fOut)
SpMapMatrixHandler * pH
Definition: strmappingext.h:56
static const doublereal a
Definition: hfluid_.h:289
Elem * ReadStructMappingExtForce(DataManager *pDM, MBDynParser &HP, unsigned int uLabel)
virtual int GetSendFlags(void) const
Definition: extforce.cc:92
virtual const Vec3 & GetXPPCurr(void) const
Definition: strnode.h:334
static doublereal buf[BUFSIZE]
Definition: discctrl.cc:333
virtual VectorHandler & MatTVecMul(VectorHandler &out, const VectorHandler &in) const
Definition: mh.cc:341
virtual void RecvFromFileDes(int infd)
double doublereal
Definition: colamd.c:52
int iCoupling
Definition: extforce.h:212
long int integer
Definition: colamd.c:51
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
Definition: mbc.h:73
unsigned int GetLabel(void) const
Definition: withlab.cc:62
const StructNode * pRefNode
Definition: strmappingext.h:48
Node * ReadNode(MBDynParser &HP, Node::Type type) const
Definition: dataman3.cc:2309
virtual VectorHandler & MatVecMul(VectorHandler &out, const VectorHandler &in) const
Definition: mh.cc:332
STLVectorHandler m_xP
Definition: strmappingext.h:85
void ReadExtForce(DataManager *pDM, MBDynParser &HP, unsigned int uLabel, ExtFileHandlerBase *&pEFH, bool &bSendAfterPredict, int &iCoupling)
Definition: extforce.cc:1141
virtual void Resize(integer iNewSize)=0
std::ostream & Forces(void) const
Definition: output.h:450
Mat3x3 R
virtual void SendToFileDes(int outfd, ExtFileHandlerBase::SendWhen when)
bool UseText(int out) const
Definition: output.cc:446
STLVectorHandler m_qP
Definition: strmappingext.h:88
SpMapMatrixHandler * ReadSparseMappingMatrix(MBDynParser &HP, integer &nRows, integer &nCols)
STLVectorHandler m_f
Definition: strmappingext.h:91
virtual doublereal GetReal(const doublereal &dDefval=0.0)
Definition: parser.cc:1056