Commit 07aee60f authored by Nathan Coppin's avatar Nathan Coppin
Browse files

Fluid damping particle omega

parent 9c9e2a7d
Pipeline #4890 passed with stage
in 26 seconds
......@@ -163,6 +163,38 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
}
}
#if D==2
void fluid_problem_compute_node_particle_torque(FluidProblem *problem, double dt, double *particle_torque) {
for (int i = 0; i < problem->n_particles; ++i) {
particle_torque[i] = problem->particle_torque[i];
}
}
static void particle_torque_f(FluidProblem *problem, double *f0, double *sol, double *dsol, const double *solold, const double c, const double *dc, const double cold, double dt, int iel, int ip, int in_derivative, int reduced_gravity) {
double a, mu, rho;
if (problem->n_fluids == 1){
rho = problem->rho[0];
mu = problem->mu[0];
}
else{
a = sol[A];
a = fmin(1.,fmax(0.,a));
rho = problem->rho[0]*a + problem->rho[1]*(1-a);
mu = problem->mu[0]*a + problem->mu[1]*(1-a);
}
double omega_res = problem->particle_omega[ip]+0.5*( dsol[2] - dsol[1] );
double r = fabs(cbrt((3*problem->particle_volume[ip]/(4*M_PI))));
double reynolds = r*r*fabs(omega_res)*rho/mu;
if(reynolds<11.0){
problem->particle_torque[ip] = -6.0*problem->particle_volume[ip]*mu*(omega_res)*(1+reynolds*reynolds/1200);
}
else{
problem->particle_torque[ip] = -problem->particle_volume[ip]*rho*r*r*(omega_res*omega_res)/(sqrt(reynolds));
}
}
#endif
static void particle_force_f(FluidProblem *problem, double *f0, double *sol, double *dsol, const double *solold, const double c, const double *dc, const double cold, double dt, int iel, int ip, int in_derivative, int reduced_gravity) {
double p = sol[P];
......@@ -284,6 +316,9 @@ static void node_force_volume(FluidProblem *problem, const double *solution_old,
dc[k] += problem->porosity[el[i]]*dphi[i][k];
}
}
#if D==2
particle_torque_f(problem,f0,s,ds,sold,c,dc,cold,dt,iel,ip,0,reduced_gravity);
#endif
particle_force_f(problem,f0,s,ds,sold,c,dc,cold,dt,iel,ip,0,reduced_gravity);
for (int ifield = 0; ifield < n_fields; ++ifield) {
for (int iphi = 0; iphi < N_SF; ++iphi){
......@@ -293,11 +328,17 @@ static void node_force_volume(FluidProblem *problem, const double *solution_old,
for (int jfield = 0; jfield < n_fields; ++jfield) {
double store = s[jfield];
s[jfield] += deps;
#if D==2
particle_torque_f(problem,g0,s,ds,sold,c,dc,cold,dt,iel,ip,1,reduced_gravity);
#endif
particle_force_f(problem,g0,s,ds,sold,c,dc,cold,dt,iel,ip,1,reduced_gravity);
s[jfield] = store;
for (int jd =0; jd < D; ++jd){
store = ds[jfield*D+jd];
ds[jfield*D+jd] += deps;
#if D==2
particle_torque_f(problem,h0+jd*n_fields,s,ds,sold,c,dc,cold,dt,iel,ip,1,reduced_gravity);
#endif
particle_force_f(problem,h0+jd*n_fields,s,ds,sold,c,dc,cold,dt,iel,ip,1,reduced_gravity);
ds[jfield*D+jd] = store;
}
......@@ -1129,6 +1170,10 @@ FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho,
problem->particle_force = NULL;
problem->particle_volume = NULL;
problem->particle_velocity = NULL;
#if DIMENSION==2
problem->particle_omega = NULL;
problem->particle_torque = NULL;
#endif
return problem;
}
......@@ -1146,8 +1191,11 @@ void fluid_problem_free(FluidProblem *problem)
free(problem->particle_position);
free(problem->particle_mass);
free(problem->particle_force);
free(problem->particle_volume);
free(problem->particle_velocity);
#if DIMENSION==2
free(problem->particle_omega);
free(problem->particle_torque);
#endif
free(problem->taup);
free(problem->tauc);
free(problem->taua);
......@@ -1478,7 +1526,60 @@ void fluid_problem_after_import(FluidProblem *problem)
#endif
}
#if DIMENSION==2
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, double *omega, long int *elid, int reload)
{
struct timespec tic, toc;
int n_fields = fluid_problem_n_fields(problem);
if(problem->n_particles != n) {
problem->n_particles = n;
free(problem->particle_uvw);
free(problem->particle_element_id);
free(problem->particle_position);
free(problem->particle_mass);
free(problem->particle_force);
free(problem->particle_torque);
free(problem->particle_volume);
free(problem->particle_velocity);
free(problem->particle_omega);
problem->particle_uvw = (double*)malloc(sizeof(double)*D*n);
problem->particle_element_id = (int*)malloc(sizeof(int)*n);
for (int i = 0; i < n; ++i) {
problem->particle_element_id[i]=-1;
}
problem->particle_position = (double*)malloc(sizeof(double)*D*n);
problem->particle_mass = (double*)malloc(sizeof(double)*n);
problem->particle_force = (double*)malloc(sizeof(double)*n*D);
problem->particle_volume = (double*)malloc(sizeof(double)*n);
problem->particle_velocity = (double*)malloc(sizeof(double)*n*D);
problem->particle_omega = (double*)malloc(sizeof(double)*n);
problem->particle_torque = (double*)malloc(sizeof(double)*n);
}
if (elid) {
for (int i = 0; i < n; ++i) {
problem->particle_element_id[i]=elid[i];
}
}
mesh_tree_particle_to_mesh(problem->mesh_tree, n, position, problem->particle_element_id, problem->particle_uvw);
for (int i = 0; i < n; ++i) {
problem->particle_mass[i] = mass[i];
problem->particle_volume[i] = volume[i];
problem->particle_omega[i] = omega[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];
}
}
if (!reload){
for (int i=0; i<problem->mesh->n_nodes; ++i){
problem->oldporosity[i] = problem->porosity[i];
}
}
fluid_problem_compute_porosity(problem);
}
#else
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid, int reload)
{
......@@ -1526,7 +1627,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
}
fluid_problem_compute_porosity(problem);
}
#endif
int *fluid_problem_particle_element_id(FluidProblem *problem)
{
return problem->particle_element_id;
......
/*
/*
* MigFlow - Copyright (C) <2010-2018>
* <Universite catholique de Louvain (UCL), Belgium
* Universite de Montpellier, France>
......@@ -76,6 +76,10 @@ struct FluidProblem {
double *particle_velocity;
double *particle_mass;
double *particle_force;
#if DIMENSION == 2
double *particle_omega;
double *particle_torque;
#endif
int n_fluids;
//stabilisation coefficients
double *taup;
......@@ -87,8 +91,15 @@ struct FluidProblem {
// complete force on the particle (including gravity)
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force);
#if DIMENSION==2
void fluid_problem_compute_node_particle_torque(FluidProblem *problem, double dt, double *particle_torque);
#endif
int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton_atol, double newton_rtol, int newton_max_it, int reduced_gravity, double stab_param);
#if DIMENSION==2
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, double *omega, long int *elid, int reload);
#else
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid, int reload);
#endif
void fluid_problem_free(FluidProblem *problem);
FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho, double coeffStab, const char *petsc_solver_type);
void fluid_problem_load_msh(FluidProblem *problem, const char *mesh_file_name);
......
......@@ -262,6 +262,17 @@ class FluidProblem :
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) :
"""Compute the angular drags to apply on each grain resulting from the fluid motion.
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, newton_atol=5e-6, newton_rtol=1e-5, newton_max_it=10, reduced_gravity=0, stab_param=0.) :
"""Solve the fluid equations.
......@@ -273,7 +284,7 @@ class FluidProblem :
"""
return self._lib.fluid_problem_implicit_euler(self._fp, c_double(dt), c_double(newton_atol), c_double(newton_rtol), c_int(newton_max_it), c_int(reduced_gravity), c_double(stab_param))
def set_particles(self, mass, volume, position, velocity, reload = False):
def set_particles(self, mass, volume, position, velocity, omega = None, reload = False):
"""Set location of the grains in the mesh and compute the porosity in each cell.
Keyword arguments:
......@@ -281,6 +292,7 @@ class FluidProblem :
volume -- list of particles volume
position -- list of particles centre positions
velocity -- list of particles velocity
omega -- Only in 2D -- list of particles angular velocity
reload -- optional arguments specifying if the simulation is reloaded or if it is the first iteration
"""
def np2c(a) :
......@@ -289,7 +301,11 @@ class FluidProblem :
r.tmp = tmp
return r
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),None,c_int(reload))
if omega is None:
self._lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),np2c(mass),np2c(volume),np2c(position),np2c(velocity),None,c_int(reload))
else:
self._lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),np2c(mass),np2c(volume),np2c(position),np2c(velocity),np2c(omega),None,c_int(reload))
def _get_matrix(self, f_name, nrow, ncol,typec=c_double) :
f = getattr(self._lib,"fluid_problem_"+f_name)
......
......@@ -172,15 +172,17 @@ class ParticleProblem :
arg_type = c_double*len(static)
self._lib.particleProblemAddFriction(self._p, arg_type(*static),len(static))
def iterate(self, dt, forces,tol=1e-8) :
def iterate(self, dt, forces, torques=None, tol=1e-8) :
"""Compute iteratively the set of velocities that respect the non-interpenetration constrain.
Keyword arguments:
dt -- numerical time step
forces -- list of vectors containing the forces applied on the particles
torques -- vector of torques on the particles
tol -- optional argument defining the interpenetration tolerance to stop the NLGS iterations of the NSCD
"""
self._get_matrix("ExternalForces",self._dim)[:] = forces
self._get_matrix("ExternalTorques",1)[:] = torques
self._lib.particleProblemIterate(self._p, c_double(np.min(self.r())), c_double(dt), c_double(tol), c_int(-1))
def write_vtk(self, odir, i, t) :
......@@ -346,3 +348,4 @@ class ParticleProblem :
(v1[0]+shift[0],v1[1]+shift[1],v1[2]+shift[2]),
(v2[0]+shift[0],v2[1]+shift[1],v2[2]+shift[2]),
stag)
......@@ -42,7 +42,7 @@ struct _ParticleProblem{
double friction_relaxation;
Particle *particles;
Contact *contacts;
double *position, *velocity, *externalForces, *omega, *mu;
double *position, *velocity, *externalForces, *externalTorques, *omega, *mu;
Disk *disks;
char **tagname;
int *diskTag, *segmentTag, *particleTag;
......@@ -347,40 +347,41 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
(p->velocity[c->o0*2+1]-p->velocity[c->o1*2+1])*t[1]+
p->omega[c->o0]*p->particles[c->o0].r +
p->omega[c->o1]*p->particles[c->o1].r;
//printf("vt avant : %.8f c->ct : %.8f \n",vt,c->ct);
*dp = fmax(0., vn - c->D/dt);
int index = p->particleTag[c->o0]*vectorSize(p->mu) + p->particleTag[c->o1];
double mu = p->mu[index];
if(p->frictionFlag && mu!=0){
if(vt>(7./2.)*mu*(*dp)){
*ct = p->friction_relaxation*mu*(*dp);//printf("Slip\n");
}
else if(vt<(-7./2.)*mu*(*dp)){
*ct = -p->friction_relaxation*mu*(*dp);//printf("Slip\n");
if(p->frictionFlag){
int index = p->particleTag[c->o0]*vectorSize(p->mu) + p->particleTag[c->o1];
double mu = p->mu[index];
if(mu>1e-12){
if(vt>(7./2.)*mu*(*dp)){
*ct = p->friction_relaxation*mu*(*dp);
}
else if(vt<(-7./2.)*mu*(*dp)){
*ct = -p->friction_relaxation*mu*(*dp);
}
else{
*ct = p->friction_relaxation*(2./7.)*vt;
}
*ct -= c->ct;
coordAdd(&p->velocity[c->o0 * DIMENSION], -*ct*c->a0, t);
coordAdd(&p->velocity[c->o1 * DIMENSION], *ct*c->a1, t);
p->omega[c->o0] -= c->In0*(*ct);
p->omega[c->o1] -= c->In1*(*ct);
}
else{
*ct = p->friction_relaxation*(2./7.)*vt;//printf("Stick\n");
*ct = 0;
}
*ct -= c->ct;
coordAdd(&p->velocity[c->o0 * DIMENSION], -*ct*c->a0, t);
coordAdd(&p->velocity[c->o1 * DIMENSION], *ct*c->a1, t);
p->omega[c->o0] -= c->In0*(*ct);
p->omega[c->o1] -= c->In1*(*ct);
}
else{
*ct = 0;
}
*dp -= c->dv;
//printf("ct : %.8f dp : %.8f\n",*ct,*dp);
coordAdd(&p->velocity[c->o0 * DIMENSION], -*dp*c->a0, c->n);
coordAdd(&p->velocity[c->o1 * DIMENSION], *dp*c->a1, c->n);
#endif
if(fabs(*dp) < tol/dt && fabs(*ct) <tol/dt){
//printf("b : %d b : %d ok\n",c->o0,c->o1);
if(fabs(*dp/p->friction_relaxation) < tol/dt && fabs(*ct/p->friction_relaxation) <tol/dt){
return 1;
}
else{
//printf("b : %d b : %d not ok\n",c->o0,c->o1);
return 0;
}
}
......@@ -409,29 +410,32 @@ static int contactParticleDiskSolve(Contact *c, ParticleProblem *p, double *ct,
p->velocity[c->o1 * DIMENSION+1]*t[1]+
p->omega[c->o1]*p->particles[c->o1].r+
p->disks[c->o0].vt;
//printf("vt avant : %.8f\n",vt);
int index = (p->particleTag[c->o1]+1)*vectorSize(p->mu) -1;
double mu = p->mu[index];
if(p->frictionFlag && mu!=0){
if(vt>(7./2.)*mu*(*dp)){
*ct = p->friction_relaxation*mu*(*dp);//printf("Slip\n");
}
else if(vt<(-7./2.)*mu*(*dp)){
*ct = -p->friction_relaxation*mu*(*dp);//printf("Slip\n");
if(p->frictionFlag){
int index = (p->particleTag[c->o1]+1)*vectorSize(p->mu) -1;
double mu = p->mu[index];
if(mu>1e-12){
if(vt>(7./2.)*mu*(*dp)){
*ct = p->friction_relaxation*mu*(*dp);
}
else if(vt<(-7./2.)*mu*(*dp)){
*ct = -p->friction_relaxation*mu*(*dp);
}
else{
*ct = p->friction_relaxation*(2./7.)*vt;
}
*ct -= c->ct;
coordAdd(&p->velocity[c->o1 * DIMENSION], -*ct*c->a1, t);
p->omega[c->o1] -= c->In1*(*ct);
}
else{
*ct = p->friction_relaxation*(2./7.)*vt;//printf("Stick\n");
*ct = 0;
}
*ct -= c->ct;
coordAdd(&p->velocity[c->o1 * DIMENSION], -*ct*c->a1, t);
p->omega[c->o1] -= c->In1*(*ct);
}
else{
*ct = 0;
}
*dp -=c->dv;
coordAdd(&p->velocity[c->o1 * DIMENSION], -*dp, c->n);
//printf("dis : %d ct : %f dp : %f\n",c->o0,*ct,*dp);
#else
......@@ -441,12 +445,10 @@ static int contactParticleDiskSolve(Contact *c, ParticleProblem *p, double *ct,
*ct = 0;
#endif
if(fabs(*dp) < tol/dt &&fabs(*ct) <tol/dt){
//printf("dis : %d b : %d ok\n",c->o0,c->o1);
if(fabs(*dp/p->friction_relaxation) < tol/dt &&fabs(*ct/p->friction_relaxation) <tol/dt){
return 1;
}
else{
//printf("dis : %d b : %d not ok\n",c->o0,c->o1);
return 0;
}
}
......@@ -472,28 +474,31 @@ static int contactParticleSegmentSolve(Contact *c, ParticleProblem *p, double *d
p->velocity[c->o1*DIMENSION +1]*t[1]+
p->omega[c->o1]*p->particles[c->o1].r+
p->segments[c->o0].vt;
//printf("ball : %d vt avant : %.8f c->ct : %.8f \n",c->o1,vt,c->ct);
int index = (p->particleTag[c->o1]+1)*vectorSize(p->mu) -1;
double mu = p->mu[index];
if(p->frictionFlag && mu!=0){
if(vt>(7./2.)*mu*(*dp)){
*ct = p->friction_relaxation*mu*(*dp);//printf("Slip\n");
}
else if(vt<(-7./2.)*mu*(*dp)){
*ct = -p->friction_relaxation*mu*(*dp);//printf("Slip\n");
if(p->frictionFlag){
int index = (p->particleTag[c->o1]+1)*vectorSize(p->mu) -1;
double mu = p->mu[index];
if(mu>1e-12){
if(vt>(7./2.)*mu*(*dp)){
*ct = p->friction_relaxation*mu*(*dp);
}
else if(vt<(-7./2.)*mu*(*dp)){
*ct = -p->friction_relaxation*mu*(*dp);
}
else{
*ct = p->friction_relaxation*(2./7.)*vt;
}
*ct -= c->ct;
coordAdd(&p->velocity[c->o1 * DIMENSION], -*ct*c->a1, t);
p->omega[c->o1] -= c->In1*(*ct);
}
else{
*ct = p->friction_relaxation*(2./7.)*vt;//printf("Stick\n");
*ct = 0;
}
*ct -= c->ct;
coordAdd(&p->velocity[c->o1 * DIMENSION], -*ct*c->a1, t);
p->omega[c->o1] -= c->In1*(*ct);
}
else{
*ct = 0;
}
*dp -=c->dv;
//printf("seg : %d ct : %f dp : %f\n",c->o0,*ct,*dp);
coordAdd(&p->velocity[c->o1 * DIMENSION], -*dp, c->n);
#else
......@@ -503,12 +508,10 @@ static int contactParticleSegmentSolve(Contact *c, ParticleProblem *p, double *d
coordAdd(&p->velocity[c->o1 * DIMENSION], -*dp, c->n);
#endif
if(fabs(*dp) < tol/dt &&fabs(*ct) <tol/dt){
//printf("seg : %d b : %d ok\n",c->o0,c->o1);
if(fabs(*dp/p->friction_relaxation) < tol/dt &&fabs(*ct/p->friction_relaxation) <tol/dt){
return 1;
}
else{
//printf("seg : %d b : %d not ok\n",c->o0,c->o1);
return 0;
}
}
......@@ -846,10 +849,10 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
int *found = NULL;
vectorPushN(&found, 100);
vectorClear(found);
// Particles
Contact *oldContacts = vectorDup(p->contacts), *cold;
vectorClear(p->contacts);
int ntot = 0;
// Particles
for (size_t i = 0; i < vectorSize(p->particles); ++i) {
particleBoundingBox(&p->particles[i], &p->position[i * DIMENSION], amin, amax);
addAlert(alert/2, amin, amax);
......@@ -869,12 +872,15 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
c->dv = cold->dv;
c->ct = cold->ct;
#if DIMENSION == 2
int index = p->particleTag[c->o0]*vectorSize(p->mu) + p->particleTag[c->o1];
if(p->frictionFlag && p->mu[index]!=0){
coordAdd(&p->velocity[c->o0 * DIMENSION], -c->ct*c->a0, t);
coordAdd(&p->velocity[c->o1 * DIMENSION], c->ct*c->a1, t);
p->omega[c->o0] -= c->In0*(c->ct);
p->omega[c->o1] -= c->In1*(c->ct);
if(p->frictionFlag){
int index = p->particleTag[c->o0]*vectorSize(p->mu) + p->particleTag[c->o1];
double mu = p->mu[index];
if(mu>1e-12){
coordAdd(&p->velocity[c->o0 * DIMENSION], -c->ct*c->a0, t);
coordAdd(&p->velocity[c->o1 * DIMENSION], c->ct*c->a1, t);
p->omega[c->o0] -= c->In0*(c->ct);
p->omega[c->o1] -= c->In1*(c->ct);
}
}
#endif
}
......@@ -900,10 +906,13 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
c->dv = cold->dv;
c->ct = cold->ct;
#if DIMENSION == 2
int index = (p->particleTag[c->o1]+1)*vectorSize(p->mu) -1;
if(p->frictionFlag && p->mu[index]!=0){
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*c->a1, t);
p->omega[c->o1] -= c->In1*(c->ct);
if(p->frictionFlag){
int index = (p->particleTag[c->o1]+1)*vectorSize(p->mu) -1;
double mu = p->mu[index];
if(mu<1e-12){
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*c->a1, t);
p->omega[c->o1] -= c->In1*(c->ct);
}
}
#endif
}
......@@ -927,15 +936,19 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
c->dv = cold->dv;
c->ct = cold->ct;
#if DIMENSION == 2
int index = (p->particleTag[c->o1]+1)*vectorSize(p->mu) -1;
if(p->frictionFlag && p->mu[index]!=0){
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*c->a1, t);
p->omega[c->o1] -= c->In1*(c->ct);
if(p->frictionFlag){
int index = (p->particleTag[c->o1]+1)*vectorSize(p->mu) -1;
double mu = p->mu[index];
if(mu>1e-12){
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*c->a1, t);
p->omega[c->o1] -= c->In1*(c->ct);
}
}
#endif
}
}
}
// Triangles
#if DIMENSION == 3
for (size_t i = 0; i < vectorSize(p->triangles); ++i) {
......@@ -1063,9 +1076,15 @@ static void _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, do
fifoPush(queue, i);
activeContact[i] = 1;
}
int iter, iterMax, trial;
iter, trial = 0;
double oldRelaxation = p->friction_relaxation;
iterMax = 250*vectorSize(p->particles);
for (size_t ic = fifoPop(queue); ic != FIFO_EMPTY; ic=fifoPop(queue)){
if(iter%iterMax==1 && iter>1 && trial<3 ){//Idée : recharger toute la queue car les contacts ont convergés sur un critère différent
p->friction_relaxation *= 0.95;
trial ++;
}
Contact *c = &p->contacts[ic];
int converged;
double dp, ct;
......@@ -1087,6 +1106,7 @@ static void _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, do
}
if (!converged) {
if(c->type == PARTICLE_PARTICLE) {
//if(iter>iterMax)printf("dp : %f\n",dp);
for (size_t j = contactByParticleP[c->o0]; j < contactByParticleP[c->o0+1]; ++j) {
if (!activeContact[contactByParticle[j]]){
activeContact[contactByParticle[j]] = 1;
......@@ -1104,8 +1124,9 @@ static void _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, do
c->dv += dp;
c->ct += ct;
activeContact[ic] = 0;
iter ++;
}
p->friction_relaxation = oldRelaxation;
fifoFree(queue);
free(activeContact);
free(contactByParticleP);
......@@ -1136,9 +1157,13 @@ void particleProblemSolve(ParticleProblem *p, double alert, double dt, double to
void particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit)
{
for (size_t j = 0; j < vectorSize(p->particles); ++j)
for (size_t j = 0; j < vectorSize(p->particles); ++j){
#if DIMENSION == 2
p->omega[j] += p->externalTorques[j]*dt/( (2./5.)*p->particles[j].m*p->particles[j].r*p->particles[j].r );
#endif
for (size_t i = 0; i < DIMENSION; ++i)
p->velocity[j * DIMENSION + i] += p->externalForces[j * DIMENSION + i] * dt / p->particles[j].m;
}
particleProblemSolve(p, alert, dt, tol, maxit);
for (size_t i = 0; i < vectorSize(p->position); ++i){
p->position[i] += p->velocity[i] * dt;
......@@ -1186,6 +1211,7 @@ void particleProblemAddParticle(ParticleProblem *p, const double x[DIMENSION], d
vectorPushN(&p->position, DIMENSION);
vectorPushN(&p->velocity, DIMENSION);
vectorPushN(&p->externalForces, DIMENSION);
vectorPushN(&p->externalTorques,1);
vectorPushN(&p->omega,1);