Commit 2d1b1eb7 authored by Cassio Kirch's avatar Cassio Kirch

first commit

parents
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "define.h"
#include "boid.h"
unsigned int
getBoxID (struct Boid boid)
{
unsigned int tmpBoxID;
boid.position[X] += 0.50000*RANGE;
boid.position[Y] += 0.50000*RANGE;
tmpBoxID = (unsigned int) (boid.position[X] * \
((double) BOXES_IN_EDGE/RANGE)) + \
BOXES_IN_EDGE * (unsigned int) (boid.position[Y] * \
(double) BOXES_IN_EDGE/RANGE);
if (boid.position[X] == RANGE)
tmpBoxID -= 1;
if (boid.position[Y] == RANGE)
tmpBoxID -= BOXES_IN_EDGE;
return tmpBoxID;
}
void
initializeCircle (struct Boid* const boid)
{
double angle, radius;
angle = RANDOM_0(PI*2.0);
radius = sqrt(RANDOM_0(1.0)) * RANGE * PROP * 0.500;
boid -> position[X] = radius*cos(angle);
boid -> position[Y] = radius*sin(angle);
}
void
initializeRectangle (struct Boid* const boid)
{
boid -> position[X] = -PROP*RANGE*0.5 + RANDOM_0(PROP * RANGE);
boid -> position[Y] = -PROP*RANGE*0.5 + RANDOM_0(PROP * RANGE);
}
void
initializeBoid (struct Boid* const boid)
{
static unsigned int boidCount = 0;
static unsigned int endoBoids = (unsigned)(ENDOPROP*N);
initializeCircle (boid);
boid -> velocity[X] = 0.0;
boid -> velocity[Y] = 0.0;
boid -> boxID = getBoxID (*boid);
boid -> next = NULL;
boid -> previous = NULL;
boid -> gamma = 0.0;
if (boidCount < endoBoids)
boid -> type = ENDODERM;
else
boid -> type = ECTODERM;
++boidCount;
}
#ifndef DERM_TYPE_ENUM
typedef enum {ENDODERM, ECTODERM} cellType;
#define DERM_TYPE_ENUM
#endif
#ifndef BOID_STRUCT
struct Boid
{
struct Boid* next; /* Considering on the same box... */
struct Boid* previous; /* ... (same list). */
unsigned int boxID;
unsigned int neighbors;
unsigned int ectoNeighbors;
cellType type; /* Endo or ectoderm. */
double gamma;
double position[2];
double velocity[2];
double newVelocity[2];
};
#define BOID_STRUCT
#endif
unsigned int getBoxID (struct Boid boid);
void initializeCircle (struct Boid* const boid);
void initializeRectangle (struct Boid* const boid);
void initializeBoid (struct Boid* const boid);
#include <stdio.h>
#include <stdbool.h>
//#include <err.h>
#include "define.h"
#include "boid.h"
#include "box.h"
void
initializeBoxes (struct Box* const box) /* The input is an array. */
{
unsigned int boxCount;
for (boxCount=0; boxCount<BOXES; ++boxCount)
{
box[boxCount].boxID = boxCount;
box[boxCount].first = NULL;
box[boxCount].last = NULL;
}
}
bool
isInEdge (unsigned int boxID)
{
if (boxID < BOXES_IN_EDGE)
return true;
else if (boxID >= BOXES - BOXES_IN_EDGE)
return true;
else if (boxID % BOXES_IN_EDGE == 0)
return true;
else if (boxID % BOXES_IN_EDGE == BOXES_IN_EDGE - 1)
return true;
else
return false;
}
/* Append boid to the end of its list. */
void
appendBoid (struct Boid* const boid, struct Box* const box)
{
struct Boid* conductor;
boid -> boxID = getBoxID(*boid);
conductor = box[boid -> boxID].first;
if (conductor == NULL)
box[boid -> boxID].first = boid;
else
{
while (conductor -> next != 0)
conductor = conductor -> next;
conductor -> next = boid;
}
boid -> previous = conductor;
boid -> next = NULL;
box[boid -> boxID].last = boid;
}
void
removeBoid (struct Boid* const boid, struct Box* const box)
{
struct Boid* boidPointer;
if (boid == box[boid -> boxID].first)
box[boid -> boxID].first = boid -> next;
if (boid == box[boid -> boxID].last)
box[boid -> boxID].last = boid -> previous;
boidPointer = boid -> previous;
if (boidPointer != 0)
boidPointer -> next = boid -> next;
boidPointer = boid -> next;
if (boidPointer != 0)
boidPointer -> previous = boid -> previous;
}
unsigned int
getSouthBoxID (const unsigned int centralBoxID)
{
if (centralBoxID >= BOXES_IN_EDGE)
return (centralBoxID - BOXES_IN_EDGE);
else
return (BOXES - BOXES_IN_EDGE + centralBoxID);
}
unsigned int
getNorthBoxID (const unsigned int centralBoxID)
{
if (centralBoxID < BOXES - BOXES_IN_EDGE)
return (centralBoxID + BOXES_IN_EDGE);
else
return (centralBoxID + BOXES_IN_EDGE - BOXES);
}
unsigned int
getWestBoxID (const unsigned int centralBoxID)
{
if (centralBoxID % BOXES_IN_EDGE != 0)
return (centralBoxID - 1);
else
return (centralBoxID + BOXES_IN_EDGE - 1);
}
unsigned int
getEastBoxID (const unsigned int centralBoxID)
{
if (centralBoxID % BOXES_IN_EDGE != BOXES_IN_EDGE - 1)
return (centralBoxID + 1);
else
return (centralBoxID - BOXES_IN_EDGE + 1);
}
#include <stdbool.h>
#ifndef BOX_STRUCT
struct Box
{
unsigned int boxID;
struct Boid* first;
struct Boid* last;
};
#define BOX_STRUCT
#endif
void initializeBoxes (struct Box* const box);
void appendBoid (struct Boid* boid, struct Box* const box);
void removeBoid (struct Boid* boid, struct Box* const box);
unsigned int getSouthBoxID (const unsigned int centralBoxID);
unsigned int getNorthBoxID (const unsigned int centralBoxID);
unsigned int getWestBoxID (const unsigned int centralBoxID);
unsigned int getEastBoxID (const unsigned int centralBoxID);
bool isInEdge (unsigned int boxID);
#if (BOXES_IN_EDGE*BOXES_IN_EDGE)!=BOXES
#error : BOXES_IN_EDGE**2 must be equal to BOXES.
#endif
/*#if (BOXES_IN_EGDE*NEIGHBOUR_DISTANCE)>(2*RANGE)
#error : Only the 8 neighbour boxes can contain neighbour boids.
#endif
*/
int
checkParameters()
{
if (DIMENSIONS != 2)
return 10;
if (NEIGHBOR_DISTANCE < ELASTIC_DISTANCE)
return -10;
if (ELASTIC_DISTANCE < CORE_RADIUS)
return -11;
if (V0*DT > CORE_RADIUS)
return -19;
else
return 0;
}
#define X 0
#define Y 1
#define Z 2
#define PI 3.14159265359
//#define PLOT_EXIT_FILES
#define GAMMA_FILE
#define ENDOPROP 0.25
#define N 2000
#define RANGE 200.0
#define PROP 0.35
#define SEED 13049 /* Must be unsigned int. */
#define DIMENSIONS 2
#define BOXES 22500
#define BOXES_IN_EDGE 150
#define NEIGHBOR_DISTANCE 1.32 /* r0 */
#define CORE_RADIUS 0.4 /* rc */
#define ELASTIC_DISTANCE 1.0 /* re */
#define INFINITE_FORCE 10000.0
#define ETA 1.0
#define V0 0.1
#define ALPHA11 0.01
#define ALPHA12 0.01
#define ALPHA22 0.01
#define BETA11 3.83
#define BETA12 2.53
#define BETA22 2.50
#define DT 1.0
#define DATE_SIZE 24 /* Size of string */
#define FILENAME_SIZE 128 /* Size of string */
#define STEPS 10000000000
#define EXIT_INTERVAL 1000
#define NUM_THREADS 8
#define RANDOM_0(max) (((double)rand()/RAND_MAX) * max)
#define absDistance(x1, y1, x2, y2) sqrt(pow ((x2) - (x1), 2) + \
pow ((y2) - (y1), 2))
#include <math.h>
#include <stdbool.h>
#include "define.h"
#include "box.h"
#include "boid.h"
#include "distance.h"
struct Distance
getDistance (const struct Boid* const boid2, const struct Boid* const boid1)
{
struct Distance distance;
double temp;
/* Distância normal. */
distance.module = absDistance (boid1 -> position[X], \
boid1 -> position[Y], \
boid2 -> position[X], \
boid2 -> position[Y]);
if (distance.module < NEIGHBOR_DISTANCE)
{
distance.sine = (boid2 -> position[Y] - boid1 -> position[Y]) / \
distance.module;
distance.cosine = (boid2 -> position[X] - boid1 -> position[X]) / \
distance.module;
return distance;
}
else if ((isInEdge (boid1 -> boxID) == true) || \
(isInEdge (boid2 -> boxID) == true))
{
/* 1o boid à esquerda. */
temp = boid1 -> position[X] + RANGE;
distance.module = absDistance (temp, \
boid1 -> position[Y], \
boid2 -> position[X], \
boid2 -> position[Y]);
if (distance.module < NEIGHBOR_DISTANCE)
{
distance.sine = (boid2 -> position[Y] - boid1 -> position[Y]) / \
distance.module;
distance.cosine = (boid2 -> position[X] - temp) / distance.module;
return distance;
}
/* 1o boid à direita */
temp = boid1 -> position[X] - RANGE;
distance.module = absDistance (temp, \
boid1 -> position[Y], \
boid2 -> position[X], \
boid2 -> position[Y]);
if (distance.module < NEIGHBOR_DISTANCE)
{
distance.sine = (boid2 -> position[Y] - boid1 -> position[Y]) / \
distance.module;
distance.cosine = (boid2 -> position[X] - temp) / distance.module;
return distance;
}
/* 1o boid abaixo */
temp = boid1 -> position[Y] + RANGE;
distance.module = absDistance (boid1 -> position[X], \
temp, \
boid2 -> position[X], \
boid2 -> position[Y]);
if (distance.module < NEIGHBOR_DISTANCE)
{
distance.sine = (boid2 -> position[Y] - temp) / distance.module;
distance.cosine = (boid2 -> position[X] - boid1 -> position[X]) / \
distance.module;
return distance;
}
/* 1o boid acima */
temp = boid1 -> position[Y] - RANGE;
distance.module = absDistance (boid1 -> position[X], \
temp, \
boid2 -> position[X], \
boid2 -> position[Y]);
if (distance.module < NEIGHBOR_DISTANCE)
{
distance.sine = (boid2 -> position[Y] - temp) / distance.module;
distance.cosine = (boid2 -> position[X] - boid1 -> position[X]) / \
distance.module;
return distance;
}
/* não vizinhos */
else
{
distance.sine = 0.0;
distance.cosine = 0.0;
distance.module = NEIGHBOR_DISTANCE * 2.0;
return distance;
}
}
/* não vizinhos */
else
{
distance.sine = 0.0;
distance.cosine = 0.0;
distance.module = NEIGHBOR_DISTANCE * 2.0;
return distance;
}
}
/* Supõe-se que a função só será chamada para boids vizinhos. */
double
getForce (const double distance)
{
if (distance < CORE_RADIUS)
return INFINITE_FORCE;
//else if (distance > NEIGHBOR_DISTANCE)
//return 0.0;
else
return (1.000000 - (distance / ELASTIC_DISTANCE));
}
#ifndef DISTANCE_STRUCT
struct Distance
{
double module;
double sine;
double cosine;
};
#define DISTANCE_STRUCT
#endif
double getForce (const double distance);
struct Distance getDistance (const struct Boid* const boid2, \
const struct Boid* const boid1);
#include <stdio.h>
#include <time.h>
#include <sys/stat.h>
#include "define.h"
#include "files.h"
char dateString[DATE_SIZE];
char*
getDate()
{
return dateString;
}
void
setDate ()
{
time_t now = time(NULL);
struct tm* tmP = localtime(&now);
strftime(dateString, DATE_SIZE, "%Y%m%d_%Hh%Mmin%ssec", tmP);
}
void
putParametersToFile(FILE* file)
{
fprintf(file, "#Modo 'um sistema' com alpha e beta fixos.\n" \
"#%s\n" \
"#ALPHA11 = %lf\n" \
"#ALPHA12 = %lf\n" \
"#ALPHA22 = %lf\n" \
"#BETA11 = %lf\n" \
"#BETA12 = %lf\n" \
"#BETA22 = %lf\n" \
"#N = %d\n" \
"#ETA = %lf\n" \
"#V0 = %lf\n" \
"#R0 = %lf\n" \
"#DT = %lf\n" \
"#CORE_RADIUS = %lf\n" \
"#ELASTIC_DISTANCE = %lf\n" \
"#INFINITE FORCE = %lf\n" \
"#RANGE = %lf\n" \
"#PROP = %lf\n" \
"#SEED = %u\n" \
"#BOXES = %u\n" \
"#BOXES_IN_EDGE = %u\n" \
"#STEPS = %u\n" \
"#EXIT_INTERVAL = %d\n#\n" \
"#THREADS: %u\n", dateString, ALPHA11, \
ALPHA12, ALPHA22, \
BETA11, BETA12, BETA22, \
(unsigned int)N, ETA, V0, \
NEIGHBOR_DISTANCE, DT, \
CORE_RADIUS, ELASTIC_DISTANCE, \
INFINITE_FORCE, RANGE, \
PROP, (unsigned int)SEED, (unsigned int)BOXES, \
(unsigned int)BOXES_IN_EDGE, (unsigned int)STEPS, \
(unsigned int)EXIT_INTERVAL, (unsigned int)NUM_THREADS);
}
FILE*
initializeSingleFile ()
{
FILE* dat;
char fileName[FILENAME_SIZE];
sprintf(fileName, "one_system_%s.dat", dateString);
dat = fopen (fileName, "w");
putParametersToFile(dat);
return dat;
}
FILE*
initializeGammaFile()
{
FILE* gammaFile;
char dateString[DATE_SIZE+1], fileName[FILENAME_SIZE];
sprintf(fileName, "gamma_%s.dat", dateString);
gammaFile = fopen(fileName, "w");
putParametersToFile(gammaFile);
return gammaFile;
}
FILE*
initializeStepAndTypeFile (const unsigned long long int step, \
const cellType type)
{
static FILE* dat;
static char fileName[FILENAME_SIZE];
sprintf(fileName, "%llu_%u.dat", step, type);
dat = fopen(fileName, "w");
return dat;
}
#include "boid.h"
FILE* initializeSingleFile ();
FILE* initializeStepAndTypeFile (const unsigned long long int step, \
const cellType type);
FILE* initializeGammaFile();
void setDate ();
char* getDate();
#include <stdio.h>
#include <stdlib.h>
#include "define.h"
#include "comp_errors.h"
#include "files.h"
#include "boid.h"
#include "box.h"
#include "nextstep.h"
#include "main.h"
double
getGamma (const struct Boid* const boids)
{
double gamma = 0.0;
unsigned int endoBoids = (unsigned)(ENDOPROP*N);
unsigned int boidCount = endoBoids;
do
{
boidCount--;
gamma += boids[boidCount].gamma/endoBoids;
}
while(boidCount != 0);
return gamma;
}
void
one_system ()
{
struct Boid boid[N];
struct Box box[BOXES];
unsigned int boidCount, boxID, threadCount;
unsigned long long int step, continuousStep = 0;
FILE* dat = initializeSingleFile();
FILE* finalConfigurationFile;
char finalConfigurationFileName[FILENAME_SIZE];
#ifdef PLOT_EXIT_FILES
FILE* endoFile;
FILE* ectoFile;
#endif
FILE* gammaFile = initializeGammaFile();
/* Set the pthread_create parameters. */
struct Parameters parametersStruct[NUM_THREADS];
void* parameters[NUM_THREADS];
for (threadCount=0; threadCount<NUM_THREADS; ++threadCount)
{
parametersStruct[threadCount].threadID = threadCount;
parametersStruct[threadCount].left = threadCount * \
((int)N/(int)NUM_THREADS);
parametersStruct[threadCount].right = (threadCount+1) * \
((int)N/(int)NUM_THREADS);
if (threadCount == NUM_THREADS-1)
parametersStruct[threadCount].right += N%NUM_THREADS;
parametersStruct[threadCount].boid = boid;
parametersStruct[threadCount].box = box;
parameters[threadCount] = &(parametersStruct[threadCount]);
}
for (boidCount=0; boidCount<N; boidCount++)
{
initializeBoid (&boid[boidCount]);
checkLimits(&(boid[boidCount]));
}
initializeBoxes(box);
for (boidCount=0; boidCount<N; boidCount++)
appendBoid(&(boid[boidCount]), box);
for (step=0; step<STEPS; ++step)
{
nextStep(parameters);
for (boidCount=0; boidCount<N; boidCount++)
{
checkLimits(&(boid[boidCount]));
boxID = getBoxID(boid[boidCount]);
if (boid[boidCount].boxID != boxID)
{
removeBoid(&boid[boidCount], box);
appendBoid(&boid[boidCount], box);
}
}
if (step%EXIT_INTERVAL == 0)
{
#ifdef PLOT_EXIT_FILES
endoFile = initializeStepAndTypeFile (continuousStep, ENDODERM);
ectoFile = initializeStepAndTypeFile (continuousStep, ECTODERM);
for (boidCount=0; boidCount<N; boidCount++)
{
/*fprintf(dat, "%llu\t" "%llu\t" "%u\t" "%u\t" "%lf\t" "%lf\n", \
continuousStep, step, boid[boidCount].type, \
boidCount, boid[boidCount].position[X], \
boid[boidCount].position[Y]);
*/
if (boid[boidCount].type == ENDODERM)
fprintf(endoFile, "%u\t" "%u\t" "%lf\t" "%lf\n", \
boid[boidCount].type, \
boidCount, boid[boidCount].position[X], \
boid[boidCount].position[Y]);
else
fprintf(ectoFile, "%u\t" "%u\t" "%lf\t" "%lf\n", \
boid[boidCount].type, \
boidCount, boid[boidCount].position[X], \
boid[boidCount].position[Y]);
}
fclose(endoFile);
fclose(ectoFile);
#endif
fprintf(gammaFile, "%llu\t%lf\n", step, getGamma(boid));
++continuousStep;
}
}
fclose (gammaFile);
fclose (dat);
sprintf(finalConfigurationFileName, "final_%s.dat", getDate());
finalConfigurationFile = fopen(finalConfigurationFileName, "w");
for (boidCount=0; boidCount<N; boidCount++)
{
fprintf(finalConfigurationFile, \
"%u\t" "%u\t" "%lf\t" "%lf\n", \
boid[boidCount].type, boidCount, \
boid[boidCount].position[X], \
boid[boidCount].position[Y]);
}
fclose (finalConfigurationFile);
}
int
main ()
{
int returned = checkParameters();
if (returned != 0)
return returned;
setDate();
srand(time(NULL));
one_system();
return 0;
}
#include "boid.h"
#include "box.h"
#ifndef PARAMETERS_STRUCT
struct Parameters
{
unsigned int threadID;
unsigned int left;
unsigned int right;
struct Boid* boid;
struct Boid* firstEndoBoid;
struct Box* box;
};
#define PARAMETERS_STRUCT
#endif
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>
#include <err.h>
#include "define.h"
#include "boid.h"
#include "box.h"
#include "distance.h"
#include "nextstep.h"
void
checkLimits (struct Boid* const boid)
{
while (boid -> position[X] < -(RANGE / 2.0))
boid -> position[X] += RANGE;
while (boid -> position[X] > (RANGE / 2.0))
boid -> position[X] -= RANGE;
while (boid -> position[Y] < -(RANGE / 2.0))
boid -> position[Y] += RANGE;
while (boid -> position[Y] > (RANGE / 2.0))
boid -> position[Y] -= RANGE;
}
void
makeSum (struct Boid* const boid, const struct Boid* conductor, \
double* const sumFX, double* const sumFY, \
double* const sumVX, double* const sumVY)
{
struct Distance distance;
double force;
while (conductor != NULL)
{
if (boid != conductor)
{
distance = getDistance (boid, conductor);
if (distance.module < NEIGHBOR_DISTANCE)
{
boid -> neighbors++;
force = getForce (distance.module);
if (conductor -> type == ENDODERM && \
boid -> type == ENDODERM)
{
*sumFX += BETA11 * force * distance.cosine;
*sumFY += BETA11 * force * distance.sine;
*sumVX += ALPHA11 * conductor -> velocity[X] / V0;
*sumVY += ALPHA11 * conductor -> velocity[Y] / V0;
}
else if (conductor -> type == ECTODERM && \
boid -> type == ECTODERM)
{
*sumFX += BETA22 * force * distance.cosine;
*sumFY += BETA22 * force * distance.sine;
*sumVX += ALPHA22 * conductor -> velocity[X] / V0;
*sumVY += ALPHA22 * conductor -> velocity[Y] / V0;
}
else
{
*sumFX += BETA12 * force * distance.cosine;
*sumFY += BETA12 * force * distance.sine;
*sumVX += ALPHA12 * conductor -> velocity[X] / V0;
*sumVY += ALPHA12 * conductor -> velocity[Y] / V0;