rk4 vectors one allocation

MBsysC/mbs_common/mbs_module/mbs_dirdyn.c

@@ -81,6 +81,9 @@ MbsDirdyn* mbs_new_dirdyn_aux(MbsData* mbs_data, MbsAux* mbs_aux)
     // initial the saving counter
     dirdyn->savePeriodCounter = 0;
 
+    // Runge Kutta 4
+    dirdyn->rk4 = (MbsDirdynRK4*) malloc(sizeof(MbsDirdynRK4));
+
     // adaptive time step
     dopr = (MbsDirdynDopri5*) malloc(sizeof(MbsDirdynDopri5));
     dopr->flag_save = 0;
@@ -95,6 +98,7 @@ MbsDirdyn* mbs_new_dirdyn_aux(MbsData* mbs_data, MbsAux* mbs_aux)
 void mbs_delete_dirdyn(MbsDirdyn* dirdyn, MbsData* mbs_data)
 {
     free(dirdyn->dopri5);
+    free(dirdyn->rk4);
     free(dirdyn->options->dopri5);
     free(dirdyn->options);
     freeMbsAux(dirdyn->mbs_aux, mbs_data);
@@ -166,6 +170,12 @@ void mbs_dirdyn_init(MbsDirdyn* dd, MbsData* mbs_data)
         dd->dopri5->dydt_save = (double*) malloc(dd->nState*sizeof(double));
         dd->dopri5->solout_last_t = mbs_data->t0 - mbs_data->dt0;
     }
+    else
+    {
+        dd->rk4->yt  = (double*) malloc(dd->nState*sizeof(double));
+        dd->rk4->dyt = (double*) malloc(dd->nState*sizeof(double));
+        dd->rk4->dym = (double*) malloc(dd->nState*sizeof(double));
+    }
 
     // real-time modules activated
     #ifdef REAL_TIME
@@ -384,6 +394,12 @@ void mbs_dirdyn_finish(MbsDirdyn* dd, MbsData* mbs_data){
         }
         free(dd->dopri5->alloc_tab);
     }
+    else
+    {
+        free(dd->rk4->yt);
+        free(dd->rk4->dyt);
+        free(dd->rk4->dym);
+    }
 
     free(dd->y);
     free(dd->yd);

MBsysC/mbs_common/mbs_numerics/rk4.c

@@ -17,9 +17,9 @@ void rk4(double y[], double dydx[], int n, double x, double h, double yout[],
     int i;
     double xh,hh,h6,*dym,*dyt,*yt;
 
-    dym=dvector(1,n);
-    dyt=dvector(1,n);
-    yt=dvector(1,n);
+    dym = dd->rk4->dym;
+    dyt = dd->rk4->dyt;
+    yt  = dd->rk4->yt;
 
     hh=h*0.5;
     h6=h/6.0;
@@ -49,8 +49,4 @@ void rk4(double y[], double dydx[], int n, double x, double h, double yout[],
     // Accumulate increments with proper weights
     for (i=0;i<n;i++)
         yout[i]=y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);
-
-    free_dvector(yt,1,n);
-    free_dvector(dyt,1,n);
-    free_dvector(dym,1,n);
 }

MBsysC/mbs_common/mbs_struct/mbs_dirdyn_struct.h

@@ -102,6 +102,19 @@ typedef struct MbsDirdynDopri5
 
 } MbsDirdynDopri5;
 
+ /**
+  * Runge-Kutta 4 structure for dirdyn
+  */
+typedef struct MbsDirdynRK4
+{
+    double *yt; ///< allocated state vector
+
+    double *dyt; ///< allocated derivative vector
+
+    double *dym; ///< allocated saved derivative vector
+
+} MbsDirdynRK4;
+
  /**
   * General structure of the direct dynamic module
   */
@@ -133,6 +146,7 @@ typedef struct MbsDirdyn
     /** counter for checking wheter results must be saved */
     int savePeriodCounter;
 
+    MbsDirdynRK4 *rk4; ///< Runge-Kutta 4 structure for dirdyn
     MbsDirdynDopri5 *dopri5; ///< Dopri5 structure for dirdyn
 
 } MbsDirdyn;