tVarious bugfixes - simple_DEM - a simple 2D Discrete Element Method code for educational purposes
(HTM) git clone git://src.adamsgaard.dk/simple_DEM
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
(DIR) commit 6890307b5b56a2d9ef0a09395342480f78a2e002
(DIR) parent 7e6917b00c4f2b71ff79a17de701b808fee18442
(HTM) Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Sun, 14 Oct 2012 21:21:29 +0200
Various bugfixes
Diffstat:
M Makefile | 28 ++++++++++++++--------------
M global_properties.h | 29 ++++++++++++++---------------
M grains.c | 17 +++++++++++------
M header.h | 25 ++++++++++++-------------
M initialize.c | 11 ++++++-----
M main.c | 25 +++++++++----------------
M vtk_export.c | 29 ++++++++++++++++-------------
M walls.c | 2 +-
8 files changed, 83 insertions(+), 83 deletions(-)
---
(DIR) diff --git a/Makefile b/Makefile
t@@ -1,19 +1,19 @@
-CC=g++
-CFLAGS=-c -Wall -O2
-LDFLAGS=
-#CFLAGS=-c -Wall -O2 -fopenmp -g -G
+CC=gcc
+CFLAGS=-Wall -O2
+LDFLAGS=-lm
+#CFLAGS=-Wall -O2 -fopenmp -g -G
#LDFLAGS=-fopenmp
-SOURCES=main.c
-OBJECTS=$(SOURCES:.c=.o)
-EXECUTABLE=simple_DEM
+SRC=grains.c initialize.c main.c vtk_export.c walls.c
+OBJ=$(SRC:.c=.o)
+BIN=simple_DEM
+DEPS=header.h global_properties.h
-all: $(SOURCES) $(EXECUTABLE)
-
-$(EXECUTABLE): $(OBJECTS)
- $(CC) $(LDFLAGS) $(OBJECTS) -o $@
+$(BIN): $(OBJ)
+ $(CC) $(LDFLAGS) $^ -o $@
-.c.o:
- $(CC) $(CFLAGS) $< -o $@
+%.o: %.c $(DEPS)
+ $(CC) -c -o $@ $< $(CFLAGS)
clean:
- rm -f $(EXECUTABLE) *.o output/*
+ $(RM) $(BIN) $(OBJ) output/*
+
(DIR) diff --git a/global_properties.h b/global_properties.h
t@@ -2,29 +2,28 @@
#define GLOBAL_PROPERTIES_H_
/* Properties of sample */
-const int ng = 5000; /* Number of grains */
+static const int np = 1000; /* Number of particles */
/* Size distribution */
-const double rmin = 1e-3; /* Min. radius */
-const double rmax = 2e-3; /* Max. radius */
+static const double rmin = 1.0e-3; /* Min. radius */
+static const double rmax = 2.0e-3; /* Max. radius */
/* Properties of grains */
-const double kn = 1e5; /* Normal stiffness */
-const double nu = 30; /* Normal damping */
-const double rho = 1000; /* Density of the grains */
-const double mu = 0.5; /* Sliding friction */
-const double kt = kn; /* Tangential stiffness */
+static const double kn = 1.0e5; /* Normal stiffness */
+static const double nu = 30.0; /* Normal damping */
+static const double rho = 1000.0; /* Density of the grains */
+static const double mu = 0.5; /* Sliding friction */
+static const double kt = 1.0e5; /* Tangential stiffness */
/* Temporal variables */
-const double dt = 1e-4; /* Time step length */
-const int maxStep = 3000; /* Number of steps */
-const int fileInterval = 20; /* No. of steps between output */
+static const double dt = 1.0e-4; /* Time step lenpth */
+static const int maxStep = 3000; /* Number of steps */
+static const int fileInterval = 20; /* No. of steps between output */
/* Physical constants */
-const double grav = 9.80; /* Gravity */
+static const double grav = 9.80; /* Gravity magnitude */
-/* Number of grains along the width in the initial configuration */
-const int ngw = 100;
+/* Number of particles along the width in the initial configuration */
+static const int npw = 100;
#endif
-
(DIR) diff --git a/grains.c b/grains.c
t@@ -1,9 +1,14 @@
+#include <math.h>
+#include "header.h"
+#include "global_properties.h"
+
+
void prediction(grain* g)
{
int i;
#pragma omp parallel for shared(g) private (i)
- for (i = 0; i < ng; i++) {
+ for (i = 0; i < np; i++) {
/* Predict positions and velocities */
g[i].x += dt * g[i].vx + 0.5 * dt * dt * g[i].ax;
g[i].y += dt * g[i].vy + 0.5 * dt * dt * g[i].ay;
t@@ -37,7 +42,7 @@ void interparticle_force(grain* g, int a, int b)
if (dn < 0.0) { /* Contact */
double xn, yn, vn, fn; /* Normal components */
- double xt, yt, vt, ft; /* Tangential components */
+ double xt, yt, vt, ft; /* Tanpential components */
/* Local axes */
xn = x_ab / dist;
yn = y_ab / dist;
t@@ -84,10 +89,10 @@ void interact_grains(grain* g)
int a, b;
#pragma omp parallel for shared(g) private (a,b)
/* Loop through particle a */
- for (a = 0; a < ng; a++) {
+ for (a = 0; a < np; a++) {
/* Loop through particle b */
- for (b = 0; b < ng; b++) {
+ for (b = 0; b < np; b++) {
interparticle_force(g, a, b);
}
t@@ -99,7 +104,7 @@ void update_acc(grain* g)
{
int i;
#pragma omp parallel for shared(g) private (i)
- for (i = 0; i < ng; i++) {
+ for (i = 0; i < np; i++) {
g[i].ax = g[i].fx / g[i].m;
g[i].ay = g[i].fy / g[i].m - grav;
g[i].ath = g[i].fth / g[i].I;
t@@ -110,7 +115,7 @@ void correction(grain* g)
{
int i;
#pragma omp parallel for shared(g) private (i)
- for (i = 0; i < ng; i++) {
+ for (i = 0; i < np; i++) {
g[i].vx += 0.5 * dt * g[i].ax;
g[i].vy += 0.5 * dt * g[i].ay;
g[i].vth += 0.5 * dt * g[i].ath;
(DIR) diff --git a/header.h b/header.h
t@@ -2,8 +2,7 @@
#define HEADER_H_
/* Structure declarations */
-
-struct grain
+typedef struct
{
double m; /* Mass */
double R; /* Radius */
t@@ -13,33 +12,33 @@ struct grain
double ax, ay, ath; /* Acceleration */
double fx, fy, fth; /* Sum of forces, decomposed */
double p; /* Pressure */
-};
+} grain;
/* Prototype functions */
-/* initialize.cpp */
+/* initialize.c */
void triangular_grid(grain* g);
-/* grains.cpp */
+/* grains.c */
void prediction(grain* g);
void interparticle_force(grain* g, int a, int b);
void interact_grains(grain* g);
void update_acc(grain* g);
void correction(grain* g);
-/* walls.cpp */
-void compute_force_upper_wall(int i, grain* g, double wfy, double wupy);
-void compute_force_lower_wall(int i, grain* g, double wfy, double wloy);
-void compute_force_left_wall(int i, grain* g, double wfy, double wlex);
-void compute_force_right_wall(int i, grain* g, double wfy, double wrix);
-
+/* walls.c */
+double compute_force_upper_wall(int i, grain* g, double wup);
+double compute_force_lower_wall(int i, grain* g, double wdown);
+double compute_force_left_wall(int i, grain* g, double wleft);
+double compute_force_right_wall(int i, grain* g, double wright);
+void interact_walls(grain* g, double wleft, double wright, double wup, double wdown,
+ double* wp_up, double* wp_down, double* wp_left, double* wp_right);
-/* vtk_export.cpp */
+/* vtk_export.c */
int vtk_export_grains(grain* g, int numfile);
/*int vtk_export_walls(int numfile, double wlex, double wrix, double wloy, double wupy);*/
int vtk_export_forces(grain* g, int numfile);
-
#endif
(DIR) diff --git a/initialize.c b/initialize.c
t@@ -5,18 +5,19 @@
void triangular_grid(grain* g)
{
+ int i;
/* Initialize grain properties */
- for (int i = 0; i < ng; i++) {
+ for (i = 0; i < np; ++i) {
g[i].R = (rand() / (double)RAND_MAX) * (rmax - rmin) + rmin;
g[i].m = M_PI * rho * g[i].R * g[i].R;
g[i].I = 0.5 * g[i].m * g[i].R * g[i].R;
g[i].p = 0.0;
}
- /* Initialize grain positions in a triangular grid */
- for (int i = 0; i < ng; i++) {
- int column = i%ngw;
- int row = i/ngw;
+ /* Initialize grain positions in a trianpular grid */
+ for (i = 0; i < np; ++i) {
+ int column = i%npw;
+ int row = i/npw;
if (row%2 == 0) /* Even row */
g[i].x = rmax + 2*column*rmax;
(DIR) diff --git a/main.c b/main.c
t@@ -7,34 +7,26 @@
/* Global and constant simulation properties */
#include "global_properties.h"
-/* Functions for exporting data to VTK formats */
-#include "vtk_export.c"
-#include "initialize.c"
-#include "grains.c"
-#include "walls.c"
-
-
-
-int main()
+int main(int argc, char* argv[])
{
printf("\n## simple_DEM ##\n");
- printf("Particles: %d\n");
- printf("maxStep: %d\n");
+ printf("Particles: %d\n", np);
+ printf("maxStep: %d\n", maxStep);
double time = 0.0; /* Time at simulation start */
/* Allocate memory */
- grain* g = new grain[ng]; /* Grain structure */
+ grain g[np]; /* Grain structure */
/* Compute simulation domain dimensions */
double wleft = 0.0; /* Left wall */
- double wright = (ngw+1)*2*rmax; /* Right wall */
+ double wright = (npw+1)*2*rmax; /* Right wall */
double wdown = 0.0; /* Lower wall */
- double wup = (ng/ngw+1)*2*rmax; /* Upper wall */
+ double wup = (np/npw+1)*2*rmax; /* Upper wall */
/* Variables for pressures on walls */
double wp_up, wp_down, wp_left, wp_right;
t@@ -45,7 +37,8 @@ int main()
/* Main time loop */
- for (int step = 0; step < maxStep; step++) {
+ int step;
+ for (step = 0; step < maxStep; ++step) {
time += dt; /* Update current time */
t@@ -74,7 +67,7 @@ int main()
/* Free dynamically allocated memory */
- delete[] g;
+ /*free(g);*/
printf("\nSimulation ended without errors.\n");
(DIR) diff --git a/vtk_export.c b/vtk_export.c
t@@ -5,14 +5,15 @@
int vtk_export_grains(grain* g, int numfile)
{
FILE* fout;
+ int i;
char filename[25]; /* File name */
sprintf(filename, "output/grains%04d.vtk", numfile);
if ((fout = fopen(filename, "wt")) == NULL) {
- printf("vtk_export error, cannot open ")
- printf(filename);
- printf("!\n");
+ fprintf(stderr, "vtk_export error, cannot open ");
+ fprintf(stderr, filename);
+ fprintf(stderr, "!\n");
return 1;
}
t@@ -26,31 +27,31 @@ int vtk_export_grains(grain* g, int numfile)
fprintf(fout, "ASCII\nDATASET UNSTRUCTURED_GRID\n");
/* Grain coordinates */
- fprintf(fout, "POINTS %d FLOAT\n", ng);
- for (int i = 0; i < ng; i++)
+ fprintf(fout, "POINTS %d FLOAT\n", np);
+ for (i = 0; i < np; i++)
fprintf(fout, "%f %f 0.0\n", g[i].x, g[i].y);
- fprintf(fout, "POINT_DATA %d\n", ng);
+ fprintf(fout, "POINT_DATA %d\n", np);
/* Grain radii */
fprintf(fout, "VECTORS Radius float\n");
- for (int i = 0; i < ng; i++)
+ for (i = 0; i < np; i++)
fprintf(fout, "%f 0.0 0.0\n", g[i].R);
/* Grain velocities */
fprintf(fout, "VECTORS Velocity float\n");
- for (int i = 0; i < ng; i++)
+ for (i = 0; i < np; i++)
fprintf(fout, "%f %f 0.0\n", g[i].vx, g[i].vy);
/* Pressure */
fprintf(fout, "SCALARS Pressure float 1\n");
fprintf(fout, "LOOKUP_TABLE default\n");
- for (int i = 0; i < ng; i++)
+ for (i = 0; i < np; i++)
fprintf(fout, "%e\n", g[i].p);
- /* Angular velocity */
- fprintf(fout, "SCALARS Angvel float 1\n");
+ /* Anpular velocity */
+ fprintf(fout, "SCALARS Anpvel float 1\n");
fprintf(fout, "LOOKUP_TABLE default\n");
- for (int i = 0; i < ng; i++)
+ for (i = 0; i < np; i++)
fprintf(fout, "%e\n", g[i].vth);
fclose(fout);
t@@ -65,7 +66,9 @@ int vtk_export_forces(grain* g, int numfile)
sprintf(filename, "output/forces%04d.vtk", numfile);
if ((fout = fopen(filename, "wt")) == NULL) {
- std::cout << "vtk_export error, cannot open " << filename << "!\n";
+ fprintf(stderr,"vtk_export error, cannot open ");
+ fprintf(stderr, filename);
+ fprintf(stderr, "!\n");
return 1;
}
(DIR) diff --git a/walls.c b/walls.c
t@@ -96,7 +96,7 @@ void interact_walls(grain* g, double wleft, double wright, double wup, double wd
int i;
#pragma omp parallel for shared(g, u, d, r, l) private (i)
- for (i = 0; i < ng; i++) {
+ for (i = 0; i < np; i++) {
u += compute_force_upper_wall(i, g, wup);
d += compute_force_lower_wall(i, g, wdown);
r += compute_force_right_wall(i, g, wright);