tadd simulation structure, more vector constructors, initialize slider values - slidergrid - grid of elastic sliders on a frictional surface
(HTM) git clone git://src.adamsgaard.dk/slidergrid
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
(DIR) commit d71a34ad2042f7d67856a76cc130c204a69a3550
(DIR) parent de268460b46629fdcfcfd155a1699fd3bd69b5c0
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Wed, 16 Mar 2016 12:34:47 -0700
add simulation structure, more vector constructors, initialize slider values
Diffstat:
M Makefile | 2 +-
M main.c | 55 +++++++++++++++++++++----------
A simulation.c | 16 ++++++++++++++++
A simulation.h | 21 +++++++++++++++++++++
M slider.c | 43 +++++++++++++++++++++++++------
M slider.h | 7 +++++--
M vector_math.c | 12 +++++++++++-
M vector_math.h | 2 ++
8 files changed, 129 insertions(+), 29 deletions(-)
---
(DIR) diff --git a/Makefile b/Makefile
t@@ -3,7 +3,7 @@ CFLAGS=-Wall -O3 -g -pg
LDLIBS=-lm
BIN=slidergrid
-$(BIN): main.o slider.o grid.o vector_math.o
+$(BIN): main.o slider.o grid.o vector_math.o simulation.o
$(CC) $(LDLIBS) $^ -o $@
profile: $(BIN)
(DIR) diff --git a/main.c b/main.c
t@@ -2,40 +2,61 @@
#include <stdlib.h>
#include "slider.h"
#include "grid.h"
+#include "simulation.h"
-int main(int argc, char** argv)
+
+simulation setup_simulation()
{
+ // create empty simulation structure
+ simulation sim;
// initialize grid of sliders
int nx = 4;
int ny = 4;
int nz = 4;
- const int N = nx*ny*nz;
- slider* sliders = create_regular_slider_grid(nx, ny, nz, 1.0, 1.0, 1.0);
+ sim.N = nx*ny*nz;
+ sim.sliders = create_regular_slider_grid(nx, ny, nz, 1.0, 1.0, 1.0);
// set slider masses and moments of inertia
int i;
- for (i=0; i<N; i++) {
- sliders[i].mass = 1.;
- sliders[i].moment_of_inertia = 1.;
+ for (i=0; i<sim.N; i++) {
+ sim.sliders[i].mass = 1.;
+ sim.sliders[i].moment_of_inertia = 1.;
}
+}
+
+int main(int argc, char** argv)
+{
+ int i;
+
+ // external function which defines the simulation setup and parameters
+ simulation sim = setup_simulation();
+
+ // set initial values for the sliders
+ for (i=0; i<sim.N; i++)
+ initialize_slider_values(sim.sliders[i]);
- find_and_bond_to_neighbors_n2(sliders, N, 1.5);
+ // create bonds between sliders and update neighbor lists
+ find_and_bond_to_neighbors_n2(sim.sliders, sim.N, sim.bond_length_limit);
+
+ // determine time step length
+ sim.dt = find_time_step(sim.sliders, sim.N);
+
+ // check simulation parameters for missing or wrong parameter values
+ if (check_simulation(sim)) {
+ fprintf(stderr, "Error: Simulation configuration invalid.\n");
+ return EXIT_FAILURE;
+ }
// temporal loop
- Float t;
- const Float t_end = 1.;
- //Float dt = calculate_time_step();
- Float dt = 1.0;
- long int iteration = 0;
- for (t=0.0; t<=t_end; t+=dt) {
- for (i=0; i<N; i++) {
- update_kinematics(sliders[i], dt, iteration);
+ for (sim.t=0.0; sim.t<=sim.t_end; sim.t+=sim.dt) {
+ for (i=0; i<sim.N; i++) {
+ update_kinematics(sim.sliders[i], sim.dt, sim.iteration);
}
- iteration++;
+ sim.iteration++;
}
// end program
- free(sliders);
+ free(sim.sliders);
return EXIT_SUCCESS;
}
(DIR) diff --git a/simulation.c b/simulation.c
t@@ -0,0 +1,16 @@
+#include "slider.h"
+#include "simulation.h"
+
+int check_simulation_values(const simulation sim)
+{
+
+
+}
+
+Float find_time_step(const slider* sliders, const int N)
+{
+ Float largest_dt = 0.;
+
+
+ return largest_dt;
+}
(DIR) diff --git a/simulation.h b/simulation.h
t@@ -0,0 +1,21 @@
+#ifndef SIMULATION_H_
+#define SIMULATION_H_
+
+typedef struct {
+ slider* sliders;
+ int N;
+ Float t;
+ Float t_end;
+ Float dt;
+ long int iteration;
+ Float bond_length_limit;
+
+} simulation;
+
+Float find_time_step(const slider* sim.sliders, const int sim.N);
+int check_simulation(const simulation sim);
+
+// user-defined function which sets up the simulation
+simulation setup_simulation();
+
+#endif
(DIR) diff --git a/slider.c b/slider.c
t@@ -8,6 +8,28 @@ void print_slider_position(slider s)
printf("%f\t%f\t%f\n", s.pos.x, s.pos.y, s.pos.z);
}
+void initialize_slider_values(slider s)
+{
+ s.force = zeroes_float3();
+ s.torque = zeroes_float3();
+
+ // define all entries in neighbor list as empty
+ int i;
+ for (i=0; i<MAX_NEIGHBORS; i++) {
+ s.neighbors[i] = -1;
+ s.neighbor_distance[i] = zeros_float3();
+ s.neighbor_relative_distance_displacement[i] = zeros_float3();
+ s.neighbor_relative_distance_velocity[i] = zeros_float3();
+ }
+
+}
+
+void check_slider_values(slider s)
+{
+
+
+}
+
/* Explicit temporal integration scheme based on three-term Taylor expansion.
* Truncation error O(dt^4) for positions, O(dt^3) for velocities. Acceleration
* change is approximated by backwards differences. */
t@@ -18,8 +40,8 @@ void update_kinematics(slider s, Float dt, long int iteration)
Float3 acc0, angacc0;
if (iteration == 0) {
- acc0 = make_float3(0., 0., 0.);
- angacc0 = make_float3(0., 0., 0.);
+ acc0 = zeroes_float3();
+ angacc0 = zeroes_float3();
} else {
acc0 = s.acc;
angacc0 = s.angacc;
t@@ -97,7 +119,7 @@ void slider_displacement(slider s1, const slider s2,
// read previous inter-slider distance vector
Float3 dist0;
if (iteration == 0)
- dist0 = make_float3(0., 0., 0.);
+ dist0 = zeroes_float3();
else
dist0 = s1.neighbor_distance[idx_neighbor];
t@@ -128,19 +150,24 @@ void slider_interaction(slider s1, const slider s2, const int idx_neighbor)
// determine the bond stiffness from the harmonic mean.
// differs from Potyondy & Stack 2004, eq. 6.
- const Float spring_stiffness =
- 2.*s1.spring_stiffness*s2.spring_stiffness/
- (s1.spring_stiffness + s2.spring_stiffness);
+ const Float bond_parallel_stiffness =
+ 2.*s1.bond_parallel_stiffness*s2.bond_parallel_stiffness/
+ (s1.bond_parallel_stiffness + s2.bond_parallel_stiffness);
+
+ // determine the bond viscosity from the harmonic mean.
+ const Float bond_parallel_viscosity =
+ 2.*s1.bond_parallel_viscosity*s2.bond_parallel_viscosity/
+ (s1.bond_parallel_viscosity + s2.bond_parallel_viscosity);
// bond-parallel elasticity
const Float3 f_n_elastic =
- multiply_scalar_float3(bond_normal_stiffness,
+ multiply_scalar_float3(bond_parallel_stiffness,
s1.neighbor_relative_distance_displacement[idx_neighbor]);
// bond-parallel viscosity
const Float3 f_n_viscous =
- multiply_scalar_float3(bond_normal_viscosity,
+ multiply_scalar_float3(bond_parallel_viscosity,
s1.neighbor_relative_distance_displacement[idx_neighbor]);
// bond-parallel force, counteracts tension and compression
(DIR) diff --git a/slider.h b/slider.h
t@@ -30,10 +30,11 @@ typedef struct {
Float moment_of_inertia;
// inter-slider contact model parameters
- Float bond_normal_stiffness; // Hookean elastic stiffness [N/m]
- Float bond_normal_viscosity; // viscosity [N/(m*s)]
+ Float bond_parallel_stiffness; // Hookean elastic stiffness [N/m]
+ Float bond_parallel_viscosity; // viscosity [N/(m*s)]
// slider-surface contact model parameters
+ // should this value belong to the lower surface instead?
Float coulomb_friction; // Coulomb frictional value [-]
// The uniquely identifying indexes of the slider neighbors which are bonded
t@@ -50,6 +51,8 @@ typedef struct {
void print_slider_position(slider s);
+void initialize_slider_values(slider s);
+//void check_slider_values(slider s);
void update_kinematics(slider s, Float dt, long int iteration);
(DIR) diff --git a/vector_math.c b/vector_math.c
t@@ -2,13 +2,23 @@
#include "typedefs.h"
-// constructor
+// constructors
inline Float3 make_float3(Float x, Float y, Float z)
{
Float3 v = {.x = x, .y = y, .z = z};
return v;
}
+inline Float3 zeroes_float3()
+{
+ return make_float3(0., 0., 0.);
+}
+
+inline Float3 ones_float3()
+{
+ return make_float3(1., 1., 1.);
+}
+
// single-vector operations
inline Float3 copy_float3(Float3 v)
(DIR) diff --git a/vector_math.h b/vector_math.h
t@@ -2,6 +2,8 @@
// constructor
Float3 make_float3(Float x, Float y, Float z);
+Float3 zeroes_float3();
+Float3 ones_float3();
// single-vector operations
Float3 copy_float3(Float3 v);