tFirst commit - simple_DEM - a simple 2D Discrete Element Method code for educational purposes
HTML git clone git://src.adamsgaard.dk/simple_DEM
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
DIR commit a17213846a2df2ba52051d8f9f9effcd8b58636d
HTML Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Sun, 14 Oct 2012 08:42:39 +0200
First commit
Diffstat:
A Makefile | 19 +++++++++++++++++++
A README.rst | 17 +++++++++++++++++
A global_properties.h | 30 ++++++++++++++++++++++++++++++
A grains.c | 119 +++++++++++++++++++++++++++++++
A header.h | 44 +++++++++++++++++++++++++++++++
A initialize.c | 28 ++++++++++++++++++++++++++++
A main.c | 86 ++++++++++++++++++++++++++++++
A output/blank.txt | 0
A vtk_export.c | 134 +++++++++++++++++++++++++++++++
A walls.c | 109 +++++++++++++++++++++++++++++++
10 files changed, 586 insertions(+), 0 deletions(-)
---
DIR diff --git a/Makefile b/Makefile
t@@ -0,0 +1,19 @@
+CC=g++
+CFLAGS=-c -Wall -O2
+LDFLAGS=
+#CFLAGS=-c -Wall -O2 -fopenmp -g -G
+#LDFLAGS=-fopenmp
+SOURCES=main.c
+OBJECTS=$(SOURCES:.c=.o)
+EXECUTABLE=simple_DEM
+
+all: $(SOURCES) $(EXECUTABLE)
+
+$(EXECUTABLE): $(OBJECTS)
+ $(CC) $(LDFLAGS) $(OBJECTS) -o $@
+
+.c.o:
+ $(CC) $(CFLAGS) $< -o $@
+
+clean:
+ rm -f $(EXECUTABLE) *.o output/*
DIR diff --git a/README.rst b/README.rst
t@@ -0,0 +1,17 @@
+simple_DEM
+==========
+
+A basic CPU-based discrete element method algorithm, created for comparison against more sophisticated codes.
+
+Requirements
+------------
+- GNU Make
+- GCC
+
+Obtaining
+---------
+`git clone https://github.com/anders-dc/simple_DEM.git`
+
+Building
+--------
+`make`
DIR diff --git a/global_properties.h b/global_properties.h
t@@ -0,0 +1,30 @@
+#ifndef GLOBAL_CONSTANTS_H_
+#define GLOBAL_CONSTANTS_H_
+
+// Properties of sample
+const int ng = 5000; // Number of grains
+
+// Size distribution
+const double rmin = 1e-3; // Min. radius
+const double rmax = 2e-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
+
+// 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
+
+// Physical constants
+const double grav = 9.80; // Gravity
+
+// Number of grains along the width in the initial configuration
+const int ngw = 100;
+
+#endif
+
DIR diff --git a/grains.c b/grains.c
t@@ -0,0 +1,119 @@
+void prediction(grain* g)
+{
+ int i;
+
+ #pragma omp parallel for shared(g) private (i)
+ for (i = 0; i < ng; 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;
+ g[i].vx += 0.5 * dt * g[i].ax;
+ g[i].vy += 0.5 * dt * g[i].ay;
+
+ g[i].th += dt * g[i].vth + 0.5 * dt * dt * g[i].ath;
+ g[i].vth += 0.5 * dt * g[i].ath;
+
+ // Zero forces
+ g[i].fx = 0.0;
+ g[i].fy = 0.0;
+ g[i].fth = 0.0;
+ g[i].p = 0.0;
+ }
+}
+
+void interparticle_force(grain* g, int a, int b)
+{
+ if (a > b) { // Use Newtons 3rd law to find both forces at once
+
+ // Particle center coordinate component differences
+ double x_ab = g[a].x - g[b].x;
+ double y_ab = g[a].y - g[b].y;
+
+ // Particle center distance
+ double dist = sqrt(x_ab*x_ab + y_ab*y_ab);
+
+ // Size of overlap
+ double dn = dist - (g[a].R + g[b].R);
+
+ if (dn < 0.0) { // Contact
+ double xn, yn, vn, fn; // Normal components
+ double xt, yt, vt, ft; // Tangential components
+ // Local axes
+ xn = x_ab / dist;
+ yn = y_ab / dist;
+ xt = -yn;
+ yt = xn;
+
+ // Compute the velocity of the contact
+ double vx_ab = g[a].vx - g[b].vy;
+ double vy_ab = g[a].vy - g[b].vy;
+ vn = vx_ab*xn + vy_ab*yn;
+ vt = vx_ab*xt + vy_ab*yt - (g[a].R*g[a].vth + g[b].R*g[b].vth);
+
+ // Compute force in local axes
+ fn = -kn * dn - nu * vn;
+
+ // Rotation
+ if (fn < 0)
+ fn = 0.0;
+ ft = fabs(kt * vt);
+ if (ft > mu*fn) // Coefficient of friction
+ ft = mu*fn;
+ if (vt > 0)
+ ft = -ft;
+
+ // Calculate sum of forces on a and b in global coordinates
+ g[a].fx += fn * xn;
+ g[a].fy += fn * yn;
+ g[a].fth += -ft*g[a].R;
+ g[a].p += fn;
+ g[b].fx -= fn * xn;
+ g[b].fy -= fn * yn;
+ g[b].p += fn;
+ g[b].fth += -ft*g[b].R;
+
+ }
+
+ } else {
+ return;
+ }
+}
+
+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++) {
+
+ // Loop through particle b
+ for (b = 0; b < ng; b++) {
+ interparticle_force(g, a, b);
+ }
+
+ }
+
+}
+
+void update_acc(grain* g)
+{
+ int i;
+ #pragma omp parallel for shared(g) private (i)
+ for (i = 0; i < ng; 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;
+ }
+}
+
+void correction(grain* g)
+{
+ int i;
+ #pragma omp parallel for shared(g) private (i)
+ for (i = 0; i < ng; 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@@ -0,0 +1,44 @@
+#ifndef HEADER_H_
+#define HEADER_H_
+
+//// Structure declarations
+struct grain
+{
+ double m; // Mass
+ double R; // Radius
+ double I; // Inertia
+ double x, y, th; // Position
+ double vx, vy, vth; // Velocities
+ double ax, ay, ath; // Acceleration
+ double fx, fy, fth; // Sum of forces, decomposed
+ double p; // Pressure
+};
+
+
+
+//// Prototype functions
+
+// initialize.cpp
+void triangular_grid(grain* g);
+
+// grains.cpp
+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);
+
+
+// vtk_export.cpp
+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@@ -0,0 +1,28 @@
+#include <stdlib.h>
+#include <math.h>
+#include "header.h"
+#include "global_properties.h"
+
+void triangular_grid(grain* g)
+{
+ // Initialize grain properties
+ for (int i = 0; i < ng; 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;
+
+ if (row%2 == 0) // Even row
+ g[i].x = rmax + 2*column*rmax;
+ else // Odd row
+ g[i].x = 2*rmax + 2*column*rmax;
+ g[i].y = rmax + 2*row*rmax;
+
+ }
+}
DIR diff --git a/main.c b/main.c
t@@ -0,0 +1,86 @@
+#include <iostream>
+#include <cmath>
+
+// Structure declarations and function prototypes
+#include "header.h"
+
+// 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()
+{
+ using std::cout;
+ using std::endl;
+
+ cout << "\n## simple_DEM ##\n"
+ << "Particles: " << ng << endl
+ << "maxStep: " << maxStep << endl;
+
+
+ double time = 0.0; // Time at simulation start
+
+ // Allocate memory
+ grain* g = new grain[ng]; // Grain structure
+
+
+ // Compute simulation domain dimensions
+ double wleft = 0.0; // Left wall
+ double wright = (ngw+1)*2*rmax; // Right wall
+ double wdown = 0.0; // Lower wall
+ double wup = (ng/ngw+1)*2*rmax; // Upper wall
+
+ // Variables for pressures on walls
+ double wp_up, wp_down, wp_left, wp_right;
+
+ // Initiailze grains inside the simulation domain
+ triangular_grid(g);
+
+
+
+ // Main time loop
+ for (int step = 0; step < maxStep; step++) {
+
+ time += dt; // Update current time
+
+ // Predict new kinematics
+ prediction(g);
+
+ // Calculate contact forces between grains
+ interact_grains(g);
+
+ // Calculate contact forces between grains and walls
+ interact_walls(g, wleft, wright, wup, wdown, &wp_up, &wp_down, &wp_left, &wp_right);
+
+ // Update acceleration from forces
+ update_acc(g);
+
+ // Correct velocities
+ correction(g);
+
+ // Write output files if the fileInterval is reached
+ if (step % fileInterval == 0) {
+ (void)vtk_export_grains(g, step);
+ (void)vtk_export_walls(step, wleft, wright, wdown, wup, wp_up, wp_down, wp_left, wp_right);
+ (void)vtk_export_forces(g, step);
+ }
+
+ } // End of main time loop
+
+
+ // Free dynamically allocated memory
+ delete[] g;
+
+
+ cout << "\nSimulation ended without errors.\n";
+ return 0; // Terminate successfully
+}
+
DIR diff --git a/output/blank.txt b/output/blank.txt
DIR diff --git a/vtk_export.c b/vtk_export.c
t@@ -0,0 +1,134 @@
+#include <iostream>
+#include <cstdio>
+#include <fstream>
+#include "header.h"
+#include "global_properties.h"
+
+int vtk_export_grains(grain* g, int numfile)
+{
+ FILE* fout;
+
+ char filename[25]; // File name
+ sprintf(filename, "output/grains%04d.vtk", numfile);
+
+ if ((fout = fopen(filename, "wt")) == NULL) {
+ std::cout << "vtk_export error, cannot open " << filename << "!\n";
+ return 1;
+ }
+
+ // Write header
+ fprintf(fout, "# vtk DataFile Version 3.0\n");
+
+ // Write title
+ fprintf(fout, "'Time: %f s'\n", numfile*dt);
+
+ // Write data type
+ 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, "%f %f 0.0\n", g[i].x, g[i].y);
+ fprintf(fout, "POINT_DATA %d\n", ng);
+
+ // Grain radii
+ fprintf(fout, "VECTORS Radius float\n");
+ for (int i = 0; i < ng; 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++)
+ 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++)
+ fprintf(fout, "%e\n", g[i].p);
+
+ // Angular velocity
+ fprintf(fout, "SCALARS Angvel float 1\n");
+ fprintf(fout, "LOOKUP_TABLE default\n");
+ for (int i = 0; i < ng; i++)
+ fprintf(fout, "%e\n", g[i].vth);
+
+ fclose(fout);
+ return 0;
+}
+
+// Save walls vtk file
+int vtk_export_walls(int numfile, double wlex, double wrix, double wloy, double wupy,
+ double wp_up, double wp_down, double wp_left, double wp_right)
+{
+ char fname[25]; // file name
+ sprintf(fname,"output/walls%04d.vtk",numfile);
+
+ using std::endl;
+
+ std::ofstream fow(fname, std::ios::out);
+ if (fow)
+ {
+ fow.precision(3); fow << std::scientific;
+ fow << "# vtk DataFile Version 3.0" << endl;
+ fow << "My walls" << endl;
+ fow << "ASCII" << endl;
+ fow << "DATASET POLYDATA" << endl;
+ fow << "POINTS 8 float" << endl;
+ // lower wall
+ fow << wlex << " " << wloy << " 0" << endl;
+ fow << wrix << " " << wloy << " 0" << endl;
+ // upper wall
+ fow << wlex << " " << wupy << " 0" << endl;
+ fow << wrix << " " << wupy << " 0" << endl;
+ // left wall
+ fow << wlex << " " << wloy << " 0" << endl;
+ fow << wlex << " " << wupy << " 0" << endl;
+ // right wall
+ fow << wrix << " " << wloy << " 0" << endl;
+ fow << wrix << " " << wupy << " 0" << endl;
+ fow << "LINES 4 12" << endl;
+ fow << "2 0 1" << endl;
+ fow << "2 2 3" << endl;
+ fow << "2 4 5" << endl;
+ fow << "2 6 7" << endl;
+ fow << "POINT_DATA 8" << endl;
+ fow << "SCALARS Rayon float" << endl;
+ fow << "LOOKUP_TABLE default" << endl;
+ fow << wp_up << endl; // No idea about the following sequence, i only know that there have to be 8 values
+ fow << wp_up << endl;
+ fow << wp_left << endl;
+ fow << wp_left << endl;
+ fow << wp_right << endl;
+ fow << wp_right << endl;
+ fow << wp_down << endl;
+ fow << wp_down << endl; // ??
+ }
+ return 0;
+}
+
+int vtk_export_forces(grain* g, int numfile)
+{
+ FILE* fout;
+
+ char filename[25]; // File name
+ sprintf(filename, "output/forces%04d.vtk", numfile);
+
+ if ((fout = fopen(filename, "wt")) == NULL) {
+ std::cout << "vtk_export error, cannot open " << filename << "!\n";
+ return 1;
+ }
+
+ // Write header
+ fprintf(fout, "# vtk DataFile Version 3.0\n");
+
+ // Write title
+ fprintf(fout, "'Time: %f s'\n", numfile*dt);
+
+ // Write data type
+ fprintf(fout, "ASCII\nDATASET POLYDATA\n");
+
+
+ fclose(fout);
+ return 0;
+}
DIR diff --git a/walls.c b/walls.c
t@@ -0,0 +1,109 @@
+#include <cmath>
+#include "header.h"
+#include "global_properties.h"
+
+// Compute force between i and upper wall
+double compute_force_upper_wall(int i, grain* g, double wup)
+{
+ double dn = wup - (g[i].y + g[i].R); // Overlap
+
+ if(dn<0.0) {
+ double vn,fn;
+ // velocity (wall velocity = 0)
+ vn = g[i].vy;
+ // force
+ fn = -kn * dn - nu * vn;
+ // Update sum of forces on grains
+ g[i].fy -= fn;
+ // Add fn to pressure on grains i
+ g[i].p += fn;
+ // Update stress to the wall
+ return fn;
+ }
+ return 0.0;
+}
+
+double compute_force_lower_wall(int i, grain* g, double wdown)
+{
+
+ double dn = g[i].y - g[i].R - wdown; // Overlap
+
+ if(dn<0.0) {
+ double vn,fn;
+ // velocity (wall velocity = 0)
+ vn = g[i].vy;
+ // force
+ fn = - kn * dn - nu * vn;
+ // Update sum of forces on grains
+ g[i].fy += fn;
+ // Add fn to pressure on grains i
+ g[i].p += fn;
+ // Update stress to the wall
+ return fn;
+ }
+ return 0.0;
+}
+
+double compute_force_left_wall(int i, grain* g, double wleft)
+{
+ double dn = wleft + g[i].x - g[i].R; // Overlap
+
+ if(dn<0.0) {
+ double vn,fn;
+ // velocity (wall velocity = 0)
+ vn = g[i].vx;
+ // force
+ fn = -kn * dn - nu * vn;
+ // Update sum of forces on grains
+ g[i].fx += fn;
+ // Add fn to pressure on grains i
+ g[i].p += fn;
+ // Update stress to the wall
+ return fn;
+ }
+ return 0.0;
+}
+
+
+double compute_force_right_wall(int i, grain* g, double wright)
+{
+ double dn = wright - (g[i].x + g[i].R); // Overlap
+
+ if(dn<0.0) {
+ double vn,fn;
+ // velocity (wall velocity = 0)
+ vn = -g[i].vx;
+ // force
+ fn = -kn * dn - nu * vn;
+ // Update sum of forces on grains
+ g[i].fx -= fn;
+ // Add fn to pressure on grains i
+ g[i].p += fn;
+ // Update stress to the wall
+ return fn;
+ }
+ return 0.0;
+}
+
+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)
+{
+ double u = 0.0;
+ double d = 0.0;
+ double r = 0.0;
+ double l = 0.0;
+
+ int i;
+
+ #pragma omp parallel for shared(g, u, d, r, l) private (i)
+ for (i = 0; i < ng; 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);
+ l += compute_force_left_wall(i, g, wleft);
+ }
+ *wp_up = u;
+ *wp_down = d;
+ *wp_left = l;
+ *wp_right = r;
+}