MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
matvecass.h
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbmath/matvecass.h,v 1.6 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 /*
33  AUTHOR: Reinhard Resch <r.resch@secop.com>
34  Copyright (C) 2013(-2017) all rights reserved.
35 
36  The copyright of this code is transferred
37  to Pierangelo Masarati and Paolo Mantegazza
38  for use in the software MBDyn as described
39  in the GNU Public License version 2.1
40 */
41 
42 #ifndef ___MATVECASS_H__INCLUDED___
43 #define ___MATVECASS_H__INCLUDED___
44 
45 #if GRADIENT_DEBUG > 0 && HAVE_FINITE == 1
46 #include <cmath>
47 #endif
48 #include "gradient.h"
49 #include "matvec.h"
50 #include "submat.h"
51 
52 namespace grad {
53 
54 template <typename T>
56 
57 template <>
59 public:
61  :vh(vh) {
62 
63  }
64 
65  void dGetCoef(integer iRow, doublereal& dVal, doublereal dCoef, LocalDofMap* pDofMap) const {
66  dVal = vh.dGetCoef(iRow);
67  }
68 
69  template <index_type N_rows>
70  void GetVec(integer iRow, Vector<double, N_rows>& v, doublereal dCoef, LocalDofMap* pDofMap) const {
71  for (integer i = 0; i < v.iGetNumRows(); ++i) {
72  v(i + 1) = vh.dGetCoef(iRow + i);
73  }
74  }
75 
76 private:
77  const VectorHandler& vh;
78 };
79 
80 template <index_type N_SIZE>
82 public:
84  :vh(vh) {
85 
86  }
87 
88  void dGetCoef(integer iRow, Gradient<N_SIZE>& gVal, doublereal dCoef, LocalDofMap* pDofMap) const {
89  gVal.SetValuePreserve(vh.dGetCoef(iRow));
90  gVal.DerivativeResizeReset(pDofMap, iRow, MapVectorBase::GLOBAL, -dCoef);
91  }
92 
93  template <index_type N_rows>
94  void GetVec(integer iRow, Vector<Gradient<N_SIZE>, N_rows>& v, doublereal dCoef, LocalDofMap* pDofMap) const {
95  for (integer i = 0; i < v.iGetNumRows(); ++i) {
96  Gradient<N_SIZE>& v_i = v(i + 1);
97  v_i.SetValuePreserve(vh.dGetCoef(iRow + i));
98  v_i.DerivativeResizeReset(pDofMap, iRow, iRow + v.iGetNumRows(), MapVectorBase::GLOBAL, 0.);
99  v_i.SetDerivativeGlobal(iRow + i, -dCoef);
100  }
101  }
102 private:
104 };
105 
107 public:
108  enum mode_t { RESET, APPEND };
109 };
110 
111 template <typename T>
113 
114 template <>
116 public:
117  GradientAssVec(SubVectorHandler& vh, enum mode_t mode = RESET)
118  :WorkVec(vh) {
119 
120  switch (mode) {
121  case APPEND:
122  iSubRow = WorkVec.iGetSize();
123  break;
124 
125  case RESET:
126  iSubRow = 0;
127  WorkVec.Resize(iSubRow);
128  break;
129 
130  default:
131  GRADIENT_ASSERT(0);
132  }
133  }
134 
135  template <typename T>
136  static void
137  AssRes(T* pElem,
138  SubVectorHandler& WorkVec,
139  doublereal dCoef,
140  const VectorHandler& XCurr,
141  const VectorHandler& XPrimeCurr,
142  enum FunctionCall func,
143  enum mode_t mode = RESET) {
144  const GradientVectorHandler<doublereal> XCurr_grad(XCurr);
145  const GradientVectorHandler<doublereal> XPrimeCurr_grad(XPrimeCurr);
146  GradientAssVec WorkVec_grad(WorkVec, mode);
147 
148  pElem->AssRes(WorkVec_grad, dCoef, XCurr_grad, XPrimeCurr_grad, func);
149  }
150 
151  template <typename T>
152  static void
153  InitialAssRes(T* pElem,
154  SubVectorHandler& WorkVec,
155  const VectorHandler& XCurr,
156  enum FunctionCall func,
157  enum mode_t mode = RESET) {
158  const GradientVectorHandler<doublereal> XCurr_grad(XCurr);
159  GradientAssVec WorkVec_grad(WorkVec, mode);
160 
161  pElem->InitialAssRes(WorkVec_grad, XCurr_grad, func);
162  }
163 
164  void AddItem(const integer iRow, const double dCoef) {
165  WorkVec.Resize(++iSubRow);
166 
167  GRADIENT_ASSERT(iSubRow <= WorkVec.iGetSize());
168 #if GRADIENT_DEBUG > 0 && HAVE_FINITE == 1
169  if (!std::isfinite(dCoef)) {
170  silent_cerr("PutItem(" << iSubRow << ", " << iRow << ", " << dCoef << ")" << std::endl);
171  }
172  GRADIENT_ASSERT(std::isfinite(dCoef));
173 #endif
174  WorkVec.PutItem(iSubRow, iRow, dCoef);
175  }
176 
177  template <index_type N_rows>
178  void AddItem(const integer iFirstRow, const Vector<doublereal, N_rows>& v) {
179  // zero based index according to VectorHandler::Put(integer iRow, const Vec3& v)
180  WorkVec.Resize(iSubRow + v.iGetNumRows());
181 
182  for (integer i = 0; i < v.iGetNumRows(); ++i) {
183  GRADIENT_ASSERT(iSubRow + 1 <= WorkVec.iGetSize());
184 #if GRADIENT_DEBUG > 0 && HAVE_FINITE == 1
185  if (!std::isfinite(v(i + 1))) {
186  silent_cerr("PutItem(" << iSubRow + 1 << ", " << iFirstRow + i << ", " << v(i + 1) << ")" << std::endl);
187  }
188  GRADIENT_ASSERT(std::isfinite(v(i + 1)));
189 #endif
190  WorkVec.PutItem(++iSubRow, iFirstRow + i, v(i + 1));
191  }
192  }
193 
194 private:
197 };
198 
199 template <index_type N_SIZE>
200 class GradientAssVec<Gradient<N_SIZE> >: public GradientAssVecBase {
201 public:
203  :WorkMat(mh) {
204 
205  switch (mode) {
206  case RESET:
207  iNextItem = 1;
208  WorkMat.ResizeReset(iNextItem, 0);
209  WorkMat.PutItem(1, 1, 1, 0.); //FIXME: avoid SIGSEGV if the matrix is empty
210  break;
211 
212  case APPEND:
213  iNextItem = WorkMat.iGetNumRows() + 1;
214  break;
215 
216  default:
217  GRADIENT_ASSERT(0);
218  }
219 
220  }
221 
222  template <typename T>
223  static void AssJac(T* pElem,
224  SparseSubMatrixHandler& WorkMat,
225  doublereal dCoef,
226  const VectorHandler& XCurr,
227  const VectorHandler& XPrimeCurr,
228  enum FunctionCall func,
229  LocalDofMap* pDofMap,
230  enum mode_t mode = RESET) {
231 
232  const GradientVectorHandler<Gradient<N_SIZE> > XCurr_grad(XCurr);
233  const GradientVectorHandler<Gradient<N_SIZE> > XPrimeCurr_grad(XPrimeCurr);
234 
235  GradientAssVec WorkMat_grad(WorkMat, mode);
236 
237  if (pDofMap) {
238  pDofMap->Reset(func);
239  }
240 
241  pElem->AssRes(WorkMat_grad, dCoef, XCurr_grad, XPrimeCurr_grad, func);
242  }
243 
244  template <typename T>
245  static void InitialAssJac(T* pElem,
246  SparseSubMatrixHandler& WorkMat,
247  const VectorHandler& XCurr,
248  enum FunctionCall func,
249  LocalDofMap* pDofMap,
250  enum mode_t mode = RESET) {
251 
252  const GradientVectorHandler<Gradient<N_SIZE> > XCurr_grad(XCurr);
253 
254  GradientAssVec WorkMat_grad(WorkMat, mode);
255 
256  if (pDofMap) {
257  pDofMap->Reset(func);
258  }
259 
260  pElem->InitialAssRes(WorkMat_grad, XCurr_grad, func);
261  }
262 
263  void AddItem(integer iRow, const Gradient<N_SIZE>& g) {
264  for (index_type i = g.iGetStartIndexLocal(); i < g.iGetEndIndexLocal(); ++i) {
265  const double dCoef = g.dGetDerivativeLocal(i);
266  if (dCoef) {
267  const index_type iDofIndex = g.iGetGlobalDof(i);
268  WorkMat.Resize(iNextItem, 0);
269  GRADIENT_ASSERT(iNextItem <= WorkMat.iGetNumRows());
270 #if GRADIENT_DEBUG > 0 && HAVE_FINITE == 1
271  if (!std::isfinite(dCoef)) {
272  silent_cerr("PutItem(" << iNextItem << ", " << iRow << ", " << iDofIndex << ", " << dCoef << ")" << std::endl);
273  }
274 
275  GRADIENT_ASSERT(std::isfinite(dCoef));
276 #endif
277  WorkMat.PutItem(iNextItem++, iRow, iDofIndex, dCoef);
278  }
279  }
280  }
281 
282  template <index_type N_rows>
283  void AddItem(integer iFirstRow, const Vector<Gradient<N_SIZE>, N_rows>& v) {
284  // zero based index according to VectorHandler::Put(integer iRow, const Vec3& v)
285 
286  for (integer i = 0; i < N_rows; ++i) {
287  AddItem(iFirstRow + i, v(i + 1));
288  }
289  }
290 private:
293 };
294 
295 } // namespace
296 
297 #endif /* ___MATVECASS_H__INCLUDED___ */
index_type iGetStartIndexLocal() const
Definition: gradient.h:2576
void GetVec(integer iRow, Vector< Gradient< N_SIZE >, N_rows > &v, doublereal dCoef, LocalDofMap *pDofMap) const
Definition: matvecass.h:94
GradientVectorHandler(const VectorHandler &vh)
Definition: matvecass.h:60
GradientAssVec(SubVectorHandler &vh, enum mode_t mode=RESET)
Definition: matvecass.h:117
void AddItem(integer iRow, const Gradient< N_SIZE > &g)
Definition: matvecass.h:263
void AddItem(const integer iRow, const double dCoef)
Definition: matvecass.h:164
integer index_type
Definition: gradient.h:104
static void InitialAssRes(T *pElem, SubVectorHandler &WorkVec, const VectorHandler &XCurr, enum FunctionCall func, enum mode_t mode=RESET)
Definition: matvecass.h:153
void GetVec(integer iRow, Vector< double, N_rows > &v, doublereal dCoef, LocalDofMap *pDofMap) const
Definition: matvecass.h:70
static void InitialAssJac(T *pElem, SparseSubMatrixHandler &WorkMat, const VectorHandler &XCurr, enum FunctionCall func, LocalDofMap *pDofMap, enum mode_t mode=RESET)
Definition: matvecass.h:245
void dGetCoef(integer iRow, doublereal &dVal, doublereal dCoef, LocalDofMap *pDofMap) const
Definition: matvecass.h:65
void func(const T &u, const T &v, const T &w, doublereal e, T &f)
SparseSubMatrixHandler & WorkMat
Definition: matvecass.h:291
void SetDerivativeGlobal(index_type iGlobalDof, scalar_deriv_type dCoef)
Definition: gradient.h:2566
static void AssRes(T *pElem, SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr, enum FunctionCall func, enum mode_t mode=RESET)
Definition: matvecass.h:137
FunctionCall
Definition: gradient.h:1018
GradientVectorHandler(const VectorHandler &vh)
Definition: matvecass.h:83
void SetValuePreserve(scalar_func_type dVal)
Definition: gradient.h:2511
static void AssJac(T *pElem, SparseSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr, enum FunctionCall func, LocalDofMap *pDofMap, enum mode_t mode=RESET)
Definition: matvecass.h:223
void Reset(enum FunctionCall func=UNKNOWN_FUNC)
Definition: gradient.h:1110
GradientAssVec(SparseSubMatrixHandler &mh, enum mode_t mode=RESET)
Definition: matvecass.h:202
void DerivativeResizeReset(LocalDofMap *pMap, index_type iStartGlobal, index_type iEndGlobal, MapVectorBase::GlobalScope s, scalar_deriv_type dVal)
Definition: gradient.h:2538
index_type iGetGlobalDof(index_type iLocalDof) const
Definition: gradient.h:2572
index_type iGetNumRows() const
Definition: matvec.h:2480
void AddItem(integer iFirstRow, const Vector< Gradient< N_SIZE >, N_rows > &v)
Definition: matvecass.h:283
void dGetCoef(integer iRow, Gradient< N_SIZE > &gVal, doublereal dCoef, LocalDofMap *pDofMap) const
Definition: matvecass.h:88
scalar_deriv_type dGetDerivativeLocal(index_type iLocalDof) const
Definition: gradient.h:2516
void AddItem(const integer iFirstRow, const Vector< doublereal, N_rows > &v)
Definition: matvecass.h:178
double doublereal
Definition: colamd.c:52
long int integer
Definition: colamd.c:51
index_type iGetEndIndexLocal() const
Definition: gradient.h:2580
#define GRADIENT_ASSERT(expr)
Definition: gradient.h:74