MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
beamslider.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/struct/beamslider.cc,v 1.50 2017/01/12 14:46:43 masarati Exp $ */
2 /*
3  * MBDyn (C) is a multibody analysis code.
4  * http://www.mbdyn.org
5  *
6  * Copyright (C) 1996-2017
7  *
8  * Pierangelo Masarati <masarati@aero.polimi.it>
9  * Paolo Mantegazza <mantegazza@aero.polimi.it>
10  *
11  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12  * via La Masa, 34 - 20156 Milano, Italy
13  * http://www.aero.polimi.it
14  *
15  * Changing this copyright notice is forbidden.
16  *
17  * This program is free software; you can redistribute it and/or modify
18  * it under the terms of the GNU General Public License as published by
19  * the Free Software Foundation (version 2 of the License).
20  *
21  *
22  * This program is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with this program; if not, write to the Free Software
29  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30  */
31 
32 /*
33  * Trave a volumi finiti, con predisposizione per forze di inerzia consistenti
34  * e legame cositutivo piezoelettico
35  */
36 
37 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
38 
39 #include "dataman.h"
40 #include "constltp.h"
41 #include "shapefnc.h"
42 #include "beamslider.h"
43 
44 
45 /* BeamConn - begin */
47  const Vec3& f1, const Vec3& f2, const Vec3& f3,
48  const Mat3x3& R1, const Mat3x3& R2, const Mat3x3& R3)
49 : m_pBeam(pB)
50 {
51  m_f[0] = f1;
52  m_f[1] = f2;
53  m_f[2] = f3;
54 
55  m_R[0] = R1;
56  m_R[1] = R2;
57  m_R[2] = R3;
58 }
59 
61 {
62  NO_OP;
63 }
64 
65 /* BeamConn - end */
66 
67 
68 /* BeamSliderJoint - begin */
69 
70 /* Punto di valutazione */
71 const doublereal dS = 1./sqrt(3.);
72 
73 /* Costruttore non banale */
74 BeamSliderJoint::BeamSliderJoint(unsigned int uL, const DofOwner* pDO,
75  const StructNode* pN, enum Type iT,
76  unsigned int nB, const BeamConn *const * ppB,
77  unsigned int uIB, unsigned int uIN,
78  doublereal dl,
79  const Vec3& fTmp, const Mat3x3& RTmp, flag fOut)
80 : Elem(uL, fOut),
81 Joint(uL, pDO, fOut),
82 nRotConstr(0), nBeams(nB), iCurrBeam(0), iType(iT),
83 pNode(pN), ppBeam(ppB),
84 f(fTmp), R(RTmp),
85 F(Zero3), M(Zero3),
86 sRef(0.), s(0.),
87 dL(dl),
88 x(Zero3), l(Zero3)
89 {
90  ASSERT(pNode != NULL);
91  ASSERT(nBeams > 0);
92  ASSERT(ppBeam != NULL);
93  ASSERT(uIB > 0 && uIB <= nBeams);
94  ASSERT(uIN > 0 && uIN <= 3);
95 
96  iCurrBeam = uIB-1;
97 
98  switch (uIN) {
99  case 1:
101  break;
102 
103  case 2:
105  break;
106 
107  case 3:
109  break;
110  }
111 
112  sRef = 2.*iCurrBeam + activeNode - 2.;
113  s = sRef;
114 
115  switch (iType) {
116  case CLASSIC:
117  nRotConstr = 2;
118  break;
119 
120  case SPLINE:
121  nRotConstr = 3;
122  break;
123 
124  default:
125  break;
126  }
127 
128 #ifdef DEBUG
129  for (unsigned int i = 0; i < nBeams; i++) {
130  ASSERT(ppBeam[i] != NULL);
131  }
132 #endif /* DEBUG */
133 }
134 
136 {
137  NO_OP;
138 }
139 
140 std::ostream&
141 BeamSliderJoint::Restart(std::ostream& out) const
142 {
143  return out << "# beam slider not implemented yet" << std::endl
144  << "beam slider;" << std::endl;
145 }
146 
147 void
149 {
150  if (bToBeOutput()) {
151  Mat3x3 RTmp(pNode->GetRCurr()*R);
152  Mat3x3 RTmpT(RTmp.Transpose());
153 
154  Joint::Output(OH.Joints(), "BeamSlider", GetLabel(),
155  RTmpT*F, M, F, RTmp*M)
156  << " " << ppBeam[iCurrBeam]->pGetBeam()->GetLabel()
157  << " " << sRef << " " << l << std::endl;
158  }
159 }
160 
161 #if 0
162 unsigned int
164 {
165  return 0;
166 }
167 
168 unsigned int
169 BeamSliderJoint::iGetPrivDataIdx(const char *s) const
170 {
171  return 0;
172 }
173 
174 doublereal
175 BeamSliderJoint::dGetPrivData(unsigned int i) const
176 {
177  return 0.;
178 }
179 #endif
180 
181 /* Assemblaggio jacobiano */
184  doublereal dCoef,
185  const VectorHandler& XCurr,
186  const VectorHandler& XPrimeCurr)
187 {
188  DEBUGCOUT("Entering BeamSliderJoint::AssJac()" << std::endl);
189 
190  /* Dimensiona e resetta la matrice di lavoro */
191  integer iNumRows = 0;
192  integer iNumCols = 0;
193  WorkSpaceDim(&iNumRows, &iNumCols);
194 
195  FullSubMatrixHandler& WM = WorkMat.SetFull();
196  WM.ResizeReset(iNumRows, iNumCols);
197 
198  /* Indici */
199  integer iNodeFirstMomIndex = pNode->iGetFirstMomentumIndex();
200  integer iNodeFirstPosIndex = pNode->iGetFirstPositionIndex();
201  integer iFirstReactionIndex = iGetFirstIndex();
202  const StructNode *pBeamNode[Beam::NUMNODES];
203 
204  /*
205  * 1 => 6: nodo corpo
206  * 7 => 12: nodo 1 trave
207  * 13 => 18: nodo 2 trave
208  * 19 => 24: nodo 3 trave
209  * 25: l^T F = 0 (s)
210  * 26 => 28: vincolo posizione (F)
211  * 29 => 31: vincoli rotazione, se presenti
212  */
213 
214  /* Indici dei nodi */
215  for (int iCnt = 1; iCnt <= 6; iCnt++) {
216  WM.PutRowIndex(iCnt, iNodeFirstMomIndex+iCnt);
217  WM.PutColIndex(iCnt, iNodeFirstPosIndex+iCnt);
218  }
219 
220  for (int nCnt = 1; nCnt <= Beam::NUMNODES; nCnt++) {
221  pBeamNode[nCnt-1] = ppBeam[iCurrBeam]->pGetNode(nCnt);
222  integer iBeamFirstMomIndex =
223  pBeamNode[nCnt-1]->iGetFirstMomentumIndex();
224  integer iBeamFirstPosIndex =
225  pBeamNode[nCnt-1]->iGetFirstPositionIndex();
226 
227  for (int iCnt = 1; iCnt <= 6; iCnt++) {
228  WM.PutRowIndex(6*nCnt+iCnt,
229  iBeamFirstMomIndex+iCnt);
230  WM.PutColIndex(6*nCnt+iCnt,
231  iBeamFirstPosIndex+iCnt);
232  }
233  }
234 
235  /* Indici del vincolo */
236  for (unsigned int iCnt = 1; iCnt <= iGetNumDof(); iCnt++) {
237  WM.PutRowIndex(6*(1+Beam::NUMNODES)+iCnt,
238  iFirstReactionIndex+iCnt);
239  WM.PutColIndex(6*(1+Beam::NUMNODES)+iCnt,
240  iFirstReactionIndex+iCnt);
241  }
242 
243  /* vincolo in posizione */
244  for (unsigned int i = 1; i <= 3; i++) {
245  /* l^T F = 0 : Delta F */
246  doublereal d = l.dGet(i)/dCoef;
247  WM.DecCoef(6*(1+Beam::NUMNODES)+1,
248  6*(1+Beam::NUMNODES)+1+i, d);
249 
250  /* xc - x = 0: l Delta s */
251  WM.IncCoef(6*(1+Beam::NUMNODES)+1+i,
252  6*(1+Beam::NUMNODES)+1, d);
253 
254  /* xc - x = 0: Delta x_b */
255  WM.DecCoef(6*(1+Beam::NUMNODES)+1+i, i, 1.);
256  }
257 
258  Vec3 lp(Zero3);
259  for (unsigned int iN = 0; iN < Beam::NUMNODES; iN++) {
260  Vec3 Tmp(fTmp[iN].Cross(F));
261 
262  /* l^T F = 0 : Delta s */
263  lp += xTmp[iN]*dNpp[iN];
264 
265  for (unsigned int i = 1; i <= 3; i++) {
266  /* l^T F = 0 : Delta x */
267  doublereal d = F.dGet(i)*dNp[iN];
268  WM.DecCoef(6*(1+Beam::NUMNODES)+1,
269  6*(1+iN)+i, d);
270 
271  /* l^T F = 0 : Delta g */
272  d = Tmp.dGet(i)*dNp[iN];
273  WM.DecCoef(6*(1+Beam::NUMNODES)+1,
274  6*(1+iN)+3+i, d);
275 
276  /* xc - x = 0: Delta x */
277  WM.IncCoef(6*(1+Beam::NUMNODES)+1+i,
278  6*(1+iN)+i, dN[iN]);
279  }
280 
281  /* xc - x = 0: Delta g */
282  WM.Sub(6*(1+Beam::NUMNODES)+1+1,
283  6*(1+iN)+3+1, Mat3x3(MatCross, fTmp[iN]*dN[iN]));
284  }
285 
286  /* l^T F = 0 : Delta s */
287  WM.DecCoef(6*(1+Beam::NUMNODES)+1,
288  6*(1+Beam::NUMNODES)+1, (F*lp)/dCoef);
289 
290  /* reazioni vincolari */
291  for (unsigned int i = 1; i <= 3; i++) {
292  /* corpo: Delta F */
293  WM.DecCoef(i, 6*(1+Beam::NUMNODES)+1+i, 1.);
294 
295  /* trave: Delta F */
296  WM.IncCoef(6*activeNode+i, 6*(1+Beam::NUMNODES)+1+i, dW[0]);
297  }
298 
299  /* corpo: Delta F (momento) */
300  Mat3x3 MTmp(MatCross, fb);
301  WM.Sub(3+1, 6*(1+Beam::NUMNODES)+1+1, MTmp);
302 
303  /* vincolo posizione: Delta gb */
304  WM.Add(6*(1+Beam::NUMNODES)+1+1, 3+1, MTmp);
305 
306  /* corpo: Delta gb (momento) */
307  Mat3x3 Ffb(MatCrossCross, F, fb*dCoef);
308  WM.Sub(3+1, 3+1, Ffb);
309 
310  /* trave: Delta gb (momento) */
311  WM.Add(6*activeNode+3+1, 3+1, Ffb*dW[0]);
312 
313  /* trave: Delta F (momento) */
314  MTmp = Mat3x3(MatCross, F*(dCoef*dW[0]));
315  WM.Add(6*activeNode+3+1, 6*(1+Beam::NUMNODES)+1+1,
316  Mat3x3(MatCross, (xc - xNod[activeNode-1])*dW[0]));
317  WM.Sub(6*activeNode+3+1, 1, MTmp);
318  WM.Add(6*activeNode+3+1, 6*activeNode+1, MTmp);
319 
320  if (dW[1] != 0.) {
321  /* trave: Delta gb (momento) */
322  WM.Add(6*(activeNode+1)+3+1, 3+1, Ffb*dW[1]);
323 
324 /* NOTE: these terms are questionable: the perturbation
325  * related to the amplitude of the smearing should go
326  * in the jacobian, but there is no evidence of improvements
327  * in convergence */
328 #define DELTADW
329 #ifdef DELTADW
330  Vec3 m1(M+(xc-xNod[activeNode-1]).Cross(F));
331  Vec3 m2(M+(xc-xNod[activeNode]).Cross(F));
332 #endif /* DELTADW */
333 
334  /* reazioni vincolari */
335  for (unsigned int i = 1; i <= 3; i++) {
336  /* trave: Delta F */
337  WM.IncCoef(6*(activeNode+1)+i,
338  6*(1+Beam::NUMNODES)+1+i, dW[1]);
339 
340 #ifdef DELTADW
341  /* trave: Delta s (Delta dW forza) */
342  doublereal d = F(i)/(2.*dL);
343  WM.DecCoef(6*activeNode+i,
344  6*(1+Beam::NUMNODES)+1, d);
345  WM.IncCoef(6*(activeNode+1)+i,
346  6*(1+Beam::NUMNODES)+1, d);
347 
348  /* trave: Delta s (Delta dW momento) */
349  d = m1(i)/(2.*dL);
350  WM.DecCoef(6*activeNode+3+i,
351  6*(1+Beam::NUMNODES)+1, d);
352  d = m2(i)/(2.*dL);
353  WM.IncCoef(6*(activeNode+1)+3+i,
354  6*(1+Beam::NUMNODES)+1, d);
355 #endif /* DELTADW */
356  }
357 
358  /* trave: Delta F (momento) */
359  Mat3x3 MTmp(MatCross, F*(dCoef*dW[1]));
360  WM.Add(6*(activeNode+1)+3+1, 6*(1+Beam::NUMNODES)+1+1,
361  Mat3x3(MatCross, (xc-xNod[activeNode]))*dW[1]);
362  WM.Sub(6*(activeNode+1)+3+1, 1, MTmp);
363  WM.Add(6*(activeNode+1)+3+1, 6*(activeNode+1)+1, MTmp);
364  }
365 
366  /* Vincolo in rotazione */
367  if (iType != SPHERICAL) {
368  Vec3 eb2 = Rb.GetVec(2);
369  Vec3 eb3 = Rb.GetVec(3);
370 
371  Vec3 mm(eb2*m(2) + eb3*m(3));
372 
373  doublereal d;
374 
375  for (unsigned int iN = 0; iN < Beam::NUMNODES; iN++) {
376  Vec3 Tmpf2(fTmp[iN].Cross(eb2));
377  Vec3 Tmpf3(fTmp[iN].Cross(eb3));
378 
379  for (unsigned int i = 1; i <= 3; i++) {
380  /* Vincolo in rotazione: Delta x */
381  d = eb2.dGet(i)*dNp[iN];
382  WM.DecCoef(6*(1+Beam::NUMNODES)+1+3+1,
383  6*(1+iN)+i, d);
384 
385  d = eb3.dGet(i)*dNp[iN];
386  WM.DecCoef(6*(1+Beam::NUMNODES)+1+3+2,
387  6*(1+iN)+i, d);
388 
389  /* Vincolo in rotazione: Delta g */
390  d = Tmpf2.dGet(i)*dNp[iN];
391  WM.DecCoef(6*(1+Beam::NUMNODES)+1+3+1,
392  6*(1+iN)+3+i, d);
393 
394  d = Tmpf3.dGet(i)*dNp[iN];
395  WM.DecCoef(6*(1+Beam::NUMNODES)+1+3+2,
396  6*(1+iN)+3+i, d);
397 
398  }
399 
400  Vec3 mmTmp(mm*(dNp[iN]*dCoef));
401  Mat3x3 mmTmp2(MatCross, mmTmp);
402  Mat3x3 mmTmp3(MatCrossCross, mmTmp, fTmp[iN]);
403 
404 #if 0
405  Vec3 MTmp(M*(dNp[iN]*dCoef));
406  Mat3x3 MTmp2(MTmp);
407  Mat3x3 MTmp3(MTmp, fTmp[iN]);
408 #endif
409 
410  /* Reazione vincolare: Delta x */
411  WM.Sub(3+1, 6*(1+iN)+1, mmTmp2);
412 
413  /* Reazione vincolare: Delta g */
414  WM.Add(3+1, 6*(1+iN)+3+1, mmTmp3);
415 
416  if (dW[1] == 0.) {
417  /* Reazione vincolare: Delta x */
418  WM.Add(6*activeNode+3+1, 6*(1+iN)+1, mmTmp2);
419 
420  /* Reazione vincolare: Delta g */
421  WM.Sub(6*activeNode+3+1, 6*(1+iN)+3+1, mmTmp3);
422  } else {
423  /* Reazione vincolare: Delta x */
424  WM.Add(6*activeNode+3+1, 6*(1+iN)+1,
425  mmTmp2*dW[0]);
426  WM.Add(6*(activeNode+1)+3+1, 6*(1+iN)+1,
427  mmTmp2*dW[1]);
428 
429  /* Reazione vincolare: Delta g */
430  WM.Sub(6*activeNode+3+1, 6*(1+iN)+3+1,
431  mmTmp3*dW[0]);
432  WM.Sub(6*(activeNode+1)+3+1, 6*(1+iN)+3+1,
433  mmTmp3*dW[1]);
434  }
435  }
436 
437  Vec3 Tmpl2(eb2.Cross(l));
438  Vec3 Tmpl3(eb3.Cross(l));
439  Vec3 Tmpmmlp(mm.Cross(lp));
440 
441  for (unsigned int i = 1; i <= 3; i++) {
442 
443  /* Vincolo in rotazione: Delta gb */
444  d = Tmpl2.dGet(i);
445  WM.DecCoef(6*(1+Beam::NUMNODES)+1+3+1, 3+i, d);
446 
447  /* Reazione vincolare: Delta M */
448  WM.DecCoef(3+i, 6*(1+Beam::NUMNODES)+1+3+1, d);
449  WM.IncCoef(6*activeNode+3+i,
450  6*(1+Beam::NUMNODES)+1+3+1, d*dW[0]);
451  if (dW[1] != 0) {
452  WM.IncCoef(6*(activeNode+1)+3+i,
453  6*(1+Beam::NUMNODES)+1+3+1,
454  d*dW[1]);
455  }
456 
457  /* Vincolo in rotazione: Delta gb */
458  d = Tmpl3.dGet(i);
459  WM.DecCoef(6*(1+Beam::NUMNODES)+1+3+2, 3+i, d);
460 
461  /* Reazione vincolare: Delta M */
462  WM.DecCoef(3+i, 6*(1+Beam::NUMNODES)+1+3+2, d);
463  WM.IncCoef(6*activeNode+3+i,
464  6*(1+Beam::NUMNODES)+1+3+2, d*dW[0]);
465  if (dW[1] != 0) {
466  WM.IncCoef(6*(activeNode+1)+3+i,
467  6*(1+Beam::NUMNODES)+1+3+2,
468  d*dW[1]);
469  }
470 
471  /* Reazione vincolare: Delta s */
472  d = Tmpmmlp(i);
473  WM.DecCoef(3+i, 6*(1+Beam::NUMNODES)+1, d);
474  WM.IncCoef(6*activeNode+3+i,
475  6*(1+Beam::NUMNODES)+1, d*dW[0]);
476  if (dW[1] != 0) {
477  WM.IncCoef(6*(activeNode+1)+3+i,
478  6*(1+Beam::NUMNODES)+1,
479  d*dW[1]);
480  }
481  }
482 
483  /* Vincolo in rotazione: Delta s */
484  d = (eb2*lp)/dCoef;
485  WM.DecCoef(6*(1+Beam::NUMNODES)+1+3+1,
486  6*(1+Beam::NUMNODES)+1, d);
487 
488  d = (eb3*lp)/dCoef;
489  WM.DecCoef(6*(1+Beam::NUMNODES)+1+3+2,
490  6*(1+Beam::NUMNODES)+1, d);
491 
492  /* Reazione vincolare: Delta gb */
493  Mat3x3 mmTmp(MatCrossCross, l, mm*dCoef);
494  WM.Sub(3+1, 3+1, mmTmp);
495  if (dW[1] == 0) {
496  WM.Add(6*activeNode+3+1, 3+1, mmTmp);
497  } else {
498  WM.Add(6*activeNode+3+1, 3+1, mmTmp*dW[0]);
499  WM.Add(6*(activeNode+1)+3+1, 3+1, mmTmp*dW[1]);
500  }
501 
502  }
503 
504  return WorkMat;
505 }
506 
507 /* Assemblaggio residuo */
510  doublereal dCoef,
511  const VectorHandler& XCurr,
512  const VectorHandler& /* XPrimeCurr */ )
513 {
514  DEBUGCOUT("Entering BeamSliderJoint::AssRes()" << std::endl);
515 
516  /*
517  * Nota: posso risparmiare tutte le righe dei nodi della trave
518  * che non sono attivi ...
519  */
520 
521  /* Dimensiona e resetta la matrice di lavoro */
522  integer iNumRows = 0;
523  integer iNumCols = 0;
524  WorkSpaceDim(&iNumRows, &iNumCols);
525  WorkVec.ResizeReset(iNumRows);
526 
527  /* Indici */
528  integer iNodeFirstMomIndex = pNode->iGetFirstMomentumIndex();
529  integer iFirstReactionIndex = iGetFirstIndex();
530  const StructNode *pBeamNode[Beam::NUMNODES];
531 
532  /* Aggiorna i dati propri */
533  sRef = XCurr(iFirstReactionIndex+1);
534  F = Vec3(XCurr, iFirstReactionIndex+1+1);
535  switch (iType) {
536  /*
537  * m(2), m(3) are the moments about the axes
538  * orthogonal to the tangent to the reference line
539  */
540  case CLASSIC:
541  m.Put(2, XCurr(iFirstReactionIndex+1+3+1));
542  m.Put(3, XCurr(iFirstReactionIndex+1+3+2));
543  break;
544 
545  /*
546  * m(1) is the moment about the tangent to the
547  * reference line;
548  * m(2), m(3) are the moments about the axes
549  * orthogonal to the tangent to the reference line
550  */
551  case SPLINE:
552  m = Vec3(XCurr, iFirstReactionIndex+1+3+1);
553  break;
554 
555  /*
556  * No moment
557  */
558  default:
559  /* M is set to zero by someone else ... */
560  break;
561  }
562 
563  /*
564  * in base al valore di s decide su quale trave sta operando
565  * (da studiare e implementare ...)
566  *
567  * Nota: passando da una trave all'altra non e' detto che la
568  * metrica sia la stessa (se hanno lunghezze diverse o i nodi
569  * non sono equispaziati, cambia).
570  * In prima approssimazione faccio finta che sia la stessa;
571  * un raffinamento si potra' avere considerando il rapporto
572  * tra le metriche.
573  */
574  s = sRef - 2*iCurrBeam;
575  if (s < -1.) {
576  /* passo alla trave precedente */
577  if (iCurrBeam > 0) {
578  s += 2.;
579  iCurrBeam--;
580  }
581 
582  } else if (s > 1.) {
583  /* passo alla trave successiva */
584  if (iCurrBeam < nBeams-1) {
585  s -= 2.;
586  iCurrBeam++;
587  }
588  }
589 
590  /* Cerco il tratto di trave a cui le forze si applicano ... */
591  /* Primo tratto */
592  if (s < -dS - dL) {
593  activeNode = 1;
594 
595  dW[0] = 1.;
596  dW[1] = 0.;
597 
598  } else if ( s < -dS + dL) {
599  activeNode = 1;
600 
601  doublereal d = .5*(dS + s)/dL;
602  dW[0] = .5 - d;
603  dW[1] = .5 + d;
604 
605  /* Ultimo tratto */
606  } else if (s > dS + dL) {
607  activeNode = 3;
608 
609  dW[0] = 1.;
610  dW[1] = 0.;
611 
612  } else if (s > dS - dL) {
613  activeNode = 2;
614 
615  doublereal d = .5*(dS - s)/dL;
616  dW[0] = .5 + d;
617  dW[1] = .5 - d;
618 
619  /* Tratto centrale */
620  } else {
621  activeNode = 2;
622 
623  dW[0] = 1.;
624  dW[1] = 0.;
625  }
626 
627  /* Indici dei nodi */
628  for (int iCnt = 1; iCnt <= 6; iCnt++) {
629  WorkVec.PutRowIndex(iCnt, iNodeFirstMomIndex+iCnt);
630  }
631  for (int nCnt = 1; nCnt <= Beam::NUMNODES; nCnt++) {
632  pBeamNode[nCnt-1] = ppBeam[iCurrBeam]->pGetNode(nCnt);
633  integer iBeamFirstMomIndex =
634  pBeamNode[nCnt-1]->iGetFirstMomentumIndex();
635 
636  for (int iCnt = 1; iCnt <= 6; iCnt++) {
637  WorkVec.PutRowIndex(6*nCnt+iCnt,
638  iBeamFirstMomIndex+iCnt);
639  }
640  }
641 
642  /* Indici del vincolo */
643  for (unsigned int iCnt = 1; iCnt <= iGetNumDof(); iCnt++) {
644  WorkVec.PutRowIndex(6*(1+Beam::NUMNODES)+iCnt,
645  iFirstReactionIndex+iCnt);
646  }
647 
648  /*
649  * Recupero dati
650  */
651  x = Vec3(Zero3);
652  l = Vec3(Zero3);
653  for (unsigned int i = 0; i < Beam::NUMNODES; i++) {
654  xNod[i] = pBeamNode[i]->GetXCurr();
655  fTmp[i] = pBeamNode[i]->GetRCurr()*ppBeam[iCurrBeam]->Getf(i+1);
656  xTmp[i] = xNod[i]+fTmp[i];
657 
658  dN[i] = ShapeFunc3N(s, i+1);
659  dNp[i] = ShapeFunc3N(s, i+1, ORD_D1);
660  dNpp[i] = ShapeFunc3N(s, i+1, ORD_D2);
661  x += xTmp[i]*dN[i];
662  l += xTmp[i]*dNp[i];
663  }
664 
665  Rb = pNode->GetRCurr()*R;
666  fb = pNode->GetRCurr()*f;
667  xc = pNode->GetXCurr()+fb;
668 
669  Vec3 eb2 = Rb.GetVec(2);
670  Vec3 eb3 = Rb.GetVec(3);
671 
672  /*
673  * vincoli di posizione
674  *
675  * FIXME: togliere la scalatura da F*l ?????
676  */
677  WorkVec.PutCoef(6*(1+Beam::NUMNODES)+1, (F*l)/dCoef);
678  WorkVec.Add(6*(1+Beam::NUMNODES)+1+1, (xc-x)/dCoef);
679 
680  /*
681  * Vincoli di rotazione
682  */
683  if (iType != SPHERICAL) {
684  /* 2 vincoli di rotazione */
685  WorkVec.PutCoef(6*(1+Beam::NUMNODES)+1+3+1, (eb2*l)/dCoef);
686  WorkVec.PutCoef(6*(1+Beam::NUMNODES)+1+3+2, (eb3*l)/dCoef);
687 
688  /* calcolo momento */
689  M = eb2.Cross(l*m(2))+eb3.Cross(l*m(3));
690 
691  if (iType == SPLINE) {
692  /* FIXME: vincolo spline */
693  NO_OP;
694  }
695  }
696 
697  /*
698  * reazioni vincolari
699  */
700  WorkVec.Add(1, F);
701  WorkVec.Add(3+1, M+fb.Cross(F));
702 
703  WorkVec.Sub(6*activeNode+1, F*dW[0]);
704  WorkVec.Sub(6*activeNode+3+1,
705  (M+(xc-xNod[activeNode-1]).Cross(F))*dW[0]);
706 
707  if (dW[1] != 0.) {
708  WorkVec.Sub(6*(activeNode+1)+1, F*dW[1]);
709  WorkVec.Sub(6*(activeNode+1)+3+1,
710  (M+(xc-xNod[activeNode]).Cross(F))*dW[1]);
711  }
712 
713  return WorkVec;
714 }
715 
716 /* Contributo allo jacobiano durante l'assemblaggio iniziale */
719  VariableSubMatrixHandler& WorkMat,
720  const VectorHandler& XCurr
721 )
722 {
723  WorkMat.SetNullMatrix();
724 
725  return WorkMat;
726 }
727 
728 /* Contributo al residuo durante l'assemblaggio iniziale */
731  SubVectorHandler& WorkVec,
732  const VectorHandler& XCurr
733 )
734 {
735  WorkVec.Resize(0);
736 
737  return WorkVec;
738 }
739 
740 /* BeamSliderJoint - end */
741 
virtual std::ostream & Restart(std::ostream &out) const
Definition: beamslider.cc:141
const Vec3 & Getf(unsigned int i) const
Definition: beamslider.h:78
void PutColIndex(integer iSubCol, integer iCol)
Definition: submat.h:325
const Vec3 Zero3(0., 0., 0.)
Vec3 Cross(const Vec3 &v) const
Definition: matvec3.h:218
Vec3 m_f[3]
Definition: beamslider.h:56
unsigned int nBeams
Definition: beamslider.h:100
long int flag
Definition: mbdyn.h:43
virtual bool bToBeOutput(void) const
Definition: output.cc:890
SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: beamslider.cc:509
Definition: matvec3.h:98
Vec3 xNod[Beam::NUMNODES]
Definition: beamslider.h:118
virtual void ResizeReset(integer)
Definition: vh.cc:55
const MatCross_Manip MatCross
Definition: matvec3.cc:639
FullSubMatrixHandler & SetFull(void)
Definition: submat.h:1168
virtual const Mat3x3 & GetRCurr(void) const
Definition: strnode.h:1012
void Add(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:209
unsigned int iCurrBeam
Definition: beamslider.h:101
virtual unsigned int iGetPrivDataIdx(const char *s) const
Definition: simentity.cc:142
virtual void Sub(integer iRow, const Vec3 &v)
Definition: vh.cc:78
doublereal ShapeFunc3N(doublereal d, integer iNode, enum Order Ord)
Definition: shapefnc.cc:173
#define NO_OP
Definition: myassert.h:74
doublereal dNp[Beam::NUMNODES]
Definition: beamslider.h:122
doublereal dNpp[Beam::NUMNODES]
Definition: beamslider.h:123
void IncCoef(integer iRow, integer iCol, const doublereal &dCoef)
Definition: submat.h:683
Vec3 GetVec(unsigned short int i) const
Definition: matvec3.h:893
BeamConn(const Beam *pB, const Vec3 &f1, const Vec3 &f2, const Vec3 &f3, const Mat3x3 &R1=Eye3, const Mat3x3 &R2=Eye3, const Mat3x3 &R3=Eye3)
Definition: beamslider.cc:46
void DecCoef(integer iRow, integer iCol, const doublereal &dCoef)
Definition: submat.h:694
doublereal dN[Beam::NUMNODES]
Definition: beamslider.h:121
const StructNode * pNode
Definition: beamslider.h:105
virtual void PutRowIndex(integer iSubRow, integer iRow)=0
const StructNode * pGetNode(unsigned int i) const
Definition: beamslider.h:73
void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
Definition: beamslider.h:160
VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Definition: beamslider.cc:183
void SetNullMatrix(void)
Definition: submat.h:1159
const BeamConn *const * ppBeam
Definition: beamslider.h:106
const doublereal & dGet(unsigned short int iRow) const
Definition: matvec3.h:285
Definition: beam.h:56
#define DEBUGCOUT(msg)
Definition: myassert.h:232
virtual integer iGetFirstMomentumIndex(void) const =0
virtual integer iGetFirstPositionIndex(void) const
Definition: strnode.h:452
Mat3x3 m_R[3]
Definition: beamslider.h:57
std::ostream & Joints(void) const
Definition: output.h:443
#define ASSERT(expression)
Definition: colamd.c:977
doublereal dL
Definition: beamslider.h:115
void Output(OutputHandler &OH) const
Definition: beamslider.cc:148
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
virtual void PutCoef(integer iRow, const doublereal &dCoef)=0
VectorExpression< VectorCrossExpr< VectorLhsExpr, VectorRhsExpr >, 3 > Cross(const VectorExpression< VectorLhsExpr, 3 > &u, const VectorExpression< VectorRhsExpr, 3 > &v)
Definition: matvec.h:3248
doublereal sRef
Definition: beamslider.h:112
virtual const Vec3 & GetXCurr(void) const
Definition: strnode.h:310
virtual void Add(integer iRow, const Vec3 &v)
Definition: vh.cc:63
enum Type iType
Definition: beamslider.h:103
BeamSliderJoint(unsigned int uL, const DofOwner *pDO, const StructNode *pN, enum Type iT, unsigned int nB, const BeamConn *const *pB, unsigned int uIB, unsigned int uIN, doublereal dl, const Vec3 &fTmp, const Mat3x3 &RTmp, flag fOut)
Definition: beamslider.cc:74
Vec3 fTmp[Beam::NUMNODES]
Definition: beamslider.h:119
virtual void ResizeReset(integer, integer)
Definition: submat.cc:182
virtual doublereal dGetPrivData(unsigned int i) const
Definition: simentity.cc:149
virtual unsigned int iGetNumPrivData(void) const
Definition: simentity.cc:136
doublereal s
Definition: beamslider.h:113
~BeamSliderJoint(void)
Definition: beamslider.cc:135
const MatCrossCross_Manip MatCrossCross
Definition: matvec3.cc:640
const doublereal dS
Definition: beamslider.cc:71
Definition: elem.h:75
void PutRowIndex(integer iSubRow, integer iRow)
Definition: submat.h:311
void Sub(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:215
doublereal dW[2]
Definition: beamslider.h:116
virtual unsigned int iGetNumDof(void) const
Definition: beamslider.h:151
Definition: joint.h:50
SubVectorHandler & InitialAssRes(SubVectorHandler &WorkVec, const VectorHandler &XCurr)
Definition: beamslider.cc:730
virtual integer iGetFirstIndex(void) const
Definition: dofown.h:127
double doublereal
Definition: colamd.c:52
std::ostream & Output(std::ostream &out, const char *sJointName, unsigned int uLabel, const Vec3 &FLocal, const Vec3 &MLocal, const Vec3 &FGlobal, const Vec3 &MGlobal) const
Definition: joint.cc:138
long int integer
Definition: colamd.c:51
Vec3 xTmp[Beam::NUMNODES]
Definition: beamslider.h:120
unsigned int GetLabel(void) const
Definition: withlab.cc:62
unsigned int nRotConstr
Definition: beamslider.h:99
virtual ~BeamConn(void)
Definition: beamslider.cc:60
void Put(unsigned short int iRow, const doublereal &dCoef)
Definition: matvec3.h:276
VariableSubMatrixHandler & InitialAssJac(VariableSubMatrixHandler &WorkMat, const VectorHandler &XCurr)
Definition: beamslider.cc:718
virtual void Resize(integer iNewSize)=0
Mat3x3 R