void printBasiliskDataStructure (struct BasiliskDataStructure * b) {
printf("Grains nparticles = %lu\n", b->particlessize);
/* Geometric aspect */
printf("Grains radius = %f\n", b->rayon);
printf("Grains center.x = %f\n", b->centreX);
printf("Grains center.y = %f\n", b->centreY);
printf("Grains center.z = %f\n", b->centreZ);
printf("Grains ncorners = %d\n", b->ncorners);
/* Coordinates of the corners, number of faces and their numerotation */
printf("Grains allPoints = %lu\n", b->allPoints);
printf("Grains allFaces = %lu\n", b->allFaces);
for (int i = 0; i < (b->allPoints); i ++) {
# if dimension == 2
printf ("Grains Coordinates of point %d is (%f,%f)\n", i, b->PointsCoord[i][0],b->PointsCoord[i][1]);
# elif dimension ==3
printf ("Grains Coordinates of point %d is (%f,%f,%f)\n", i, b->PointsCoord[i][0], b->PointsCoord[i][1], b->PointsCoord[i][2]);
# endif
}
/* Velocities */
printf("Grains U.x = %f\n", b->vitesseTX);
printf("Grains U.y = %f\n", b->vitesseTY);
printf("Grains U.z = %f\n", b->vitesseTZ);
printf("Grains w.x = %f\n", b->vitesseRX);
printf("Grains w.y = %f\n", b->vitesseRY);
printf("Grains w.z = %f\n", b->vitesseRZ);
/* Physical properties */
printf ("Grains mass = %f\n", b->masse);
printf ("Grains rho_s = %f\n", b->masseVol);
for (int k = 0; k < 6; k++) {
printf ("Grains inertia = %g\n", b->inertie[k]);
}
}
void UpdateParticlesBasilisk (struct BasiliskDataStructure * b, particle * p, const int m) {
int r;
for (int k = 0; k < m; k++) {
#if _MPI
MPI_Bcast (&(b[k].particlessize), 1, MPI_INT, 0, MPI_COMM_WORLD);
/* Geometric aspect */
MPI_Bcast (&(b[k].centreX), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast (&(b[k].centreY), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast (&(b[k].centreZ), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast (&(b[k].rayon), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast (&(b[k].ncorners), 1, MPI_INT, 0, MPI_COMM_WORLD);
/* Coordinates of the corners, number of faces and their numerotation */
MPI_Bcast (&(b[k].allPoints), 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Bcast (&(b[k].allFaces), 1, MPI_INT, 0, MPI_COMM_WORLD);
r = b[k].allPoints;
/* Allocate the structure in the other threads (Grains is serial) */
if (pid() > 0) {
b[k].PointsCoord = (double **) malloc (r*sizeof(double *));
#if dimension == 2
for (int i = 0; i < r; i++)
b[k].PointsCoord[i] = (double *) malloc (2 * sizeof(double));
#elif dimension == 3
for (int i = 0; i < r; i++)
b[k].PointsCoord[i] = (double *) malloc (3 * sizeof(double));
#endif
}
for (int i = 0; i < r; i++) {
MPI_Bcast (&(b[k].PointsCoord[i][0]), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast (&(b[k].PointsCoord[i][1]), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
#if dimension ==3
MPI_Bcast (&(b[k].PointsCoord[i][2]), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
#endif
}
/* Allocate the structures in the other threads (Grains is serial) */
r = b[k].allFaces;
if (pid() > 0) {
b[k].FacesIndex = (long int **) malloc (r * sizeof(long int *));
b[k].numPointsOnFaces = (long int *) malloc (r * sizeof(long int));
}
/* Broadcast-allocate-broadcast */
for (int i = 0; i < r; i++) {
MPI_Bcast (&(b[k].numPointsOnFaces[i]), 1, MPI_LONG, 0, MPI_COMM_WORLD);
int rr = b[k].numPointsOnFaces[i];
if (pid() > 0)
b[k].FacesIndex[i] = (long int *)malloc(rr * sizeof(long int));
MPI_Barrier (MPI_COMM_WORLD);
for (int j = 0; j < rr; j++)
MPI_Bcast (&(b[k].FacesIndex[i][j]), 1, MPI_LONG, 0, MPI_COMM_WORLD);
}
/* Velocities */
MPI_Bcast (&(b[k].vitesseTX), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast (&(b[k].vitesseTY), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast (&(b[k].vitesseTZ), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast (&(b[k].vitesseRX), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast (&(b[k].vitesseRY), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast (&(b[k].vitesseRZ), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
/* Physical properties */
MPI_Bcast (&(b[k].masseVol), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast (&(b[k].masse), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
for (int i = 0; i < 6; i++)
MPI_Bcast (&(b[k].inertie[i]), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
#endif
/* Geometric aspect */
GeomParameter gg;
coord cc = {0., 0., 0.};
cc.x = b[k].centreX;
cc.y = b[k].centreY;
cc.z = b[k].centreZ;
gg.center = cc;
gg.radius = b[k].rayon;
gg.ncorners = b[k].ncorners;
gg.allPoints = b[k].allPoints;
gg.allFaces = b[k].allFaces;
/* We need the coordinates of the corners for a cube */
/* Allocate on all threads */
r = gg.ncorners;
gg.cornersCoord = (double **) malloc (r*sizeof(double *));
#if dimension == 2
for (int i = 0; i < r; i++)
gg.cornersCoord[i] = (double *) malloc (2 * sizeof(double));
#elif dimension == 3
for (int i = 0; i < r; i++)
gg.cornersCoord[i] = (double *) malloc (3 * sizeof(double));
#endif
int ndim = 0;
for (int i = 0; i < r; i++) {
#if dimension == 2
ndim = 2;
#elif dimension == 3
ndim = 3;
#endif
for (int j = 0; j < ndim; j++) {
gg.cornersCoord[i][j] = b[k].PointsCoord[i][j];
}
}
/* We need the indices of the corners for a cube */
/* Allocate */
r = gg.allFaces;
gg.cornersIndex = (long int **) malloc (r * sizeof(long int *));
gg.numPointsOnFaces = (long int *) malloc (r * sizeof(long int));
for (int i = 0; i < r; i++) {
int rr = b[k].numPointsOnFaces[i];
gg.cornersIndex[i] = (long int *) malloc (rr * sizeof(long int));
gg.numPointsOnFaces[i] = rr;
for (int j = 0; j < rr; j++) {
gg.cornersIndex[i][j] = b[k].FacesIndex[i][j];
}
}
p[k].g = gg;
#if DLM_Moving_particle
/* Velocities */
cc.x = b[k].vitesseTX;
cc.y = b[k].vitesseTY;
cc.z = b[k].vitesseTZ;
#if TRANSLATION
/* Save the previous translational velocity before updating */
p[k].Unm1 = p[k].U;
p[k].U = cc;
#endif
cc.x = b[k].vitesseRX;
cc.y = b[k].vitesseRY;
cc.z = b[k].vitesseRZ;
#if ROTATION
/* Save the previous angular velocity before updating */
p[k].wnm1 = p[k].w;
p[k].w = cc;
#endif
#endif
/* Physical properties */
p[k].M = b[k].masse;
p[k].rho_s = b[k].masseVol;
p[k].Vp = (p[k].M)/(p[k].rho_s);
/* Inertia tensor: Ggrains store them as */
/* inertie[0] = Ixx; */
/* inertie[1] = Ixy; */
/* inertie[2] = Ixz; */
/* inertie[3] = Iyy; */
/* inertie[4] = Iyz; */
/* inertie[5] = Izz; */
/* Basilisk stores these as */
/* Ip[0] = Ixx */
/* Ip[1] = Iyy */
/* Ip[2] = Izz */
/* Ip[3] = Ixy */
/* Ip[4] = Ixz */
/* Ip[5] = Iyz */
/* Ixx */
p[k].Ip[0] = b[k].inertie[0];
/* Iyy */
p[k].Ip[1] = b[k].inertie[3];
/* Izz */
p[k].Ip[2] = b[k].inertie[5];
/* Ixy */
p[k].Ip[3] = b[k].inertie[1];
/* Ixz */
p[k].Ip[4] = b[k].inertie[2];
/* Iyz */
p[k].Ip[5] = b[k].inertie[4];
}
}
void UpdateBasiliskStructure (struct BasiliskDataStructure * b, particle * p, const int m) {
coord U = {0., 0., 0.};
coord w = {0., 0., 0.};
for (int k = 0; k < m; k++) {
#if DLM_Moving_particle
#if TRANSLATION
U = p[k].U;
#endif
#if ROTATION
w = p[k].w;
#endif
#endif
b[k].vitesseTX = U.x;
b[k].vitesseTY = U.y;
b[k].vitesseTZ = U.z;
b[k].vitesseRX = w.x;
b[k].vitesseRY = w.y;
b[k].vitesseRZ = w.z;
}
}
void unallocateBasiliskDataStructure(struct BasiliskDataStructure * b, const int m) {
double * c;
for (int k = 0; k < m; k++) {
/* Unallocate fields of BasiliskDataStructure that are no longer needed (we transferred everything on the particles structure */
/* printf("addresse of b[k].PointsCoord = %p\n", (void *) b[k].PointsCoord); */
for (signed long int i = 0; i < b[k].allPoints; i++) {
c = b[k].PointsCoord[i];
/* printf("addresse of b[k].PointsCoord[i] = %p\n",(void *) b[k].PointsCoord[i]); */
/* printf("c = %p\n",(void *) c); */
if(c)
free(c);
}
/* printf("ok free(b[k].PointsCoord[i]) on thread %d\n",pid()); */
free(b[k].PointsCoord);
/* printf("ok free(b[k].PointsCoord) on thread %d\n",pid()); */
for (signed long int i = 0; i < b[k].allFaces; i++) {
free(b[k].FacesIndex[i]);
}
/* printf("ok free(b[k].FaceIndex[i]) on thread %d\n",pid()); */
free(b[k].FacesIndex);
free(b[k].numPointsOnFaces);
}
}
