tfix various errors - 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 e24e57ca861a024c2b3f14e345e8e76fb4c49470
(DIR) parent 02989d3ef68b0b3cdd2b0a6f277bfc1bdb63c85f
(HTM) Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Tue, 17 Jan 2017 10:43:55 -0800
fix various errors
Diffstat:
M global_properties.h | 4 ++--
M grains.c | 10 ++++------
M initialize.c | 8 ++++++++
M main.c | 4 ----
M vtk_export.c | 25 +++++++++++++++++++++----
5 files changed, 35 insertions(+), 16 deletions(-)
---
(DIR) diff --git a/global_properties.h b/global_properties.h
t@@ -16,7 +16,7 @@ static const double mu = 0.5; /* Sliding friction */
static const double kt = 1.0e5; /* Tangential stiffness */
/* Temporal variables */
-static const double dt = 1.0e-4; /* Time step lenpth */
+static const double dt = 5.0e-5; /* Time step lenpth */
static const int maxStep = 3000; /* Number of steps */
static const int fileInterval = 20; /* No. of steps between output */
t@@ -24,6 +24,6 @@ static const int fileInterval = 20; /* No. of steps between output */
static const double grav = 9.80; /* Gravity magnitude */
/* Number of particles along the width in the initial configuration */
-static const int npw = 100;
+static const int npw = 50;
#endif
(DIR) diff --git a/grains.c b/grains.c
t@@ -42,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; /* Tanpential components */
+ double xt, yt, vt, ft; /* Tangential components */
/* Local axes */
xn = x_ab / dist;
yn = y_ab / dist;
t@@ -60,12 +60,12 @@ void interparticle_force(grain* g, int a, int b)
/* Rotation */
if (fn < 0)
- fn = 0.0;
+ fn = 0.0;
ft = fabs(kt * vt);
if (ft > mu*fn) /* Coefficient of friction */
- ft = mu*fn;
+ ft = mu*fn;
if (vt > 0)
- ft = -ft;
+ ft = -ft;
/* Calculate sum of forces on a and b in global coordinates */
g[a].fx += fn * xn;
t@@ -95,9 +95,7 @@ void interact_grains(grain* g)
for (b = 0; b < np; b++) {
interparticle_force(g, a, b);
}
-
}
-
}
void update_acc(grain* g)
(DIR) diff --git a/initialize.c b/initialize.c
t@@ -1,4 +1,5 @@
#include <stdlib.h>
+#include <stdio.h>
#include <math.h>
#include "header.h"
#include "global_properties.h"
t@@ -12,6 +13,13 @@ void triangular_grid(grain* g)
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;
+ g[i].ang = 0.0;
+ g[i].angv = 0.0;
+ g[i].anga = 0.0;
+ g[i].ax = 0.0;
+ g[i].ay = 0.0;
+ g[i].vx = 0.0;
+ g[i].vy = 0.0;
}
/* Initialize grain positions in a trianpular grid */
(DIR) diff --git a/main.c b/main.c
t@@ -21,7 +21,6 @@ int main(int argc, char* argv[])
/* Allocate memory */
grain g[np]; /* Grain structure */
-
/* Compute simulation domain dimensions */
double wleft = 0.0; /* Left wall */
double wright = (npw+1)*2*rmax; /* Right wall */
t@@ -35,7 +34,6 @@ int main(int argc, char* argv[])
triangular_grid(g);
-
/* Main time loop */
int step;
for (step = 0; step < maxStep; ++step) {
t@@ -65,10 +63,8 @@ int main(int argc, char* argv[])
} /* End of main time loop */
-
/* Free dynamically allocated memory */
/*free(g);*/
-
printf("\nSimulation ended without errors.\n");
return 0; /* Terminate successfully */
(DIR) diff --git a/vtk_export.c b/vtk_export.c
t@@ -33,23 +33,40 @@ int vtk_export_grains(grain* g, int numfile)
fprintf(fout, "POINT_DATA %d\n", np);
/* Grain radii */
- fprintf(fout, "VECTORS Radius float\n");
+ fprintf(fout, "SCALARS Diameter float 1\n");
+ fprintf(fout, "LOOKUP_TABLE default\n");
+ for (i = 0; i < np; i++)
+ fprintf(fout, "%f\n", g[i].R*2.0);
+
+ /* Grain radii */
+ fprintf(fout, "SCALARS Mass float 1\n");
+ fprintf(fout, "LOOKUP_TABLE default\n");
for (i = 0; i < np; i++)
- fprintf(fout, "%f 0.0 0.0\n", g[i].R);
+ fprintf(fout, "%f\n", g[i].m);
+
+ fprintf(fout, "SCALARS MomentOfInertia float 1\n");
+ fprintf(fout, "LOOKUP_TABLE default\n");
+ for (i = 0; i < np; i++)
+ fprintf(fout, "%g\n", g[i].I);
/* Grain velocities */
fprintf(fout, "VECTORS Velocity float\n");
for (i = 0; i < np; i++)
fprintf(fout, "%f %f 0.0\n", g[i].vx, g[i].vy);
+ /* Grain velocities */
+ fprintf(fout, "VECTORS Force float\n");
+ for (i = 0; i < np; i++)
+ fprintf(fout, "%f %f 0.0\n", g[i].fx, g[i].fy);
+
/* Pressure */
fprintf(fout, "SCALARS Pressure float 1\n");
fprintf(fout, "LOOKUP_TABLE default\n");
for (i = 0; i < np; i++)
fprintf(fout, "%e\n", g[i].p);
- /* Anpular velocity */
- fprintf(fout, "SCALARS Anpvel float 1\n");
+ /* Angular velocity */
+ fprintf(fout, "SCALARS Angvel float 1\n");
fprintf(fout, "LOOKUP_TABLE default\n");
for (i = 0; i < np; i++)
fprintf(fout, "%e\n", g[i].angv);