Skip to content

GitLab

Commit df657eef authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

wip

parent e6e1505a
Pipeline #9371 failed with stages
in 1 minute
Hide whitespace changes
Inline Side-by-side
fluid/fluid_problem.c
View file @ df657eef
......@@ -35,8 +35,7 @@
#define LOCAL_MATRIX(i,j,a,b) local_matrix[((i)*n_fields*N_SF+N_SF*(a)+(j))*n_fields+(b)]
static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, double *dfddp, double *sol, double *dsol, const double *solold, const double c, const double *dc, double a, double dt, int iel, int ip);
void fluid_problem_interpolate_rho_and_nu(const FluidProblem *problem, double a, double *rho, double *mu) {
void fluid_problem_interpolate_rho_and_mu(const FluidProblem *problem, double a, double *rho, double *mu) {
if (problem->n_fluids == 1){
*rho = problem->rho[0];
*mu = problem->mu[0];
......@@ -133,12 +132,13 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
dc[k] += problem->porosity[el[i]]*dphi[i][k];
}
}
double f[D],dfdu,dfddp;
particle_force_f(problem,f,&dfdu,&dfddp,s,ds,sold,c,dc,a,dt,iel,ip);
//double f[D],dfdu,dfddp;
//particle_force_f(problem,f,&dfdu,&dfddp,s,ds,sold,c,dc,a,dt,iel,ip);
const double *f = problem->particle_force + ip*DIMENSION;
if (!local_vector)
continue;
double rho,nu;
fluid_problem_interpolate_rho_and_nu(problem,a, &rho,&nu);
double rho,mu;
fluid_problem_interpolate_rho_and_mu(problem,a, &rho,&mu);
double pspg = drag_in_stab?problem->taup[iel]/rho:0;
for (int iD = 0; iD < D; ++iD) {
for (int iphi = 0; iphi < N_SF; ++iphi){
......@@ -150,7 +150,7 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
}
}
}
if (!local_matrix)
/*if (!local_matrix)
continue;
for (int iD = 0; iD < D; ++iD){
for (int iphi = 0; iphi < N_SF; ++iphi){
......@@ -166,7 +166,7 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
}
}
}
}
}*/
}
free(s);
free(ds);
......@@ -287,7 +287,7 @@ void fluid_problem_compute_stab_parameters(FluidProblem *problem, double dt) {
double rho, mu;
a /= N_N;
a = fmin(1.,fmax(0.,a));
fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
fluid_problem_interpolate_rho_and_mu(problem,a, &rho, &mu);
double nu = mu/rho;
problem->taup[iel] = 1/sqrt(
(problem->temporal?pow2(2/dt):0.)
......@@ -330,6 +330,7 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
}
}
#if 0
static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, double *dfddp, double *sol, double *dsol, const double *solold, const double c, const double *dc, double a, double dt, int iel, int ip) {
const double *contact = &problem->particle_contact[ip*D];
double u[D], uold[D], dp[D];
......@@ -339,7 +340,7 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
dp[iD] = dsol[P*D+iD];
}
double mu, rho;
fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
fluid_problem_interpolate_rho_and_mu(problem,a, &rho, &mu);
double Due[D];
const double *up = &problem->particle_velocity[ip*D];
for (int j = 0; j < D; ++j) {
......@@ -373,6 +374,7 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
*dfdu = gammastar/c;
*dfddp = gammastar*dt/mass*vol;//-vol;
}
#endif
static void fluid_problem_f(const FluidProblem *problem, const double *sol, double *dsol, const double *solold, const double *dsolold, const double* mesh_velocity, const double c, const double *dc, const double *bf, const double cold, const double rho, const double mu, double dt, int iel, double *f0, double *f1, double *f00, double *f10, double *f01, double *f11) {
double p = sol[P];
......@@ -563,7 +565,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
}
}
double rho, mu;
fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
fluid_problem_interpolate_rho_and_mu(problem,a, &rho, &mu);
const double jw = LQW[iqp]*detbnd;
if (wbnd->type != BND_WALL){
for (int i = 0; i < D; ++i){
......@@ -848,7 +850,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
}
double mu;
double rho;
fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
fluid_problem_interpolate_rho_and_mu(problem,a, &rho, &mu);
const double jw = QW[iqp]*detj;
problem->kinetic_energy += rho*sqrt(s[U]*s[U]+s[U+1]*s[U+1])*sqrt(s[U]*s[U]+s[U+1]*s[U+1])/2*jw;
......@@ -1019,7 +1021,7 @@ static void fluid_problem_surface_tension(FluidProblem *problem, const double *s
}
double mu;
double rho;
fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);const double jw = QW[iqp]*det;
fluid_problem_interpolate_rho_and_mu(problem,a, &rho, &mu);const double jw = QW[iqp]*det;
for (int iphi = 0; iphi < N_SF; ++iphi) {
for (int id = 0; id < D; ++id) {
//local_vector[(U+id)*N_N+iphi] += c*sigma*jw*(divnda[id]-divext_da[id])*phi[iphi];
......@@ -1755,7 +1757,7 @@ void fluid_problem_set_coordinates(FluidProblem *problem, double *x) {
compute_porosity(problem->mesh, problem->node_volume, problem->porosity, problem->n_particles, problem->particle_position, problem->particle_volume, problem->particle_element_id, problem->particle_uvw, NULL);
}
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, double *contact) {
void fluid_problem_set_particles(FluidProblem *problem, int n, double *volume, double *position) {
if (problem->n_particles != n) {
problem->n_particles = n;
free(problem->particle_uvw);
......@@ -1776,12 +1778,9 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
problem->particle_contact = (double*)malloc(sizeof(double)*n*D);
}
for (int i = 0; i < n; ++i) {
problem->particle_mass[i] = mass[i];
problem->particle_volume[i] = volume[i];
for (int k = 0; k < D; ++k) {
problem->particle_position[i*D+k] = position[i*D+k];
problem->particle_velocity[i*D+k] = velocity[i*D+k];
problem->particle_contact[i*D+k] = contact[i*D+k];
}
problem->particle_element_id[i]=-1;
}
......@@ -2054,3 +2053,72 @@ void fluid_problem_reset_boundary_force(FluidProblem *problem) {
for (int i=0;i<D*problem->mesh->n_boundaries;i++)
problem->boundary_force[i] = 0;
}
double *fluid_problem_particle_force(FluidProblem *p) {
return p->particle_force;
}
typedef struct{
double v[DIMENSION];
double dp[DIMENSION];
double mu;
double rho;
double porosity;
}FluidData;
void fluid_problem_get_particle_data(const FluidProblem *problem, FluidData *data)
{
const Mesh *mesh = problem->mesh;
const double *porosity = problem->porosity;
const double *solution = problem->solution;
int n_fields = fluid_problem_n_fields(problem);
double *s = malloc(sizeof(double)*n_fields);
double *ds = malloc(sizeof(double)*n_fields*D);
for (int ip=0; ip < problem->n_particles; ++ip) {
int iel = problem->particle_element_id[ip];
if(iel < 0){
for (int d = 0; d < D; ++d)
problem->particle_force[ip*D+d] = 0;
continue;
}
double *xi = problem->particle_uvw + D*ip;
double phi[N_SF];
shape_functions(xi,phi);
const int *el = &mesh->elements[iel*N_N];
double dxidx[D][D], dphi[N_SF][D];
double detj = mesh_dxidx(mesh, iel, dxidx);
grad_shape_functions(dxidx, dphi);
for (int i = 0; i < n_fields; ++i) {
s[i] = 0;
for (int j = 0; j < D; ++j) {
ds[i*D+j] = 0;
}
}
double c = 0;
double a = 0;
double cold = 0;
double dc[D] = {0};
for (int i = 0; i < N_SF; ++i) {
c += problem->porosity[el[i]]*phi[i];
cold += problem->oldporosity[el[i]]*phi[i];
if (problem->n_fluids == 2) {
a += problem->concentration[iel*N_N+i]*phi[i];
}
for (int j = 0; j < n_fields; ++j) {
double dof = solution[el[i]*n_fields+j];
s[j] += phi[i]*dof;
for (int k = 0; k < D; ++k)
ds[j*D+k] += dphi[i][k]*dof;
}
}
FluidData *datap = data + ip;
fluid_problem_interpolate_rho_and_mu(problem,a, &datap->rho,&datap->mu);
for (int d = 0; d < D; ++d) {
datap->v[d] = s[U+d];
datap->dp[d] = ds[P*D+d];
}
datap->porosity = c;
}
free(s);
free(ds);
}
fluid/fluid_problem.h
View file @ df657eef
......@@ -105,7 +105,7 @@ struct FluidProblem {
// complete force on the particle (including gravity)
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, double *contact);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *volume, double *position);
void fluid_problem_move_particles(FluidProblem *problem, int n, double *position, double *velocity, double *contact);
void fluid_problem_set_bulk_force(FluidProblem *problem, double *force);
void fluid_problem_set_concentration_cg(FluidProblem *problem, double *concentration);
......
python/fluid.py
View file @ df657eef
......@@ -165,8 +165,9 @@ class FluidProblem :
raise ValueError("dim should be 2 or 3.")
self._lib.fluid_problem_new.restype = c_void_p
n_fluids = np.require(rho,"float","C").reshape([-1]).shape[0]
self._n_fluids = n_fluids
self._dim = dim
self._fluiddatatype = np.dtype([('v', np.float64,(self._dim,)), ('dp', np.float64,(self._dim,)), ('mu',np.float64),('rho', np.float64) , ('porosity', np.float64)],align=True)
self._n_fluids = n_fluids
self._fp = c_void_p(self._lib.fluid_problem_new(_np2c(g), n_fluids, _np2c(mu), _np2c(rho), c_double(sigma), c_double(coeff_stab), c_double(volume_drag), c_double(quadratic_drag), c_int(drag_in_stab),c_double(drag_coefficient_factor),c_double(ip_factor),c_int(temporal), c_int(advection)))
if self._fp == None :
raise NameError("Cannot create fluid problem.")
......@@ -425,27 +426,6 @@ class FluidProblem :
self.old_porosity()[:] = data["old_porosity"]
self.sys = None
def compute_node_force(self, dt) :
"""Computes the forces to apply on each grain resulting from the fluid motion.
Keyword arguments:
dt -- Computation time step
"""
forces = np.ndarray([self.n_particles,self._dim],"d",order="C")
self._lib.fluid_problem_compute_node_particle_force(self._fp, c_double(dt), c_void_p(forces.ctypes.data))
return forces
def compute_node_torque(self, dt) :
"""Computes the angular drags to apply on each grain resulting from the fluid motion.
Only in 2D
Keyword arguments:
dt -- Computation time step
"""
torques = np.ndarray([self.n_particles,1],"d",order="C")
self._lib.fluid_problem_compute_node_particle_torque(self._fp, c_double(dt), c_void_p(torques.ctypes.data))
return torques
def implicit_euler(self, dt, check_residual_norm=-1, reduced_gravity=0, stab_param=0.) :
"""Solves the fluid equations.
......@@ -509,28 +489,24 @@ class FluidProblem :
for i in range(nsub) :
self._lib.fluid_problem_advance_concentration(self._fp,c_double(dt/nsub))
def set_particles(self, mass, volume, position, velocity, contact):
def set_particles(self, volume, position):
"""Set location of the grains in the mesh and compute the porosity in each cell.
Keyword arguments:
mass -- List of particles mass
volume -- List of particles volume
position -- List of particles centre positions
velocity -- List of particles velocity
contact -- List of particles contact resultant forces
"""
self.n_particles = mass.shape[0]
self._lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),_np2c(mass),_np2c(volume),_np2c(position),_np2c(velocity),_np2c(contact))
self.n_particles = volume.shape[0]
self._lib.fluid_problem_set_particles(self._fp,c_int(volume.shape[0]),_np2c(volume),_np2c(position))
def move_particles(self, position, velocity, contact):
def move_particles(self, position, velocity):
"""Set location of the grains in the mesh and compute the porosity in each cell.
Keyword arguments:
position -- List of particles centre positions
velocity -- List of particles velocity
contact -- List of particles contact resultant forces
"""
self._lib.fluid_problem_move_particles(self._fp,c_int(velocity.shape[0]),_np2c(position),_np2c(velocity),_np2c(contact))
self._lib.fluid_problem_move_particles(self._fp,c_int(velocity.shape[0]),_np2c(position),_np2c(velocity))
def _get_matrix(self, f_name, nrow, ncol,typec=c_double) :
f = getattr(self._lib,"fluid_problem_"+f_name)
......@@ -644,6 +620,15 @@ class FluidProblem :
def particle_volume(self) :
"""Gives access to the fluid field value at the mesh nodes."""
return self._get_matrix("particle_volume", self.n_particles, 1)
def particle_force(self) :
"""Gives access to the force apply by each particle on the fluid """
return self._get_matrix("particle_force", self.n_particles, 2)
def get_data_at_particles(self):
data = np.ndarray(self.n_particles, dtype=self._fluiddatatype)
self._lib.fluid_problem_get_particle_data(self._fp, c_void_p(data.ctypes.data))
return data
def bulk_force(self) :
"""Gives access to the volume force at fluid nodes."""
......
python/scontact.py
View file @ df657eef
......@@ -109,6 +109,7 @@ class ParticleProblem :
self._periodicTriangletype = np.dtype([('entity_id', np.int64),('p', np.float64,(3,dim))])
self._bodytype = np.dtype([('inverseI', np.float64),('inverseM', np.float64), ('theta', np.float64),('omega',np.float64),('position',np.float64,dim),('velocity',np.float64,dim)], align=True)
self._contacttype = np.dtype([('o', np.uint64,(2,)), ('type', np.int32), ('space',np.float64),('r', np.float64,(2,self._dim)) , ('reaction', np.float64,(self._dim,))],align=True)
self._fluiddatatype = np.dtype([('v', np.float64,(self._dim,)), ('dp', np.float64,(self._dim,)), ('mu',np.float64),('rho', np.float64) , ('porosity', np.float64)],align=True)
self.material2id = {}
self.tag2id = {}
self.id2tag = []
......@@ -163,6 +164,13 @@ class ParticleProblem :
else:
return self._get_array("body", self._bodytype)["omega"]
def body_theta(self) :
d = 1 if self._dim==2 else 3
if self._dim == 2:
return self._get_array("body", self._bodytype)["theta"][:,None]
else:
return self._get_array("body", self._bodytype)["theta"]
def save_state(self) :
self._saved_velocity = np.copy(self.body_velocity())
self._saved_position = np.copy(self.body_position())
......@@ -170,13 +178,21 @@ class ParticleProblem :
if self.dim() == 3 :
self._saved_triangles = np.copy(self.triangles())
self._saved_omega = np.copy(self.body_omega())
self._saved_theta = np.copy(self._bodies()["theta"])
self._saved_theta = np.copy(self.body_theta())
def set_fluid_data(self, fluiddata):
self._lib.particle_problem_set_fluid_data(self._p, _np2c(fluiddata,self._fluiddatatype))
def get_fluid_force(self):
f = np.zeros((self.n_particles(), self._dim))
self._lib.particle_problem_get_fluid_force(self._p, _np2c(f))
return f
def restore_state(self) :
self.body_velocity()[:] = self._saved_velocity
self.body_position()[:] = self._saved_position
self.body_omega()[:] = self._saved_omega
self._bodies()["theta"] = _self.saved_theta
self.body_theta()[:] = _self.saved_theta
self.segments()[:] = self._saved_segments
if self.dim() == 3 :
self.triangles()[:] = self._saved_triangles
......
python/time_integration.py
View file @ df657eef
......@@ -89,19 +89,13 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
external_particles_forces = np.zeros_like(particles.velocity())
# Copy the fluid solution of the previous time step before computing the prediction
sf0 = np.copy(fluid.solution())
# Copy the solid solution of the previous time step before computing the prediction
x0 = np.copy(particles.position())
cp0 = np.copy(particles.contact_forces())
seg0 = np.copy(particles.segments())
disk0 = np.copy(particles.disks())
if particles.dim() == 3 :
tri0 = np.copy(particles.triangles())
# Predictor
#
# Compute the fluid solution and keep a copy of this solution and the forces applied by the fluid on the grains
fluid.particle_force()[:] = particles.get_fluid_force()
fluid.implicit_euler(dt)
sf1 = np.copy(fluid.solution())
f0 = np.copy(fluid.compute_node_force(dt))+external_particles_forces
particles.set_fluid_data(fluid.get_data_at_particles())
_advance_particles(particles,f0,dt,min_nsub,contact_tol,max_nsub=max_nsub)
# Keep same contact forces and positions while giving to the fluid the predicted solid velocities to compute the correction
......@@ -111,6 +105,7 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
#
# Compute the fluid solution from the previous time-step one using the predicted solid velocities in the fluid-grain interaction force
fluid.solution()[:,:] = sf0
fluid.particle_force()[:] = particles.get_fluid_force()
fluid.implicit_euler(dt)
# Fluid solution is a weighted average of the predicted one and the corrected one.
# The alpha parametre is the weight
......@@ -122,9 +117,10 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
fluid.advance_concentration(dt)
# Reset solid velocities
particles.restore_state()
particles.set_fluid_data(fluid.get_data_at_particles())
_advance_particles(particles,(alpha*f0+(1-alpha)*f1),dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter,max_nsub=max_nsub)
# Give to the fluid the solid information
fluid.move_particles(particles.position(), particles.velocity(), particles.contact_forces())
fluid.move_particles(particles.position(), particles.velocity())
def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, fixed_grains=False, after_sub_iter=None,max_nsub=None) :
......
scontact/scontact.c
View file @ df657eef
......@@ -40,6 +40,7 @@ typedef struct _Triangle Triangle;
typedef struct _PeriodicEntity PeriodicEntity;
typedef struct _PeriodicSegment PeriodicSegment;
typedef struct _PeriodicTriangle PeriodicTriangle;
typedef struct _FluidData FluidData;
struct _ParticleProblem{
Particle *particles;
......@@ -55,6 +56,7 @@ struct _ParticleProblem{
PeriodicTriangle *periodicTriangles;
#endif
int use_queue;
FluidData *fluid_data;
};
struct _Particle{
......@@ -95,10 +97,44 @@ struct _Triangle {
double v[DIMENSION];
};
struct _FluidData{
double v[DIMENSION];
double dp[DIMENSION];
double mu;
double rho;
double porosity;
};
#if DIMENSION == 2
static double particle_drag_coeff(double u[2], double r, const FluidData *fluid)
{
double d = 2*r;
double normu = hypot(u[0],u[1]);
//Reynolds/|u_p-u_f|
double Re_O_u = d*fluid->porosity/fluid->mu*fluid->rho;
double Cd_u = 0.63*sqrt(normu)+4.8/sqrt(Re_O_u);
Cd_u = Cd_u*Cd_u;
double f = pow(fluid->porosity,-(1.8-0.65*exp(-(1.5-log(Re_O_u*normu))*(1.5-log(Re_O_u*normu))/2.)));
return f*Cd_u*fluid->rho/2*d;
}
#else
static double particle_drag_coeff(double u[3], double r, const FluidData *fluid)
{
double d = 2*r;
double normu = sqrt(u[0]*u[0]+u[1]*u[1]+u[2]*u[2]);
//Reynolds/|u_p-u_f|
double Re_O_u = d*c/fluid->mu*fluid->rho;
double Cd_u = 0.63*sqrt(normu)+4.8/sqrt(Re_O_u);
Cd_u = Cd_u*Cd_u;
double f = pow(fluid->porosity,-1.8);
return f*Cd_u*fluid->rho/2*(M_PI*r*r);
}
#endif
#if DIMENSION == 3
typedef enum {PARTICLE_PARTICLE=0, PARTICLE_DISK, PARTICLE_SEGMENT, PARTICLE_TRIANGLE} ContactType;
typedef enum {PARTICLE_PARTICLE=0, PARTICLE_FLUID, PARTICLE_SEGMENT, PARTICLE_TRIANGLE} ContactType;
#else
typedef enum {PARTICLE_PARTICLE=0, PARTICLE_DISK, PARTICLE_SEGMENT} ContactType;
typedef enum {PARTICLE_PARTICLE=0, PARTICLE_FLUID, PARTICLE_SEGMENT} ContactType;
#endif
struct _Contact {
......@@ -469,6 +505,20 @@ static void contact_get_omega_pointers(const Contact *c, ParticleProblem *p, dou
}
}
static int contact_init_fluid_particle(Contact *c, ParticleProblem *p, size_t particle_id) {
const Particle *particle = &p->particles[particle_id];
const Body *body = &p->bodies[particle->body];
if (body->invertI == 0 && body->invertM == 0)
return 0;
c->o0 = -1;
c->o1 = particle_id;
for (int i = 0; i < DIMENSION; ++i) {
c->reaction[i] = 0;
}
c->type = PARTICLE_FLUID;
return 1;
}
static void contact_get_velocity_pointers(const Contact *c, ParticleProblem *p, double **v0, double **v1) {
*v1 = p->bodies[p->particles[c->o1].body].velocity;
......@@ -795,6 +845,17 @@ static void particle_problem_gen_contacts(ParticleProblem *p, const double alert
contact_tree_gen_boundary_contact(tree, PARTICLE_TRIANGLE, i, old_contacts, &iold);
}
#endif
if (vectorSize(p->particles) == vectorSize(p->fluid_data)) {
for (size_t i = 0; i < vectorSize(p->particles); ++i){
const FluidData *data = &p->fluid_data[i];
const Body *body = &p->bodies[p->particles[i].body];
if (body->invertI == 0 && body->invertM == 0) continue;
Contact *c = vectorPush(&p->contacts);
int r = contact_init_fluid_particle(c, p, i);
if(!r)
vectorPop(p->contacts);
}
}
vectorFree(old_contacts);
vectorFree(found);
contact_tree_free(tree);
......@@ -843,10 +904,66 @@ static void body_advance(Body *b, double dt) {
}
}
void fluid_particle_force(double v[DIMENSION], double r, const FluidData *fluid, double f[DIMENSION]) {
#if DIMENSION == 2
double vol = M_PI*r*r;
double u[2] = {v[0]-fluid->v[0], v[1]-fluid->v[1]};
#else
double vol = 4./3*M_PI*r*r*r;
double u[3] = {v[0]-fluid->v[0], v[1]-fluid->v[1], v[2]-fluid->v[2]};
#endif
double coeff = particle_drag_coeff(u, r, fluid);
for (int d = 0; d < DIMENSION; ++d)
f[d] = -fluid->dp[d]*vol -coeff*u[d];
}
static int contact_solve_fluid(ParticleProblem *p, Contact *c, double dt, double tol)
{
double r1[DIMENSION];
const Particle *p1 = &p->particles[c->o1];
Body *body1 = &p->bodies[p1->body];
body_advance(body1, dt);
double x1[DIMENSION];
body_get_point(body1, p1->x, x1);
r1[0] = x1[0]-body1->position[0];
r1[1] = x1[1]-body1->position[1];
// compute the drag
double oldv[2]= {
body1->velocity[0]-body1->omega*c->r1[1], body1->velocity[1]+body1->omega*c->r1[0]
};
double coeff = particle_drag_coeff(oldv, p1->r, &p->fluid_data[c->o1]);
body_advance(body1, -dt);
// un-apply the old contact
double rold[DIMENSION];
for (int d = 0; d < DIMENSION; ++d) {
rold[d] = c->reaction[d];
}
for (int i = 0; i < DIMENSION; ++i) {
body1->velocity[i] += body1->invertM*(-dt*c->reaction[i]);
}
body1->omega += body1->invertI*(c->r1[0]*(-dt*c->reaction[1])-c->r1[1]*(-dt*c->reaction[0]));
// apply the new contact
for (int d=0; d<DIMENSION; ++d) {
c->r1[d] = r1[d];
}
for (int i = 0; i < DIMENSION; ++i) {
body1->velocity[i] += body1->invertM*(dt*c->reaction[i]);
}
body1->omega += body1->invertI*(c->r1[0]*(dt*c->reaction[1])-c->r1[1]*(dt*c->reaction[0]));
double newv[2]= {
body1->velocity[0]-body1->omega*c->r1[1], body1->velocity[1]+body1->omega*c->r1[0]
};
double dv[2] = {newv[0]-oldv[0], newv[1]-oldv[1]};
return (dv[0]*dv[0]+dv[1]*dv[1])*dt*dt < tol*tol;
}
static int contact_solve(ParticleProblem *p, Contact *c, double dt, double tol) {
// position prediction with the velocity taking into account the previous contact
// then initialize the contact with those positions and revert to the position in n
if (c->type == PARTICLE_FLUID)
return contact_solve_fluid(p, c, dt, tol);
double D;
double nx, ny;
int contact_avoided = 0;
......@@ -873,10 +990,10 @@ static int contact_solve(ParticleProblem *p, Contact *c, double dt, double tol)
nnorm = sqrt(nx*nx+ny*ny);//hypot(nx,ny);
nx /= nnorm;
ny /= nnorm;
r0[0] = x0[0]-body0->position[0] + nx*p0->r;
r1[0] = x1[0]-body1->position[0] - nx*p1->r;
r0[1] = x0[1]-body0->position[1] + ny*p0->r;
r1[1] = x1[1]-body1->position[1] - ny*p1->r;
r0[0] = x0[0]-body0->position[0] + nx*p0->r;
r1[0] = x1[0]-body1->position[0] - nx*p1->r;
r0[1] = x0[1]-body0->position[1] + ny*p0->r;
r1[1] = x1[1]-body1->position[1] - ny*p1->r;
D = nnorm - (p0->r + p1->r);
}
}
......@@ -1464,6 +1581,7 @@ ParticleProblem *particle_problem_new() {
p->segments = NULL;
p->periodicEntities = NULL;
p->periodicSegments = NULL;
p->fluid_data = NULL;
p->mu = NULL;
p->n_material = 0;
#if DIMENSION == 3
......@@ -1492,3 +1610,40 @@ void particle_problem_allocate_material(ParticleProblem *p, size_t n) {
vectorFree(oldmu);
}
void particle_problem_set_fluid_data(ParticleProblem *p, FluidData *data) {
vectorClear(p->fluid_data);
vectorPushN(&p->fluid_data, vectorSize(p->particles));
memcpy(p->fluid_data, data, vectorSize(p->particles)*sizeof(FluidData));
printf("set fluid data with size %lu\n", vectorSize(p->fluid_data));
}
void particle_problem_get_fluid_force(ParticleProblem *p, double *force) {
for (size_t i = 0; i < vectorSize(p->particles); ++i) {
const Particle *particle = &p->particles[i];
if (particle->r == 0) {
for (int d = 0; d < DIMENSION; ++d) {
force[i*DIMENSION+d] = 0;
}
continue;