tgather headers in granular.h, rename np to ng, split operations on grain arrays - granular - granular dynamics simulation
 (HTM) git clone git://src.adamsgaard.dk/granular
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) README
 (DIR) LICENSE
       ---
 (DIR) commit df197a3f20247716945dd8fa4ce272bbf797d422
 (DIR) parent fef9eb22a304c90e4aed9a94978a8efe6c13fb26
 (HTM) Author: Anders Damsgaard <anders@adamsgaard.dk>
       Date:   Sun, 28 Mar 2021 09:22:21 +0200
       
       gather headers in granular.h, rename np to ng, split operations on grain arrays
       
       Diffstat:
         M arrays.h                            |       2 ++
         M grain.c                             |       1 -
         A grains.c                            |     128 +++++++++++++++++++++++++++++++
         A grains.h                            |      12 ++++++++++++
         M granular.c                          |       4 ++--
         M granular.h                          |       7 +++++++
         M granular2vtu.c                      |       2 +-
         M granularpacking.c                   |       9 +++------
         M packing.c                           |      38 ++++++++++++++++----------------
         M simulation.c                        |     132 +++++++------------------------
         M simulation.h                        |       9 +++++----
       
       11 files changed, 208 insertions(+), 136 deletions(-)
       ---
 (DIR) diff --git a/arrays.h b/arrays.h
       t@@ -1,6 +1,8 @@
        #ifndef ARRAYS_
        #define ARRAYS_
        
       +#include <stdio.h>
       +
        unsigned int idx3(
                const unsigned int i, const unsigned int j, const unsigned int k,
                const unsigned int nx, const unsigned int ny);
 (DIR) diff --git a/grain.c b/grain.c
       t@@ -279,4 +279,3 @@ grain_kinetic_energy(const struct grain *g)
                return grain_translational_kinetic_energy(g) +
                        grain_rotational_kinetic_energy(g);
        }
       -
 (DIR) diff --git a/grains.c b/grains.c
       t@@ -0,0 +1,128 @@
       +#include <stdio.h>
       +#include <math.h>
       +#include "grain.h"
       +
       +#define VTK_FLOAT_FMT "%.17g "
       +
       +#define VTK_XML_SCALAR(M, N, T, F) \
       +        fprintf(stream,\
       +                "\t\t\t\t<DataArray type=\"" T "\" Name=\"" N "\" "\
       +                "NumberOfComponents=\"1\" format=\"ascii\">\n");\
       +        for (i = 0; i < ng; i++)\
       +                fprintf(stream, F, grains[i].M);\
       +        fprintf(stream, "\n\t\t\t\t</DataArray>\n");
       +
       +#define VTK_XML_VECTOR(M, N, T, F) \
       +        fprintf(stream,\
       +                "\t\t\t\t<DataArray type=\"" T "\" Name=\"" N "\" "\
       +                "NumberOfComponents=\"3\" format=\"ascii\">\n");\
       +        for (i = 0; i < ng; i++)\
       +                for (d = 0; d < 3; d++)\
       +                        fprintf(stream, F, grains[i].M[d]);\
       +        fprintf(stream, "\n\t\t\t\t</DataArray>\n");
       +
       +void
       +grains_print(FILE *stream, const struct grain *grains, size_t ng)
       +{
       +        size_t i;
       +
       +        for (i = 0; i < ng; i++)
       +                grain_print(stream, &grains[i]);
       +}
       +
       +void
       +grains_print_vtk(FILE *stream, const struct grain *grains, size_t ng)
       +{
       +        size_t i, d;
       +
       +        fprintf(stream,
       +                "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n"
       +                "<VTKFile type=\"UnstructuredGrid\" version=\"1.0\" "
       +                "byte_order=\"LittleEndian\">\n"
       +                "\t<UnstructuredGrid>\n"
       +                "\t\t<Piece NumberOfPoints=\"%zu\" NumberOfCells=\"0\">\n", ng);
       +        fprintf(stream, "\t\t\t<Points>\n");
       +        VTK_XML_VECTOR(pos, "Position [m]", "Float64", VTK_FLOAT_FMT);
       +        fprintf(stream, "\t\t\t</Points>\n");
       +
       +        fprintf(stream,
       +                "\t\t\t<Cells>\n"
       +                "\t\t\t\t<DataArray type=\"Int32\" Name=\"connectivity\" "
       +                "NumberOfComponents=\"1\" format=\"ascii\"/>\n"
       +                "\t\t\t\t<DataArray type=\"Int32\" Name=\"offsets\" "
       +                "NumberOfComponents=\"1\" format=\"ascii\"/>\n"
       +                "\t\t\t\t<DataArray type=\"UInt8\" Name=\"types\" "
       +                "NumberOfComponents=\"1\" format=\"ascii\"/>\n"
       +                "\t\t\t</Cells>\n");
       +
       +        fprintf(stream,
       +                "\t\t\t<PointData Scalars=\"Diameter [m]\" "
       +                "Vectors=\"Angular position [-]\">\n");
       +
       +        fprintf(stream,
       +                "\t\t\t\t<DataArray type=\"Float64\" Name=\"Diameter [m]\" "
       +                "NumberOfComponents=\"1\" format=\"ascii\">\n");
       +        for (i = 0; i < ng; i++)
       +                fprintf(stream, VTK_FLOAT_FMT, grains[i].radius * 2.0);
       +        fprintf(stream,
       +                "\n"
       +                "\t\t\t\t</DataArray>\n");
       +
       +        VTK_XML_VECTOR(vel, "Velocity [m/s]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_VECTOR(acc, "Acceleration [m/s^2]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_VECTOR(force, "Force [N]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_VECTOR(angpos, "Angular position [rad]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_VECTOR(angvel, "Angular velocity [rad/s]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_VECTOR(angacc, "Angular acceleration [rad/s^2]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_VECTOR(torque, "Torque [N/m]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_VECTOR(disp, "Displacement [m]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_VECTOR(forceext, "External body force [N]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_SCALAR(density, "Density [kg/m^3]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_SCALAR(fixed, "Fixed [-]", "Int64", "%d ");
       +        VTK_XML_SCALAR(rotating, "Rotating [-]", "Int64", "%d ");
       +        VTK_XML_SCALAR(enabled, "Enabled [-]", "Int64", "%d ");
       +        VTK_XML_SCALAR(youngs_modulus, "Young's modulus [Pa]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_SCALAR(poissons_ratio, "Poisson's ratio [-]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_SCALAR(friction_coeff, "Friction coefficient [-]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_SCALAR(tensile_strength, "Tensile strength [Pa]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_SCALAR(shear_strength, "Shear strength [Pa]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_SCALAR(fracture_toughness, "Fracture toughness [Pa]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_VECTOR(gridpos, "Grid position [-]", "UInt64", "%zu ");
       +        VTK_XML_SCALAR(ncontacts, "Number of contacts [-]", "UInt64", "%zu ");
       +        VTK_XML_VECTOR(contact_stress, "Contact stress [Pa]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_SCALAR(thermal_energy, "Thermal energy [J]", "Float64", VTK_FLOAT_FMT);
       +        VTK_XML_SCALAR(color, "Color [-]", "Int64", "%d ");
       +
       +        fprintf(stream, "\t\t\t</PointData>\n");
       +
       +        fprintf(stream,
       +                "\t\t</Piece>\n"
       +                "\t</UnstructuredGrid>\n"
       +                "</VTKFile>\n");
       +}
       +
       +double
       +grains_minimum_grain_edge(const struct grain *grains, size_t ng, int d)
       +{
       +        size_t i;
       +        double res = INFINITY;
       +
       +        for (i = 0; i < ng; i++)
       +                if (res > grains[i].pos[d] - grains[i].radius)
       +                        res = grains[i].pos[d] - grains[i].radius;
       +
       +        return res;
       +}
       +
       +double
       +grains_maximum_grain_edge(const struct grain *grains, size_t ng, int d)
       +{
       +        size_t i;
       +        double res = -INFINITY;
       +
       +        for (i = 0; i < ng; i++)
       +                if (res < grains[i].pos[d] + grains[i].radius)
       +                        res = grains[i].pos[d] + grains[i].radius;
       +
       +        return res;
       +}
 (DIR) diff --git a/grains.h b/grains.h
       t@@ -0,0 +1,12 @@
       +#ifndef GRANULAR_GRAINS_
       +#define GRANULAR_GRAINS_
       +
       +#include "grain.h"
       +
       +void grains_print(FILE *stream, const struct grain *grains, size_t ng);
       +void grains_print_vtk(FILE *stream, const struct grain *grains, size_t ng);
       +
       +double sim_minimum_grain_edge(const struct grain *grains, size_t ng, int d);
       +double sim_maximum_grain_edge(const struct grain *grains, size_t ng, int d);
       +
       +#endif
 (DIR) diff --git a/granular.c b/granular.c
       t@@ -1,7 +1,7 @@
        #include <stdlib.h>
        #include <err.h>
        #include <unistd.h>
       -#include "simulation.h"
       +#include "granular.h"
        #include "arg.h"
        
        char *argv0;
       t@@ -42,7 +42,7 @@ main(int argc, char *argv[])
                        usage();
                        break;
                case 'I':
       -                sim.file_dt = atof(EARGF(usage()));
       +                sim.dt_file = atof(EARGF(usage()));
                        break;
                case 'j':
                        sim.dt = atof(EARGF(usage()));
 (DIR) diff --git a/granular.h b/granular.h
       t@@ -1,6 +1,13 @@
        #ifndef GRANULAR_
        #define GRANULAR_
        
       +#include "arrays.h"
       +#include "grain.h"
       +#include "grains.h"
       +#include "packing.h"
       +#include "simulation.h"
       +#include "util.h"
       +
        #define PI 3.14159265358979323846
        #define DEG2RAD(x) (x*PI/180.0)
        
 (DIR) diff --git a/granular2vtu.c b/granular2vtu.c
       t@@ -1,7 +1,7 @@
        #include <stdlib.h>
        #include <unistd.h>
        #include <err.h>
       -#include "simulation.h"
       +#include "granular.h"
        
        int
        main(void)
 (DIR) diff --git a/granularpacking.c b/granularpacking.c
       t@@ -1,10 +1,7 @@
        #include <stdlib.h>
        #include <stdlib.h>
        #include <err.h>
       -#include "packing.h"
       -#include "simulation.h"
       -#include "arrays.h"
       -#include "util.h"
       +#include "granular.h"
        #include "arg.h"
        
        char *argv0;
       t@@ -83,12 +80,12 @@ main(int argc, char *argv[])
                
                switch (packing) {
                case 0:
       -                sim.np += rectangular_packing(&sim.grains, n,
       +                sim.ng += rectangular_packing(&sim.grains, n,
                                                      d_min / 2.0, d_max / 2.0,
                                                      sizefunc, padding, origo);
                        break;
                case 1:
       -                sim.np += triangular_packing(&sim.grains, n,
       +                sim.ng += triangular_packing(&sim.grains, n,
                                                     d_min / 2.0, d_max / 2.0,
                                                     sizefunc, padding, origo);
                        break;
 (DIR) diff --git a/packing.c b/packing.c
       t@@ -15,7 +15,7 @@ rectangular_packing(struct grain **grains,
                            double padding_factor,
                            double origo[3])
        {
       -        size_t i, j, k, l, np, N = 0;
       +        size_t i, j, k, l, ng, N = 0;
                double dx_padding = radius_max * 2.0 * padding_factor;
                double dx = radius_max * 2.0 + dx_padding;
        
       t@@ -30,18 +30,18 @@ rectangular_packing(struct grain **grains,
                        for (j = 0; j < n[1]; j++)
                                for (i = 0; i < n[0]; i++) {
                                        N++;
       -                                np = idx3(i, j, k, n[0], n[1]);
       -                                grain_defaults(&(*grains)[np]);
       -                                (*grains)[np].radius = gsd(radius_min, radius_max);
       -                                (*grains)[np].pos[0] = i * dx + 0.5 * dx + origo[0];
       -                                (*grains)[np].pos[1] = j * dx + 0.5 * dx + origo[1];
       -                                (*grains)[np].pos[2] = k * dx + 0.5 * dx + origo[2];
       +                                ng = idx3(i, j, k, n[0], n[1]);
       +                                grain_defaults(&(*grains)[ng]);
       +                                (*grains)[ng].radius = gsd(radius_min, radius_max);
       +                                (*grains)[ng].pos[0] = i * dx + 0.5 * dx + origo[0];
       +                                (*grains)[ng].pos[1] = j * dx + 0.5 * dx + origo[1];
       +                                (*grains)[ng].pos[2] = k * dx + 0.5 * dx + origo[2];
                                        for (l = 0; l < 3; l++)
       -                                        (*grains)[np].pos[l] +=
       +                                        (*grains)[ng].pos[l] +=
                                                        random_value_uniform(-0.5 * dx_padding,
                                                                             0.5 * dx_padding);
        #ifdef VERBOSE
       -                                printf("added grain %zu\n", np);
       +                                printf("added grain %zu\n", ng);
        #endif
                                }
        
       t@@ -56,7 +56,7 @@ triangular_packing(struct grain **grains,
                           double padding_factor,
                           double origo[3])
        {
       -        size_t i, j, k, l, np, N = 0;
       +        size_t i, j, k, l, ng, N = 0;
                double dx_padding = radius_max * 2.0 * padding_factor;
                double dx = radius_max * 2.0 + dx_padding;
                double dy = dx * sin(PI/3.0);
       t@@ -72,20 +72,20 @@ triangular_packing(struct grain **grains,
                        for (j = 0; j < n[1]; j++)
                                for (i = 0; i < n[0]; i++) {
                                        N++;
       -                                np = idx3(i, j, k, n[0], n[1]);
       -                                grain_defaults(&(*grains)[np]);
       -                                (*grains)[np].radius = gsd(radius_min, radius_max);
       -                                (*grains)[np].pos[0] = i * dx + 0.5 * dx + origo[0];
       +                                ng = idx3(i, j, k, n[0], n[1]);
       +                                grain_defaults(&(*grains)[ng]);
       +                                (*grains)[ng].radius = gsd(radius_min, radius_max);
       +                                (*grains)[ng].pos[0] = i * dx + 0.5 * dx + origo[0];
                                        if (j % 2 == 0)
       -                                        (*grains)[np].pos[0] += 0.5 * dx;
       -                                (*grains)[np].pos[1] = j * dy + 0.5 * dx + origo[1];
       -                                (*grains)[np].pos[2] = k * dx + 0.5 * dx + origo[2];
       +                                        (*grains)[ng].pos[0] += 0.5 * dx;
       +                                (*grains)[ng].pos[1] = j * dy + 0.5 * dx + origo[1];
       +                                (*grains)[ng].pos[2] = k * dx + 0.5 * dx + origo[2];
                                        for (l = 0; l < 2; l++)
       -                                        (*grains)[np].pos[l] +=
       +                                        (*grains)[ng].pos[l] +=
                                                        random_value_uniform(-0.5 * dx_padding,
                                                                                                 0.5 * dx_padding);
        #ifdef VERBOSE
       -                                printf("added grain %zu\n", np);
       +                                printf("added grain %zu\n", ng);
        #endif
                                }
        
 (DIR) diff --git a/simulation.c b/simulation.c
       t@@ -7,29 +7,11 @@
        #include <err.h>
        #include <string.h>
        #include "grain.h"
       +#include "grains.h"
        #include "simulation.h"
        #include "arrays.h"
        #include "util.h"
        
       -#define VTK_FLOAT_FMT "%.17g "
       -
       -#define VTK_XML_SCALAR(M, N, T, F) \
       -        fprintf(stream,\
       -                "\t\t\t\t<DataArray type=\"" T "\" Name=\"" N "\" "\
       -                "NumberOfComponents=\"1\" format=\"ascii\">\n");\
       -        for (i = 0; i < n; i++)\
       -                fprintf(stream, F, grains[i].M);\
       -        fprintf(stream, "\n\t\t\t\t</DataArray>\n");
       -
       -#define VTK_XML_VECTOR(M, N, T, F) \
       -        fprintf(stream,\
       -                "\t\t\t\t<DataArray type=\"" T "\" Name=\"" N "\" "\
       -                "NumberOfComponents=\"3\" format=\"ascii\">\n");\
       -        for (i = 0; i < n; i++)\
       -                for (d = 0; d < 3; d++)\
       -                        fprintf(stream, F, grains[i].M[d]);\
       -        fprintf(stream, "\n\t\t\t\t</DataArray>\n");
       -
        struct simulation
        sim_new(void)
        {
       t@@ -44,7 +26,7 @@ sim_defaults(struct simulation *sim)
                size_t i;
        
                snprintf(sim->name, sizeof(sim->name), DEFAULT_SIMULATION_NAME);
       -        sim->np = 0;
       +        sim->ng = 0;
                for (i = 0; i < 3; i++) {
                        sim->constacc[i] = 0.0;
                        sim->nd[i] = 1;
       t@@ -54,9 +36,11 @@ sim_defaults(struct simulation *sim)
                sim->t = 0.0;
                sim->t_end = 0.0;
                sim->dt = 0.0;
       -        sim->file_dt = 1.0;
       +        sim->dt_file = 1.0;
                sim->n_file = 0;
       -        sim->np = 0;
       +        sim->iter = 0;
       +        sim->t_file = 0;
       +        sim->ng = 0;
                sim->grains = NULL;
        }
        
       t@@ -65,15 +49,15 @@ sim_free(struct simulation *sim)
        {
                free(sim->grains);
                sim->grains = NULL;
       -        sim->np = 0;
       +        sim->ng = 0;
        }
        
        void
        sim_add_grain(struct simulation *sim, struct grain *g)
        {
       -        if (!(sim->grains = xreallocarray(sim->grains, sim->np + 1, sizeof(*g))))
       +        if (!(sim->grains = xreallocarray(sim->grains, sim->ng + 1, sizeof(*g))))
                        err(1, "%s: sim.grains reallocarray", __func__);
       -        sim->grains[sim->np++] = *g;
       +        sim->grains[sim->ng++] = *g;
                free(g);
        }
        
       t@@ -91,98 +75,40 @@ sim_read_grains(struct simulation *sim, FILE *stream)
        }
        
        void
       -print_grains(FILE *stream, const struct grain *grains, size_t n)
       +sim_print_grains(FILE *stream, const struct simulation *sim)
        {
       -        size_t i;
       -
       -        for (i = 0; i < n; i++)
       -                grain_print(stream, &grains[i]);
       +        grains_print(stream, sim->grains, sim->ng);
        }
        
        void
       -sim_print_grains(FILE *stream, const struct simulation *sim)
       +sim_print_grains_vtk(FILE *stream, const struct simulation *sim)
        {
       -        print_grains(stream, sim->grains, sim->np);
       +        grains_print_vtk(stream, sim->grains, sim->ng);
        }
        
        void
       -print_grains_vtk(FILE *stream, const struct grain *grains, size_t n)
       +sim_step_time(struct simulation *sim)
        {
       -        size_t i, d;
       -
       -        fprintf(stream,
       -                "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n"
       -                "<VTKFile type=\"UnstructuredGrid\" version=\"1.0\" "
       -                "byte_order=\"LittleEndian\">\n"
       -                "\t<UnstructuredGrid>\n"
       -                "\t\t<Piece NumberOfPoints=\"%zu\" NumberOfCells=\"0\">\n", n);
       -        fprintf(stream, "\t\t\t<Points>\n");
       -        VTK_XML_VECTOR(pos, "Position [m]", "Float64", VTK_FLOAT_FMT);
       -        fprintf(stream, "\t\t\t</Points>\n");
       -
       -        fprintf(stream,
       -                "\t\t\t<Cells>\n"
       -                "\t\t\t\t<DataArray type=\"Int32\" Name=\"connectivity\" "
       -                "NumberOfComponents=\"1\" format=\"ascii\"/>\n"
       -                "\t\t\t\t<DataArray type=\"Int32\" Name=\"offsets\" "
       -                "NumberOfComponents=\"1\" format=\"ascii\"/>\n"
       -                "\t\t\t\t<DataArray type=\"UInt8\" Name=\"types\" "
       -                "NumberOfComponents=\"1\" format=\"ascii\"/>\n"
       -                "\t\t\t</Cells>\n");
       -
       -        fprintf(stream,
       -                "\t\t\t<PointData Scalars=\"Diameter [m]\" "
       -                "Vectors=\"Angular position [-]\">\n");
       -
       -        fprintf(stream,
       -                "\t\t\t\t<DataArray type=\"Float64\" Name=\"Diameter [m]\" "
       -                "NumberOfComponents=\"1\" format=\"ascii\">\n");
       -        for (i = 0; i < n; i++)
       -                fprintf(stream, VTK_FLOAT_FMT, grains[i].radius * 2.0);
       -        fprintf(stream,
       -                "\n"
       -                "\t\t\t\t</DataArray>\n");
       -
       -        VTK_XML_VECTOR(vel, "Velocity [m/s]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_VECTOR(acc, "Acceleration [m/s^2]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_VECTOR(force, "Force [N]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_VECTOR(angpos, "Angular position [rad]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_VECTOR(angvel, "Angular velocity [rad/s]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_VECTOR(angacc, "Angular acceleration [rad/s^2]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_VECTOR(torque, "Torque [N/m]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_VECTOR(disp, "Displacement [m]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_VECTOR(forceext, "External body force [N]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_SCALAR(density, "Density [kg/m^3]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_SCALAR(fixed, "Fixed [-]", "Int64", "%d ");
       -        VTK_XML_SCALAR(rotating, "Rotating [-]", "Int64", "%d ");
       -        VTK_XML_SCALAR(enabled, "Enabled [-]", "Int64", "%d ");
       -        VTK_XML_SCALAR(youngs_modulus, "Young's modulus [Pa]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_SCALAR(poissons_ratio, "Poisson's ratio [-]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_SCALAR(friction_coeff, "Friction coefficient [-]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_SCALAR(tensile_strength, "Tensile strength [Pa]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_SCALAR(shear_strength, "Shear strength [Pa]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_SCALAR(fracture_toughness, "Fracture toughness [Pa]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_VECTOR(gridpos, "Grid position [-]", "UInt64", "%zu ");
       -        VTK_XML_SCALAR(ncontacts, "Number of contacts [-]", "UInt64", "%zu ");
       -        VTK_XML_VECTOR(contact_stress, "Contact stress [Pa]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_SCALAR(thermal_energy, "Thermal energy [J]", "Float64", VTK_FLOAT_FMT);
       -        VTK_XML_SCALAR(color, "Color [-]", "Int64", "%d ");
       -
       -        fprintf(stream, "\t\t\t</PointData>\n");
       +        size_t i;
        
       -        fprintf(stream,
       -                "\t\t</Piece>\n"
       -                "\t</UnstructuredGrid>\n"
       -                "</VTKFile>\n");
       -}
       +        for (i = 0; i < sim->ng; i++) {
       +                /* sim_detect_contacts(sim); */
       +                /* sim_resolve_interactions(sim); */
       +                /* grain_temporal_integration(sim->g[i], sim->dt); */
       +        }
        
       -void
       -sim_print_grains_vtk(FILE *stream, const struct simulation *sim)
       -{
       -        print_grains_vtk(stream, sim->grains, sim->np);
       +        sim->t += sim->dt;
       +        sim->t_file += sim->dt;
       +        sim->iter++;
        }
        
        void
        sim_run_time_loop(struct simulation *sim)
        {
       +        if (sim->dt > sim->dt_file)
       +                sim->dt = sim->dt_file;
       +
       +        do {
       +                sim_step_time(sim);
       +        } while (sim->t - sim->dt < sim->t_end);
        }
 (DIR) diff --git a/simulation.h b/simulation.h
       t@@ -24,11 +24,13 @@ struct simulation {
                double t;
                double t_end;
                double dt;
       -        double file_dt;
       +        double dt_file;
                size_t n_file;
       +        size_t iter;
       +        double t_file;
        
                /* grains */
       -        size_t np;
       +        size_t ng;
                struct grain *grains;
        };
        
       t@@ -39,11 +41,10 @@ void sim_free(struct simulation *sim);
        void sim_add_grain(struct simulation *sim, struct grain *g);
        void sim_read_grains(struct simulation *sim, FILE *stream);
        
       -void print_grains(FILE *stream, const struct grain *grains, size_t n);
        void sim_print_grains(FILE *stream, const struct simulation *sim);
       -void print_grains_vtk(FILE *stream, const struct grain *grains, size_t n);
        void sim_print_grains_vtk(FILE *stream, const struct simulation *sim);
        
       +void sim_step_time(struct simulation *sim);
        void sim_run_time_loop(struct simulation *sim);
        
        #endif