tRemoved legacy c++ syntax - 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 c718a3c6307caad2ff4b083b9abfe77a1360bdf9
(DIR) parent 303835daee186810677bd22cd133d18386f643f7
(HTM) Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Sun, 14 Oct 2012 09:02:13 +0200
Removed legacy c++ syntax
Diffstat:
M global_properties.h | 40 ++++++++++++++++----------------
M grains.c | 34 ++++++++++++++++----------------
M header.h | 31 ++++++++++++++++---------------
M initialize.c | 8 ++++----
M main.c | 63 +++++++++++++++----------------
M vtk_export.c | 84 +++++++------------------------
M walls.c | 52 ++++++++++++++++----------------
7 files changed, 130 insertions(+), 182 deletions(-)
---
(DIR) diff --git a/global_properties.h b/global_properties.h
t@@ -1,29 +1,29 @@
-#ifndef GLOBAL_CONSTANTS_H_
-#define GLOBAL_CONSTANTS_H_
+#ifndef GLOBAL_PROPERTIES_H_
+#define GLOBAL_PROPERTIES_H_
-// Properties of sample
-const int ng = 5000; // Number of grains
+/* 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
+/* 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
+/* 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
+/* 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
+/* Physical constants */
+const double grav = 9.80; /* Gravity */
-// Number of grains along the width in the initial configuration
+/* 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@@ -4,7 +4,7 @@ void prediction(grain* g)
#pragma omp parallel for shared(g) private (i)
for (i = 0; i < ng; i++) {
- // Predict positions and velocities
+ /* 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;
t@@ -13,7 +13,7 @@ void prediction(grain* g)
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
+ /* Zero forces */
g[i].fx = 0.0;
g[i].fy = 0.0;
g[i].fth = 0.0;
t@@ -23,46 +23,46 @@ void prediction(grain* g)
void interparticle_force(grain* g, int a, int b)
{
- if (a > b) { // Use Newtons 3rd law to find both forces at once
+ if (a > b) { /* Use Newtons 3rd law to find both forces at once */
- // Particle center coordinate component differences
+ /* 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
+ /* Particle center distance */
double dist = sqrt(x_ab*x_ab + y_ab*y_ab);
- // Size of overlap
+ /* 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
+ 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
+ /* 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
+ /* Compute force in local axes */
fn = -kn * dn - nu * vn;
- // Rotation
+ /* Rotation */
if (fn < 0)
fn = 0.0;
ft = fabs(kt * vt);
- if (ft > mu*fn) // Coefficient of friction
+ 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
+ /* 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;
t@@ -83,10 +83,10 @@ void interact_grains(grain* g)
{
int a, b;
#pragma omp parallel for shared(g) private (a,b)
- // Loop through particle a
+ /* Loop through particle a */
for (a = 0; a < ng; a++) {
- // Loop through particle b
+ /* Loop through particle b */
for (b = 0; b < ng; b++) {
interparticle_force(g, a, b);
}
(DIR) diff --git a/header.h b/header.h
t@@ -1,43 +1,44 @@
#ifndef HEADER_H_
#define HEADER_H_
-//// Structure declarations
+/* 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
+ 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
+/* Prototype functions */
-// initialize.cpp
+/* initialize.cpp */
void triangular_grid(grain* g);
-// grains.cpp
+/* 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
+/* 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
+/* 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_walls(int numfile, double wlex, double wrix, double wloy, double wupy);*/
int vtk_export_forces(grain* g, int numfile);
(DIR) diff --git a/initialize.c b/initialize.c
t@@ -5,7 +5,7 @@
void triangular_grid(grain* g)
{
- // Initialize grain properties
+ /* 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;
t@@ -13,14 +13,14 @@ void triangular_grid(grain* g)
g[i].p = 0.0;
}
- // Initialize grain positions in a triangular grid
+ /* 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
+ if (row%2 == 0) /* Even row */
g[i].x = rmax + 2*column*rmax;
- else // Odd row
+ 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@@ -1,13 +1,13 @@
-#include <iostream>
-#include <cmath>
+#include <stdio.h>
+#include <math.h>
-// Structure declarations and function prototypes
+/* Structure declarations and function prototypes */
#include "header.h"
-// Global and constant simulation properties
+/* Global and constant simulation properties */
#include "global_properties.h"
-// Functions for exporting data to VTK formats
+/* Functions for exporting data to VTK formats */
#include "vtk_export.c"
#include "initialize.c"
t@@ -18,69 +18,66 @@
int main()
{
- using std::cout;
- using std::endl;
- cout << "\n## simple_DEM ##\n"
- << "Particles: " << ng << endl
- << "maxStep: " << maxStep << endl;
+ printf("\n## simple_DEM ##\n");
+ printf("Particles: %d\n");
+ printf("maxStep: %d\n");
- double time = 0.0; // Time at simulation start
+ double time = 0.0; /* Time at simulation start */
- // Allocate memory
- grain* g = new grain[ng]; // Grain structure
+ /* 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
+ /* 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
+ /* Variables for pressures on walls */
double wp_up, wp_down, wp_left, wp_right;
- // Initiailze grains inside the simulation domain
+ /* Initiailze grains inside the simulation domain */
triangular_grid(g);
- // Main time loop
+ /* Main time loop */
for (int step = 0; step < maxStep; step++) {
- time += dt; // Update current time
+ time += dt; /* Update current time */
- // Predict new kinematics
+ /* Predict new kinematics */
prediction(g);
- // Calculate contact forces between grains
+ /* Calculate contact forces between grains */
interact_grains(g);
- // Calculate contact forces between grains and walls
+ /* 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 acceleration from forces */
update_acc(g);
- // Correct velocities
+ /* Correct velocities */
correction(g);
- // Write output files if the fileInterval is reached
+ /* 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
+ } /* End of main time loop */
- // Free dynamically allocated memory
+ /* Free dynamically allocated memory */
delete[] g;
- cout << "\nSimulation ended without errors.\n";
- return 0; // Terminate successfully
+ printf("\nSimulation ended without errors.\n");
+ return 0; /* Terminate successfully */
}
(DIR) diff --git a/vtk_export.c b/vtk_export.c
t@@ -1,6 +1,4 @@
-#include <iostream>
-#include <cstdio>
-#include <fstream>
+#include <stdio.h>
#include "header.h"
#include "global_properties.h"
t@@ -8,46 +6,48 @@ int vtk_export_grains(grain* g, int numfile)
{
FILE* fout;
- char filename[25]; // File name
+ 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";
+ printf("vtk_export error, cannot open ")
+ printf(filename);
+ printf("!\n");
return 1;
}
- // Write header
+ /* Write header */
fprintf(fout, "# vtk DataFile Version 3.0\n");
- // Write title
+ /* Write title */
fprintf(fout, "'Time: %f s'\n", numfile*dt);
- // Write data type
+ /* Write data type */
fprintf(fout, "ASCII\nDATASET UNSTRUCTURED_GRID\n");
- // Grain coordinates
+ /* 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
+ /* 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
+ /* 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
+ /* 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
+ /* Angular velocity */
fprintf(fout, "SCALARS Angvel float 1\n");
fprintf(fout, "LOOKUP_TABLE default\n");
for (int i = 0; i < ng; i++)
t@@ -57,61 +57,11 @@ int vtk_export_grains(grain* g, int numfile)
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
+ char filename[25]; /* File name */
sprintf(filename, "output/forces%04d.vtk", numfile);
if ((fout = fopen(filename, "wt")) == NULL) {
t@@ -119,13 +69,13 @@ int vtk_export_forces(grain* g, int numfile)
return 1;
}
- // Write header
+ /* Write header */
fprintf(fout, "# vtk DataFile Version 3.0\n");
- // Write title
+ /* Write title */
fprintf(fout, "'Time: %f s'\n", numfile*dt);
- // Write data type
+ /* Write data type */
fprintf(fout, "ASCII\nDATASET POLYDATA\n");
(DIR) diff --git a/walls.c b/walls.c
t@@ -1,23 +1,23 @@
-#include <cmath>
+#include <math.h>
#include "header.h"
#include "global_properties.h"
-// Compute force between i and upper wall
+/* 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
+ double dn = wup - (g[i].y + g[i].R); /* Overlap */
if(dn<0.0) {
double vn,fn;
- // velocity (wall velocity = 0)
+ /* velocity (wall velocity = 0) */
vn = g[i].vy;
- // force
+ /* force */
fn = -kn * dn - nu * vn;
- // Update sum of forces on grains
+ /* Update sum of forces on grains */
g[i].fy -= fn;
- // Add fn to pressure on grains i
+ /* Add fn to pressure on grains i */
g[i].p += fn;
- // Update stress to the wall
+ /* Update stress to the wall */
return fn;
}
return 0.0;
t@@ -26,19 +26,19 @@ double compute_force_upper_wall(int i, grain* g, double wup)
double compute_force_lower_wall(int i, grain* g, double wdown)
{
- double dn = g[i].y - g[i].R - wdown; // Overlap
+ double dn = g[i].y - g[i].R - wdown; /* Overlap */
if(dn<0.0) {
double vn,fn;
- // velocity (wall velocity = 0)
+ /* velocity (wall velocity = 0) */
vn = g[i].vy;
- // force
+ /* force */
fn = - kn * dn - nu * vn;
- // Update sum of forces on grains
+ /* Update sum of forces on grains */
g[i].fy += fn;
- // Add fn to pressure on grains i
+ /* Add fn to pressure on grains i */
g[i].p += fn;
- // Update stress to the wall
+ /* Update stress to the wall */
return fn;
}
return 0.0;
t@@ -46,19 +46,19 @@ double compute_force_lower_wall(int i, grain* g, double wdown)
double compute_force_left_wall(int i, grain* g, double wleft)
{
- double dn = wleft + g[i].x - g[i].R; // Overlap
+ double dn = wleft + g[i].x - g[i].R; /* Overlap */
if(dn<0.0) {
double vn,fn;
- // velocity (wall velocity = 0)
+ /* velocity (wall velocity = 0) */
vn = g[i].vx;
- // force
+ /* force */
fn = -kn * dn - nu * vn;
- // Update sum of forces on grains
+ /* Update sum of forces on grains */
g[i].fx += fn;
- // Add fn to pressure on grains i
+ /* Add fn to pressure on grains i */
g[i].p += fn;
- // Update stress to the wall
+ /* Update stress to the wall */
return fn;
}
return 0.0;
t@@ -67,19 +67,19 @@ double compute_force_left_wall(int i, grain* g, double wleft)
double compute_force_right_wall(int i, grain* g, double wright)
{
- double dn = wright - (g[i].x + g[i].R); // Overlap
+ double dn = wright - (g[i].x + g[i].R); /* Overlap */
if(dn<0.0) {
double vn,fn;
- // velocity (wall velocity = 0)
+ /* velocity (wall velocity = 0) */
vn = -g[i].vx;
- // force
+ /* force */
fn = -kn * dn - nu * vn;
- // Update sum of forces on grains
+ /* Update sum of forces on grains */
g[i].fx -= fn;
- // Add fn to pressure on grains i
+ /* Add fn to pressure on grains i */
g[i].p += fn;
- // Update stress to the wall
+ /* Update stress to the wall */
return fn;
}
return 0.0;