/** # Helper functions for the fictitious-domain method with distributed Lagrangian multipliers. */ /** Structure for the coordinates of the fictitious domain's boundary. */ typedef struct { double *x; double *y; double *z; int m; int nm; } SolidBodyBoundary; /** Geometric parameters of the particle. */ typedef struct { coord center; double radius; int ncorners, allPoints, allFaces; double ** cornersCoord; long int ** cornersIndex; long int * numPointsOnFaces; /*3 principals vectors of the cube */ coord u1, v1, w1; coord mins, maxs; } GeomParameter; typedef struct { SolidBodyBoundary s; GeomParameter g; GeomParameter gnm1; #if DLM_Moving_particle double kn, en, vzero, wished_ratio; coord normalvector; coord gravity; double M, Ip[6], rho_s, Vp, DLMFD_couplingfactor; #if TRANSLATION coord U, Unm1, qU, tU; #endif #if ROTATION coord w, wnm1, qw, tw; #endif #endif size_t pnum, iswall, iscube; coord wallmax, wallmin, wallpos; coord imposedU, imposedw; Cache Interior; Cache reduced_domain; coord adforce; long tcells, tmultipliers; } particle; /** Function that allocates in memory the SolidBodyBoundary structure. */ void allocate_SolidBodyBoundary (SolidBodyBoundary *sbm, const int m) { sbm->x = (double*) calloc( m, sizeof(double)); sbm->y = (double*) calloc( m, sizeof(double)); if ( dimension == 3 ) sbm->z = (double*) calloc( m, sizeof(double)); else sbm->z = NULL; sbm->m = m; } void reallocate_SolidBodyBoundary (SolidBodyBoundary *sbm, const int m) { sbm->x = (double*) realloc (sbm->x, m*sizeof(double)); sbm->y = (double*) realloc (sbm->y, m*sizeof(double)); if (dimension == 3) sbm->z = (double*) realloc (sbm->z, m*sizeof(double)); sbm->m = m; } /** Function that de-allocates in memory the SolidBodyBoundary structure. */ void free_SolidBodyBoundary (SolidBodyBoundary *sbm) { free(sbm->x); sbm->x = NULL; free(sbm->y); sbm->y = NULL; if ( dimension == 3 ){free(sbm->z); sbm->z = NULL;} } /** Function that allocates an initial Cache structure */ void allocate_Cache (Cache * p) { p->n = 0; p->nm = 0; if (p->n >= p->nm) { p->p = (Index *) calloc (5, sizeof(int)); } } /** Set of functions for the sphere/circle as fictitious domain */ void compute_nboundary (const GeomParameter gcp, int * nb) { coord pos[6]; Cache poscache = {0}; Point lpoint; pos[0].x = gcp.center.x + gcp.radius + X0; pos[0].y = gcp.center.y + Y0; pos[1].x = gcp.center.x - gcp.radius + X0; pos[1].y = gcp.center.y + Y0; pos[2].x = gcp.center.x + X0; pos[2].y = gcp.center.y + gcp.radius + Y0; pos[3].x = gcp.center.x + X0; pos[3].y = gcp.center.y - gcp.radius + Y0; #if dimension == 3 pos[0].z = gcp.center.z + Z0; pos[1].z = gcp.center.z + Z0; pos[2].z = gcp.center.z + Z0; pos[3].z = gcp.center.z + Z0; pos[4].x = gcp.center.x + X0; pos[4].y = gcp.center.y + Y0; pos[4].z = gcp.center.z + gcp.radius + Z0; pos[5].x = gcp.center.x + X0; pos[5].y = gcp.center.y + Y0; pos[5].z = gcp.center.z - gcp.radius + Z0; #endif *nb = 0; lpoint.level = -1; int ip = 0; /* We assume that around the sphere at least one point has to be at maximum level of refinement. */ for (ip = 0; ip < dimension*2; ip++) { #if dimension == 3 lpoint = locate (pos[ip].x, pos[ip].y, pos[ip].z); #elif dimension ==2 lpoint = locate (pos[ip].x, pos[ip].y); #endif if (lpoint.level == depth()) break; } /** Multiple threads can have a different point at the same level of refinement. This is fine, they will do the same computation for *nb. This problem does not exist in serial but the code below still works. */ if (lpoint.level == depth()) { /** Only this thread creates the Cache ... */ cache_append(&poscache, lpoint, 0); /** and only this thread computes the number of boundary points ... */ foreach_cache(poscache) { #if dimension == 2 *nb += (int){floor(gcp.radius*2*pi/(sqrt(2)*Delta))}; /* Debug periodic boundaries, we use only one multiplier */ /* *nb =1.; */ #elif dimension == 3 *nb += (int){floor(pow( 3.809 * gcp.radius / (sqrt(3)*Delta), 2.))}; #endif } /** and finally, this thread destroys the cache. */ free(poscache.p); } #if _MPI MPI_Barrier(MPI_COMM_WORLD); mpi_all_reduce (*nb, MPI_INT, MPI_MAX); #endif if(*nb == 0) printf("nboundary = 0: No boundary points !!!\n"); } void create_FD_Interior_Particle (particle * p, vector Index_lambda, vector shift, scalar flag) { GeomParameter gci = p->g; Cache * fd = &(p->Interior); /** Create the cache for the interior points */ foreach() { #if dimension == 2 if ((sq(x - (gci.center.x + X0)) + sq(y - (gci.center.y + Y0))) < sq(gci.radius)) { cache_append (fd, point, 0); /* tagg cell with the number of the particle */ if ((int)Index_lambda.y[] == -1) Index_lambda.y[] = p->pnum; } if (Period.x) { if ((sq(x - L0 - (gci.center.x + X0)) + sq(y - (gci.center.y + Y0))) < sq(gci.radius)) { cache_append (fd, point, 0); /* tagg cell with the number of the particle */ if ((int)Index_lambda.y[] == -1) Index_lambda.y[] = p->pnum; if ((int)Index_lambda.y[] == p->pnum) shift.x[] = L0; } if ((sq(x + L0 - (gci.center.x + X0)) + sq(y - (gci.center.y + Y0))) < sq(gci.radius)) { cache_append (fd, point, 0); /* tagg cell with the number of the particle */ if ((int)Index_lambda.y[] == -1) Index_lambda.y[] = p->pnum; if ((int)Index_lambda.y[] == p->pnum) shift.x[] = -L0; } } if (Period.y) { if ((sq(x - (gci.center.x + X0)) + sq(y - L0 - (gci.center.y + Y0))) < sq(gci.radius)) { cache_append (fd, point, 0); /* tagg cell with the number of the particle */ if ((int)Index_lambda.y[] == -1) { Index_lambda.y[] = p->pnum; if (flag[] < 1) shift.y[] = L0; } } if ((sq(x - (gci.center.x + X0)) + sq(y + L0 - (gci.center.y + Y0))) < sq(gci.radius)) { cache_append (fd, point, 0); /* tagg cell with the number of the particle */ if ((int)Index_lambda.y[] == -1) { Index_lambda.y[] = p->pnum; if (flag[] < 1) shift.y[] = -L0; } } } #endif #if dimension == 3 if ((sq(x - (gci.center.x + X0)) + sq(y - (gci.center.y + Y0)) + sq(z - (gci.center.z + Z0))) < sq(gci.radius)) { cache_append (fd, point, 0); /* tagg cell with the number of the particle */ if ((int)Index_lambda.y[] == -1) Index_lambda.y[] = p->pnum; } if (Period.x) { if ((sq(x - L0 - (gci.center.x + X0) ) + sq(y - (gci.center.y + Y0)) + sq(z - (gci.center.z + Z0))) < sq(gci.radius)) { cache_append (fd, point, 0); /* tagg cell with the number of the particle */ if ((int)Index_lambda.y[] == -1) Index_lambda.y[] = p->pnum; if ((int)Index_lambda.y[] == p->pnum) shift.x[] = L0; } if ((sq(x + L0 - (gci.center.x + X0) ) + sq(y - (gci.center.y + Y0)) + sq(z - (gci.center.z + Z0))) < sq(gci.radius)) { cache_append (fd, point, 0); /* tagg cell with the number of the particle */ if ((int)Index_lambda.y[] == -1) Index_lambda.y[] = p->pnum; if ((int)Index_lambda.y[] == p->pnum) shift.x[] = -L0; } } if (Period.y) { if ((sq(x - (gci.center.x + X0)) + sq(y - L0 - (gci.center.y + Y0)) + sq(z - (gci.center.z + Z0))) < sq(gci.radius)) { cache_append (fd, point, 0); /* tagg cell with the number of the particle */ if ((int)Index_lambda.y[] == -1) { Index_lambda.y[] = p->pnum; if (flag[] < 1) shift.y[] = L0; } } if ((sq(x - (gci.center.x + X0)) + sq(y + L0 - (gci.center.y + Y0)) + sq(z - (gci.center.z + Z0))) < sq(gci.radius)) { cache_append (fd, point, 0); /* tagg cell with the number of the particle */ if ((int)Index_lambda.y[] == -1) { Index_lambda.y[] = p->pnum; if (flag[] < 1) shift.y[] = -L0; } } } if (Period.z) { if ((sq(x - (gci.center.x + X0)) + sq(y - (gci.center.y + Y0)) + sq(z - L0 - (gci.center.z + Z0))) < sq(gci.radius)) { cache_append (fd, point, 0); /* tagg cell with the number of the particle */ if ((int)Index_lambda.y[] == -1) { Index_lambda.y[] = p->pnum; if (flag[] < 1) shift.z[] = L0; } } if ((sq(x - (gci.center.x + X0)) + sq(y - (gci.center.y + Y0)) + sq(z + L0 - (gci.center.z + Z0))) < sq(gci.radius)) { cache_append (fd, point, 0); /* tagg cell with the number of the particle */ if ((int)Index_lambda.y[] == -1) { Index_lambda.y[] = p->pnum; if (flag[] < 1) shift.z[] = -L0; } } } #endif } cache_shrink (fd); } void create_FD_Boundary_Particle (GeomParameter gcp, SolidBodyBoundary * dlm_bd, const int nsphere, vector pshift) { int m = nsphere; coord pos; Point lpoint; coord shift = {0., 0., 0.}; coord ori = {X0, Y0, Z0}; #if dimension == 2 int i; double theta[m]; double radius = gcp.radius; coord fact_theta[m]; for (i = 0; i < m; i++) { theta[i] = i*2*pi/m; fact_theta[i].x = cos (theta[i]); fact_theta[i].y = sin (theta[i]); /* Assign positions x,y on a circle's boundary */ pos.x = radius*fact_theta[i].x + gcp.center.x + X0; pos.y = radius*fact_theta[i].y + gcp.center.y + Y0; /* Check if the point falls outside of the domain */ foreach_dimension() { if (Period.x) { shift.x = 0.; if (pos.x > (L0 + ori.x)) { pos.x -= L0; shift.x = -L0; } if (pos.x < (0. + ori.x)) { pos.x += L0; shift.x = L0; } } dlm_bd->x[i] = pos.x; } Cache poscache = {0}; lpoint = locate (pos.x, pos.y); if (lpoint.level > -1) { cache_append (&poscache, lpoint, 0); foreach_cache (poscache) { foreach_dimension() { pshift.x[] = shift.x; } } free (poscache.p); } } #elif dimension == 3 /* Lay out points on a sphere with a Z oriented spiral scheme */ /* This portion of code is recovered from Peligriff and adapted */ /* More information of this method can be found at: */ /* Saff, E. & Kuijlaars, A. Distributing many points on a sphere The Mathematical Intelligencer, 1997 */ double spiral_spacing_correction = 0; foreach_level(depth()) spiral_spacing_correction = 2.*Delta/sqrt(3.); double hydro_radius = gcp.radius; size_t k; /* Number of points */ size_t NSpiral = m; /* Spiral points construction */ double hk, thetak, phik, phikm1 = 0., TwoPi = 2. * pi, Cspiral = 3.6 / sqrt( NSpiral ), dphi = 0.; for (k = 0; k < NSpiral; ++k) { hk = - 1. + 2. * (double)(k) / ( NSpiral - 1. ); thetak = acos( hk ) ; if ( k == 0 ) { phik = 0.; thetak = pi - 0.5 * spiral_spacing_correction / hydro_radius ; } else if ( k == NSpiral - 1 ) { phik = phikm1 + 1. * dphi; if ( phik > TwoPi ) phik -= TwoPi ; thetak = 0.5 * spiral_spacing_correction / hydro_radius ; } else { dphi = Cspiral / sqrt( 1. - hk * hk ) ; phik = phikm1 + dphi ; if ( phik > TwoPi ) phik -= TwoPi ; } phikm1 = phik ; if ( k == 1 ) thetak -= 0.4 * spiral_spacing_correction / hydro_radius ; if ( k == NSpiral - 2 ) { phik -= 0.1 * dphi ; thetak += 0.25 * spiral_spacing_correction / hydro_radius ; } pos.x = hydro_radius * cos( phik ) * sin( thetak ) + gcp.center.x + X0; pos.y = hydro_radius * sin( phik ) * sin( thetak ) + gcp.center.y + Y0; pos.z = hydro_radius * cos( thetak ) + gcp.center.z + Z0; foreach_dimension() { if (Period.x) { shift.x = 0.; if (pos.x > (L0 + ori.x)) { pos.x -= L0; shift.x = -L0; } if (pos.x < (0. + ori.x)) { pos.x += L0; shift.x = L0; } } } Cache poscache = {0}; /* find the cell of this multiplier */ lpoint = locate (pos.x, pos.y, pos.z); if (lpoint.level > -1) { cache_append(&poscache, lpoint, 0); foreach_cache(poscache) { foreach_dimension() { pshift.x[] = shift.x; } } free(poscache.p); } dlm_bd->x[k] = pos.x; dlm_bd->y[k] = pos.y; dlm_bd->z[k] = pos.z; } #endif boundary((scalar*){pshift}); } /** Set of functions for the walls as fictitious domain */ void create_FD_Interior_Wall (particle *p, vector index, scalar flag) { coord wallmin, wallmax; Cache * c; /** Create the cache of the interior points for a wall*/ c = &(p->Interior); wallmin = p->wallmin; wallmax = p->wallmax; foreach() { #if dimension == 2 if ((x >= wallmin.x) && (y >= wallmin.y) && (x <= wallmax.x) && (y <= wallmax.y)) { index.y[] = p->pnum; cache_append (c, point, 0); } #elif dimension == 3 if ((x > wallmin.x) && (y > wallmin.y) && (z > wallmin.z) && (x < wallmax.x) && (y < wallmax.y) && (z < wallmax.z)) { if ((int)index.y[] == -1 && flag[] < 1) index.y[] = p->pnum; cache_append (c, point, 0); } #endif } cache_shrink (c); } void create_FD_Boundary_Wall (particle *p) { SolidBodyBoundary * sbm; int m = 0; coord wallpos; double h = 0; foreach_level(depth()) h = Delta; sbm = &(p->s); m = sbm->m; wallpos = p->wallpos; for (int i = 0; i < m; i++) { sbm->x[i] = (3.*h/4.) + i*h; sbm->y[i] = wallpos.y; } for (int i = m/2; i < m; i++) { sbm->x[i] -= h/2.; } } void compute_nboundary_Wall (coord wallpos, int * nb) { double mindelta = 0.; foreach_level(depth()) mindelta = Delta; *nb = floor(wallpos.x/mindelta); } /** Set of functions for the cube as fictitious domain */ void compute_nboundary_Cube_v2 (GeomParameter * gcp, int * nb, int * lN) { Cache poscache = {0}; Point lpoint; *nb = 0; coord pos = {0., 0., 0.}; int ip = 0; while (*nb == 0) { pos.x = gcp->cornersCoord[ip][0]; pos.y = gcp->cornersCoord[ip][1]; #if dimension == 3 pos.z = gcp->cornersCoord[ip][2]; lpoint = locate(pos.x, pos.y, pos.z); #elif dimension ==2 lpoint = locate(pos.x, pos.y); #endif /** Only one thread has the point in its domain (works in serial too). */ if (lpoint.level > -1) { /** Only this thread creates the Cache ... */ cache_append(&poscache, lpoint, 0); /** and only this thread computes the number of boundary points ... */ /* Grains sends an equivalent radius that is multiplied by sqrt(3)/2 */ double lengthedge = 2.*(gcp->radius)/sqrt(3.); /* printf ("lengthedge = %f on thread %d\n", lengthedge, pid() ); */ foreach_cache (poscache) { *lN = floor (sqrt(3.)*lengthedge/(2.*Delta)); } #if dimension == 2 /* number of points required for the 4 edges of the square */ *nb += (*lN-2)*4; /* number of points required for the face */ *nb += (*lN-2)*(*lN-2); /* number of points required for the 4 corners */ *nb += 4; #elif dimension == 3 /* number of points required for the 12 edges of the cube */ *nb += (*lN-2)*12; /* number of points required for the 6 faces of the cube */ *nb += 6*(*lN-2)*(*lN-2); /* number of points required for the 8 corners */ *nb += 8; #endif /** and finally, this thread destroys the cache. */ free(poscache.p); } #if _MPI MPI_Barrier(MPI_COMM_WORLD); mpi_all_reduce(*nb, MPI_INT, MPI_MAX); mpi_all_reduce(*lN, MPI_INT, MPI_MAX); #endif if (ip < gcp->ncorners) ip++; else break; } /* printf("lN number of points for a complete edge:%d\n",*lN); */ if (*nb == 0) fprintf(stderr,"nboundary = 0: No boundary points for the cube/square !!!\n"); } void distribute_points_edge_Cube_v2 (coord const corner1, coord const corner2, SolidBodyBoundary * dlm_bd, int const lN, int const istart) { # if dimension == 3 if (lN > 0) { coord dinc; foreach_dimension() { dinc.x = (corner2.x - corner1.x)/(lN-1); } /* printf( "length edge corner to corner = %f\n", sqrt( sq(corner2.x - corner1.x) + sq(corner2.y - corner1.y)+ sq(corner2.z - corner1.z))); */ for (int i = 1; i <= lN-2; i++) { dlm_bd->x[istart + i -1] = corner1.x + (double)i * dinc.x; dlm_bd->y[istart + i -1] = corner1.y + (double)i * dinc.y; dlm_bd->z[istart + i -1] = corner1.z + (double)i * dinc.z; } } #endif } bool is_it_in_cube_v2 (coord * u, coord * v, coord * w, coord * mins, coord * maxs, coord * checkpt) { /* a point x with coord (x,y,z) lies in the rectangle if the 3 following conditions are satisfied */ /* 1- u.p0 <= u.x <= u.p3 */ /* 2- v.p0 <= v.x <= v.p1 */ /* 3- w.p0 <= w.w <= w.p4 */ double x = checkpt->x; double y = checkpt->y; double z = checkpt->z; double checkval = 0.; coord u1 = *u; coord v1 = *v; coord w1 = *w; coord cpt = {0., 0., 0.}; bool isin = false; checkval = u1.x*x + u1.y*y + u1.z*z; cpt.x = checkval; checkval = v1.x*x + v1.y*y + v1.z*z; cpt.y = checkval; checkval = w1.x*x + w1.y*y + w1.z*z; cpt.z = checkval; if ((cpt.x >= maxs->x) && (cpt.x <= mins->x) && (cpt.y >= maxs->y) && (cpt.y <= mins->y) && (cpt.z >= maxs->z) && (cpt.z <= mins->z) ) { isin = true; } return isin; } void compute_principal_vectors_Cubes (particle * p) { GeomParameter * gcp = &(p->g); int nfaces = gcp->allFaces; int npoints; /* get the 3 directions u1, u2, u3 of the cube with such that */ /* u1 = corner0 - corner3; */ /* u2 = corner0 - corner1; */ /* u3 = corner0 - corner4; */ coord corner0 = {0, 0, 0}, corner1 = {0, 0, 0}, corner3 = {0, 0, 0}, corner4 = {0, 0, 0}; coord u1, v1, w1; coord mins, maxs; /* find the coordinates of these 4 corners */ for (int i = 0; i < nfaces; i++) { npoints = gcp->numPointsOnFaces[i]; for (int j = 0; j < npoints; j++) { /* printf("index = %lu\n", gcp->cornersIndex[i][j]); */ if (gcp->cornersIndex[i][j] == 0) { long int ii = gcp->cornersIndex[i][j]; corner0.x = gcp->cornersCoord[ii][0]; corner0.y = gcp->cornersCoord[ii][1]; corner0.z = gcp->cornersCoord[ii][2]; /* printf("index = %lu with coordinate (%f,%f,%f)\n", gcp->cornersIndex[i][j], corner0.x, corner0.y, corner0.z); */ } if (gcp->cornersIndex[i][j] == 1) { long int ii = gcp->cornersIndex[i][j]; corner1.x = gcp->cornersCoord[ii][0]; corner1.y = gcp->cornersCoord[ii][1]; corner1.z = gcp->cornersCoord[ii][2]; /* printf("index = %lu with coordinate (%f,%f,%f)\n", gcp->cornersIndex[i][j], corner1.x, corner1.y, corner1.z); */ } if (gcp->cornersIndex[i][j] == 3) { long int ii = gcp->cornersIndex[i][j]; corner3.x = gcp->cornersCoord[ii][0]; corner3.y = gcp->cornersCoord[ii][1]; corner3.z = gcp->cornersCoord[ii][2]; /* printf("index = %lu with coordinate (%f,%f,%f)\n", gcp->cornersIndex[i][j], corner3.x, corner3.y, corner3.z); */ } if (gcp->cornersIndex[i][j] == 4) { long int ii = gcp->cornersIndex[i][j]; corner4.x = gcp->cornersCoord[ii][0]; corner4.y = gcp->cornersCoord[ii][1]; corner4.z = gcp->cornersCoord[ii][2]; /* printf("index = %lu with coordinate (%f,%f,%f)\n", gcp->cornersIndex[i][j], corner4.x, corner4.y, corner4.z); */ } } } foreach_dimension() { u1.x = corner0.x - corner3.x; v1.x = corner0.x - corner1.x; w1.x = corner0.x - corner4.x; mins.x = 0.; maxs.x = 0.; } gcp->u1 = u1; gcp->v1 = v1; gcp->w1 = w1; double minval = 0., maxval = 0.; foreach_dimension() { minval += u1.x*corner0.x; maxval += u1.x*corner3.x; } mins.x = (minval); maxs.x = (maxval); minval = 0; maxval = 0; foreach_dimension() { minval += v1.x*corner0.x; maxval += v1.x*corner1.x; } mins.y = (minval); maxs.y = (maxval); minval = 0; maxval = 0; foreach_dimension() { minval += w1.x*corner0.x; maxval += w1.x*corner4.x; } mins.z = (minval); maxs.z = (maxval); gcp->mins = mins; gcp->maxs = maxs; } /* void compute_principal_mins_maxs_Cubes (particle * p) { */ /* GeomParameter * gcp = &(p->g); */ /* coord u1 = gcp->u1; */ /* coord v1 = gcp->v1; */ /* coord w1 = gcp->w1; */ /* } */ void create_FD_Boundary_Cube_v2 (GeomParameter * gcp, SolidBodyBoundary * dlm_bd, const int m, const int lN, vector pshift) { int nfaces = gcp->allFaces; int iref, i1, i2, ichoice; ichoice = 0; int isb = 0; int npoints; /* Add first interrior points on surfaces */ for (int i = 0; i < nfaces; i++) { npoints = gcp->numPointsOnFaces[i]; iref = gcp->cornersIndex[i][ichoice]; i1 = gcp->cornersIndex[i][ichoice + 1]; i2 = gcp->cornersIndex[i][npoints-1]; /* printf("on face %d, iref = %d, i1 = %d, i2 = %d\n", i, iref, i1, i2); */ coord refcorner = {gcp->cornersCoord[iref][0], gcp->cornersCoord[iref][1],gcp->cornersCoord[iref][2]} ; coord dir1 = {gcp->cornersCoord[i1][0], gcp->cornersCoord[i1][1],gcp->cornersCoord[i1][2]}; coord dir2 = {gcp->cornersCoord[i2][0], gcp->cornersCoord[i2][1],gcp->cornersCoord[i2][2]}; /* printf ("corresponding coord for ref = (%f,%f,%f)\n",refcorner.x, refcorner.y, refcorner.z); */ /* printf ("corresponding coord for i1 = (%f,%f,%f)\n", dir1.x, dir1.y, dir1.z); */ /* printf ("corresponding coord for i2 = (%f,%f,%f)\n", dir2.x, dir2.y, dir2.z); */ foreach_dimension() { dir1.x -= refcorner.x; dir2.x -= refcorner.x; dir1.x /= (lN-1); dir2.x /= (lN-1); } for (int ii = 1; ii <= lN-2; ii++) { for (int jj = 1; jj <= lN-2; jj++) { dlm_bd->x[isb] = refcorner.x + (double) ii * dir1.x + (double) jj * dir2.x; dlm_bd->y[isb] = refcorner.y + (double) ii * dir1.y + (double) jj * dir2.y; dlm_bd->z[isb] = refcorner.z + (double) ii * dir1.z + (double) jj * dir2.z; isb++; } } } int allindextable[8][8] = {{0}}; int j1,jm1; /* Add points on the edges without the corners*/ for (int i = 0; i < nfaces; i++) { npoints = gcp->numPointsOnFaces[i]; i1 = gcp->cornersIndex[i][1]; for (int j = 1; j < npoints; j++) { jm1 = gcp->cornersIndex[i][j-1]; j1 = gcp->cornersIndex[i][j]; if (jm1 > j1) { if (allindextable[jm1][j1] == 0) { coord c1 = {gcp->cornersCoord[jm1][0], gcp->cornersCoord[jm1][1], gcp->cornersCoord[jm1][2]}; coord c2 = {gcp->cornersCoord[j1][0], gcp->cornersCoord[j1][1], gcp->cornersCoord[j1][2]}; distribute_points_edge_Cube_v2 (c1, c2, dlm_bd, lN, isb); allindextable[jm1][j1] = 1; isb +=lN-2; } } else { if (allindextable[j1][jm1] == 0) { coord c1 = {gcp->cornersCoord[j1][0], gcp->cornersCoord[j1][1], gcp->cornersCoord[j1][2]}; coord c2 = {gcp->cornersCoord[jm1][0], gcp->cornersCoord[jm1][1], gcp->cornersCoord[jm1][2]}; distribute_points_edge_Cube_v2 (c1, c2, dlm_bd, lN, isb); allindextable[j1][jm1] = 1; isb +=lN-2; } } /* printf ("isb after face %d = %d\n", i, isb); */ } } /* printf ("isb after all edges = %d\n", isb); */ /* Add the final 8 corners points */ for (int i = 0; i < gcp->ncorners; i++) { dlm_bd->x[isb] = gcp->cornersCoord[i][0]; dlm_bd->y[isb] = gcp->cornersCoord[i][1]; dlm_bd->z[isb] = gcp->cornersCoord[i][2]; isb++; } /* printf ("isb after 8 corners = %d\n", isb); */ } void create_FD_Interior_Cube_v2 (particle *p, vector Index_lambda, vector pshift) { Cache * c; /** Create the cache of the interior points for a cube*/ c = &(p->Interior); /* compute the 3 principal vector of the cube */ compute_principal_vectors_Cubes (p); /* a point x with coord (x,y,z) lies in the cube if the 3 following conditions are satisfied */ /* 1- u.p0 <= u.x <= u.p3 */ /* 2- v.p0 <= v.x <= v.p1 */ /* 3- w.p0 <= w.w <= w.p4 */ coord checkpt; coord u1 = p->g.u1; coord v1 = p->g.v1; coord w1 = p->g.w1; coord mins = p->g.mins; coord maxs = p->g.maxs; /* Min/Max coordinates for the AABB (Axed-Aligned-Bounding-Box) */ coord mincoord = {HUGE, HUGE, HUGE}; coord maxcoord = {-HUGE, -HUGE, -HUGE}; double ** table = p->g.cornersCoord; for (int ii = 0; ii < p->g.ncorners; ii++) { if (mincoord.x > table[ii][0]) mincoord.x = table[ii][0]; if (mincoord.y > table[ii][1]) mincoord.y = table[ii][1]; if (mincoord.z > table[ii][2]) mincoord.z = table[ii][2]; if (maxcoord.x < table[ii][0]) maxcoord.x = table[ii][0]; if (maxcoord.y < table[ii][1]) maxcoord.y = table[ii][1]; if (maxcoord.z < table[ii][2]) maxcoord.z = table[ii][2]; } /* if (pid() == 0) { */ /* printf ("mincoord = (%f,%f,%f)\n",mincoord.x, mincoord.y, mincoord.z); */ /* printf ("maxcoord = (%f,%f,%f)\n",maxcoord.x, maxcoord.y, maxcoord.z); */ /* } */ foreach() { checkpt.x = x; checkpt.y = y; checkpt.z = z; /* Check only if the point is in the AABB (Axed-Aligned-Bounding-Box) */ if ((x > mincoord.x) && (x < maxcoord.x)) if ((y > mincoord.y) && (y < maxcoord.y)) if ((z > mincoord.z) && (z < maxcoord.z)) /* If yes: check if it is inside the cube now */ if (is_it_in_cube_v2 (&u1, &v1, &w1, &mins, &maxs, &checkpt)) { cache_append (c, point, 0); /* tagg cell with the number of the particle */ if ((int)Index_lambda.y[] == -1) Index_lambda.y[] = p->pnum; } /* /\* For periodic boundaries, we have to check for the image-cubes as well *\/ */ /* if (Period.x) */ /* if (mincoord.x < 0 + X0) */ } /* foreach_dimension() { */ /* if (Period.x) { */ /* } */ /* } */ cache_shrink (c); } void free_particles (particle * pp, const int n) { for (int k = 0; k < n; k++) { SolidBodyBoundary * sbm = &(pp[k].s); free_SolidBodyBoundary(sbm); Cache * c = &(pp[k].Interior); free(c->p); c = &(pp[k].reduced_domain); free(c->p); /* free corner coordinates tables */ double * cc; GeomParameter * gg = &(pp[k].g); for (int i = 0; i < gg->ncorners; i++) { cc = &(gg->cornersCoord[i][0]); if(cc) { free(cc); } } free(gg->cornersCoord); /* free corner's indices tables */ for (int i = 0; i < gg->allFaces; i++) { free(gg->cornersIndex[i]); } free(gg->cornersIndex); free(gg->numPointsOnFaces); } } /** Function that creates the scalar field where the cells contain the indices of the dlmfd points. If there is no dlmfd point it is tagged -1. */ void create_index_lambda_scalar (const SolidBodyBoundary dlm_bd, vector Index_lambda, const int kk) { Point lpoint; int i; Cache *fdlocal; fdlocal = (Cache *){calloc(dlm_bd.m, sizeof(Cache))}; for (i = 0; i < dlm_bd.m; i++) { #if dimension == 2 lpoint = locate(dlm_bd.x[i], dlm_bd.y[i]); #elif dimension == 3 lpoint = locate(dlm_bd.x[i], dlm_bd.y[i], dlm_bd.z[i]); #endif /* */ if ((lpoint.level) == depth()) { /* Create a cache for each fictitious domain's boundary point */ cache_append(&fdlocal[i], lpoint, 0); foreach_cache(fdlocal[i]) { /* Tag cell only if it was not tagged by another particle */ if (Index_lambda.x[] < 0){ Index_lambda.x[] = i; Index_lambda.y[] = kk; if (level != depth()) { printf("On thread %d, point dlmfd %d at (%f, %f, %f) is in a cell that has not the maximum level of refinnement %d, it is on level %d \n", pid(), i, x, y, z, depth(), level); } } } } } for (i = 0; i < dlm_bd.m; i++) { free(fdlocal[i].p); } free (fdlocal); boundary ((scalar*) {Index_lambda}); } #include "weights_functions_backup.h" # ifndef STENCIL_EXTERIOR # define STENCIL_EXTERIOR 0 # endif # ifndef STENCIL_INTERIOR # define STENCIL_INTERIOR 0 # endif /** Function that returns the weight associated to a Lagrange multipliers for a cell within a foreach() loop (ie for a "real" or local cell) associated to a Lagrange multiplier in it's neighborhood. The neighboor can be in another process (i.e be a ghost cell for the current thread). */ double reversed_weight (particle * pp, const coord weightcellpos, const coord lambdacellpos, const coord lambdapos, const double delta, const coord pshift, const double a, const double b) { /* this function has to be embedded within a double foreach_cache() and foreach_neighbor() loop */ coord rel = {0., 0., 0.}; coord relnl = {0., 0., 0.}; int NCX = 0, CX = 0, weight_id = 0; size_t goflag = 0; double weight = 0.; GeomParameter gcb = pp->g; GeomParameter gcbdum = gcb; coord lambdaposdum = lambdapos; /* Shift artificially the center of the particle and position of the Lagrange multipliers for periodic boundary cases */ foreach_dimension() { gcbdum.center.x = gcb.center.x + a*pshift.x; lambdaposdum.x = lambdapos.x + b*pshift.x; } /* compute relative vector X_boundary - X_local = rel from the cell (containning the boundary) position to the boundary's (analytical) position .i.e vector goes from first argument to the second */ compute_relative_vector (lambdacellpos, lambdaposdum, &rel); /* reset dials integers */ NCX = 0; CX = 0; weight_id = 0; goflag = 0; /* assign first normal (gives the quadrant) */ assign_dial (rel, &CX); /* assign fictitious-boundary's normal (use boundary's analytical position) */ assign_dial_fd_boundary (pp, lambdaposdum, gcbdum, delta, &NCX); if (pp->iswall) { #if STENCIL_EXTERIOR /* Stencils oriented toward the wall */ if (lambdapos.y > 0) { if (lambdapos.x < L0/2) NCX = 1; if (lambdapos.x > L0/2) NCX = 2; } if (lambdapos.y < 0) { if (lambdapos.x < L0/2) NCX = 4; if (lambdapos.x > L0/2) NCX = 3; } #endif #if STENCIL_INTERIOR /* Second normal oriented toward the fluid */ if (lambdapos.y > 0) { if (lambdapos.x < L0/2) { NCX = 4; } if (lambdapos.x > L0/2) { NCX = 3; } } if (lambdapos.y < 0) { if (lambdapos.x < L0/2) { NCX = 1; } if (lambdapos.x > L0/2) { NCX = 2; } } #endif } /* compute relative vector from the neighbor cell to the local cell */ compute_relative_vector (weightcellpos, lambdacellpos , &relnl); /* Assign weight ids. */ assign_weight_id_quad_outward (NCX, CX, relnl, delta, &weight_id, &goflag); if (goflag == 1) { weight = compute_weight_Quad (weight_id, lambdaposdum, lambdacellpos, NCX, CX, delta); } return weight; } void remove_too_close_multipliers (particle * p, vector index_lambda) { for (int k = 0; k < NPARTICLES; k++) { SolidBodyBoundary dlm_lambda_to_desactivate; int allocated = 1; allocate_SolidBodyBoundary (&dlm_lambda_to_desactivate, allocated*BSIZE); int countalloc = 0; int other_part; int is_in_other_particle; foreach_level(depth()) { int direct_neigh = 0; is_in_other_particle = 0; other_part = -1; if (((int)index_lambda.x[] > -1) && level == depth() && is_leaf(cell) && (p[k].pnum == (int)index_lambda.y[])) { /* check here if this multiplier is not in another particle's domain, if yes desactive it */ for (int l = 0; l < NPARTICLES; l++) { if (l != p[k].pnum) { particle * other_particle = &(p[l]); /* printf("thread %d, this is particle %zu checkin on particle %zu iteration %d\n", pid(), p[k].pnum, other_particle->pnum, l); */ /* Check particle's type */ if (other_particle->iscube) { compute_principal_vectors_Cubes (other_particle); coord u = other_particle->g.u1; coord v = other_particle->g.v1; coord w = other_particle->g.w1; coord mins = other_particle->g.mins; coord maxs = other_particle->g.maxs; /* current cell's position */ coord checkpt = {x, y, z}; if (is_it_in_cube_v2 (&u, &v, &w, &mins, &maxs, &checkpt)) { is_in_other_particle = 1; other_part = other_particle->pnum; break; } } /* if sphere */ else { GeomParameter gp = other_particle->g; if (is_it_in_sphere (x, y, z, gp)) { is_in_other_particle = 1; other_part = other_particle->pnum; break; } } } } if (is_in_other_particle) { /* printf("thread %d, this is particle %zu which has a cell in particle %d\n", pid(), p[k].pnum, other_part); */ index_lambda.x[] = -1; index_lambda.y[] = other_part; } /* Check if two (or more) Lagrange multipliers (from different particles) are in the same neighborhoods (i.e a 5x5 stencil in 2D), if yes desactivate them. Otherwise the cells of their stencils may be located in each other's domain */ direct_neigh = 0; foreach_neighbor() { if (((int)index_lambda.x[] > -1) && level == depth() && is_leaf(cell) && (p[k].pnum != (int)index_lambda.y[])) { direct_neigh = 1; /* We may catch and identical multiplier more than once here */ /* Add the multiplier(s) indices to a list because MPI imposes that we can't modify the cell within a foreach_neighbor loop */ if (countalloc >= allocated*BSIZE) { allocated ++; reallocate_SolidBodyBoundary (&dlm_lambda_to_desactivate, allocated*BSIZE); } dlm_lambda_to_desactivate.x[countalloc] = x; dlm_lambda_to_desactivate.y[countalloc] = y; #if dimension == 3 dlm_lambda_to_desactivate.z[countalloc] = z; #endif countalloc ++; } } /* Desactivate the local multiplier here */ if (direct_neigh) { index_lambda.x[] = -1; index_lambda.y[] = -1; } } } /* if (countalloc > 0) */ /* printf("there is %d points to desactivate on thread %d\n", countalloc, pid()); */ Cache c = {0}; Point lpoint; #if _MPI int size = npe(); int counts[size]; /* Each thread tells the root how many multipliers it holds */ MPI_Gather(&countalloc, 1, MPI_INT, &counts, 1, MPI_INT, 0, MPI_COMM_WORLD); /* if (pid() == 0) */ /* for (int i = 0; i < size; i++) */ /* printf("particle %d, this is root receiving %d elements by thread %d\n",k, counts[i], i); */ /* Displacements in the receive buffer for MPI_GATHERV */ int disps[size]; /* Displacement for the first chunk of data - 0 */ for (int i = 0; i < size; i++) disps[i] = (i > 0) ? (disps[i-1] + counts[i-1]) : 0; double * alldatax = NULL; double * alldatay = NULL; #if dimension == 3 double * alldataz = NULL; #endif int m = 0; /* Threads 0 allocates and compute the total number of multipliers */ if (pid() == 0) { /* disps[size-1] + counts[size-1] == total number of elements */ /* printf("thread %d : disps[size-1]+counts[size-1] = %d\n", pid(), disps[size-1]+counts[size-1]); */ m = disps[size-1] + counts[size-1]; alldatax = (double*) calloc (m, sizeof(double)); alldatay = (double*) calloc (m, sizeof(double)); #if dimension == 3 alldataz = (double*) calloc (m, sizeof(double)); #endif } /* Gather on thread 0 the coordinates of the multipliers */ MPI_Gatherv (&dlm_lambda_to_desactivate.x[0], countalloc, MPI_DOUBLE, &alldatax[0], counts, disps, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Gatherv (&dlm_lambda_to_desactivate.y[0], countalloc, MPI_DOUBLE, &alldatay[0], counts, disps, MPI_DOUBLE, 0, MPI_COMM_WORLD); #if dimension == 3 MPI_Gatherv (&dlm_lambda_to_desactivate.z[0], countalloc, MPI_DOUBLE, &alldataz[0], counts, disps, MPI_DOUBLE, 0, MPI_COMM_WORLD); #endif /* send the total number of multipliers to all threads */ mpi_all_reduce (m, MPI_INT, MPI_MAX); /* allocate now alldatax, alldatay, alldataz on threads !=0 */ if (pid() != 0) { alldatax = (double*) calloc (m, sizeof(double)); alldatay = (double*) calloc (m, sizeof(double)); #if dimension == 3 alldataz = (double*) calloc (m, sizeof(double)); #endif } /* Now broadcast alldatax,alldatay,alldataz from thread 0 to other threads */ MPI_Bcast (alldatax, m, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast (alldatay, m, MPI_DOUBLE, 0, MPI_COMM_WORLD); #if dimension == 3 MPI_Bcast (alldataz, m, MPI_DOUBLE, 0, MPI_COMM_WORLD); #endif /* for (int i = 0; i < m; i++) { */ /* #if dimension == 2 */ /* printf("thread %d: particle %d, coord of the multiplier %d to be removed (%g,%g,%g)\n", pid(), k, i, alldatax[i], alldatay[i], 0.); */ /* #elif dimension ==3 */ /* printf("thread %d: particle %d, coord of the multiplier %d to be removed (%g,%g,%g)\n", pid(), k, i, alldatax[i], alldatay[i], alldataz[i]); */ /* #endif */ /* } */ for (int i = 0; i < m; i++) { #if dimension == 2 lpoint = locate(alldatax[i], alldatay[i]); #elif dimension == 3 lpoint = locate(alldatax[i], alldatay[i], alldataz[i]); #endif if (lpoint.level > -1) cache_append(&c, lpoint, 0); } free(alldatax); free(alldatay); #if dimension == 3 free(alldataz); #endif #elif _MPI == 0 for (int i = 0; i < countalloc; i++) { #if dimension == 2 lpoint = locate(dlm_lambda_to_desactivate.x[i], dlm_lambda_to_desactivate.y[i]); #elif dimension == 3 lpoint = locate(dlm_lambda_to_desactivate.x[i], dlm_lambda_to_desactivate.y[i], dlm_lambda_to_desactivate.z[i]); #endif if (lpoint.level > -1) cache_append(&c, lpoint, 0); } #endif /* End if _MPI */ /* Finally desactivate the multiplicators */ int counter = 0; foreach_cache(c) { if (index_lambda.x[] > -1 && index_lambda.y[] > -1) { index_lambda.x[] = -1; index_lambda.y[] = -1; counter ++; } } /* printf("thread %d: removed %d multipliers \n", pid(), counter); */ free(c.p); free_SolidBodyBoundary(&dlm_lambda_to_desactivate); } /* End loop NPARTICLES */ boundary((scalar*) {index_lambda}); } /** Function that tags cells associated to a Lagrange multiplier on the boundary of the fictitious domain. */ void reverse_fill_flagfield (particle * p, scalar f, vector index_lambda, const int cacheflag) { for (int k = 0; k < NPARTICLES; k++) { coord rel = {0., 0., 0.}; coord relnl = {0., 0., 0.}; int NCX, CX, weight_id; size_t goflag = 0; coord lambdacellpos = {0., 0., 0.}; coord lambdapos = {0., 0., 0.}; coord localcellpos = {0., 0., 0.}; SolidBodyBoundary dlm_lambda = p[k].s; GeomParameter gcb = p[k].g; Cache * reduced_domain = &(p[k].reduced_domain); GeomParameter gcbdum = gcb; coord lambdaposdum = lambdapos; foreach_level(depth()) { localcellpos.x = x; localcellpos.y = y; #if dimension == 3 localcellpos.z = z; #endif goflag = 0; // Check if there is a Lagrange multiplier in the neigborhood of this cell; foreach_neighbor() { if (((int)index_lambda.x[] > -1) && level == depth() && is_leaf(cell) && (p[k].pnum == (int)index_lambda.y[])) { lambdacellpos.x = x; lambdacellpos.y = y; lambdapos.x = dlm_lambda.x[(int)index_lambda.x[]]; lambdapos.y = dlm_lambda.y[(int)index_lambda.x[]]; #if dimension == 3 lambdacellpos.z = z; lambdapos.z = dlm_lambda.z[(int)index_lambda.x[]]; #endif /* Shift artificially the center of the particle and position of the Lagrange multipliers for periodic boundary cases */ coord ppshift = {0., 0., 0.}; foreach_dimension() { ppshift.x = 0.; if (lambdacellpos.x < 0.) ppshift.x = -L0; if (lambdacellpos.x > L0) ppshift.x = +L0; gcbdum.center.x = gcb.center.x + ppshift.x; lambdaposdum.x = lambdapos.x + ppshift.x; } /* compute relative vector X_boundary - X_local = rel from the cell (containning the boundary) position to the boundary's (analytical) position */ compute_relative_vector (lambdacellpos, lambdaposdum, &rel); /* reset dials integers */ NCX = 0; CX = 0; weight_id = 0; /* assign first normal (gives the quadrant) */ assign_dial (rel, &CX); /* assign fictitious-boundary's normal (use boundary's position) */ assign_dial_fd_boundary (&p[k], lambdaposdum, gcbdum, Delta, &NCX); if (p->iswall) { #if STENCIL_EXTERIOR /* Stencils oriented toward the wall */ if (lambdapos.y > 0) { if (lambdapos.x < L0/2) NCX = 1; if (lambdapos.x > L0/2) NCX = 2; } if (lambdapos.y <0) { if (lambdapos.x < L0/2) NCX = 4; if (lambdapos.x > L0/2) NCX = 3; } #endif #if STENCIL_INTERIOR /* Stencils oriented toward the fluid */ if ((lambdapos.y) > 0) { if (lambdapos.x < L0/2) NCX = 4; if (lambdapos.x > L0/2) NCX = 3; } if (lambdapos.y <0) { if (lambdapos.x < L0/2) { NCX = 1; } if (lambdapos.x > L0/2) NCX = 2; } #endif } /* compute relative vector neighbor to the local cell */ compute_relative_vector (localcellpos, lambdacellpos , &relnl); /* Assign weight ids. */ assign_weight_id_quad_outward (NCX, CX, relnl, Delta, &weight_id, &goflag); } // end if (lambda.x[] > -1) } // end foreach_neigboor loop /* If the cell belongs to the stencil of a Lagrange multiplier tagg it and create a reduced domain to optimize the stencil-traversal cost */ if (goflag == 1) { f[] = 1; if (cacheflag == 1) cache_append (reduced_domain, point, 0); } } boundary ({f, index_lambda.x, index_lambda.y}); }// loop on particles id } /** Function that initialize the particles and the scalar/vector fields needed to the method */ void allocate_and_init_particles (particle * p, const int n, vector e, scalar g, scalar h, vector pshift) { foreach() { e.x[] = -1;// index_lambda.x => lambda indices e.y[] = -1;// particle structure/fd number; e.z[] = -1;// constrained cells (-1:not constrained) g[] = 0; // flagfield h[] = 0; // flagfield_mailleur /* Reset the shift vector for the periodic case */ if (Period.x || Period.y || Period.z) { foreach_dimension() { pshift.x[] = 0.; } } } boundary ({g, h}); boundary ((scalar *){e, pshift}); for (int k = 0; k < n; k++) { Cache * c = NULL; #if debugBD == 0 GeomParameter gci = p[k].g; int m = 0; /* default case, sphere */ if (!(p[k].iswall) && !(p[k].iscube)) { compute_nboundary (gci, &m); allocate_SolidBodyBoundary(&(p[k].s), m); create_FD_Boundary_Particle (gci, &(p[k].s), m, pshift); } if (p[k].iswall) { coord wallpos = p[k].wallmax; compute_nboundary_Wall (wallpos, &m); allocate_SolidBodyBoundary(&(p[k].s), m); create_FD_Boundary_Wall (&p[k]); } if (p[k].iscube) { int lN = 0; compute_nboundary_Cube_v2 (&gci, &m, &lN); allocate_SolidBodyBoundary(&(p[k].s), m); create_FD_Boundary_Cube_v2 (&gci, &(p[k].s), m, lN, pshift); } create_index_lambda_scalar ((p[k].s), e, k); c = &(p[k].reduced_domain); allocate_Cache(c); #endif #if debugInterior == 0 c = &(p[k].Interior); allocate_Cache(c); #endif } } #if DLM_Moving_particle #if TRANSLATION void compute_contact_distance (particle * p, const coord gci, double * delta) { GeomParameter gc = p->g; double radius = gc.radius; *delta = gci.y - radius; } double compute_fwo (const double en, const double v, const double wo) { double k = log(en)/sqrt(sq(pi) + sq(log(en))); double g = atan(-sqrt(1 - sq(k))/k); double val; return val = v*exp(k*g/sqrt(1- sq(k)))*sin(g)/(wo*sqrt(1 - sq(k))); } double derivative_fwo (const double en, const double v, const double wo) { double k = log(en)/sqrt(sq(pi) + sq(log(en))); double g = atan(-sqrt(1 - sq(k))/k); double val; return val = -v*exp(k*g/sqrt(1 - sq(k)))*sin(g)/(sqrt(1 - sq(k))*sq(wo)); } void compute_wo (particle * p) { /* compute an estimate of (kn,gamman) derived from (deltamax,en) */ /* Guess value for wo */ double wo = 100.; double epsilo = 1e-2; int maxiter = 200; double aa = 0.; double bb = 0.; GeomParameter g = p->g; double r = g.radius; double deltamax = (p->wished_ratio)*r; double en = p->en; double v = p->vzero; /* Newton iteration to find the root of delta(w0) */ if (en < 1) { for (int pp = 1; pp <= maxiter; pp++) { aa = compute_fwo (wo, v, en) - deltamax; bb = derivative_fwo (wo, v, en); if (abs(aa-deltamax)/deltamax < epsilo) { break; } wo += -aa/bb; } p->kn = (p->M)*sq(wo); } else { /* if en = 1 the formula simplifies as such */ p->kn = (p->M)*sq(p->vzero)/(sq(deltamax)); } } void compute_Fontact (coord * Fc, particle * p, coord * gci, coord * U, const double gamman) { double delta_colision = 0.; foreach_dimension() { (*Fc).x = 0.; } compute_contact_distance(p, *gci, &delta_colision); if (delta_colision < 0.) { double kn = p->kn; double M = p->M; coord vrel = *U; coord Fel = {0., 0., 0.}; coord Fdm = {0., 0., 0.}; coord normalvec = p->normalvector; foreach_dimension() { /* compute Hookean elastic restoring force */ Fel.x = -kn*delta_colision*normalvec.x; /* compute viscous dynamic force */ Fdm.x = -2.*gamman*M*fabs(vrel.x); /* Fc is the sum of these two forces */ (*Fc).x = Fel.x + Fdm.x; } } } void granular_subproblem (particle * p, const int gravity_flag, const double dt, const double rho_f) { // Mini Granular solver, which solves (1 - rho_f/rho_s)MdU/dt = (1-rho_f/rho_s)Mg + F_c // (1 - rho_f/rho_s) Ip dw/dt = sum_all_particles (r x F_c) - (1 - rho_f/rho_s)w x Ip w // where F_c is the contact force, g the gravity acceleration and M the particle's mass and Ip the inertia tensor double wo = 0. ; double gamman = 0.; double T_c = 0., dtg = 0., M = 0., kn = 0., en = 0.; int miter, gi; coord Uold, Xold; coord k1, k2, k3, k4, Xtemp, Utemp; double fsf; /* particle's structure pointers */ GeomParameter * gci; GeomParameter * gcinm1; coord * U; coord * Unm1; coord * gravity; coord decal = {X0, Y0, Z0}; for (int k = 0; k < NPARTICLES; k++) { fsf = (1. - (rho_f)/(p[k].rho_s)); miter = 0; M = p[k].M; kn = p[k].kn; en = p[k].en; /* compute wo and gamman */ wo = sqrt(kn/M); // wo = sqrt(2kn/M) for sphere-sphere contact and wo = sqrt(kn/M) for sphere-wall contact (Powder tech. Rakotonirina 2018) gamman = -wo*log(en)/sqrt(sq(pi) + sq(log(en))); /* compute the contact time Tc */ T_c = pi/(sqrt(sq(wo) - sq(gamman))); /* set granular timestep */ dtg = T_c/20.; miter = ceil(dt/dtg); dtg = dt/miter; if (k == 0) fprintf (stderr,"Tc = %g, dt = %20.18f, dtg = %20.18f, Tc/dtg = %g, miter = %d, particle = %d\n",T_c, dt, dtg, T_c/dtg, miter, k); /* get particle structure pointers */ /* position */ gcinm1 = &(p[k].gnm1); gci = &(p[k].g); /* translational velocity */ U = &(p[k].U); Unm1 = &(p[k].Unm1); gravity = &(p[k].gravity); /* fprintf (stderr,"gravity = (%f,%f,%f)\n",(*gravity).x,(*gravity).y,(*gravity).z); */ if (gravity_flag) { /* Before solving the granular problem: save previous particle's position and velocity (predictor step with gravity, first subproblem after N-S) */ foreach_dimension() { (*Unm1).x = (*U).x; /* Check if the domain is periodic, if yes shift the particle's position when the end of the domain is reached */ if (Period.x) { if ((*gci).center.x > (L0 + decal.x)) { (*gci).center.x -= L0; } if ((*gci).center.x < (0. + decal.x)) { (*gci).center.x += L0; } } (*gcinm1).center.x = (*gci).center.x; } } else { /* corrector step without gravity: fourth subproblem after the fictitious domain problem (correction step) */ foreach_dimension() { (*gci).center.x = (*gcinm1).center.x; } } /* integrating in time with rk4 */ for (gi = 1; gi <= miter; gi++) { Uold = (*U); Xold = (*gci).center; /* compute k1 */ compute_Fontact (&k1, &p[k], &Xold, &Uold, gamman); foreach_dimension() { k1.x /= (fsf*M); k1.x += gravity_flag*(*gravity).x; Xtemp.x = Xold.x + 0.5*dtg*Uold.x; Utemp.x = Uold.x + 0.5*dtg*k1.x; } /* compute k2 */ compute_Fontact (&k2, &p[k], &Xtemp, &Utemp, gamman); foreach_dimension() { k2.x /= (fsf*M); k2.x += gravity_flag*(*gravity).x; Xtemp.x = Xold.x + 0.5*dtg*Uold.x + 0.25*sq(dtg)*k1.x; Utemp.x = Uold.x + 0.5*dtg*k2.x; } /* compute k3 */ compute_Fontact (&k3, &p[k], &Xtemp, &Utemp, gamman); foreach_dimension() { k3.x /= (fsf*M); k3.x += gravity_flag*(*gravity).x; Xtemp.x = Xold.x + dtg*Uold.x + 0.5*sq(dtg)*k2.x; Utemp.x = Uold.x + dtg*k3.x; } /* compute k4 */ compute_Fontact (&k4, &p[k], &Xtemp, &Utemp, gamman); foreach_dimension() { k4.x /= (fsf*M); k4.x += gravity_flag*(*gravity).x; } foreach_dimension () { (*gci).center.x = Xold.x + dtg*Uold.x + sq(dtg)*(k1.x + k2.x + k3.x)/6.; /* Check if the domain is periodic, if yes shift the particle's position when the end of the domain is reached */ if (Period.x) { if ((*gci).center.x > (L0 + decal.x)) { (*gci).center.x -= L0; } if ((*gci).center.x < (0. + decal.x)) { (*gci).center.x += L0; } } (*U).x = Uold.x + dtg*(k1.x + 2*k2.x + 2*k3.x + k4.x)/6.; } } if (k == 0) { if (gravity_flag) { fprintf (stderr,"Prediction: particle-0's velocity on thread %d is (%20.18f, %20.18f, %20.18f)\n",pid(), (*U).x, (*U).y, (*U).z); fprintf (stderr,"Prediction: particle-0's position on thread %d is (%20.18f, %20.18f, %20.18f)\n",pid(), (*gci).center.x, (*gci).center.y, (*gci).center.z); } else { fprintf(stderr,"Correction: particle-0's velocity on thread %d is (%20.18f, %20.18f, %20.18f)\n",pid(), (*U).x, (*U).y, (*U).z); fprintf(stderr,"Correction: particle-0's position on thread %d is (%20.18f, %20.18f, %20.18f)\n",pid(), (*gci).center.x, (*gci).center.y, (*gci).center.z); } } } } #endif #endif /** Diagnostic functions. */ void writer_headers (FILE * pdata, FILE * sl) { #if _MPI if(pid() == 0) #endif { /* Write header for particles positions and velocities */ #if dimension == 2 fprintf(pdata, "#t\t position.x\t position.y\t U.x\t U.y\t w.z\n"); fprintf(sl, "#t\t F.x\t F.y\t T.z\n"); #elif dimension == 3 fprintf (pdata, "#t\t position.x\t position.y\t position.z\t U.x\t U.y\t U.z\t w.x\t w.y\t w.z\n"); fprintf (sl, "#t\t F.x\t F.y\t F.z\t T.x\t T.y\t T.z\n"); #endif fflush(pdata); fflush(sl); } } void particle_data (particle * p, const double t, const int i, FILE ** pdata) { for (int k = 0; k < NPARTICLES; k++) { GeomParameter * GCi = &(p[k].g); #if DLM_Moving_particle #if TRANSLATION coord * U = &(p[k].U); #endif #if ROTATION coord * w = &(p[k].w); #endif #endif #if dimension == 2 if(pid() == 0) { fprintf(pdata[k], "%20.18f\t %20.18f\t %20.18f %20.18f\t %20.18f\t %20.18f \n", t, (*GCi).center.x, (*GCi).center.y #if DLM_Moving_particle #if TRANSLATION , (*U).x, (*U).y #elif TRANSLATION == 0 , 0., 0. #endif #if ROTATION , (*w).z #elif ROTATION == 0 , 0. #endif #elif DLM_Moving_particle == 0 , 0., 0. , 0. #endif ); fflush(pdata[k]); } #elif dimension == 3 if(pid() == 0) { fprintf(pdata[k], "%20.18f\t %20.18f\t %20.18f\t %20.18f\t %20.18f\t %20.18f %20.18f\t %20.18f\t %20.18f\t %20.18f\n", t, (*GCi).center.x, (*GCi).center.y, (*GCi).center.z #if DLM_Moving_particle #if TRANSLATION , (*U).x, (*U).y, (*U).z #elif TRANSLATION == 0 , 0., 0., 0. #endif #if ROTATION ,(*w).x, (*w).y, (*w).z #elif ROTATION == 0 , 0., 0., 0. #endif #elif DLM_Moving_particle == 0 , 0., 0., 0. , 0., 0., 0. #endif ); fflush(pdata[k]); } #endif } } coord sumLambda (particle * p, FILE ** sl, const double t, const double dt, scalar Flagfield, vector Lambda, vector Index_lambda, const double rho_f, vector pshift) { /* Compute hydrodynamic force & torque and write to a file */ /* Fh = _P + (rho_f/rho_s)MdU/dt */ /* Mh = _P + (rho_f/rho_s).( Idom/dt + om x (I.om)) */ /* The inertia tensor is: */ /* Ip[0] = Ixx */ /* Ip[1] = Iyy */ /* Ip[2] = Izz */ /* Ip[3] = Ixy */ /* Ip[4] = Ixz */ /* Ip[5] = Iyz */ coord lambdasumint; coord lambdasumboundary; coord lambdasum; coord crossLambdaSumInt; coord crossLambdaSumBoundary; coord crossLambdaSum; Cache * Interior[NPARTICLES]; Cache * Boundary[NPARTICLES]; GeomParameter * gci; SolidBodyBoundary * sbm; /* Loop over all particles */ for (int k = 0; k < NPARTICLES; k++) { foreach_dimension() { lambdasumint.x = 0.; lambdasumboundary.x = 0.; lambdasum.x = 0; crossLambdaSumInt.x = 0.; crossLambdaSumBoundary.x = 0.; crossLambdaSum.x = 0.; } #if dimension ==2 lambdasumint.z = 0.; lambdasumboundary.z = 0.; lambdasum.z = 0; crossLambdaSumInt.z = 0.; crossLambdaSumBoundary.z = 0.; crossLambdaSum.z = 0.; #endif /* Interior domain of the particle */ Interior[k] = &(p[k].Interior); /* Boundary domain of the particle */ Boundary[k] = &(p[k].reduced_domain); gci = &(p[k].g); sbm = &(p[k].s); coord lambdapos; #if DLM_Moving_particle double rho_s = p[k].rho_s; #if TRANSLATION double M = p[k].M; coord Unm1 = p[k].Unm1; coord U = p[k].U; #endif #if ROTATION coord wnm1 = p[k].wnm1; coord w = p[k].w; coord Iom, omIom, Idomdt; #endif #endif // Particule's interior multipliers foreach_cache((*Interior[k])) { if (Flagfield[] < 1 && k == Index_lambda.y[]) { /* Compute Fh = _P */ /* For moving particles the additional term + (rho_f/rho_s).M.dU/dt is added after the mpi_calls below */ foreach_dimension() { lambdasumint.x += Lambda.x[]; } /* Compute Mh = _P */ /* For moving particles, the additional term + (rho_f/rho_s).(I.dom/dt + om ^ (I.om)) is added after mpi calls below */ /* Modify temporerly the particle center position for periodic boundary condition */ foreach_dimension() (*gci).center.x += pshift.x[]; crossLambdaSumInt.x += (Lambda.z[]*(y - (*gci).center.y) - Lambda.y[]*(z - (*gci).center.z)); crossLambdaSumInt.y += (Lambda.x[]*(z - (*gci).center.z) - Lambda.z[]*(x - (*gci).center.x)); crossLambdaSumInt.z += (Lambda.y[]*(x - (*gci).center.x) - Lambda.x[]*(y - (*gci).center.y)); foreach_dimension() (*gci).center.x -= pshift.x[]; } }//end foreach_cache() // Particle's boundary multipliers foreach_cache((*Boundary[k])) { if ((Index_lambda.x[] > -1) && (p[k].pnum == (int)Index_lambda.y[])) { lambdapos.x = (*sbm).x[(int)Index_lambda.x[]]; lambdapos.y = (*sbm).y[(int)Index_lambda.x[]]; lambdapos.z = 0.; #if dimension == 3 lambdapos.z = (*sbm).z[(int)Index_lambda.x[]]; #endif /* Compute Fh = _P */ foreach_dimension() { lambdasumboundary.x += Lambda.x[]; } /* Modify temporerly the particle center position for periodic boundary condition */ foreach_dimension() (*gci).center.x += pshift.x[]; /* Compute Mh = _P */ crossLambdaSumBoundary.x += (Lambda.z[]*(lambdapos.y - (*gci).center.y) - Lambda.y[]*(lambdapos.z - (*gci).center.z)); crossLambdaSumBoundary.y += (Lambda.x[]*(lambdapos.z - (*gci).center.z) - Lambda.z[]*(lambdapos.x - (*gci).center.x)); crossLambdaSumBoundary.z += (Lambda.y[]*(lambdapos.x - (*gci).center.x) - Lambda.x[]*(lambdapos.y - (*gci).center.y)); foreach_dimension() (*gci).center.x -= pshift.x[]; } } #if _MPI foreach_dimension() { mpi_all_reduce (lambdasumint.x, MPI_DOUBLE, MPI_SUM); mpi_all_reduce (lambdasumboundary.x, MPI_DOUBLE, MPI_SUM); mpi_all_reduce (crossLambdaSumInt.x, MPI_DOUBLE, MPI_SUM); mpi_all_reduce (crossLambdaSumBoundary.x, MPI_DOUBLE, MPI_SUM); } #if dimension == 2 mpi_all_reduce (crossLambdaSumInt.z, MPI_DOUBLE, MPI_SUM); mpi_all_reduce (crossLambdaSumBoundary.z, MPI_DOUBLE, MPI_SUM); #endif #endif /* The Force term (rho_f/rho_s)MdU/dt is added here for translating particles */ #if DLM_Moving_particle #if TRANSLATION foreach_dimension() lambdasumint.x += (rho_f/rho_s)*M*(U.x - Unm1.x)/dt; #endif /* The Torque term (rho_f/rho_s).(Idom/dt + om ^ (I.om)) is added here for rotating particles */ #if ROTATION /* I.om term */ Iom.x = p[k].Ip[0]*w.x + p[k].Ip[3]*w.y + p[k].Ip[4]*w.z; Iom.y = p[k].Ip[3]*w.x + p[k].Ip[1]*w.y + p[k].Ip[5]*w.z; Iom.z = p[k].Ip[4]*w.x + p[k].Ip[5]*w.y + p[k].Ip[2]*w.z; /* om^(I.om) term */ omIom.x = w.y*Iom.z - w.z*Iom.y; omIom.y = w.z*Iom.x - w.x*Iom.z; omIom.z = w.x*Iom.y - w.y*Iom.x; /* Idom/dt term; */ Idomdt.x = p[k].Ip[0]*(w.x - wnm1.x) + p[k].Ip[3]*(w.y - wnm1.y) + p[k].Ip[4]*(w.z - wnm1.z); Idomdt.y = p[k].Ip[3]*(w.x - wnm1.x) + p[k].Ip[1]*(w.y - wnm1.y) + p[k].Ip[5]*(w.z - wnm1.z); Idomdt.z = p[k].Ip[4]*(w.x - wnm1.x) + p[k].Ip[5]*(w.y - wnm1.y) + p[k].Ip[2]*(w.z - wnm1.z); crossLambdaSumInt.x += (rho_f/rho_s)*((Idomdt.x)/dt + omIom.x); crossLambdaSumInt.y += (rho_f/rho_s)*((Idomdt.y)/dt + omIom.y); crossLambdaSumInt.z += (rho_f/rho_s)*((Idomdt.z)/dt + omIom.z); #endif #endif foreach_dimension() { lambdasum.x = lambdasumint.x + lambdasumboundary.x; crossLambdaSum.x = crossLambdaSumInt.x + crossLambdaSumBoundary.x; } #if dimension == 2 crossLambdaSum.z = crossLambdaSumInt.z + crossLambdaSumBoundary.z; if(pid() == 0) { fprintf(sl[k], "%f\t %20.18f\t %20.18f\t %20.18f\n", t, lambdasum.x, lambdasum.y, crossLambdaSum.z); fflush(sl[k]); } #elif dimension == 3 if(pid() == 0) { fprintf(sl[k], "%20.18f\t %20.18f\t %20.18f\t %20.18f\t %20.18f\t %20.18f\t %20.18f\n", t, lambdasum.x, lambdasum.y, lambdasum.z, crossLambdaSum.x, crossLambdaSum.y, crossLambdaSum.z); fflush(sl[k]); } #endif } return lambdasum; } void init_file_pointers (FILE ** p, FILE ** d, const size_t rflag) { char name[800]; char name2[800]; #if _MPI if (pid() == 0) #endif { for (int k = 0; k < NPARTICLES; k++) { sprintf (name, "particle-data-%d", k); sprintf (name2, "sum_lambda-%d", k); if (!rflag) { p[k] = fopen (name, "w"); d[k] = fopen (name2, "w"); /* if (k == (NPARTICLES - 1)) { */ /* printf ("pointer sum lambda= %p\n", (void *)d[k]); */ /* printf ("pointer particle data = %p\n", (void *)p[k]); */ /* } */ /* Write headers in these files */ writer_headers (p[k], d[k]); } else { p[k] = fopen (name, "a"); d[k] = fopen (name2, "a"); } } } } void vorticity_3D (const vector u, vector omega) { foreach() foreach_dimension() omega.x[] = ( (u.z[0,1,0] - u.z[0,-1,0]) - (u.y[0,0,1] - u.y[0,0,-1]) )/2.*Delta; boundary((scalar *){omega}); } void compute_flowrate_right (FILE * pf, const vector u, const int level) { coord velo = {0., 0., 0.}; foreach_dimension() velo.x = 0.; #if !adaptive Cache intDomain = {0}; Point lpoint; foreach() { if (x > (L0 - Delta)) { lpoint = locate(x, y, z); cache_append(&intDomain, lpoint, 0); } } cache_shrink(&intDomain); foreach_cache(intDomain){ foreach_dimension() velo.x += sq(Delta)*u.x[]; } free(intDomain.p); #endif #if adaptive double hh = L0/pow(2,level); double xi = 0., yj = 0., zval = 0.; int ii = 0, jj = 0; foreach_level(level) { zval = L0 - 0.5*Delta + Z0; } /* printf("zz = %f, h = %f\n", zval, hh); */ for (ii = 0; ii < pow(2,level); ii++) { xi = 0.5*hh + ii*hh + X0; for (jj = 0; jj < pow(2,level); jj++) { yj = 0.5*hh + jj*hh + Y0; coord uu = {0., 0., 0.}; foreach_dimension() { uu.x = interpolate (u.x, xi, yj, zval); if (uu.x < nodata) velo.x += sq(hh)*uu.x; } /* printf("thread %d, u = (%f,%f,%f)\n",pid(), interpolate(u.x, xi, yj, zz), interpolate(u.y, xi, yj, zz), interpolate(u.z, xi, yj, zz)); */ } } #endif /* printf("thread %d, velo = (%f,%f,%f)\n",pid(), velo.x, velo.y, velo.z); */ #if _MPI foreach_dimension () mpi_all_reduce (velo.x, MPI_DOUBLE, MPI_SUM); #endif fprintf(pf, "%f\t %20.18f\t %20.18f\t %20.18f\n", t, velo.x, velo.y, velo.z); fflush(pf); } void pressure_law (scalar pres, FILE * ppf, const coord hv, const double t) { double plaw = 0.; plaw = interpolate (pres, hv.x, hv.y, hv.z); #if _MPI MPI_Barrier(MPI_COMM_WORLD); if (plaw == nodata) plaw = 0.; MPI_Barrier(MPI_COMM_WORLD); if (plaw != nodata) mpi_all_reduce(plaw, MPI_DOUBLE, MPI_SUM); #endif if (pid() == 0) { fprintf(ppf, "%20.18f\t %20.18f\n", t, plaw); fflush(ppf); } } /** dump/resto particle's structure functions */ void dump_particle (particle * p) { FILE * fp; char * file = "dump_particle"; if (file && (fp = fopen (file, "w")) == NULL) { perror (file); exit (1); } assert (fp); fwrite(p, sizeof(*p), NPARTICLES, fp); fclose(fp); } bool restore_particle (particle * p) { FILE * fp; char * file = "dump_particle"; if (file && (fp = fopen (file, "r")) == NULL) return false; assert (fp); if (fread (p, sizeof(*p), NPARTICLES, fp) < 1) { fprintf (ferr, "restore_particle(): error reading particle data\n"); exit (1); } fclose(fp); return true; } void inverse3by3matrix(double oldMatrix[3][3], double ** inversedMatrix) { double block1 = oldMatrix[1][1]*oldMatrix[2][2] - oldMatrix[1][2]*oldMatrix[2][1]; double block2 = oldMatrix[1][2]*oldMatrix[2][0] - oldMatrix[1][0]*oldMatrix[2][2]; double block3 = oldMatrix[1][0]*oldMatrix[2][1] - oldMatrix[1][1]*oldMatrix[2][0]; double determinat = oldMatrix[0][0]*block1 + oldMatrix[0][1]*block2 + oldMatrix[0][2]*block3; /* m[i][j], *(*(m + i) + j) */ *(*(inversedMatrix + 0) + 0) = block1/determinat; *(*(inversedMatrix + 0) + 1) = (oldMatrix[0][2]*oldMatrix[2][1] - oldMatrix[0][1]*oldMatrix[2][2])/determinat; *(*(inversedMatrix + 0) + 2) = (oldMatrix[0][1]*oldMatrix[1][2] - oldMatrix[0][2]*oldMatrix[1][1])/determinat; *(*(inversedMatrix + 1) + 0) = block2 / determinat; *(*(inversedMatrix + 1) + 1) = (oldMatrix[0][0]*oldMatrix[2][2] - oldMatrix[0][2]*oldMatrix[2][0])/determinat; *(*(inversedMatrix + 1) + 2) = (oldMatrix[0][2]*oldMatrix[1][0] - oldMatrix[0][0]*oldMatrix[1][2])/determinat; *(*(inversedMatrix + 2) + 0) = block3/determinat; *(*(inversedMatrix + 2) + 1) = (oldMatrix[0][1]*oldMatrix[2][0] - oldMatrix[0][0]*oldMatrix[2][1])/determinat; *(*(inversedMatrix + 2) + 2) = (oldMatrix[0][0]*oldMatrix[1][1] - oldMatrix[0][1]*oldMatrix[1][0])/determinat ; } int totalcells() { int t = 0; foreach() { t++; } #if _MPI mpi_all_reduce (t, MPI_INTEGER, MPI_SUM); #endif return t; } int total_dlmfd_cells (particle * p, const int np) { int apts = 0; for (int k = 0; k < np; k++) { #if debugInterior == 0 apts += p[k].Interior.n; #endif #if debugBD == 0 apts += p[k].reduced_domain.n; #endif } #if _MPI mpi_all_reduce (apts, MPI_INTEGER, MPI_SUM); #endif return apts; } int total_dlmfd_multipliers (particle * p, const int np) { int apts = 0; for (int k = 0; k < np; k++) { #if debugInterior == 0 apts += p[k].Interior.n; #endif } #if _MPI mpi_all_reduce (apts, MPI_INTEGER, MPI_SUM); #endif for (int k = 0; k < np; k++) { #if debugBD == 0 SolidBodyBoundary * bla; bla = &(p[k].s); apts += bla->m; #endif } return apts; } double getDelta() { /* Return the smallest grid size */ /* ----------------------------- */ int ii = 0; double dd = 0.; foreach_level(depth()) { ii++; dd = Delta; if (ii > 0) break; } #if _MPI mpi_all_reduce (dd, MPI_DOUBLE, MPI_MIN); #endif return dd; }