MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
vh.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbmath/vh.cc,v 1.31 2017/01/12 14:43:54 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 /* solution manager */
33 
34 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
35 
36 #include <cstring> /* for memset() */
37 #include <cmath>
38 
39 #include <iostream>
40 #include <iomanip>
41 
42 #include "matvec3.h"
43 #include "submat.h"
44 
45 class SubVectorHandler;
46 
47 /* VectorHandler - begin */
48 
50 {
51  NO_OP;
52 }
53 
54 void
56 {
57  Resize(iNewRow);
58  Reset();
59 }
60 
61 /* Somma un Vec3 nella posizione desiderata */
62 void
64 {
65 #ifdef DEBUG
66  IsValid();
67  ASSERT(iRow > 0);
68  ASSERT(iGetSize() >= iRow + 2);
69 #endif /* DEBUG */
70 
71  operator()(iRow) += v(1);
72  operator()(++iRow) += v(2);
73  operator()(++iRow) += v(3);
74 }
75 
76 /* Sottrae un Vec3 nella posizione desiderata */
77 void
79 {
80 #ifdef DEBUG
81  IsValid();
82  ASSERT(iRow > 0);
83  ASSERT(iGetSize() >= iRow + 2);
84 #endif /* DEBUG */
85 
86  operator()(iRow) -= v(1);
87  operator()(++iRow) -= v(2);
88  operator()(++iRow) -= v(3);
89 }
90 
91 /* Scrive un Vec3 nella posizione desiderata */
92 void
94 {
95 #ifdef DEBUG
96  IsValid();
97  ASSERT(iRow > 0);
98  ASSERT(iGetSize() >= iRow + 2);
99 #endif /* DEBUG */
100 
101  operator()(iRow) = v(1);
102  operator()(++iRow) = v(2);
103  operator()(++iRow) = v(3);
104 }
105 
106 /* Somma e moltiplica per uno scalare */
109 {
110 #ifdef DEBUG
111  IsValid();
112  VH.IsValid();
113  ASSERT(iGetSize() == VH.iGetSize());
114 #endif /* DEBUG */
115 
116  for (integer i = iGetSize(); i > 0; i--) {
117  operator()(i) += d*VH(i);
118  }
119 
120  return *this;
121 }
122 
123 /* Somma e moltiplica per uno scalare this = VH + d * VH1 */
126  const doublereal& d)
127 {
128 #ifdef DEBUG
129  IsValid();
130  VH.IsValid();
131  ASSERT(iGetSize() == VH.iGetSize());
132  VH1.IsValid();
133  ASSERT(iGetSize() == VH1.iGetSize());
134 #endif /* DEBUG */
135 
136  for (integer i = iGetSize(); i > 0; i--) {
137  operator()(i) = VH(i) + d*VH1(i);
138  }
139 
140  return *this;
141 }
142 
143 /* Moltiplica per uno scalare */
146 {
147 #ifdef DEBUG
148  IsValid();
149  VH.IsValid();
150  ASSERT(iGetSize() == VH.iGetSize());
151 #endif /* DEBUG */
152 
153  for (integer i = iGetSize(); i > 0; i--) {
154  operator()(i) = d*VH(i);
155  }
156 
157  return *this;
158 }
159 
160 
161 /* Overload di += usato per la correzione della soluzione */
164 {
165 #ifdef DEBUG
166  IsValid();
167  VH.IsValid();
168  ASSERT(iGetSize() == VH.iGetSize());
169 #endif /* DEBUG */
170 
171  for (integer i = iGetSize(); i > 0; i--) {
172  operator()(i) += VH(i);
173  }
174 
175  return *this;
176 }
177 
178 /* Overload di += usato per l'assemblaggio del residuo */
181 {
182 #ifdef DEBUG
183  SubVH.IsValid();
184 #endif /* DEBUG */
185 
186  if (SubVH.iGetSize() == 0) {
187  return *this;
188  }
189 
190  return SubVH.AddTo(*this);
191 }
192 
193 /* Overload di -= */
196 {
197 #ifdef DEBUG
198  IsValid();
199  VH.IsValid();
200  ASSERT(iGetSize() == VH.iGetSize());
201 #endif /* DEBUG */
202 
203  for (integer i = iGetSize(); i > 0; i--) {
204  operator()(i) -= VH(i);
205  }
206 
207  return *this;
208 }
209 
210 /* Overload di *= */
213 {
214 #ifdef DEBUG
215  IsValid();
216 #endif /* DEBUG */
217 
218  for (integer i = iGetSize(); i > 0; i--) {
219  operator()(i) *= d;
220  }
221 
222  return *this;
223 }
224 
225 /* Assegnazione che copia il contenuto della memoria di due handlers */
228 {
229 #ifdef DEBUG
230  VH.IsValid();
231 #endif /* DEBUG */
232 
233  integer nr = VH.iGetSize();
234  Resize(nr);
235  for (integer i = nr; i > 0; i--) {
236  operator()(i) = VH(i);
237  }
238 
239  return *this;
240 }
241 
242 /* Norma 2 del vettore */
245 {
246 #ifdef DEBUG
247  IsValid();
248 #endif /* DEBUG */
249 
250  doublereal d2 = 0.;
251 
252  for (integer i = iGetSize(); i > 0; i--) {
253  doublereal d = operator()(i);
254  d2 += d*d;
255  }
256 
257  return d2;
258 }
259 
260 /* Norma del vettore */
263 {
264  return std::sqrt(Dot());
265 }
266 
267 /* Prodotto Scalare fra due Vettori */
270 {
271 #ifdef DEBUG
272  IsValid();
273  VH.IsValid();
274  ASSERT(iGetSize() == VH.iGetSize());
275 #endif /* DEBUG */
276 
277  doublereal d2 = 0.;
278 
279  for (integer i = iGetSize(); i > 0; i--) {
280  d2 += operator()(i)*VH(i);
281  }
282 
283  return d2;
284 }
285 
286 std::ostream&
287 operator << (std::ostream& out, const VectorHandler& VH)
288 {
289  for (integer i = 1; i <= VH.iGetSize(); i++) {
290  out << std::setw(16) << VH(i) << std::endl;
291  }
292  return out;
293 }
294 
295 /* VectorHandler - end */
296 
297 
298 /* MyVectorHandler - begin */
299 
301 : bOwnsMemory(false),
302 iMaxSize(iSize), iCurSize(iSize), pdVecm1(0)
303 {
304  if (iSize == 0) {
305  ASSERT(pdVecm1 == NULL);
306 
307  } else {
308  if (pdTmpVec == NULL) {
309  bOwnsMemory = true;
310  Resize(iSize);
311  Reset();
312  } else {
313  pdVecm1 = pdTmpVec - 1;
314  }
315 #ifdef DEBUG
316  IsValid();
317 #endif /* DEBUG */
318  }
319 }
320 
322 : bOwnsMemory(false),
323 iMaxSize(VH.iCurSize), iCurSize(VH.iCurSize), pdVecm1(0)
324 {
325  if (iCurSize == 0) {
326  ASSERT(VH.pdVecm1 == 0);
327 
328  } else {
329  bOwnsMemory = true;
330  Resize(iCurSize);
331 #ifdef DEBUG
332  IsValid();
333 #endif /* DEBUG */
334 
335  for (integer i = 1; i <= iCurSize; i++) {
336  pdVecm1[i] = VH.pdVecm1[i];
337  }
338  }
339 }
340 
342 {
343  Detach();
344 }
345 
346 void
348 {
349  if (iSize < 0) {
350  silent_cerr("Negative size!" << std::endl);
352  }
353 
354  if (pdVecm1 == NULL || bOwnsMemory) {
355  if (pdVecm1 != NULL) {
356  if (iSize > iMaxSize) {
357  doublereal* pd = NULL;
358 
359  SAFENEWARR(pd, doublereal, iSize);
360  pd--;
361 #ifdef HAVE_MEMMOVE
362  memmove(pd + 1, pdVecm1 + 1, iCurSize*sizeof(doublereal));
363 #else /* ! HAVE_MEMMOVE */
364  for (integer i = iCurSize; i > 0; i--) {
365  pd[i] = pdVecm1[i];
366  }
367 #endif /* ! HAVE_MEMMOVE */
368  doublereal *pdv = pdVecm1 + 1;
369  SAFEDELETEARR(pdv);
370  pdVecm1 = pd;
371  iMaxSize = iCurSize = iSize;
372 
373  } else {
374  iCurSize = iSize;
375  }
376 
377  } else {
378  SAFENEWARR(pdVecm1, doublereal, iSize);
379  pdVecm1--;
380  iMaxSize = iCurSize = iSize;
381  bOwnsMemory = true;
382  }
383 
384  } else {
385  if (pdVecm1 != NULL) {
386  if (iSize > iMaxSize) {
387  silent_cerr("Can't resize to " << iSize
388  << ": larger than "
389  "max size " << iMaxSize << std::endl);
391 
392  } else {
393  iCurSize = iSize;
394  }
395  } else {
396  silent_cerr("internal error!" << std::endl);
398  }
399  }
400 }
401 
402 void
404 {
405  if (bOwnsMemory) {
406  if (pdVecm1 != NULL) {
407  doublereal *pd = pdVecm1 + 1;
408  SAFEDELETEARR(pd);
409  }
410  bOwnsMemory = false;
411  }
412 
413  iMaxSize = iCurSize = 0;
414  pdVecm1 = NULL;
415 }
416 
417 void
419 {
420  if (bOwnsMemory || pdVecm1 != NULL) {
421  Detach();
422  bOwnsMemory = false;
423  }
424 
425  iMaxSize = iCurSize = iSize;
426  if (iMSize > 0) {
427  if (iMSize >= iSize) {
428  iMaxSize = iMSize;
429 
430  } else if (iMSize < iSize) {
432  }
433  }
434 
435  pdVecm1 = pd - 1;
436 }
437 
438 #ifdef DEBUG
439 void
440 MyVectorHandler::IsValid(void) const
441 {
442  ASSERT(iMaxSize > 0);
443  ASSERT(iCurSize >= 0 && iCurSize <= iMaxSize);
444  ASSERT(pdVecm1 != NULL);
445 
446 #ifdef DEBUG_MEMMANAGER
447  if (bOwnsMemory) {
448  ASSERT(defaultMemoryManager.fIsBlock(pdVecm1 + 1,
449  iMaxSize*sizeof(doublereal)));
450  } else {
451  ASSERT(defaultMemoryManager.fIsValid(pdVecm1 + 1,
452  iMaxSize*sizeof(doublereal)));
453  }
454 #endif /* DEBUG_MEMMANAGER */
455 }
456 #endif /* DEBUG */
457 
458 void
460 {
461 #ifdef DEBUG
462  IsValid();
463 #endif /* DEBUG */
464 
465  ASSERT(iCurSize > 0);
466 
467 #if defined HAVE_MEMSET
468  memset(pdVecm1 + 1, 0, iGetSize()*sizeof(doublereal));
469 #else /* !HAVE_MEMSET */
470  for (integer i = iGetSize(); i > 0; i--) {
471  pdVecm1[i] = 0.;
472  }
473 #endif /* HAVE_MEMSET */
474 }
475 
476 /* Somma e moltiplica per uno scalare */
479 {
480 #ifdef DEBUG
481  IsValid();
482  VH.IsValid();
483  ASSERT(iCurSize > 0);
484  ASSERT(VH.iGetSize() > 0);
485  ASSERT(iCurSize == VH.iGetSize());
486 #endif /* DEBUG */
487 
488  if (d != 0.) {
489  for (integer i = iGetSize(); i > 0; i--) {
490  pdVecm1[i] += d*VH(i);
491  }
492  }
493 
494  return *this;
495 }
496 
497 /* Somma e moltiplica per uno scalare this = VH + d * VH1 */
500  const VectorHandler& VH1, const doublereal& d)
501 {
502 #ifdef DEBUG
503  IsValid();
504  VH.IsValid();
505  ASSERT(iGetSize() == VH.iGetSize());
506  VH1.IsValid();
507  ASSERT(iGetSize() == VH1.iGetSize());
508 #endif /* DEBUG */
509 
510  for (integer i = iGetSize(); i > 0; i--) {
511  pdVecm1[i] = VH(i) + d*VH1(i);
512  }
513 
514  return *this;
515 }
516 
517 /* Moltiplica per uno scalare */
520 {
521 #ifdef DEBUG
522  IsValid();
523  VH.IsValid();
524  ASSERT(iCurSize > 0);
525  ASSERT(VH.iGetSize() > 0);
526  ASSERT(iCurSize == VH.iGetSize());
527 #endif /* DEBUG */
528 
529  if (d == 0.) {
530  Reset();
531  } else {
532  for (integer i = iGetSize(); i > 0; i--) {
533  pdVecm1[i] = d*VH(i);
534  }
535  }
536 
537  return *this;
538 }
539 
540 /* Overload di += usato per la correzione della soluzione */
543 {
544 #ifdef DEBUG
545  IsValid();
546  VH.IsValid();
547 #endif /* DEBUG */
548 
549  ASSERT(VH.iGetSize() > 0);
550  ASSERT(iCurSize == VH.iGetSize());
551 
552  for (integer i = iGetSize(); i > 0; i--) {
553  pdVecm1[i] += VH(i);
554  }
555 
556  return *this;
557 }
558 
559 /* Overload di += usato per la correzione della soluzione */
562 {
563 #ifdef DEBUG
564  IsValid();
565  VH.IsValid();
566 #endif /* DEBUG */
567 
568  ASSERT(VH.iGetSize() > 0);
569  ASSERT(iCurSize == VH.iGetSize());
570 
571  doublereal* pdFrom = VH.pdGetVec() - 1;
572  for (integer i = iGetSize(); i > 0; i--) {
573  pdVecm1[i] += pdFrom[i];
574  }
575 
576  return *this;
577 }
578 
579 /* Overload di -= usato per la correzione della soluzione */
582 {
583 #ifdef DEBUG
584  IsValid();
585  VH.IsValid();
586 #endif /* DEBUG */
587 
588  ASSERT(VH.iGetSize() > 0);
589  ASSERT(iCurSize == VH.iGetSize());
590 
591  for (integer i = iGetSize(); i > 0; i--) {
592  pdVecm1[i] -= VH(i);
593  }
594 
595  return *this;
596 }
597 
598 /* Overload di *= */
601 {
602 #ifdef DEBUG
603  IsValid();
604 #endif /* DEBUG */
605 
606  for (integer i = iGetSize(); i > 0; i--) {
607  pdVecm1[i] *= d;
608  }
609 
610  return *this;
611 }
612 
613 /* Overload di -= usato per la correzione della soluzione */
616 {
617 #ifdef DEBUG
618  IsValid();
619  VH.IsValid();
620 #endif /* DEBUG */
621 
622  ASSERT(VH.iGetSize() > 0);
623  ASSERT(iCurSize == VH.iGetSize());
624 
625  doublereal* pdFrom = VH.pdGetVec() - 1;
626  for (integer i = iGetSize(); i > 0; i--) {
627  pdVecm1[i] -= pdFrom[i];
628  }
629 
630  return *this;
631 }
632 
633 /* Overload di = con copia termine per termine */
636 {
637 #ifdef DEBUG
638  IsValid();
639  VH.IsValid();
640 #endif /* DEBUG */
641 
642  ASSERT(VH.iGetSize() > 0);
643  ASSERT(iCurSize == VH.iGetSize());
644 
645  for (integer i = iGetSize(); i > 0; i--) {
646  pdVecm1[i] = VH(i);
647  }
648 
649  return *this;
650 }
651 
652 /* Overload di = con copia termine per termine */
655 {
656 #ifdef DEBUG
657  IsValid();
658  VH.IsValid();
659 #endif /* DEBUG */
660 
661  ASSERT(VH.iGetSize() > 0);
662  ASSERT(iCurSize == VH.iGetSize());
663 
664  doublereal* pdFrom = VH.pdGetVec() - 1;
665  for (integer i = iGetSize(); i > 0; i--) {
666  pdVecm1[i] = pdFrom[i];
667  }
668 
669  return *this;
670 }
671 
672 void
674 {
675 #ifdef DEBUG
676  IsValid();
677  ASSERT((iRow > 0) && (iRow <= iCurSize-2));
678 #endif /* DEBUG */
679 
680  pdVecm1[iRow] += v.dGet(1);
681  pdVecm1[++iRow] += v.dGet(2);
682  pdVecm1[++iRow] += v.dGet(3);
683 }
684 
685 void
687 {
688 #ifdef DEBUG
689  IsValid();
690  ASSERT((iRow > 0) && (iRow <= iCurSize-2));
691 #endif /* DEBUG */
692 
693  pdVecm1[iRow] -= v.dGet(1);
694  pdVecm1[++iRow] -= v.dGet(2);
695  pdVecm1[++iRow] -= v.dGet(3);
696 }
697 
698 void
700 {
701 #ifdef DEBUG
702  IsValid();
703  ASSERT((iRow > 0) && (iRow <= iCurSize-2));
704 #endif /* DEBUG */
705 
706  pdVecm1[iRow] = v.dGet(1);
707  pdVecm1[++iRow] = v.dGet(2);
708  pdVecm1[++iRow] = v.dGet(3);
709 }
710 
711 /* Norma 2 del vettore */
714 {
715 #ifdef DEBUG
716  IsValid();
717 #endif /* DEBUG */
718 
719  doublereal d2 = 0.;
720 
721  for (integer i = iCurSize; i > 0; i--) {
722  d2 += pdVecm1[i]*pdVecm1[i];
723  }
724 
725  return d2;
726 }
727 
728 /* MyVectorHandler - end */
virtual void Reset(void)=0
integer iCurSize
Definition: vh.h:156
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Definition: matvec3.h:98
virtual void ResizeReset(integer)
Definition: vh.cc:55
virtual VectorHandler & operator*=(const doublereal &d)
Definition: vh.cc:212
virtual const doublereal & operator()(integer iRow) const =0
virtual VectorHandler & operator=(const VectorHandler &VH)
Definition: vh.cc:635
virtual void Sub(integer iRow, const Vec3 &v)
Definition: vh.cc:78
#define SAFEDELETEARR(pnt)
Definition: mynewmem.h:713
virtual doublereal InnerProd(const VectorHandler &VH) const
Definition: vh.cc:269
void Detach(void)
Definition: vh.cc:403
virtual VectorHandler & ScalarAddMul(const VectorHandler &VH, const VectorHandler &VH1, const doublereal &d)
Definition: vh.cc:499
virtual void Put(integer iRow, const Vec3 &v)
Definition: vh.cc:699
#define NO_OP
Definition: myassert.h:74
virtual doublereal Dot(void) const
Definition: vh.cc:244
virtual integer iGetSize(void) const =0
virtual doublereal Norm(void) const
Definition: vh.cc:262
doublereal * pdVecm1
Definition: vh.h:158
virtual VectorHandler & operator-=(const VectorHandler &VH)
Definition: vh.cc:581
virtual ~VectorHandler(void)
Definition: vh.cc:49
virtual VectorHandler & operator+=(const VectorHandler &VH)
Definition: vh.cc:163
virtual VectorHandler & ScalarMul(const VectorHandler &VH, const doublereal &d)
Definition: vh.cc:519
virtual ~MyVectorHandler(void)
Definition: vh.cc:341
virtual void Resize(integer iSize)
Definition: vh.cc:347
integer iMaxSize
Definition: vh.h:155
virtual integer iGetSize(void) const
Definition: vh.h:255
doublereal Dot(void) const
Definition: vh.cc:713
virtual VectorHandler & operator-=(const VectorHandler &VH)
Definition: vh.cc:195
virtual void Add(integer iRow, const Vec3 &v)
Definition: vh.cc:673
Definition: mbdyn.h:76
virtual VectorHandler & operator=(const VectorHandler &VH)
Definition: vh.cc:227
const doublereal & dGet(unsigned short int iRow) const
Definition: matvec3.h:285
virtual VectorHandler & ScalarMul(const VectorHandler &VH, const doublereal &d)
Definition: vh.cc:145
MyVectorHandler(const MyVectorHandler &)
Definition: vh.cc:321
std::ostream & operator<<(std::ostream &out, const VectorHandler &VH)
Definition: vh.cc:287
virtual VectorHandler & AddTo(VectorHandler &VH) const =0
#define ASSERT(expression)
Definition: colamd.c:977
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
virtual void Reset(void)
Definition: vh.cc:459
virtual void Add(integer iRow, const Vec3 &v)
Definition: vh.cc:63
#define defaultMemoryManager
Definition: mynewmem.h:691
virtual VectorHandler & ScalarAddMul(const VectorHandler &VH, const doublereal &d)
Definition: vh.cc:108
virtual void Put(integer iRow, const Vec3 &v)
Definition: vh.cc:93
virtual VectorHandler & operator*=(const doublereal &d)
Definition: vh.cc:600
#define SAFENEWARR(pnt, item, sz)
Definition: mynewmem.h:701
virtual VectorHandler & operator+=(const VectorHandler &VH)
Definition: vh.cc:542
void Attach(integer iSize, doublereal *pd, integer iMSize=0)
Definition: vh.cc:418
double doublereal
Definition: colamd.c:52
virtual doublereal * pdGetVec(void) const
Definition: vh.h:245
long int integer
Definition: colamd.c:51
virtual void Resize(integer iNewSize)=0
bool bOwnsMemory
Definition: vh.h:152
virtual void Sub(integer iRow, const Vec3 &v)
Definition: vh.cc:686