/******************************************************************************\ Copyright 1995 The University of North Carolina at Chapel Hill. All Rights Reserved. Permission to use, copy, modify and distribute this software and its documentation for educational, research and non-profit purposes, without fee, and without a written agreement is hereby granted, provided that the above copyright notice and the following three paragraphs appear in all copies. IN NO EVENT SHALL THE UNIVERSITY OF NORTH CAROLINA AT CHAPEL HILL BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF NORTH CAROLINA HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES. Permission to use, copy, modify and distribute this software and its documentation for educational, research and non-profit purposes, without fee, and without a written agreement is hereby granted, provided that the above copyright notice and the following three paragraphs appear in all copies. THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION S TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. The authors may be contacted via: US Mail: J. Cohen/M. Lin/D. Manocha/K. Ponamgi Department of Computer Science Sitterson Hall, CB #3175 University of N. Carolina Chapel Hill, NC 27599-3175 Phone: (919)962-1749 EMail: {cohenj,lin,manocha,ponamgi}@cs.unc.edu \*****************************************************************************/ /*****************************************************************************\ nbody.c -- Description : Nbody application to exercise libcollide. This simulation may be run with or without graphics support. \*****************************************************************************/ /*----------------------------- Local Includes -----------------------------*/ #include #include #include #include #include #include "polytope.h" #include "graphics.h" /*----------------------------- Local Constants -----------------------------*/ #if 0 #define PRINT_COLLISIONS #endif /*------------------------------ Local Macros -------------------------------*/ #define ABS(a) (((a) > 0.0) ? (a) : (-a)) /*------------------------------- Local Types -------------------------------*/ /*------------------------ Local Function Prototypes ------------------------*/ /*------------------------------ Local Globals ------------------------------*/ Status status; Polytope *poly=0; int frames=0; int num_objects=0; int total_collisions=0; double cube_bound=0.0; extern __col_Polyhedron __col_polyhedronLibrary[1000]; /*---------------------------------Functions-------------------------------- */ /*****************************************************************************\ @ polygon_area() ----------------------------------------------------------------------------- description : Compute the signed area of the 2D projection of a face input : A face, and which two dimensions to use for the area computation output : area of the projection of the face notes : \*****************************************************************************/ double polygon_area(__col_Face *face, int axis1, int axis2) { double area; col_Fnode *feature; __col_Vertex *v1, *v2; area = 0.0; for (feature=face->verts, v1 = feature->f.v; feature; feature = feature->next, v1 = v2) { v2 = (feature->next) ? feature->next->f.v : face->verts->f.v; area += v1->coords[axis1]*v2->coords[axis2] - v1->coords[axis2]*v2->coords[axis1]; } area /= 2.0; return area; } /** End of polygon_area() **/ /*****************************************************************************\ @ face_area() ----------------------------------------------------------------------------- description : Compute the unsigned area of a face input : output : notes : \*****************************************************************************/ double face_area(__col_Face *face) { double area, xy_area, yz_area, xz_area; /* project the polygon onto the three coordinate planes and calculate its area there */ xy_area = polygon_area(face, 0, 1); yz_area = polygon_area(face, 1, 2); xz_area = polygon_area(face, 0, 2); /* use the projected areas to calculate the area of the face */ area = sqrt(xy_area * xy_area + yz_area * yz_area + xz_area * xz_area); return area; } /** End of face_area **/ /*****************************************************************************\ @ polytope_volume() ----------------------------------------------------------------------------- description : Compute the volume of a polyhedron input : output : notes : \*****************************************************************************/ double polytope_volume(__col_Polyhedron *poly) { double volume, area; col_Fnode *feature; volume = 0.0; for(feature=poly->faces; feature; feature = feature->next) { area = face_area(feature->f.f); volume -= area * feature->f.f->plane[3]; } volume /= 3.0; return volume; } /** End of polytope_volume() **/ /*****************************************************************************\ @ init_movement() ----------------------------------------------------------------------------- description : Initializes the simulation volume and the velocity/rotation of each polytope. input : output : notes : \*****************************************************************************/ void init_movement(Status *status, Polytope *poly) { int i; int coord; int x_ratio, y_ratio, z_ratio, tot_ratio; double v_squared; double delta_vel, real_velocity, percent_off; double x_percent, y_percent, z_percent; double poly_volume; int x_sign, y_sign, z_sign; double x_rot, y_rot, z_rot, delta_rot; int num_div_lib, num_mod_lib; double lib_volume, mod_volume, vol; float initial_orientation[3]; double max_radius; int *unused_cells; int num_unused_cells, cells_across, num_cells; double cell_radius; int cell_index, cell_number; int cell_x_number, cell_y_number, cell_z_number; num_div_lib = status->num_polytopes / status->num_lib_polytopes; num_mod_lib = status->num_polytopes % status->num_lib_polytopes; lib_volume = mod_volume = 0.0; for (i = 0; i < status->num_lib_polytopes; i++) { vol = polytope_volume(&(__col_polyhedronLibrary[i])); lib_volume += vol; if (i <= num_mod_lib) mod_volume += vol; } poly_volume = num_div_lib * lib_volume + num_mod_lib * mod_volume; status->working_volume = ((double)poly_volume/((double)status->density/100.0)); cube_bound = (double)(pow((status->working_volume), (1.0/3.0)) / 2.0); printf("Density: %.3f%%\n", status->density); printf("Polytope volume: %.3f\n", poly_volume); printf("Working Volume: %.3f\n", status->working_volume); printf("Cube length: %.3f\n", cube_bound * 2.0); printf("Velocity: %d%% object radius per frame\n", status->trans_velocity); printf("Rotational Velocity: %d degrees per frame\n", status->rot_velocity); fflush(stdout); printf("Initializing polytope positions and velocities..."); fflush(stdout); max_radius = 0.0; for (i=0; i < status->num_polytopes; i++) { col_calc_cuboid(poly[i].id, poly[i].object_center, &poly[i].max_radius); poly[i].old_colliding = poly[i].colliding = -1; if (poly[i].max_radius > max_radius) max_radius = poly[i].max_radius; } cells_across = ceil(pow(status->num_polytopes, (1.0/3.0))); num_cells = cells_across * cells_across * cells_across; cell_radius = cube_bound/cells_across; if (cell_radius < max_radius) { #if 0 fprintf(stderr, "init_movement: Some objects may initially colliding.\n"); fflush(stderr); #endif } ALLOCATE(unused_cells, int, num_cells); for (i=0; i < num_cells; i++) unused_cells[i] = i; num_unused_cells = num_cells; for (i= 0; i < status->num_polytopes; ++i) { x_rot = (double) ((rand() % 100) - 50); y_rot = (double) ((rand() % 100) - 50); z_rot = (double) ((rand() % 100) - 50); delta_rot = status->rot_velocity * Q_PI/180.; q_make(poly[i].delta_quat, x_rot, y_rot, z_rot, delta_rot); /* V = sqrt(Vx*Vx + Vy*Vy + Vz*Vz); choose some ratio of Vx:Vy:Vz, and use this to find Vx, Vy, and Vz that come to the proper total velocity Actually, pick ratio of the velocities squared, and calculate percentage of the squared total velocity */ v_squared = status->trans_velocity * poly[i].max_radius / 100.0; v_squared *= v_squared; x_ratio = rand() % 100; y_ratio = rand() % 100; z_ratio = rand() % 100; tot_ratio = x_ratio + y_ratio + z_ratio; x_percent = (double)x_ratio/(double)tot_ratio; y_percent = (double)y_ratio/(double)tot_ratio; z_percent = (double)z_ratio/(double)tot_ratio; poly[i].trans_vel[X] = (float)sqrt(x_percent*v_squared); poly[i].trans_vel[Y] = (float)sqrt(y_percent*v_squared); poly[i].trans_vel[Z] = (float)sqrt(z_percent*v_squared); x_sign = rand() % 100; y_sign = rand() % 100; z_sign = rand() % 100; if (x_sign < 50) poly[i].trans_vel[X] *= -1.0; if (y_sign < 50) poly[i].trans_vel[Y] *= -1.0; if (z_sign < 50) poly[i].trans_vel[Z] *= -1.0; initial_orientation[X] = rand() % 360; initial_orientation[Y] = rand() % 360; initial_orientation[Z] = rand() % 360; q_from_euler(poly[i].tot_quat, initial_orientation[X], initial_orientation[Y], initial_orientation[Z]); /* try to choose initial positions so objects aren't colliding */ cell_index = rand() % num_unused_cells; cell_number = unused_cells[cell_index]; unused_cells[cell_index] = unused_cells[num_unused_cells-1]; num_unused_cells--; cell_x_number = cell_number % cells_across; cell_y_number = (cell_number % (cells_across*cells_across)) / cells_across; cell_z_number = cell_number / (cells_across*cells_across); poly[i].trans_tot[X] = -cube_bound + (cell_x_number*2 + 1)*cell_radius; poly[i].trans_tot[Y] = -cube_bound + (cell_y_number*2 + 1)*cell_radius; poly[i].trans_tot[Z] = -cube_bound + (cell_z_number*2 + 1)*cell_radius; } free(unused_cells); printf("done.\n"); fflush(stdout); } /** End of init_movement() **/ /*****************************************************************************\ @ mat_identity() ----------------------------------------------------------------------------- description : input : output : notes : \*****************************************************************************/ static void mat_identity(col_Mat4 m) { register int i,j; for (i = 0; i < 4; i++) for(j = 0; j < 4; j++) m[j][i] = (i==j) ? 1.0 : 0.0; } /** End of mat_identity() **/ /*****************************************************************************\ @ calc_translation() ----------------------------------------------------------------------------- description : input : output : notes : \*****************************************************************************/ static void calc_translation(col_Mat4 m, float x, float y, float z) { mat_identity(m); m[0][3] = x; m[1][3] = y; m[2][3] = z; } /** End of calc_translation() **/ /*****************************************************************************\ @ load_all_lib_objects() ----------------------------------------------------------------------------- description : Loads ALL objects in collision library for collision detection. input : output : notes : \*****************************************************************************/ void load_all_lib_objects(Status *status, Polytope **poly) { int i; graphics_per_library_object_init(status); ALLOCATE(*poly, Polytope, status->num_polytopes); printf("Adding %d polytopes...", status->num_polytopes); fflush(stdout); for (i=0; i < status->num_polytopes; i++) { col_instance_polytope_by_num(i%(status->num_lib_polytopes), &((*poly)[i].id)); } printf("done.\n"); fflush(stdout); graphics_per_object_instance_init(status); } /** End of load_all_lib_objects() **/ /*****************************************************************************\ @ init_collision() ----------------------------------------------------------------------------- description : Read in the polytopes from a file. 1. Read in the object. 2. Pass the objects to collision library so it can create Voronoi cones. 3. Create initial positions for the objects. 4. Give initial velocities for the objects. 5. Give initial O(n^2) matrix the relationship between objects (bounce, stick, etc). input : output : notes : \*****************************************************************************/ void init_collision(Status *status, Polytope **poly, char *objFile) { int i; printf("Loading library file: %s ...", objFile); fflush(stdout); if (col_open(status->num_polytopes, objFile, 1.0, &(status->num_lib_polytopes)) != COL_OK) { fprintf(stderr, "init_collision: couldn't open collision library\n"); exit(1); } printf("done.\n"); fflush(stdout); load_all_lib_objects(status, poly); /* different libcollide matrix allocation methods might have different effects on the random number generator, so seed the generator again after initializing the library */ if (status->seed > 0) { srand(status->seed); } init_movement(status, *poly); printf("Activating all pairs (slow for COL_SPARSE_MATRIX representation)..."); fflush(stdout); col_activate_all(); col_enable_nbody(); printf("done.\n"); fflush(stdout); if (status->use_cuboid) { printf("Using fixed size cube for bounding box\n"); fflush(stdout); for (i=0; i < status->num_polytopes; i++) { col_use_cuboid((*poly)[i].id); } } } /** End of init_collision() **/ /*****************************************************************************\ @ usage() ----------------------------------------------------------------------------- description : input : output : notes : \*****************************************************************************/ void usage() { printf("usage: nbody [-l ] [-f #frames] [-n #polytopes]\n"); printf(" [-v velocity] [-r rotation] [-d density] [-s seed] [-b]\n"); printf(" velocity is %% of object radius per frame\n"); printf(" rotation is degrees of rotation per frame\n"); printf(" density is total object volume/simulation volume\n"); printf(" seed is an integer to seed the random number generator\n"); printf(" -b specifies using dynamic bounding boxes rather than\n"); printf(" fixed-sized bounding cubes\n"); printf("\n"); exit(1); } /** End of usage() **/ /*****************************************************************************\ @ initialize() ----------------------------------------------------------------------------- description : Control initialization. 1. Decide which file the polytopes will be read from. 2. Give initial signal to the collision library. 3. Read in polytopes. 4. Update State. state Updates the state of the simulation (whether do n-body and etc.). argc/argv Variables passed from the environment input : output : notes : \*****************************************************************************/ void initialize(int argc, char **argv, Status *status, Polytope **poly) { char objFile[256]; extern char *optarg; int c; int args_found; status->num_polytopes = DEFAULT_NUM_OBJECTS; status->trans_velocity = DEFAULT_VELOCITY; status->rot_velocity = DEFAULT_ROTATION; status->density = DEFAULT_DENSITY; status->seed = 0; status->use_cuboid = 1; strcpy(objFile, DEFAULT_LIBRARY); status->num_frames = DEFAULT_FRAMES; args_found = 1; while ((c = getopt(argc, argv, "bBhHl:L:f:F:n:N:v:V:r:R:d:D:s:S:")) != EOF) /* Parse the various options. */ switch (c) { case 'l': case 'L': strcpy(objFile, optarg); args_found += 2; break; case 'n': case 'N': status->num_polytopes = atoi(optarg); args_found += 2; break; case 'v': case 'V': status->trans_velocity = atoi(optarg); args_found += 2; break; case 'r': case 'R': status->rot_velocity = atoi(optarg); args_found += 2; break; case 'b': case 'B': status->use_cuboid = 0; args_found++; break; case 'f': case 'F': status->num_frames = atoi(optarg); args_found += 2; break; case 'd': case 'D': status->density = atof(optarg); args_found += 2; break; case 's': case 'S': status->seed = atoi(optarg); args_found += 2; break; case 'h': case 'H': args_found++; default: usage(); break; } if (args_found != argc) usage(); if (status->seed > 0) { srand(status->seed); printf("Random number seed: %d\n", status->seed); fflush(stdout); } /* allocate space for reporting all n^2 possible collisions, there aren't likely to be anywhere near this many */ ALLOCATE(status->all_collisions, col_Report, status->num_polytopes * status->num_polytopes); graphics_init(argv[0]); init_collision(status, poly, objFile); } /** End of initialize() **/ /*****************************************************************************\ @ compute_transformation() ----------------------------------------------------------------------------- description : Compute a polytope's current matrix transformation from its total translation and rotation input : output : notes : \*****************************************************************************/ static void compute_transformation(Polytope *poly) { col_Mat4 xrot, yrot, zrot, trans; col_Mat4 trans_to_center, xrot_cent, zy, zyxt2c; col_Mat4 rot, rot_cent; calc_translation(trans_to_center, -poly->object_center[X], -poly->object_center[Y], -poly->object_center[Z]); calc_translation(trans, poly->trans_tot[X] + poly->object_center[X], poly->trans_tot[Y] + poly->object_center[Y], poly->trans_tot[Z] + poly->object_center[Z]); q_to_row_matrix(rot, poly->tot_quat); __col_matMultXform(rot, trans_to_center, rot_cent); __col_matMultXform(trans, rot_cent, poly->tot); } /** End of compute_transformation() **/ /*****************************************************************************\ @ compute_collision() ----------------------------------------------------------------------------- description : Compute all object transformations and perform collision test. Swap translations and rotations of colliding pairs input : output : notes : \*****************************************************************************/ void compute_collision(Status *status, Polytope *poly) { col_Mat4 T; double angle[3]; int i, j, obj1, obj2; float temp; q_type quat_temp; /* save previous status */ for (i=0; i < status->num_polytopes; i++) { poly[i].old_colliding = poly[i].colliding; poly[i].colliding = -1; } /* inform collision library of new object positions */ for (i= 0; i < status->num_polytopes; ++i) { compute_transformation(&(poly[i])); col_set_transform(poly[i].id, poly[i].tot); } /* perform collision test */ if (col_test() != COL_OK) /* The biggie */ { fprintf(stderr, "col_test didn't return COL_OK\n"); exit(1); } /* inquire results of collision test */ status->num_collisions = col_report_collision_all(status->all_collisions, status->num_polytopes*status->num_polytopes); total_collisions += status->num_collisions; #ifdef PRINT_COLLISIONS if (status->num_collisions) printf("FRAME %d: %d collisions\n", frames, status->num_collisions); #endif /* exchange velocities of colliding objects, providing they weren't already colliding last frame (i.e. don't keep swapping until objects have seperated) */ for (i= 0; i < status->num_collisions; ++i) { obj1 = status->all_collisions[i].ids[0]; obj2 = status->all_collisions[i].ids[1]; #ifdef PRINT_COLLISIONS printf("ids: %d %d\n", obj1, obj2); #endif poly[obj1].colliding = obj2; poly[obj2].colliding = obj1; if ((poly[obj1].old_colliding != obj2) && (poly[obj2].old_colliding != obj1)) { temp = poly[obj1].trans_vel[X]; poly[obj1].trans_vel[X] = poly[obj2].trans_vel[X]; poly[obj2].trans_vel[X] = temp; temp = poly[obj1].trans_vel[Y]; poly[obj1].trans_vel[Y] = poly[obj2].trans_vel[Y]; poly[obj2].trans_vel[Y] = temp; temp = poly[obj1].trans_vel[Z]; poly[obj1].trans_vel[Z] = poly[obj2].trans_vel[Z]; poly[obj2].trans_vel[Z] = temp; for (j=0; j<3; j++) quat_temp[j] = poly[obj1].delta_quat[j]; for (j=0; j<3; j++) poly[obj1].delta_quat[j] = poly[obj2].delta_quat[j]; for (j=0; j<3; j++) poly[obj2].delta_quat[j] = quat_temp[j]; } } #ifdef PRINT_COLLISIONS if (status->num_collisions) printf("\n"); #endif } /** End of compute_collision() **/ /*****************************************************************************\ @ move_polytopes() ----------------------------------------------------------------------------- description : Moves each polytope by amount described in its rotational and translational velocities. Checks to see if move exceeds bounding box. input : output : notes : \*****************************************************************************/ void move_polytopes(Status *status, Polytope *poly) { int i; float trans; float center[3]; /* accumulate the changes to the translation and rotation -- this has some error accumulation over a long period of time, but we don't really care for this application */ for (i= 0; i < status->num_polytopes; ++i) { q_mult(poly[i].tot_quat, poly[i].tot_quat, poly[i].delta_quat); poly[i].trans_tot[X] += poly[i].trans_vel[X]; poly[i].trans_tot[Y] += poly[i].trans_vel[Y]; poly[i].trans_tot[Z] += poly[i].trans_vel[Z]; center[X] = poly[i].object_center[X] + poly[i].trans_tot[X]; center[Y] = poly[i].object_center[Y] + poly[i].trans_tot[Y]; center[Z] = poly[i].object_center[Z] + poly[i].trans_tot[Z]; /* Check for penetration of the walls of the simulation volume, and rebound off the walls if necessary. This is a bit of a hack, since it only keeps the object center within the walls. */ if (center[X] > cube_bound) { trans = center[X] - cube_bound; poly[i].trans_tot[X] -= trans; poly[i].trans_vel[X] = -poly[i].trans_vel[X]; } if (center[X] < -cube_bound) { trans = center[X] - (-cube_bound); poly[i].trans_tot[X] -= trans; poly[i].trans_vel[X] = -poly[i].trans_vel[X]; } if (center[Y] > cube_bound) { trans = center[Y] - cube_bound; poly[i].trans_tot[Y] -= trans; poly[i].trans_vel[Y] = -poly[i].trans_vel[Y]; } if (center[Y] < -cube_bound) { trans = center[Y] - (-cube_bound); poly[i].trans_tot[Y] -= trans; poly[i].trans_vel[Y] = -poly[i].trans_vel[Y]; } if (center[Z] > cube_bound) { trans = center[Z] - cube_bound; poly[i].trans_tot[Z] -= trans; poly[i].trans_vel[Z] = -poly[i].trans_vel[Z]; } if (center[Z] < -cube_bound) { trans = center[Z] - (-cube_bound); poly[i].trans_tot[Z] -= trans; poly[i].trans_vel[Z] = -poly[i].trans_vel[Z]; } } } /** End of move_polytopes() **/ /*****************************************************************************\ @ main() ----------------------------------------------------------------------------- description : main program input : output : notes : \*****************************************************************************/ int main(int argc, char **argv) { initialize(argc, argv, &status, &poly); printf("Simulating %d frames...", status.num_frames); fflush(stdout); for (frames=0; frames < status.num_frames; frames++) { move_polytopes(&status, poly); graphics_draw_scene(&status, poly); compute_collision(&status, poly); } printf("done.\n"); printf("Total collisions: %d\n", total_collisions); printf("Average collisions per frame: %.3f\n", (double)total_collisions / (double)status.num_frames); return 0; } /** End of main() **/ /*****************************************************************************\ \*****************************************************************************/