MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
module-octave.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/modules/module-octave/module-octave.cc,v 1.21 2017/01/12 14:56:18 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) 2011(-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 #include <mbconfig.h> /* This goes first in every *.c,*.cc file */
43 
44 #ifdef USE_OCTAVE
45 
46 // FIXME: there is a conflict between the MBDyn and octave real type
47 #define real mbdyn_real_type
48 #include <mbdefs.h>
49 #include <matvec3.h>
50 #include <matvec6.h>
51 #include <dataman.h>
52 #include <drive.h>
53 #include <tpldrive.h>
54 #include <userelem.h>
55 #include <constltp.h>
56 #undef real
57 
58 #include <cmath>
59 #include <cfloat>
60 #include <string>
61 #include <cassert>
62 #include <limits>
63 #include <map>
64 #include <sstream>
65 #include <stdint.h>
66 
67 // note: undef MPI_Comm and MPI_Info because of a conflict
68 // between some mpi implementations and hdf5 (used by NetCDF)
69 #undef MPI_Comm
70 #undef MPI_Info
71 #include <octave/oct.h>
72 #include <octave/parse.h>
73 #include <octave/toplev.h>
74 #include <octave/octave.h>
75 
76 #include "module-octave.h"
77 #include "octave_object.h"
78 
79 //#define DEBUG
80 
81 #ifdef DEBUG
82 #define TRACE(msg) ((void)(std::cerr << __FILE__ << ":" << __LINE__ << ":" << __PRETTY_FUNCTION__ << ":" << msg << std::endl))
83 #define ASSERT(expr) assert(expr)
84 #else
85 #define TRACE(msg) ((void)0)
86 #endif
87 
88 namespace oct {
89 
90 class OctaveInterface {
91 public:
92  enum OctaveCallFlags_t {
93  DEFAULT_CALL_FLAGS = 0x0,
94  PASS_DATA_MANAGER = 0x1,
95  UPDATE_OCTAVE_VARIABLES = 0x2,
96  UPDATE_MBDYN_VARIABLES = 0x4,
97  UPDATE_GLOBAL_VARIABLES = UPDATE_OCTAVE_VARIABLES | UPDATE_MBDYN_VARIABLES,
98  OPTIONAL_OUTPUT_ARGS = 0x8
99  };
100 
101 private:
102  explicit OctaveInterface(const DataManager* pDM, MBDynParser* pHP);
103  virtual ~OctaveInterface(void);
104 
105 public:
106  static OctaveInterface* CreateInterface(const DataManager* pDM, MBDynParser* pHP);
107  void Destroy();
108  static OctaveInterface* GetInterface(void) { return pOctaveInterface; };
109  void UpdateOctaveVariables(void);
110  void UpdateMBDynVariables(void);
111  const DataManager* GetDataManager(void) const{ return pDM; }
112  MBDynParser* GetMBDynParser(void) const { return pHP; }
113  const octave_value& GetDataManagerInterface(void) const { return octDM; }
114  const octave_value& GetMBDynParserInterface(void) const { return octHP; }
115  bool AddOctaveSearchPath(const std::string& path);
116  static bool ConvertMBDynToOctave(const TypedValue& mbValue, octave_value& octValue);
117  static bool ConvertMBDynToOctave(doublereal mbValue, octave_value& octValue); // useful mainly for templates
118  static bool ConvertMBDynToOctave(const Vec3& mbValue, octave_value& octValue);
119  static bool ConvertMBDynToOctave(const Vec6& mbValue, octave_value& octValue);
120  static bool ConvertMBDynToOctave(const Mat3x3& mbValue, octave_value& octValue);
121  static bool ConvertMBDynToOctave(const Mat6x6& mbValue, octave_value& octValue);
122  static bool ConvertOctaveToMBDyn(const ColumnVector& octValue, doublereal& mbValue);
123  static bool ConvertOctaveToMBDyn(const ColumnVector& octValue, Vec3& mbValue);
124  static bool ConvertOctaveToMBDyn(const ColumnVector& octValue, Vec6& mbValue);
125  static bool ConvertOctaveToMBDyn(const Matrix& octValue, doublereal& mbValue);
126  static bool ConvertOctaveToMBDyn(const Matrix& octValue, Vec3& mbValue);
127  static bool ConvertOctaveToMBDyn(const Matrix& octValue, Vec6& mbValue);
128  static bool ConvertOctaveToMBDyn(const Matrix& octValue, Mat3x3& mbValue);
129  static bool ConvertOctaveToMBDyn(const Matrix& octValue, Mat6x6& mbValue);
130  static bool ConvertOctaveToMBDyn(const octave_value& octValue, TypedValue& mbValue);
131  static bool ConvertOctaveToMBDyn(const octave_value& octValue, doublereal& mbValue);
132  inline octave_value_list EvalFunction(const std::string& func, const octave_value_list& args, int nargout = 1, int flags = DEFAULT_CALL_FLAGS);
133  inline octave_value_list EvalFunctionDerivative(const std::string& func, const octave_value_list& args, int flags = DEFAULT_CALL_FLAGS);
134  inline octave_value_list EvalFunctionDerivative(const std::string& func, doublereal dVar, int flags = DEFAULT_CALL_FLAGS);
135  inline octave_value_list EvalFunction(const std::string& func, doublereal dVar, int nargout = 1, int flags = DEFAULT_CALL_FLAGS);
136  inline doublereal EvalScalarFunction(const std::string& func, doublereal dVar, int flags = DEFAULT_CALL_FLAGS);
137  inline doublereal EvalScalarFunction(const std::string& func, const octave_value_list& args, int flags = DEFAULT_CALL_FLAGS);
138  inline doublereal EvalScalarFunctionDerivative(const std::string& func, doublereal dVar, int flags = DEFAULT_CALL_FLAGS);
139  inline doublereal EvalScalarFunctionDerivative(const std::string& func, const octave_value_list& args, int flags = DEFAULT_CALL_FLAGS);
140  template<typename T>
141  inline void EvalMatrixFunction(const std::string& func, T& res, doublereal dVar, int flags = DEFAULT_CALL_FLAGS);
142  template <typename T>
143  inline void EvalMatrixFunction(const std::string& func, T& res, const octave_value_list& args, int flags = DEFAULT_CALL_FLAGS);
144  template<typename T>
145  inline void EvalMatrixFunctionDerivative(const std::string& func, T& res, doublereal dVar, int flags = DEFAULT_CALL_FLAGS);
146  template<typename T>
147  inline void EvalMatrixFunctionDerivative(const std::string& func, T& res, const octave_value_list& args, int flags = DEFAULT_CALL_FLAGS);
148  static bool HaveADPackage(void) { return bHaveADPackage; };
149  void AddEmbedFileName(const std::string& strFile);
150  bool bHaveMethod(const octave_value& octObject, const std::string& strClass, const std::string& strName);
151  inline octave_value_list MakeArgList(doublereal dVar, const octave_value_list& args, int iFlags) const;
152 private:
153  template<typename T>
154  inline static bool ConvertOctaveToMBDyn(const ColumnVector& octValue, T& mbValue, int rows);
155  template<typename T>
156  inline static bool ConvertOctaveToMBDyn(const Matrix& octValue, T& mbValue, int rows);
157  template<typename T>
158  inline static bool ConvertOctaveToMBDyn(const Matrix& octValue, T& mbValue, int rows, int cols);
159  bool LoadADPackage(void);
160  template<typename T>
161  inline static bool ConvertMBDynToOctave(const T& mbValue, octave_value& octValue, int rows);
162  template<typename T>
163  inline static bool ConvertMBDynToOctave(const T& mbValue, octave_value& octValue, int rows, int cols);
164  static void exit(int);
165 
166 private:
167  typedef std::map<std::string, bool> EmbedFileNameMap_t;
168  typedef EmbedFileNameMap_t::iterator EmbedFileNameIter_t;
169  bool bFirstCall;
170  bool bEmbedFileDirty;
171  const octave_value octDM;
172  const octave_value octHP;
173  const DataManager* const pDM;
174  MBDynParser* const pHP;
175  EmbedFileNameMap_t strEmbedFileNames;
176  static OctaveInterface* pOctaveInterface;
177  static int iRefCount;
178  static const std::string strADFunc;
179  static const std::string strIsMethod;
180  static bool bHaveADPackage;
181 };
182 
183 class MBDynInterface : public octave_object {
184 public:
185  // An instance of this class might be created directly from within octave
186  // In this case pInterface must be defined
187  explicit MBDynInterface(OctaveInterface* pInterface = OctaveInterface::GetInterface());
188  virtual ~MBDynInterface(void);
189  virtual void print(std::ostream& os, bool pr_as_read_syntax = false) const;
190 protected:
191  BEGIN_METHOD_TABLE_DECLARE()
192  METHOD_DECLARE(GetVersion)
193  END_METHOD_TABLE_DECLARE()
194 
195 private:
196  DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA
197  DECLARE_OCTAVE_ALLOCATOR
198 
199 protected:
200  OctaveInterface* GetInterface(void) const { return pInterface; };
201 
202 private:
203  OctaveInterface* const pInterface;
204 };
205 
206 class ConstVectorHandlerInterface : public MBDynInterface {
207 public:
208  explicit ConstVectorHandlerInterface(OctaveInterface* pInterface = OctaveInterface::GetInterface(), const VectorHandler* pX=0);
209  virtual ~ConstVectorHandlerInterface();
210  void Set(const VectorHandler* pX){ this->pX = const_cast<VectorHandler*>(pX); }
211  virtual void print(std::ostream& os, bool pr_as_read_syntax = false) const;
212  virtual octave_value operator()(const octave_value_list& idx) const;
213  virtual dim_vector dims (void) const;
214 protected:
215  BEGIN_METHOD_TABLE_DECLARE()
216  METHOD_DECLARE(dGetCoef)
217  METHOD_DECLARE(GetVec)
218  METHOD_DECLARE(iGetSize)
219  END_METHOD_TABLE_DECLARE()
220 
221 private:
222  DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA
223  DECLARE_OCTAVE_ALLOCATOR
224 
225 protected:
226  VectorHandler* pX;
227 };
228 
229 class VectorHandlerInterface : public ConstVectorHandlerInterface {
230 public:
231  explicit VectorHandlerInterface(OctaveInterface* pInterface = OctaveInterface::GetInterface(), VectorHandler* X=0);
232  virtual ~VectorHandlerInterface();
233 protected:
234  BEGIN_METHOD_TABLE_DECLARE()
235  METHOD_DECLARE(PutCoef)
236  METHOD_DECLARE(PutVec)
237  END_METHOD_TABLE_DECLARE()
238 
239 private:
240  DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA
241  DECLARE_OCTAVE_ALLOCATOR
242 
243 };
244 
245 class OStreamInterface : public MBDynInterface {
246 public:
247  explicit OStreamInterface(OctaveInterface* pInterface = OctaveInterface::GetInterface(), std::ostream* pOS = 0);
248  virtual ~OStreamInterface();
249  void Set(std::ostream* pOS){ this->pOS = pOS; }
250  std::ostream* Get()const{ return pOS; }
251 protected:
252  BEGIN_METHOD_TABLE_DECLARE()
253  METHOD_DECLARE(printf)
254  END_METHOD_TABLE_DECLARE()
255 
256 private:
257  DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA
258  DECLARE_OCTAVE_ALLOCATOR
259 
260 private:
261  std::ostream* pOS;
262  static const std::string strsprintf;
263 };
264 
265 class SimulationEntityInterface: public MBDynInterface {
266 public:
267  explicit SimulationEntityInterface(OctaveInterface* pInterface = OctaveInterface::GetInterface());
268  virtual ~SimulationEntityInterface();
269  virtual const SimulationEntity* Get()const=0;
270 protected:
271  BEGIN_METHOD_TABLE_DECLARE()
272  METHOD_DECLARE(iGetNumPrivData)
273  METHOD_DECLARE(iGetPrivDataIdx)
274  METHOD_DECLARE(dGetPrivData)
275  END_METHOD_TABLE_DECLARE()
276 };
277 
278 class NodeInterface: public SimulationEntityInterface {
279 public:
280  explicit NodeInterface(OctaveInterface* pInterface = OctaveInterface::GetInterface());
281  virtual ~NodeInterface();
282  virtual void print(std::ostream& os, bool pr_as_read_syntax = false) const;
283  virtual const Node* Get()const=0;
284 protected:
285  BEGIN_METHOD_TABLE_DECLARE()
286  METHOD_DECLARE(GetLabel)
287  METHOD_DECLARE(GetName)
288  METHOD_DECLARE(iGetFirstIndex)
289  METHOD_DECLARE(iGetFirstRowIndex)
290  METHOD_DECLARE(iGetFirstColIndex)
291  METHOD_DECLARE(dGetDofValue)
292  METHOD_DECLARE(dGetDofValuePrev)
293  END_METHOD_TABLE_DECLARE()
294 };
295 
296 class ScalarNodeInterface: public NodeInterface {
297 public:
298  explicit ScalarNodeInterface(OctaveInterface* pInterface = OctaveInterface::GetInterface(), ScalarNode* pNode = 0);
299  virtual ~ScalarNodeInterface();
300  virtual const ScalarNode* Get()const;
301 
302 protected:
303  BEGIN_METHOD_TABLE_DECLARE()
304  METHOD_DECLARE(SetX)
305  METHOD_DECLARE(dGetX)
306  METHOD_DECLARE(SetXPrime)
307  METHOD_DECLARE(dGetXPrime)
308  END_METHOD_TABLE_DECLARE()
309 
310 private:
311  DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA
312  DECLARE_OCTAVE_ALLOCATOR
313 
314 private:
315  ScalarNode* const pNode;
316 };
317 
318 class StructDispNodeBaseInterface : public NodeInterface {
319 public:
320  explicit StructDispNodeBaseInterface(OctaveInterface* pInterface = OctaveInterface::GetInterface());
321  virtual ~StructDispNodeBaseInterface();
322  virtual const StructDispNode* Get()const=0;
323 
324 protected:
325  BEGIN_METHOD_TABLE_DECLARE()
326  METHOD_DECLARE(iGetFirstPositionIndex)
327  METHOD_DECLARE(iGetFirstMomentumIndex)
328  METHOD_DECLARE(GetLabel)
329  METHOD_DECLARE(GetXCurr)
330  METHOD_DECLARE(GetXPrev)
331  METHOD_DECLARE(GetVCurr)
332  METHOD_DECLARE(GetVPrev)
333  METHOD_DECLARE(GetXPPCurr)
334  METHOD_DECLARE(GetXPPPrev)
335  END_METHOD_TABLE_DECLARE()
336 
337 protected:
338  octave_value GetVec3(const Vec3& v, const octave_value_list& args) const;
339  octave_value GetMat3x3(const Mat3x3& m, const octave_value_list& args) const;
340 };
341 
342 class StructDispNodeInterface: public StructDispNodeBaseInterface {
343 public:
344  explicit StructDispNodeInterface(OctaveInterface* pInterface = OctaveInterface::GetInterface(), const StructDispNode* pNode = 0);
345  virtual ~StructDispNodeInterface();
346  virtual const StructDispNode* Get()const;
347 
348 private:
349  DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA
350  DECLARE_OCTAVE_ALLOCATOR
351 
352 private:
353  const StructDispNode* const pNode;
354 };
355 
356 class StructNodeInterface : public StructDispNodeBaseInterface {
357 public:
358  explicit StructNodeInterface(OctaveInterface* pInterface = OctaveInterface::GetInterface(), const StructNode* pNode = 0);
359  virtual ~StructNodeInterface();
360  virtual const StructNode* Get()const;
361 protected:
362  BEGIN_METHOD_TABLE_DECLARE()
363  METHOD_DECLARE(GetgCurr)
364  METHOD_DECLARE(GetgRef)
365  METHOD_DECLARE(GetgPCurr)
366  METHOD_DECLARE(GetgPRef)
367  METHOD_DECLARE(GetRCurr)
368  METHOD_DECLARE(GetRPrev)
369  METHOD_DECLARE(GetRRef)
370  METHOD_DECLARE(GetWCurr)
371  METHOD_DECLARE(GetWPrev)
372  METHOD_DECLARE(GetWRef)
373  METHOD_DECLARE(GetWPCurr)
374  METHOD_DECLARE(GetWPPrev)
375  END_METHOD_TABLE_DECLARE()
376 
377 private:
378  DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA
379  DECLARE_OCTAVE_ALLOCATOR
380 
381 private:
382  const StructNode* const pNode;
383 };
384 
385 class DataManagerInterface : public MBDynInterface {
386 public:
387  // An instance of this class might be created directly from within octave
388  // In this case pInterface must be defined
389  explicit DataManagerInterface(OctaveInterface* pInterface = OctaveInterface::GetInterface(), const DataManager* pDM = 0);
390  virtual ~DataManagerInterface();
391  virtual void print(std::ostream& os, bool pr_as_read_syntax = false) const;
392  inline const DataManager* GetDataManager(void)const;
393  inline const Table& GetSymbolTable(void)const;
394 
395 protected:
396  BEGIN_METHOD_TABLE_DECLARE()
397  METHOD_DECLARE(GetVariable)
398  METHOD_DECLARE(GetStructNodePos)
399  METHOD_DECLARE(GetStructNode)
400  METHOD_DECLARE(pFindNode)
401  METHOD_DECLARE(ReadNode)
402  METHOD_DECLARE(dGetTime)
403  END_METHOD_TABLE_DECLARE()
404 
405 private:
406  DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA
407  DECLARE_OCTAVE_ALLOCATOR
408 
409 private:
410  Node* pFindNode_(const octave_value_list& args, Node::Type& Type) const;
411  Node::Type GetNodeType(const std::string& strType) const;
412  NodeInterface* CreateNodeInterface(Node* pNode, Node::Type Type) const;
413 
414 private:
415  const DataManager* const pDM;
416 };
417 
418 class MBDynParserInterface : public MBDynInterface {
419 public:
420  // An instance of this class might be created directly from within octave
421  // In this case pInterface must be defined
422  explicit MBDynParserInterface(OctaveInterface* pInterface = OctaveInterface::GetInterface(), MBDynParser* pHP=0);
423  virtual ~MBDynParserInterface();
424  virtual void print(std::ostream& os, bool pr_as_read_syntax = false) const;
425  MBDynParser* Get()const{ return pHP; }
426 
427 protected:
428  BEGIN_METHOD_TABLE_DECLARE()
429  METHOD_DECLARE(IsKeyWord)
430  METHOD_DECLARE(IsArg)
431  METHOD_DECLARE(IsStringWithDelims)
432  METHOD_DECLARE(GetReal)
433  METHOD_DECLARE(GetInt)
434  METHOD_DECLARE(GetBool)
435  METHOD_DECLARE(GetYesNoOrBool)
436  METHOD_DECLARE(GetString)
437  METHOD_DECLARE(GetStringWithDelims)
438  METHOD_DECLARE(GetValue)
439  METHOD_DECLARE(GetPosRel)
440  METHOD_DECLARE(GetRotRel)
441  METHOD_DECLARE(GetLineData)
442  METHOD_DECLARE(GetDriveCaller)
443  END_METHOD_TABLE_DECLARE()
444 
445 private:
446  DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA
447  DECLARE_OCTAVE_ALLOCATOR
448 
449 private:
450  const StructDispNode* GetStructNode(const octave_value& arg)const;
451  static bool GetDelimsFromString(const std::string& strDelims, HighParser::Delims& delims);
452 private:
453  MBDynParser* pHP;
454  static const struct MBDynStringDelims {
455  HighParser::Delims value;
456  char name[15];
457  } mbStringDelims[6];
458 };
459 
460 class DriveCallerInterface : public MBDynInterface {
461 public:
462  explicit DriveCallerInterface(OctaveInterface* pInterface = OctaveInterface::GetInterface(), const DriveCaller* pDC=0);
463  virtual ~DriveCallerInterface();
464  void Set(const DriveCaller* pDC){ DC.Set(pDC); }
465 protected:
466  BEGIN_METHOD_TABLE_DECLARE()
467  METHOD_DECLARE(dGet)
468  METHOD_DECLARE(dGetP)
469  END_METHOD_TABLE_DECLARE()
470 
471 private:
472  DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA
473  DECLARE_OCTAVE_ALLOCATOR
474 
475 protected:
476  DriveOwner DC;
477 };
478 
479 class OctaveDriveCaller : public DriveCaller {
480 public:
481  explicit
482  OctaveDriveCaller(const std::string& strFunc,
483  OctaveInterface* pInterface,
484  int iFlags,
485  const octave_value_list& args);
486  virtual ~OctaveDriveCaller(void);
487 
488  /* Copia */
489  virtual DriveCaller* pCopy(void) const;
490 
491  /* Scrive il contributo del DriveCaller al file di restart */
492  virtual std::ostream& Restart(std::ostream& out) const;
493 
494  inline doublereal dGet(const doublereal& dVar) const;
495  inline doublereal dGet(void) const;
496 
497  /* this is about drives that are differentiable */
498  virtual bool bIsDifferentiable(void) const;
499  virtual doublereal dGetP(const doublereal& dVar) const;
500  virtual inline doublereal dGetP(void) const;
501 
502 private:
503  inline octave_value_list MakeArgList(doublereal dVar) const;
504 
505 private:
506  const int iFlags;
507  const std::string strFunc;
508  OctaveInterface* const pInterface;
509  const octave_value_list args;
510 };
511 
512 template <class T>
513 class OctaveTplDriveCaller : public TplDriveCaller<T> {
514 public:
515  OctaveTplDriveCaller(const std::string& strFunction, OctaveInterface* pInterface, int iFlags, const octave_value_list& args);
516  ~OctaveTplDriveCaller(void);
517  virtual TplDriveCaller<T>* pCopy(void) const;
518  virtual std::ostream& Restart(std::ostream& out) const;
519  virtual std::ostream& Restart_int(std::ostream& out) const;
520  virtual inline T Get(const doublereal& dVar) const;
521  virtual inline T Get(void) const;
522  virtual inline bool bIsDifferentiable(void) const;
523  virtual inline T GetP(void) const;
524  virtual inline int getNDrives(void) const;
525 
526 private:
527  inline octave_value_list MakeArgList(doublereal dVar) const;
528 
529 private:
530  const std::string strFunction;
531  OctaveInterface* const pInterface;
532  const int iFlags;
533  const octave_value_list args;
534 };
535 
536 class DerivativeDriveCaller : public DriveCaller {
537 public:
538  explicit DerivativeDriveCaller(DriveCaller* pDriveCaller);
539  virtual ~DerivativeDriveCaller(void);
540 
541  /* Copia */
542  virtual DriveCaller* pCopy(void) const;
543 
544  /* Scrive il contributo del DriveCaller al file di restart */
545  virtual std::ostream& Restart(std::ostream& out) const;
546 
547  inline doublereal dGet(const doublereal& dVar) const;
548  inline doublereal dGet(void) const;
549 
550  /* this is about drives that are differentiable */
551  virtual bool bIsDifferentiable(void) const;
552  virtual doublereal dGetP(const doublereal& dVar) const;
553  virtual inline doublereal dGetP(void) const;
554 
555 private:
556  DriveOwner pDriveCaller;
557 };
558 
559 class OctaveScalarFunction : public DifferentiableScalarFunction {
560 public:
561  OctaveScalarFunction(const std::string& strFunc, OctaveInterface* pInterface, int iFlags, const octave_value_list& args);
562  virtual ~OctaveScalarFunction(void);
563  virtual doublereal operator()(const doublereal x) const;
564  virtual doublereal ComputeDiff(const doublereal t, const integer order = 1) const;
565 
566 private:
567  inline octave_value_list MakeArgList(doublereal dVar) const;
568 
569 private:
570  const std::string strFunc;
571  OctaveInterface* const pInterface;
572  const int iFlags;
573  const octave_value_list args;
574 };
575 
576 class OctaveConstitutiveLawBase
577 {
578 public:
579  explicit inline OctaveConstitutiveLawBase(const std::string& strClass, OctaveInterface* pInterface, int iFlags);
580  inline bool bHaveMethod(const std::string& strName) const;
581  OctaveInterface* GetInterface()const{ return pInterface; }
582  const std::string& GetClass()const{ return strClass; }
583  int GetFlags()const{ return iFlags; }
584 
585 protected:
586  octave_value octObject;
587 
588 private:
589  const std::string strClass;
590  OctaveInterface* const pInterface;
591  const int iFlags;
592  static const std::string strGetConstLawType;
593  static const std::string strUpdate;
594 };
595 
596 template <class T, class Tder>
597 class OctaveConstitutiveLaw
598 : public ConstitutiveLaw<T, Tder>, private OctaveConstitutiveLawBase {
599  typedef ConstitutiveLaw<T, Tder> Base_t;
600 public:
601  OctaveConstitutiveLaw(const std::string& strClass, OctaveInterface* pInterface, int iFlags);
602  virtual ~OctaveConstitutiveLaw(void);
603  ConstLawType::Type GetConstLawType(void) const;
604  virtual ConstitutiveLaw<T, Tder>* pCopy(void) const;
605  virtual std::ostream& Restart(std::ostream& out) const;
606  virtual void Update(const T& mbEps, const T& mbEpsPrime);
607 
608 private:
609  mutable ConstLawType::Type clType;
610 };
611 
612 class OctaveBaseDCR {
613 public:
614  OctaveBaseDCR(void);
615  virtual ~OctaveBaseDCR(void);
616  void Read(const DataManager* pDM, MBDynParser& HP, bool bDeferred = false) const;
617  OctaveInterface* GetInterface(void) const { return pInterface; };
618  const std::string& GetFunction(void) const { return strFunction; };
619  int GetFlags(void) const { return iFlags; };
620 
621 private:
622  // FIXME: mutable needed for ScalarFunctionRead
623  mutable OctaveInterface* pInterface;
624  mutable std::string strFunction;
625  mutable int iFlags;
626 };
627 
628 class OctaveFunctionDCR: public OctaveBaseDCR {
629 public:
630  void Read(const DataManager* pDM, MBDynParser& HP, bool bDeferred = false) const;
631  const octave_value_list& GetArgs()const{ return args; }
632 
633 private:
634  mutable octave_value_list args;
635 };
636 
637 class OctaveDCR : public DriveCallerRead, public OctaveFunctionDCR {
638 public:
639  virtual DriveCaller *
640  Read(const DataManager* pDM, MBDynParser& HP, bool bDeferred);
641 };
642 
643 template <class T>
644 class OctaveTDCR : public TplDriveCallerRead<T>, public OctaveFunctionDCR {
645 public:
646  virtual TplDriveCaller<T> *
647  Read(const DataManager* pDM, MBDynParser& HP);
648 };
649 
650 class DerivativeDCR : public DriveCallerRead {
651 public:
652  virtual DriveCaller *
653  Read(const DataManager* pDM, MBDynParser& HP, bool bDeferred);
654 };
655 
656 class OctaveSFR : public ScalarFunctionRead, public OctaveFunctionDCR {
657 public:
658  const BasicScalarFunction *
659  Read(DataManager* pDM, MBDynParser& HP) const;
660 };
661 
662 template <class T, class Tder>
663 struct OctaveCLR : public ConstitutiveLawRead<T, Tder>, public OctaveBaseDCR {
664  virtual ConstitutiveLaw<T, Tder> *
665  Read(const DataManager* pDM, MBDynParser& HP, ConstLawType::Type& CLType) {
666  ConstitutiveLaw<T, Tder>* pCL = 0;
667 
668  OctaveBaseDCR::Read(pDM, HP);
669 
670  typedef OctaveConstitutiveLaw<T, Tder> L;
671  SAFENEWWITHCONSTRUCTOR(pCL, L, L(GetFunction(), GetInterface(), GetFlags()));
672 
673  CLType = pCL->GetConstLawType();
674 
675  return pCL;
676  };
677 };
678 
679 class OctaveElement: virtual public Elem, public UserDefinedElem {
680 public:
681  OctaveElement(unsigned uLabel, const DofOwner *pDO,
682  DataManager* pDM, MBDynParser& HP);
683  virtual ~OctaveElement(void);
684 
685  virtual void Output(OutputHandler& OH) const;
686  virtual void WorkSpaceDim(integer* piNumRows, integer* piNumCols) const;
688  AssJac(VariableSubMatrixHandler& WorkMat,
689  doublereal dCoef,
690  const VectorHandler& XCurr,
691  const VectorHandler& XPrimeCurr);
693  AssRes(SubVectorHandler& WorkVec,
694  doublereal dCoef,
695  const VectorHandler& XCurr,
696  const VectorHandler& XPrimeCurr);
697  unsigned int iGetNumPrivData(void) const;
698  virtual unsigned int iGetPrivDataIdx(const char *s) const;
699  virtual doublereal dGetPrivData(unsigned int i) const;
700  int iGetNumConnectedNodes(void) const;
701  void GetConnectedNodes(std::vector<const Node *>& connectedNodes) const;
702  void SetValue(DataManager *pDM, VectorHandler& X, VectorHandler& XP,
704  virtual unsigned int iGetNumDof(void) const;
705  virtual DofOrder::Order GetDofType(unsigned int i) const;
706  virtual DofOrder::Order GetEqType(unsigned int i) const;
707  virtual std::ostream& DescribeDof(std::ostream& out, const char *prefix, bool bInitial) const;
708  virtual std::ostream& DescribeEq(std::ostream& out, const char *prefix, bool bInitial) const;
709  virtual void AfterPredict(VectorHandler& X, VectorHandler& XP);
710  virtual void Update(const VectorHandler& XCurr,const VectorHandler& XPrimeCurr);
711  virtual void AfterConvergence(const VectorHandler& X,
712  const VectorHandler& XP);
713  std::ostream& Restart(std::ostream& out) const;
714  virtual unsigned int iGetInitialNumDof(void) const;
715  virtual void
716  InitialWorkSpaceDim(integer* piNumRows, integer* piNumCols) const;
718  InitialAssJac(VariableSubMatrixHandler& WorkMat,
719  const VectorHandler& XCurr);
721  InitialAssRes(SubVectorHandler& WorkVec, const VectorHandler& XCurr);
722  virtual void SetInitialValue(VectorHandler& X);
723  OctaveInterface* GetInterface(void) const { return dcr.GetInterface(); };
724  int GetFlags(void) const { return dcr.GetFlags(); };
725  const std::string& GetClass(void) const { return dcr.GetFunction(); };
726 
727 private:
729  AssMatrix(VariableSubMatrixHandler& WorkMatVar, const octave_value_list& ans, bool bInitial);
731  AssMatrix(VariableSubMatrixHandler& WorkMatVar, const Matrix& Jac, const int32NDArray& ridx, const int32NDArray& cidx, bool bSparse, bool bInitial);
732  bool bHaveMethod(const std::string& strName)const;
733 private:
734  octave_value octObject;
735  octave_value mbdObject;
736  OctaveBaseDCR dcr;
737  enum OctaveMethods_t {
738  HAVE_DEFAULT = 0x0000,
739  HAVE_JACOBIAN = 0x0001,
740  HAVE_UPDATE = 0x0002,
741  HAVE_SET_VALUE = 0x0004,
742  HAVE_PRIVATE_DOF = 0x0008,
743  HAVE_INITIAL_ASSEMBLY = 0x0010,
744  HAVE_SET_INITIAL_VALUE = 0x0020,
745  HAVE_AFTER_CONVERGENCE = 0x0040,
746  HAVE_PRIVATE_DATA = 0x0080,
747  HAVE_CONNECTED_NODES = 0x0100,
748  HAVE_OUTPUT = 0x0200,
749  HAVE_DESCRIBE_DOF = 0x0400,
750  HAVE_DESCRIBE_EQ = 0x0800,
751  HAVE_AFTER_PREDICT = 0x1000,
752  HAVE_RESTART = 0x2000
753  };
754  int haveMethod;
755  octave_object_ptr<ConstVectorHandlerInterface> X;
756  octave_object_ptr<ConstVectorHandlerInterface> XP;
757  octave_object_ptr<OStreamInterface> OS;
758  static const std::string strWorkSpaceDim;
759  static const std::string striGetNumDof;
760  static const std::string strAssRes;
761  static const std::string strAssJac;
762  static const std::string strUpdate;
763  static const std::string strSetValue;
764  static const std::string striGetInitialNumDof;
765  static const std::string strSetInitialValue;
766  static const std::string strInitialAssRes;
767  static const std::string strInitialAssJac;
768  static const std::string strInitialWorkSpaceDim;
769  static const std::string strGetDofType;
770  static const std::string strGetEqType;
771  static const std::string strAfterConvergence;
772  static const std::string striGetNumPrivData;
773  static const std::string striGetPrivDataIdx;
774  static const std::string strdGetPrivData;
775  static const std::string striGetNumConnectedNodes;
776  static const std::string strGetConnectedNodes;
777  static const std::string strOutput;
778  static const std::string strDescribeDof;
779  static const std::string strDescribeEq;
780  static const std::string strAfterPredict;
781  static const std::string strRestart;
782 };
783 
784 class OctaveElementInterface: public MBDynInterface {
785 public:
786  explicit OctaveElementInterface(OctaveInterface* pInterface = OctaveInterface::GetInterface(), OctaveElement* pElem = 0);
787  virtual ~OctaveElementInterface(void);
788  virtual void print(std::ostream& os, bool pr_as_read_syntax = false) const;
789 
790 protected:
791  BEGIN_METHOD_TABLE_DECLARE()
792  METHOD_DECLARE(GetLabel)
793  METHOD_DECLARE(iGetFirstIndex)
794  END_METHOD_TABLE_DECLARE()
795 
796 private:
797  DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA
798  DECLARE_OCTAVE_ALLOCATOR
799 
800 private:
801  OctaveElement* const pElem;
802 };
803 
804 OctaveInterface* OctaveInterface::pOctaveInterface = 0;
805 int OctaveInterface::iRefCount = 0;
806 const std::string OctaveInterface::strADFunc("mbdyn_derivative");
807 const std::string OctaveInterface::strIsMethod("ismethod");
808 bool OctaveInterface::bHaveADPackage = false;
809 
810 OctaveInterface::OctaveInterface(const DataManager* pDM, MBDynParser* pHP)
811 : bFirstCall(true),
812 bEmbedFileDirty(false),
813 octDM(new DataManagerInterface(this, pDM)),
814 octHP(new MBDynParserInterface(this, pHP)),
815 pDM(pDM),
816 pHP(pHP)
817 {
818  TRACE("constructor");
819 
820  ASSERT(pDM != 0);
821  ASSERT(pHP != 0);
822  ASSERT(pOctaveInterface == 0);
823 
824  pOctaveInterface = this;
825 
826  const int nmax_args = 4;
827  const int nmax_char = 14;
828  char args[nmax_char] = "octave"; //FIXME: this might be unsafe
829  char *argv[nmax_args] = { &args[0] } ;
830  int argc = 1;
831 
832  char* p = strchr(args, '\0') + 1;
833 
834  if (silent_output) {
835  strncat(p, "-q", &args[nmax_char] - p - 2);
836  if (argc < nmax_args - 1)
837  argv[argc++] = p;
838  p = strchr(p, '\0') + 1;
839  }
840 
841  if (pedantic_output) {
842  strncat(p, "-V", &args[nmax_char] - p - 2);
843  if (argc < nmax_args - 1)
844  argv[argc++] = p;
845  p = strchr(p, '\0') + 1;
846  }
847 
848 #if defined(DEBUG)
849  int i;
850  char** pp;
851  for (i = 0, pp = argv; *pp; ++pp, ++i) {
852  silent_cerr("octave>argv[" << i << "]=\"" << *pp << "\"" << std::endl);
853  }
854 #endif
855 
856  octave_main(argc, argv, 1);
857 
858  LoadADPackage();
859 }
860 
861 void OctaveInterface::exit(int)
862 {
863  silent_cerr("octave interface has been uninitialized" << std::endl);
864 }
865 
866 OctaveInterface::~OctaveInterface(void)
867 {
868  TRACE("destructor");
869 
870  ASSERT(this == pOctaveInterface);
871 
872  octave_exit = &OctaveInterface::exit;
873 
874 #if defined(HAVE_DO_OCTAVE_ATEXIT)
875  do_octave_atexit();
876 #elif defined(HAVE_CLEAN_UP_AND_EXIT)
877  clean_up_and_exit(0, true);
878 #else
879  #warning "do_octave_atexit() and clean_up_and_exit() are not defined"
880 #endif
881 
882  pOctaveInterface = 0;
883 }
884 
885 bool
886 OctaveInterface::LoadADPackage(void)
887 {
888  octave_value_list args;
889  args.append(octave_value("load"));
890  args.append(octave_value("ad"));
891 
892  feval("pkg", args, 0);
893 
894  bHaveADPackage = true; // Use finite differences in var/mbdyn_derivative.m if AD is not available
895 
896  if ((error_state != 0)) {
897  silent_cerr("warning: octave package for automatic forward differentiation is not available" << std::endl);
898  error_state = 0; // ignore error
899  }
900 
901  return bHaveADPackage;
902 }
903 
904 OctaveInterface *
905 OctaveInterface::CreateInterface(const DataManager* pDM, MBDynParser* pHP)
906 {
907  ++iRefCount;
908 
909  if (pOctaveInterface) {
910  return pOctaveInterface;
911  }
912 
913  return new OctaveInterface(pDM, pHP);
914 }
915 
916 void
917 OctaveInterface::Destroy(void)
918 {
919  ASSERT(iRefCount >= 1);
920 
921  if (0 == --iRefCount) {
922  delete this;
923  }
924 }
925 
926 void
927 OctaveInterface::UpdateOctaveVariables(void)
928 {
929  typedef Table::VM::const_iterator iterator;
930 
931  const Table& symbolTable = GetDataManager()->GetMathParser().GetSymbolTable();
932 
933  for (iterator it = symbolTable.begin(); it != symbolTable.end(); ++it)
934  {
935  const std::string& mbName(it->first);
936  const NamedValue*const namedValue = it->second;
937  const TypedValue mbValue(namedValue->GetVal());
938  octave_value octValue;
939 
940  if (!ConvertMBDynToOctave(mbValue, octValue)) {
941  silent_cerr("octave error: data type \"" << mbValue.GetTypeName() << "\" of variable \"" << mbName << "\": not handled in switch statement " << std::endl);
942  ASSERT(0);
944  }
945 
946  set_global_value(mbName, octValue);
947  }
948 }
949 
950 void
951 OctaveInterface::UpdateMBDynVariables(void)
952 {
953  typedef Table::VM::const_iterator iterator;
954 
955  Table& symbolTable = GetDataManager()->GetMathParser().GetSymbolTable();
956 
957  for (iterator it = symbolTable.begin(); it != symbolTable.end(); ++it)
958  {
959  Var* const varValue = dynamic_cast<Var*>(it->second);
960 
961  if (!varValue || varValue->Const()) {
962  continue;
963  }
964 
965  const std::string& mbName(it->first);
966  //static const std::string type("any");
967 
968  const octave_value octValue(get_global_value(mbName));
969 
970  if (!octValue.is_defined()) {
971  pedantic_cerr("octave warning: global variable " << mbName << " is not defined in octave" << std::endl);
972  continue;
973  }
974 
975  TypedValue mbValue;
976 
977  if (!ConvertOctaveToMBDyn(octValue, mbValue)) {
978  silent_cerr("octave error: data type \"" << octValue.type_name() << "\" of variable \"" << mbName << "\": can not be converted into MBDyn format " << std::endl);
980  }
981 
982  varValue->SetVal(mbValue);
983  }
984 }
985 
986 bool
987 OctaveInterface::AddOctaveSearchPath(const std::string& path)
988 {
989  feval("addpath", octave_value(path));
990 
991  return error_state == 0;
992 }
993 
994 bool
995 OctaveInterface::ConvertMBDynToOctave(const TypedValue& mbValue, octave_value& octValue)
996 {
997  switch (mbValue.GetType()) {
999  octValue = mbValue.GetBool();
1000  break;
1001  case TypedValue::VAR_INT:
1002  octValue = mbValue.GetInt();
1003  break;
1004  case TypedValue::VAR_REAL:
1005  octValue = mbValue.GetReal();
1006  break;
1008  octValue = mbValue.GetString();
1009  break;
1010  default:
1011  return false;
1012  }
1013 
1014  return true;
1015 }
1016 
1017 bool
1018 OctaveInterface::ConvertMBDynToOctave(doublereal mbValue, octave_value& octValue)
1019 {
1020  // This function is especially useful for template based constitutive laws
1021  octValue = mbValue;
1022 
1023  return true;
1024 }
1025 
1026 bool
1027 OctaveInterface::ConvertMBDynToOctave(const Vec3& mbValue, octave_value& octValue)
1028 {
1029  return ConvertMBDynToOctave(mbValue, octValue, 3);
1030 }
1031 
1032 bool
1033 OctaveInterface::ConvertMBDynToOctave(const Vec6& mbValue, octave_value& octValue)
1034 {
1035  return ConvertMBDynToOctave(mbValue, octValue, 6);
1036 }
1037 
1038 bool
1039 OctaveInterface::ConvertMBDynToOctave(const Mat3x3& mbValue, octave_value& octValue)
1040 {
1041  return ConvertMBDynToOctave(mbValue, octValue, 3, 3);
1042 }
1043 
1044 bool
1045 OctaveInterface::ConvertMBDynToOctave(const Mat6x6& mbValue, octave_value& octValue)
1046 {
1047  return ConvertMBDynToOctave(mbValue, octValue, 6, 6);
1048 }
1049 
1050 bool
1051 OctaveInterface::ConvertOctaveToMBDyn(const ColumnVector& octValue, doublereal& mbValue)
1052 {
1053  if (octValue.length() != 1) {
1054  silent_cerr("octave error: invalid vector size " << octValue.length() << " expected 1" << std::endl);
1055  return false;
1056  }
1057 
1058  mbValue = octValue(0);
1059 
1060  return true;
1061 }
1062 
1063 bool
1064 OctaveInterface::ConvertOctaveToMBDyn(const ColumnVector& octValue, Vec3& mbValue)
1065 {
1066  return ConvertOctaveToMBDyn(octValue, mbValue, 3);
1067 }
1068 
1069 bool
1070 OctaveInterface::ConvertOctaveToMBDyn(const ColumnVector& octValue, Vec6& mbValue)
1071 {
1072  return ConvertOctaveToMBDyn(octValue, mbValue, 6);
1073 }
1074 
1075 template<typename T>
1076 bool
1077 OctaveInterface::ConvertOctaveToMBDyn(const ColumnVector& octValue, T& mbValue, int rows)
1078 {
1079  if (octValue.length() != rows) {
1080  silent_cerr("octave error: invalid vector size " << octValue.length() << " expected " << rows << std::endl);
1081  return false;
1082  }
1083 
1084  for (int i = 0; i < rows; ++i) {
1085  mbValue(i + 1) = octValue(i);
1086  }
1087 
1088  return true;
1089 }
1090 
1091 template<typename T>
1092 bool
1093 OctaveInterface::ConvertOctaveToMBDyn(const Matrix& octValue, T& mbValue, int rows)
1094 {
1095  if (octValue.rows() != rows || octValue.columns() != 1) {
1096  silent_cerr("octave error: invalid matrix size " << octValue.rows() << "x" << octValue.columns() << " expected " << rows << "x1" << std::endl);
1097  return false;
1098  }
1099 
1100  for (int i = 0; i < rows; ++i) {
1101  mbValue(i + 1) = octValue(i, 0);
1102  }
1103 
1104  return true;
1105 }
1106 
1107 template<typename T>
1108 bool
1109 OctaveInterface::ConvertOctaveToMBDyn(const Matrix& octValue, T& mbValue, int rows, int cols)
1110 {
1111  if (octValue.rows() != rows || octValue.columns() != cols) {
1112  silent_cerr("octave error: invalid matrix size " << octValue.rows() << "x" << octValue.columns() << " expected " << rows << "x" << cols << std::endl);
1113  return false;
1114  }
1115 
1116  for (int i = 0; i < rows; ++i) {
1117  for (int j = 0; j < cols; ++j) {
1118  mbValue(i + 1, j + 1) = octValue(i, j);
1119  }
1120  }
1121 
1122  return true;
1123 }
1124 
1125 template<typename T>
1126 bool
1127 OctaveInterface::ConvertMBDynToOctave(const T& mbValue, octave_value& octValue, int rows)
1128 {
1129  ColumnVector V(rows);
1130 
1131  for (int i = 0; i < rows; ++i) {
1132  V(i) = mbValue(i + 1);
1133  }
1134 
1135  octValue = V;
1136 
1137  return true;
1138 }
1139 
1140 template<typename T>
1141 bool
1142 OctaveInterface::ConvertMBDynToOctave(const T& mbValue, octave_value& octValue, int rows, int cols)
1143 {
1144  Matrix M(rows, cols);
1145 
1146  for (int i = 0; i < rows; ++i) {
1147  for (int j = 0; j < cols; ++j) {
1148  M(i, j) = mbValue(i + 1, j + 1);
1149  }
1150  }
1151 
1152  octValue = M;
1153 
1154  return true;
1155 }
1156 
1157 bool
1158 OctaveInterface::ConvertOctaveToMBDyn(const Matrix& octValue, doublereal& mbValue)
1159 {
1160  if (octValue.rows() != 1 || octValue.columns() != 1) {
1161  silent_cerr("octave error: invalid matrix size " << octValue.rows() << "x" << octValue.columns() << " expected 1x1" << std::endl);
1162  return false;
1163  }
1164 
1165  mbValue = octValue(0, 0);
1166 
1167  return true;
1168 }
1169 
1170 bool
1171 OctaveInterface::ConvertOctaveToMBDyn(const Matrix& octValue, Vec3& mbValue)
1172 {
1173  return ConvertOctaveToMBDyn(octValue, mbValue, 3);
1174 }
1175 
1176 bool
1177 OctaveInterface::ConvertOctaveToMBDyn(const Matrix& octValue, Vec6& mbValue)
1178 {
1179  return ConvertOctaveToMBDyn(octValue, mbValue, 6);
1180 }
1181 
1182 bool
1183 OctaveInterface::ConvertOctaveToMBDyn(const Matrix& octValue, Mat3x3& mbValue)
1184 {
1185  return ConvertOctaveToMBDyn(octValue, mbValue, 3, 3);
1186 }
1187 
1188 bool
1189 OctaveInterface::ConvertOctaveToMBDyn(const Matrix& octValue, Mat6x6& mbValue)
1190 {
1191  return ConvertOctaveToMBDyn(octValue, mbValue, 6, 6);
1192 }
1193 
1194 bool
1195 OctaveInterface::ConvertOctaveToMBDyn(const octave_value& octValue, TypedValue& mbValue)
1196 {
1197  if (!octValue.is_scalar_type()) {
1198  return false;
1199  }
1200 
1201  if (octValue.is_bool_scalar()) {
1202  mbValue = TypedValue(octValue.bool_value());
1203  } else if (octValue.is_int32_type()) {
1204  mbValue = TypedValue(static_cast<int32_t>(octValue.int32_scalar_value()));
1205  } else if (octValue.is_real_scalar()) {
1206  mbValue = TypedValue(octValue.scalar_value());
1207  } else if (octValue.is_string()) {
1208  mbValue = TypedValue(octValue.string_value());
1209  } else {
1210  return false;
1211  }
1212 
1213  return true;
1214 }
1215 
1216 bool OctaveInterface::ConvertOctaveToMBDyn(const octave_value& octValue, doublereal& mbValue)
1217 {
1218  if (!octValue.is_real_scalar()) {
1219  return false;
1220  }
1221 
1222  mbValue = octValue.scalar_value();
1223 
1224  return true;
1225 }
1226 
1227 octave_value_list
1228 OctaveInterface::EvalFunction(const std::string& func, doublereal dVar, int nargout, int flags)
1229 {
1230  octave_value_list args = octave_value(dVar);
1231 
1232  if (flags & PASS_DATA_MANAGER) {
1233  args.append(octDM);
1234  }
1235 
1236  return EvalFunction(func, args, nargout, flags);
1237 }
1238 
1239 octave_value_list
1240 OctaveInterface::EvalFunction(const std::string& func, const octave_value_list& args, int nargout, int flags)
1241 {
1242  if ((flags & UPDATE_OCTAVE_VARIABLES) || bFirstCall) {
1243  UpdateOctaveVariables();
1244  }
1245 
1246  if (bFirstCall) {
1247  // octave .oct files by default are installed here
1248  if (!AddOctaveSearchPath(OCTAVEBINPATH)) {
1249  silent_cerr("OctaveInterface error: addpath(\"" << BINPATH << "\") failed" << std::endl);
1251  }
1252 
1253  // octave .m files by default are installed here
1254  if (!AddOctaveSearchPath(OCTAVEPATH)) {
1255  silent_cerr("OctaveInterface error: addpath(\"" << OCTAVEPATH << "\") failed" << std::endl);
1257  }
1258  }
1259 
1260  bFirstCall = false;
1261 
1262  if ( bEmbedFileDirty ) {
1263  for (EmbedFileNameIter_t pFile = strEmbedFileNames.begin();
1264  pFile != strEmbedFileNames.end(); ++pFile )
1265  {
1266  if (!pFile->second) {
1267  feval("source", octave_value(pFile->first));
1268 
1269  if (error_state) {
1271  }
1272 
1273  pFile->second = true;
1274  }
1275  }
1276  bEmbedFileDirty = false;
1277  }
1278 
1279  octave_value_list ans = feval(func, args, nargout);
1280 
1281  if (error_state) {
1282  // An error message has been displayed by octave
1283  // There is no need to output further error information
1285  }
1286 
1287  if (flags & OPTIONAL_OUTPUT_ARGS) {
1288  if (ans.length() > nargout) {
1289  silent_cerr("octave error: function \"" << func << "\" returned " << ans.length() << " output parameters\n"
1290  "expected maximum " << nargout << " output parameters" << std::endl);
1292  }
1293  }
1294  else if (ans.length() != nargout) {
1295  silent_cerr("octave error: function \"" << func << "\" returned " << ans.length() << " output parameters\n"
1296  "expected " << nargout << " output parameters" << std::endl);
1298  }
1299 
1300  for (int i = 0; i < ans.length(); ++i) {
1301  if (!ans(i).is_defined()) {
1302  silent_cerr("octave error: result " << i + 1 << " of function \"" << func << "\" is undefined" << std::endl);
1304  }
1305  }
1306 
1307  if (flags & UPDATE_MBDYN_VARIABLES) {
1308  UpdateMBDynVariables();
1309  }
1310 
1311  return ans;
1312 }
1313 
1314 octave_value_list
1315 OctaveInterface::EvalFunctionDerivative(const std::string& func, const octave_value_list& args, int flags)
1316 {
1317  TRACE("func=" << func);
1318  TRACE("flags=" << flags);
1319 
1320  octave_value_list derivArgs = octave_value(func);
1321  derivArgs.append(args);
1322 
1323  return EvalFunction(strADFunc, derivArgs, 2, flags);
1324 }
1325 
1326 octave_value_list
1327 OctaveInterface::EvalFunctionDerivative(const std::string& func, doublereal dVar, int flags)
1328 {
1329  TRACE("func=" << func);
1330  TRACE("flags=" << flags);
1331 
1332  octave_value_list args = octave_value(dVar);
1333 
1334  if (flags & PASS_DATA_MANAGER) {
1335  args.append(octDM);
1336  }
1337 
1338  return EvalFunctionDerivative(func, args, flags);
1339 }
1340 
1341 inline doublereal
1342 OctaveInterface::EvalScalarFunction(const std::string& func, doublereal dVar, int flags)
1343 {
1344  return EvalScalarFunction(func, octave_value(dVar), flags);
1345 }
1346 
1347 inline doublereal
1348 OctaveInterface::EvalScalarFunction(const std::string& func, const octave_value_list& args, int flags)
1349 {
1350  const octave_value_list ans = EvalFunction(func, args, 1, flags);
1351 
1352  if (!ans(0).is_real_scalar()) {
1353  silent_cerr("octave error: result of function \"" << func << "\" is not a scalar value" << std::endl);
1355  }
1356 
1357  return ans(0).scalar_value();
1358 }
1359 
1360 inline doublereal
1361 OctaveInterface::EvalScalarFunctionDerivative(const std::string& func, doublereal dVar, int flags)
1362 {
1363  return EvalScalarFunctionDerivative(func, octave_value(dVar), flags);
1364 }
1365 
1366 inline doublereal
1367 OctaveInterface::EvalScalarFunctionDerivative(const std::string& func, const octave_value_list& args, int flags)
1368 {
1369  TRACE("func=" << func);
1370  TRACE("flags=" << flags);
1371 
1372  octave_value_list ans = EvalFunctionDerivative(func, args, flags);
1373 
1374  ASSERT(ans.length() == 2);
1375 
1376  if (!ans(1).is_real_scalar()) {
1377  silent_cerr("octave error: derivative of function \"" << func << "\" is not a scalar value" << std::endl);
1379  }
1380 
1381  return ans(1).scalar_value();
1382 }
1383 
1384 template <typename T>
1385 inline void
1386 OctaveInterface::EvalMatrixFunction(const std::string& func, T& res, doublereal dVar, int flags)
1387 {
1388  EvalMatrixFunction(func, res, octave_value(dVar), flags);
1389 }
1390 
1391 template <typename T>
1392 inline void
1393 OctaveInterface::EvalMatrixFunction(const std::string& func, T& res, const octave_value_list& args, int flags)
1394 {
1395  const octave_value_list ans = EvalFunction(func, args, 1, flags);
1396 
1397  ASSERT(ans.length() == 1);
1398  // accept also scalar values for scalar template drive caller
1399  if (!(ans(0).is_real_matrix() || ans(0).is_real_scalar())) {
1400  silent_cerr("octave error: result of function \"" << func << "\" is not a matrix value" << std::endl);
1402  }
1403 
1404  const Matrix A(ans(0).matrix_value());
1405 
1406  if (!ConvertOctaveToMBDyn(A, res)) {
1407  silent_cerr("octave error: conversion of octave matrix failed!" << std::endl);
1409  }
1410 }
1411 
1412 template <typename T>
1413 inline void
1414 OctaveInterface::EvalMatrixFunctionDerivative(const std::string& func, T& res, doublereal dVar, int flags)
1415 {
1416  EvalMatrixFunctionDerivative(func, res, octave_value(dVar), flags);
1417 }
1418 
1419 template<typename T>
1420 inline void
1421 OctaveInterface::EvalMatrixFunctionDerivative(const std::string& func, T& res, const octave_value_list& args, int flags)
1422 {
1423  const octave_value_list ans = EvalFunctionDerivative(func, args, flags);
1424 
1425  ASSERT(ans.length() == 2);
1426 
1427  if (!(ans(0).is_real_matrix() || ans(0).is_real_scalar())) {
1428  silent_cerr("octave error: result of function \"" << func << "\" is not a matrix value" << std::endl);
1430  }
1431 
1432  if (!(ans(1).is_real_matrix() || ans(1).is_real_scalar())) {
1433  silent_cerr("octave error: result of derivative of function \"" << func << "\" is not a matrix value" << std::endl);
1435  }
1436 
1437  const Matrix A(ans(1).matrix_value());
1438 
1439  if (!ConvertOctaveToMBDyn(A, res)) {
1440  silent_cerr("octave error: conversion of octave matrix failed!" << std::endl);
1442  }
1443 }
1444 
1445 void
1446 OctaveInterface::AddEmbedFileName(const std::string& strFile)
1447 {
1448  strEmbedFileNames.insert(std::pair<std::string, bool>(strFile,false));
1449  bEmbedFileDirty = true;
1450 }
1451 
1452 bool OctaveInterface::bHaveMethod(const octave_value& octObject, const std::string& strClass, const std::string& strName)
1453 {
1454  octave_value_list args(octObject);
1455  args.append(octave_value(strName));
1456 
1457  octave_value_list ans = EvalFunction(strIsMethod, args, 1);
1458 
1459  ASSERT(ans.length() == 1);
1460 
1461  if (!ans(0).is_bool_scalar()) {
1462  silent_cerr("octave error: unexpected error in function " << strIsMethod << std::endl);
1464  }
1465 
1466  return ans(0).bool_value(true);
1467 }
1468 
1469 octave_value_list OctaveInterface::MakeArgList(doublereal dVar, const octave_value_list& args, int iFlags) const
1470 {
1471  octave_value_list fargs = octave_value(dVar);
1472  fargs.append(args);
1473 
1474  if (iFlags & PASS_DATA_MANAGER) {
1475  fargs.append(GetDataManagerInterface());
1476  }
1477 
1478  return fargs;
1479 }
1480 
1481 MBDynInterface::MBDynInterface(OctaveInterface* pInterface)
1482 : pInterface(pInterface)
1483 {
1484  TRACE("constructor");
1485  ASSERT(pInterface != 0);
1486 }
1487 
1488 MBDynInterface::~MBDynInterface(void)
1489 {
1490  TRACE("destructor");
1491 }
1492 
1493 void
1494 MBDynInterface::print(std::ostream& os, bool pr_as_read_syntax) const
1495 {
1496  os << type_name() << ": MBDyn version" << VERSION << std::endl;
1497 }
1498 
1499 METHOD_DEFINE(MBDynInterface, GetVersion, args, nargout)
1500 {
1501  octave_value octValue;
1502 
1503  if (args.length() != 0) {
1504  error("%s: invalid number of arguments: %ld",
1505  type_name().c_str(),
1506  long(args.length()));
1507  return octValue;
1508  }
1509 
1510  octValue = VERSION;
1511 
1512  return octValue;
1513 }
1514 
1515 BEGIN_METHOD_TABLE(MBDynInterface, octave_object)
1516  METHOD_DISPATCH(MBDynInterface, GetVersion)
1517 END_METHOD_TABLE()
1518 
1519 DEFINE_OCTAVE_ALLOCATOR(MBDynInterface);
1520 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA(MBDynInterface, "MBDyn", "MBDyn");
1521 
1522 ConstVectorHandlerInterface::ConstVectorHandlerInterface(OctaveInterface* pInterface, const VectorHandler* X)
1523  :MBDynInterface(pInterface),
1524  pX(const_cast<VectorHandler*>(X))
1525 {
1526  TRACE("constructor");
1527 }
1528 
1529 ConstVectorHandlerInterface::~ConstVectorHandlerInterface()
1530 {
1531  TRACE("destructor");
1532 }
1533 
1534 void
1535 ConstVectorHandlerInterface::print(std::ostream& os, bool pr_as_read_syntax) const
1536 {
1537  if ( pX == 0 ) {
1538  error("%s: not connected", type_name().c_str());
1539  return;
1540  }
1541 
1542  for ( int i = 1; i <= pX->iGetSize(); ++i ) {
1543  os << '\t' << pX->dGetCoef(i) << std::endl;
1544  }
1545 }
1546 
1547 octave_value ConstVectorHandlerInterface::operator()(const octave_value_list& idx) const
1548 {
1549  if ( pX == 0 ) {
1550  error("%s: not connected", type_name().c_str());
1551  return octave_value();
1552  }
1553 
1554  if (idx.length() != 1) {
1555  error("%s: invalid number of indices %ld", type_name().c_str(), long(idx.length()));
1556  return octave_value();
1557  }
1558 
1559  const integer iSize = pX->iGetSize();
1560 
1561  ColumnVector X;
1562 
1563  if (idx(0).is_magic_colon()) {
1564  X.resize(iSize);
1565 
1566  for (int i = 0; i < X.length(); ++i) {
1567  X(i) = (*pX)(i + 1);
1568  }
1569  } else {
1570  if (!(idx(0).is_range()
1571  || (idx(0).is_integer_type()
1572  && (idx(0).rows() == 1 || idx(0).columns() == 1)))) {
1573  error("%s: invalid index type %dx%d (%s)\n"
1574  "expected integer vector",
1575  type_name().c_str(),
1576  idx(0).rows(),
1577  idx(0).columns(),
1578  idx(0).type_name().c_str());
1579  return octave_value();
1580  }
1581 
1582  const int32NDArray iRow(idx(0).int32_array_value());
1583 
1584  if (error_state) {
1585  return octave_value();
1586  }
1587 
1588  X.resize(iRow.length());
1589 
1590  for (int i = 0; i < iRow.length(); ++i) {
1591  if (int32_t(iRow(i)) < 1 || int32_t(iRow(i)) > iSize) {
1592  error("%s: index %d out of range [%d-%d]",
1593  type_name().c_str(),
1594  int32_t(iRow(i)),
1595  1,
1596  int32_t(iSize));
1597 
1598  return octave_value();
1599  }
1600 
1601  X(i) = (*pX)(int32_t(iRow(i)));
1602  }
1603  }
1604 
1605  return octave_value(X);
1606 }
1607 
1608 dim_vector ConstVectorHandlerInterface::dims() const
1609 {
1610  if ( pX == 0 ) {
1611  error("%s: not connected", type_name().c_str());
1612  return dim_vector(-1, -1);
1613  }
1614 
1615  return dim_vector(pX->iGetSize(), 1);
1616 }
1617 
1618 METHOD_DEFINE(ConstVectorHandlerInterface, dGetCoef, args, nargout)
1619 {
1620  if ( pX == 0 ) {
1621  error("%s: not connected", type_name().c_str());
1622  return octave_value();
1623  }
1624 
1625  if ( args.length() != 1 ) {
1626  error("%s: invalid number of arguments %ld\n"
1627  "one argument expected",
1628  type_name().c_str(),
1629  long(args.length()));
1630  return octave_value();
1631  }
1632 
1633  if ( !(args(0).is_integer_type() && args(0).is_scalar_type()) ) {
1634  error("%s: invalid argument type \"%s\"\n"
1635  "iRow must be an integer",
1636  type_name().c_str(),
1637  args(0).type_name().c_str());
1638  return octave_value();
1639  }
1640 
1641  const integer iRow = args(0).int32_scalar_value();
1642 
1643  if ( iRow < 1 || iRow > pX->iGetSize() ) {
1644  error("VectorHandler: index %ld out of range [%ld-%ld]", long(iRow), 1L, long(pX->iGetSize()));
1645  return octave_value();
1646  }
1647 
1648  return octave_value(pX->dGetCoef(iRow));
1649 }
1650 
1651 METHOD_DEFINE(ConstVectorHandlerInterface, GetVec, args, nargout)
1652 {
1653  if (pX == 0) {
1654  error("%s: not connected", type_name().c_str());
1655  return octave_value();
1656  }
1657 
1658  if (args.length() != 1) {
1659  error("%s: invalid number of arguments %ld\n"
1660  "one argument expected",
1661  type_name().c_str(),
1662  long(args.length()));
1663  return octave_value();
1664  }
1665 
1666  if (!(args(0).is_integer_type() && args(0).columns() == 1)) {
1667  error("%s: invalid argument type \"%s\"\n"
1668  "iRow must be an integer column vector",
1669  type_name().c_str(),
1670  args(0).type_name().c_str());
1671  return octave_value();
1672  }
1673 
1674  const int32NDArray iRow(args(0).int32_array_value());
1675 
1676  ColumnVector X(iRow.length());
1677 
1678  const int32_t iSize = pX->iGetSize();
1679 
1680  for (int i = 0; i < iRow.length(); ++i) {
1681  if (int32_t(iRow(i)) < 1 || int32_t(iRow(i)) > iSize) {
1682  error("%s: index %d out of range [%d-%d]",
1683  type_name().c_str(),
1684  int32_t(iRow(i)),
1685  1,
1686  iSize);
1687 
1688  return octave_value();
1689  }
1690 
1691  X(i) = pX->dGetCoef(iRow(i));
1692  }
1693 
1694  return octave_value(X);
1695 }
1696 
1697 METHOD_DEFINE(ConstVectorHandlerInterface, iGetSize, args, nargout)
1698 {
1699  if ( pX == 0 ) {
1700  error("%s: not connected", type_name().c_str());
1701  return octave_value();
1702  }
1703 
1704  if ( args.length() != 0 ) {
1705  error("%s: invalid number of arguments %ld\n"
1706  "no arguments expected",
1707  type_name().c_str(),
1708  long(args.length()));
1709  return octave_value();
1710  }
1711 
1712  return octave_value(octave_int<integer>(pX->iGetSize()));
1713 }
1714 
1715 BEGIN_METHOD_TABLE(ConstVectorHandlerInterface, MBDynInterface)
1716  METHOD_DISPATCH(ConstVectorHandlerInterface, dGetCoef)
1717  METHOD_DISPATCH(ConstVectorHandlerInterface, GetVec)
1718  METHOD_DISPATCH(ConstVectorHandlerInterface, iGetSize)
1719 END_METHOD_TABLE()
1720 
1721 DEFINE_OCTAVE_ALLOCATOR(ConstVectorHandlerInterface);
1722 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA(ConstVectorHandlerInterface, "ConstVectorHandler", "ConstVectorHandler");
1723 
1724 VectorHandlerInterface::VectorHandlerInterface(OctaveInterface* pInterface, VectorHandler* X)
1725  :ConstVectorHandlerInterface(pInterface, X)
1726 {
1727  TRACE("constructor");
1728 }
1729 
1730 VectorHandlerInterface::~VectorHandlerInterface()
1731 {
1732  TRACE("destructor");
1733 }
1734 
1735 METHOD_DEFINE(VectorHandlerInterface, PutCoef, args, nargout)
1736 {
1737  if ( pX == 0 ) {
1738  error("%s: not connected", type_name().c_str());
1739  return octave_value();
1740  }
1741 
1742  if ( args.length() != 2 ) {
1743  error("%s: invalid number of arguments %ld\n"
1744  "two arguments expected",
1745  type_name().c_str(),
1746  long(args.length()));
1747 
1748  return octave_value();
1749  }
1750 
1751  if (!(args(0).is_scalar_type() && args(0).is_integer_type())) {
1752  error("%s: invalid argument type (%s)\n"
1753  "iRow must be an integer",
1754  type_name().c_str(),
1755  args(0).type_name().c_str());
1756 
1757  return octave_value();
1758  }
1759 
1760  const integer iRow = args(0).int32_scalar_value();
1761 
1762  if ( iRow < 1 || iRow > pX->iGetSize() ) {
1763  error("%s: index %ld out of range [%ld-%ld]",
1764  type_name().c_str(),
1765  long(iRow),
1766  1L,
1767  long(pX->iGetSize()));
1768  return octave_value();
1769  }
1770 
1771  if ( !args(1).is_scalar_type() ) {
1772  error("%s: invalid argument type \"%s\"\n"
1773  "dCoef must be a scalar value",
1774  type_name().c_str(),
1775  args(1).type_name().c_str());
1776  return octave_value();
1777  }
1778 
1779  const doublereal dCoef = args(1).scalar_value();
1780 
1781  pX->PutCoef(iRow, dCoef);
1782 
1783  return octave_value();
1784 }
1785 
1786 METHOD_DEFINE(VectorHandlerInterface, PutVec, args, nargout)
1787 {
1788  if ( pX == 0 ) {
1789  error("%s: not connected", type_name().c_str());
1790  return octave_value();
1791  }
1792 
1793  if ( args.length() != 2 ) {
1794  error("%s: invalid number of arguments %ld\n"
1795  "two arguments expected",
1796  type_name().c_str(),
1797  long(args.length()));
1798 
1799  return octave_value();
1800  }
1801 
1802  if (!(args(0).is_integer_type() && args(0).columns() == 1)) {
1803  error("%s: invalid argument type %ldx%ld (%s)\n"
1804  "iRow must be an integer column vector",
1805  type_name().c_str(),
1806  long(args(0).rows()),
1807  long(args(0).columns()),
1808  args(0).type_name().c_str());
1809 
1810  return octave_value();
1811  }
1812 
1813  const int32NDArray iRow(args(0).int32_array_value());
1814 
1815  if (!(args(1).is_real_matrix()
1816  && args(1).columns() == 1
1817  && args(1).rows() == iRow.length())) {
1818  error("%s: invalid argument type %ldx%ld (%s)\n"
1819  "X must be a real column vector "
1820  "with the same number of rows like iRow",
1821  type_name().c_str(),
1822  long(args(1).rows()),
1823  long(args(1).columns()),
1824  args(1).type_name().c_str());
1825 
1826  return octave_value();
1827  }
1828 
1829  const ColumnVector X(args(1).column_vector_value());
1830 
1831  ASSERT(X.length() == iRow.length());
1832 
1833  const int32_t iSize = pX->iGetSize();
1834 
1835  for (int i = 0; i < iRow.length(); ++i) {
1836  if (int32_t(iRow(i)) < 1 || int32_t(iRow(i)) > iSize) {
1837  error("%s: index %d out of range [%d-%d]",
1838  type_name().c_str(),
1839  int32_t(iRow(i)),
1840  1,
1841  iSize);
1842  return octave_value();
1843  }
1844 
1845  pX->PutCoef(iRow(i), X(i));
1846  }
1847 
1848  return octave_value();
1849 }
1850 
1851 BEGIN_METHOD_TABLE(VectorHandlerInterface, ConstVectorHandlerInterface)
1852  METHOD_DISPATCH(VectorHandlerInterface, PutCoef)
1853  METHOD_DISPATCH(VectorHandlerInterface, PutVec)
1854 END_METHOD_TABLE()
1855 
1856 DEFINE_OCTAVE_ALLOCATOR(VectorHandlerInterface);
1857 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA(VectorHandlerInterface, "VectorHandler", "VectorHandler");
1858 
1859 const std::string OStreamInterface::strsprintf("sprintf");
1860 
1861 OStreamInterface::OStreamInterface(OctaveInterface* pInterface, std::ostream* pOS)
1862 : MBDynInterface(pInterface),
1863  pOS(pOS)
1864 {
1865 
1866 }
1867 
1868 OStreamInterface::~OStreamInterface()
1869 {
1870 
1871 }
1872 
1873 METHOD_DEFINE(OStreamInterface, printf, args, nargout)
1874 {
1875  if (!pOS) {
1876  error("ostream: not connected");
1877  return octave_value();
1878  }
1879 
1880  const octave_value_list ans = feval(strsprintf, args, 1);
1881 
1882  if (error_state) {
1883  return octave_value();
1884  }
1885 
1886  if (!(ans.length() >= 1 && ans(0).is_string())) { // sprintf returns more than one output argument
1887  error("ostream: %s failed", strsprintf.c_str());
1888  return octave_value();
1889  }
1890 
1891  const std::string str(ans(0).string_value());
1892 
1893  for (std::string::const_iterator p = str.begin(); p != str.end(); ++p) {
1894  *pOS << *p;
1895 
1896  if (*p == '\n') // conform to default behavior in C and octave
1897  pOS->flush();
1898  }
1899 
1900  return octave_value(octave_int<size_t>(str.length()));
1901 }
1902 
1903 BEGIN_METHOD_TABLE(OStreamInterface, MBDynInterface)
1904  METHOD_DISPATCH(OStreamInterface, printf)
1905 END_METHOD_TABLE()
1906 
1907 DEFINE_OCTAVE_ALLOCATOR(OStreamInterface);
1908 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA(OStreamInterface, "ostream", "ostream");
1909 
1910 
1911 SimulationEntityInterface::SimulationEntityInterface(OctaveInterface* pInterface)
1912 : MBDynInterface(pInterface)
1913 {
1914 
1915 }
1916 
1917 SimulationEntityInterface::~SimulationEntityInterface()
1918 {
1919 
1920 }
1921 
1922 METHOD_DEFINE(SimulationEntityInterface, iGetNumPrivData, args, nargout)
1923 {
1924  const SimulationEntity* const pSimulationEntity = Get();
1925 
1926  if (!pSimulationEntity) {
1927  error("%s: not connected", type_name().c_str());
1928  return octave_value();
1929  }
1930 
1931  if (args.length() != 0) {
1932  error("%s: invalid number of arguments\n"
1933  "expected no arguments", type_name().c_str());
1934 
1935  return octave_value();
1936  }
1937 
1938  return octave_value(octave_int<unsigned int>(pSimulationEntity->iGetNumPrivData()));
1939 }
1940 
1941 METHOD_DEFINE(SimulationEntityInterface, iGetPrivDataIdx, args, nargout)
1942 {
1943  const SimulationEntity* pSimulationEntity = Get();
1944 
1945  if (!pSimulationEntity) {
1946  error("%s: not connected",
1947  type_name().c_str());
1948  return octave_value();
1949  }
1950 
1951  if (args.length() != 1) {
1952  error("%s: invalid number of arguments\n"
1953  "expected one argument",
1954  type_name().c_str());
1955 
1956  return octave_value();
1957  }
1958 
1959  if ( !args(0).is_string() ) {
1960  error("%s: invalid argument type (%s)\n"
1961  "expected string",
1962  type_name().c_str(),
1963  args(0).type_name().c_str());
1964 
1965  return octave_value();
1966  }
1967 
1968  const std::string name = args(0).string_value();
1969 
1970  return octave_value(octave_int<unsigned int>(pSimulationEntity->iGetPrivDataIdx(name.c_str())));
1971 }
1972 
1973 METHOD_DEFINE(SimulationEntityInterface, dGetPrivData, args, nargout)
1974 {
1975  const SimulationEntity* pSimulationEntity = Get();
1976 
1977  if (!pSimulationEntity) {
1978  error("%s: not connected", type_name().c_str());
1979  return octave_value();
1980  }
1981 
1982  if (args.length() != 1) {
1983  error("%s: invalid number of arguments\n"
1984  "expected one arguments", type_name().c_str());
1985 
1986  return octave_value();
1987  }
1988 
1989  if (!(args(0).is_integer_type() && args(0).is_scalar_type())) {
1990  error("%s: invalid argument type (%s)\n"
1991  "expected integer scalar",
1992  type_name().c_str(),
1993  args(0).type_name().c_str());
1994  return octave_value();
1995  }
1996 
1997  const int32_t idx = args(0).int32_scalar_value();
1998 
1999  return octave_value(pSimulationEntity->dGetPrivData(idx));
2000 }
2001 
2002 BEGIN_METHOD_TABLE(SimulationEntityInterface, MBDynInterface)
2003  METHOD_DISPATCH(SimulationEntityInterface, iGetNumPrivData)
2004  METHOD_DISPATCH(SimulationEntityInterface, iGetPrivDataIdx)
2005  METHOD_DISPATCH(SimulationEntityInterface, dGetPrivData)
2006 END_METHOD_TABLE()
2007 
2008 NodeInterface::NodeInterface(OctaveInterface* pInterface)
2009 : SimulationEntityInterface(pInterface)
2010 {
2011 
2012 }
2013 
2014 NodeInterface::~NodeInterface()
2015 {
2016 
2017 }
2018 
2019 void NodeInterface::print(std::ostream& os, bool pr_as_read_syntax) const
2020 {
2021  const Node* const pNode = Get();
2022 
2023  if (!pNode) {
2024  error("%s: not connected", type_name().c_str());
2025  return;
2026  }
2027 
2028  os << type_name() << "(" << pNode->GetLabel() << "):" << pNode->GetName() << std::endl;
2029 }
2030 
2031 METHOD_DEFINE(NodeInterface, GetLabel, args, nargout)
2032 {
2033  const Node* const pNode = Get();
2034 
2035  if (pNode == 0) {
2036  error("%s: not connected", type_name().c_str());
2037  return octave_value();
2038  }
2039 
2040  if (args.length() != 0) {
2041  error("%s: invalid number of arguments %ld\n"
2042  "no arguments expected",
2043  type_name().c_str(),
2044  long(args.length()));
2045 
2046  return octave_value();
2047  }
2048 
2049  return octave_value(octave_int<unsigned int>(pNode->GetLabel()));
2050 }
2051 
2052 METHOD_DEFINE(NodeInterface, GetName, args, nargout)
2053 {
2054  const Node* const pNode = Get();
2055 
2056  if (pNode == 0) {
2057  error("%s: not connected", type_name().c_str());
2058  return octave_value();
2059  }
2060 
2061  if (args.length() != 0) {
2062  error("%s: invalid number of arguments %ld\n"
2063  "no arguments expected",
2064  type_name().c_str(),
2065  long(args.length()));
2066 
2067  return octave_value();
2068  }
2069 
2070  return octave_value(pNode->GetName());
2071 }
2072 
2073 METHOD_DEFINE(NodeInterface, iGetFirstIndex, args, nargout)
2074 {
2075  const Node* const pNode = Get();
2076 
2077  if (pNode == 0) {
2078  error("%s: not connected", type_name().c_str());
2079  return octave_value();
2080  }
2081 
2082  if (args.length() != 0) {
2083  error("%s: invalid number of arguments %ld\n"
2084  "no arguments expected",
2085  type_name().c_str(),
2086  long(args.length()));
2087 
2088  return octave_value();
2089  }
2090 
2091  return octave_value(octave_int<integer>(pNode->iGetFirstIndex()));
2092 }
2093 
2094 METHOD_DEFINE(NodeInterface, iGetFirstRowIndex, args, nargout)
2095 {
2096  const Node* const pNode = Get();
2097 
2098  if (pNode == 0) {
2099  error("%s: not connected", type_name().c_str());
2100  return octave_value();
2101  }
2102 
2103  if (args.length() != 0) {
2104  error("%s: invalid number of arguments %ld\n"
2105  "no arguments expected",
2106  type_name().c_str(),
2107  long(args.length()));
2108 
2109  return octave_value();
2110  }
2111 
2112  return octave_value(octave_int<integer>(pNode->iGetFirstRowIndex()));
2113 }
2114 
2115 METHOD_DEFINE(NodeInterface, iGetFirstColIndex, args, nargout)
2116 {
2117  const Node* const pNode = Get();
2118 
2119  if (pNode == 0) {
2120  error("%s: not connected", type_name().c_str());
2121  return octave_value();
2122  }
2123 
2124  if (args.length() != 0) {
2125  error("%s: invalid number of arguments %ld\n"
2126  "no arguments expected",
2127  type_name().c_str(),
2128  long(args.length()));
2129 
2130  return octave_value();
2131  }
2132 
2133  return octave_value(octave_int<integer>(pNode->iGetFirstColIndex()));
2134 }
2135 
2136 METHOD_DEFINE(NodeInterface, dGetDofValue, args, nargout)
2137 {
2138  const Node* const pNode = Get();
2139 
2140  if (pNode == 0) {
2141  error("%s: not connected", type_name().c_str());
2142  return octave_value();
2143  }
2144 
2145  if (args.length() != 2) {
2146  error("%s: invalid number of arguments %ld\n"
2147  "%ld arguments expected",
2148  type_name().c_str(),
2149  long(args.length()),
2150  2L);
2151 
2152  return octave_value();
2153  }
2154 
2155  if (!(args(0).is_integer_type() && args(0).is_scalar_type())) {
2156  error("%s: invalid argument type (%s) for argument 1\n"
2157  "integer scalar expected",
2158  type_name().c_str(),
2159  args(0).type_name().c_str());
2160 
2161  return octave_value();
2162  }
2163 
2164  const int32_t iDof = args(0).int32_scalar_value();
2165 
2166  if (!(args(1).is_integer_type() && args(1).is_scalar_type())) {
2167  error("%s: invalid argument type (%s) for argument 2\n"
2168  "integer scalar expected",
2169  type_name().c_str(),
2170  args(1).type_name().c_str());
2171 
2172  return octave_value();
2173  }
2174 
2175  const int32_t iOrder = args(1).int32_scalar_value();
2176 
2177  return octave_value(pNode->dGetDofValue(iDof, iOrder));
2178 }
2179 
2180 METHOD_DEFINE(NodeInterface, dGetDofValuePrev, args, nargout)
2181 {
2182  const Node* const pNode = Get();
2183 
2184  if (pNode == 0) {
2185  error("%s: not connected", type_name().c_str());
2186  return octave_value();
2187  }
2188 
2189  if (args.length() != 2) {
2190  error("%s: invalid number of arguments %ld\n"
2191  "%ld arguments expected",
2192  type_name().c_str(),
2193  long(args.length()),
2194  2L);
2195 
2196  return octave_value();
2197  }
2198 
2199  if (!(args(0).is_integer_type() && args(0).is_scalar_type())) {
2200  error("%s: invalid argument type (%s) for argument 1\n"
2201  "integer scalar expected",
2202  type_name().c_str(),
2203  args(0).type_name().c_str());
2204 
2205  return octave_value();
2206  }
2207 
2208  const int32_t iDof = args(0).int32_scalar_value();
2209 
2210  if (!(args(1).is_integer_type() && args(1).is_scalar_type())) {
2211  error("%s: invalid argument type (%s) for argument 2\n"
2212  "integer scalar expected",
2213  type_name().c_str(),
2214  args(1).type_name().c_str());
2215 
2216  return octave_value();
2217  }
2218 
2219  const int32_t iOrder = args(1).int32_scalar_value();
2220 
2221  return octave_value(pNode->dGetDofValuePrev(iDof, iOrder));
2222 }
2223 
2224 BEGIN_METHOD_TABLE(NodeInterface, SimulationEntityInterface)
2225  METHOD_DISPATCH(NodeInterface, GetLabel)
2226  METHOD_DISPATCH(NodeInterface, GetName)
2227  METHOD_DISPATCH(NodeInterface, iGetFirstIndex)
2228  METHOD_DISPATCH(NodeInterface, iGetFirstRowIndex)
2229  METHOD_DISPATCH(NodeInterface, iGetFirstColIndex)
2230  METHOD_DISPATCH(NodeInterface, dGetDofValue)
2231  METHOD_DISPATCH(NodeInterface, dGetDofValuePrev)
2232 END_METHOD_TABLE()
2233 
2234 ScalarNodeInterface::ScalarNodeInterface(OctaveInterface* pInterface, ScalarNode* pNode)
2235 : NodeInterface(pInterface),
2236 pNode(pNode)
2237 {
2238 
2239 }
2240 
2241 ScalarNodeInterface::~ScalarNodeInterface()
2242 {
2243 
2244 }
2245 
2246 const ScalarNode* ScalarNodeInterface::Get()const
2247 {
2248  return pNode;
2249 }
2250 
2251 METHOD_DEFINE(ScalarNodeInterface, SetX, args, nargout)
2252 {
2253  if (!pNode) {
2254  error("%s: not connected", type_name().c_str());
2255  return octave_value();
2256  }
2257 
2258  if (args.length() != 1) {
2259  error("%s: invalid number of arguments %ld\n"
2260  "expected one argument",
2261  type_name().c_str(),
2262  long(args.length()));
2263 
2264  return octave_value();
2265  }
2266 
2267  if (!args(0).is_real_scalar()) {
2268  error("%s: invalid argument type (%s)\n"
2269  "expected real scalar",
2270  type_name().c_str(),
2271  args(0).type_name().c_str());
2272 
2273  return octave_value();
2274  }
2275 
2276  pNode->SetX(args(0).scalar_value());
2277 
2278  return octave_value();
2279 }
2280 
2281 METHOD_DEFINE(ScalarNodeInterface, dGetX, args, nargout)
2282 {
2283  if (!pNode) {
2284  error("%s: not connected", type_name().c_str());
2285  return octave_value();
2286  }
2287 
2288  if (args.length() != 0) {
2289  error("%s: invalid number of arguments %ld\n"
2290  "expected no argument",
2291  type_name().c_str(),
2292  long(args.length()));
2293 
2294  return octave_value();
2295  }
2296 
2297  return octave_value(pNode->dGetX());
2298 }
2299 
2300 METHOD_DEFINE(ScalarNodeInterface, SetXPrime, args, nargout)
2301 {
2302  if (!pNode) {
2303  error("%s: not connected", type_name().c_str());
2304  return octave_value();
2305  }
2306 
2307  if (args.length() != 1) {
2308  error("%s: invalid number of arguments %ld\n"
2309  "expected one argument",
2310  type_name().c_str(),
2311  long(args.length()));
2312 
2313  return octave_value();
2314  }
2315 
2316  if (!args(0).is_real_scalar()) {
2317  error("%s: invalid argument type (%s)\n"
2318  "expected real scalar",
2319  type_name().c_str(),
2320  args(0).type_name().c_str());
2321 
2322  return octave_value();
2323  }
2324 
2325  pNode->SetXPrime(args(0).scalar_value());
2326 
2327  return octave_value();
2328 }
2329 
2330 METHOD_DEFINE(ScalarNodeInterface, dGetXPrime, args, nargout)
2331 {
2332  if (!pNode) {
2333  error("%s: not connected", type_name().c_str());
2334  return octave_value();
2335  }
2336 
2337  if (args.length() != 0) {
2338  error("%s: invalid number of arguments %ld\n"
2339  "expected no argument",
2340  type_name().c_str(),
2341  long(args.length()));
2342 
2343  return octave_value();
2344  }
2345 
2346  return octave_value(pNode->dGetXPrime());
2347 }
2348 
2349 BEGIN_METHOD_TABLE(ScalarNodeInterface, NodeInterface)
2350  METHOD_DISPATCH(ScalarNodeInterface, SetX)
2351  METHOD_DISPATCH(ScalarNodeInterface, dGetX)
2352  METHOD_DISPATCH(ScalarNodeInterface, SetXPrime)
2353  METHOD_DISPATCH(ScalarNodeInterface, dGetXPrime)
2354 END_METHOD_TABLE()
2355 
2356 DEFINE_OCTAVE_ALLOCATOR(ScalarNodeInterface);
2357 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA(ScalarNodeInterface, "ScalarNode", "ScalarNode");
2358 
2359 StructDispNodeBaseInterface::StructDispNodeBaseInterface(OctaveInterface* pInterface)
2360 : NodeInterface(pInterface)
2361 {
2362  TRACE("constructor");
2363 }
2364 
2365 StructDispNodeBaseInterface::~StructDispNodeBaseInterface()
2366 {
2367  TRACE("destructor");
2368 }
2369 
2370 octave_value
2371 StructDispNodeBaseInterface::GetVec3(const Vec3& V, const octave_value_list& args) const
2372 {
2373  octave_value octValue;
2374 
2375  if (args.length() != 0) {
2376  error("StructNode: invalid number of arguments %ld"
2377  "no arguments expected", long(args.length()));
2378  return octValue;
2379  }
2380 
2381  if (!GetInterface()->ConvertMBDynToOctave(V, octValue)) {
2382  error("StructNode: could not convert data");
2383  }
2384 
2385  return octValue;
2386 }
2387 
2388 octave_value
2389 StructDispNodeBaseInterface::GetMat3x3(const Mat3x3& M, const octave_value_list& args) const
2390 {
2391  octave_value octValue;
2392 
2393  if (args.length() != 0) {
2394  error("StructNode: invalid number of arguments %ld"
2395  "no arguments expected", long(args.length()));
2396  return octValue;
2397  }
2398 
2399  if (!GetInterface()->ConvertMBDynToOctave(M, octValue)) {
2400  error("StructNode: could not convert data");
2401  }
2402 
2403  return octValue;
2404 }
2405 
2406 METHOD_DEFINE(StructDispNodeBaseInterface, iGetFirstPositionIndex, args, nargout)
2407 {
2408  const StructDispNode* const pNode = Get();
2409 
2410  if (pNode == 0) {
2411  error("StructNode: not connected");
2412  return octave_value();
2413  }
2414 
2415  if (args.length() != 0) {
2416  error("StructNode: invalid number of arguments %ld"
2417  "no arguments expected", long(args.length()));
2418  return octave_value();
2419  }
2420 
2421  return octave_value(octave_int<integer>(pNode->iGetFirstPositionIndex()));
2422 }
2423 
2424 METHOD_DEFINE(StructDispNodeBaseInterface, iGetFirstMomentumIndex, args, nargout)
2425 {
2426  const StructDispNode* const pNode = Get();
2427 
2428  if (pNode == 0) {
2429  error("StructNode: not connected");
2430  return octave_value();
2431  }
2432 
2433  if (args.length() != 0) {
2434  error("StructNode: invalid number of arguments %ld"
2435  "no arguments expected", long(args.length()));
2436  return octave_value();
2437  }
2438 
2439  return octave_value(octave_int<integer>(pNode->iGetFirstMomentumIndex()));
2440 }
2441 
2442 METHOD_DEFINE(StructDispNodeBaseInterface, GetLabel, args, nargout)
2443 {
2444  const StructDispNode* const pNode = Get();
2445 
2446  if (pNode == 0) {
2447  error("StructNode: not connected");
2448  return octave_value();
2449  }
2450 
2451  if (args.length() != 0) {
2452  error("StructNode: invalid number of arguments %ld\n"
2453  "no arguments expected", long(args.length()));
2454  return octave_value();
2455  }
2456 
2457  return octave_value(octave_int<integer>(pNode->GetLabel()));
2458 }
2459 
2460 METHOD_DEFINE(StructDispNodeBaseInterface, GetXCurr, args, nargout)
2461 {
2462  const StructDispNode* const pNode = Get();
2463 
2464  if (pNode == 0) {
2465  error("StructNode: not connected");
2466  return octave_value();
2467  }
2468 
2469  return GetVec3(pNode->GetXCurr(), args);
2470 }
2471 
2472 METHOD_DEFINE(StructDispNodeBaseInterface, GetXPrev, args, nargout)
2473 {
2474  const StructDispNode* const pNode = Get();
2475 
2476  if (pNode == 0) {
2477  error("StructNode: not connected");
2478  return octave_value();
2479  }
2480 
2481  return GetVec3(pNode->GetXPrev(), args);
2482 }
2483 
2484 METHOD_DEFINE(StructDispNodeBaseInterface, GetVCurr, args, nargout)
2485 {
2486  const StructDispNode* const pNode = Get();
2487 
2488  if (pNode == 0) {
2489  error("StructNode: not connected");
2490  return octave_value();
2491  }
2492 
2493  return GetVec3(pNode->GetVCurr(), args);
2494 }
2495 
2496 METHOD_DEFINE(StructDispNodeBaseInterface, GetVPrev, args, nargout)
2497 {
2498  const StructDispNode* const pNode = Get();
2499 
2500  if (pNode == 0) {
2501  error("StructNode: not connected");
2502  return octave_value();
2503  }
2504 
2505  return GetVec3(pNode->GetVPrev(), args);
2506 }
2507 
2508 METHOD_DEFINE(StructDispNodeBaseInterface, GetXPPCurr, args, nargout)
2509 {
2510  const StructDispNode* const pNode = Get();
2511 
2512  if (pNode == 0) {
2513  error("StructNode: not connected");
2514  return octave_value();
2515  }
2516 
2517  if (!pNode->bComputeAccelerations()) {
2518  error("StructNode: accelerations are not available for node %u", pNode->GetLabel());
2519  return octave_value();
2520  }
2521 
2522  return GetVec3(pNode->GetXPPCurr(), args);
2523 }
2524 
2525 METHOD_DEFINE(StructDispNodeBaseInterface, GetXPPPrev, args, nargout)
2526 {
2527  const StructDispNode* const pNode = Get();
2528 
2529  if (pNode == 0) {
2530  error("StructNode: not connected");
2531  return octave_value();
2532  }
2533 
2534  if (!pNode->bComputeAccelerations()) {
2535  error("StructNode: accelerations are not available for node %u", pNode->GetLabel());
2536  return octave_value();
2537  }
2538 
2539  return GetVec3(pNode->GetXPPPrev(), args);
2540 }
2541 
2542 BEGIN_METHOD_TABLE(StructDispNodeBaseInterface, NodeInterface)
2543  METHOD_DISPATCH(StructDispNodeBaseInterface, iGetFirstIndex)
2544  METHOD_DISPATCH(StructDispNodeBaseInterface, iGetFirstPositionIndex)
2545  METHOD_DISPATCH(StructDispNodeBaseInterface, iGetFirstMomentumIndex)
2546  METHOD_DISPATCH(StructDispNodeBaseInterface, GetLabel)
2547  METHOD_DISPATCH(StructDispNodeBaseInterface, GetXCurr)
2548  METHOD_DISPATCH(StructDispNodeBaseInterface, GetXPrev)
2549  METHOD_DISPATCH(StructNodeInterface, GetVCurr)
2550  METHOD_DISPATCH(StructNodeInterface, GetVPrev)
2551  METHOD_DISPATCH(StructNodeInterface, GetXPPCurr)
2552  METHOD_DISPATCH(StructNodeInterface, GetXPPPrev)
2553 END_METHOD_TABLE()
2554 
2555 StructDispNodeInterface::StructDispNodeInterface(OctaveInterface* pInterface, const StructDispNode* pNode)
2556 : StructDispNodeBaseInterface(pInterface),
2557  pNode(pNode)
2558 {
2559 
2560 }
2561 
2562 StructDispNodeInterface::~StructDispNodeInterface()
2563 {
2564 
2565 }
2566 
2567 const StructDispNode* StructDispNodeInterface::Get()const
2568 {
2569  return pNode;
2570 }
2571 
2572 DEFINE_OCTAVE_ALLOCATOR(StructDispNodeInterface);
2573 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA(StructDispNodeInterface, "StructDispNode", "StructDispNode");
2574 
2575 StructNodeInterface::StructNodeInterface(OctaveInterface* pInterface, const StructNode* pNode)
2576 : StructDispNodeBaseInterface(pInterface),
2577 pNode(pNode)
2578 {
2579  TRACE("constructor");
2580 }
2581 
2582 StructNodeInterface::~StructNodeInterface(void)
2583 {
2584  TRACE("destructor");
2585 }
2586 
2587 const StructNode*
2588 StructNodeInterface::Get()const
2589 {
2590  return pNode;
2591 }
2592 
2593 METHOD_DEFINE(StructNodeInterface, GetgCurr, args, nargout)
2594 {
2595  if (pNode == 0) {
2596  error("StructNode: not connected");
2597  return octave_value();
2598  }
2599 
2600  return GetVec3(pNode->GetgCurr(), args);
2601 }
2602 
2603 METHOD_DEFINE(StructNodeInterface, GetgRef, args, nargout)
2604 {
2605  if (pNode == 0) {
2606  error("StructNode: not connected");
2607  return octave_value();
2608  }
2609 
2610  return GetVec3(pNode->GetgRef(), args);
2611 }
2612 
2613 METHOD_DEFINE(StructNodeInterface, GetgPCurr, args, nargout)
2614 {
2615  if (pNode == 0) {
2616  error("StructNode: not connected");
2617  return octave_value();
2618  }
2619 
2620  return GetVec3(pNode->GetgPCurr(), args);
2621 }
2622 
2623 METHOD_DEFINE(StructNodeInterface, GetgPRef, args, nargout)
2624 {
2625  if (pNode == 0) {
2626  error("StructNode: not connected");
2627  return octave_value();
2628  }
2629 
2630  return GetVec3(pNode->GetgPRef(), args);
2631 }
2632 
2633 METHOD_DEFINE(StructNodeInterface, GetRCurr, args, nargout)
2634 {
2635  if (pNode == 0) {
2636  error("StructNode: not connected");
2637  return octave_value();
2638  }
2639 
2640  return GetMat3x3(pNode->GetRCurr(), args);
2641 }
2642 
2643 METHOD_DEFINE(StructNodeInterface, GetRPrev, args, nargout)
2644 {
2645  if (pNode == 0) {
2646  error("StructNode: not connected");
2647  return octave_value();
2648  }
2649 
2650  return GetMat3x3(pNode->GetRPrev(), args);
2651 }
2652 
2653 METHOD_DEFINE(StructNodeInterface, GetRRef, args, nargout)
2654 {
2655  if (pNode == 0) {
2656  error("StructNode: not connected");
2657  return octave_value();
2658  }
2659 
2660  return GetMat3x3(pNode->GetRRef(), args);
2661 }
2662 
2663 METHOD_DEFINE(StructNodeInterface, GetWCurr, args, nargout)
2664 {
2665  if (pNode == 0) {
2666  error("StructNode: not connected");
2667  return octave_value();
2668  }
2669 
2670  return GetVec3(pNode->GetWCurr(), args);
2671 }
2672 
2673 METHOD_DEFINE(StructNodeInterface, GetWPrev, args, nargout)
2674 {
2675  if (pNode == 0) {
2676  error("StructNode: not connected");
2677  return octave_value();
2678  }
2679 
2680  return GetVec3(pNode->GetWPrev(), args);
2681 }
2682 
2683 METHOD_DEFINE(StructNodeInterface, GetWRef, args, nargout)
2684 {
2685  if (pNode == 0) {
2686  error("StructNode: not connected");
2687  return octave_value();
2688  }
2689 
2690  return GetVec3(pNode->GetWRef(), args);
2691 }
2692 
2693 METHOD_DEFINE(StructNodeInterface, GetWPCurr, args, nargout)
2694 {
2695  if (pNode == 0) {
2696  error("StructNode: not connected");
2697  return octave_value();
2698  }
2699 
2700  if (!pNode->bComputeAccelerations()) {
2701  error("StructNode: accelerations are not available for node %u", pNode->GetLabel());
2702  return octave_value();
2703  }
2704 
2705  return GetVec3(pNode->GetWPCurr(), args);
2706 }
2707 
2708 METHOD_DEFINE(StructNodeInterface, GetWPPrev, args, nargout)
2709 {
2710  if (pNode == 0) {
2711  error("StructNode: not connected");
2712  return octave_value();
2713  }
2714 
2715  if (!pNode->bComputeAccelerations()) {
2716  error("StructNode: accelerations are not available for node %u", pNode->GetLabel());
2717  return octave_value();
2718  }
2719 
2720  return GetVec3(pNode->GetWPPrev(), args);
2721 }
2722 
2723 BEGIN_METHOD_TABLE(StructNodeInterface, StructDispNodeBaseInterface)
2724  METHOD_DISPATCH(StructNodeInterface, GetgCurr)
2725  METHOD_DISPATCH(StructNodeInterface, GetgRef)
2726  METHOD_DISPATCH(StructNodeInterface, GetgPCurr)
2727  METHOD_DISPATCH(StructNodeInterface, GetgPRef)
2728  METHOD_DISPATCH(StructNodeInterface, GetRCurr)
2729  METHOD_DISPATCH(StructNodeInterface, GetRPrev)
2730  METHOD_DISPATCH(StructNodeInterface, GetRRef)
2731  METHOD_DISPATCH(StructNodeInterface, GetWCurr)
2732  METHOD_DISPATCH(StructNodeInterface, GetWPrev)
2733  METHOD_DISPATCH(StructNodeInterface, GetWRef)
2734  METHOD_DISPATCH(StructNodeInterface, GetWPCurr)
2735  METHOD_DISPATCH(StructNodeInterface, GetWPPrev)
2736 END_METHOD_TABLE()
2737 
2738 DEFINE_OCTAVE_ALLOCATOR(StructNodeInterface);
2739 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA(StructNodeInterface, "StructNode", "StructNode");
2740 
2741 DataManagerInterface::DataManagerInterface(OctaveInterface* pInterface, const DataManager* pDM)
2742  :MBDynInterface(pInterface),
2743  pDM(pDM)
2744 {
2745  TRACE("constructor");
2746 
2747  ASSERT(pInterface != 0);
2748  ASSERT(pDM != 0);
2749 }
2750 
2751 DataManagerInterface::~DataManagerInterface()
2752 {
2753  TRACE("destructor");
2754 }
2755 
2756 void DataManagerInterface::print(std::ostream& os, bool pr_as_read_syntax)const
2757 {
2758  os << "DataManager" << std::endl
2759  << "iTotDofs=" << GetDataManager()->iGetNumDofs() << std::endl;
2760 }
2761 
2762 METHOD_DEFINE(DataManagerInterface, GetVariable, args, nargout)
2763 {
2764  octave_value octValue;
2765 
2766  if (args.length() != 1) {
2767  error("DataManager: invalid number of arguments %ld\n"
2768  "one argument expected", long(args.length()));
2769  return octValue;
2770  }
2771 
2772  if (!args(0).is_string()) {
2773  error("DataManager: invalid argument type \"%s\"\n"
2774  "varName must be a string", args(0).type_name().c_str());
2775  return octValue;
2776  }
2777 
2778  const std::string varName = args(0).string_value();
2779 
2780  const NamedValue* const pVar = GetSymbolTable().Get(varName.c_str());
2781 
2782  if (pVar == 0) {
2783  error("DataManager: variable \"%s\" not found in symbol table", varName.c_str());
2784  return octValue;
2785  }
2786 
2787  const TypedValue mbValue = pVar->GetVal();
2788 
2789  if (!GetInterface()->ConvertMBDynToOctave(mbValue, octValue)) {
2790  error("%s: failed to convert MBDyn variable %s"
2791  " of type %s value into octave value",
2792  type_name().c_str(),
2793  varName.c_str(),
2794  mbValue.GetTypeName());
2795  }
2796 
2797  return octValue;
2798 }
2799 
2800 METHOD_DEFINE(DataManagerInterface, GetStructNodePos, args, nargout)
2801 {
2802  octave_value octValue;
2803 
2804  if (args.length() != 1) {
2805  error("DataManager: invalid number of arguments\n"
2806  "one argument expected");
2807  return octValue;
2808  }
2809 
2810  if (!(args(0).is_scalar_type() && args(0).is_integer_type())) {
2811  error("DataManager: invalid argument type \"%s\"\n"
2812  "iNodeLabel must be an integer", args(0).type_name().c_str());
2813  return octValue;
2814  }
2815 
2816  integer iNodeLabel = args(0).int32_scalar_value();
2817 
2818  Node* pNode = GetDataManager()->pFindNode(Node::STRUCTURAL, iNodeLabel);
2819 
2820  if (pNode == 0) {
2821  error("DataManager: invalid node label\n"
2822  "could not find node %ld", long(iNodeLabel));
2823  return octValue;
2824  }
2825 
2826  StructDispNode* pStructNode = dynamic_cast<StructDispNode*>(pNode);
2827 
2828  if (pStructNode == 0) {
2829  error("DataManager: invalid node type\n"
2830  "node %ld is not a structural displacement node", long(iNodeLabel));
2831  return octValue;
2832  }
2833 
2834  const Vec3& X = pStructNode->GetXCurr();
2835 
2836  ColumnVector octX(3);
2837 
2838  for (int i = 0; i < 3; ++i) {
2839  octX(i) = X(i + 1);
2840  }
2841 
2842  octValue = octX;
2843 
2844  return octValue;
2845 }
2846 
2847 METHOD_DEFINE(DataManagerInterface, GetStructNode, args, nargout)
2848 {
2849  octave_value octValue;
2850 
2851  if (args.length() != 1) {
2852  error("DataManager: invalid number of arguments\n"
2853  "one argument expected");
2854  return octValue;
2855  }
2856 
2857  if (!(args(0).is_scalar_type() && args(0).is_integer_type())) {
2858  error("DataManager: invalid argument type \"%s\"\n"
2859  "iNodeLabel must be an integer", args(0).type_name().c_str());
2860  return octValue;
2861  }
2862 
2863  integer iNodeLabel = static_cast<int32_t>(args(0).int32_scalar_value());
2864 
2865  Node* pNode = GetDataManager()->pFindNode(Node::STRUCTURAL, iNodeLabel);
2866 
2867  if (pNode == 0) {
2868  error("DataManager: could not find node %ld", long(iNodeLabel));
2869  return octValue;
2870  }
2871 
2872  StructNode* pStructNode = dynamic_cast<StructNode*>(pNode);
2873 
2874  if (pStructNode) {
2875  octValue = new StructNodeInterface(GetInterface(), pStructNode);
2876  return octValue;
2877  }
2878 
2879  StructDispNode* pStructDispNode = dynamic_cast<StructDispNode*>(pNode);
2880 
2881  if (!pStructDispNode) {
2882  error("%s: node %ld is not a structural node", type_name().c_str(), long(iNodeLabel));
2883  return octValue;
2884  }
2885 
2886  octValue = new StructDispNodeInterface(GetInterface(), pStructDispNode);
2887 
2888  return octValue;
2889 }
2890 
2891 METHOD_DEFINE(DataManagerInterface, pFindNode, args, nargout)
2892 {
2893  Node::Type Type;
2894 
2895  Node* const pNode = pFindNode_(args, Type);
2896 
2897  if (!pNode) {
2898  return octave_value(); // error has been called
2899  }
2900 
2901  NodeInterface* const pNodeInterface = CreateNodeInterface(pNode, Type);
2902 
2903  if (pNodeInterface == 0) {
2904  error("%s: unknown node type %ld", type_name().c_str(), long(Type));
2905  return octave_value();
2906  }
2907 
2908  return octave_value(pNodeInterface);
2909 }
2910 
2911 METHOD_DEFINE(DataManagerInterface, ReadNode, args, nargout)
2912 {
2913  if (args.length() != 2) {
2914  error("%s: invalid number of arguments (%ld)\n"
2915  "arguments expected: %ld",
2916  type_name().c_str(),
2917  long(args.length()),
2918  2L);
2919 
2920  return octave_value();
2921  }
2922 
2923  const MBDynParserInterface* const pHPI =
2924  dynamic_cast<const MBDynParserInterface*>(&args(0).get_rep());
2925 
2926  if (pHPI == 0) {
2927  error("%s: invalid argument type (%s)\n"
2928  "expected MBDynParser",
2929  type_name().c_str(),
2930  args(0).type_name().c_str());
2931 
2932  return octave_value();
2933  }
2934 
2935  MBDynParser* const pHP = pHPI->Get();
2936 
2937  if (pHP == 0) {
2938  error("%s: not connected", pHPI->type_name().c_str());
2939  return octave_value();
2940  }
2941 
2942  if (!args(1).is_string()) {
2943  error("%s: invalid argument type (%s)\n"
2944  "expected string",
2945  type_name().c_str(),
2946  args(1).type_name().c_str());
2947 
2948  return octave_value();
2949  }
2950 
2951  const std::string strType(args(1).string_value());
2952 
2953  const Node::Type Type = GetNodeType(strType);
2954 
2955  if (Type == Node::UNKNOWN) {
2956  error("%s: invalid node type (%s)", type_name().c_str(), strType.c_str());
2957  return octave_value();
2958  }
2959 
2960  Node* pNode = 0;
2961 
2962  try {
2963  // FIXME: This function should not be called from a drive caller
2964  // FIXME: because DataManager is assumed to be constant in this case?
2965  pNode = const_cast<DataManager*>(pDM)->ReadNode(*pHP, Type);
2966  } catch(...) { // do not throw across call back boundaries
2967  error("%s: ReadNode failed", type_name().c_str());
2968  return octave_value();
2969  }
2970 
2971  NodeInterface* const pNodeInterface = CreateNodeInterface(pNode, Type);
2972 
2973  if (pNodeInterface == 0) {
2974  error("%s: unknown node type %ld", type_name().c_str(), long(Type));
2975  return octave_value();
2976  }
2977 
2978  return octave_value(pNodeInterface);
2979 }
2980 
2981 METHOD_DEFINE(DataManagerInterface, dGetTime, args, nargout)
2982 {
2983  if (args.length() != 0) {
2984  error("DataManager: invalid number of arguments\n"
2985  "no argument expected");
2986  return octave_value();
2987  }
2988 
2989  return octave_value(GetDataManager()->dGetTime());
2990 }
2991 
2992 const DataManager* DataManagerInterface::GetDataManager()const
2993 {
2994  ASSERT(GetInterface() != 0);
2995  ASSERT(pDM != 0);
2996 
2997  return pDM;
2998 }
2999 
3000 const Table& DataManagerInterface::GetSymbolTable()const
3001 {
3002  return GetDataManager()->GetMathParser().GetSymbolTable();
3003 }
3004 
3005 Node* DataManagerInterface::pFindNode_(const octave_value_list& args, Node::Type& Type) const
3006 {
3007  if (args.length() != 2) {
3008  error("%s: invalid number of arguments\n"
3009  "arguments expected 2", type_name().c_str());
3010 
3011  return 0;
3012  }
3013 
3014  if (!(args(0).is_string())) {
3015  error("%s: invalid argument type (%s)\n"
3016  "Type must be a string",
3017  type_name().c_str(),
3018  args(0).type_name().c_str());
3019 
3020  return 0;
3021  }
3022 
3023  const std::string strType(args(0).string_value());
3024 
3025  Type = GetNodeType(strType);
3026 
3027  if (Type == Node::UNKNOWN) {
3028  error("%s: unknown node type (%s)", type_name().c_str(), strType.c_str());
3029  return 0;
3030  }
3031 
3032  if (!(args(1).is_scalar_type() && args(1).is_integer_type())) {
3033  error("%s: invalid argument type (%s)\n"
3034  "iNodeLabel must be an integer",
3035  type_name().c_str(),
3036  args(1).type_name().c_str());
3037 
3038  return 0;
3039  }
3040 
3041  const int32_t iNodeLabel = args(1).int32_scalar_value();
3042 
3043  Node* pNode = GetDataManager()->pFindNode(Type, iNodeLabel);
3044 
3045  if (pNode == 0) {
3046  error("%s: could not find node %ld", type_name().c_str(), long(iNodeLabel));
3047  return 0;
3048  }
3049 
3050  return pNode;
3051 }
3052 
3053 Node::Type DataManagerInterface::GetNodeType(const std::string& strType) const
3054 {
3055  static const struct {
3056  char name[11];
3057  Node::Type type;
3058  }
3059  nodeTypes[] = {
3060  { "ABSTRACT", Node::ABSTRACT },
3061  { "STRUCTURAL", Node::STRUCTURAL },
3062  { "ELECTRIC", Node::ELECTRIC },
3063  { "THERMAL", Node::THERMAL },
3064  { "PARAMETER", Node::PARAMETER },
3065  { "HYDRAULIC", Node::HYDRAULIC }
3066  };
3067 
3068  static const int count = sizeof(nodeTypes) / sizeof(nodeTypes[0]);
3069 
3070  for (int i = 0; i < count; ++i) {
3071  if (strType == nodeTypes[i].name) {
3072  return nodeTypes[i].type;
3073  }
3074  }
3075 
3076  return Node::UNKNOWN;
3077 }
3078 
3079 NodeInterface* DataManagerInterface::CreateNodeInterface(Node* pNode, Node::Type Type) const
3080 {
3081  switch (Type) {
3082  case Node::STRUCTURAL:
3083  {
3084  StructNode* const pStructNode = dynamic_cast<StructNode*>(pNode);
3085 
3086  if (pStructNode == 0) {
3087  StructDispNode* pStructDispNode = dynamic_cast<StructDispNode*>(pNode);
3088  if (pStructDispNode == 0) {
3089  ASSERT(0);
3090  return 0;
3091  }
3092 
3093  return new StructDispNodeInterface(GetInterface(), pStructDispNode);
3094  }
3095  return new StructNodeInterface(GetInterface(), pStructNode);
3096  }
3097  case Node::ABSTRACT:
3098  case Node::ELECTRIC:
3099  case Node::THERMAL:
3100  case Node::PARAMETER:
3101  case Node::HYDRAULIC:
3102  {
3103  ScalarNode* const pScalarNode = dynamic_cast<ScalarNode*>(pNode);
3104  ASSERT(pScalarNode != 0);
3105  return new ScalarNodeInterface(GetInterface(), pScalarNode);
3106  }
3107  default:
3108  ASSERT(0); // FIXME: add StructuralDisplacementNode
3109  return 0;
3110  }
3111 }
3112 
3113 BEGIN_METHOD_TABLE(DataManagerInterface, MBDynInterface)
3114  METHOD_DISPATCH(DataManagerInterface, GetVariable)
3115  METHOD_DISPATCH(DataManagerInterface, GetStructNodePos)
3116  METHOD_DISPATCH(DataManagerInterface, GetStructNode)
3117  METHOD_DISPATCH(DataManagerInterface, pFindNode)
3118  METHOD_DISPATCH(DataManagerInterface, ReadNode)
3119  METHOD_DISPATCH(DataManagerInterface, dGetTime)
3120 END_METHOD_TABLE()
3121 
3122 DEFINE_OCTAVE_ALLOCATOR(DataManagerInterface);
3123 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA(DataManagerInterface, "DataManager", "DataManager");
3124 
3125 const MBDynParserInterface::MBDynStringDelims
3126 MBDynParserInterface::mbStringDelims[6] = {
3127  { HighParser::PLAINBRACKETS, "PLAINBRACKETS" },
3128  { HighParser::SQUAREBRACKETS, "SQUAREBRACKETS" },
3129  { HighParser::CURLYBRACKETS, "CURLYBRACKETS" },
3130  { HighParser::SINGLEQUOTE, "SINGLEQUOTE" },
3131  { HighParser::DOUBLEQUOTE, "DOUBLEQUOTE" },
3132  { HighParser::DEFAULTDELIM, "DEFAULTDELIM" }
3133 };
3134 
3135 MBDynParserInterface::MBDynParserInterface(OctaveInterface* pInterface, MBDynParser* pHP)
3136  :MBDynInterface(pInterface),
3137  pHP(pHP)
3138 {
3139  if (!pHP) {
3140  if (pInterface) {
3141  ASSERT(0); // should not happen
3142  this->pHP = &pInterface->GetDataManager()->GetMBDynParser();
3143  }
3144  }
3145 }
3146 
3147 MBDynParserInterface::~MBDynParserInterface()
3148 {
3149 
3150 }
3151 
3152 void MBDynParserInterface::print(std::ostream& os, bool pr_as_read_syntax) const
3153 {
3154  os << "MBDynParser";
3155 
3156  if (pHP) {
3157  const IncludeParser::ErrOut lineData(pHP->GetLineData());
3158  os << " at " << lineData.sFileName << ":" << lineData.iLineNumber;
3159  }
3160 
3161  os << std::endl;
3162 }
3163 
3164 METHOD_DEFINE(MBDynParserInterface, IsKeyWord, args, nargout)
3165 {
3166  if (!pHP) {
3167  error("MBDynParser: not connected");
3168  return octave_value();
3169  }
3170 
3171  if (args.length() != 1) {
3172  error("MBDynParser: invalid number of arguments %ld\n"
3173  "expected %ld arguments", long(args.length()), 1L);
3174  return octave_value();
3175  }
3176 
3177  if (!args(0).is_string()) {
3178  error("MBDynParser: invalid argument type (%s)\n"
3179  "expected string", args(0).type_name().c_str());
3180  return octave_value();
3181  }
3182 
3183  const std::string strKeyWord = args(0).string_value();
3184 
3185  bool bVal;
3186 
3187  try {
3188  bVal = pHP->IsKeyWord(strKeyWord.c_str());
3189  } catch(...) {
3190  error("%s: IsKeyWord failed", type_name().c_str());
3191  return octave_value();
3192  }
3193 
3194  return octave_value(bVal);
3195 }
3196 
3197 METHOD_DEFINE(MBDynParserInterface, IsArg, args, nargout)
3198 {
3199  if (!pHP) {
3200  error("%s: not connected", type_name().c_str());
3201  return octave_value();
3202  }
3203 
3204  if (args.length() != 0) {
3205  error("%s: invalid number of arguments %ld\n"
3206  "no arguments expected",
3207  type_name().c_str(),
3208  long(args.length()));
3209  return octave_value();
3210  }
3211 
3212  return octave_value(pHP->IsArg());
3213 }
3214 
3215 METHOD_DEFINE(MBDynParserInterface, IsStringWithDelims, args, nargout)
3216 {
3217  if (!pHP) {
3218  error("%s: not connected", type_name().c_str());
3219  return octave_value();
3220  }
3221 
3222  if (args.length() > 1) {
3223  error("%s: invalid number of arguments %ld\n"
3224  "one argument expected",
3225  type_name().c_str(),
3226  long(args.length()));
3227  return octave_value();
3228  }
3229 
3231 
3232  if (args.length() >= 1) {
3233  if (!args(0).is_string()) {
3234  error("%s: invalid argument type %s\n"
3235  "expected string",
3236  type_name().c_str(),
3237  args(0).type_name().c_str());
3238  return octave_value();
3239  }
3240 
3241  const std::string strDelims(args(0).string_value());
3242 
3243  if (!GetDelimsFromString(strDelims, delims)) {
3244  error("%s: invalid argument %s",
3245  type_name().c_str(),
3246  strDelims.c_str());
3247  return octave_value();
3248  }
3249  }
3250 
3251  return octave_value(pHP->IsStringWithDelims(delims));
3252 }
3253 
3254 METHOD_DEFINE(MBDynParserInterface, GetReal, args, nargout)
3255 {
3256  if (!pHP) {
3257  error("MBDynParser: not connected");
3258  return octave_value();
3259  }
3260 
3261  doublereal dDefVal = 0.0;
3262 
3263  if (args.length() > 1) {
3264  error("MBDynParser: invalid number of arguments %ld\n"
3265  "expected 0-1 arguments", long(args.length()));
3266  return octave_value();
3267  } else if (args.length() >= 1) {
3268  if (!args(0).is_real_scalar()) {
3269  error("MBDynParser: invalid argument type (%s)\n"
3270  "expected real scalar", args(0).type_name().c_str());
3271  return octave_value();
3272  }
3273  dDefVal = args(0).scalar_value();
3274  }
3275 
3276  doublereal dVal;
3277 
3278  try {
3279  dVal = pHP->GetReal(dDefVal);
3280  } catch(...) {
3281  error("%s: GetReal failed", type_name().c_str());
3282  return octave_value();
3283  }
3284 
3285  return octave_value(dVal);
3286 }
3287 
3288 METHOD_DEFINE(MBDynParserInterface, GetInt, args, nargout)
3289 {
3290  if (!pHP) {
3291  error("%s: not connected", type_name().c_str());
3292  return octave_value();
3293  }
3294 
3295  integer iDefVal = 0;
3296 
3297  if (args.length() > 1) {
3298  error("%s: invalid number of arguments %ld\n"
3299  "expected 0-1 arguments",
3300  type_name().c_str(),
3301  long(args.length()));
3302  return octave_value();
3303  } else if (args.length() >= 1) {
3304  if (!args(0).is_integer_type() && args(0).is_scalar_type()) {
3305  error("%s: invalid argument type (%s)\n"
3306  "expected integer scalar",
3307  type_name().c_str(),
3308  args(0).type_name().c_str());
3309  return octave_value();
3310  }
3311  iDefVal = args(0).int_value();
3312  }
3313 
3314  integer iVal;
3315 
3316  try {
3317  iVal = pHP->GetInt(iDefVal);
3318  } catch(...) {
3319  error("%s: GetInt failed", type_name().c_str());
3320  return octave_value();
3321  }
3322 
3323  return octave_value(iVal);
3324 }
3325 
3326 METHOD_DEFINE(MBDynParserInterface, GetBool, args, nargout)
3327 {
3328  if (!pHP) {
3329  error("%s: not connected", type_name().c_str());
3330  return octave_value();
3331  }
3332 
3333  bool bDefVal = 0;
3334 
3335  if (args.length() > 1) {
3336  error("%s: invalid number of arguments %ld\n"
3337  "expected 0-1 arguments",
3338  type_name().c_str(),
3339  long(args.length()));
3340  return octave_value();
3341  } else if (args.length() >= 1) {
3342  if (!args(0).is_bool_scalar()) {
3343  error("%s: invalid argument type (%s)\n"
3344  "expected bool scalar",
3345  type_name().c_str(),
3346  args(0).type_name().c_str());
3347  return octave_value();
3348  }
3349  bDefVal = args(0).bool_value();
3350  }
3351 
3352  bool bVal;
3353 
3354  try {
3355  bVal = pHP->GetBool(bDefVal);
3356  } catch(...) {
3357  error("%s: GetBool failed", type_name().c_str());
3358  return octave_value();
3359  }
3360 
3361  return octave_value(bVal);
3362 }
3363 
3364 METHOD_DEFINE(MBDynParserInterface, GetYesNoOrBool, args, nargout)
3365 {
3366  if (!pHP) {
3367  error("%s: not connected", type_name().c_str());
3368  return octave_value();
3369  }
3370 
3371  bool bDefVal = 0;
3372 
3373  if (args.length() > 1) {
3374  error("%s: invalid number of arguments %ld\n"
3375  "expected 0-1 arguments",
3376  type_name().c_str(),
3377  long(args.length()));
3378  return octave_value();
3379  } else if (args.length() >= 1) {
3380  if (!args(0).is_bool_scalar()) {
3381  error("%s: invalid argument type (%s)\n"
3382  "expected bool scalar",
3383  type_name().c_str(),
3384  args(0).type_name().c_str());
3385  return octave_value();
3386  }
3387  bDefVal = args(0).bool_value();
3388  }
3389 
3390  bool bVal;
3391 
3392  try {
3393  bVal = pHP->GetYesNoOrBool(bDefVal);
3394  } catch(...) {
3395  error("%s: GetYesNoOrBool failed", type_name().c_str());
3396  return octave_value();
3397  }
3398 
3399  return octave_value(bVal);
3400 }
3401 
3402 METHOD_DEFINE(MBDynParserInterface, GetString, args, nargout)
3403 {
3404  if (!pHP) {
3405  error("%s: not connected", type_name().c_str());
3406  return octave_value();
3407  }
3408 
3409  std::string strDefVal;
3410 
3411  if (args.length() > 1) {
3412  error("%s: invalid number of arguments %ld\n"
3413  "expected 0-1 arguments",
3414  type_name().c_str(),
3415  long(args.length()));
3416  return octave_value();
3417  } else if (args.length() >= 1) {
3418  if (!args(0).is_string()) {
3419  error("%s: invalid argument type (%s)\n"
3420  "expected string",
3421  type_name().c_str(),
3422  args(0).type_name().c_str());
3423  return octave_value();
3424  }
3425  strDefVal = args(0).string_value();
3426  }
3427 
3428  std::string strVal;
3429 
3430  try {
3431  strVal = pHP->GetString(strDefVal);
3432  } catch(...) {
3433  error("%s: GetString failed", type_name().c_str());
3434  return octave_value();
3435  }
3436 
3437  return octave_value(strVal);
3438 }
3439 
3440 METHOD_DEFINE(MBDynParserInterface, GetStringWithDelims, args, nargout)
3441 {
3442  if (!pHP) {
3443  error("%s: not connected", type_name().c_str());
3444  return octave_value();
3445  }
3446 
3447  std::string strDelims;
3448 
3449  if (args.length() > 2) {
3450  error("%s: invalid number of arguments %ld\n"
3451  "expected 0-2 arguments",
3452  type_name().c_str(),
3453  long(args.length()));
3454  return octave_value();
3455  }
3456 
3458 
3459  if (args.length() >= 1) {
3460  if (!args(0).is_string()) {
3461  error("%s: invalid argument type (%s)\n"
3462  "expected string",
3463  type_name().c_str(),
3464  args(0).type_name().c_str());
3465  return octave_value();
3466  }
3467 
3468  std::string strDelims(args(0).string_value());
3469 
3470  if (!GetDelimsFromString(strDelims, delims)) {
3471  error("%s: invalid argument %s",
3472  type_name().c_str(),
3473  strDelims.c_str());
3474  return octave_value();
3475  }
3476  }
3477 
3478  bool bEscape = true;
3479 
3480  if (args.length() >= 2) {
3481  if (!args(1).is_bool_scalar()) {
3482  error("%s: invalid argument type (%s)\n"
3483  "expected bool",
3484  type_name().c_str(),
3485  args(1).type_name().c_str());
3486  return octave_value();
3487  }
3488 
3489  bEscape = args(1).bool_value();
3490  }
3491 
3492  std::string strVal;
3493 
3494  try {
3495  strVal = pHP->GetStringWithDelims(delims, bEscape);
3496  } catch(...) {
3497  error("%s: GetStringWithDelims failed", type_name().c_str());
3498  return octave_value();
3499  }
3500 
3501  return octave_value(strVal);
3502 }
3503 
3504 METHOD_DEFINE(MBDynParserInterface, GetPosRel, args, nargout)
3505 {
3506  if (!pHP) {
3507  error("%s: not connected", type_name().c_str());
3508  return octave_value();
3509  }
3510 
3511  if (args.length() != 1) {
3512  error("%s: invalid number of arguments\n"
3513  "expected one argument", type_name().c_str());
3514 
3515  return octave_value();
3516  }
3517 
3518  const StructDispNode* const pNode = GetStructNode(args(0));
3519 
3520  if (pNode == 0) {
3521  return octave_value();
3522  }
3523 
3524  Vec3 mbX;
3525 
3526  try {
3527  mbX = pHP->GetPosRel(ReferenceFrame(pNode));
3528  } catch(...) { // do not throw across call back boundaries
3529  error("%s: GetPosRel failed", type_name().c_str());
3530  return octave_value();
3531  }
3532 
3533  octave_value octX;
3534 
3535  if (!GetInterface()->ConvertMBDynToOctave(mbX, octX)) {
3536  error("%s: could not convert data", type_name().c_str());
3537  return octave_value();
3538  }
3539 
3540  return octX;
3541 }
3542 
3543 METHOD_DEFINE(MBDynParserInterface, GetRotRel, args, nargout)
3544 {
3545  if (!pHP) {
3546  error("%s: not connected", type_name().c_str());
3547  return octave_value();
3548  }
3549 
3550  if (args.length() != 1) {
3551  error("%s: invalid number of arguments\n"
3552  "expected one argument", type_name().c_str());
3553 
3554  return octave_value();
3555  }
3556 
3557  const StructDispNode* const pNode = GetStructNode(args(0));
3558 
3559  if (pNode == 0) {
3560  return octave_value();
3561  }
3562 
3563  Mat3x3 mbR;
3564 
3565  try {
3566  mbR = pHP->GetRotRel(ReferenceFrame(pNode));
3567  } catch(...) { // do not throw across call back boundaries
3568  error("%s: GetRotRel failed", type_name().c_str());
3569  return octave_value();
3570  }
3571 
3572  octave_value octR;
3573 
3574  if (!GetInterface()->ConvertMBDynToOctave(mbR, octR)) {
3575  error("%s: could not convert data", type_name().c_str());
3576  return octave_value();
3577  }
3578 
3579  return octR;
3580 }
3581 
3582 METHOD_DEFINE(MBDynParserInterface, GetValue, args, nargout)
3583 {
3584  if (!pHP) {
3585  error("%s: not connected", type_name().c_str());
3586  return octave_value();
3587  }
3588 
3589  if (args.length() > 1) {
3590  error("%s: invalid number of arguments\n"
3591  "expected 0-1 arguments", type_name().c_str());
3592  return octave_value();
3593  }
3594 
3595  TypedValue mbDefValue;
3596 
3597  if (args.length() >= 1) {
3598  if (!GetInterface()->ConvertOctaveToMBDyn(args(0), mbDefValue)) {
3599  error("%s: could not convert octave data type %s to MBDyn",
3600  type_name().c_str(),
3601  args(0).type_name().c_str());
3602  return octave_value();
3603  }
3604  }
3605 
3606  TypedValue mbValue;
3607 
3608  try {
3609  mbValue = pHP->GetValue(mbDefValue);
3610  } catch (...) {
3611  error("%s: GetValue failed", type_name().c_str());
3612  return octave_value();
3613  }
3614 
3615  octave_value octValue;
3616 
3617  if (!GetInterface()->ConvertMBDynToOctave(mbValue, octValue)) {
3618  error("%s: could not convert data type %s into octave value",
3619  type_name().c_str(),
3620  mbValue.GetTypeName());
3621  }
3622 
3623  return octValue;
3624 }
3625 
3626 METHOD_DEFINE(MBDynParserInterface, GetLineData, args, nargout)
3627 {
3628  if (!pHP) {
3629  error("%s: not connected", type_name().c_str());
3630  return octave_value();
3631  }
3632 
3633  if (args.length() != 0) {
3634  error("%s: invalid number of arguments\n"
3635  "expected no argument", type_name().c_str());
3636 
3637  return octave_value();
3638  }
3639 
3640  std::ostringstream os;
3641 
3642  os << pHP->GetLineData();
3643 
3644  return octave_value(os.str());
3645 }
3646 
3647 const StructDispNode* MBDynParserInterface::GetStructNode(const octave_value& arg) const
3648 {
3649  const StructDispNodeBaseInterface* const pNodeInterface =
3650  dynamic_cast<const StructDispNodeBaseInterface*>(&arg.get_rep());
3651 
3652  if (pNodeInterface == 0) {
3653  error("%s: invalid argument type (%s)\n"
3654  "expected StructNode or StructDispNode",
3655  type_name().c_str(),
3656  arg.type_name().c_str());
3657 
3658  return 0;
3659  }
3660 
3661  const StructDispNode* const pNode = pNodeInterface->Get();
3662 
3663  if (pNode == 0) {
3664  error("%s: not connected", pNodeInterface->type_name().c_str());
3665  return 0;
3666  }
3667 
3668  return pNode;
3669 }
3670 
3671 bool MBDynParserInterface::GetDelimsFromString(const std::string& strDelims, HighParser::Delims& delims)
3672 {
3673  const int count = sizeof(mbStringDelims) / sizeof(mbStringDelims[0]);
3674 
3675  for (int i = 0; i < count; ++i) {
3676  if (strDelims == mbStringDelims[i].name) {
3677  delims = mbStringDelims[i].value;
3678  return true;
3679  }
3680  }
3681 
3682  delims = HighParser::DEFAULTDELIM;
3683 
3684  return false;
3685 }
3686 
3687 METHOD_DEFINE(MBDynParserInterface, GetDriveCaller, args, nargout)
3688 {
3689  if (!pHP) {
3690  error("%s: not connected", type_name().c_str());
3691  return octave_value();
3692  }
3693 
3694  if (args.length() != 0) {
3695  error("%s: invalid number of arguments\n"
3696  "expected no argument", type_name().c_str());
3697 
3698  return octave_value();
3699  }
3700 
3701  DriveCaller* pDC = 0;
3702 
3703  try {
3704  pDC = pHP->GetDriveCaller();
3705  } catch (...) {
3706  error("%s: GetDriveCaller failed", type_name().c_str());
3707  return octave_value();
3708  }
3709 
3710  DriveCallerInterface* pInterface = new DriveCallerInterface(GetInterface(), pDC);
3711 
3712  return octave_value(pInterface);
3713 }
3714 
3715 BEGIN_METHOD_TABLE(MBDynParserInterface, MBDynInterface)
3716  METHOD_DISPATCH(MBDynParserInterface, IsKeyWord)
3717  METHOD_DISPATCH(MBDynParserInterface, IsArg)
3718  METHOD_DISPATCH(MBDynParserInterface, IsStringWithDelims)
3719  METHOD_DISPATCH(MBDynParserInterface, GetReal)
3720  METHOD_DISPATCH(MBDynParserInterface, GetInt)
3721  METHOD_DISPATCH(MBDynParserInterface, GetBool)
3722  METHOD_DISPATCH(MBDynParserInterface, GetYesNoOrBool)
3723  METHOD_DISPATCH(MBDynParserInterface, GetString)
3724  METHOD_DISPATCH(MBDynParserInterface, GetStringWithDelims)
3725  METHOD_DISPATCH(MBDynParserInterface, GetValue)
3726  METHOD_DISPATCH(MBDynParserInterface, GetPosRel)
3727  METHOD_DISPATCH(MBDynParserInterface, GetRotRel)
3728  METHOD_DISPATCH(MBDynParserInterface, GetLineData)
3729  METHOD_DISPATCH(MBDynParserInterface, GetDriveCaller)
3730 END_METHOD_TABLE()
3731 
3732 DEFINE_OCTAVE_ALLOCATOR(MBDynParserInterface);
3733 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA(MBDynParserInterface, "MBDynParser", "MBDynParser");
3734 
3735 DriveCallerInterface::DriveCallerInterface(OctaveInterface* pInterface, const DriveCaller* pDC)
3736 :MBDynInterface(pInterface), DC(pDC)
3737 {
3738 
3739 }
3740 
3741 DriveCallerInterface::~DriveCallerInterface()
3742 {
3743 
3744 }
3745 
3746 
3747 METHOD_DEFINE(DriveCallerInterface, dGet, args, nargout)
3748 {
3749  octave_value y;
3750 
3751  if (!DC.pGetDriveCaller()) {
3752  error("%s: not connected", type_name().c_str());
3753  return y;
3754  }
3755 
3756  switch (args.length())
3757  {
3758  case 0:
3759  y = DC.dGet();
3760  break;
3761 
3762  case 1:
3763  {
3764  doublereal x;
3765 
3766  if (!GetInterface()->ConvertOctaveToMBDyn(args(0), x)) {
3767  error("%s: invalid argument type %s", type_name().c_str(), args(0).type_name().c_str());
3768  break;
3769  }
3770 
3771  y = DC.dGet(x);
3772  break;
3773  }
3774  default:
3775  error("%s: invalid number of arguments\n"
3776  "expected zero or one argument", type_name().c_str());
3777  };
3778 
3779  return y;
3780 }
3781 
3782 METHOD_DEFINE(DriveCallerInterface, dGetP, args, nargout)
3783 {
3784  octave_value y;
3785 
3786  if (!DC.pGetDriveCaller()) {
3787  error("%s: not connected", type_name().c_str());
3788  return y;
3789  }
3790 
3791  if (!DC.pGetDriveCaller()->bIsDifferentiable()) {
3792  error("%s: drive caller(%d %s) is not differentiable", type_name().c_str(), DC.pGetDriveCaller()->GetLabel(), DC.pGetDriveCaller()->GetName().c_str());
3793  return y;
3794  }
3795 
3796  switch (args.length())
3797  {
3798  case 0:
3799  y = DC.dGetP();
3800  break;
3801 
3802  case 1:
3803  {
3804  doublereal x;
3805 
3806  if (!GetInterface()->ConvertOctaveToMBDyn(args(0), x)) {
3807  error("%s: invalid argument type %s", type_name().c_str(), args(0).type_name().c_str());
3808  break;
3809  }
3810 
3811  y = DC.dGetP(x);
3812  break;
3813  }
3814  default:
3815  error("%s: invalid number of arguments\n"
3816  "expected zero or one argument", type_name().c_str());
3817  };
3818 
3819  return y;
3820 }
3821 
3822 BEGIN_METHOD_TABLE(DriveCallerInterface, MBDynInterface)
3823  METHOD_DISPATCH(DriveCallerInterface, dGet)
3824  METHOD_DISPATCH(DriveCallerInterface, dGetP)
3825 END_METHOD_TABLE()
3826 
3827 DEFINE_OCTAVE_ALLOCATOR(DriveCallerInterface);
3828 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA(DriveCallerInterface, "DriveCaller", "DriveCaller");
3829 
3830 OctaveDriveCaller::OctaveDriveCaller(const std::string& strFunc, OctaveInterface* pInterface, int iFlags, const octave_value_list& args)
3831 : DriveCaller(0),
3832  iFlags(iFlags),
3833  strFunc(strFunc),
3834  pInterface(pInterface),
3835  args(args)
3836 {
3837  TRACE("constructor");
3838  TRACE("strFunc=" << strFunc);
3839  TRACE("iFlags=" << iFlags);
3840 }
3841 
3842 OctaveDriveCaller::~OctaveDriveCaller(void)
3843 {
3844  TRACE("destructor");
3845 }
3846 
3847 DriveCaller *
3848 OctaveDriveCaller::pCopy(void) const
3849 {
3850  return new OctaveDriveCaller(strFunc, pInterface, iFlags, args);
3851 }
3852 
3853 std::ostream&
3854 OctaveDriveCaller::Restart(std::ostream& out) const
3855 {
3856  return out << "octave, \"" << strFunc << "\"";
3857 }
3858 
3859 inline doublereal
3860 OctaveDriveCaller::dGet(const doublereal& dVar) const
3861 {
3862  return pInterface->EvalScalarFunction(strFunc, MakeArgList(dVar), iFlags);
3863 }
3864 
3865 inline doublereal
3866 OctaveDriveCaller::dGet(void) const
3867 {
3868  return dGet(pInterface->GetDataManager()->dGetTime());
3869 }
3870 
3871 inline bool
3872 OctaveDriveCaller::bIsDifferentiable(void) const
3873 {
3874  return pInterface->HaveADPackage();
3875 }
3876 
3877 inline doublereal
3878 OctaveDriveCaller::dGetP(const doublereal& dVar) const
3879 {
3880  return pInterface->EvalScalarFunctionDerivative(strFunc, MakeArgList(dVar), iFlags);
3881 }
3882 
3883 inline doublereal
3884 OctaveDriveCaller::dGetP(void) const
3885 {
3886  if (!bIsDifferentiable()) {
3887  return 0.;
3888  }
3889 
3890  return dGetP(pInterface->GetDataManager()->dGetTime());
3891 }
3892 
3893 octave_value_list OctaveDriveCaller::MakeArgList(doublereal dVar) const
3894 {
3895  return pInterface->MakeArgList(dVar, args, iFlags);
3896 }
3897 
3898 template <class T>
3899 OctaveTplDriveCaller<T>::OctaveTplDriveCaller(const std::string& strFunction, OctaveInterface* pInterface, int iFlags, const octave_value_list& args)
3900  :strFunction(strFunction),
3901  pInterface(pInterface),
3902  iFlags(iFlags),
3903  args(args)
3904 {
3905  TRACE("constructor");
3906 };
3907 
3908 template <class T>
3909 OctaveTplDriveCaller<T>::~OctaveTplDriveCaller(void)
3910 {
3911  TRACE("destructor");
3912 };
3913 
3914 
3915 template <class T>
3916 TplDriveCaller<T>* OctaveTplDriveCaller<T>::pCopy(void) const
3917 {
3918  return new OctaveTplDriveCaller(strFunction, pInterface, iFlags, args);
3919 };
3920 
3921 template <class T>
3922 std::ostream& OctaveTplDriveCaller<T>::Restart(std::ostream& out) const
3923 {
3924  return out;
3925 };
3926 
3927 template <class T>
3928 std::ostream& OctaveTplDriveCaller<T>::Restart_int(std::ostream& out) const
3929 {
3930  return out;
3931 };
3932 
3933 template <class T>
3934 T OctaveTplDriveCaller<T>::Get(const doublereal& dVar) const
3935 {
3936  T X;
3937  pInterface->EvalMatrixFunction(strFunction, X, MakeArgList(dVar), iFlags);
3938  return X;
3939 };
3940 
3941 template <class T>
3942 T OctaveTplDriveCaller<T>::Get(void) const
3943 {
3944  return Get(pInterface->GetDataManager()->dGetTime());
3945 };
3946 
3947 template <class T>
3948 bool OctaveTplDriveCaller<T>::bIsDifferentiable(void) const
3949 {
3950  return pInterface->HaveADPackage();
3951 };
3952 
3953 template <class T>
3954 T OctaveTplDriveCaller<T>::GetP(void) const
3955 {
3956  const doublereal t = pInterface->GetDataManager()->dGetTime();
3957  T XP;
3958  pInterface->EvalMatrixFunctionDerivative(strFunction, XP, MakeArgList(t), iFlags);
3959  return XP;
3960 };
3961 
3962 template <class T>
3963 int OctaveTplDriveCaller<T>::getNDrives(void) const
3964 {
3965  return 0;
3966 };
3967 
3968 template <class T>
3969 octave_value_list OctaveTplDriveCaller<T>::MakeArgList(doublereal dVar) const
3970 {
3971  return pInterface->MakeArgList(dVar, args, iFlags);
3972 }
3973 
3974 template <class T>
3976 OctaveTDCR<T>::Read(const DataManager* pDM, MBDynParser& HP)
3977 {
3978  OctaveFunctionDCR::Read(pDM, HP);
3979 
3980  return new OctaveTplDriveCaller<T>(GetFunction(), GetInterface(), GetFlags(), GetArgs());
3981 };
3982 
3983 DerivativeDriveCaller::DerivativeDriveCaller(DriveCaller* pDriveCaller)
3984  :DriveCaller(0), pDriveCaller(pDriveCaller)
3985 {
3986  TRACE("constructor");
3987  ASSERT(pDriveCaller->bIsDifferentiable());
3988 }
3989 
3990 DerivativeDriveCaller::~DerivativeDriveCaller(void)
3991 {
3992  TRACE("destructor");
3993 }
3994 
3995 DriveCaller* DerivativeDriveCaller::pCopy(void) const
3996 {
3997  return new DerivativeDriveCaller(pDriveCaller.pGetDriveCaller()->pCopy());
3998 }
3999 
4000 std::ostream& DerivativeDriveCaller::Restart(std::ostream& out) const
4001 {
4002  return out;
4003 }
4004 
4005 doublereal DerivativeDriveCaller::dGet(const doublereal& dVar) const
4006 {
4007  return pDriveCaller.dGetP(dVar);
4008 }
4009 
4010 doublereal DerivativeDriveCaller::dGet(void) const
4011 {
4012  return pDriveCaller.dGetP();
4013 }
4014 
4015 bool DerivativeDriveCaller::bIsDifferentiable(void) const
4016 {
4017  return false;
4018 }
4019 
4020 doublereal DerivativeDriveCaller::dGetP(const doublereal& dVar) const
4021 {
4022  return 0.;
4023 }
4024 
4025 doublereal DerivativeDriveCaller::dGetP(void) const
4026 {
4027  return 0.;
4028 }
4029 
4030 OctaveScalarFunction::OctaveScalarFunction(const std::string& strFunc, OctaveInterface* pInterface, int iFlags, const octave_value_list& args)
4031  :strFunc(strFunc),
4032  pInterface(pInterface),
4033  iFlags(iFlags),
4034  args(args)
4035 {
4036 
4037 }
4038 
4039 OctaveScalarFunction::~OctaveScalarFunction(void)
4040 {
4041 
4042 }
4043 
4044 doublereal OctaveScalarFunction::operator()(const doublereal x) const
4045 {
4046  return pInterface->EvalScalarFunction(strFunc, MakeArgList(x), iFlags);
4047 }
4048 
4049 doublereal OctaveScalarFunction::ComputeDiff(const doublereal t, const integer order) const
4050 {
4051  if (order != 1) {
4052  silent_cerr("octave scalar function \"" << strFunc << "\" derivative of order " << order << " not supported" << std::endl);
4054  }
4055 
4056  return pInterface->EvalScalarFunctionDerivative(strFunc, MakeArgList(t), iFlags);
4057 }
4058 
4059 octave_value_list OctaveScalarFunction::MakeArgList(doublereal dVar) const
4060 {
4061  return pInterface->MakeArgList(dVar, args, iFlags);
4062 }
4063 
4064 OctaveConstitutiveLawBase::OctaveConstitutiveLawBase(const std::string& strClass, OctaveInterface* pInterface, int iFlags)
4065  : strClass(strClass),
4066  pInterface(pInterface),
4067  iFlags(iFlags)
4068 {
4069 
4070 }
4071 
4072 bool OctaveConstitutiveLawBase::bHaveMethod(const std::string& strName) const
4073 {
4074  return GetInterface()->bHaveMethod(octObject, strClass, strName);
4075 }
4076 
4077 const std::string OctaveConstitutiveLawBase::strGetConstLawType("GetConstLawType");
4078 const std::string OctaveConstitutiveLawBase::strUpdate("Update");
4079 
4080 template <class T, class Tder>
4081 OctaveConstitutiveLaw<T, Tder>::OctaveConstitutiveLaw(const std::string& strClass, OctaveInterface* pInterface, int iFlags)
4082 : OctaveConstitutiveLawBase(strClass, pInterface, iFlags),
4083  clType(ConstLawType::UNKNOWN)
4084 {
4085  octave_value_list args(GetInterface()->GetDataManagerInterface());
4086  args.append(GetInterface()->GetMBDynParserInterface());
4087 
4088  const int flags = GetFlags() | OctaveInterface::UPDATE_OCTAVE_VARIABLES;
4089 
4090  octave_value_list ans = GetInterface()->EvalFunction(GetClass(), args, 1, flags);
4091 
4092  ASSERT(ans.length() == 1);
4093 
4094  octObject = ans(0);
4095 
4096  if (!octObject.is_object()) {
4097  silent_cerr("octave constitutive law(" << Base_t::GetLabel() << "): result of constructor ("
4098  << ans(0).type_name() << ") is not an object at line "
4099  << GetInterface()->GetMBDynParser()->GetLineData() << std::endl);
4101  }
4102 
4103  if (!bHaveMethod(strGetConstLawType)) {
4104  silent_cerr("octave constitutive law(" << Base_t::GetLabel()
4105  << "): method " << GetClass() << "." << strGetConstLawType
4106  << " is undefined" << std::endl);
4108  }
4109 
4110  if (!bHaveMethod(strUpdate)) {
4111  silent_cerr("octave constitutive law(" << Base_t::GetLabel()
4112  << "): method " << GetClass() << "." << strUpdate
4113  << " is undefined" << std::endl);
4115  }
4116 };
4117 
4118 template <class T, class Tder>
4119 OctaveConstitutiveLaw<T, Tder>::~OctaveConstitutiveLaw(void) {
4120  NO_OP;
4121 };
4122 
4123 template <class T, class Tder>
4124 ConstLawType::Type OctaveConstitutiveLaw<T, Tder>::GetConstLawType(void) const {
4125 
4126  switch (clType) {
4127  case ConstLawType::UNKNOWN:
4128  break;
4129  default:
4130  return clType;
4131  }
4132 
4133  static const struct {
4134  ConstLawType::Type value;
4135  char name[13];
4136  } types[] = {
4137  { ConstLawType::ELASTIC, "ELASTIC" },
4138  { ConstLawType::VISCOUS, "VISCOUS" },
4139  { ConstLawType::VISCOELASTIC, "VISCOELASTIC" }
4140  };
4141 
4142  static const int count = sizeof(types)/sizeof(types[0]);
4143 
4144  octave_value_list args(octObject);
4145  const octave_value_list ans = GetInterface()->EvalFunction(strGetConstLawType, args, 1, GetFlags());
4146 
4147  ASSERT(ans.length() == 1);
4148 
4149  if (!ans(0).is_string()) {
4150  silent_cerr("octave constitutive law(" << Base_t::GetLabel()
4151  << "): " << GetClass() << "." << strGetConstLawType
4152  << " returned a invalid data type (" << ans(0).type_name() << std::endl
4153  << "expected string" << std::endl);
4154 
4156  }
4157 
4158  const std::string strCLType(ans(0).string_value());
4159 
4160  for (int i = 0; i < count; ++i) {
4161  if (strCLType == types[i].name) {
4162  clType = types[i].value;
4163  return clType;
4164  }
4165  }
4166 
4167  silent_cerr("octave constitutive law(" << Base_t::GetLabel()
4168  << "): " << GetClass() << "." << strGetConstLawType
4169  << " returned a invalid value (" << strCLType << ")" << std::endl);
4170 
4172 };
4173 
4174 template <class T, class Tder>
4175 ConstitutiveLaw<T, Tder>* OctaveConstitutiveLaw<T, Tder>::pCopy(void) const {
4176  Base_t* pCL = NULL;
4177 
4178  typedef OctaveConstitutiveLaw<T, Tder> cl;
4179  SAFENEWWITHCONSTRUCTOR(pCL, cl, cl(GetClass(), GetInterface(), GetFlags()));
4180  return pCL;
4181 };
4182 
4183 template <class T, class Tder>
4184 std::ostream& OctaveConstitutiveLaw<T, Tder>::Restart(std::ostream& out) const {
4185  return out << "# octave constitutive law not implemented" << std::endl;
4186 };
4187 
4188 template <class T, class Tder>
4189 void OctaveConstitutiveLaw<T, Tder>::Update(const T& mbEps, const T& mbEpsPrime) {
4190 
4191  octave_value octEps, octEpsPrime;
4192 
4193  if (!GetInterface()->ConvertMBDynToOctave(mbEps, octEps)) {
4194  silent_cerr("octave constitutive law(" << Base_t::GetLabel()
4195  << "): failed to convert MBDyn data type to octave_value" << std::endl);
4196 
4197  ASSERT(0);
4199  }
4200 
4201  if (!GetInterface()->ConvertMBDynToOctave(mbEpsPrime, octEpsPrime)) {
4202  silent_cerr("octave constitutive law(" << Base_t::GetLabel()
4203  << "): failed to convert MBDyn data type to octave_value" << std::endl);
4204 
4205  ASSERT(0);
4207  }
4208 
4209  octave_value_list args(octObject);
4210  args.append(octEps);
4211  args.append(octEpsPrime);
4212 
4213  const octave_value_list ans = GetInterface()->EvalFunction(strUpdate, args, 4, GetFlags());
4214 
4215  ASSERT(ans.length() == 4);
4216 
4217  if (!(ans(0).is_object() && ans(0).type_id() == octObject.type_id())) {
4218  silent_cerr("octave constitutive law(" << Base_t::GetLabel()
4219  << "): output argument 1 is not an instance of " << octObject.type_name() << std::endl);
4220 
4222  }
4223 
4224  octObject = ans(0);
4225 
4226  if (!((ans(1).is_real_matrix() || ans(1).is_real_scalar()) && ans(1).columns() == 1)) {
4227  silent_cerr("octave constitutive law(" << Base_t::GetLabel()
4228  << "): data type of output argument F is invalid (" << ans(1).type_name() << ")" << std::endl
4229  << "expected real vector" << std::endl);
4230 
4232  }
4233 
4234  const ColumnVector octF(ans(1).column_vector_value());
4235 
4236  if (!(ans(2).is_real_matrix() || ans(2).is_real_scalar())) {
4237  silent_cerr("octave constitutive law(" << Base_t::GetLabel()
4238  << "): data type of output argument FDE is invalid (" << ans(2).type_name() << ")" << std::endl
4239  << "expected real matrix" << std::endl);
4240 
4242  }
4243 
4244  const Matrix octFDE(ans(2).matrix_value());
4245 
4246  if (!(ans(3).is_real_matrix() || ans(3).is_real_scalar())) {
4247  silent_cerr("octave constitutive law(" << Base_t::GetLabel()
4248  << "): data type of output argument FDEPrime is invalid (" << ans(3).type_name() << ")" << std::endl
4249  << "expected real matrix" << std::endl);
4250 
4252  }
4253 
4254  const Matrix octFDEPrime(ans(3).matrix_value());
4255 
4256  Base_t::Epsilon = mbEps;
4257  Base_t::EpsilonPrime = mbEpsPrime;
4258 
4259  if (!GetInterface()->ConvertOctaveToMBDyn(octF, Base_t::F)) {
4260  silent_cerr("octave constitutive law(" << Base_t::GetLabel()
4261  << "): data type of output argument F is invalid "
4262  << octF.rows() << "x" << octF.columns() << " (" << ans(1).type_name() << ")" << std::endl);
4263 
4265  }
4266 
4267  if (!GetInterface()->ConvertOctaveToMBDyn(octFDE, Base_t::FDE)) {
4268  silent_cerr("octave constitutive law(" << Base_t::GetLabel()
4269  << "): data type of output argument FDE is invalid "
4270  << octFDE.rows() << "x" << octFDE.columns() << " (" << ans(2).type_name() << ")" << std::endl);
4271 
4273  }
4274 
4275  if (!GetInterface()->ConvertOctaveToMBDyn(octFDEPrime, Base_t::FDEPrime)) {
4276  silent_cerr("octave constitutive law(" << Base_t::GetLabel()
4277  << "): data type of output argument FDEPrime is invalid "
4278  << octFDEPrime.rows() << "x" << octFDEPrime.columns() << " (" << ans(3).type_name() << ")" << std::endl);
4279 
4281  }
4282 };
4283 
4284 OctaveBaseDCR::OctaveBaseDCR()
4285  :pInterface(0), iFlags(OctaveInterface::DEFAULT_CALL_FLAGS)
4286 {
4287  TRACE("constructor");
4288  TRACE("iFlags=" << iFlags);
4289 }
4290 
4291 OctaveBaseDCR::~OctaveBaseDCR()
4292 {
4293  TRACE("destructor");
4294 
4295  if (pInterface) {
4296  pInterface->Destroy();
4297  }
4298 }
4299 
4300 void OctaveBaseDCR::Read(const DataManager* pDM, MBDynParser& HP, bool bDeferred)const
4301 {
4302  if (pInterface == 0) {
4303  pInterface = OctaveInterface::CreateInterface(pDM, &HP);
4304  }
4305 
4306  strFunction = HP.GetStringWithDelims(HighParser::DOUBLEQUOTE);
4307 
4308  iFlags = OctaveInterface::DEFAULT_CALL_FLAGS;
4309 
4310  if (HP.IsKeyWord("update" "octave" "variables") && HP.GetYesNoOrBool()) {
4311  iFlags |= OctaveInterface::UPDATE_OCTAVE_VARIABLES;
4312  }
4313 
4314  if (HP.IsKeyWord("update" "mbdyn" "variables") && HP.GetYesNoOrBool()) {
4315  iFlags |= OctaveInterface::UPDATE_MBDYN_VARIABLES;
4316  }
4317 
4318  if (HP.IsKeyWord("update" "global" "variables") && HP.GetYesNoOrBool()) {
4319  iFlags |= OctaveInterface::UPDATE_GLOBAL_VARIABLES;
4320  }
4321 
4322  if (HP.IsKeyWord("pass" "data" "manager") && HP.GetYesNoOrBool()) {
4323  iFlags |= OctaveInterface::PASS_DATA_MANAGER;
4324  }
4325 
4326  if (HP.IsKeyWord("embed" "octave") && HP.GetYesNoOrBool()) {
4327  IncludeParser::ErrOut sFileInfo = HP.GetLineData();
4328  pInterface->AddEmbedFileName(sFileInfo.sFileName);
4329  }
4330 
4331  if (HP.IsKeyWord("octave" "search" "path")) {
4333  std::string path = HP.GetFileName(HighParser::DOUBLEQUOTE);
4334  pedantic_cout("octave: adding \"" << path << "\" to octave search path at line " << HP.GetLineData() << std::endl);
4335  if (!pInterface->AddOctaveSearchPath(path)) {
4336  silent_cerr("octave error: addpath(\"" << path << "\") failed at line " << HP.GetLineData() << std::endl);
4338  }
4339  }
4340  }
4341 
4342  TRACE("strFunction=" << strFunction);
4343  TRACE("iFlags=" << iFlags);
4344 }
4345 
4346 void OctaveFunctionDCR::Read(const DataManager* pDM, MBDynParser& HP, bool bDeferred) const
4347 {
4348  args.resize(0);
4349 
4350  OctaveBaseDCR::Read(pDM, HP, bDeferred);
4351 
4352  if (HP.IsKeyWord("arguments")) {
4353  integer iNumArgs = HP.GetInt();
4354  while (iNumArgs-- > 0) {
4355  TypedValue mbVal(HP.GetValue(TypedValue()));
4356  octave_value octVal;
4357 
4358  if (!GetInterface()->ConvertMBDynToOctave(mbVal, octVal)) {
4359  silent_cerr("octave error: could not convert MBDyn data type " << mbVal.GetTypeName() << " to octave" << std::endl);
4361  }
4362 
4363  args.append(octVal);
4364  }
4365  }
4366 }
4367 
4368 DriveCaller *
4369 OctaveDCR::Read(const DataManager* pDM, MBDynParser& HP, bool bDeferred)
4370 {
4371  OctaveFunctionDCR::Read(pDM, HP, bDeferred);
4372 
4373  return new OctaveDriveCaller(GetFunction(), GetInterface(), GetFlags(), GetArgs());
4374 };
4375 
4376 DriveCaller *
4377 DerivativeDCR::Read(const DataManager* pDM, MBDynParser& HP, bool bDeferred)
4378 {
4379  DriveCaller* pDriveCaller = HP.GetDriveCaller(bDeferred);
4380 
4381  if (!pDriveCaller->bIsDifferentiable()) {
4382  silent_cerr("DriveCaller " << pDriveCaller->GetLabel() << " " << pDriveCaller->GetName() << " is not differentiable at line " << HP.GetLineData() << std::endl);
4383 
4384  delete pDriveCaller;
4385 
4387  }
4388 
4389  return new DerivativeDriveCaller(pDriveCaller);
4390 };
4391 
4392 
4393 const BasicScalarFunction *
4394 OctaveSFR::Read(DataManager* pDM, MBDynParser& HP) const
4395 {
4396  OctaveFunctionDCR::Read(pDM, HP);
4397  return new OctaveScalarFunction(GetFunction(), GetInterface(), GetFlags(), GetArgs());
4398 };
4399 
4400 
4401 const std::string OctaveElement::strWorkSpaceDim("WorkSpaceDim");
4402 const std::string OctaveElement::striGetNumDof("iGetNumDof");
4403 const std::string OctaveElement::strAssRes("AssRes");
4404 const std::string OctaveElement::strAssJac("AssJac");
4405 const std::string OctaveElement::strUpdate("Update");
4406 const std::string OctaveElement::strSetValue("SetValue");
4407 const std::string OctaveElement::striGetInitialNumDof("iGetInitialNumDof");
4408 const std::string OctaveElement::strSetInitialValue("SetInitialValue");
4409 const std::string OctaveElement::strInitialAssRes("InitialAssRes");
4410 const std::string OctaveElement::strInitialAssJac("InitialAssJac");
4411 const std::string OctaveElement::strInitialWorkSpaceDim("InitialWorkSpaceDim");
4412 const std::string OctaveElement::strGetDofType("GetDofType");
4413 const std::string OctaveElement::strGetEqType("GetEqType");
4414 const std::string OctaveElement::strAfterConvergence("AfterConvergence");
4415 const std::string OctaveElement::striGetNumPrivData("iGetNumPrivData");
4416 const std::string OctaveElement::striGetPrivDataIdx("iGetPrivDataIdx");
4417 const std::string OctaveElement::strdGetPrivData("dGetPrivData");
4418 const std::string OctaveElement::striGetNumConnectedNodes("iGetNumConnectedNodes");
4419 const std::string OctaveElement::strGetConnectedNodes("GetConnectedNodes");
4420 const std::string OctaveElement::strOutput("Output");
4421 const std::string OctaveElement::strDescribeDof("DescribeDof");
4422 const std::string OctaveElement::strDescribeEq("DescribeEq");
4423 const std::string OctaveElement::strAfterPredict("AfterPredict");
4424 const std::string OctaveElement::strRestart("Restart");
4425 
4426 OctaveElement::OctaveElement(
4427  unsigned uLabel, const DofOwner *pDO,
4428  DataManager* pDM, MBDynParser& HP)
4429 : Elem(uLabel, flag(0)),
4430  UserDefinedElem(uLabel, pDO),
4431  haveMethod(HAVE_DEFAULT)
4432 {
4433  // help
4434  if (HP.IsKeyWord("help")) {
4435  silent_cout("Module: octave\n"
4436  << std::endl);
4437 
4438  if (!HP.IsArg()) {
4439  /*
4440  * Exit quietly if nothing else is provided
4441  */
4442  throw NoErr(MBDYN_EXCEPT_ARGS);
4443  }
4444  }
4445 
4446  dcr.Read(pDM, HP);
4447 
4448  X = new ConstVectorHandlerInterface(GetInterface());
4449  XP = new ConstVectorHandlerInterface(GetInterface());
4450  mbdObject = new OctaveElementInterface(GetInterface(), this);
4451  OS = new OStreamInterface(GetInterface());
4452 
4453  octave_value_list args(mbdObject);
4454 
4455  args.append(GetInterface()->GetDataManagerInterface());
4456  args.append(GetInterface()->GetMBDynParserInterface());
4457  args.append(OS);
4458  const int flags = OctaveInterface::UPDATE_OCTAVE_VARIABLES;
4459 
4460  OS->Set(&pDM->GetLogFile());
4461 
4462  octave_value_list ans = GetInterface()->EvalFunction(GetClass(), args, 1, flags);
4463 
4464  OS->Set(0);
4465 
4466  ASSERT(ans.length() == 1);
4467 
4468  octObject = ans(0);
4469 
4470  if (!octObject.is_object()) {
4471  silent_cerr("octave(" << GetLabel() << "): result of constructor ("
4472  << ans(0).type_name() << ") is not an object at line "
4473  << HP.GetLineData() << std::endl);
4475  }
4476 
4477  SetOutputFlag(pDM->fReadOutput(HP, Elem::LOADABLE));
4478 
4479  if (bHaveMethod(strAssJac)) {
4480  haveMethod |= HAVE_JACOBIAN;
4481  } else {
4482  pedantic_cerr("octave waring: method " << strAssJac << " not found in class " << GetClass() << std::endl);
4483  }
4484 
4485  if (bHaveMethod(strUpdate)) {
4486  haveMethod |= HAVE_UPDATE;
4487  } else {
4488  pedantic_cerr("octave warning: method " << strUpdate << " not found in class " << GetClass() << std::endl);
4489  }
4490 
4491  if (bHaveMethod(strAfterConvergence)) {
4492  haveMethod |= HAVE_AFTER_CONVERGENCE;
4493  } else {
4494  pedantic_cerr("octave warning: method " << strAfterConvergence << " not found in class " << GetClass() << std::endl);
4495  }
4496 
4497  if (bHaveMethod(strSetValue)) {
4498  haveMethod |= HAVE_SET_VALUE;
4499  } else {
4500  pedantic_cerr("octave warning: method " << strSetValue << " not found in class " << GetClass() << std::endl);
4501  }
4502 
4503  if (bHaveMethod(striGetNumDof)) {
4504  if (!bHaveMethod(strGetDofType)) {
4505  silent_cerr("octave error: method " << strGetDofType << " not found in class " << GetClass() << std::endl);
4507  }
4508 
4509  if (!bHaveMethod(strGetEqType)) {
4510  silent_cerr("octave error: method " << strGetEqType << " not found in class " << GetClass() << std::endl);
4512  }
4513 
4514  haveMethod |= HAVE_PRIVATE_DOF;
4515 
4516  if (bHaveMethod(strDescribeDof)) {
4517  haveMethod |= HAVE_DESCRIBE_DOF;
4518  }
4519 
4520  if (bHaveMethod(strDescribeEq)) {
4521  haveMethod |= HAVE_DESCRIBE_EQ;
4522  }
4523  } else {
4524  pedantic_cerr("octave warning: method " << striGetNumDof << " not found in class " << GetClass() << std::endl);
4525  }
4526 
4527  if (bHaveMethod(striGetNumPrivData)) {
4528  if (!bHaveMethod(striGetPrivDataIdx)) {
4529  silent_cerr("octave error: method " << striGetPrivDataIdx << " not found in class " << GetClass() << std::endl);
4531  }
4532 
4533  if (!bHaveMethod(strdGetPrivData)) {
4534  silent_cerr("octave error: method " << strdGetPrivData << " not found in class " << GetClass() << std::endl);
4536  }
4537 
4538  haveMethod |= HAVE_PRIVATE_DATA;
4539  } else {
4540  pedantic_cerr("octave warning: method " << striGetNumPrivData << " not found in class " << GetClass() << std::endl);
4541  }
4542 
4543  if (bHaveMethod(striGetNumConnectedNodes)) {
4544  if (!bHaveMethod(strGetConnectedNodes)) {
4545  silent_cerr("octave error: method " << strGetConnectedNodes << " not found in class " << GetClass() << std::endl);
4547  }
4548 
4549  haveMethod |= HAVE_CONNECTED_NODES;
4550  } else {
4551  pedantic_cerr("octave warning: method " << striGetNumConnectedNodes << " not found in class " << GetClass() << std::endl);
4552  }
4553 
4554  if (bHaveMethod(strInitialAssRes)) {
4555  if (!bHaveMethod(strInitialAssJac)) {
4556  silent_cerr("octave error: method " << strInitialAssJac << " not found in class " << GetClass() << std::endl);
4558  }
4559 
4560  if (!bHaveMethod(strInitialWorkSpaceDim)) {
4561  silent_cerr("octave error: method " << strInitialWorkSpaceDim << " not found in class " << GetClass() << std::endl);
4563  }
4564 
4565  if (!bHaveMethod(striGetInitialNumDof)) {
4566  silent_cerr("octave error: method " << striGetInitialNumDof << " not found in class " << GetClass() << std::endl);
4568  }
4569 
4570  haveMethod |= HAVE_INITIAL_ASSEMBLY;
4571 
4572  if (bHaveMethod(strSetInitialValue)) {
4573  haveMethod |= HAVE_SET_INITIAL_VALUE;
4574  } else {
4575  pedantic_cerr("octave warning: method " << strSetInitialValue << " not found in class " << GetClass() << std::endl);
4576  }
4577  } else {
4578  pedantic_cerr("octave warning: method " << strInitialAssRes
4579  << " not found in class " << GetClass() << std::endl
4580  << " initial assembly has been disabled" << std::endl);
4581  }
4582 
4583  if (bHaveMethod(strAfterPredict)) {
4584  haveMethod |= HAVE_AFTER_PREDICT;
4585  } else {
4586  pedantic_cerr("octave warning: method " << strAfterPredict << " not found in class " << GetClass() << std::endl);
4587  }
4588 
4589  if (bHaveMethod(strRestart)) {
4590  haveMethod |= HAVE_RESTART;
4591  } else {
4592  pedantic_cerr("octave warning: method " << strRestart << " not found in class " << GetClass() << std::endl);
4593  }
4594 
4595  if (bHaveMethod(strOutput)) {
4596  haveMethod |= HAVE_OUTPUT;
4597  } else {
4598  pedantic_cerr("octave warning: method " << strOutput << " not found in class " << GetClass() << std::endl);
4599  }
4600 }
4601 
4602 OctaveElement::~OctaveElement(void)
4603 {
4604 
4605 }
4606 
4607 void
4608 OctaveElement::Output(OutputHandler& OH) const
4609 {
4610  if (!(bToBeOutput() && OH.UseText(OutputHandler::LOADABLE))) {
4611  return;
4612  }
4613 
4614  if (!(haveMethod & HAVE_OUTPUT)) {
4615  return;
4616  }
4617 
4618  OS->Set(&OH.Loadable());
4619 
4620  octave_value_list args(octObject);
4621  args.append(OS);
4622 
4623  GetInterface()->EvalFunction(strOutput, args, 0, GetFlags());
4624 
4625  OS->Set(0);
4626 }
4627 
4628 void
4629 OctaveElement::WorkSpaceDim(integer* piNumRows, integer* piNumCols) const
4630 {
4631  octave_value_list args(octObject);
4632 
4633  octave_value_list ans = GetInterface()->EvalFunction(strWorkSpaceDim, args, 2, GetFlags());
4634 
4635  ASSERT(ans.length() == 2);
4636 
4637  if (!(ans(0).is_scalar_type()
4638  && ans(0).is_integer_type()
4639  && ans(1).is_scalar_type()
4640  && ans(1).is_integer_type())) {
4641  silent_cerr("octave(" << GetLabel() << "): function " << GetClass() << "." << strWorkSpaceDim
4642  << " must return two integer values\nreturned(" << ans(0).type_name() << ", " << ans(1).type_name() << ")" << std::endl);
4644  }
4645 
4646  *piNumRows = static_cast<int32_t>(ans(0).int32_scalar_value());
4647  *piNumCols = static_cast<int32_t>(ans(1).int32_scalar_value());
4648 }
4649 
4651 OctaveElement::AssJac(VariableSubMatrixHandler& WorkMatVar,
4652  doublereal dCoef,
4653  const VectorHandler& XCurr,
4654  const VectorHandler& XPrimeCurr)
4655 {
4656  if ( !(haveMethod & HAVE_JACOBIAN) ) {
4657  WorkMatVar.SetNullMatrix();
4658  return WorkMatVar;
4659  }
4660 
4661  integer iNumRows, iNumCols;
4662 
4663  WorkSpaceDim(&iNumRows, &iNumCols);
4664 
4665  octave_value_list args(octObject);
4666 
4667  X->Set(&XCurr);
4668  XP->Set(&XPrimeCurr);
4669 
4670  args.append(octave_value(dCoef));
4671  args.append(X);
4672  args.append(XP);
4673 
4674  const octave_value_list ans = GetInterface()->EvalFunction(strAssJac, args, 5, GetFlags() | OctaveInterface::OPTIONAL_OUTPUT_ARGS);
4675 
4676  ASSERT(ans.length() <= 5);
4677 
4678  X->Set(0);
4679  XP->Set(0);
4680 
4681  if (ans.length() >= 5) {
4682  if (!(ans(4).is_object() && ans(4).type_id() == octObject.type_id())) {
4683  silent_cerr("octave(" << GetLabel() << "): function " << GetClass() << "." << strAssRes << ": output argument 5 is not an instance of " << octObject.type_name() << std::endl);
4685  }
4686  octObject = ans(4);
4687  }
4688 
4689  return AssMatrix(WorkMatVar, ans, false);
4690 }
4691 
4693 OctaveElement::AssRes(SubVectorHandler& WorkVec,
4694  doublereal dCoef,
4695  const VectorHandler& XCurr,
4696  const VectorHandler& XPrimeCurr)
4697 {
4698  int iNumRows, iNumCols;
4699 
4700  WorkSpaceDim(&iNumRows, &iNumCols);
4701 
4702  WorkVec.ResizeReset(iNumRows);
4703 
4704  X->Set(&XCurr);
4705  XP->Set(&XPrimeCurr);
4706 
4707  octave_value_list args(octObject);
4708  args.append(octave_value(dCoef));
4709  args.append(X);
4710  args.append(XP);
4711 
4712  octave_value_list ans = GetInterface()->EvalFunction(strAssRes, args, 3, GetFlags() | OctaveInterface::OPTIONAL_OUTPUT_ARGS);
4713 
4714  X->Set(0);
4715  XP->Set(0);
4716 
4717  ASSERT(ans.length() <= 3);
4718 
4719  if (ans.length() < 2) {
4720  silent_cerr("octave(" << GetLabel() << "): function " << GetClass() << "." << strAssRes << ": expected 2-3 output arguments" << std::endl);
4722  }
4723 
4724  ASSERT(ans.length() >= 2);
4725  if (!(ans(0).is_matrix_type()
4726  && ans(0).is_real_type()
4727  && ans(0).columns() == 1)) {
4728  silent_cerr("octave(" << GetLabel()
4729  << "): function " << GetClass() << "." << strAssRes
4730  << ": output argument f (" << ans(0).type_name()
4731  << ") must be a real column vector" << std::endl);
4733  }
4734 
4735  const ColumnVector f = ans(0).column_vector_value();
4736 
4737  if (!(ans(1).is_integer_type()
4738  && ans(1).is_matrix_type()
4739  && ans(1).columns() == 1)) {
4740  silent_cerr("octave(" << GetLabel()
4741  << "): function " << GetClass() << "." << strAssRes
4742  << ": output argument ridx (" << ans(1).type_name()
4743  << ") must be an integer column vector" << std::endl);
4745  }
4746 
4747  const int32NDArray ridx = ans(1).int32_array_value();
4748 
4749  if (ridx.length() != iNumRows || f.length() != iNumRows) {
4750  silent_cerr("octave(" << GetLabel() << "): function " << GetClass() << "." << strAssRes
4751  << "\nlength(f)=" << f.length()
4752  << " length(ridx)=" << ridx.length()
4753  << " is not equal to iNumRows=" << iNumRows << std::endl);
4755  }
4756 
4757  const integer iNumDof = GetInterface()->GetDataManager()->iGetNumDofs();
4758 
4759  for (octave_idx_type i = 0; i < f.length(); ++i) {
4760  if (int(ridx(i)) <= 0 || int(ridx(i)) > iNumDof) {
4761  silent_cerr("octave(" << GetLabel() << "): function " << GetClass() << "." << strAssRes
4762  << ": row index " << ridx(i)
4763  << " out of range [1:" << iNumDof << "]" << std::endl);
4765  }
4766 
4767  WorkVec.PutItem(i + 1, ridx(i), f(i));
4768  }
4769 
4770  if (ans.length() >= 3) {
4771  if (!(ans(2).is_object() && ans(2).type_id() == octObject.type_id())) {
4772  silent_cerr("octave(" << GetLabel() << "): function " << GetClass() << "." << strAssRes << ": output argument 3 is not an instance of " << octObject.type_name() << std::endl);
4774  }
4775  octObject = ans(2);
4776  }
4777 
4778  return WorkVec;
4779 }
4780 
4781 unsigned int
4782 OctaveElement::iGetNumPrivData(void) const
4783 {
4784  if ( !(haveMethod & HAVE_PRIVATE_DATA) ) {
4785  return 0u;
4786  }
4787 
4788  octave_value_list args(octObject);
4789 
4790  octave_value_list ans = GetInterface()->EvalFunction(striGetNumPrivData, args, 1, GetFlags());
4791 
4792  ASSERT(ans.length() == 1);
4793 
4794  if ( !(ans(0).is_integer_type() && ans(0).is_scalar_type()) ) {
4795  silent_cerr("octave error: method " << GetClass() << "." << striGetNumPrivData
4796  << " returned (" << ans(0).type_name() << ")\n"
4797  << "expected integer scalar" << std::endl);
4799  }
4800 
4801  const int32_t iNumPrivData = ans(0).int32_scalar_value();
4802 
4803  if ( iNumPrivData < 0 ) {
4804  silent_cerr("octave error: method " << GetClass() << "." << striGetNumPrivData
4805  << " returned a negative value" << std::endl);
4807  }
4808 
4809  return iNumPrivData;
4810 }
4811 
4812 unsigned int
4813 OctaveElement::iGetPrivDataIdx(const char *s) const
4814 {
4815  if ( !(haveMethod & HAVE_PRIVATE_DATA) ) {
4816  ASSERT(0); // We should not get here!
4818  }
4819 
4820  octave_value_list args(octObject);
4821  args.append(octave_value(s));
4822 
4823  octave_value_list ans = GetInterface()->EvalFunction(striGetPrivDataIdx, args, 1, GetFlags());
4824 
4825  ASSERT(ans.length() == 1);
4826 
4827  if ( !(ans(0).is_integer_type() && ans(0).is_scalar_type()) ) {
4828  silent_cerr("octave error: method " << GetClass() << "." << striGetPrivDataIdx
4829  << " returned (" << ans(0).type_name() << ")\n"
4830  << "expected integer scalar" << std::endl);
4832  }
4833 
4834  const int32_t iPrivDataIdx = ans(0).int32_scalar_value();
4835 
4836  if ( iPrivDataIdx <= 0 ) {
4837  silent_cerr("octave error: method " << GetClass() << "." << striGetPrivDataIdx
4838  << " returned a negative or zero value" << std::endl);
4840  }
4841 
4842  return iPrivDataIdx;
4843 }
4844 
4845 doublereal
4846 OctaveElement::dGetPrivData(unsigned int i) const
4847 {
4848  if ( !(haveMethod & HAVE_PRIVATE_DATA) ) {
4849  ASSERT(0); // We should not get here!
4851  }
4852 
4853  octave_value_list args(octObject);
4854  args.append(octave_int32(i));
4855 
4856  octave_value_list ans = GetInterface()->EvalFunction(strdGetPrivData, args, 1, GetFlags());
4857 
4858  ASSERT(ans.length() == 1);
4859 
4860  if ( !(ans(0).is_real_scalar()) ) {
4861  silent_cerr("octave error: method " << GetClass() << "." << strdGetPrivData
4862  << " returned (" << ans(0).type_name() << ")\n"
4863  << "expected real scalar" << std::endl);
4865  }
4866 
4867  return ans(0).scalar_value();
4868 }
4869 
4870 int
4871 OctaveElement::iGetNumConnectedNodes(void) const
4872 {
4873  if (!(haveMethod & HAVE_CONNECTED_NODES)) {
4874  return 0u;
4875  }
4876 
4877  octave_value_list args(octObject);
4878 
4879  octave_value_list ans = GetInterface()->EvalFunction(striGetNumConnectedNodes, args, 1, GetFlags());
4880 
4881  ASSERT(ans.length() == 1);
4882 
4883  if ( !(ans(0).is_integer_type() && ans(0).is_scalar_type()) ) {
4884  silent_cerr("octave error: method " << GetClass() << "." << striGetNumConnectedNodes
4885  << " returned (" << ans(0).type_name() << ")\n"
4886  << "expected integer scalar" << std::endl);
4888  }
4889 
4890  const int32_t iNumConnectedNodes = ans(0).int32_scalar_value();
4891 
4892  if ( iNumConnectedNodes < 0 ) {
4893  silent_cerr("octave error: method " << GetClass() << "." << striGetNumConnectedNodes
4894  << " returned a negative value" << std::endl);
4896  }
4897 
4898  return iNumConnectedNodes;
4899 }
4900 
4901 void
4902 OctaveElement::GetConnectedNodes(std::vector<const Node *>& connectedNodes) const
4903 {
4904  if (!(haveMethod & HAVE_CONNECTED_NODES)) {
4905  connectedNodes.resize(0);
4906  return;
4907  }
4908 
4909  connectedNodes.resize(iGetNumConnectedNodes());
4910 
4911  octave_value_list args(octObject);
4912 
4913  octave_value_list ans = GetInterface()->EvalFunction(strGetConnectedNodes, args, 1, GetFlags());
4914 
4915  ASSERT(ans.length() == 1);
4916 
4917  if (!ans(0).is_cell()) {
4918  silent_cerr("octave error: method " << GetClass() << "." << strGetConnectedNodes
4919  << " returned (" << ans(0).type_name() << ")\n"
4920  << "expected cell" << std::endl);
4922  }
4923 
4924  const Cell nodes(ans(0).cell_value());
4925 
4926  if (size_t(nodes.length()) != connectedNodes.size()) {
4927  silent_cerr("octave error: method " << GetClass() << "." << strGetConnectedNodes
4928  << " returned an object of size " << nodes.length()
4929  << " whereas " << GetClass() << "." << striGetNumConnectedNodes
4930  << " returned "<< connectedNodes.size() << std::endl);
4932  }
4933 
4934  for (size_t i = 0; i < connectedNodes.size(); ++i) {
4935  const NodeInterface* pNode = dynamic_cast<const NodeInterface*>(&nodes(i).get_rep());
4936 
4937  if (pNode == 0) {
4938  silent_cerr("octave error: cell(" << i + 1 << ") (data type " << nodes(i).type_name() << ") is not a node" << std::endl);
4940  }
4941  connectedNodes[i] = pNode->Get();
4942  }
4943 }
4944 
4945 void
4946 OctaveElement::SetValue(DataManager *pDM,
4947  VectorHandler& X, VectorHandler& XP,
4949 {
4950  if (!(haveMethod & HAVE_SET_VALUE)) {
4951  return;
4952  }
4953 
4954  octave_value_list args(octObject);
4955  args.append(octave_value(new VectorHandlerInterface(GetInterface(), &X)));
4956  args.append(octave_value(new VectorHandlerInterface(GetInterface(), &XP)));
4957 
4958  GetInterface()->EvalFunction(strSetValue, args, 0, GetFlags());
4959 }
4960 
4961 unsigned int OctaveElement::iGetNumDof(void) const
4962 {
4963  if (!(haveMethod & HAVE_PRIVATE_DOF)) {
4964  return 0u;
4965  }
4966 
4967  octave_value_list args(octObject);
4968 
4969  octave_value_list ans = GetInterface()->EvalFunction(striGetNumDof, args, 1, GetFlags());
4970 
4971  ASSERT(ans.length() == 1);
4972 
4973  if (!(ans(0).is_integer_type() && ans(0).is_scalar_type())){
4974  silent_cerr("octave error: method " << GetClass() << "." << striGetNumDof
4975  << " returned (" << ans(0).type_name()
4976  << ")\nmethod must return an integer" << std::endl);
4978  }
4979 
4980  return static_cast<int32_t>(ans(0).int32_scalar_value());
4981 }
4982 
4983 DofOrder::Order OctaveElement::GetDofType(unsigned int i) const
4984 {
4985  if (!(haveMethod & HAVE_PRIVATE_DOF)) {
4986  ASSERT(0);
4988  }
4989 
4990  octave_value_list args(octObject);
4991  args.append(octave_int32(i + 1)); // FIXME: Why is this index zero based?
4992 
4993  octave_value_list ans = GetInterface()->EvalFunction(strGetDofType, args, 1, GetFlags());
4994 
4995  ASSERT(ans.length() == 1);
4996 
4997  if (ans(0).is_integer_type() && ans(0).is_scalar_type()) {
4998  const int32_t order = static_cast<int32_t>(ans(0).int32_scalar_value());
4999 
5000  switch ( order ) {
5001  case DofOrder::ALGEBRAIC:
5003  return static_cast<DofOrder::Order>(order);
5004  default:
5005  silent_cerr("octave error: method " << GetClass() << "." << strGetDofType
5006  << " returned an illegal value" << std::endl);
5008  }
5009  }
5010 
5011  if ( !ans(0).is_string() ) {
5012  silent_cerr("octave error: method " << GetClass() << "." << strGetDofType
5013  << " returned (" << ans(0).type_name()
5014  << ")\nmethod must return integer or string values" << std::endl);
5016  }
5017 
5018  const std::string order = ans(0).string_value();
5019 
5020  if ( order == "ALGEBRAIC" ) {
5021  return DofOrder::ALGEBRAIC;
5022  } else if ( order == "DIFFERENTIAL" ) {
5023  return DofOrder::DIFFERENTIAL;
5024  } else {
5025  silent_cerr("octave error: method " << GetClass() << "." << strGetDofType
5026  << " returned an illegal value (" << order << ")" << std::endl);
5028  }
5029 }
5030 
5031 DofOrder::Order OctaveElement::GetEqType(unsigned int i) const
5032 {
5033  if (!(haveMethod & HAVE_PRIVATE_DOF)) {
5034  ASSERT(0);
5036  }
5037 
5038  octave_value_list args(octObject);
5039  args.append(octave_int32(i + 1)); // FIXME: Why is this index zero based?
5040 
5041  octave_value_list ans = GetInterface()->EvalFunction(strGetEqType, args, 1, GetFlags());
5042 
5043  ASSERT(ans.length() == 1);
5044 
5045  if (ans(0).is_integer_type() && ans(0).is_scalar_type()) {
5046  const int32_t order = static_cast<int32_t>(ans(0).int32_scalar_value());
5047 
5048  switch ( order ) {
5049  case DofOrder::ALGEBRAIC:
5051  return static_cast<DofOrder::Order>(order);
5052  default:
5053  silent_cerr("octave error: method " << GetClass() << "." << strGetEqType
5054  << " returned an illegal value" << std::endl);
5056  }
5057  }
5058 
5059  if ( !ans(0).is_string() ) {
5060  silent_cerr("octave error: method " << GetClass() << "." << strGetEqType
5061  << " returned (" << ans(0).type_name()
5062  << ")\nmethod must return integer or string values" << std::endl);
5064  }
5065 
5066  const std::string order = ans(0).string_value();
5067 
5068  if ( order == "ALGEBRAIC" ) {
5069  return DofOrder::ALGEBRAIC;
5070  } else if ( order == "DIFFERENTIAL" ) {
5071  return DofOrder::DIFFERENTIAL;
5072  } else {
5073  silent_cerr("octave error: method " << GetClass() << "." << strGetEqType
5074  << " returned an illegal value (" << order << ")" << std::endl);
5076  }
5077 }
5078 
5079 std::ostream& OctaveElement::DescribeDof(std::ostream& out, const char *prefix, bool bInitial) const
5080 {
5081  if (!(haveMethod & HAVE_DESCRIBE_DOF)) {
5082  return out;
5083  }
5084 
5085  OS->Set(&out);
5086 
5087  octave_value_list args(octObject);
5088  args.append(OS);
5089  args.append(octave_value(prefix));
5090  args.append(octave_value(bInitial));
5091 
5092  GetInterface()->EvalFunction(strDescribeDof, args, 0, GetFlags());
5093 
5094  OS->Set(0);
5095 
5096  return out;
5097 }
5098 
5099 std::ostream& OctaveElement::DescribeEq(std::ostream& out, const char *prefix, bool bInitial) const
5100 {
5101  if (!(haveMethod & HAVE_DESCRIBE_EQ)) {
5102  return out;
5103  }
5104 
5105  OS->Set(&out);
5106 
5107  octave_value_list args(octObject);
5108  args.append(OS);
5109  args.append(octave_value(prefix));
5110  args.append(octave_value(bInitial));
5111 
5112  GetInterface()->EvalFunction(strDescribeEq, args, 0, GetFlags());
5113 
5114  OS->Set(0);
5115 
5116  return out;
5117 }
5118 
5119 void OctaveElement::Update(const VectorHandler& XCurr,const VectorHandler& XPrimeCurr)
5120 {
5121  if (!(haveMethod & HAVE_UPDATE)) {
5122  return;
5123  }
5124 
5125  X->Set(&XCurr);
5126  XP->Set(&XPrimeCurr);
5127 
5128  octave_value_list args(octObject);
5129  args.append(X);
5130  args.append(XP);
5131 
5132  octave_value_list ans = GetInterface()->EvalFunction(strUpdate, args, 1, GetFlags());
5133 
5134  ASSERT(ans.length() == 1);
5135 
5136  if (!(ans(0).is_object() && ans(0).type_id() == octObject.type_id())) {
5137  silent_cerr("octave(" << GetLabel() << "): function " << GetClass() << "." << strUpdate << ": output argument 1 is not an instance of " << octObject.type_name() << std::endl);
5139  }
5140 
5141  octObject = ans(0);
5142 
5143  X->Set(0);
5144  XP->Set(0);
5145 }
5146 
5147 void OctaveElement::AfterPredict(VectorHandler& XCurr, VectorHandler& XPrimeCurr)
5148 {
5149  if (!(haveMethod & HAVE_AFTER_PREDICT)) {
5150  return;
5151  }
5152 
5153  X->Set(&XCurr);
5154  XP->Set(&XPrimeCurr);
5155 
5156  octave_value_list args(octObject);
5157  args.append(X);
5158  args.append(XP);
5159 
5160  octave_value_list ans = GetInterface()->EvalFunction(strAfterPredict, args, 1, GetFlags());
5161 
5162  ASSERT(ans.length() == 1);
5163 
5164  if (!(ans(0).is_object() && ans(0).type_id() == octObject.type_id())) {
5165  silent_cerr("octave(" << GetLabel() << "): function " << GetClass() << "." << strAfterPredict << ": output argument 1 is not an instance of " << octObject.type_name() << std::endl);
5167  }
5168 
5169  octObject = ans(0);
5170 
5171  X->Set(0);
5172  XP->Set(0);
5173 }
5174 
5175 void OctaveElement::AfterConvergence(const VectorHandler& X,
5176  const VectorHandler& XP)
5177 {
5178  if ( !(haveMethod & HAVE_AFTER_CONVERGENCE) ) {
5179  return;
5180  }
5181 
5182  this->X->Set(&X);
5183  this->XP->Set(&XP);
5184 
5185  octave_value_list args(octObject);
5186  args.append(this->X);
5187  args.append(this->XP);
5188 
5189  octave_value_list ans = GetInterface()->EvalFunction(strAfterConvergence, args, 1, GetFlags());
5190 
5191  this->X->Set(0);
5192  this->XP->Set(0);
5193 
5194  ASSERT(ans.length() == 1);
5195 
5196  if ( !ans(0).is_object() ) {
5197  silent_cerr("octave error: return value of " << strAfterConvergence << " ("
5198  << ans(0).type_name() << ") is not an object" << std::endl);
5200  }
5201 
5202  octObject = ans(0);
5203 }
5204 
5205 std::ostream&
5206 OctaveElement::Restart(std::ostream& out) const
5207 {
5208  if (!(haveMethod & HAVE_RESTART)) {
5209  out << "# OctaveElement(" << GetClass() << "): not implemented" << std::endl;
5210  return out;
5211  }
5212 
5213  OS->Set(&out);
5214 
5215  octave_value_list args(octObject);
5216  args.append(OS);
5217 
5218  GetInterface()->EvalFunction(strRestart, args, 0, GetFlags());
5219 
5220  OS->Set(0);
5221 
5222  return out;
5223 }
5224 
5225 unsigned int
5226 OctaveElement::iGetInitialNumDof(void) const
5227 {
5228  if (!(haveMethod & HAVE_INITIAL_ASSEMBLY)) {
5229  return 0u;
5230  }
5231 
5232  octave_value_list args(octObject);
5233 
5234  octave_value_list ans = GetInterface()->EvalFunction(striGetInitialNumDof, args, 1, GetFlags());
5235 
5236  ASSERT(ans.length() == 1);
5237 
5238  if (!(ans(0).is_integer_type() && ans(0).is_scalar_type())){
5239  silent_cerr("octave error: method " << GetClass() << "." << striGetInitialNumDof
5240  << " returned (" << ans(0).type_name()
5241  << ")\nmethod must return an integer" << std::endl);
5243  }
5244 
5245  return static_cast<int32_t>(ans(0).int32_scalar_value());
5246 }
5247 
5248 void
5249 OctaveElement::InitialWorkSpaceDim(
5250  integer* piNumRows,
5251  integer* piNumCols) const
5252 {
5253  if (!(haveMethod & HAVE_INITIAL_ASSEMBLY)) {
5254  *piNumRows = 0;
5255  *piNumCols = 0;
5256  return;
5257  }
5258 
5259  octave_value_list args(octObject);
5260 
5261  octave_value_list ans = GetInterface()->EvalFunction(strInitialWorkSpaceDim, args, 2, GetFlags());
5262 
5263  ASSERT(ans.length() == 2);
5264 
5265  if (!(ans(0).is_integer_type()
5266  && ans(0).is_scalar_type()
5267  && ans(1).is_integer_type()
5268  && ans(1).is_scalar_type())) {
5269  silent_cerr("octave error: method " << GetClass() << "." << strInitialWorkSpaceDim
5270  << " must return two integer values\n"
5271  << "returned (" << ans(0).type_name() << ", " << ans(1).type_name() << ")" << std::endl);
5273  }
5274 
5275  *piNumRows = static_cast<int32_t>(ans(0).int32_scalar_value());
5276  *piNumCols = static_cast<int32_t>(ans(1).int32_scalar_value());
5277 }
5278 
5280 OctaveElement::InitialAssJac(
5281  VariableSubMatrixHandler& WorkMatVar,
5282  const VectorHandler& XCurr)
5283 {
5284  if (!(haveMethod & HAVE_INITIAL_ASSEMBLY)) {
5285  WorkMatVar.SetNullMatrix();
5286  return WorkMatVar;
5287  }
5288 
5289  octave_value_list args(octObject);
5290 
5291  X->Set(&XCurr);
5292 
5293  args.append(X);
5294 
5295  const octave_value_list ans = GetInterface()->EvalFunction(strInitialAssJac, args, 4, GetFlags() | OctaveInterface::OPTIONAL_OUTPUT_ARGS);
5296 
5297  ASSERT(ans.length() <= 4);
5298 
5299  X->Set(0);
5300 
5301  return AssMatrix(WorkMatVar, ans, true);
5302 }
5303 
5305 OctaveElement::InitialAssRes(
5306  SubVectorHandler& WorkVec,
5307  const VectorHandler& XCurr)
5308 {
5309  if ( !(haveMethod & HAVE_INITIAL_ASSEMBLY) ) {
5310  WorkVec.ResizeReset(0);
5311  return WorkVec;
5312  }
5313 
5314  int iNumRows, iNumCols;
5315 
5316  InitialWorkSpaceDim(&iNumRows, &iNumCols);
5317 
5318  WorkVec.ResizeReset(iNumRows);
5319 
5320  X->Set(&XCurr);
5321 
5322  octave_value_list args(octObject);
5323  args.append(X);
5324 
5325  octave_value_list ans = GetInterface()->EvalFunction(strInitialAssRes, args, 2, GetFlags());
5326 
5327  X->Set(0);
5328 
5329  ASSERT(ans.length() == 2);
5330 
5331  if (!(ans(0).is_matrix_type()
5332  && ans(0).is_real_type()
5333  && ans(0).columns() == 1)) {
5334  silent_cerr("octave(" << GetLabel() << "): function " << GetClass() << "." << strInitialAssRes
5335  << " output argument f (" << ans(0).type_name() << ") must be a real column vector" << std::endl);
5337  }
5338 
5339  const ColumnVector f = ans(0).column_vector_value();
5340 
5341  if (!(ans(1).is_integer_type()
5342  && ans(1).is_matrix_type()
5343  && ans(1).columns() == 1)) {
5344  silent_cerr("octave(" << GetLabel() << "): function " << GetClass() << "." << strInitialAssRes
5345  << ": output argument ridx (" << ans(1).type_name() << ") must be an integer column vector" << std::endl);
5347  }
5348 
5349  const int32NDArray ridx = ans(1).int32_array_value();
5350 
5351  if (ridx.length() != iNumRows || f.length() != iNumRows) {
5352  silent_cerr("octave(" << GetLabel()
5353  << "): function " << GetClass() << "." << strInitialAssRes
5354  << ": length(f)=" << f.length()
5355  << " length(ridx)=" << ridx.length()
5356  << " is not equal to iNumRows=" << iNumRows << std::endl);
5358  }
5359 
5360  const integer iNumDof = GetInterface()->GetDataManager()->iGetNumDofs();
5361 
5362  for (octave_idx_type i = 0; i < f.length(); ++i) {
5363  if (int(ridx(i)) <= 0 || int(ridx(i)) > iNumDof) {
5364  silent_cerr("octave(" << GetLabel() << "): function " << GetClass() << "." << strAssRes << ": row index " << ridx(i) << " out of range [1:" << iNumDof << "]" << std::endl);
5366  }
5367 
5368  WorkVec.PutItem(i + 1, ridx(i), f(i));
5369  }
5370 
5371  return WorkVec;
5372 }
5373 
5374 void OctaveElement::SetInitialValue(VectorHandler& X)
5375 {
5376  if (!(haveMethod & HAVE_SET_INITIAL_VALUE)) {
5377  return;
5378  }
5379 
5380  octave_value_list args(octObject);
5381  args.append(octave_value(new VectorHandlerInterface(GetInterface(), &X)));
5382 
5383  GetInterface()->EvalFunction(strSetInitialValue, args, 0, GetFlags());
5384 }
5385 
5386 bool OctaveElement::bHaveMethod(const std::string& strName)const
5387 {
5388  return GetInterface()->bHaveMethod(octObject, GetClass(), strName);
5389 }
5390 
5392 OctaveElement::AssMatrix(VariableSubMatrixHandler& WorkMatVar, const octave_value_list& ans, bool bInitial)
5393 {
5394  const std::string& strFunction = bInitial ? strInitialAssJac : strAssJac;
5395 
5396  ASSERT(ans.length() <= 5);
5397 
5398  if (ans.length() < 3) {
5399  silent_cerr("octave(" << GetLabel()
5400  << "): function " << GetClass() << "." << strFunction
5401  << " returned " << ans.length() << " output arguments\n"
5402  << "expected 3-4 output arguments" << std::endl);
5403 
5405  }
5406 
5407  ASSERT(ans.length() >= 3);
5408 
5409  if (!ans(0).is_real_matrix()) {
5410  silent_cerr("octave(" << GetLabel()
5411  << "): function " << GetClass() << "." << strFunction
5412  << " output argument Jac (" << ans(0).type_name()
5413  << ") must be a real matrix" << std::endl);
5415  }
5416 
5417  const Matrix Jac = ans(0).matrix_value();
5418 
5419  if (!(ans(1).is_integer_type() && ans(1).is_matrix_type())) {
5420  silent_cerr("octave(" << GetLabel()
5421  << "): function " << GetClass() << "." << strFunction
5422  << " output argument ridx (" << ans(1).type_name()
5423  << ") must be a vector of integer values" << std::endl);
5425  }
5426 
5427  const int32NDArray ridx = ans(1).int32_array_value();
5428 
5429  if (!(ans(2).is_integer_type() && ans(2).is_matrix_type())) {
5430  silent_cerr("octave(" << GetLabel()
5431  << "): function " << GetClass() << "." << strFunction
5432  << " output argument cidx (" << ans(2).type_name()
5433  << ") must be a vector of integer values" << std::endl);
5435  }
5436 
5437  const int32NDArray cidx = ans(2).int32_array_value();
5438 
5439  bool bSparse = false;
5440 
5441  if (ans.length() >= 4) {
5442  if (!ans(3).is_bool_scalar()) {
5443  silent_cerr("octave(" << GetLabel()
5444  << "): invalid data type (" << ans(3).type_name() << ")\n"
5445  << "expected bool scalar" << std::endl);
5447  }
5448 
5449  bSparse = ans(3).bool_value();
5450  }
5451 
5452  return AssMatrix(WorkMatVar, Jac, ridx, cidx, bSparse, bInitial);
5453 }
5454 
5456 OctaveElement::AssMatrix(VariableSubMatrixHandler& WorkMatVar, const Matrix& Jac, const int32NDArray& ridx, const int32NDArray& cidx, bool bSparse, bool bInitial)
5457 {
5458  const std::string& strFunction = bInitial ? strInitialAssJac : strAssJac;
5459 
5460  const integer iNumRows = Jac.rows();
5461  const integer iNumCols = Jac.cols();
5462 
5463  const integer iNumDof = GetInterface()->GetDataManager()->iGetNumDofs();
5464 
5465  if (bSparse) {
5466  if (iNumCols != 1
5467  || ridx.length() != iNumRows
5468  || cidx.length() != iNumRows) {
5469  silent_cerr("octave(" << GetLabel() << "):"
5470  << " rows(Jac)=" << Jac.rows()
5471  << " columns(Jac)= " << Jac.columns()
5472  << " length(ridx)= " << ridx.length()
5473  << " length(cidx)=" << cidx.length()
5474  << " are not consistent for a sparse submatrix" << std::endl);
5476  }
5477 
5478  SparseSubMatrixHandler& WorkMat = WorkMatVar.SetSparse();
5479  WorkMat.ResizeReset(iNumRows, 1);
5480 
5481  for (int i = 0; i < iNumRows; ++i) {
5482  if (int32_t(ridx(i)) <= 0 || int32_t(ridx(i)) > iNumDof) {
5483  silent_cerr("octave(" << GetLabel() << "):"
5484  << " function " << GetClass() << "." << strFunction
5485  << ": row index " << ridx(i)
5486  << " out of range [1:" << iNumDof << "]" << std::endl);
5488  }
5489 
5490  if (int32_t(cidx(i)) <= 0 || int32_t(cidx(i)) > iNumDof) {
5491  silent_cerr("octave(" << GetLabel() << "):"
5492  << " function " << GetClass() << "." << strFunction
5493  << ": column index " << cidx(i)
5494  << " out of range [1:" << iNumDof << "]" << std::endl);
5496  }
5497 
5498  WorkMat.PutItem(i + 1, ridx(i), cidx(i), Jac(i, 0));
5499  }
5500  } else {
5501  FullSubMatrixHandler& WorkMat = WorkMatVar.SetFull();
5502 
5503  WorkMat.ResizeReset(iNumRows, iNumCols);
5504 
5505  if (ridx.length() != iNumRows
5506  || cidx.length() != iNumCols) {
5507  silent_cerr("octave(" << GetLabel() << "):"
5508  << " rows(Jac)=" << Jac.rows()
5509  << " columns(Jac)= " << Jac.columns()
5510  << " length(ridx)= " << ridx.length()
5511  << " length(cidx)=" << cidx.length()
5512  << " are not consistent for a full submatrix" << std::endl);
5514  }
5515 
5516  for (int i = 0; i < iNumRows; ++i) {
5517  if (int32_t(ridx(i)) <= 0 || int32_t(ridx(i)) > iNumDof) {
5518  silent_cerr("octave(" << GetLabel() << "):"
5519  << " function " << GetClass() << "." << strFunction
5520  << ": row index " << ridx(i)
5521  << " out of range [1:" << iNumDof << "]" << std::endl);
5523  }
5524 
5525  WorkMat.PutRowIndex(i + 1, ridx(i));
5526  }
5527 
5528  for (int j = 0; j < iNumCols; ++j) {
5529  if (int32_t(cidx(j)) <= 0 || int32_t(cidx(j)) > iNumDof) {
5530  silent_cerr("octave(" << GetLabel() << "):"
5531  << " function " << GetClass() << "." << strFunction
5532  << ": column index " << cidx(j)
5533  << " out of range [1:" << iNumDof << "]" << std::endl);
5535  }
5536 
5537  WorkMat.PutColIndex(j + 1, cidx(j));
5538  }
5539 
5540  for (int i = 0; i < iNumRows; ++i) {
5541  for (int j = 0; j < iNumCols; ++j) {
5542  WorkMat.PutCoef(i + 1, j + 1, Jac(i, j));
5543  }
5544  }
5545  }
5546 
5547  return WorkMatVar;
5548 }
5549 
5550 OctaveElementInterface::OctaveElementInterface(OctaveInterface* pInterface, OctaveElement* pElem)
5551 : MBDynInterface(pInterface),
5552 pElem(pElem)
5553 {
5554 
5555 }
5556 
5557 OctaveElementInterface::~OctaveElementInterface()
5558 {
5559 
5560 }
5561 
5562 void
5563 OctaveElementInterface::print(std::ostream& os, bool pr_as_read_syntax) const
5564 {
5565  os << "MBDynElement(" << (pElem ? pElem->GetLabel() : (unsigned)(-1)) << ")" << std::endl;
5566 }
5567 
5568 METHOD_DEFINE(OctaveElementInterface, GetLabel, args, nargout)
5569 {
5570  if (args.length() != 0) {
5571  error("OctaveElement: invalid number of arguments %ld\n"
5572  "expected no arguments", long(args.length()));
5573  return octave_value();
5574  }
5575 
5576  return octave_value(octave_int<unsigned int>(pElem->GetLabel()));
5577 }
5578 
5579 METHOD_DEFINE(OctaveElementInterface, iGetFirstIndex, args, nargout)
5580 {
5581  if (args.length() != 0) {
5582  error("OctaveElement: invalid number of arguments %ld\n"
5583  "expected no arguments", long(args.length()));
5584  return octave_value();
5585  }
5586 
5587  return octave_value(octave_int<integer>(pElem->iGetFirstIndex()));
5588 }
5589 
5590 BEGIN_METHOD_TABLE(OctaveElementInterface, MBDynInterface)
5591  METHOD_DISPATCH(OctaveElementInterface, GetLabel)
5592  METHOD_DISPATCH(OctaveElementInterface, iGetFirstIndex)
5593 END_METHOD_TABLE()
5594 
5595 DEFINE_OCTAVE_ALLOCATOR(OctaveElementInterface);
5596 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA(OctaveElementInterface, "MBDynElement", "MBDynElement");
5597 
5598 }
5599 
5600 #endif // USE_OCTAVE
5601 
5602 bool
5604 {
5605 #ifdef USE_OCTAVE
5606  using namespace oct;
5607 
5608  DriveCallerRead *rf = new OctaveDCR;
5609  if (!SetDriveCallerData("octave", rf)) {
5610  delete rf;
5611  return false;
5612  }
5613 
5614  OctaveTDCR<doublereal>* const rf1D = new OctaveTDCR<doublereal>;
5615  if (!SetDC1D("octave", rf1D)) {
5616  delete rf1D;
5617  return false;
5618  }
5619 
5620  OctaveTDCR<Vec3>* const rf3D = new OctaveTDCR<Vec3>;
5621  if (!SetDC3D("octave", rf3D)) {
5622  delete rf3D;
5623  return false;
5624  }
5625 
5626  OctaveTDCR<Vec6>* const rf6D = new OctaveTDCR<Vec6>;
5627  if (!SetDC6D("octave", rf6D)) {
5628  delete rf6D;
5629  return false;
5630  }
5631 
5632  OctaveTDCR<Mat3x3>* const rf3x3D = new OctaveTDCR<Mat3x3>;
5633  if (!SetDC3x3D("octave", rf3x3D)) {
5634  delete rf3x3D;
5635  return false;
5636  }
5637 
5638  OctaveTDCR<Mat6x6>* const rf6x6D = new OctaveTDCR<Mat6x6>;
5639  if (!SetDC6x6D("octave", rf6x6D)) {
5640  delete rf6x6D;
5641  return false;
5642  }
5643 
5644  rf = new DerivativeDCR;
5645  if (!SetDriveCallerData("derivative", rf)) {
5646  delete rf;
5647  return false;
5648  }
5649 
5650  ScalarFunctionRead* rsf = new OctaveSFR;
5651  if (!SetSF("octave", rsf)) {
5652  delete rsf;
5653  return false;
5654  }
5655 
5657  if (!SetUDE("octave", urf)) {
5658  delete urf;
5659  return false;
5660  }
5661 
5662  ConstitutiveLawRead<doublereal, doublereal> *rfcl1D = new OctaveCLR<doublereal, doublereal>;
5663  if (!SetCL1D("octave", rfcl1D)) {
5664  delete rfcl1D;
5665  return false;
5666  }
5667 
5668  ConstitutiveLawRead<Vec3, Mat3x3> *rfcl3D = new OctaveCLR<Vec3, Mat3x3>;
5669  if (!SetCL3D("octave", rfcl3D)) {
5670  delete rfcl3D;
5671  return false;
5672  }
5673 
5674  ConstitutiveLawRead<Vec6, Mat6x6> *rfcl6D = new OctaveCLR<Vec6, Mat6x6>;
5675  if (!SetCL6D("octave", rfcl6D)) {
5676  delete rfcl6D;
5677  return false;
5678  }
5679 #else
5680  pedantic_cerr("warning: MBDyn has been configured without octave support\n"
5681  "warning: module-octave is not available in this version of MBDyn" << std::endl);
5682 #endif
5683 
5684  // Return true also if USE_OCTAVE is not enabled
5685  // This prevents one assertion to fail in userelem.cc if built as a static module
5686  return true;
5687 }
5688 
5689 #ifndef STATIC_MODULES
5690 
5691 extern "C" int
5692 module_init(const char *module_name, void *pdm, void *php)
5693 {
5694  if (!mbdyn_octave_set()) {
5695  silent_cerr("octave: "
5696  "module_init(" << module_name << ") "
5697  "failed" << std::endl);
5698  return -1;
5699  }
5700 
5701  return 0;
5702 }
5703 
5704 #endif // ! STATIC_MODULES
flag fReadOutput(MBDynParser &HP, const T &t) const
Definition: dataman.h:1064
Mat3x3 GetRotRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1795
const char * sFileName
Definition: parser.h:304
virtual const doublereal & dGetDofValue(int iDof, int iOrder=0) const =0
void PutColIndex(integer iSubCol, integer iCol)
Definition: submat.h:325
bool SetDriveCallerData(const char *name, DriveCallerRead *rf)
Definition: drive_.cc:1324
Definition: mathtyp.h:175
Real GetReal(void) const
Definition: mathp.cc:1228
long int flag
Definition: mbdyn.h:43
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Definition: matvec3.h:98
Definition: node.h:67
TypedValue::Type GetType(void) const
Definition: mathp.cc:1155
virtual void ResizeReset(integer)
Definition: vh.cc:55
#define TRACE(expr)
Definition: linesearch.cc:77
virtual integer GetInt(integer iDefval=0)
Definition: parser.cc:1050
bool Const(void) const
Definition: mathp.cc:1811
FullSubMatrixHandler & SetFull(void)
Definition: submat.h:1168
bool SetDC3x3D(const char *name, TplDriveCallerRead< Mat3x3 > *rf)
virtual TypedValue GetVal(void) const =0
int error(const char *test, int value)
virtual const Vec3 & GetXPrev(void) const
Definition: strnode.h:304
virtual bool GetBool(bool bDefval=false)
Definition: parser.cc:1010
bool SetDC6D(const char *name, TplDriveCallerRead< Vec6 > *rf)
virtual unsigned int iGetPrivDataIdx(const char *s) const
Definition: simentity.cc:142
virtual const char * GetFileName(enum Delims Del=DEFAULTDELIM)
Definition: parsinc.cc:673
void PutCoef(integer iRow, integer iCol, const doublereal &dCoef)
Definition: submat.h:672
static const unsigned int iNumPrivData
Definition: beam.cc:249
virtual const doublereal & dGetDofValuePrev(int iDof, int iOrder=0) const =0
void ResizeReset(integer iNewRow, integer iNewCol)
Definition: submat.cc:1084
#define NO_OP
Definition: myassert.h:74
std::vector< Hint * > Hints
Definition: simentity.h:89
const char *const GetTypeName(void) const
Definition: mathp.cc:1176
virtual void PutItem(integer iSubRow, integer iRow, const doublereal &dCoef)
Definition: submat.h:1445
virtual bool GetYesNoOrBool(bool bDefval=false)
Definition: parser.cc:1038
virtual TypedValue GetValue(const TypedValue &v)
Definition: parser.cc:1004
virtual doublereal dGetP(const doublereal &dVar) const
Definition: drive.cc:499
VM::const_iterator end(void) const
Definition: table.h:65
void func(const T &u, const T &v, const T &w, doublereal e, T &f)
enum @55 order
void PutItem(integer iSubIt, integer iRow, integer iCol, const doublereal &dCoef)
Definition: submat.h:997
Definition: matvec6.h:37
virtual bool bIsDifferentiable(void) const
Definition: drive.h:495
Vec3 GetPosRel(const ReferenceFrame &rf)
Definition: mbpar.cc:1331
virtual ConstLawType::Type GetConstLawType(void) const =0
void SetNullMatrix(void)
Definition: submat.h:1159
Definition: mbdyn.h:76
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
bool SetDC3D(const char *name, TplDriveCallerRead< Vec3 > *rf)
std::ostream & Loadable(void) const
Definition: output.h:506
virtual const Vec3 & GetVPrev(void) const
Definition: strnode.h:316
static int count
Definition: modules.cc:41
Type
Definition: node.h:71
virtual const char * GetStringWithDelims(enum Delims Del=DEFAULTDELIM, bool escape=true)
Definition: parser.cc:1228
Definition: mbdyn.h:77
bool SetDC1D(const char *name, TplDriveCallerRead< doublereal > *rf)
bool SetCL3D(const char *name, ConstitutiveLawRead< Vec3, Mat3x3 > *rf)
virtual integer iGetFirstRowIndex(void) const
Definition: node.cc:82
virtual integer iGetFirstMomentumIndex(void) const =0
virtual integer iGetFirstPositionIndex(void) const
Definition: strnode.h:452
const std::string & GetString(void) const
Definition: mathp.cc:1247
#define ASSERT(expression)
Definition: colamd.c:977
bool SetCL1D(const char *name, ConstitutiveLawRead< doublereal, doublereal > *rf)
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
virtual const Vec3 & GetXCurr(void) const
Definition: strnode.h:310
virtual void ResizeReset(integer, integer)
Definition: submat.cc:182
Definition: except.h:79
virtual DriveCaller * pCopy(void) const =0
static int nodes
virtual doublereal dGetPrivData(unsigned int i) const
Definition: simentity.cc:149
virtual unsigned int iGetNumPrivData(void) const
Definition: simentity.cc:136
bool SetCL6D(const char *name, ConstitutiveLawRead< Vec6, Mat6x6 > *rf)
virtual std::string GetString(const std::string &sDefVal)
Definition: parser.cc:1074
bool SetSF(const std::string &s, const ScalarFunctionRead *rf)
virtual bool IsArg(void)
Definition: parser.cc:807
Definition: elem.h:75
virtual bool IsStringWithDelims(enum Delims Del=DEFAULTDELIM)
Definition: parser.cc:1210
void PutRowIndex(integer iSubRow, integer iRow)
Definition: submat.h:311
bool mbdyn_octave_set(void)
std::ostream & GetLogFile(void) const
Definition: dataman.h:326
bool SetUDE(const std::string &s, UserDefinedElemRead *rude)
Definition: userelem.cc:97
virtual const Vec3 & GetVCurr(void) const
Definition: strnode.h:322
DriveCaller * GetDriveCaller(bool bDeferred=false)
Definition: mbpar.cc:2033
Definition: table.h:43
bool SetDC6x6D(const char *name, TplDriveCallerRead< Mat6x6 > *rf)
virtual const Vec3 & GetXPPCurr(void) const
Definition: strnode.h:334
virtual const Vec3 & GetXPPPrev(void) const
Definition: strnode.h:328
SparseSubMatrixHandler & SetSparse(void)
Definition: submat.h:1178
virtual integer iGetFirstIndex(void) const
Definition: dofown.h:127
const std::string & GetName(void) const
Definition: withlab.cc:68
virtual doublereal Get(const doublereal &d)
Definition: mbpar.cc:2213
Int GetInt(void) const
Definition: mathp.cc:1209
double doublereal
Definition: colamd.c:52
VM::const_iterator begin(void) const
Definition: table.h:64
long int integer
Definition: colamd.c:51
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
void SetVal(const bool &b)
Definition: mathp.cc:1829
static std::stack< const HighParser * > pHP
Definition: parser.cc:598
bool GetBool(void) const
Definition: mathp.cc:1188
unsigned int GetLabel(void) const
Definition: withlab.cc:62
virtual bool bComputeAccelerations(void) const
Definition: strnode.h:433
const char * path
Definition: autopilot.c:141
int module_init(const char *module_name, void *pdm, void *php)
This function registers our user defined element for the math parser.
bool UseText(int out) const
Definition: output.cc:446
virtual doublereal GetReal(const doublereal &dDefval=0.0)
Definition: parser.cc:1056
virtual integer iGetFirstColIndex(void) const
Definition: node.cc:88