Hi,
I have a function called Verlet_integration(), that requires a counter as an input variable, along with some other variables. This counter is an integer called 'count_nodes' and the value of it is 96.
Within that function, I call for five other functions in this order: Positions(), Int_velocities(), Forces(), Accelerations() and Velocities(). All of these functions, except for Forces(), require the counter as an input argument.
However, in between the functions Forces() and Accelerations(), the variable 'count_nodes' for some reason becomes 1178657792 - always the same number.
If I disable Forces() everything works fine, so it has to be somewhere within that function. But the weird thing is if I print the value at the end of Forces(), just before its return, the counter is still 96 - see code below.
I cannot find the mistake I have made. I tried comparing the current version with an older but working version using UltraCompare, but I could not find the mistake and I don't know how to find the problem.
Hopefully someone can give me some directions to track down the problem.
/*******************************************************************
*
* Function : Verlet_integration()
* Purpose : Output node-coördinates to XML-file
* External file : lib_time.c
*
*******************************************************************/
int Verlet_integration (
Node* Nodes ,
Link* Links ,
Head* Paul ,
int count_nodes ,
int count_links )
{
float iv_x[count_nodes];
float iv_y[count_nodes];
// Step 1 : Calculate new positions
Positions( Nodes, Paul, count_nodes );
// Step 2 : Intermediate velocities
Int_velocities( Nodes, Paul, count_nodes, iv_x, iv_y );
// Step 3 : Forces (pressure, springs, drag)
if ( Forces( Nodes, Links, Paul, count_nodes, count_links ) != 0 )
{
return 1;
}
printf("%i\n", count_nodes); // Returns 1178657792
// Step 4 : Accelerations
Accelerations( Nodes, Paul, count_nodes );
// Step 5 : Update velocities
Velocities( Nodes, Paul, count_nodes, iv_x, iv_y );
return 0;
}
/*******************************************************************
*
* Function : Forces
* Purpose : Forces on nodes
*
*******************************************************************/
int Forces (
Node* Nodes ,
Link* Links ,
Head* Paul ,
int count_nodes ,
int count_links )
{
float rho = RHO;
float vsound = V_SOUND;
int node1;
int node2;
float k; // Stiffness
float perm; // Permeability
float length; // Length of link
float strain; // Strain
float rx; // Length-vector
float ry;
float tx; // Tangential vector
float ty;
float nx; // Normal vector
float ny;
float fp_x; // Force-vector due to pressure
float fp_y;
float fs_x; // Springforce vector
float fs_y;
float avx; // Average velocity of link
float avy;
float fd_x; // Drag force due to pressure drag
float fd_y;
float sum_x1; // Sumvector for Node 1
float sum_y1;
float sum_x2; // Sumvector for Node 2
float sum_y2;
Force Fc[count_nodes+1]; // Forces due to impact
int j;
Fc[count_nodes+1].x = 12345; // Set check-number
// Turn this off for testing
CheckImpact(Nodes, Paul, Fc, count_nodes);
// Safety check
if ( Fc[count_nodes+1].x != 12345 )
{
printf("Bounds exceeded for Fc.\n");
return 1;
}
for ( j = 0; j < count_links; j++ )
{
node1 = Links[j].node1;
node2 = Links[j].node2;
k = Links[j].k;
perm = Links[j].perm;
rx = Nodes[node2].r_x - Nodes[node1].r_x; // Vector in link-direction
ry = Nodes[node2].r_y - Nodes[node1].r_y;
length = sqrt( pow(rx,2) + pow(ry,2) ); // Length of link
tx = rx / length; // Tangential vector
ty = ry / length;
nx = ty; // Normal vector n = R * t with R = [0 1; -1 0]
ny = -tx;
// Return 1 if vectors > 1 (impossible)
if ( sqrt( pow(tx,2) + pow(ty, 2) ) > 1.001 || sqrt( pow(nx,2) + pow(ny, 2) ) > 1.001)
{
printf("[%f %f] [%f %f]\n", tx, ty, nx, ny);
printf("%f %f\n", sqrt( pow(tx,2) + pow(ty, 2) ), sqrt( pow(nx,2) + pow(ny, 2) ) );
return 1;
}
// Pressure
fp_x = 0.5 * 4000 * length * nx;
fp_y = 0.5 * 4000 * length * ny;
//fp_x = 0;
//fp_y = 0;
// Springs
// F = k * strain * t
strain = ( length / Links[j].l0 ) - 1;
fs_x = k * strain * tx;
fs_y = k * strain * ty;
//fs_x = 0;
//fs_y = 0;
// Drag forces
// Fd = rho * vsound * v * L * n
// Turn off for testing
avx = 0.5 * (Nodes[node1].v_x + Nodes[node2].v_x) * nx;
avy = 0.5 * (Nodes[node1].v_y + Nodes[node2].v_y) * ny;
fd_x = - 0.5 * rho * vsound * length * avx * nx;
fd_y = - 0.5 * rho * vsound * length * avy * ny;
//fd_x = 0;
//fd_y = 0;
// Sum of forces
sum_x1 = (1 - perm)*(fp_x + fd_x) + fs_x + Fc[node1].x;
sum_y1 = (1 - perm)*(fp_y + fd_y) + fs_y + Fc[node1].y;
sum_x2 = (1 - perm)*(fp_x + fd_x) - fs_x + Fc[node2].x;
sum_y2 = (1 - perm)*(fp_y + fd_y) - fs_y + Fc[node2].y;
// Store sum of forces to linked nodes
Nodes[node1].f_x = Nodes[node1].f_x + sum_x1;
Nodes[node1].f_y = Nodes[node1].f_y + sum_y1;
Nodes[node2].f_x = Nodes[node2].f_x + sum_x2;
Nodes[node2].f_y = Nodes[node2].f_y + sum_y2;
}
printf("%i\n", count_nodes); // Returns 96, as it is supposed to
return 0;
}
